#include <G4RPGStrangeProduction.hh>
Inheritance diagram for G4RPGStrangeProduction:
Public Member Functions | |
G4RPGStrangeProduction () | |
G4bool | ReactionStage (const G4HadProjectile *, G4ReactionProduct &, G4bool &, const G4DynamicParticle *, G4ReactionProduct &, G4bool &, const G4Nucleus &, G4ReactionProduct &, G4FastVector< G4ReactionProduct, 256 > &, G4int &, G4bool, G4ReactionProduct &) |
Definition at line 49 of file G4RPGStrangeProduction.hh.
G4RPGStrangeProduction::G4RPGStrangeProduction | ( | ) |
G4bool G4RPGStrangeProduction::ReactionStage | ( | const G4HadProjectile * | , | |
G4ReactionProduct & | , | |||
G4bool & | , | |||
const G4DynamicParticle * | , | |||
G4ReactionProduct & | , | |||
G4bool & | , | |||
const G4Nucleus & | , | |||
G4ReactionProduct & | , | |||
G4FastVector< G4ReactionProduct, 256 > & | , | |||
G4int & | , | |||
G4bool | , | |||
G4ReactionProduct & | ||||
) |
Reimplemented from G4RPGReaction.
Definition at line 42 of file G4RPGStrangeProduction.cc.
References G4AntiLambda::AntiLambda(), G4AntiNeutron::AntiNeutron(), G4AntiProton::AntiProton(), G4AntiSigmaMinus::AntiSigmaMinus(), G4AntiSigmaPlus::AntiSigmaPlus(), G4AntiSigmaZero::AntiSigmaZero(), G4UniformRand, G4DynamicParticle::GetDefinition(), G4ReactionProduct::GetDefinition(), G4ReactionProduct::GetMass(), G4ParticleDefinition::GetPDGMass(), G4ReactionProduct::GetTotalEnergy(), G4KaonMinus::KaonMinus(), G4KaonPlus::KaonPlus(), G4KaonZeroLong::KaonZeroLong(), G4KaonZeroShort::KaonZeroShort(), G4Lambda::Lambda(), G4Neutron::Neutron(), G4Proton::Proton(), G4ReactionProduct::SetDefinition(), G4ReactionProduct::SetDefinitionAndUpdateE(), G4FastVector< Type, N >::SetElement(), G4ReactionProduct::SetMayBeKilled(), G4ReactionProduct::SetSide(), G4SigmaMinus::SigmaMinus(), G4SigmaPlus::SigmaPlus(), and G4SigmaZero::SigmaZero().
00054 { 00055 // Derived from H. Fesefeldt's original FORTRAN code STPAIR 00056 // 00057 // Choose charge combinations K+ K-, K+ K0B, K0 K0B, K0 K-, 00058 // K+ Y0, K0 Y+, K0 Y- 00059 // For antibaryon induced reactions half of the cross sections KB YB 00060 // pairs are produced. Charge is not conserved, no experimental data available 00061 // for exclusive reactions, therefore some average behaviour assumed. 00062 // The ratio L/SIGMA is taken as 3:1 (from experimental low energy) 00063 // 00064 00065 if( vecLen == 0 )return true; 00066 // 00067 // the following protects against annihilation processes 00068 // 00069 if( currentParticle.GetMass() == 0.0 || targetParticle.GetMass() == 0.0 )return true; 00070 00071 const G4double etOriginal = modifiedOriginal.GetTotalEnergy()/GeV; 00072 const G4double mOriginal = modifiedOriginal.GetDefinition()->GetPDGMass()/GeV; 00073 G4double targetMass = originalTarget->GetDefinition()->GetPDGMass()/GeV; 00074 G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal + 00075 targetMass*targetMass + 00076 2.0*targetMass*etOriginal ); // GeV 00077 G4double currentMass = currentParticle.GetMass()/GeV; 00078 G4double availableEnergy = centerofmassEnergy-(targetMass+currentMass); 00079 if( availableEnergy <= 1.0 )return true; 00080 00081 G4ParticleDefinition *aProton = G4Proton::Proton(); 00082 G4ParticleDefinition *anAntiProton = G4AntiProton::AntiProton(); 00083 G4ParticleDefinition *aNeutron = G4Neutron::Neutron(); 00084 G4ParticleDefinition *anAntiNeutron = G4AntiNeutron::AntiNeutron(); 00085 G4ParticleDefinition *aSigmaMinus = G4SigmaMinus::SigmaMinus(); 00086 G4ParticleDefinition *aSigmaPlus = G4SigmaPlus::SigmaPlus(); 00087 G4ParticleDefinition *aSigmaZero = G4SigmaZero::SigmaZero(); 00088 G4ParticleDefinition *anAntiSigmaMinus = G4AntiSigmaMinus::AntiSigmaMinus(); 00089 G4ParticleDefinition *anAntiSigmaPlus = G4AntiSigmaPlus::AntiSigmaPlus(); 00090 G4ParticleDefinition *anAntiSigmaZero = G4AntiSigmaZero::AntiSigmaZero(); 00091 G4ParticleDefinition *aKaonMinus = G4KaonMinus::KaonMinus(); 00092 G4ParticleDefinition *aKaonPlus = G4KaonPlus::KaonPlus(); 00093 G4ParticleDefinition *aKaonZL = G4KaonZeroLong::KaonZeroLong(); 00094 G4ParticleDefinition *aKaonZS = G4KaonZeroShort::KaonZeroShort(); 00095 G4ParticleDefinition *aLambda = G4Lambda::Lambda(); 00096 G4ParticleDefinition *anAntiLambda = G4AntiLambda::AntiLambda(); 00097 00098 const G4double protonMass = aProton->GetPDGMass()/GeV; 00099 const G4double sigmaMinusMass = aSigmaMinus->GetPDGMass()/GeV; 00100 // 00101 // determine the center of mass energy bin 00102 // 00103 const G4double avrs[] = {3.,4.,5.,6.,7.,8.,9.,10.,20.,30.,40.,50.}; 00104 00105 G4int ibin, i3, i4; 00106 G4double avk, avy, avn, ran; 00107 G4int i = 1; 00108 while( (i<12) && (centerofmassEnergy>avrs[i]) )++i; 00109 if( i == 12 ) 00110 ibin = 11; 00111 else 00112 ibin = i; 00113 // 00114 // the fortran code chooses a random replacement of produced kaons 00115 // but does not take into account charge conservation 00116 // 00117 if (vecLen == 1) { // we know that vecLen > 0 00118 i3 = 0; 00119 i4 = 1; // note that we will be adding a new secondary particle in this case only 00120 } else { // otherwise 0 <= i3,i4 < vecLen 00121 G4double rnd = G4UniformRand(); 00122 while (rnd == 1.0) rnd = G4UniformRand(); 00123 i4 = i3 = G4int(vecLen*rnd); 00124 while(i3 == i4) { 00125 rnd = G4UniformRand(); 00126 while( rnd == 1.0 ) rnd = G4UniformRand(); 00127 i4 = G4int(vecLen*rnd); 00128 } 00129 } 00130 00131 // use linear interpolation or extrapolation by y=centerofmassEnergy*x+b 00132 // 00133 const G4double avkkb[] = { 0.0015, 0.005, 0.012, 0.0285, 0.0525, 0.075, 00134 0.0975, 0.123, 0.28, 0.398, 0.495, 0.573 }; 00135 const G4double avky[] = { 0.005, 0.03, 0.064, 0.095, 0.115, 0.13, 00136 0.145, 0.155, 0.20, 0.205, 0.210, 0.212 }; 00137 const G4double avnnb[] = { 0.00001, 0.0001, 0.0006, 0.0025, 0.01, 0.02, 00138 0.04, 0.05, 0.12, 0.15, 0.18, 0.20 }; 00139 00140 avk = (std::log(avkkb[ibin])-std::log(avkkb[ibin-1]))*(centerofmassEnergy-avrs[ibin-1]) 00141 /(avrs[ibin]-avrs[ibin-1]) + std::log(avkkb[ibin-1]); 00142 avk = std::exp(avk); 00143 00144 avy = (std::log(avky[ibin])-std::log(avky[ibin-1]))*(centerofmassEnergy-avrs[ibin-1]) 00145 /(avrs[ibin]-avrs[ibin-1]) + std::log(avky[ibin-1]); 00146 avy = std::exp(avy); 00147 00148 avn = (std::log(avnnb[ibin])-std::log(avnnb[ibin-1]))*(centerofmassEnergy-avrs[ibin-1]) 00149 /(avrs[ibin]-avrs[ibin-1]) + std::log(avnnb[ibin-1]); 00150 avn = std::exp(avn); 00151 00152 if( avk+avy+avn <= 0.0 )return true; 00153 00154 if( currentMass < protonMass )avy /= 2.0; 00155 if( targetMass < protonMass )avy = 0.0; 00156 avy += avk+avn; 00157 avk += avn; 00158 ran = G4UniformRand(); 00159 if( ran < avn ) 00160 { 00161 if( availableEnergy < 2.0 )return true; 00162 if( vecLen == 1 ) // add a new secondary 00163 { 00164 G4ReactionProduct *p1 = new G4ReactionProduct; 00165 if( G4UniformRand() < 0.5 ) 00166 { 00167 vec[0]->SetDefinition( aNeutron ); 00168 p1->SetDefinition( anAntiNeutron ); 00169 (G4UniformRand() < 0.5) ? p1->SetSide( -1 ) : p1->SetSide( 1 ); 00170 vec[0]->SetMayBeKilled(false); 00171 p1->SetMayBeKilled(false); 00172 } 00173 else 00174 { 00175 vec[0]->SetDefinition( aProton ); 00176 p1->SetDefinition( anAntiProton ); 00177 (G4UniformRand() < 0.5) ? p1->SetSide( -1 ) : p1->SetSide( 1 ); 00178 vec[0]->SetMayBeKilled(false); 00179 p1->SetMayBeKilled(false); 00180 } 00181 vec.SetElement( vecLen++, p1 ); 00182 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen); 00183 } 00184 else 00185 { // replace two secondaries 00186 if( G4UniformRand() < 0.5 ) 00187 { 00188 vec[i3]->SetDefinition( aNeutron ); 00189 vec[i4]->SetDefinition( anAntiNeutron ); 00190 vec[i3]->SetMayBeKilled(false); 00191 vec[i4]->SetMayBeKilled(false); 00192 } 00193 else 00194 { 00195 vec[i3]->SetDefinition( aProton ); 00196 vec[i4]->SetDefinition( anAntiProton ); 00197 vec[i3]->SetMayBeKilled(false); 00198 vec[i4]->SetMayBeKilled(false); 00199 } 00200 } 00201 } 00202 else if( ran < avk ) 00203 { 00204 if( availableEnergy < 1.0 )return true; 00205 00206 const G4double kkb[] = { 0.2500, 0.3750, 0.5000, 0.5625, 0.6250, 00207 0.6875, 0.7500, 0.8750, 1.000 }; 00208 const G4int ipakkb1[] = { 10, 10, 10, 11, 11, 12, 12, 11, 12 }; 00209 const G4int ipakkb2[] = { 13, 11, 12, 11, 12, 11, 12, 13, 13 }; 00210 ran = G4UniformRand(); 00211 i = 0; 00212 while( (i<9) && (ran>=kkb[i]) )++i; 00213 if( i == 9 )return true; 00214 // 00215 // ipakkb[] = { 10,13, 10,11, 10,12, 11,11, 11,12, 12,11, 12,12, 11,13, 12,13 }; 00216 // charge + - + 0 + 0 0 0 0 0 0 0 0 0 0 - 0 - 00217 // 00218 switch( ipakkb1[i] ) 00219 { 00220 case 10: 00221 vec[i3]->SetDefinition( aKaonPlus ); 00222 vec[i3]->SetMayBeKilled(false); 00223 break; 00224 case 11: 00225 vec[i3]->SetDefinition( aKaonZS ); 00226 vec[i3]->SetMayBeKilled(false); 00227 break; 00228 case 12: 00229 vec[i3]->SetDefinition( aKaonZL ); 00230 vec[i3]->SetMayBeKilled(false); 00231 break; 00232 } 00233 00234 if( vecLen == 1 ) // add a secondary 00235 { 00236 G4ReactionProduct *p1 = new G4ReactionProduct; 00237 switch( ipakkb2[i] ) 00238 { 00239 case 11: 00240 p1->SetDefinition( aKaonZS ); 00241 p1->SetMayBeKilled(false); 00242 break; 00243 case 12: 00244 p1->SetDefinition( aKaonZL ); 00245 p1->SetMayBeKilled(false); 00246 break; 00247 case 13: 00248 p1->SetDefinition( aKaonMinus ); 00249 p1->SetMayBeKilled(false); 00250 break; 00251 } 00252 (G4UniformRand() < 0.5) ? p1->SetSide( -1 ) : p1->SetSide( 1 ); 00253 vec.SetElement( vecLen++, p1 ); 00254 00255 } 00256 else // replace 00257 { 00258 switch( ipakkb2[i] ) 00259 { 00260 case 11: 00261 vec[i4]->SetDefinition( aKaonZS ); 00262 vec[i4]->SetMayBeKilled(false); 00263 break; 00264 case 12: 00265 vec[i4]->SetDefinition( aKaonZL ); 00266 vec[i4]->SetMayBeKilled(false); 00267 break; 00268 case 13: 00269 vec[i4]->SetDefinition( aKaonMinus ); 00270 vec[i4]->SetMayBeKilled(false); 00271 break; 00272 } 00273 } 00274 } 00275 else if( ran < avy ) 00276 { 00277 if( availableEnergy < 1.6 )return true; 00278 00279 const G4double ky[] = { 0.200, 0.300, 0.400, 0.550, 0.625, 0.700, 00280 0.800, 0.850, 0.900, 0.950, 0.975, 1.000 }; 00281 const G4int ipaky1[] = { 18, 18, 18, 20, 20, 20, 21, 21, 21, 22, 22, 22 }; 00282 const G4int ipaky2[] = { 10, 11, 12, 10, 11, 12, 10, 11, 12, 10, 11, 12 }; 00283 const G4int ipakyb1[] = { 19, 19, 19, 23, 23, 23, 24, 24, 24, 25, 25, 25 }; 00284 const G4int ipakyb2[] = { 13, 12, 11, 13, 12, 11, 13, 12, 11, 13, 12, 11 }; 00285 ran = G4UniformRand(); 00286 i = 0; 00287 while( (i<12) && (ran>ky[i]) )++i; 00288 if( i == 12 )return true; 00289 if( (currentMass<protonMass) || (G4UniformRand()<0.5) ) 00290 { 00291 // ipaky[] = { 18,10, 18,11, 18,12, 20,10, 20,11, 20,12, 00292 // 0 + 0 0 0 0 + + + 0 + 0 00293 // 00294 // 21,10, 21,11, 21,12, 22,10, 22,11, 22,12 } 00295 // 0 + 0 0 0 0 - + - 0 - 0 00296 switch( ipaky1[i] ) 00297 { 00298 case 18: 00299 targetParticle.SetDefinition( aLambda ); 00300 break; 00301 case 20: 00302 targetParticle.SetDefinition( aSigmaPlus ); 00303 break; 00304 case 21: 00305 targetParticle.SetDefinition( aSigmaZero ); 00306 break; 00307 case 22: 00308 targetParticle.SetDefinition( aSigmaMinus ); 00309 break; 00310 } 00311 targetHasChanged = true; 00312 switch( ipaky2[i] ) 00313 { 00314 case 10: 00315 vec[i3]->SetDefinition( aKaonPlus ); 00316 vec[i3]->SetMayBeKilled(false); 00317 break; 00318 case 11: 00319 vec[i3]->SetDefinition( aKaonZS ); 00320 vec[i3]->SetMayBeKilled(false); 00321 break; 00322 case 12: 00323 vec[i3]->SetDefinition( aKaonZL ); 00324 vec[i3]->SetMayBeKilled(false); 00325 break; 00326 } 00327 } 00328 else // (currentMass >= protonMass) && (G4UniformRand() >= 0.5) 00329 { 00330 // ipakyb[] = { 19,13, 19,12, 19,11, 23,13, 23,12, 23,11, 00331 // 24,13, 24,12, 24,11, 25,13, 25,12, 25,11 }; 00332 if( (currentParticle.GetDefinition() == anAntiProton) || 00333 (currentParticle.GetDefinition() == anAntiNeutron) || 00334 (currentParticle.GetDefinition() == anAntiLambda) || 00335 (currentMass > sigmaMinusMass) ) 00336 { 00337 switch( ipakyb1[i] ) 00338 { 00339 case 19: 00340 currentParticle.SetDefinitionAndUpdateE( anAntiLambda ); 00341 break; 00342 case 23: 00343 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaPlus ); 00344 break; 00345 case 24: 00346 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaZero ); 00347 break; 00348 case 25: 00349 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaMinus ); 00350 break; 00351 } 00352 incidentHasChanged = true; 00353 switch( ipakyb2[i] ) 00354 { 00355 case 11: 00356 vec[i3]->SetDefinition( aKaonZS ); 00357 vec[i3]->SetMayBeKilled(false); 00358 break; 00359 case 12: 00360 vec[i3]->SetDefinition( aKaonZL ); 00361 vec[i3]->SetMayBeKilled(false); 00362 break; 00363 case 13: 00364 vec[i3]->SetDefinition( aKaonMinus ); 00365 vec[i3]->SetMayBeKilled(false); 00366 break; 00367 } 00368 } 00369 else 00370 { 00371 switch( ipaky1[i] ) 00372 { 00373 case 18: 00374 currentParticle.SetDefinitionAndUpdateE( aLambda ); 00375 break; 00376 case 20: 00377 currentParticle.SetDefinitionAndUpdateE( aSigmaPlus ); 00378 break; 00379 case 21: 00380 currentParticle.SetDefinitionAndUpdateE( aSigmaZero ); 00381 break; 00382 case 22: 00383 currentParticle.SetDefinitionAndUpdateE( aSigmaMinus ); 00384 break; 00385 } 00386 incidentHasChanged = true; 00387 switch( ipaky2[i] ) 00388 { 00389 case 10: 00390 vec[i3]->SetDefinition( aKaonPlus ); 00391 vec[i3]->SetMayBeKilled(false); 00392 break; 00393 case 11: 00394 vec[i3]->SetDefinition( aKaonZS ); 00395 vec[i3]->SetMayBeKilled(false); 00396 break; 00397 case 12: 00398 vec[i3]->SetDefinition( aKaonZL ); 00399 vec[i3]->SetMayBeKilled(false); 00400 break; 00401 } 00402 } 00403 } 00404 } 00405 else return true; 00406 00407 // 00408 // check the available energy 00409 // if there is not enough energy for kkb/ky pair production 00410 // then reduce the number of secondary particles 00411 // NOTE: 00412 // the number of secondaries may have been changed 00413 // the incident and/or target particles may have changed 00414 // charge conservation is ignored (as well as strangness conservation) 00415 // 00416 currentMass = currentParticle.GetMass()/GeV; 00417 targetMass = targetParticle.GetMass()/GeV; 00418 00419 G4double energyCheck = centerofmassEnergy-(currentMass+targetMass); 00420 for( i=0; i<vecLen; ++i ) 00421 { 00422 energyCheck -= vec[i]->GetMass()/GeV; 00423 if( energyCheck < 0.0 ) // chop off the secondary List 00424 { 00425 vecLen = std::max( 0, --i ); // looks like a memory leak @@@@@@@@@@@@ 00426 G4int j; 00427 for(j=i; j<vecLen; j++) delete vec[j]; 00428 break; 00429 } 00430 } 00431 00432 return true; 00433 }