45 for (
G4int i = 0; i < 20; i++)
dndl[i] = 0.0;
58 for (
G4int i = 1; i < 20; i++) {
60 term1 = 1. + parMass*parMass*x*
x;
61 term2 = pt*x*et*pt*x*et + pt*pt + secMass*secMass;
62 dndl[i] = dx / std::sqrt( term1*term1*term1*term2 )
71 G4bool& incidentHasChanged,
96 if (vecLen == 0)
return false;
105 G4bool veryForward =
false;
112 G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
113 targetMass*targetMass +
114 2.0*targetMass*etOriginal );
121 for (i=0; i<vecLen; ++i) {
124 *vec[itemp] = *vec[i];
128 if (currentMass == 0.0 && targetMass == 0.0) {
134 currentParticle = *vec[0];
136 targetParticle = *vec[1];
138 for( i=0; i<(vecLen-2); ++i )*vec[i] = *vec[i+2];
141 temp = vec[vecLen-2];
146 incidentHasChanged =
true;
147 targetHasChanged =
true;
159 currentParticle = targetParticle;
160 targetParticle = temp;
161 incidentHasChanged =
true;
162 targetHasChanged =
true;
167 0.312+0.200*
G4Log(
G4Log(centerofmassEnergy*centerofmassEnergy))+
168 std::pow(centerofmassEnergy*centerofmassEnergy,1.5)/6000.0 );
170 G4double freeEnergy = centerofmassEnergy-currentMass-targetMass;
171 G4double forwardEnergy = freeEnergy/2.;
172 G4int forwardCount = 1;
174 G4double backwardEnergy = freeEnergy/2.;
175 G4int backwardCount = 1;
178 if(currentParticle.
GetSide()==-1)
180 forwardEnergy += currentMass;
182 backwardEnergy -= currentMass;
185 if(targetParticle.
GetSide()!=-1)
187 backwardEnergy += targetMass;
189 forwardEnergy -= targetMass;
194 for (i=0; i<vecLen; ++i) {
195 if( vec[i]->GetSide() == -1 )
198 backwardEnergy -= vec[i]->GetMass()/
GeV;
201 forwardEnergy -= vec[i]->GetMass()/
GeV;
209 if (backwardEnergy < 0.0) {
210 for (i = 0; i < vecLen; ++i) {
211 if (vec[i]->GetSide() == -1) {
212 backwardEnergy += vec[i]->GetMass()/
GeV;
215 forwardEnergy -= vec[i]->GetMass()/
GeV;
217 if (backwardEnergy > 0.0)
break;
222 if (forwardEnergy < 0.0) {
223 for (i = 0; i < vecLen; ++i) {
224 if (vec[i]->GetSide() == 1) {
225 forwardEnergy += vec[i]->GetMass()/
GeV;
228 backwardEnergy -= vec[i]->GetMass()/
GeV;
230 if (forwardEnergy > 0.0)
break;
238 if (forwardEnergy > 0.0 && backwardEnergy < 0.0) {
239 forwardEnergy += backwardEnergy;
244 if (forwardEnergy + backwardEnergy < 0.0)
return false;
262 xtarg = afc * (a13-1.0) * (2.0*backwardCount+vecLen+2)/2.0;
264 xtarg = afc * (a13-1.0) * (2.0*backwardCount);
266 if( xtarg <= 0.0 )xtarg = 0.01;
270 if(atomicWeight<1.0001) nuclearExcitationCount = 0;
271 G4int extraNucleonCount = 0;
274 if (nuclearExcitationCount > 0) {
275 const G4double nucsup[] = { 1.00, 0.7, 0.5, 0.4, 0.35, 0.3 };
276 const G4double psup[] = { 3., 6., 20., 50., 100., 1000. };
277 G4int momentumBin = 0;
281 ed <<
" While count exceeded " <<
G4endl;
282 while( (momentumBin < 6) &&
292 momentumBin =
std::min( 5, momentumBin );
298 for (i = 0; i < nuclearExcitationCount; ++i) {
322 else if( ran < 0.6819 )
345 for (i = 0; i < 8; ++i) pseudoParticle[i].SetZero();
348 pseudoParticle[0].
SetMomentum( 0.0, 0.0, pOriginal*GeV );
350 std::sqrt( pOriginal*pOriginal + mOriginal*mOriginal )*GeV );
352 pseudoParticle[1].
SetMass(protonMass);
355 pseudoParticle[3].
SetMass(protonMass*(1+extraNucleonCount) );
356 pseudoParticle[3].
SetTotalEnergy(protonMass*(1+extraNucleonCount) );
358 pseudoParticle[2] = pseudoParticle[0] + pseudoParticle[1];
359 pseudoParticle[3] = pseudoParticle[3] + pseudoParticle[0];
361 pseudoParticle[0].
Lorentz( pseudoParticle[0], pseudoParticle[2] );
362 pseudoParticle[1].Lorentz( pseudoParticle[1], pseudoParticle[2] );
368 G4int innerCounter, outerCounter;
369 G4bool eliminateThisParticle, resetEnergies, constantCrossSection;
379 G4int backwardNucleonCount = 0;
380 G4double totalEnergy, kineticEnergy, vecMass;
383 for (i = vecLen-1; i >= 0; --i) {
385 if (vec[i]->GetNewlyAdded()) {
386 if (vec[i]->GetSide() == -2) {
387 if (backwardNucleonCount < 18) {
388 if (vec[i]->GetDefinition()->GetParticleSubType() ==
"pi") {
389 for (
G4int j = 0; j < vecLen; j++)
delete vec[j];
392 "G4RPGFragmentation::ReactionStage : a pion has been counted as a backward nucleon");
395 ++backwardNucleonCount;
406 vecMass = vec[i]->GetMass()/
GeV;
409 if (vec[i]->GetSide() == -2) {
411 pt = std::sqrt( std::pow( ran, 1.2 ) );
414 if (vec[i]->GetDefinition()->GetParticleSubType() ==
"pi") {
416 pt = std::sqrt( std::pow( ran, 1.7 ) );
417 }
else if (vec[i]->GetDefinition()->GetParticleSubType() ==
"kaon") {
419 pt = std::sqrt( std::pow( ran, 1.7 ) );
422 pt = std::sqrt( std::pow( ran, 1.5 ) );
428 vec[i]->SetMomentum( pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV );
429 if (vec[i]->GetSide() > 0)
430 et = pseudoParticle[0].GetTotalEnergy()/
GeV;
432 et = pseudoParticle[1].GetTotalEnergy()/
GeV;
438 eliminateThisParticle =
true;
439 resetEnergies =
true;
442 while (++outerCounter < 3) {
446 vec[i]->SetMomentum( pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV );
450 while (++innerCounter < 7) {
457 ed <<
" While count exceeded " <<
G4endl;
458 while( ( ran >
dndl[l] ) && ( l < 19 ) ) {
468 if (vec[i]->GetSide() < 0) x *= -1.;
469 vec[i]->SetMomentum( x*et*GeV );
470 totalEnergy = std::sqrt( x*et*x*et + pt*pt + vecMass*vecMass );
471 vec[i]->SetTotalEnergy( totalEnergy*GeV );
472 kineticEnergy = vec[i]->GetKineticEnergy()/
GeV;
474 if (vec[i]->GetSide() > 0) {
475 if( (forwardKinetic+kineticEnergy) < 0.95*forwardEnergy ) {
478 pseudoParticle[4] = pseudoParticle[4] + (*vec[i]);
479 forwardKinetic += kineticEnergy;
481 eliminateThisParticle =
false;
482 resetEnergies =
false;
485 if( innerCounter > 5 )
break;
486 if( backwardEnergy >= vecMass )
489 forwardEnergy += vecMass;
490 backwardEnergy -= vecMass;
498 if (extraNucleonCount < 20) xxx = 0.95+0.05*extraNucleonCount/20.0;
500 if ((backwardKinetic+kineticEnergy) < xxx*backwardEnergy) {
501 pseudoParticle[5] = pseudoParticle[5] + (*vec[i]);
502 backwardKinetic += kineticEnergy;
504 eliminateThisParticle =
false;
505 resetEnergies =
false;
508 if (innerCounter > 5)
break;
509 if (forwardEnergy >= vecMass) {
511 forwardEnergy -= vecMass;
512 backwardEnergy += vecMass;
517 vec[i]->SetMomentum( momentum.
x() * 0.9, momentum.
y() * 0.9 );
530 pseudoParticle[4], pseudoParticle[5],
535 if (eliminateThisParticle && vec[i]->GetMayBeKilled()) {
538 if (vec[i]->GetSide() > 0) {
540 forwardEnergy += vecMass;
543 if (vec[i]->GetSide() == -2) {
545 extraNucleonMass -= vecMass;
547 backwardEnergy += vecMass;
551 for(
G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];
557 G4cout <<
" FALSE RETURN DUE TO ENERGY BALANCE " <<
G4endl;
565 G4double forwardKEDiff = forwardEnergy - forwardKinetic;
566 G4double backwardKEDiff = backwardEnergy - backwardKinetic;
568 if (forwardKEDiff < 0.0 || backwardKEDiff < 0.0) {
571 pseudoParticle[4], pseudoParticle[5],
574 forwardKEDiff = forwardEnergy - forwardKinetic;
575 backwardKEDiff = backwardEnergy - backwardKinetic;
576 if (backwardKEDiff < 0.0) {
577 if (forwardKEDiff + backwardKEDiff > 0.0) {
578 backwardEnergy = backwardKinetic;
579 forwardEnergy += backwardKEDiff;
580 forwardKEDiff = forwardEnergy - forwardKinetic;
581 backwardKEDiff = 0.0;
583 G4cout <<
" False return due to insufficient backward energy " <<
G4endl;
588 if (forwardKEDiff < 0.0) {
589 if (forwardKEDiff + backwardKEDiff > 0.0) {
590 forwardEnergy = forwardKinetic;
591 backwardEnergy += forwardKEDiff;
592 backwardKEDiff = backwardEnergy - backwardKinetic;
595 G4cout <<
" False return due to insufficient forward energy " <<
G4endl;
609 pt = std::sqrt( std::pow( ran/6.0, 1.7 ) );
612 pt = std::sqrt( std::pow( ran/5.0, 1.4 ) );
615 pt = std::sqrt( std::pow( ran/4.0, 1.2 ) );
619 currentParticle.
SetMomentum(pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV);
620 et = pseudoParticle[0].GetTotalEnergy()/
GeV;
631 ed <<
" While count exceeded " <<
G4endl;
632 while( ( ran >
dndl[l] ) && ( l < 19 ) ) {
644 if (forwardEnergy < forwardKinetic) {
645 totalEnergy = vecMass + 0.04*std::fabs(
normal());
646 G4cout <<
" Not enough forward energy: forwardEnergy = "
647 << forwardEnergy <<
" forwardKinetic = "
648 << forwardKinetic <<
" total energy left = "
649 << backwardKEDiff + forwardKEDiff <<
G4endl;
651 totalEnergy = vecMass + forwardEnergy - forwardKinetic;
652 forwardKinetic = forwardEnergy;
655 pp = std::sqrt(
std::abs( totalEnergy*totalEnergy - vecMass*vecMass) )*
GeV;
658 if (pp1 < 1.0
e-6*GeV) {
664 pseudoParticle[4] = pseudoParticle[4] + currentParticle;
670 if (backwardNucleonCount < 18) {
672 ++backwardNucleonCount;
681 pt =
std::max( 0.001, std::sqrt( std::pow( ran/4.0, 1.2 ) ) );
683 targetParticle.
SetMomentum(pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV);
684 et = pseudoParticle[1].GetTotalEnergy()/
GeV;
687 G4bool marginalEnergy =
true;
690 if( extraNucleonCount < 20 ) xxx = 0.95+0.05*extraNucleonCount/20.0;
693 while (++outerCounter < 4) {
696 for (innerCounter = 0; innerCounter < 6; innerCounter++) {
702 eda <<
" While count exceeded " <<
G4endl;
703 while( ( ran >
dndl[l] ) && ( l < 19 ) ) {
714 totalEnergy = std::sqrt(x*et*x*et + pt*pt + vecMass*vecMass);
717 if ((backwardKinetic+totalEnergy-vecMass) < xxx*backwardEnergy) {
718 pseudoParticle[5] = pseudoParticle[5] + targetParticle;
719 backwardKinetic += totalEnergy - vecMass;
721 marginalEnergy =
false;
725 targetParticle.
SetMomentum(momentum.
x() * 0.9, momentum.
y() * 0.9);
731 if (marginalEnergy) {
732 G4cout <<
" Extra backward kinetic energy = "
733 << 0.999*backwardEnergy - backwardKinetic <<
G4endl;
734 totalEnergy = vecMass + 0.999*backwardEnergy - backwardKinetic;
736 pp = std::sqrt(
std::abs(totalEnergy*totalEnergy - vecMass*vecMass) )*
GeV;
737 targetParticle.
SetMomentum(momentum.
x()/0.9, momentum.
y()/0.9);
740 pseudoParticle[5] = pseudoParticle[5] + targetParticle;
741 backwardKinetic = 0.999*backwardEnergy;
746 if (backwardEnergy < backwardKinetic)
747 G4cout <<
" Backward Edif = " << backwardEnergy - backwardKinetic <<
G4endl;
748 if (forwardEnergy != forwardKinetic)
749 G4cout <<
" Forward Edif = " << forwardEnergy - forwardKinetic <<
G4endl;
759 pseudoParticle[6].Lorentz( pseudoParticle[3], pseudoParticle[2] );
760 pseudoParticle[6] = pseudoParticle[6] - pseudoParticle[4];
761 pseudoParticle[6] = pseudoParticle[6] - pseudoParticle[5];
763 if (backwardNucleonCount == 1) {
768 std::min(backwardEnergy-backwardKinetic, centerofmassEnergy/2.0-protonMass/GeV);
770 if( ekin < 0.04 )ekin = 0.04 * std::fabs(
normal() );
772 totalEnergy = ekin + vecMass;
774 pp = std::sqrt(
std::abs(totalEnergy*totalEnergy - vecMass*vecMass) )*
GeV;
775 pp1 = pseudoParticle[6].GetMomentum().mag();
776 if (pp1 < 1.0
e-6*GeV) {
780 targetParticle.
SetMomentum( pseudoParticle[6].GetMomentum() * (pp/pp1));
782 pseudoParticle[5] = pseudoParticle[5] + targetParticle;
784 }
else if (backwardNucleonCount > 1) {
788 if (backwardNucleonCount < 5) tempCount = backwardNucleonCount;
792 if (targetParticle.
GetSide() == -3)
794 for (i = 0; i < vecLen; ++i)
795 if (vec[i]->GetSide() == -3) clusterMass += vec[i]->GetMass()/
GeV;
796 clusterMass += backwardEnergy - backwardKinetic;
798 totalEnergy = pseudoParticle[6].GetTotalEnergy()/
GeV;
799 pseudoParticle[6].SetMass(clusterMass*GeV);
801 pp = std::sqrt(
std::abs(totalEnergy*totalEnergy -
802 clusterMass*clusterMass) )*
GeV;
803 pp1 = pseudoParticle[6].GetMomentum().mag();
804 if (pp1 < 1.0
e-6*GeV) {
806 pseudoParticle[6].SetMomentum(iso.
x(), iso.
y(), iso.
z());
808 pseudoParticle[6].SetMomentum(pseudoParticle[6].GetMomentum() * (-pp/pp1));
811 std::vector<G4ReactionProduct*> tempList;
812 if (targetParticle.
GetSide() == -3) tempList.push_back(&targetParticle);
813 for (i = 0; i < vecLen; ++i)
814 if (vec[i]->GetSide() == -3) tempList.push_back(vec[i]);
816 constantCrossSection =
true;
818 if (tempList.size() > 1) {
821 constantCrossSection, tempList);
823 if (targetParticle.
GetSide() == -3) {
824 targetParticle = *tempList[0];
825 targetParticle.
Lorentz(targetParticle, pseudoParticle[6]);
829 for (i = 0; i < vecLen; ++i) {
830 if (vec[i]->GetSide() == -3) {
831 *vec[i] = *tempList[n_entry];
832 vec[i]->Lorentz(*vec[i], pseudoParticle[6]);
839 if (vecLen == 0)
return false;
843 currentParticle.Lorentz( currentParticle, pseudoParticle[1] );
844 targetParticle.
Lorentz( targetParticle, pseudoParticle[1] );
845 for (i = 0; i < vecLen; ++i) vec[i]->Lorentz(*vec[i], pseudoParticle[1]);
855 G4bool leadingStrangeParticleHasChanged =
true;
858 if (currentParticle.GetDefinition() == leadingStrangeParticle.
GetDefinition())
859 leadingStrangeParticleHasChanged =
false;
860 if (leadingStrangeParticleHasChanged &&
862 leadingStrangeParticleHasChanged =
false;
863 if( leadingStrangeParticleHasChanged )
865 for( i=0; i<vecLen; i++ )
867 if( vec[i]->GetDefinition() == leadingStrangeParticle.
GetDefinition() )
869 leadingStrangeParticleHasChanged =
false;
874 if( leadingStrangeParticleHasChanged )
885 if( (leadTest&&targetTest) || !(leadTest||targetTest) )
888 targetHasChanged =
true;
892 currentParticle.SetDefinitionAndUpdateE( leadingStrangeParticle.
GetDefinition() );
893 incidentHasChanged =
false;
901 std::pair<G4int, G4int> finalStateNucleons =
904 G4int protonsInFinalState = finalStateNucleons.first;
905 G4int neutronsInFinalState = finalStateNucleons.second;
907 G4int numberofFinalStateNucleons =
908 protonsInFinalState + neutronsInFinalState;
910 if (currentParticle.GetDefinition()->GetBaryonNumber() == 1 &&
914 numberofFinalStateNucleons++;
916 numberofFinalStateNucleons =
std::max(1, numberofFinalStateNucleons);
923 pseudoParticle[3].SetMomentum( 0.0, 0.0, pOriginal*GeV );
924 pseudoParticle[3].SetMass( mOriginal*GeV );
925 pseudoParticle[3].SetTotalEnergy(
926 std::sqrt( pOriginal*pOriginal + mOriginal*mOriginal )*GeV );
931 if(numberofFinalStateNucleons == 1) diff = 0;
932 pseudoParticle[4].SetMomentum( 0.0, 0.0, 0.0 );
933 pseudoParticle[4].SetMass( protonMass*(numberofFinalStateNucleons-diff) );
934 pseudoParticle[4].SetTotalEnergy( protonMass*(numberofFinalStateNucleons-diff) );
937 pseudoParticle[3].GetTotalEnergy() + pseudoParticle[4].GetTotalEnergy() -
938 currentParticle.GetMass() - targetParticle.
GetMass();
939 for (i = 0; i < vecLen; ++i) theoreticalKinetic -= vec[i]->GetMass();
943 for (i = 0; i < vecLen; ++i)
944 simulatedKinetic += vec[i]->GetKineticEnergy();
946 pseudoParticle[5] = pseudoParticle[3] + pseudoParticle[4];
947 pseudoParticle[3].Lorentz( pseudoParticle[3], pseudoParticle[5] );
948 pseudoParticle[4].Lorentz( pseudoParticle[4], pseudoParticle[5] );
950 pseudoParticle[7].SetZero();
951 pseudoParticle[7] = pseudoParticle[7] + currentParticle;
952 pseudoParticle[7] = pseudoParticle[7] + targetParticle;
953 for (i = 0; i < vecLen; ++i)
954 pseudoParticle[7] = pseudoParticle[7] + *vec[i];
999 if (simulatedKinetic != 0.0) {
1000 wgt = theoreticalKinetic/simulatedKinetic;
1001 theoreticalKinetic = currentParticle.GetKineticEnergy() * wgt;
1002 simulatedKinetic = theoreticalKinetic;
1003 currentParticle.SetKineticEnergy(theoreticalKinetic);
1004 pp = currentParticle.GetTotalMomentum();
1005 pp1 = currentParticle.GetMomentum().mag();
1006 if (pp1 < 1.0
e-6*GeV) {
1008 currentParticle.SetMomentum( iso.
x(), iso.
y(), iso.
z() );
1010 currentParticle.SetMomentum(currentParticle.GetMomentum() * (pp/pp1));
1013 theoreticalKinetic = targetParticle.GetKineticEnergy() * wgt;
1014 targetParticle.SetKineticEnergy(theoreticalKinetic);
1015 simulatedKinetic += theoreticalKinetic;
1016 pp = targetParticle.GetTotalMomentum();
1017 pp1 = targetParticle.GetMomentum().mag();
1019 if (pp1 < 1.0
e-6*GeV) {
1021 targetParticle.SetMomentum(iso.
x(), iso.
y(), iso.
z() );
1023 targetParticle.SetMomentum(targetParticle.GetMomentum() * (pp/pp1) );
1026 for (i = 0; i < vecLen; ++i ) {
1027 theoreticalKinetic = vec[i]->GetKineticEnergy() * wgt;
1028 simulatedKinetic += theoreticalKinetic;
1029 vec[i]->SetKineticEnergy(theoreticalKinetic);
1030 pp = vec[i]->GetTotalMomentum();
1031 pp1 = vec[i]->GetMomentum().mag();
1032 if( pp1 < 1.0
e-6*GeV ) {
1034 vec[i]->SetMomentum(iso.
x(), iso.
y(), iso.
z() );
1036 vec[i]->SetMomentum(vec[i]->GetMomentum() * (pp/pp1) );
1049 if( atomicWeight >= 1.5 )
1089 if (epnb > pnCutOff)
1091 npnb =
G4Poisson((1.5+1.25*numberofFinalStateNucleons)*epnb/(epnb+edta));
1092 if (numberofFinalStateNucleons + npnb > atomicWeight)
1093 npnb =
G4int(atomicWeight+0.00001 - numberofFinalStateNucleons);
1094 npnb =
std::min( npnb, 127-vecLen );
1096 if( edta >= dtaCutOff )
1098 ndta =
G4Poisson((1.5+1.25*numberofFinalStateNucleons)*edta/(epnb+edta));
1099 ndta =
std::min( ndta, 127-vecLen );
1101 if (npnb == 0 && ndta == 0) npnb = 1;
1104 PinNucleus, NinNucleus, targetNucleus,
1115 if( (atomicWeight >= 1.5) && (atomicWeight <= 230.0) && (ekOriginal <= 0.2) )
1116 currentParticle.SetTOF(
1119 currentParticle.SetTOF( 1.0 );
1147 forwardKinetic = 0.0;
1148 backwardKinetic = 0.0;
1149 forwardPseudoParticle.
SetZero();
1150 backwardPseudoParticle.
SetZero();
1152 for (i = startingIndex; i < vecLen; i++) {
1158 pp = std::sqrt(
std::abs( totalEnergy*totalEnergy - mass*mass ) );
1160 if (pp1 < 1.0
e-6*
GeV) {
1169 pt =
std::max(1.0, std::sqrt( px*px + py*py ) )/
GeV;
1172 forwardPseudoParticle = forwardPseudoParticle + (*pVec);
1175 backwardPseudoParticle = backwardPseudoParticle + (*pVec);