diff --git a/app/PrepareNEUT.cxx b/app/PrepareNEUT.cxx index 3b37e40..dd1f8bf 100644 --- a/app/PrepareNEUT.cxx +++ b/app/PrepareNEUT.cxx @@ -1,238 +1,418 @@ #include #include +#include "FitLogger.h" +#include "PlotUtils.h" +#include "StatUtils.h" #include "TFile.h" #include "TH1D.h" #include "TTree.h" -#include "PlotUtils.h" -#include "StatUtils.h" -#include "FitLogger.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 = ""; +std::string fFluxFile = ""; bool fFluxInGeV = false; +bool fIsMonoEFlux = false; +double fMonoEEnergy = 0xdeadbeef; void PrintOptions(); void ParseOptions(int argc, char* argv[]); -void CreateRateHistogram(std::string inputList, std::string flux, std::string output); +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[]){ -//******************************* +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; - CreateRateHistogram(fInputFiles, fFluxFile, fOutputFile); + 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"); + // 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){ -//******************************* +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,","); + std::vector inputs = GeneralUtils::ParseToStr(inputList, ","); for (std::vector::iterator it = inputs.begin(); - it != inputs.end(); ++it){ + 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; + 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..." <SetBranchAddress("vectorbranch", &fNeutVect); // Get Flux Hist - std::vector fluxvect = GeneralUtils::ParseToStr(flux,","); + 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()); + 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; + 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(); + TH1D* entryHist = (TH1D*)xsecHist->Clone(); - for (int i = 0; i < nevts; ++i){ + for (int i = 0; i < nevts; ++i) { tn->GetEntry(i); - NeutPart* part = fNeutVect->PartInfo(0); + 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); + 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; + if (i % (nevts / 20) == 0) { + LOG(FIT) << "Processed " << i << "/" << nevts << " NEUT events." + << std::endl; } } LOG(FIT) << "Processed all events" << std::endl; - - xsecHist ->Divide(entryHist); + + 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; + + // 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())); + evtHist->Scale(fluxHist->Integral() / float(evtHist->Integral())); } else { evtHist = (TH1D*)xsecHist->Clone(); - evtHist ->Multiply(fluxHist); + 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; + // 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()){ + 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()){ + 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(); + 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; +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 << " -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 << " 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 + << " 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::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[]){ +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; } - if (i+1 != argc){ - + 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 { - ERR(FTL) << "ERROR: unknown command line option given! - '" - <. *******************************************************************************/ #include "InputHandler.h" #include "InputUtils.h" InputHandlerBase::InputHandlerBase() { fName = ""; fFluxHist = NULL; fEventHist = NULL; fNEvents = 0; fNUISANCEEvent = NULL; fBaseEvent = NULL; kRemoveUndefParticles = FitPar::Config().GetParB("RemoveUndefParticles"); kRemoveFSIParticles = FitPar::Config().GetParB("RemoveFSIParticles"); kRemoveNuclearParticles = FitPar::Config().GetParB("RemoveNuclearParticles"); fMaxEvents = FitPar::Config().GetParI("MAXEVENTS"); fTTreePerformance = NULL; }; InputHandlerBase::~InputHandlerBase() { if (fFluxHist) delete fFluxHist; if (fEventHist) delete fEventHist; // if (fXSecHist) delete fXSecHist; // if (fNUISANCEEvent) delete fNUISANCEEvent; jointfluxinputs.clear(); jointeventinputs.clear(); jointindexlow.clear(); jointindexhigh.clear(); jointindexallowed.clear(); jointindexscale.clear(); // if (fTTreePerformance) { // fTTreePerformance->SaveAs(("ttreeperfstats_" + fName + // ".root").c_str()); // } } void InputHandlerBase::Print(){}; TH1D* InputHandlerBase::GetXSecHistogram(void) { fXSecHist = (TH1D*)fFluxHist->Clone(); fXSecHist->Divide(fEventHist); return fXSecHist; }; double InputHandlerBase::PredictedEventRate(double low, double high, std::string intOpt) { int minBin = fFluxHist->GetXaxis()->FindBin(low); int maxBin = fFluxHist->GetXaxis()->FindBin(high); return fEventHist->Integral(minBin, maxBin + 1, intOpt.c_str()); }; double InputHandlerBase::TotalIntegratedFlux(double low, double high, std::string intOpt) { Int_t minBin = fFluxHist->GetXaxis()->FindFixBin(low); Int_t maxBin = fFluxHist->GetXaxis()->FindFixBin(high); if ((fFluxHist->IsBinOverflow(minBin) && (low != -9999.9))) { minBin = 1; } if ((fFluxHist->IsBinOverflow(maxBin) && (high != -9999.9))) { maxBin = fFluxHist->GetXaxis()->GetNbins() + 1; } // If we are within a single bin if (minBin == maxBin) { // Get the contained fraction of the single bin's width return ((high - low) / fFluxHist->GetXaxis()->GetBinWidth(minBin)) * fFluxHist->Integral(minBin, minBin, intOpt.c_str()); } double lowBinUpEdge = fFluxHist->GetXaxis()->GetBinUpEdge(minBin); double highBinLowEdge = fFluxHist->GetXaxis()->GetBinLowEdge(maxBin); double lowBinfracIntegral = ((lowBinUpEdge - low) / fFluxHist->GetXaxis()->GetBinWidth(minBin)) * fFluxHist->Integral(minBin, minBin, intOpt.c_str()); double highBinfracIntegral = ((high - highBinLowEdge) / fFluxHist->GetXaxis()->GetBinWidth(maxBin)) * fFluxHist->Integral(maxBin, maxBin, intOpt.c_str()); // If they are neighbouring bins if ((minBin + 1) == maxBin) { std::cout << "Get lowfrac + highfrac" << std::endl; // Get the contained fraction of the two bin's width return lowBinfracIntegral + highBinfracIntegral; } + double ContainedIntegral = + fFluxHist->Integral(minBin + 1, maxBin - 1, intOpt.c_str()); // If there are filled bins between them - return lowBinfracIntegral + highBinfracIntegral + - fFluxHist->Integral(minBin + 1, maxBin - 1, intOpt.c_str()); + return lowBinfracIntegral + highBinfracIntegral + ContainedIntegral; // return fFluxHist->Integral(minBin + 1, maxBin - 1, intOpt.c_str()); } std::vector InputHandlerBase::GetFluxList(void) { return std::vector(1, fFluxHist); }; std::vector InputHandlerBase::GetEventList(void) { return std::vector(1, fEventHist); }; std::vector InputHandlerBase::GetXSecList(void) { return std::vector(1, GetXSecHistogram()); }; FitEvent* InputHandlerBase::FirstNuisanceEvent() { fCurrentIndex = 0; return GetNuisanceEvent(fCurrentIndex); }; FitEvent* InputHandlerBase::NextNuisanceEvent() { fCurrentIndex++; if ((fMaxEvents != -1) && (fCurrentIndex > fMaxEvents)) { return NULL; } return GetNuisanceEvent(fCurrentIndex); }; BaseFitEvt* InputHandlerBase::FirstBaseEvent() { fCurrentIndex = 0; return GetBaseEvent(fCurrentIndex); }; BaseFitEvt* InputHandlerBase::NextBaseEvent() { fCurrentIndex++; if (jointinput and fMaxEvents != -1) { while (fCurrentIndex < jointindexlow[jointindexswitch] || fCurrentIndex >= jointindexhigh[jointindexswitch]) { jointindexswitch++; // Loop Around if (jointindexswitch == jointindexlow.size()) { jointindexswitch = 0; } } if (fCurrentIndex > jointindexlow[jointindexswitch] + jointindexallowed[jointindexswitch]) { fCurrentIndex = jointindexlow[jointindexswitch]; } } return GetBaseEvent(fCurrentIndex); }; void InputHandlerBase::RegisterJointInput(std::string input, int n, TH1D* f, TH1D* e) { if (jointfluxinputs.size() == 0) { jointindexswitch = 0; fNEvents = 0; } // Push into individual input vectors jointfluxinputs.push_back((TH1D*)f->Clone()); jointeventinputs.push_back((TH1D*)e->Clone()); jointindexlow.push_back(fNEvents); jointindexhigh.push_back(fNEvents + n); fNEvents += n; // Add to the total flux/event hist if (!fFluxHist) fFluxHist = (TH1D*)f->Clone(); else fFluxHist->Add(f); if (!fEventHist) fEventHist = (TH1D*)e->Clone(); else fEventHist->Add(e); } void InputHandlerBase::SetupJointInputs() { if (jointeventinputs.size() <= 1) { jointinput = false; } else if (jointeventinputs.size() > 1) { jointinput = true; jointindexswitch = 0; } fMaxEvents = FitPar::Config().GetParI("MAXEVENTS"); if (fMaxEvents != -1 and jointeventinputs.size() > 1) { THROW("Can only handle joint inputs when config MAXEVENTS = -1!"); } if (jointeventinputs.size() > 1) { ERROR(WRN, "GiBUU sample contains multiple inputs. This will only work for " "samples that expect multi-species inputs. If this sample does, you " "can ignore this warning."); } for (size_t i = 0; i < jointeventinputs.size(); i++) { double scale = double(fNEvents) / fEventHist->Integral("width"); scale *= jointeventinputs.at(i)->Integral("width"); scale /= double(jointindexhigh[i] - jointindexlow[i]); jointindexscale.push_back(scale); } fEventHist->SetNameTitle((fName + "_EVT").c_str(), (fName + "_EVT").c_str()); fFluxHist->SetNameTitle((fName + "_FLUX").c_str(), (fName + "_FLUX").c_str()); // Setup Max Events if (fMaxEvents > 1 && fMaxEvents < fNEvents) { if (LOG_LEVEL(SAM)) { std::cout << "\t\t|-> Read Max Entries : " << fMaxEvents << std::endl; } fNEvents = fMaxEvents; } // Print out Status if (LOG_LEVEL(SAM)) { std::cout << "\t\t|-> Total Entries : " << fNEvents << std::endl << "\t\t|-> Event Integral : " << fEventHist->Integral("width") * 1.E-38 << " events/nucleon" << std::endl << "\t\t|-> Flux Integral : " << fFluxHist->Integral("width") << " /cm2" << std::endl << "\t\t|-> Event/Flux : " << fEventHist->Integral("width") * 1.E-38 / fFluxHist->Integral("width") << " cm2/nucleon" << std::endl; } } BaseFitEvt* InputHandlerBase::GetBaseEvent(const UInt_t entry) { return static_cast(GetNuisanceEvent(entry, true)); } double InputHandlerBase::GetInputWeight(int entry) { if (!jointinput) return 1.0; // Find Switch Scale while (entry < jointindexlow[jointindexswitch] || entry >= jointindexhigh[jointindexswitch]) { jointindexswitch++; // Loop Around if (jointindexswitch >= jointindexlow.size()) { jointindexswitch = 0; } } return jointindexscale[jointindexswitch]; }; diff --git a/src/MCStudies/GenericFlux_Tester.cxx b/src/MCStudies/GenericFlux_Tester.cxx index e3383b5..5834274 100644 --- a/src/MCStudies/GenericFlux_Tester.cxx +++ b/src/MCStudies/GenericFlux_Tester.cxx @@ -1,559 +1,564 @@ // 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 "GenericFlux_Tester.h" //******************************************************************** /// @brief Class to perform MC Studies on a custom measurement GenericFlux_Tester::GenericFlux_Tester(std::string name, std::string inputfile, FitWeight *rw, std::string type, std::string fakeDataFile) { //******************************************************************** // Measurement Details fName = name; eventVariables = NULL; // Define our energy range for flux calcs EnuMin = 0.; EnuMax = 100.; // Arbritrarily high energy limit // Set default fitter flags fIsDiag = true; fIsShape = false; fIsRawEvents = false; nu_4mom = new TLorentzVector(0, 0, 0, 0); pmu = new TLorentzVector(0, 0, 0, 0); ppip = new TLorentzVector(0, 0, 0, 0); ppim = new TLorentzVector(0, 0, 0, 0); ppi0 = new TLorentzVector(0, 0, 0, 0); pprot = new TLorentzVector(0, 0, 0, 0); pneut = new TLorentzVector(0, 0, 0, 0); // This function will sort out the input files automatically and parse all the // inputs,flags,etc. // There may be complex cases where you have to do this by hand, but usually // this will do. Measurement1D::SetupMeasurement(inputfile, type, rw, fakeDataFile); eventVariables = NULL; liteMode = FitPar::Config().GetParB("isLiteMode"); // Setup fDataHist as a placeholder this->fDataHist = new TH1D(("empty_data"), ("empty-data"), 1, 0, 1); this->SetupDefaultHist(); fFullCovar = StatUtils::MakeDiagonalCovarMatrix(fDataHist); covar = StatUtils::GetInvert(fFullCovar); // 1. The generator is organised in SetupMeasurement so it gives the // cross-section in "per nucleon" units. // So some extra scaling for a specific measurement may be required. For // Example to get a "per neutron" measurement on carbon // which we do here, we have to multiple by the number of nucleons 12 and // divide by the number of neutrons 6. this->fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38 / (fNEvents + 0.)) / this->TotalIntegratedFlux(); - LOG(SAM) << " Generic Flux Scaling Factor = " << fScaleFactor << std::endl; + LOG(SAM) << " Generic Flux Scaling Factor = " << fScaleFactor + << " [= " << (GetEventHistogram()->Integral("width") * 1E-38) << "/(" + << (fNEvents + 0.) << "*" << this->TotalIntegratedFlux() << ")]" + << std::endl; if (fScaleFactor <= 0.0) { ERR(WRN) << "SCALE FACTOR TOO LOW " << std::endl; sleep(20); } // Setup our TTrees this->AddEventVariablesToTree(); this->AddSignalFlagsToTree(); } void GenericFlux_Tester::AddEventVariablesToTree() { // Setup the TTree to save everything if (!eventVariables) { Config::Get().out->cd(); eventVariables = new TTree((this->fName + "_VARS").c_str(), (this->fName + "_VARS").c_str()); } LOG(SAM) << "Adding Event Variables" << std::endl; eventVariables->Branch("Mode", &Mode, "Mode/I"); eventVariables->Branch("PDGnu", &PDGnu, "PDGnu/I"); eventVariables->Branch("Enu_true", &Enu_true, "Enu_true/F"); eventVariables->Branch("Nleptons", &Nleptons, "Nleptons/I"); eventVariables->Branch("MLep", &MLep, "MLep/F"); eventVariables->Branch("ELep", &ELep, "ELep/F"); eventVariables->Branch("TLep", &TLep, "TLep/F"); eventVariables->Branch("CosLep", &CosLep, "CosLep/F"); eventVariables->Branch("CosPmuPpip", &CosPmuPpip, "CosPmuPpip/F"); eventVariables->Branch("CosPmuPpim", &CosPmuPpim, "CosPmuPpim/F"); eventVariables->Branch("CosPmuPpi0", &CosPmuPpi0, "CosPmuPpi0/F"); eventVariables->Branch("CosPmuPprot", &CosPmuPprot, "CosPmuPprot/F"); eventVariables->Branch("CosPmuPneut", &CosPmuPneut, "CosPmuPneut/F"); eventVariables->Branch("Nprotons", &Nprotons, "Nprotons/I"); eventVariables->Branch("MPr", &MPr, "MPr/F"); eventVariables->Branch("EPr", &EPr, "EPr/F"); eventVariables->Branch("TPr", &TPr, "TPr/F"); eventVariables->Branch("CosPr", &CosPr, "CosPr/F"); eventVariables->Branch("CosPprotPneut", &CosPprotPneut, "CosPprotPneut/F"); eventVariables->Branch("Nneutrons", &Nneutrons, "Nneutrons/I"); eventVariables->Branch("MNe", &MNe, "MNe/F"); eventVariables->Branch("ENe", &ENe, "ENe/F"); eventVariables->Branch("TNe", &TNe, "TNe/F"); eventVariables->Branch("CosNe", &CosNe, "CosNe/F"); eventVariables->Branch("Npiplus", &Npiplus, "Npiplus/I"); eventVariables->Branch("MPiP", &MPiP, "MPiP/F"); eventVariables->Branch("EPiP", &EPiP, "EPiP/F"); eventVariables->Branch("TPiP", &TPiP, "TPiP/F"); eventVariables->Branch("CosPiP", &CosPiP, "CosPiP/F"); eventVariables->Branch("CosPpipPprot", &CosPpipPprot, "CosPpipProt/F"); eventVariables->Branch("CosPpipPneut", &CosPpipPneut, "CosPpipPneut/F"); eventVariables->Branch("CosPpipPpim", &CosPpipPpim, "CosPpipPpim/F"); eventVariables->Branch("CosPpipPpi0", &CosPpipPpi0, "CosPpipPpi0/F"); eventVariables->Branch("Npineg", &Npineg, "Npineg/I"); eventVariables->Branch("MPiN", &MPiN, "MPiN/F"); eventVariables->Branch("EPiN", &EPiN, "EPiN/F"); eventVariables->Branch("TPiN", &TPiN, "TPiN/F"); eventVariables->Branch("CosPiN", &CosPiN, "CosPiN/F"); eventVariables->Branch("CosPpimPprot", &CosPpimPprot, "CosPpimPprot/F"); eventVariables->Branch("CosPpimPneut", &CosPpimPneut, "CosPpimPneut/F"); eventVariables->Branch("CosPpimPpi0", &CosPpimPpi0, "CosPpimPpi0/F"); eventVariables->Branch("Npi0", &Npi0, "Npi0/I"); eventVariables->Branch("MPi0", &MPi0, "MPi0/F"); eventVariables->Branch("EPi0", &EPi0, "EPi0/F"); eventVariables->Branch("TPi0", &TPi0, "TPi0/F"); eventVariables->Branch("CosPi0", &CosPi0, "CosPi0/F"); eventVariables->Branch("CosPi0Pprot", &CosPi0Pprot, "CosPi0Pprot/F"); eventVariables->Branch("CosPi0Pneut", &CosPi0Pneut, "CosPi0Pneut/F"); eventVariables->Branch("Nother", &Nother, "Nother/I"); eventVariables->Branch("Q2_true", &Q2_true, "Q2_true/F"); eventVariables->Branch("q0_true", &q0_true, "q0_true/F"); eventVariables->Branch("q3_true", &q3_true, "q3_true/F"); eventVariables->Branch("Enu_QE", &Enu_QE, "Enu_QE/F"); eventVariables->Branch("Q2_QE", &Q2_QE, "Q2_QE/F"); eventVariables->Branch("W_nuc_rest", &W_nuc_rest, "W_nuc_rest/F"); eventVariables->Branch("bjorken_x", &bjorken_x, "bjorken_x/F"); eventVariables->Branch("bjorken_y", &bjorken_y, "bjorken_y/F"); eventVariables->Branch("Erecoil_true", &Erecoil_true, "Erecoil_true/F"); eventVariables->Branch("Erecoil_charged", &Erecoil_charged, "Erecoil_charged/F"); eventVariables->Branch("Erecoil_minerva", &Erecoil_minerva, "Erecoil_minerva/F"); if (!liteMode) { eventVariables->Branch("nu_4mom", &nu_4mom); eventVariables->Branch("pmu_4mom", &pmu); eventVariables->Branch("hm_ppip_4mom", &ppip); eventVariables->Branch("hm_ppim_4mom", &ppim); eventVariables->Branch("hm_ppi0_4mom", &ppi0); eventVariables->Branch("hm_pprot_4mom", &pprot); eventVariables->Branch("hm_pneut_4mom", &pneut); } // Event Scaling Information eventVariables->Branch("Weight", &Weight, "Weight/F"); eventVariables->Branch("InputWeight", &InputWeight, "InputWeight/F"); eventVariables->Branch("RWWeight", &RWWeight, "RWWeight/F"); eventVariables->Branch("FluxWeight", &FluxWeight, "FluxWeight/F"); eventVariables->Branch("fScaleFactor", &fScaleFactor, "fScaleFactor/D"); return; } void GenericFlux_Tester::AddSignalFlagsToTree() { if (!eventVariables) { Config::Get().out->cd(); eventVariables = new TTree((this->fName + "_VARS").c_str(), (this->fName + "_VARS").c_str()); } LOG(SAM) << "Adding Samples" << std::endl; // Signal Definitions from SignalDef.cxx eventVariables->Branch("flagCCINC", &flagCCINC, "flagCCINC/O"); eventVariables->Branch("flagNCINC", &flagNCINC, "flagNCINC/O"); eventVariables->Branch("flagCCQE", &flagCCQE, "flagCCQE/O"); eventVariables->Branch("flagCC0pi", &flagCC0pi, "flagCC0pi/O"); eventVariables->Branch("flagCCQELike", &flagCCQELike, "flagCCQELike/O"); eventVariables->Branch("flagNCEL", &flagNCEL, "flagNCEL/O"); eventVariables->Branch("flagNC0pi", &flagNC0pi, "flagNC0pi/O"); eventVariables->Branch("flagCCcoh", &flagCCcoh, "flagCCcoh/O"); eventVariables->Branch("flagNCcoh", &flagNCcoh, "flagNCcoh/O"); eventVariables->Branch("flagCC1pip", &flagCC1pip, "flagCC1pip/O"); eventVariables->Branch("flagNC1pip", &flagNC1pip, "flagNC1pip/O"); eventVariables->Branch("flagCC1pim", &flagCC1pim, "flagCC1pim/O"); eventVariables->Branch("flagNC1pim", &flagNC1pim, "flagNC1pim/O"); eventVariables->Branch("flagCC1pi0", &flagCC1pi0, "flagCC1pi0/O"); eventVariables->Branch("flagNC1pi0", &flagNC1pi0, "flagNC1pi0/O"); }; //******************************************************************** void GenericFlux_Tester::FillEventVariables(FitEvent *event) { //******************************************************************** // Fill Signal Variables FillSignalFlags(event); LOG(DEB) << "Filling signal" << std::endl; // Function used to extract any variables of interest to the event Mode = event->Mode; Nleptons = 0; Nparticles = 0; PDGnu = 0; PDGLep = 0; Enu_true = Enu_QE = Q2_true = Q2_QE = TLep = TPr = TNe = TPiP = TPiN = TPi0 = -999.9; Nprotons = 0; PPr = EPr = MPr = CosPr = -999.9; Nneutrons = 0; PNe = ENe = MNe = CosNe = -999.9; Npiplus = 0; PPiP = EPiP = MPiP = CosPiP = -999.9; Npineg = 0; PPiN = EPiN = MPiN = CosPiN = -999.9; Npi0 = 0; PPi0 = EPi0 = MPi0 = CosPi0 = -999.9; // All of the angles Clarence added CosPmuPpip = CosPmuPpim = CosPmuPpi0 = CosPmuPprot = CosPmuPneut = CosPpipPprot = CosPpipPneut = CosPpipPpim = CosPpipPpi0 = CosPpimPprot = CosPpimPneut = CosPpimPpi0 = CosPi0Pprot = CosPi0Pneut = CosPprotPneut = -999.9; float proton_highmom = -999.9; float neutron_highmom = -999.9; float piplus_highmom = -999.9; float pineg_highmom = -999.9; float pi0_highmom = -999.9; (*nu_4mom) = event->PartInfo(0)->fP; if (!liteMode) { (*pmu) = TLorentzVector(0, 0, 0, 0); (*ppip) = TLorentzVector(0, 0, 0, 0); (*ppim) = TLorentzVector(0, 0, 0, 0); (*ppi0) = TLorentzVector(0, 0, 0, 0); (*pprot) = TLorentzVector(0, 0, 0, 0); (*pneut) = TLorentzVector(0, 0, 0, 0); } Enu_true = nu_4mom->E(); PDGnu = event->PartInfo(0)->fPID; bool cc = (abs(event->Mode) < 30); (void)cc; // Add all pion distributions for the event. // Add classifier for CC0pi or CC1pi or CCOther // Save Modes Properly // Save low recoil measurements // Start Particle Loop UInt_t npart = event->Npart(); for (UInt_t i = 0; i < npart; i++) { // Skip particles that weren't in the final state - bool part_alive = event->PartInfo(i)->fIsAlive and event->PartInfo(i)->Status() == kFinalState; + bool part_alive = event->PartInfo(i)->fIsAlive and + event->PartInfo(i)->Status() == kFinalState; if (!part_alive) continue; // PDG Particle int PDGpart = event->PartInfo(i)->fPID; TLorentzVector part_4mom = event->PartInfo(i)->fP; Nparticles++; // Get Charged Lepton if (abs(PDGpart) == abs(PDGnu) - 1) { Nleptons++; PDGLep = PDGpart; TLep = FitUtils::T(part_4mom) * 1000.0; PLep = (part_4mom.Vect().Mag()); ELep = (part_4mom.E()); MLep = (part_4mom.Mag()); CosLep = cos(part_4mom.Vect().Angle(nu_4mom->Vect())); (*pmu) = part_4mom; Q2_true = -1 * (part_4mom - (*nu_4mom)).Mag2(); float ThetaLep = (event->PartInfo(0)) ->fP.Vect() .Angle((event->PartInfo(i))->fP.Vect()); q0_true = (part_4mom - (*nu_4mom)).E(); q3_true = (part_4mom - (*nu_4mom)).Vect().Mag(); // Get W_true with assumption of initial state nucleon at rest float m_n = (float)PhysConst::mass_proton * 1000.; W_nuc_rest = sqrt(-Q2_true + 2 * m_n * (Enu_true - ELep) + m_n * m_n); // Get the Bjorken x and y variables // Assume that E_had = Enu - Emu as in MINERvA bjorken_x = Q2_true / (2 * m_n * (Enu_true - ELep)); bjorken_y = 1 - ELep / Enu_true; // Quasi-elastic ---------------------- // ------------------------------------ // Q2 QE Assuming Carbon Input. Should change this to be dynamic soon. Q2_QE = FitUtils::Q2QErec(part_4mom, cos(ThetaLep), 34., true) * 1000000.0; Enu_QE = FitUtils::EnuQErec(part_4mom, cos(ThetaLep), 34., true) * 1000.0; // Pion Production ---------------------- // -------------------------------------- } else if (PDGpart == 2212) { Nprotons++; if (part_4mom.Vect().Mag() > proton_highmom) { proton_highmom = part_4mom.Vect().Mag(); PPr = (part_4mom.Vect().Mag()); EPr = (part_4mom.E()); TPr = FitUtils::T(part_4mom) * 1000.; MPr = (part_4mom.Mag()); CosPr = cos(part_4mom.Vect().Angle(nu_4mom->Vect())); (*pprot) = part_4mom; } } else if (PDGpart == 2112) { Nneutrons++; if (part_4mom.Vect().Mag() > neutron_highmom) { neutron_highmom = part_4mom.Vect().Mag(); PNe = (part_4mom.Vect().Mag()); ENe = (part_4mom.E()); TNe = FitUtils::T(part_4mom) * 1000.; MNe = (part_4mom.Mag()); CosNe = cos(part_4mom.Vect().Angle(nu_4mom->Vect())); (*pneut) = part_4mom; } } else if (PDGpart == 211) { Npiplus++; if (part_4mom.Vect().Mag() > piplus_highmom) { piplus_highmom = part_4mom.Vect().Mag(); PPiP = (part_4mom.Vect().Mag()); EPiP = (part_4mom.E()); TPiP = FitUtils::T(part_4mom) * 1000.; MPiP = (part_4mom.Mag()); CosPiP = cos(part_4mom.Vect().Angle(nu_4mom->Vect())); (*ppip) = part_4mom; } } else if (PDGpart == -211) { Npineg++; if (part_4mom.Vect().Mag() > pineg_highmom) { pineg_highmom = part_4mom.Vect().Mag(); PPiN = (part_4mom.Vect().Mag()); EPiN = (part_4mom.E()); TPiN = FitUtils::T(part_4mom) * 1000.; MPiN = (part_4mom.Mag()); CosPiN = cos(part_4mom.Vect().Angle(nu_4mom->Vect())); (*ppim) = part_4mom; } } else if (PDGpart == 111) { Npi0++; if (part_4mom.Vect().Mag() > pi0_highmom) { pi0_highmom = part_4mom.Vect().Mag(); PPi0 = (part_4mom.Vect().Mag()); EPi0 = (part_4mom.E()); TPi0 = FitUtils::T(part_4mom) * 1000.; MPi0 = (part_4mom.Mag()); CosPi0 = cos(part_4mom.Vect().Angle(nu_4mom->Vect())); (*ppi0) = part_4mom; } } else { Nother++; } } // Get Recoil Definitions ------ // ----------------------------- Erecoil_true = FitUtils::GetErecoil_TRUE(event); Erecoil_charged = FitUtils::GetErecoil_CHARGED(event); Erecoil_minerva = FitUtils::GetErecoil_MINERvA_LowRecoil(event); // Do the angles between final state particles if (Nleptons > 0 && Npiplus > 0) CosPmuPpip = cos(pmu->Vect().Angle(ppip->Vect())); if (Nleptons > 0 && Npineg > 0) CosPmuPpim = cos(pmu->Vect().Angle(ppim->Vect())); if (Nleptons > 0 && Npi0 > 0) CosPmuPpi0 = cos(pmu->Vect().Angle(ppi0->Vect())); if (Nleptons > 0 && Nprotons > 0) CosPmuPprot = cos(pmu->Vect().Angle(pprot->Vect())); if (Nleptons > 0 && Nneutrons > 0) CosPmuPneut = cos(pmu->Vect().Angle(pneut->Vect())); if (Npiplus > 0 && Nprotons > 0) CosPpipPprot = cos(ppip->Vect().Angle(pprot->Vect())); if (Npiplus > 0 && Nneutrons > 0) CosPpipPneut = cos(ppip->Vect().Angle(pneut->Vect())); if (Npiplus > 0 && Npineg > 0) CosPpipPpim = cos(ppip->Vect().Angle(ppim->Vect())); if (Npiplus > 0 && Npi0 > 0) CosPpipPpi0 = cos(ppip->Vect().Angle(ppi0->Vect())); if (Npineg > 0 && Nprotons > 0) CosPpimPprot = cos(ppim->Vect().Angle(pprot->Vect())); if (Npineg > 0 && Nneutrons > 0) CosPpimPneut = cos(ppim->Vect().Angle(pneut->Vect())); if (Npineg > 0 && Npi0 > 0) CosPpimPpi0 = cos(ppim->Vect().Angle(ppi0->Vect())); if (Npi0 > 0 && Nprotons > 0) CosPi0Pprot = cos(ppi0->Vect().Angle(pprot->Vect())); if (Npi0 > 0 && Nneutrons > 0) CosPi0Pneut = cos(ppi0->Vect().Angle(pneut->Vect())); if (Nprotons > 0 && Nneutrons > 0) CosPprotPneut = cos(pprot->Vect().Angle(pneut->Vect())); // Event Weights ---- // ------------------ Weight = event->RWWeight * event->InputWeight; RWWeight = event->RWWeight; InputWeight = event->InputWeight; FluxWeight = - GetFluxHistogram()->GetBinContent(GetFluxHistogram()->FindBin(Enu)) / GetFluxHistogram()->Integral(); + GetFluxHistogram()->GetBinContent(GetFluxHistogram()->FindBin(Enu)) / + GetFluxHistogram()->Integral(); xsecScaling = fScaleFactor; if (fScaleFactor <= 0.0) { ERR(WRN) << "SCALE FACTOR TOO LOW " << std::endl; sleep(20); } // Fill the eventVariables Tree eventVariables->Fill(); return; }; //******************************************************************** void GenericFlux_Tester::Write(std::string drawOpt) { //******************************************************************** // First save the TTree eventVariables->Write(); // Save Flux and Event Histograms too GetInput()->GetFluxHistogram()->Write(); GetInput()->GetEventHistogram()->Write(); return; } //******************************************************************** void GenericFlux_Tester::FillSignalFlags(FitEvent *event) { //******************************************************************** // Some example flags are given from SignalDef. // See src/Utils/SignalDef.cxx for more. int nuPDG = event->PartInfo(0)->fPID; // Generic signal flags flagCCINC = SignalDef::isCCINC(event, nuPDG); flagNCINC = SignalDef::isNCINC(event, nuPDG); flagCCQE = SignalDef::isCCQE(event, nuPDG); flagCCQELike = SignalDef::isCCQELike(event, nuPDG); flagCC0pi = SignalDef::isCC0pi(event, nuPDG); flagNCEL = SignalDef::isNCEL(event, nuPDG); flagNC0pi = SignalDef::isNC0pi(event, nuPDG); flagCCcoh = SignalDef::isCCCOH(event, nuPDG, 211); flagNCcoh = SignalDef::isNCCOH(event, nuPDG, 111); flagCC1pip = SignalDef::isCC1pi(event, nuPDG, 211); flagNC1pip = SignalDef::isNC1pi(event, nuPDG, 211); flagCC1pim = SignalDef::isCC1pi(event, nuPDG, -211); flagNC1pim = SignalDef::isNC1pi(event, nuPDG, -211); flagCC1pi0 = SignalDef::isCC1pi(event, nuPDG, 111); flagNC1pi0 = SignalDef::isNC1pi(event, nuPDG, 111); } // ------------------------------------------------------------------- // Purely MC Plot // Following functions are just overrides to handle this // ------------------------------------------------------------------- //******************************************************************** /// Everything is classed as signal... bool GenericFlux_Tester::isSignal(FitEvent *event) { //******************************************************************** (void)event; return true; }; //******************************************************************** void GenericFlux_Tester::ScaleEvents() { //******************************************************************** // Saving everything to a TTree so no scaling required return; } //******************************************************************** void GenericFlux_Tester::ApplyNormScale(float norm) { //******************************************************************** // Saving everything to a TTree so no scaling required this->fCurrentNorm = norm; return; } //******************************************************************** void GenericFlux_Tester::FillHistograms() { //******************************************************************** // No Histograms need filling........ return; } //******************************************************************** void GenericFlux_Tester::ResetAll() { //******************************************************************** eventVariables->Reset(); return; } //******************************************************************** float GenericFlux_Tester::GetChi2() { //******************************************************************** // No Likelihood to test, purely MC return 0.0; }