diff --git a/app/PrepareGENIE.cxx b/app/PrepareGENIE.cxx index 9521eaf..ca21ce0 100644 --- a/app/PrepareGENIE.cxx +++ b/app/PrepareGENIE.cxx @@ -1,860 +1,866 @@ #include #include #include "FitLogger.h" #include "PlotUtils.h" #include "TFile.h" #include "TH1D.h" #include "TTree.h" #ifdef __GENIE_ENABLED__ +#ifdef GENIE_PRE_R3 #include "Conventions/Units.h" #include "GHEP/GHepParticle.h" #include "PDG/PDGUtils.h" +#else +#include "Framework/Conventions/Units.h" +#include "Framework/GHEP/GHepParticle.h" +#include "Framework/ParticleData/PDGUtils.h" +#endif #endif std::string gInputFiles = ""; std::string gOutputFile = ""; std::string gFluxFile = ""; std::string gTarget = ""; double MonoEnergy; int gNEvents = -999; bool IsMonoE = false; void PrintOptions(); void ParseOptions(int argc, char* argv[]); void RunGENIEPrepareMono(std::string input, std::string target, std::string output); void RunGENIEPrepare(std::string input, std::string flux, std::string target, std::string output); int main(int argc, char* argv[]) { ParseOptions(argc, argv); if (IsMonoE) { RunGENIEPrepareMono(gInputFiles, gTarget, gOutputFile); } else { RunGENIEPrepare(gInputFiles, gFluxFile, gTarget, gOutputFile); } } void RunGENIEPrepareMono(std::string input, std::string target, std::string output) { LOG(FIT) << "Running GENIE Prepare in mono energetic with E = " << MonoEnergy << " GeV" << std::endl; // Setup TTree TChain* tn = new TChain("gtree"); tn->AddFile(input.c_str()); if (tn->GetFile() == NULL) { tn->Print(); ERROR(FTL, "gtree not located in GENIE file: " << input); THROW("Check your inputs, they may need to be completely regenerated!"); throw; } int nevt = tn->GetEntries(); if (gNEvents != -999) { LOG(FIT) << "Overriding number of events by user from " << nevt << " to " << gNEvents << std::endl; nevt = gNEvents; } NtpMCEventRecord* genientpl = NULL; tn->SetBranchAddress("gmcrec", &genientpl); // Have the TH1D go from MonoEnergy/2 to MonoEnergy/2 TH1D* fluxhist = new TH1D("flux", "flux", 1000, MonoEnergy/2., MonoEnergy*2.); fluxhist->Fill(MonoEnergy); fluxhist->Scale(1, "width"); // Make Event Hist TH1D* eventhist = (TH1D*)fluxhist->Clone(); eventhist->Reset(); TH1D* xsechist = (TH1D*)eventhist->Clone(); // Create maps std::map modexsec; std::map modecount; std::vector genieids; std::vector targetids; std::vector interids; // Loop over all events for (int i = 0; i < nevt; i++) { tn->GetEntry(i); StopTalking(); EventRecord& event = *(genientpl->event); GHepParticle* neu = event.Probe(); StartTalking(); // Get XSec From Spline GHepRecord genie_record = static_cast(event); double xsec = (genie_record.XSec() / (1E-38 * genie::units::cm2)); // Parse Interaction String std::string mode = genie_record.Summary()->AsString(); std::vector modevec = GeneralUtils::ParseToStr(mode, ";"); std::string targ = (modevec[0] + ";" + modevec[1]); std::string inter = mode; // Fill lists of Unique IDS if (std::find(targetids.begin(), targetids.end(), targ) == targetids.end()) { targetids.push_back(targ); } if (std::find(interids.begin(), interids.end(), inter) == interids.end()) { interids.push_back(inter); } // Create entries Mode Maps if (modexsec.find(mode) == modexsec.end()) { genieids.push_back(mode); modexsec[mode] = (TH1D*)xsechist->Clone(); modecount[mode] = (TH1D*)xsechist->Clone(); modexsec[mode]->GetYaxis()->SetTitle("d#sigma/dE_{#nu} #times 10^{-38} (events weighted by #sigma)"); modecount[mode]->GetYaxis()->SetTitle("Number of events in file"); } // Fill XSec Histograms modexsec[mode]->Fill(neu->E(), xsec); modecount[mode]->Fill(neu->E()); // Fill total event hist eventhist->Fill(neu->E()); // Clear Event genientpl->Clear(); size_t freq = nevt / 20; if (freq && !(i % freq)) { LOG(FIT) << "Processed " << i << "/" << nevt << " GENIE events (E: " << neu->E() << " GeV, xsec: " << xsec << " E-38 cm^2/nucleon)" << std::endl; } } LOG(FIT) << "Processed all events" << std::endl; TFile* outputfile; // If no output is specified just append to the file if (!gOutputFile.length()) { tn->GetEntry(0); outputfile = tn->GetFile(); outputfile->cd(); } else { outputfile = new TFile(gOutputFile.c_str(), "RECREATE"); outputfile->cd(); QLOG(FIT, "Cloning input vector to output file: " << gOutputFile); TTree* cloneTree = tn->CloneTree(); cloneTree->SetDirectory(outputfile); cloneTree->Write(); QLOG(FIT, "Cloning input nova_wgts to output file: " << gOutputFile); // *********************************** // *********************************** // FUDGE FOR NOVA MINERVA WORKSHOP // Also check for the nova_wgts tree from Jeremy TChain *nova_chain = new TChain("nova_wgts"); nova_chain->AddFile(input.c_str()); TTree* nova_tree = nova_chain->GetTree(); if (!nova_tree) { QLOG(FIT, "Could not find nova_wgts tree in " << gOutputFile); } else { QLOG(FIT, "Found nova_wgts tree in " << gOutputFile); } if (nova_tree) { nova_tree->SetDirectory(outputfile); nova_tree->Write(); } QLOG(FIT, "Done cloning tree."); } LOG(FIT) << "Getting splines in mono-energetic..." << std::endl; // Save each of the reconstructed splines to file std::map modeavg; TDirectory* inddir = (TDirectory*)outputfile->Get("IndividualGENIESplines"); if (!inddir) inddir = (TDirectory*)outputfile->mkdir("IndividualGENIESplines"); inddir->cd(); // Loop over GENIE ID's and get MEC count int MECcount = 0; bool MECcorrect = FitPar::Config().GetParB("CorrectGENIEMECNorm"); for (UInt_t i = 0; i < genieids.size(); i++) { if (genieids[i].find("MEC") != std::string::npos) { MECcount++; } } LOG(FIT) << "Found " << MECcount << " repeated MEC instances." << std::endl; for (UInt_t i = 0; i < genieids.size(); i++) { std::string mode = genieids[i]; modexsec[mode]->Write((mode + "_summed_xsec").c_str(), TObject::kOverwrite); modecount[mode]->Write((mode + "_summed_evt").c_str(), TObject::kOverwrite); // Form extra avg xsec map -> Reconstructed spline modeavg[mode] = (TH1D*)modexsec[mode]->Clone(); modeavg[mode]->GetYaxis()->SetTitle("#sigma (E_{#nu}) #times 10^{-38} (cm^{2}/target)"); modeavg[mode]->Divide(modecount[mode]); if (MECcorrect && (mode.find("MEC") != std::string::npos)) { modeavg[mode]->Scale(1.0 / double(MECcount)); } modeavg[mode]->Write((mode + "_rec_spline").c_str(), TObject::kOverwrite); } TDirectory* targdir = (TDirectory*)outputfile->Get("TargetGENIESplines"); if (!targdir) targdir = (TDirectory*)outputfile->mkdir("TargetGENIESplines"); targdir->cd(); LOG(FIT) << "Getting Target Splines" << std::endl; // For each target save a total spline std::map targetsplines; for (uint i = 0; i < targetids.size(); i++) { std::string targ = targetids[i]; LOG(FIT) << "Getting target " << i << ": " << targ << std::endl; targetsplines[targ] = (TH1D*)xsechist->Clone(); targetsplines[targ]->GetYaxis()->SetTitle("#sigma (E_{#nu}) #times 10^{-38} (cm^{2}/target)"); LOG(FIT) << "Created target spline for " << targ << std::endl; for (uint j = 0; j < genieids.size(); j++) { std::string mode = genieids[j]; if (mode.find(targ) != std::string::npos) { LOG(FIT) << " Mode " << mode << " contains " << targ << " target" << std::endl; targetsplines[targ]->Add(modeavg[mode]); LOG(FIT) << "Finished with Mode " << mode << " " << modeavg[mode]->Integral() << std::endl; } } LOG(FIT) << "Saving target spline:" << targ << std::endl; targetsplines[targ]->Write(("Total_" + targ).c_str(), TObject::kOverwrite); } LOG(FIT) << "Getting total splines" << std::endl; // Now we have each of the targets we need to create a total cross-section. int totalnucl = 0; // Get the targets specified by the user, separated by commas // This has structure target1[fraction1], target2[fraction2] std::vector targprs = GeneralUtils::ParseToStr(target, ","); std::vector targ_list; std::vector frac_list; // Chop up the target string which has format TARGET1[fraction1],TARGET2[fraction2] //std::cout << "Targets: " << std::endl; // Loop over the vector of strings "TARGET1[fraction1]" "TARGET2[fraction2]" for (std::vector::iterator it = targprs.begin(); it != targprs.end(); ++it) { // Cut into "TARGET1" and "fraction1]" std::vector targind = GeneralUtils::ParseToStr(*it, "["); //std::cout << " " << *it << std::endl; // Cut into "TARGET1" and "fraction1" for (std::vector::iterator jt = targind.begin(); jt != targind.end(); ++jt) { if ((*jt).find("]") != std::string::npos) { (*jt) = (*jt).substr(0, (*jt).find("]")); //*jt = "hello"; frac_list.push_back(*jt); // Won't find bracket for target } else { targ_list.push_back(*jt); } } } targprs = targ_list; std::vector targ_fractions; double minimum = 1.0; for (std::vector::iterator it = frac_list.begin(); it != frac_list.end(); it++) { //std::cout << " " << *it << std::endl; double frac = std::atof((*it).c_str()); targ_fractions.push_back(frac); if (frac < minimum) minimum = frac; } std::vector::iterator it = targ_fractions.begin(); std::vector::iterator jt = targ_list.begin(); double scaling = 0; for (; it != targ_fractions.end(); it++, jt++) { // First get the mass number from the targ_list int nucl = atoi((*jt).c_str()); nucl = (nucl%10000)/10; // Gets the relative portions right *it = (*it)/minimum; // Scale relative the atomic mass //(*it) *= (double(nucl)/(*it)); double tempscaling = double(nucl)/(*it); if (tempscaling > scaling) scaling=tempscaling; } it = targ_fractions.begin(); for (; it != targ_fractions.end(); it++) { // Round the scaling to nearest integer and multiply *it *= int(scaling+0.5); // Round to nearest integer *it = int(*it+0.5); totalnucl += *it; } if (totalnucl == 0) { THROW("Didn't find any nucleons in input file. Did you really specify the target ratios?\ne.g. TARGET1[fraction1],TARGET2[fraction2]" << std::endl); } TH1D* totalxsec = (TH1D*)xsechist->Clone(); for (uint i = 0; i < targprs.size(); i++) { std::string targpdg = targprs[i]; // Check that we found the user requested target in GENIE bool FoundTarget = false; for (std::map::iterator iter = targetsplines.begin(); iter != targetsplines.end(); iter++) { std::string targstr = iter->first; TH1D* xsec = iter->second; // Match the user targets to the targets found in GENIE if (targstr.find(targpdg) != std::string::npos) { FoundTarget = true; LOG(FIT) << "Adding target spline " << targstr << " Integral = " << xsec->Integral("width") << std::endl; totalxsec->Add(xsec); //int nucl = atoi(targpdg.c_str()); //totalnucl += int((nucl % 10000) / 10); } } // Check that targets were all found if (!FoundTarget) { ERR(WRN) << "Didn't find target " << targpdg << " in the list of targets recorded by GENIE" << std::endl; ERR(WRN) << " The list of targets you requested is: " << std::endl; for (uint i = 0; i < targprs.size(); ++i) ERR(WRN) << " " << targprs[i] << std::endl; ERR(WRN) << " The list of targets found in GENIE is: " << std::endl; for (std::map::iterator iter = targetsplines.begin(); iter != targetsplines.end(); iter++) ERR(WRN) << " " << iter->first<< std::endl; } } outputfile->cd(); totalxsec->GetYaxis()->SetTitle("#sigma (E_{#nu}) #times 10^{-38} (cm^{2}/nucleon)"); totalxsec->Write("nuisance_xsec", TObject::kOverwrite); eventhist = (TH1D*)fluxhist->Clone(); eventhist->Multiply(totalxsec); eventhist->GetYaxis()->SetTitle((std::string("Event rate (N = #sigma #times #Phi) #times 10^{-38} (cm^{2}/nucleon) #times ")+eventhist->GetYaxis()->GetTitle()).c_str()); LOG(FIT) << "Dividing by Total Nucl = " << totalnucl << std::endl; eventhist->Scale(1.0 / double(totalnucl)); eventhist->Write("nuisance_events", TObject::kOverwrite); fluxhist->Write("nuisance_flux", TObject::kOverwrite); LOG(FIT) << "Inclusive XSec Per Nucleon = " << eventhist->Integral("width") * 1E-38 / fluxhist->Integral("width") << std::endl; LOG(FIT) << "XSec Hist Integral = " << totalxsec->Integral("width") << std::endl; outputfile->Close(); return; } void RunGENIEPrepare(std::string input, std::string flux, std::string target, std::string output) { LOG(FIT) << "Running GENIE Prepare with flux..." << std::endl; // Get Flux Hist std::vector fluxvect = GeneralUtils::ParseToStr(flux, ","); TH1* fluxhist = NULL; if (fluxvect.size() == 3) { double from = GeneralUtils::StrToDbl(fluxvect[0]); double to = GeneralUtils::StrToDbl(fluxvect[1]); double step = GeneralUtils::StrToDbl(fluxvect[2]); int nstep = ceil((to - from) / step); to = from + step * nstep; QLOG(FIT, "Generating flat flux histogram from " << from << " to " << to << " with bins " << step << " wide (NBins = " << nstep << ")."); fluxhist = new TH1D("spectrum", ";E_{#nu} (GeV);Count (A.U.)", nstep, from, to); for (Int_t bi_it = 1; bi_it < fluxhist->GetXaxis()->GetNbins(); ++bi_it) { fluxhist->SetBinContent(bi_it, 1.0 / double(step * nstep)); } fluxhist->SetDirectory(0); } else if (fluxvect.size() == 2) { TFile* fluxfile = new TFile(fluxvect[0].c_str(), "READ"); if (!fluxfile->IsZombie()) { fluxhist = dynamic_cast(fluxfile->Get(fluxvect[1].c_str())); if (!fluxhist) { ERR(FTL) << "Couldn't find histogram named: \"" << fluxvect[1] << "\" in file: \"" << fluxvect[0] << std::endl; throw; } fluxhist->SetDirectory(0); } } else if (fluxvect.size() == 1) { MonoEnergy = GeneralUtils::StrToDbl(fluxvect[0]); RunGENIEPrepareMono(input, target, output); return; } else { LOG(FTL) << "Bad flux specification: \"" << flux << "\"." << std::endl; throw; } // Setup TTree TChain* tn = new TChain("gtree"); if (input.find_first_of(',') != std::string::npos) { std::vector inputvect = GeneralUtils::ParseToStr(input, ","); for (size_t iv_it = 0; iv_it < inputvect.size(); ++iv_it) { tn->AddFile(inputvect[iv_it].c_str()); QLOG(FIT, "Added input file: " << inputvect[iv_it]); } } else { // The Add form can accept wildcards. tn->Add(input.c_str()); } if (tn->GetFile() == NULL) { tn->Print(); ERROR(FTL, "gtree not located in GENIE file: " << input); THROW("Check your inputs, they may need to be completely regenerated!"); throw; } int nevt = tn->GetEntries(); if (gNEvents != -999) { LOG(FIT) << "Overriding number of events by user from " << nevt << " to " << gNEvents << std::endl; nevt = gNEvents; } if (!nevt) { THROW("Couldn't load any events from input specification: \"" << input.c_str() << "\""); } else { QLOG(FIT, "Found " << nevt << " input entries in " << input); } NtpMCEventRecord* genientpl = NULL; tn->SetBranchAddress("gmcrec", &genientpl); // Make Event and xsec Hist TH1D* eventhist = (TH1D*)fluxhist->Clone(); eventhist->Reset(); TH1D* xsechist = (TH1D*)eventhist->Clone(); // Create maps std::map modexsec; std::map modecount; std::vector genieids; std::vector targetids; std::vector interids; // Loop over all events for (int i = 0; i < nevt; i++) { tn->GetEntry(i); // Hussssch GENIE StopTalking(); // Get the event EventRecord& event = *(genientpl->event); // Get the neutrino GHepParticle* neu = event.Probe(); StartTalking(); // Get XSec From Spline // Get the GHepRecord GHepRecord genie_record = static_cast(event); double xsec = (genie_record.XSec() / (1E-38 * genie::units::cm2)); // Parse Interaction String std::string mode = genie_record.Summary()->AsString(); std::vector modevec = GeneralUtils::ParseToStr(mode, ";"); std::string targ = (modevec[0] + ";" + modevec[1]); std::string inter = mode; // Get target nucleus // Alternative ways of getting the summaries //GHepParticle *target = genie_record.TargetNucleus(); //int pdg = target->Pdg(); // Fill lists of Unique IDS (neutrino and target) if (std::find(targetids.begin(), targetids.end(), targ) == targetids.end()) { targetids.push_back(targ); } // The full interaction list if (std::find(interids.begin(), interids.end(), inter) == interids.end()) { interids.push_back(inter); } // Create entries Mode Maps if (modexsec.find(mode) == modexsec.end()) { genieids.push_back(mode); modexsec[mode] = (TH1D*)xsechist->Clone(); modecount[mode] = (TH1D*)xsechist->Clone(); modexsec[mode]->GetYaxis()->SetTitle("d#sigma/dE_{#nu} #times 10^{-38} (events weighted by #sigma)"); modecount[mode]->GetYaxis()->SetTitle("Number of events in file"); } // Fill XSec Histograms modexsec[mode]->Fill(neu->E(), xsec); modecount[mode]->Fill(neu->E()); // Fill total event hist eventhist->Fill(neu->E()); if (i % (nevt / 20) == 0) { LOG(FIT) << "Processed " << i << "/" << nevt << " GENIE events (E: " << neu->E() << " GeV, xsec: " << xsec << " E-38 cm^2/nucleon)" << std::endl; } // Clear Event genientpl->Clear(); } LOG(FIT) << "Processed all events" << std::endl; // Once event loop is done we can start saving stuff into the file TFile* outputfile; if (!gOutputFile.length()) { tn->GetEntry(0); outputfile = tn->GetFile(); outputfile->cd(); } else { outputfile = new TFile(gOutputFile.c_str(), "RECREATE"); outputfile->cd(); QLOG(FIT, "Cloning input vector to output file: " << gOutputFile); TTree* cloneTree = tn->CloneTree(); cloneTree->SetDirectory(outputfile); cloneTree->Write(); // ******************************** // CLUDGE KLUDGE KLUDGE FOR NOVA QLOG(FIT, "Cloning input nova_wgts to output file: " << gOutputFile); // Also check for the nova_wgts tree from Jeremy TChain *nova_chain = new TChain("nova_wgts"); nova_chain->AddFile(input.c_str()); TTree* nova_tree = nova_chain->CloneTree(); if (!nova_tree) { QLOG(FIT, "Could not find nova_wgts tree in " << input); } else { QLOG(FIT, "Found nova_wgts tree in " << input); nova_tree->SetDirectory(outputfile); nova_tree->Write(); } QLOG(FIT, "Done cloning tree."); } LOG(FIT) << "Getting splines..." << std::endl; // Save each of the reconstructed splines to file std::map modeavg; TDirectory* inddir = (TDirectory*)outputfile->Get("IndividualGENIESplines"); if (!inddir) inddir = (TDirectory*)outputfile->mkdir("IndividualGENIESplines"); inddir->cd(); // Loop over GENIE ID's and get MEC count int MECcount = 0; bool MECcorrect = FitPar::Config().GetParB("CorrectGENIEMECNorm"); for (UInt_t i = 0; i < genieids.size(); i++) { if (genieids[i].find("MEC") != std::string::npos) { MECcount++; } } LOG(FIT) << "Found " << MECcount << " repeated MEC instances." << std::endl; for (UInt_t i = 0; i < genieids.size(); i++) { std::string mode = genieids[i]; modexsec[mode]->Write((mode + "_summed_xsec").c_str(), TObject::kOverwrite); modecount[mode]->Write((mode + "_summed_evt").c_str(), TObject::kOverwrite); // Form extra avg xsec map -> Reconstructed spline modeavg[mode] = (TH1D*)modexsec[mode]->Clone(); modeavg[mode]->GetYaxis()->SetTitle("#sigma (E_{#nu}) #times 10^{-38} (cm^{2}/target)"); modeavg[mode]->Divide(modecount[mode]); if (MECcorrect && (mode.find("MEC") != std::string::npos)) { modeavg[mode]->Scale(1.0 / double(MECcount)); } modeavg[mode]->Write((mode + "_rec_spline").c_str(), TObject::kOverwrite); } TDirectory* targdir = (TDirectory*)outputfile->Get("TargetGENIESplines"); if (!targdir) targdir = (TDirectory*)outputfile->mkdir("TargetGENIESplines"); targdir->cd(); LOG(FIT) << "Getting Target Splines" << std::endl; // For each target save a total spline std::map targetsplines; for (uint i = 0; i < targetids.size(); i++) { std::string targ = targetids[i]; LOG(FIT) << "Getting target " << i << ": " << targ << std::endl; targetsplines[targ] = (TH1D*)xsechist->Clone(); targetsplines[targ]->GetYaxis()->SetTitle("#sigma (E_{#nu}) #times 10^{-38} (cm^{2}/target)"); LOG(FIT) << "Created target spline for " << targ << std::endl; for (uint j = 0; j < genieids.size(); j++) { std::string mode = genieids[j]; // Look at all matching modes/targets if (mode.find(targ) != std::string::npos) { LOG(FIT) << " Mode " << mode << " contains " << targ << " target" << std::endl; targetsplines[targ]->Add(modeavg[mode]); LOG(FIT) << "Finished with Mode " << mode << " " << modeavg[mode]->Integral() << std::endl; } } LOG(FIT) << "Saving target spline: " << targ << std::endl; targetsplines[targ]->Write(("Total_" + targ).c_str(), TObject::kOverwrite); } LOG(FIT) << "Getting total splines" << std::endl; // Now we have each of the targets we need to create a total cross-section. int totalnucl = 0; // This has structure target1[fraction1], target2[fraction2] std::vector targprs = GeneralUtils::ParseToStr(target, ","); std::vector targ_list; std::vector frac_list; // Chop up the target string which has format TARGET1[fraction1],TARGET2[fraction2] //std::cout << "Targets: " << std::endl; // Loop over the vector of strings "TARGET1[fraction1]" "TARGET2[fraction2]" for (std::vector::iterator it = targprs.begin(); it != targprs.end(); ++it) { // Cut into "TARGET1" and "fraction1]" std::vector targind = GeneralUtils::ParseToStr(*it, "["); //std::cout << " " << *it << std::endl; // Cut into "TARGET1" and "fraction1" for (std::vector::iterator jt = targind.begin(); jt != targind.end(); ++jt) { if ((*jt).find("]") != std::string::npos) { (*jt) = (*jt).substr(0, (*jt).find("]")); //*jt = "hello"; frac_list.push_back(*jt); // Won't find bracket for target } else { targ_list.push_back(*jt); } } } targprs = targ_list; std::vector targ_fractions; double minimum = 1.0; for (std::vector::iterator it = frac_list.begin(); it != frac_list.end(); it++) { //std::cout << " " << *it << std::endl; double frac = std::atof((*it).c_str()); targ_fractions.push_back(frac); if (frac < minimum) minimum = frac; } std::vector::iterator it = targ_fractions.begin(); std::vector::iterator jt = targ_list.begin(); double scaling = 0; for (; it != targ_fractions.end(); it++, jt++) { // First get the mass number from the targ_list int nucl = atoi((*jt).c_str()); nucl = (nucl%10000)/10; // Gets the relative portions right *it = (*it)/minimum; // Scale relative the atomic mass //(*it) *= (double(nucl)/(*it)); double tempscaling = double(nucl)/(*it); if (tempscaling > scaling) scaling=tempscaling; } it = targ_fractions.begin(); for (; it != targ_fractions.end(); it++) { // Round the scaling to nearest integer and multiply *it *= int(scaling+0.5); // Round to nearest integer *it = int(*it+0.5); totalnucl += *it; } if (totalnucl == 0) { THROW("Didn't find any nucleons in input file. Did you really specify the target ratios?\ne.g. TARGET1[fraction1],TARGET2[fraction2]" << std::endl); } TH1D* totalxsec = (TH1D*)xsechist->Clone(); // Loop over the specified targets by the user for (uint i = 0; i < targprs.size(); i++) { std::string targpdg = targprs[i]; // Check that we found the user requested target in GENIE bool FoundTarget = false; for (std::map::iterator iter = targetsplines.begin(); iter != targetsplines.end(); iter++) { std::string targstr = iter->first; TH1D* xsec = iter->second; // Match the user targets to the targets found in GENIE if (targstr.find(targpdg) != std::string::npos) { FoundTarget = true; LOG(FIT) << "Adding target spline " << targstr << " Integral = " << xsec->Integral("width") << std::endl; totalxsec->Add(xsec); //int nucl = atoi(targpdg.c_str()); //totalnucl += int((nucl % 10000) / 10); } } // Looped over target splines // Check that targets were all found if (!FoundTarget) { ERR(WRN) << "Didn't find target " << targpdg << " in the list of targets recorded by GENIE" << std::endl; ERR(WRN) << " The list of targets you requested is: " << std::endl; for (uint i = 0; i < targprs.size(); ++i) ERR(WRN) << " " << targprs[i] << std::endl; ERR(WRN) << " The list of targets found in GENIE is: " << std::endl; for (std::map::iterator iter = targetsplines.begin(); iter != targetsplines.end(); iter++) ERR(WRN) << " " << iter->first<< std::endl; } } outputfile->cd(); totalxsec->GetYaxis()->SetTitle("#sigma (E_{#nu}) #times 10^{-38} (cm^{2}/nucleon)"); totalxsec->Write("nuisance_xsec", TObject::kOverwrite); eventhist = (TH1D*)fluxhist->Clone(); eventhist->Multiply(totalxsec); eventhist->GetYaxis()->SetTitle((std::string("Event rate (N = #sigma #times #Phi) #times 10^{-38} (cm^{2}/nucleon) #times ")+eventhist->GetYaxis()->GetTitle()).c_str()); LOG(FIT) << "Dividing by Total Nucl = " << totalnucl << std::endl; eventhist->Scale(1.0 / double(totalnucl)); eventhist->Write("nuisance_events", TObject::kOverwrite); fluxhist->Write("nuisance_flux", TObject::kOverwrite); LOG(FIT) << "Inclusive XSec Per Nucleon = " << eventhist->Integral("width") * 1E-38 / fluxhist->Integral("width") << std::endl; LOG(FIT) << "XSec Hist Integral = " << totalxsec->Integral() << std::endl; outputfile->Close(); return; }; void PrintOptions() { std::cout << "PrepareGENIEEvents NUISANCE app. " << std::endl << "Takes GHep Outputs and prepares events for NUISANCE." << std::endl << std::endl << "PrepareGENIE [-h,-help,--h,--help] [-i " "inputfile1.root,inputfile2.root,inputfile3.root,...] " << "[-f flux_root_file.root,flux_hist_name] [-t " "target1[frac1],target2[frac2],...]" << "[-n number_of_events (experimental)]" << std::endl << std::endl; std::cout << "Prepare Mode [Default] : Takes a single GHep file, " "reconstructs the original GENIE splines, " << " and creates a duplicate file that also contains the flux, " "event rate, and xsec predictions that NUISANCE needs. " << std::endl; std::cout << "Following options are required for Prepare Mode:" << std::endl; std::cout << " [ -i inputfile.root ] : Reads in a single GHep input file " "that needs the xsec calculation ran on it. " << std::endl; std::cout << " [ -f flux_file.root,hist_name ] : Path to root file " "containing the flux histogram the GHep records were generated " "with." << " A simple method is to point this to the flux histogram genie " "generatrs '-f /path/to/events/input-flux.root,spectrum'. " << std::endl; std::cout << " [ -f elow,ehigh,estep ] : Energy range specification when no " "flux file was used." << std::endl; std::cout << " [ -t target ] : Target that GHepRecords were generated with. " "Comma seperated list. E.g. for CH2 " "target=1000060120,1000010010,1000010010" << std::endl; std::cout << " [ -o outputfile.root ] : File to write prepared input file to." << std::endl; std::cout << " [ -m Mono_E_nu_GeV ] : Run in mono-energetic mode with m GeV neutrino energy." << std::endl; std::cout << " [ -n number_of_evt ] : Run with a reduced number of events for debugging purposes" << std::endl; } void ParseOptions(int argc, char* argv[]) { bool flagopt = false; // If No Arguments print commands for (int i = 1; i < argc; ++i) { if (!std::strcmp(argv[i], "-h")) { flagopt = true; break; } if (i + 1 != argc) { // Cardfile if (!std::strcmp(argv[i], "-h")) { flagopt = true; break; } else if (!std::strcmp(argv[i], "-i")) { gInputFiles = argv[i + 1]; ++i; } else if (!std::strcmp(argv[i], "-o")) { gOutputFile = argv[i + 1]; ++i; } else if (!std::strcmp(argv[i], "-f")) { gFluxFile = argv[i + 1]; ++i; } else if (!std::strcmp(argv[i], "-t")) { gTarget = argv[i + 1]; ++i; } else if (!std::strcmp(argv[i], "-n")) { gNEvents = GeneralUtils::StrToInt(argv[i + 1]); ++i; } else if (!std::strcmp(argv[i], "-m")) { MonoEnergy = GeneralUtils::StrToDbl(argv[i + 1]); IsMonoE = true; ++i; } else { ERR(FTL) << "ERROR: unknown command line option given! - '" << argv[i] << " " << argv[i + 1] << "'" << std::endl; PrintOptions(); break; } } } if (gInputFiles == "" && !flagopt) { ERR(FTL) << "No input file(s) specified!" << std::endl; flagopt = true; } if (gFluxFile == "" && !flagopt && !IsMonoE) { ERR(FTL) << "No flux input specified for Prepare Mode" << std::endl; flagopt = true; } if (gTarget == "" && !flagopt) { ERR(FTL) << "No target specified for Prepare Mode" << std::endl; flagopt = true; } if (gTarget.find("[") == std::string::npos || gTarget.find("]") == std::string::npos) { ERR(FTL) << "Didn't specify target ratios in Prepare Mode" << std::endl; ERR(FTL) << "Are you sure you gave it as -t \"TARGET1[fraction1],TARGET2[fraction]\"?" << std::endl; flagopt = true; } if (argc < 1 || flagopt) { PrintOptions(); exit(-1); } return; } diff --git a/app/nuismin.cxx b/app/nuismin.cxx index 5bb48d8..afe4224 100644 --- a/app/nuismin.cxx +++ b/app/nuismin.cxx @@ -1,113 +1,113 @@ // Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ // Author: Callum Wilkinson 01/2014 // Patrick Stowell 09/2015 /** Usage: ./GSLminimizerWithReWeight.exe -c card file, where samples and parameters are defined -o output file, where the results of the fit are stored where: */ #include "MinimizerRoutines.h" //******************************* void printInputCommands(){ //******************************* std::cout<<"ExtFit_minimizer.exe -c cardFile -o outFile [-f fitStategy] [-d fakeDataFile] [-i inputFile] [-q config_name=config_val] \n"; std::cout<SaveNominal(); if (FitPar::Config().GetParB("saveprefit")) min->SavePrefit(); - // Run the fit rotines + // Run the fit routines min->Run(); // Save by default min->SaveResults(); // Get Status status = min->GetStatus(); // Show Final Status LOG(FIT)<<"-------------------------------------"<. ################################################################################ # TODO # check system for libxml2 # check whether we need the includes # check if we can use a subset of the GENIE libraries ################################################################################ # Check Dependencies ################################################################################ ################################# GENIE ###################################### if(GENIE STREQUAL "") cmessage(FATAL_ERROR "Variable GENIE is not defined. " "The location of a pre-built GENIE install must be defined either as" " $ cmake -DGENIE=/path/to/GENIE or as and environment vairable" " $ export GENIE=/path/to/GENIE") endif() if (BUILD_GEVGEN) cmessage(STATUS "Building custom gevgen") LIST(APPEND EXTRA_CXX_FLAGS -D__GEVGEN_ENABLED__) endif() if(GENIE_EMPMEC_REWEIGHT) cmessage(STATUS "Enable EMPMEC dials") LIST(APPEND EXTRA_CXX_FLAGS -D__GENIE_EMP_MECRW_ENABLED) endif() # Extract GENIE VERSION if (GENIE_VERSION STREQUAL "AUTO") execute_process (COMMAND ${CMAKE_SOURCE_DIR}/cmake/getgenieversion.sh ${GENIE} OUTPUT_VARIABLE GENIE_VERSION OUTPUT_STRIP_TRAILING_WHITESPACE) endif() +execute_process(COMMAND genie-config --version +OUTPUT_VARIABLE GENIE_VER OUTPUT_STRIP_TRAILING_WHITESPACE) +cmessage(STATUS "genie_ver: ${GENIE_VER}") +if(GENIE_VER VERSION_GREATER 3.0.0) + set(GENIE_POST_R3 1) + string(REPLACE "." "" GENIE_VERSION ${GENIE_VER}) + cmessage(STATUS "set genie_post_r3") +endif() + +if(NOT GENIE_POST_R3) +LIST(APPEND EXTRA_CXX_FLAGS -DGENIE_PRE_R3) +cmessage(STATUS "setting genie_pre_r3 ${EXTRA_CXX_FLAGS}") +endif() + execute_process (COMMAND genie-config --libs OUTPUT_VARIABLE GENIE_LD_FLAGS_STR OUTPUT_STRIP_TRAILING_WHITESPACE) execute_process (COMMAND genie-config --topsrcdir OUTPUT_VARIABLE GENIE_INCLUDES_DIR OUTPUT_STRIP_TRAILING_WHITESPACE) string(REGEX MATCH "-L\([^ ]+\) \(.*\)$" PARSE_GENIE_LIBS_MATCH ${GENIE_LD_FLAGS_STR}) cmessage(DEBUG "genie-config --libs: ${GENIE_LD_FLAGS_STR}") if(NOT PARSE_GENIE_LIBS_MATCH) cmessage(FATAL_ERROR "Expected to be able to parse the result of genie-config --libs to a lib directory and a list of libraries to include, but got: \"${GENIE_LD_FLAGS_STR}\"") endif() -set(GENIE_LIB_DIR ${CMAKE_MATCH_1}) +if(${CMAKE_MATCH_1}) + set(GENIE_LIB_DIR ${CMAKE_MATCH_1}) +else() + set(GENIE_LIB_DIR $ENV{GENIE_LIB}) +endif() + set(GENIE_LIBS_RAW ${CMAKE_MATCH_2}) string(REPLACE "-l" "" GENIE_LIBS_STRIPED "${GENIE_LIBS_RAW}") cmessage(STATUS "GENIE version : ${GENIE_VERSION}") cmessage(STATUS "GENIE libdir : ${GENIE_LIB_DIR}") cmessage(STATUS "GENIE libs : ${GENIE_LIBS_STRIPED}") string(REGEX MATCH "ReinSeghal" WASMATCHED ${GENIE_LIBS_STRIPED}) if(WASMATCHED AND GENIE_VERSION STREQUAL "210") set(GENIE_SEHGAL ${GENIE_LIBS_STRIPED}) STRING(REPLACE "ReinSeghal" "ReinSehgal" GENIE_LIBS_STRIPED ${GENIE_SEHGAL}) cmessage(DEBUG "Fixed inconsistency in library naming: ${GENIE_LIBS_STRIPED}") endif() -string(REGEX MATCH "ReWeight" WASMATCHED ${GENIE_LIBS_STRIPED}) -if(NOT WASMATCHED) - set(GENIE_LIBS_STRIPED "GReWeight ${GENIE_LIBS_STRIPED}") - cmessage(DEBUG "Force added ReWeight library: ${GENIE_LIBS_STRIPED}") +if(NOT GENIE_POST_R3) + string(REGEX MATCH "ReWeight" WASMATCHED ${GENIE_LIBS_STRIPED}) + if(NOT WASMATCHED) + set(GENIE_LIBS_STRIPED "GReWeight ${GENIE_LIBS_STRIPED}") + cmessage(DEBUG "Force added ReWeight library: ${GENIE_LIBS_STRIPED}") + endif() +else() + string(REGEX MATCH "RwFwk" WASMATCHED ${GENIE_LIBS_STRIPED}) + if(NOT WASMATCHED) + set(GENIE_LIBS_STRIPED "GRwClc GRwFwk GRwIO ${GENIE_LIBS_STRIPED}") + cmessage(DEBUG "Force added ReWeight library: ${GENIE_LIBS_STRIPED}") + endif() endif() string(REPLACE " " ";" GENIE_LIBS_LIST "-Wl,--no-as-needed -Wl,--start-group ${GENIE_LIBS_STRIPED} -Wl,--end-group") cmessage(DEBUG "genie-config --libs -- MATCH1: ${CMAKE_MATCH_1}") cmessage(DEBUG "genie-config --libs -- MATCH2: ${CMAKE_MATCH_2}") cmessage(DEBUG "genie-config --libs -- libs stripped: ${GENIE_LIBS_STRIPED}") cmessage(DEBUG "genie-config --libs -- libs list: ${GENIE_LIBS_LIST}") ################################ LHAPDF ###################################### if(LHAPDF_LIB STREQUAL "") cmessage(FATAL_ERROR "Variable LHAPDF_LIB is not defined. The location of a pre-built lhapdf install must be defined either as $ cmake -DLHAPDF_LIB=/path/to/LHAPDF_libraries or as and environment vairable $ export LHAPDF_LIB=/path/to/LHAPDF_libraries") endif() if(LHAPDF_INC STREQUAL "") cmessage(FATAL_ERROR "Variable LHAPDF_INC is not defined. The location of a pre-built lhapdf install must be defined either as $ cmake -DLHAPDF_INC=/path/to/LHAPDF_includes or as and environment vairable $ export LHAPDF_INC=/path/to/LHAPDF_includes") endif() if(LHAPATH STREQUAL "") cmessage(FATAL_ERROR "Variable LHAPATH is not defined. The location of a the LHAPATH directory must be defined either as $ cmake -DLHAPATH=/path/to/LHAPATH or as and environment variable $ export LHAPATH=/path/to/LHAPATH") endif() ################################ LIBXML ###################################### if(LIBXML2_LIB STREQUAL "") # Check for xml2-config find_program(LIBXMLCFGLIBS xml2-config) if(NOT LIBXMLCFGLIBS STREQUAL "LIBXMLLIBS-NOTFOUND") execute_process (COMMAND ${LIBXMLCFGLIBS} --libs OUTPUT_VARIABLE LIBXML2_LIB OUTPUT_STRIP_TRAILING_WHITESPACE) else() message(FATAL_ERROR "Variable LIBXML2_LIB is not defined and could not find xml2-config. The location of a pre-built libxml2 install must be defined either as $ cmake -DLIBXML2_LIB=/path/to/LIBXML2_libraries or as and environment vairable $ export LIBXML2_LIB=/path/to/LIBXML2_libraries") endif() endif() if(LIBXML2_INC STREQUAL "") # Check for xml2-config find_program(LIBXMLCFGINCS xml2-config) if(NOT LIBXMLCFGINCS STREQUAL "LIBXMLINCS-NOTFOUND") execute_process (COMMAND ${LIBXMLCFGINCS} --cflags OUTPUT_VARIABLE LIBXML2_INC OUTPUT_STRIP_TRAILING_WHITESPACE) else() message(FATAL_ERROR "Variable LIBXML2_INC is not defined and could not find xml2-config. The location of a pre-built libxml2 install must be defined either as $ cmake -DLIBXML2_INC=/path/to/LIBXML2_libraries or as and environment vairable $ export LIBXML2_INC=/path/to/LIBXML2_libraries") endif() endif() +string(REPLACE "-llzma" "" LIBXML2_LIB ${LIBXML2_LIB}) + ############################### log4cpp ###################################### if(LOG4CPP_LIB STREQUAL "") find_program(LOG4CPPCFG log4cpp-config) if(NOT LOG4CPPCFG STREQUAL "LOG4CPPCFG-NOTFOUND") execute_process (COMMAND ${LOG4CPPCFG} --pkglibdir OUTPUT_VARIABLE LOG4CPP_LIB OUTPUT_STRIP_TRAILING_WHITESPACE) else() message(FATAL_ERROR "Variable LOG4CPP_LIB is not defined. The location of a pre-built log4cpp install must be defined either as $ cmake -DLOG4CPP_LIB=/path/to/LOG4CPP_libraries or as and environment vairable $ export LOG4CPP_LIB=/path/to/LOG4CPP_libraries") endif() endif() if(LOG4CPP_INC STREQUAL "") find_program(LOG4CPPCFG log4cpp-config) if(NOT LOG4CPPCFG STREQUAL "LOG4CPPCFG-NOTFOUND") execute_process (COMMAND ${LOG4CPPCFG} --pkgincludedir OUTPUT_VARIABLE LOG4CPP_INC OUTPUT_STRIP_TRAILING_WHITESPACE) else() message(FATAL_ERROR "Variable LOG4CPP_INC is not defined. The location of a pre-built log4cpp install must be defined either as $ cmake -DGENIE_LOG4CPP_INC=/path/to/LOG4CPP_includes or as and environment vairable $ export LOG4CPP_INC=/path/to/LOG4CPP_includes") endif() endif() ################################################################################ LIST(APPEND EXTRA_CXX_FLAGS -D__GENIE_ENABLED__ -D__GENIE_VERSION__=${GENIE_VERSION}) -LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES - ${GENIE_INCLUDES_DIR} - ${GENIE_INCLUDES_DIR}/GHEP - ${GENIE_INCLUDES_DIR}/Ntuple - ${GENIE_INCLUDES_DIR}/ReWeight - ${GENIE_INCLUDES_DIR}/Apps - ${GENIE_INCLUDES_DIR}/FluxDrivers - ${GENIE_INCLUDES_DIR}/EVGDrivers - ${LHAPDF_INC} - ${LIBXML2_INC} - ${LOG4CPP_INC}) +LIST(APPEND EXTRA_LIBS ${GENIE_LIBS_LIST}) -SAYVARS() +LIST(APPEND EXTRA_LIBS LHAPDF xml2 log4cpp) LIST(APPEND EXTRA_LINK_DIRS ${GENIE_LIB_DIR} ${LHAPDF_LIB} ${LIBXML2_LIB} ${LOG4CPP_LIB}) -LIST(APPEND EXTRA_LIBS ${GENIE_LIBS_LIST}) +if(NOT GENIE_POST_R3) + LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES + ${GENIE_INCLUDES_DIR} + ${GENIE_INCLUDES_DIR}/GHEP + ${GENIE_INCLUDES_DIR}/Ntuple + ${GENIE_INCLUDES_DIR}/ReWeight + ${GENIE_INCLUDES_DIR}/Apps + ${GENIE_INCLUDES_DIR}/FluxDrivers + ${GENIE_INCLUDES_DIR}/EVGDrivers + ${LHAPDF_INC} + ${LIBXML2_INC} + ${LOG4CPP_INC}) +else() + LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES + ${GENIE_INCLUDES_DIR} + $ENV{GENIE_REWEIGHT}/src + ${GSL_INC} + ${LHAPDF_INC} + ${LIBXML2_INC} + ${LOG4CPP_INC}) + LIST(APPEND EXTRA_LINK_DIRS $ENV{GSL_LIB}) + LIST(APPEND EXTRA_LIBS gsl gslcblas) +endif() -LIST(APPEND EXTRA_LIBS LHAPDF xml2 log4cpp) +SAYVARS() if(USE_PYTHIA8) set(NEED_PYTHIA8 TRUE) set(NEED_ROOTPYTHIA8 TRUE) else() set(NEED_PYTHIA6 TRUE) set(NEED_ROOTPYTHIA6 TRUE) endif() set(NEED_ROOTEVEGEN TRUE) SET(USE_GENIE TRUE CACHE BOOL "Whether to enable GENIE (reweight) support. Requires external libraries. " FORCE) diff --git a/parameters/config.xml b/parameters/config.xml index b4a01c0..f093aee 100644 --- a/parameters/config.xml +++ b/parameters/config.xml @@ -1,215 +1,219 @@ + + + + diff --git a/src/FitBase/MeasurementBase.h b/src/FitBase/MeasurementBase.h index 00dddc7..87eae43 100644 --- a/src/FitBase/MeasurementBase.h +++ b/src/FitBase/MeasurementBase.h @@ -1,362 +1,366 @@ // Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ #ifndef INPUTHANDLER_H_SEEN #define INPUTHANDLER_H_SEEN /*! * \addtogroup FitBase * @{ */ // C Includes #include #include #include #include #include #include #include #include #include #include // ROOT includes #include #include #include #include #include #include #include #include #include #include // External data fit includes #include "FitEvent.h" #include "FitUtils.h" #include "GeneralUtils.h" #include "PlotUtils.h" #include "StatUtils.h" #include "InputFactory.h" #include "FitWeight.h" #include "TMultiDimFit.h" #ifdef __GENIE_ENABLED__ +#ifdef GENIE_PRE_R3 #include "Conventions/Units.h" +#else +#include "Framework/Conventions/Units.h" +#endif #endif #include "EventManager.h" #include "TObject.h" #include "InputHandler.h" #include "NuisConfig.h" #include "NuisKey.h" #include "SampleSettings.h" #include "StackBase.h" #include "StandardStacks.h" /// Enumerations to help with extra plot functions enum extraplotflags { kCMD_Reset = 0, kCMD_Fill, kCMD_Scale, kCMD_Norm, kCMD_Write, kCMD_Error, kCMD_extraplotflags }; enum MeasurementSpeciesClass { kSingleSpeciesMeasurement = 0, kNumuWithWrongSignMeasurement, kNueWithWrongSignMeasurement, kFourSpeciesMeasurement, }; /// InputHandler Class /// /// Inherits from Measurement base to handle whatever input is throwna t the /// fitter automatically. /// All functions here handle how files are read in, converted to custom formats /// and reconfigures are called. /// Used generally for the MC inputs. //! 2nd level experiment class that handles converting MC into a common format //! and calling reconfigure class MeasurementBase { public: /* Constructor/Destructors */ //! Default Constructor. Set everything to NULL MeasurementBase(); //! Default virtual destructor virtual ~MeasurementBase(void); virtual void InitialSetup(void) {}; /* Reconfigure Functions */ //! Function called if MC tuning dials haven't been changed and all we want to //! do is update the normalisation. virtual void Renormalise(void); //! Call reconfigure only looping over signal events to save time. virtual void ReconfigureFast(void); virtual void FillHistograms(double weight); //! Call reconfigure looping over all MC events including background virtual void Reconfigure(void); // virtual TH2D GetCovarMatrix(void) = 0; virtual double GetLikelihood(void) { return 0.0; }; virtual int GetNDOF(void) { return 0; }; virtual void ThrowCovariance(void) = 0; virtual void ThrowDataToy(void) = 0; virtual void SetFakeDataValues(std::string fkdt) = 0; //! Get the total integrated flux between this samples energy range virtual double TotalIntegratedFlux(std::string intOpt = "width", double low = -9999.9, double high = -9999.9); //! Get the predicted event rate for this sample virtual double PredictedEventRate(std::string intOpt = "width", double low = -9999.9, double high = -9999.9); virtual SampleSettings LoadSampleSettings(nuiskey samplekey); virtual SampleSettings LoadSampleSettings(std::string name, std::string input, std::string type); virtual void FinaliseSampleSettings(); virtual void FinaliseMeasurement(); virtual void ProcessExtraHistograms(int cmd, MeasurementVariableBox* vars, double weight = 1.0); virtual void FillExtraHistograms(MeasurementVariableBox* vars, double weight = 1.0); virtual void ScaleExtraHistograms(MeasurementVariableBox* vars); virtual void ResetExtraHistograms(); virtual void NormExtraHistograms(MeasurementVariableBox* vars, double norm = 1.0); virtual void WriteExtraHistograms(); virtual MeasurementVariableBox* CreateBox() {return new MeasurementVariableBox();}; int GetPassed() { int signalSize = 0; return signalSize; } int GetTotal() { return fNEvents; } /* Reconfigure LOOP */ // All these should be virtual ///! Reset Histograms (Handled at Measurement Stage) virtual void ResetAll(void) = 0; ///! Fill the event variables for this sample (Handled in each inherited /// sample) virtual void FillEventVariables(FitEvent* event) { (void)event; }; ///! Check whether this event is signle (Handled in each inherited sample) virtual bool isSignal(FitEvent* event) { (void)event; return false; }; ///! Fill the histogram for this event using fXVar and fYVar (Handled in each /// inherited sample) virtual void FillHistograms(void) {}; ///! Convert event rates to whatever distributions you need. virtual void ConvertEventRates(void); ///! Call scale events after the plots have been filled at the end of /// reconfigure. virtual void ScaleEvents(void) {}; ///! Apply the scale factor at the end of reconfigure. virtual void ApplyNormScale(double norm) { (void)norm; }; ///! Save Histograms virtual void Write(std::string drawOpt = "") = 0; virtual MeasurementVariableBox* FillVariableBox(FitEvent* event); virtual MeasurementVariableBox* GetBox(); void FillHistogramsFromBox(MeasurementVariableBox* var, double weight); /* Histogram Access Functions */ ///! Virtual function to get data histogram virtual std::vector GetDataList(void) = 0; ///! Virtual function to get MC histogram virtual std::vector GetMCList(void) = 0; virtual std::vector GetFineList(void) = 0; virtual std::vector GetMaskList(void) = 0; ///! Return flux histograms in a vector virtual std::vector GetFluxList(void); virtual std::vector GetEventRateList(void); virtual std::vector GetXSecList(void); virtual TH1D* GetEventHistogram() { return fInput->GetEventHistogram(); }; virtual TH1D* GetXSecHistogram() { return fInput->GetXSecHistogram(); }; virtual TH1D* GetFluxHistogram() { return fInput->GetFluxHistogram(); }; ///! Return input for this sample InputHandlerBase* GetInput(void); std::string GetName(void) { return fName; }; double GetScaleFactor(void) { return fScaleFactor; }; double GetXVar(void) { return fXVar; }; double GetYVar(void) { return fYVar; }; double GetZVar(void) { return fZVar; }; double GetMode(void) { return this->Mode; }; double GetEnu(void) { return this->Enu; }; void SetupInputs(std::string inputfile); int GetInputID(void); std::string GetInputFileName() { return fInputFileName; }; void SetSignal(bool sig); void SetSignal(FitEvent* evt); void SetWeight(double wght); void SetMode(int md); void SetNoData(bool isTrue = true) { fNoData = isTrue; }; inline void SetXVar(double xvar) { fXVar = xvar; }; inline void SetYVar(double yvar) { fYVar = yvar; }; inline void SetZVar(double zvar) { fZVar = zvar; }; virtual std::vector GetSubSamples() { return std::vector(1, this); } void SetAutoProcessTH1(TH1* hist, int c1 = -1, int c2 = -1, int c3 = -1, int c4 = -1, int c5 = -1); void SetAutoProcess(TH1* hist, int c1 = -1, int c2 = -1, int c3 = -1, int c4 = -1, int c5 = -1); void SetAutoProcess(TGraph* g, int c1 = -1, int c2 = -1, int c3 = -1, int c4 = -1, int c5 = -1); void SetAutoProcess(TF1* f, int c1 = -1, int c2 = -1, int c3 = -1, int c4 = -1, int c5 = -1); void SetAutoProcess(StackBase* hist, int c1 = -1, int c2 = -1, int c3 = -1, int c4 = -1, int c5 = -1); void SetAutoProcessTH1(StackBase* hist, int c1 = -1, int c2 = -1, int c3 = -1, int c4 = -1, int c5 = -1); void AutoFillExtraTH1(); void AutoResetExtraTH1(); void AutoScaleExtraTH1(); void AutoNormExtraTH1(double norm); void AutoWriteExtraTH1(); // functions that need to be added. // - Initial Check // - Check Target/Beam loop. // - Check flux shape if suggested one given. // - Return MeasurementList (returns ) protected: // Minimum and maximum energies double Enu; //!< Neutrino Energy double EnuMin; //!< Minimum incoming particle energy of events to include double EnuMax; //!< Maximum incoming particle energy of events to include BaseFitEvt* signal_event; FitEvent* cust_event; FitWeight* fRW; //!< Pointer to the rw engine InputHandlerBase* fInput; //!< Instance of the input handler std::string fName; //!< Name of the sample int fEventType; double fBeamDistance; //!< Incoming Particle flight distance (for oscillation //! analysis) double fScaleFactor; //!< fScaleFactor applied to events to convert from //! eventrate to final distribution double fCurrentNorm; //!< current normalisation factor applied if fit is "FREE" bool fMCFilled; //!< flag whether MC plots have been filled (For //! ApplyNormalisation) bool fNoData; //!< flag whether data plots do not exist (for ratios) bool fIsNoWidth; ///< Flag : Don't scale by bin width // TEMP OBJECTS TO HANDLE MERGE double fXVar, fYVar, fZVar, Mode, Weight; bool Signal; int ievt; int fNEvents; double Enu_rec, ThetaMu, CosThetaMu; InputUtils::InputType fInputType; std::string fInputFileName; TH1D* fFluxHist; TH1D* fEventHist; MeasurementSpeciesClass fMeasurementSpeciesType; SampleSettings fSettings; MeasurementVariableBox* fEventVariables; std::map > fExtraTH1s; int NSignal; // std::map fExtaStacks; bool fIsJoint; double fNPOT, fFluxIntegralOverride, fTargetVolume, fTargetMaterialDensity; double fEvtRateScaleFactor; }; // Class TypeDefs typedef std::list::const_iterator MeasListConstIter; typedef std::list::iterator MeasListIter; typedef std::vector::const_iterator MeasVectConstIter; typedef std::vector::iterator MeasVectIter; /*! @} */ #endif diff --git a/src/InputHandler/BaseFitEvt.h b/src/InputHandler/BaseFitEvt.h index 84316aa..84e3d88 100644 --- a/src/InputHandler/BaseFitEvt.h +++ b/src/InputHandler/BaseFitEvt.h @@ -1,141 +1,147 @@ // Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ #ifndef FITEVENTBASE_H_SEEN #define FITEVENTBASE_H_SEEN /*! * \addtogroup InputHandler * @{ */ #ifdef __NEUT_ENABLED__ #include "neutpart.h" #include "neutvect.h" #endif #ifdef __NUWRO_ENABLED__ #ifdef __USE_NUWRO_SRW_EVENTS__ #include "NuwroReWeightSimpleEvent.h" #else #include "event1.h" #endif #endif #ifdef __GENIE_ENABLED__ +#ifdef GENIE_PRE_R3 #include "EVGCore/EventRecord.h" #include "GHEP/GHepRecord.h" #include "Ntuple/NtpMCEventRecord.h" +#else +#include "Framework/EventGen/EventRecord.h" +#include "Framework/GHEP/GHepRecord.h" +#include "Framework/Ntuple/NtpMCEventRecord.h" +#endif using namespace genie; #endif #ifdef __NUANCE_ENABLED__ #include "NuanceEvent.h" #endif #include "StdHepEvt.h" #include "SplineReader.h" #include "InputTypes.h" #include "GeneratorInfoBase.h" /// Base Event Class used to store just the generator event pointers class BaseFitEvt { public: /// Base Constructor BaseFitEvt(void); ~BaseFitEvt(); /// Copy base fit event pointers BaseFitEvt(const BaseFitEvt* obj); BaseFitEvt(BaseFitEvt const &); BaseFitEvt operator=(BaseFitEvt const &); /// Reset weight to 1.0 void ResetWeight(); /// Return combined weight for this event double GetWeight(); /// Manually set event type inline void SetType(int type){fType = type;}; // Global Event Variables/Weights int Mode; ///< True interaction mode double probe_E; ///< True probe energy double probe_pdg; // Weighting Info double Weight; ///< Total Weight For Event double InputWeight; ///< Input Starting Weight (used for GiBUU) double RWWeight; ///< Saved RW from FitWeight double CustomWeight; ///< Extra custom weight that samples can set double SavedRWWeight; ///< Saved RW value for FitEvents double CustomWeightArray[6]; ///< For custom tuning using arrays, e.g. NOvA MINERvA WS // Spline Info Coefficients and Readers float* fSplineCoeff; ///< ND Array of Spline Coefficients SplineReader* fSplineRead; ///< Spline Interpretter // Generator Info GeneratorInfoBase* fGenInfo; ///< Generator Variable Box UInt_t fType; ///< Generator Event Type #ifdef __NEUT_ENABLED__ /// Setup Event Reading from NEUT Event void SetNeutVect(NeutVect* v); NeutVect* fNeutVect; ///< Pointer to Neut Vector #endif #ifdef __NUWRO_ENABLED__ #ifdef __USE_NUWRO_SRW_EVENTS__ SRW::SRWEvent fNuwroSRWEvent; ///< Pointer to Nuwro event params * fNuwroParams; #endif event* fNuwroEvent; ///< Pointer to Nuwro event #endif #ifdef __GENIE_ENABLED__ /// Setup Event Reading from GENIE Event void SetGenieEvent(NtpMCEventRecord* ntpl); NtpMCEventRecord* genie_event; ///< Pointer to NTuple Genie Event GHepRecord* genie_record; ///< Pointer to actually accessible Genie Record #endif #ifdef __NUANCE_ENABLED__ /// Setup Event Reading from NUANCE Event void SetNuanceEvent(NuanceEvent* e); NuanceEvent* nuance_event; ///< Pointer to Nuance reader #endif #ifdef __GiBUU_ENABLED__ GiBUUStdHepReader* GiRead; ///< Pointer to GiBUU reader #endif /// Setup Event Type to FitEvent void SetInputFitEvent(); /// Setup Event Type to FitSpline void SetInputFitSpline(); /// Setup Event Type to HepMC void SetInputHepMC(); }; /*! @} */ #endif diff --git a/src/InputHandler/GENIEInputHandler.h b/src/InputHandler/GENIEInputHandler.h index d2a1d18..16a7092 100644 --- a/src/InputHandler/GENIEInputHandler.h +++ b/src/InputHandler/GENIEInputHandler.h @@ -1,119 +1,130 @@ // Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ #ifndef GENIEINPUTHANDLER_H #define GENIEINPUTHANDLER_H /*! * \addtogroup InputHandler * @{ */ #ifdef __GENIE_ENABLED__ #include "InputHandler.h" #include "InputUtils.h" #include "PlotUtils.h" +#ifdef GENIE_PRE_R3 #include "GHEP/GHepParticle.h" #include "PDG/PDGUtils.h" #include "GHEP/GHepUtils.h" #include "Conventions/Units.h" #include "EVGCore/EventRecord.h" #include "GHEP/GHepRecord.h" #include "Ntuple/NtpMCEventRecord.h" +#else +#include "Framework/GHEP/GHepParticle.h" +#include "Framework/ParticleData/PDGUtils.h" +#include "Framework/GHEP/GHepUtils.h" +#include "Framework/Conventions/Units.h" +#include "Framework/EventGen/EventRecord.h" +#include "Framework/GHEP/GHepRecord.h" +#include "Framework/Ntuple/NtpMCEventRecord.h" +#endif + using namespace genie; /// GENIE Generator Container to save extra particle status codes. class GENIEGeneratorInfo : public GeneratorInfoBase { public: GENIEGeneratorInfo() {}; virtual ~GENIEGeneratorInfo(); /// Assigns information to branches void AddBranchesToTree(TTree* tn); /// Setup reading information from branches void SetBranchesFromTree(TTree* tn); /// Allocate any dynamic arrays for a new particle stack size void AllocateParticleStack(int stacksize); /// Clear any dynamic arrays void DeallocateParticleStack(); /// Read extra genie information from the event void FillGeneratorInfo(NtpMCEventRecord* ntpl); /// Reset extra information to default/empty values void Reset(); int kMaxParticles; ///< Number of particles in stack int* fGenieParticlePDGs; ///< GENIE Particle PDGs (example) }; /// Main GENIE InputHandler class GENIEInputHandler : public InputHandlerBase { public: /// Standard constructor given a name and input files GENIEInputHandler(std::string const& handle, std::string const& rawinputs); virtual ~GENIEInputHandler(); /// Create a TTree Cache to speed up file read void CreateCache(); /// Remove TTree Cache to save memory void RemoveCache(); /// Returns a NUISANCE format event from the GENIE TTree. If !lightweight /// then CalcNUISANCEKinematics() is called to convert the GENIE event into /// a standard NUISANCE format. FitEvent* GetNuisanceEvent(const UInt_t entry, const bool lightweight = false); /// Converts GENIE event into standard NUISANCE FitEvent by looping over all /// particles in the event and adding them to stack in fNUISANCEEvent. void CalcNUISANCEKinematics(); /// Placeholder for GENIE related event printing. void Print(); /// Converts GENIE particle status codes into NUISANCE status codes. int GetGENIEParticleStatus(genie::GHepParticle* part, int mode = 0); /// Converts GENIE event reaction codes into NUISANCE reaction codes. int ConvertGENIEReactionCode(GHepRecord* gheprec); GHepRecord* fGenieGHep; ///< Pointer to actual event record NtpMCEventRecord* fGenieNtpl; ///< Ntpl Wrapper Class TChain* fGENIETree; ///< Main GENIE Event TTree bool fSaveExtra; ///< Flag to save Extra GENIE info into Nuisance Event GENIEGeneratorInfo* fGenieInfo; ///< Extra GENIE Generator Info Writer bool fNOvAWeights; ///< Flag to save nova weights or not // Extra weights from Jeremy for NOvA weights double MAQEw; double NonResw; double RPAQEw; double RPARESw; double MECw; double DISw; double NOVAw; }; /*! @} */ #endif #endif diff --git a/src/Reweight/GENIEWeightEngine.cxx b/src/Reweight/GENIEWeightEngine.cxx index 21eb1b0..31afebb 100644 --- a/src/Reweight/GENIEWeightEngine.cxx +++ b/src/Reweight/GENIEWeightEngine.cxx @@ -1,267 +1,277 @@ #include "GENIEWeightEngine.h" #ifdef __GENIE_EMP_MECRW_ENABLED #include "ReWeight/GReWeightXSecEmpiricalMEC.h" #endif GENIEWeightEngine::GENIEWeightEngine(std::string name) { #ifdef __GENIE_ENABLED__ // Setup the NEUT Reweight engien fCalcName = name; LOG(FIT) << "Setting up GENIE RW : " << fCalcName << std::endl; // Create RW Engine suppressing cout StopTalking(); fGenieRW = new genie::rew::GReWeight(); // Get List of Vetos (Just for debugging) std::string rw_engine_list = FitPar::Config().GetParS("FitWeight_fGenieRW_veto"); bool xsec_ncel = rw_engine_list.find("xsec_ncel") == std::string::npos; bool xsec_ccqe = rw_engine_list.find("xsec_ccqe") == std::string::npos; bool xsec_coh = rw_engine_list.find("xsec_coh") == std::string::npos; bool xsec_nnres = rw_engine_list.find("xsec_nonresbkg") == std::string::npos; bool xsec_nudis = rw_engine_list.find("nuclear_dis") == std::string::npos; bool xsec_resdec = rw_engine_list.find("hadro_res_decay") == std::string::npos; bool xsec_fzone = rw_engine_list.find("hadro_intranuke") == std::string::npos; bool xsec_intra = rw_engine_list.find("hadro_fzone") == std::string::npos; bool xsec_agky = rw_engine_list.find("hadro_agky") == std::string::npos; bool xsec_qevec = rw_engine_list.find("xsec_ccqe_vec") == std::string::npos; bool xsec_dis = rw_engine_list.find("xsec_dis") == std::string::npos; bool xsec_nc = rw_engine_list.find("xsec_nc") == std::string::npos; bool xsec_ccres = rw_engine_list.find("xsec_ccres") == std::string::npos; bool xsec_ncres = rw_engine_list.find("xsec_ncres") == std::string::npos; bool xsec_nucqe = rw_engine_list.find("nuclear_qe") == std::string::npos; bool xsec_qeaxial = rw_engine_list.find("xsec_ccqe_axial") == std::string::npos; #ifdef __GENIE_EMP_MECRW_ENABLED bool xsec_empMEC = rw_engine_list.find("xsec_empMEC") == std::string::npos; #endif +#ifndef GENIE_PRE_R3 + genie::RunOpt* grunopt = genie::RunOpt::Instance(); + grunopt->EnableBareXSecPreCalc(true); + grunopt->SetEventGeneratorList(Config::GetParS("GENIEEventGeneratorList")); + grunopt->SetTuneName(Config::GetParS("GENIETune")); + grunopt->BuildTune(); + std::string genv = std::string(getenv("GENIE")) + "/config/Messenger_laconic.xml"; + genie::utils::app_init::MesgThresholds(genv); +#endif + // Now actually add the RW Calcs if (xsec_ncel) fGenieRW->AdoptWghtCalc("xsec_ncel", new genie::rew::GReWeightNuXSecNCEL); if (xsec_ccqe) { fGenieRW->AdoptWghtCalc("xsec_ccqe", new genie::rew::GReWeightNuXSecCCQE); // (dynamic_cast (fGenieRW->WghtCalc("xsec_ccqe"))) // ->SetXSecModel( FitPar::Config().GetParS("GENIEXSecModelCCQE") ); } #ifdef __GENIE_EMP_MECRW_ENABLED if (xsec_empMEC) { fGenieRW->AdoptWghtCalc("xsec_empMEC", new genie::rew::GReWeightXSecEmpiricalMEC); } #endif if (xsec_coh) { fGenieRW->AdoptWghtCalc("xsec_coh", new genie::rew::GReWeightNuXSecCOH()); // (dynamic_cast (fGenieRW->WghtCalc("xsec_coh"))) // ->SetXSecModel( FitPar::Config().GetParS("GENIEXSecModelCOH") ); } if (xsec_nnres) fGenieRW->AdoptWghtCalc("xsec_nonresbkg", new genie::rew::GReWeightNonResonanceBkg); if (xsec_nudis) fGenieRW->AdoptWghtCalc("nuclear_dis", new genie::rew::GReWeightDISNuclMod); if (xsec_resdec) fGenieRW->AdoptWghtCalc("hadro_res_decay", new genie::rew::GReWeightResonanceDecay); if (xsec_fzone) fGenieRW->AdoptWghtCalc("hadro_fzone", new genie::rew::GReWeightFZone); if (xsec_intra) fGenieRW->AdoptWghtCalc("hadro_intranuke", new genie::rew::GReWeightINuke); if (xsec_agky) fGenieRW->AdoptWghtCalc("hadro_agky", new genie::rew::GReWeightAGKY); if (xsec_qevec) fGenieRW->AdoptWghtCalc("xsec_ccqe_vec", new genie::rew::GReWeightNuXSecCCQEvec); #if __GENIE_VERSION__ >= 212 if (xsec_qeaxial) fGenieRW->AdoptWghtCalc("xsec_ccqe_axial", new genie::rew::GReWeightNuXSecCCQEaxial); #endif if (xsec_dis) fGenieRW->AdoptWghtCalc("xsec_dis", new genie::rew::GReWeightNuXSecDIS); if (xsec_nc) fGenieRW->AdoptWghtCalc("xsec_nc", new genie::rew::GReWeightNuXSecNC); if (xsec_ccres) { #if __GENIE_VERSION__ < 213 fGenieRW->AdoptWghtCalc("xsec_ccres", new genie::rew::GReWeightNuXSecCCRES); #else fGenieRW->AdoptWghtCalc( "xsec_ccres", new genie::rew::GReWeightNuXSecCCRES( FitPar::Config().GetParS("GENIEXSecModelCCRES"), "Default")); #endif } if (xsec_ncres) fGenieRW->AdoptWghtCalc("xsec_ncres", new genie::rew::GReWeightNuXSecNCRES); if (xsec_nucqe) fGenieRW->AdoptWghtCalc("nuclear_qe", new genie::rew::GReWeightFGM); if (xsec_ccqe) { GReWeightNuXSecCCQE *rwccqe = dynamic_cast(fGenieRW->WghtCalc("xsec_ccqe")); rwccqe->SetMode(GReWeightNuXSecCCQE::kModeMa); } if (xsec_ccres) { // Default to include shape and normalization changes for CCRES (can be // changed downstream if desired) GReWeightNuXSecCCRES *rwccres = dynamic_cast(fGenieRW->WghtCalc("xsec_ccres")); std::string marestype = FitPar::Config().GetParS("GENIEWeightEngine_CCRESMode"); if (!marestype.compare("kModeNormAndMaMvShape")) { rwccres->SetMode(GReWeightNuXSecCCRES::kModeNormAndMaMvShape); } else if (!marestype.compare("kModeMaMv")) { rwccres->SetMode(GReWeightNuXSecCCRES::kModeMaMv); } else { THROW("Unkown MARES Mode in GENIE Weight Engine : " << marestype); } } if (xsec_ncres) { // Default to include shape and normalization changes for NCRES (can be // changed downstream if desired) GReWeightNuXSecNCRES *rwncres = dynamic_cast(fGenieRW->WghtCalc("xsec_ncres")); rwncres->SetMode(GReWeightNuXSecNCRES::kModeMaMv); } if (xsec_dis) { // Default to include shape and normalization changes for DIS (can be // changed downstream if desired) GReWeightNuXSecDIS *rwdis = dynamic_cast(fGenieRW->WghtCalc("xsec_dis")); rwdis->SetMode(GReWeightNuXSecDIS::kModeABCV12u); // Set Abs Twk Config fIsAbsTwk = (FitPar::Config().GetParB("setabstwk")); } // allow cout again StartTalking(); #else ERR(FTL) << "GENIE RW NOT ENABLED" << std::endl; #endif }; void GENIEWeightEngine::IncludeDial(std::string name, double startval) { #ifdef __GENIE_ENABLED__ // Get First enum int nuisenum = Reweight::ConvDial(name, kGENIE); // Setup Maps fEnumIndex[nuisenum]; // = std::vector(0); fNameIndex[name]; // = std::vector(0); // Split by commas std::vector allnames = GeneralUtils::ParseToStr(name, ","); for (uint i = 0; i < allnames.size(); i++) { std::string singlename = allnames[i]; // Get RW genie::rew::GSyst_t rwsyst = GSyst::FromString(singlename); // Fill Maps int index = fValues.size(); fValues.push_back(0.0); fGENIESysts.push_back(rwsyst); // Initialize dial std::cout << "Registering " << singlename << " from " << name << std::endl; fGenieRW->Systematics().Init(fGENIESysts[index]); // If Absolute if (fIsAbsTwk) { GSystUncertainty::Instance()->SetUncertainty(rwsyst, 1.0, 1.0); } // Setup index fEnumIndex[nuisenum].push_back(index); fNameIndex[name].push_back(index); } // Set Value if given if (startval != -999.9) { SetDialValue(nuisenum, startval); } #endif }; void GENIEWeightEngine::SetDialValue(int nuisenum, double val) { #ifdef __GENIE_ENABLED__ std::vector indices = fEnumIndex[nuisenum]; for (uint i = 0; i < indices.size(); i++) { fValues[indices[i]] = val; fGenieRW->Systematics().Set(fGENIESysts[indices[i]], val); } #endif } void GENIEWeightEngine::SetDialValue(std::string name, double val) { #ifdef __GENIE_ENABLED__ std::vector indices = fNameIndex[name]; for (uint i = 0; i < indices.size(); i++) { fValues[indices[i]] = val; fGenieRW->Systematics().Set(fGENIESysts[indices[i]], val); } #endif } void GENIEWeightEngine::Reconfigure(bool silent) { #ifdef __GENIE_ENABLED__ // Hush now... if (silent) StopTalking(); // Reconf fGenieRW->Reconfigure(); fGenieRW->Print(); // Shout again if (silent) StartTalking(); #endif } double GENIEWeightEngine::CalcWeight(BaseFitEvt *evt) { double rw_weight = 1.0; #ifdef __GENIE_ENABLED__ // Make nom weight if (!evt) { THROW("evt not found : " << evt); } // Skip Non GENIE if (evt->fType != kGENIE) return 1.0; if (!(evt->genie_event)) { THROW("evt->genie_event not found!" << evt->genie_event); } if (!(evt->genie_event->event)) { THROW("evt->genie_event->event GHepRecord not found!" << (evt->genie_event->event)); } if (!fGenieRW) { THROW("GENIE RW Not Found!" << fGenieRW); } rw_weight = fGenieRW->CalcWeight(*(evt->genie_event->event)); // std::cout << "Returning GENIE Weight for electron scattering = " << // rw_weight << std::endl; #endif // Return rw_weight return rw_weight; } diff --git a/src/Reweight/GENIEWeightEngine.h b/src/Reweight/GENIEWeightEngine.h index b0d3d03..e0a7f7e 100644 --- a/src/Reweight/GENIEWeightEngine.h +++ b/src/Reweight/GENIEWeightEngine.h @@ -1,64 +1,88 @@ #ifndef WEIGHT_ENGINE_GENIE_H #define WEIGHT_ENGINE_GENIE_H #include "FitLogger.h" #ifdef __GENIE_ENABLED__ -#include "EVGCore/EventRecord.h" +#ifdef GENIE_PRE_R3 #include "EVGCore/EventRecord.h" #include "GHEP/GHepRecord.h" -#include "GSyst.h" -#include "GSystUncertainty.h" +#include "ReWeight/GSyst.h" +#include "ReWeight/GSystUncertainty.h" #include "Ntuple/NtpMCEventRecord.h" #include "ReWeight/GReWeight.h" #include "ReWeight/GReWeightAGKY.h" #include "ReWeight/GReWeightDISNuclMod.h" #include "ReWeight/GReWeightFGM.h" #include "ReWeight/GReWeightFZone.h" #include "ReWeight/GReWeightINuke.h" #include "ReWeight/GReWeightNonResonanceBkg.h" #include "ReWeight/GReWeightNuXSecCCQE.h" #include "ReWeight/GReWeightNuXSecCCQEvec.h" #include "ReWeight/GReWeightNuXSecCCRES.h" #include "ReWeight/GReWeightNuXSecCOH.h" #include "ReWeight/GReWeightNuXSecDIS.h" #include "ReWeight/GReWeightNuXSecNC.h" #include "ReWeight/GReWeightNuXSecNCEL.h" #include "ReWeight/GReWeightNuXSecNCRES.h" #include "ReWeight/GReWeightResonanceDecay.h" - #if __GENIE_VERSION__ >= 212 #include "ReWeight/GReWeightNuXSecCCQEaxial.h" #endif - +#else +#include "Framework/EventGen/EventRecord.h" +#include "Framework/GHEP/GHepRecord.h" +#include "Framework/GHEP/GHepParticle.h" +#include "Framework/Ntuple/NtpMCEventRecord.h" +#include "Framework/Utils/RunOpt.h" +#include "Framework/Utils/AppInit.h" +#include "RwFramework/GSyst.h" +#include "RwFramework/GSystUncertainty.h" +#include "RwFramework/GReWeight.h" +#include "RwCalculators/GReWeightAGKY.h" +#include "RwCalculators/GReWeightDISNuclMod.h" +#include "RwCalculators/GReWeightFGM.h" +#include "RwCalculators/GReWeightFZone.h" +#include "RwCalculators/GReWeightINuke.h" +#include "RwCalculators/GReWeightNonResonanceBkg.h" +#include "RwCalculators/GReWeightNuXSecCCQE.h" +#include "RwCalculators/GReWeightNuXSecCCQEvec.h" +#include "RwCalculators/GReWeightNuXSecCCRES.h" +#include "RwCalculators/GReWeightNuXSecCOH.h" +#include "RwCalculators/GReWeightNuXSecDIS.h" +#include "RwCalculators/GReWeightNuXSecNC.h" +#include "RwCalculators/GReWeightNuXSecNCEL.h" +#include "RwCalculators/GReWeightNuXSecNCRES.h" +#include "RwCalculators/GReWeightResonanceDecay.h" +#include "RwCalculators/GReWeightNuXSecCCQEaxial.h" +#endif using namespace genie; using namespace genie::rew; #endif - #include "GeneratorUtils.h" #include "WeightEngineBase.h" #include "FitWeight.h" class GENIEWeightEngine : public WeightEngineBase { public: GENIEWeightEngine(std::string name); ~GENIEWeightEngine() {}; void IncludeDial(std::string name, double startval); void SetDialValue(int rwenum, double val); void SetDialValue(std::string name, double val); void Reconfigure(bool silent = false); double CalcWeight(BaseFitEvt* evt); inline bool NeedsEventReWeight() { return true; }; #ifdef __GENIE_ENABLED__ std::vector fGENIESysts; genie::rew::GReWeight* fGenieRW; //!< Genie RW Object #endif }; #endif diff --git a/src/Reweight/MINERvAWeightCalcs.h b/src/Reweight/MINERvAWeightCalcs.h index 79e4654..965010d 100644 --- a/src/Reweight/MINERvAWeightCalcs.h +++ b/src/Reweight/MINERvAWeightCalcs.h @@ -1,148 +1,158 @@ #ifndef MINERVA_WEIGHT_CALCS #define MINERVA_WEIGHT_CALCS #include #ifdef __MINERVA_RW_ENABLED__ #ifdef __GENIE_ENABLED__ +#ifdef GENIE_PRE_R3 #include "Conventions/Units.h" #include "EVGCore/EventRecord.h" -#include "FitEvent.h" - #include "GHEP/GHepParticle.h" #include "GHEP/GHepRecord.h" #include "GHEP/GHepUtils.h" -#include "GeneralUtils.h" -#include "NUISANCESyst.h" -#include "NUISANCEWeightCalcs.h" #include "Ntuple/NtpMCEventRecord.h" #include "PDG/PDGUtils.h" +#else +#include "Framework/Conventions/Units.h" +#include "Framework/EventGen/EventRecord.h" +#include "Framework/GHEP/GHepParticle.h" +#include "Framework/GHEP/GHepRecord.h" +#include "Framework/GHEP/GHepUtils.h" +#include "Framework/Ntuple/NtpMCEventRecord.h" +#include "Framework/ParticleData/PDGUtils.h" +#endif + +#include "NUISANCEWeightCalcs.h" +#include "GeneralUtils.h" +#include "NUISANCESyst.h" +#include "FitEvent.h" #include "WeightUtils.h" #include "weightRPA.h" using namespace genie; class BaseFitEvt; namespace nuisance { namespace reweight { // MEC Dials class MINERvAReWeight_QE : public NUISANCEWeightCalc { public: MINERvAReWeight_QE(); virtual ~MINERvAReWeight_QE(); double CalcWeight(BaseFitEvt* evt); void SetDialValue(std::string name, double val); void SetDialValue(int rwenum, double val); bool IsHandled(int rwenum); double fTwk_NormCCQE; double fCur_NormCCQE; double fDef_NormCCQE; bool fTweaked; }; // MEC Dials class MINERvAReWeight_MEC : public NUISANCEWeightCalc { public: MINERvAReWeight_MEC(); virtual ~MINERvAReWeight_MEC(); double CalcWeight(BaseFitEvt* evt); void SetDialValue(std::string name, double val); void SetDialValue(int rwenum, double val); bool IsHandled(int rwenum); double fTwk_NormCCMEC; double fCur_NormCCMEC; double fDef_NormCCMEC; bool fTweaked; }; // RES Dials class MINERvAReWeight_RES : public NUISANCEWeightCalc { public: MINERvAReWeight_RES(); virtual ~MINERvAReWeight_RES(); double CalcWeight(BaseFitEvt* evt); void SetDialValue(std::string name, double val); void SetDialValue(int rwenum, double val); bool IsHandled(int rwenum); double fTwk_NormCCRES; double fCur_NormCCRES; double fDef_NormCCRES; bool fTweaked; }; /// RPA Weight Calculator that applies RPA systematics /// to GENIE events. GENIE EVENTS ONLY! class RikRPA : public NUISANCEWeightCalc { public: RikRPA(); ~RikRPA(); double CalcWeight(BaseFitEvt* evt); void SetDialValue(std::string name, double val); void SetDialValue(int rwenum, double val); bool IsHandled(int rwenum); void SetupRPACalculator(int calcenum); int GetRPACalcEnum(int bpdg, int tpdg); bool fApplyDial_RPACorrection; double fTwkDial_RPALowQ2; double fDefDial_RPALowQ2; double fCurDial_RPALowQ2; double fErrDial_RPALowQ2; double fTwkDial_RPAHighQ2; double fDefDial_RPAHighQ2; double fCurDial_RPAHighQ2; double fErrDial_RPAHighQ2; bool fApplyDial_RESRPACorrection; double fTwkDial_RESRPALowQ2; double fDefDial_RESRPALowQ2; double fCurDial_RESRPALowQ2; double fErrDial_RESRPALowQ2; double fTwkDial_RESRPAHighQ2; double fDefDial_RESRPAHighQ2; double fCurDial_RESRPAHighQ2; double fErrDial_RESRPAHighQ2; double* fEventWeights; bool fTweaked; const static int kMaxCalculators = 10; enum rpacalcenums { kNuMuC12, kNuMuO16, kNuMuAr40, kNuMuCa40, kNuMuFe56, kNuMuBarC12, kNuMuBarO16, kNuMuBarAr40, kNuMuBarCa40, kNuMuBarFe56 }; weightRPA* fRPACalculators[kMaxCalculators]; }; }; // namespace reweight }; // namespace nuisance #endif // __GENIE_ENABLED__ #endif //__MINERVA_RW_ENABLED__ #endif diff --git a/src/Reweight/WeightUtils.cxx b/src/Reweight/WeightUtils.cxx index 1a1bd77..6beb1dc 100644 --- a/src/Reweight/WeightUtils.cxx +++ b/src/Reweight/WeightUtils.cxx @@ -1,614 +1,637 @@ #include "WeightUtils.h" #include "FitLogger.h" #ifdef __T2KREW_ENABLED__ #include "T2KGenieReWeight.h" #include "T2KNIWGReWeight.h" #include "T2KNIWGUtils.h" #include "T2KNeutReWeight.h" #include "T2KNeutUtils.h" #include "T2KReWeight.h" using namespace t2krew; #endif #ifdef __NIWG_ENABLED__ #include "NIWGReWeight.h" #include "NIWGReWeight1piAngle.h" #include "NIWGReWeight2010a.h" #include "NIWGReWeight2012a.h" #include "NIWGReWeight2014a.h" #include "NIWGReWeightDeltaMass.h" #include "NIWGReWeightEffectiveRPA.h" #include "NIWGReWeightHadronMultSwitch.h" #include "NIWGReWeightMEC.h" #include "NIWGReWeightPiMult.h" #include "NIWGReWeightProtonFSIbug.h" #include "NIWGReWeightRPA.h" #include "NIWGReWeightSpectralFunc.h" #include "NIWGReWeightSplineEnu.h" #include "NIWGSyst.h" #include "NIWGSystUncertainty.h" #endif #ifdef __NEUT_ENABLED__ #include "NReWeight.h" #include "NReWeightCasc.h" #include "NReWeightNuXSecCCQE.h" #include "NReWeightNuXSecCCRES.h" #include "NReWeightNuXSecCOH.h" #include "NReWeightNuXSecDIS.h" #include "NReWeightNuXSecNC.h" #include "NReWeightNuXSecNCEL.h" #include "NReWeightNuXSecNCRES.h" #include "NReWeightNuXSecRES.h" #include "NReWeightNuclPiless.h" #include "NSyst.h" #include "NSystUncertainty.h" #include "neutpart.h" #include "neutvect.h" #endif #ifdef __NUWRO_ENABLED__ #include "event1.h" #endif #ifdef __NUWRO_REWEIGHT_ENABLED__ #include "NuwroReWeight.h" #include "NuwroReWeight_FlagNorm.h" #include "NuwroReWeight_QEL.h" #include "NuwroReWeight_SPP.h" #include "NuwroSyst.h" #include "NuwroSystUncertainty.h" #endif #ifdef __GENIE_ENABLED__ -#include "EVGCore/EventRecord.h" +#ifdef GENIE_PRE_R3 #include "EVGCore/EventRecord.h" #include "GHEP/GHepRecord.h" -#include "GSyst.h" -#include "GSystUncertainty.h" +#include "ReWeight/GSyst.h" +#include "ReWeight/GSystUncertainty.h" #include "Ntuple/NtpMCEventRecord.h" #include "ReWeight/GReWeight.h" #include "ReWeight/GReWeightAGKY.h" #include "ReWeight/GReWeightDISNuclMod.h" #include "ReWeight/GReWeightFGM.h" #include "ReWeight/GReWeightFZone.h" #include "ReWeight/GReWeightINuke.h" #include "ReWeight/GReWeightNonResonanceBkg.h" #include "ReWeight/GReWeightNuXSecCCQE.h" #include "ReWeight/GReWeightNuXSecCCQEvec.h" #include "ReWeight/GReWeightNuXSecCCRES.h" #include "ReWeight/GReWeightNuXSecCOH.h" #include "ReWeight/GReWeightNuXSecDIS.h" #include "ReWeight/GReWeightNuXSecNC.h" #include "ReWeight/GReWeightNuXSecNCEL.h" #include "ReWeight/GReWeightNuXSecNCRES.h" #include "ReWeight/GReWeightResonanceDecay.h" +#else +#include "Framework/EventGen/EventRecord.h" +#include "Framework/GHEP/GHepRecord.h" +#include "RwFramework/GSyst.h" +#include "RwFramework/GSystUncertainty.h" +#include "Framework/Ntuple/NtpMCEventRecord.h" +#include "RwFramework/GReWeight.h" +#include "RwCalculators/GReWeightAGKY.h" +#include "RwCalculators/GReWeightDISNuclMod.h" +#include "RwCalculators/GReWeightFGM.h" +#include "RwCalculators/GReWeightFZone.h" +#include "RwCalculators/GReWeightINuke.h" +#include "RwCalculators/GReWeightNonResonanceBkg.h" +#include "RwCalculators/GReWeightNuXSecCCQE.h" +#include "RwCalculators/GReWeightNuXSecCCQEvec.h" +#include "RwCalculators/GReWeightNuXSecCCRES.h" +#include "RwCalculators/GReWeightNuXSecCOH.h" +#include "RwCalculators/GReWeightNuXSecDIS.h" +#include "RwCalculators/GReWeightNuXSecNC.h" +#include "RwCalculators/GReWeightNuXSecNCEL.h" +#include "RwCalculators/GReWeightNuXSecNCRES.h" +#include "RwCalculators/GReWeightResonanceDecay.h" +#endif using namespace genie; using namespace genie::rew; #endif #include "GlobalDialList.h" #include "ModeNormEngine.h" #include "NUISANCESyst.h" #include "OscWeightEngine.h" //******************************************************************** TF1 FitBase::GetRWConvFunction(std::string const &type, std::string const &name) { //******************************************************************** std::string dialfunc = "x"; std::string parType = type; double low = -10000.0; double high = 10000.0; if (parType.find("parameter") == std::string::npos) parType += "_parameter"; std::string line; std::ifstream card( (GeneralUtils::GetTopLevelDir() + "/parameters/dial_conversion.card") .c_str(), std::ifstream::in); while (std::getline(card >> std::ws, line, '\n')) { std::vector inputlist = GeneralUtils::ParseToStr(line, " "); // Check the line length if (inputlist.size() < 4) continue; // Check whether this is a comment if (inputlist[0].c_str()[0] == '#') continue; // Check whether this is the correct parameter type if (inputlist[0].compare(parType) != 0) continue; // Check the parameter name if (inputlist[1].compare(name) != 0) continue; // inputlist[2] should be the units... ignore for now dialfunc = inputlist[3]; // High and low are optional, check whether they exist if (inputlist.size() > 4) low = GeneralUtils::StrToDbl(inputlist[4]); if (inputlist.size() > 5) high = GeneralUtils::StrToDbl(inputlist[5]); } TF1 convfunc = TF1((name + "_convfunc").c_str(), dialfunc.c_str(), low, high); return convfunc; } //******************************************************************** std::string FitBase::GetRWUnits(std::string const &type, std::string const &name) { //******************************************************************** std::string unit = "sig."; std::string parType = type; if (parType.find("parameter") == std::string::npos) { parType += "_parameter"; } std::string line; std::ifstream card( (GeneralUtils::GetTopLevelDir() + "/parameters/dial_conversion.card") .c_str(), std::ifstream::in); while (std::getline(card >> std::ws, line, '\n')) { std::vector inputlist = GeneralUtils::ParseToStr(line, " "); // Check the line length if (inputlist.size() < 3) continue; // Check whether this is a comment if (inputlist[0].c_str()[0] == '#') continue; // Check whether this is the correct parameter type if (inputlist[0].compare(parType) != 0) continue; // Check the parameter name if (inputlist[1].compare(name) != 0) continue; unit = inputlist[2]; break; } return unit; } //******************************************************************** double FitBase::RWAbsToSigma(std::string const &type, std::string const &name, double val) { //******************************************************************** TF1 f1 = GetRWConvFunction(type, name); double conv_val = f1.GetX(val); if (fabs(conv_val) < 1E-10) conv_val = 0.0; QLOG(FIT, "AbsToSigma(" << name << ") = " << val << " -> " << conv_val); return conv_val; } //******************************************************************** double FitBase::RWSigmaToAbs(std::string const &type, std::string const &name, double val) { //******************************************************************** TF1 f1 = GetRWConvFunction(type, name); double conv_val = f1.Eval(val); return conv_val; } //******************************************************************** double FitBase::RWFracToSigma(std::string const &type, std::string const &name, double val) { //******************************************************************** TF1 f1 = GetRWConvFunction(type, name); double conv_val = f1.GetX((val * f1.Eval(0.0))); if (fabs(conv_val) < 1E-10) conv_val = 0.0; return conv_val; } //******************************************************************** double FitBase::RWSigmaToFrac(std::string const &type, std::string const &name, double val) { //******************************************************************** TF1 f1 = GetRWConvFunction(type, name); double conv_val = f1.Eval(val) / f1.Eval(0.0); return conv_val; } int FitBase::ConvDialType(std::string const &type) { if (!type.compare("neut_parameter")) return kNEUT; else if (!type.compare("niwg_parameter")) return kNIWG; else if (!type.compare("nuwro_parameter")) return kNUWRO; else if (!type.compare("t2k_parameter")) return kT2K; else if (!type.compare("genie_parameter")) return kGENIE; else if (!type.compare("custom_parameter")) return kCUSTOM; else if (!type.compare("norm_parameter")) return kNORM; else if (!type.compare("likeweight_parameter")) return kLIKEWEIGHT; else if (!type.compare("spline_parameter")) return kSPLINEPARAMETER; else if (!type.compare("osc_parameter")) return kOSCILLATION; else if (!type.compare("modenorm_parameter")) return kMODENORM; else return kUNKNOWN; } std::string FitBase::ConvDialType(int type) { switch (type) { case kNEUT: { return "neut_parameter"; } case kNIWG: { return "niwg_parameter"; } case kNUWRO: { return "nuwro_parameter"; } case kT2K: { return "t2k_parameter"; } case kGENIE: { return "genie_parameter"; } case kNORM: { return "norm_parameter"; } case kCUSTOM: { return "custom_parameter"; } case kLIKEWEIGHT: { return "likeweight_parameter"; } case kSPLINEPARAMETER: { return "spline_parameter"; } case kOSCILLATION: { return "osc_parameter"; } case kMODENORM: { return "modenorm_parameter"; } default: return "unknown_parameter"; } } int FitBase::GetDialEnum(std::string const &type, std::string const &name) { return FitBase::GetDialEnum(FitBase::ConvDialType(type), name); } int FitBase::GetDialEnum(int type, std::string const &name) { int offset = type * 1000; int this_enum = Reweight::kNoDialFound; // Not Found QLOG(DEB, "Getting dial enum " << type << " " << name); // Select Types switch (type) { // NEUT DIAL TYPE case kNEUT: { #ifdef __NEUT_ENABLED__ int neut_enum = (int)neut::rew::NSyst::FromString(name); if (neut_enum != 0) { this_enum = neut_enum + offset; } #else this_enum = Reweight::kNoTypeFound; // Not enabled #endif break; } // NIWG DIAL TYPE case kNIWG: { #ifdef __NIWG_ENABLED__ int niwg_enum = (int)niwg::rew::NIWGSyst::FromString(name); if (niwg_enum != 0) { this_enum = niwg_enum + offset; } #else this_enum = Reweight::kNoTypeFound; #endif break; } // NUWRO DIAL TYPE case kNUWRO: { #ifdef __NUWRO_REWEIGHT_ENABLED__ int nuwro_enum = (int)nuwro::rew::NuwroSyst::FromString(name); if (nuwro_enum > 0) { this_enum = nuwro_enum + offset; } #else this_enum = Reweight::kNoTypeFound; #endif } // GENIE DIAL TYPE case kGENIE: { #ifdef __GENIE_ENABLED__ int genie_enum = (int)genie::rew::GSyst::FromString(name); if (genie_enum > 0) { this_enum = genie_enum + offset; } #else this_enum = Reweight::kNoTypeFound; #endif break; } case kCUSTOM: { int custom_enum = 0; // PLACEHOLDER this_enum = custom_enum + offset; break; } // T2K DIAL TYPE case kT2K: { #ifdef __T2KREW_ENABLED__ int t2k_enum = (int)t2krew::T2KSyst::FromString(name); if (t2k_enum > 0) { this_enum = t2k_enum + offset; } #else this_enum = Reweight::kNoTypeFound; #endif break; } case kNORM: { if (gNormEnums.find(name) == gNormEnums.end()) { gNormEnums[name] = gNormEnums.size() + 1 + offset; } this_enum = gNormEnums[name]; break; } case kLIKEWEIGHT: { if (gLikeWeightEnums.find(name) == gLikeWeightEnums.end()) { gLikeWeightEnums[name] = gLikeWeightEnums.size() + 1 + offset; } this_enum = gLikeWeightEnums[name]; break; } case kSPLINEPARAMETER: { if (gSplineParameterEnums.find(name) == gSplineParameterEnums.end()) { gSplineParameterEnums[name] = gSplineParameterEnums.size() + 1 + offset; } this_enum = gSplineParameterEnums[name]; break; } case kOSCILLATION: { #ifdef __PROB3PP_ENABLED__ int oscEnum = OscWeightEngine::SystEnumFromString(name); if (oscEnum != 0) { this_enum = oscEnum + offset; } #else this_enum = Reweight::kNoTypeFound; // Not enabled #endif } case kMODENORM: { size_t us_pos = name.find_first_of('_'); std::string numstr = name.substr(us_pos + 1); int mode_num = std::atoi(numstr.c_str()); LOG(FTL) << "Getting mode num " << mode_num << std::endl; if (!mode_num) { ERR(FTL) << "Attempting to parse dial name: \"" << name << "\" as a mode norm dial but failed." << std::endl; throw; } this_enum = 60 + mode_num + offset; break; } } // If Not Enabled if (this_enum == Reweight::kNoTypeFound) { ERR(FTL) << "RW Engine not supported for " << FitBase::ConvDialType(type) << std::endl; ERR(FTL) << "Check dial " << name << std::endl; } // If Not Found if (this_enum == Reweight::kNoDialFound) { ERR(FTL) << "Dial " << name << " not found." << std::endl; } return this_enum; } int Reweight::ConvDialType(std::string const &type) { return FitBase::ConvDialType(type); } std::string Reweight::ConvDialType(int type) { return FitBase::ConvDialType(type); } int Reweight::GetDialType(int type) { int t = (type / 1000); return t > kMODENORM ? Reweight::kNoDialFound : t; } int Reweight::RemoveDialType(int type) { return (type % 1000); } int Reweight::NEUTEnumFromName(std::string const &name) { #ifdef __NEUT_ENABLED__ int neutenum = (int)neut::rew::NSyst::FromString(name); return (neutenum > 0) ? neutenum : Reweight::kNoDialFound; #else return Reweight::kGeneratorNotBuilt; #endif } int Reweight::NIWGEnumFromName(std::string const &name) { #ifdef __NIWG_ENABLED__ int niwgenum = (int)niwg::rew::NIWGSyst::FromString(name); return (niwgenum != 0) ? niwgenum : Reweight::kNoDialFound; #else return Reweight::kGeneratorNotBuilt; #endif } int Reweight::NUWROEnumFromName(std::string const &name) { #ifdef __NUWRO_REWEIGHT_ENABLED__ int nuwroenum = (int)nuwro::rew::NuwroSyst::FromString(name); return (nuwroenum > 0) ? nuwroenum : Reweight::kNoDialFound; #else return Reweight::kGeneratorNotBuilt; #endif } int Reweight::GENIEEnumFromName(std::string const &name) { #ifdef __GENIE_ENABLED__ int genieenum = (int)genie::rew::GSyst::FromString(name); return (genieenum > 0) ? genieenum : Reweight::kNoDialFound; #else return Reweight::kGeneratorNotBuilt; #endif } int Reweight::T2KEnumFromName(std::string const &name) { #ifdef __T2KREW_ENABLED__ int t2kenum = (int)t2krew::T2KSyst::FromString(name); return (t2kenum > 0) ? t2kenum : Reweight::kNoDialFound; #else return Reweight::kGeneratorNotBuilt; #endif } int Reweight::OscillationEnumFromName(std::string const &name) { #ifdef __PROB3PP_ENABLED__ int oscEnum = OscWeightEngine::SystEnumFromString(name); return (oscEnum > 0) ? oscEnum : Reweight::kNoDialFound; #else return Reweight::kGeneratorNotBuilt; #endif } int Reweight::NUISANCEEnumFromName(std::string const &name, int type) { int nuisenum = Reweight::DialList().EnumFromNameAndType(name, type); return nuisenum; } int Reweight::CustomEnumFromName(std::string const &name) { int custenum = Reweight::ConvertNUISANCEDial(name); return custenum; } int Reweight::ConvDial(std::string const &name, std::string const &type, bool exceptions) { return Reweight::ConvDial(name, Reweight::ConvDialType(type), exceptions); } int Reweight::ConvDial(std::string const &fullname, int type, bool exceptions) { std::string name = GeneralUtils::ParseToStr(fullname, ",")[0]; // Only use first dial given // Produce offset seperating each type. int offset = type * 1000; int genenum = Reweight::kNoDialFound; switch (type) { case kNEUT: genenum = NEUTEnumFromName(name); break; case kNIWG: genenum = NIWGEnumFromName(name); break; case kNUWRO: genenum = NUWROEnumFromName(name); break; case kGENIE: genenum = GENIEEnumFromName(name); break; case kT2K: genenum = T2KEnumFromName(name); break; case kCUSTOM: genenum = CustomEnumFromName(name); break; case kNORM: case kLIKEWEIGHT: case kSPLINEPARAMETER: case kNEWSPLINE: genenum = NUISANCEEnumFromName(name, type); break; case kOSCILLATION: genenum = OscillationEnumFromName(name); break; case kMODENORM: genenum = ModeNormEngine::SystEnumFromString(name); break; default: genenum = Reweight::kNoTypeFound; break; } // Throw if required. if (exceptions) { // If Not Enabled if (genenum == Reweight::kGeneratorNotBuilt) { ERR(FTL) << "RW Engine not supported for " << FitBase::ConvDialType(type) << std::endl; ERR(FTL) << "Check dial " << name << std::endl; throw; } // If no type enabled if (genenum == Reweight::kNoTypeFound) { ERR(FTL) << "Type mismatch inside ConvDialEnum" << std::endl; throw; } // If Not Found if (genenum == Reweight::kNoDialFound) { ERR(FTL) << "Dial " << name << " not found." << std::endl; throw; } } // Add offset if no issue int nuisenum = genenum; if ((genenum != Reweight::kGeneratorNotBuilt) && (genenum != Reweight::kNoTypeFound) && (genenum != Reweight::kNoDialFound)) { nuisenum += offset; } // Now register dial Reweight::DialList().RegisterDialEnum(name, type, nuisenum); return nuisenum; } std::string Reweight::ConvDial(int nuisenum) { for (size_t i = 0; i < Reweight::DialList().fAllDialEnums.size(); i++) { if (Reweight::DialList().fAllDialEnums[i] == nuisenum) { return Reweight::DialList().fAllDialNames[i]; } } LOG(FIT) << "Cannot find dial with enum = " << nuisenum << std::endl; return ""; }