From 7ac0caa05470143f6433daad57d18b73abacd6c7 Mon Sep 17 00:00:00 2001 From: Gray Putnam Date: Thu, 20 Feb 2025 16:47:57 -0600 Subject: [PATCH 1/3] Properly sample lorentz-invariant-phase-space for pion production kinematics. Add in optional bias towards forward nucleon momentum. --- .../HadronTransport/INukeUtils2018.cxx | 213 ++++++++++-------- src/Physics/HadronTransport/INukeUtils2018.h | 2 +- 2 files changed, 121 insertions(+), 94 deletions(-) diff --git a/src/Physics/HadronTransport/INukeUtils2018.cxx b/src/Physics/HadronTransport/INukeUtils2018.cxx index 773e1c7224..fe116bce01 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 ) ) diff --git a/src/Physics/HadronTransport/INukeUtils2018.h b/src/Physics/HadronTransport/INukeUtils2018.h index dbbb00c5eb..874080bb88 100644 --- a/src/Physics/HadronTransport/INukeUtils2018.h +++ b/src/Physics/HadronTransport/INukeUtils2018.h @@ -89,7 +89,7 @@ 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=6); bool PionProduction( GHepRecord* ev, GHepParticle* p, GHepParticle* s1, GHepParticle* s2, GHepParticle* s3, int &RemnA, int &RemnZ, From 2a8e515ded2b1e44dc03a6c4403f7896011541a4 Mon Sep 17 00:00:00 2001 From: Steven Gardiner Date: Mon, 15 Sep 2025 17:41:04 +0100 Subject: [PATCH 2/3] Adjust biasing of the three-body phase space used in hA/hN to allow XML-based configuration of the bias parameter. --- config/HAIntranuke2018.xml | 4 +- .../HadronTransport/HAIntranuke2018.cxx | 7 +- .../HadronTransport/HNIntranuke2018.cxx | 104 +++++++++--------- .../HadronTransport/INukeUtils2018.cxx | 9 +- src/Physics/HadronTransport/INukeUtils2018.h | 5 +- src/Physics/HadronTransport/Intranuke2018.h | 3 + 6 files changed, 72 insertions(+), 60 deletions(-) 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 682470dde1..108bf7b819 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/HNIntranuke2018.cxx b/src/Physics/HadronTransport/HNIntranuke2018.cxx index 3a1b6d796a..2a6657f196 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/INukeUtils2018.cxx b/src/Physics/HadronTransport/INukeUtils2018.cxx index fe116bce01..cb4b40535d 100644 --- a/src/Physics/HadronTransport/INukeUtils2018.cxx +++ b/src/Physics/HadronTransport/INukeUtils2018.cxx @@ -1158,9 +1158,9 @@ bool genie::utils::intranuke2018::ThreeBodyKinematics( // G.P. 2/20/2025 // Sample uniformly in lorentz invariany phase space // - // Just because you sample uniformly does not mean you are + // 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 @@ -1327,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) // @@ -1647,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 874080bb88..00399f062c 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, double bias=6); + 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 ad067cb558..f2d4bede01 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 From dfd8a7b0e0eae2c6e476034afa08d1250a90d6b8 Mon Sep 17 00:00:00 2001 From: Steven Gardiner Date: Tue, 2 Jun 2026 11:12:17 -0500 Subject: [PATCH 3/3] Apply Gray's fixes for three-body phase space sampling in pion production FSI to the hA/hN 2025 models. The adjustments are intended to match what was previously done for hA/hN 2018. --- .../HadronTransport/HAIntranuke2025.cxx | 29 +-- .../HadronTransport/HNIntranuke2025.cxx | 101 ++++---- src/Physics/HadronTransport/HNIntranuke2025.h | 8 +- .../HadronTransport/INukeUtils2025.cxx | 221 ++++++++++-------- src/Physics/HadronTransport/INukeUtils2025.h | 4 +- src/Physics/HadronTransport/Intranuke2025.h | 4 + 6 files changed, 203 insertions(+), 164 deletions(-) diff --git a/src/Physics/HadronTransport/HAIntranuke2025.cxx b/src/Physics/HadronTransport/HAIntranuke2025.cxx index cb5027f69f..7462116455 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/HNIntranuke2025.cxx b/src/Physics/HadronTransport/HNIntranuke2025.cxx index b00ffbcc67..81640fe16d 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 c629b90746..928c948b29 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/INukeUtils2025.cxx b/src/Physics/HadronTransport/INukeUtils2025.cxx index b1951bbbf8..f927ff7e89 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 514b58d0bf..7eb3033870 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/Intranuke2025.h b/src/Physics/HadronTransport/Intranuke2025.h index 65a79d1522..a683b21bf0 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