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 "G4RPGFragmentation.hh"
00033 #include "G4PhysicalConstants.hh"
00034 #include "G4SystemOfUnits.hh"
00035 #include "G4AntiProton.hh"
00036 #include "G4AntiNeutron.hh"
00037 #include "Randomize.hh"
00038 #include "G4Poisson.hh"
00039 #include "G4HadReentrentException.hh"
00040
00041 G4RPGFragmentation::G4RPGFragmentation()
00042 : G4RPGReaction()
00043 {
00044 for (G4int i = 0; i < 20; i++) dndl[i] = 0.0;
00045 }
00046
00047
00048 void G4RPGFragmentation::
00049 FragmentationIntegral(G4double pt, G4double et, G4double parMass, G4double secMass)
00050 {
00051 pt = std::max( 0.001, pt );
00052 G4double dx = 1./(19.*pt);
00053 G4double x;
00054 G4double term1;
00055 G4double term2;
00056
00057 for (G4int i = 1; i < 20; i++) {
00058 x = (G4double(i) - 0.5)*dx;
00059 term1 = 1. + parMass*parMass*x*x;
00060 term2 = pt*x*et*pt*x*et + pt*pt + secMass*secMass;
00061 dndl[i] = dx / std::sqrt( term1*term1*term1*term2 )
00062 + dndl[i-1];
00063 }
00064 }
00065
00066
00067 G4bool G4RPGFragmentation::
00068 ReactionStage(const G4HadProjectile* originalIncident,
00069 G4ReactionProduct& modifiedOriginal,
00070 G4bool& incidentHasChanged,
00071 const G4DynamicParticle* originalTarget,
00072 G4ReactionProduct& targetParticle,
00073 G4bool& targetHasChanged,
00074 const G4Nucleus& targetNucleus,
00075 G4ReactionProduct& currentParticle,
00076 G4FastVector<G4ReactionProduct,256>& vec,
00077 G4int& vecLen,
00078 G4bool leadFlag,
00079 G4ReactionProduct& leadingStrangeParticle)
00080 {
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095 if (vecLen == 0) return false;
00096
00097 G4ParticleDefinition* aPiMinus = G4PionMinus::PionMinus();
00098 G4ParticleDefinition* aProton = G4Proton::Proton();
00099 G4ParticleDefinition* aNeutron = G4Neutron::Neutron();
00100 G4ParticleDefinition* aPiPlus = G4PionPlus::PionPlus();
00101 G4ParticleDefinition* aPiZero = G4PionZero::PionZero();
00102
00103 G4int i, l;
00104 G4bool veryForward = false;
00105
00106 const G4double ekOriginal = modifiedOriginal.GetKineticEnergy()/GeV;
00107 const G4double etOriginal = modifiedOriginal.GetTotalEnergy()/GeV;
00108 const G4double mOriginal = modifiedOriginal.GetMass()/GeV;
00109 const G4double pOriginal = modifiedOriginal.GetMomentum().mag()/GeV;
00110 G4double targetMass = targetParticle.GetDefinition()->GetPDGMass()/GeV;
00111 G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
00112 targetMass*targetMass +
00113 2.0*targetMass*etOriginal );
00114 G4double currentMass = currentParticle.GetMass()/GeV;
00115 targetMass = targetParticle.GetMass()/GeV;
00116
00117
00118
00119
00120 for (i=0; i<vecLen; ++i) {
00121 G4int itemp = G4int( G4UniformRand()*vecLen );
00122 G4ReactionProduct pTemp = *vec[itemp];
00123 *vec[itemp] = *vec[i];
00124 *vec[i] = pTemp;
00125 }
00126
00127 if (currentMass == 0.0 && targetMass == 0.0) {
00128
00129
00130
00131 G4double ek = currentParticle.GetKineticEnergy();
00132 G4ThreeVector mom = currentParticle.GetMomentum();
00133 currentParticle = *vec[0];
00134 currentParticle.SetSide(1);
00135 targetParticle = *vec[1];
00136 targetParticle.SetSide(-1);
00137 for( i=0; i<(vecLen-2); ++i )*vec[i] = *vec[i+2];
00138 G4ReactionProduct *temp = vec[vecLen-1];
00139 delete temp;
00140 temp = vec[vecLen-2];
00141 delete temp;
00142 vecLen -= 2;
00143 currentMass = currentParticle.GetMass()/GeV;
00144 targetMass = targetParticle.GetMass()/GeV;
00145 incidentHasChanged = true;
00146 targetHasChanged = true;
00147 currentParticle.SetKineticEnergy( ek );
00148 currentParticle.SetMomentum(mom);
00149 veryForward = true;
00150 }
00151 const G4double atomicWeight = targetNucleus.GetA_asInt();
00152 const G4double atomicNumber = targetNucleus.GetZ_asInt();
00153 const G4double protonMass = aProton->GetPDGMass();
00154
00155 if (originalIncident->GetDefinition()->GetParticleSubType() == "kaon"
00156 && G4UniformRand() >= 0.7) {
00157 G4ReactionProduct temp = currentParticle;
00158 currentParticle = targetParticle;
00159 targetParticle = temp;
00160 incidentHasChanged = true;
00161 targetHasChanged = true;
00162 currentMass = currentParticle.GetMass()/GeV;
00163 targetMass = targetParticle.GetMass()/GeV;
00164 }
00165 const G4double afc = std::min( 0.75,
00166 0.312+0.200*std::log(std::log(centerofmassEnergy*centerofmassEnergy))+
00167 std::pow(centerofmassEnergy*centerofmassEnergy,1.5)/6000.0 );
00168
00169 G4double freeEnergy = centerofmassEnergy-currentMass-targetMass;
00170 G4double forwardEnergy = freeEnergy/2.;
00171 G4int forwardCount = 1;
00172
00173 G4double backwardEnergy = freeEnergy/2.;
00174 G4int backwardCount = 1;
00175
00176 if(veryForward) {
00177 if(currentParticle.GetSide()==-1)
00178 {
00179 forwardEnergy += currentMass;
00180 forwardCount --;
00181 backwardEnergy -= currentMass;
00182 backwardCount ++;
00183 }
00184 if(targetParticle.GetSide()!=-1)
00185 {
00186 backwardEnergy += targetMass;
00187 backwardCount --;
00188 forwardEnergy -= targetMass;
00189 forwardCount ++;
00190 }
00191 }
00192
00193 for (i=0; i<vecLen; ++i) {
00194 if( vec[i]->GetSide() == -1 )
00195 {
00196 ++backwardCount;
00197 backwardEnergy -= vec[i]->GetMass()/GeV;
00198 } else {
00199 ++forwardCount;
00200 forwardEnergy -= vec[i]->GetMass()/GeV;
00201 }
00202 }
00203
00204
00205
00206
00207
00208 if (backwardEnergy < 0.0) {
00209 for (i = 0; i < vecLen; ++i) {
00210 if (vec[i]->GetSide() == -1) {
00211 backwardEnergy += vec[i]->GetMass()/GeV;
00212 --backwardCount;
00213 vec[i]->SetSide(1);
00214 forwardEnergy -= vec[i]->GetMass()/GeV;
00215 ++forwardCount;
00216 if (backwardEnergy > 0.0) break;
00217 }
00218 }
00219 }
00220
00221 if (forwardEnergy < 0.0) {
00222 for (i = 0; i < vecLen; ++i) {
00223 if (vec[i]->GetSide() == 1) {
00224 forwardEnergy += vec[i]->GetMass()/GeV;
00225 --forwardCount;
00226 vec[i]->SetSide(-1);
00227 backwardEnergy -= vec[i]->GetMass()/GeV;
00228 ++backwardCount;
00229 if (forwardEnergy > 0.0) break;
00230 }
00231 }
00232 }
00233
00234
00235
00236
00237 if (forwardEnergy > 0.0 && backwardEnergy < 0.0) {
00238 forwardEnergy += backwardEnergy;
00239 backwardEnergy = 0;
00240 }
00241
00242
00243 if (forwardEnergy + backwardEnergy < 0.0) return false;
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258 G4double xtarg;
00259 if( centerofmassEnergy < (2.0+G4UniformRand()) )
00260 xtarg = afc * (std::pow(atomicWeight,0.33)-1.0) * (2.0*backwardCount+vecLen+2)/2.0;
00261 else
00262 xtarg = afc * (std::pow(atomicWeight,0.33)-1.0) * (2.0*backwardCount);
00263 if( xtarg <= 0.0 )xtarg = 0.01;
00264 G4int nuclearExcitationCount = G4Poisson( xtarg );
00265
00266
00267 if(atomicWeight<1.0001) nuclearExcitationCount = 0;
00268 G4int extraNucleonCount = 0;
00269 G4double extraNucleonMass = 0.0;
00270
00271 if (nuclearExcitationCount > 0) {
00272 const G4double nucsup[] = { 1.00, 0.7, 0.5, 0.4, 0.35, 0.3 };
00273 const G4double psup[] = { 3., 6., 20., 50., 100., 1000. };
00274 G4int momentumBin = 0;
00275 while( (momentumBin < 6) &&
00276 (modifiedOriginal.GetTotalMomentum()/GeV > psup[momentumBin]) )
00277 ++momentumBin;
00278 momentumBin = std::min( 5, momentumBin );
00279
00280
00281
00282
00283
00284 for (i = 0; i < nuclearExcitationCount; ++i) {
00285 G4ReactionProduct * pVec = new G4ReactionProduct();
00286 if (G4UniformRand() < nucsup[momentumBin]) {
00287
00288 if (G4UniformRand() > 1.0-atomicNumber/atomicWeight)
00289 pVec->SetDefinition( aProton );
00290 else
00291 pVec->SetDefinition( aNeutron );
00292
00293
00294
00295 pVec->SetSide( -2 );
00296 ++extraNucleonCount;
00297 extraNucleonMass += pVec->GetMass()/GeV;
00298
00299 pVec->SetNewlyAdded( true );
00300 vec.SetElement( vecLen++, pVec );
00301 ++backwardCount;
00302
00303 } else {
00304
00305 G4double ran = G4UniformRand();
00306 if( ran < 0.3181 )
00307 pVec->SetDefinition( aPiPlus );
00308 else if( ran < 0.6819 )
00309 pVec->SetDefinition( aPiZero );
00310 else
00311 pVec->SetDefinition( aPiMinus );
00312
00313 if (backwardEnergy > pVec->GetMass()/GeV) {
00314 backwardEnergy -= pVec->GetMass()/GeV;
00315 ++backwardCount;
00316 pVec->SetSide( -1 );
00317 pVec->SetNewlyAdded( true );
00318 vec.SetElement( vecLen++, pVec );
00319 }
00320
00321
00322
00323 }
00324 }
00325 }
00326
00327
00328
00329
00330 G4ReactionProduct pseudoParticle[8];
00331 for (i = 0; i < 8; ++i) pseudoParticle[i].SetZero();
00332
00333 pseudoParticle[0].SetMass( mOriginal*GeV );
00334 pseudoParticle[0].SetMomentum( 0.0, 0.0, pOriginal*GeV );
00335 pseudoParticle[0].SetTotalEnergy(
00336 std::sqrt( pOriginal*pOriginal + mOriginal*mOriginal )*GeV );
00337
00338 pseudoParticle[1].SetMass(protonMass);
00339 pseudoParticle[1].SetTotalEnergy(protonMass);
00340
00341 pseudoParticle[3].SetMass(protonMass*(1+extraNucleonCount) );
00342 pseudoParticle[3].SetTotalEnergy(protonMass*(1+extraNucleonCount) );
00343
00344 pseudoParticle[2] = pseudoParticle[0] + pseudoParticle[1];
00345 pseudoParticle[3] = pseudoParticle[3] + pseudoParticle[0];
00346
00347 pseudoParticle[0].Lorentz( pseudoParticle[0], pseudoParticle[2] );
00348 pseudoParticle[1].Lorentz( pseudoParticle[1], pseudoParticle[2] );
00349
00350
00351
00352
00353 G4double aspar, pt, et, x, pp, pp1, wgt;
00354 G4int innerCounter, outerCounter;
00355 G4bool eliminateThisParticle, resetEnergies, constantCrossSection;
00356
00357 G4double forwardKinetic = 0.0;
00358 G4double backwardKinetic = 0.0;
00359
00360
00361
00362
00363
00364
00365 G4int backwardNucleonCount = 0;
00366 G4double totalEnergy, kineticEnergy, vecMass;
00367 G4double phi;
00368
00369 for (i = vecLen-1; i >= 0; --i) {
00370
00371 if (vec[i]->GetNewlyAdded()) {
00372 if (vec[i]->GetSide() == -2) {
00373 if (backwardNucleonCount < 18) {
00374 if (vec[i]->GetDefinition()->GetParticleSubType() == "pi") {
00375 for (G4int j = 0; j < vecLen; j++) delete vec[j];
00376 vecLen = 0;
00377 throw G4HadReentrentException(__FILE__, __LINE__,
00378 "G4RPGFragmentation::ReactionStage : a pion has been counted as a backward nucleon");
00379 }
00380 vec[i]->SetSide(-3);
00381 ++backwardNucleonCount;
00382 continue;
00383
00384
00385 }
00386 }
00387 }
00388
00389
00390
00391
00392 vecMass = vec[i]->GetMass()/GeV;
00393 G4double ran = -std::log(1.0-G4UniformRand())/3.5;
00394
00395 if (vec[i]->GetSide() == -2) {
00396 aspar = 0.20;
00397 pt = std::sqrt( std::pow( ran, 1.2 ) );
00398
00399 } else {
00400 if (vec[i]->GetDefinition()->GetParticleSubType() == "pi") {
00401 aspar = 0.75;
00402 pt = std::sqrt( std::pow( ran, 1.7 ) );
00403 } else if (vec[i]->GetDefinition()->GetParticleSubType() == "kaon") {
00404 aspar = 0.70;
00405 pt = std::sqrt( std::pow( ran, 1.7 ) );
00406 } else {
00407 aspar = 0.65;
00408 pt = std::sqrt( std::pow( ran, 1.5 ) );
00409 }
00410 }
00411
00412 pt = std::max( 0.001, pt );
00413 phi = G4UniformRand()*twopi;
00414 vec[i]->SetMomentum( pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV );
00415 if (vec[i]->GetSide() > 0)
00416 et = pseudoParticle[0].GetTotalEnergy()/GeV;
00417 else
00418 et = pseudoParticle[1].GetTotalEnergy()/GeV;
00419
00420
00421
00422
00423 outerCounter = 0;
00424 eliminateThisParticle = true;
00425 resetEnergies = true;
00426 dndl[0] = 0.0;
00427
00428 while (++outerCounter < 3) {
00429
00430 FragmentationIntegral(pt, et, aspar, vecMass);
00431 innerCounter = 0;
00432 vec[i]->SetMomentum( pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV );
00433
00434
00435
00436 while (++innerCounter < 7) {
00437
00438 ran = G4UniformRand()*dndl[19];
00439 l = 1;
00440 while( ( ran > dndl[l] ) && ( l < 19 ) ) l++;
00441 x = (G4double(l-1) + G4UniformRand())/19.;
00442 if (vec[i]->GetSide() < 0) x *= -1.;
00443 vec[i]->SetMomentum( x*et*GeV );
00444 totalEnergy = std::sqrt( x*et*x*et + pt*pt + vecMass*vecMass );
00445 vec[i]->SetTotalEnergy( totalEnergy*GeV );
00446 kineticEnergy = vec[i]->GetKineticEnergy()/GeV;
00447
00448 if (vec[i]->GetSide() > 0) {
00449 if( (forwardKinetic+kineticEnergy) < 0.95*forwardEnergy ) {
00450
00451
00452 pseudoParticle[4] = pseudoParticle[4] + (*vec[i]);
00453 forwardKinetic += kineticEnergy;
00454 outerCounter = 2;
00455 eliminateThisParticle = false;
00456 resetEnergies = false;
00457 break;
00458 }
00459 if( innerCounter > 5 )break;
00460 if( backwardEnergy >= vecMass )
00461 {
00462 vec[i]->SetSide(-1);
00463 forwardEnergy += vecMass;
00464 backwardEnergy -= vecMass;
00465 ++backwardCount;
00466 }
00467 } else {
00468
00469
00470
00471 G4double xxx = 0.999;
00472 if (extraNucleonCount < 20) xxx = 0.95+0.05*extraNucleonCount/20.0;
00473
00474 if ((backwardKinetic+kineticEnergy) < xxx*backwardEnergy) {
00475 pseudoParticle[5] = pseudoParticle[5] + (*vec[i]);
00476 backwardKinetic += kineticEnergy;
00477 outerCounter = 2;
00478 eliminateThisParticle = false;
00479 resetEnergies = false;
00480 break;
00481 }
00482 if (innerCounter > 5) break;
00483 if (forwardEnergy >= vecMass) {
00484 vec[i]->SetSide(1);
00485 forwardEnergy -= vecMass;
00486 backwardEnergy += vecMass;
00487 backwardCount--;
00488 }
00489 }
00490 G4ThreeVector momentum = vec[i]->GetMomentum();
00491 vec[i]->SetMomentum( momentum.x() * 0.9, momentum.y() * 0.9 );
00492 pt *= 0.9;
00493 dndl[19] *= 0.9;
00494 }
00495
00496
00497
00498
00499
00500
00501 if (resetEnergies)
00502 ReduceEnergiesOfSecondaries(i+1, forwardKinetic, backwardKinetic,
00503 vec, vecLen,
00504 pseudoParticle[4], pseudoParticle[5],
00505 pt);
00506
00507 }
00508
00509 if (eliminateThisParticle && vec[i]->GetMayBeKilled()) {
00510
00511
00512 if (vec[i]->GetSide() > 0) {
00513 --forwardCount;
00514 forwardEnergy += vecMass;
00515 } else {
00516 --backwardCount;
00517 if (vec[i]->GetSide() == -2) {
00518 --extraNucleonCount;
00519 extraNucleonMass -= vecMass;
00520 } else {
00521 backwardEnergy += vecMass;
00522 }
00523 }
00524
00525 for( G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];
00526 G4ReactionProduct* temp = vec[vecLen-1];
00527 delete temp;
00528
00529
00530 if( --vecLen == 0 ){
00531 G4cout << " FALSE RETURN DUE TO ENERGY BALANCE " << G4endl;
00532 return false;
00533 }
00534 }
00535 }
00536
00537
00538
00539 G4double forwardKEDiff = forwardEnergy - forwardKinetic;
00540 G4double backwardKEDiff = backwardEnergy - backwardKinetic;
00541
00542 if (forwardKEDiff < 0.0 || backwardKEDiff < 0.0) {
00543 ReduceEnergiesOfSecondaries(0, forwardKinetic, backwardKinetic,
00544 vec, vecLen,
00545 pseudoParticle[4], pseudoParticle[5],
00546 pt);
00547
00548 forwardKEDiff = forwardEnergy - forwardKinetic;
00549 backwardKEDiff = backwardEnergy - backwardKinetic;
00550 if (backwardKEDiff < 0.0) {
00551 if (forwardKEDiff + backwardKEDiff > 0.0) {
00552 backwardEnergy = backwardKinetic;
00553 forwardEnergy += backwardKEDiff;
00554 forwardKEDiff = forwardEnergy - forwardKinetic;
00555 backwardKEDiff = 0.0;
00556 } else {
00557 G4cout << " False return due to insufficient backward energy " << G4endl;
00558 return false;
00559 }
00560 }
00561
00562 if (forwardKEDiff < 0.0) {
00563 if (forwardKEDiff + backwardKEDiff > 0.0) {
00564 forwardEnergy = forwardKinetic;
00565 backwardEnergy += forwardKEDiff;
00566 backwardKEDiff = backwardEnergy - backwardKinetic;
00567 forwardKEDiff = 0.0;
00568 } else {
00569 G4cout << " False return due to insufficient forward energy " << G4endl;
00570 return false;
00571 }
00572 }
00573 }
00574
00575
00576
00577
00578
00579
00580 G4double ran = -std::log(1.0-G4UniformRand());
00581 if (currentParticle.GetDefinition()->GetParticleSubType() == "pi") {
00582 aspar = 0.60;
00583 pt = std::sqrt( std::pow( ran/6.0, 1.7 ) );
00584 } else if (currentParticle.GetDefinition()->GetParticleSubType() == "kaon") {
00585 aspar = 0.50;
00586 pt = std::sqrt( std::pow( ran/5.0, 1.4 ) );
00587 } else {
00588 aspar = 0.40;
00589 pt = std::sqrt( std::pow( ran/4.0, 1.2 ) );
00590 }
00591
00592 phi = G4UniformRand()*twopi;
00593 currentParticle.SetMomentum(pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV);
00594 et = pseudoParticle[0].GetTotalEnergy()/GeV;
00595 dndl[0] = 0.0;
00596 vecMass = currentParticle.GetMass()/GeV;
00597
00598 FragmentationIntegral(pt, et, aspar, vecMass);
00599
00600 ran = G4UniformRand()*dndl[19];
00601 l = 1;
00602 while( ( ran > dndl[l] ) && ( l < 19 ) ) l++;
00603 x = (G4double(l-1) + G4UniformRand())/19.;
00604 currentParticle.SetMomentum( x*et*GeV );
00605
00606 if (forwardEnergy < forwardKinetic) {
00607 totalEnergy = vecMass + 0.04*std::fabs(normal());
00608 G4cout << " Not enough forward energy: forwardEnergy = "
00609 << forwardEnergy << " forwardKinetic = "
00610 << forwardKinetic << " total energy left = "
00611 << backwardKEDiff + forwardKEDiff << G4endl;
00612 } else {
00613 totalEnergy = vecMass + forwardEnergy - forwardKinetic;
00614 forwardKinetic = forwardEnergy;
00615 }
00616 currentParticle.SetTotalEnergy( totalEnergy*GeV );
00617 pp = std::sqrt(std::abs( totalEnergy*totalEnergy - vecMass*vecMass) )*GeV;
00618 pp1 = currentParticle.GetMomentum().mag();
00619
00620 if (pp1 < 1.0e-6*GeV) {
00621 G4ThreeVector iso = Isotropic(pp);
00622 currentParticle.SetMomentum( iso.x(), iso.y(), iso.z() );
00623 } else {
00624 currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
00625 }
00626 pseudoParticle[4] = pseudoParticle[4] + currentParticle;
00627
00628
00629
00630
00631
00632 if (backwardNucleonCount < 18) {
00633 targetParticle.SetSide(-3);
00634 ++backwardNucleonCount;
00635
00636 } else {
00637
00638
00639
00640 vecMass = targetParticle.GetMass()/GeV;
00641 ran = -std::log(1.0-G4UniformRand());
00642 aspar = 0.40;
00643 pt = std::max( 0.001, std::sqrt( std::pow( ran/4.0, 1.2 ) ) );
00644 phi = G4UniformRand()*twopi;
00645 targetParticle.SetMomentum(pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV);
00646 et = pseudoParticle[1].GetTotalEnergy()/GeV;
00647 outerCounter = 0;
00648 innerCounter = 0;
00649 G4bool marginalEnergy = true;
00650 dndl[0] = 0.0;
00651 G4double xxx = 0.999;
00652 if( extraNucleonCount < 20 ) xxx = 0.95+0.05*extraNucleonCount/20.0;
00653 G4ThreeVector momentum;
00654
00655 while (++outerCounter < 4) {
00656 FragmentationIntegral(pt, et, aspar, vecMass);
00657
00658 for (innerCounter = 0; innerCounter < 6; innerCounter++) {
00659 ran = G4UniformRand()*dndl[19];
00660 l = 1;
00661 while( ( ran > dndl[l] ) && ( l < 19 ) ) l++;
00662 x = -(G4double(l-1) + G4UniformRand())/19.;
00663 targetParticle.SetMomentum( x*et*GeV );
00664 totalEnergy = std::sqrt(x*et*x*et + pt*pt + vecMass*vecMass);
00665 targetParticle.SetTotalEnergy( totalEnergy*GeV );
00666
00667 if ((backwardKinetic+totalEnergy-vecMass) < xxx*backwardEnergy) {
00668 pseudoParticle[5] = pseudoParticle[5] + targetParticle;
00669 backwardKinetic += totalEnergy - vecMass;
00670 outerCounter = 3;
00671 marginalEnergy = false;
00672 break;
00673 }
00674 momentum = targetParticle.GetMomentum();
00675 targetParticle.SetMomentum(momentum.x() * 0.9, momentum.y() * 0.9);
00676 pt *= 0.9;
00677 dndl[19] *= 0.9;
00678 }
00679 }
00680
00681 if (marginalEnergy) {
00682 G4cout << " Extra backward kinetic energy = "
00683 << 0.999*backwardEnergy - backwardKinetic << G4endl;
00684 totalEnergy = vecMass + 0.999*backwardEnergy - backwardKinetic;
00685 targetParticle.SetTotalEnergy(totalEnergy*GeV);
00686 pp = std::sqrt(std::abs(totalEnergy*totalEnergy - vecMass*vecMass) )*GeV;
00687 targetParticle.SetMomentum(momentum.x()/0.9, momentum.y()/0.9);
00688 pp1 = targetParticle.GetMomentum().mag();
00689 targetParticle.SetMomentum(targetParticle.GetMomentum() * pp/pp1 );
00690 pseudoParticle[5] = pseudoParticle[5] + targetParticle;
00691 backwardKinetic = 0.999*backwardEnergy;
00692 }
00693
00694 }
00695
00696 if (backwardEnergy < backwardKinetic)
00697 G4cout << " Backward Edif = " << backwardEnergy - backwardKinetic << G4endl;
00698 if (forwardEnergy != forwardKinetic)
00699 G4cout << " Forward Edif = " << forwardEnergy - forwardKinetic << G4endl;
00700
00701
00702
00703
00704
00705
00706
00707
00708
00709 pseudoParticle[6].Lorentz( pseudoParticle[3], pseudoParticle[2] );
00710 pseudoParticle[6] = pseudoParticle[6] - pseudoParticle[4];
00711 pseudoParticle[6] = pseudoParticle[6] - pseudoParticle[5];
00712
00713 if (backwardNucleonCount == 1) {
00714
00715
00716
00717 G4double ekin =
00718 std::min(backwardEnergy-backwardKinetic, centerofmassEnergy/2.0-protonMass/GeV);
00719
00720 if( ekin < 0.04 )ekin = 0.04 * std::fabs( normal() );
00721 vecMass = targetParticle.GetMass()/GeV;
00722 totalEnergy = ekin + vecMass;
00723 targetParticle.SetTotalEnergy( totalEnergy*GeV );
00724 pp = std::sqrt(std::abs(totalEnergy*totalEnergy - vecMass*vecMass) )*GeV;
00725 pp1 = pseudoParticle[6].GetMomentum().mag();
00726 if (pp1 < 1.0e-6*GeV) {
00727 G4ThreeVector iso = Isotropic(pp);
00728 targetParticle.SetMomentum( iso.x(), iso.y(), iso.z() );
00729 } else {
00730 targetParticle.SetMomentum( pseudoParticle[6].GetMomentum() * (pp/pp1));
00731 }
00732 pseudoParticle[5] = pseudoParticle[5] + targetParticle;
00733
00734 } else if (backwardNucleonCount > 1) {
00735
00736
00737 G4int tempCount = 5;
00738 if (backwardNucleonCount < 5) tempCount = backwardNucleonCount;
00739 tempCount -= 2;
00740
00741 G4double clusterMass = 0.;
00742 if (targetParticle.GetSide() == -3)
00743 clusterMass = targetParticle.GetMass()/GeV;
00744 for (i = 0; i < vecLen; ++i)
00745 if (vec[i]->GetSide() == -3) clusterMass += vec[i]->GetMass()/GeV;
00746 clusterMass += backwardEnergy - backwardKinetic;
00747
00748 totalEnergy = pseudoParticle[6].GetTotalEnergy()/GeV;
00749 pseudoParticle[6].SetMass(clusterMass*GeV);
00750
00751 pp = std::sqrt(std::abs(totalEnergy*totalEnergy -
00752 clusterMass*clusterMass) )*GeV;
00753 pp1 = pseudoParticle[6].GetMomentum().mag();
00754 if (pp1 < 1.0e-6*GeV) {
00755 G4ThreeVector iso = Isotropic(pp);
00756 pseudoParticle[6].SetMomentum(iso.x(), iso.y(), iso.z());
00757 } else {
00758 pseudoParticle[6].SetMomentum(pseudoParticle[6].GetMomentum() * (-pp/pp1));
00759 }
00760
00761 std::vector<G4ReactionProduct*> tempList;
00762 if (targetParticle.GetSide() == -3) tempList.push_back(&targetParticle);
00763 for (i = 0; i < vecLen; ++i)
00764 if (vec[i]->GetSide() == -3) tempList.push_back(vec[i]);
00765
00766 constantCrossSection = true;
00767
00768 if (tempList.size() > 1) {
00769 G4int n_entry = 0;
00770 wgt = GenerateNBodyEventT(pseudoParticle[6].GetMass(),
00771 constantCrossSection, tempList);
00772
00773 if (targetParticle.GetSide() == -3) {
00774 targetParticle = *tempList[0];
00775 targetParticle.Lorentz(targetParticle, pseudoParticle[6]);
00776 n_entry++;
00777 }
00778
00779 for (i = 0; i < vecLen; ++i) {
00780 if (vec[i]->GetSide() == -3) {
00781 *vec[i] = *tempList[n_entry];
00782 vec[i]->Lorentz(*vec[i], pseudoParticle[6]);
00783 n_entry++;
00784 }
00785 }
00786 }
00787 } else return false;
00788
00789 if (vecLen == 0) return false;
00790
00791
00792
00793 currentParticle.Lorentz( currentParticle, pseudoParticle[1] );
00794 targetParticle.Lorentz( targetParticle, pseudoParticle[1] );
00795 for (i = 0; i < vecLen; ++i) vec[i]->Lorentz(*vec[i], pseudoParticle[1]);
00796
00797
00798
00799
00800
00801
00802
00803
00804
00805 G4bool leadingStrangeParticleHasChanged = true;
00806 if (leadFlag)
00807 {
00808 if (currentParticle.GetDefinition() == leadingStrangeParticle.GetDefinition())
00809 leadingStrangeParticleHasChanged = false;
00810 if (leadingStrangeParticleHasChanged &&
00811 (targetParticle.GetDefinition() == leadingStrangeParticle.GetDefinition()) )
00812 leadingStrangeParticleHasChanged = false;
00813 if( leadingStrangeParticleHasChanged )
00814 {
00815 for( i=0; i<vecLen; i++ )
00816 {
00817 if( vec[i]->GetDefinition() == leadingStrangeParticle.GetDefinition() )
00818 {
00819 leadingStrangeParticleHasChanged = false;
00820 break;
00821 }
00822 }
00823 }
00824 if( leadingStrangeParticleHasChanged )
00825 {
00826 G4bool leadTest =
00827 (leadingStrangeParticle.GetDefinition()->GetParticleSubType() == "kaon" ||
00828 leadingStrangeParticle.GetDefinition()->GetParticleSubType() == "pi");
00829 G4bool targetTest =
00830 (targetParticle.GetDefinition()->GetParticleSubType() == "kaon" ||
00831 targetParticle.GetDefinition()->GetParticleSubType() == "pi");
00832
00833
00834
00835 if( (leadTest&&targetTest) || !(leadTest||targetTest) )
00836 {
00837 targetParticle.SetDefinitionAndUpdateE( leadingStrangeParticle.GetDefinition() );
00838 targetHasChanged = true;
00839 }
00840 else
00841 {
00842 currentParticle.SetDefinitionAndUpdateE( leadingStrangeParticle.GetDefinition() );
00843 incidentHasChanged = false;
00844 }
00845 }
00846 }
00847
00848
00849
00850
00851 std::pair<G4int, G4int> finalStateNucleons =
00852 GetFinalStateNucleons(originalTarget, vec, vecLen);
00853
00854 G4int protonsInFinalState = finalStateNucleons.first;
00855 G4int neutronsInFinalState = finalStateNucleons.second;
00856
00857 G4int numberofFinalStateNucleons =
00858 protonsInFinalState + neutronsInFinalState;
00859
00860 if (currentParticle.GetDefinition()->GetBaryonNumber() == 1 &&
00861 targetParticle.GetDefinition()->GetBaryonNumber() == 1 &&
00862 originalIncident->GetDefinition()->GetPDGMass() <
00863 G4Lambda::Lambda()->GetPDGMass())
00864 numberofFinalStateNucleons++;
00865
00866 numberofFinalStateNucleons = std::max(1, numberofFinalStateNucleons);
00867
00868 G4int PinNucleus = std::max(0,
00869 G4int(targetNucleus.GetZ_asInt()) - protonsInFinalState);
00870 G4int NinNucleus = std::max(0,
00871 G4int(targetNucleus.GetA_asInt()-targetNucleus.GetZ_asInt()) - neutronsInFinalState);
00872
00873 pseudoParticle[3].SetMomentum( 0.0, 0.0, pOriginal*GeV );
00874 pseudoParticle[3].SetMass( mOriginal*GeV );
00875 pseudoParticle[3].SetTotalEnergy(
00876 std::sqrt( pOriginal*pOriginal + mOriginal*mOriginal )*GeV );
00877
00878 G4ParticleDefinition * aOrgDef = modifiedOriginal.GetDefinition();
00879 G4int diff = 0;
00880 if(aOrgDef == G4Proton::Proton() || aOrgDef == G4Neutron::Neutron() ) diff = 1;
00881 if(numberofFinalStateNucleons == 1) diff = 0;
00882 pseudoParticle[4].SetMomentum( 0.0, 0.0, 0.0 );
00883 pseudoParticle[4].SetMass( protonMass*(numberofFinalStateNucleons-diff) );
00884 pseudoParticle[4].SetTotalEnergy( protonMass*(numberofFinalStateNucleons-diff) );
00885
00886 G4double theoreticalKinetic =
00887 pseudoParticle[3].GetTotalEnergy() + pseudoParticle[4].GetTotalEnergy() -
00888 currentParticle.GetMass() - targetParticle.GetMass();
00889 for (i = 0; i < vecLen; ++i) theoreticalKinetic -= vec[i]->GetMass();
00890
00891 G4double simulatedKinetic =
00892 currentParticle.GetKineticEnergy() + targetParticle.GetKineticEnergy();
00893 for (i = 0; i < vecLen; ++i)
00894 simulatedKinetic += vec[i]->GetKineticEnergy();
00895
00896 pseudoParticle[5] = pseudoParticle[3] + pseudoParticle[4];
00897 pseudoParticle[3].Lorentz( pseudoParticle[3], pseudoParticle[5] );
00898 pseudoParticle[4].Lorentz( pseudoParticle[4], pseudoParticle[5] );
00899
00900 pseudoParticle[7].SetZero();
00901 pseudoParticle[7] = pseudoParticle[7] + currentParticle;
00902 pseudoParticle[7] = pseudoParticle[7] + targetParticle;
00903 for (i = 0; i < vecLen; ++i)
00904 pseudoParticle[7] = pseudoParticle[7] + *vec[i];
00905
00906
00907
00908
00909
00910
00911
00912
00913
00914
00915
00916
00917
00918
00919
00920
00921
00922
00923
00924
00925
00926
00927
00928
00929
00930
00931
00932
00933
00934
00935
00936
00937
00938
00939
00940
00941
00942
00943
00944
00945
00946
00947
00948
00949 if (simulatedKinetic != 0.0) {
00950 wgt = theoreticalKinetic/simulatedKinetic;
00951 theoreticalKinetic = currentParticle.GetKineticEnergy() * wgt;
00952 simulatedKinetic = theoreticalKinetic;
00953 currentParticle.SetKineticEnergy(theoreticalKinetic);
00954 pp = currentParticle.GetTotalMomentum();
00955 pp1 = currentParticle.GetMomentum().mag();
00956 if (pp1 < 1.0e-6*GeV) {
00957 G4ThreeVector iso = Isotropic(pp);
00958 currentParticle.SetMomentum( iso.x(), iso.y(), iso.z() );
00959 } else {
00960 currentParticle.SetMomentum(currentParticle.GetMomentum() * (pp/pp1));
00961 }
00962
00963 theoreticalKinetic = targetParticle.GetKineticEnergy() * wgt;
00964 targetParticle.SetKineticEnergy(theoreticalKinetic);
00965 simulatedKinetic += theoreticalKinetic;
00966 pp = targetParticle.GetTotalMomentum();
00967 pp1 = targetParticle.GetMomentum().mag();
00968
00969 if (pp1 < 1.0e-6*GeV) {
00970 G4ThreeVector iso = Isotropic(pp);
00971 targetParticle.SetMomentum(iso.x(), iso.y(), iso.z() );
00972 } else {
00973 targetParticle.SetMomentum(targetParticle.GetMomentum() * (pp/pp1) );
00974 }
00975
00976 for (i = 0; i < vecLen; ++i ) {
00977 theoreticalKinetic = vec[i]->GetKineticEnergy() * wgt;
00978 simulatedKinetic += theoreticalKinetic;
00979 vec[i]->SetKineticEnergy(theoreticalKinetic);
00980 pp = vec[i]->GetTotalMomentum();
00981 pp1 = vec[i]->GetMomentum().mag();
00982 if( pp1 < 1.0e-6*GeV ) {
00983 G4ThreeVector iso = Isotropic(pp);
00984 vec[i]->SetMomentum(iso.x(), iso.y(), iso.z() );
00985 } else {
00986 vec[i]->SetMomentum(vec[i]->GetMomentum() * (pp/pp1) );
00987 }
00988 }
00989 }
00990
00991
00992
00993
00994
00995
00996
00997
00998
00999 if( atomicWeight >= 1.5 )
01000 {
01001
01002
01003
01004
01005
01006
01007 G4int npnb = 0;
01008 G4int ndta = 0;
01009
01010 G4double epnb, edta;
01011 if (veryForward) {
01012 epnb = targetNucleus.GetAnnihilationPNBlackTrackEnergy();
01013 edta = targetNucleus.GetAnnihilationDTABlackTrackEnergy();
01014 } else {
01015 epnb = targetNucleus.GetPNBlackTrackEnergy();
01016 edta = targetNucleus.GetDTABlackTrackEnergy();
01017 }
01018
01019
01020
01021
01022
01023
01024
01025
01026
01027
01028
01029
01030
01031 const G4double pnCutOff = 0.001;
01032 const G4double dtaCutOff = 0.001;
01033
01034
01035
01036
01037
01038
01039 if (epnb > pnCutOff)
01040 {
01041 npnb = G4Poisson((1.5+1.25*numberofFinalStateNucleons)*epnb/(epnb+edta));
01042 if (numberofFinalStateNucleons + npnb > atomicWeight)
01043 npnb = G4int(atomicWeight+0.00001 - numberofFinalStateNucleons);
01044 npnb = std::min( npnb, 127-vecLen );
01045 }
01046 if( edta >= dtaCutOff )
01047 {
01048 ndta = G4Poisson((1.5+1.25*numberofFinalStateNucleons)*edta/(epnb+edta));
01049 ndta = std::min( ndta, 127-vecLen );
01050 }
01051 if (npnb == 0 && ndta == 0) npnb = 1;
01052
01053 AddBlackTrackParticles(epnb, npnb, edta, ndta, modifiedOriginal,
01054 PinNucleus, NinNucleus, targetNucleus,
01055 vec, vecLen);
01056 }
01057
01058
01059
01060
01061
01062
01063
01064
01065 if( (atomicWeight >= 1.5) && (atomicWeight <= 230.0) && (ekOriginal <= 0.2) )
01066 currentParticle.SetTOF(
01067 1.0-500.0*std::exp(-ekOriginal/0.04)*std::log(G4UniformRand()) );
01068 else
01069 currentParticle.SetTOF( 1.0 );
01070 return true;
01071
01072 }
01073
01074
01075 void G4RPGFragmentation::
01076 ReduceEnergiesOfSecondaries(G4int startingIndex,
01077 G4double& forwardKinetic,
01078 G4double& backwardKinetic,
01079 G4FastVector<G4ReactionProduct,256>& vec,
01080 G4int& vecLen,
01081 G4ReactionProduct& forwardPseudoParticle,
01082 G4ReactionProduct& backwardPseudoParticle,
01083 G4double& pt)
01084 {
01085
01086
01087
01088 G4double totalEnergy;
01089 G4double pp;
01090 G4double pp1;
01091 G4double px;
01092 G4double py;
01093 G4double mass;
01094 G4ReactionProduct* pVec;
01095 G4int i;
01096
01097 forwardKinetic = 0.0;
01098 backwardKinetic = 0.0;
01099 forwardPseudoParticle.SetZero();
01100 backwardPseudoParticle.SetZero();
01101
01102 for (i = startingIndex; i < vecLen; i++) {
01103 pVec = vec[i];
01104 if (pVec->GetSide() != -3) {
01105 mass = pVec->GetMass();
01106 totalEnergy = 0.95*pVec->GetTotalEnergy() + 0.05*mass;
01107 pVec->SetTotalEnergy(totalEnergy);
01108 pp = std::sqrt( std::abs( totalEnergy*totalEnergy - mass*mass ) );
01109 pp1 = pVec->GetMomentum().mag();
01110 if (pp1 < 1.0e-6*GeV) {
01111 G4ThreeVector iso = Isotropic(pp);
01112 pVec->SetMomentum( iso.x(), iso.y(), iso.z() );
01113 } else {
01114 pVec->SetMomentum(pVec->GetMomentum() * (pp/pp1) );
01115 }
01116
01117 px = pVec->GetMomentum().x();
01118 py = pVec->GetMomentum().y();
01119 pt = std::max(1.0, std::sqrt( px*px + py*py ) )/GeV;
01120 if (pVec->GetSide() > 0) {
01121 forwardKinetic += pVec->GetKineticEnergy()/GeV;
01122 forwardPseudoParticle = forwardPseudoParticle + (*pVec);
01123 } else {
01124 backwardKinetic += pVec->GetKineticEnergy()/GeV;
01125 backwardPseudoParticle = backwardPseudoParticle + (*pVec);
01126 }
01127 }
01128 }
01129 }
01130
01131