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 "G4RPGInelastic.hh"
00030 #include "Randomize.hh"
00031 #include "G4PhysicalConstants.hh"
00032 #include "G4SystemOfUnits.hh"
00033 #include "G4HadReentrentException.hh"
00034 #include "G4RPGStrangeProduction.hh"
00035 #include "G4RPGTwoBody.hh"
00036
00037
00038 G4RPGInelastic::G4RPGInelastic(const G4String& modelName)
00039 : G4HadronicInteraction(modelName)
00040 {
00041 cache = 0.0;
00042 particleDef[0] = G4PionZero::PionZero();
00043 particleDef[1] = G4PionPlus::PionPlus();
00044 particleDef[2] = G4PionMinus::PionMinus();
00045 particleDef[3] = G4KaonPlus::KaonPlus();
00046 particleDef[4] = G4KaonMinus::KaonMinus();
00047 particleDef[5] = G4KaonZero::KaonZero();
00048 particleDef[6] = G4AntiKaonZero::AntiKaonZero();
00049 particleDef[7] = G4Proton::Proton();
00050 particleDef[8] = G4Neutron::Neutron();
00051 particleDef[9] = G4Lambda::Lambda();
00052 particleDef[10] = G4SigmaPlus::SigmaPlus();
00053 particleDef[11] = G4SigmaZero::SigmaZero();
00054 particleDef[12] = G4SigmaMinus::SigmaMinus();
00055 particleDef[13] = G4XiZero::XiZero();
00056 particleDef[14] = G4XiMinus::XiMinus();
00057 particleDef[15] = G4OmegaMinus::OmegaMinus();
00058 particleDef[16] = G4AntiProton::AntiProton();
00059 particleDef[17] = G4AntiNeutron::AntiNeutron();
00060
00061 G4cout << " **************************************************** " << G4endl;
00062 G4cout << " * The RPG model is currently under development and * " << G4endl;
00063 G4cout << " * should not be used. * " << G4endl;
00064 G4cout << " **************************************************** " << G4endl;
00065 }
00066
00067
00068 G4double G4RPGInelastic::Pmltpc(G4int np, G4int nneg, G4int nz,
00069 G4int n, G4double b, G4double c)
00070 {
00071 const G4double expxu = 82.;
00072 const G4double expxl = -expxu;
00073 G4double npf = 0.0;
00074 G4double nmf = 0.0;
00075 G4double nzf = 0.0;
00076 G4int i;
00077 for( i=2; i<=np; i++ )npf += std::log((double)i);
00078 for( i=2; i<=nneg; i++ )nmf += std::log((double)i);
00079 for( i=2; i<=nz; i++ )nzf += std::log((double)i);
00080 G4double r;
00081 r = std::min( expxu, std::max( expxl, -(np-nneg+nz+b)*(np-nneg+nz+b)/(2*c*c*n*n)-npf-nmf-nzf ) );
00082 return std::exp(r);
00083 }
00084
00085
00086 G4int G4RPGInelastic::Factorial( G4int n )
00087 {
00088 G4int j = std::min(n,10);
00089 G4int result = 1;
00090 if (j <= 1) return result;
00091 for (G4int i = 2; i <= j; ++i) result *= i;
00092 return result;
00093 }
00094
00095
00096 G4bool G4RPGInelastic::MarkLeadingStrangeParticle(
00097 const G4ReactionProduct ¤tParticle,
00098 const G4ReactionProduct &targetParticle,
00099 G4ReactionProduct &leadParticle )
00100 {
00101
00102
00103
00104
00105
00106
00107 G4bool lead = false;
00108 if( (currentParticle.GetMass() >= G4KaonPlus::KaonPlus()->GetPDGMass()) &&
00109 (currentParticle.GetDefinition() != G4Proton::Proton()) &&
00110 (currentParticle.GetDefinition() != G4Neutron::Neutron()) )
00111 {
00112 lead = true;
00113 leadParticle = currentParticle;
00114 }
00115 else if( (targetParticle.GetMass() >= G4KaonPlus::KaonPlus()->GetPDGMass()) &&
00116 (targetParticle.GetDefinition() != G4Proton::Proton()) &&
00117 (targetParticle.GetDefinition() != G4Neutron::Neutron()) )
00118 {
00119 lead = true;
00120 leadParticle = targetParticle;
00121 }
00122 return lead;
00123 }
00124
00125
00126 void G4RPGInelastic::SetUpPions(const G4int np, const G4int nneg,
00127 const G4int nz,
00128 G4FastVector<G4ReactionProduct,256> &vec,
00129 G4int &vecLen)
00130 {
00131 if( np+nneg+nz == 0 )return;
00132 G4int i;
00133 G4ReactionProduct *p;
00134 for( i=0; i<np; ++i )
00135 {
00136 p = new G4ReactionProduct;
00137 p->SetDefinition( G4PionPlus::PionPlus() );
00138 (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
00139 vec.SetElement( vecLen++, p );
00140 }
00141 for( i=np; i<np+nneg; ++i )
00142 {
00143 p = new G4ReactionProduct;
00144 p->SetDefinition( G4PionMinus::PionMinus() );
00145 (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
00146 vec.SetElement( vecLen++, p );
00147 }
00148 for( i=np+nneg; i<np+nneg+nz; ++i )
00149 {
00150 p = new G4ReactionProduct;
00151 p->SetDefinition( G4PionZero::PionZero() );
00152 (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
00153 vec.SetElement( vecLen++, p );
00154 }
00155 }
00156
00157
00158 void G4RPGInelastic::GetNormalizationConstant(
00159 const G4double energy,
00160 G4double &n,
00161 G4double &anpn )
00162 {
00163 const G4double expxu = 82.;
00164 const G4double expxl = -expxu;
00165 const G4int numSec = 60;
00166
00167
00168
00169
00170 G4int iBegin = 1;
00171 G4double en = energy;
00172 if( energy < 0.0 )
00173 {
00174 iBegin = 2;
00175 en *= -1.0;
00176 }
00177
00178
00179
00180 G4double aleab = std::log(en/GeV);
00181 n = 3.62567 + aleab*(0.665843 + aleab*(0.336514 + aleab*(0.117712 + 0.0136912*aleab)));
00182 n -= 2.0;
00183
00184
00185
00186 anpn = 0.0;
00187 G4double test, temp;
00188 for( G4int i=iBegin; i<=numSec; ++i )
00189 {
00190 temp = pi*i/(2.0*n*n);
00191 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(i*i)/(n*n) ) ) );
00192 if( temp < 1.0 )
00193 {
00194 if( test >= 1.0e-10 )anpn += temp*test;
00195 }
00196 else
00197 anpn += temp*test;
00198 }
00199 }
00200
00201 void
00202 G4RPGInelastic::CalculateMomenta(G4FastVector<G4ReactionProduct,256>& vec,
00203 G4int& vecLen,
00204 const G4HadProjectile* originalIncident,
00205 const G4DynamicParticle* originalTarget,
00206 G4ReactionProduct& modifiedOriginal,
00207 G4Nucleus& targetNucleus,
00208 G4ReactionProduct& currentParticle,
00209 G4ReactionProduct& targetParticle,
00210 G4bool& incidentHasChanged,
00211 G4bool& targetHasChanged,
00212 G4bool quasiElastic)
00213 {
00214 cache = 0;
00215 what = originalIncident->Get4Momentum().vect();
00216
00217 G4ReactionProduct leadingStrangeParticle;
00218
00219
00220
00221
00222
00223
00224
00225
00226 if( quasiElastic )
00227 {
00228 twoBody.ReactionStage(originalIncident, modifiedOriginal,
00229 incidentHasChanged, originalTarget,
00230 targetParticle, targetHasChanged,
00231 targetNucleus, currentParticle,
00232 vec, vecLen,
00233 false, leadingStrangeParticle);
00234 return;
00235 }
00236
00237 G4bool leadFlag = MarkLeadingStrangeParticle(currentParticle,
00238 targetParticle,
00239 leadingStrangeParticle );
00240
00241
00242
00243
00244 G4bool finishedGenXPt = false;
00245 G4bool annihilation = false;
00246 if( originalIncident->GetDefinition()->GetPDGEncoding() < 0 &&
00247 currentParticle.GetMass() == 0.0 && targetParticle.GetMass() == 0.0 )
00248 {
00249
00250 annihilation = true;
00251 G4double ekcor = 1.0;
00252 G4double ek = originalIncident->GetKineticEnergy();
00253 G4double ekOrg = ek;
00254
00255 const G4double tarmas = originalTarget->GetDefinition()->GetPDGMass();
00256 if( ek > 1.0*GeV )ekcor = 1./(ek/GeV);
00257 const G4double atomicWeight = targetNucleus.GetA_asInt();
00258 ek = 2*tarmas + ek*(1.+ekcor/atomicWeight);
00259 G4double tkin = targetNucleus.Cinema(ek);
00260 ek += tkin;
00261 ekOrg += tkin;
00262
00263
00264
00265
00266
00267 tkin = targetNucleus.AnnihilationEvaporationEffects(ek, ekOrg);
00268 ekOrg -= tkin;
00269 ekOrg = std::max( 0.0001*GeV, ekOrg );
00270 modifiedOriginal.SetKineticEnergy( ekOrg );
00271 G4double amas = originalIncident->GetDefinition()->GetPDGMass();
00272 G4double et = ekOrg + amas;
00273 G4double p = std::sqrt( std::abs(et*et-amas*amas) );
00274 G4double pp = modifiedOriginal.GetMomentum().mag();
00275 if( pp > 0.0 )
00276 {
00277 G4ThreeVector momentum = modifiedOriginal.GetMomentum();
00278 modifiedOriginal.SetMomentum( momentum * (p/pp) );
00279 }
00280 if( ekOrg <= 0.0001 )
00281 {
00282 modifiedOriginal.SetKineticEnergy( 0.0 );
00283 modifiedOriginal.SetMomentum( 0.0, 0.0, 0.0 );
00284 }
00285 }
00286
00287
00288
00289 const G4double twsup[] = { 1.0, 0.7, 0.5, 0.3, 0.2, 0.1 };
00290 G4double rand1 = G4UniformRand();
00291 G4double rand2 = G4UniformRand();
00292
00293
00294 G4ReactionProduct saveCurrent = currentParticle;
00295 G4ReactionProduct saveTarget = targetParticle;
00296 std::vector<G4ReactionProduct> savevec;
00297 for (G4int i = 0; i < vecLen; i++) savevec.push_back(*vec[i]);
00298
00299
00300
00301
00302
00303
00304
00305
00306 if( annihilation || vecLen > 5 ||
00307 ( modifiedOriginal.GetKineticEnergy()/GeV >= 1.0 &&
00308
00309 (((originalIncident->GetDefinition() == G4KaonPlus::KaonPlus() ||
00310 originalIncident->GetDefinition() == G4KaonMinus::KaonMinus() ||
00311 originalIncident->GetDefinition() == G4KaonZeroLong::KaonZeroLong() ||
00312 originalIncident->GetDefinition() == G4KaonZeroShort::KaonZeroShort()) &&
00313 rand1 < 0.5)
00314 || rand2 > twsup[vecLen]) ) )
00315
00316 finishedGenXPt =
00317 fragmentation.ReactionStage(originalIncident, modifiedOriginal,
00318 incidentHasChanged, originalTarget,
00319 targetParticle, targetHasChanged,
00320 targetNucleus, currentParticle,
00321 vec, vecLen,
00322 leadFlag, leadingStrangeParticle);
00323
00324 if (finishedGenXPt) return;
00325
00326 G4bool finishedTwoClu = false;
00327
00328 if (modifiedOriginal.GetTotalMomentum() < 1.0) {
00329 for (G4int i = 0; i < vecLen; i++) delete vec[i];
00330 vecLen = 0;
00331
00332 } else {
00333
00334
00335
00336
00337 if (!finishedGenXPt && annihilation) {
00338 currentParticle = saveCurrent;
00339 targetParticle = saveTarget;
00340 for (G4int i = 0; i < vecLen; i++) delete vec[i];
00341 vecLen = 0;
00342 vec.Initialize( 0 );
00343 for (G4int i = 0; i < G4int(savevec.size()); i++) {
00344 G4ReactionProduct* p = new G4ReactionProduct;
00345 *p = savevec[i];
00346 vec.SetElement( vecLen++, p );
00347 }
00348 }
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359 try
00360 {
00361 finishedTwoClu =
00362 twoCluster.ReactionStage(originalIncident, modifiedOriginal,
00363 incidentHasChanged, originalTarget,
00364 targetParticle, targetHasChanged,
00365 targetNucleus, currentParticle,
00366 vec, vecLen,
00367 leadFlag, leadingStrangeParticle);
00368 }
00369 catch(G4HadReentrentException aC)
00370 {
00371 aC.Report(G4cout);
00372 throw G4HadReentrentException(__FILE__, __LINE__, "Failing to calculate momenta");
00373 }
00374 }
00375
00376 if (finishedTwoClu) return;
00377
00378 twoBody.ReactionStage(originalIncident, modifiedOriginal,
00379 incidentHasChanged, originalTarget,
00380 targetParticle, targetHasChanged,
00381 targetNucleus, currentParticle,
00382 vec, vecLen,
00383 false, leadingStrangeParticle);
00384 }
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402 void
00403 G4RPGInelastic::SetUpChange(G4FastVector<G4ReactionProduct,256>& vec,
00404 G4int& vecLen,
00405 G4ReactionProduct& currentParticle,
00406 G4ReactionProduct& targetParticle,
00407 G4bool& incidentHasChanged )
00408 {
00409 theParticleChange.Clear();
00410 G4ParticleDefinition* aKaonZL = G4KaonZeroLong::KaonZeroLong();
00411 G4ParticleDefinition* aKaonZS = G4KaonZeroShort::KaonZeroShort();
00412 G4int i;
00413
00414 if (currentParticle.GetDefinition() == particleDef[k0]) {
00415 if (G4UniformRand() < 0.5) {
00416 currentParticle.SetDefinitionAndUpdateE(aKaonZL);
00417 incidentHasChanged = true;
00418 } else {
00419 currentParticle.SetDefinitionAndUpdateE(aKaonZS);
00420 }
00421 } else if (currentParticle.GetDefinition() == particleDef[k0b]) {
00422 if (G4UniformRand() < 0.5) {
00423 currentParticle.SetDefinitionAndUpdateE(aKaonZL);
00424 } else {
00425 currentParticle.SetDefinitionAndUpdateE(aKaonZS);
00426 incidentHasChanged = true;
00427 }
00428 }
00429
00430 if (targetParticle.GetDefinition() == particleDef[k0] ||
00431 targetParticle.GetDefinition() == particleDef[k0b] ) {
00432 if (G4UniformRand() < 0.5) {
00433 targetParticle.SetDefinitionAndUpdateE(aKaonZL);
00434 } else {
00435 targetParticle.SetDefinitionAndUpdateE(aKaonZS);
00436 }
00437 }
00438
00439 for (i = 0; i < vecLen; ++i) {
00440 if (vec[i]->GetDefinition() == particleDef[k0] ||
00441 vec[i]->GetDefinition() == particleDef[k0b] ) {
00442 if (G4UniformRand() < 0.5) {
00443 vec[i]->SetDefinitionAndUpdateE(aKaonZL);
00444 } else {
00445 vec[i]->SetDefinitionAndUpdateE(aKaonZS);
00446 }
00447 }
00448 }
00449
00450 if (incidentHasChanged) {
00451 G4DynamicParticle* p0 = new G4DynamicParticle;
00452 p0->SetDefinition(currentParticle.GetDefinition() );
00453 p0->SetMomentum(currentParticle.GetMomentum() );
00454 theParticleChange.AddSecondary( p0 );
00455 theParticleChange.SetStatusChange( stopAndKill );
00456 theParticleChange.SetEnergyChange( 0.0 );
00457
00458 } else {
00459 G4double p = currentParticle.GetMomentum().mag()/MeV;
00460 G4ThreeVector mom = currentParticle.GetMomentum();
00461 if (p > DBL_MIN)
00462 theParticleChange.SetMomentumChange(mom.x()/p, mom.y()/p, mom.z()/p );
00463 else
00464 theParticleChange.SetMomentumChange(0.0, 0.0, 1.0);
00465
00466 G4double aE = currentParticle.GetKineticEnergy();
00467 if (std::fabs(aE)<.1*eV) aE=.1*eV;
00468 theParticleChange.SetEnergyChange( aE );
00469 }
00470
00471 if (targetParticle.GetMass() > 0.0)
00472 {
00473 G4ThreeVector momentum = targetParticle.GetMomentum();
00474 momentum = momentum.rotate(cache, what);
00475 G4double targKE = targetParticle.GetKineticEnergy();
00476 G4ThreeVector dir(0.0, 0.0, 1.0);
00477 if (targKE < DBL_MIN)
00478 targKE = DBL_MIN;
00479 else
00480 dir = momentum/momentum.mag();
00481
00482 G4DynamicParticle* p1 =
00483 new G4DynamicParticle(targetParticle.GetDefinition(), dir, targKE);
00484
00485 theParticleChange.AddSecondary( p1 );
00486 }
00487
00488 G4DynamicParticle* p;
00489 for (i = 0; i < vecLen; ++i) {
00490 G4double secKE = vec[i]->GetKineticEnergy();
00491 G4ThreeVector momentum = vec[i]->GetMomentum();
00492 G4ThreeVector dir(0.0, 0.0, 1.0);
00493 if (secKE < DBL_MIN)
00494 secKE = DBL_MIN;
00495 else
00496 dir = momentum/momentum.mag();
00497
00498 p = new G4DynamicParticle(vec[i]->GetDefinition(), dir, secKE);
00499 theParticleChange.AddSecondary( p );
00500 delete vec[i];
00501 }
00502 }
00503
00504
00505 std::pair<G4int, G4double>
00506 G4RPGInelastic::interpolateEnergy(G4double e) const
00507 {
00508 G4int index = 29;
00509 G4double fraction = 0.0;
00510
00511 for (G4int i = 1; i < 30; i++) {
00512 if (e < energyScale[i]) {
00513 index = i-1;
00514 fraction = (e - energyScale[index]) / (energyScale[i] - energyScale[index]);
00515 break;
00516 }
00517 }
00518 return std::pair<G4int, G4double>(index, fraction);
00519 }
00520
00521
00522 G4int
00523 G4RPGInelastic::sampleFlat(std::vector<G4double> sigma) const
00524 {
00525 G4int i;
00526 G4double sum(0.);
00527 for (i = 0; i < G4int(sigma.size()); i++) sum += sigma[i];
00528
00529 G4double fsum = sum*G4UniformRand();
00530 G4double partialSum = 0.0;
00531 G4int channel = 0;
00532
00533 for (i = 0; i < G4int(sigma.size()); i++) {
00534 partialSum += sigma[i];
00535 if (fsum < partialSum) {
00536 channel = i;
00537 break;
00538 }
00539 }
00540
00541 return channel;
00542 }
00543
00544
00545 void G4RPGInelastic::CheckQnums(G4FastVector<G4ReactionProduct,256> &vec,
00546 G4int &vecLen,
00547 G4ReactionProduct ¤tParticle,
00548 G4ReactionProduct &targetParticle,
00549 G4double Q, G4double B, G4double S)
00550 {
00551 G4ParticleDefinition* projDef = currentParticle.GetDefinition();
00552 G4ParticleDefinition* targDef = targetParticle.GetDefinition();
00553 G4double chargeSum = projDef->GetPDGCharge() + targDef->GetPDGCharge();
00554 G4double baryonSum = projDef->GetBaryonNumber() + targDef->GetBaryonNumber();
00555 G4double strangenessSum = projDef->GetQuarkContent(3) -
00556 projDef->GetAntiQuarkContent(3) +
00557 targDef->GetQuarkContent(3) -
00558 targDef->GetAntiQuarkContent(3);
00559
00560 G4ParticleDefinition* secDef = 0;
00561 for (G4int i = 0; i < vecLen; i++) {
00562 secDef = vec[i]->GetDefinition();
00563 chargeSum += secDef->GetPDGCharge();
00564 baryonSum += secDef->GetBaryonNumber();
00565 strangenessSum += secDef->GetQuarkContent(3)
00566 - secDef->GetAntiQuarkContent(3);
00567 }
00568
00569 G4bool OK = true;
00570 if (chargeSum != Q) {
00571 G4cout << " Charge not conserved " << G4endl;
00572 OK = false;
00573 }
00574 if (baryonSum != B) {
00575 G4cout << " Baryon number not conserved " << G4endl;
00576 OK = false;
00577 }
00578 if (strangenessSum != S) {
00579 G4cout << " Strangeness not conserved " << G4endl;
00580 OK = false;
00581 }
00582
00583 if (!OK) {
00584 G4cout << " projectile: " << projDef->GetParticleName()
00585 << " target: " << targDef->GetParticleName() << G4endl;
00586 for (G4int i = 0; i < vecLen; i++) {
00587 secDef = vec[i]->GetDefinition();
00588 G4cout << secDef->GetParticleName() << " " ;
00589 }
00590 G4cout << G4endl;
00591 }
00592
00593 }
00594
00595
00596 const G4double G4RPGInelastic::energyScale[30] = {
00597 0.0, 0.01, 0.013, 0.018, 0.024, 0.032, 0.042, 0.056, 0.075, 0.1,
00598 0.13, 0.18, 0.24, 0.32, 0.42, 0.56, 0.75, 1.0, 1.3, 1.8,
00599 2.4, 3.2, 4.2, 5.6, 7.5, 10.0, 13.0, 18.0, 24.0, 32.0 };
00600
00601
00602