82 {2190., 1920., 1700., 1600., 1440., 1232. };
96115.603, 133.424, 153.991, 177.729, 205.126, 236.746, 273.24, 315.361, 363.973, 420.08, 484.836, 559.573, 645.832,
97745.387, 860.289, 992.903, 1145.96, 1322.61, 1526.49, 1761.8, 2033.38, 2346.83, 2708.59, 3126.12, 3608.02, 4164.19,
984806.1, 5546.97, 6402.04, 7388.91, 8527.92, 9842.5, 11359.7, 13110.8, 15131.9, 17464.5, 20156.6, 23263.8, 26849.9,
9930988.8, 35765.7, 41279, 47642.2, 54986.3, 63462.4, 73245.2, 84536, 97567.2, 112607, 129966 };
126 fMu = 105.6583745*MeV;
127 fMpi = 139.57018*MeV;
128 fM1 = 939.5654133*MeV;
129 fM2 = 938.2720813*MeV;
198 outFile <<
"G4NeutrinoNucleusModel is a neutrino-nucleus general\n"
199 <<
"model which uses the standard model \n"
200 <<
"transfer parameterization. The model is fully relativistic\n";
227 G4int i(0), nBin(50);
230 for( i = 0; i < nBin; ++i )
242 xx =
GetXkr( nBin-1, prob);
255 else xx = x1 + (e-e1)*(x2-x1)/(e2-e1);
269 for( i = 0; i < nBin; ++i )
297 else xx = x1 + (prob-p1)*(x2-x1)/(p2-p1);
316 qq1 =
GetQkr( 0, jX, prob);
318 else if ( iE >= nBin-1)
320 qq1 =
GetQkr( nBin-1, jX, prob);
332 else qq1 = q1 + (e-e1)*(q2-q1)/(e2-e1);
339 qq2 =
GetQkr( iE, 0, prob);
341 else if ( jX >= nBin)
343 qq2 =
GetQkr( iE, nBin, prob);
355 else qq2 = q1 + (e-e1)*(q2-q1)/(e2-e1);
371 for( i = 0; i < nBin; ++i )
399 else qq = q1 + (prob-p1)*(q2-q1)/(p2-p1);
415 if( pdg == 211 || pdg == -211 || pdg == 111)
430 for(
unsigned int i = 0; i < ddktv->size(); i++ )
434 ddktv->operator[](i)->Get4Momentum());
439 delete ddktv->operator[](i);
458 if( pdg == 2212 || pdg == 2112)
469 G4double det(0.), det2(0.), rM(0.), mX = lvB.
m();
484 rM = electron_mass_c2;
490 G4double a = 4.*(sumE*sumE - pX*pX);
494 if( det2 <= 0. ) det = 0.;
495 else det = sqrt(det2);
506 eX = sqrt( pX*pX +
fMr*
fMr );
510 if( pdg == 2212 || pdg == 2112)
525 for(
unsigned int i = 0; i < ddktv->size(); i++ )
529 ddktv->operator[](i)->Get4Momentum() );
534 delete ddktv->operator[](i);
540 G4double eRecoil = sqrt( rM*rM + dP*dP );
578 if( products !=
nullptr )
580 for(
auto & prod : *products )
583 prod->GetTotalEnergy(),
584 prod->GetMomentum() ),
fSecID );
600 G4double rM(0.), mN(938.), det(0.), det2(0.);
651 G4double a = 4.*(sumE*sumE - pX*pX);
655 if(det2 > 0.) det = sqrt(det2);
663 if( pX < 0. ) pX = 0.;
665 eX = sqrt( dP*dP +
fMr*
fMr );
668 if(
A >= 1 ) lvN.
boost(bst);
678 G4double eRecoil = sqrt( rM*rM + pX*pX );
723 G4int pdgB(0), i(0), qM(0), qB(0);
724 G4double mM(0.), mB(0.), eM(0.), eB(0.), pM(0.), pB(0.);
725 G4double mm1(0.), mm22(0.), M1(0.), M2(0.), mX(0.);
733 G4double deltaMr[4] = { 0.*MeV, 0.*MeV, 100.*MeV, 0.*MeV};
751 if ( qX == 2 || qX == 0) { pdgB = 2212; qB = 1;}
753 else if( qX == 1 || qX == -1) { pdgB = 2112; qB = 0;}
757 else if( mX <
fBarMass[i] + deltaMr[i] || mX < mN + mPi )
761 if ( qX == 1 && pdgB != 2212) pdgB = pdgB - 10;
762 else if( qX == 0 && pdgB != 2212) pdgB = pdgB - 110;
763 else if( qX == 0 && pdgB == 2212) pdgB = pdgB - 100;
776 if( qX == 2 || qX == 0) { pdgB = 2212; qB = 1;}
777 else if( qX == 1 || qX == -1) { pdgB = 2112; qB = 0;}
789 mM = mm1 + rand*(mm22-mm1);
806 if ( qX == 2 || qX == 0) { pdgB = 2212; qB = 1;}
807 else if( qX == 1 || qX == -1) { pdgB = 2112; qB = 0;}
824 eM = 0.5*(mX*mX + mM*mM - mB*mB)/mX;
825 pM = sqrt(eM*eM - mM*mM);
829 eB = 0.5*(mX*mX + mB*mB - mM*mM)/mX;
830 pB = sqrt(eB*eB - mB*mB);
838 if ( qX == 2 ) { qM = 1; qB = 1;}
839 else if( qX == 1 ) { qM = 0; qB = 1;}
840 else if( qX == 0 ) { qM = 0; qB = 0;}
841 else if( qX == -1 ) { qM = -1; qB = 0;}
862 G4int pdgM(0), pdgB(0), i(0), qM(0), qB(0);
863 G4double mM(0.), mB(0.), eM(0.), eB(0.), pM(0.), pB(0.);
864 G4double mm1(0.), mm22(0.), M1(0.), M2(0.), mX(0.), Tkin(0.);
870 if( mX < 1080*MeV && mX > 990*MeV && Tkin < 600*MeV )
876 G4double deltaMr[4] = { 0.*MeV, 0.*MeV, 100.*MeV, 0.*MeV};
894 if ( qX == 1) { pdgB = 211; qB = 1;}
895 else if( qX == 0 ) { pdgB = 111; qB = 0;}
896 else if( qX == -1) { pdgB = -211; qB = -1;}
900 else if( mX <
fMesMass[i] + deltaMr[i] )
907 if( qX == 0 ) pdgB = pdgB - 100;
908 else if( qX == -1 ) pdgB = -pdgB;
910 if( finB )
return FinalMeson( lvX, qX, pdgB );
917 mm22 = mX - mPi - 1.*MeV;
921 if ( qX == 1) { pdgB = 211; qB = 1;}
922 else if( qX == 0 ) { pdgB = 111; qB = 0;}
923 else if( qX == -1) { pdgB = -211; qB = -1;}
932 if ( qX == 1 ) { qM = 1; qB = 0;}
933 else if( qX == 0 ) { qM = -1; qB = 1;}
934 else if( qX == -1 ) { qM = -1; qB = 0;}
943 mM = mm1 + rand*(mm22-mm1);
958 if ( qX == 1) { pdgB = 211; qB = 1;}
959 else if( qX == 0 ) { pdgB = 111; qB = 0;}
960 else if( qX == -1) { pdgB = -211; qB = -1;}
964 else if( mX <
fMesMass[i] + deltaMr[i] )
971 if( qX == 0 ) pdgB = pdgB - 100;
972 else if( qX == -1 ) pdgB = -pdgB;
974 if( finB )
return FinalMeson( lvX, qX, pdgB );
982 if ( qX == 1) { pdgB = 211; qB = 1;}
983 else if( qX == 0 ) { pdgB = 111; qB = 0;}
984 else if( qX == -1) { pdgB = -211; qB = -1;}
997 eM = 0.5*(mX*mX + mM*mM - mB*mB)/mX;
998 pM = sqrt(eM*eM - mM*mM);
1002 eB = 0.5*(mX*mX + mB*mB - mM*mM)/mX;
1003 pB = sqrt(eB*eB - mB*mB);
1013 if ( qM == 0 ) pdgM = pdgM - 100;
1014 else if( qM == -1 ) pdgM = -pdgM;
1045 if( delta2 >= 0. ) delta = sqrt(delta2);
1047 result = 0.5*(-b-delta)/a;
1070 if (
Z == 1 &&
A == 1 ) { kF = 0.; }
1071 else if (
Z == 1 &&
A == 2 ) { kF = 87.*MeV; }
1072 else if (
Z == 2 &&
A == 3 ) { kF = 134.*MeV; }
1073 else if (
Z == 6 &&
A == 12 ) { kF = 221.*MeV; }
1074 else if (
Z == 14 &&
A == 28 ) { kF = 239.*MeV; }
1075 else if (
Z == 26 &&
A == 56 ) { kF = 257.*MeV; }
1076 else if (
Z == 82 &&
A == 208 ) { kF = 265.*MeV; }
1079 kF = kp*ZpA*( 1 - pow(
G4double(
A), -t1 ) ) + kn*NpA*( 1 - pow(
G4double(
A), -t2 ) );
1119 const G4int maxBin = 12;
1121 G4double refA[maxBin] = { 2., 6., 12., 16., 27., 28., 40., 50., 56., 58., 197., 208. };
1123 G4double pEx[maxBin] = { 0., 12.2, 10.1, 10.9, 21.6, 12.4, 17.8, 17., 19., 16.8, 19.5, 14.7 };
1125 G4double nEx[maxBin] = { 0., 12.2, 10., 10.2, 21.6, 12.4, 21.8, 17., 19., 16.8, 19.5, 16.9 };
1129 if(fP)
for( i = 0; i < maxBin; ++i ) dE[i] = pEx[i];
1130 else dE[i] = nEx[i];
1132 for( i = 0; i < maxBin; ++i )
1134 if( aa <= refA[i] )
break;
1136 if( i >= maxBin ) eX = dE[maxBin-1];
1137 else if( i <= 0 ) eX = dE[0];
1144 if (a1 == a2 || e1 == e2 ) eX = dE[i];
1145 else eX = e1 + (e2-e1)*(aa-a1)/(a2-a1);
1159 G4double shift(1.), xx(1.), mom(0.), th(0.1);
1167 if(
A <= 12) th = 0.1;
1182 if(
A <= 12 ) ll = 6.0;
1201 if( mom > 2.*kF )
f2p2h =
true;
1215 G4int i, eIndex = 0;
1217 for( i = 0; i <
fIndex; i++)
1238 if( index <= 0 || energy <
fNuMuEnergy[0] ) ratio = 0.;
1251 ratio = y1 + (energy-x1)*angle;
1261 0.112103, 0.117359, 0.123119, 0.129443, 0.136404,
1262 0.144084, 0.152576, 0.161991, 0.172458, 0.184126,
1263 0.197171, 0.211801, 0.228261, 0.24684, 0.267887,
1264 0.291816, 0.319125, 0.350417, 0.386422, 0.428032,
1265 0.47634, 0.532692, 0.598756, 0.676612, 0.768868,
1266 0.878812, 1.01062, 1.16963, 1.36271, 1.59876,
1267 1.88943, 2.25002, 2.70086, 3.26916, 3.99166,
1268 4.91843, 6.11836, 7.6872, 9.75942, 12.5259,
1269 16.2605, 21.3615, 28.4141, 38.2903, 52.3062,
1270 72.4763, 101.93, 145.6, 211.39, 312.172
1280 0.98, 0.98, 0.98, 0.98, 0.98, 0.98, 0.98, 0.98, 0.98, 0.98,
1281 0.98, 0.98, 0.98, 0.98, 0.98, 0.98, 0.98, 0.98, 0.98, 0.98,
1282 0.97, 0.96, 0.95, 0.93,
1283 0.917794, 0.850239, 0.780412, 0.709339, 0.638134, 0.568165,
1284 0.500236, 0.435528, 0.375015, 0.319157, 0.268463, 0.2232, 0.183284,
1285 0.148627, 0.119008, 0.0940699, 0.0733255, 0.0563819, 0.0427312, 0.0319274,
1286 0.0235026, 0.0170486, 0.0122149, 0.00857825, 0.00594018, 0.00405037
1295 G4int i, eIndex = 0;
1331 ratio = y1 + (energy-x1)*angle;
1342 0.275314, 0.293652, 0.31729, 0.33409, 0.351746, 0.365629, 0.380041, 0.400165, 0.437941, 0.479237,
1343 0.504391, 0.537803, 0.588487, 0.627532, 0.686839, 0.791905, 0.878332, 0.987405, 1.08162, 1.16971,
1344 1.2982, 1.40393, 1.49854, 1.64168, 1.7524, 1.87058, 2.02273, 2.15894, 2.3654, 2.55792, 2.73017,
1345 3.03005, 3.40733, 3.88128, 4.53725, 5.16786, 5.73439, 6.53106, 7.43879, 8.36214, 9.39965, 10.296,
1346 11.5735, 13.1801, 15.2052, 17.5414, 19.7178, 22.7462, 25.9026, 29.4955, 33.5867, 39.2516, 46.4716,
1347 53.6065, 63.4668, 73.2147, 85.5593, 99.9854
1355 0.0019357, 0.0189361, 0.0378722, 0.0502758, 0.0662559, 0.0754581, 0.0865008, 0.0987275, 0.124112,
1356 0.153787, 0.18308, 0.213996, 0.245358, 0.274425, 0.301536, 0.326612, 0.338208, 0.337806, 0.335948,
1357 0.328092, 0.313557, 0.304965, 0.292169, 0.28481, 0.269474, 0.254138, 0.247499, 0.236249, 0.221654,
1358 0.205492, 0.198781, 0.182216, 0.162251, 0.142878, 0.128631, 0.116001, 0.108435, 0.0974843, 0.082092,
1359 0.0755204, 0.0703121, 0.0607066, 0.0554278, 0.0480401, 0.0427023, 0.0377123, 0.0323248, 0.0298584,
1360 0.0244296, 0.0218526, 0.019121, 0.016477, 0.0137309, 0.0137963, 0.0110371, 0.00834028, 0.00686127, 0.00538226
1368 G4double qerata(0.5), rr(0.), x1(0.), x2(0.), y1(0.), y2(0.), aa(0.);
1371 if(
A >
Z )
N =
A-
Z;
1373 for( i = 0; i < 50; i++)
1377 if(i <= 0)
return 1.;
1378 else if (i >= 49)
return 0.;
1384 if( nepdg == 12 || nepdg == 14 )
1398 aa = (y2-y1)/(x2-x1);
1399 rr = y1 + (energy-x1)*aa;
1401 if( nepdg == 12 || nepdg == 14 ) qerata =
N*rr/(
N*rr +
A*( 1 - rr ) );
1402 else qerata =
Z*rr/(
Z*rr +
A*( 1 - rr ) );
1411 0.12, 0.1416, 0.167088, 0.197164, 0.232653, 0.274531, 0.323946, 0.382257, 0.451063, 0.532254,
1412 0.62806, 0.741111, 0.874511, 1.03192, 1.21767, 1.43685, 1.69548, 2.00067, 2.36079, 2.78573,
1413 3.28716, 3.87885, 4.57705, 5.40092, 6.37308, 7.52024, 8.87388, 10.4712, 12.356, 14.5801,
1414 17.2045, 20.3013, 23.9555, 28.2675, 33.3557, 39.3597, 46.4444, 54.8044, 64.6692, 76.3097,
1415 90.0454, 106.254, 125.379, 147.947, 174.578, 206.002, 243.082, 286.837, 338.468, 399.392
1420 1, 1, 1, 1, 1, 1, 1, 0.97506, 0.920938, 0.847671, 0.762973, 0.677684, 0.597685,
1421 0.52538, 0.461466, 0.405329, 0.356154, 0.312944, 0.274984, 0.241341, 0.211654, 0.185322,
1422 0.161991, 0.141339, 0.123078, 0.106952, 0.0927909, 0.0803262, 0.0693698, 0.0598207, 0.0514545,
1423 0.044193, 0.0378696, 0.0324138, 0.0276955, 0.0236343, 0.0201497, 0.0171592, 0.014602, 0.0124182,
1424 0.0105536, 0.00896322, 0.00761004, 0.00645821, 0.00547859, 0.00464595, 0.00393928,
1425 0.00333961, 0.00283086, 0.00239927
1430 1, 1, 1, 1, 1, 1, 1, 0.977592, 0.926073, 0.858783, 0.783874, 0.706868, 0.63113, 0.558681,
1431 0.490818, 0.428384, 0.371865, 0.321413, 0.276892, 0.237959, 0.204139, 0.1749, 0.149706, 0.128047,
1432 0.109456, 0.093514, 0.0798548, 0.0681575, 0.0581455, 0.0495804, 0.0422578, 0.036002, 0.0306614,
1433 0.0261061, 0.0222231, 0.0189152, 0.0160987, 0.0137011, 0.0116604, 0.00992366, 0.00844558, 0.00718766,
1434 0.00611714, 0.00520618, 0.00443105, 0.00377158, 0.00321062, 0.0027335, 0.00232774, 0.00198258
G4double B(G4double temperature)
G4double G4Log(G4double x)
CLHEP::HepLorentzVector G4LorentzVector
std::vector< G4ReactionProduct * > G4ReactionProductVector
CLHEP::Hep3Vector G4ThreeVector
Hep3Vector orthogonal() const
Hep3Vector boostVector() const
HepLorentzVector & boost(double, double, double)
void SetNumberOfHoles(G4int valueTot, G4int valueP=0)
void SetExcEnergyAndMomentum(G4double eexc, const G4LorentzVector &)
void AddSecondary(G4DynamicParticle *aP, G4int mod=-1)
void SetLocalEnergyDeposit(G4double aE)
const G4ParticleDefinition * GetDefinition() const
G4double GetTotalEnergy() const
G4HadFinalState theParticleChange
void SetMinEnergy(G4double anEnergy)
const G4String & GetModelName() const
void SetMaxEnergy(const G4double anEnergy)
G4KineticTrackVector * Decay()
static G4MuonMinus * MuonMinus()
static G4MuonPlus * MuonPlus()
static const G4double fBarMass[4]
void CoherentPion(G4LorentzVector &lvP, G4int pdgP, G4Nucleus &targetNucleus)
static G4double fNuMuQarrayKR[50][51][51]
virtual void ModelDescription(std::ostream &) const
G4GeneratorPrecompoundInterface * fPrecoInterface
G4double FermiMomentum(G4Nucleus &targetNucleus)
static G4double fNuMuXarrayKR[50][51]
G4double NucleonMomentum(G4Nucleus &targetNucleus)
G4int GetOnePionIndex(G4double energy)
static const G4double fResMass[6]
static const G4int fMesPDG[4]
static const G4double fNuMuEnergy[50]
G4ExcitationHandler * fDeExcitation
G4double SampleXkr(G4double energy)
G4ParticleDefinition * theMuonMinus
static const G4double fQEnergy[50]
static const G4int fResNumber
G4double SampleQkr(G4double energy, G4double xx)
G4double GetNuMuQeTotRat(G4int index, G4double energy)
G4int GetEnergyIndex(G4double energy)
G4double GetXkr(G4int iEnergy, G4double prob)
void MesonDecay(G4LorentzVector &lvX, G4int qX)
G4double GgSampleNM(G4Nucleus &nucl)
G4double GetNuMuOnePionProb(G4int index, G4double energy)
virtual ~G4NeutrinoNucleusModel()
static const G4double fNeMuQEratio[50]
void FinalMeson(G4LorentzVector &lvM, G4int qM, G4int pdgM)
G4double GetQkr(G4int iE, G4int jX, G4double prob)
static const G4double fMesMass[4]
G4PreCompoundModel * fPreCompound
static const G4int fBarPDG[4]
static G4double fNuMuXdistrKR[50][50]
static G4double fNuMuQdistrKR[50][51][50]
G4double CalculateQEratioA(G4int Z, G4int A, G4double energy, G4int nepdg)
virtual G4bool IsApplicable(const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
static const G4double fNuMuQeTotRat[50]
G4NeutrinoNucleusModel(const G4String &name="neutrino-nucleus")
G4double FinalMomentum(G4double mI, G4double mF, G4double mP, G4LorentzVector lvX)
G4ParticleDefinition * theMuonPlus
G4double GetMinNuMuEnergy()
void ClusterDecay(G4LorentzVector &lvX, G4int qX)
void FinalBarion(G4LorentzVector &lvB, G4int qB, G4int pdgB)
static const G4double fOnePionProb[58]
static const G4double fANeMuQEratio[50]
static const G4int fClustNumber
void RecoilDeexcitation(G4Fragment &fragment)
static const G4double fNuMuEnergyLogVector[50]
static const G4double fOnePionEnergy[58]
static G4double GetNuclearMass(const G4double A, const G4double Z)
G4double AtomicMass(const G4double A, const G4double Z, const G4int numberOfLambdas=0) const
G4double GetPDGMass() const
const G4String & GetParticleName() const
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
static G4ParticleTable * GetParticleTable()
static G4int GetModelID(const G4int modelIndex)
virtual G4ReactionProductVector * DeExcite(G4Fragment &aFragment) final
static G4Proton * Proton()
void SetDeExcitation(G4VPreCompoundModel *ptr)