diff --git a/app/PrepareGENIE.cxx b/app/PrepareGENIE.cxx
index 3839119..0e74407 100644
--- a/app/PrepareGENIE.cxx
+++ b/app/PrepareGENIE.cxx
@@ -1,614 +1,615 @@
#include <stdio.h>
#include <stdlib.h>
#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;
// Setup TTree
TChain* tn = new TChain("gtree");
tn->AddFile(input.c_str());
int nevt = tn->GetEntries();
NtpMCEventRecord* genientpl = NULL;
tn->SetBranchAddress("gmcrec", &genientpl);
TH1D* fluxhist = new TH1D("flux", "flux", 1000, 0, 10);
fluxhist->Fill(MonoEnergy);
fluxhist->Scale(1, "width");
// Make Event Hist
TH1D* eventhist = (TH1D*)fluxhist->Clone();
eventhist->Reset();
TH1D* xsechist = (TH1D*)eventhist->Clone();
// Create maps
std::map<std::string, TH1D*> modexsec;
std::map<std::string, TH1D*> modecount;
std::vector<std::string> genieids;
std::vector<std::string> targetids;
std::vector<std::string> interids;
// Loop over all events
for (int i = 0; i < nevt; i++) {
tn->GetEntry(i);
StopTalking();
EventRecord& event = *(genientpl->event);
GHepParticle* neu = event.Probe();
StartTalking();
// Get XSec From Spline
GHepRecord genie_record = static_cast<GHepRecord>(event);
double xsec = (genie_record.XSec() / (1E-38 * genie::units::cm2));
// Parse Interaction String
std::string mode = genie_record.Summary()->AsString();
std::vector<std::string> modevec = GeneralUtils::ParseToStr(mode, ";");
std::string targ = (modevec[0] + ";" + modevec[1]);
std::string inter = mode;
// Fill lists of Unique IDS
if (std::find(targetids.begin(), targetids.end(), targ) ==
targetids.end()) {
targetids.push_back(targ);
}
if (std::find(interids.begin(), interids.end(), inter) == interids.end()) {
interids.push_back(inter);
}
// Create entries Mode Maps
if (modexsec.find(mode) == modexsec.end()) {
genieids.push_back(mode);
modexsec[mode] = (TH1D*)xsechist->Clone();
modecount[mode] = (TH1D*)xsechist->Clone();
}
// Fill XSec Histograms
modexsec[mode]->Fill(neu->E(), xsec);
modecount[mode]->Fill(neu->E());
// Fill total event hist
eventhist->Fill(neu->E());
// Clear Event
genientpl->Clear();
- if (i % (nevt / 20) == 0) {
+ 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();
LOG(FIT) << "Getting splines in mono" << std::endl;
// Save each of the reconstructed splines to file
std::map<std::string, TH1D*> modeavg;
TDirectory* inddir = (TDirectory*)outputfile->Get("IndividualGENIESplines");
if (!inddir)
inddir = (TDirectory*)outputfile->mkdir("IndividualGENIESplines");
inddir->cd();
// Loop over GENIE ID's and get MEC count
int MECcount = 0;
bool MECcorrect = FitPar::Config().GetParB("CorrectGENIEMECNorm");
for (UInt_t i = 0; i < genieids.size(); i++) {
if (genieids[i].find("MEC") != std::string::npos) {
MECcount++;
}
}
LOG(FIT) << "Found " << MECcount << " repeated MEC instances." << std::endl;
for (UInt_t i = 0; i < genieids.size(); i++) {
std::string mode = genieids[i];
modexsec[mode]->Write((mode + "_summed_xsec").c_str(), TObject::kOverwrite);
modecount[mode]->Write((mode + "_summed_evt").c_str(), TObject::kOverwrite);
// Form extra avg xsec map -> Reconstructed spline
modeavg[mode] = (TH1D*)modexsec[mode]->Clone();
modeavg[mode]->Divide(modecount[mode]);
if (MECcorrect && (mode.find("MEC") != std::string::npos)) {
modeavg[mode]->Scale(1.0 / double(MECcount));
}
modeavg[mode]->Write((mode + "_rec_spline").c_str(), TObject::kOverwrite);
}
TDirectory* targdir = (TDirectory*)outputfile->Get("TargetGENIESplines");
if (!targdir) targdir = (TDirectory*)outputfile->mkdir("TargetGENIESplines");
targdir->cd();
LOG(FIT) << "Getting Target Splines" << std::endl;
// For each target save a total spline
std::map<std::string, TH1D*> targetsplines;
for (uint i = 0; i < targetids.size(); i++) {
LOG(FIT) << "Getting target " << i << std::endl;
std::string targ = targetids[i];
targetsplines[targ] = (TH1D*)xsechist->Clone();
LOG(FIT) << "Created target spline for " << targ << std::endl;
for (uint j = 0; j < genieids.size(); j++) {
std::string mode = genieids[j];
if (mode.find(targ) != std::string::npos) {
LOG(FIT) << "Mode " << mode << " contains " << targ << " target!"
<< std::endl;
targetsplines[targ]->Add(modeavg[mode]);
LOG(FIT) << "Finished with Mode " << mode << " "
<< modeavg[mode]->Integral() << std::endl;
}
}
LOG(FIT) << "Saving target spline:" << targ << std::endl;
targetsplines[targ]->Write(("Total" + targ).c_str(), TObject::kOverwrite);
}
LOG(FIT) << "Getting total splines" << std::endl;
// Now we have each of the targets we need to create a total cross-section.
int totalnucl = 0;
std::vector<std::string> targprs = GeneralUtils::ParseToStr(target, ",");
TH1D* totalxsec = (TH1D*)xsechist->Clone();
for (uint i = 0; i < targprs.size(); i++) {
std::string targpdg = targprs[i];
for (std::map<std::string, TH1D*>::iterator iter = targetsplines.begin();
iter != targetsplines.end(); iter++) {
std::string targstr = iter->first;
TH1D* xsec = iter->second;
if (targstr.find(targpdg) != std::string::npos) {
LOG(FIT) << "Adding target spline " << targstr
<< " Integral = " << xsec->Integral("width") << std::endl;
totalxsec->Add(xsec);
int nucl = atoi(targpdg.c_str());
totalnucl += int((nucl % 10000) / 10);
}
}
}
outputfile->cd();
totalxsec->Write("nuisance_Xsec", TObject::kOverwrite);
eventhist = (TH1D*)totalxsec->Clone();
eventhist->Multiply(fluxhist);
eventhist->Write("nuisance_events", TObject::kOverwrite);
fluxhist->Write("nuisance_flux", TObject::kOverwrite);
LOG(FIT) << "Inclusive XSec Per Nucleon = "
<< eventhist->Integral("width") * 1E-38 / fluxhist->Integral("width")
<< std::endl;
std::cout << "XSec Hist Integral = " << xsechist->Integral("width")
<< std::endl;
return;
}
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;
// Get Flux Hist
std::vector<std::string> 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<TH1*>(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<std::string> 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<std::string, TH1D*> modexsec;
std::map<std::string, TH1D*> modecount;
std::vector<std::string> genieids;
std::vector<std::string> targetids;
std::vector<std::string> interids;
// Loop over all events
for (int i = 0; i < nevt; i++) {
tn->GetEntry(i);
StopTalking();
EventRecord& event = *(genientpl->event);
GHepParticle* neu = event.Probe();
StartTalking();
// Get XSec From Spline
GHepRecord genie_record = static_cast<GHepRecord>(event);
double xsec = (genie_record.XSec() / (1E-38 * genie::units::cm2));
// Parse Interaction String
std::string mode = genie_record.Summary()->AsString();
std::vector<std::string> modevec = GeneralUtils::ParseToStr(mode, ";");
std::string targ = (modevec[0] + ";" + modevec[1]);
std::string inter = mode;
// Fill lists of Unique IDS
if (std::find(targetids.begin(), targetids.end(), targ) ==
targetids.end()) {
targetids.push_back(targ);
}
if (std::find(interids.begin(), interids.end(), inter) == interids.end()) {
interids.push_back(inter);
}
// Create entries Mode Maps
if (modexsec.find(mode) == modexsec.end()) {
genieids.push_back(mode);
modexsec[mode] = (TH1D*)xsechist->Clone();
modecount[mode] = (TH1D*)xsechist->Clone();
}
// Fill XSec Histograms
modexsec[mode]->Fill(neu->E(), xsec);
modecount[mode]->Fill(neu->E());
// Fill total event hist
eventhist->Fill(neu->E());
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<std::string, TH1D*> modeavg;
TDirectory* inddir = (TDirectory*)outputfile->Get("IndividualGENIESplines");
if (!inddir)
inddir = (TDirectory*)outputfile->mkdir("IndividualGENIESplines");
inddir->cd();
// Loop over GENIE ID's and get MEC count
int MECcount = 0;
bool MECcorrect = FitPar::Config().GetParB("CorrectGENIEMECNorm");
for (UInt_t i = 0; i < genieids.size(); i++) {
if (genieids[i].find("MEC") != std::string::npos) {
MECcount++;
}
}
LOG(FIT) << "Found " << MECcount << " repeated MEC instances." << std::endl;
for (UInt_t i = 0; i < genieids.size(); i++) {
std::string mode = genieids[i];
modexsec[mode]->Write((mode + "_summed_xsec").c_str(), TObject::kOverwrite);
modecount[mode]->Write((mode + "_summed_evt").c_str(), TObject::kOverwrite);
// Form extra avg xsec map -> Reconstructed spline
modeavg[mode] = (TH1D*)modexsec[mode]->Clone();
modeavg[mode]->Divide(modecount[mode]);
if (MECcorrect && (mode.find("MEC") != std::string::npos)) {
modeavg[mode]->Scale(1.0 / double(MECcount));
}
modeavg[mode]->Write((mode + "_rec_spline").c_str(), TObject::kOverwrite);
}
TDirectory* targdir = (TDirectory*)outputfile->Get("TargetGENIESplines");
if (!targdir) targdir = (TDirectory*)outputfile->mkdir("TargetGENIESplines");
targdir->cd();
LOG(FIT) << "Getting Target Splines" << std::endl;
// For each target save a total spline
std::map<std::string, TH1D*> targetsplines;
for (uint i = 0; i < targetids.size(); i++) {
LOG(FIT) << "Getting target " << i << std::endl;
std::string targ = targetids[i];
targetsplines[targ] = (TH1D*)xsechist->Clone();
LOG(FIT) << "Created target spline for " << targ << std::endl;
for (uint j = 0; j < genieids.size(); j++) {
std::string mode = genieids[j];
// Look at all matching modes/targets
if (mode.find(targ) != std::string::npos) {
LOG(FIT) << "Mode " << mode << " contains " << targ << " target!"
<< std::endl;
// modeavg[mode]->Write( (mode + "_cont_" + targ).c_str() ,
// TObject::kOverwrite);
targetsplines[targ]->Add(modeavg[mode]);
LOG(FIT) << "Finished with Mode " << mode << " "
<< modeavg[mode]->Integral() << std::endl;
}
}
LOG(FIT) << "Saving target spline:" << targ << std::endl;
targetsplines[targ]->Write(("Total" + targ).c_str(), TObject::kOverwrite);
}
LOG(FIT) << "Getting total splines" << std::endl;
// Now we have each of the targets we need to create a total cross-section.
int totalnucl = 0;
std::vector<std::string> targprs = GeneralUtils::ParseToStr(target, ",");
TH1D* totalxsec = (TH1D*)xsechist->Clone();
for (uint i = 0; i < targprs.size(); i++) {
std::string targpdg = targprs[i];
for (std::map<std::string, TH1D*>::iterator iter = targetsplines.begin();
iter != targetsplines.end(); iter++) {
std::string targstr = iter->first;
TH1D* xsec = iter->second;
if (targstr.find(targpdg) != std::string::npos) {
LOG(FIT) << "Adding target spline " << targstr
<< " Integral = " << xsec->Integral("width") << std::endl;
totalxsec->Add(xsec);
int nucl = atoi(targpdg.c_str());
totalnucl += int((nucl % 10000) / 10);
}
}
}
LOG(FIT) << "Total XSec Integral = " << totalxsec->Integral("width")
<< std::endl;
outputfile->cd();
totalxsec->Write("nuisance_xsec", TObject::kOverwrite);
eventhist = (TH1D*)fluxhist->Clone();
eventhist->Multiply(totalxsec);
LOG(FIT) << "Dividing by Total Nucl = " << totalnucl << std::endl;
eventhist->Scale(1.0 / double(totalnucl));
eventhist->Write("nuisance_events", TObject::kOverwrite);
fluxhist->Write("nuisance_flux", TObject::kOverwrite);
LOG(FIT) << "Inclusive XSec Per Nucleon = "
<< eventhist->Integral("width") * 1E-38 / fluxhist->Integral("width")
<< std::endl;
std::cout << "XSec Hist Integral = " << xsechist->Integral("width")
<< std::endl;
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/cmake/GENIESetup.cmake b/cmake/GENIESetup.cmake
index 7c0b2c5..b4d12f0 100644
--- a/cmake/GENIESetup.cmake
+++ b/cmake/GENIESetup.cmake
@@ -1,157 +1,160 @@
# Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
################################################################################
# This file is part of NUISANCE.
#
# NUISANCE is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# NUISANCE is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with NUISANCE. If not, see <http://www.gnu.org/licenses/>.
################################################################################
# TODO
# check system for libxml2
# check whether we need the includes
# check if we can use a subset of the GENIE libraries
################################################################################
# Check Dependencies
################################################################################
################################# GENIE ######################################
if(GENIE STREQUAL "")
cmessage(FATAL_ERROR "Variable GENIE is not defined. "
"The location of a pre-built GENIE install must be defined either as"
" $ cmake -DGENIE=/path/to/GENIE or as and environment vairable"
" $ export GENIE=/path/to/GENIE")
endif()
if (BUILD_GEVGEN)
cmessage(STATUS "Building custom gevgen")
LIST(APPEND EXTRA_CXX_FLAGS -D__GEVGEN_ENABLED__)
endif()
# Extract GENIE VERSION
if (GENIE_VERSION STREQUAL "AUTO")
execute_process (COMMAND ${CMAKE_SOURCE_DIR}/cmake/getgenieversion.sh ${GENIE}
OUTPUT_VARIABLE GENIE_VERSION OUTPUT_STRIP_TRAILING_WHITESPACE)
endif()
execute_process (COMMAND genie-config
--libs OUTPUT_VARIABLE GENIE_LD_FLAGS_STR OUTPUT_STRIP_TRAILING_WHITESPACE)
execute_process (COMMAND genie-config
--topsrcdir OUTPUT_VARIABLE GENIE_INCLUDES_DIR OUTPUT_STRIP_TRAILING_WHITESPACE)
string(REGEX MATCH "-L\([^ ]+\) \(.*\)$" PARSE_GENIE_LIBS_MATCH ${GENIE_LD_FLAGS_STR})
cmessage(DEBUG "genie-config --libs: ${GENIE_LD_FLAGS_STR}")
if(NOT PARSE_GENIE_LIBS_MATCH)
cmessage(FATAL_ERROR "Expected to be able to parse the result of genie-config --libs to a lib directory and a list of libraries to include, but got: \"${GENIE_LD_FLAGS_STR}\"")
endif()
set(GENIE_LIB_DIR ${CMAKE_MATCH_1})
set(GENIE_LIBS_RAW ${CMAKE_MATCH_2})
string(REPLACE "-l" "" GENIE_LIBS_STRIPED "${GENIE_LIBS_RAW}")
cmessage(STATUS "GENIE version : ${GENIE_VERSION}")
cmessage(STATUS "GENIE libdir : ${GENIE_LIB_DIR}")
cmessage(STATUS "GENIE libs : ${GENIE_LIBS_STRIPED}")
string(REGEX MATCH "ReinSeghal" WASMATCHED ${GENIE_LIBS_STRIPED})
if(WASMATCHED AND GENIE_VERSION STREQUAL "210")
set(GENIE_SEHGAL ${GENIE_LIBS_STRIPED})
STRING(REPLACE "ReinSeghal" "ReinSehgal" GENIE_LIBS_STRIPED ${GENIE_SEHGAL})
cmessage(DEBUG "Fixed inconsistency in library naming: ${GENIE_LIBS_STRIPED}")
endif()
string(REGEX MATCH "ReWeight" WASMATCHED ${GENIE_LIBS_STRIPED})
if(NOT WASMATCHED)
set(GENIE_LIBS_STRIPED "GReWeight ${GENIE_LIBS_STRIPED}")
cmessage(DEBUG "Force added ReWeight library: ${GENIE_LIBS_STRIPED}")
endif()
string(REPLACE " " ";" GENIE_LIBS_LIST "${GENIE_LIBS_STRIPED}")
cmessage(DEBUG "genie-config --libs -- MATCH1: ${CMAKE_MATCH_1}")
cmessage(DEBUG "genie-config --libs -- MATCH2: ${CMAKE_MATCH_2}")
cmessage(DEBUG "genie-config --libs -- libs stripped: ${GENIE_LIBS_STRIPED}")
cmessage(DEBUG "genie-config --libs -- libs list: ${GENIE_LIBS_LIST}")
################################ LHAPDF ######################################
if(LHAPDF_LIB STREQUAL "")
cmessage(FATAL_ERROR "Variable LHAPDF_LIB is not defined. The location of a pre-built lhapdf install must be defined either as $ cmake -DLHAPDF_LIB=/path/to/LHAPDF_libraries or as and environment vairable $ export LHAPDF_LIB=/path/to/LHAPDF_libraries")
endif()
if(LHAPDF_INC STREQUAL "")
cmessage(FATAL_ERROR "Variable LHAPDF_INC is not defined. The location of a pre-built lhapdf install must be defined either as $ cmake -DLHAPDF_INC=/path/to/LHAPDF_includes or as and environment vairable $ export LHAPDF_INC=/path/to/LHAPDF_includes")
endif()
if(LHAPATH STREQUAL "")
cmessage(FATAL_ERROR "Variable LHAPATH is not defined. The location of a the LHAPATH directory must be defined either as $ cmake -DLHAPATH=/path/to/LHAPATH or as and environment variable $ export LHAPATH=/path/to/LHAPATH")
endif()
################################ LIBXML ######################################
if(LIBXML2_LIB STREQUAL "")
cmessage(FATAL_ERROR "Variable LIBXML2_LIB is not defined. The location of a pre-built libxml2 install must be defined either as $ cmake -DLIBXML2_LIB=/path/to/LIBXML2_libraries or as and environment vairable $ export LIBXML2_LIB=/path/to/LIBXML2_libraries")
endif()
if(LIBXML2_INC STREQUAL "")
cmessage(FATAL_ERROR "Variable LIBXML2_INC is not defined. The location of a pre-built libxml2 install must be defined either as $ cmake -DLIBXML2_INC=/path/to/LIBXML2_includes or as and environment vairable $ export LIBXML2_INC=/path/to/LIBXML2_includes")
endif()
############################### log4cpp ######################################
if(LOG4CPP_LIB STREQUAL "")
cmessage(FATAL_ERROR "Variable LOG4CPP_LIB is not defined. The location of a pre-built log4cpp install must be defined either as $ cmake -DLOG4CPP_LIB=/path/to/LOG4CPP_libraries or as and environment vairable $ export LOG4CPP_LIB=/path/to/LOG4CPP_libraries")
endif()
if(LOG4CPP_INC STREQUAL "")
cmessage(FATAL_ERROR "Variable LOG4CPP_INC is not defined. The location of a pre-built log4cpp install must be defined either as $ cmake -DGENIE_LOG4CPP_INC=/path/to/LOG4CPP_includes or as and environment vairable $ export LOG4CPP_INC=/path/to/LOG4CPP_includes")
endif()
################################################################################
LIST(APPEND EXTRA_CXX_FLAGS -D__GENIE_ENABLED__ -D__GENIE_VERSION__=${GENIE_VERSION})
LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES
${GENIE_INCLUDES_DIR}
${GENIE_INCLUDES_DIR}/GHEP
${GENIE_INCLUDES_DIR}/Ntuple
${GENIE_INCLUDES_DIR}/ReWeight
${GENIE_INCLUDES_DIR}/Apps
${GENIE_INCLUDES_DIR}/FluxDrivers
${GENIE_INCLUDES_DIR}/EVGDrivers
${LHAPDF_INC}
${LIBXML2_INC}
${LOG4CPP_INC})
SAYVARS()
LIST(APPEND EXTRA_LINK_DIRS
${GENIE_LIB_DIR}
${LHAPDF_LIB}
${LIBXML2_LIB}
${LOG4CPP_LIB})
#LIST(REVERSE EXTRA_LIBS)
#LIST(REVERSE GENIE_LIBS_LIST)
+
+LIST(APPEND EXTRA_LIBS -Wl,--start-group)
LIST(APPEND EXTRA_LIBS ${GENIE_LIBS_LIST})
+LIST(APPEND EXTRA_LIBS -Wl,--end-group)
#LIST(REVERSE EXTRA_LIBS)
LIST(APPEND EXTRA_LIBS LHAPDF xml2 log4cpp)
if(USE_PYTHIA8)
set(NEED_PYTHIA8 TRUE)
set(NEED_ROOTPYTHIA8 TRUE)
else()
set(NEED_PYTHIA6 TRUE)
set(NEED_ROOTPYTHIA6 TRUE)
endif()
set(NEED_ROOTEVEGEN TRUE)
SET(USE_GENIE TRUE CACHE BOOL "Whether to enable GENIE (reweight) support. Requires external libraries. <FALSE>" FORCE)
diff --git a/cmake/c++CompilerSetup.cmake b/cmake/c++CompilerSetup.cmake
index a59c304..bd9adcb 100644
--- a/cmake/c++CompilerSetup.cmake
+++ b/cmake/c++CompilerSetup.cmake
@@ -1,126 +1,130 @@
# 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 <http://www.gnu.org/licenses/>.
################################################################################
if(USE_OMP)
LIST(APPEND EXTRA_CXX_FLAGS -fopenmp)
endif()
if(USE_DYNSAMPLES)
LIST(APPEND EXTRA_LIBS dl)
LIST(APPEND EXTRA_CXX_FLAGS -D__USE_DYNSAMPLES__)
endif()
set(CXX_WARNINGS -Wall )
cmessage(DEBUG "EXTRA_CXX_FLAGS: ${EXTRA_CXX_FLAGS}")
string(REPLACE ";" " " STR_EXTRA_CXX_FLAGS "${EXTRA_CXX_FLAGS}")
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${STR_EXTRA_CXX_FLAGS} ${CXX_WARNINGS}")
set(CMAKE_Fortran_FLAGS_RELEASE "-fPIC")
set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -O0")
if(USE_DYNSAMPLES)
set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -fPIC")
set(CMAKE_Fortran_FLAGS_DEBUG "-fPIC")
endif()
set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -fPIC -O3")
if(CMAKE_BUILD_TYPE MATCHES DEBUG)
set(CURRENT_CMAKE_CXX_FLAGS ${CMAKE_CXX_FLAGS_DEBUG})
elseif(CMAKE_BUILD_TYPE MATCHES RELEASE)
set(CURRENT_CMAKE_CXX_FLAGS ${CMAKE_CXX_FLAGS_RELEASE})
else()
cmessage(FATAL_ERROR "[ERROR]: Unknown CMAKE_BUILD_TYPE (\"${CMAKE_BUILD_TYPE}\"): Should be \"DEBUG\" or \"RELEASE\".")
endif()
SET(STR_EXTRA_LINK_DIRS)
if(NOT EXTRA_LINK_DIRS STREQUAL "")
string(REPLACE ";" " -L" STR_EXTRA_LINK_DIRS "-L${EXTRA_LINK_DIRS}")
endif()
+
SET(STR_EXTRA_LIBS)
if(NOT EXTRA_LIBS STREQUAL "")
- string(REPLACE ";" " -l" STR_EXTRA_LIBS "-l${EXTRA_LIBS}")
+ SET(STR_EXTRA_LIBS_NO_SCRUB_LINKOPTS)
+ string(REPLACE ";" " -l" STR_EXTRA_LIBS_NO_SCRUB_LINKOPTS "-l${EXTRA_LIBS}")
+ string(REPLACE "-l-" "-" STR_EXTRA_LIBS ${STR_EXTRA_LIBS_NO_SCRUB_LINKOPTS})
endif()
+
SET(STR_EXTRA_SHAREDOBJS)
if(NOT EXTRA_SHAREDOBJS STREQUAL "")
string(REPLACE ";" " " STR_EXTRA_SHAREDOBJS "${EXTRA_SHAREDOBJS}")
endif()
SET(STR_EXTRA_LINK_FLAGS)
if(NOT EXTRA_LINK_FLAGS STREQUAL "")
string(REPLACE ";" " " STR_EXTRA_LINK_FLAGS "${EXTRA_LINK_FLAGS}")
endif()
cmessage(DEBUG "EXTRA_LINK_DIRS: ${STR_EXTRA_LINK_DIRS}")
cmessage(DEBUG "EXTRA_LIBS: ${STR_EXTRA_LIBS}")
cmessage(DEBUG "EXTRA_SHAREDOBJS: ${STR_EXTRA_SHAREDOBJS}")
cmessage(DEBUG "EXTRA_LINK_FLAGS: ${STR_EXTRA_LINK_FLAGS}")
if(NOT STR_EXTRA_LINK_DIRS STREQUAL "" AND NOT STR_EXTRA_LIBS STREQUAL "")
SET(CMAKE_DEPENDLIB_FLAGS "${STR_EXTRA_LINK_DIRS} ${STR_EXTRA_LIBS}")
endif()
if(USE_NEUT)
foreach(OBJ ${NEUT_ROOT_LIBS})
if(NOT CMAKE_DEPENDLIB_FLAGS STREQUAL "")
SET(CMAKE_DEPENDLIB_FLAGS "${CMAKE_DEPENDLIB_FLAGS} ${OBJ}")
else()
SET(CMAKE_DEPENDLIB_FLAGS "${OBJ}")
endif()
endforeach()
foreach(OBJ ${NEUT_ROOT_LIBS})
if(NOT CMAKE_DEPENDLIB_FLAGS STREQUAL "")
SET(CMAKE_DEPENDLIB_FLAGS "${CMAKE_DEPENDLIB_FLAGS} ${OBJ}")
else()
SET(CMAKE_DEPENDLIB_FLAGS "${OBJ}")
endif()
endforeach()
endif()
if(NOT EXTRA_SHAREDOBJS STREQUAL "")
if(NOT STR_EXTRA_LINK_FLAGS STREQUAL "")
SET(STR_EXTRA_LINK_FLAGS "${STR_EXTRA_SHAREDOBJS} ${STR_EXTRA_LINK_FLAGS}")
else()
SET(STR_EXTRA_LINK_FLAGS "${STR_EXTRA_SHAREDOBJS}")
endif()
endif()
if(NOT EXTRA_LINK_FLAGS STREQUAL "")
if(NOT CMAKE_LINK_FLAGS STREQUAL "")
SET(CMAKE_LINK_FLAGS "${CMAKE_LINK_FLAGS} ${STR_EXTRA_LINK_FLAGS}")
else()
SET(CMAKE_LINK_FLAGS "${STR_EXTRA_LINK_FLAGS}")
endif()
endif()
if(USE_OMP)
cmessage(FATAL_ERROR "No OMP features currently enabled so this is a FATAL_ERROR to let you know that you don't gain anything with this declaration.")
endif()
if (VERBOSE)
cmessage (STATUS "C++ Compiler : ${CXX_COMPILER_NAME}")
cmessage (STATUS " flags : ${CMAKE_CXX_FLAGS}")
cmessage (STATUS " Release flags : ${CMAKE_CXX_FLAGS_RELEASE}")
cmessage (STATUS " Debug flags : ${CMAKE_CXX_FLAGS_DEBUG}")
cmessage (STATUS " Link Flags : ${CMAKE_LINK_FLAGS}")
cmessage (STATUS " Lib Flags : ${CMAKE_DEPENDLIB_FLAGS}")
endif()
diff --git a/src/InputHandler/InputHandler.cxx b/src/InputHandler/InputHandler.cxx
index b400590..80e9f65 100644
--- a/src/InputHandler/InputHandler.cxx
+++ b/src/InputHandler/InputHandler.cxx
@@ -1,267 +1,301 @@
// 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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
#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);
+ Int_t minBin = fEventHist->GetXaxis()->FindFixBin(low);
+ Int_t maxBin = fEventHist->GetXaxis()->FindFixBin(high);
- return fEventHist->Integral(minBin, maxBin + 1, intOpt.c_str());
+ if ((fEventHist->IsBinOverflow(minBin) && (low != -9999.9))) {
+ minBin = 1;
+ }
+
+ if ((fEventHist->IsBinOverflow(maxBin) && (high != -9999.9))) {
+ maxBin = fEventHist->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) / fEventHist->GetXaxis()->GetBinWidth(minBin)) *
+ fEventHist->Integral(minBin, minBin, intOpt.c_str());
+ }
+
+ double lowBinUpEdge = fEventHist->GetXaxis()->GetBinUpEdge(minBin);
+ double highBinLowEdge = fEventHist->GetXaxis()->GetBinLowEdge(maxBin);
+
+ double lowBinfracIntegral =
+ ((lowBinUpEdge - low) / fEventHist->GetXaxis()->GetBinWidth(minBin)) *
+ fEventHist->Integral(minBin, minBin, intOpt.c_str());
+ double highBinfracIntegral =
+ ((high - highBinLowEdge) / fEventHist->GetXaxis()->GetBinWidth(maxBin)) *
+ fEventHist->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 =
+ fEventHist->Integral(minBin + 1, maxBin - 1, intOpt.c_str());
+ // If there are filled bins between them
+ return lowBinfracIntegral + highBinfracIntegral + ContainedIntegral;
};
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 + ContainedIntegral;
- // return fFluxHist->Integral(minBin + 1, maxBin - 1, intOpt.c_str());
}
std::vector<TH1*> InputHandlerBase::GetFluxList(void) {
return std::vector<TH1*>(1, fFluxHist);
};
std::vector<TH1*> InputHandlerBase::GetEventList(void) {
return std::vector<TH1*>(1, fEventHist);
};
std::vector<TH1*> InputHandlerBase::GetXSecList(void) {
return std::vector<TH1*>(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<BaseFitEvt*>(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/InputHandler/InputHandler.h b/src/InputHandler/InputHandler.h
index 3377403..1110091 100644
--- a/src/InputHandler/InputHandler.h
+++ b/src/InputHandler/InputHandler.h
@@ -1,145 +1,145 @@
// 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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
#ifndef INPUTHANDLER2_H
#define INPUTHANDLER2_H
/*!
* \addtogroup InputHandler
* @{
*/
#include "TH1D.h"
#include "FitEvent.h"
#include "BaseFitEvt.h"
#include "TTreePerfStats.h"
/// Base InputHandler class defining how events are requested and setup.
class InputHandlerBase {
public:
/// Base constructor resets everything to default
InputHandlerBase();
/// Removes flux/event rate histograms
virtual ~InputHandlerBase();
/// Return NUISANCE FitEvent Class from given event entry.
/// Must be overriden by GeneratorInputHandler. Lightweight allows a faster option
/// to be given where only RW information is needed.
virtual FitEvent* GetNuisanceEvent(const UInt_t entry, const bool lightweight=false) = 0;
/// Calls GetNuisanceEvent(entry, TRUE);
virtual BaseFitEvt* GetBaseEvent(const UInt_t entry);
/// Print current event information
virtual void Print();
/// Return handler ID
inline std::string GetName (void) {return fName; };
/// Return Handler Event Type Index
inline int GetType (void) {return fEventType;};
/// Get Total Number of Events being Handled
inline virtual int GetNEvents (void) {return fNEvents; };
/// Get the Total Flux Histogram these events were generated with
inline virtual TH1D* GetFluxHistogram (void) {return fFluxHist; };
/// Get the Total Event Histogram these events were generated with
inline virtual TH1D* GetEventHistogram (void) {return fEventHist;};
/// Get the Total Cross-section Histogram (EventHist/FluxHist)
virtual TH1D* GetXSecHistogram(void);
/// Return all Flux Histograms for all InputFiles.
virtual std::vector<TH1*> GetFluxList(void);
/// Return all Event Histograms for all InputFiles
virtual std::vector<TH1*> GetEventList(void);
/// Return all Xsec Histograms for all InputFiles
virtual std::vector<TH1*> GetXSecList(void);
/// Placeholder to create a cache to speed up reads in GeneratorInputHandler
inline virtual void CreateCache(){};
/// Placeholder to remove optional cache to free up memory
inline virtual void RemoveCache(){};
/// Return starting NUISANCE event pointer (entry=0)
FitEvent* FirstNuisanceEvent();
/// Iterate to next NUISANCE event. Returns NULL when entry > fNEvents.
FitEvent* NextNuisanceEvent();
/// Returns starting Base Event Pointer (entry=0)
BaseFitEvt* FirstBaseEvent();
/// Iterate to next NUISANCE Base Event. Returns NULL when entry > fNEvents.
BaseFitEvt* NextBaseEvent();
/// Register an input file and update event/flux information
virtual void RegisterJointInput(std::string input, int n, TH1D* f, TH1D* e);
/// Finalise setup of Input event/flux information and calculate
/// joint input weights if joint input is provided.
virtual void SetupJointInputs();
/// Calculate a weight for the event given the joint input information.
/// Used to scale the relative proportion of multiple inputs correctly
/// with respect to one another.
virtual double GetInputWeight(int entry);
/// Returns the total predicted event rate for this input given the
/// low and high energy ranges. intOpt specifies the option the ROOT
/// TH1D integral should use. e.g. "" or "width"
- double PredictedEventRate(double low, double high,
- std::string intOpt);
+ double PredictedEventRate(double low = -9999.9, double high = -9999.9,
+ std::string intOpt = "");
/// Returns the total generated flux for this input given the
/// low and high energy ranges. intOpt specifies the option the ROOT
/// TH1D integral should use. e.g. "" or "width"
double TotalIntegratedFlux(double low = -9999.9, double high = -9999.9,
std::string intOpt = "");
/// Actual data members.
std::vector<TH1D*> jointfluxinputs;
std::vector<TH1D*> jointeventinputs;
std::vector<int> jointindexlow;
std::vector<int> jointindexhigh;
std::vector<int> jointindexallowed;
size_t jointindexswitch;
bool jointinput;
std::vector<double> jointindexscale;
std::string fName;
TH1D* fFluxHist;
TH1D* fEventHist;
TH1D* fXSecHist;
int fNEvents;
int fMaxEvents;
FitEvent* fNUISANCEEvent;
BaseFitEvt* fBaseEvent;
int fEventType;
int fCurrentIndex;
int fCacheSize;
bool kRemoveUndefParticles;
bool kRemoveFSIParticles;
bool kRemoveNuclearParticles;
TTreePerfStats* fTTreePerformance;
};
/*! @} */
#endif
diff --git a/src/MCStudies/GenericFlux_Tester.cxx b/src/MCStudies/GenericFlux_Tester.cxx
index 424a3c1..1bd7f2b 100644
--- a/src/MCStudies/GenericFlux_Tester.cxx
+++ b/src/MCStudies/GenericFlux_Tester.cxx
@@ -1,582 +1,590 @@
// 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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
#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");
+ 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.
this->fScaleFactor =
- (GetEventHistogram()->Integral("width") * 1E-38 / (fNEvents + 0.)) /
+ (this->PredictedEventRate("width") * 1E-38 / (fNEvents + 0.)) /
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);
}
// 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;
// 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;
// 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);
}
// 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;
}