From 29823701348f6f23cf962a025a3c18e702d032c1 Mon Sep 17 00:00:00 2001 From: mking Date: Wed, 17 Jun 2026 15:30:48 -0500 Subject: [PATCH 1/2] Bugfix: For HA Pion Absorption, adding in the pion mass to the final state energy budget for the multinucleon absorption mechanism --- .../HadronTransport/HAIntranuke2018.cxx | 17 +++++++++------ .../HadronTransport/HAIntranuke2025.cxx | 21 ++++++++++++------- 2 files changed, 24 insertions(+), 14 deletions(-) diff --git a/src/Physics/HadronTransport/HAIntranuke2018.cxx b/src/Physics/HadronTransport/HAIntranuke2018.cxx index 682470dde1..0bf2c8405a 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. @@ -1241,10 +1241,13 @@ void HAIntranuke2018::Inelastic( // simple for now, each (of 5) in hadron cluster has 1/5 of mom and KE double probM = pLib->Find(pdgc) ->Mass(); - probM -= .025; // BE correction + double BE_correction_per_nucleon = .005; //GeV + probM -= BE_correction_per_nucleon*5; // BE correction TVector3 pP3 = p->P4()->Vect() * (1./5.); - double probKE = p->P4()->E() -probM; - double clusKE = probKE * (1./5.); + double probEAvail; + if ( pdgc==kPdgPiP || pdgc==kPdgPi0 || pdgc==kPdgPiM) probEAvail = p->P4()->E() + BE_correction_per_nucleon*5; //pion probe -- includes mass, BE correction for 4 nucleons + else probEAvail = p->P4()->E() -probM; //proton or neutron probe -- does not include mass + double clusKE = probEAvail * (1./5.); TLorentzVector clusP4(pP3,clusKE); //no mass LOG("HAIntranuke2018",pINFO) << "probM = " << probM << " ;clusKE= " << clusKE; TLorentzVector X4(*p->X4()); @@ -1391,8 +1394,10 @@ void HAIntranuke2018::Inelastic( double probM = pLib->Find(pdgc) ->Mass(); double probBE = (np+nn)*.005; // BE correction TVector3 pP3 = p->P4()->Vect(); - double probKE = p->P4()->E() - (probM - probBE); - double clusKE = probKE; // + np*0.9383 + nn*.9396; + double probEAvail; + if ( pdgc==kPdgPiP || pdgc==kPdgPi0 || pdgc==kPdgPiM) probEAvail = p->P4()->E() + probBE; //pion probe -- includes mass, BE correction for 4 nucleons + else probEAvail = p->P4()->E() - (probM - probBE); //proton or neutron probe -- does not include mass + double clusKE = probEAvail; // + np*0.9383 + nn*.9396; TLorentzVector clusP4(pP3,clusKE); //no mass is correct LOG("HAIntranuke2018",pINFO) << "probM = " << probM << " ;clusKE= " << clusKE; TLorentzVector X4(*p->X4()); diff --git a/src/Physics/HadronTransport/HAIntranuke2025.cxx b/src/Physics/HadronTransport/HAIntranuke2025.cxx index cb5027f69f..bf936b0ce5 100644 --- a/src/Physics/HadronTransport/HAIntranuke2025.cxx +++ b/src/Physics/HadronTransport/HAIntranuke2025.cxx @@ -1258,11 +1258,14 @@ void HAIntranuke2025::Inelastic( // simple for now, each (of 5) in hadron cluster has 1/5 of mom and KE double probM = pLib->Find(pdgc) ->Mass(); - probM -= .025; // BE correction - TVector3 pP3 = p->P4()->Vect() * (1./5.); - double probKE = p->P4()->E() -probM; - double clusKE = probKE * (1./5.); - TLorentzVector clusP4(pP3,clusKE); //no mass + double BE_correction_per_nucleon = .005; //GeV + probM -= BE_correction_per_nucleon*5; // BE correction + TVector3 pP3 = p->P4()->Vect() * (1./5.); + double probEAvail; + if ( pdgc==kPdgPiP || pdgc==kPdgPi0 || pdgc==kPdgPiM) probEAvail = p->P4()->E() + BE_correction_per_nucleon*5; //pion probe -- includes mass, BE correction for 4 nucleons + else probEAvail = p->P4()->E() -probM; //proton or neutron probe -- does not include mass + double clusKE = probEAvail * (1./5.); + TLorentzVector clusP4(pP3,clusKE); //no mass LOG("HAIntranuke2025",pINFO) << "probM = " << probM << " ;clusKE= " << clusKE; TLorentzVector X4(*p->X4()); GHepStatus_t ist = kIStNucleonClusterTarget; @@ -1408,9 +1411,11 @@ void HAIntranuke2025::Inelastic( double probM = pLib->Find(pdgc) ->Mass(); double probBE = (np+nn)*.005; // BE correction TVector3 pP3 = p->P4()->Vect(); - double probKE = p->P4()->E() - (probM - probBE); - double clusKE = probKE; // + np*0.9383 + nn*.9396; - TLorentzVector clusP4(pP3,clusKE); //no mass is correct + double probEAvail; + if ( pdgc==kPdgPiP || pdgc==kPdgPi0 || pdgc==kPdgPiM) probEAvail = p->P4()->E() + probBE; //pion probe -- includes mass, BE correction for 4 nucleons + else probEAvail = p->P4()->E() - (probM - probBE); //proton or neutron probe -- does not include mass + double clusKE = probEAvail; // + np*0.9383 + nn*.9396; + TLorentzVector clusP4(pP3,clusKE); //no mass is correct LOG("HAIntranuke2025",pINFO) << "probM = " << probM << " ;clusKE= " << clusKE; TLorentzVector X4(*p->X4()); GHepStatus_t ist = kIStNucleonClusterTarget; From 781694849a95fd52151572831bd0a2ba3ce4489c Mon Sep 17 00:00:00 2001 From: mking Date: Wed, 17 Jun 2026 16:08:51 -0500 Subject: [PATCH 2/2] Updating bugfix code, fixing indentation and utilizing IsPion() function --- src/Physics/HadronTransport/HAIntranuke2018.cxx | 8 ++++---- src/Physics/HadronTransport/HAIntranuke2025.cxx | 8 ++++---- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/Physics/HadronTransport/HAIntranuke2018.cxx b/src/Physics/HadronTransport/HAIntranuke2018.cxx index 0bf2c8405a..47dcc7865a 100644 --- a/src/Physics/HadronTransport/HAIntranuke2018.cxx +++ b/src/Physics/HadronTransport/HAIntranuke2018.cxx @@ -1245,10 +1245,10 @@ void HAIntranuke2018::Inelastic( probM -= BE_correction_per_nucleon*5; // BE correction TVector3 pP3 = p->P4()->Vect() * (1./5.); double probEAvail; - if ( pdgc==kPdgPiP || pdgc==kPdgPi0 || pdgc==kPdgPiM) probEAvail = p->P4()->E() + BE_correction_per_nucleon*5; //pion probe -- includes mass, BE correction for 4 nucleons + if ( genie::pdg::IsPion(pgdc) ) probEAvail = p->P4()->E() + BE_correction_per_nucleon*5; //pion probe -- includes mass, BE correction for 4 nucleons else probEAvail = p->P4()->E() -probM; //proton or neutron probe -- does not include mass double clusKE = probEAvail * (1./5.); - TLorentzVector clusP4(pP3,clusKE); //no mass + TLorentzVector clusP4(pP3,clusKE); //using no mass here is correct LOG("HAIntranuke2018",pINFO) << "probM = " << probM << " ;clusKE= " << clusKE; TLorentzVector X4(*p->X4()); GHepStatus_t ist = kIStNucleonClusterTarget; @@ -1395,10 +1395,10 @@ void HAIntranuke2018::Inelastic( double probBE = (np+nn)*.005; // BE correction TVector3 pP3 = p->P4()->Vect(); double probEAvail; - if ( pdgc==kPdgPiP || pdgc==kPdgPi0 || pdgc==kPdgPiM) probEAvail = p->P4()->E() + probBE; //pion probe -- includes mass, BE correction for 4 nucleons + if ( genie::pdg::IsPion(pgdc) ) probEAvail = p->P4()->E() + probBE; //pion probe -- includes mass, BE correction for 4 nucleons else probEAvail = p->P4()->E() - (probM - probBE); //proton or neutron probe -- does not include mass double clusKE = probEAvail; // + np*0.9383 + nn*.9396; - TLorentzVector clusP4(pP3,clusKE); //no mass is correct + TLorentzVector clusP4(pP3,clusKE); //using no mass here is correct LOG("HAIntranuke2018",pINFO) << "probM = " << probM << " ;clusKE= " << clusKE; TLorentzVector X4(*p->X4()); GHepStatus_t ist = kIStNucleonClusterTarget; diff --git a/src/Physics/HadronTransport/HAIntranuke2025.cxx b/src/Physics/HadronTransport/HAIntranuke2025.cxx index bf936b0ce5..5885a98d49 100644 --- a/src/Physics/HadronTransport/HAIntranuke2025.cxx +++ b/src/Physics/HadronTransport/HAIntranuke2025.cxx @@ -1262,10 +1262,10 @@ void HAIntranuke2025::Inelastic( probM -= BE_correction_per_nucleon*5; // BE correction TVector3 pP3 = p->P4()->Vect() * (1./5.); double probEAvail; - if ( pdgc==kPdgPiP || pdgc==kPdgPi0 || pdgc==kPdgPiM) probEAvail = p->P4()->E() + BE_correction_per_nucleon*5; //pion probe -- includes mass, BE correction for 4 nucleons + if ( genie::pdg::IsPion(pgdc) ) probEAvail = p->P4()->E() + BE_correction_per_nucleon*5; //pion probe -- includes mass, BE correction for 4 nucleons else probEAvail = p->P4()->E() -probM; //proton or neutron probe -- does not include mass double clusKE = probEAvail * (1./5.); - TLorentzVector clusP4(pP3,clusKE); //no mass + TLorentzVector clusP4(pP3,clusKE); //using no mass here is correct LOG("HAIntranuke2025",pINFO) << "probM = " << probM << " ;clusKE= " << clusKE; TLorentzVector X4(*p->X4()); GHepStatus_t ist = kIStNucleonClusterTarget; @@ -1412,10 +1412,10 @@ void HAIntranuke2025::Inelastic( double probBE = (np+nn)*.005; // BE correction TVector3 pP3 = p->P4()->Vect(); double probEAvail; - if ( pdgc==kPdgPiP || pdgc==kPdgPi0 || pdgc==kPdgPiM) probEAvail = p->P4()->E() + probBE; //pion probe -- includes mass, BE correction for 4 nucleons + if ( genie::pdg::IsPion(pgdc) ) probEAvail = p->P4()->E() + probBE; //pion probe -- includes mass, BE correction for 4 nucleons else probEAvail = p->P4()->E() - (probM - probBE); //proton or neutron probe -- does not include mass double clusKE = probEAvail; // + np*0.9383 + nn*.9396; - TLorentzVector clusP4(pP3,clusKE); //no mass is correct + TLorentzVector clusP4(pP3,clusKE); //using no mass here is correct LOG("HAIntranuke2025",pINFO) << "probM = " << probM << " ;clusKE= " << clusKE; TLorentzVector X4(*p->X4()); GHepStatus_t ist = kIStNucleonClusterTarget;