diff --git a/app/PrepareGENIE.cxx b/app/PrepareGENIE.cxx index 0e74407..4452750 100644 --- a/app/PrepareGENIE.cxx +++ b/app/PrepareGENIE.cxx @@ -1,615 +1,635 @@ #include #include #include "FitLogger.h" #include "PlotUtils.h" #include "TFile.h" #include "TH1D.h" #include "TTree.h" #ifdef __GENIE_ENABLED__ #include "Conventions/Units.h" #include "GHEP/GHepParticle.h" #include "PDG/PDGUtils.h" #endif std::string gInputFiles = ""; std::string gOutputFile = ""; std::string gFluxFile = ""; std::string gTarget = ""; double MonoEnergy; bool IsMonoE = false; void PrintOptions(); void ParseOptions(int argc, char* argv[]); void RunGENIEPrepareMono(std::string input, std::string target, std::string output); void RunGENIEPrepare(std::string input, std::string flux, std::string target, std::string output); int main(int argc, char* argv[]) { ParseOptions(argc, argv); if (IsMonoE) { RunGENIEPrepareMono(gInputFiles, gTarget, gOutputFile); } else { RunGENIEPrepare(gInputFiles, gFluxFile, gTarget, gOutputFile); } } void RunGENIEPrepareMono(std::string input, std::string target, std::string output) { - std::cout << "Running in mono" << std::endl; + LOG(FIT) << "Running in mono energetic with E = " << MonoEnergy << " GeV" << std::endl; // 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); + // Have the TH1D go from MonoEnergy/2 to MonoEnergy/2 + TH1D* fluxhist = new TH1D("flux", "flux", 1000, MonoEnergy/2., MonoEnergy*2.); fluxhist->Fill(MonoEnergy); fluxhist->Scale(1, "width"); // Make Event Hist TH1D* eventhist = (TH1D*)fluxhist->Clone(); eventhist->Reset(); TH1D* xsechist = (TH1D*)eventhist->Clone(); // Create maps std::map 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(); size_t freq = nevt / 20; if (freq && !(i % freq)) { 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(); + TFile* outputfile; - LOG(FIT) << "Getting splines in mono" << std::endl; + if (!gOutputFile.length()) { + tn->GetEntry(0); + outputfile = tn->GetFile(); + outputfile->cd(); + } else { + outputfile = new TFile(gOutputFile.c_str(), "RECREATE"); + outputfile->cd(); + + QLOG(FIT, "Cloning input vector to output file: " << gOutputFile); + TTree* cloneTree = tn->CloneTree(); + cloneTree->SetDirectory(outputfile); + cloneTree->Write(); + QLOG(FIT, "Done."); + } + + LOG(FIT) << "Getting splines in mono-energetic..." << 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); } } } + LOG(FIT) << "Total XSec Integral = " << totalxsec->Integral("width") << std::endl; outputfile->cd(); - totalxsec->Write("nuisance_Xsec", TObject::kOverwrite); + totalxsec->Write("nuisance_xsec", TObject::kOverwrite); eventhist = (TH1D*)totalxsec->Clone(); eventhist->Multiply(fluxhist); + 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") + LOG(FIT) << "XSec Hist Integral = " << xsechist->Integral("width") << std::endl; + outputfile->Write(); + outputfile->Close(); + return; } void RunGENIEPrepare(std::string input, std::string flux, std::string target, std::string output) { LOG(FIT) << "Running GENIE Prepare" << std::endl; - std::cout << "Running in prepare" << std::endl; + LOG(FIT) << "Running in prepare" << std::endl; // Get Flux Hist std::vector fluxvect = GeneralUtils::ParseToStr(flux, ","); TH1* fluxhist = NULL; if (fluxvect.size() == 3) { double from = GeneralUtils::StrToDbl(fluxvect[0]); double to = GeneralUtils::StrToDbl(fluxvect[1]); double step = GeneralUtils::StrToDbl(fluxvect[2]); int nstep = ceil((to - from) / step); to = from + step * nstep; QLOG(FIT, "Generating flat flux histogram from " << from << " to " << to << " with bins " << step << " wide (NBins = " << nstep << ")."); fluxhist = new TH1D("spectrum", ";E_{#nu} (GeV);Count (A.U.)", nstep, from, to); for (Int_t bi_it = 1; bi_it < fluxhist->GetXaxis()->GetNbins(); ++bi_it) { fluxhist->SetBinContent(bi_it, 1.0 / double(step * nstep)); } fluxhist->SetDirectory(0); } else if (fluxvect.size() == 2) { TFile* fluxfile = new TFile(fluxvect[0].c_str(), "READ"); if (!fluxfile->IsZombie()) { fluxhist = dynamic_cast(fluxfile->Get(fluxvect[1].c_str())); if (!fluxhist) { ERR(FTL) << "Couldn't find histogram named: \"" << fluxvect[1] << "\" in file: \"" << fluxvect[0] << std::endl; throw; } fluxhist->SetDirectory(0); } } else if (fluxvect.size() == 1) { MonoEnergy = GeneralUtils::StrToDbl(fluxvect[0]); RunGENIEPrepareMono(input, target, output); return; } else { LOG(FTL) << "Bad flux specification: \"" << flux << "\"." << std::endl; throw; } // Setup TTree TChain* tn = new TChain("gtree"); if (input.find_first_of(',') != std::string::npos) { std::vector inputvect = GeneralUtils::ParseToStr(input, ","); for (size_t iv_it = 0; iv_it < inputvect.size(); ++iv_it) { tn->AddFile(inputvect[iv_it].c_str()); QLOG(FIT, "Added input file: " << inputvect[iv_it]); } } else { // The Add form can accept wildcards. tn->Add(input.c_str()); } int nevt = tn->GetEntries(); if (!nevt) { THROW("Couldn't load any events from input specification: \"" << input.c_str() << "\""); } else { QLOG(FIT, "Found " << nevt << " input entries."); } 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()); if (i % (nevt / 20) == 0) { LOG(FIT) << "Processed " << i << "/" << nevt << " GENIE events (Last event: { E: " << neu->E() << ", xsec: " << xsec << " }." << std::endl; } // Clear Event genientpl->Clear(); } LOG(FIT) << "Processed all events" << std::endl; // Once event loop is done we can start saving stuff into the file TFile* outputfile; if (!gOutputFile.length()) { tn->GetEntry(0); outputfile = tn->GetFile(); outputfile->cd(); } else { outputfile = new TFile(gOutputFile.c_str(), "RECREATE"); outputfile->cd(); QLOG(FIT, "Cloning input vector to output file: " << gOutputFile); TTree* cloneTree = tn->CloneTree(); cloneTree->SetDirectory(outputfile); cloneTree->Write(); QLOG(FIT, "Done."); } 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; + 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") + LOG(FIT) << "XSec Hist Integral = " << xsechist->Integral("width") << std::endl; outputfile->Write(); outputfile->Close(); - delete outputfile; 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 << " [ -f elow,ehigh,estep ] : Energy range specification when no " "flux file was used." << std::endl; std::cout << " [ -t target ] : Target that GHepRecords were generated with. " "Comma seperated list. E.g. for CH2 " "target=1000060120,1000010010,1000010010" << std::endl; std::cout << " [ -o outputfile.root ] : File to write prepared input file to." << std::endl; std::cout << " [ -m Mono_E_nu_GeV ] : Run in mono-energetic mode." << std::endl; } void ParseOptions(int argc, char* argv[]) { bool flagopt = false; // If No Arguments print commands for (int i = 1; i < argc; ++i) { if (!std::strcmp(argv[i], "-h")) { flagopt = true; break; } if (i + 1 != argc) { // Cardfile if (!std::strcmp(argv[i], "-h")) { flagopt = true; break; } else if (!std::strcmp(argv[i], "-i")) { gInputFiles = argv[i + 1]; ++i; } else if (!std::strcmp(argv[i], "-o")) { gOutputFile = argv[i + 1]; ++i; } else if (!std::strcmp(argv[i], "-f")) { gFluxFile = argv[i + 1]; ++i; } else if (!std::strcmp(argv[i], "-t")) { gTarget = argv[i + 1]; ++i; } else if (!std::strcmp(argv[i], "-m")) { MonoEnergy = GeneralUtils::StrToDbl(argv[i + 1]); IsMonoE = true; ++i; } else { ERR(FTL) << "ERROR: unknown command line option given! - '" << argv[i] << " " << argv[i + 1] << "'" << std::endl; PrintOptions(); break; } } } if (gInputFiles == "" && !flagopt) { ERR(FTL) << "No input file(s) specified!" << std::endl; flagopt = true; } if (gFluxFile == "" && !flagopt && !IsMonoE) { ERR(FTL) << "No flux input specified for Prepare Mode" << std::endl; flagopt = true; } if (gTarget == "" && !flagopt) { ERR(FTL) << "No target specified for Prepare Mode" << std::endl; flagopt = true; } if (argc < 1 || flagopt) { PrintOptions(); exit(-1); } return; } diff --git a/freshbuild.sh b/freshbuild.sh index 8018e18..1ee2e60 100644 --- a/freshbuild.sh +++ b/freshbuild.sh @@ -1,5 +1,11 @@ +#!/bin/bash + + mkdir build cd ./build -cmake -DUSE_NuWro=0 -DUSE_NIWG=0 -DUSE_NEUT=1 -DUSE_GENIE=0 -DUSE_T2K=0 -DUSE_GiBUU=0 -DUSE_NUANCE=0 -DUSE_MINIMIZER=false ../ +if [[ -e CMakeCache.txt ]]; then + rm -v CMakeCache.txt +fi +cmake -DUSE_NuWro=0 -DUSE_NIWG=0 -DUSE_NEUT=1 -DUSE_GENIE=1 -DUSE_T2K=0 -DUSE_GiBUU=0 -DUSE_NUANCE=0 -DUSE_MINIMIZER=false ../ make clean && make -j8 cd - diff --git a/parameters/config.xml b/parameters/config.xml index 3feb65d..c11de34 100644 --- a/parameters/config.xml +++ b/parameters/config.xml @@ -1,202 +1,201 @@ - diff --git a/src/MCStudies/GenericFlux_Tester.cxx b/src/MCStudies/GenericFlux_Tester.cxx index f81dc1a..1d28b4c 100644 --- a/src/MCStudies/GenericFlux_Tester.cxx +++ b/src/MCStudies/GenericFlux_Tester.cxx @@ -1,599 +1,596 @@ // 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 + EnuMax = 1E10; // 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 = Config::Get().GetParB("isLiteMode"); if (Config::HasPar("EnuMin")) { EnuMin = Config::GetParD("EnuMin"); } if (Config::HasPar("EnuMax")) { EnuMax = Config::GetParD("EnuMax"); } // 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. // N.B. MeasurementBase::PredictedEventRate includes the 1E-38 factor that is // often included here in other classes that directly integrate the event // histogram. This method is used here as it now respects EnuMin and EnuMax // correctly. this->fScaleFactor = (this->PredictedEventRate("width", 0, 1000) / double(fNEvents)) / this->TotalIntegratedFlux(); 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); + throw; } // 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"); // all sensible eventVariables->Branch("MLep", &MLep, "MLep/F"); eventVariables->Branch("ELep", &ELep, "ELep/F"); // negative -999 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; + LOG(SAM) << "Adding signal flags" << 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::ResetVariables() { //******************************************************************** // Reset neutrino PDG PDGnu = 0; // Reset energies Enu_true = Enu_QE = __BAD_FLOAT__; // Reset auxillaries Q2_true = Q2_QE = W_nuc_rest = bjorken_x = bjorken_y = q0_true = q3_true = Erecoil_true = Erecoil_charged = Erecoil_minerva = __BAD_FLOAT__; // Reset particle counters Nparticles = Nleptons = Nother = Nprotons = Nneutrons = Npiplus = Npineg = Npi0 = 0; // Reset Lepton PDG PDGLep = 0; // Reset Lepton variables TLep = CosLep = ELep = PLep = MLep = __BAD_FLOAT__; // Rset proton variables PPr = CosPr = EPr = TPr = MPr = __BAD_FLOAT__; // Reset neutron variables PNe = CosNe = ENe = TNe = MNe = __BAD_FLOAT__; // Reset pi+ variables PPiP = CosPiP = EPiP = TPiP = MPiP = __BAD_FLOAT__; // Reset pi- variables PPiN = CosPiN = EPiN = TPiN = MPiN = __BAD_FLOAT__; // Reset pi0 variables PPi0 = CosPi0 = EPi0 = TPi0 = MPi0 = __BAD_FLOAT__; // Reset the cos angles CosPmuPpip = CosPmuPpim = CosPmuPpi0 = CosPmuPprot = CosPmuPneut = CosPpipPprot = CosPpipPneut = CosPpipPpim = CosPpipPpi0 = CosPpimPprot = CosPpimPneut = CosPpimPpi0 = CosPi0Pprot = CosPi0Pneut = CosPprotPneut = __BAD_FLOAT__; } //******************************************************************** void GenericFlux_Tester::FillEventVariables(FitEvent *event) { //******************************************************************** // Fill Signal Variables FillSignalFlags(event); LOG(DEB) << "Filling signal" << std::endl; // Reset the private variables (see header) ResetVariables(); // Function used to extract any variables of interest to the event Mode = event->Mode; // Reset the highest momentum variables float proton_highmom = __BAD_FLOAT__; float neutron_highmom = __BAD_FLOAT__; float piplus_highmom = __BAD_FLOAT__; float pineg_highmom = __BAD_FLOAT__; float pi0_highmom = __BAD_FLOAT__; (*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; - if (!part_alive) - continue; + 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(); - xsecScaling = fScaleFactor; - if (fScaleFactor <= 0.0) { ERR(WRN) << "SCALE FACTOR TOO LOW " << std::endl; - sleep(20); + throw; } // 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; } diff --git a/src/MCStudies/GenericFlux_Tester.h b/src/MCStudies/GenericFlux_Tester.h index d050aa5..a14cc0e 100644 --- a/src/MCStudies/GenericFlux_Tester.h +++ b/src/MCStudies/GenericFlux_Tester.h @@ -1,201 +1,199 @@ // 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 . *******************************************************************************/ #ifndef GenericFlux_Tester_H_SEEN #define GenericFlux_Tester_H_SEEN #include "Measurement1D.h" #ifndef __BAD__FLOAT__ #define __BAD_FLOAT__ -999.99 #endif //******************************************************************** class GenericFlux_Tester : public Measurement1D { //******************************************************************** public: GenericFlux_Tester(std::string name, std::string inputfile, FitWeight *rw, std::string type, std::string fakeDataFile); virtual ~GenericFlux_Tester() {}; //! Clear private variables inline void ResetVariables(); //! Grab info from event void FillEventVariables(FitEvent *event); //! Fill Custom Histograms void FillHistograms(); //! ResetAll void ResetAll(); //! Scale void ScaleEvents(); //! Norm void ApplyNormScale(float norm); //! Define this samples signal bool isSignal(FitEvent *nvect); //! Write Files void Write(std::string drawOpt); //! Get Chi2 float GetChi2(); //! Fill all signal flags we currently have void FillSignalFlags(FitEvent *event); void AddEventVariablesToTree(); void AddSignalFlagsToTree(); private: // Lighter flat trees that don't include vectors bool liteMode; TTree* eventVariables; TLorentzVector *nu_4mom; TLorentzVector *pmu; TLorentzVector *ppip; TLorentzVector *ppim; TLorentzVector *ppi0; TLorentzVector *pprot; TLorentzVector *pneut; // Saved Variables float Enu_true; float Enu_QE; int PDGnu; // Auxillairies float Q2_true; float Q2_QE; float W_nuc_rest; float bjorken_x; float bjorken_y; float q0_true; float q3_true; float Erecoil_true; float Erecoil_charged; float Erecoil_minerva; // Interaction mode int Mode; // Particle counters int Nparticles; int Nleptons; int Nother; int Nprotons; int Nneutrons; int Npiplus; int Npineg; int Npi0; // Lepton variables int PDGLep; float TLep; float CosLep; float ELep; float PLep; float MLep; // Proton variables float PPr; //!< Highest Mom Proton float CosPr; //!< Highest Mom Proton float EPr; float TPr; float MPr; // Neutron variables float PNe; float CosNe; float ENe; float TNe; float MNe; // Pi+ variables float PPiP; float CosPiP; float EPiP; float TPiP; float MPiP; // Pi- variables float PPiN; float CosPiN; float EPiN; float TPiN; float MPiN; float PPi0; float CosPi0; float EPi0; float TPi0; float MPi0; // Angular variables float CosPmuPpip; float CosPmuPpim; float CosPmuPpi0; float CosPmuPprot; float CosPmuPneut; float CosPpipPprot; float CosPpipPneut; float CosPpipPpim; float CosPpipPpi0; float CosPpimPprot; float CosPpimPneut; float CosPpimPpi0; float CosPi0Pprot; float CosPi0Pneut; float CosPprotPneut; // Weights float Weight; float RWWeight; float InputWeight; float FluxWeight; // Generic signal flags bool flagCCINC; bool flagNCINC; bool flagCCQE; bool flagCC0pi; bool flagCCQELike; bool flagNCEL; bool flagNC0pi; bool flagCCcoh; bool flagNCcoh; bool flagCC1pip; bool flagNC1pip; bool flagCC1pim; bool flagNC1pim; bool flagCC1pi0; bool flagNC1pi0; - float xsecScaling; - }; #endif diff --git a/src/MCStudies/GenericFlux_Vectors.cxx b/src/MCStudies/GenericFlux_Vectors.cxx index 7e131b2..d0f8613 100644 --- a/src/MCStudies/GenericFlux_Vectors.cxx +++ b/src/MCStudies/GenericFlux_Vectors.cxx @@ -1,247 +1,280 @@ // 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_Vectors.h" GenericFlux_Vectors::GenericFlux_Vectors(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 + EnuMax = 1E10; // Arbritrarily high energy limit if (Config::HasPar("EnuMin")) { EnuMin = Config::GetParD("EnuMin"); } if (Config::HasPar("EnuMax")) { EnuMax = Config::GetParD("EnuMax"); } // Set default fitter flags fIsDiag = true; fIsShape = false; fIsRawEvents = false; // 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; // 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. // N.B. MeasurementBase::PredictedEventRate includes the 1E-38 factor that is // often included here in other classes that directly integrate the event // histogram. This method is used here as it now respects EnuMin and EnuMax // correctly. this->fScaleFactor = - (this->PredictedEventRate("width", 0, 1000) / double(fNEvents)) / + (this->PredictedEventRate("width", 0, EnuMax) / double(fNEvents)) / 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; + throw; } // Setup our TTrees this->AddEventVariablesToTree(); + this->AddSignalFlagsToTree(); } void GenericFlux_Vectors::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("cc", &cc, "cc/B"); eventVariables->Branch("PDGnu", &PDGnu, "PDGnu/I"); eventVariables->Branch("Enu_true", &Enu_true, "Enu_true/F"); eventVariables->Branch("tgt", &tgt, "tgt/I"); eventVariables->Branch("PDGLep", &PDGLep, "PDGLep/I"); eventVariables->Branch("ELep", &ELep, "ELep/F"); eventVariables->Branch("CosLep", &CosLep, "CosLep/F"); // Basic interaction kinematics eventVariables->Branch("Q2", &Q2, "Q2/F"); eventVariables->Branch("q0", &q0, "q0/F"); eventVariables->Branch("q3", &q3, "q3/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("W", &W, "W/F"); eventVariables->Branch("x", &x, "x/F"); eventVariables->Branch("y", &y, "y/F"); // Save outgoing particle vectors eventVariables->Branch("nfsp", &nfsp, "nfsp/I"); eventVariables->Branch("px", px, "px[nfsp]/F"); eventVariables->Branch("py", py, "py[nfsp]/F"); eventVariables->Branch("pz", pz, "pz[nfsp]/F"); eventVariables->Branch("E", E, "E[nfsp]/F"); eventVariables->Branch("pdg", pdg, "pdg[nfsp]/I"); // Event Scaling Information - eventVariables->Branch("Weight", &Weight, "Weight/D"); - eventVariables->Branch("InputWeight", &InputWeight, "InputWeight/D"); - eventVariables->Branch("RWWeight", &RWWeight, "RWWeight/D"); + eventVariables->Branch("Weight", &Weight, "Weight/F"); + eventVariables->Branch("InputWeight", &InputWeight, "InputWeight/F"); + eventVariables->Branch("RWWeight", &RWWeight, "RWWeight/F"); + // Should be a double because may be 1E-39 and less eventVariables->Branch("fScaleFactor", &fScaleFactor, "fScaleFactor/D"); return; } void GenericFlux_Vectors::FillEventVariables(FitEvent *event) { // Reset all Function used to extract any variables of interest to the event Mode = PDGnu = tgt = PDGLep = 0; Enu_true = ELep = CosLep = Q2 = q0 = q3 = Enu_QE = Q2_QE = W_nuc_rest = W = x = y = -999.9; nfsp = 0; for (int i = 0; i < kMAX; ++i) { px[i] = py[i] = pz[i] = E[i] = -999; pdg[i] = 0; } Weight = InputWeight = RWWeight = 0; partList.clear(); // Now fill the information Mode = event->Mode; cc = (abs(event->Mode) < 30); // Get the incoming neutrino and outgoing lepton FitParticle *nu = event->GetNeutrinoIn(); FitParticle *lep = event->GetHMFSAnyLepton(); PDGnu = nu->fPID; Enu_true = nu->fP.E() / 1E3; tgt = event->fTargetPDG; if (lep != NULL) { PDGLep = lep->fPID; ELep = lep->fP.E() / 1E3; CosLep = cos(nu->fP.Vect().Angle(lep->fP.Vect())); // Basic interaction kinematics Q2 = -1 * (nu->fP - lep->fP).Mag2() / 1E6; q0 = (nu->fP - lep->fP).E() / 1E3; q3 = (nu->fP - lep->fP).Vect().Mag() / 1E3; // These assume C12 binding from MINERvA... not ideal Enu_QE = FitUtils::EnuQErec(lep->fP, CosLep, 34., true); Q2_QE = FitUtils::Q2QErec(lep->fP, CosLep, 34., true); // Get W_true with assumption of initial state nucleon at rest float m_n = (float)PhysConst::mass_proton; W_nuc_rest = sqrt(-Q2 + 2 * m_n * q0 + m_n * m_n); W = sqrt(-Q2 + 2 * m_n * q0 + m_n * m_n); x = Q2 / (2 * m_n * q0); y = 1 - ELep / Enu_true; } // Loop over the particles and store all the final state particles in a vector for (UInt_t i = 0; i < event->Npart(); ++i) { - bool part_alive = event->PartInfo(i)->fIsAlive and + bool part_alive = event->PartInfo(i)->fIsAlive && event->PartInfo(i)->Status() == kFinalState; - if (!part_alive) - continue; + if (!part_alive) continue; partList.push_back(event->PartInfo(i)); } // Save outgoing particle vectors nfsp = (int)partList.size(); for (int i = 0; i < nfsp; ++i) { px[i] = partList[i]->fP.X() / 1E3; py[i] = partList[i]->fP.Y() / 1E3; pz[i] = partList[i]->fP.Z() / 1E3; E[i] = partList[i]->fP.E() / 1E3; pdg[i] = partList[i]->fPID; } // Fill event weights Weight = event->RWWeight * event->InputWeight; RWWeight = event->RWWeight; InputWeight = event->InputWeight; // Fill the eventVariables Tree eventVariables->Fill(); return; }; +void GenericFlux_Vectors::AddSignalFlagsToTree() { + if (!eventVariables) { + Config::Get().out->cd(); + eventVariables = new TTree((this->fName + "_VARS").c_str(), + (this->fName + "_VARS").c_str()); + } + + LOG(SAM) << "Adding signal flags" << 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_Vectors::Write(std::string drawOpt) { // First save the TTree eventVariables->Write(); // Save Flux and Event Histograms too GetInput()->GetFluxHistogram()->Write(); GetInput()->GetEventHistogram()->Write(); return; } // Override functions which aren't really necessary bool GenericFlux_Vectors::isSignal(FitEvent *event) { (void)event; return true; }; void GenericFlux_Vectors::ScaleEvents() { return; } void GenericFlux_Vectors::ApplyNormScale(float norm) { this->fCurrentNorm = norm; return; } void GenericFlux_Vectors::FillHistograms() { return; } void GenericFlux_Vectors::ResetAll() { eventVariables->Reset(); return; } float GenericFlux_Vectors::GetChi2() { return 0.0; } diff --git a/src/MCStudies/GenericFlux_Vectors.h b/src/MCStudies/GenericFlux_Vectors.h index 237e467..ec618cd 100644 --- a/src/MCStudies/GenericFlux_Vectors.h +++ b/src/MCStudies/GenericFlux_Vectors.h @@ -1,98 +1,115 @@ // 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 . *******************************************************************************/ #ifndef GenericFlux_Vectors_H_SEEN #define GenericFlux_Vectors_H_SEEN #include "Measurement1D.h" class GenericFlux_Vectors : public Measurement1D { public: GenericFlux_Vectors(std::string name, std::string inputfile, FitWeight *rw, std::string type, std::string fakeDataFile); virtual ~GenericFlux_Vectors() {}; //! Grab info from event void FillEventVariables(FitEvent *event); //! Fill Custom Histograms void FillHistograms(); //! ResetAll void ResetAll(); //! Scale void ScaleEvents(); //! Norm void ApplyNormScale(float norm); //! Define this samples signal bool isSignal(FitEvent *nvect); //! Write Files void Write(std::string drawOpt); //! Get Chi2 float GetChi2(); void AddEventVariablesToTree(); + void AddSignalFlagsToTree(); private: TTree* eventVariables; std::vector partList; int Mode; bool cc; int PDGnu; int tgt; int PDGLep; float ELep; float CosLep; // Basic interaction kinematics float Q2; float q0; float q3; float Enu_QE; float Enu_true; float Q2_QE; float W_nuc_rest; float W; float x; float y; // Save outgoing particle vectors int nfsp; static const int kMAX = 200; float px[kMAX]; float py[kMAX]; float pz[kMAX]; float E[kMAX]; int pdg[kMAX]; // Basic event info - double Weight; - double InputWeight; - double RWWeight; + float Weight; + float InputWeight; + float RWWeight; double fScaleFactor; + // Generic signal flags + bool flagCCINC; + bool flagNCINC; + bool flagCCQE; + bool flagCC0pi; + bool flagCCQELike; + bool flagNCEL; + bool flagNC0pi; + bool flagCCcoh; + bool flagNCcoh; + bool flagCC1pip; + bool flagNC1pip; + bool flagCC1pim; + bool flagNC1pim; + bool flagCC1pi0; + bool flagNC1pi0; }; #endif