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 <iostream>
00030 #include <signal.h>
00031
00032 #include "G4RPGStrangeProduction.hh"
00033 #include "Randomize.hh"
00034 #include "G4SystemOfUnits.hh"
00035 #include "G4HadReentrentException.hh"
00036
00037 G4RPGStrangeProduction::G4RPGStrangeProduction()
00038 : G4RPGReaction() {}
00039
00040
00041 G4bool G4RPGStrangeProduction::
00042 ReactionStage(const G4HadProjectile* ,
00043 G4ReactionProduct& modifiedOriginal,
00044 G4bool& incidentHasChanged,
00045 const G4DynamicParticle* originalTarget,
00046 G4ReactionProduct& targetParticle,
00047 G4bool& targetHasChanged,
00048 const G4Nucleus& ,
00049 G4ReactionProduct& currentParticle,
00050 G4FastVector<G4ReactionProduct,256>& vec,
00051 G4int& vecLen,
00052 G4bool ,
00053 G4ReactionProduct& )
00054 {
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065 if( vecLen == 0 )return true;
00066
00067
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 );
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
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
00115
00116
00117 if (vecLen == 1) {
00118 i3 = 0;
00119 i4 = 1;
00120 } else {
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
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 )
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
00183 }
00184 else
00185 {
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
00216
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 )
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
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
00292
00293
00294
00295
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
00329 {
00330
00331
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
00409
00410
00411
00412
00413
00414
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 )
00424 {
00425 vecLen = std::max( 0, --i );
00426 G4int j;
00427 for(j=i; j<vecLen; j++) delete vec[j];
00428 break;
00429 }
00430 }
00431
00432 return true;
00433 }
00434
00435
00436