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 <stdio.h>
#include <stdlib.h>
#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<std::string, TH1D*> modexsec;
std::map<std::string, TH1D*> modecount;
std::vector<std::string> genieids;
std::vector<std::string> targetids;
std::vector<std::string> 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<GHepRecord>(event);
double xsec = (genie_record.XSec() / (1E-38 * genie::units::cm2));
// Parse Interaction String
std::string mode = genie_record.Summary()->AsString();
std::vector<std::string> 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<std::string, TH1D*> 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<std::string, TH1D*> 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<std::string> targprs = GeneralUtils::ParseToStr(target, ",");
std::vector<std::string> targ_list;
std::vector<std::string> 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<std::string>::iterator it = targprs.begin(); it != targprs.end(); ++it) {
// Cut into "TARGET1" and "fraction1]"
std::vector<std::string> targind = GeneralUtils::ParseToStr(*it, "[");
//std::cout << " " << *it << std::endl;
// Cut into "TARGET1" and "fraction1"
for (std::vector<std::string>::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<double> targ_fractions;
double minimum = 1.0;
for (std::vector<std::string>::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<double>::iterator it = targ_fractions.begin();
std::vector<std::string>::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<std::string, TH1D*>::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<std::string, TH1D*>::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<std::string> 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<TH1*>(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<std::string> 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<std::string, TH1D*> modexsec;
std::map<std::string, TH1D*> modecount;
std::vector<std::string> genieids;
std::vector<std::string> targetids;
std::vector<std::string> 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<GHepRecord>(event);
double xsec = (genie_record.XSec() / (1E-38 * genie::units::cm2));
// Parse Interaction String
std::string mode = genie_record.Summary()->AsString();
std::vector<std::string> 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<std::string, TH1D*> 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<std::string, TH1D*> 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<std::string> targprs = GeneralUtils::ParseToStr(target, ",");
std::vector<std::string> targ_list;
std::vector<std::string> 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<std::string>::iterator it = targprs.begin(); it != targprs.end(); ++it) {
// Cut into "TARGET1" and "fraction1]"
std::vector<std::string> targind = GeneralUtils::ParseToStr(*it, "[");
//std::cout << " " << *it << std::endl;
// Cut into "TARGET1" and "fraction1"
for (std::vector<std::string>::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<double> targ_fractions;
double minimum = 1.0;
for (std::vector<std::string>::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<double>::iterator it = targ_fractions.begin();
std::vector<std::string>::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<std::string, TH1D*>::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<std::string, TH1D*>::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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
// 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<<std::endl;
std::cout<<"Arguments:"<<std::endl;
std::cout<<" -c cardFile: Path to card file that defines fit samples, free parameters, and config overrides \n";
std::cout<<" -o outFile: Path to root file that will be created to save output file.\n";
std::cout<<" To turn automatically overwrite outFile if one exists turn off use 'config overwrite_output 1'\n";
std::cout<<" To automatically use previous outFile as an inputFile if it exists so that the fit can be continued\n";
std::cout<<" use the flag 'config use_previous_output 1'. (A warning will be printed when doing this). \n";
std::cout<<" -f fitStategy: Pass a comma seperated list of fit routines to run in order. Default is Migrad,FixAtLim \n";
std::cout<<" Possible Options: \n";
std::cout<<" 1. Migrad - Minuit2 Migrad Minimizer \n";
std::cout<<" 2. Simplex - Simplex Minimizer \n";
std::cout<<" 3. Scan - Brute force scan of parameter space \n";
std::cout<<" 4. FixAtLim - Takes any free parameters close to a limit and fixes them \n";
std::cout<<" 5. Scan1D - Make 1D Scans and save them in a folder \n";
std::cout<<" 6. Contours - Make Contour Scans \n";
std::cout<<" 7. Save - Will save the state of the fitter (Always done by default at the end) \n";
std::cout<<" Extra option LowStatFit will perform each of these options with a lower number \n";
std::cout<<" of fit events (config lowstat). Example: LowStatMigrad, LowStatScan \n";
std::cout<<" -d fakeDataFile: Uses the MC generated from a previous fit as a fake data set for these fits \n";
std::cout<<" Can also specify MC to set the fake-data to the Monte-Carlo prediction\n";
std::cout<<" In this case, you can specify fake_parameter PARAM_NAME PARAM_VALUE in the card\n";
std::cout<<" to reweight the MC parameter PARAM_NAME to some PARAM_VALUE. The minimiser will start\n";
std::cout<<" at whatever neut_parameter or genie_parameter is set to in the cardfile.\n";
std::cout<<" -i inputFile: Uses the results from a previous fit file as starting input for these fits \n";
std::cout<<" -q config_name=config_val : Allows any config parameter to be overridden from the command line.\n";
std::cout<<" This will take priority over those given in the default, or cardFile. \n";
std::cout<<" example: -q verbosity=6 -q maxevents=10000 \n";
exit(-1);
};
//*******************************
int main(int argc, char* argv[]){
//*******************************
// Program status;
int status = 0;
// If No Arguments print commands
if (argc == 1) printInputCommands();
for (int i = 1; i< argc; ++i){
// Cardfile
if (!std::strcmp(argv[i], "-h")) printInputCommands();
else break;
}
// Read input arguments such as card file, parameter arguments, and fit routines
LOG(FIT)<<"Starting ExtFit_minimizer.exe"<<std::endl;
// Make minimizer class and run fit
MinimizerRoutines* min = new MinimizerRoutines(argc, argv);
// Save Starting States
if (FitPar::Config().GetParB("savenominal")) min->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)<<"-------------------------------------"<<std::endl;
if (status == 0) LOG(FIT)<<"Minimizer Complete."<<std::endl;
else ERR(WRN)<<"Minimizer Failed (error state = "<<status<<")"<<std::endl;
LOG(FIT)<<"-------------------------------------"<<std::endl;
return status;
}
diff --git a/cmake/GENIESetup.cmake b/cmake/GENIESetup.cmake
index fe49732..09d6572 100644
--- a/cmake/GENIESetup.cmake
+++ b/cmake/GENIESetup.cmake
@@ -1,186 +1,227 @@
# 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 <http://www.gnu.org/licenses/>.
################################################################################
# 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. <FALSE>" 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 @@
<nuisance>
<!-- # ###################################################### -->
<!-- # NUISANCE CONFIGURATION OPTIONS -->
<!-- # This file is read in by default at runtime -->
<!-- # If you want to override on a case by case bases use -q at runtime -->
<!-- # ###################################################### -->
<!-- # MAIN Configs -->
<!-- # ###################################################### -->
<!-- # Logger goes from -->
<!-- # 1 Quiet -->
<!-- # 2 Fitter -->
<!-- # 3 Samples -->
<!-- # 4 Reconfigure Loops -->
<!-- # 5 Every Event print out (SHOUT) -->
<!-- # -1 DEBUGGING -->
<config VERBOSITY='4'/>
<!-- # ERROR goes from -->
<!-- # 0 NONE -->
<!-- # 1 FATAL -->
<!-- # 2 WARN -->
<config ERROR='2'/>
<config TRACE='0'/>
<config cores='1' />
<config spline_test_throws='50' />
<config spline_cores='1' />
<config spline_chunks='20' />
<config spline_procchunk='-1' />
<config Electron_NThetaBins='4' />
<config Electron_NEnergyBins='4' />
<config Electron_ThetaWidth='1.0' />
<config Electron_EnergyWidth='0.10' />
<config RemoveFSIParticles='0' />
<config RemoveUndefParticles='0' />
<config RemoveNuclearParticles='0'/>
<config MINERvASaveExtraCCQE='0' />
<!-- # Input Configs -->
<!-- # ###################################################### -->
<!-- # Default Requirements file for the externalDataFitter Package -->
<!-- # MAX Events : -1 is no cut. Events will be scaled automatically to give good xsec predictions. -->
<config MAXEVENTS='-1'/>
<!-- Include empty stacks in the THStack -->
<config includeemptystackhists='0'/>
<!-- # Turn on/off event manager -->
<!-- # EventManager enables us to only loop number of events once for multiple projections of the same measurements -->
<!-- # e.g. MiniBooNE CC1pi+ Q2 and MiniBooNE CC1pi+ Tmu would ordinarily require 2 reconfigures, but with this enabled it requires only one -->
<config EventManager='1'/>
<!-- # Event Directories -->
<!-- # Can setup default directories and use @EVENT_DIR/path to link to it -->
<config EVENT_DIR='/data2/stowell/NIWG/'/>
<config NEUT_EVENT_DIR='/data2/stowell/NIWG/neut/fit_samples_neut5.3.3/'/>
<config GENIE_EVENT_DIR='/data2/stowell/NIWG/genie/fit_samples_R.2.10.0/'/>
<config NUWRO_EVENT_DIR='/data2/stowell/NIWG/nuwro/fit_samples/'/>
<config GIBUU_EVENT_DIR='/data/GIBUU/DIR/'/>
<config SaveNuWroExtra='0' />
<!-- # In PrepareGENIE the reconstructed splines can be saved into the file -->
<config save_genie_splines='1'/>
<!-- # In InputHandler the option to regenerate NuWro flux/xsec plots is available -->
<!-- # Going to move this to its own app soon -->
<!-- # DEVEL CONFIG OPTION, don't touch! -->
<config CacheSize='0'/>
<!-- # ReWeighting Configuration Options -->
<!-- # ###################################################### -->
<!-- # Convert Dials in output statements using dial conversion card -->
<config convert_dials='0'/>
<!-- # Vetos can be used to specify RW dials NOT to be loaded in -->
<!-- # Useful if one specific one has an issue -->
<config FitWeight_fNIWGRW_veto=''/>
<config FitWeight_fNuwroRW_veto=''/>
<config FitWeight_fNeutRW_veto=''/>
<config FitWeight_fGenieRW_veto=''/>
<!-- # Output Options -->
<!-- # ###################################################### -->
<!-- # Save Nominal prediction with all rw engines at default -->
<config savenominal='0'/>
<!-- # Save prefit with values at starting values -->
<config saveprefit='0'/>
<!-- # Here's the full list of drawing options -->
<!-- DATA/MC/EVT/FINE/RATIO/MODES/SHAPE/WGHT/WEIGHTS/FLUX/XSEC/MASK/COV/INCOV/DECOMP/CANVPDG/CANVMC/SETTINGS'/>
<!-- #config drawopts DATA/MC/EVT/FINE/RATIO/MODES/SHAPE/RESIDUAL/MATRIX/FLUX/MASK/MAP -->
<!-- #config drawopts DATA/MC -->
<config drawopts='DATA/MC/EVT/FINE/RATIO/MODES/SHAPE/WEIGHTS/FLUX/XSEC/MASK/COV/INCOV/DECOMP/CANVPDG/CANVMC/SETTINGS'/>
<config InterpolateSigmaQ0Histogram='1' />
<config InterpolateSigmaQ0HistogramRes='100' />
<config InterpolateSigmaQ0HistogramThrow='1' />
<config InterpolateSigmaQ0HistogramNTHROWS='100000' />
<!-- # Save the shape scaling applied with option SHAPE into the main MC hist -->
<config saveshapescaling='0'/>
<config CorrectGENIEMECNorm='1'/>
<!-- # Set style of 1D output histograms -->
<config linecolour='1'/>
<config linestyle='1'/>
<config linewidth='1'/>
<!-- # For GenericFlux -->
<config isLiteMode='0'/>
<!-- # Statistical Options -->
<!-- # ###################################################### -->
<!-- # Add MC Statistical error to likelihoods -->
<config addmcerror='0'/>
<!-- # NUISMIN Configurations -->
<!-- # ###################################################### -->
<config MAXCALLS='1000000'/>
<config MAXITERATIONS='1000000'/>
<config TOLERANCE='0.001'/>
<!-- # Number of events required in low stats routines -->
<config LOWSTATEVENTS='25000'/>
<!-- # Error band generator configs -->
<!-- # ###################################################### -->
<!-- # For -f ErrorBands creates error bands for given measurements -->
<!-- # How many throws do we want (The higher the better precision) -->
<config error_throws='250'/>
<!-- # Are we throwing uniform or according to Gaussian? -->
<!-- # Only use uniform if wanting to study the limits of a dial. -->
<config error_uniform='0'/>
<config WriteSeperateStacks='1'/>
<!-- # Other Individual Case Configs -->
<!-- # ###################################################### -->
<!-- # Covariance throw options for fake data studies with MiniBooNE data. -->
<config thrown_covariance='FULL'/>
<config throw_mc_stat='0.0'/>
<config throw_diag_syst='0'/>
<config throw_corr_syst='0'/>
<config throw_mc_stat='0'/>
<!-- # Apply a shift to the muon momentum before calculation of Q2 -->
<config muon_momentum_shift='0.0'/>
<config muon_momentum_throw='0'/>
<!-- # MINERvA Specific Configs -->
<config MINERvA_CCinc_XSec_2DEavq3_nu_hadron_cut='0'/>
<config MINERvA_CCinc_XSec_2DEavq3_nu_useq3true='0'/>
<config Modes_split_PN_NN='0'/>
<config SignalReconfigures='0'/>
<!-- # SciBooNE specific -->
<config SciBarDensity='1.04'/>
<config SciBarRecoDist='12.0'/>
<config PenetratingMuonEnergy='1.4'/>
<config FlatEfficiency='1.0'/>
<config NumRangeSteps='50'/>
<config UseProton='true'/>
<config UseZackEff='false'/>
<config MINERvADensity='1.04'/>
<config MINERvARecoDist='10.0'/>
<!-- Different way of reweighting in GENIE -->
<config GENIEWeightEngine_CCRESMode="kModeMaMv"/>
<!--
<config GENIEWeightEngine_CCQEMode="kModeMa"/>
-->
<!-- CCQE/2p2h/1pi Gaussian enhancement method -->
<!-- Apply tilt-shift weights or normal Gaussian parameters-->
<config Gaussian_Enhancement="Normal" />
<!--
<config Gaussian_Enhancement="Tilt-Shift" />
-->
<config NToyThrows='100000' />
<!-- Use NOvA Weights or not -->
<config NOvA_Weights="false" />
+<!-- Tune name, for GENIE v3+ -->
+<config GENIETune="G18_02a_00_000" />
+<config GENIEEventGeneratorList="Default" />
+
<!--
<config GENIEXSecModelCCRES="genie::ReinSehgalRESPXSec" />
<config GENIEXSecModelCOH="genie::ReinSehgalCOHPiPXSec" />
<config GENIEXSecModelCCQE="genie::LwlynSmithQELCCPXSec" />
-->
<!--
<config GENIEXSecModelCCQE="genie::NievesQELCCPXSec" />
<config GENIEXSecModelCCRES="genie::BergerSehgalRESPXSec2014" />
<config GENIEXSecModelCOH="genie::BergerSehgalCOHPiPXSec2015" />
-->
<config UseShapeCovar="0" />
<config CC0piNBINS="156" />
</nuisance>
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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
#ifndef INPUTHANDLER_H_SEEN
#define INPUTHANDLER_H_SEEN
/*!
* \addtogroup FitBase
* @{
*/
// C Includes
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <deque>
#include <iomanip>
#include <iostream>
#include <list>
#include <numeric>
#include <sstream>
#include <string>
// ROOT includes
#include <TArrayF.h>
#include <TDecompChol.h>
#include <TDecompSVD.h>
#include <TGraph.h>
#include <TGraphErrors.h>
#include <TH1D.h>
#include <TH2D.h>
#include <TMatrixDSym.h>
#include <TROOT.h>
#include <TSystem.h>
// 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<TH1*> GetDataList(void) = 0;
///! Virtual function to get MC histogram
virtual std::vector<TH1*> GetMCList(void) = 0;
virtual std::vector<TH1*> GetFineList(void) = 0;
virtual std::vector<TH1*> GetMaskList(void) = 0;
///! Return flux histograms in a vector
virtual std::vector<TH1*> GetFluxList(void);
virtual std::vector<TH1*> GetEventRateList(void);
virtual std::vector<TH1*> 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<MeasurementBase*> GetSubSamples() {
return std::vector<MeasurementBase*>(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<StackBase*, std::vector<int> > fExtraTH1s;
int NSignal;
// std::map<TH1*, bool[6] > fExtaStacks;
bool fIsJoint;
double fNPOT, fFluxIntegralOverride, fTargetVolume, fTargetMaterialDensity;
double fEvtRateScaleFactor;
};
// Class TypeDefs
typedef std::list<MeasurementBase*>::const_iterator MeasListConstIter;
typedef std::list<MeasurementBase*>::iterator MeasListIter;
typedef std::vector<MeasurementBase*>::const_iterator MeasVectConstIter;
typedef std::vector<MeasurementBase*>::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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
#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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
#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<GReWeightNuXSecCCQE*> (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<GReWeightNuXSecCOH*> (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<GReWeightNuXSecCCQE *>(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<GReWeightNuXSecCCRES *>(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<GReWeightNuXSecNCRES *>(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<GReWeightNuXSecDIS *>(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<size_t>(0);
fNameIndex[name]; // = std::vector<size_t>(0);
// Split by commas
std::vector<std::string> 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<size_t> 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<size_t> 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<genie::rew::GSyst_t> 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 <string>
#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<std::string> 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<std::string> 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 "";
}