From a8d06c2b96cb1d0503b0615f37fdf4e2dd7d8e15 Mon Sep 17 00:00:00 2001 From: Marco Del Tutto Date: Tue, 9 Dec 2025 14:10:45 -0600 Subject: [PATCH 1/3] Make flux calculator compatible with g4bnb flux --- .../Calculators/BNBFlux/FluxCalcPrep.cxx | 47 ++++++++++++++++--- .../Calculators/BNBFlux/FluxCalcPrep.h | 4 ++ 2 files changed, 44 insertions(+), 7 deletions(-) diff --git a/sbncode/SBNEventWeight/Calculators/BNBFlux/FluxCalcPrep.cxx b/sbncode/SBNEventWeight/Calculators/BNBFlux/FluxCalcPrep.cxx index 640a84769..b1677ce35 100644 --- a/sbncode/SBNEventWeight/Calculators/BNBFlux/FluxCalcPrep.cxx +++ b/sbncode/SBNEventWeight/Calculators/BNBFlux/FluxCalcPrep.cxx @@ -213,8 +213,16 @@ namespace sbn { // std::cout<<"SBNEventWeight : getweight for the "< > mcFluxHandle; - e.getByLabel(fGeneratorModuleLabel,mcFluxHandle); + e.getByLabel(fGeneratorModuleLabel, mcFluxHandle); std::vector const& fluxlist = *mcFluxHandle; + + art::Handle< std::vector > dk2nuHandle; + std::vector const* dk2nu_v = nullptr; + e.getByLabel(fGeneratorModuleLabel, dk2nuHandle); + if (dk2nuHandle.isValid()) { + dk2nu_v = dk2nuHandle.product(); + } + //or do the above 3 lines in one line // auto const& mclist = *e.getValidHandle>(fGeneratorModuleLabel); @@ -287,9 +295,35 @@ namespace sbn { // First let's check that the parent of the neutrino we are looking for is // the particle we intended it to be, if not set all weights to 1 - if (std::find(fprimaryHad.begin(), fprimaryHad.end(),(fluxlist[inu].ftptype)) == fprimaryHad.end() ){//if it does not contain any particles we need get 1 + + simb::MCFlux flux; + flux.ftptype = fluxlist[inu].ftptype; + flux.ftpx = fluxlist[inu].ftpx; + flux.ftpy = fluxlist[inu].ftpy; + flux.ftpz = fluxlist[inu].ftpz; + + // int tptype = fluxlist[inu].ftptype; + + // If Dk2Nu flux, use ancestors to evaluate tptype + if (fluxlist[inu].fFluxType == simb::kDk2Nu) { + + // Find the first inelastic interaction + int firstInelastic = 0; + while (dk2nu_v->at(inu).ancestor[firstInelastic].proc.find("HadronInelastic")==std::string::npos||dk2nu_v->at(inu).ancestor[firstInelastic].proc.find("QEBooNE")!=std::string::npos) firstInelastic++; + + // Get the parent/ancestor (this should be the secondary in the p+Be interaction) + flux.ftptype = dk2nu_v->at(inu).ancestor[firstInelastic].pdg; + flux.ftpx = dk2nu_v->at(inu).ancestor[firstInelastic].startpx; + flux.ftpy = dk2nu_v->at(inu).ancestor[firstInelastic].startpy; + flux.ftpz = dk2nu_v->at(inu).ancestor[firstInelastic].startpz; + + } + + + if (std::find(fprimaryHad.begin(), fprimaryHad.end(),(flux.ftptype)) == fprimaryHad.end() ){//if it does not contain any particles we need get 1 weights.resize( NUni); std::fill(weights.begin(), weights.end(), 1); + std::cout << "We don't need this parent, returning." << std::endl; return weights;//done, all 1 }// Hadronic parent check @@ -301,19 +335,18 @@ namespace sbn { std::vector< float > Vrandom = (fParameterSet.fParameterMap.begin())->second;//vector of random # std::vector< float > subVrandom;//sub-vector of random numbers; if( CalcType == "PrimaryHadronNormalization"){//Normalization - test_weight = PHNWeightCalc(fluxlist[inu], Vrandom[i]); + test_weight = PHNWeightCalc(flux, Vrandom[i]); } else { subVrandom = {Vrandom.begin()+i*FitCov->GetNcols(), Vrandom.begin()+(i+1)*FitCov->GetNcols()}; if( CalcType == "PrimaryHadronFeynmanScaling"){//FeynmanScaling - test_weight = PHFSWeightCalc(fluxlist[inu], subVrandom); + test_weight = PHFSWeightCalc(flux, subVrandom); } else if( CalcType == "PrimaryHadronSanfordWang"){//SanfordWang - test_weight = PHSWWeightCalc(fluxlist[inu], subVrandom); + test_weight = PHSWWeightCalc(flux, subVrandom); } else if( CalcType == "PrimaryHadronSWCentralSplineVariation"){//SWCentaralSplineVariation - test_weight = PHSWCSVWeightCalc(fluxlist[inu], subVrandom); - + test_weight = PHSWCSVWeightCalc(flux, subVrandom); } else throw cet::exception(__PRETTY_FUNCTION__) << GetName() << ": this shouldnt happen.."< Date: Tue, 9 Dec 2025 14:28:33 -0600 Subject: [PATCH 2/3] Remove spurious cout --- sbncode/SBNEventWeight/Calculators/BNBFlux/FluxCalcPrep.cxx | 1 - 1 file changed, 1 deletion(-) diff --git a/sbncode/SBNEventWeight/Calculators/BNBFlux/FluxCalcPrep.cxx b/sbncode/SBNEventWeight/Calculators/BNBFlux/FluxCalcPrep.cxx index b1677ce35..7e0e9265e 100644 --- a/sbncode/SBNEventWeight/Calculators/BNBFlux/FluxCalcPrep.cxx +++ b/sbncode/SBNEventWeight/Calculators/BNBFlux/FluxCalcPrep.cxx @@ -323,7 +323,6 @@ namespace sbn { if (std::find(fprimaryHad.begin(), fprimaryHad.end(),(flux.ftptype)) == fprimaryHad.end() ){//if it does not contain any particles we need get 1 weights.resize( NUni); std::fill(weights.begin(), weights.end(), 1); - std::cout << "We don't need this parent, returning." << std::endl; return weights;//done, all 1 }// Hadronic parent check From 2a55bc48996f6361c33fbd4a167f77580bb22e73 Mon Sep 17 00:00:00 2001 From: Marco Del Tutto Date: Tue, 16 Dec 2025 09:48:13 -0600 Subject: [PATCH 3/3] Use BooNEpBeInteraction to identify pBe interaction, also better loop Co-authored-by: John Plows kom.plows@gmail.com --- .../Calculators/BNBFlux/FluxCalcPrep.cxx | 22 ++++++++----------- 1 file changed, 9 insertions(+), 13 deletions(-) diff --git a/sbncode/SBNEventWeight/Calculators/BNBFlux/FluxCalcPrep.cxx b/sbncode/SBNEventWeight/Calculators/BNBFlux/FluxCalcPrep.cxx index 7e0e9265e..ff461c562 100644 --- a/sbncode/SBNEventWeight/Calculators/BNBFlux/FluxCalcPrep.cxx +++ b/sbncode/SBNEventWeight/Calculators/BNBFlux/FluxCalcPrep.cxx @@ -302,24 +302,20 @@ namespace sbn { flux.ftpy = fluxlist[inu].ftpy; flux.ftpz = fluxlist[inu].ftpz; - // int tptype = fluxlist[inu].ftptype; - // If Dk2Nu flux, use ancestors to evaluate tptype if (fluxlist[inu].fFluxType == simb::kDk2Nu) { - // Find the first inelastic interaction - int firstInelastic = 0; - while (dk2nu_v->at(inu).ancestor[firstInelastic].proc.find("HadronInelastic")==std::string::npos||dk2nu_v->at(inu).ancestor[firstInelastic].proc.find("QEBooNE")!=std::string::npos) firstInelastic++; - - // Get the parent/ancestor (this should be the secondary in the p+Be interaction) - flux.ftptype = dk2nu_v->at(inu).ancestor[firstInelastic].pdg; - flux.ftpx = dk2nu_v->at(inu).ancestor[firstInelastic].startpx; - flux.ftpy = dk2nu_v->at(inu).ancestor[firstInelastic].startpy; - flux.ftpz = dk2nu_v->at(inu).ancestor[firstInelastic].startpz; - + for( const bsim::Ancestor & ancestor : dk2nu_v->at(inu).ancestor ) { + std::string aproc = ancestor.proc; + if( (aproc.find("BooNEHadronInelastic:BooNEpBeInteraction") != std::string::npos) && (aproc.find("QEBooNE") == std::string::npos) ) { + flux.ftptype = ancestor.pdg; + flux.ftpx = ancestor.startpx; + flux.ftpy = ancestor.startpy; + flux.ftpz = ancestor.startpz; + } // found it + } // find first inelastic } - if (std::find(fprimaryHad.begin(), fprimaryHad.end(),(flux.ftptype)) == fprimaryHad.end() ){//if it does not contain any particles we need get 1 weights.resize( NUni); std::fill(weights.begin(), weights.end(), 1);