51 G4cout <<
" G4RPGReactionStage must be overridden in a derived class "
100 G4int local_npnb = npnb;
102 local_npnb =
std::min(PinNucleus + NinNucleus , local_npnb);
104 if (ndta == 0) local_epnb += edta;
107 remainingE = local_epnb;
108 for (i = 0; i < local_npnb; ++i)
127 if (PinNucleus > 0) {
141 if (NinNucleus > 0) {
153 if (i < local_npnb - 1) {
155 remainingE -= kinetic;
157 kinetic = remainingE;
162 G4double sint = std::sqrt(std::fabs(1.0-cost*cost));
164 vec[vecLen]->SetNewlyAdded(
true );
165 vec[vecLen]->SetKineticEnergy( kinetic*
GeV );
167 pp = vec[vecLen]->GetTotalMomentum();
168 vec[vecLen]->SetMomentum(pp*sint*std::sin(phi),
169 pp*sint*std::cos(phi),
174 if (NinNucleus > 0) {
175 if( (atomicWeight >= 10.0) && (ekOriginal <= 2.0*
GeV) )
179 if( ekw > 1.0 )ekw *= ekw;
181 ika =
G4int(ika1*
G4Exp((atomicNumber*atomicNumber/
182 atomicWeight-ika2)/ika3)/ekw);
185 for( i=(vecLen-1); i>=0; --i )
187 if( (vec[i]->GetDefinition() == aProton) && vec[i]->GetNewlyAdded() )
189 vec[i]->SetDefinitionAndUpdateE( aNeutron );
192 if( ++kk > ika )
break;
205 G4int local_ndta=ndta;
208 if (npnb == 0) local_edta += epnb;
211 remainingE = local_edta;
212 for (i = 0; i < local_ndta; ++i) {
228 if (PinNucleus > 0 && NinNucleus > 0) {
232 }
else if (NinNucleus > 0) {
235 }
else if (PinNucleus > 0) {
242 }
else if (ran < 0.90) {
243 if (PinNucleus > 0 && NinNucleus > 1) {
247 }
else if (PinNucleus > 0 && NinNucleus > 0) {
251 }
else if (NinNucleus > 0) {
254 }
else if (PinNucleus > 0) {
262 if (PinNucleus > 1 && NinNucleus > 1) {
266 }
else if (PinNucleus > 0 && NinNucleus > 1) {
270 }
else if (PinNucleus > 0 && NinNucleus > 0) {
274 }
else if (NinNucleus > 0) {
277 }
else if (PinNucleus > 0) {
286 if (i < local_ndta - 1) {
288 remainingE -= kinetic;
290 kinetic = remainingE;
297 vec[vecLen]->SetNewlyAdded(
true );
298 vec[vecLen]->SetKineticEnergy( kinetic*
GeV );
301 pp = vec[vecLen]->GetTotalMomentum();
302 vec[vecLen]->SetMomentum( pp*sint*std::sin(phi),
303 pp*sint*std::cos(phi),
313 const G4bool constantCrossSection,
322 G4cerr <<
"*** Error in G4RPGReaction::GenerateNBodyEvent" <<
G4endl;
324 G4cerr <<
"totalEnergy = " << totalEnergy <<
"MeV, vecLen = " << vecLen <<
G4endl;
336 for (i=0; i<vecLen; ++i) {
337 mass[i] = vec[i]->GetMass()/
GeV;
338 if(vec[i]->GetSide() == -2) extraMass+=vec[i]->GetMass()/
GeV;
339 vec[i]->SetMomentum( 0.0, 0.0, 0.0 );
344 totalMass += mass[i];
349 if (totalMass > totalE) {
357 G4double kineticEnergy = totalE - totalMass;
360 emm[vecLen-1] = totalE;
365 for (i=0; i<vecLen-2; ++i) {
366 for (
G4int j=vecLen-2; j>i; --j) {
367 if (ran[i] > ran[j]) {
374 for( i=1; i<vecLen-1; ++i )emm[i] = ran[i-1]*kineticEnergy + sm[i];
381 if (constantCrossSection) {
382 G4double emmax = kineticEnergy + mass[0];
384 for( i=1; i<vecLen; ++i )
389 if( emmax*emmax > 0.0 )
392 + (emmin*emmin-mass[i]*mass[i])*(emmin*emmin-mass[i]*mass[i])/(emmax*emmax)
393 - 2.0*(emmin*emmin+mass[i]*mass[i]);
394 if( arg > 0.0 )wtfc = 0.5*std::sqrt( arg );
401 wtmax +=
G4Log( wtfc );
409 const G4double ffq[18] = { 0., 3.141592, 19.73921, 62.01255, 129.8788, 204.0131,
410 256.3704, 268.4705, 240.9780, 189.2637,
411 132.1308, 83.0202, 47.4210, 24.8295,
412 12.0006, 5.3858, 2.2560, 0.8859 };
413 wtmax =
G4Log( std::pow( kineticEnergy, vecLen-2 ) * ffq[vecLen-1] / totalE );
421 for( i=0; i<vecLen-1; ++i )
424 if( emm[i+1]*emm[i+1] > 0.0 )
427 + (emm[i]*emm[i]-mass[i+1]*mass[i+1])*(emm[i]*emm[i]-mass[i+1]*mass[i+1])
429 - 2.0*(emm[i]*emm[i]+mass[i+1]*mass[i+1]);
430 if( arg > 0.0 )pd[i] = 0.5*std::sqrt( arg );
435 wtmax +=
G4Log( pd[i] );
440 G4double bang, cb, sb,
s0, s1, s2, c, esys,
a, b, gama, beta;
445 for( i=1; i<vecLen; ++i )
448 pcm[1][i] = -pd[i-1];
454 ss = std::sqrt( std::fabs( 1.0-c*c ) );
457 esys = std::sqrt(pd[i]*pd[i] + emm[i]*emm[i]);
460 for(
G4int j=0; j<=i; ++j )
465 energy[j] = std::sqrt( s0*s0 + s1*s1 + s2*s2 + mass[j]*mass[j] );
467 pcm[1][j] = s0*ss + s1*c;
469 pcm[0][j] = a*cb - b*sb;
470 pcm[2][j] = a*sb + b*cb;
471 pcm[1][j] = gama*(pcm[1][j] + beta*energy[j]);
476 for(
G4int j=0; j<=i; ++j )
481 energy[j] = std::sqrt( s0*s0 + s1*s1 + s2*s2 + mass[j]*mass[j] );
483 pcm[1][j] = s0*ss + s1*c;
485 pcm[0][j] = a*cb - b*sb;
486 pcm[2][j] = a*sb + b*cb;
491 for (i=0; i<vecLen; ++i) {
492 vec[i]->SetMomentum( pcm[0][i]*
GeV, pcm[1][i]*GeV, pcm[2][i]*GeV );
493 vec[i]->SetTotalEnergy( energy[i]*GeV );
502 const G4bool constantCrossSection,
503 std::vector<G4ReactionProduct*>& tempList)
508 G4int listLen = tempList.size();
511 G4cerr <<
"*** Error in G4RPGReaction::GenerateNBodyEvent" <<
G4endl;
513 G4cerr <<
"totalEnergy = " << totalEnergy <<
"MeV, listLen = " << listLen <<
G4endl;
525 for (i=0; i<listLen; ++i) {
526 mass[i] = tempList[i]->GetMass()/
GeV;
527 if(tempList[i]->GetSide() == -2) extraMass+=tempList[i]->GetMass()/
GeV;
528 tempList[i]->SetMomentum( 0.0, 0.0, 0.0 );
533 totalMass += mass[i];
538 if (totalMass > totalE) {
543 G4double kineticEnergy = totalE - totalMass;
546 emm[listLen-1] = totalE;
551 for (i=0; i<listLen-2; ++i) {
552 for (
G4int j=listLen-2; j>i; --j) {
553 if (ran[i] > ran[j]) {
560 for( i=1; i<listLen-1; ++i )emm[i] = ran[i-1]*kineticEnergy + sm[i];
567 if (constantCrossSection) {
568 G4double emmax = kineticEnergy + mass[0];
570 for( i=1; i<listLen; ++i )
575 if( emmax*emmax > 0.0 )
578 + (emmin*emmin-mass[i]*mass[i])*(emmin*emmin-mass[i]*mass[i])/(emmax*emmax)
579 - 2.0*(emmin*emmin+mass[i]*mass[i]);
580 if( arg > 0.0 )wtfc = 0.5*std::sqrt( arg );
587 wtmax +=
G4Log( wtfc );
595 const G4double ffq[18] = { 0., 3.141592, 19.73921, 62.01255, 129.8788, 204.0131,
596 256.3704, 268.4705, 240.9780, 189.2637,
597 132.1308, 83.0202, 47.4210, 24.8295,
598 12.0006, 5.3858, 2.2560, 0.8859 };
599 wtmax =
G4Log( std::pow( kineticEnergy, listLen-2 ) * ffq[listLen-1] / totalE );
607 for( i=0; i<listLen-1; ++i )
610 if( emm[i+1]*emm[i+1] > 0.0 )
613 + (emm[i]*emm[i]-mass[i+1]*mass[i+1])*(emm[i]*emm[i]-mass[i+1]*mass[i+1])
615 - 2.0*(emm[i]*emm[i]+mass[i+1]*mass[i+1]);
616 if( arg > 0.0 )pd[i] = 0.5*std::sqrt( arg );
621 wtmax +=
G4Log( pd[i] );
626 G4double bang, cb, sb,
s0, s1, s2, c, esys,
a, b, gama, beta;
631 for( i=1; i<listLen; ++i )
634 pcm[1][i] = -pd[i-1];
640 ss = std::sqrt( std::fabs( 1.0-c*c ) );
643 esys = std::sqrt(pd[i]*pd[i] + emm[i]*emm[i]);
646 for(
G4int j=0; j<=i; ++j )
651 energy[j] = std::sqrt( s0*s0 + s1*s1 + s2*s2 + mass[j]*mass[j] );
653 pcm[1][j] = s0*ss + s1*c;
655 pcm[0][j] = a*cb - b*sb;
656 pcm[2][j] = a*sb + b*cb;
657 pcm[1][j] = gama*(pcm[1][j] + beta*energy[j]);
662 for(
G4int j=0; j<=i; ++j )
667 energy[j] = std::sqrt( s0*s0 + s1*s1 + s2*s2 + mass[j]*mass[j] );
669 pcm[1][j] = s0*ss + s1*c;
671 pcm[0][j] = a*cb - b*sb;
672 pcm[2][j] = a*sb + b*cb;
677 for (i=0; i<listLen; ++i) {
678 tempList[i]->SetMomentum(pcm[0][i]*
GeV, pcm[1][i]*GeV, pcm[2][i]*GeV);
679 tempList[i]->SetTotalEnergy(energy[i]*GeV);
707 if (pjx*pjx+pjy*pjy > 0.0) {
708 G4double cost, sint, ph, cosp, sinp, pix, piy, piz;
710 sint = std::sqrt(std::abs((1.0-cost)*(1.0+cost)));
715 if( std::abs( pjx ) > 0.001*
MeV )ph = std::atan2(pjy,pjx);
721 currentParticle.
SetMomentum((cost*cosp*pix - sinp*piy + sint*cosp*piz)*
MeV,
722 (cost*sinp*pix + cosp*piy + sint*sinp*piz)*MeV,
723 (-sint*pix + cost*piz)*MeV);
727 targetParticle.
SetMomentum((cost*cosp*pix - sinp*piy + sint*cosp*piz)*MeV,
728 (cost*sinp*pix + cosp*piy + sint*sinp*piz)*MeV,
729 (-sint*pix + cost*piz)*MeV);
731 for (
G4int i=0; i<vecLen; ++i) {
732 pix = vec[i]->GetMomentum().x()/
MeV;
733 piy = vec[i]->GetMomentum().y()/
MeV;
734 piz = vec[i]->GetMomentum().z()/
MeV;
735 vec[i]->SetMomentum((cost*cosp*pix - sinp*piy + sint*cosp*piz)*MeV,
736 (cost*sinp*pix + cosp*piy + sint*sinp*piz)*MeV,
737 (-sint*pix + cost*piz)*MeV);
744 for (
G4int i=0; i<vecLen; ++i) vec[i]->SetMomentum( -vec[i]->GetMomentum().z() );
751 const G4double numberofFinalStateNucleons,
775 for( i=0; i<4; ++i )pseudoParticle[i].set(0,0,0);
778 for( i=0; i<vecLen; ++i )
779 pseudoParticle[0] = pseudoParticle[0] + (vec[i]->GetMomentum());
785 G4double r1, r2,
a1, ran1, ran2, xxh, exh, pxTemp, pyTemp, pzTemp;
789 a1 = std::sqrt(-2.0*
G4Log(r2));
790 ran1 = a1*std::sin(r1)*0.020*numberofFinalStateNucleons*
GeV;
791 ran2 = a1*std::cos(r1)*0.020*numberofFinalStateNucleons*
GeV;
794 pseudoParticle[0] = pseudoParticle[0]+fermir;
795 pseudoParticle[2] = temp;
796 pseudoParticle[3] = pseudoParticle[0];
798 pseudoParticle[1] = pseudoParticle[2].cross(pseudoParticle[3]);
800 pseudoParticle[1] = pseudoParticle[1].rotate(rotation, pseudoParticle[3]);
801 pseudoParticle[2] = pseudoParticle[3].cross(pseudoParticle[1]);
802 for(
G4int ii=1; ii<=3; ii++)
804 p = pseudoParticle[ii].mag();
808 pseudoParticle[ii]= pseudoParticle[ii] * (1./p);
811 pxTemp = pseudoParticle[1].dot(currentParticle.
GetMomentum());
812 pyTemp = pseudoParticle[2].dot(currentParticle.
GetMomentum());
813 pzTemp = pseudoParticle[3].dot(currentParticle.
GetMomentum());
814 currentParticle.
SetMomentum( pxTemp, pyTemp, pzTemp );
816 pxTemp = pseudoParticle[1].dot(targetParticle.
GetMomentum());
817 pyTemp = pseudoParticle[2].dot(targetParticle.
GetMomentum());
818 pzTemp = pseudoParticle[3].dot(targetParticle.
GetMomentum());
819 targetParticle.
SetMomentum( pxTemp, pyTemp, pzTemp );
821 for( i=0; i<vecLen; ++i )
823 pxTemp = pseudoParticle[1].dot(vec[i]->GetMomentum());
824 pyTemp = pseudoParticle[2].dot(vec[i]->GetMomentum());
825 pzTemp = pseudoParticle[3].dot(vec[i]->GetMomentum());
826 vec[i]->SetMomentum( pxTemp, pyTemp, pzTemp );
832 Defs1( modifiedOriginal, currentParticle, targetParticle, vec, vecLen );
837 if( atomicWeight >= 1.5 )
841 const G4double alem[] = { 1.40, 2.30, 2.70, 3.00, 3.40, 4.60, 7.00 };
842 const G4double val0[] = { 0.00, 0.40, 0.48, 0.51, 0.54, 0.60, 0.65 };
845 if( alekw > alem[0] )
848 for(
G4int j=1; j<7; ++j )
850 if( alekw < alem[j] )
852 G4double rcnve = (val0[j] - val0[j-1]) / (alem[j] - alem[j-1]);
853 exh = rcnve * alekw + val0[j-1] - rcnve * alem[j-1];
859 const G4double cfa = 0.025*((atomicWeight-1.)/120.)*
G4Exp(-(atomicWeight-1.)/120.);
867 dekin += ekin*(1.0-xxh);
876 if( pp1 < 0.001*
MeV )
879 G4double sintheta = std::sqrt(1. - costheta*costheta);
882 pp*sintheta*std::sin(phi)*MeV,
894 dekin += ekin*(1.0-xxh);
903 if( pp1 < 0.001*
MeV )
906 G4double sintheta = std::sqrt(1. - costheta*costheta);
909 pp*sintheta*std::sin(phi)*MeV,
914 for( i=0; i<vecLen; ++i )
916 ekin = vec[i]->GetKineticEnergy()/GeV - cfa*(1+
normal()/2.0);
921 vec[i]->GetDefinition() == aPiZero &&
923 dekin += ekin*(1.0-xxh);
925 if (vec[i]->GetDefinition()->GetParticleSubType() ==
"pi") {
929 vec[i]->SetKineticEnergy( ekin*GeV );
930 pp = vec[i]->GetTotalMomentum()/
MeV;
931 pp1 = vec[i]->GetMomentum().mag()/
MeV;
932 if( pp1 < 0.001*
MeV )
935 G4double sintheta = std::sqrt(1. - costheta*costheta);
937 vec[i]->SetMomentum( pp*sintheta*std::cos(phi)*
MeV,
938 pp*sintheta*std::sin(phi)*MeV,
942 vec[i]->SetMomentum( vec[i]->GetMomentum() * (pp/pp1) );
945 if( (ek1 != 0.0) && (npions > 0) )
947 dekin = 1.0 + dekin/ek1;
960 G4double sintheta = std::sqrt(1. - costheta*costheta);
963 pp*sintheta*std::sin(phi)*MeV,
979 G4double sintheta = std::sqrt(1. - costheta*costheta);
982 pp*sintheta*std::sin(phi)*MeV,
989 for( i=0; i<vecLen; ++i )
991 if (vec[i]->GetDefinition()->GetParticleSubType() ==
"pi")
993 vec[i]->SetKineticEnergy(
std::max( 0.001*
MeV, dekin*vec[i]->GetKineticEnergy() ) );
994 pp = vec[i]->GetTotalMomentum()/
MeV;
995 pp1 = vec[i]->GetMomentum().mag()/
MeV;
999 G4double sintheta = std::sqrt(1. - costheta*costheta);
1001 vec[i]->SetMomentum( pp*sintheta*std::cos(phi)*
MeV,
1002 pp*sintheta*std::sin(phi)*MeV,
1005 vec[i]->SetMomentum( vec[i]->GetMomentum() * (pp/pp1) );
1015 const G4int& vecLen)
1019 G4int protonsRemoved = 0;
1020 G4int neutronsRemoved = 0;
1027 for (
G4int i = 0; i < vecLen; i++) {
1028 secName = vec[i]->GetDefinition()->GetParticleName();
1029 if (secName ==
"proton") {
1031 }
else if (secName ==
"neutron") {
1033 }
else if (secName ==
"anti_proton") {
1035 }
else if (secName ==
"anti_neutron") {
1040 return std::pair<G4int, G4int>(protonsRemoved, neutronsRemoved);
1047 G4double sintheta = std::sqrt(1. - costheta*costheta);
1050 pp*sintheta*std::sin(phi),
1065 if( testMomentum >= pOriginal )
1069 std::sqrt( pMass*pMass + pOriginal*pOriginal )*
MeV );
1071 currentParticle.
GetMomentum() * (pOriginal/testMomentum) );
1074 if( testMomentum >= pOriginal )
1078 std::sqrt( pMass*pMass + pOriginal*pOriginal )*
MeV );
1080 targetParticle.
GetMomentum() * (pOriginal/testMomentum) );
1082 for(
G4int i=0; i<vecLen; ++i )
1084 testMomentum = vec[i]->GetMomentum().mag()/
MeV;
1085 if( testMomentum >= pOriginal )
1087 pMass = vec[i]->GetMass()/
MeV;
1088 vec[i]->SetTotalEnergy(
1089 std::sqrt( pMass*pMass + pOriginal*pOriginal )*
MeV );
1090 vec[i]->SetMomentum( vec[i]->GetMomentum() * (pOriginal/testMomentum) );
1121 currentParticle = *originalIncident;
1129 if( pp <= 0.001*
MeV )
1133 currentParticle.
SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
1134 p*std::sin(rthnve)*std::sin(phinve),
1135 p*std::cos(rthnve) );
1144 G4double qv = currentKinetic + theAtomicMass + currentMass;
1147 qval[0] = qv - mass[0];
1148 qval[1] = qv - mass[1] - aNeutronMass;
1149 qval[2] = qv - mass[2] - aProtonMass;
1150 qval[3] = qv - mass[3] - aDeuteronMass;
1151 qval[4] = qv - mass[4] - aTritonMass;
1152 qval[5] = qv - mass[5] - anAlphaMass;
1153 qval[6] = qv - mass[6] - aNeutronMass - aNeutronMass;
1154 qval[7] = qv - mass[7] - aNeutronMass - aProtonMass;
1155 qval[8] = qv - mass[8] - aProtonMass - aProtonMass;
1163 qval[2] = qval[3] = qval[4] = qval[5] = qval[8] = 0.0;
1170 for( i=0; i<9; ++i )
1172 if( mass[i] < 500.0*
MeV )qval[i] = 0.0;
1173 if( qval[i] < 0.0 )qval[i] = 0.0;
1179 for( index=0; index<9; ++index )
1181 if( qval[index] > 0.0 )
1183 qv1 += qval[index]/qv;
1184 if( ran <= qv1 )
break;
1190 "G4RPGReaction::NuclearReaction: inelastic reaction kinematically not possible");
1245 pseudo1.
SetMass( theAtomicMass*MeV );
1257 if( nt == 3 )tempV.
SetElement( tempLen++, v[2] );
1258 G4bool constantCrossSection =
true;
1260 v[0]->
Lorentz( *v[0], pseudo2 );
1261 v[1]->
Lorentz( *v[1], pseudo2 );
1262 if( nt == 3 )v[2]->
Lorentz( *v[2], pseudo2 );
1264 G4bool particleIsDefined =
false;
1265 if( v[0]->GetMass()/MeV - aProtonMass < 0.1 )
1268 particleIsDefined =
true;
1270 else if( v[0]->GetMass()/MeV - aNeutronMass < 0.1 )
1273 particleIsDefined =
true;
1275 else if( v[0]->GetMass()/MeV - aDeuteronMass < 0.1 )
1278 particleIsDefined =
true;
1280 else if( v[0]->GetMass()/MeV - aTritonMass < 0.1 )
1283 particleIsDefined =
true;
1285 else if( v[0]->GetMass()/MeV - anAlphaMass < 0.1 )
1288 particleIsDefined =
true;
1294 if( pp <= 0.001*MeV )
1298 currentParticle.
SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
1299 p*std::sin(rthnve)*std::sin(phinve),
1300 p*std::cos(rthnve) );
1305 if( particleIsDefined )
1311 if( pp <= 0.001*MeV )
1315 v[0]->
SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
1316 p*std::sin(rthnve)*std::sin(phinve),
1317 p*std::cos(rthnve) );
1320 v[0]->
SetMomentum( v[0]->GetMomentum() * (p/pp) );
1322 if( (v[1]->GetDefinition() == aDeuteron) ||
1323 (v[1]->GetDefinition() == aTriton) ||
1324 (v[1]->GetDefinition() == anAlpha) )
1332 if( pp <= 0.001*MeV )
1336 v[1]->
SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
1337 p*std::sin(rthnve)*std::sin(phinve),
1338 p*std::cos(rthnve) );
1341 v[1]->
SetMomentum( v[1]->GetMomentum() * (p/pp) );
1345 if( (v[2]->GetDefinition() == aDeuteron) ||
1346 (v[2]->GetDefinition() == aTriton) ||
1347 (v[2]->GetDefinition() == anAlpha) )
1355 if( pp <= 0.001*MeV )
1359 v[2]->
SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
1360 p*std::sin(rthnve)*std::sin(phinve),
1361 p*std::cos(rthnve) );
1364 v[2]->
SetMomentum( v[2]->GetMomentum() * (p/pp) );
1367 for(del=0; del<vecLen; del++)
delete vec[del];
1369 if( particleIsDefined )
void SetElement(G4int anIndex, Type *anElement)
void Rotate(const G4double numberofFinalStateNucleons, const G4ThreeVector &temp, const G4ReactionProduct &modifiedOriginal, const G4HadProjectile *originalIncident, const G4Nucleus &targetNucleus, G4ReactionProduct ¤tParticle, G4ReactionProduct &targetParticle, G4FastVector< G4ReactionProduct, 256 > &vec, G4int &vecLen)
void MomentumCheck(const G4ReactionProduct &modifiedOriginal, G4ReactionProduct ¤tParticle, G4ReactionProduct &targetParticle, G4FastVector< G4ReactionProduct, 256 > &vec, G4int &vecLen)
static const double halfpi
G4double GetTotalMomentum() const
void Lorentz(const G4ReactionProduct &p1, const G4ReactionProduct &p2)
void Defs1(const G4ReactionProduct &modifiedOriginal, G4ReactionProduct ¤tParticle, G4ReactionProduct &targetParticle, G4FastVector< G4ReactionProduct, 256 > &vec, G4int &vecLen)
CLHEP::Hep3Vector G4ThreeVector
void SetKineticEnergy(const G4double en)
void SetMomentum(const G4double x, const G4double y, const G4double z)
G4bool ReactionStage(const G4HadProjectile *, G4ReactionProduct &, G4bool &, const G4DynamicParticle *, G4ReactionProduct &, G4bool &, const G4Nucleus &, G4ReactionProduct &, G4FastVector< G4ReactionProduct, 256 > &, G4int &, G4bool, G4ReactionProduct &)
G4ParticleDefinition * GetDefinition() const
const G4String & GetParticleSubType() const
void Initialize(G4int items)
const G4String & GetParticleName() const
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
std::pair< G4int, G4int > GetFinalStateNucleons(const G4DynamicParticle *originalTarget, const G4FastVector< G4ReactionProduct, 256 > &vec, const G4int &vecLen)
const G4ParticleDefinition * GetDefinition() const
void SetMass(const G4double mas)
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
G4double GetKineticEnergy() const
void SetTotalEnergy(const G4double en)
static const double twopi
static G4Triton * Triton()
static G4Proton * Proton()
static G4PionPlus * PionPlus()
static G4Neutron * Neutron()
static G4PionZero * PionZero()
static G4Deuteron * Deuteron()
G4double GetKineticEnergy() const
G4double G4Log(G4double x)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
G4double GenerateNBodyEventT(const G4double totalEnergy, const G4bool constantCrossSection, std::vector< G4ReactionProduct * > &list)
G4double GetPDGMass() const
void NuclearReaction(G4FastVector< G4ReactionProduct, 4 > &vec, G4int &vecLen, const G4HadProjectile *originalIncident, const G4Nucleus &aNucleus, const G4double theAtomicMass, const G4double *massVec)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double energy(const ThreeVector &p, const G4double m)
static G4PionMinus * PionMinus()
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4ThreeVector GetMomentum() const
void AddBlackTrackParticles(const G4double, const G4int, const G4double, const G4int, const G4ReactionProduct &, G4int, G4int, const G4Nucleus &, G4FastVector< G4ReactionProduct, 256 > &, G4int &)
G4double GenerateNBodyEvent(const G4double totalEnergy, const G4bool constantCrossSection, G4FastVector< G4ReactionProduct, 256 > &vec, G4int &vecLen)
G4ThreeVector Isotropic(const G4double &)
G4GLOB_DLL std::ostream G4cerr