diff --git a/src/Apps/gRwght1Param.cxx b/src/Apps/gRwght1Param.cxx index d975d7d90..95ccd5a7c 100644 --- a/src/Apps/gRwght1Param.cxx +++ b/src/Apps/gRwght1Param.cxx @@ -123,6 +123,9 @@ #include "RwCalculators/GReWeightNuXSecCCQEvec.h" #include "RwCalculators/GReWeightNuXSecNCRES.h" #include "RwCalculators/GReWeightNuXSecDIS.h" +#include "RwCalculators/GReWeighthA2025.h" + + #include "RwCalculators/GReWeightINukeParams.h" #include "RwCalculators/GReWeightNuXSecNC.h" @@ -192,6 +195,15 @@ int main(int argc, char ** argv) const float twk_dial_min = gOptMinTwk; const float twk_dial_max = gOptMaxTwk; const float twk_dial_step = (twk_dial_max - twk_dial_min) / (n_points-1); + + LOG("grwght1scan", pNOTICE) + << "\nSummary of tweak dial parameters:" + << "\n - n_points = " << n_points + << "\n - twk_dial_min = " << twk_dial_min + << "\n - twk_dial_max = " << twk_dial_max + << "\n - twk_dial_step = " << twk_dial_step + << "\n"; + // Work-out the range of events to process Long64_t nfirst = 0; @@ -248,7 +260,8 @@ int main(int argc, char ** argv) rw.AdoptWghtCalc( "hadro_fzone", new GReWeightFZone ); rw.AdoptWghtCalc( "hadro_intranuke", new GReWeightINuke ); rw.AdoptWghtCalc( "hadro_agky", new GReWeightAGKY ); - + rw.AdoptWghtCalc( "hA_2025_cex_weight", new GReWeighthA2025 ); + // GReWeightDISNuclMod::CalcWeight() not implemented - don't try to use it .. // will return 1 if tweak dial = 0, hard fail otherwise rw.AdoptWghtCalc( "nuclear_dis", new GReWeightDISNuclMod ); diff --git a/src/RwCalculators/GReWeighthA2025.cxx b/src/RwCalculators/GReWeighthA2025.cxx new file mode 100644 index 000000000..3bf8b3980 --- /dev/null +++ b/src/RwCalculators/GReWeighthA2025.cxx @@ -0,0 +1,297 @@ +//____________________________________________________________________________ +/* + Copyright (c) 2003-2025, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + + Authors: Mohamed Ismail + University of Pittsburgh + +*/ +//____________________________________________________________________________ + +#include +#include + +#include +#include +#include +#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" +#include "Physics/HadronTransport/HAIntranuke2018.h" +#include "Physics/HadronTransport/Intranuke2018.h" +#include "Physics/HadronTransport/INukeHadroData2018.h" +#include "Physics/HadronTransport/INukeHadroFates2018.h" +#include "Physics/HadronTransport/INukeUtils2018.h" + +#include "Physics/HadronTransport/INukeHadroData2025.h" +#include "Physics/HadronTransport/INukeHadroFates2025.h" +#include "Physics/HadronTransport/INukeUtils2025.h" +#include "Physics/HadronTransport/HAIntranuke2025.h" +#include "Physics/HadronTransport/Intranuke2025.h" + +// GENIE/Reweight includes +#include "RwCalculators/GReWeighthA2025.h" +#include "RwCalculators/GReWeightUtils.h" +#include "RwFramework/GSystUncertainty.h" + +using namespace genie; +using namespace genie::rew; + +//_______________________________________________________________________________________ +GReWeighthA2025::GReWeighthA2025() : +GReWeightModel("IntraNuke") +{ +#ifdef _G_REWEIGHT_INUKE_DEBUG_NTP_ + fTestFile = new TFile("./intranuke_reweight_test.root","recreate"); + fTestNtp = new TNtuple("testntp","","pdg:E:mfp_twk_dial:d:d_mfp:fate:interact:w_mfp:w_fate"); +#endif + + // Look up the FSI model for the current tune. Also check whether FSIs are + // actually enabled. + AlgConfigPool* conf_pool = AlgConfigPool::Instance(); + Registry* gpl = conf_pool->GlobalParameterList(); + RgAlg fsi_alg = gpl->GetAlg( "HadronTransp-Model" ); + bool fsi_enabled = gpl->GetBool( "HadronTransp-Enable" ); + + if ( !fsi_enabled ) { + LOG( "ReW", pFATAL ) << "FSIs are not enabled for the current tune." + << " Refusing to reweight FSIs."; + std::exit( 1 ); + } + + AlgId id( fsi_alg ); + + AlgFactory* algf = AlgFactory::Instance(); + + Algorithm* alg = algf->AdoptAlgorithm( id ); + fFSIModel = dynamic_cast< HAIntranuke2018* >( alg ); + + if ( !fFSIModel ) { + LOG( "ReW", pFATAL ) << "Reweighting events produced with the FSI model " + << fsi_alg << " is not currently supported."; + std::exit( 1 ); + } + + fFSIModel->AdoptSubstructure(); +} +//_______________________________________________________________________________________ +GReWeighthA2025::~GReWeighthA2025() +{ +#ifdef _G_REWEIGHT_INUKE_DEBUG_NTP_ + assert(fTestFile); + assert(fTestNtp); + fTestFile->cd(); + fTestNtp->Write(); + fTestFile->Close(); + delete fTestFile; + //delete fTestNtp; +#endif +} +//_______________________________________________________________________________________ +bool GReWeighthA2025::IsHandled(GSyst_t syst) const +{ + bool handle; + + switch(syst) { + case ( kINukehA2025_cex ) : + handle = true; + break; + + default: + handle = false; + } + + return handle; +} +//_______________________________________________________________________________________ +bool GReWeighthA2025::AppliesTo(ScatteringType_t type, bool /*is_cc*/) const +{ + if (type != kScCoherentProduction ) { + return true; + } + return false; +} +//_______________________________________________________________________________________ + +void GReWeighthA2025::SetSystematic(GSyst_t syst, double val) +{ + // if(this->IsHandled(syst)) { + // fINukeRwParams.SetTwkDial(syst, val); + // } +} +//_______________________________________________________________________________________ +void GReWeighthA2025::Reset(void) +{ + // fINukeRwParams.Reset(); + // this->Reconfigure(); +} +//_______________________________________________________________________________________ +void GReWeighthA2025::Reconfigure(void) +{ + // fINukeRwParams.Reconfigure(); +} + +//_______________________________________________________________________________________ +double GReWeighthA2025::CalcWeight(const EventRecord & event) +{ + // 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; + + //fINukeRwParams.SetTargetA( A ); + + // Get the pre-FSI nuclear remnant. The A and Z values for this particle + // (distinct from both the post-FSI remnant and the target nucleus) + // are used to compute mean free paths in Intranuke2018. + GHepParticle* pre_fsi_remnant = 0; + TObjArrayIter piter( &event ); + // First loop over all particles in the event record and find the + // final-state nuclear remnant (i.e., the post-FSI nuclear remnant) + while( GHepParticle* p = dynamic_cast(piter.Next()) ) { + if ( p->Status() == genie::kIStFinalStateNuclearRemnant ) { + // The pre-FSI nuclear remnant is set to be the first mother for + // the post-FSI one by Intranuke2018 + int mother_idx = p->FirstMother(); + pre_fsi_remnant = event.Particle( mother_idx ); + break; + } + } + + // Skip this event if we couldn't find a pre-FSI nuclear remnant that + // is an ion (something went wrong in the search above) + bool pre_fsi_remnant_ok = false; + if ( pre_fsi_remnant ) { + pre_fsi_remnant_ok = genie::pdg::IsIon( pre_fsi_remnant->Pdg() ); + } + if ( !pre_fsi_remnant_ok ) { + LOG( "ReW", pWARN ) << "Could not find a suitable pre-FSI remnant" + << " in GReWeighthA2025::CalcWeight()"; + return 1.; + } + + // Store the nucleon and proton numbers for the pre-FSI remnant. + // These will be used for mean free path calculations in + // genie::utils::rew::MeanFreePathWeight + int remnA = pre_fsi_remnant->A(); + int remnZ = pre_fsi_remnant->Z(); + LOG( "ReW", pDEBUG ) << "Found pre-FSI remnant with A = " << remnA + << ", Z = " << remnZ << ". Target had A = " << A << ", Z = " << Z; + + + double event_weight = 1.0; + + // 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_pion && !is_nucleon && !is_kaon) + { + continue; + } + + + // Get 4-momentum and 4-position + TLorentzVector x4 (p->Vx(), p->Vy(), p->Vz(), 0. ); + TLorentzVector p4 (p->Px(), p->Py(), p->Pz(), p->E()); + + int fsi_code = p->RescatterCode(); + bool escaped = (fsi_code == (int)kIHAFtNoInteraction); + bool interacted = !escaped; + + // new + double fate_frac2018 = 0.0; + double fate_frac2025 = 0.0; + if (is_pion and fsi_code > 1 ){ + INukeHadroData2018 * hd2018 = INukeHadroData2018::Instance(); + INukeHadroData2025 * hd2025 = INukeHadroData2025::Instance(); + + // convert to MeV and + + double KE = p4.Energy() - p4.M(); // kinetic energy + double ke = KE / units::MeV; + + // get particle fate + auto fate_rescatter = p->RescatterCode(); + + // get fate frac + + if (fate_rescatter == kIHAFtAbs){ + fate_frac2018 = hd2018->FracADep(pdgc, kIHAFtAbs, ke, remnA); + fate_frac2025 = hd2025->FracADep(pdgc, kIHAFtAbs, ke, remnA); + } + + else if (fate_rescatter == kIHAFtInelas){ + fate_frac2018 = hd2018->FracADep(pdgc, kIHAFtInelas, ke, remnA); + fate_frac2025 = hd2025->FracADep(pdgc, kIHAFtInelas, ke, remnA); + } + + else if (fate_rescatter == kIHAFtCEx){ + fate_frac2018 = hd2018->FracADep(pdgc, kIHAFtCEx, ke, remnA); + fate_frac2025 = hd2025->FracADep(pdgc, kIHAFtCEx, ke, remnA); + } + + else if (fate_rescatter == kIHAFtPiProd){ + fate_frac2018 = hd2018->FracADep(pdgc, kIHAFtPiProd, ke, remnA); + fate_frac2025 = hd2025->FracADep(pdgc, kIHAFtPiProd, ke, remnA); + } + + LOG("ReW", pNOTICE) + << "New:::::Reweighted hadron at position = " << ip + << " with PDG code = " << pdgc + << ", FSI code = " << fsi_code + << ", KE= " << ke + << ", A= " << remnA + << " (" << INukeHadroFates::AsString((INukeFateHA_t)fsi_code) << ") :" + << " frac2018 = " << fate_frac2018 + <<", frac2025 = " << fate_frac2025; + + + + } + double frac_ratio = 0.0; + if (fate_frac2018 != 0.0) { + frac_ratio = fate_frac2025 / fate_frac2018; + } else { + frac_ratio =1.0 ; + } + + + + + + + // Update the current event weight + event_weight *= frac_ratio; + + + + }//particle loop + + return event_weight; +} + diff --git a/src/RwCalculators/GReWeighthA2025.h b/src/RwCalculators/GReWeighthA2025.h new file mode 100644 index 000000000..7da3e44c3 --- /dev/null +++ b/src/RwCalculators/GReWeighthA2025.h @@ -0,0 +1,87 @@ +//____________________________________________________________________________ +/*! + +\class genie::rew::GReWeightINuke + +\brief Reweighting GENIE INTRANUKE/hA hadron transport model. + + The reweighting code considers two sets of physics changes: + * Change in the hadron mean free path \lambda, i.e. change in + the total rescattering probability P_{rescat}. + * Changes in probabilty for rescattering mode X, given a fixed + total rescattering probability P(X | \lambda). + X = {elastic, inelastic, charge exchange, pion production, + absorption + multi-nucleon knockout}. + + Physics changes are considered separately for pions and nucleons. + Unitarity is explicitly conserved. + +\author Mohamed Ismail + Univeirsty of Pittsburgh + + +\created Sep 2025 + +\cpright Copyright (c) 2003-2025, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org +*/ +//____________________________________________________________________________ + +#ifndef GREWEIGHTHA2025_H +#define GREWEIGHTHA2025_H + +//#define _G_REWEIGHT_INUKE_DEBUG_NTP_ + +// GENIE/Reweight includes +#include "RwCalculators/GReWeightModel.h" +//#include "RwCalculators/GReWeightINukeParams.h" + +using namespace genie::rew; +using namespace genie; + +class TFile; +class TNtuple; +class TLorentzVector; + +namespace genie { + + class HAIntranuke2018; + +namespace rew { + + class GReWeighthA2025 : public GReWeightModel + { + public: + GReWeighthA2025(); + ~GReWeighthA2025(); + + // implement the GReWeightI interface + bool AppliesTo (ScatteringType_t type, bool is_cc) 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: + +// void CalcDeltaAZ( const EventRecord& event, const GHepParticle& p, + // int& deltaA, int& deltaZ ); + + // void UpdateRemnantAZ( int deltaA, int deltaZ ); + + // GReWeightINukeParams fINukeRwParams; + + HAIntranuke2018* fFSIModel; + +#ifdef _G_REWEIGHT_INUKE_DEBUG_NTP_ + TFile * fTestFile; + TNtuple * fTestNtp; +#endif + + }; + +} // rew +} // genie + +#endif diff --git a/src/RwFramework/GSyst.cxx b/src/RwFramework/GSyst.cxx index d2cfd0cf4..45b25fd6d 100644 --- a/src/RwFramework/GSyst.cxx +++ b/src/RwFramework/GSyst.cxx @@ -138,6 +138,9 @@ std::map GSyst::BuildGSystToStringMap() { temp_map[ kXSecTwkDial_ZExpELFF_BN2 ] = "ZExpELFFBN2CCQE"; temp_map[ kXSecTwkDial_ZExpELFF_BN3 ] = "ZExpELFFBN3CCQE"; temp_map[ kXSecTwkDial_ZExpELFF_BN4 ] = "ZExpELFFBN4CCQE"; + temp_map[ kINukehA2025_cex ] = "rw_hA2025_cex"; + + return temp_map; } diff --git a/src/RwFramework/GSyst.h b/src/RwFramework/GSyst.h index bec7ddda5..53c2fe107 100644 --- a/src/RwFramework/GSyst.h +++ b/src/RwFramework/GSyst.h @@ -132,7 +132,7 @@ typedef enum EGSyst { kINukeTwkDial_FrInel_N, ///< tweak inelastic probability for nucleons, for given total rescattering probability 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 - + // // Nuclear model // @@ -240,6 +240,8 @@ typedef enum EGSyst { kXSecTwkDial_ZExpELFF_BN3, kXSecTwkDial_ZExpELFF_BN4, + // hA2018 FSI reweight dial + kINukehA2025_cex, // // Misc //