diff --git a/.gitignore b/.gitignore index 27ad5fe5..61ec2b74 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,5 @@ bin/ lib/ -puml/ build-env.log src/make/Make.config src/make/Make.config_no_paths diff --git a/src/Apps/Makefile b/src/Apps/Makefile index 235afc3a..32c4b4d5 100644 --- a/src/Apps/Makefile +++ b/src/Apps/Makefile @@ -21,7 +21,8 @@ LIBRARIES := $(LIBRARIES) $(CERN_LIBRARIES) $(GENIE_LIBS) $(GENIE_REWEIGHT_LIBS TGT_BASE = grwght1p \ grwghtnp \ - grwproftest + grwproftest \ + grwproftestHad TGT = $(addprefix $(GENIE_REWEIGHT_BIN_PATH)/,$(TGT_BASE)) @@ -49,9 +50,16 @@ $(GENIE_REWEIGHT_BIN_PATH)/grwproftest: gRwProfTest.o $(call find_libs,grwghtnp) @echo "** Building grwghtnp" $(LD) $(LDFLAGS) gRwProfTest.o $(LIBRARIES) -o $(GENIE_REWEIGHT_BIN_PATH)/grwproftest +$(GENIE_REWEIGHT_BIN_PATH)/grwproftestHad: gRwProfTestHadronization.o $(call find_libs,grwghtnp) + @echo "** Building gRwProfTestHadronization" + $(LD) $(LDFLAGS) gRwProfTestHadronization.o $(LIBRARIES) -o $(GENIE_REWEIGHT_BIN_PATH)/grwproftestHad + +$(GENIE_REWEIGHT_BIN_PATH)/grwprofn: gRwProfNCorr.o $(call find_libs,grwghtnp) + @echo "** Building gRwProfNCorr" + $(LD) $(LDFLAGS) gRwProfNCorr.o $(LIBRARIES) -o $(GENIE_REWEIGHT_BIN_PATH)/grwprofn %.o : %.cxx - $(CXX) $(CXXFLAGS) -MMD -MP -c $(CPP_INCLUDES) $< -o $@ + $(CXX) $(CXXFLAGS) -MMD -MP -c $(CPP_INCLUDES) $< -o $@ -std=c++23 # CLEANING-UP diff --git a/src/Apps/gRwProfNCorr.cxx b/src/Apps/gRwProfNCorr.cxx new file mode 100644 index 00000000..a1b482cf --- /dev/null +++ b/src/Apps/gRwProfNCorr.cxx @@ -0,0 +1,345 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// GENIE/Generator includes +#include "Framework/EventGen/EventRecord.h" +#include "Framework/GHEP/GHepParticle.h" +#include "Framework/GHEP/GHepStatus.h" +#include "Framework/Messenger/Messenger.h" +#include "Framework/Ntuple/NtpMCEventRecord.h" +#include "Framework/Ntuple/NtpMCTreeHeader.h" +#include "Framework/Numerical/MathUtils.h" +#include "Framework/ParticleData/PDGUtils.h" +#include "Framework/Utils/AppInit.h" +#include "Framework/Utils/CmdLnArgParser.h" +#include "Framework/Utils/RunOpt.h" +#include "Framework/Utils/XSecSplineList.h" + +// GENIE/Reweight includes + +#include "Rtypes.h" +#include "RwCalculators/GReWeightDeltaradAngle.h" +#include "RwCalculators/GReWeightINukeParams.h" +#include "RwCalculators/GReWeightNuXSecNC.h" +#include "RwCalculators/GReWeightXSecEmpiricalMEC.h" +#include "RwCalculators/GReWeightXSecMEC.h" + +#include "RwCalculators/GReWeightProfessor.h" +#include "RwFramework/GReWeight.h" +#include "TArrayD.h" +#include "TAttLine.h" +#include "TF1.h" +#include "TLorentzVector.h" +#include "TMatrixDfwd.h" +#include "TPad.h" +#include "TROOT.h" +#include "TSpline.h" +#include "TVectorDfwd.h" + +#include +#include +#include +#include + +#include +#include + +using namespace genie; +using namespace genie::rew; +using std::stringstream; + +void PrintSyntax(); +void GetEventRange(Long64_t nev_in_file, Long64_t &nfirst, Long64_t &nlast); +void GetCommandLineArgs(int argc, char **argv); +void GetCorrelationMatrix(string fname, TMatrixD *&cmat); +void AdoptWeightCalcs(vector lsyst, GReWeight &rw); +bool FindIncompatibleSystematics(vector lsyst); + +// vector gOptVSyst; // not used as the syst are specifed in scanning +// file +static string gOptInpFilename; +static string gOptOutFilename; +static string gOptInpCovariance; +static Long64_t gOptNEvt1; +static Long64_t gOptNEvt2; +static int gOptRunKey = 0; +static int gOptNSyst = 0; +static int gOptNTwk = 0; +static long int gOptRanSeed; ///< random number seed + +static std::string binningxml, ipol_path; + +void GetCommandLineArgs(int argc, char **argv) { + LOG("grwghtnp", pINFO) << "*** Parsing command line arguments"; + + // Common run options. Set defaults and read. + RunOpt::Instance()->ReadFromCommandLine(argc, argv); + + CmdLnArgParser parser(argc, argv); + + // get GENIE event sample + if (parser.OptionExists('f')) { + LOG("grwghtnp", pINFO) << "Reading event sample filename"; + gOptInpFilename = parser.ArgAsString('f'); + } else { + LOG("grwghtnp", pFATAL) << "Unspecified input filename - Exiting"; + PrintSyntax(); + exit(1); + } + + // output weight file + if (parser.OptionExists('o')) { + LOG("grwghtnp", pINFO) << "Reading requested output filename"; + gOptOutFilename = parser.ArgAsString('o'); + } else { + LOG("grwghtnp", pINFO) << "Setting default output filename"; + gOptOutFilename = "rwt_cov.root"; + } + + // get ROOT covariance matrix binary file + if (parser.OptionExists('c')) { + LOG("grwghtnp", pINFO) << "Reading covariance matrix filename"; + gOptInpCovariance = parser.ArgAsString('c'); + } else { + LOG("grwghtnp", pFATAL) << "Unspecified covariance filename - Exiting"; + PrintSyntax(); + exit(1); + } + + if (parser.OptionExists('n')) { + // + LOG("grwghtnp", pINFO) << "Reading number of events to analyze"; + string nev = parser.ArgAsString('n'); + if (nev.find(",") != string::npos) { + vector vecn = parser.ArgAsLongTokens('n', ","); + if (vecn.size() != 2) { + LOG("grwghtnp", pFATAL) << "Invalid syntax"; + gAbortingInErr = true; + PrintSyntax(); + exit(1); + } + // User specified a comma-separated set of values n1,n2. + // Use [n1,n2] as the event range to process. + gOptNEvt1 = vecn[0]; + gOptNEvt2 = vecn[1]; + } else { + // User specified a single number n. + // Use [0,n] as the event range to process. + gOptNEvt1 = -1; + gOptNEvt2 = parser.ArgAsLong('n'); + } + } else { + LOG("grwghtnp", pINFO) + << "Unspecified number of events to analyze - Use all"; + gOptNEvt1 = -1; + gOptNEvt2 = -1; + } + LOG("grwghtnp", pDEBUG) << "Input event range: " << gOptNEvt1 << ", " + << gOptNEvt2; + + // number of tweaks: + if (parser.OptionExists('t')) { + LOG("grwghtnp", pINFO) << "Reading number of tweaks"; + gOptNTwk = parser.ArgAsInt('t'); + + if (gOptNTwk < 1) { + LOG("grwghtnp", pFATAL) << "Must have at least 1 tweak - Exiting"; + PrintSyntax(); + exit(1); + } + LOG("grwghtnp", pINFO) << "Number of tweaks : " << gOptNTwk; + } else { + LOG("grwghtnp", pFATAL) << "Unspecified tweaks for parameters - Exiting"; + PrintSyntax(); + exit(1); + } + + // run key: + if (parser.OptionExists('r')) { + LOG("grwghtnp", pINFO) << "Reading run key"; + gOptRunKey = parser.ArgAsInt('r'); + + LOG("grwghtnp", pINFO) << "Run key set to " << gOptRunKey; + } + + // binning_xml + if (parser.OptionExists('b')) { + LOG("grwghtnp", pINFO) << "Reading binning XML file"; + binningxml = parser.ArgAsString('b'); + } else { + LOG("grwghtnp", pFATAL) << "Unspecified binning XML file - Exiting"; + PrintSyntax(); + exit(1); + } + + // ipol_path + if (parser.OptionExists('p')) { + LOG("grwghtnp", pINFO) << "Reading interpolation path"; + ipol_path = parser.ArgAsString('p'); + } else { + LOG("grwghtnp", pFATAL) << "Unspecified interpolation path - Exiting"; + PrintSyntax(); + exit(1); + } +} + +TMatrixD cov_to_corr(const TMatrixD &cov) { + // Convert covariance matrix to correlation matrix + TMatrixD corr(cov); + for (int i = 0; i < cov.GetNrows(); ++i) { + for (int j = 0; j < cov.GetNcols(); ++j) { + if (i == j) { + corr(i, j) = 1.0; + } else { + corr(i, j) = cov(i, j) / std::sqrt(cov(i, i) * cov(j, j)); + } + } + } + return corr; +} + +std::vector cov_to_err(const TMatrixD &cov) { + // Convert covariance matrix to a vector of uncertainties + std::vector err(cov.GetNrows()); + for (int i = 0; i < cov.GetNrows(); ++i) { + err[i] = std::sqrt(cov(i, i)); + } + return err; +} + +int main(int argc, char **argv) { + GetCommandLineArgs(argc, argv); + utils::app_init::MesgThresholds(RunOpt::Instance()->MesgThresholdFiles()); + utils::app_init::RandGen(gOptRanSeed); + GHepRecord::SetPrintLevel(RunOpt::Instance()->EventRecordPrintLevel()); + + if (!RunOpt::Instance()->Tune()) { + LOG("greweight", pFATAL) << " No TuneId in RunOption"; + exit(-1); + } + RunOpt::Instance()->BuildTune(); + + // open the ROOT file and get the TTree & its header + TTree *tree = 0; + NtpMCTreeHeader *thdr = 0; + TFile file(gOptInpFilename.c_str(), "READ"); + tree = dynamic_cast(file.Get("gtree")); + thdr = dynamic_cast(file.Get("header")); + if (!tree) { + LOG("grwghtnp", pFATAL) + << "Can't find a GHEP tree in input file: " << file.GetName(); + gAbortingInErr = true; + PrintSyntax(); + exit(1); + } + LOG("grwghtnp", pNOTICE) << "Input tree header: " << *thdr; + + // LOG("grwghtnp", pNOTICE) << "Correlation matrix:"; + // cmat->Print(); + // LOG("grwghtnp", pNOTICE) << "Lower triangular matrix:"; + // lTri.Print(); + + NtpMCEventRecord *mcrec = 0; + tree->SetBranchAddress("gmcrec", &mcrec); + + Long64_t nev_in_file = tree->GetEntries(); + Long64_t nfirst = 0; + Long64_t nlast = 0; + GetEventRange(nev_in_file, nfirst, nlast); + int nev = int(nlast - nfirst + 1); + + auto observable_splines = std::make_unique(""); + observable_splines->ReadComparisonXML(binningxml, ipol_path); + + std::vector> systematics_values{}; + + // Gets Cor, which is needed in decompositions + // Assumed errors from covariance are stored in one sigma errors for + // parameters + // GetCorrelationMatrix(gOptInpCovariance, cmat); + TFile input_file = TFile(gOptInpCovariance.c_str(), "READ"); + auto cov_matrix = input_file.Get("param_covariance"); + if (!cov_matrix) { + LOG("grwghtnp", pFATAL) + << "Cannot find covariance matrix in file " << gOptInpCovariance; + exit(1); + } + auto corr_matrix = cov_to_corr(*cov_matrix); + auto std_err = cov_to_err(*cov_matrix); + + auto h_param_results = + input_file.Get("h_param_result"); // central values + if (!h_param_results) { + LOG("grwghtnp", pFATAL) + << "Cannot find h_param_result in file " << gOptInpCovariance; + exit(1); + } + + std::vector central_values{h_param_results->GetArray(), + h_param_results->GetArray() + + h_param_results->GetSize()}; + + using namespace genie::utils::math; + TMatrixD lTri = CholeskyDecomposition(corr_matrix); + // todo: read the values from tunning file + + // + + for (int itk = 0; itk < gOptNTwk; itk++) { + auto twkvals = CholeskyGenerateCorrelatedParamVariations(lTri); + assert(twkvals.GetNrows() == central_values.size()); + auto vec = + std::vector(twkvals.GetMatrixArray(), + twkvals.GetMatrixArray() + twkvals.GetNrows()); + for (int i = 0; i < vec.size(); ++i) { + // Apply the tweak to the central value + vec[i] = central_values[i] + vec[i] * std_err[i]; + } + systematics_values.push_back(vec); + } + + auto *wght_file = new TFile(gOptOutFilename.c_str(), "RECREATE"); + auto *wght_tree = new TTree("covrwt", "GENIE covariant reweighting tree"); + wght_tree->SetDirectory(wght_file); + TArrayD branch_weight(gOptNTwk); + wght_tree->Branch("weights", &branch_weight); + + for (auto id_event{nfirst}; id_event <= nlast; ++id_event) { + // Process each event with the corresponding systematic variations + tree->GetEntry(id_event); + EventRecord &event = *(mcrec->event); + for (size_t idx = 0; idx < systematics_values.size(); ++idx) { + const auto &sys_var = systematics_values[idx]; + // Apply each systematic variation to the event + observable_splines->SetSystematic(sys_var, central_values); + auto weight = observable_splines->CalcWeight(event); + branch_weight[idx] = weight; + } + wght_tree->Fill(); + mcrec->Clear(); + } + + wght_file->Write(); + wght_file->Close(); + delete wght_file; + + return 0; +} \ No newline at end of file diff --git a/src/Apps/gRwProfTestHadronization.cxx b/src/Apps/gRwProfTestHadronization.cxx new file mode 100644 index 00000000..c57bd237 --- /dev/null +++ b/src/Apps/gRwProfTestHadronization.cxx @@ -0,0 +1,708 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +// GENIE/Generator includes +#include "Framework/EventGen/EventRecord.h" +#include "Framework/GHEP/GHepParticle.h" +#include "Framework/GHEP/GHepStatus.h" +#include "Framework/Messenger/Messenger.h" +#include "Framework/Ntuple/NtpMCEventRecord.h" +#include "Framework/Utils/XSecSplineList.h" + +// GENIE/Reweight includes + +#include "RwCalculators/GReWeightDeltaradAngle.h" +#include "RwCalculators/GReWeightINukeParams.h" +#include "RwCalculators/GReWeightNuXSecNC.h" +#include "RwCalculators/GReWeightXSecEmpiricalMEC.h" +#include "RwCalculators/GReWeightXSecMEC.h" + +#include "RwCalculators/GReWeightProfessor.h" +#include "TAttLine.h" +#include "TF1.h" +#include "TLorentzVector.h" +#include "TPad.h" +#include "TROOT.h" +#include "TSpline.h" + +#include +#include +#include +#include + +#include +#include + +using std::ostringstream; +using std::string; + +using namespace genie; +using namespace genie::rew; + +std::vector fromdat(std::string path) { + std::vector data{}; + std::ifstream infile(path); + for (std::string line; std::getline(infile, line);) { + std::istringstream iss(line); + std::string _{}; + double a; + iss >> _ >> a; + data.push_back(a); + } + return data; +} + +template +std::unique_ptr get_object(std::string file_path, std::string obj_path) { + TFile root_file{file_path.c_str(), "READ"}; + auto objptr = static_cast(root_file.Get(obj_path.c_str())->Clone()); + assert(objptr); + return std::unique_ptr{objptr}; +} + +std::pair get_xsec(TH1 *h_rate, TGraph *spline) { + double fluxint{}; + // spline->SaveAs("wrong.root"); + TSpline3 sp("sp", spline); + TF1 func( + "spline", [&](double *x, double *) { return sp.Eval(*x); }, 0, + h_rate->GetXaxis()->GetXmax(), 0); + for (int ii = 1; ii <= h_rate->GetNbinsX(); ii++) { + double bin_c = h_rate->GetBinContent(ii); + double bin_up = h_rate->GetXaxis()->GetBinUpEdge(ii); + double bin_low = h_rate->GetXaxis()->GetBinLowEdge(ii); + double bin_width = bin_up - bin_low; + if (bin_c < 1 || func.Integral(bin_low, bin_up) == 0) { + continue; + } + fluxint += bin_c / func.Integral(bin_low, bin_up) * bin_width; + } + double event_rate = h_rate->Integral(); + return {event_rate, event_rate / fluxint}; +} + +template auto common_def(T &&df) { + return df + .Filter( + [](const NtpMCEventRecord &event) { + return event.event->Summary()->ProcInfo().IsWeakCC() && + (event.event->Summary()->ProcInfo().IsResonant() || + event.event->Summary()->ProcInfo().IsDeepInelastic()) + + ; + }, + {"gmcrec"}) + .Define("boost", + [](const NtpMCEventRecord &event) { + GHepParticle *target_nucleus_p = event.event->HitNucleon(); + + auto target_nucleus = target_nucleus_p + ? *target_nucleus_p + : *(event.event->TargetNucleus()); + return -target_nucleus.P4()->BoostVector(); + }, + {"gmcrec"}) + .Define("muonp", + [](const NtpMCEventRecord &event) { + return event.event->FinalStatePrimaryLepton()->P4()->P(); + }, + {"gmcrec"}) + .Define("muon_theta", + [](const NtpMCEventRecord &event) { + return event.event->FinalStatePrimaryLepton()->P4()->Theta(); + }, + {"gmcrec"}) + .Define("q", + [](const NtpMCEventRecord &event, const TVector3 &boost) { + // return event.event.event->Summary()->KinePtr()->Q2(); + auto &p4_lep = *(event.event->FinalStatePrimaryLepton()->P4()); + auto &nu_lep = *(event.event->Probe()->P4()); + // GHepParticle *target_nucleus_p = + // event.event.event->HitNucleon(); + + // auto target_nucleus = target_nucleus_p + // ? *target_nucleus_p + // : + // *(event.event.event->TargetNucleus()); + // auto boost_vec = -target_nucleus.P4()->BoostVector(); + auto vec = nu_lep - p4_lep; + vec.Boost(boost); + return vec; + }, + {"gmcrec", "boost"}) + .Define("q0", + [](const TLorentzVector &q) { + // return event.event.event->Summary()->KinePtr()->Q2(); + return q.T(); + }, + {"q"}) + .Define("Q2", + [](const TLorentzVector &q) { + // return event.event.event->Summary()->KinePtr()->Q2(); + return -q.M2(); + }, + {"q"}) + .Define("visE", + [](const NtpMCEventRecord &event) { + double visE{}; + auto nentries = event.event->GetEntriesFast(); + for (int i = 0; i < nentries; i++) { + auto part = event.event->Particle(i); + if (part->Status() == genie::kIStStableFinalState && + part->Charge() != 0) { + // remove mass part for p/n + visE += part->E(); + if (pdg::IsNeutron(part->Pdg()) || + pdg::IsProton(part->Pdg())) { + visE -= part->Mass(); + } + } + } + return visE; + }, + {"gmcrec"}) + .Define("nch", + [](const NtpMCEventRecord &event) { + int nch{}; + auto nentries = event.event->GetEntriesFast(); + for (int i = 0; i < nentries; i++) { + auto part = event.event->Particle(i); + if (part->Status() == genie::kIStStableFinalState && + part->Charge() != 0) { + nch++; + } + } + return (double)nch; + }, + {"gmcrec"}) + // todo: add angle / mag of p_had_system p_max_hadron + .Define("p_had_system", + [](const NtpMCEventRecord &event) { + auto np = event.event->GetEntriesFast(); + TLorentzVector p{}; + for (int i{}; i < np; i++) { + auto particle = event.event->Particle(i); + if (particle->Status() == kIStStableFinalState && + pdg::IsHadron(particle->Pdg())) { + p += *particle->P4(); + } + } + return p; + }, + {"gmcrec"}) + .Define("angle_p_had_system", + [](const TLorentzVector &v) { return v.Theta(); }, + {"p_had_system"}) + .Define("mag_p_had_system", [](const TLorentzVector &v) { return v.P(); }, + {"p_had_system"}) + .Define("p_max_hadron", + [](const NtpMCEventRecord &event) { + auto np = event.event->GetEntriesFast(); + TLorentzVector p{}; + for (int i{}; i < np; i++) { + auto particle = event.event->Particle(i); + if (particle->Status() == kIStStableFinalState && + pdg::IsHadron(particle->Pdg())) { + if (particle->P4()->P() > p.P()) { + p = *particle->P4(); + } + } + } + return p; + }, + {"gmcrec"}) + .Define("angle_p_max_hadron", + [](const TLorentzVector &v) { return v.Theta(); }, + {"p_max_hadron"}) + .Define("mag_p_max_hadron", [](const TLorentzVector &v) { return v.P(); }, + {"p_max_hadron"}) + .Define("mag_pl_max_hadron", + [](const TLorentzVector &v) { return v.Pz(); }, {"p_max_hadron"}) + .Define("pt", + [](const NtpMCEventRecord &event) { + auto np = event.event->GetEntriesFast(); + TLorentzVector p{}; + for (int i{}; i < np; i++) { + auto particle = event.event->Particle(i); + if (particle->Status() == kIStStableFinalState && + pdg::IsHadron(particle->Pdg())) { + p += *particle->P4(); + } + } + auto pt = p.Pt(); + return pt; + }, + {"gmcrec"}) + .Define("pl", + [](const NtpMCEventRecord &event) { + auto np = event.event->GetEntriesFast(); + TLorentzVector p{}; + for (int i{}; i < np; i++) { + auto particle = event.event->Particle(i); + if (particle->Status() == kIStStableFinalState && + pdg::IsHadron(particle->Pdg())) { + p += *particle->P4(); + } + } + auto pl = p.Pz(); + return pl; + }, + {"gmcrec"}) + .Define("W", + [](const NtpMCEventRecord &event) { + return event.event->Summary()->Kine().W(true); + }, + {"gmcrec"}); + ; + // for now we filter out event outside of the binning range +} + +// we don't have a proper configuration mechanism yet, so we need to +// hardcode sth. +int main(int argc, char **argv) { + // ROOT::EnableImplicitMT(4); + + TH1::AddDirectory(false); + std::string from_id = argc > 1 ? argv[1] : "0000"; + std::string to_id = argc > 2 ? argv[2] : "0001"; + std::string basedir = + argc > 3 ? argv[3] : "/var/home/yan/neutrino/had_scan_test/had_scan"; + if (basedir.back() != '/') + basedir += "/"; + + auto gevgen_filename = "numu_on_p.ghep.root"; + // auto gevgen_filename = "gntp.70000.ghep.root"; + + std::string input_from_file = + basedir + "scan/" + from_id + + "/3.04.02-routine_validation_01-hadronization/" + gevgen_filename; + std::string input_to_file = basedir + "scan/" + to_id + + "/3.04.02-routine_validation_01-hadronization/" + + gevgen_filename; + auto spline_path_from = + basedir + "scan/" + from_id + + "/3.04.02-routine_validation_01-xsec_vA/total_xsec.root"; + auto spline_path_to = + basedir + "scan/" + to_id + + "/3.04.02-routine_validation_01-xsec_vA/total_xsec.root"; + auto ipol_path = basedir + "ipol_test.dat"; + auto binningxml = basedir + "binning.xml"; + + auto input_from_tree = + ROOT::RDataFrame{"gtree", input_from_file}; // read the tree from the file + ROOT::RDF::Experimental::AddProgressBar(input_from_tree); + auto observable_splines = std::make_unique(""); + observable_splines->ReadComparisonXML(binningxml, ipol_path); + LOG("Test", pNOTICE) << "Read Professor Reweight success"; + + auto from = fromdat(basedir + "scan/" + from_id + "/params.dat"); + auto to = fromdat(basedir + "scan/" + to_id + "/params.dat"); + + observable_splines->SetSystematic(to, from); + LOG("Test", pNOTICE) << "Set Systematic success"; + + LOG("Test", pNOTICE) << "Testing RDataFrame::Define to calculate weight\n"; + auto tree_to = ROOT::RDataFrame{"gtree", input_to_file}; + ROOT::RDF::Experimental::AddProgressBar(tree_to); + + auto ds_ref = common_def(tree_to); + + auto tree_rew = + common_def(input_from_tree) + .Define( + "weight", + [&](NtpMCEventRecord &event) { + auto weight = observable_splines->CalcWeight(*(event.event)); + // if (abs(weight - 1.) > 1e-3) + // std::cout << "At lease we have " << weight + // << std::endl; + // event.Clear(); + return (weight < 0 || std::isnan(weight) || std::isinf(weight)) + ? 0. + : weight; + // return weight; + }, + {"gmcrec"}); + + auto modelmap = std::to_array>( + {{"muonp", {"muonp", ";p_{#mu};", 10, 0, 10.}}, + {"muon_theta", {"muon_theta", ";#theta_{#mu};", 20, 0., 2.}}, + {"Q2", {"Q2", ";Q^{2};", 30, 0., 6.}}, + {"q0", {"q0", ";q^{0};", 30, 0., 6.}}, + {"visE", {"visE", ";E_{vis};", 10, 0., 100.}}, + {"visE", {"visE_low", ";E_{vis};", 10, 0., 10.}}, + {"nch", {"nch", ";nch;", 19, -0.5, 18.5}}, + {"pt", {"pt", ";pt;", 11, 0., 10.}}, + {"pl", {"pl", ";pl;", 11, 0., 10.}}, + {"W", {"W", ";W;", 20, 0., 30.}}, + {"W", {"W_low", ";W;", 10, 0., 6.}}, + {"W", {"W_3low", ";W;", 10, 0., 3.}}, + {"angle_p_max_hadron", + {"angle_p_max_hadron", ";#theta_{p_{max}};", 10, 0, 1.5}}, + {"mag_p_max_hadron", {"mag_p_max_hadron", ";|p_{max}|;", 10, 0., 10.}}, + {"angle_p_had_system", + {"angle_p_had_system", ";#theta_{p_{had}};", 10, 0, 1.5}}, + {"mag_p_had_system", {"mag_p_had_system", ";|p_{had}|;", 10, 0., 10.}}, + {"mag_pl_max_hadron", + {"mag_pl_max_hadron", ";|p_{had}|;", 10, 0., 10.}}}); + + // auto + + auto merge_cut_lists = [merge_cut_func = [](auto &&...args) + -> std::function { + return [args...](NtpMCEventRecord &event) { return (args(event) && ...); }; + }](auto &&lista, auto &&listb) { + return std::views::transform( + std::views::cartesian_product(lista, listb), + [merge_cut_func](auto &&cuts) { + auto [cuta, cutb] = cuts; + auto [namea, cutfunc_a] = cuta; + auto [nameb, cutfunc_b] = cutb; + auto newcut = merge_cut_func(cutfunc_a, cutfunc_b); + return std::make_tuple(namea + nameb, newcut); + }) | + std::ranges::to(); + }; + + auto W_cuts = std::to_array< + std::tuple>>( + {{"LowW", + [](const NtpMCEventRecord &event) -> bool { + double W = event.event->Summary()->Kine().W(true); + return W < 2.3; + }}, + {"MidW", + [](const NtpMCEventRecord &event) -> bool { + double W = event.event->Summary()->Kine().W(true); + return W > 2.3 && W < 3; + }}, + {"HiW", [](const NtpMCEventRecord &event) -> bool { + double W = event.event->Summary()->Kine().W(true); + return W > 3; + }}}); + + auto to_npi = [](const NtpMCEventRecord &event) { + size_t npi{}; + auto np = event.event->GetEntriesFast(); + for (int i{}; i < np; i++) { + auto particle = event.event->Particle(i); + if (particle->Status() == kIStStableFinalState && + pdg::IsPion(particle->Pdg())) { + npi++; + } + } + return npi; + }; + + auto multi_cuts = std::to_array< + std::tuple>>( + {{"0pi", + [&](const NtpMCEventRecord &event) -> bool { + return to_npi(event) == 0; + }}, + {"1pi", + [&](const NtpMCEventRecord &event) -> bool { + return to_npi(event) == 1; + }}, + {"2pi", + [&](const NtpMCEventRecord &event) -> bool { + return to_npi(event) == 2; + }}, + {"Mpi", [&](const NtpMCEventRecord &event) -> bool { + return to_npi(event) >= 3; + }}}); + + auto channelcuts = merge_cut_lists(W_cuts, multi_cuts); + channelcuts.insert(channelcuts.end(), W_cuts.begin(), W_cuts.end()); + + // , + // {"LowW_0pi", + // [](const NtpMCEventRecord &event) -> bool { + // double W = event.event->Summary()->Kine().W(true); + // size_t npi{}; + // auto np = event.event->GetEntriesFast(); + // for (int i{}; i < np; i++) { + // auto particle = event.event->Particle(i); + // if (particle->Status() == kIStStableFinalState && + // pdg::IsPion(particle->Pdg())) { + // npi++; + // } + // } + // return W < 3 && npi == 0; + // }}, + // {"LowW_1pi", + // [](const NtpMCEventRecord &event) -> bool { + // double W = event.event->Summary()->Kine().W(true); + // size_t npi{}; + // auto np = event.event->GetEntriesFast(); + // for (int i{}; i < np; i++) { + // auto particle = event.event->Particle(i); + // if (particle->Status() == kIStStableFinalState && + // pdg::IsPion(particle->Pdg())) { + // npi++; + // } + // } + // return W < 3 && npi == 1; + // }}, + // {"LowW_2pi", + // [](const NtpMCEventRecord &event) -> bool { + // double W = event.event->Summary()->Kine().W(true); + // size_t npi{}; + // auto np = event.event->GetEntriesFast(); + // for (int i{}; i < np; i++) { + // auto particle = event.event->Particle(i); + // if (particle->Status() == kIStStableFinalState && + // pdg::IsPion(particle->Pdg())) { + // npi++; + // } + // } + // return W < 3 && npi == 2; + // }}, + // {"LowW_Mpi", [](const NtpMCEventRecord &event) -> bool { + // double W = event.event->Summary()->Kine().W(true); + // size_t npi{}; + // auto np = event.event->GetEntriesFast(); + // for (int i{}; i < np; i++) { + // auto particle = event.event->Particle(i); + // if (particle->Status() == kIStStableFinalState && + // pdg::IsPion(particle->Pdg())) { + // npi++; + // } + // } + // return W < 3 && npi >= 3; + // }}} + + auto hists_ref = + modelmap | + std::views::transform( + [&](const std::tuple &obj) { + auto [var, model] = obj; + auto plot_cut = + channelcuts | + std::views::transform( + [&](const std::tuple< + std::string, std::function> + &obj) { + auto &[cutname, cut] = obj; + auto thismodel = model; + thismodel.fName += cutname; + return ds_ref.Filter(cut, {"gmcrec"}) + .Histo1D(thismodel, var); + }) | + std::ranges::to(); + return std::make_tuple(ds_ref.Histo1D(model, var), + plot_cut); + }) | + std::ranges::to(); + + auto hists_rew = + modelmap | + std::views::transform( + [&](const std::tuple &obj) { + auto [var, model] = obj; + auto plot_cut = + channelcuts | + std::views::transform( + [&](const std::tuple< + std::string, std::function> + &obj) { + auto &[cutname, cut] = obj; + auto thismodel = model; + thismodel.fName += cutname; + return tree_rew.Filter(cut, {"gmcrec"}) + .Histo1D(thismodel, var, "weight"); + }) | + std::ranges::to(); + return std::make_tuple( + tree_rew.Histo1D(model, var, "weight"), + plot_cut); + }) | + std::ranges::to(); + + auto hists_unrew = + modelmap | + std::views::transform( + [&](const std::tuple &obj) { + auto [var, model] = obj; + auto plot_cut = + channelcuts | + std::views::transform( + [&](const std::tuple< + std::string, std::function> + &obj) { + auto &[cutname, cut] = obj; + auto thismodel = model; + thismodel.fName += cutname; + return tree_rew.Filter(cut, {"gmcrec"}) + .Histo1D(thismodel, var); + }) | + std::ranges::to(); + return std::make_tuple(tree_rew.Histo1D(model, var), + plot_cut); + }) | + std::ranges::to(); + + std::cout << "Running RDF\n\n" << std::flush; + ROOT::RDF::RunGraphs({ + std::get<0>(hists_ref[0]), + std::get<0>(hists_rew[0]), + std::get<0>(hists_unrew[0]), + }); + + auto do_plot = [&](TH1D *hist_rew, TH1D *hist_reference, TH1D *hist_unrew) { + std::string ytitle = "Count (Count * Weight)"; + + auto chi2 = hist_reference->Chi2Test(hist_rew, "CHI2/NDF"); + + auto c1 = std::make_unique(); + c1->Draw(); + + constexpr double fair_share = 0.7; + auto pad1 = std::make_unique("pad1", "", 0, 1 - fair_share, 1, 1); + auto pad2 = std::make_unique("pad2", "", 0, 0, 1, 1 - fair_share); + pad1->SetBottomMargin(0); + pad1->SetTopMargin(pad1->GetTopMargin() / 2.); + pad2->SetTopMargin(0); + pad2->SetBottomMargin(pad2->GetBottomMargin() * 2.); + + // c1->cd(1); + pad1->cd(); + auto max = std::max({hist_reference->GetMaximum(), hist_rew->GetMaximum(), + hist_unrew->GetMaximum()}); + gStyle->SetOptStat(0); + auto legend = std::make_unique(0.7, 0.7, 0.9, 0.9); + double ytitleoffset = 1.2; + hist_rew->SetLineColor(kBlue); + hist_rew->SetMaximum(max * 1.1); + hist_rew->GetYaxis()->SetTitle(ytitle.c_str()); + hist_rew->GetYaxis()->SetTitleOffset(ytitleoffset); + hist_rew->Draw("hist E"); + legend->AddEntry(hist_rew, "reweighted", "l"); + hist_unrew->SetLineColor(kGreen); + hist_unrew->SetMaximum(max * 1.1); + hist_unrew->SetLineStyle(kDashed); + // hist_unrew->SetLineStyle(kDotted); + hist_unrew->GetYaxis()->SetTitle(ytitle.c_str()); + hist_unrew->GetYaxis()->SetTitleOffset(ytitleoffset); + hist_unrew->Draw("hist E same"); + legend->AddEntry(hist_unrew, "unreweighted", "l"); + hist_reference->SetMaximum(max * 1.1); + hist_reference->SetLineColor(kRed); + // hist_reference->SetLineStyle(kDotted); + hist_reference->GetYaxis()->SetTitle(ytitle.c_str()); + hist_reference->GetYaxis()->SetTitleOffset(ytitleoffset); + hist_reference->SetLineStyle(kDashed); + hist_reference->Draw("hist E SAME"); + legend->AddEntry(hist_reference, "reference", "l"); + legend->SetHeader(("chi2/ndf: " + std::to_string(chi2)).c_str()); + legend->Draw("SAME"); + // c1->cd(2); + pad2->cd(); + auto hist_diff_rew_ref = + std::unique_ptr{static_cast(hist_rew->Clone())}; + hist_diff_rew_ref->Add(hist_reference, -1); + // hist_diff_rew_ref: Rew - Ref + + // auto hist_diff_unrew_ref = + // std::unique_ptr{static_cast(hist_unrew->Clone())}; + // hist_diff_unrew_ref->Add(hist_reference, -1); + // (hist_rew - hist_reference) / (hist_unrew - hist_reference); + hist_diff_rew_ref->Divide(hist_reference); + // hist_diff_rew_ref->SetMaximum(); + // hist_diff_rew_ref->Set + // for (int bin_id = 1; bin_id < hist_diff_rew_ref->GetNbinsX(); bin_id++) { + // double sigma_rew = hist_rew->GetBinError(bin_id); + // double sigma_unrew = hist_unrew->GetBinError(bin_id); + // double sigma_ref = hist_reference->GetBinError(bin_id); + // double rew = hist_rew->GetBinContent(bin_id); + // double unrew = hist_unrew->GetBinContent(bin_id); + // double ref = hist_reference->GetBinContent(bin_id); + // using std::sqrt, std::pow; + // double new_err = sqrt( + // (pow(1 / (ref - unrew), 2) * pow(sigma_rew, 2)) + + // (pow((rew - ref) / pow(ref - unrew, 2), 2) * pow(sigma_unrew, 2)) + + // (pow(-(rew - unrew) / pow(ref - unrew, 2), 2) * pow(sigma_ref, + // 2))); + // if (rew == 0 || unrew == 0 || ref == 0) { + // new_err = 0; + // } + // hist_diff_rew_ref->SetBinError(bin_id, new_err); + // } + auto text_size = hist_diff_rew_ref->GetYaxis()->GetLabelSize(); + auto title_offset = hist_reference->GetYaxis()->GetTitleOffset(); + + auto text_size_scaled = text_size * fair_share / (1 - fair_share); + auto title_offset_scaled = title_offset / (fair_share / (1 - fair_share)); + hist_diff_rew_ref->GetXaxis()->SetLabelSize(text_size_scaled); + hist_diff_rew_ref->GetXaxis()->SetTitleSize(text_size_scaled); + hist_diff_rew_ref->GetYaxis()->SetLabelSize(text_size_scaled); + hist_diff_rew_ref->GetYaxis()->SetTitleSize(text_size_scaled); + + hist_diff_rew_ref->GetYaxis()->SetTitleOffset(title_offset_scaled); + hist_diff_rew_ref->GetYaxis()->SetTitle("rew - ref / ref"); + hist_diff_rew_ref->GetYaxis()->CenterTitle(true); + + // hist_diff_rew_ref->SetMaximum(hist_diff_rew_ref->GetMaximum() + 0.5); + // hist_diff_rew_ref->SetMinimum(hist_diff_rew_ref->GetMinimum() - 0.5); + + hist_diff_rew_ref->Draw("histE1"); + c1->cd(); + pad1->Draw(); + pad2->Draw(); + c1->Draw(); + // std::string filename = + // varname + "_" + (do_normalize ? "normalized" : "unnormalized") + cut; + std::string filename = hist_rew->GetName(); + c1->SaveAs((filename + ".pdf").c_str()); + c1->SaveAs((filename + ".eps").c_str()); + }; + + for (auto &&[var, ent_ref, ent_rew, ent_unrew] : std::views::zip( + modelmap | std::views::keys, hists_ref, hists_rew, hists_unrew)) { + auto [hist_ref, cut_plot_ref] = ent_ref; + auto [hist_rew, cut_plot_rew] = ent_rew; + auto [hist_unrew, cut_plot_unrew] = ent_unrew; + std::cout << "Plotting " << var << '\n'; + do_plot(hist_rew.GetPtr(), hist_ref.GetPtr(), hist_unrew.GetPtr()); + for (auto &&[cutname, hist_ref, hist_rew, hist_unrew] : + std::views::zip(channelcuts | std::views::keys, cut_plot_ref, + cut_plot_rew, cut_plot_unrew)) { + std::cout << "Plotting " << var << " with cut: " << cutname << '\n'; + do_plot(hist_rew.GetPtr(), hist_ref.GetPtr(), hist_unrew.GetPtr()); + } + } + + // tree_rew.Snapshot("out", "out.root", {"weight", "Q2"}); + + // // do_plot("muonp", true); + // // do_plot("muonp", false); + // // do_plot("muon_theta", true); + // // do_plot("muon_theta", false); + // // do_plot("Q2", true); + // // do_plot("Q2", false); + // for (auto &&var : varnames) { + // // do_plot(var, false); + // do_plot(var, false); + // do_plot(var, false, "LowW"); + // do_plot(var, false, "HiW"); + // } + + return 0; +} \ No newline at end of file diff --git a/src/ProfSpline/KinematicVariables.h b/src/ProfSpline/KinematicVariables.h index fbb70496..3f6720b3 100644 --- a/src/ProfSpline/KinematicVariables.h +++ b/src/ProfSpline/KinematicVariables.h @@ -39,6 +39,15 @@ class ChannelIDs : public std::vector { } push_back(std::stoi(str)); } + + ChannelIDs subset(size_t size = 2) const { + ChannelIDs sub; + sub.resize(size); + for (size_t i = 0; i < size; ++i) { + sub[i] = (*this)[i]; + } + return sub; + } }; using KinematicVariables = std::vector; diff --git a/src/ProfSpline/LinkDef.h b/src/ProfSpline/LinkDef.h index cffe2218..4cff64f9 100644 --- a/src/ProfSpline/LinkDef.h +++ b/src/ProfSpline/LinkDef.h @@ -13,6 +13,9 @@ #pragma link C++ class genie::rew::ObservableMuonMomentum; #pragma link C++ class genie::rew::ObservablePMuEnu; #pragma link C++ class genie::rew::ObservablePMuEnuW; +// #pragma link C++ class genie::rew::ObservableHadronization; +#pragma link C++ class genie::rew::ObservableHadronizationLowW; +#pragma link C++ class genie::rew::ObservableHadronizationHighW; #pragma link C++ class genie::rew::KinematicVariables; #pragma link C++ class genie::rew::ChannelIDs; diff --git a/src/ProfSpline/ObservableHadronization.cxx b/src/ProfSpline/ObservableHadronization.cxx new file mode 100644 index 00000000..45e2880a --- /dev/null +++ b/src/ProfSpline/ObservableHadronization.cxx @@ -0,0 +1,167 @@ +#include + +#include "Framework/GHEP/GHepParticle.h" +#include "Framework/GHEP/GHepStatus.h" +#include "Framework/ParticleData/PDGUtils.h" +#include "Framework/Utils/KineUtils.h" +#include "ObservableHadronization.h" +#include "TLorentzVector.h" + +namespace genie::rew { + +// ObservableHadronization::ObservableHadronization() +// : RwgKineSpace("genie::rew::ObservableHadronization") {} + +// ObservableHadronization::ObservableHadronization(std::string config) +// : RwgKineSpace("genie::rew::ObservableHadronization", config) {} + +KinematicVariables ObservableHadronization::CalcKinematicVariables( + const EventRecord &event) const { + // This is a dummy implementation + std::vector ret; + + double W = event.Summary()->Kine().W(true); + ret.push_back(W); + return ret; +} + +ChannelIDs ObservableHadronization::ChannelID(const EventRecord &event) const { + ChannelIDs ret{}; + + if (event.TargetNucleus()) { + ret.push_back(event.Probe()->Pdg()); + ret.push_back(event.TargetNucleus()->Pdg()); + } else if (event.HitNucleon()) { + ret.push_back(event.Probe()->Pdg()); + ret.push_back(event.HitNucleon()->Pdg()); + } else { + ret.push_back(event.Probe()->Pdg()); + ret.push_back(0); + } + + const auto nucleus = ret[1]; + + size_t nch{}, nn{}; + auto ne = event.GetEntriesFast(); + for (int i{}; i < ne; ++i) { + auto particle = event.Particle(i); + if (particle->Status() == ((nucleus == 2212 || nucleus == 2112) + ? kIStStableFinalState + : kIStHadronInTheNucleus)) { + nch += (particle->Charge() != 0); + nn += (particle->Charge() == 0); + } + } + nch = std::min(nch, 17); + ret.push_back(nch); + + nn = std::min(nn, 3); + ret.push_back(nn); + + bool cc = event.Summary()->ProcInfo().IsWeakCC(); + ret.push_back(cc); + + return ret; +} + +bool ObservableHadronization::IsHandled(const EventRecord &event) const { + if (!(event.Summary()->ProcInfo().IsResonant() || + event.Summary()->ProcInfo().IsDeepInelastic())) { + return false; + } + + // exclude event with Kaons + for (int i{}; i < event.GetEntriesFast(); ++i) { + auto particle = event.Particle(i); + if (particle->Status() == kIStStableFinalState && + pdg::IsKaon(particle->Pdg())) { + return false; + } + } + + return true; +} + +void ObservableHadronization::LoadConfig(void) {} + +ObservableHadronizationLowW::ObservableHadronizationLowW() + : ObservableHadronization{"genie::rew::ObservableHadronizationLowW"} {} + +ObservableHadronizationLowW::ObservableHadronizationLowW(std::string config) + : ObservableHadronization{"genie::rew::ObservableHadronizationLowW", + config} {} + +KinematicVariables ObservableHadronizationLowW::CalcKinematicVariables( + const EventRecord &event) const { + auto ret = ObservableHadronization::CalcKinematicVariables(event); + double visE{}; + auto nentries = event.GetEntriesFast(); + for (int i = 0; i < nentries; i++) { + auto part = event.Particle(i); + if (part->Status() == genie::kIStStableFinalState && part->Charge() != 0) { + // remove mass part for p/n + visE += part->E(); + if (pdg::IsNeutron(part->Pdg()) || pdg::IsProton(part->Pdg())) { + visE -= part->Mass(); + } + } + } + ret.push_back(visE); + return ret; +} + +bool ObservableHadronizationLowW::IsHandled(const EventRecord &event) const { + if (!ObservableHadronization::IsHandled(event)) { + return false; + } + + return event.NEntries(92) == 0; +} + +ObservableHadronizationHighW::ObservableHadronizationHighW() + : ObservableHadronization{"genie::rew::ObservableHadronizationHighW"} {} + +ObservableHadronizationHighW::ObservableHadronizationHighW(std::string config) + : ObservableHadronization{"genie::rew::ObservableHadronizationHighW", + config} {} + +KinematicVariables ObservableHadronizationHighW::CalcKinematicVariables( + const EventRecord &event) const { + auto ret = ObservableHadronization::CalcKinematicVariables(event); + auto np = event.GetEntriesFast(); + TLorentzVector had_system{}; + for (int i{}; i < np; i++) { + auto particle = event.Particle(i); + if (particle->Status() == kIStStableFinalState && + pdg::IsHadron(particle->Pdg())) { + had_system += *particle->P4(); + } + } + auto had_system_boost = had_system.BoostVector(); + // auto had_system_dir = had_system.Vect().Unit(); + + // double sum_of_transverse_momentum{}; + double p_leading_pion{}; + for (int i{}; i < np; i++) { + auto particle = event.Particle(i); + if (particle->Status() == kIStStableFinalState && + pdg::IsPion(particle->Pdg())) { + auto p = *(particle->P4()); + p.Boost(-had_system_boost); + p_leading_pion = std::max(p_leading_pion, p.P()); + } + } + ret.push_back(p_leading_pion); + return ret; +} + +bool ObservableHadronizationHighW::IsHandled(const EventRecord &event) const { + if (!ObservableHadronization::IsHandled(event)) { + return false; + } + + // return event.Summary()->Kine().W(true) >= 3; + return event.NEntries(92) != 0; +} + +} // namespace genie::rew diff --git a/src/ProfSpline/ObservableHadronization.h b/src/ProfSpline/ObservableHadronization.h new file mode 100644 index 00000000..e47a6e63 --- /dev/null +++ b/src/ProfSpline/ObservableHadronization.h @@ -0,0 +1,45 @@ +#ifndef OBSERVABLEHADRONIZATION_H +#define OBSERVABLEHADRONIZATION_H +#include "ProfSpline/RwgKineSpace.h" + +namespace genie::rew { +class ObservableHadronization : public RwgKineSpace { +public: + using RwgKineSpace::RwgKineSpace; + + virtual KinematicVariables + CalcKinematicVariables(const EventRecord &event) const override; + + virtual ChannelIDs ChannelID(const EventRecord &event) const override; + + virtual ~ObservableHadronization() = default; + + virtual bool IsHandled(const EventRecord &event) const override; + +private: + virtual void LoadConfig(void) override; +}; + +class ObservableHadronizationLowW final : public ObservableHadronization { +public: + ObservableHadronizationLowW(); + ObservableHadronizationLowW(std::string config); + virtual KinematicVariables + CalcKinematicVariables(const EventRecord &event) const override; + virtual ~ObservableHadronizationLowW() = default; + virtual bool IsHandled(const EventRecord &event) const override; +}; + +class ObservableHadronizationHighW final : public ObservableHadronization { +public: + ObservableHadronizationHighW(); + ObservableHadronizationHighW(std::string config); + virtual KinematicVariables + CalcKinematicVariables(const EventRecord &event) const override; + virtual ~ObservableHadronizationHighW() = default; + virtual bool IsHandled(const EventRecord &event) const override; +}; + +} // namespace genie::rew + +#endif diff --git a/src/ProfSpline/ObservableSplines.cxx b/src/ProfSpline/ObservableSplines.cxx index 120696e1..8be2b12a 100644 --- a/src/ProfSpline/ObservableSplines.cxx +++ b/src/ProfSpline/ObservableSplines.cxx @@ -65,6 +65,11 @@ void ObservableSplines::InitializeIpols(const std::vector &lines) { LOG("ObservableSplines", pNOTICE) << "Initializing " << lines.size() << " bins"; for (const auto &line : lines) { + if (line.empty()) { + LOG("ObservableSplines", pERROR) + << "Empty line in spline definition, this could lead to errors"; + exit(1); + } bins.emplace_back(line); } } diff --git a/src/RwCalculators/GReWeightProfessor.cxx b/src/RwCalculators/GReWeightProfessor.cxx index 1ac3b7fa..fd6833c2 100644 --- a/src/RwCalculators/GReWeightProfessor.cxx +++ b/src/RwCalculators/GReWeightProfessor.cxx @@ -81,17 +81,25 @@ GReWeightProfessor::ReadProf2Spline(std::string filepath) { } else if (name == "Dimension") { std::stringstream ss{var}; } - } else if (line.find("#") != std::string::npos) { - // "/${orbname}/${orbname}__#..." - auto path = line.substr(0, line.find("#")); - auto name = path.substr(path.find_last_of("/") + 1); - auto orbname = path.substr(1, path.find_last_of("/") - 1); + } else if (auto loc_pound = line.find('#'); + loc_pound != std::string::npos) { + // "/${orbname}/${orbname}__#: ..." + auto path = line.substr(0, loc_pound); + auto loc_colon = path.find(':'); + auto id_str = path.substr(loc_pound + 1, loc_colon - loc_pound - 1); + auto id = std::stoul(id_str); + auto name = path.substr(path.find_last_of('/') + 1); + auto orbname = path.substr(1, path.find_last_of('/') - 1); auto sep_loc = name.find("__"); auto channelstr = name.substr(sep_loc + 2); ChannelIDs channel(channelstr, "_"); std::string errline{}; // not used now - // auto &varline = var_lines.emplace_back(); - auto &varline = ret[std::make_tuple(orbname, channel)].emplace_back(); + auto lines_vec = ret[std::make_tuple(orbname, channel)]; + if (lines_vec.size() <= id) { + lines_vec.resize(id + 1); + } + // auto &varline = ret[std::make_tuple(orbname, channel)].emplace_back(); + auto &varline = lines_vec[id]; std::getline(spline_file, varline); std::getline(spline_file, errline); } else {