From c28d8246e99842f23c4a6c902d1a8a05228d2ee5 Mon Sep 17 00:00:00 2001 From: Gray Putnam Date: Fri, 4 Apr 2025 17:38:04 -0500 Subject: [PATCH 1/7] First draft of adding N-N kinematics reweighting parameters. --- .../GReWeightINukeKinematics.cxx | 164 ++++++++++++++++ src/RwCalculators/GReWeightINukeKinematics.h | 68 +++++++ .../GReWeightINukeKinematicsParams.cxx | 180 ++++++++++++++++++ .../GReWeightINukeKinematicsParams.h | 48 +++++ 4 files changed, 460 insertions(+) create mode 100644 src/RwCalculators/GReWeightINukeKinematics.cxx create mode 100644 src/RwCalculators/GReWeightINukeKinematics.h create mode 100644 src/RwCalculators/GReWeightINukeKinematicsParams.cxx create mode 100644 src/RwCalculators/GReWeightINukeKinematicsParams.h diff --git a/src/RwCalculators/GReWeightINukeKinematics.cxx b/src/RwCalculators/GReWeightINukeKinematics.cxx new file mode 100644 index 00000000..6eed93d4 --- /dev/null +++ b/src/RwCalculators/GReWeightINukeKinematics.cxx @@ -0,0 +1,164 @@ +#include +#include + +// GENIE/Generator includes +#include "Framework/Algorithm/AlgFactory.h" +#include "Framework/Algorithm/AlgConfigPool.h" +#include "Framework/Conventions/Units.h" +#include "Framework/EventGen/EventRecord.h" +#include "Framework/GHEP/GHepParticle.h" +#include "Framework/Messenger/Messenger.h" +#include "Framework/Numerical/Spline.h" +#include "Framework/ParticleData/PDGUtils.h" +#include "Framework/Registry/Registry.h" +#include "Physics/NuclearState/NuclearUtils.h" + +// GENIE/Reweight includes +#include "RwCalculators/GReWeightINukeKinematics.h" +#include "RwCalculators/GReWeightUtils.h" +#include "RwFramework/GSystUncertainty.h" + +using namespace genie; +using namespace genie::rew; + +//_______________________________________________________________________________________ +bool GReWeightINukeKinematics::AppliesTo(const EventRecord & event) const +{ + auto type = event.Summary()->ProcInfo().ScatteringTypeId(); + switch (type) { + case kScCoherentProduction: + case kScDiffractive: + case kScNuElectronElastic: + case kScAMNuGamma: + case kScCoherentElastic: + return false; + default: + return true; + } +} + + +//_______________________________________________________________________________________ +bool GReWeightINukeKinematics::IsHandled(GSyst_t syst) const +{ + switch (syst) { + case ( kINukeKinematicsTwkDial_NP_N ): + default: + case ( kINukeKinematicsTwkDial_PP_N ): + return true; + } + + return false; + +} + +//_______________________________________________________________________________________ +void GReWeightINukeKinematics::SetSystematic(GSyst_t syst, double val) +{ + if (IsHandled(syst)) { + fRwParams.SetTwkDial(syst, val); + } +} + +//_______________________________________________________________________________________ +void GReWeightINukeKinematics::Reset() { + fRwParams.Reset(); +} + +//_______________________________________________________________________________________ +void GReWeightINukeKinematics::Reconfigure() {} + +// "Lab" energy of particle p in p + t -> f1 + f2 scattering process. +// That is -- the kinetic energy of p in the frame where t is at rest +double TLab(const TLorentzVector &p, const TLorentzVector &f1, const TLorentzVector &f2) { + TLorentzVector t = f1 + f2 - p; // E.M. conservation + double E_p = ((p + t).Mag2() - t.M2() - p.M2())/(2.0*t.M()); + + return E_p - p.M(); +} + +// Scattering angle in p + t -> f1 + f2 scattering process, +// in center-of-momentum (C.O.M.) frame +double costhcm(const TLorentzVector &p, const TLorentzVector &f1, const TLorentzVector &f2) { + TLorentzVector t = f1 + f2 - p; // E.M. conservation + + double Ecm2 = (f1 + f2).Mag2(); + + double E_p_cm = (Ecm2 + p.M2() - t.M2()) / (2*sqrt(Ecm2)); + double P_p_cm = sqrt(E_p_cm*E_p_cm - p.M2()); + double E_f1_cm = (Ecm2 + f1.M2() - f2.M2()) / (2*sqrt(Ecm2)); + double P_f1_cm = sqrt(E_f1_cm*E_f1_cm - f1.M2()); + + return ((p + f1).M2() - p.M2() - f1.M2() - 2*E_p_cm*E_f1_cm) / (2*P_f1_cm*P_p_cm); +} + +double GReWeightINukeKinematics::CalcWeight(const EventRecord &event) { + // to return + double event_weight = 1.; + + // get the atomic mass number for the hit nucleus + GHepParticle * tgt = event.TargetNucleus(); + if (!tgt) return 1.0; + double A = tgt->A(); + double Z = tgt->Z(); + if (A<=1) return 1.0; + if (Z<=1) return 1.0; + + fRwParams.SetTargetA( A ); + + // Loop over stdhep entries and only calculate weights for particles. + // All particles that are not hadrons generated inside the nucleus are given weights of 1.0 + int ip=-1; + GHepParticle * p = 0; + TIter event_iter(&event); + while ( (p = dynamic_cast(event_iter.Next())) ) { + ip++; + + // Skip particles not rescattered by the actual hadron transport code + int pdgc = p->Pdg(); + // bool is_pion = pdg::IsPion (pdgc); + bool is_nucleon = pdg::IsNucleon(pdgc); + // bool is_kaon = pdg::IsKaon( pdgc ); + if (!is_nucleon) + { + continue; + } + + // Skip particles with code other than 'hadron in the nucleus' + GHepStatus_t ist = p->Status(); + if(ist != kIStHadronInTheNucleus) + { + continue; + } + + // Determine the interaction type for current hadron in nucleus, if any + int fsi_code = p->RescatterCode(); + LOG("RWINukeKin", pDEBUG) + << "Attempting to reweight hadron at position = " << ip + << " with PDG code = " << pdgc + << " and FSI code = " << fsi_code + << " (" << INukeHadroFates::AsString((INukeFateHA_t)fsi_code) << ")\n"; + + // Inelastic, Charge Exchange + if (fsi_code == (int)kIHAFtInelas || fsi_code == (int)kIHAFtCEx) { + const GHepParticle &f1 = *event.Particle(p->FirstDaughter()); + const GHepParticle &f2 = *event.Particle(p->LastDaughter()); + + int np_pp = 0; // == 0 if nucleons in scatter are different, == 1 if they are the same + if (fsi_code == (int)kIHAFtInelas) np_pp = (f1.Pdg() != f2.Pdg()) ? 0 : 1; + else if (fsi_code == (int)kIHAFtCEx) np_pp = (f1.Pdg() == f2.Pdg()) ? 0 : 1; + + double tlab_v = TLab(*p->P4(), *f1.P4(), *f2.P4()); + double costhcm_v = costhcm(*p->P4(), *f1.P4(), *f2.P4()); + + double weight = fRwParams.CalcWeight(np_pp, tlab_v, costhcm_v); + + // std::cerr << "TLab: " << tlab_v << " costhcm: " << costhcm_v << " np/pp " << np_pp << " weight: " << weight << std::endl; + + event_weight *= weight; + } + } + + return event_weight; +} + diff --git a/src/RwCalculators/GReWeightINukeKinematics.h b/src/RwCalculators/GReWeightINukeKinematics.h new file mode 100644 index 00000000..ddf1cfd5 --- /dev/null +++ b/src/RwCalculators/GReWeightINukeKinematics.h @@ -0,0 +1,68 @@ +//____________________________________________________________________________ +/*! + +\class genie::rew::GReWeightINukeKinematics + +\brief Reweighting GENIE INTRANUKE/hA hadron transport model. + +\author Gray Putnam + Fermi National Accelerator Laboratory + +\created Feb 14, 2025 + +\cpright Copyright (c) 2003-2024, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org +*/ +//____________________________________________________________________________ + +#ifndef _G_REWEIGHT_INUKEKINEMATICS_H_ +#define _G_REWEIGHT_INUKEKINEMATICS_H_ + +//#define _G_REWEIGHT_INUKEKINEMATICS_DEBUG_NTP_ + +// GENIE/Reweight includes +#include "RwCalculators/GReWeightModel.h" +#include "RwCalculators/GReWeightINukeKinematicsParams.h" + +using namespace genie::rew; +using namespace genie; + +class TFile; +class TNtuple; +class TLorentzVector; + +namespace genie { + + class HAIntranuke2018; + class GHepParticle; + +namespace rew { + + class GReWeightINukeKinematics : public GReWeightModel + { + public: + GReWeightINukeKinematics(): GReWeightModel("IntraNukeKin") {} + ~GReWeightINukeKinematics() {} + + // implement the GReWeightI interface + bool AppliesTo (const EventRecord & event) const; + bool IsHandled (GSyst_t syst) const; + void SetSystematic (GSyst_t syst, double val); + void Reset (void); + void Reconfigure (void); + double CalcWeight (const EventRecord & event); + + private: + GReWeightINukeKinematicsParams fRwParams; + +#ifdef _G_REWEIGHT_INUKEKINEMATICS_DEBUG_NTP_ + TFile * fTestFile; + TNtuple * fTestNtp; +#endif + + }; + +} // rew +} // genie + +#endif diff --git a/src/RwCalculators/GReWeightINukeKinematicsParams.cxx b/src/RwCalculators/GReWeightINukeKinematicsParams.cxx new file mode 100644 index 00000000..180af3cb --- /dev/null +++ b/src/RwCalculators/GReWeightINukeKinematicsParams.cxx @@ -0,0 +1,180 @@ +#include +#include + +#include + +// GENIE/Generator includes +#include "Framework/Conventions/Controls.h" +#include "Framework/Conventions/Constants.h" +#include "Framework/Conventions/Units.h" +#include "Framework/Messenger/Messenger.h" +#include "Framework/Messenger/Messenger.h" +#include "Framework/Numerical/Spline.h" + +// GENIE/Reweight includes +#include "RwCalculators/GReWeightINukeKinematicsParams.h" +#include "RwCalculators/GReWeightUtils.h" +#include "RwFramework/GSystUncertainty.h" + +using namespace genie; +using namespace genie::rew; +using namespace genie::constants; + +//___________________________________________________________________________ +GReWeightINukeKinematicsParams::GReWeightINukeKinematicsParams(): + fNPwgt("pn"), + fPPwgt("pp") {} + +//___________________________________________________________________________ +void GReWeightINukeKinematicsParams::SetTargetA(int A) {(void) A;} // nothing, for now + +double GReWeightINukeKinematicsParams::CalcWeight(int np_pp, float tlab, float costhcm) { + if (np_pp == 0) + return fNPwgt.CalcWeight(tlab, costhcm); + + return fPPwgt.CalcWeight(tlab, costhcm); +} + +//___________________________________________________________________________ +void GReWeightINukeKinematicsParams::SetTwkDial(GSyst_t syst, double val) +{ + if (syst == kINukeKinematicsTwkDial_NP_N) { + fNPwgt.SetUniverse((int)val); + } + else if (syst == kINukeKinematicsTwkDial_PP_N) { + fPPwgt.SetUniverse((int)val); + } +} + +//___________________________________________________________________________ +void GReWeightINukeKinematicsParams::Reset() { + fNPwgt.Reset(); + fPPwgt.Reset(); +} + +void GReWeightINukeKinematicsParams::ReBounce::Reset() { + fRewt.reset(); + fCache.clear(); +} + +//___________________________________________________________________________ +void GReWeightINukeKinematicsParams::ReBounce::SetUniverse(int u) { + fUniverse = u; + + if (u == 0) { // set to CV + fRewt.reset(); + return; + } + + LoadUniverse(u); + fRewt = fCache[u]; +} + +double GReWeightINukeKinematicsParams::ReBounce::CalcWeight(float tlab, float costhcm) { + // std::cerr << "UNIVERSE: " << fUniverse << " VALID: " << ((int)!!fRewt) << " dw: " << ((!fRewt) ? 0. : fRewt->Evaluate(tlab, costhcm)) << std::endl; + if (!fRewt) return 1; + + return 1 + fRewt->Evaluate(tlab, costhcm); +} + +// Helper function for reading hN cross section file +//____________________________________________________________________________ +void ReadhNFile( + string filename, double ke, int npoints, int & curr_point, + double * costh_array, double * wgt_array, int cols) +{ + // open + std::ifstream hN_stream(filename.c_str(), std::ios::in); + if(!hN_stream.good()) { + LOG("RWINukeKin", pERROR) + << "Error reading INTRANUKE/hN data from: " << filename; + return; + } + if(cols<2) { + LOG("RWINukeKin", pERROR) + << "Error reading INTRANUKE/hN data from: " << filename; + LOG("RWINukeKin", pERROR) + << "Too few columns: " << cols; + return; + } + + LOG("RWINukeKin", pINFO) + << "Reading INTRANUKE/hN data from: " << filename; + + // skip initial comments + char cbuf[501]; + hN_stream.getline(cbuf,400); + hN_stream.getline(cbuf,400); + hN_stream.getline(cbuf,400); + + // read + double angle = 0; + double wgt = 0; + double trash = 0; + + for(int ip = 0; ip < npoints; ip++) { + hN_stream >> angle >> wgt; + + for(int ic = 0; ic < (cols-2); ic++) { + hN_stream >> trash; + } + + LOG("RWINukeKin", pDEBUG) + << "Adding data point: (KE = " << ke << " MeV, angle = " + << angle << ", sigma = " << wgt << ")"; + costh_array[ip] = TMath::Cos(angle*kPi/180.); + wgt_array [curr_point] = wgt; + curr_point++; + } +} + +std::shared_ptr LoadWgts(std::string fname) { + string data_dir = (gSystem->Getenv("GINUKEHADRONDATA")) ? + string(gSystem->Getenv("GINUKEHADRONDATA")) : + string(gSystem->Getenv("GENIE")) + string("/data/syst/diff_ang_variations/"); + + // Hard-coded parameters for hadron data, taken from + // Physics/HadronTransport/INukeHadroData2018.cxx + const int hN_wgt_nfiles = 20; + const int hN_wgt_points_per_file = 21; + const int hN_wgt_npoints = hN_wgt_points_per_file * hN_wgt_nfiles; + + double hN_wgt_energies[hN_wgt_nfiles] = { + 50, 100, 150, 200, 250, 300, 350, 400, 450, 500, + 550, 600, 650, 700, 750, 800, 850, 900, 950, 1000 + }; + + double hN_wgt_costh [hN_wgt_points_per_file] = {}; + double hN_wgt_wgt [hN_wgt_npoints] = {}; + + int ipoint=0; + + for(int ifile = 0; ifile < hN_wgt_nfiles; ifile++) { + // build filename + std::ostringstream hN_datafile; + double ke = hN_wgt_energies[ifile]; + hN_datafile << data_dir << "/" << fname << ke << ".txt"; + // read data + ReadhNFile( + hN_datafile.str(), ke, hN_wgt_points_per_file, + ipoint, hN_wgt_costh, hN_wgt_wgt,2); + }//loop over files + + return std::shared_ptr(new BLI2DNonUnifGrid(hN_wgt_nfiles,hN_wgt_points_per_file, + hN_wgt_energies,hN_wgt_costh,hN_wgt_wgt)); + +} + +GReWeightINukeKinematicsParams::ReBounce::~ReBounce() {} + +GReWeightINukeKinematicsParams::ReBounce::ReBounce(std::string n): + fUniverse(0), fName(n) {} + +void GReWeightINukeKinematicsParams::ReBounce::LoadUniverse(int u) { + // Do nothing if we already have the wgt + if (fCache.count(u) > 0) return; + + std::string file_name = "diff_ang_univ" + std::to_string(u) + "/" + fName + "/" + fName; + fCache[u] = LoadWgts(file_name); +} + diff --git a/src/RwCalculators/GReWeightINukeKinematicsParams.h b/src/RwCalculators/GReWeightINukeKinematicsParams.h new file mode 100644 index 00000000..2c634935 --- /dev/null +++ b/src/RwCalculators/GReWeightINukeKinematicsParams.h @@ -0,0 +1,48 @@ +#ifndef _G_REWEIGHT_INTRANUKEKINEMATICS_PARAMS_H_ +#define _G_REWEIGHT_INTRANUKEKINEMATICS_PARAMS_H_ + +#include + +// GENIE/Reweight includes +#include "RwFramework/GSyst.h" +#include "Framework/Numerical/BLI2D.h" + +namespace genie { +namespace rew { + + class GReWeightINukeKinematicsParams { + + public: + GReWeightINukeKinematicsParams(); + ~GReWeightINukeKinematicsParams() {} + + void Reset (void); ///< + void SetTwkDial (GSyst_t s, double val); ///< + void SetTargetA (int target_A); ///< Set the mass number of the hit nucleus + double CalcWeight (int np_pp, float tlab, float costhcm); + + struct ReBounce { + void SetUniverse(int u); + void LoadUniverse(int u); + void Reset(); + double CalcWeight(float tlab, float costhcm); + + ReBounce(std::string n); + ~ReBounce(); + + int fUniverse; + std::string fName; + std::shared_ptr fRewt; + std::map> fCache; + }; + + private: + ReBounce fNPwgt; + ReBounce fPPwgt; + + }; //GReWeightINukeKinematicsParams + +} // rew namespace +} // genie namespace + +#endif // _G_REWEIGHT_INTRANUKEKINEMATICS_PARAMS_H_ From 9eea238d1808ce864d2bb2be3a78ef6901005c57 Mon Sep 17 00:00:00 2001 From: Gray Putnam Date: Fri, 4 Apr 2025 17:39:44 -0500 Subject: [PATCH 2/7] Integrate N-N kinematics parameters into framework. --- src/Apps/gRwght1Param.cxx | 2 ++ src/RwFramework/GSyst.cxx | 2 ++ src/RwFramework/GSyst.h | 8 ++++++-- 3 files changed, 10 insertions(+), 2 deletions(-) diff --git a/src/Apps/gRwght1Param.cxx b/src/Apps/gRwght1Param.cxx index f1cd77bb..8836faa9 100644 --- a/src/Apps/gRwght1Param.cxx +++ b/src/Apps/gRwght1Param.cxx @@ -118,6 +118,7 @@ #include "RwCalculators/GReWeightResonanceDecay.h" #include "RwCalculators/GReWeightFZone.h" #include "RwCalculators/GReWeightINuke.h" +#include "RwCalculators/GReWeightINukeKinematics.h" #include "RwCalculators/GReWeightAGKY.h" #include "RwCalculators/GReWeightNuXSecCCQEaxial.h" #include "RwCalculators/GReWeightNuXSecCCQEvec.h" @@ -247,6 +248,7 @@ int main(int argc, char ** argv) rw.AdoptWghtCalc( "hadro_res_decay", new GReWeightResonanceDecay ); rw.AdoptWghtCalc( "hadro_fzone", new GReWeightFZone ); rw.AdoptWghtCalc( "hadro_intranuke", new GReWeightINuke ); + rw.AdoptWghtCalc( "hadro_inuke_kin", new GReWeightINukeKinematics ); rw.AdoptWghtCalc( "hadro_agky", new GReWeightAGKY ); // GReWeightDISNuclMod::CalcWeight() not implemented - don't try to use it .. diff --git a/src/RwFramework/GSyst.cxx b/src/RwFramework/GSyst.cxx index c5c3030d..fba65ede 100644 --- a/src/RwFramework/GSyst.cxx +++ b/src/RwFramework/GSyst.cxx @@ -95,6 +95,8 @@ std::map GSyst::BuildGSystToStringMap() { temp_map[ kINukeTwkDial_FrInel_N ] = "FrInel_N"; temp_map[ kINukeTwkDial_FrAbs_N ] = "FrAbs_N"; temp_map[ kINukeTwkDial_FrPiProd_N ] = "FrPiProd_N"; + temp_map[ kINukeKinematicsTwkDial_NP_N ] = "FrKin_NP_N"; + temp_map[ kINukeKinematicsTwkDial_PP_N ] = "FrKin_PP_N"; temp_map[ kSystNucl_CCQEPauliSupViaKF ] = "CCQEPauliSupViaKF"; temp_map[ kSystNucl_CCQEMomDistroFGtoSF ] = "CCQEMomDistroFGtoSF"; temp_map[ kRDcyTwkDial_BR1gamma ] = "RDecBR1gamma"; diff --git a/src/RwFramework/GSyst.h b/src/RwFramework/GSyst.h index d911a677..a93bf389 100644 --- a/src/RwFramework/GSyst.h +++ b/src/RwFramework/GSyst.h @@ -113,11 +113,12 @@ typedef enum EGSyst { // // Intranuclear rescattering systematics. - // There are 2 sets of parameters: + // There are 3 sets of parameters: // - parameters that control the total rescattering probability, P(total) // - parameters that control the fraction of each process (`fate'), given a total rescat. prob., P(fate|total) // These parameters are considered separately for pions and nucleons. - // + // - parameters that control the kinematics of the scattering process + // kINukeTwkDial_MFP_pi, ///< tweak mean free path for pions kINukeTwkDial_MFP_N, ///< tweak mean free path for nucleons @@ -133,6 +134,9 @@ typedef enum EGSyst { kINukeTwkDial_FrAbs_N, ///< tweak absorption probability for nucleons, for given total rescattering probability kINukeTwkDial_FrPiProd_N, ///< tweak pion production probability for nucleons, for given total rescattering probability + kINukeKinematicsTwkDial_NP_N, ///< tweak scattering angle of NP scatters for nucleons + kINukeKinematicsTwkDial_PP_N, ///< tweak scattering angle of PP and NN scatters for nucleons + // // Nuclear model // From 9d9f9b55fbd111d5fea26e6680db4d09093d71f7 Mon Sep 17 00:00:00 2001 From: Gray Putnam Date: Fri, 4 Apr 2025 17:43:01 -0500 Subject: [PATCH 3/7] bugfix --- src/RwCalculators/GReWeightINukeKinematics.cxx | 1 - 1 file changed, 1 deletion(-) diff --git a/src/RwCalculators/GReWeightINukeKinematics.cxx b/src/RwCalculators/GReWeightINukeKinematics.cxx index 6eed93d4..abcb51b2 100644 --- a/src/RwCalculators/GReWeightINukeKinematics.cxx +++ b/src/RwCalculators/GReWeightINukeKinematics.cxx @@ -43,7 +43,6 @@ bool GReWeightINukeKinematics::IsHandled(GSyst_t syst) const { switch (syst) { case ( kINukeKinematicsTwkDial_NP_N ): - default: case ( kINukeKinematicsTwkDial_PP_N ): return true; } From 01a41c01da7cfb6f309b3a5e8ecc926c5caabd4a Mon Sep 17 00:00:00 2001 From: Gray Putnam Date: Fri, 4 Apr 2025 17:43:25 -0500 Subject: [PATCH 4/7] add default to switch. --- src/RwCalculators/GReWeightINukeKinematics.cxx | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/RwCalculators/GReWeightINukeKinematics.cxx b/src/RwCalculators/GReWeightINukeKinematics.cxx index abcb51b2..e3491d15 100644 --- a/src/RwCalculators/GReWeightINukeKinematics.cxx +++ b/src/RwCalculators/GReWeightINukeKinematics.cxx @@ -45,6 +45,8 @@ bool GReWeightINukeKinematics::IsHandled(GSyst_t syst) const case ( kINukeKinematicsTwkDial_NP_N ): case ( kINukeKinematicsTwkDial_PP_N ): return true; + default: + return false; } return false; From 21923aeba2d5a91fa94329a90c95516cbc972f1b Mon Sep 17 00:00:00 2001 From: Gray Putnam Date: Wed, 11 Jun 2025 11:08:33 -0500 Subject: [PATCH 5/7] Update logging to be consistent with other reweight code. --- src/RwCalculators/GReWeightINukeKinematics.cxx | 2 +- src/RwCalculators/GReWeightINukeKinematicsParams.cxx | 12 ++++++------ 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/RwCalculators/GReWeightINukeKinematics.cxx b/src/RwCalculators/GReWeightINukeKinematics.cxx index e3491d15..c339c092 100644 --- a/src/RwCalculators/GReWeightINukeKinematics.cxx +++ b/src/RwCalculators/GReWeightINukeKinematics.cxx @@ -134,7 +134,7 @@ double GReWeightINukeKinematics::CalcWeight(const EventRecord &event) { // Determine the interaction type for current hadron in nucleus, if any int fsi_code = p->RescatterCode(); - LOG("RWINukeKin", pDEBUG) + LOG("ReW", pDEBUG) << "Attempting to reweight hadron at position = " << ip << " with PDG code = " << pdgc << " and FSI code = " << fsi_code diff --git a/src/RwCalculators/GReWeightINukeKinematicsParams.cxx b/src/RwCalculators/GReWeightINukeKinematicsParams.cxx index 180af3cb..49ae978a 100644 --- a/src/RwCalculators/GReWeightINukeKinematicsParams.cxx +++ b/src/RwCalculators/GReWeightINukeKinematicsParams.cxx @@ -71,7 +71,7 @@ void GReWeightINukeKinematicsParams::ReBounce::SetUniverse(int u) { } double GReWeightINukeKinematicsParams::ReBounce::CalcWeight(float tlab, float costhcm) { - // std::cerr << "UNIVERSE: " << fUniverse << " VALID: " << ((int)!!fRewt) << " dw: " << ((!fRewt) ? 0. : fRewt->Evaluate(tlab, costhcm)) << std::endl; + // std::cout << "UNIVERSE: " << fUniverse << " VALID: " << ((int)!!fRewt) << " dw: " << ((!fRewt) ? 0. : fRewt->Evaluate(tlab, costhcm)) << std::endl; if (!fRewt) return 1; return 1 + fRewt->Evaluate(tlab, costhcm); @@ -86,19 +86,19 @@ void ReadhNFile( // open std::ifstream hN_stream(filename.c_str(), std::ios::in); if(!hN_stream.good()) { - LOG("RWINukeKin", pERROR) + LOG("ReW", pERROR) << "Error reading INTRANUKE/hN data from: " << filename; return; } if(cols<2) { - LOG("RWINukeKin", pERROR) + LOG("ReW", pERROR) << "Error reading INTRANUKE/hN data from: " << filename; - LOG("RWINukeKin", pERROR) + LOG("ReW", pERROR) << "Too few columns: " << cols; return; } - LOG("RWINukeKin", pINFO) + LOG("ReW", pINFO) << "Reading INTRANUKE/hN data from: " << filename; // skip initial comments @@ -119,7 +119,7 @@ void ReadhNFile( hN_stream >> trash; } - LOG("RWINukeKin", pDEBUG) + LOG("ReW", pDEBUG) << "Adding data point: (KE = " << ke << " MeV, angle = " << angle << ", sigma = " << wgt << ")"; costh_array[ip] = TMath::Cos(angle*kPi/180.); From 65a1c9e7ed6b281df10cee0407a479dd46734f3b Mon Sep 17 00:00:00 2001 From: Gray Putnam Date: Wed, 11 Jun 2025 11:09:48 -0500 Subject: [PATCH 6/7] Add in pion production FSI fix, and kinematic biasing, to weight calculators. --- .../GReWeightINukeKinematics.cxx | 13 ++ .../GReWeightINukeKinematicsParams.cxx | 212 +++++++++++++++++- .../GReWeightINukeKinematicsParams.h | 5 + src/RwFramework/GSyst.cxx | 3 + src/RwFramework/GSyst.h | 3 + 5 files changed, 235 insertions(+), 1 deletion(-) diff --git a/src/RwCalculators/GReWeightINukeKinematics.cxx b/src/RwCalculators/GReWeightINukeKinematics.cxx index c339c092..ea94c357 100644 --- a/src/RwCalculators/GReWeightINukeKinematics.cxx +++ b/src/RwCalculators/GReWeightINukeKinematics.cxx @@ -44,6 +44,9 @@ bool GReWeightINukeKinematics::IsHandled(GSyst_t syst) const switch (syst) { case ( kINukeKinematicsTwkDial_NP_N ): case ( kINukeKinematicsTwkDial_PP_N ): + case ( kINukeKinematicsFixPiPro ): + case ( kINukeKinematicsPiProBias ): + case ( kINukeKinematicsPiProBiaswFix): return true; default: return false; @@ -156,6 +159,16 @@ double GReWeightINukeKinematics::CalcWeight(const EventRecord &event) { // std::cerr << "TLab: " << tlab_v << " costhcm: " << costhcm_v << " np/pp " << np_pp << " weight: " << weight << std::endl; + event_weight *= weight; + } + // pion production + else if (fsi_code == (int)kIHAFtPiProd) { + const GHepParticle &f1 = *event.Particle(p->FirstDaughter()); + const GHepParticle &f2 = *event.Particle(p->FirstDaughter()+1); + const GHepParticle &f3 = *event.Particle(p->LastDaughter()); + + double weight = fRwParams.CalcPionWeight(*p->P4(), *f1.P4(), *f2.P4(), *f3.P4()); + event_weight *= weight; } } diff --git a/src/RwCalculators/GReWeightINukeKinematicsParams.cxx b/src/RwCalculators/GReWeightINukeKinematicsParams.cxx index 49ae978a..3eb11aa3 100644 --- a/src/RwCalculators/GReWeightINukeKinematicsParams.cxx +++ b/src/RwCalculators/GReWeightINukeKinematicsParams.cxx @@ -16,6 +16,8 @@ #include "RwCalculators/GReWeightUtils.h" #include "RwFramework/GSystUncertainty.h" +#include + using namespace genie; using namespace genie::rew; using namespace genie::constants; @@ -23,11 +25,205 @@ using namespace genie::constants; //___________________________________________________________________________ GReWeightINukeKinematicsParams::GReWeightINukeKinematicsParams(): fNPwgt("pn"), - fPPwgt("pp") {} + fPPwgt("pp"), + fFixPiPro(false), + fBiasPiPro(0.) {} //___________________________________________________________________________ void GReWeightINukeKinematicsParams::SetTargetA(int A) {(void) A;} // nothing, for now +double ThreeBodyLorentzWeight(double E, double m1, double m2, double m3, double m23) { + double p1cm = sqrt((E*E - (m1 + m23)*(m1 + m23))*(E*E - (m1 - m23)*(m1 - m23))) / (2*E); + double p2cm2 = sqrt((m23*m23 - (m2+m3)*(m2+m3))*(m23*m23 - (m2-m3)*(m2-m3))) / (2*m23); + return p1cm*p2cm2; +} + +double ThreeBodyBias(double t1p, double bias) { + return exp(bias*t1p); +} + +double ThreeBodyLorentzWeightBias(double E, double m1, double m2, double m3, double m23, double t1p, double bias) { + return ThreeBodyLorentzWeight(E, m1, m2, m3, m23)*ThreeBodyBias(t1p, bias); +} + +// param should point to a list of 4 doubles: E, m1, m2, m3 +double ParamsThreeBodyLorentzWeight(double m23, void *param) { + double *M = (double *)param; + return ThreeBodyLorentzWeight(M[0], M[1], M[2], M[3], m23); +} + +double AvgThreeBodyLorentzWeight(double E, double m1, double m2, double m3) { + double M[4] = {E, m1, m2, m3}; + + gsl_integration_workspace * w = gsl_integration_workspace_alloc(1000); + gsl_function F; + F.function = &ParamsThreeBodyLorentzWeight; + F.params = (void*)M; + double result, error; + gsl_integration_qags(&F, m2+m3, E-m1, 0, 1.0e-7, 1000, w, &result, &error); + gsl_integration_workspace_free(w); + return result / (E - m1 - m2 - m3); +} + +// param should point to a list of 7 doubles: E, m1, m2, m3, m23, mp, bias +double ParamsThreeBodyLorentzWeightBiasInner(double costh1p, void *param) { + double *M = (double *)param; + double E = M[0]; + double m1 = M[1]; + double m2 = M[2]; + double m3 = M[3]; + double m23 = M[4]; + double mp = M[5]; + double bias = M[6]; + + // calculate the energy transfer from the projectile to particle 1 + double p1 = sqrt((E*E - (m1 + m23)*(m1 + m23))*(E*E - (m1 - m23)*(m1 - m23))) / (2*E); + double e1 = sqrt(p1*p1 + m1*m1); + // assume equal masses in collision + double ep = E / 2; + double pp = sqrt(ep*ep - mp*mp); + double t1p = m1*m1 + mp*mp - 2*(ep*e1 - pp*p1*costh1p); // == (pp - p1)**2 + + return ThreeBodyLorentzWeightBias(E, m1, m2, m3, m23, t1p, bias); +} + +// param should point to a list of 6 doubles: E, m1, m2, m3, mp, bias +double ParamsThreeBodyLorentzWeightBias(double m23, void *param) { + double *M = (double *)param; + double E = M[0]; + double m1 = M[1]; + double m2 = M[2]; + double m3 = M[3]; + double mp = M[4]; + double bias = M[5]; + + double Minner[7] = {E, m1, m2, m3, m23, mp, bias}; + + gsl_integration_workspace * w = gsl_integration_workspace_alloc(1000); + gsl_function F; + F.function = &ParamsThreeBodyLorentzWeightBiasInner; + F.params = (void*)Minner; + double result, error; + gsl_integration_qags(&F, -1, 1, 0, 1.0e-7, 1000, w, &result, &error); + gsl_integration_workspace_free(w); + return result / 2; // divide out integration range +} + +double AvgThreeBodyLorentzWeightBias(double E, double m1, double m2, double m3, double mp, double bias) { + double M[6] = {E, m1, m2, m3, mp, bias}; + + gsl_integration_workspace * w = gsl_integration_workspace_alloc(1000); + gsl_function F; + F.function = &ParamsThreeBodyLorentzWeightBias; + F.params = (void*)M; + double result, error; + gsl_integration_qags(&F, m2+m3, E-m1, 0, 1.0e-7, 1000, w, &result, &error); + gsl_integration_workspace_free(w); + return result / (E - m1 - m2 - m3); +} + +// param should point to a list of 5 doubles: E, m1, m23, mp, bias +double ParamsThreeBodyBiasInner(double costh1p, void *param) { + double *M = (double *)param; + double E = M[0]; + double m1 = M[1]; + double m23 = M[2]; + double mp = M[3]; + double bias = M[4]; + + // calculate the energy transfer from the projectile to particle 1 + double p1 = sqrt((E*E - (m1 + m23)*(m1 + m23))*(E*E - (m1 - m23)*(m1 - m23))) / (2*E); + double e1 = sqrt(p1*p1 + m1*m1); + // assume equal masses in collision + double ep = E / 2; + double pp = sqrt(ep*ep - mp*mp); + double t1p = m1*m1 + mp*mp - 2*(ep*e1 - pp*p1*costh1p); // == (pp - p1)**2 + + return ThreeBodyBias(t1p, bias); +} + +// param should point to a list of 4 doubles: E, m1, mp, bias +double ParamsThreeBodyBias(double m23, void *param) { + double *M = (double *)param; + double E = M[0]; + double m1 = M[1]; + double mp = M[2]; + double bias = M[3]; + + double Minner[5] = {E, m1, m23, mp, bias}; + + gsl_integration_workspace * w = gsl_integration_workspace_alloc(1000); + gsl_function F; + F.function = &ParamsThreeBodyBiasInner; + F.params = (void*)Minner; + double result, error; + gsl_integration_qags(&F, -1, 1, 0, 1.0e-7, 1000, w, &result, &error); + gsl_integration_workspace_free(w); + return result / 2; // divide out integration range +} + +double AvgThreeBodyBias(double E, double m1, double m2, double m3, double mp, double bias) { + double M[4] = {E, m1, mp, bias}; + + gsl_integration_workspace * w = gsl_integration_workspace_alloc(1000); + gsl_function F; + F.function = &ParamsThreeBodyBias; + F.params = (void*)M; + double result, error; + gsl_integration_qags(&F, m2+m3, E-m1, 0, 1.0e-7, 1000, w, &result, &error); + gsl_integration_workspace_free(w); + return result / (E - m1 - m2 - m3); +} + +double GReWeightINukeKinematicsParams::CalcPionWeight(TLorentzVector p, TLorentzVector f1, TLorentzVector f2, TLorentzVector f3) { + double weight = 1.; + + double E = (f1 + f2 + f3).M(); + if (fFixPiPro && fBiasPiPro != 0.) { + // Momentum transfer from parent to child nulceon + double t1p = (p - f1).M2(); + + // mass of 2-3 system + double m23 = (f2 + f3).M(); + + double bias_fix_weight = ThreeBodyLorentzWeightBias(E, f1.M(), f2.M(), f3.M(), m23, t1p, fBiasPiPro); + double ave_weight = AvgThreeBodyLorentzWeightBias(E, f1.M(), f2.M(), f3.M(), p.M(), fBiasPiPro); + weight = bias_fix_weight / ave_weight; + } + else if (fBiasPiPro != 0.) { + // Momentum transfer from parent to child nulceon + double t1p = (p - f1).M2(); + double bias_weight = ThreeBodyBias(t1p, fBiasPiPro); + double ave_weight = AvgThreeBodyBias(E, f1.M(), f2.M(), f3.M(), p.M(), fBiasPiPro); + weight = bias_weight / ave_weight; + } + else if (fFixPiPro) { + // mass of 2-3 system + double m23 = (f2 + f3).M(); + + // fix angles + TVector3 boost_to_COM = -(f1+f2+f3).BoostVector(); + TLorentzVector f1_COM(f1); + f1_COM.Boost(boost_to_COM); + TVector3 z_dir = (f1+f2+f3).Vect().Unit(); + double costh1 = f1_COM.Vect().Unit().Dot(z_dir); + double sinth1 = sqrt(1 - costh1*costh1); + + TVector3 boost_to_COM2 = -(f2+f3).BoostVector(); + TLorentzVector f2_COM2(f2); + f2_COM2.Boost(boost_to_COM2); + TVector3 z_dir2 = (f2+f3).Vect().Unit(); + double costh2 = f2_COM2.Vect().Unit().Dot(z_dir2); + double sinth2 = sqrt(1 - costh2*costh2); + + double fix_weight = ThreeBodyLorentzWeight(E, f1.M(), f2.M(), f3.M(), m23)*sinth1*sinth2; + double ave_weight = AvgThreeBodyLorentzWeight(E, f1.M(), f2.M(), f3.M()) * (2./M_PI) * (2./M_PI); + weight = fix_weight / ave_weight; + } + + return weight; +} + double GReWeightINukeKinematicsParams::CalcWeight(int np_pp, float tlab, float costhcm) { if (np_pp == 0) return fNPwgt.CalcWeight(tlab, costhcm); @@ -44,12 +240,26 @@ void GReWeightINukeKinematicsParams::SetTwkDial(GSyst_t syst, double val) else if (syst == kINukeKinematicsTwkDial_PP_N) { fPPwgt.SetUniverse((int)val); } + else if (syst == kINukeKinematicsFixPiPro) { + if (val != 0.) fFixPiPro = true; + } + else if (syst == kINukeKinematicsPiProBias) { + GSystUncertainty * fracerr = GSystUncertainty::Instance(); + fBiasPiPro = val*fracerr->OneSigmaErr(kINukeKinematicsPiProBias); + } + else if (syst == kINukeKinematicsPiProBiaswFix) { + GSystUncertainty * fracerr = GSystUncertainty::Instance(); + fBiasPiPro = val*fracerr->OneSigmaErr(kINukeKinematicsPiProBiaswFix); + fFixPiPro = true; + } } //___________________________________________________________________________ void GReWeightINukeKinematicsParams::Reset() { fNPwgt.Reset(); fPPwgt.Reset(); + fFixPiPro = false; + fBiasPiPro = 0.; } void GReWeightINukeKinematicsParams::ReBounce::Reset() { diff --git a/src/RwCalculators/GReWeightINukeKinematicsParams.h b/src/RwCalculators/GReWeightINukeKinematicsParams.h index 2c634935..f876dcc9 100644 --- a/src/RwCalculators/GReWeightINukeKinematicsParams.h +++ b/src/RwCalculators/GReWeightINukeKinematicsParams.h @@ -2,6 +2,7 @@ #define _G_REWEIGHT_INTRANUKEKINEMATICS_PARAMS_H_ #include +#include // GENIE/Reweight includes #include "RwFramework/GSyst.h" @@ -20,6 +21,8 @@ namespace rew { void SetTwkDial (GSyst_t s, double val); ///< void SetTargetA (int target_A); ///< Set the mass number of the hit nucleus double CalcWeight (int np_pp, float tlab, float costhcm); + double CalcPionWeight (TLorentzVector p, TLorentzVector f1, TLorentzVector f2, TLorentzVector f3); + struct ReBounce { void SetUniverse(int u); @@ -39,6 +42,8 @@ namespace rew { private: ReBounce fNPwgt; ReBounce fPPwgt; + bool fFixPiPro; + double fBiasPiPro; }; //GReWeightINukeKinematicsParams diff --git a/src/RwFramework/GSyst.cxx b/src/RwFramework/GSyst.cxx index fba65ede..588414da 100644 --- a/src/RwFramework/GSyst.cxx +++ b/src/RwFramework/GSyst.cxx @@ -97,6 +97,9 @@ std::map GSyst::BuildGSystToStringMap() { temp_map[ kINukeTwkDial_FrPiProd_N ] = "FrPiProd_N"; temp_map[ kINukeKinematicsTwkDial_NP_N ] = "FrKin_NP_N"; temp_map[ kINukeKinematicsTwkDial_PP_N ] = "FrKin_PP_N"; + temp_map[ kINukeKinematicsFixPiPro ] = "FrKin_PiProFix_N"; + temp_map[ kINukeKinematicsPiProBiaswFix ] = "FrKin_PiProBiaswFix_N"; + temp_map[ kINukeKinematicsPiProBias ] = "FrKin_PiProBias_N"; temp_map[ kSystNucl_CCQEPauliSupViaKF ] = "CCQEPauliSupViaKF"; temp_map[ kSystNucl_CCQEMomDistroFGtoSF ] = "CCQEMomDistroFGtoSF"; temp_map[ kRDcyTwkDial_BR1gamma ] = "RDecBR1gamma"; diff --git a/src/RwFramework/GSyst.h b/src/RwFramework/GSyst.h index a93bf389..374e87df 100644 --- a/src/RwFramework/GSyst.h +++ b/src/RwFramework/GSyst.h @@ -136,6 +136,9 @@ typedef enum EGSyst { kINukeKinematicsTwkDial_NP_N, ///< tweak scattering angle of NP scatters for nucleons kINukeKinematicsTwkDial_PP_N, ///< tweak scattering angle of PP and NN scatters for nucleons + kINukeKinematicsFixPiPro, ///< tweak to fix pion production kinematics + kINukeKinematicsPiProBias, ///< tweak to introduce bias to pion production momenta + kINukeKinematicsPiProBiaswFix, ///< tweak to introduce bias to pion production momenta, including lorentz weight fixing // // Nuclear model From 9deb658a56c3263130e85472e35f2f9cf59db021 Mon Sep 17 00:00:00 2001 From: Gray Putnam Date: Wed, 6 Aug 2025 19:39:00 -0500 Subject: [PATCH 7/7] Print error message if user configures bad pion production kinematic bias. --- src/RwCalculators/GReWeightINukeKinematicsParams.cxx | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/RwCalculators/GReWeightINukeKinematicsParams.cxx b/src/RwCalculators/GReWeightINukeKinematicsParams.cxx index 3eb11aa3..3d6c5b40 100644 --- a/src/RwCalculators/GReWeightINukeKinematicsParams.cxx +++ b/src/RwCalculators/GReWeightINukeKinematicsParams.cxx @@ -244,6 +244,8 @@ void GReWeightINukeKinematicsParams::SetTwkDial(GSyst_t syst, double val) if (val != 0.) fFixPiPro = true; } else if (syst == kINukeKinematicsPiProBias) { + if (val < 0) LOG("ReW", pERROR) << "It is unphysical to bias pion production with a negative 'B' value. Therefore, only positive sigma values are physical."; + GSystUncertainty * fracerr = GSystUncertainty::Instance(); fBiasPiPro = val*fracerr->OneSigmaErr(kINukeKinematicsPiProBias); }