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 #include +#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 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 modexsec; - std::map modecount; - std::vector genieids; - std::vector targetids; - std::vector interids; - - - // Loop over all events - for (int i = 0; i < nevt; i++) { - tn->GetEntry(i); - - StopTalking(); - EventRecord& event = *(genientpl->event); - GHepParticle* neu = event.Probe(); - StartTalking(); - - // Get XSec From Spline - GHepRecord genie_record = static_cast(event); - double xsec = (genie_record.XSec() / (1E-38 * genie::units::cm2)); - - // Parse Interaction String - std::string mode = genie_record.Summary()->AsString(); - std::vector modevec = GeneralUtils::ParseToStr(mode, ";"); - std::string targ = ( modevec[0] + ";" + modevec[1] ); - std::string inter = mode; - - // Fill lists of Unique IDS - if (std::find(targetids.begin(), targetids.end(), targ) - == targetids.end()) { - targetids.push_back(targ); - } - - if (std::find(interids.begin(), interids.end(), inter) - == interids.end()) { - interids.push_back(inter); - } - - // Create entries Mode Maps - if (modexsec.find(mode) == modexsec.end()) { - genieids.push_back(mode); - - modexsec[mode] = (TH1D*)xsechist->Clone(); - modecount[mode] = (TH1D*)xsechist->Clone(); - } - - // 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 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 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 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::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 modexsec; + std::map modecount; + std::vector genieids; + std::vector targetids; + std::vector interids; + + // Loop over all events + for (int i = 0; i < nevt; i++) { + tn->GetEntry(i); + + StopTalking(); + EventRecord& event = *(genientpl->event); + GHepParticle* neu = event.Probe(); + StartTalking(); + + // Get XSec From Spline + GHepRecord genie_record = static_cast(event); + double xsec = (genie_record.XSec() / (1E-38 * genie::units::cm2)); + + // Parse Interaction String + std::string mode = genie_record.Summary()->AsString(); + std::vector modevec = GeneralUtils::ParseToStr(mode, ";"); + std::string targ = (modevec[0] + ";" + modevec[1]); + std::string inter = mode; + + // Fill lists of Unique IDS + if (std::find(targetids.begin(), targetids.end(), targ) == + targetids.end()) { + targetids.push_back(targ); + } + + if (std::find(interids.begin(), interids.end(), inter) == interids.end()) { + interids.push_back(inter); + } + + // Create entries Mode Maps + if (modexsec.find(mode) == modexsec.end()) { + genieids.push_back(mode); + + modexsec[mode] = (TH1D*)xsechist->Clone(); + modecount[mode] = (TH1D*)xsechist->Clone(); + } + + // 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 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 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 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::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 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 modexsec; + std::map modecount; + std::vector genieids; + std::vector targetids; + std::vector interids; + + // Loop over all events + for (int i = 0; i < nevt; i++) { + tn->GetEntry(i); + + StopTalking(); + EventRecord& event = *(genientpl->event); + GHepParticle* neu = event.Probe(); + StartTalking(); + + // Get XSec From Spline + GHepRecord genie_record = static_cast(event); + double xsec = (genie_record.XSec() / (1E-38 * genie::units::cm2)); + + // Parse Interaction String + std::string mode = genie_record.Summary()->AsString(); + std::vector modevec = GeneralUtils::ParseToStr(mode, ";"); + std::string targ = (modevec[0] + ";" + modevec[1]); + std::string inter = mode; + + // Fill lists of Unique IDS + if (std::find(targetids.begin(), targetids.end(), targ) == + targetids.end()) { + targetids.push_back(targ); + } + + if (std::find(interids.begin(), interids.end(), inter) == interids.end()) { + interids.push_back(inter); + } + + // Create entries Mode Maps + if (modexsec.find(mode) == modexsec.end()) { + genieids.push_back(mode); + + modexsec[mode] = (TH1D*)xsechist->Clone(); + modecount[mode] = (TH1D*)xsechist->Clone(); + } + + // 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 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 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 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::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 #include #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 inputs = GeneralUtils::ParseToStr(inputList, ","); for (std::vector::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 inputs = GeneralUtils::ParseToStr(inputList, ","); for (std::vector::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 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; }