diff --git a/app/PrepareGENIE.cxx b/app/PrepareGENIE.cxx index ca21ce0..b67d856 100644 --- a/app/PrepareGENIE.cxx +++ b/app/PrepareGENIE.cxx @@ -1,866 +1,942 @@ -#include -#include #include "FitLogger.h" #include "PlotUtils.h" #include "TFile.h" #include "TH1D.h" #include "TTree.h" +#include +#include #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 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[]) { +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; + LOG(FIT) << "Running GENIE Prepare in mono energetic with E = " << MonoEnergy + << " GeV" << std::endl; // Setup TTree - TChain* tn = new TChain("gtree"); + 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; + LOG(FIT) << "Overriding number of events by user from " << nevt << " to " + << gNEvents << std::endl; nevt = gNEvents; } - NtpMCEventRecord* genientpl = NULL; + 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.); + 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(); + TH1D *eventhist = (TH1D *)fluxhist->Clone(); eventhist->Reset(); - TH1D* xsechist = (TH1D*)eventhist->Clone(); + TH1D *xsechist = (TH1D *)eventhist->Clone(); // Create maps - std::map modexsec; - std::map modecount; + std::map modexsec; + std::map modecount; std::vector genieids; std::vector targetids; std::vector interids; // Loop over all events for (int i = 0; i < nevt; i++) { tn->GetEntry(i); StopTalking(); - EventRecord& event = *(genientpl->event); - GHepParticle* neu = event.Probe(); + EventRecord &event = *(genientpl->event); + GHepParticle *neu = event.Probe(); StartTalking(); // Get XSec From Spline GHepRecord genie_record = static_cast(event); double xsec = (genie_record.XSec() / (1E-38 * genie::units::cm2)); // Parse Interaction String std::string mode = genie_record.Summary()->AsString(); std::vector modevec = GeneralUtils::ParseToStr(mode, ";"); std::string targ = (modevec[0] + ";" + modevec[1]); std::string inter = mode; // Fill lists of Unique IDS if (std::find(targetids.begin(), targetids.end(), targ) == targetids.end()) { targetids.push_back(targ); } if (std::find(interids.begin(), interids.end(), inter) == interids.end()) { interids.push_back(inter); } // Create entries Mode Maps if (modexsec.find(mode) == modexsec.end()) { genieids.push_back(mode); - modexsec[mode] = (TH1D*)xsechist->Clone(); - modecount[mode] = (TH1D*)xsechist->Clone(); + modexsec[mode] = (TH1D *)xsechist->Clone(); + modecount[mode] = (TH1D *)xsechist->Clone(); - modexsec[mode]->GetYaxis()->SetTitle("d#sigma/dE_{#nu} #times 10^{-38} (events weighted by #sigma)"); + 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; + << " GENIE events (E: " << neu->E() << " GeV, xsec: " << xsec + << " E-38 cm^2/nucleon)" << std::endl; } } LOG(FIT) << "Processed all events" << std::endl; - TFile* outputfile; + 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(); + 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(); + TTree *nova_tree = nova_chain->GetTree(); if (!nova_tree) { QLOG(FIT, "Could not find nova_wgts tree in " << gOutputFile); } else { QLOG(FIT, "Found nova_wgts tree in " << gOutputFile); } if (nova_tree) { nova_tree->SetDirectory(outputfile); nova_tree->Write(); } QLOG(FIT, "Done cloning tree."); } LOG(FIT) << "Getting splines in mono-energetic..." << std::endl; // Save each of the reconstructed splines to file - std::map modeavg; + std::map modeavg; - TDirectory* inddir = (TDirectory*)outputfile->Get("IndividualGENIESplines"); - if (!inddir) inddir = (TDirectory*)outputfile->mkdir("IndividualGENIESplines"); + 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] = (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"); + TDirectory *targdir = (TDirectory *)outputfile->Get("TargetGENIESplines"); + if (!targdir) + targdir = (TDirectory *)outputfile->mkdir("TargetGENIESplines"); targdir->cd(); LOG(FIT) << "Getting Target Splines" << std::endl; // For each target save a total spline - std::map targetsplines; + std::map targetsplines; for (uint i = 0; i < targetids.size(); i++) { std::string targ = targetids[i]; LOG(FIT) << "Getting target " << i << ": " << targ << std::endl; - targetsplines[targ] = (TH1D*)xsechist->Clone(); - targetsplines[targ]->GetYaxis()->SetTitle("#sigma (E_{#nu}) #times 10^{-38} (cm^{2}/target)"); + 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; + 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) << "Finished with Mode " << mode << " " + << modeavg[mode]->Integral() << std::endl; } } LOG(FIT) << "Saving target spline:" << targ << std::endl; targetsplines[targ]->Write(("Total_" + targ).c_str(), TObject::kOverwrite); } LOG(FIT) << "Getting total splines" << std::endl; // Now we have each of the targets we need to create a total cross-section. int totalnucl = 0; // Get the targets specified by the user, separated by commas // This has structure target1[fraction1], target2[fraction2] std::vector targprs = GeneralUtils::ParseToStr(target, ","); std::vector targ_list; std::vector frac_list; - // Chop up the target string which has format TARGET1[fraction1],TARGET2[fraction2] + // Chop up the target string which has format + // TARGET1[fraction1],TARGET2[fraction2] - //std::cout << "Targets: " << std::endl; + // std::cout << "Targets: " << std::endl; // Loop over the vector of strings "TARGET1[fraction1]" "TARGET2[fraction2]" - for (std::vector::iterator it = targprs.begin(); it != targprs.end(); ++it) { + for (std::vector::iterator it = targprs.begin(); + it != targprs.end(); ++it) { // Cut into "TARGET1" and "fraction1]" std::vector targind = GeneralUtils::ParseToStr(*it, "["); - //std::cout << " " << *it << std::endl; + // std::cout << " " << *it << std::endl; // Cut into "TARGET1" and "fraction1" - for (std::vector::iterator jt = targind.begin(); jt != targind.end(); ++jt) { + for (std::vector::iterator jt = targind.begin(); + jt != targind.end(); ++jt) { if ((*jt).find("]") != std::string::npos) { (*jt) = (*jt).substr(0, (*jt).find("]")); //*jt = "hello"; frac_list.push_back(*jt); // Won't find bracket for target } else { targ_list.push_back(*jt); } } } targprs = targ_list; std::vector targ_fractions; double minimum = 1.0; - for (std::vector::iterator it = frac_list.begin(); it != frac_list.end(); it++) { - //std::cout << " " << *it << std::endl; + for (std::vector::iterator it = frac_list.begin(); + it != frac_list.end(); it++) { + // std::cout << " " << *it << std::endl; double frac = std::atof((*it).c_str()); targ_fractions.push_back(frac); - if (frac < minimum) minimum = frac; + if (frac < minimum) + minimum = frac; } std::vector::iterator it = targ_fractions.begin(); std::vector::iterator jt = targ_list.begin(); double scaling = 0; for (; it != targ_fractions.end(); it++, jt++) { // First get the mass number from the targ_list int nucl = atoi((*jt).c_str()); - nucl = (nucl%10000)/10; + nucl = (nucl % 10000) / 10; // Gets the relative portions right - *it = (*it)/minimum; + *it = (*it) / minimum; // Scale relative the atomic mass //(*it) *= (double(nucl)/(*it)); - double tempscaling = double(nucl)/(*it); - if (tempscaling > scaling) scaling=tempscaling; + 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); + *it *= int(scaling + 0.5); // Round to nearest integer - *it = int(*it+0.5); + *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); + 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(); + TH1D *totalxsec = (TH1D *)xsechist->Clone(); for (uint i = 0; i < targprs.size(); i++) { std::string targpdg = targprs[i]; // Check that we found the user requested target in GENIE bool FoundTarget = false; - for (std::map::iterator iter = targetsplines.begin(); - iter != targetsplines.end(); iter++) { + for (std::map::iterator iter = targetsplines.begin(); + iter != targetsplines.end(); iter++) { std::string targstr = iter->first; - TH1D* xsec = iter->second; + 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; + 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); - } + // 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) << "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; + for (uint i = 0; i < targprs.size(); ++i) + ERR(WRN) << " " << targprs[i] << std::endl; ERR(WRN) << " The list of targets found in GENIE is: " << std::endl; - for (std::map::iterator iter = targetsplines.begin(); iter != targetsplines.end(); iter++) ERR(WRN) << " " << iter->first<< std::endl; + for (std::map::iterator iter = targetsplines.begin(); + iter != targetsplines.end(); iter++) + ERR(WRN) << " " << iter->first << std::endl; } } outputfile->cd(); - totalxsec->GetYaxis()->SetTitle("#sigma (E_{#nu}) #times 10^{-38} (cm^{2}/nucleon)"); + totalxsec->GetYaxis()->SetTitle( + "#sigma (E_{#nu}) #times 10^{-38} (cm^{2}/nucleon)"); totalxsec->Write("nuisance_xsec", TObject::kOverwrite); - eventhist = (TH1D*)fluxhist->Clone(); + 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()); + 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; + 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) { + std::string output) { LOG(FIT) << "Running GENIE Prepare with flux..." << std::endl; // Get Flux Hist std::vector fluxvect = GeneralUtils::ParseToStr(flux, ","); - TH1* fluxhist = NULL; + 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 << ")."); + << from << " to " << to << " with bins " << step + << " wide (NBins = " << nstep << ")."); - fluxhist = new TH1D("spectrum", ";E_{#nu} (GeV);Count (A.U.)", nstep, from, to); + 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"); + TFile *fluxfile = new TFile(fluxvect[0].c_str(), "READ"); if (!fluxfile->IsZombie()) { - fluxhist = dynamic_cast(fluxfile->Get(fluxvect[1].c_str())); + fluxhist = dynamic_cast(fluxfile->Get(fluxvect[1].c_str())); if (!fluxhist) { ERR(FTL) << "Couldn't find histogram named: \"" << fluxvect[1] - << "\" in file: \"" << fluxvect[0] << std::endl; + << "\" 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"); + TChain *tn = new TChain("gtree"); + std::string first_file = ""; if (input.find_first_of(',') != std::string::npos) { std::vector inputvect = GeneralUtils::ParseToStr(input, ","); for (size_t iv_it = 0; iv_it < inputvect.size(); ++iv_it) { tn->AddFile(inputvect[iv_it].c_str()); QLOG(FIT, "Added input file: " << inputvect[iv_it]); + if (!first_file.length()) { + first_file = inputvect[iv_it]; + } } - } else { // The Add form can accept wildcards. + } else { // The Add form can accept wildcards. tn->Add(input.c_str()); + first_file = input; } 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; + 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() << "\""); + << input.c_str() << "\""); } else { QLOG(FIT, "Found " << nevt << " input entries in " << input); } - NtpMCEventRecord* genientpl = NULL; + NtpMCEventRecord *genientpl = NULL; tn->SetBranchAddress("gmcrec", &genientpl); // Make Event and xsec Hist - TH1D* eventhist = (TH1D*)fluxhist->Clone(); + TH1D *eventhist = (TH1D *)fluxhist->Clone(); + eventhist->SetDirectory(NULL); eventhist->Reset(); - TH1D* xsechist = (TH1D*)eventhist->Clone(); + TH1D *xsechist = (TH1D *)eventhist->Clone(); + xsechist->SetDirectory(NULL); // Create maps - std::map modexsec; - std::map modecount; + std::map modexsec; + std::map modecount; std::vector genieids; std::vector targetids; std::vector interids; // Loop over all events for (int i = 0; i < nevt; i++) { tn->GetEntry(i); // Hussssch GENIE StopTalking(); // Get the event - EventRecord& event = *(genientpl->event); + EventRecord &event = *(genientpl->event); // Get the neutrino - GHepParticle* neu = event.Probe(); + GHepParticle *neu = event.Probe(); StartTalking(); // Get XSec From Spline // Get the GHepRecord GHepRecord genie_record = static_cast(event); double xsec = (genie_record.XSec() / (1E-38 * genie::units::cm2)); // Parse Interaction String std::string mode = genie_record.Summary()->AsString(); std::vector modevec = GeneralUtils::ParseToStr(mode, ";"); std::string targ = (modevec[0] + ";" + modevec[1]); std::string inter = mode; // Get target nucleus // Alternative ways of getting the summaries - //GHepParticle *target = genie_record.TargetNucleus(); - //int pdg = target->Pdg(); + // 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()) { + 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] = (TH1D *)xsechist->Clone(); + modecount[mode] = (TH1D *)xsechist->Clone(); + + modexsec[mode]->SetDirectory(NULL); + modecount[mode]->SetDirectory(NULL); - modexsec[mode]->GetYaxis()->SetTitle("d#sigma/dE_{#nu} #times 10^{-38} (events weighted by #sigma)"); + 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; + << " 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; + TFile *outputfile; if (!gOutputFile.length()) { - tn->GetEntry(0); - outputfile = tn->GetFile(); - outputfile->cd(); + // Shut the chain; + delete tn; + outputfile = new TFile(first_file.c_str(), "UPDATE"); } else { outputfile = new TFile(gOutputFile.c_str(), "RECREATE"); outputfile->cd(); QLOG(FIT, "Cloning input vector to output file: " << gOutputFile); - TTree* cloneTree = tn->CloneTree(); + 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(); + TTree *nova_tree = nova_chain->CloneTree(); if (!nova_tree) { QLOG(FIT, "Could not find nova_wgts tree in " << input); } else { QLOG(FIT, "Found nova_wgts tree in " << input); nova_tree->SetDirectory(outputfile); nova_tree->Write(); } QLOG(FIT, "Done cloning tree."); } LOG(FIT) << "Getting splines..." << std::endl; // Save each of the reconstructed splines to file - std::map modeavg; + std::map modeavg; - TDirectory* inddir = (TDirectory*)outputfile->Get("IndividualGENIESplines"); - if (!inddir) inddir = (TDirectory*)outputfile->mkdir("IndividualGENIESplines"); + 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] = (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"); + TDirectory *targdir = (TDirectory *)outputfile->Get("TargetGENIESplines"); + if (!targdir) + targdir = (TDirectory *)outputfile->mkdir("TargetGENIESplines"); targdir->cd(); LOG(FIT) << "Getting Target Splines" << std::endl; // For each target save a total spline - std::map targetsplines; + std::map targetsplines; for (uint i = 0; i < targetids.size(); i++) { std::string targ = targetids[i]; LOG(FIT) << "Getting target " << i << ": " << targ << std::endl; - targetsplines[targ] = (TH1D*)xsechist->Clone(); - targetsplines[targ]->GetYaxis()->SetTitle("#sigma (E_{#nu}) #times 10^{-38} (cm^{2}/target)"); + 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; + 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) << "Finished with Mode " << mode << " " + << modeavg[mode]->Integral() << std::endl; } } LOG(FIT) << "Saving target spline: " << targ << std::endl; targetsplines[targ]->Write(("Total_" + targ).c_str(), TObject::kOverwrite); } LOG(FIT) << "Getting total splines" << std::endl; // Now we have each of the targets we need to create a total cross-section. int totalnucl = 0; // This has structure target1[fraction1], target2[fraction2] std::vector targprs = GeneralUtils::ParseToStr(target, ","); std::vector targ_list; std::vector frac_list; - // Chop up the target string which has format TARGET1[fraction1],TARGET2[fraction2] + // Chop up the target string which has format + // TARGET1[fraction1],TARGET2[fraction2] - //std::cout << "Targets: " << std::endl; + // std::cout << "Targets: " << std::endl; // Loop over the vector of strings "TARGET1[fraction1]" "TARGET2[fraction2]" - for (std::vector::iterator it = targprs.begin(); it != targprs.end(); ++it) { + for (std::vector::iterator it = targprs.begin(); + it != targprs.end(); ++it) { // Cut into "TARGET1" and "fraction1]" std::vector targind = GeneralUtils::ParseToStr(*it, "["); - //std::cout << " " << *it << std::endl; + // std::cout << " " << *it << std::endl; // Cut into "TARGET1" and "fraction1" - for (std::vector::iterator jt = targind.begin(); jt != targind.end(); ++jt) { + for (std::vector::iterator jt = targind.begin(); + jt != targind.end(); ++jt) { if ((*jt).find("]") != std::string::npos) { (*jt) = (*jt).substr(0, (*jt).find("]")); //*jt = "hello"; frac_list.push_back(*jt); // Won't find bracket for target } else { targ_list.push_back(*jt); } } } targprs = targ_list; std::vector targ_fractions; double minimum = 1.0; - for (std::vector::iterator it = frac_list.begin(); it != frac_list.end(); it++) { - //std::cout << " " << *it << std::endl; + for (std::vector::iterator it = frac_list.begin(); + it != frac_list.end(); it++) { + // std::cout << " " << *it << std::endl; double frac = std::atof((*it).c_str()); targ_fractions.push_back(frac); - if (frac < minimum) minimum = frac; + if (frac < minimum) + minimum = frac; } std::vector::iterator it = targ_fractions.begin(); std::vector::iterator jt = targ_list.begin(); double scaling = 0; for (; it != targ_fractions.end(); it++, jt++) { // First get the mass number from the targ_list int nucl = atoi((*jt).c_str()); - nucl = (nucl%10000)/10; + nucl = (nucl % 10000) / 10; // Gets the relative portions right - *it = (*it)/minimum; + *it = (*it) / minimum; // Scale relative the atomic mass //(*it) *= (double(nucl)/(*it)); - double tempscaling = double(nucl)/(*it); - if (tempscaling > scaling) scaling=tempscaling; + 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); + *it *= int(scaling + 0.5); // Round to nearest integer - *it = int(*it+0.5); + *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); + 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(); + TH1D *totalxsec = (TH1D *)xsechist->Clone(); // Loop over the specified targets by the user for (uint i = 0; i < targprs.size(); i++) { std::string targpdg = targprs[i]; // Check that we found the user requested target in GENIE bool FoundTarget = false; - for (std::map::iterator iter = targetsplines.begin(); iter != targetsplines.end(); iter++) { + for (std::map::iterator iter = targetsplines.begin(); + iter != targetsplines.end(); iter++) { std::string targstr = iter->first; - TH1D* xsec = iter->second; + 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; + 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); - } + // 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) << "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; + for (uint i = 0; i < targprs.size(); ++i) + ERR(WRN) << " " << targprs[i] << std::endl; ERR(WRN) << " The list of targets found in GENIE is: " << std::endl; - for (std::map::iterator iter = targetsplines.begin(); iter != targetsplines.end(); iter++) ERR(WRN) << " " << iter->first<< std::endl; + for (std::map::iterator iter = targetsplines.begin(); + iter != targetsplines.end(); iter++) + ERR(WRN) << " " << iter->first << std::endl; } } outputfile->cd(); - totalxsec->GetYaxis()->SetTitle("#sigma (E_{#nu}) #times 10^{-38} (cm^{2}/nucleon)"); + totalxsec->GetYaxis()->SetTitle( + "#sigma (E_{#nu}) #times 10^{-38} (cm^{2}/nucleon)"); totalxsec->Write("nuisance_xsec", TObject::kOverwrite); - eventhist = (TH1D*)fluxhist->Clone(); + 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()); + 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) << "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; + << "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; + "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; + "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; + "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; + "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; + "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; + << 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[]) { +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; + << " " << 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) { + 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; + 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/nuiscomp.cxx b/app/nuiscomp.cxx index bea4c32..8d53d93 100644 --- a/app/nuiscomp.cxx +++ b/app/nuiscomp.cxx @@ -1,228 +1,226 @@ // 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 . *******************************************************************************/ #include "ComparisonRoutines.h" //******************************* void printInputCommands() { //******************************* std::cout << "nuiscomp : NUISANCE Data Comparison App \n" << std::endl; std::cout << "# Running nuiscomp with a card file #\n" << "######################################\n" << "nuiscomp -c cardfile.xml [ -o outputfile.root ] [ -f routines ] [ -n maxevents ] \n" << " [ -i 'cardstructure' ] [ -d fakedata ] [ -q config=val ] [ +e/-e ] [ +v/-v ] \n" << "\n" << "# Running nuiscomp with structures at cmd line #\n" << "################################################\n" << "nuiscomp -i 'cardstructure' -o outputfile.root [ -c cardfile.xml [ -f routines ] [ -n maxevents ] \n" << " [ -d fakedata ] [ -q config=val ] [ +e/-e ] [ +v/-v ] \n" << std::endl; sleep(4); std::cout << "" << "\n" << " \n" << " -c cardfile.xml : NUISANCE format card file defining comparisons. \n" << " \n" << " -o outputfile.root : Output file that histograms will be saved in. If a card file is \n" << " given but no output file this will default to cardfile.xml.root \n" << " \n" << " -f routines : Comma separated list of comparison routines to run in order. \n" << " Allowed Routines : \n" << " Compare : Fixed comparison at nominal dial values. \n" << " \n" << " -n maxevents : Set limit on the number of event entries to process. \n" << " \n" << " -i \'cardstructure\' : Define card structure like those available in the standard NUISANCE \n" << " card format, but on the command line at runtime. MUST be enclosed \n" << " in single quotation marks. See examples below for usage. \n" << " \n" << " It is possible to entirely define the comparison using \' -i\' commands \n" << " without the need to write a card file explicitly. If you do this, \n" << " make sure to also use the \' -o\' flag to tell it where to go. \n" << " \n" << " -d fakedata : Define a fake data set to be used. All data in NUISANCE will be set \n" << " to the values defined in this fake data before comparisons are made. \n" << " \n" << " There are two possible methods. Fake data from MC or a previous file.\n" << " fakedata = \'MC\' : Sets the MC to the values defined by \'fake_parameters\' \n" << " shown in the examples below, and then sets the data \n" << " to be equal to this MC prediction. \n" << " fakedata = \'file.root\' : Reads in the ROOT file at the specified path \n" << " assuming its a standard NUISANCE file. Takes \n" << " MC predictions in this file and uses them as \n" << " fake data. \n" << " \n" << " -q config=val : Overrides default configurations provided in the cardfile and in \n" << " '$NUISANCE/parameters/config.xml\'. Any config parameter can be set. \n" << " Examples : \n" << " \'-q VERBOSITY=4\' \n" << " \'-q EventManager=1\' \n" << " \'-q drawOpts=DATA/MC\' \n" << " \n" << " +e/-e : Increase/Decrease the default error verbosity by 1. \n" << " \n" << " +v/-v : Increase/Decrease the default logging verbosity by 1.\n" << " \n\n" << std::endl; sleep(4); std::cout << "# nuiscomp Running Examples #" << "############################# \n" << " \n" << " 1. Generate cardfile comparisons with increased verbosity and only 50000 events \n\n" << " nuiscomp -c cardfile.card -o mycomp.root -n 50000 +v +v \n" << " \n\n" << " 2. Generate a comparison to MiniBooNE data using simple structure, saving it to outfile.root \n\n" << " nuiscomp -o outfile.root -i \'sample MiniBooNE_CCQE_XSec_1DQ2_nu NEUT:neutevents.root\' \n" << " \n\n" << " 3. Generate a comparison to MiniBooNE data using xml structure, reweight MaCCQE, and save the prediction to outfile.root \n\n" << " nuiscomp -o outfile.root -i \'sample name=\"MiniBooNE_CCQE_XSec_1DQ2_nu\" input=\"NEUT:neutevents.root\"\' \\ \n" << " -i \'sample name=\"MiniBooNE_CC1pip_XSec_1DQ2_nu\" input=\"NEUT:neutevents.root\"\' \\ \n" << " -i \'parameter name=\"MaCCQE\" nominal=\"1.0\" type=\"neut_parameter\"\' \n " << " \n\n" << " 4. Generate a comparison, using fake data from the MC predictions inside the fakedata.root \n\n" << " nuiscomp -c cardfile.card -o myfakecomp.root -d fakedata.root \n" << " \n\n" << " 5. Generate a comparison using fake data defined on the command line use fake parameters \n\n" << " nuiscomp -c cardfile.card -d MC -i \'fakeparameter name=\"MaCCQE\" nominal=\"1.0\"\' \n " << " -i \'parameter name=\"MaCCQE\" nominal=\"1.0\" type=\"neut_parameter\"' " << " \n\n" << std::endl; sleep(4); std::cout << "# NUISANCE Card Format Structure Examples # \n" << "########################################### \n" << "\n" << "The NUISANCE card can be defined as a simple text file, or an xml file. \n" << "Examples for both with relevant structures are given below. \n" << std::endl; std::cout << "# XML Card File Example # \n" << "cardfile.xml: \n" << "" << "\n" << " \n" << " \n" << " \n" << "\n" << " \n" << " \n" << " \n" << " \n" << " \n" << " \n" << " \n" << "\n" << " \n" << " \n" << "\n" << "\n" << " \n" << " \n" << " \n" << " \n" << " \n" << "\n" << " \n" << "\n" << " \n" << " \n" << " \n" << " \n" << " \n" << "\n" << " \n" << "\n" << " \n\n" << std::endl; std::cout << "# Simple Card File Example # \n" << "cardfile.card: \n" << "\n" << "# CONFIG STRUCTURE \n" << "# config name val \n" << "config VERBOSITY 4 \n" << "\n" << "# Sample Structure \n" << "# ID Corresponds to names given in src/FCN/SampleList.cxx \n" << "# TYPE is the generator type (NEUT,NUWRO,GENIE,GIBUU). \n" << "# FILE is the input generator events file. \n" << "# TYPE is optional and used to define options for a class. e.g. FREE \n" << "# NORM is optional and sets sample normalisations. \n" << "# sample ID TYPE:FILE TYPE 1.0 \n" << "\n" << "sample MiniBooNE_CCQE_XSec_1DQ2_nu GENIE:genieevents.root \n" << "sample MiniBooNE_CC1pip_XSec_1DQ2_nu GENIE:genieevents.root SHAPE \n" << "\n" << "\n" << "# Parameter Structure \n" << "# ID is the name of the dial in each generator RW engine \n" << "# TYPE is the dial type (neut,newer,genie,niwg,t2k,custom,norm) \n" << "# VAL is the nominal value in 1-sigma variation for the comparison \n" << "# TYPE_parameter ID VAL \n" << "\n" << "neut_parameter MaCCQE 0.5 \n" << "\n" << "# Fake Parameter Structure \n" << "# Sets values for fake data defined using the ‘MC’ flag. \n" << "# ID is the dial name, it MUST be specified before hand using a normal parameter structure \n" << "# VAL is the value to use for the fake data \n" << "# fake_parameter ID VAL \n" << "\n" << "fake_parameter MaCCQE 1.0 \n" << "\n" << std::endl; 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 nuiscomp.exe" << std::endl; // Make minimizer class and run fit ComparisonRoutines* comp = new ComparisonRoutines(argc, argv); comp->Run(); + delete comp; // Show Final Status LOG(FIT) << "------------------------------------ -" << std::endl; LOG(FIT) << "Comparison Complete." << std::endl; LOG(FIT) << "------------------------------------ -" << std::endl; return status; } - - - diff --git a/cmake/GENIESetup.cmake b/cmake/GENIESetup.cmake index 74dfd3e..c7f9eb3 100644 --- a/cmake/GENIESetup.cmake +++ b/cmake/GENIESetup.cmake @@ -1,236 +1,232 @@ # 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 . ################################################################################ # TODO # check system for libxml2 # check whether we need the includes # check if we can use a subset of the GENIE libraries include(${CMAKE_SOURCE_DIR}/cmake/parseConfigApp.cmake) ################################################################################ # 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 an environment variable" " $ 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 --topsrcdir OUTPUT_VARIABLE GENIE_INCLUDES_DIR OUTPUT_STRIP_TRAILING_WHITESPACE) #Allows for external override in the case where genie-config lies. if(NOT DEFINED GENIE_LIB_DIR OR GENIE_LIB_DIR STREQUAL "") GETLIBDIRS(genie-config --libs GENIE_LIB_DIR) endif() GETLIBS(genie-config --libs GENIE_LIBS) cmessage(STATUS "GENIE version : ${GENIE_VERSION}") cmessage(STATUS "GENIE libdir : ${GENIE_LIB_DIR}") cmessage(STATUS "GENIE libs : ${GENIE_LIBS}") string(REGEX MATCH "ReinSeghal" WASMATCHED ${GENIE_LIBS}) if(WASMATCHED AND GENIE_VERSION STREQUAL "210") set(GENIE_SEHGAL ${GENIE_LIBS}) STRING(REPLACE "ReinSeghal" "ReinSehgal" GENIE_LIBS ${GENIE_SEHGAL}) cmessage(DEBUG "Fixed inconsistency in library naming: ${GENIE_LIBS}") endif() if(NOT GENIE_POST_R3) LIST(FIND GENIE_LIBS GReWeight WAS_FOUND) if(WAS_FOUND STREQUAL "-1") LIST(APPEND GENIE_LIBS GReWeight) cmessage(DEBUG "Force added ReWeight library: ${GENIE_LIBS}") endif() else() LIST(FIND GENIE_LIBS GRwFwk WAS_FOUND) if(WAS_FOUND STREQUAL "-1") LIST(APPEND GENIE_LIBS GRwClc GRwFwk GRwIO) cmessage(DEBUG "Force added ReWeight library: ${GENIE_LIBS}") endif() endif() LIST(APPEND GENIE_LIBS -Wl,--end-group ) LIST(REVERSE GENIE_LIBS) LIST(APPEND GENIE_LIBS -Wl,--start-group -Wl,--no-as-needed ) LIST(REVERSE GENIE_LIBS) cmessage(DEBUG "GENIE_LIBS: ${GENIE_LIBS}") ################################ 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 an environment variable $ 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 an environment variable $ 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 an environment variable $ export LHAPATH=/path/to/LHAPATH") endif() ################################ LIBXML ###################################### if(LIBXML2_LIB STREQUAL "") GETLIBDIR(xml2-config --libs LIBXML2_LIB IGNORE_EMPTY_RESPONSE) if(LIBXML2_LIB STREQUAL "") message(WARNING "Variable LIBXML2_LIB is not defined, as xml2-config was found and didn't report a library include path, it is likely that libxml2.so can be found in the standard system location, lets hope so.") endif() endif() if(LIBXML2_INC STREQUAL "") GETINCDIR(xml2-config --cflags LIBXML2_INC) if(LIBXML2_INC STREQUAL "") message(FATAL_ERROR "Variable LIBXML2_INC is not defined and could not be found with xml2-config. The location of a pre-built libxml2 install must be defined either as $ cmake -DLIBXML2_INC=/path/to/LIBXML2_includes or as an environment variable $ export LIBXML2_INC=/path/to/LIBXML2_includes") endif() endif() ############################### log4cpp ###################################### if(LOG4CPP_LIB STREQUAL "") GETLIBDIR(log4cpp-config --libs LOG4CPP_LIB) if(LOG4CPP_LIB STREQUAL "") message(FATAL_ERROR "Variable LOG4CPP_LIB is not defined and could not be found with log4cpp-config. The location of a pre-built log4cpp install must be defined either as $ cmake -DLOG4CPP_LIB=/path/to/LOG4CPP_libraries or as an environment variable $ export LOG4CPP_LIB=/path/to/LOG4CPP_libraries") endif() endif() if(LOG4CPP_INC STREQUAL "") GETINCDIR(log4cpp-config --cflags LOG4CPP_INC) if(LOG4CPP_INC STREQUAL "") message(FATAL_ERROR "Variable LOG4CPP_INC is not defined and could not be found with log4cpp-config. The location of a pre-built log4cpp install must be defined either as $ cmake -DLOG4CPP_INC=/path/to/LOG4CPP_includes or as an environment variable $ export LOG4CPP_INC=/path/to/LOG4CPP_includes") endif() endif() ################################################################################ LIST(APPEND EXTRA_CXX_FLAGS -D__GENIE_ENABLED__ -D__GENIE_VERSION__=${GENIE_VERSION}) LIST(APPEND EXTRA_LIBS ${GENIE_LIBS}) ############################### GSL ###################################### if(GENIE_POST_R3) if(GSL_LIB STREQUAL "") GETLIBDIR(gsl-config --libs GSL_LIB) if(GSL_LIB STREQUAL "") message(FATAL_ERROR "Variable GSL_LIB is not defined and could not be found with gsl-config. The location of a pre-built gsl install must be defined either as $ cmake -DGSL_LIB=/path/to/GSL_libraries or as an environment variable $ export GSL_LIB=/path/to/GSL_libraries") endif() endif() if(GSL_INC STREQUAL "") GETINCDIR(gsl-config --cflags GSL_INC) if(GSL_INC STREQUAL "") message(FATAL_ERROR "Variable GSL_INC is not defined and could not be found with gsl-config. The location of a pre-built gsl install must be defined either as $ cmake -DGSL_INC=/path/to/GSL_includes or as an environment variable $ export GSL_INC=/path/to/GSL_includes") endif() endif() - #cmessage(WARNING "GSL_LIB: ${GSL_LIB}, GSL_LIB_LIST: ${GSL_LIB_LIST}") GETLIBS(gsl-config --libs GSL_LIB_LIST) - #cmessage(AUTHOR_WARNING "GSL_LIB: ${GSL_LIB}, GSL_LIB_LIST: ${GSL_LIB_LIST}") if(GENIE_REWEIGHT STREQUAL "") message(FATAL_ERROR "Variable GENIE_REWEIGHT is not defined. When using GENIE v3+, we require the reweight product to be built and accessible via the environment variable GENIE_REWEIGHT") endif() endif() ################################################################################ -#cmessage(FATAL_ERROR "GSL_LIB: ${GSL_LIB}, GSL_LIB_LIST: ${GSL_LIB_LIST}") - LIST(APPEND EXTRA_LIBS LHAPDF xml2 log4cpp) LIST(APPEND EXTRA_LINK_DIRS ${GENIE_LIB_DIR} ${LHAPDF_LIB} ${LOG4CPP_LIB}) 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} ${GENIE_REWEIGHT}/src ${GSL_INC} ${LHAPDF_INC} ${LIBXML2_INC} ${LOG4CPP_INC}) LIST(APPEND EXTRA_LINK_DIRS ${GENIE_REWEIGHT}/lib ${GSL_LIB} ) LIST(APPEND EXTRA_LIBS ${GSL_LIB_LIST}) endif() cmessage(WARNING ${EXTRA_LINK_DIRS}) cmessage(WARNING ${EXTRA_LIBS}) SAYVARS() if(USE_PYTHIA8) set(NEED_PYTHIA8 TRUE) set(NEED_ROOTPYTHIA8 TRUE) else() set(NEED_PYTHIA6 TRUE) set(NEED_ROOTPYTHIA6 TRUE) endif() set(NEED_ROOTEVEGEN TRUE) SET(USE_GENIE TRUE CACHE BOOL "Whether to enable GENIE (reweight) support. Requires external libraries. " FORCE) diff --git a/cmake/setup.sh.in b/cmake/setup.sh.in index 123aaef..9956d8a 100644 --- a/cmake/setup.sh.in +++ b/cmake/setup.sh.in @@ -1,155 +1,160 @@ # 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 . ################################################################################ #!/bin/sh ### Adapted from https://unix.stackexchange.com/questions/4965/keep-duplicates-out-of-path-on-source function add_to_PATH () { for d; do d=$(cd -- "$d" && { pwd -P || pwd; }) 2>/dev/null # canonicalize symbolic links if [ -z "$d" ]; then continue; fi # skip nonexistent directory if [ "$d" == "/usr/bin" ] || [ "$d" == "/usr/bin64" ] || [ "$d" == "/usr/local/bin" ] || [ "$d" == "/usr/local/bin64" ]; then case ":$PATH:" in *":$d:"*) :;; *) export PATH=$PATH:$d;; esac else case ":$PATH:" in *":$d:"*) :;; *) export PATH=$d:$PATH;; esac fi done } function add_to_LD_LIBRARY_PATH () { for d; do d=$(cd -- "$d" && { pwd -P || pwd; }) 2>/dev/null # canonicalize symbolic links if [ -z "$d" ]; then continue; fi # skip nonexistent directory if [ "$d" == "/usr/lib" ] || [ "$d" == "/usr/lib64" ] || [ "$d" == "/usr/local/lib" ] || [ "$d" == "/usr/local/lib64" ]; then case ":$LD_LIBRARY_PATH:" in *":$d:"*) :;; *) export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$d;; esac else case ":$LD_LIBRARY_PATH:" in *":$d:"*) :;; *) export LD_LIBRARY_PATH=$d:$LD_LIBRARY_PATH;; esac fi done } if [ "@EXTRA_SETUP_SCRIPT@" ]; then if [ ! -e @EXTRA_SETUP_SCRIPT@ ]; then echo "[WARN]: Extra setup script \"@EXTRA_SETUP_SCRIPT@\" requested, but could not be found. Skipping..." else echo "[INFO]: Sourcing extra setup from \"@EXTRA_SETUP_SCRIPT@\"." . @EXTRA_SETUP_SCRIPT@ fi fi add_to_PATH "@CMAKE_INSTALL_PREFIX@/bin" add_to_LD_LIBRARY_PATH "@CMAKE_INSTALL_PREFIX@/lib" if [ ! "${ROOTSYS}" ]; then echo "[INFO]: Sourcing ROOT from: @CMAKE_ROOTSYS@" source "@CMAKE_ROOTSYS@/bin/thisroot.sh" fi if [ "@USE_T2K@" != "FALSE" ]; then echo "[INFO]: Adding T2K paths to the environment." export T2KREWEIGHT=@T2KREWEIGHT@ add_to_LD_LIBRARY_PATH "@T2KREWEIGHT@/lib" fi if [ "@USE_NIWG@" != "FALSE" ]; then echo "[INFO]: Adding NIWG paths to the environment." export NIWG=@NIWG_ROOT@ export NIWGREWEIGHT_INPUTS=@NIWG_ROOT@/inputs add_to_LD_LIBRARY_PATH "@NIWG_ROOT@" fi if [ "@USE_NEUT@" != "FALSE" ]; then echo "[INFO]: Adding NEUT library paths to the environment." export NEUT_ROOT=@NEUT_ROOT@ export CERN=@CERN@ export CERN_LEVEL=@CERN_LEVEL@ add_to_LD_LIBRARY_PATH "${NEUT_LIB_DIR}" "${NEUT_ROOT}/src/reweight" fi if [ "@USE_NuWro@" != "FALSE" ]; then if [ "@NUWRO_BUILT_FROM_FILE@" == "FALSE" ]; then echo "[INFO]: Adding NuWro library paths to the environment." export NUWRO="@NUWRO@" add_to_PATH "@NUWRO@/bin" add_to_LD_LIBRARY_PATH "@NUWRO@/build/@CMAKE_SYSTEM_NAME@/lib" if [ "@NUWRO_INC@" ]; then export NUWRO_INC=@NUWRO_INC@ fi else echo "[INFO]: NuWro support included from input event file." fi fi if [ "@NEED_PYTHIA6@" != "FALSE" ]; then echo "[INFO]: Adding PYTHIA6 library paths to the environment." export PYTHIA6="@PYTHIA6@" add_to_LD_LIBRARY_PATH "@PYTHIA6@" fi if [ "@USE_GENIE@" != "FALSE" ]; then echo "[INFO]: Adding GENIE paths to the environment." export GENIE="@GENIE@" export LHAPDF_LIB="@LHAPDF_LIB@" export LHAPDF_INC="@LHAPDF_INC@" export LIBXML2_LIB="@LIBXML2_LIB@" export LIBXML2_INC="@LIBXML2_INC@" export LOG4CPP_LIB="@LOG4CPP_LIB@" export LOG4CPP_INC="@LOG4CPP_INC@" if [ "@LHAPATH@" ]; then export LHAPATH="@LHAPATH@" fi add_to_PATH "@GENIE@/bin" add_to_LD_LIBRARY_PATH "@GENIE@/lib" "@LHAPDF_LIB@" "@LIBXML2_LIB@" "@LOG4CPP_LIB@" + if [ "@GENIE_REWEIGHT@" ]; then + export GENIE_REWEIGHT="@GENIE_REWEIGHT@" + add_to_LD_LIBRARY_PATH "@GENIE_REWEIGHT@/lib" + fi + fi if [ "@BUILD_GiBUU@" != "FALSE" ]; then echo "[INFO]: Sourcing GiBUU tools." source @CMAKE_BINARY_DIR@/GiBUUTools/src/GiBUUTools-build/Linux/setup.sh fi export NUISANCE="@CMAKE_SOURCE_DIR@" diff --git a/src/InputHandler/GENIEInputHandler.cxx b/src/InputHandler/GENIEInputHandler.cxx index 78ff0e3..600d076 100644 --- a/src/InputHandler/GENIEInputHandler.cxx +++ b/src/InputHandler/GENIEInputHandler.cxx @@ -1,508 +1,523 @@ // 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 . *******************************************************************************/ #ifdef __GENIE_ENABLED__ #include "GENIEInputHandler.h" +#pragma push_macro("LOG") +#undef LOG +#pragma push_macro("ERROR") +#undef ERROR +#ifdef GENIE_PRE_R3 +#include "Messenger/Messenger.h" +#else +#include "Framework/Messenger/Messenger.h" +#endif +#pragma pop_macro("LOG") +#pragma pop_macro("ERROR") + + #include "InputUtils.h" GENIEGeneratorInfo::~GENIEGeneratorInfo() { DeallocateParticleStack(); } void GENIEGeneratorInfo::AddBranchesToTree(TTree* tn) { tn->Branch("GenieParticlePDGs", &fGenieParticlePDGs, "GenieParticlePDGs/I"); } void GENIEGeneratorInfo::SetBranchesFromTree(TTree* tn) { tn->SetBranchAddress("GenieParticlePDGs", &fGenieParticlePDGs); } void GENIEGeneratorInfo::AllocateParticleStack(int stacksize) { fGenieParticlePDGs = new int[stacksize]; } void GENIEGeneratorInfo::DeallocateParticleStack() { delete fGenieParticlePDGs; } void GENIEGeneratorInfo::FillGeneratorInfo(NtpMCEventRecord* ntpl) { Reset(); // Check for GENIE Event if (!ntpl) return; if (!ntpl->event) return; // Cast Event Record GHepRecord* ghep = static_cast(ntpl->event); if (!ghep) return; // Fill Particle Stack GHepParticle* p = 0; TObjArrayIter iter(ghep); // Loop over all particles int i = 0; while ((p = (dynamic_cast((iter).Next())))) { if (!p) continue; // Get PDG fGenieParticlePDGs[i] = p->Pdg(); i++; } } void GENIEGeneratorInfo::Reset() { for (int i = 0; i < kMaxParticles; i++) { fGenieParticlePDGs[i] = 0; } } GENIEInputHandler::GENIEInputHandler(std::string const& handle, std::string const& rawinputs) { LOG(SAM) << "Creating GENIEInputHandler : " << handle << std::endl; + genie::Messenger::Instance()->SetPriorityLevel("GHepUtils", pFATAL); + // Run a joint input handling fName = handle; // Setup the TChain fGENIETree = new TChain("gtree"); fSaveExtra = FitPar::Config().GetParB("SaveExtraGenie"); fCacheSize = FitPar::Config().GetParI("CacheSize"); fMaxEvents = FitPar::Config().GetParI("MAXEVENTS"); // Are we running with NOvA weights fNOvAWeights = FitPar::Config().GetParB("NOvA_Weights"); MAQEw = 1.0; NonResw = 1.0; RPAQEw = 1.0; RPARESw = 1.0; MECw = 1.0; DISw = 1.0; NOVAw = 1.0; // Loop over all inputs and grab flux, eventhist, and nevents std::vector inputs = InputUtils::ParseInputFileList(rawinputs); for (size_t inp_it = 0; inp_it < inputs.size(); ++inp_it) { // Open File for histogram access TFile* inp_file = new TFile( InputUtils::ExpandInputDirectories(inputs[inp_it]).c_str(), "READ"); if (!inp_file or inp_file->IsZombie()) { THROW("GENIE File IsZombie() at : '" << inputs[inp_it] << "'" << std::endl << "Check that your file paths are correct and the file exists!" << std::endl << "$ ls -lh " << inputs[inp_it]); } // Get Flux/Event hist TH1D* fluxhist = (TH1D*)inp_file->Get("nuisance_flux"); TH1D* eventhist = (TH1D*)inp_file->Get("nuisance_events"); if (!fluxhist or !eventhist) { ERROR(FTL, "Input File Contents: " << inputs[inp_it]); inp_file->ls(); THROW("GENIE FILE doesn't contain flux/xsec info." << std::endl << "Try running the app PrepareGENIE first on :" << inputs[inp_it] << std::endl << "$ PrepareGENIE -h"); } // Get N Events TTree* genietree = (TTree*)inp_file->Get("gtree"); if (!genietree) { ERROR(FTL, "gtree not located in GENIE file: " << inputs[inp_it]); THROW("Check your inputs, they may need to be completely regenerated!"); throw; } - + int nevents = genietree->GetEntries(); if (nevents <= 0) { THROW("Trying to a TTree with " << nevents << " to TChain from : " << inputs[inp_it]); } // Check for precomputed weights TTree *weighttree = (TTree*)inp_file->Get("nova_wgts"); if (fNOvAWeights) { if (!weighttree) { THROW("Did not find nova_wgts tree in file " << inputs[inp_it] << " but you specified it" << std::endl); } else { LOG(FIT) << "Found nova_wgts tree in file " << inputs[inp_it] << std::endl; } } // Register input to form flux/event rate hists RegisterJointInput(inputs[inp_it], nevents, fluxhist, eventhist); // Add To TChain fGENIETree->AddFile(inputs[inp_it].c_str()); if (weighttree != NULL) fGENIETree->AddFriend(weighttree); } // Registor all our file inputs SetupJointInputs(); // Assign to tree fEventType = kGENIE; fGenieNtpl = NULL; fGENIETree->SetBranchAddress("gmcrec", &fGenieNtpl); // Set up the custom weights if (fNOvAWeights) { fGENIETree->SetBranchAddress("MAQEwgt", &MAQEw); fGENIETree->SetBranchAddress("nonResNormWgt", &NonResw); fGENIETree->SetBranchAddress("RPAQEWgt", &RPAQEw); fGENIETree->SetBranchAddress("RPARESWgt", &RPARESw); fGENIETree->SetBranchAddress("MECWgt", &MECw); fGENIETree->SetBranchAddress("DISWgt", &DISw); fGENIETree->SetBranchAddress("nova2018CVWgt", &NOVAw); } // Libraries should be seen but not heard... StopTalking(); fGENIETree->GetEntry(0); StartTalking(); // Create Fit Event fNUISANCEEvent = new FitEvent(); fNUISANCEEvent->SetGenieEvent(fGenieNtpl); if (fSaveExtra) { fGenieInfo = new GENIEGeneratorInfo(); fNUISANCEEvent->AddGeneratorInfo(fGenieInfo); } fNUISANCEEvent->HardReset(); }; GENIEInputHandler::~GENIEInputHandler() { // if (fGenieGHep) delete fGenieGHep; // if (fGenieNtpl) delete fGenieNtpl; // if (fGENIETree) delete fGENIETree; // if (fGenieInfo) delete fGenieInfo; } void GENIEInputHandler::CreateCache() { if (fCacheSize > 0) { // fGENIETree->SetCacheEntryRange(0, fNEvents); fGENIETree->AddBranchToCache("*", 1); fGENIETree->SetCacheSize(fCacheSize); } } void GENIEInputHandler::RemoveCache() { // fGENIETree->SetCacheEntryRange(0, fNEvents); fGENIETree->AddBranchToCache("*", 0); fGENIETree->SetCacheSize(0); } FitEvent* GENIEInputHandler::GetNuisanceEvent(const UInt_t entry, const bool lightweight) { if (entry >= (UInt_t)fNEvents) return NULL; // Clear the previous event (See Note 1 in ROOT TClonesArray documentation) if (fGenieNtpl) { fGenieNtpl->Clear(); } // Read Entry from TTree to fill NEUT Vect in BaseFitEvt; fGENIETree->GetEntry(entry); // Run NUISANCE Vector Filler if (!lightweight) { CalcNUISANCEKinematics(); } #ifdef __PROB3PP_ENABLED__ else { // Check for GENIE Event if (!fGenieNtpl) return NULL; if (!fGenieNtpl->event) return NULL; // Cast Event Record fGenieGHep = static_cast(fGenieNtpl->event); if (!fGenieGHep) return NULL; TObjArrayIter iter(fGenieGHep); genie::GHepParticle* p; while ((p = (dynamic_cast((iter).Next())))) { if (!p) { continue; } // Get Status int state = GetGENIEParticleStatus(p, fNUISANCEEvent->Mode); if (state != genie::kIStInitialState) { continue; } fNUISANCEEvent->probe_E = p->E() * 1.E3; fNUISANCEEvent->probe_pdg = p->Pdg(); break; } } #endif // Setup Input scaling for joint inputs fNUISANCEEvent->InputWeight = GetInputWeight(entry); return fNUISANCEEvent; } int GENIEInputHandler::GetGENIEParticleStatus(genie::GHepParticle* p, int mode) { /* kIStUndefined = -1, kIStInitialState = 0, / generator-level initial state / kIStStableFinalState = 1, / generator-level final state: particles to be tracked by detector-level MC / kIStIntermediateState = 2, kIStDecayedState = 3, kIStCorrelatedNucleon = 10, kIStNucleonTarget = 11, kIStDISPreFragmHadronicState = 12, kIStPreDecayResonantState = 13, kIStHadronInTheNucleus = 14, / hadrons inside the nucleus: marked for hadron transport modules to act on / kIStFinalStateNuclearRemnant = 15, / low energy nuclear fragments entering the record collectively as a 'hadronic blob' pseudo-particle / kIStNucleonClusterTarget = 16, // for composite nucleons before phase space decay */ int state = kUndefinedState; switch (p->Status()) { case genie::kIStNucleonTarget: case genie::kIStInitialState: case genie::kIStCorrelatedNucleon: case genie::kIStNucleonClusterTarget: state = kInitialState; break; case genie::kIStStableFinalState: state = kFinalState; break; case genie::kIStHadronInTheNucleus: if (abs(mode) == 2) state = kInitialState; else state = kFSIState; break; case genie::kIStPreDecayResonantState: case genie::kIStDISPreFragmHadronicState: case genie::kIStIntermediateState: state = kFSIState; break; case genie::kIStFinalStateNuclearRemnant: case genie::kIStUndefined: case genie::kIStDecayedState: default: break; } // Flag to remove nuclear part in genie if (p->Pdg() > 1000000) { if (state == kInitialState) state = kNuclearInitial; else if (state == kFinalState) state = kNuclearRemnant; } return state; } #endif #ifdef __GENIE_ENABLED__ int GENIEInputHandler::ConvertGENIEReactionCode(GHepRecord* gheprec) { // Electron Scattering if (gheprec->Summary()->ProcInfo().IsEM()) { if (gheprec->Summary()->InitState().ProbePdg() == 11) { if (gheprec->Summary()->ProcInfo().IsQuasiElastic()) return 1; else if (gheprec->Summary()->ProcInfo().IsMEC()) return 2; else if (gheprec->Summary()->ProcInfo().IsResonant()) return 13; else if (gheprec->Summary()->ProcInfo().IsDeepInelastic()) return 26; else { ERROR(WRN, "Unknown GENIE Electron Scattering Mode!" << std::endl << "ScatteringTypeId = " << gheprec->Summary()->ProcInfo().ScatteringTypeId() << " " << "InteractionTypeId = " << gheprec->Summary()->ProcInfo().InteractionTypeId() << std::endl << genie::ScatteringType::AsString( gheprec->Summary()->ProcInfo().ScatteringTypeId()) << " " << genie::InteractionType::AsString( gheprec->Summary()->ProcInfo().InteractionTypeId()) << " " << gheprec->Summary()->ProcInfo().IsMEC()); return 0; } } // Weak CC } else if (gheprec->Summary()->ProcInfo().IsWeakCC()) { // CC MEC if (gheprec->Summary()->ProcInfo().IsMEC()) { if (pdg::IsNeutrino(gheprec->Summary()->InitState().ProbePdg())) return 2; else if (pdg::IsAntiNeutrino(gheprec->Summary()->InitState().ProbePdg())) return -2; // CC OTHER } else { return utils::ghep::NeutReactionCode(gheprec); } // Weak NC } else if (gheprec->Summary()->ProcInfo().IsWeakNC()) { // NC MEC if (gheprec->Summary()->ProcInfo().IsMEC()) { if (pdg::IsNeutrino(gheprec->Summary()->InitState().ProbePdg())) return 32; else if (pdg::IsAntiNeutrino(gheprec->Summary()->InitState().ProbePdg())) return -32; // NC OTHER } else { return utils::ghep::NeutReactionCode(gheprec); } } return 0; } void GENIEInputHandler::CalcNUISANCEKinematics() { // Reset all variables fNUISANCEEvent->ResetEvent(); // Check for GENIE Event if (!fGenieNtpl) return; if (!fGenieNtpl->event) return; // Cast Event Record fGenieGHep = static_cast(fGenieNtpl->event); if (!fGenieGHep) return; // Convert GENIE Reaction Code fNUISANCEEvent->Mode = ConvertGENIEReactionCode(fGenieGHep); // Set Event Info fNUISANCEEvent->fEventNo = 0.0; fNUISANCEEvent->fTotCrs = fGenieGHep->XSec(); fNUISANCEEvent->fTargetA = 0.0; fNUISANCEEvent->fTargetZ = 0.0; fNUISANCEEvent->fTargetH = 0; fNUISANCEEvent->fBound = 0.0; fNUISANCEEvent->InputWeight = 1.0; //(1E+38 / genie::units::cm2) * fGenieGHep->XSec(); // And the custom weights if (fNOvAWeights) { fNUISANCEEvent->CustomWeight = NOVAw; fNUISANCEEvent->CustomWeightArray[0] = MAQEw; fNUISANCEEvent->CustomWeightArray[1] = NonResw; fNUISANCEEvent->CustomWeightArray[2] = RPAQEw; fNUISANCEEvent->CustomWeightArray[3] = RPARESw; fNUISANCEEvent->CustomWeightArray[4] = MECw; fNUISANCEEvent->CustomWeightArray[5] = NOVAw; } else { fNUISANCEEvent->CustomWeight = 1.0; fNUISANCEEvent->CustomWeightArray[0] = 1.0; fNUISANCEEvent->CustomWeightArray[1] = 1.0; fNUISANCEEvent->CustomWeightArray[2] = 1.0; fNUISANCEEvent->CustomWeightArray[3] = 1.0; fNUISANCEEvent->CustomWeightArray[4] = 1.0; fNUISANCEEvent->CustomWeightArray[5] = 1.0; } // Get N Particle Stack unsigned int npart = fGenieGHep->GetEntries(); unsigned int kmax = fNUISANCEEvent->kMaxParticles; if (npart > kmax) { ERR(WRN) << "GENIE has too many particles, expanding stack." << std::endl; fNUISANCEEvent->ExpandParticleStack(npart); } // Fill Particle Stack GHepParticle* p = 0; TObjArrayIter iter(fGenieGHep); fNUISANCEEvent->fNParticles = 0; // Loop over all particles while ((p = (dynamic_cast((iter).Next())))) { if (!p) continue; // Get Status int state = GetGENIEParticleStatus(p, fNUISANCEEvent->Mode); // Remove Undefined if (kRemoveUndefParticles && state == kUndefinedState) continue; // Remove FSI if (kRemoveFSIParticles && state == kFSIState) continue; if (kRemoveNuclearParticles && (state == kNuclearInitial || state == kNuclearRemnant)) continue; // Fill Vectors int curpart = fNUISANCEEvent->fNParticles; fNUISANCEEvent->fParticleState[curpart] = state; // Mom fNUISANCEEvent->fParticleMom[curpart][0] = p->Px() * 1.E3; fNUISANCEEvent->fParticleMom[curpart][1] = p->Py() * 1.E3; fNUISANCEEvent->fParticleMom[curpart][2] = p->Pz() * 1.E3; fNUISANCEEvent->fParticleMom[curpart][3] = p->E() * 1.E3; // PDG fNUISANCEEvent->fParticlePDG[curpart] = p->Pdg(); // Add to N particle count fNUISANCEEvent->fNParticles++; // Extra Check incase GENIE fails. if ((UInt_t)fNUISANCEEvent->fNParticles == kmax) { ERR(WRN) << "Number of GENIE Particles exceeds maximum!" << std::endl; ERR(WRN) << "Extend kMax, or run without including FSI particles!" << std::endl; break; } } // Fill Extra Stack if (fSaveExtra) fGenieInfo->FillGeneratorInfo(fGenieNtpl); // Run Initial, FSI, Final, Other ordering. fNUISANCEEvent->OrderStack(); FitParticle* ISNeutralLepton = fNUISANCEEvent->GetHMISParticle(PhysConst::pdg_neutrinos); if (ISNeutralLepton) { fNUISANCEEvent->probe_E = ISNeutralLepton->E(); fNUISANCEEvent->probe_pdg = ISNeutralLepton->PDG(); } return; } void GENIEInputHandler::Print() {} #endif diff --git a/src/Routines/ComparisonRoutines.cxx b/src/Routines/ComparisonRoutines.cxx index 1f4f7bb..a4adcd6 100755 --- a/src/Routines/ComparisonRoutines.cxx +++ b/src/Routines/ComparisonRoutines.cxx @@ -1,532 +1,531 @@ // 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 . *******************************************************************************/ #include "ComparisonRoutines.h" /* Constructor/Destructor */ //************************ void ComparisonRoutines::Init() { //************************ fOutputFile = ""; fOutputRootFile = NULL; fStrategy = "Compare"; fRoutines.clear(); fCardFile = ""; fFakeDataInput = ""; fSampleFCN = NULL; fAllowedRoutines = ("Compare"); }; //************************************* ComparisonRoutines::~ComparisonRoutines() { //************************************* +delete fOutputRootFile; }; /* Input Functions */ //************************************* ComparisonRoutines::ComparisonRoutines(int argc, char* argv[]) { //************************************* // Initialise Defaults Init(); nuisconfig configuration = Config::Get(); // Default containers std::string cardfile = ""; std::string maxevents = "-1"; int errorcount = 0; int verbocount = 0; std::vector xmlcmds; std::vector configargs; // Make easier to handle arguments. std::vector args = GeneralUtils::LoadCharToVectStr(argc, argv); ParserUtils::ParseArgument(args, "-c", fCardFile, true); ParserUtils::ParseArgument(args, "-o", fOutputFile, false, false); ParserUtils::ParseArgument(args, "-n", maxevents, false, false); ParserUtils::ParseArgument(args, "-f", fStrategy, false, false); ParserUtils::ParseArgument(args, "-d", fFakeDataInput, false, false); ParserUtils::ParseArgument(args, "-i", xmlcmds); ParserUtils::ParseArgument(args, "-q", configargs); ParserUtils::ParseCounter(args, "e", errorcount); ParserUtils::ParseCounter(args, "v", verbocount); ParserUtils::CheckBadArguments(args); // Add extra defaults if none given if (fCardFile.empty() and xmlcmds.empty()) { ERR(FTL) << "No input supplied!" << std::endl; throw; } if (fOutputFile.empty() and !fCardFile.empty()) { fOutputFile = fCardFile + ".root"; ERR(WRN) << "No output supplied so saving it to: " << fOutputFile << std::endl; } else if (fOutputFile.empty()) { ERR(FTL) << "No output file or cardfile supplied!" << std::endl; throw; } // Configuration Setup ============================= // Check no comp key is available nuiskey fCompKey; if (Config::Get().GetNodes("nuiscomp").empty()) { fCompKey = Config::Get().CreateNode("nuiscomp"); } else { fCompKey = Config::Get().GetNodes("nuiscomp")[0]; } if (!fCardFile.empty()) fCompKey.Set("cardfile", fCardFile); if (!fOutputFile.empty()) fCompKey.Set("outputfile", fOutputFile); if (!fStrategy.empty()) fCompKey.Set("strategy", fStrategy); // Load XML Cardfile configuration.LoadSettings( fCompKey.GetS("cardfile"), ""); // Add Config Args for (size_t i = 0; i < configargs.size(); i++) { configuration.OverrideConfig(configargs[i]); } if (maxevents.compare("-1")){ configuration.OverrideConfig("MAXEVENTS=" + maxevents); } // Finish configuration XML configuration.FinaliseSettings(fCompKey.GetS("outputfile") + ".xml"); // Add Error Verbo Lines verbocount += Config::GetParI("VERBOSITY"); errorcount += Config::GetParI("ERROR"); bool trace = Config::GetParB("TRACE"); std::cout << "[ NUISANCE ]: Setting VERBOSITY=" << verbocount << std::endl; std::cout << "[ NUISANCE ]: Setting ERROR=" << errorcount << std::endl; SETVERBOSITY(verbocount); SETTRACE(trace); // Comparison Setup ======================================== // Proper Setup fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE"); SetupComparisonsFromXML(); SetupRWEngine(); SetupFCN(); return; }; //************************************* void ComparisonRoutines::SetupComparisonsFromXML() { //************************************* LOG(FIT) << "Setting up nuiscomp" << std::endl; // Setup Parameters ------------------------------------------ std::vector parkeys = Config::QueryKeys("parameter"); if (!parkeys.empty()) { LOG(FIT) << "Number of parameters : " << parkeys.size() << std::endl; } for (size_t i = 0; i < parkeys.size(); i++) { nuiskey key = parkeys.at(i); // Check for type,name,nom if (!key.Has("type")) { ERR(FTL) << "No type given for parameter " << i << std::endl; ERR(FTL) << "type='PARAMETER_TYPE'" << std::endl; throw; } else if (!key.Has("name")) { ERR(FTL) << "No name given for parameter " << i << std::endl; ERR(FTL) << "name='SAMPLE_NAME'" << std::endl; throw; } else if (!key.Has("nominal")) { ERR(FTL) << "No nominal given for parameter " << i << std::endl; ERR(FTL) << "nominal='NOMINAL_VALUE'" << std::endl; throw; } // Get Inputs std::string partype = key.GetS("type"); std::string parname = key.GetS("name"); double parnom = key.GetD("nominal"); double parlow = parnom - 1; double parhigh = parnom + 1; double parstep = 1; // override if state not given if (!key.Has("state")){ key.SetS("state","FIX"); } std::string parstate = key.GetS("state"); // Check for incomplete limtis int limdef = ((int)key.Has("low") + (int)key.Has("high") + (int)key.Has("step")); if (limdef > 0 and limdef < 3){ ERR(FTL) << "Incomplete limit set given for parameter : " << parname << std::endl; ERR(FTL) << "Requires: low='LOWER_LIMIT' high='UPPER_LIMIT' step='STEP_SIZE' " << std::endl; throw; } // Extra limits if (key.Has("low")) { parlow = key.GetD("low"); parhigh = key.GetD("high"); parstep = key.GetD("step"); LOG(FIT) << "Read " << partype << " : " << parname << " = " << parnom << " : " << parlow << " < p < " << parhigh << " : " << parstate << std::endl; } else { LOG(FIT) << "Read " << partype << " : " << parname << " = " << parnom << " : " << parstate << std::endl; } // Convert if required if (parstate.find("ABS") != std::string::npos) { parnom = FitBase::RWAbsToSigma( partype, parname, parnom ); parlow = FitBase::RWAbsToSigma( partype, parname, parlow ); parhigh = FitBase::RWAbsToSigma( partype, parname, parhigh ); parstep = FitBase::RWAbsToSigma( partype, parname, parstep ); } else if (parstate.find("FRAC") != std::string::npos) { parnom = FitBase::RWFracToSigma( partype, parname, parnom ); parlow = FitBase::RWFracToSigma( partype, parname, parlow ); parhigh = FitBase::RWFracToSigma( partype, parname, parhigh ); parstep = FitBase::RWFracToSigma( partype, parname, parstep ); } // Push into vectors fParams.push_back(parname); fTypeVals[parname] = FitBase::ConvDialType(partype);; fCurVals[parname] = parnom; fStateVals[parname] = parstate; } // Setup Samples ---------------------------------------------- std::vector samplekeys = Config::QueryKeys("sample"); if (!samplekeys.empty()) { LOG(FIT) << "Number of samples : " << samplekeys.size() << std::endl; } for (size_t i = 0; i < samplekeys.size(); i++) { nuiskey key = samplekeys.at(i); // Get Sample Options std::string samplename = key.GetS("name"); std::string samplefile = key.GetS("input"); std::string sampletype = key.Has("type") ? key.GetS("type") : "DEFAULT"; double samplenorm = key.Has("norm") ? key.GetD("norm") : 1.0; // Print out LOG(FIT) << "Read Sample " << i << ". : " << samplename << " (" << sampletype << ") [Norm=" << samplenorm<<"]"<< std::endl << " -> input='" << samplefile << "'" << std::endl; // If FREE add to parameters otherwise continue if (sampletype.find("FREE") == std::string::npos) { continue; } // Form norm dial from samplename + sampletype + "_norm"; std::string normname = samplename + "_norm"; // Check normname not already present if (fTypeVals.find("normname") != fTypeVals.end()) { continue; } // Add new norm dial to list if its passed above checks fParams.push_back(normname); fTypeVals[normname] = kNORM; fStateVals[normname] = sampletype; fCurVals[normname] = samplenorm; } // Setup Fake Parameters ----------------------------- std::vector fakekeys = Config::QueryKeys("fakeparameter"); if (!fakekeys.empty()) { LOG(FIT) << "Number of fake parameters : " << fakekeys.size() << std::endl; } for (size_t i = 0; i < fakekeys.size(); i++) { nuiskey key = fakekeys.at(i); // Check for type,name,nom if (!key.Has("name")) { ERR(FTL) << "No name given for fakeparameter " << i << std::endl; throw; } else if (!key.Has("nominal")) { ERR(FTL) << "No nominal given for fakeparameter " << i << std::endl; throw; } // Get Inputs std::string parname = key.GetS("name"); double parnom = key.GetD("nominal"); // Push into vectors fFakeVals[parname] = parnom; } } //************************************* void ComparisonRoutines::SetupRWEngine() { //************************************* LOG(FIT) << "Setting up FitWeight Engine" << std::endl; for (UInt_t i = 0; i < fParams.size(); i++) { std::string name = fParams[i]; FitBase::GetRW()->IncludeDial(name, fTypeVals.at(name)); } return; } //************************************* void ComparisonRoutines::SetupFCN() { //************************************* LOG(FIT) << "Building the SampleFCN" << std::endl; if (fSampleFCN) delete fSampleFCN; Config::Get().out = fOutputRootFile; fOutputRootFile->cd(); fSampleFCN = new JointFCN(fOutputRootFile); SetFakeData(); return; } //************************************* void ComparisonRoutines::SetFakeData() { //************************************* if (fFakeDataInput.empty()) return; if (fFakeDataInput.compare("MC") == 0) { LOG(FIT) << "Setting fake data from MC starting prediction." << std::endl; UpdateRWEngine(fFakeVals); FitBase::GetRW()->Reconfigure(); fSampleFCN->ReconfigureAllEvents(); fSampleFCN->SetFakeData("MC"); UpdateRWEngine(fCurVals); LOG(FIT) << "Set all data to fake MC predictions." << std::endl; } else { LOG(FIT) << "Setting fake data from: " << fFakeDataInput << std::endl; fSampleFCN->SetFakeData(fFakeDataInput); } return; } /* Fitting Functions */ //************************************* void ComparisonRoutines::UpdateRWEngine( std::map& updateVals) { //************************************* for (UInt_t i = 0; i < fParams.size(); i++) { std::string name = fParams[i]; if (updateVals.find(name) == updateVals.end()) continue; FitBase::GetRW()->SetDialValue(name, updateVals.at(name)); } FitBase::GetRW()->Reconfigure(); return; } //************************************* void ComparisonRoutines::Run() { //************************************* LOG(FIT) << "Running ComparisonRoutines : " << fStrategy << std::endl; if (FitPar::Config().GetParB("save_nominal")) { SaveNominal(); } // Parse given routines fRoutines = GeneralUtils::ParseToStr(fStrategy, ","); if (fRoutines.empty()) { ERR(FTL) << "Trying to run ComparisonRoutines with no routines given!" << std::endl; throw; } for (UInt_t i = 0; i < fRoutines.size(); i++) { std::string routine = fRoutines.at(i); LOG(FIT) << "Routine: " << routine << std::endl; if (!routine.compare("Compare")) { UpdateRWEngine(fCurVals); GenerateComparison(); PrintState(); SaveCurrentState(); } } return; } //************************************* void ComparisonRoutines::GenerateComparison() { //************************************* LOG(FIT) << "Generating Comparison." << std::endl; // Main Event Loop from event Manager fSampleFCN->ReconfigureAllEvents(); return; } //************************************* void ComparisonRoutines::PrintState() { //************************************* LOG(FIT) << "------------" << std::endl; // Count max size int maxcount = 0; for (UInt_t i = 0; i < fParams.size(); i++) { maxcount = max(int(fParams[i].size()), maxcount); } // Header LOG(FIT) << " # " << left << setw(maxcount) << "Parameter " << " = " << setw(10) << "Value" << " +- " << setw(10) << "Error" << " " << setw(8) << "(Units)" << " " << setw(10) << "Conv. Val" << " +- " << setw(10) << "Conv. Err" << " " << setw(8) << "(Units)" << std::endl; // Parameters for (UInt_t i = 0; i < fParams.size(); i++) { std::string syst = fParams.at(i); std::string typestr = FitBase::ConvDialType(fTypeVals[syst]); std::string curunits = "(sig.)"; double curval = fCurVals[syst]; double curerr = 0.0; if (fStateVals[syst].find("ABS") != std::string::npos) { curval = FitBase::RWSigmaToAbs(typestr, syst, curval); curerr = (FitBase::RWSigmaToAbs(typestr, syst, curerr) - FitBase::RWSigmaToAbs(typestr, syst, 0.0)); curunits = "(Abs.)"; } else if (fStateVals[syst].find("FRAC") != std::string::npos) { curval = FitBase::RWSigmaToFrac(typestr, syst, curval); curerr = (FitBase::RWSigmaToFrac(typestr, syst, curerr) - FitBase::RWSigmaToFrac(typestr, syst, 0.0)); curunits = "(Frac)"; } std::string convunits = "(" + FitBase::GetRWUnits(typestr, syst) + ")"; double convval = FitBase::RWSigmaToAbs(typestr, syst, curval); double converr = (FitBase::RWSigmaToAbs(typestr, syst, curerr) - FitBase::RWSigmaToAbs(typestr, syst, 0.0)); std::ostringstream curparstring; curparstring << " " << setw(3) << left << i << ". " << setw(maxcount) << syst << " = " << setw(10) << curval << " +- " << setw(10) << curerr << " " << setw(8) << curunits << " " << setw(10) << convval << " +- " << setw(10) << converr << " " << setw(8) << convunits; LOG(FIT) << curparstring.str() << std::endl; } LOG(FIT) << "------------" << std::endl; double like = fSampleFCN->GetLikelihood(); LOG(FIT) << std::left << std::setw(46) << "Likelihood for JointFCN: " << like << std::endl; LOG(FIT) << "------------" << std::endl; } /* Write Functions */ //************************************* void ComparisonRoutines::SaveCurrentState(std::string subdir) { //************************************* LOG(FIT) << "Saving current full FCN predictions" << std::endl; // Setup DIRS TDirectory* curdir = gDirectory; if (!subdir.empty()) { TDirectory* newdir = (TDirectory*)gDirectory->mkdir(subdir.c_str()); newdir->cd(); } fSampleFCN->Write(); // Change back to current DIR curdir->cd(); return; } //************************************* void ComparisonRoutines::SaveNominal() { //************************************* fOutputRootFile->cd(); LOG(FIT) << "Saving Nominal Predictions (be cautious with this)" << std::endl; FitBase::GetRW()->Reconfigure(); GenerateComparison(); SaveCurrentState("nominal"); }; - -