diff --git a/app/PrepareGENIE.cxx b/app/PrepareGENIE.cxx
index 60f02ca..374d270 100644
--- a/app/PrepareGENIE.cxx
+++ b/app/PrepareGENIE.cxx
@@ -1,386 +1,554 @@
#include <stdio.h>
#include <stdlib.h>
+#include "FitLogger.h"
+#include "PlotUtils.h"
#include "TFile.h"
#include "TH1D.h"
#include "TTree.h"
-#include "PlotUtils.h"
-#include "FitLogger.h"
#ifdef __GENIE_ENABLED__
#include "Conventions/Units.h"
#include "GHEP/GHepParticle.h"
#include "PDG/PDGUtils.h"
#endif
-bool gFlagMerge = false;
std::string gInputFiles = "";
std::string gOutputFile = "";
-std::string gFluxFile = "";
-std::string gTarget = "";
+std::string gFluxFile = "";
+std::string gTarget = "";
+double MonoEnergy;
+bool IsMonoE = false;
void PrintOptions();
void ParseOptions(int argc, char* argv[]);
-void RunGENIEMerger(std::string inputs, std::string output);
-void RunGENIEPrepare(std::string input, std::string flux, std::string target, std::string output);
+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 (gFlagMerge) RunGENIEMerger(gInputFiles, gOutputFile);
- else RunGENIEPrepare(gInputFiles, gFluxFile, gTarget, gOutputFile);
-
+ ParseOptions(argc, argv);
+ if (IsMonoE) {
+ RunGENIEPrepareMono(gInputFiles, gTarget, gOutputFile);
+ } else {
+ RunGENIEPrepare(gInputFiles, gFluxFile, gTarget, gOutputFile);
+ }
}
-void RunGENIEMerger(std::string inputs, std::string output) {
-
-};
-
-void RunGENIEPrepare(std::string input, std::string flux, std::string target, std::string output) {
-
- LOG(FIT) << "Running GENIE Prepare" << std::endl;
-
- // Setup TTree
- TChain* tn = new TChain("gtree");
- tn->AddFile(input.c_str());
-
- int nevt = tn->GetEntries();
- NtpMCEventRecord * genientpl = NULL;
- bool ismono = false;
- tn->SetBranchAddress("gmcrec", &genientpl);
-
- // Get Flux Hist
- std::vector<std::string> fluxvect = GeneralUtils::ParseToStr(flux, ",");
- TH1D* fluxhist = NULL;
- if (fluxvect.size() > 1) {
- TFile* fluxfile = new TFile(fluxvect[0].c_str(), "READ");
- if (!fluxfile->IsZombie()) {
- fluxhist = (TH1D*) fluxfile->Get(fluxvect[1].c_str());
- fluxhist->SetDirectory(0);
- } else {
-
- // Function with EnuRange
- if (fluxvect.size() == 3) {
-
- ERR(FTL) << "FUNCTION WITH ENU RANGE NOT SUPPORTED SORRY!" << std::endl;
- throw;
-
- // Single Enu
- }
- }
-
- } else if (fluxvect.size() == 1) {
-
- // double E = GeneralUtils::StrToDbl(fluxvect[0]);
- // fluxhist = new TH1D("fluxhist","fluxhist",1, E-0.00001, E+0.00001);
- // fluxhist->SetBinContent(1, 1.0);
- double E = GeneralUtils::StrToDbl(fluxvect[0]);
- LOG(FIT) << "Adding mono-energetic genie flux, E: " << E << std::endl;
-
- fluxhist = new TH1D("fluxhist", "fluxhist", 100, 0, E * 2);
- fluxhist->SetBinContent(fluxhist->FindBin(E), 1);
- ismono = true;
-
- // if (fluxvect[0] == '1.0') {
- // fluxhist = new TH1D("fluxhist", "fluxhist", 40, GeneralUtils::StrToDbl(fluxvect[1]), GeneralUtils::StrToDbl(fluxvect[2]));
- // for (int i = 0; i < fluxhist->GetNbinsX(); i++) {
- // fluxhist->SetBinContent(i + 1, 1.0);
- // }
- // } else {
- // // TF1 f1 = TF1("1.0", GeneralUtils::StrToDbl(fluxvect[1]), GeneralUtils::StrToDbl(fluxvect[2]), 1);
- // // std::cout << "Created TF1 with '" << fluxvect[0].c_str()<<"' " << GeneralUtils::StrToDbl(fluxvect[1]) << " " << GeneralUtils::StrToDbl(fluxvect[2]) << std::endl;
- // fluxhist = new TH1D("fluxhist", "fluxhist", 40, GeneralUtils::StrToDbl(fluxvect[1]), GeneralUtils::StrToDbl(fluxvect[2]));
- // for (int i = 0; i < fluxhist->GetNbinsX(); i++) {
- // fluxhist->SetBinContent(i + 1, 1.0); //f1.Eval(fluxhist->GetXaxis()->GetBinCenter(i+1)));
- // // std::cout << "Filling Flux Hist " << f1.Eval(fluxhist->GetXaxis()->GetBinCenter(i+1)) << std::endl;
- // }
- // sleep(10);
-
- } else {
- LOG(FTL) << "NO FLUX SPECIFIED" << std::endl;
- throw;
- }
-
- // 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();
- }
-
- // 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();
-
- if (i % (nevt / 20) == 0) {
- LOG(FIT) << "Processed " << i << "/" << nevt << " GENIE events." << std::endl;
- }
- }
- LOG(FIT) << "Processed all events" << std::endl;
-
- // If Mono// You have filled the total xsec for each event
- // Then divided by the total number of events, and then multiplied by the flux integral * AvgXSec )
- // eventlist[pdg] -> Scale(1.0 / zeroevents->Integral());
- // eventlist[pdg] -> Scale( fluxlist[0]->Integral()*AvgXSec );
-
- // Once event loop is done we can start saving stuff into the file
- // bool savesplines = FitPar::Config().GetParB("save_genie_splines"); // Currently not implemented
-
- TFile* outputfile = new TFile(input.c_str(), "UPDATE");
- outputfile->cd();
-
- 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]->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++) {
- LOG(FIT) << "Getting target " << i << std::endl;
- std::string targ = targetids[i];
- targetsplines[targ] = (TH1D*) xsechist->Clone();
- 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;
- // modeavg[mode]->Write( (mode + "_cont_" + targ).c_str() , TObject::kOverwrite);
- 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;
- std::vector<std::string> targprs = GeneralUtils::ParseToStr(target, ",");
- TH1D* totalxsec = (TH1D*) xsechist->Clone();
-
- for (uint i = 0; i < targprs.size(); i++) {
- std::string targpdg = targprs[i];
-
- for (std::map<std::string, TH1D*>::iterator iter = targetsplines.begin();
- iter != targetsplines.end(); iter++) {
- std::string targstr = iter->first;
- TH1D* xsec = iter->second;
-
- if (targstr.find(targpdg) != std::string::npos) {
- 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);
-
- }
- }
- }
-
- if (!ismono){
- LOG(FIT) << "Total XSec Integral = " << totalxsec->Integral("width") << std::endl;
-
- outputfile->cd();
- totalxsec->Write("nuisance_xsec", TObject::kOverwrite);
- eventhist = (TH1D*)fluxhist->Clone();
- eventhist->Multiply(totalxsec);
-
- LOG(FIT) << "Dividing by Total Nucl = " << totalnucl << std::endl;
- eventhist->Scale(1.0 / double(totalnucl) );
- } else {
-
- outputfile->cd();
- totalxsec->Write("nuisance_Xsec", TObject::kOverwrite);
- eventhist = (TH1D*) fluxhist->Clone();
- eventhist->Scale(totalxsec->Integral());
-
- }
-
- 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;
- std::cout << "XSec Hist Integral = " << xsechist->Integral("width") << std::endl;
+void RunGENIEPrepareMono(std::string input, std::string target,
+ std::string output) {
+ // Setup TTree
+ TChain* tn = new TChain("gtree");
+ tn->AddFile(input.c_str());
+
+ int nevt = tn->GetEntries();
+ NtpMCEventRecord* genientpl = NULL;
+ tn->SetBranchAddress("gmcrec", &genientpl);
+
+ TH1D *fluxhist = new TH1D("flux", "flux", 1000, 0, 10);
+ 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();
+ }
+
+ // 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();
+
+ if (i % (nevt / 20) == 0) {
+ LOG(FIT) << "Processed " << i << "/" << nevt << " GENIE events."
+ << std::endl;
+ }
+ }
+ LOG(FIT) << "Processed all events" << std::endl;
+
+ TFile* outputfile = new TFile(input.c_str(), "UPDATE");
+ outputfile->cd();
+
+ 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]->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++) {
+ LOG(FIT) << "Getting target " << i << std::endl;
+ std::string targ = targetids[i];
+ targetsplines[targ] = (TH1D*)xsechist->Clone();
+ 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;
+ std::vector<std::string> targprs = GeneralUtils::ParseToStr(target, ",");
+ TH1D* totalxsec = (TH1D*)xsechist->Clone();
+
+ for (uint i = 0; i < targprs.size(); i++) {
+ std::string targpdg = targprs[i];
+
+ for (std::map<std::string, TH1D*>::iterator iter = targetsplines.begin();
+ iter != targetsplines.end(); iter++) {
+ std::string targstr = iter->first;
+ TH1D* xsec = iter->second;
+
+ if (targstr.find(targpdg) != std::string::npos) {
+ 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);
+ }
+ }
+ }
+
+ outputfile->cd();
+ totalxsec->Write("nuisance_Xsec", TObject::kOverwrite);
+ eventhist = (TH1D*)totalxsec->Clone();
+ eventhist->Multiply(fluxHist);
+
+ 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;
+ std::cout << "XSec Hist Integral = " << xsechist->Integral("width")
+ << std::endl;
+
+ return;
+}
- return;
+void RunGENIEPrepare(std::string input, std::string flux, std::string target,
+ std::string output) {
+ LOG(FIT) << "Running GENIE Prepare" << std::endl;
+
+ // Get Flux Hist
+ std::vector<std::string> fluxvect = GeneralUtils::ParseToStr(flux, ",");
+ TH1D* fluxhist = NULL;
+ if (fluxvect.size() > 1) {
+ TFile* fluxfile = new TFile(fluxvect[0].c_str(), "READ");
+ if (!fluxfile->IsZombie()) {
+ fluxhist = (TH1D*)fluxfile->Get(fluxvect[1].c_str());
+ fluxhist->SetDirectory(0);
+ } else {
+ // Function with EnuRange
+ if (fluxvect.size() == 3) {
+ ERR(FTL) << "FUNCTION WITH ENU RANGE NOT SUPPORTED SORRY!" << std::endl;
+ throw;
+ }
+ }
+
+ } else if (fluxvect.size() == 1) {
+ MonoEnergy = GeneralUtils::StrToDbl(fluxvect[0]);
+ RunGENIEPrepareMono(input, target, output);
+ return;
+ } else {
+ LOG(FTL) << "NO FLUX SPECIFIED" << std::endl;
+ throw;
+ }
+
+ // Setup TTree
+ TChain* tn = new TChain("gtree");
+ tn->AddFile(input.c_str());
+
+ int nevt = tn->GetEntries();
+ NtpMCEventRecord* genientpl = NULL;
+ tn->SetBranchAddress("gmcrec", &genientpl);
+
+ // 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();
+ }
+
+ // 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();
+
+ if (i % (nevt / 20) == 0) {
+ LOG(FIT) << "Processed " << i << "/" << nevt << " GENIE events."
+ << std::endl;
+ }
+ }
+ LOG(FIT) << "Processed all events" << std::endl;
+
+ // Once event loop is done we can start saving stuff into the file
+
+ TFile* outputfile = new TFile(input.c_str(), "UPDATE");
+ outputfile->cd();
+
+ 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]->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++) {
+ LOG(FIT) << "Getting target " << i << std::endl;
+ std::string targ = targetids[i];
+ targetsplines[targ] = (TH1D*)xsechist->Clone();
+ 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;
+ // modeavg[mode]->Write( (mode + "_cont_" + targ).c_str() ,
+ // TObject::kOverwrite);
+ 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;
+ std::vector<std::string> targprs = GeneralUtils::ParseToStr(target, ",");
+ TH1D* totalxsec = (TH1D*)xsechist->Clone();
+
+ for (uint i = 0; i < targprs.size(); i++) {
+ std::string targpdg = targprs[i];
+
+ for (std::map<std::string, TH1D*>::iterator iter = targetsplines.begin();
+ iter != targetsplines.end(); iter++) {
+ std::string targstr = iter->first;
+ TH1D* xsec = iter->second;
+
+ if (targstr.find(targpdg) != std::string::npos) {
+ 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);
+ }
+ }
+ }
+ LOG(FIT) << "Total XSec Integral = " << totalxsec->Integral("width")
+ << std::endl;
+
+ outputfile->cd();
+ totalxsec->Write("nuisance_xsec", TObject::kOverwrite);
+ eventhist = (TH1D*)fluxhist->Clone();
+ eventhist->Multiply(totalxsec);
+
+ 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;
+ std::cout << "XSec Hist Integral = " << xsechist->Integral("width")
+ << std::endl;
+
+ return;
};
-
-
void PrintOptions() {
-
- std::cout << "PrepareGENIEEvents NUISANCE app. " << std::endl
- << "Takes GHep Outputs and prepares events for NUISANCE." << std::endl << std::endl
- << "PrepareGENIEEvents [-h,-help,--h,--help] [-i inputfile1.root,inputfile2.root,inputfile3.root,...] "
- << "[-f flux_root_file.root,flux_hist_name] [-t target1[frac1],target2[frac2],...]"
- << 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 << " [ -t target ] : Target that GHepRecords were generated with. Comma seperated list. E.g. for CH2 target=1000060120,1000010010,1000010010" << std::endl;
-
- /*
- std::cout << "Merger Mode [activate with -m] : Takes the list of input files assuming 'Prepare Mode' has already been ran on them and merges them "
- << "into a single file with associated Friend Tree to help with conserving ratios of events (e.g. adding nue and nueb beams together into a single file" << std::endl;
- std::cout << "Following optoins are required for Merger Mode:" << std::endl;
- std::cout << " [ -i inputfile1.root,inputfile2.root ] : Comma Seperated list of files to be merged. " << std::endl;
- std::cout << " [ -o outputfile.root ] : Output file for merger." << std::endl;
- */
-
+ std::cout << "PrepareGENIEEvents NUISANCE app. " << std::endl
+ << "Takes GHep Outputs and prepares events for NUISANCE."
+ << std::endl
+ << std::endl
+ << "PrepareGENIEEvents [-h,-help,--h,--help] [-i "
+ "inputfile1.root,inputfile2.root,inputfile3.root,...] "
+ << "[-f flux_root_file.root,flux_hist_name] [-t "
+ "target1[frac1],target2[frac2],...]"
+ << 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 << " [ -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."
+ << 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 (!std::strcmp(argv[i], "-m")) { gFlagMerge = 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 {
- ERR(FTL) << "ERROR: unknown command line option given! - '"
- << argv[i] << " " << argv[i + 1] << "'" << std::endl;
- PrintOptions();
- break;
- }
- }
- }
-
- /*
- if (gOutputFile == "" && !flagopt){
- ERR(FTL) << "No output file specificed!" << std::endl;
- flagopt = true;
- }
- */
-
- if (gInputFiles == "" && !flagopt) {
- ERR(FTL) << "No input file(s) specified!" << std::endl;
- flagopt = true;
- }
-
- if (!gFlagMerge && gFluxFile == "" && !flagopt) {
- ERR(FTL) << "No flux input specified for Prepare Mode" << std::endl;
- flagopt = true;
- }
-
- if (!gFlagMerge && gTarget == "" && !flagopt) {
- ERR(FTL) << "No target specified for Prepare Mode" << std::endl;
- flagopt = true;
- }
-
- if (argc < 1 || flagopt) {
- PrintOptions();
- exit(-1);
- }
-
- return;
+ 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], "-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) {
+ 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 (argc < 1 || flagopt) {
+ PrintOptions();
+ exit(-1);
+ }
+
+ return;
}
diff --git a/app/PrepareNEUT.cxx b/app/PrepareNEUT.cxx
index dd1f8bf..e465855 100644
--- a/app/PrepareNEUT.cxx
+++ b/app/PrepareNEUT.cxx
@@ -1,418 +1,416 @@
#include <stdio.h>
#include <stdlib.h>
#include "FitLogger.h"
#include "PlotUtils.h"
#include "StatUtils.h"
#include "TFile.h"
#include "TH1D.h"
#include "TTree.h"
// If you don't have NEUT enabled, you shouldn't compile this...
#include "neutpart.h"
#include "neutvect.h"
std::string fInputFiles = "";
std::string fOutputFile = "";
std::string fFluxFile = "";
bool fFluxInGeV = false;
bool fIsMonoEFlux = false;
double fMonoEEnergy = 0xdeadbeef;
void PrintOptions();
void ParseOptions(int argc, char* argv[]);
void AddMonoRateHistogram(std::string inputList, double MonoE,
std::string output);
void CreateRateHistogram(std::string inputList, std::string flux,
std::string output);
//*******************************
int main(int argc, char* argv[]) {
//*******************************
LOG_VERB(FitPar::Config().GetParI("VERBOSITY"));
ERR_VERB(FitPar::Config().GetParI("ERROR"));
ParseOptions(argc, argv);
LOG(FIT) << "Running PrepareNEUT" << std::endl;
if (fIsMonoEFlux) {
AddMonoRateHistogram(fInputFiles, fMonoEEnergy, fOutputFile);
} else {
CreateRateHistogram(fInputFiles, fFluxFile, fOutputFile);
}
};
void AddMonoRateHistogram(std::string inputList, double MonoE,
std::string output) {
// Need to allow for more than one file... will do soon
TChain* tn = new TChain("neuttree");
std::vector<std::string> inputs = GeneralUtils::ParseToStr(inputList, ",");
for (std::vector<std::string>::iterator it = inputs.begin();
it != inputs.end(); ++it) {
LOG(FIT) << "Adding " << *it << " to the output" << std::endl;
tn->AddFile((*it).c_str());
}
if (inputs.size() > 1 && output.empty()) {
ERR(FTL) << "You must provide a new output file name if you want to have "
"more than 1 input file!"
<< std::endl;
throw;
}
int nevts = tn->GetEntries();
if (!nevts) {
ERR(FTL) << "Either the input file is not from NEUT, or it's empty..."
<< std::endl;
throw;
}
NeutVect* fNeutVect = NULL;
tn->SetBranchAddress("vectorbranch", &fNeutVect);
TH1D* fluxHist = new TH1D("flux", "flux", 1000, 0, fFluxInGeV ? 10 : 10000);
fluxHist->Fill(MonoE);
- fluxHist->Scale(1,"width");
+ fluxHist->Scale(1, "width");
// Make Event Hist
TH1D* xsecHist = (TH1D*)fluxHist->Clone();
xsecHist->Reset();
// Make a total cross section hist for shits and giggles
TH1D* entryHist = (TH1D*)xsecHist->Clone();
for (int i = 0; i < nevts; ++i) {
tn->GetEntry(i);
NeutPart* part = fNeutVect->PartInfo(0);
double E = part->fP.E();
double xsec = fNeutVect->Totcrs;
// Unit conversion
if (fFluxInGeV) E *= 1E-3;
xsecHist->Fill(E, xsec);
entryHist->Fill(E);
if (i % (nevts / 20) == 0) {
LOG(FIT) << "Processed " << i << "/" << nevts << " NEUT events."
<< std::endl;
}
}
LOG(FIT) << "Processed all events" << std::endl;
xsecHist->Divide(entryHist);
// This will be the evtrt histogram
TH1D* evtHist = (TH1D*)xsecHist->Clone();
evtHist->Multiply(fluxHist);
// Check whether the overflow is empty. If not, advise that either the wrong
// flux histogram or units were used...
// If the events were generated with a limited range of the flux histogram,
// this may be benign
if (evtHist->Integral(0, -1) != evtHist->Integral() ||
evtHist->Integral(0, -1) == 0) {
ERR(WRN) << "The input file and flux histogram provided do not match... "
<< std::endl;
ERR(WRN) << "Are the units correct? Did you provide the correct flux file?"
<< std::endl;
ERR(WRN) << "Use output with caution..." << std::endl;
}
// Pick where the output should go
TFile* outFile = NULL;
if (!output.empty()) {
LOG(FIT) << "Saving histograms in " << output << std::endl;
outFile = new TFile(output.c_str(), "RECREATE");
} else {
LOG(FIT) << "Saving histograms in " << inputs[0] << std::endl;
outFile = new TFile(inputs[0].c_str(), "UPDATE");
}
outFile->cd();
std::string xsec_name = "xsec_PrepareNeut";
std::string flux_name = "flux_PrepareNeut";
std::string rate_name = "evtrt_PrepareNeut";
if (output.empty()) {
// Check whether we should overwrite existing histograms
std::string input_xsec = PlotUtils::GetObjectWithName(outFile, "xsec");
std::string input_flux = PlotUtils::GetObjectWithName(outFile, "flux");
std::string input_rate = PlotUtils::GetObjectWithName(outFile, "evtrt");
if (!input_xsec.empty()) {
LOG(FIT) << "Updating histogram: " << input_xsec << std::endl;
xsec_name = input_xsec;
}
if (!input_flux.empty()) {
LOG(FIT) << "Updating histogram: " << input_flux << std::endl;
flux_name = input_flux;
}
if (!input_rate.empty()) {
LOG(FIT) << "Updating histogram: " << input_rate << std::endl;
rate_name = input_rate;
}
} else {
LOG(FIT) << "Cloning neuttree into output file." << std::endl;
StopTalking();
TTree* newtree = (TTree*)tn->CloneTree(-1, "fast");
StartTalking();
newtree->Write();
}
xsecHist->Write(xsec_name.c_str(), TObject::kOverwrite);
fluxHist->Write(flux_name.c_str(), TObject::kOverwrite);
evtHist->Write(rate_name.c_str(), TObject::kOverwrite);
outFile->Close();
}
//*******************************
void CreateRateHistogram(std::string inputList, std::string flux,
std::string output) {
//*******************************
// Need to allow for more than one file... will do soon
TChain* tn = new TChain("neuttree");
std::vector<std::string> inputs = GeneralUtils::ParseToStr(inputList, ",");
for (std::vector<std::string>::iterator it = inputs.begin();
it != inputs.end(); ++it) {
LOG(FIT) << "Adding " << *it << " to the output" << std::endl;
tn->AddFile((*it).c_str());
}
if (inputs.size() > 1 && output.empty()) {
ERR(FTL) << "You must provide a new output file name if you want to have "
"more than 1 input file!"
<< std::endl;
throw;
}
int nevts = tn->GetEntries();
if (!nevts) {
ERR(FTL) << "Either the input file is not from NEUT, or it's empty..."
<< std::endl;
throw;
}
NeutVect* fNeutVect = NULL;
tn->SetBranchAddress("vectorbranch", &fNeutVect);
// Get Flux Hist
std::vector<std::string> fluxvect = GeneralUtils::ParseToStr(flux, ",");
TH1D* fluxHist = NULL;
if (fluxvect.size() > 1) {
TFile* fluxfile = new TFile(fluxvect[0].c_str(), "READ");
fluxHist = (TH1D*)fluxfile->Get(fluxvect[1].c_str());
fluxHist->SetDirectory(0);
} else {
ERR(FTL) << "NO FLUX SPECIFIED" << std::endl;
throw;
}
// Decide what type of flux was given
if (fFluxInGeV)
LOG(FIT) << "Assuming flux histogram is in GeV" << std::endl;
else
LOG(FIT) << "Assuming flux histogram is in MeV" << std::endl;
// Make Event Hist
TH1D* xsecHist = (TH1D*)fluxHist->Clone();
xsecHist->Reset();
// Make a total cross section hist for shits and giggles
TH1D* entryHist = (TH1D*)xsecHist->Clone();
for (int i = 0; i < nevts; ++i) {
tn->GetEntry(i);
NeutPart* part = fNeutVect->PartInfo(0);
double E = part->fP.E();
double xsec = fNeutVect->Totcrs;
// Unit conversion
if (fFluxInGeV) E *= 1E-3;
xsecHist->Fill(E, xsec);
entryHist->Fill(E);
if (i % (nevts / 20) == 0) {
LOG(FIT) << "Processed " << i << "/" << nevts << " NEUT events."
<< std::endl;
}
}
LOG(FIT) << "Processed all events" << std::endl;
xsecHist->Divide(entryHist);
// This will be the evtrt histogram
TH1D* evtHist = NULL;
// If the integral of xsecHist is 0 the input file used a really old version
// of NEUT without Totcrs
if (!xsecHist->Integral(0, -1)) {
ERR(WRN) << "Old NEUT input file: events will not be correctly normalized"
<< std::endl;
evtHist = (TH1D*)entryHist->Clone();
if (evtHist->Integral() != 0)
evtHist->Scale(fluxHist->Integral() / float(evtHist->Integral()));
} else {
evtHist = (TH1D*)xsecHist->Clone();
evtHist->Multiply(fluxHist);
}
// Check whether the overflow is empty. If not, advise that either the wrong
// flux histogram or units were used...
// If the events were generated with a limited range of the flux histogram,
// this may be benign
if (evtHist->Integral(0, -1) != evtHist->Integral() ||
evtHist->Integral(0, -1) == 0) {
ERR(WRN) << "The input file and flux histogram provided do not match... "
<< std::endl;
ERR(WRN) << "Are the units correct? Did you provide the correct flux file?"
<< std::endl;
ERR(WRN) << "Use output with caution..." << std::endl;
}
// Pick where the output should go
TFile* outFile = NULL;
if (!output.empty()) {
LOG(FIT) << "Saving histograms in " << output << std::endl;
outFile = new TFile(output.c_str(), "RECREATE");
} else {
LOG(FIT) << "Saving histograms in " << inputs[0] << std::endl;
outFile = new TFile(inputs[0].c_str(), "UPDATE");
}
outFile->cd();
std::string xsec_name = "xsec_PrepareNeut";
std::string flux_name = "flux_PrepareNeut";
std::string rate_name = "evtrt_PrepareNeut";
if (output.empty()) {
// Check whether we should overwrite existing histograms
std::string input_xsec = PlotUtils::GetObjectWithName(outFile, "xsec");
std::string input_flux = PlotUtils::GetObjectWithName(outFile, "flux");
std::string input_rate = PlotUtils::GetObjectWithName(outFile, "evtrt");
if (!input_xsec.empty()) {
LOG(FIT) << "Updating histogram: " << input_xsec << std::endl;
xsec_name = input_xsec;
}
if (!input_flux.empty()) {
LOG(FIT) << "Updating histogram: " << input_flux << std::endl;
flux_name = input_flux;
}
if (!input_rate.empty()) {
LOG(FIT) << "Updating histogram: " << input_rate << std::endl;
rate_name = input_rate;
}
} else {
LOG(FIT) << "Cloning neuttree into output file." << std::endl;
StopTalking();
TTree* newtree = (TTree*)tn->CloneTree(-1, "fast");
StartTalking();
newtree->Write();
}
xsecHist->Write(xsec_name.c_str(), TObject::kOverwrite);
fluxHist->Write(flux_name.c_str(), TObject::kOverwrite);
evtHist->Write(rate_name.c_str(), TObject::kOverwrite);
outFile->Close();
return;
}
void PrintOptions() {
std::cout << "PrepareNEUT NUISANCE app. " << std::endl
<< "Produces or recalculates evtrt and flux histograms necessary "
"for NUISANCE normalization."
<< std::endl;
std::cout << "PrepareNEUT: " << std::endl;
std::cout << " [-h,-help,--h,--help]" << std::endl;
std::cout << " -i inputfile1.root,inputfile2.root,inputfile3.root,..."
<< std::endl;
std::cout << " Takes any number of files, but assumes all are "
"produced with a single flux"
<< std::endl;
std::cout << " -f flux_root_file.root,flux_hist_name" << std::endl;
std::cout << " Path to root file containing the flux histogram used "
"when generating the NEUT files"
<< std::endl;
std::cout << " [-o outputfile.root] " << std::endl;
std::cout
<< " If an output file is not given, the input file will be used"
<< std::endl;
std::cout << " If more than one input file is given, an output file "
"must be given"
<< std::endl;
std::cout << " [-G]" << std::endl;
std::cout << " Flux is assumed to be in MeV. This switch indicates "
"the input flux is in GeV"
<< std::endl;
std::cout << " [-m E_nu]" << std::endl;
std::cout << " Used to add dummy flux and evt rate histograms to "
"mono-energetic vectors. Adheres to the -G flag."
<< 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;
- } else if (!std::strcmp(argv[i], "-G")) {
- fFluxInGeV = true;
- continue;
- } else if (!std::strcmp(argv[i], "-m")) {
- fIsMonoEFlux = true;
- fMonoEEnergy = GeneralUtils::StrToDbl(argv[i + 1]);
- ++i;
- continue;
}
if (i + 1 != argc) {
// Cardfile
if (!std::strcmp(argv[i], "-h")) {
flagopt = true;
break;
} else if (!std::strcmp(argv[i], "-i")) {
fInputFiles = argv[i + 1];
++i;
} else if (!std::strcmp(argv[i], "-o")) {
fOutputFile = argv[i + 1];
++i;
} else if (!std::strcmp(argv[i], "-f")) {
fFluxFile = argv[i + 1];
++i;
+ } else if (!std::strcmp(argv[i], "-G")) {
+ fFluxInGeV = true;
+ } else if (!std::strcmp(argv[i], "-m")) {
+ fIsMonoEFlux = true;
+ fMonoEEnergy = GeneralUtils::StrToDbl(argv[i + 1]);
+ ++i;
} else {
ERR(FTL) << "ERROR: unknown command line option given! - '" << argv[i]
<< " " << argv[i + 1] << "'" << std::endl;
PrintOptions();
break;
}
}
}
if (fInputFiles == "" && !flagopt) {
ERR(FTL) << "No input file(s) specified!" << std::endl;
flagopt = true;
}
if (fFluxFile == "" && (!flagopt) && (!fIsMonoEFlux)) {
ERR(FTL) << "No flux input specified!" << std::endl;
flagopt = true;
}
if (argc < 1 || flagopt) {
PrintOptions();
exit(-1);
}
return;
}