diff --git a/config/HAIntranuke2018.xml b/config/HAIntranuke2018.xml
index 9db9957e6..7e47d22d1 100644
--- a/config/HAIntranuke2018.xml
+++ b/config/HAIntranuke2018.xml
@@ -52,7 +52,7 @@ DelRNucleon double Yes mult. factor for nucleon de-Broglie wavelength
true
true
-
+
0.0
0.0
@@ -74,7 +74,7 @@ DelRNucleon double Yes mult. factor for nucleon de-Broglie wavelength
1.000
1.000
-->
-
+
diff --git a/src/Physics/HadronTransport/HAIntranuke2018.cxx b/src/Physics/HadronTransport/HAIntranuke2018.cxx
index 682470dde..108bf7b81 100644
--- a/src/Physics/HadronTransport/HAIntranuke2018.cxx
+++ b/src/Physics/HadronTransport/HAIntranuke2018.cxx
@@ -2,7 +2,7 @@
/*
Copyright (c) 2003-2025, The GENIE Collaboration
For the full text of the license visit http://copyright.genie-mc.org
-
+
Author: Steve Dytman , Pittsburgh Univ.
Aaron Meyer , Pittsburgh Univ.
@@ -781,7 +781,7 @@ void HAIntranuke2018::Inelastic(
GHepParticle s3(*p);
bool success = utils::intranuke2018::PionProduction(
- ev,p,&s1,&s2,&s3,fRemnA,fRemnZ,fRemnP4, fDoFermi,fFermiFac,fFermiMomentum,fNuclmodel);
+ ev,p,&s1,&s2,&s3,fRemnA,fRemnZ,fRemnP4, fDoFermi,fFermiFac,fFermiMomentum,fNuclmodel,fPiProdThreeBodyBias);
if (success){
LOG ("HAIntranuke2018",pINFO) << " successful pion production fate";
@@ -1562,6 +1562,8 @@ void HAIntranuke2018::LoadConfig(void)
GetParamDef( "FSI-Nucleon-FracAbsScale", fNucleonFracAbsScale, 1.0 ) ;
GetParamDef( "FSI-Nucleon-FracPiProdScale", fNucleonFracPiProdScale, 1.0 ) ;
+ GetParamDef( "FSI-PiProd-ThreeBodyBias", fPiProdThreeBodyBias, 0.0 ) ;
+
// report
LOG("HAIntranuke2018", pINFO) << "Settings for INTRANUKE mode: " << INukeMode::AsString(kIMdHA);
LOG("HAIntranuke2018", pINFO) << "R0 = " << fR0 << " fermi";
@@ -1577,6 +1579,7 @@ void HAIntranuke2018::LoadConfig(void)
LOG("HAIntranuke2018", pINFO) << "DoFermi? = " << ((fDoFermi)?(true):(false));
LOG("HAIntranuke2018", pINFO) << "DoCmpndNuc? = " << ((fDoCompoundNucleus)?(true):(false));
LOG("HAIntranuke2018", pINFO) << "XsecNNCorr? = " << ((fXsecNNCorr)?(true):(false));
+ LOG("HAIntranuke2018", pINFO) << "PiProdBias = " << fPiProdThreeBodyBias;
}
//___________________________________________________________________________
/*
diff --git a/src/Physics/HadronTransport/HAIntranuke2025.cxx b/src/Physics/HadronTransport/HAIntranuke2025.cxx
index cb5027f69..746211645 100644
--- a/src/Physics/HadronTransport/HAIntranuke2025.cxx
+++ b/src/Physics/HadronTransport/HAIntranuke2025.cxx
@@ -2,7 +2,7 @@
/*
Copyright (c) 2003-2025, The GENIE Collaboration
For the full text of the license visit http://copyright.genie-mc.org
-
+
Author: Steve Dytman , Pittsburgh Univ.
Aaron Meyer , Pittsburgh Univ.
@@ -15,11 +15,11 @@
Important revisions after version 2.0.0 :
@ sep 2025 - Mohamed Ismail, SD
- new class for hA2025. no changes from the 2018 class version, major changes made in INukeHadroData2025
-for new hA pion splines Add data, use hN for high pion KE, use INCL for low energy. Use splines for channel
+ new class for hA2025. no changes from the 2018 class version, major changes made in INukeHadroData2025
+for new hA pion splines Add data, use hN for high pion KE, use INCL for low energy. Use splines for channel
and total reac xs to improve accuracy. Also, smooth results to avoid discontinuities.
-
+
*/
//____________________________________________________________________________
@@ -242,8 +242,8 @@ INukeFateHA_t HAIntranuke2025::HadronFateHA(const GHepParticle * p) const
if (pdgc==kPdgPiP || pdgc==kPdgPiM) frac_abs *= fChPionFracAbsScale;
if (pdgc==kPdgPi0) frac_abs *= fNeutralPionFracAbsScale;
frac_piprod *= fPionFracPiProdScale;
-
-
+
+
// Flag to enable or disable π0/π+ ratio corrections
bool apply_pi0_ratio_correction = true;
@@ -261,7 +261,7 @@ INukeFateHA_t HAIntranuke2025::HadronFateHA(const GHepParticle * p) const
double ratio_piprod = 0.0004402 * ke_ratio + 0.47418;
// Apply the corrections only if kinetic energy is below 1000 MeV
-
+
frac_cex *= ratio_cex;
frac_abs *= ratio_abs;
frac_inel *= ratio_inel;
@@ -270,11 +270,11 @@ INukeFateHA_t HAIntranuke2025::HadronFateHA(const GHepParticle * p) const
if (ke > 400.0) {
frac_piprod *= ratio_piprod;
}
-
+
}
-
-
+
+
double frac_rescale = 1./(frac_cex + frac_inel + frac_abs + frac_piprod);
@@ -793,7 +793,7 @@ void HAIntranuke2025::Inelastic(
GHepParticle s3(*p);
bool success = utils::intranuke2025::PionProduction(
- ev,p,&s1,&s2,&s3,fRemnA,fRemnZ,fRemnP4, fDoFermi,fFermiFac,fFermiMomentum,fNuclmodel);
+ ev,p,&s1,&s2,&s3,fRemnA,fRemnZ,fRemnP4, fDoFermi,fFermiFac,fFermiMomentum,fNuclmodel,fPiProdThreeBodyBias);
if (success){
LOG ("HAIntranuke2025",pINFO) << " successful pion production fate";
@@ -1131,10 +1131,10 @@ void HAIntranuke2025::Inelastic(
if ( pdg::IsNeutronOrProton (pdgc) ) //nucleon probe
{
- // work done in 2011 by Meyer and Dytman. Make fits to hN simulation output for p and n.
+ // work done in 2011 by Meyer and Dytman. Make fits to hN simulation output for p and n.
ns = -TMath::Log(rnd->RndFsi().Rndm())/gam_ns; // exponential random variable
}
- if ( pdg::IsKaon (pdgc) ) //charged kaon probe - either 2 or 3 nucleons to stay simple.
+ if ( pdg::IsKaon (pdgc) ) //charged kaon probe - either 2 or 3 nucleons to stay simple.
// No data available, use hN simulation (sd)
{
ns = (rnd->RndFsi().Rndm()<0.5?2:3);
@@ -1579,6 +1579,8 @@ void HAIntranuke2025::LoadConfig(void)
GetParamDef( "FSI-Nucleon-FracAbsScale", fNucleonFracAbsScale, 1.0 ) ;
GetParamDef( "FSI-Nucleon-FracPiProdScale", fNucleonFracPiProdScale, 1.0 ) ;
+ GetParamDef( "FSI-PiProd-ThreeBodyBias", fPiProdThreeBodyBias, 0.0 );
+
// report
LOG("HAIntranuke2025", pINFO) << "Settings for INTRANUKE mode: " << INukeMode::AsString(kIMdHA);
LOG("HAIntranuke2025", pINFO) << "R0 = " << fR0 << " fermi";
@@ -1594,6 +1596,7 @@ void HAIntranuke2025::LoadConfig(void)
LOG("HAIntranuke2025", pINFO) << "DoFermi? = " << ((fDoFermi)?(true):(false));
LOG("HAIntranuke2025", pINFO) << "DoCmpndNuc? = " << ((fDoCompoundNucleus)?(true):(false));
LOG("HAIntranuke2025", pINFO) << "XsecNNCorr? = " << ((fXsecNNCorr)?(true):(false));
+ LOG("HAIntranuke2025", pINFO) << "PiProdBias = " << fPiProdThreeBodyBias;
}
//___________________________________________________________________________
/*
diff --git a/src/Physics/HadronTransport/HNIntranuke2018.cxx b/src/Physics/HadronTransport/HNIntranuke2018.cxx
index 3a1b6d796..2a6657f19 100644
--- a/src/Physics/HadronTransport/HNIntranuke2018.cxx
+++ b/src/Physics/HadronTransport/HNIntranuke2018.cxx
@@ -3,7 +3,7 @@
/*
Copyright (c) 2003-2025, The GENIE Collaboration
For the full text of the license visit http://copyright.genie-mc.org
-
+
Author: Steve Dytman , Pittsburgh Univ.
Aaron Meyer , Pittsburgh Univ.
@@ -18,22 +18,22 @@
@ Nov 30, 2007 - SD
Changed the hadron tracking algorithm to take into account the radial
nuclear density dependence. Using the somewhat empirical approach of
- increasing the nuclear radius by a const (tunable) number times the tracked
- particle's de Broglie wavelength as this helps getting the hadron+nucleus
+ increasing the nuclear radius by a const (tunable) number times the tracked
+ particle's de Broglie wavelength as this helps getting the hadron+nucleus
cross sections right.
@ Mar 08, 2008 - CA
Fixed code retrieving the remnant nucleus which stopped working as soon as
simulation of nuclear de-excitation started pushing photons in the target
nucleus daughter list.
@ Jun 20, 2008 - CA
- Fix a mem leak: The (clone of the) GHepParticle being re-scattered was not
+ Fix a mem leak: The (clone of the) GHepParticle being re-scattered was not
deleted after it was added at the GHEP event record.
@ Jul 15, 2010 - AM
The hN mode is now implemented in Intranuke. Similar to hA mode, but particles
produced by reactions are stepped through the nucleus like probe particles.
Particles react with nucleons instead of the entire nucleus, and final states
are determined after reactions are finished, not before.
- @ Dec 15, 2014 - SD, Nick Geary
+ @ Dec 15, 2014 - SD, Nick Geary
Update fates to include Compound Nucleus final state correctly.
@ Jan 9, 2015 - SD, NG, Tomek Golan
Added 2014 version of INTRANUKE codes (new class) for independent development.
@@ -112,10 +112,10 @@ HNIntranuke2018::~HNIntranuke2018()
//___________________________________________________________________________
void HNIntranuke2018::ProcessEventRecord(GHepRecord * evrec) const
{
- LOG("HNIntranuke2018", pNOTICE)
+ LOG("HNIntranuke2018", pNOTICE)
<< "************ Running hN2018 MODE INTRANUKE ************";
-
- /* LOG("HNIntranuke2018", pWARN)
+
+ /* LOG("HNIntranuke2018", pWARN)
<< print::PrintFramedMesg(
"Experimental code (INTRANUKE/hN model) - Run at your own risk");
*/
@@ -140,7 +140,7 @@ void HNIntranuke2018::SimulateHadronicFinalState(GHepRecord* ev, GHepParticle* p
bool is_pion = (pdgc==kPdgPiP || pdgc==kPdgPiM || pdgc==kPdgPi0);
bool is_kaon = (pdgc==kPdgKP);
bool is_baryon = (pdgc==kPdgProton || pdgc==kPdgNeutron);
- bool is_gamma = (pdgc==kPdgGamma);
+ bool is_gamma = (pdgc==kPdgGamma);
if(!(is_pion || is_baryon || is_gamma || is_kaon))
{
LOG("HNIntranuke2018", pERROR) << "** Cannot handle particle: " << p->Name();
@@ -157,7 +157,7 @@ void HNIntranuke2018::SimulateHadronicFinalState(GHepRecord* ev, GHepParticle* p
if(fate == kIHNFtUndefined)
{
LOG("HNIntranuke2018", pERROR) << "** Couldn't select a fate";
- LOG("HNIntranuke2018", pERROR) << "** Num Protons: " << fRemnZ
+ LOG("HNIntranuke2018", pERROR) << "** Num Protons: " << fRemnZ
<< ", Num Neutrons: "<<(fRemnA-fRemnZ);
LOG("HNIntranuke2018", pERROR) << "** Particle: " << "\n" << (*p);
//LOG("HNIntranuke2018", pERROR) << "** Event Record: " << "\n" << (*ev);
@@ -176,7 +176,7 @@ void HNIntranuke2018::SimulateHadronicFinalState(GHepRecord* ev, GHepParticle* p
this->ElasHN(ev,p,fate);
}
else if(fate == kIHNFtAbs) {this-> AbsorbHN(ev,p,fate);}
- else if(fate == kIHNFtInelas && pdgc != kPdgGamma)
+ else if(fate == kIHNFtInelas && pdgc != kPdgGamma)
{
#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
LOG("HNIntranuke2018", pDEBUG)
@@ -198,7 +198,7 @@ void HNIntranuke2018::SimulateHadronicFinalState(GHepRecord* ev, GHepParticle* p
catch(exceptions::INukeException exception)
{
this->SimulateHadronicFinalState(ev,p);
- LOG("HNIntranuke2018", pNOTICE)
+ LOG("HNIntranuke2018", pNOTICE)
<< "retry call to SimulateHadronicFinalState ";
LOG("HNIntranuke2018", pNOTICE) << exception;
@@ -218,8 +218,8 @@ INukeFateHN_t HNIntranuke2018::HadronFateHN(const GHepParticle * p) const
bool isPion = (pdgc == kPdgPiP or pdgc == kPdgPi0 or pdgc == kPdgPiM);
if (isPion and fUseOset and ke < 350.0) return HadronFateOset ();
-
- LOG("HNIntranuke2018", pNOTICE)
+
+ LOG("HNIntranuke2018", pNOTICE)
<< "Selecting hN fate for " << p->Name() << " with KE = " << ke << " MeV";
// try to generate a hadron fate
@@ -244,7 +244,7 @@ INukeFateHN_t HNIntranuke2018::HadronFateHN(const GHepParticle * p) const
frac_elas *= fNucQEFac;
if(pdgc==kPdgPi0) frac_abs*= 0.665; //isospin factor
- LOG("HNIntranuke2018", pNOTICE)
+ LOG("HNIntranuke2018", pNOTICE)
<< "\n frac{" << INukeHadroFates::AsString(kIHNFtCEx) << "} = " << frac_cex
<< "\n frac{" << INukeHadroFates::AsString(kIHNFtElas) << "} = " << frac_elas
<< "\n frac{" << INukeHadroFates::AsString(kIHNFtInelas) << "} = " << frac_inel
@@ -253,7 +253,7 @@ INukeFateHN_t HNIntranuke2018::HadronFateHN(const GHepParticle * p) const
// compute total fraction (can be <1 if fates have been switched off)
double tf = frac_cex +
frac_elas +
- frac_inel +
+ frac_inel +
frac_abs;
double r = tf * rnd->RndFsi().Rndm();
@@ -266,10 +266,10 @@ INukeFateHN_t HNIntranuke2018::HadronFateHN(const GHepParticle * p) const
if(r < (cf += frac_cex )) return kIHNFtCEx; //cex
if(r < (cf += frac_elas )) return kIHNFtElas; //elas
if(r < (cf += frac_inel )) return kIHNFtInelas; //inelas
- if(r < (cf += frac_abs )) return kIHNFtAbs; //abs
+ if(r < (cf += frac_abs )) return kIHNFtAbs; //abs
- LOG("HNIntranuke2018", pWARN)
- << "No selection after going through all fates! "
+ LOG("HNIntranuke2018", pWARN)
+ << "No selection after going through all fates! "
<< "Total fraction = " << tf << " (r = " << r << ")";
////////////////////////////
return kIHNFtUndefined;
@@ -285,7 +285,7 @@ INukeFateHN_t HNIntranuke2018::HadronFateHN(const GHepParticle * p) const
double frac_cmp = this->FateWeight(pdgc, kIHNFtCmp)
* fHadroData2018->Frac(pdgc, kIHNFtCmp, ke, fRemnA , fRemnZ);
- LOG("HNIntranuke2018", pINFO)
+ LOG("HNIntranuke2018", pINFO)
<< "\n frac{" << INukeHadroFates::AsString(kIHNFtElas) << "} = " << frac_elas
<< "\n frac{" << INukeHadroFates::AsString(kIHNFtInelas) << "} = " << frac_inel;
@@ -305,14 +305,14 @@ INukeFateHN_t HNIntranuke2018::HadronFateHN(const GHepParticle * p) const
if(r < (cf += frac_inel )) return kIHNFtInelas; // inelas
if(r < (cf += frac_cmp )) return kIHNFtCmp; // cmp
- LOG("HNIntranuke2018", pWARN)
+ LOG("HNIntranuke2018", pWARN)
<< "No selection after going through all fates! "
<< "Total fraction = " << tf << " (r = " << r << ")";
//////////////////////////
return kIHNFtUndefined;
}
- // handle gamma -- does not currently consider the elastic case
+ // handle gamma -- does not currently consider the elastic case
else if (pdgc==kPdgGamma) return kIHNFtInelas;
// Handle kaon -- elastic + charge exchange
else if (pdgc==kPdgKP){
@@ -324,7 +324,7 @@ INukeFateHN_t HNIntranuke2018::HadronFateHN(const GHepParticle * p) const
// frac_cex *= fNucCEXFac; // scaling factors
// frac_elas *= fNucQEFac; // Flor - Correct scaling factors?
- LOG("HNIntranuke", pINFO)
+ LOG("HNIntranuke", pINFO)
<< "\n frac{" << INukeHadroFates::AsString(kIHNFtCEx) << "} = " << frac_cex
<< "\n frac{" << INukeHadroFates::AsString(kIHNFtElas) << "} = " << frac_elas;
@@ -340,10 +340,10 @@ INukeFateHN_t HNIntranuke2018::HadronFateHN(const GHepParticle * p) const
double cf=0; // current fraction
if(r < (cf += frac_cex )) return kIHNFtCEx; //cex
- if(r < (cf += frac_elas )) return kIHNFtElas; //elas
+ if(r < (cf += frac_elas )) return kIHNFtElas; //elas
- LOG("HNIntranuke", pWARN)
- << "No selection after going through all fates! "
+ LOG("HNIntranuke", pWARN)
+ << "No selection after going through all fates! "
<< "Total fraction = " << tf << " (r = " << r << ")";
////////////////////////////
return kIHNFtUndefined;
@@ -367,7 +367,7 @@ double HNIntranuke2018::FateWeight(int pdgc, INukeFateHN_t fate) const
int np = fRemnZ;
int nn = fRemnA - fRemnZ;
-
+
if (np < 1 && nn < 1)
{
LOG("HNIntranuke2018", pERROR) << "** Nothing left in nucleus!!! **";
@@ -389,7 +389,7 @@ void HNIntranuke2018::AbsorbHN(
GHepRecord * ev, GHepParticle * p, INukeFateHN_t fate) const
{
// handles pi+d->2p, pi-d->nn, pi0 d->pn absorbtion, all using pi+d values
-
+
int pdgc = p->Pdg();
#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
@@ -437,7 +437,7 @@ void HNIntranuke2018::AbsorbHN(
// Library instance for reference
PDGLibrary * pLib = PDGLibrary::Instance();
-
+
// Handle fermi target
Target target(ev->TargetNucleus()->Pdg());
@@ -473,7 +473,7 @@ void HNIntranuke2018::AbsorbHN(
<< "AbsorbHN() cannot handle probe: " << pdgc;
return;
}
-
+
// assign proper masses
M1 = pLib->Find(pcode) ->Mass();
M2_1 = pLib->Find(t1code)->Mass();
@@ -481,14 +481,14 @@ void HNIntranuke2018::AbsorbHN(
M3 = pLib->Find(scode) ->Mass();
M4 = pLib->Find(s2code)->Mass();
- // handle fermi momentum
+ // handle fermi momentum
if(fDoFermi)
{
target.SetHitNucPdg(t1code);
fNuclmodel->GenerateNucleon(target);
tP2_1L=fFermiFac * fNuclmodel->Momentum3();
E2_1L = TMath::Sqrt(tP2_1L.Mag2() + M2_1*M2_1);
-
+
target.SetHitNucPdg(t2code);
fNuclmodel->GenerateNucleon(target);
tP2_2L=fFermiFac * fNuclmodel->Momentum3();
@@ -508,7 +508,7 @@ void HNIntranuke2018::AbsorbHN(
// adjust p to reflect scattering
// get random scattering angle
C3CM = fHadroData2018->IntBounce(p,t1code,scode,fate);
- if (C3CM<-1.)
+ if (C3CM<-1.)
{
p->SetStatus(kIStStableFinalState);
ev->AddParticle(*p);
@@ -598,7 +598,7 @@ void HNIntranuke2018::AbsorbHN(
p->SetStatus(kIStHadronInTheNucleus);
//disable until needed
// utils::intranuke2018::StepParticle(p,fFreeStep,fTrackingRadius);
- ev->AddParticle(*p);
+ ev->AddParticle(*p);
return;
*/
// new attempt at error handling:
@@ -616,12 +616,12 @@ void HNIntranuke2018::AbsorbHN(
// get random phi angle, distributed uniformally in 360 deg
PHI3 = 2 * kPi * rnd->RndFsi().Rndm();
-
+
tP3L = P3zL*bDir + P3tL*tTrans;
tP4L = P4zL*bDir + P4tL*tTrans;
tP3L.Rotate(PHI3,bDir); // randomize transverse components
- tP4L.Rotate(PHI3,bDir);
+ tP4L.Rotate(PHI3,bDir);
E3L = TMath::Sqrt(P3L*P3L + M3*M3);
E4L = TMath::Sqrt(P4L*P4L + M4*M4);
@@ -708,7 +708,7 @@ void HNIntranuke2018::ElasHN(
// get random scattering angle
double C3CM = fHadroData2018->IntBounce(p,tcode,scode,fate);
- if (C3CM<-1.)
+ if (C3CM<-1.)
{
p->SetStatus(kIStStableFinalState);
ev->AddParticle(*p);
@@ -721,7 +721,7 @@ void HNIntranuke2018::ElasHN(
double Mt = t->Mass();
//t->SetMomentum(TLorentzVector(0,0,0,Mt));
t->SetRemovalEnergy(0);
- // handle fermi momentum
+ // handle fermi momentum
if(fDoFermi)
{
// Handle fermi target
@@ -768,24 +768,25 @@ void HNIntranuke2018::ElasHN(
void HNIntranuke2018::InelasticHN(GHepRecord* ev, GHepParticle* p) const
{
// Aaron Meyer (Jan 2010)
- // Updated version of InelasticHN
+ // Updated version of InelasticHN
- GHepParticle s1(*p);
+ GHepParticle s1(*p);
GHepParticle s2(*p);
GHepParticle s3(*p);
s2.SetRemovalEnergy(0);
s3.SetRemovalEnergy(0);
-
-
-
- if (utils::intranuke2018::PionProduction(ev,p,&s1,&s2,&s3,fRemnA,fRemnZ,fRemnP4,fDoFermi,fFermiFac,fFermiMomentum,fNuclmodel))
+
+
+
+ if (utils::intranuke2018::PionProduction(ev,p,&s1,&s2,&s3,fRemnA,fRemnZ,fRemnP4,fDoFermi,fFermiFac,fFermiMomentum,
+ fNuclmodel,fPiProdThreeBodyBias))
{
// set status of particles and return
-
+
s1.SetStatus(kIStHadronInTheNucleus);
s2.SetStatus(kIStHadronInTheNucleus);
s3.SetStatus(kIStHadronInTheNucleus);
-
+
ev->AddParticle(s1);
ev->AddParticle(s2);
ev->AddParticle(s3);
@@ -801,7 +802,7 @@ void HNIntranuke2018::InelasticHN(GHepRecord* ev, GHepParticle* p) const
}
//___________________________________________________________________________
-void HNIntranuke2018::GammaInelasticHN(GHepRecord* ev, GHepParticle* p, INukeFateHN_t fate) const
+void HNIntranuke2018::GammaInelasticHN(GHepRecord* ev, GHepParticle* p, INukeFateHN_t fate) const
{
// This function handles pion photoproduction reactions
@@ -850,7 +851,7 @@ void HNIntranuke2018::GammaInelasticHN(GHepRecord* ev, GHepParticle* p, INukeFat
<< "Error: could not determine particle final states";
ev->AddParticle(*p);
return;
- }
+ }
LOG("HNIntranuke2018", pNOTICE)
<< "GammaInelastic fate: " << INukeHadroFates::AsString(fate);
@@ -863,7 +864,7 @@ void HNIntranuke2018::GammaInelasticHN(GHepRecord* ev, GHepParticle* p, INukeFat
t->SetPdgCode(tcode);
double Mt = t->Mass();
- // handle fermi momentum
+ // handle fermi momentum
if(fDoFermi)
{
// Handle fermi target
@@ -910,7 +911,7 @@ int HNIntranuke2018::HandleCompoundNucleus(GHepRecord* ev, GHepParticle* p, int
// handle compound nucleus option
// -- Call the PreEquilibrium function
- if( fDoCompoundNucleus && IsInNucleus(p) && pdg::IsNeutronOrProton(p->Pdg()))
+ if( fDoCompoundNucleus && IsInNucleus(p) && pdg::IsNeutronOrProton(p->Pdg()))
{ // random number generator
//unused var - quiet compiler warning//RandomGen * rnd = RandomGen::Instance();
@@ -986,6 +987,8 @@ void HNIntranuke2018::LoadConfig(void)
GetParamDef( "FSI-NeutralPion-MFPScale", fNeutralPionMFPScale, 1.0 ) ;
GetParamDef( "FSI-Nucleon-MFPScale", fNucleonMFPScale, 1.0 ) ;
+ GetParamDef( "FSI-PiProd-ThreeBodyBias", fPiProdThreeBodyBias, 0.0 ) ;
+
// report
LOG("HNIntranuke2018", pINFO) << "Settings for Intranuke2018 mode: " << INukeMode::AsString(kIMdHN);
LOG("HNIntranuke2018", pWARN) << "R0 = " << fR0 << " fermi";
@@ -1006,6 +1009,7 @@ void HNIntranuke2018::LoadConfig(void)
LOG("HNIntranuke2018", pWARN) << "XsecNNCorr? = " << ((fXsecNNCorr)?(true):(false));
LOG("HNIntranuke2018", pWARN) << "FSI-ChargedPion-MFPScale = " << fChPionMFPScale;
LOG("HNIntranuke2018", pWARN) << "FSI-NeutralPion-MFPScale = " << fNeutralPionMFPScale;
+ LOG("HAIntranuke2018", pWARN) << "PiProdBias = " << fPiProdThreeBodyBias;
}
//___________________________________________________________________________
diff --git a/src/Physics/HadronTransport/HNIntranuke2025.cxx b/src/Physics/HadronTransport/HNIntranuke2025.cxx
index b00ffbcc6..81640fe16 100644
--- a/src/Physics/HadronTransport/HNIntranuke2025.cxx
+++ b/src/Physics/HadronTransport/HNIntranuke2025.cxx
@@ -3,7 +3,7 @@
/*
Copyright (c) 2003-2025, The GENIE Collaboration
For the full text of the license visit http://copyright.genie-mc.org
-
+
Author: Steve Dytman , Pittsburgh Univ.
Aaron Meyer , Pittsburgh Univ.
@@ -18,22 +18,22 @@
@ Nov 30, 2007 - SD
Changed the hadron tracking algorithm to take into account the radial
nuclear density dependence. Using the somewhat empirical approach of
- increasing the nuclear radius by a const (tunable) number times the tracked
- particle's de Broglie wavelength as this helps getting the hadron+nucleus
+ increasing the nuclear radius by a const (tunable) number times the tracked
+ particle's de Broglie wavelength as this helps getting the hadron+nucleus
cross sections right.
@ Mar 08, 2008 - CA
Fixed code retrieving the remnant nucleus which stopped working as soon as
simulation of nuclear de-excitation started pushing photons in the target
nucleus daughter list.
@ Jun 20, 2008 - CA
- Fix a mem leak: The (clone of the) GHepParticle being re-scattered was not
+ Fix a mem leak: The (clone of the) GHepParticle being re-scattered was not
deleted after it was added at the GHEP event record.
@ Jul 15, 2010 - AM
The hN mode is now implemented in Intranuke. Similar to hA mode, but particles
produced by reactions are stepped through the nucleus like probe particles.
Particles react with nucleons instead of the entire nucleus, and final states
are determined after reactions are finished, not before.
- @ Dec 15, 2014 - SD, Nick Geary
+ @ Dec 15, 2014 - SD, Nick Geary
Update fates to include Compound Nucleus final state correctly.
@ Jan 9, 2015 - SD, NG, Tomek Golan
Added 2014 version of INTRANUKE codes (new class) for independent development.
@@ -113,9 +113,9 @@ HNIntranuke2025::~HNIntranuke2025()
//___________________________________________________________________________
void HNIntranuke2025::ProcessEventRecord(GHepRecord * evrec) const
{
- LOG("HNIntranuke2025", pNOTICE)
+ LOG("HNIntranuke2025", pNOTICE)
<< "************ Running hN2025 MODE INTRANUKE ************";
-
+
Intranuke2025::ProcessEventRecord(evrec);
LOG("HNIntranuke2025", pINFO) << "Done with this event";
@@ -136,7 +136,7 @@ void HNIntranuke2025::SimulateHadronicFinalState(GHepRecord* ev, GHepParticle* p
bool is_pion = (pdgc==kPdgPiP || pdgc==kPdgPiM || pdgc==kPdgPi0);
bool is_kaon = (pdgc==kPdgKP);
bool is_baryon = (pdgc==kPdgProton || pdgc==kPdgNeutron);
- bool is_gamma = (pdgc==kPdgGamma);
+ bool is_gamma = (pdgc==kPdgGamma);
if(!(is_pion || is_baryon || is_gamma || is_kaon))
{
LOG("HNIntranuke2025", pERROR) << "** Cannot handle particle: " << p->Name();
@@ -153,7 +153,7 @@ void HNIntranuke2025::SimulateHadronicFinalState(GHepRecord* ev, GHepParticle* p
if(fate == kIHNFtUndefined)
{
LOG("HNIntranuke2025", pERROR) << "** Couldn't select a fate";
- LOG("HNIntranuke2025", pERROR) << "** Num Protons: " << fRemnZ
+ LOG("HNIntranuke2025", pERROR) << "** Num Protons: " << fRemnZ
<< ", Num Neutrons: "<<(fRemnA-fRemnZ);
LOG("HNIntranuke2025", pERROR) << "** Particle: " << "\n" << (*p);
//LOG("HNIntranuke2025", pERROR) << "** Event Record: " << "\n" << (*ev);
@@ -172,7 +172,7 @@ void HNIntranuke2025::SimulateHadronicFinalState(GHepRecord* ev, GHepParticle* p
this->ElasHN(ev,p,fate);
}
else if(fate == kIHNFtAbs) {this-> AbsorbHN(ev,p,fate);}
- else if(fate == kIHNFtInelas && pdgc != kPdgGamma)
+ else if(fate == kIHNFtInelas && pdgc != kPdgGamma)
{
#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
LOG("HNIntranuke2025", pDEBUG)
@@ -195,7 +195,7 @@ void HNIntranuke2025::SimulateHadronicFinalState(GHepRecord* ev, GHepParticle* p
catch(exceptions::INukeException exception)
{
this->SimulateHadronicFinalState(ev,p);
- LOG("HNIntranuke2025", pNOTICE)
+ LOG("HNIntranuke2025", pNOTICE)
<< "retry call to SimulateHadronicFinalState ";
LOG("HNIntranuke2025", pNOTICE) << exception;
@@ -215,8 +215,8 @@ INukeFateHN_t HNIntranuke2025::HadronFateHN(const GHepParticle * p) const
bool isPion = (pdgc == kPdgPiP or pdgc == kPdgPi0 or pdgc == kPdgPiM);
if (isPion and fUseOset and ke < 350.0) return HadronFateOset ();
-
- LOG("HNIntranuke2025", pNOTICE)
+
+ LOG("HNIntranuke2025", pNOTICE)
<< "Selecting hN fate for " << p->Name() << " with KE = " << ke << " MeV";
// try to generate a hadron fate
@@ -241,7 +241,7 @@ INukeFateHN_t HNIntranuke2025::HadronFateHN(const GHepParticle * p) const
frac_elas *= fNucQEFac;
if(pdgc==kPdgPi0) frac_abs*= 0.665; //isospin factor
- LOG("HNIntranuke2025", pNOTICE)
+ LOG("HNIntranuke2025", pNOTICE)
<< "\n frac{" << INukeHadroFates2025::AsString(kIHNFtCEx) << "} = " << frac_cex
<< "\n frac{" << INukeHadroFates2025::AsString(kIHNFtElas) << "} = " << frac_elas
<< "\n frac{" << INukeHadroFates2025::AsString(kIHNFtInelas) << "} = " << frac_inel
@@ -250,7 +250,7 @@ INukeFateHN_t HNIntranuke2025::HadronFateHN(const GHepParticle * p) const
// compute total fraction (can be <1 if fates have been switched off)
double tf = frac_cex +
frac_elas +
- frac_inel +
+ frac_inel +
frac_abs;
double r = tf * rnd->RndFsi().Rndm();
@@ -263,10 +263,10 @@ INukeFateHN_t HNIntranuke2025::HadronFateHN(const GHepParticle * p) const
if(r < (cf += frac_cex )) return kIHNFtCEx; //cex
if(r < (cf += frac_elas )) return kIHNFtElas; //elas
if(r < (cf += frac_inel )) return kIHNFtInelas; //inelas
- if(r < (cf += frac_abs )) return kIHNFtAbs; //abs
+ if(r < (cf += frac_abs )) return kIHNFtAbs; //abs
- LOG("HNIntranuke2025", pWARN)
- << "No selection after going through all fates! "
+ LOG("HNIntranuke2025", pWARN)
+ << "No selection after going through all fates! "
<< "Total fraction = " << tf << " (r = " << r << ")";
////////////////////////////
return kIHNFtUndefined;
@@ -282,7 +282,7 @@ INukeFateHN_t HNIntranuke2025::HadronFateHN(const GHepParticle * p) const
double frac_cmp = this->FateWeight(pdgc, kIHNFtCmp)
* fHadroData2025->Frac(pdgc, kIHNFtCmp, ke, fRemnA , fRemnZ);
- LOG("HNIntranuke2025", pINFO)
+ LOG("HNIntranuke2025", pINFO)
<< "\n frac{" << INukeHadroFates2025::AsString(kIHNFtElas) << "} = " << frac_elas
<< "\n frac{" << INukeHadroFates2025::AsString(kIHNFtInelas) << "} = " << frac_inel;
@@ -302,14 +302,14 @@ INukeFateHN_t HNIntranuke2025::HadronFateHN(const GHepParticle * p) const
if(r < (cf += frac_inel )) return kIHNFtInelas; // inelas
if(r < (cf += frac_cmp )) return kIHNFtCmp; // cmp
- LOG("HNIntranuke2025", pWARN)
+ LOG("HNIntranuke2025", pWARN)
<< "No selection after going through all fates! "
<< "Total fraction = " << tf << " (r = " << r << ")";
//////////////////////////
return kIHNFtUndefined;
}
- // handle gamma -- does not currently consider the elastic case
+ // handle gamma -- does not currently consider the elastic case
else if (pdgc==kPdgGamma) return kIHNFtInelas;
// Handle kaon -- elastic + charge exchange
else if (pdgc==kPdgKP){
@@ -321,7 +321,7 @@ INukeFateHN_t HNIntranuke2025::HadronFateHN(const GHepParticle * p) const
// frac_cex *= fNucCEXFac; // scaling factors
// frac_elas *= fNucQEFac; // Flor - Correct scaling factors?
- LOG("HNIntranuke", pINFO)
+ LOG("HNIntranuke", pINFO)
<< "\n frac{" << INukeHadroFates2025::AsString(kIHNFtCEx) << "} = " << frac_cex
<< "\n frac{" << INukeHadroFates2025::AsString(kIHNFtElas) << "} = " << frac_elas;
@@ -337,10 +337,10 @@ INukeFateHN_t HNIntranuke2025::HadronFateHN(const GHepParticle * p) const
double cf=0; // current fraction
if(r < (cf += frac_cex )) return kIHNFtCEx; //cex
- if(r < (cf += frac_elas )) return kIHNFtElas; //elas
+ if(r < (cf += frac_elas )) return kIHNFtElas; //elas
- LOG("HNIntranuke", pWARN)
- << "No selection after going through all fates! "
+ LOG("HNIntranuke", pWARN)
+ << "No selection after going through all fates! "
<< "Total fraction = " << tf << " (r = " << r << ")";
////////////////////////////
return kIHNFtUndefined;
@@ -364,7 +364,7 @@ double HNIntranuke2025::FateWeight(int pdgc, INukeFateHN_t fate) const
int np = fRemnZ;
int nn = fRemnA - fRemnZ;
-
+
if (np < 1 && nn < 1)
{
LOG("HNIntranuke2025", pERROR) << "** Nothing left in nucleus!!! **";
@@ -386,7 +386,7 @@ void HNIntranuke2025::AbsorbHN(
GHepRecord * ev, GHepParticle * p, INukeFateHN_t fate) const
{
// handles pi+d->2p, pi-d->nn, pi0 d->pn absorbtion, all using pi+d values
-
+
int pdgc = p->Pdg();
#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
@@ -434,7 +434,7 @@ void HNIntranuke2025::AbsorbHN(
// Library instance for reference
PDGLibrary * pLib = PDGLibrary::Instance();
-
+
// Handle fermi target
Target target(ev->TargetNucleus()->Pdg());
@@ -470,7 +470,7 @@ void HNIntranuke2025::AbsorbHN(
<< "AbsorbHN() cannot handle probe: " << pdgc;
return;
}
-
+
// assign proper masses
M1 = pLib->Find(pcode) ->Mass();
M2_1 = pLib->Find(t1code)->Mass();
@@ -478,14 +478,14 @@ void HNIntranuke2025::AbsorbHN(
M3 = pLib->Find(scode) ->Mass();
M4 = pLib->Find(s2code)->Mass();
- // handle fermi momentum
+ // handle fermi momentum
if(fDoFermi)
{
target.SetHitNucPdg(t1code);
fNuclmodel->GenerateNucleon(target);
tP2_1L=fFermiFac * fNuclmodel->Momentum3();
E2_1L = TMath::Sqrt(tP2_1L.Mag2() + M2_1*M2_1);
-
+
target.SetHitNucPdg(t2code);
fNuclmodel->GenerateNucleon(target);
tP2_2L=fFermiFac * fNuclmodel->Momentum3();
@@ -505,7 +505,7 @@ void HNIntranuke2025::AbsorbHN(
// adjust p to reflect scattering
// get random scattering angle
C3CM = fHadroData2025->IntBounce(p,t1code,scode,fate);
- if (C3CM<-1.)
+ if (C3CM<-1.)
{
p->SetStatus(kIStStableFinalState);
ev->AddParticle(*p);
@@ -595,7 +595,7 @@ void HNIntranuke2025::AbsorbHN(
p->SetStatus(kIStHadronInTheNucleus);
//disable until needed
// utils::intranuke2025::StepParticle(p,fFreeStep,fTrackingRadius);
- ev->AddParticle(*p);
+ ev->AddParticle(*p);
return;
*/
// new attempt at error handling:
@@ -613,12 +613,12 @@ void HNIntranuke2025::AbsorbHN(
// get random phi angle, distributed uniformally in 360 deg
PHI3 = 2 * kPi * rnd->RndFsi().Rndm();
-
+
tP3L = P3zL*bDir + P3tL*tTrans;
tP4L = P4zL*bDir + P4tL*tTrans;
tP3L.Rotate(PHI3,bDir); // randomize transverse components
- tP4L.Rotate(PHI3,bDir);
+ tP4L.Rotate(PHI3,bDir);
E3L = TMath::Sqrt(P3L*P3L + M3*M3);
E4L = TMath::Sqrt(P4L*P4L + M4*M4);
@@ -705,7 +705,7 @@ void HNIntranuke2025::ElasHN(
// get random scattering angle
double C3CM = fHadroData2025->IntBounce(p,tcode,scode,fate);
- if (C3CM<-1.)
+ if (C3CM<-1.)
{
p->SetStatus(kIStStableFinalState);
ev->AddParticle(*p);
@@ -718,7 +718,7 @@ void HNIntranuke2025::ElasHN(
double Mt = t->Mass();
//t->SetMomentum(TLorentzVector(0,0,0,Mt));
t->SetRemovalEnergy(0);
- // handle fermi momentum
+ // handle fermi momentum
if(fDoFermi)
{
// Handle fermi target
@@ -765,24 +765,24 @@ void HNIntranuke2025::ElasHN(
void HNIntranuke2025::InelasticHN(GHepRecord* ev, GHepParticle* p) const
{
// Aaron Meyer (Jan 2010)
- // Updated version of InelasticHN
+ // Updated version of InelasticHN
- GHepParticle s1(*p);
+ GHepParticle s1(*p);
GHepParticle s2(*p);
GHepParticle s3(*p);
s2.SetRemovalEnergy(0);
s3.SetRemovalEnergy(0);
-
-
-
- if (utils::intranuke2025::PionProduction(ev,p,&s1,&s2,&s3,fRemnA,fRemnZ,fRemnP4,fDoFermi,fFermiFac,fFermiMomentum,fNuclmodel))
+
+
+
+ if (utils::intranuke2025::PionProduction(ev,p,&s1,&s2,&s3,fRemnA,fRemnZ,fRemnP4,fDoFermi,fFermiFac,fFermiMomentum,fNuclmodel,fPiProdThreeBodyBias))
{
// set status of particles and return
-
+
s1.SetStatus(kIStHadronInTheNucleus);
s2.SetStatus(kIStHadronInTheNucleus);
s3.SetStatus(kIStHadronInTheNucleus);
-
+
ev->AddParticle(s1);
ev->AddParticle(s2);
ev->AddParticle(s3);
@@ -798,7 +798,7 @@ void HNIntranuke2025::InelasticHN(GHepRecord* ev, GHepParticle* p) const
}
//___________________________________________________________________________
-void HNIntranuke2025::GammaInelasticHN(GHepRecord* ev, GHepParticle* p, INukeFateHN_t fate) const
+void HNIntranuke2025::GammaInelasticHN(GHepRecord* ev, GHepParticle* p, INukeFateHN_t fate) const
{
// This function handles pion photoproduction reactions
@@ -847,7 +847,7 @@ void HNIntranuke2025::GammaInelasticHN(GHepRecord* ev, GHepParticle* p, INukeFat
<< "Error: could not determine particle final states";
ev->AddParticle(*p);
return;
- }
+ }
LOG("HNIntranuke2025", pNOTICE)
<< "GammaInelastic fate: " << INukeHadroFates2025::AsString(fate);
@@ -860,7 +860,7 @@ void HNIntranuke2025::GammaInelasticHN(GHepRecord* ev, GHepParticle* p, INukeFat
t->SetPdgCode(tcode);
double Mt = t->Mass();
- // handle fermi momentum
+ // handle fermi momentum
if(fDoFermi)
{
// Handle fermi target
@@ -907,7 +907,7 @@ int HNIntranuke2025::HandleCompoundNucleus(GHepRecord* ev, GHepParticle* p, int
// handle compound nucleus option
// -- Call the PreEquilibrium function
- if( fDoCompoundNucleus && IsInNucleus(p) && pdg::IsNeutronOrProton(p->Pdg()))
+ if( fDoCompoundNucleus && IsInNucleus(p) && pdg::IsNeutronOrProton(p->Pdg()))
{ // random number generator
//unused var - quiet compiler warning//RandomGen * rnd = RandomGen::Instance();
@@ -983,6 +983,8 @@ void HNIntranuke2025::LoadConfig(void)
GetParamDef( "FSI-NeutralPion-MFPScale", fNeutralPionMFPScale, 1.0 ) ;
GetParamDef( "FSI-Nucleon-MFPScale", fNucleonMFPScale, 1.0 ) ;
+ GetParamDef( "FSI-PiProd-ThreeBodyBias", fPiProdThreeBodyBias, 0.0 );
+
// report
LOG("HNIntranuke2025", pINFO) << "Settings for Intranuke2025 mode: " << INukeMode::AsString(kIMdHN);
LOG("HNIntranuke2025", pWARN) << "R0 = " << fR0 << " fermi";
@@ -1003,6 +1005,7 @@ void HNIntranuke2025::LoadConfig(void)
LOG("HNIntranuke2025", pWARN) << "XsecNNCorr? = " << ((fXsecNNCorr)?(true):(false));
LOG("HNIntranuke2025", pWARN) << "FSI-ChargedPion-MFPScale = " << fChPionMFPScale;
LOG("HNIntranuke2025", pWARN) << "FSI-NeutralPion-MFPScale = " << fNeutralPionMFPScale;
+ LOG("HNIntranuke2025", pWARN) << "PiProdBias = " << fPiProdThreeBodyBias;
}
//___________________________________________________________________________
diff --git a/src/Physics/HadronTransport/HNIntranuke2025.h b/src/Physics/HadronTransport/HNIntranuke2025.h
index c629b9074..928c948b2 100644
--- a/src/Physics/HadronTransport/HNIntranuke2025.h
+++ b/src/Physics/HadronTransport/HNIntranuke2025.h
@@ -23,7 +23,7 @@
\cpright Copyright (c) 2003-2025, The GENIE Collaboration
For the full text of the license visit http://copyright.genie-mc.org
-
+
*/
//____________________________________________________________________________
@@ -48,7 +48,7 @@ class GHepParticle;
class INukeHadroData;
class PDGCodeList;
-class HNIntranuke2025 : public Intranuke2025 {
+class HNIntranuke2025 : public Intranuke2025 {
friend class IntranukeTester;
@@ -74,9 +74,9 @@ public :
void ElasHN (GHepRecord* ev, GHepParticle* p, INukeFateHN_t fate) const;
void AbsorbHN (GHepRecord* ev, GHepParticle* p, INukeFateHN_t fate) const;
void InelasticHN (GHepRecord* ev, GHepParticle* p) const;
- void GammaInelasticHN (GHepRecord* ev, GHepParticle* p, INukeFateHN_t fate) const;
+ void GammaInelasticHN (GHepRecord* ev, GHepParticle* p, INukeFateHN_t fate) const;
bool HandleCompoundNucleusHN (GHepRecord* ev, GHepParticle* p) const;
- int HandleCompoundNucleus(GHepRecord* ev, GHepParticle* p, int mom) const;
+ int HandleCompoundNucleus(GHepRecord* ev, GHepParticle* p, int mom) const;
mutable int nuclA; ///< value of A for the target nucleus in hA mode
diff --git a/src/Physics/HadronTransport/INukeUtils2018.cxx b/src/Physics/HadronTransport/INukeUtils2018.cxx
index 773e1c722..cb4b40535 100644
--- a/src/Physics/HadronTransport/INukeUtils2018.cxx
+++ b/src/Physics/HadronTransport/INukeUtils2018.cxx
@@ -1071,7 +1071,7 @@ bool genie::utils::intranuke2018::TwoBodyKinematics(
//___________________________________________________________________________
bool genie::utils::intranuke2018::ThreeBodyKinematics(
GHepRecord* ev, GHepParticle* p, int tcode, GHepParticle* s1, GHepParticle* s2, GHepParticle* s3,
- bool DoFermi, double FermiFac, double FermiMomentum, const NuclearModelI* Nuclmodel)
+ bool DoFermi, double FermiFac, double FermiMomentum, const NuclearModelI* Nuclmodel, double bias)
{
// Aaron Meyer (7/15/10)
@@ -1093,7 +1093,7 @@ bool genie::utils::intranuke2018::ThreeBodyKinematics(
double P3zL, P4zL, P4tL, P5zL, P5tL;
double Et, M, theta1, theta2;
double P1zL, P2zL;
- double theta3, theta4, phi3, phi4, theta5;
+ double costheta3, costheta4, phi3, phi4, theta5;
TVector3 tP2L, tP1L, tPtot, tbeta, tbetadir, tTrans, tP4L, tP5L;
TVector3 tP1zCM, tP2zCM, tP3L, tPiL, tbeta2, tbetadir2, tVect, tTrans2;
@@ -1154,105 +1154,132 @@ bool genie::utils::intranuke2018::ThreeBodyKinematics(
E2CM = gm*E2L - gm*beta*P2zL;
tP2zCM = gm*P2zL*tbetadir - gm*tbeta*E2L;
Et = E1CM + E2CM;
- M = (rnd->RndFsi().Rndm()*(Et - M3 - M4 - M5)) + (M4 + M5);
- E3CM = (Et*Et + M3*M3 - M*M)/(2*Et);
- EiCM = Et - E3CM;
- if(E3CM*E3CM - M3*M3<0)
- {
- LOG("INukeUtils",pNOTICE)
- << "PionProduction P3 has non-real momentum - retry kinematics";
- LOG("INukeUtils",pNOTICE) << "Energy, masses of 3 fs particales:"
- << E3CM << " " << M3 << " " << " " << M4 << " " << M5;
- exceptions::INukeException exception;
- exception.SetReason("PionProduction particle 3 has non-real momentum");
- throw exception;
- return false;
- }
- P3CM = TMath::Sqrt(E3CM*E3CM - M3*M3);
-
- theta3 = kPi * rnd->RndFsi().Rndm();
- theta4 = kPi * rnd->RndFsi().Rndm();
- phi3 = 2*kPi * rnd->RndFsi().Rndm();
- phi4 = 2*kPi * rnd->RndFsi().Rndm();
-
- P3zL = gm*beta*E3CM + gm*P3CM*TMath::Cos(theta3);
- P3tL = P3CM*TMath::Sin(theta3);
- PizL = gm*beta*EiCM - gm*P3CM*TMath::Cos(theta3);
- PitL = -P3CM*TMath::Sin(theta3);
- P3L = TMath::Sqrt(P3zL*P3zL + P3tL*P3tL);
- PiL = TMath::Sqrt(PizL*PizL + PitL*PitL);
- E3L = TMath::Sqrt(P3L*P3L + M3*M3);
- EiL = TMath::Sqrt(PiL*PiL + M*M);
-
- // handle very low momentum particles
- if(!(TMath::Finite(P3L)) || P3L < .001)
+ // G.P. 2/20/2025
+ // Sample uniformly in lorentz invariany phase space
+ //
+ // Just because you sample uniformly does not mean you are
+ // uniformly sampling Lorentz invariant phase space.
+ //
+ // D-LIPS taken from PDG: https://pdg.lbl.gov/2018/reviews/rpp2018-rev-kinematics.pdf
+ // Algorithm is morally similar to Raubdo-Lynch method, hard-coded to 3 bodies
+
+ // compute naive max-weight
+ double E3CM_max = (Et*Et + M3*M3 - (M4+M5)*(M4+M5))/(2*Et);
+ double P3CM_max = TMath::Sqrt(E3CM_max*E3CM_max - M3*M3);
+ double E4CM2_max = ((Et-M3)*(Et-M3) + M4*M4 - M5*M5) / (2*(Et-M3));
+ double P4CM2_max = TMath::Sqrt(E4CM2_max*E4CM2_max - M4*M4);
+ double max_weight = P3CM_max * P4CM2_max;
+ double weight = 0.;
+
+ do {
+ M = (rnd->RndFsi().Rndm()*(Et - M3 - M4 - M5)) + (M4 + M5);
+ E3CM = (Et*Et + M3*M3 - M*M)/(2*Et);
+ EiCM = Et - E3CM;
+ if(E3CM*E3CM - M3*M3<0)
{
- LOG("INukeUtils",pINFO)
- << "Particle 3 " << M3 << " momentum small or non-finite: " << P3L
- << "\n" << "--> Assigning .001 as new momentum";
- P3tL = 0;
- P3zL = .001;
- P3L = .001;
- E3L = TMath::Sqrt(P3L*P3L + M3*M3);
+ LOG("INukeUtils",pNOTICE)
+ << "PionProduction P3 has non-real momentum - retry kinematics";
+ LOG("INukeUtils",pNOTICE) << "Energy, masses of 3 fs particales:"
+ << E3CM << " " << M3 << " " << " " << M4 << " " << M5;
+ exceptions::INukeException exception;
+ exception.SetReason("PionProduction particle 3 has non-real momentum");
+ throw exception;
+ return false;
}
+ P3CM = TMath::Sqrt(E3CM*E3CM - M3*M3);
+
+ costheta3 = 2*rnd->RndFsi().Rndm()-1;
+ costheta4 = 2*rnd->RndFsi().Rndm()-1;
+ double sintheta3 = TMath::Sqrt(1 - costheta3*costheta3);
+ double sintheta4 = TMath::Sqrt(1 - costheta4*costheta4);
+ phi3 = 2*kPi * rnd->RndFsi().Rndm();
+ phi4 = 2*kPi * rnd->RndFsi().Rndm();
+
+ P3zL = gm*beta*E3CM + gm*P3CM*costheta3;
+ P3tL = P3CM*sintheta3;
+ PizL = gm*beta*EiCM - gm*P3CM*costheta3;
+ PitL = -P3CM*sintheta3;
+
+ P3L = TMath::Sqrt(P3zL*P3zL + P3tL*P3tL);
+ PiL = TMath::Sqrt(PizL*PizL + PitL*PitL);
+ E3L = TMath::Sqrt(P3L*P3L + M3*M3);
+ EiL = TMath::Sqrt(PiL*PiL + M*M);
+
+ // handle very low momentum particles
+ if(!(TMath::Finite(P3L)) || P3L < .001)
+ {
+ LOG("INukeUtils",pINFO)
+ << "Particle 3 " << M3 << " momentum small or non-finite: " << P3L
+ << "\n" << "--> Assigning .001 as new momentum";
+ P3tL = 0;
+ P3zL = .001;
+ P3L = .001;
+ E3L = TMath::Sqrt(P3L*P3L + M3*M3);
+ }
- tP3L = P3zL*tbetadir + P3tL*tTrans;
- tPiL = PizL*tbetadir + PitL*tTrans;
- tP3L.Rotate(phi3,tbetadir);
- tPiL.Rotate(phi3,tbetadir);
-
- // second sequence, handle formally composite particles 4 and 5
- tbeta2 = tPiL * (1.0 / EiL);
- tbetadir2 = tbeta2.Unit();
- beta2 = tbeta2.Mag();
- gm2 = 1.0 / TMath::Sqrt(1.0 - beta2*beta2);
-
- E4CM2 = (M*M + M4*M4 - M5*M5) / (2*M);
- E5CM2 = M - E4CM2;
- P4CM2 = TMath::Sqrt(E4CM2*E4CM2 - M4*M4);
-
- tVect.SetXYZ(1,0,0);
- if(TMath::Abs((tVect - tbetadir2).Mag())<.01) tVect.SetXYZ(0,1,0);
- theta5 = tVect.Angle(tbetadir2);
- tTrans2 = (tVect - TMath::Cos(theta5)*tbetadir2).Unit();
+ tP3L = P3zL*tbetadir + P3tL*tTrans;
+ tPiL = PizL*tbetadir + PitL*tTrans;
+ tP3L.Rotate(phi3,tbetadir);
+ tPiL.Rotate(phi3,tbetadir);
+
+ // second sequence, handle formally composite particles 4 and 5
+ tbeta2 = tPiL * (1.0 / EiL);
+ tbetadir2 = tbeta2.Unit();
+ beta2 = tbeta2.Mag();
+ gm2 = 1.0 / TMath::Sqrt(1.0 - beta2*beta2);
+
+ E4CM2 = (M*M + M4*M4 - M5*M5) / (2*M);
+ E5CM2 = M - E4CM2;
+ P4CM2 = TMath::Sqrt(E4CM2*E4CM2 - M4*M4);
+
+ tVect.SetXYZ(1,0,0);
+ if(TMath::Abs((tVect - tbetadir2).Mag())<.01) tVect.SetXYZ(0,1,0);
+ theta5 = tVect.Angle(tbetadir2);
+ tTrans2 = (tVect - TMath::Cos(theta5)*tbetadir2).Unit();
+
+ P4zL = gm2*beta2*E4CM2 + gm2*P4CM2*costheta4;
+ P4tL = P4CM2*sintheta4;
+ P5zL = gm2*beta2*E5CM2 - gm2*P4CM2*costheta4;
+ P5tL = - P4tL;
+
+ P4L = TMath::Sqrt(P4zL*P4zL + P4tL*P4tL);
+ P5L = TMath::Sqrt(P5zL*P5zL + P5tL*P5tL);
+ E4L = TMath::Sqrt(P4L*P4L + M4*M4);
+ E5L = TMath::Sqrt(P5L*P5L + M5*M5);
+
+ // handle very low momentum particles
+ if(!(TMath::Finite(P4L)) || P4L < .001)
+ {
+ LOG("INukeUtils",pINFO)
+ << "Particle 4 " << M4 << " momentum small or non-finite: " << P4L
+ << "\n" << "--> Assigning .001 as new momentum";
+ P4tL = 0;
+ P4zL = .001;
+ P4L = .001;
+ E4L = TMath::Sqrt(P4L*P4L + M4*M4);
+ }
+ if(!(TMath::Finite(P5L)) || P5L < .001)
+ {
+ LOG("INukeUtils",pINFO)
+ << "Particle 5 " << M5 << " momentum small or non-finite: " << P5L
+ << "\n" << "--> Assigning .001 as new momentum";
+ P5tL = 0;
+ P5zL = .001;
+ P5L = .001;
+ E5L = TMath::Sqrt(P5L*P5L + M5*M5);
+ }
- P4zL = gm2*beta2*E4CM2 + gm2*P4CM2*TMath::Cos(theta4);
- P4tL = P4CM2*TMath::Sin(theta4);
- P5zL = gm2*beta2*E5CM2 - gm2*P4CM2*TMath::Cos(theta4);
- P5tL = - P4tL;
+ tP4L = P4zL*tbetadir2 + P4tL*tTrans2;
+ tP5L = P5zL*tbetadir2 + P5tL*tTrans2;
+ tP4L.Rotate(phi4,tbetadir2);
+ tP5L.Rotate(phi4,tbetadir2);
- P4L = TMath::Sqrt(P4zL*P4zL + P4tL*P4tL);
- P5L = TMath::Sqrt(P5zL*P5zL + P5tL*P5tL);
- E4L = TMath::Sqrt(P4L*P4L + M4*M4);
- E5L = TMath::Sqrt(P5L*P5L + M5*M5);
+ weight = P3CM*P4CM2;
- // handle very low momentum particles
- if(!(TMath::Finite(P4L)) || P4L < .001)
- {
- LOG("INukeUtils",pINFO)
- << "Particle 4 " << M4 << " momentum small or non-finite: " << P4L
- << "\n" << "--> Assigning .001 as new momentum";
- P4tL = 0;
- P4zL = .001;
- P4L = .001;
- E4L = TMath::Sqrt(P4L*P4L + M4*M4);
- }
- if(!(TMath::Finite(P5L)) || P5L < .001)
- {
- LOG("INukeUtils",pINFO)
- << "Particle 5 " << M5 << " momentum small or non-finite: " << P5L
- << "\n" << "--> Assigning .001 as new momentum";
- P5tL = 0;
- P5zL = .001;
- P5L = .001;
- E5L = TMath::Sqrt(P5L*P5L + M5*M5);
- }
+ if (bias != 0) weight *= TMath::Exp(bias*(TLorentzVector(tP3L,E3L) - *p->P4()).M2());
- tP4L = P4zL*tbetadir2 + P4tL*tTrans2;
- tP5L = P5zL*tbetadir2 + P5tL*tTrans2;
- tP4L.Rotate(phi4,tbetadir2);
- tP5L.Rotate(phi4,tbetadir2);
+ } while (rnd->RndFsi().Rndm() > weight/max_weight);
// pauli blocking
if(P3L < FermiMomentum || ( pdg::IsNeutronOrProton(s2->Pdg()) && P4L < FermiMomentum ) )
@@ -1300,7 +1327,8 @@ bool genie::utils::intranuke2018::ThreeBodyKinematics(
//___________________________________________________________________________
bool genie::utils::intranuke2018::PionProduction(
GHepRecord* ev, GHepParticle* p, GHepParticle* s1, GHepParticle* s2, GHepParticle* s3, int &RemnA, int &RemnZ,
- TLorentzVector &RemnP4, bool DoFermi, double FermiFac, double FermiMomentum, const NuclearModelI* Nuclmodel)
+ TLorentzVector &RemnP4, bool DoFermi, double FermiFac, double FermiMomentum, const NuclearModelI* Nuclmodel,
+ double bias)
{
// Aaron Meyer (7/15/2010)
//
@@ -1620,7 +1648,7 @@ bool genie::utils::intranuke2018::PionProduction(
s3->SetPdgCode(p5code);
if(genie::utils::intranuke2018::ThreeBodyKinematics(
- ev,p,(ptarg?kPdgProton:kPdgNeutron),s1,s2,s3,DoFermi,FermiFac,FermiMomentum,Nuclmodel))
+ ev,p,(ptarg?kPdgProton:kPdgNeutron),s1,s2,s3,DoFermi,FermiFac,FermiMomentum,Nuclmodel,bias))
{
// okay, handle remnants and return true
// assumes first particle is always the nucleon,
diff --git a/src/Physics/HadronTransport/INukeUtils2018.h b/src/Physics/HadronTransport/INukeUtils2018.h
index dbbb00c5e..00399f062 100644
--- a/src/Physics/HadronTransport/INukeUtils2018.h
+++ b/src/Physics/HadronTransport/INukeUtils2018.h
@@ -89,11 +89,12 @@ namespace intranuke2018
bool ThreeBodyKinematics(
GHepRecord* ev, GHepParticle* p, int tcode, GHepParticle* s1, GHepParticle* s2, GHepParticle* s3,
- bool DoFermi=false, double FermiFac=0, double FermiMomentum=0, const NuclearModelI* Nuclmodel=(const NuclearModelI*)0);
+ bool DoFermi=false, double FermiFac=0, double FermiMomentum=0, const NuclearModelI* Nuclmodel=(const NuclearModelI*)0, double bias=0);
bool PionProduction(
GHepRecord* ev, GHepParticle* p, GHepParticle* s1, GHepParticle* s2, GHepParticle* s3, int &RemnA, int &RemnZ,
- TLorentzVector &RemnP4,bool DoFermi, double FermiFac, double FermiMomentum, const NuclearModelI* Nuclmodel);
+ TLorentzVector &RemnP4,bool DoFermi, double FermiFac, double FermiMomentum, const NuclearModelI* Nuclmodel,
+ double bias);
double CalculateEta(
double Minc, double ke, double Mtarg, double Mtwopart, double Mpi);
diff --git a/src/Physics/HadronTransport/INukeUtils2025.cxx b/src/Physics/HadronTransport/INukeUtils2025.cxx
index b1951bbbf..f927ff7e8 100644
--- a/src/Physics/HadronTransport/INukeUtils2025.cxx
+++ b/src/Physics/HadronTransport/INukeUtils2025.cxx
@@ -62,7 +62,7 @@ using namespace genie::controls;
//____________________________________________________________________________
double genie::utils::intranuke2025::MeanFreePath(
int pdgc, const TLorentzVector & x4, const TLorentzVector & p4,
- double A, double Z, double nRpi, double nRnuc, const bool useOset,
+ double A, double Z, double nRpi, double nRnuc, const bool useOset,
const bool altOset, const bool xsecNNCorr, string INukeMode)
{
// Calculate the mean free path (in fm) for a pions and nucleons in a nucleus
@@ -1056,7 +1056,7 @@ bool genie::utils::intranuke2025::TwoBodyKinematics(
//___________________________________________________________________________
bool genie::utils::intranuke2025::ThreeBodyKinematics(
GHepRecord* ev, GHepParticle* p, int tcode, GHepParticle* s1, GHepParticle* s2, GHepParticle* s3,
- bool DoFermi, double FermiFac, double FermiMomentum, const NuclearModelI* Nuclmodel)
+ bool DoFermi, double FermiFac, double FermiMomentum, const NuclearModelI* Nuclmodel, double bias)
{
// Aaron Meyer (7/15/10)
@@ -1078,7 +1078,7 @@ bool genie::utils::intranuke2025::ThreeBodyKinematics(
double P3zL, P4zL, P4tL, P5zL, P5tL;
double Et, M, theta1, theta2;
double P1zL, P2zL;
- double theta3, theta4, phi3, phi4, theta5;
+ double costheta3, costheta4, phi3, phi4, theta5;
TVector3 tP2L, tP1L, tPtot, tbeta, tbetadir, tTrans, tP4L, tP5L;
TVector3 tP1zCM, tP2zCM, tP3L, tPiL, tbeta2, tbetadir2, tVect, tTrans2;
@@ -1139,105 +1139,134 @@ bool genie::utils::intranuke2025::ThreeBodyKinematics(
E2CM = gm*E2L - gm*beta*P2zL;
tP2zCM = gm*P2zL*tbetadir - gm*tbeta*E2L;
Et = E1CM + E2CM;
- M = (rnd->RndFsi().Rndm()*(Et - M3 - M4 - M5)) + (M4 + M5);
- E3CM = (Et*Et + M3*M3 - M*M)/(2*Et);
- EiCM = Et - E3CM;
- if(E3CM*E3CM - M3*M3<0)
- {
- LOG("INukeUtils",pNOTICE)
- << "PionProduction P3 has non-real momentum - retry kinematics";
- LOG("INukeUtils",pNOTICE) << "Energy, masses of 3 fs particales:"
- << E3CM << " " << M3 << " " << " " << M4 << " " << M5;
- exceptions::INukeException exception;
- exception.SetReason("PionProduction particle 3 has non-real momentum");
- throw exception;
- return false;
- }
- P3CM = TMath::Sqrt(E3CM*E3CM - M3*M3);
-
- theta3 = kPi * rnd->RndFsi().Rndm();
- theta4 = kPi * rnd->RndFsi().Rndm();
- phi3 = 2*kPi * rnd->RndFsi().Rndm();
- phi4 = 2*kPi * rnd->RndFsi().Rndm();
-
- P3zL = gm*beta*E3CM + gm*P3CM*TMath::Cos(theta3);
- P3tL = P3CM*TMath::Sin(theta3);
- PizL = gm*beta*EiCM - gm*P3CM*TMath::Cos(theta3);
- PitL = -P3CM*TMath::Sin(theta3);
- P3L = TMath::Sqrt(P3zL*P3zL + P3tL*P3tL);
- PiL = TMath::Sqrt(PizL*PizL + PitL*PitL);
- E3L = TMath::Sqrt(P3L*P3L + M3*M3);
- EiL = TMath::Sqrt(PiL*PiL + M*M);
-
- // handle very low momentum particles
- if(!(TMath::Finite(P3L)) || P3L < .001)
+ // G.P. 2/20/2025
+ // Sample uniformly in lorentz invariant phase space
+ //
+ // Just because you sample uniformly does not mean you are
+ // uniformly sampling Lorentz invariant phase space.
+ //
+ // D-LIPS taken from PDG:
+ // https://pdg.lbl.gov/2018/reviews/rpp2018-rev-kinematics.pdf
+ // Algorithm is morally similar to Raubdo-Lynch method, hard-coded to 3 bodies
+
+ // compute naive max-weight
+ double E3CM_max = (Et*Et + M3*M3 - (M4+M5)*(M4+M5))/(2*Et);
+ double P3CM_max = TMath::Sqrt(E3CM_max*E3CM_max - M3*M3);
+ double E4CM2_max = ((Et-M3)*(Et-M3) + M4*M4 - M5*M5) / (2*(Et-M3));
+ double P4CM2_max = TMath::Sqrt(E4CM2_max*E4CM2_max - M4*M4);
+ double max_weight = P3CM_max * P4CM2_max;
+ double weight = 0.;
+
+ do {
+
+ M = (rnd->RndFsi().Rndm()*(Et - M3 - M4 - M5)) + (M4 + M5);
+ E3CM = (Et*Et + M3*M3 - M*M)/(2*Et);
+ EiCM = Et - E3CM;
+ if(E3CM*E3CM - M3*M3<0)
{
- LOG("INukeUtils",pINFO)
- << "Particle 3 " << M3 << " momentum small or non-finite: " << P3L
- << "\n" << "--> Assigning .001 as new momentum";
- P3tL = 0;
- P3zL = .001;
- P3L = .001;
- E3L = TMath::Sqrt(P3L*P3L + M3*M3);
+ LOG("INukeUtils",pNOTICE)
+ << "PionProduction P3 has non-real momentum - retry kinematics";
+ LOG("INukeUtils",pNOTICE) << "Energy, masses of 3 fs particales:"
+ << E3CM << " " << M3 << " " << " " << M4 << " " << M5;
+ exceptions::INukeException exception;
+ exception.SetReason("PionProduction particle 3 has non-real momentum");
+ throw exception;
+ return false;
}
+ P3CM = TMath::Sqrt(E3CM*E3CM - M3*M3);
+
+ costheta3 = 2.*rnd->RndFsi().Rndm() - 1.;
+ costheta4 = 2.*rnd->RndFsi().Rndm() - 1.;
+ double sintheta3 = TMath::Sqrt( 1. - costheta3*costheta3 );
+ double sintheta4 = TMath::Sqrt( 1. - costheta4*costheta4 );
+ phi3 = 2*kPi * rnd->RndFsi().Rndm();
+ phi4 = 2*kPi * rnd->RndFsi().Rndm();
+
+ P3zL = gm*beta*E3CM + gm*P3CM*costheta3;
+ P3tL = P3CM*sintheta3;
+ PizL = gm*beta*EiCM - gm*P3CM*costheta3;
+ PitL = -P3CM*sintheta3;
+
+ P3L = TMath::Sqrt(P3zL*P3zL + P3tL*P3tL);
+ PiL = TMath::Sqrt(PizL*PizL + PitL*PitL);
+ E3L = TMath::Sqrt(P3L*P3L + M3*M3);
+ EiL = TMath::Sqrt(PiL*PiL + M*M);
+
+ // handle very low momentum particles
+ if(!(TMath::Finite(P3L)) || P3L < .001)
+ {
+ LOG("INukeUtils",pINFO)
+ << "Particle 3 " << M3 << " momentum small or non-finite: " << P3L
+ << "\n" << "--> Assigning .001 as new momentum";
+ P3tL = 0;
+ P3zL = .001;
+ P3L = .001;
+ E3L = TMath::Sqrt(P3L*P3L + M3*M3);
+ }
- tP3L = P3zL*tbetadir + P3tL*tTrans;
- tPiL = PizL*tbetadir + PitL*tTrans;
- tP3L.Rotate(phi3,tbetadir);
- tPiL.Rotate(phi3,tbetadir);
-
- // second sequence, handle formally composite particles 4 and 5
- tbeta2 = tPiL * (1.0 / EiL);
- tbetadir2 = tbeta2.Unit();
- beta2 = tbeta2.Mag();
- gm2 = 1.0 / TMath::Sqrt(1.0 - beta2*beta2);
-
- E4CM2 = (M*M + M4*M4 - M5*M5) / (2*M);
- E5CM2 = M - E4CM2;
- P4CM2 = TMath::Sqrt(E4CM2*E4CM2 - M4*M4);
-
- tVect.SetXYZ(1,0,0);
- if(TMath::Abs((tVect - tbetadir2).Mag())<.01) tVect.SetXYZ(0,1,0);
- theta5 = tVect.Angle(tbetadir2);
- tTrans2 = (tVect - TMath::Cos(theta5)*tbetadir2).Unit();
+ tP3L = P3zL*tbetadir + P3tL*tTrans;
+ tPiL = PizL*tbetadir + PitL*tTrans;
+ tP3L.Rotate(phi3,tbetadir);
+ tPiL.Rotate(phi3,tbetadir);
+
+ // second sequence, handle formally composite particles 4 and 5
+ tbeta2 = tPiL * (1.0 / EiL);
+ tbetadir2 = tbeta2.Unit();
+ beta2 = tbeta2.Mag();
+ gm2 = 1.0 / TMath::Sqrt(1.0 - beta2*beta2);
+
+ E4CM2 = (M*M + M4*M4 - M5*M5) / (2*M);
+ E5CM2 = M - E4CM2;
+ P4CM2 = TMath::Sqrt(E4CM2*E4CM2 - M4*M4);
+
+ tVect.SetXYZ(1,0,0);
+ if(TMath::Abs((tVect - tbetadir2).Mag())<.01) tVect.SetXYZ(0,1,0);
+ theta5 = tVect.Angle(tbetadir2);
+ tTrans2 = (tVect - TMath::Cos(theta5)*tbetadir2).Unit();
+
+ P4zL = gm2*beta2*E4CM2 + gm2*P4CM2*costheta4;
+ P4tL = P4CM2*sintheta4;
+ P5zL = gm2*beta2*E5CM2 - gm2*P4CM2*costheta4;
+ P5tL = - P4tL;
+
+ P4L = TMath::Sqrt(P4zL*P4zL + P4tL*P4tL);
+ P5L = TMath::Sqrt(P5zL*P5zL + P5tL*P5tL);
+ E4L = TMath::Sqrt(P4L*P4L + M4*M4);
+ E5L = TMath::Sqrt(P5L*P5L + M5*M5);
+
+ // handle very low momentum particles
+ if(!(TMath::Finite(P4L)) || P4L < .001)
+ {
+ LOG("INukeUtils",pINFO)
+ << "Particle 4 " << M4 << " momentum small or non-finite: " << P4L
+ << "\n" << "--> Assigning .001 as new momentum";
+ P4tL = 0;
+ P4zL = .001;
+ P4L = .001;
+ E4L = TMath::Sqrt(P4L*P4L + M4*M4);
+ }
+ if(!(TMath::Finite(P5L)) || P5L < .001)
+ {
+ LOG("INukeUtils",pINFO)
+ << "Particle 5 " << M5 << " momentum small or non-finite: " << P5L
+ << "\n" << "--> Assigning .001 as new momentum";
+ P5tL = 0;
+ P5zL = .001;
+ P5L = .001;
+ E5L = TMath::Sqrt(P5L*P5L + M5*M5);
+ }
- P4zL = gm2*beta2*E4CM2 + gm2*P4CM2*TMath::Cos(theta4);
- P4tL = P4CM2*TMath::Sin(theta4);
- P5zL = gm2*beta2*E5CM2 - gm2*P4CM2*TMath::Cos(theta4);
- P5tL = - P4tL;
+ tP4L = P4zL*tbetadir2 + P4tL*tTrans2;
+ tP5L = P5zL*tbetadir2 + P5tL*tTrans2;
+ tP4L.Rotate(phi4,tbetadir2);
+ tP5L.Rotate(phi4,tbetadir2);
- P4L = TMath::Sqrt(P4zL*P4zL + P4tL*P4tL);
- P5L = TMath::Sqrt(P5zL*P5zL + P5tL*P5tL);
- E4L = TMath::Sqrt(P4L*P4L + M4*M4);
- E5L = TMath::Sqrt(P5L*P5L + M5*M5);
+ weight = P3CM*P4CM2;
- // handle very low momentum particles
- if(!(TMath::Finite(P4L)) || P4L < .001)
- {
- LOG("INukeUtils",pINFO)
- << "Particle 4 " << M4 << " momentum small or non-finite: " << P4L
- << "\n" << "--> Assigning .001 as new momentum";
- P4tL = 0;
- P4zL = .001;
- P4L = .001;
- E4L = TMath::Sqrt(P4L*P4L + M4*M4);
- }
- if(!(TMath::Finite(P5L)) || P5L < .001)
- {
- LOG("INukeUtils",pINFO)
- << "Particle 5 " << M5 << " momentum small or non-finite: " << P5L
- << "\n" << "--> Assigning .001 as new momentum";
- P5tL = 0;
- P5zL = .001;
- P5L = .001;
- E5L = TMath::Sqrt(P5L*P5L + M5*M5);
- }
+ if ( bias != 0 ) weight *= TMath::Exp( bias*(TLorentzVector(tP3L,E3L) - *p->P4()).M2() );
- tP4L = P4zL*tbetadir2 + P4tL*tTrans2;
- tP5L = P5zL*tbetadir2 + P5tL*tTrans2;
- tP4L.Rotate(phi4,tbetadir2);
- tP5L.Rotate(phi4,tbetadir2);
+ } while ( rnd->RndFsi().Rndm() > weight/max_weight );
// pauli blocking
if(P3L < FermiMomentum || ( pdg::IsNeutronOrProton(s2->Pdg()) && P4L < FermiMomentum ) )
@@ -1285,7 +1314,7 @@ bool genie::utils::intranuke2025::ThreeBodyKinematics(
//___________________________________________________________________________
bool genie::utils::intranuke2025::PionProduction(
GHepRecord* ev, GHepParticle* p, GHepParticle* s1, GHepParticle* s2, GHepParticle* s3, int &RemnA, int &RemnZ,
- TLorentzVector &RemnP4, bool DoFermi, double FermiFac, double FermiMomentum, const NuclearModelI* Nuclmodel)
+ TLorentzVector &RemnP4, bool DoFermi, double FermiFac, double FermiMomentum, const NuclearModelI* Nuclmodel, double bias)
{
// Aaron Meyer (7/15/2010)
//
@@ -1605,7 +1634,7 @@ bool genie::utils::intranuke2025::PionProduction(
s3->SetPdgCode(p5code);
if(genie::utils::intranuke2025::ThreeBodyKinematics(
- ev,p,(ptarg?kPdgProton:kPdgNeutron),s1,s2,s3,DoFermi,FermiFac,FermiMomentum,Nuclmodel))
+ ev,p,(ptarg?kPdgProton:kPdgNeutron),s1,s2,s3,DoFermi,FermiFac,FermiMomentum,Nuclmodel,bias))
{
// okay, handle remnants and return true
// assumes first particle is always the nucleon,
diff --git a/src/Physics/HadronTransport/INukeUtils2025.h b/src/Physics/HadronTransport/INukeUtils2025.h
index 514b58d0b..7eb303387 100644
--- a/src/Physics/HadronTransport/INukeUtils2025.h
+++ b/src/Physics/HadronTransport/INukeUtils2025.h
@@ -89,11 +89,11 @@ namespace intranuke2025
bool ThreeBodyKinematics(
GHepRecord* ev, GHepParticle* p, int tcode, GHepParticle* s1, GHepParticle* s2, GHepParticle* s3,
- bool DoFermi=false, double FermiFac=0, double FermiMomentum=0, const NuclearModelI* Nuclmodel=(const NuclearModelI*)0);
+ bool DoFermi=false, double FermiFac=0, double FermiMomentum=0, const NuclearModelI* Nuclmodel=(const NuclearModelI*)0, double bias=0);
bool PionProduction(
GHepRecord* ev, GHepParticle* p, GHepParticle* s1, GHepParticle* s2, GHepParticle* s3, int &RemnA, int &RemnZ,
- TLorentzVector &RemnP4,bool DoFermi, double FermiFac, double FermiMomentum, const NuclearModelI* Nuclmodel);
+ TLorentzVector &RemnP4,bool DoFermi, double FermiFac, double FermiMomentum, const NuclearModelI* Nuclmodel, double bias);
double CalculateEta(
double Minc, double ke, double Mtarg, double Mtwopart, double Mpi);
diff --git a/src/Physics/HadronTransport/Intranuke2018.h b/src/Physics/HadronTransport/Intranuke2018.h
index ad067cb55..f2d4bede0 100644
--- a/src/Physics/HadronTransport/Intranuke2018.h
+++ b/src/Physics/HadronTransport/Intranuke2018.h
@@ -153,6 +153,9 @@ public :
double fNucleonFracAbsScale;
double fNucleonFracPiProdScale;
+ /// Bias parameter used to adjust pion production distribution away from pure
+ /// three-body phase space
+ double fPiProdThreeBodyBias;
};
} // genie namespace
diff --git a/src/Physics/HadronTransport/Intranuke2025.h b/src/Physics/HadronTransport/Intranuke2025.h
index 65a79d152..a683b21bf 100644
--- a/src/Physics/HadronTransport/Intranuke2025.h
+++ b/src/Physics/HadronTransport/Intranuke2025.h
@@ -155,6 +155,10 @@ public :
double fNucleonFracAbsScale;
double fNucleonFracPiProdScale;
+ /// Bias parameter used to adjust pion production distribution away from pure
+ /// three-body phase space
+ double fPiProdThreeBodyBias;
+
};
} // genie namespace