diff --git a/config/HAIntranuke2018.xml b/config/HAIntranuke2018.xml
index 9db9957e63..7e47d22d16 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 a817050da2..4ecec030f4 100644
--- a/src/Physics/HadronTransport/HAIntranuke2018.cxx
+++ b/src/Physics/HadronTransport/HAIntranuke2018.cxx
@@ -2,7 +2,7 @@
/*
Copyright (c) 2003-2024, 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/HNIntranuke2018.cxx b/src/Physics/HadronTransport/HNIntranuke2018.cxx
index 44074e420f..910f7b8b3b 100644
--- a/src/Physics/HadronTransport/HNIntranuke2018.cxx
+++ b/src/Physics/HadronTransport/HNIntranuke2018.cxx
@@ -3,7 +3,7 @@
/*
Copyright (c) 2003-2024, 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/INukeUtils2018.cxx b/src/Physics/HadronTransport/INukeUtils2018.cxx
index a235ff7cf3..28e212bff3 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 aac37e9a29..c16e00d716 100644
--- a/src/Physics/HadronTransport/INukeUtils2018.h
+++ b/src/Physics/HadronTransport/INukeUtils2018.h
@@ -85,11 +85,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/Intranuke2018.h b/src/Physics/HadronTransport/Intranuke2018.h
index 7eaeab557c..55c399fd14 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