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
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044 #include "G4InelasticInteraction.hh"
00045 #include "G4PhysicalConstants.hh"
00046 #include "G4SystemOfUnits.hh"
00047 #include "Randomize.hh"
00048 #include "G4HadReentrentException.hh"
00049 #include "G4IsoResult.hh"
00050 #include "G4IsoParticleChange.hh"
00051
00052 G4IsoParticleChange* G4InelasticInteraction::theIsoResult = 0;
00053 G4IsoParticleChange* G4InelasticInteraction::theOldIsoResult = 0;
00054
00055 G4InelasticInteraction::G4InelasticInteraction(const G4String& name)
00056 : G4HadronicInteraction(name), isotopeProduction(false), cache(0.0)
00057 {}
00058
00059 G4InelasticInteraction::~G4InelasticInteraction()
00060 {}
00061
00062
00063 G4double
00064 G4InelasticInteraction::Pmltpc(G4int npos, G4int nneg, G4int nzero, G4int n,
00065 G4double b, G4double c)
00066 {
00067 const G4double expxu = 82.;
00068 const G4double expxl = -expxu;
00069 G4double npf = 0.0;
00070 G4double nmf = 0.0;
00071 G4double nzf = 0.0;
00072 G4int i;
00073 for (i = 2; i <= npos; i++) npf += std::log((G4double)i);
00074 for (i = 2; i <= nneg; i++) nmf += std::log((G4double)i);
00075 for (i = 2; i <= nzero; i++) nzf += std::log((G4double)i);
00076 G4double r = std::min(expxu, std::max(expxl,
00077 -(npos-nneg+nzero+b)*(npos-nneg+nzero+b)/(2*c*c*n*n)
00078 - npf - nmf - nzf ) );
00079 return std::exp(r);
00080 }
00081
00082
00083 G4bool G4InelasticInteraction::MarkLeadingStrangeParticle(
00084 const G4ReactionProduct& currentParticle,
00085 const G4ReactionProduct& targetParticle,
00086 G4ReactionProduct& leadParticle)
00087 {
00088
00089
00090
00091
00092
00093 G4bool lead = false;
00094 if ((currentParticle.GetMass() >= G4KaonPlus::KaonPlus()->GetPDGMass()) &&
00095 (currentParticle.GetDefinition() != G4Proton::Proton()) &&
00096 (currentParticle.GetDefinition() != G4Neutron::Neutron()) ) {
00097 lead = true;
00098 leadParticle = currentParticle;
00099
00100 } else if ((targetParticle.GetMass() >= G4KaonPlus::KaonPlus()->GetPDGMass()) &&
00101 (targetParticle.GetDefinition() != G4Proton::Proton()) &&
00102 (targetParticle.GetDefinition() != G4Neutron::Neutron() ) ) {
00103 lead = true;
00104 leadParticle = targetParticle;
00105 }
00106
00107 return lead;
00108 }
00109
00110
00111 void
00112 G4InelasticInteraction::SetUpPions(const G4int npos, const G4int nneg,
00113 const G4int nzero,
00114 G4FastVector<G4ReactionProduct,GHADLISTSIZE>& vec,
00115 G4int& vecLen)
00116 {
00117 if (npos + nneg + nzero == 0) return;
00118 G4int i;
00119 G4ReactionProduct* p;
00120
00121 for (i = 0; i < npos; ++i) {
00122 p = new G4ReactionProduct;
00123 p->SetDefinition( G4PionPlus::PionPlus() );
00124 (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
00125 vec.SetElement( vecLen++, p );
00126 }
00127 for (i = npos; i < npos+nneg; ++i) {
00128 p = new G4ReactionProduct;
00129 p->SetDefinition( G4PionMinus::PionMinus() );
00130 (G4UniformRand() < 0.5) ? p->SetSide(-1) : p->SetSide(1);
00131 vec.SetElement( vecLen++, p );
00132 }
00133 for (i = npos+nneg; i < npos+nneg+nzero; ++i) {
00134 p = new G4ReactionProduct;
00135 p->SetDefinition( G4PionZero::PionZero() );
00136 (G4UniformRand() < 0.5) ? p->SetSide(-1) : p->SetSide(1);
00137 vec.SetElement( vecLen++, p );
00138 }
00139 }
00140
00141
00142 void
00143 G4InelasticInteraction::GetNormalizationConstant(
00144 const G4double energy,
00145 G4double &n,
00146 G4double &anpn )
00147 {
00148 const G4double expxu = 82.;
00149 const G4double expxl = -expxu;
00150 const G4int numSec = 60;
00151
00152
00153
00154
00155 G4int iBegin = 1;
00156 G4double en = energy;
00157 if( energy < 0.0 )
00158 {
00159 iBegin = 2;
00160 en *= -1.0;
00161 }
00162
00163
00164
00165 G4double aleab = std::log(en/GeV);
00166 n = 3.62567 + aleab*(0.665843 + aleab*(0.336514 + aleab*(0.117712 + 0.0136912*aleab)));
00167 n -= 2.0;
00168
00169
00170
00171 anpn = 0.0;
00172 G4double test, temp;
00173 for( G4int i=iBegin; i<=numSec; ++i )
00174 {
00175 temp = pi*i/(2.0*n*n);
00176 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(i*i)/(n*n) ) ) );
00177 if( temp < 1.0 )
00178 {
00179 if( test >= 1.0e-10 )anpn += temp*test;
00180 }
00181 else
00182 anpn += temp*test;
00183 }
00184 }
00185
00186 void G4InelasticInteraction::CalculateMomenta(
00187 G4FastVector<G4ReactionProduct,GHADLISTSIZE>& vec,
00188 G4int& vecLen,
00189 const G4HadProjectile* originalIncident,
00190 const G4DynamicParticle* originalTarget,
00191 G4ReactionProduct& modifiedOriginal,
00192 G4Nucleus& targetNucleus,
00193 G4ReactionProduct& currentParticle,
00194 G4ReactionProduct& targetParticle,
00195 G4bool& incidentHasChanged,
00196 G4bool& targetHasChanged,
00197 G4bool quasiElastic)
00198 {
00199 cache = 0;
00200 what = originalIncident->Get4Momentum().vect();
00201
00202 theReactionDynamics.ProduceStrangeParticlePairs(vec, vecLen, modifiedOriginal,
00203 originalTarget, currentParticle,
00204 targetParticle, incidentHasChanged,
00205 targetHasChanged);
00206
00207 if (quasiElastic) {
00208 theReactionDynamics.TwoBody(vec, vecLen,
00209 modifiedOriginal, originalTarget,
00210 currentParticle, targetParticle,
00211 targetNucleus, targetHasChanged);
00212 return;
00213 }
00214 G4ReactionProduct leadingStrangeParticle;
00215 G4bool leadFlag = MarkLeadingStrangeParticle(currentParticle,
00216 targetParticle,
00217 leadingStrangeParticle);
00218
00219
00220 G4bool finishedGenXPt = false;
00221 G4bool annihilation = false;
00222 if( originalIncident->GetDefinition()->GetPDGEncoding() < 0 &&
00223 currentParticle.GetMass() == 0.0 && targetParticle.GetMass() == 0.0 )
00224 {
00225
00226 annihilation = true;
00227 G4double ekcor = 1.0;
00228 G4double ek = originalIncident->GetKineticEnergy();
00229 G4double ekOrg = ek;
00230
00231 const G4double tarmas = originalTarget->GetDefinition()->GetPDGMass();
00232 if( ek > 1.0*GeV )ekcor = 1./(ek/GeV);
00233 const G4double atomicWeight = targetNucleus.GetA_asInt();
00234 ek = 2*tarmas + ek*(1.+ekcor/atomicWeight);
00235 G4double tkin = targetNucleus.Cinema(ek);
00236 ek += tkin;
00237 ekOrg += tkin;
00238
00239
00240
00241
00242
00243 tkin = targetNucleus.AnnihilationEvaporationEffects(ek, ekOrg);
00244 ekOrg -= tkin;
00245 ekOrg = std::max( 0.0001*GeV, ekOrg );
00246 modifiedOriginal.SetKineticEnergy( ekOrg );
00247 G4double amas = originalIncident->GetDefinition()->GetPDGMass();
00248 G4double et = ekOrg + amas;
00249 G4double p = std::sqrt( std::abs(et*et-amas*amas) );
00250 G4double pp = modifiedOriginal.GetMomentum().mag();
00251 if( pp > 0.0 )
00252 {
00253 G4ThreeVector momentum = modifiedOriginal.GetMomentum();
00254 modifiedOriginal.SetMomentum( momentum * (p/pp) );
00255 }
00256 if( ekOrg <= 0.0001 )
00257 {
00258 modifiedOriginal.SetKineticEnergy( 0.0 );
00259 modifiedOriginal.SetMomentum( 0.0, 0.0, 0.0 );
00260 }
00261 }
00262 const G4double twsup[] = { 1.0, 0.7, 0.5, 0.3, 0.2, 0.1 };
00263 G4double rand1 = G4UniformRand();
00264 G4double rand2 = G4UniformRand();
00265
00266
00267 G4ReactionProduct saveCurrent = currentParticle;
00268 G4ReactionProduct saveTarget = targetParticle;
00269 std::vector<G4ReactionProduct> savevec;
00270 for (G4int i = 0; i < vecLen; i++) savevec.push_back(*vec[i]);
00271
00272 if (annihilation ||
00273 vecLen >= 6 ||
00274 ( modifiedOriginal.GetKineticEnergy()/GeV >= 1.0 &&
00275 ( ( (originalIncident->GetDefinition() == G4KaonPlus::KaonPlus() ||
00276 originalIncident->GetDefinition() == G4KaonMinus::KaonMinus() ||
00277 originalIncident->GetDefinition() == G4KaonZeroLong::KaonZeroLong() ||
00278 originalIncident->GetDefinition() == G4KaonZeroShort::KaonZeroShort() )
00279 &&
00280 rand1 < 0.5 )
00281 || rand2 > twsup[vecLen] ) ) )
00282
00283 finishedGenXPt =
00284 theReactionDynamics.GenerateXandPt( vec, vecLen,
00285 modifiedOriginal, originalIncident,
00286 currentParticle, targetParticle,
00287 originalTarget,
00288 targetNucleus, incidentHasChanged,
00289 targetHasChanged, leadFlag,
00290 leadingStrangeParticle );
00291 if( finishedGenXPt )
00292 {
00293 Rotate(vec, vecLen);
00294 return;
00295 }
00296
00297 G4bool finishedTwoClu = false;
00298 if( modifiedOriginal.GetTotalMomentum()/MeV < 1.0 )
00299 {
00300 for(G4int i=0; i<vecLen; i++) delete vec[i];
00301 vecLen = 0;
00302 }
00303 else
00304 {
00305
00306
00307
00308
00309 if (!finishedGenXPt && annihilation) {
00310 currentParticle = saveCurrent;
00311 targetParticle = saveTarget;
00312 for (G4int i = 0; i < vecLen; i++) delete vec[i];
00313 vecLen = 0;
00314 vec.Initialize( 0 );
00315 for (G4int i = 0; i < G4int(savevec.size()); i++) {
00316 G4ReactionProduct* p = new G4ReactionProduct;
00317 *p = savevec[i];
00318 vec.SetElement( vecLen++, p );
00319 }
00320 }
00321
00322 theReactionDynamics.SuppressChargedPions( vec, vecLen,
00323 modifiedOriginal, currentParticle,
00324 targetParticle, targetNucleus,
00325 incidentHasChanged, targetHasChanged );
00326 try
00327 {
00328 finishedTwoClu = theReactionDynamics.TwoCluster( vec, vecLen,
00329 modifiedOriginal, originalIncident,
00330 currentParticle, targetParticle,
00331 originalTarget,
00332 targetNucleus, incidentHasChanged,
00333 targetHasChanged, leadFlag,
00334 leadingStrangeParticle );
00335 }
00336 catch(G4HadReentrentException aC)
00337 {
00338 aC.Report(G4cout);
00339 throw G4HadReentrentException(__FILE__, __LINE__, "Failing to calculate momenta");
00340 }
00341 }
00342
00343 if (finishedTwoClu) {
00344 Rotate(vec, vecLen);
00345 return;
00346 }
00347
00348 theReactionDynamics.TwoBody(vec, vecLen, modifiedOriginal, originalTarget,
00349 currentParticle, targetParticle,
00350 targetNucleus, targetHasChanged);
00351 }
00352
00353
00354 void G4InelasticInteraction::Rotate(
00355 G4FastVector<G4ReactionProduct,GHADLISTSIZE>& vec, G4int& vecLen)
00356 {
00357 G4double rotation = 2.*pi*G4UniformRand();
00358 cache = rotation;
00359 G4int i;
00360 for (i = 0; i < vecLen; ++i) {
00361 G4ThreeVector momentum = vec[i]->GetMomentum();
00362 momentum = momentum.rotate(rotation, what);
00363 vec[i]->SetMomentum(momentum);
00364 }
00365 }
00366
00367
00368 void G4InelasticInteraction::SetUpChange(
00369 G4FastVector<G4ReactionProduct,GHADLISTSIZE>& vec,
00370 G4int& vecLen,
00371 G4ReactionProduct& currentParticle,
00372 G4ReactionProduct& targetParticle,
00373 G4bool& incidentHasChanged)
00374 {
00375 theParticleChange.Clear();
00376 G4ParticleDefinition* aKaonZL = G4KaonZeroLong::KaonZeroLong();
00377 G4ParticleDefinition* aKaonZS = G4KaonZeroShort::KaonZeroShort();
00378 G4int i;
00379 if (currentParticle.GetDefinition() == aKaonZL) {
00380 if (G4UniformRand() <= 0.5) {
00381 currentParticle.SetDefinition(aKaonZS);
00382 incidentHasChanged = true;
00383 }
00384 } else if (currentParticle.GetDefinition() == aKaonZS) {
00385 if (G4UniformRand() > 0.5) {
00386 currentParticle.SetDefinition(aKaonZL);
00387 incidentHasChanged = true;
00388 }
00389 }
00390
00391 if (targetParticle.GetDefinition() == aKaonZL) {
00392 if (G4UniformRand() <= 0.5) targetParticle.SetDefinition(aKaonZS);
00393 } else if (targetParticle.GetDefinition() == aKaonZS) {
00394 if (G4UniformRand() > 0.5 ) targetParticle.SetDefinition(aKaonZL);
00395 }
00396
00397 for (i = 0; i < vecLen; ++i) {
00398 if (vec[i]->GetDefinition() == aKaonZL) {
00399 if( G4UniformRand() <= 0.5) vec[i]->SetDefinition(aKaonZS);
00400 } else if (vec[i]->GetDefinition() == aKaonZS) {
00401 if (G4UniformRand() > 0.5) vec[i]->SetDefinition(aKaonZL);
00402 }
00403 }
00404
00405 if (incidentHasChanged) {
00406 G4DynamicParticle* p0 = new G4DynamicParticle;
00407 p0->SetDefinition( currentParticle.GetDefinition() );
00408 p0->SetMomentum( currentParticle.GetMomentum() );
00409 theParticleChange.AddSecondary( p0 );
00410 theParticleChange.SetStatusChange( stopAndKill );
00411 theParticleChange.SetEnergyChange( 0.0 );
00412 } else {
00413 G4double p = currentParticle.GetMomentum().mag()/MeV;
00414 G4ThreeVector pvec = currentParticle.GetMomentum();
00415 if (p > DBL_MIN)
00416 theParticleChange.SetMomentumChange(pvec.x()/p, pvec.y()/p, pvec.z()/p);
00417 else
00418 theParticleChange.SetMomentumChange(1.0, 0.0, 0.0);
00419
00420 G4double aE = currentParticle.GetKineticEnergy();
00421 if (std::fabs(aE) < .1*eV) aE = .1*eV;
00422 theParticleChange.SetEnergyChange(aE);
00423 }
00424
00425 if (targetParticle.GetMass() > 0.0) {
00426
00427 G4DynamicParticle* p1 = new G4DynamicParticle;
00428 p1->SetDefinition(targetParticle.GetDefinition() );
00429 G4ThreeVector momentum = targetParticle.GetMomentum();
00430 momentum = momentum.rotate(cache, what);
00431 p1->SetMomentum(momentum);
00432 theParticleChange.AddSecondary(p1);
00433 }
00434
00435 G4DynamicParticle* p;
00436 for (i = 0; i < vecLen; ++i) {
00437 p = new G4DynamicParticle(vec[i]->GetDefinition(), vec[i]->GetMomentum() );
00438 theParticleChange.AddSecondary( p );
00439 delete vec[i];
00440 }
00441 }
00442
00443
00444 G4IsoParticleChange* G4InelasticInteraction::GetIsotopeProductionInfo()
00445 {
00446 G4IsoParticleChange* anIsoResult = theIsoResult;
00447 if(theIsoResult) theOldIsoResult = theIsoResult;
00448 theIsoResult = 0;
00449 return anIsoResult;
00450 }
00451
00452
00453 void
00454 G4InelasticInteraction::DoIsotopeCounting(const G4HadProjectile* theProjectile,
00455 const G4Nucleus& aNucleus)
00456 {
00457 delete theOldIsoResult;
00458 theOldIsoResult = 0;
00459 delete theIsoResult;
00460 theIsoResult = new G4IsoParticleChange;
00461 G4IsoResult* anIsoResult = 0;
00462 G4int nModels = theProductionModels.size();
00463 if (nModels > 0) {
00464 for (G4int i = 0; i < nModels; i++) {
00465 anIsoResult = theProductionModels[i]->GetIsotope(theProjectile, aNucleus);
00466 if (anIsoResult) break;
00467 }
00468 if (!anIsoResult) anIsoResult = ExtractResidualNucleus(aNucleus);
00469 } else {
00470
00471 anIsoResult = ExtractResidualNucleus(aNucleus);
00472 }
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485 theIsoResult->SetIsotope(anIsoResult->GetIsotope());
00486 theIsoResult->SetProductionTime(theProjectile->GetGlobalTime() );
00487 theIsoResult->SetParentParticle(theProjectile->GetDefinition() );
00488 theIsoResult->SetMotherNucleus(anIsoResult->GetMotherNucleus());
00489 theIsoResult->SetProducer(this->GetModelName() );
00490
00491 delete anIsoResult;
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503 }
00504
00505
00506 G4IsoResult*
00507 G4InelasticInteraction::ExtractResidualNucleus(const G4Nucleus& aNucleus)
00508 {
00509 G4double A = aNucleus.GetA_asInt();
00510 G4double Z = aNucleus.GetZ_asInt();
00511 G4double bufferA = 0;
00512 G4double bufferZ = 0;
00513
00514
00515
00516 for (G4int i = 0; i < theParticleChange.GetNumberOfSecondaries(); ++i) {
00517 G4HadSecondary* aSecTrack = theParticleChange.GetSecondary(i);
00518 const G4ParticleDefinition* part = aSecTrack->GetParticle()->GetParticleDefinition();
00519 G4double Q = part->GetPDGCharge()/eplus;
00520 G4double BN = part->GetBaryonNumber();
00521 if (bufferA < BN) {
00522 bufferA = BN;
00523 bufferZ = Q;
00524 }
00525 Z -= Q;
00526 A -= BN;
00527 }
00528
00529
00530
00531 if (A < 0.1) {
00532 A = bufferA;
00533 Z = bufferZ;
00534 }
00535
00536
00537 std::ostringstream ost1;
00538 ost1 << Z << "_" << A;
00539 G4String biff = ost1.str();
00540 G4IsoResult* theResult = new G4IsoResult(biff, aNucleus);
00541
00542 return theResult;
00543 }
00544
00545 const std::pair<G4double, G4double> G4InelasticInteraction::GetFatalEnergyCheckLevels() const
00546 {
00547
00548 return std::pair<G4double, G4double>(5*perCent,250*GeV);
00549 }