diff --git a/sbncode/EventGenerator/MeVPrtl/Tools/Constants.cpp b/sbncode/EventGenerator/MeVPrtl/Tools/Constants.cpp index 39d23b13d..2746dc4bf 100644 --- a/sbncode/EventGenerator/MeVPrtl/Tools/Constants.cpp +++ b/sbncode/EventGenerator/MeVPrtl/Tools/Constants.cpp @@ -22,7 +22,7 @@ Constants::Constants() { eta_mass = 0.547862; // GeV rho_mass = 0.77526; // GeV https://pdg.lbl.gov/2019/listings/rpp2019-list-rho-770.pdf etap_mass = 0.95778; // GeV - + gamma_mass = 0; // Couplings fine_structure_constant = 7.2973525693e-3;// https://pdg.lbl.gov/2019/reviews/rpp2019-rev-phys-constants.pdf Gfermi = 1.166379e-5; // 1/GeV^2 https://pdg.lbl.gov/2020/reviews/rpp2020-rev-phys-constants.pdf @@ -269,6 +269,8 @@ double PDG2Mass(int pdg) { return Constants::Instance().klong_mass; case 221: return Constants::Instance().eta_mass; + case 22: + return Constants::Instance().gamma_mass; case 11: case -11: return Constants::Instance().elec_mass; @@ -280,6 +282,8 @@ double PDG2Mass(int pdg) { case 15: case -15: return Constants::Instance().tau_mass; + case 111: + return Constants::Instance().pizero_mass; case 211: case -211: return Constants::Instance().piplus_mass; diff --git a/sbncode/EventGenerator/MeVPrtl/Tools/Constants.h b/sbncode/EventGenerator/MeVPrtl/Tools/Constants.h index bcce6caa7..57421f71d 100644 --- a/sbncode/EventGenerator/MeVPrtl/Tools/Constants.h +++ b/sbncode/EventGenerator/MeVPrtl/Tools/Constants.h @@ -21,6 +21,7 @@ class Constants { public: // Masses + double gamma_mass; double elec_mass; double muon_mass; double piplus_mass; diff --git a/sbncode/FluxReader/CMakeLists.txt b/sbncode/FluxReader/CMakeLists.txt index ad70d3c22..8e1eb9098 100644 --- a/sbncode/FluxReader/CMakeLists.txt +++ b/sbncode/FluxReader/CMakeLists.txt @@ -74,6 +74,7 @@ cet_make_library( cet_build_plugin(FluxReaderAna art::module LIBRARIES sbncode_FluxReader sbncode_FluxReader_FluxInterface + sbnobj::Common_SBNEventWeight nusimdata::SimulationBase larcoreobj::SummaryData larcore::Geometry_Geometry_service diff --git a/sbncode/FluxReader/FluxReaderAna_module.cc b/sbncode/FluxReader/FluxReaderAna_module.cc index 914c71da7..11080a1a3 100644 --- a/sbncode/FluxReader/FluxReaderAna_module.cc +++ b/sbncode/FluxReader/FluxReaderAna_module.cc @@ -35,6 +35,20 @@ #include "nusimdata/SimulationBase/MCFlux.h" #include "nusimdata/SimulationBase/MCTruth.h" + +#include "canvas/Persistency/Common/FindManyP.h" +#include "canvas/Persistency/Common/Ptr.h" +#include "canvas/Persistency/Common/PtrVector.h" +#include "sbnobj/Common/SBNEventWeight/EventWeightMap.h" + + +#include +#include +#include +#include +#include +#include + class FluxReaderAna; @@ -101,6 +115,19 @@ class FluxReaderAna : public art::EDAnalyzer { float _nu_other_y; /// Y poisition of neutrino at the front face of the TPC (other) float _nu_other_z; /// Z poisition of neutrino at the front face of the TPC (other) float _nu_other_t; /// Time of the neutrino (other) + + + //For weighting: + std::string _flux_eventweight_multisim_producer; + bool _save_systs_flux; + + int _evtwgt_flux_nfunc; /// Number of flux weight functions + std::vector _evtwgt_flux_funcname; /// Names of weight functions + std::vector _evtwgt_flux_nweight; /// Number of weights per function + std::vector> _evtwgt_flux_weight; /// Raw weights by function + std::vector _evtwgt_flux_oneweight; /// Combined flux multisim weight + + }; @@ -122,6 +149,10 @@ FluxReaderAna::FluxReaderAna(fhicl::ParameterSet const& p) _x_cut = p.get("XCut"); // cm _y_cut = p.get("YCut"); // cm + _flux_eventweight_multisim_producer = + p.get("FluxEventWeightProducer", "fluxweight"); + + _save_systs_flux = p.get("SaveSystsFlux", false); art::ServiceHandle fs; _tree = fs->make("tree", ""); @@ -154,6 +185,14 @@ FluxReaderAna::FluxReaderAna(fhicl::ParameterSet const& p) _tree->Branch("nu_other_y", &_nu_other_y, "nu_other_y/F"); _tree->Branch("nu_other_z", &_nu_other_z, "nu_other_z/F"); _tree->Branch("nu_other_t", &_nu_other_t, "nu_other_t/F"); + + + _tree->Branch("evtwgt_flux_nfunc", &_evtwgt_flux_nfunc, "evtwgt_flux_nfunc/I"); + _tree->Branch("evtwgt_flux_funcname", "std::vector", &_evtwgt_flux_funcname); + _tree->Branch("evtwgt_flux_nweight", "std::vector", &_evtwgt_flux_nweight); + _tree->Branch("evtwgt_flux_weight", "std::vector>", &_evtwgt_flux_weight); + _tree->Branch("evtwgt_flux_oneweight", "std::vector", &_evtwgt_flux_oneweight); + } void FluxReaderAna::analyze(art::Event const& e) @@ -167,6 +206,12 @@ void FluxReaderAna::analyze(art::Event const& e) e.getByLabel(_flux_label, mctruthHandle); std::vector const& mclist = *mctruthHandle; + std::vector> mct_v; + art::fill_ptr_vector(mct_v, mctruthHandle); + + art::FindManyP mct_to_fluxewm(mctruthHandle, e, + _flux_eventweight_multisim_producer); + for(unsigned int inu = 0; inu < mclist.size(); inu++) { simb::MCParticle nu = mclist[inu].GetNeutrino().Nu(); simb::MCFlux flux = fluxlist[inu]; @@ -228,6 +273,76 @@ void FluxReaderAna::analyze(art::Event const& e) } } + // Reset flux systematic containers every neutrino + _evtwgt_flux_nfunc = 0; + _evtwgt_flux_funcname.clear(); + _evtwgt_flux_nweight.clear(); + _evtwgt_flux_weight.clear(); + _evtwgt_flux_oneweight.clear(); + + if (_save_systs_flux) { + + std::vector> flux_ewm_v = mct_to_fluxewm.at(inu); + + if (flux_ewm_v.size() == 1) { + + std::map> evtwgt_map = *(flux_ewm_v[0]); + + std::vector previous_weights; + std::vector final_weights; + + int countFunc = 0; + + for (auto const& it : evtwgt_map) { + std::string const& func_name = it.first; + std::vector const& weight_v = it.second; + + if (previous_weights.empty()) { + previous_weights.resize(weight_v.size(), 1.f); + final_weights.resize(weight_v.size(), 1.f); + } + + // Protect against mismatched universe counts + if (weight_v.size() != previous_weights.size()) { + mf::LogWarning("FluxReaderAna") + << "Weight function " << func_name + << " has size " << weight_v.size() + << " but previous combined weight size is " + << previous_weights.size() + << ". Skipping this function."; + continue; + } + + _evtwgt_flux_funcname.push_back(func_name); + _evtwgt_flux_weight.push_back(weight_v); + _evtwgt_flux_nweight.push_back((int)weight_v.size()); + countFunc++; + + std::transform(previous_weights.begin(), previous_weights.end(), + weight_v.begin(), + final_weights.begin(), + std::multiplies()); + + previous_weights = final_weights; + } + + _evtwgt_flux_nfunc = countFunc; + _evtwgt_flux_oneweight = final_weights; + } + else if (flux_ewm_v.empty()) { + mf::LogWarning("FluxReaderAna") + << "No associated EventWeightMap found for neutrino index " + << inu << " from producer " << _flux_eventweight_multisim_producer; + } + else { + mf::LogWarning("FluxReaderAna") + << "Found " << flux_ewm_v.size() + << " associated EventWeightMaps for neutrino index " << inu + << " from producer " << _flux_eventweight_multisim_producer + << ", expected 1."; + } + } + _tree->Fill(); } } diff --git a/sbncode/FluxReader/job/run_fluxreader_dk2nu_sbnd.fcl b/sbncode/FluxReader/job/run_fluxreader_dk2nu_sbnd.fcl new file mode 100644 index 000000000..51c91068a --- /dev/null +++ b/sbncode/FluxReader/job/run_fluxreader_dk2nu_sbnd.fcl @@ -0,0 +1,7 @@ +# Driver fcl file for reading in BooNE files +# for SBND +#include "genie_sbnd.fcl" +#include "run_fluxreader_sbnd.fcl" + +source.inputType: "dk2nu" + diff --git a/sbncode/FluxReader/job/run_fluxreader_sbnd.fcl b/sbncode/FluxReader/job/run_fluxreader_sbnd.fcl index 646ff8953..45a272c1b 100644 --- a/sbncode/FluxReader/job/run_fluxreader_sbnd.fcl +++ b/sbncode/FluxReader/job/run_fluxreader_sbnd.fcl @@ -23,11 +23,68 @@ source: module_type: FluxReader skipEvents: 0 maxEvents: -1 - inputType: "gsimple" + inputType: "dk2nu" + dk2nuConfig: "dk2nu_sbnd_v3" nBins: 200 Elow: 0 Ehigh: 10 SelfIncrementRun: true + dk2nu_sbnd: { + rotmatrix: [ 1, 0, 0, + 0, 1, 0, + 0, 0, 1] + userbeam: [ 73.78, 0, 11000] + + #x_detAV: [-200, 200] + #y_detAV: [-200, 200] + #z_detAV: [0, 500] + + windowBase: [ 500, 500, -1000 ] + window1: [ 500, -500, -1000 ] + window2: [ -500, 500, -1000 ] + + } + dk2nu_sbnd_v2: { + rotmatrix: [ 1, 0, 0, + 0, 1, 0, + 0, 0, 1] + userbeam: [ 0, 0, 0, 73.78, 0, 11000] + + + windowBase: [ 500, 500, -1000 ] + window1: [ 0, -1000, 0 ] + window2: [ -1000, 0, 0 ] + + } +dk2nu_sbnd_v3: { + userbeam: [ 73.78, 0, 11000 ] + rotmatrix: [ 1, 0, 0, + 0, 1, 0, + 0, 0, 1 ] + windowBase: [ 500, 500, -1000 ] + window1: [ -500, 500, -1000 ] + window2: [ 500, -500, -1000 ] + } + +dk2nu_sbnd_leo: { + userbeam: [ -73.78, 0, 11000 ] + rotmatrix: [ 1, 0, 0, + 0, 1, 0, + 0, 0, 1 ] + windowBase: [ 500, 500, -1000 ] + window1: [ -500, 500, -1000 ] + window2: [ 500, -500, -1000 ] + } + dk2nu_bnb_at_uboone_v1: { + userbeam: [ -130., 0., 47000 ] + #userbeam: [ -73.78, 0., -11000.] + rotmatrix: [ 1, 0, 0, + 0, 1, 0, + 0, 0, 1 ] + windowBase: [ 630, 500, -1000 ] + window1: [ 630, -500, -1000 ] + window2: [ -370, 500, -1000 ] + } } outputs: @@ -38,7 +95,7 @@ outputs: fileName: "fluxreader.root" compressionLevel: 1 dataTier: "simulated" - SelectEvents: [ filter ] + SelectEvents: [ filter ] } } diff --git a/sbncode/FluxReader/job/run_fluxreader_sbnd_SBNEventWeight.fcl b/sbncode/FluxReader/job/run_fluxreader_sbnd_SBNEventWeight.fcl index 37949eb6f..89269228d 100644 --- a/sbncode/FluxReader/job/run_fluxreader_sbnd_SBNEventWeight.fcl +++ b/sbncode/FluxReader/job/run_fluxreader_sbnd_SBNEventWeight.fcl @@ -38,6 +38,7 @@ physics.producers.fluxweight.weight_functions_flux: [horncurrent , kzero , piplus , piminus + , fluxhist ] # Need to overwrite the parameter weight_functions as weight_functions_flux; diff --git a/sbncode/FluxReader/job/run_fluxreader_sbnd_boone_SBNEventWeight.fcl b/sbncode/FluxReader/job/run_fluxreader_sbnd_boone_SBNEventWeight.fcl new file mode 100644 index 000000000..f7f447bbe --- /dev/null +++ b/sbncode/FluxReader/job/run_fluxreader_sbnd_boone_SBNEventWeight.fcl @@ -0,0 +1,50 @@ +# Driver fcl file for reading in gsimple files +# for SBND +# Add SBNEventWeight -- Keng Lin June 2021 + +#include "seedservice.fcl" +#include "services_sbnd.fcl" +#include "eventweight_flux_sbn.fcl" +#include "eventweight_genie_sbn.fcl" + +#include "run_fluxreader_sbnd.fcl" + + +source.inputType: "boone" + +physics.producers.fluxweight: @local::sbn_eventweight_flux +# physics.producers.genieweight: @local::sbn_eventweight_genie + +physics.filter: [ rns + , fluxfilter + , fluxweight + # , genieweight + ] + + +# physics.filters.fluxfilter.volumes: ["volWorld"] + +physics.producers.fluxweight.generator_module_label: flux + +# Customize what kind of weights we need, i.e. define weight_functions_flux: +physics.producers.fluxweight.weight_functions_flux: [horncurrent + , expskin + , pioninexsec + , pionqexsec + , piontotxsec + , nucleoninexsec + , nucleonqexsec + , nucleontotxsec + , kplus + , kminus + , kzero + , piplus + , piminus + , fluxhist + ] + +# Need to overwrite the parameter weight_functions as weight_functions_flux; +physics.producers.fluxweight.weight_functions: @local::physics.producers.fluxweight.weight_functions_flux + + + diff --git a/sbncode/FluxReader/job/run_fluxreader_sbnd_dk2nu_SBNEventWeight.fcl b/sbncode/FluxReader/job/run_fluxreader_sbnd_dk2nu_SBNEventWeight.fcl new file mode 100644 index 000000000..521629c6f --- /dev/null +++ b/sbncode/FluxReader/job/run_fluxreader_sbnd_dk2nu_SBNEventWeight.fcl @@ -0,0 +1,49 @@ +# Driver fcl file for reading in gsimple files +# for SBND +# Add SBNEventWeight -- Keng Lin June 2021 + +#include "seedservice.fcl" +#include "services_sbnd.fcl" +#include "eventweight_flux_sbn.fcl" +#include "eventweight_genie_sbn.fcl" + +#include "run_fluxreader_sbnd.fcl" + +source.inputType: "dk2nu" + +physics.producers.fluxweight: @local::sbn_eventweight_flux +# physics.producers.genieweight: @local::sbn_eventweight_genie + +physics.filter: [ rns + , fluxfilter + , fluxweight + # , genieweight + ] + + +# physics.filters.fluxfilter.volumes: ["volWorld"] + +physics.producers.fluxweight.generator_module_label: flux + +# Customize what kind of weights we need, i.e. define weight_functions_flux: +physics.producers.fluxweight.weight_functions_flux: [horncurrent + , expskin + , pioninexsec + , pionqexsec + , piontotxsec + , nucleoninexsec + , nucleonqexsec + , nucleontotxsec + , kplus + , kminus + , kzero + , piplus + , piminus + , fluxhist + ] + +# Need to overwrite the parameter weight_functions as weight_functions_flux; +physics.producers.fluxweight.weight_functions: @local::physics.producers.fluxweight.weight_functions_flux + + + diff --git a/sbncode/SBNEventWeight/Calculators/BNBFlux/CMakeLists.txt b/sbncode/SBNEventWeight/Calculators/BNBFlux/CMakeLists.txt index 37733ab95..ff8e10d4d 100644 --- a/sbncode/SBNEventWeight/Calculators/BNBFlux/CMakeLists.txt +++ b/sbncode/SBNEventWeight/Calculators/BNBFlux/CMakeLists.txt @@ -13,6 +13,7 @@ art_make_library( fhiclcpp::fhiclcpp cetlib::cetlib nusimdata::SimulationBase + larsim::EventWeight_Base ROOT::Hist ) diff --git a/sbncode/SBNEventWeight/Calculators/BNBFlux/FluxCalcPrep.cxx b/sbncode/SBNEventWeight/Calculators/BNBFlux/FluxCalcPrep.cxx index d832e8638..713e98d20 100644 --- a/sbncode/SBNEventWeight/Calculators/BNBFlux/FluxCalcPrep.cxx +++ b/sbncode/SBNEventWeight/Calculators/BNBFlux/FluxCalcPrep.cxx @@ -9,7 +9,7 @@ namespace sbn { //Configure everything here! void FluxWeightCalc::Configure(fhicl::ParameterSet const& p, - CLHEP::HepRandomEngine& engine) { + CLHEP::HepRandomEngine& engine) { std::cout<<"SBNEventWeight Flux: configure "<< GetName() << std::endl; fGeneratorModuleLabel = p.get("generator_module_label");//use this label to get MC*Handle @@ -22,12 +22,12 @@ namespace sbn { std::vector< float > parsigmas(pars.size(), 1.0); if (pars.size() != parsigmas.size()) { throw cet::exception(__PRETTY_FUNCTION__) << GetName() << ": " - << "parameter_list and parameter_sigma length mismatch." - << std::endl; + << "parameter_list and parameter_sigma length mismatch." + << std::endl; } if(!pset.get_if_present("parameter_sigma", parsigmas)){ - std::cout<<"SBNEventWeight Flux: `parameter_sigma` was not set; now it is set as 1"<("mode");//3 types: multisim/pmNsigma/fixed @@ -44,10 +44,13 @@ namespace sbn { CalcType = pset.get< std::string >("calc_type");//Unisim,PrimaryHadronSWCentralSplineVariation,PrimaryHadronFeynmanScaling,PrimaryHadronSanfordWang,PrimaryHadronNormalization std::cout<<"SBNEventWeight Flux: Calculator type: "<("use_flux_hist", false); + fScalePos = pset.get("scale_factor_pos"); if(!pset.get_if_present("scale_factor_neg",fScaleNeg)){ - std::cout<<"SBNEventWeight Flux: auto-assignment: scale_factor_neg = 1."<("CentralValue_hist_file"); - std::string cvfile_path = sp.find_file(dataInput1); + std::string dataInput1 = pset.get< std::string >("CentralValue_hist_file"); + std::string cvfile_path = sp.find_file(dataInput1); TFile fcv(Form("%s",cvfile_path.c_str())); - std::string dataInput2pos = pset.get< std::string >("PositiveSystematicVariation_hist_file"); - std::string rwfilepos = sp.find_file(dataInput2pos); + std::string dataInput2pos = pset.get< std::string >("PositiveSystematicVariation_hist_file"); + std::string rwfilepos = sp.find_file(dataInput2pos); TFile frwpos(Form("%s", rwfilepos.c_str())); - std::string dataInput2neg = pset.get< std::string >("NegativeSystematicVariation_hist_file"); - std::string rwfileneg = sp.find_file(dataInput2neg); + std::string dataInput2neg = pset.get< std::string >("NegativeSystematicVariation_hist_file"); + std::string rwfileneg = sp.find_file(dataInput2neg); TFile frwneg(Form("%s", rwfileneg.c_str())); if(dataInput2pos == dataInput2neg) PosOnly = true;//true - for skin depth @@ -89,12 +92,151 @@ namespace sbn { }//type of hadron parent fcv.Close(); frwpos.Close(); - frwneg.Close(); + frwneg.Close(); + + } + + // G4BNB to BooNE + + + else if (CalcType == "FluxHist") { + + // For debugging you can keep hardcoded paths for now. + // Later you can switch back to FHiCL-controlled file names. + // std::string cv_file = + //"/exp/sbnd/app/users/rohanr/larsoft_v10_14_02/srcs/sbncode/sbncode/SBNEventWeight/jobs/flux/nuMom_immParent_merged.root"; + // "/pnfs/sbnd/scratch/users/rohanr/beamweight/final_hists/nuMom_immParent_merged_200.root"; + // std::string rw_file = + //"/exp/sbnd/app/users/rohanr/larsoft_v10_14_02/srcs/sbncode/sbncode/SBNEventWeight/jobs/flux/nuMom_immParent_oldFiles_correctedName.root"; + // "/pnfs/sbnd/scratch/users/rohanr/beamweight/final_hists/nuMom_immParent_oldFiles_200_merged.root"; + + std::string dataInputCV = pset.get("cv_hist_file"); + std::string dataInputRW = pset.get("rw_hist_file"); + + cet::search_path sp("FW_SEARCH_PATH"); + + std::string cv_file = sp.find_file(dataInputCV); + std::string rw_file = sp.find_file(dataInputRW); + + TFile fcv(cv_file.c_str()); + TFile frw(rw_file.c_str()); + + if (fcv.IsZombie() || frw.IsZombie()) { + throw cet::exception("FluxHist") + << "Could not open FluxHist ROOT file(s):\n" + << " CV = " << cv_file << "\n" + << " RW = " << rw_file << "\n"; + } + + std::vector> validCombinations = { + { 211, "numu" }, + { 211, "nue" }, + {-211, "numubar" }, + {-211, "nuebar" }, + { 321, "numu" }, + { 321, "nue" }, + {-321, "numubar" }, + {-321, "nuebar" }, + { 130, "numu" }, + { 130, "numubar" }, + { 130, "nue" }, + { 130, "nuebar" }, + { 13, "numu" }, + { 13, "nuebar" }, + { -13, "numubar" }, + { -13, "nue" } + }; + + for (size_t i = 0; i < validCombinations.size(); ++i) { + int pdg = std::get<0>(validCombinations[i]); + std::string nut = std::get<1>(validCombinations[i]); + + std::string hname = Form("hSBND_%d_%s_3mom", pdg, nut.c_str()); + + TH3D* hcv = dynamic_cast(fcv.Get(hname.c_str())); + TH3D* hrw = dynamic_cast(frw.Get(hname.c_str())); + + if (!hcv || !hrw) { + throw cet::exception("FluxHist") + << "Missing TH3D histogram: " << hname << "\n"; + } + + fCVHist3D[i] = dynamic_cast(hcv->Clone(Form("%s_cv_clone", hname.c_str()))); + fRWHist3D[i] = dynamic_cast(hrw->Clone(Form("%s_rw_clone", hname.c_str()))); + + if (!fCVHist3D[i] || !fRWHist3D[i]) { + throw cet::exception("FluxHist") + << "Failed to clone TH3D histogram: " << hname << "\n"; + } + + fCVHist3D[i]->SetDirectory(nullptr); + fRWHist3D[i]->SetDirectory(nullptr); + + std::cout << "Loaded 3D histogram: " << hname + << " CV entries=" << fCVHist3D[i]->GetEntries() + << " RW entries=" << fRWHist3D[i]->GetEntries() + << std::endl; + } + + /* + std::string compNames[3] = {"px", "py", "pz"}; + + for (size_t i = 0; i < validCombinations.size(); ++i) { + int pdg = std::get<0>(validCombinations[i]); + std::string nut = std::get<1>(validCombinations[i]); + + for (int icomp = 0; icomp < 3; ++icomp) { + std::string hname = Form("hSBND_%d_%s_%s", pdg, nut.c_str(), compNames[icomp].c_str()); + + TH1* hcv = dynamic_cast(fcv.Get(hname.c_str())); + TH1* hrw = dynamic_cast(frw.Get(hname.c_str())); + + if (!hcv || !hrw) { + throw cet::exception("FluxHist") + << "Missing histogram: " << hname << "\n"; + } + + // std::cout << "Loaded histogram: " << hname + << " CV entries=" << hcv->GetEntries() + << " RW entries=" << hrw->GetEntries() + << std::endl; + // + if (hcv->GetNbinsX() < 200 || hrw->GetNbinsX() < 200) { + throw cet::exception("FluxHist") + << "Histogram " << hname + << " has fewer than 200 bins.\n"; + } + + fCVHist[i][icomp] = + dynamic_cast(hcv->Clone(Form("%s_cv_clone", hname.c_str()))); + fRWHist[i][icomp] = + dynamic_cast(hrw->Clone(Form("%s_rw_clone", hname.c_str()))); + + if (!fCVHist[i][icomp] || !fRWHist[i][icomp]) { + throw cet::exception("FluxHist") + << "Failed to clone histogram: " << hname << "\n"; + } + + fCVHist[i][icomp]->SetDirectory(nullptr); + fRWHist[i][icomp]->SetDirectory(nullptr); + + std::cout << "Loaded histogram: " << hname + << " CV entries=" << fCVHist[i][icomp]->GetEntries() + << " RW entries=" << fRWHist[i][icomp]->GetEntries() + << std::endl; + + } + } +*/ + fcv.Close(); + frw.Close(); + } + //------------- //-- Hadrons -- //------------- - } else if( CalcType.compare(0, 13,"PrimaryHadron") == 0){//Hadron Calculators + else if( CalcType.compare(0, 13,"PrimaryHadron") == 0){//Hadron Calculators fprimaryHad = pset.get< std::vector>("PrimaryHadronGeantCode"); @@ -104,11 +246,11 @@ namespace sbn { } else{//Other Hadron Calculators - std::string dataInput = pset.get< std::string >("ExternalData"); - std::string ExternalDataInput = sp.find_file(dataInput); + std::string dataInput = pset.get< std::string >("ExternalData"); + std::string ExternalDataInput = sp.find_file(dataInput); TFile* file = new TFile(Form("%s",ExternalDataInput.c_str())); - std::vector< std::string > pname; // these are keys to histograms + std::vector< std::string > pname; // these are keys to histograms if( CalcType == "PrimaryHadronFeynmanScaling" ){//k+ pname.push_back("FS/KPlus/FSKPlusFitVal"); @@ -130,9 +272,9 @@ namespace sbn { }else if( CalcType == "PrimaryHadronSWCentralSplineVariation" ){//pi+- - std::string fitInput = pset.get< std::string >("ExternalFit"); - std::string HadronName; - std::string HadronAbriviation; + std::string fitInput = pset.get< std::string >("ExternalFit"); + std::string HadronName; + std::string HadronAbriviation; if(fprimaryHad[0] == 211){ HadronName = "PiPlus"; HadronAbriviation = "PP"; @@ -175,10 +317,10 @@ namespace sbn { // //////////////// - std::string ExternalFitInput = sp.find_file(fitInput); + std::string ExternalFitInput = sp.find_file(fitInput); TFile* Fitfile = new TFile(Form("%s",ExternalFitInput.c_str())); - std::string fitname; // these are what we will extract from the file + std::string fitname; // these are what we will extract from the file fitname = Form("SW/%s/SW%sFitVal",HadronName.c_str(),HadronName.c_str()); // Sanford-Wang Fit Parameters TArrayD* SWParamArray = (TArrayD*) Fitfile->Get(fitname.c_str()); @@ -196,21 +338,21 @@ namespace sbn { }else { throw cet::exception(__PRETTY_FUNCTION__) << GetName() << ": " - <<" calculator "+CalcType + "is invalid" - < FluxWeightCalc::GetWeight(art::Event& e, size_t inu) { bool count_weights = false; -// std::cout<<"SBNEventWeight : getweight for the "< > mcFluxHandle; e.getByLabel(fGeneratorModuleLabel, mcFluxHandle); @@ -224,18 +366,18 @@ namespace sbn { } //or do the above 3 lines in one line -// auto const& mclist = *e.getValidHandle>(fGeneratorModuleLabel); + auto const& mclist = *e.getValidHandle>(fGeneratorModuleLabel); // If no neutrinos in this event, gives 0 weight; int NUni = fParameterSet.fNuniverses; std::vector weights;//( mclist.size(), 0); if(fluxlist.size() == 0){ - std::cout<<"SBNEventWeight Flux: EMPTY WEIGHTS"<second[i]; weights[i]=UnisimWeightCalc(enu, ptype, ntype, randomN , PosOnly);//AddParameter result does not work here; - if(count_weights){ - if(weights[i]<0){ - wcn++; - }else if((weights[i]-0)<1e-30){ - wc0++; - }else if(fabs(weights[i]-30) < 1e-30){ - wc30++; - } else if(fabs(weights[i]-1)<1e-30){ - wc1++; - } else { - wc++; + if(count_weights){ + if(weights[i]<0){ + wcn++; + }else if((weights[i]-0)<1e-30){ + wc0++; + }else if(fabs(weights[i]-30) < 1e-30){ + wc30++; + } else if(fabs(weights[i]-1)<1e-30){ + wc1++; + } else { + wc++; + } + } + }//Iterate through the number of universes } } - }//Iterate through the number of universes + + else if (CalcType == "FluxHist") { + + weights.resize(NUni, 1.0); + + if (inu >= fluxlist.size()) { + throw cet::exception("FluxHist") + << "inu out of range for fluxlist: inu=" << inu + << " fluxlist.size()=" << fluxlist.size() << "\n"; + } + + if (inu >= mclist.size()) { + throw cet::exception("FluxHist") + << "inu out of range for mclist: inu=" << inu + << " mclist.size()=" << mclist.size() << "\n"; + } + + int parent_pdg = fluxlist[inu].fptype; + int nu_pdg = fluxlist[inu].fntype; + + double px = mclist[inu].GetNeutrino().Nu().Px(); + double py = mclist[inu].GetNeutrino().Nu().Py(); + double pz = mclist[inu].GetNeutrino().Nu().Pz(); + + double w = FluxHistWeightCalc(parent_pdg, nu_pdg, px, py, pz); + + for (int i = 0; i < NUni; ++i) { + weights[i] = std::isfinite(w) ? w : 1.0; + } + } + + /* + else if (CalcType == "FluxHist") { + + weights.resize(NUni, 1.0); + + if (inu >= fluxlist.size()) { + throw cet::exception("FluxHist") + << "inu out of range for fluxlist: inu=" << inu + << " fluxlist.size()=" << fluxlist.size() << "\n"; } - } else{//then this must be PrimaryHadron + if (inu >= mclist.size()) { + throw cet::exception("FluxHist") + << "inu out of range for mclist: inu=" << inu + << " mclist.size()=" << mclist.size() << "\n"; + } - // 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 + int parent_pdg = fluxlist[inu].fptype; + int nu_pdg = fluxlist[inu].fntype; - simb::MCFlux flux; - flux.ftptype = fluxlist[inu].ftptype; - flux.ftpx = fluxlist[inu].ftpx; - flux.ftpy = fluxlist[inu].ftpy; - flux.ftpz = fluxlist[inu].ftpz; - - // If Dk2Nu flux, use ancestors to evaluate tptype - if (fluxlist[inu].fFluxType == simb::kDk2Nu) { - - 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 + int idx = GetFluxHistIndex(parent_pdg, nu_pdg); + if (idx < 0 || idx >= 14) { + return weights; + } + + double px = mclist[inu].GetNeutrino().Nu().Px(); + double py = mclist[inu].GetNeutrino().Nu().Py(); + double pz = mclist[inu].GetNeutrino().Nu().Pz(); + + double w = FluxHistWeightCalc(parent_pdg, nu_pdg, px, py, pz); + + // Deterministic event-level correction: same value in every universe slot + for (int i = 0; i < NUni; ++i) { + weights[i] = std::isfinite(w) ? w : 1.0; + } } - if (std::find(fprimaryHad.begin(), fprimaryHad.end(),(flux.ftptype)) == fprimaryHad.end() ){//if it does not contain any particles we need get 1 + */ + /* + else if (CalcType == "FluxHist") { + + std::cout << "Entering FluxHist GetWeight for inu = " << inu << std::endl; + + weights.resize(NUni, 1.0); + + if (inu >= fluxlist.size()) { + throw cet::exception("FluxHist") + << "inu out of range for fluxlist: inu=" << inu + << " fluxlist.size()=" << fluxlist.size() << "\n"; + } + + if (inu >= mclist.size()) { + throw cet::exception("FluxHist") + << "inu out of range for mclist: inu=" << inu + << " mclist.size()=" << mclist.size() << "\n"; + } + + int parent_pdg = fluxlist[inu].fptype; + int nu_pdg = fluxlist[inu].fntype; + + + int idx = GetFluxHistIndex(parent_pdg, nu_pdg); + + //std::cout << " FluxHist index = " << idx << std::endl; + + if (idx < 0 || idx >= 14) { + std::cout << " Unsupported parent/flavor combination, returning unit weights" + << std::endl; + return weights; + } + + double pz = mclist[inu].GetNeutrino().Nu().Pz(); + + std::cout << " neutrino pz = " << pz << std::endl; + + int bin = static_cast(pz / 0.05); + + if (bin < 0) bin = 0; + if (bin > 199) bin = 199; + + std::cout << " histogram bin = " << bin << std::endl; + + if (fParameterSet.fRWType == EventWeightParameterSet::kMultiSim) { + + auto const& randomVec = (fParameterSet.fParameterMap.begin())->second; + + std::cout << " NUni = " << NUni + << ", randomVec.size() = " << randomVec.size() << std::endl; + + for (int i = 0; i < NUni; ++i) { + + if (i >= (int)randomVec.size()) { + throw cet::exception("FluxHist") + << "Random vector too short: i=" << i + << " size=" << randomVec.size() << "\n"; + } + + // double randomN = randomVec[i]; + + double cv = fCVHist[idx][bin]; + double rw = fRWHist[idx][bin]; + + + double w = 1.0; + + if (cv > 0.0) { + double ratio = rw / cv; + if (std::isfinite(ratio)) { + w = ratio;//1.0 - (1.0 - ratio) * randomN; + } + } + + weights[i] = std::isfinite(w) ? w : 1.0; + } + } + + std::cout << "Leaving FluxHist GetWeight" << std::endl; + } + + */ + else{//then this must be PrimaryHadron + + + // 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 + + simb::MCFlux flux; + flux.ftptype = fluxlist[inu].ftptype; + flux.ftpx = fluxlist[inu].ftpx; + flux.ftpy = fluxlist[inu].ftpy; + flux.ftpz = fluxlist[inu].ftpz; + + // If Dk2Nu flux, use ancestors to evaluate tptype + if (fluxlist[inu].fFluxType == simb::kDk2Nu) { + + 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); + std::fill(weights.begin(), weights.end(), 1); return weights;//done, all 1 }// Hadronic parent check if(fParameterSet.fRWType == EventWeightParameterSet::kMultiSim){ for (unsigned int i = 0; int(weights.size()) < NUni; i++) {//if all weights are 1, no need to calculate weights; - std::pair test_weight; + std::pair test_weight; - std::vector< float > Vrandom = (fParameterSet.fParameterMap.begin())->second;//vector of random # - std::vector< float > subVrandom;//sub-vector of random numbers; + 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(flux, Vrandom[i]); @@ -348,20 +642,20 @@ namespace sbn { if(test_weight.first){ weights.push_back(test_weight.second); double tmp_weight = test_weight.second; -// std::cout<GetXaxis()->FindBin(px); + int biny = hcv->GetYaxis()->FindBin(py); + int binz = hcv->GetZaxis()->FindBin(pz); + + if (binx < 1) binx = 1; + if (binx > hcv->GetNbinsX()) binx = hcv->GetNbinsX(); + + if (biny < 1) biny = 1; + if (biny > hcv->GetNbinsY()) biny = hcv->GetNbinsY(); + + if (binz < 1) binz = 1; + if (binz > hcv->GetNbinsZ()) binz = hcv->GetNbinsZ(); + + double cv = hcv->GetBinContent(binx, biny, binz); + double rw = hrw->GetBinContent(binx, biny, binz); + + double weight = 1.0; + + if (cv > 0.0) { + weight = rw / cv; + if (!std::isfinite(weight) || weight <= 0.0) { + weight = 1.0; + } + } + + if (debug_counter < 20) { + std::cout << "\n[FluxHist 3D DEBUG]" + << " parent=" << parent_pdg + << " nu=" << nu_pdg + << "\n px=" << px << " binx=" << binx + << "\n py=" << py << " biny=" << biny + << "\n pz=" << pz << " binz=" << binz + << "\n CV=" << cv + << "\n RW=" << rw + << "\n weight=" << weight + << std::endl; + } + + debug_counter++; + + return std::isfinite(weight) ? weight : 1.0; + } + /* + double FluxWeightCalc::FluxHistWeightCalc(int parent_pdg, int nu_pdg, + double px, double py, double pz) + { + static int debug_counter = 0; + + int idx = GetFluxHistIndex(parent_pdg, nu_pdg); + if (idx < 0) return 1.0; + + double comps[3] = {px, py, pz}; + const char* compNames[3] = {"px", "py", "pz"}; + + double log_sum = 0.0; + int n_used = 0; + + if (debug_counter < 20) { + std::cout << "\n[FluxHist DEBUG] Event " << debug_counter + << " parent=" << parent_pdg + << " nu=" << nu_pdg << std::endl; + } + + for (int icomp = 0; icomp < 3; ++icomp) { + TH1* hcv = fCVHist[idx][icomp]; + TH1* hrw = fRWHist[idx][icomp]; + + if (!hcv || !hrw) continue; + + double val = comps[icomp]; + int bin = hcv->GetXaxis()->FindBin(val); + + if (bin < 1) bin = 1; + if (bin > hcv->GetNbinsX()) bin = hcv->GetNbinsX(); + + double cv = hcv->GetBinContent(bin); + double rw = hrw->GetBinContent(bin); + + double ratio = 1.0; + + if (cv > 0.0) { + ratio = rw / cv; + if (!std::isfinite(ratio) || ratio <= 0.0) ratio = 1.0; + } + + + log_sum += std::log(ratio); + n_used++; + + if (debug_counter < 20) { + std::cout << " " << compNames[icomp] + << " = " << val + << " bin=" << bin + << " CV=" << cv + << " RW=" << rw + << " ratio=" << ratio + << std::endl; + } + } + + if (n_used == 0) return 1.0; + + double weight = std::exp(log_sum / n_used); + + if (debug_counter < 20) { + std::cout << " --> FINAL WEIGHT = " << weight << std::endl; + } + + debug_counter++; + + return std::isfinite(weight) ? weight : 1.0; + } + */ + REGISTER_WEIGHTCALC(FluxWeightCalc) } // namespace evwgh } // namespace sbn - diff --git a/sbncode/SBNEventWeight/Calculators/BNBFlux/FluxCalcPrep.h b/sbncode/SBNEventWeight/Calculators/BNBFlux/FluxCalcPrep.h index 57f1d568b..45c2b3185 100644 --- a/sbncode/SBNEventWeight/Calculators/BNBFlux/FluxCalcPrep.h +++ b/sbncode/SBNEventWeight/Calculators/BNBFlux/FluxCalcPrep.h @@ -15,7 +15,8 @@ #include "sbncode/SBNEventWeight/Base/SmearingUtils.h"//MultiGaussianSmearing! //#include //for exit(0); debugging purpose - +#include "TH1.h" +#include "TH3D.h" #include "TH1F.h" #include "TFile.h" #include "TDecompChol.h"//for Choleskey Decomposition @@ -24,73 +25,117 @@ namespace sbn { namespace evwgh { class FluxWeightCalc : public WeightCalc { - public: - FluxWeightCalc() : WeightCalc() {} - - //Read FHiCL and store the settings for the reweighting environment and the calculator. - void Configure(fhicl::ParameterSet const& p, - CLHEP::HepRandomEngine& engine) override; - - - //GetWeight() returns the final weights as a vector - // each weight evaluaed by *WeightCalc() function; - //inu - the ith parameter. - std::vector GetWeight(art::Event& e, size_t inu) override; - - //UnisimWeightCalc() - Function for evaluating a specific weight - //enu - neutrino energy from simb::MCTruth; - //ptype- parent particles label: pi, k, k0, mu from simb:MCFlux - //ntype - neutrino flavor label: numu, numubar, nue, nuebar from simb:MCFlux - //randomN - input randmo number - //noNeg - determine what formulas to use for weights depending on input histograms. - - //4 *WeightCalc() functions - double UnisimWeightCalc(double enu, int ptype, int ntype, double randomN, bool noNeg);//Unisim - std::pair PHNWeightCalc (simb::MCFlux flux, float rand);//PrimaryHadronNormalizationWeightCalc - std::pair PHFSWeightCalc (simb::MCFlux flux, std::vector rand);//PrimaryHadronFeynmanScaling - std::pair PHSWWeightCalc (simb::MCFlux flux, std::vector rand);//PrimaryHadronSanfordWangWeightCalc - std::pair PHSWCSVWeightCalc(simb::MCFlux flux, std::vector rand);//PrimaryHadronSWCentralSplineVariationWeightCalc + public: + FluxWeightCalc() : WeightCalc() {} + + //Read FHiCL and store the settings for the reweighting environment and the calculator. + void Configure(fhicl::ParameterSet const& p, + CLHEP::HepRandomEngine& engine) override; + + + //GetWeight() returns the final weights as a vector + // each weight evaluaed by *WeightCalc() function; + //inu - the ith parameter. + std::vector GetWeight(art::Event& e, size_t inu) override; + + //UnisimWeightCalc() - Function for evaluating a specific weight + //enu - neutrino energy from simb::MCTruth; + //ptype- parent particles label: pi, k, k0, mu from simb:MCFlux + //ntype - neutrino flavor label: numu, numubar, nue, nuebar from simb:MCFlux + //randomN - input randmo number + //noNeg - determine what formulas to use for weights depending on input histograms. + + //4 *WeightCalc() functions + double UnisimWeightCalc(double enu, int ptype, int ntype, double randomN, bool noNeg);//Unisim + std::pair PHNWeightCalc (simb::MCFlux flux, float rand);//PrimaryHadronNormalizationWeightCalc + std::pair PHFSWeightCalc (simb::MCFlux flux, std::vector rand);//PrimaryHadronFeynmanScaling + std::pair PHSWWeightCalc (simb::MCFlux flux, std::vector rand);//PrimaryHadronSanfordWangWeightCalc + std::pair PHSWCSVWeightCalc(simb::MCFlux flux, std::vector rand);//PrimaryHadronSWCentralSplineVariationWeightCalc - //Handy tool - std::vector ConvertToVector(TArrayD const* array); - - private: - //fParameterSet was prepared in `sbncode/Base/WeightManager.h` - std::string fGeneratorModuleLabel; - std::string CalcType; - - double fScalePos{}; - double fScaleNeg = 1; //for Unisim - - //for Unisim - //load histograms values: [mu/pi/k-/k] x [nue,anue,numu,anumu] x [bin#] - double fCV[4][4][200]; - double fRWpos[4][4][200]; - double fRWneg[4][4][200]; - bool PosOnly{false}; - - //for HadronsProduction - std::vector fprimaryHad = {0}; - //-- FeynmanScaling - std::vector FitVal{};//Shared to SanfordWang - TMatrixD* FitCov{nullptr};//Shared to SanfordWang - - //-- SWCentaralSplineVariation - TMatrixD* HARPXSec{nullptr}; - std::vector HARPmomentumBounds{}; - std::vector HARPthetaBounds{}; - std::vector SWParam{}; - bool fIsDecomposed{false}; - - //Weight Counter - int wc = 0; - int wcn = 0; - int wc0 = 0; - int wc1 = 0; - int wc30 = 0; - - DECLARE_WEIGHTCALC(FluxWeightCalc) - }; + //Handy tool + std::vector ConvertToVector(TArrayD const* array); + + private: + //fParameterSet was prepared in `sbncode/Base/WeightManager.h` + std::string fGeneratorModuleLabel; + std::string CalcType; + + double fScalePos{}; + double fScaleNeg = 1; //for Unisim + + //for Unisim + //load histograms values: [mu/pi/k-/k] x [nue,anue,numu,anumu] x [bin#] + double fCV[4][4][200]; + double fRWpos[4][4][200]; + double fRWneg[4][4][200]; + bool PosOnly{false}; + + //for FluxHist + bool fUseFluxHist = false; + //double fCVHist[7][4][200]; + //double fRWHist[7][4][200]; + + //double fCVHist[14][200]{}; + //double fRWHist[14][200]{}; + + //TH1* fCVHist[14][3]{}; + //TH1* fRWHist[14][3]{}; + + TH3D* fCVHist3D[16]{}; + TH3D* fRWHist3D[16]{}; + + int GetFluxHistIndex(int parent_pdg, int nu_pdg) const { + if (parent_pdg == 211 && nu_pdg == 14) return 0; + else if (parent_pdg == 211 && nu_pdg == 12) return 1; + + else if (parent_pdg == -211 && nu_pdg == -14) return 2; + else if (parent_pdg == -211 && nu_pdg == -12) return 3; + + else if (parent_pdg == 321 && nu_pdg == 14) return 4; + else if (parent_pdg == 321 && nu_pdg == 12) return 5; + + else if (parent_pdg == -321 && nu_pdg == -14) return 6; + else if (parent_pdg == -321 && nu_pdg == -12) return 7; + + else if (parent_pdg == 130 && nu_pdg == 14) return 8; + else if (parent_pdg == 130 && nu_pdg == 12) return 9; + else if (parent_pdg == 130 && nu_pdg == -14) return 10; + else if (parent_pdg == 130 && nu_pdg == -12) return 11; + + else if (parent_pdg == 13 && nu_pdg == 14) return 12; + else if (parent_pdg == 13 && nu_pdg == -12) return 13; + + else if (parent_pdg == -13 && nu_pdg == -14) return 14; + else if (parent_pdg == -13 && nu_pdg == 12) return 15; + + else return -1; + } + + double FluxHistWeightCalc(int parent_pdg, int nu_pdg, + double px, double py, double pz); + + //for HadronsProduction + std::vector fprimaryHad = {0}; + //-- FeynmanScaling + std::vector FitVal{};//Shared to SanfordWang + TMatrixD* FitCov{nullptr};//Shared to SanfordWang + + //-- SWCentaralSplineVariation + TMatrixD* HARPXSec{nullptr}; + std::vector HARPmomentumBounds{}; + std::vector HARPthetaBounds{}; + std::vector SWParam{}; + bool fIsDecomposed{false}; + + //Weight Counter + int wc = 0; + int wcn = 0; + int wc0 = 0; + int wc1 = 0; + int wc30 = 0; + + DECLARE_WEIGHTCALC(FluxWeightCalc) + }; } // namespace evwgh } // namespace sbn diff --git a/sbncode/SBNEventWeight/Calculators/BNBFlux/FluxHistWeightCalc.cxx b/sbncode/SBNEventWeight/Calculators/BNBFlux/FluxHistWeightCalc.cxx new file mode 100644 index 000000000..f783b7ad8 --- /dev/null +++ b/sbncode/SBNEventWeight/Calculators/BNBFlux/FluxHistWeightCalc.cxx @@ -0,0 +1,241 @@ +#include "larsim/EventWeight/Base/WeightCalcCreator.h" +#include "larsim/EventWeight/Base/WeightCalc.h" + +#include + +#include "art/Framework/Core/ModuleMacros.h" +#include "art/Framework/Core/EDProducer.h" +#include "nurandom/RandomUtils/NuRandomService.h" + +#include "CLHEP/Random/RandGaussQ.h" + +#include "nusimdata/SimulationBase/MCFlux.h" +#include "nusimdata/SimulationBase/MCTruth.h" + +#include "TFile.h" +#include "TH1F.h" +#include "TH1.h" + +#include +#include + +namespace evwgh { + class FluxHistWeightCalc : public WeightCalc + { + public: + FluxHistWeightCalc() = default; + void Configure(fhicl::ParameterSet const& pset, + CLHEP::HepRandomEngine& engine); + std::vector > GetWeight(art::Event & e); + + private: + std::vector fWeightArray{}; + int fNmultisims{}; + std::string fMode{}; + std::string fGenieModuleLabel{}; + + int GetFluxHistIndex(int parent_pdg, int nu_pdg) const { + if (parent_pdg == 211 && nu_pdg == 14) return 0; // pi+ -> numu + else if (parent_pdg == -211 && nu_pdg == -14) return 1; // pi- -> numubar + else if (parent_pdg == 321 && nu_pdg == 14) return 2; // k+ -> numu + else if (parent_pdg == -321 && nu_pdg == -14) return 3; // k- -> numubar + else if (parent_pdg == 130 && nu_pdg == 12) return 4; // k0 -> nue + else if (parent_pdg == 13 && nu_pdg == 14) return 5; // mu- -> numu + else if (parent_pdg == -13 && nu_pdg == -14) return 6; // mu+ -> numubar + else return -1; + } + + + // pi+-,k+-,k0,mu+- + // | numu, numubar, nue, nuebar + // | | 50MeV bins + // | | | + // double fCV[7][4][200]; + // double fRW[7][4][200]; + + double fCV[14][200]{}; + double fRW[14][200]{}; + + + DECLARE_WEIGHTCALC(FluxHistWeightCalc) + }; + + void FluxHistWeightCalc::Configure(fhicl::ParameterSet const& p, + CLHEP::HepRandomEngine& engine) + { + //global config + fGenieModuleLabel= p.get< std::string > ("genie_module_label"); + + fhicl::ParameterSet const &pset=p.get (GetName()); + //calc config + fNmultisims = pset.get("number_of_multisims"); + fMode = pset.get("mode"); + std::string dataInput1 = pset.get< std::string >("cv_hist_file"); + std::string dataInput2 = pset.get< std::string >("rw_hist_file"); + + cet::search_path sp("FW_SEARCH_PATH"); + std::string cvfile = sp.find_file(dataInput1); + std::string rwfile = sp.find_file(dataInput2); + + std::string ptype[] = {"pi+", "pi-", "k+", "k-", "k0", "mu-", "mu+"}; + std::string ntype[] = {"numu", "numubar", "nue", "nuebar"}; + + std::vector> validCombinations = { + { 211, "numu" }, + { 211, "nue" }, + {-211, "numubar" }, + {-211, "nuebar" }, + { 321, "numu" }, + { 321, "nue" }, + {-321, "numubar" }, + {-321, "nuebar" }, + { 130, "nue" }, + { 130, "nuebar" }, + { 13, "numu" }, + { 13, "nuebar" }, + { -13, "numubar" }, + { -13, "nue" } + }; + + TFile fcv(Form("%s",cvfile.c_str())); + TFile frw(Form("%s",rwfile.c_str())); + /* for (int iptyp=0;iptyp<7;iptyp++) { + for (int intyp=0;intyp<4;intyp++) { + for (int ibin=0;ibin<200;ibin++) { + fCV[iptyp][intyp][ibin]=(dynamic_cast (fcv.Get(Form("h_%s_%s",ptype[iptyp].c_str(),ntype[intyp].c_str()))))->GetBinContent(ibin+1); + fRW[iptyp][intyp][ibin]=(dynamic_cast (frw.Get(Form("h_%s_%s",ptype[iptyp].c_str(),ntype[intyp].c_str()))))->GetBinContent(ibin+1); + } + } + }*/ + + for (size_t i = 0; i < validCombinations.size(); ++i) { + int pdg = std::get<0>(validCombinations[i]); + std::string nut = std::get<1>(validCombinations[i]); + + std::string hname = Form("hSBND_%d_%s_pz", pdg, nut.c_str()); + + TH1* hcv = dynamic_cast(fcv.Get(hname.c_str())); + TH1* hrw = dynamic_cast(frw.Get(hname.c_str())); + + if (!hcv || !hrw) { + throw cet::exception("FluxHistWeightCalc") + << "Missing histogram: " << hname << "\n"; + } + + if (hcv->GetNbinsX() < 200 || hrw->GetNbinsX() < 200) { + throw cet::exception("FluxHistWeightCalc") + << "Histogram " << hname << " has fewer than 200 bins.\n"; + } + + for (int ibin = 0; ibin < 200; ++ibin) { + fCV[i][ibin] = hcv->GetBinContent(ibin + 1); + fRW[i][ibin] = hrw->GetBinContent(ibin + 1); + } + } + + fcv.Close(); + frw.Close(); + + fWeightArray.resize(fNmultisims); + + if (fMode.find("multisim") != std::string::npos ) + for (double& weight : fWeightArray) weight = CLHEP::RandGaussQ::shoot(&engine, 0, 1.); + else + for (double& weight : fWeightArray) weight = 1.; + } + + std::vector > FluxHistWeightCalc::GetWeight(art::Event & e) + { + //calculate weight(s) here + std::vector > weight; + + // * MC flux information + auto mcfluxListHandle = e.getHandle< std::vector >(fGenieModuleLabel); + if (!mcfluxListHandle) { + return weight; + } + + // * MC truth information + auto mctruthListHandle = e.getHandle< std::vector >(fGenieModuleLabel); + if (!mctruthListHandle) { + return weight; + } + + std::vector const& fluxlist = *mcfluxListHandle; + std::vector const& mclist = *mctruthListHandle; + + weight.resize(mclist.size()); + for (unsigned int inu=0;inu(pz / 0.05); + if (bin < 0) bin = 0; + if (bin > 199) bin = 199; + + for (int i = 0; i < fNmultisims; ++i) { + double cv = fCV[idx][bin]; + double rw = fRW[idx][bin]; + + double test = 1.0; + + if (cv > 0.0) { + double ratio = rw / cv; + if (std::isfinite(ratio)) { + test = 1.0 - (1.0 - ratio) * fWeightArray[i]; + } + } + + weight[inu][i] = std::isfinite(test) ? test : 1.0; + } + + /* + int ptype=-9999; + int ntype=-9999; + int bin=-9999; + + if ( fluxlist[inu].fptype==211 ) ptype = 0; + else if ( fluxlist[inu].fptype==-211 ) ptype = 1; + else if ( fluxlist[inu].fptype==321 ) ptype = 2; + else if ( fluxlist[inu].fptype==-321 ) ptype = 3; + else if ( fluxlist[inu].fptype==130 ) ptype = 4; + else if ( fluxlist[inu].fptype==13 ) ptype = 5; + else if ( fluxlist[inu].fptype==-13 ) ptype = 6; + else { + throw cet::exception(__FUNCTION__) << GetName()<<"::Unknown ptype "<