50 G4cout <<
" G4RPGReactionStage must be overridden in a derived class "
99 G4int local_npnb = npnb;
101 local_npnb =
std::min(PinNucleus + NinNucleus , local_npnb);
103 if (ndta == 0) local_epnb += edta;
106 remainingE = local_epnb;
107 for (i = 0; i < local_npnb; ++i)
126 if (PinNucleus > 0) {
140 if (NinNucleus > 0) {
152 if (i < local_npnb - 1) {
154 remainingE -= kinetic;
156 kinetic = remainingE;
161 G4double sint = std::sqrt(std::fabs(1.0-cost*cost));
163 vec[vecLen]->SetNewlyAdded(
true );
164 vec[vecLen]->SetKineticEnergy( kinetic*
GeV );
166 pp = vec[vecLen]->GetTotalMomentum();
167 vec[vecLen]->SetMomentum(pp*sint*std::sin(phi),
168 pp*sint*std::cos(phi),
173 if (NinNucleus > 0) {
174 if( (atomicWeight >= 10.0) && (ekOriginal <= 2.0*
GeV) )
178 if( ekw > 1.0 )ekw *= ekw;
180 ika =
G4int(ika1*std::exp((atomicNumber*atomicNumber/
181 atomicWeight-ika2)/ika3)/ekw);
184 for( i=(vecLen-1); i>=0; --i )
186 if( (vec[i]->GetDefinition() == aProton) && vec[i]->GetNewlyAdded() )
188 vec[i]->SetDefinitionAndUpdateE( aNeutron );
191 if( ++kk > ika )
break;
204 G4int local_ndta=ndta;
207 if (npnb == 0) local_edta += epnb;
210 remainingE = local_edta;
211 for (i = 0; i < local_ndta; ++i) {
227 if (PinNucleus > 0 && NinNucleus > 0) {
231 }
else if (NinNucleus > 0) {
234 }
else if (PinNucleus > 0) {
241 }
else if (ran < 0.90) {
242 if (PinNucleus > 0 && NinNucleus > 1) {
246 }
else if (PinNucleus > 0 && NinNucleus > 0) {
250 }
else if (NinNucleus > 0) {
253 }
else if (PinNucleus > 0) {
261 if (PinNucleus > 1 && NinNucleus > 1) {
265 }
else if (PinNucleus > 0 && NinNucleus > 1) {
269 }
else if (PinNucleus > 0 && NinNucleus > 0) {
273 }
else if (NinNucleus > 0) {
276 }
else if (PinNucleus > 0) {
285 if (i < local_ndta - 1) {
287 remainingE -= kinetic;
289 kinetic = remainingE;
296 vec[vecLen]->SetNewlyAdded(
true );
297 vec[vecLen]->SetKineticEnergy( kinetic*
GeV );
300 pp = vec[vecLen]->GetTotalMomentum();
301 vec[vecLen]->SetMomentum( pp*sint*std::sin(phi),
302 pp*sint*std::cos(phi),
312 const G4bool constantCrossSection,
321 G4cerr <<
"*** Error in G4RPGReaction::GenerateNBodyEvent" <<
G4endl;
323 G4cerr <<
"totalEnergy = " << totalEnergy <<
"MeV, vecLen = " << vecLen <<
G4endl;
335 for (i=0; i<vecLen; ++i) {
336 mass[i] = vec[i]->GetMass()/
GeV;
337 if(vec[i]->GetSide() == -2) extraMass+=vec[i]->GetMass()/
GeV;
338 vec[i]->SetMomentum( 0.0, 0.0, 0.0 );
343 totalMass += mass[i];
348 if (totalMass > totalE) {
356 G4double kineticEnergy = totalE - totalMass;
359 emm[vecLen-1] = totalE;
364 for (i=0; i<vecLen-2; ++i) {
365 for (
G4int j=vecLen-2; j>i; --j) {
366 if (ran[i] > ran[j]) {
373 for( i=1; i<vecLen-1; ++i )emm[i] = ran[i-1]*kineticEnergy + sm[i];
380 if (constantCrossSection) {
381 G4double emmax = kineticEnergy + mass[0];
383 for( i=1; i<vecLen; ++i )
388 if( emmax*emmax > 0.0 )
391 + (emmin*emmin-mass[i]*mass[i])*(emmin*emmin-mass[i]*mass[i])/(emmax*emmax)
392 - 2.0*(emmin*emmin+mass[i]*mass[i]);
393 if( arg > 0.0 )wtfc = 0.5*std::sqrt( arg );
400 wtmax += std::log( wtfc );
408 const G4double ffq[18] = { 0., 3.141592, 19.73921, 62.01255, 129.8788, 204.0131,
409 256.3704, 268.4705, 240.9780, 189.2637,
410 132.1308, 83.0202, 47.4210, 24.8295,
411 12.0006, 5.3858, 2.2560, 0.8859 };
412 wtmax = std::log( std::pow( kineticEnergy, vecLen-2 ) * ffq[vecLen-1] / totalE );
420 for( i=0; i<vecLen-1; ++i )
423 if( emm[i+1]*emm[i+1] > 0.0 )
426 + (emm[i]*emm[i]-mass[i+1]*mass[i+1])*(emm[i]*emm[i]-mass[i+1]*mass[i+1])
428 - 2.0*(emm[i]*emm[i]+mass[i+1]*mass[i+1]);
429 if( arg > 0.0 )pd[i] = 0.5*std::sqrt( arg );
434 wtmax += std::log( pd[i] );
439 G4double bang, cb, sb,
s0, s1, s2,
c, esys,
a,
b, gama, beta;
444 for( i=1; i<vecLen; ++i )
447 pcm[1][i] = -pd[i-1];
453 ss = std::sqrt( std::fabs( 1.0-c*c ) );
456 esys = std::sqrt(pd[i]*pd[i] + emm[i]*emm[i]);
459 for(
G4int j=0; j<=i; ++j )
464 energy[j] = std::sqrt( s0*s0 + s1*s1 + s2*s2 + mass[j]*mass[j] );
466 pcm[1][j] = s0*ss + s1*
c;
468 pcm[0][j] = a*cb - b*sb;
469 pcm[2][j] = a*sb + b*cb;
470 pcm[1][j] = gama*(pcm[1][j] + beta*energy[j]);
475 for(
G4int j=0; j<=i; ++j )
480 energy[j] = std::sqrt( s0*s0 + s1*s1 + s2*s2 + mass[j]*mass[j] );
482 pcm[1][j] = s0*ss + s1*
c;
484 pcm[0][j] = a*cb - b*sb;
485 pcm[2][j] = a*sb + b*cb;
490 for (i=0; i<vecLen; ++i) {
491 vec[i]->SetMomentum( pcm[0][i]*
GeV, pcm[1][i]*GeV, pcm[2][i]*GeV );
492 vec[i]->SetTotalEnergy( energy[i]*GeV );
501 const G4bool constantCrossSection,
502 std::vector<G4ReactionProduct*>& tempList)
507 G4int listLen = tempList.size();
510 G4cerr <<
"*** Error in G4RPGReaction::GenerateNBodyEvent" <<
G4endl;
512 G4cerr <<
"totalEnergy = " << totalEnergy <<
"MeV, listLen = " << listLen <<
G4endl;
524 for (i=0; i<listLen; ++i) {
525 mass[i] = tempList[i]->GetMass()/
GeV;
526 if(tempList[i]->GetSide() == -2) extraMass+=tempList[i]->GetMass()/
GeV;
527 tempList[i]->SetMomentum( 0.0, 0.0, 0.0 );
532 totalMass += mass[i];
537 if (totalMass > totalE) {
542 G4double kineticEnergy = totalE - totalMass;
545 emm[listLen-1] = totalE;
550 for (i=0; i<listLen-2; ++i) {
551 for (
G4int j=listLen-2; j>i; --j) {
552 if (ran[i] > ran[j]) {
559 for( i=1; i<listLen-1; ++i )emm[i] = ran[i-1]*kineticEnergy + sm[i];
566 if (constantCrossSection) {
567 G4double emmax = kineticEnergy + mass[0];
569 for( i=1; i<listLen; ++i )
574 if( emmax*emmax > 0.0 )
577 + (emmin*emmin-mass[i]*mass[i])*(emmin*emmin-mass[i]*mass[i])/(emmax*emmax)
578 - 2.0*(emmin*emmin+mass[i]*mass[i]);
579 if( arg > 0.0 )wtfc = 0.5*std::sqrt( arg );
586 wtmax += std::log( wtfc );
594 const G4double ffq[18] = { 0., 3.141592, 19.73921, 62.01255, 129.8788, 204.0131,
595 256.3704, 268.4705, 240.9780, 189.2637,
596 132.1308, 83.0202, 47.4210, 24.8295,
597 12.0006, 5.3858, 2.2560, 0.8859 };
598 wtmax = std::log( std::pow( kineticEnergy, listLen-2 ) * ffq[listLen-1] / totalE );
606 for( i=0; i<listLen-1; ++i )
609 if( emm[i+1]*emm[i+1] > 0.0 )
612 + (emm[i]*emm[i]-mass[i+1]*mass[i+1])*(emm[i]*emm[i]-mass[i+1]*mass[i+1])
614 - 2.0*(emm[i]*emm[i]+mass[i+1]*mass[i+1]);
615 if( arg > 0.0 )pd[i] = 0.5*std::sqrt( arg );
620 wtmax += std::log( pd[i] );
625 G4double bang, cb, sb,
s0, s1, s2,
c, esys,
a,
b, gama, beta;
630 for( i=1; i<listLen; ++i )
633 pcm[1][i] = -pd[i-1];
639 ss = std::sqrt( std::fabs( 1.0-c*c ) );
642 esys = std::sqrt(pd[i]*pd[i] + emm[i]*emm[i]);
645 for(
G4int j=0; j<=i; ++j )
650 energy[j] = std::sqrt( s0*s0 + s1*s1 + s2*s2 + mass[j]*mass[j] );
652 pcm[1][j] = s0*ss + s1*
c;
654 pcm[0][j] = a*cb - b*sb;
655 pcm[2][j] = a*sb + b*cb;
656 pcm[1][j] = gama*(pcm[1][j] + beta*energy[j]);
661 for(
G4int j=0; j<=i; ++j )
666 energy[j] = std::sqrt( s0*s0 + s1*s1 + s2*s2 + mass[j]*mass[j] );
668 pcm[1][j] = s0*ss + s1*
c;
670 pcm[0][j] = a*cb - b*sb;
671 pcm[2][j] = a*sb + b*cb;
676 for (i=0; i<listLen; ++i) {
677 tempList[i]->SetMomentum(pcm[0][i]*
GeV, pcm[1][i]*GeV, pcm[2][i]*GeV);
678 tempList[i]->SetTotalEnergy(energy[i]*GeV);
706 if (pjx*pjx+pjy*pjy > 0.0) {
707 G4double cost, sint, ph, cosp, sinp, pix, piy, piz;
709 sint = std::sqrt(std::abs((1.0-cost)*(1.0+cost)));
714 if( std::abs( pjx ) > 0.001*
MeV )ph = std::atan2(pjy,pjx);
720 currentParticle.
SetMomentum((cost*cosp*pix - sinp*piy + sint*cosp*piz)*
MeV,
721 (cost*sinp*pix + cosp*piy + sint*sinp*piz)*MeV,
722 (-sint*pix + cost*piz)*MeV);
726 targetParticle.
SetMomentum((cost*cosp*pix - sinp*piy + sint*cosp*piz)*MeV,
727 (cost*sinp*pix + cosp*piy + sint*sinp*piz)*MeV,
728 (-sint*pix + cost*piz)*MeV);
730 for (
G4int i=0; i<vecLen; ++i) {
731 pix = vec[i]->GetMomentum().x()/
MeV;
732 piy = vec[i]->GetMomentum().y()/
MeV;
733 piz = vec[i]->GetMomentum().z()/
MeV;
734 vec[i]->SetMomentum((cost*cosp*pix - sinp*piy + sint*cosp*piz)*MeV,
735 (cost*sinp*pix + cosp*piy + sint*sinp*piz)*MeV,
736 (-sint*pix + cost*piz)*MeV);
743 for (
G4int i=0; i<vecLen; ++i) vec[i]->SetMomentum( -vec[i]->GetMomentum().z() );
750 const G4double numberofFinalStateNucleons,
766 const G4double logWeight = std::log(atomicWeight);
774 for( i=0; i<4; ++i )pseudoParticle[i].set(0,0,0);
777 for( i=0; i<vecLen; ++i )
778 pseudoParticle[0] = pseudoParticle[0] + (vec[i]->GetMomentum());
784 G4double r1, r2, a1, ran1, ran2, xxh, exh, pxTemp, pyTemp, pzTemp;
788 a1 = std::sqrt(-2.0*std::log(r2));
789 ran1 = a1*std::sin(r1)*0.020*numberofFinalStateNucleons*
GeV;
790 ran2 = a1*std::cos(r1)*0.020*numberofFinalStateNucleons*
GeV;
793 pseudoParticle[0] = pseudoParticle[0]+fermir;
794 pseudoParticle[2] = temp;
795 pseudoParticle[3] = pseudoParticle[0];
797 pseudoParticle[1] = pseudoParticle[2].
cross(pseudoParticle[3]);
799 pseudoParticle[1] = pseudoParticle[1].
rotate(rotation, pseudoParticle[3]);
800 pseudoParticle[2] = pseudoParticle[3].
cross(pseudoParticle[1]);
801 for(
G4int ii=1; ii<=3; ii++)
803 p = pseudoParticle[ii].
mag();
807 pseudoParticle[ii]= pseudoParticle[ii] * (1./
p);
813 currentParticle.
SetMomentum( pxTemp, pyTemp, pzTemp );
818 targetParticle.
SetMomentum( pxTemp, pyTemp, pzTemp );
820 for( i=0; i<vecLen; ++i )
822 pxTemp = pseudoParticle[1].
dot(vec[i]->GetMomentum());
823 pyTemp = pseudoParticle[2].
dot(vec[i]->GetMomentum());
824 pzTemp = pseudoParticle[3].
dot(vec[i]->GetMomentum());
825 vec[i]->SetMomentum( pxTemp, pyTemp, pzTemp );
831 Defs1( modifiedOriginal, currentParticle, targetParticle, vec, vecLen );
836 if( atomicWeight >= 1.5 )
840 const G4double alem[] = { 1.40, 2.30, 2.70, 3.00, 3.40, 4.60, 7.00 };
841 const G4double val0[] = { 0.00, 0.40, 0.48, 0.51, 0.54, 0.60, 0.65 };
844 if( alekw > alem[0] )
847 for(
G4int j=1; j<7; ++j )
849 if( alekw < alem[j] )
851 G4double rcnve = (val0[j] - val0[j-1]) / (alem[j] - alem[j-1]);
852 exh = rcnve * alekw + val0[j-1] - rcnve * alem[j-1];
858 const G4double cfa = 0.025*((atomicWeight-1.)/120.)*std::exp(-(atomicWeight-1.)/120.);
866 dekin += ekin*(1.0-xxh);
875 if( pp1 < 0.001*
MeV )
878 G4double sintheta = std::sqrt(1. - costheta*costheta);
881 pp*sintheta*std::sin(phi)*MeV,
893 dekin += ekin*(1.0-xxh);
902 if( pp1 < 0.001*
MeV )
905 G4double sintheta = std::sqrt(1. - costheta*costheta);
908 pp*sintheta*std::sin(phi)*MeV,
913 for( i=0; i<vecLen; ++i )
915 ekin = vec[i]->GetKineticEnergy()/GeV - cfa*(1+
normal()/2.0);
920 vec[i]->GetDefinition() == aPiZero &&
922 dekin += ekin*(1.0-xxh);
924 if (vec[i]->GetDefinition()->GetParticleSubType() ==
"pi") {
928 vec[i]->SetKineticEnergy( ekin*GeV );
929 pp = vec[i]->GetTotalMomentum()/
MeV;
930 pp1 = vec[i]->GetMomentum().mag()/
MeV;
931 if( pp1 < 0.001*
MeV )
934 G4double sintheta = std::sqrt(1. - costheta*costheta);
936 vec[i]->SetMomentum( pp*sintheta*std::cos(phi)*
MeV,
937 pp*sintheta*std::sin(phi)*MeV,
941 vec[i]->SetMomentum( vec[i]->GetMomentum() * (pp/pp1) );
944 if( (ek1 != 0.0) && (npions > 0) )
946 dekin = 1.0 + dekin/ek1;
959 G4double sintheta = std::sqrt(1. - costheta*costheta);
962 pp*sintheta*std::sin(phi)*MeV,
978 G4double sintheta = std::sqrt(1. - costheta*costheta);
981 pp*sintheta*std::sin(phi)*MeV,
988 for( i=0; i<vecLen; ++i )
990 if (vec[i]->GetDefinition()->GetParticleSubType() ==
"pi")
992 vec[i]->SetKineticEnergy(
std::max( 0.001*
MeV, dekin*vec[i]->GetKineticEnergy() ) );
993 pp = vec[i]->GetTotalMomentum()/
MeV;
994 pp1 = vec[i]->GetMomentum().mag()/
MeV;
998 G4double sintheta = std::sqrt(1. - costheta*costheta);
1000 vec[i]->SetMomentum( pp*sintheta*std::cos(phi)*
MeV,
1001 pp*sintheta*std::sin(phi)*MeV,
1004 vec[i]->SetMomentum( vec[i]->GetMomentum() * (pp/pp1) );
1014 const G4int& vecLen)
1018 G4int protonsRemoved = 0;
1019 G4int neutronsRemoved = 0;
1026 for (
G4int i = 0; i < vecLen; i++) {
1027 secName = vec[i]->GetDefinition()->GetParticleName();
1028 if (secName ==
"proton") {
1030 }
else if (secName ==
"neutron") {
1032 }
else if (secName ==
"anti_proton") {
1034 }
else if (secName ==
"anti_neutron") {
1039 return std::pair<G4int, G4int>(protonsRemoved, neutronsRemoved);
1046 G4double sintheta = std::sqrt(1. - costheta*costheta);
1049 pp*sintheta*std::sin(phi),
1064 if( testMomentum >= pOriginal )
1068 std::sqrt( pMass*pMass + pOriginal*pOriginal )*
MeV );
1070 currentParticle.
GetMomentum() * (pOriginal/testMomentum) );
1073 if( testMomentum >= pOriginal )
1077 std::sqrt( pMass*pMass + pOriginal*pOriginal )*
MeV );
1079 targetParticle.
GetMomentum() * (pOriginal/testMomentum) );
1081 for(
G4int i=0; i<vecLen; ++i )
1083 testMomentum = vec[i]->GetMomentum().mag()/
MeV;
1084 if( testMomentum >= pOriginal )
1086 pMass = vec[i]->GetMass()/
MeV;
1087 vec[i]->SetTotalEnergy(
1088 std::sqrt( pMass*pMass + pOriginal*pOriginal )*
MeV );
1089 vec[i]->SetMomentum( vec[i]->GetMomentum() * (pOriginal/testMomentum) );
1120 currentParticle = *originalIncident;
1128 if( pp <= 0.001*
MeV )
1132 currentParticle.
SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
1133 p*std::sin(rthnve)*std::sin(phinve),
1134 p*std::cos(rthnve) );
1143 G4double qv = currentKinetic + theAtomicMass + currentMass;
1146 qval[0] = qv - mass[0];
1147 qval[1] = qv - mass[1] - aNeutronMass;
1148 qval[2] = qv - mass[2] - aProtonMass;
1149 qval[3] = qv - mass[3] - aDeuteronMass;
1150 qval[4] = qv - mass[4] - aTritonMass;
1151 qval[5] = qv - mass[5] - anAlphaMass;
1152 qval[6] = qv - mass[6] - aNeutronMass - aNeutronMass;
1153 qval[7] = qv - mass[7] - aNeutronMass - aProtonMass;
1154 qval[8] = qv - mass[8] - aProtonMass - aProtonMass;
1162 qval[2] = qval[3] = qval[4] = qval[5] = qval[8] = 0.0;
1169 for( i=0; i<9; ++i )
1171 if( mass[i] < 500.0*
MeV )qval[i] = 0.0;
1172 if( qval[i] < 0.0 )qval[i] = 0.0;
1178 for( index=0; index<9; ++
index )
1180 if( qval[index] > 0.0 )
1182 qv1 += qval[
index]/qv;
1183 if( ran <= qv1 )
break;
1189 "G4RPGReaction::NuclearReaction: inelastic reaction kinematically not possible");
1244 pseudo1.
SetMass( theAtomicMass*MeV );
1256 if( nt == 3 )tempV.
SetElement( tempLen++, v[2] );
1257 G4bool constantCrossSection =
true;
1259 v[0]->
Lorentz( *v[0], pseudo2 );
1260 v[1]->
Lorentz( *v[1], pseudo2 );
1261 if( nt == 3 )v[2]->
Lorentz( *v[2], pseudo2 );
1263 G4bool particleIsDefined =
false;
1264 if( v[0]->GetMass()/MeV - aProtonMass < 0.1 )
1267 particleIsDefined =
true;
1269 else if( v[0]->GetMass()/MeV - aNeutronMass < 0.1 )
1272 particleIsDefined =
true;
1274 else if( v[0]->GetMass()/MeV - aDeuteronMass < 0.1 )
1277 particleIsDefined =
true;
1279 else if( v[0]->GetMass()/MeV - aTritonMass < 0.1 )
1282 particleIsDefined =
true;
1284 else if( v[0]->GetMass()/MeV - anAlphaMass < 0.1 )
1287 particleIsDefined =
true;
1293 if( pp <= 0.001*MeV )
1297 currentParticle.
SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
1298 p*std::sin(rthnve)*std::sin(phinve),
1299 p*std::cos(rthnve) );
1304 if( particleIsDefined )
1310 if( pp <= 0.001*MeV )
1314 v[0]->
SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
1315 p*std::sin(rthnve)*std::sin(phinve),
1316 p*std::cos(rthnve) );
1319 v[0]->
SetMomentum( v[0]->GetMomentum() * (p/pp) );
1321 if( (v[1]->GetDefinition() == aDeuteron) ||
1322 (v[1]->GetDefinition() == aTriton) ||
1323 (v[1]->GetDefinition() == anAlpha) )
1331 if( pp <= 0.001*MeV )
1335 v[1]->
SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
1336 p*std::sin(rthnve)*std::sin(phinve),
1337 p*std::cos(rthnve) );
1340 v[1]->
SetMomentum( v[1]->GetMomentum() * (p/pp) );
1344 if( (v[2]->GetDefinition() == aDeuteron) ||
1345 (v[2]->GetDefinition() == aTriton) ||
1346 (v[2]->GetDefinition() == anAlpha) )
1354 if( pp <= 0.001*MeV )
1358 v[2]->
SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
1359 p*std::sin(rthnve)*std::sin(phinve),
1360 p*std::cos(rthnve) );
1363 v[2]->
SetMomentum( v[2]->GetMomentum() * (p/pp) );
1366 for(del=0; del<vecLen; del++)
delete vec[del];
1368 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)
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
double dot(const Hep3Vector &) const
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
G4ParticleDefinition * GetDefinition() const
std::pair< G4int, G4int > GetFinalStateNucleons(const G4DynamicParticle *originalTarget, const G4FastVector< G4ReactionProduct, 256 > &vec, const G4int &vecLen)
void SetMass(const G4double mas)
G4GLOB_DLL std::ostream G4cout
G4double GetKineticEnergy() const
void SetTotalEnergy(const G4double en)
static G4Triton * Triton()
static G4Proton * Proton()
static G4PionPlus * PionPlus()
static G4Neutron * Neutron()
static G4PionZero * PionZero()
static G4Deuteron * Deuteron()
G4double GetKineticEnergy() const
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
static G4PionMinus * PionMinus()
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
void SetDefinition(G4ParticleDefinition *aParticleDefinition)
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)
Hep3Vector cross(const Hep3Vector &) const
G4ThreeVector Isotropic(const G4double &)
Hep3Vector & rotate(double, const Hep3Vector &)
G4GLOB_DLL std::ostream G4cerr