00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 #include "G4RPGXiMinusInelastic.hh"
00030 #include "G4PhysicalConstants.hh"
00031 #include "G4SystemOfUnits.hh"
00032 #include "Randomize.hh"
00033
00034 G4HadFinalState*
00035 G4RPGXiMinusInelastic::ApplyYourself( const G4HadProjectile &aTrack,
00036 G4Nucleus &targetNucleus )
00037 {
00038 const G4HadProjectile *originalIncident = &aTrack;
00039 if (originalIncident->GetKineticEnergy()<= 0.1*MeV)
00040 {
00041 theParticleChange.SetStatusChange(isAlive);
00042 theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy());
00043 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
00044 return &theParticleChange;
00045 }
00046
00047
00048
00049 G4DynamicParticle *originalTarget = targetNucleus.ReturnTargetParticle();
00050
00051 if( verboseLevel > 1 )
00052 {
00053 const G4Material *targetMaterial = aTrack.GetMaterial();
00054 G4cout << "G4RPGXiMinusInelastic::ApplyYourself called" << G4endl;
00055 G4cout << "kinetic energy = " << originalIncident->GetKineticEnergy()/MeV << "MeV, ";
00056 G4cout << "target material = " << targetMaterial->GetName() << ", ";
00057 G4cout << "target particle = " << originalTarget->GetDefinition()->GetParticleName()
00058 << G4endl;
00059 }
00060
00061
00062
00063
00064 G4double ek = originalIncident->GetKineticEnergy()/MeV;
00065 G4double amas = originalIncident->GetDefinition()->GetPDGMass()/MeV;
00066 G4ReactionProduct modifiedOriginal;
00067 modifiedOriginal = *originalIncident;
00068
00069 G4double tkin = targetNucleus.Cinema( ek );
00070 ek += tkin;
00071 modifiedOriginal.SetKineticEnergy( ek*MeV );
00072 G4double et = ek + amas;
00073 G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
00074 G4double pp = modifiedOriginal.GetMomentum().mag()/MeV;
00075 if( pp > 0.0 )
00076 {
00077 G4ThreeVector momentum = modifiedOriginal.GetMomentum();
00078 modifiedOriginal.SetMomentum( momentum * (p/pp) );
00079 }
00080
00081
00082
00083 tkin = targetNucleus.EvaporationEffects( ek );
00084 ek -= tkin;
00085 modifiedOriginal.SetKineticEnergy( ek*MeV );
00086 et = ek + amas;
00087 p = std::sqrt( std::abs((et-amas)*(et+amas)) );
00088 pp = modifiedOriginal.GetMomentum().mag()/MeV;
00089 if( pp > 0.0 )
00090 {
00091 G4ThreeVector momentum = modifiedOriginal.GetMomentum();
00092 modifiedOriginal.SetMomentum( momentum * (p/pp) );
00093 }
00094 G4ReactionProduct currentParticle = modifiedOriginal;
00095 G4ReactionProduct targetParticle;
00096 targetParticle = *originalTarget;
00097 currentParticle.SetSide( 1 );
00098 targetParticle.SetSide( -1 );
00099 G4bool incidentHasChanged = false;
00100 G4bool targetHasChanged = false;
00101 G4bool quasiElastic = false;
00102 G4FastVector<G4ReactionProduct,GHADLISTSIZE> vec;
00103 G4int vecLen = 0;
00104 vec.Initialize( 0 );
00105
00106 const G4double cutOff = 0.1;
00107 if( currentParticle.GetKineticEnergy()/MeV > cutOff )
00108 Cascade( vec, vecLen,
00109 originalIncident, currentParticle, targetParticle,
00110 incidentHasChanged, targetHasChanged, quasiElastic );
00111
00112 CalculateMomenta( vec, vecLen,
00113 originalIncident, originalTarget, modifiedOriginal,
00114 targetNucleus, currentParticle, targetParticle,
00115 incidentHasChanged, targetHasChanged, quasiElastic );
00116
00117 SetUpChange( vec, vecLen,
00118 currentParticle, targetParticle,
00119 incidentHasChanged );
00120
00121 delete originalTarget;
00122 return &theParticleChange;
00123 }
00124
00125
00126 void
00127 G4RPGXiMinusInelastic::Cascade(G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
00128 G4int& vecLen,
00129 const G4HadProjectile* originalIncident,
00130 G4ReactionProduct& currentParticle,
00131 G4ReactionProduct& targetParticle,
00132 G4bool& incidentHasChanged,
00133 G4bool& targetHasChanged,
00134 G4bool& quasiElastic)
00135 {
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146 const G4double mOriginal = originalIncident->GetDefinition()->GetPDGMass()/MeV;
00147 const G4double etOriginal = originalIncident->GetTotalEnergy()/MeV;
00148 const G4double targetMass = targetParticle.GetMass()/MeV;
00149 G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
00150 targetMass*targetMass +
00151 2.0*targetMass*etOriginal );
00152 G4double availableEnergy = centerofmassEnergy-(targetMass+mOriginal);
00153 if (availableEnergy <= G4PionPlus::PionPlus()->GetPDGMass()/MeV) {
00154 quasiElastic = true;
00155 return;
00156 }
00157 static G4bool first = true;
00158 const G4int numMul = 1200;
00159 const G4int numSec = 60;
00160 static G4double protmul[numMul], protnorm[numSec];
00161 static G4double neutmul[numMul], neutnorm[numSec];
00162
00163
00164 G4int counter, nt = 0, np = 0, nneg = 0, nz = 0;
00165 G4double test;
00166 const G4double c = 1.25;
00167 const G4double b[] = { 0.7, 0.7 };
00168 if (first) {
00169 first = false;
00170 G4int i;
00171 for (i = 0; i < numMul; ++i) protmul[i] = 0.0;
00172 for (i = 0; i < numSec; ++i) protnorm[i] = 0.0;
00173 counter = -1;
00174 for (np = 0; np < (numSec/3); ++np) {
00175 for (nneg = std::max(0,np-1); nneg <= (np+1); ++nneg) {
00176 for (nz = 0; nz < numSec/3; ++nz) {
00177 if (++counter < numMul) {
00178 nt = np + nneg + nz;
00179 if (nt > 0 && nt <= numSec) {
00180 protmul[counter] = Pmltpc(np,nneg,nz,nt,b[0],c);
00181 protnorm[nt-1] += protmul[counter];
00182 }
00183 }
00184 }
00185 }
00186 }
00187
00188 for( i=0; i<numMul; ++i )neutmul[i] = 0.0;
00189 for( i=0; i<numSec; ++i )neutnorm[i] = 0.0;
00190 counter = -1;
00191 for( np=0; np<numSec/3; ++np )
00192 {
00193 for( nneg=np; nneg<=(np+2); ++nneg )
00194 {
00195 for( nz=0; nz<numSec/3; ++nz )
00196 {
00197 if( ++counter < numMul )
00198 {
00199 nt = np+nneg+nz;
00200 if( nt>0 && nt<=numSec )
00201 {
00202 neutmul[counter] = Pmltpc(np,nneg,nz,nt,b[1],c);
00203 neutnorm[nt-1] += neutmul[counter];
00204 }
00205 }
00206 }
00207 }
00208 }
00209 for( i=0; i<numSec; ++i )
00210 {
00211 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
00212 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
00213 }
00214 }
00215
00216 const G4double expxu = 82.;
00217 const G4double expxl = -expxu;
00218 G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
00219 G4ParticleDefinition *aProton = G4Proton::Proton();
00220 G4ParticleDefinition *aKaonMinus = G4KaonMinus::KaonMinus();
00221 G4ParticleDefinition *aSigmaPlus = G4SigmaPlus::SigmaPlus();
00222 G4ParticleDefinition *aXiZero = G4XiZero::XiZero();
00223
00224
00225
00226 G4double n, anpn;
00227 GetNormalizationConstant( availableEnergy, n, anpn );
00228 G4double ran = G4UniformRand();
00229 G4double dum, excs = 0.0;
00230 if( targetParticle.GetDefinition() == aProton )
00231 {
00232 counter = -1;
00233 for( np=0; np<numSec/3 && ran>=excs; ++np )
00234 {
00235 for( nneg=std::max(0,np-1); nneg<=(np+1) && ran>=excs; ++nneg )
00236 {
00237 for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
00238 {
00239 if( ++counter < numMul )
00240 {
00241 nt = np+nneg+nz;
00242 if( nt>0 && nt<=numSec )
00243 {
00244 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00245 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
00246 if( std::fabs(dum) < 1.0 )
00247 {
00248 if( test >= 1.0e-10 )excs += dum*test;
00249 }
00250 else
00251 excs += dum*test;
00252 }
00253 }
00254 }
00255 }
00256 }
00257 if( ran >= excs )
00258 {
00259 quasiElastic = true;
00260 return;
00261 }
00262 np--; nneg--; nz--;
00263
00264
00265
00266
00267
00268
00269 if( np < nneg )
00270 {
00271 if( np+1 == nneg )
00272 {
00273 currentParticle.SetDefinitionAndUpdateE( aXiZero );
00274 incidentHasChanged = true;
00275 }
00276 else
00277 {
00278 currentParticle.SetDefinitionAndUpdateE( aSigmaPlus );
00279 incidentHasChanged = true;
00280
00281
00282
00283 vec.Initialize( 1 );
00284 G4ReactionProduct *p = new G4ReactionProduct;
00285 p->SetDefinition( aKaonMinus );
00286 (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
00287 vec.SetElement( vecLen++, p );
00288 --nneg;
00289 }
00290 }
00291 else if( np == nneg )
00292 {
00293 if( G4UniformRand() >= 0.5 )
00294 {
00295 currentParticle.SetDefinitionAndUpdateE( aXiZero );
00296 incidentHasChanged = true;
00297 targetParticle.SetDefinitionAndUpdateE( aNeutron );
00298 targetHasChanged = true;
00299 }
00300 }
00301 else
00302 {
00303 targetParticle.SetDefinitionAndUpdateE( aNeutron );
00304 targetHasChanged = true;
00305 }
00306 }
00307 else
00308 {
00309 counter = -1;
00310 for( np=0; np<numSec/3 && ran>=excs; ++np )
00311 {
00312 for( nneg=np; nneg<=(np+2) && ran>=excs; ++nneg )
00313 {
00314 for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
00315 {
00316 if( ++counter < numMul )
00317 {
00318 nt = np+nneg+nz;
00319 if( nt>0 && nt<=numSec )
00320 {
00321 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00322 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
00323 if( std::fabs(dum) < 1.0 )
00324 {
00325 if( test >= 1.0e-10 )excs += dum*test;
00326 }
00327 else
00328 excs += dum*test;
00329 }
00330 }
00331 }
00332 }
00333 }
00334 if( ran >= excs )
00335 {
00336 quasiElastic = true;
00337 return;
00338 }
00339 np--; nneg--; nz--;
00340 if( np+1 < nneg )
00341 {
00342 if( np+2 == nneg )
00343 {
00344 currentParticle.SetDefinitionAndUpdateE( aXiZero );
00345 incidentHasChanged = true;
00346 targetParticle.SetDefinitionAndUpdateE( aProton );
00347 targetHasChanged = true;
00348 }
00349 else
00350 {
00351 currentParticle.SetDefinitionAndUpdateE( aSigmaPlus );
00352 incidentHasChanged = true;
00353 targetParticle.SetDefinitionAndUpdateE( aProton );
00354 targetHasChanged = true;
00355
00356
00357
00358 vec.Initialize( 1 );
00359 G4ReactionProduct *p = new G4ReactionProduct;
00360 p->SetDefinition( aKaonMinus );
00361 (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
00362 vec.SetElement( vecLen++, p );
00363 --nneg;
00364 }
00365 }
00366 else if( np+1 == nneg )
00367 {
00368 if( G4UniformRand() < 0.5 )
00369 {
00370 currentParticle.SetDefinitionAndUpdateE( aXiZero );
00371 incidentHasChanged = true;
00372 }
00373 else
00374 {
00375 targetParticle.SetDefinitionAndUpdateE( aProton );
00376 targetHasChanged = true;
00377 }
00378 }
00379 }
00380 SetUpPions(np, nneg, nz, vec, vecLen);
00381 return;
00382 }
00383
00384
00385