diff --git a/app/PrepareGENIE.cxx b/app/PrepareGENIE.cxx
index b67d856..8af366c 100644
--- a/app/PrepareGENIE.cxx
+++ b/app/PrepareGENIE.cxx
@@ -1,942 +1,942 @@
#include "FitLogger.h"
#include "PlotUtils.h"
#include "TFile.h"
#include "TH1D.h"
#include "TTree.h"
#include <stdio.h>
#include <stdlib.h>
#ifdef __GENIE_ENABLED__
#ifdef GENIE_PRE_R3
#include "Conventions/Units.h"
#include "GHEP/GHepParticle.h"
#include "PDG/PDGUtils.h"
#else
#include "Framework/Conventions/Units.h"
#include "Framework/GHEP/GHepParticle.h"
#include "Framework/ParticleData/PDGUtils.h"
#endif
#endif
std::string gInputFiles = "";
std::string gOutputFile = "";
std::string gFluxFile = "";
std::string gTarget = "";
double MonoEnergy;
int gNEvents = -999;
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) {
LOG(FIT) << "Running GENIE Prepare in mono energetic with E = " << MonoEnergy
<< " GeV" << std::endl;
// Setup TTree
TChain *tn = new TChain("gtree");
tn->AddFile(input.c_str());
if (tn->GetFile() == NULL) {
tn->Print();
ERROR(FTL, "gtree not located in GENIE file: " << input);
THROW("Check your inputs, they may need to be completely regenerated!");
throw;
}
int nevt = tn->GetEntries();
if (gNEvents != -999) {
LOG(FIT) << "Overriding number of events by user from " << nevt << " to "
<< gNEvents << std::endl;
nevt = gNEvents;
}
NtpMCEventRecord *genientpl = NULL;
tn->SetBranchAddress("gmcrec", &genientpl);
// 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<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();
modexsec[mode]->GetYaxis()->SetTitle(
"d#sigma/dE_{#nu} #times 10^{-38} (events weighted by #sigma)");
modecount[mode]->GetYaxis()->SetTitle("Number of events in file");
}
// 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 (E: " << neu->E() << " GeV, xsec: " << xsec
<< " E-38 cm^2/nucleon)" << std::endl;
}
}
LOG(FIT) << "Processed all events" << std::endl;
TFile *outputfile;
// If no output is specified just append to the file
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();
+ TTree *cloneTree = tn->CloneTree(-1, "fast");
cloneTree->SetDirectory(outputfile);
cloneTree->Write();
QLOG(FIT, "Cloning input nova_wgts to output file: " << gOutputFile);
// ***********************************
// ***********************************
// FUDGE FOR NOVA MINERVA WORKSHOP
// Also check for the nova_wgts tree from Jeremy
TChain *nova_chain = new TChain("nova_wgts");
nova_chain->AddFile(input.c_str());
TTree *nova_tree = nova_chain->GetTree();
if (!nova_tree) {
QLOG(FIT, "Could not find nova_wgts tree in " << gOutputFile);
} else {
QLOG(FIT, "Found nova_wgts tree in " << gOutputFile);
}
if (nova_tree) {
nova_tree->SetDirectory(outputfile);
nova_tree->Write();
}
QLOG(FIT, "Done cloning tree.");
}
LOG(FIT) << "Getting splines in mono-energetic..." << 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]->GetYaxis()->SetTitle(
"#sigma (E_{#nu}) #times 10^{-38} (cm^{2}/target)");
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++) {
std::string targ = targetids[i];
LOG(FIT) << "Getting target " << i << ": " << targ << std::endl;
targetsplines[targ] = (TH1D *)xsechist->Clone();
targetsplines[targ]->GetYaxis()->SetTitle(
"#sigma (E_{#nu}) #times 10^{-38} (cm^{2}/target)");
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;
// Get the targets specified by the user, separated by commas
// This has structure target1[fraction1], target2[fraction2]
std::vector<std::string> targprs = GeneralUtils::ParseToStr(target, ",");
std::vector<std::string> targ_list;
std::vector<std::string> frac_list;
// Chop up the target string which has format
// TARGET1[fraction1],TARGET2[fraction2]
// std::cout << "Targets: " << std::endl;
// Loop over the vector of strings "TARGET1[fraction1]" "TARGET2[fraction2]"
for (std::vector<std::string>::iterator it = targprs.begin();
it != targprs.end(); ++it) {
// Cut into "TARGET1" and "fraction1]"
std::vector<std::string> targind = GeneralUtils::ParseToStr(*it, "[");
// std::cout << " " << *it << std::endl;
// Cut into "TARGET1" and "fraction1"
for (std::vector<std::string>::iterator jt = targind.begin();
jt != targind.end(); ++jt) {
if ((*jt).find("]") != std::string::npos) {
(*jt) = (*jt).substr(0, (*jt).find("]"));
//*jt = "hello";
frac_list.push_back(*jt);
// Won't find bracket for target
} else {
targ_list.push_back(*jt);
}
}
}
targprs = targ_list;
std::vector<double> targ_fractions;
double minimum = 1.0;
for (std::vector<std::string>::iterator it = frac_list.begin();
it != frac_list.end(); it++) {
// std::cout << " " << *it << std::endl;
double frac = std::atof((*it).c_str());
targ_fractions.push_back(frac);
if (frac < minimum)
minimum = frac;
}
std::vector<double>::iterator it = targ_fractions.begin();
std::vector<std::string>::iterator jt = targ_list.begin();
double scaling = 0;
for (; it != targ_fractions.end(); it++, jt++) {
// First get the mass number from the targ_list
int nucl = atoi((*jt).c_str());
nucl = (nucl % 10000) / 10;
// Gets the relative portions right
*it = (*it) / minimum;
// Scale relative the atomic mass
//(*it) *= (double(nucl)/(*it));
double tempscaling = double(nucl) / (*it);
if (tempscaling > scaling)
scaling = tempscaling;
}
it = targ_fractions.begin();
for (; it != targ_fractions.end(); it++) {
// Round the scaling to nearest integer and multiply
*it *= int(scaling + 0.5);
// Round to nearest integer
*it = int(*it + 0.5);
totalnucl += *it;
}
if (totalnucl == 0) {
THROW("Didn't find any nucleons in input file. Did you really specify the "
"target ratios?\ne.g. TARGET1[fraction1],TARGET2[fraction2]"
<< std::endl);
}
TH1D *totalxsec = (TH1D *)xsechist->Clone();
for (uint i = 0; i < targprs.size(); i++) {
std::string targpdg = targprs[i];
// Check that we found the user requested target in GENIE
bool FoundTarget = false;
for (std::map<std::string, TH1D *>::iterator iter = targetsplines.begin();
iter != targetsplines.end(); iter++) {
std::string targstr = iter->first;
TH1D *xsec = iter->second;
// Match the user targets to the targets found in GENIE
if (targstr.find(targpdg) != std::string::npos) {
FoundTarget = true;
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);
}
}
// Check that targets were all found
if (!FoundTarget) {
ERR(WRN) << "Didn't find target " << targpdg
<< " in the list of targets recorded by GENIE" << std::endl;
ERR(WRN) << " The list of targets you requested is: " << std::endl;
for (uint i = 0; i < targprs.size(); ++i)
ERR(WRN) << " " << targprs[i] << std::endl;
ERR(WRN) << " The list of targets found in GENIE is: " << std::endl;
for (std::map<std::string, TH1D *>::iterator iter = targetsplines.begin();
iter != targetsplines.end(); iter++)
ERR(WRN) << " " << iter->first << std::endl;
}
}
outputfile->cd();
totalxsec->GetYaxis()->SetTitle(
"#sigma (E_{#nu}) #times 10^{-38} (cm^{2}/nucleon)");
totalxsec->Write("nuisance_xsec", TObject::kOverwrite);
eventhist = (TH1D *)fluxhist->Clone();
eventhist->Multiply(totalxsec);
eventhist->GetYaxis()->SetTitle(
(std::string("Event rate (N = #sigma #times #Phi) #times 10^{-38} "
"(cm^{2}/nucleon) #times ") +
eventhist->GetYaxis()->GetTitle())
.c_str());
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;
LOG(FIT) << "XSec Hist Integral = " << totalxsec->Integral("width")
<< std::endl;
outputfile->Close();
return;
}
void RunGENIEPrepare(std::string input, std::string flux, std::string target,
std::string output) {
LOG(FIT) << "Running GENIE Prepare with flux..." << 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");
std::string first_file = "";
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]);
if (!first_file.length()) {
first_file = inputvect[iv_it];
}
}
} else { // The Add form can accept wildcards.
tn->Add(input.c_str());
first_file = input;
}
if (tn->GetFile() == NULL) {
tn->Print();
ERROR(FTL, "gtree not located in GENIE file: " << input);
THROW("Check your inputs, they may need to be completely regenerated!");
throw;
}
int nevt = tn->GetEntries();
if (gNEvents != -999) {
LOG(FIT) << "Overriding number of events by user from " << nevt << " to "
<< gNEvents << std::endl;
nevt = gNEvents;
}
if (!nevt) {
THROW("Couldn't load any events from input specification: \""
<< input.c_str() << "\"");
} else {
QLOG(FIT, "Found " << nevt << " input entries in " << input);
}
NtpMCEventRecord *genientpl = NULL;
tn->SetBranchAddress("gmcrec", &genientpl);
// Make Event and xsec Hist
TH1D *eventhist = (TH1D *)fluxhist->Clone();
eventhist->SetDirectory(NULL);
eventhist->Reset();
TH1D *xsechist = (TH1D *)eventhist->Clone();
xsechist->SetDirectory(NULL);
// 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);
// Hussssch GENIE
StopTalking();
// Get the event
EventRecord &event = *(genientpl->event);
// Get the neutrino
GHepParticle *neu = event.Probe();
StartTalking();
// Get XSec From Spline
// Get the GHepRecord
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;
// Get target nucleus
// Alternative ways of getting the summaries
// GHepParticle *target = genie_record.TargetNucleus();
// int pdg = target->Pdg();
// Fill lists of Unique IDS (neutrino and target)
if (std::find(targetids.begin(), targetids.end(), targ) ==
targetids.end()) {
targetids.push_back(targ);
}
// The full interaction list
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();
modexsec[mode]->SetDirectory(NULL);
modecount[mode]->SetDirectory(NULL);
modexsec[mode]->GetYaxis()->SetTitle(
"d#sigma/dE_{#nu} #times 10^{-38} (events weighted by #sigma)");
modecount[mode]->GetYaxis()->SetTitle("Number of events in file");
}
// 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 (E: " << neu->E() << " GeV, xsec: " << xsec
<< " E-38 cm^2/nucleon)" << 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()) {
// Shut the chain;
delete tn;
outputfile = new TFile(first_file.c_str(), "UPDATE");
} else {
outputfile = new TFile(gOutputFile.c_str(), "RECREATE");
outputfile->cd();
QLOG(FIT, "Cloning input vector to output file: " << gOutputFile);
- TTree *cloneTree = tn->CloneTree();
+ TTree *cloneTree = tn->CloneTree(-1, "fast");
cloneTree->SetDirectory(outputfile);
cloneTree->Write();
// ********************************
// CLUDGE KLUDGE KLUDGE FOR NOVA
QLOG(FIT, "Cloning input nova_wgts to output file: " << gOutputFile);
// Also check for the nova_wgts tree from Jeremy
TChain *nova_chain = new TChain("nova_wgts");
nova_chain->AddFile(input.c_str());
- TTree *nova_tree = nova_chain->CloneTree();
+ TTree *nova_tree = nova_chain->CloneTree(-1, "fast");
if (!nova_tree) {
QLOG(FIT, "Could not find nova_wgts tree in " << input);
} else {
QLOG(FIT, "Found nova_wgts tree in " << input);
nova_tree->SetDirectory(outputfile);
nova_tree->Write();
}
QLOG(FIT, "Done cloning tree.");
}
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]->GetYaxis()->SetTitle(
"#sigma (E_{#nu}) #times 10^{-38} (cm^{2}/target)");
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++) {
std::string targ = targetids[i];
LOG(FIT) << "Getting target " << i << ": " << targ << std::endl;
targetsplines[targ] = (TH1D *)xsechist->Clone();
targetsplines[targ]->GetYaxis()->SetTitle(
"#sigma (E_{#nu}) #times 10^{-38} (cm^{2}/target)");
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;
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;
// This has structure target1[fraction1], target2[fraction2]
std::vector<std::string> targprs = GeneralUtils::ParseToStr(target, ",");
std::vector<std::string> targ_list;
std::vector<std::string> frac_list;
// Chop up the target string which has format
// TARGET1[fraction1],TARGET2[fraction2]
// std::cout << "Targets: " << std::endl;
// Loop over the vector of strings "TARGET1[fraction1]" "TARGET2[fraction2]"
for (std::vector<std::string>::iterator it = targprs.begin();
it != targprs.end(); ++it) {
// Cut into "TARGET1" and "fraction1]"
std::vector<std::string> targind = GeneralUtils::ParseToStr(*it, "[");
// std::cout << " " << *it << std::endl;
// Cut into "TARGET1" and "fraction1"
for (std::vector<std::string>::iterator jt = targind.begin();
jt != targind.end(); ++jt) {
if ((*jt).find("]") != std::string::npos) {
(*jt) = (*jt).substr(0, (*jt).find("]"));
//*jt = "hello";
frac_list.push_back(*jt);
// Won't find bracket for target
} else {
targ_list.push_back(*jt);
}
}
}
targprs = targ_list;
std::vector<double> targ_fractions;
double minimum = 1.0;
for (std::vector<std::string>::iterator it = frac_list.begin();
it != frac_list.end(); it++) {
// std::cout << " " << *it << std::endl;
double frac = std::atof((*it).c_str());
targ_fractions.push_back(frac);
if (frac < minimum)
minimum = frac;
}
std::vector<double>::iterator it = targ_fractions.begin();
std::vector<std::string>::iterator jt = targ_list.begin();
double scaling = 0;
for (; it != targ_fractions.end(); it++, jt++) {
// First get the mass number from the targ_list
int nucl = atoi((*jt).c_str());
nucl = (nucl % 10000) / 10;
// Gets the relative portions right
*it = (*it) / minimum;
// Scale relative the atomic mass
//(*it) *= (double(nucl)/(*it));
double tempscaling = double(nucl) / (*it);
if (tempscaling > scaling)
scaling = tempscaling;
}
it = targ_fractions.begin();
for (; it != targ_fractions.end(); it++) {
// Round the scaling to nearest integer and multiply
*it *= int(scaling + 0.5);
// Round to nearest integer
*it = int(*it + 0.5);
totalnucl += *it;
}
if (totalnucl == 0) {
THROW("Didn't find any nucleons in input file. Did you really specify the "
"target ratios?\ne.g. TARGET1[fraction1],TARGET2[fraction2]"
<< std::endl);
}
TH1D *totalxsec = (TH1D *)xsechist->Clone();
// Loop over the specified targets by the user
for (uint i = 0; i < targprs.size(); i++) {
std::string targpdg = targprs[i];
// Check that we found the user requested target in GENIE
bool FoundTarget = false;
for (std::map<std::string, TH1D *>::iterator iter = targetsplines.begin();
iter != targetsplines.end(); iter++) {
std::string targstr = iter->first;
TH1D *xsec = iter->second;
// Match the user targets to the targets found in GENIE
if (targstr.find(targpdg) != std::string::npos) {
FoundTarget = true;
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);
}
} // Looped over target splines
// Check that targets were all found
if (!FoundTarget) {
ERR(WRN) << "Didn't find target " << targpdg
<< " in the list of targets recorded by GENIE" << std::endl;
ERR(WRN) << " The list of targets you requested is: " << std::endl;
for (uint i = 0; i < targprs.size(); ++i)
ERR(WRN) << " " << targprs[i] << std::endl;
ERR(WRN) << " The list of targets found in GENIE is: " << std::endl;
for (std::map<std::string, TH1D *>::iterator iter = targetsplines.begin();
iter != targetsplines.end(); iter++)
ERR(WRN) << " " << iter->first << std::endl;
}
}
outputfile->cd();
totalxsec->GetYaxis()->SetTitle(
"#sigma (E_{#nu}) #times 10^{-38} (cm^{2}/nucleon)");
totalxsec->Write("nuisance_xsec", TObject::kOverwrite);
eventhist = (TH1D *)fluxhist->Clone();
eventhist->Multiply(totalxsec);
eventhist->GetYaxis()->SetTitle(
(std::string("Event rate (N = #sigma #times #Phi) #times 10^{-38} "
"(cm^{2}/nucleon) #times ") +
eventhist->GetYaxis()->GetTitle())
.c_str());
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;
LOG(FIT) << "XSec Hist Integral = " << totalxsec->Integral() << std::endl;
outputfile->Close();
return;
};
void PrintOptions() {
std::cout << "PrepareGENIEEvents NUISANCE app. " << std::endl
<< "Takes GHep Outputs and prepares events for NUISANCE."
<< std::endl
<< std::endl
<< "PrepareGENIE [-h,-help,--h,--help] [-i "
"inputfile1.root,inputfile2.root,inputfile3.root,...] "
<< "[-f flux_root_file.root,flux_hist_name] [-t "
"target1[frac1],target2[frac2],...]"
<< "[-n number_of_events (experimental)]" << 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 with m GeV "
"neutrino energy."
<< std::endl;
std::cout << " [ -n number_of_evt ] : Run with a reduced number of events "
"for debugging purposes"
<< 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], "-n")) {
gNEvents = GeneralUtils::StrToInt(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 (gTarget.find("[") == std::string::npos ||
gTarget.find("]") == std::string::npos) {
ERR(FTL) << "Didn't specify target ratios in Prepare Mode" << std::endl;
ERR(FTL) << "Are you sure you gave it as -t "
"\"TARGET1[fraction1],TARGET2[fraction]\"?"
<< std::endl;
flagopt = true;
}
if (argc < 1 || flagopt) {
PrintOptions();
exit(-1);
}
return;
}
diff --git a/parameters/config.xml b/parameters/config.xml
index fb75bda..48492b7 100644
--- a/parameters/config.xml
+++ b/parameters/config.xml
@@ -1,229 +1,232 @@
<nuisance>
<!-- # ###################################################### -->
<!-- # NUISANCE CONFIGURATION OPTIONS -->
<!-- # This file is read in by default at runtime -->
<!-- # If you want to override on a case by case bases use -q at runtime -->
<!-- # ###################################################### -->
<!-- # MAIN Configs -->
<!-- # ###################################################### -->
<!-- # Logger goes from -->
<!-- # 1 Quiet -->
<!-- # 2 Fitter -->
<!-- # 3 Samples -->
<!-- # 4 Reconfigure Loops -->
<!-- # 5 Every Event print out (SHOUT) -->
<!-- # -1 DEBUGGING -->
<config VERBOSITY='4'/>
<!-- # ERROR goes from -->
<!-- # 0 NONE -->
<!-- # 1 FATAL -->
<!-- # 2 WARN -->
<config ERROR='2'/>
<config TRACE='0'/>
<config cores='1' />
<config spline_test_throws='50' />
<config spline_cores='1' />
<config spline_chunks='20' />
<config spline_procchunk='-1' />
<config Electron_NThetaBins='4' />
<config Electron_NEnergyBins='4' />
<config Electron_ThetaWidth='1.0' />
<config Electron_EnergyWidth='0.10' />
+<!-- Do we want to remove FSI, undefined and nuclear particles from the GENIE particle stack? -->
<config RemoveFSIParticles='0' />
<config RemoveUndefParticles='0' />
<config RemoveNuclearParticles='0'/>
<config MINERvASaveExtraCCQE='0' />
<!-- # Input Configs -->
<!-- # ###################################################### -->
<!-- # Default Requirements file for the externalDataFitter Package -->
<!-- # MAX Events : -1 is no cut. Events will be scaled automatically to give good xsec predictions. -->
<config MAXEVENTS='-1'/>
<!-- Include empty stacks in the THStack -->
<config includeemptystackhists='0'/>
<!-- # Turn on/off event manager -->
<!-- # EventManager enables us to only loop number of events once for multiple projections of the same measurements -->
<!-- # e.g. MiniBooNE CC1pi+ Q2 and MiniBooNE CC1pi+ Tmu would ordinarily require 2 reconfigures, but with this enabled it requires only one -->
<config EventManager='1'/>
<!-- # Event Directories -->
<!-- # Can setup default directories and use @EVENT_DIR/path to link to it -->
<config EVENT_DIR='/data2/stowell/NIWG/'/>
<config NEUT_EVENT_DIR='/data2/stowell/NIWG/neut/fit_samples_neut5.3.3/'/>
<config GENIE_EVENT_DIR='/data2/stowell/NIWG/genie/fit_samples_R.2.10.0/'/>
<config NUWRO_EVENT_DIR='/data2/stowell/NIWG/nuwro/fit_samples/'/>
<config GIBUU_EVENT_DIR='/data/GIBUU/DIR/'/>
<config SaveNuWroExtra='0' />
<!-- # In PrepareGENIE the reconstructed splines can be saved into the file -->
<config save_genie_splines='1'/>
<!-- # In InputHandler the option to regenerate NuWro flux/xsec plots is available -->
<!-- # Going to move this to its own app soon -->
<!-- # DEVEL CONFIG OPTION, don't touch! -->
<config CacheSize='0'/>
<!-- # ReWeighting Configuration Options -->
<!-- # ###################################################### -->
<!-- # Convert Dials in output statements using dial conversion card -->
<config convert_dials='0'/>
<!-- # Vetos can be used to specify RW dials NOT to be loaded in -->
<!-- # Useful if one specific one has an issue -->
<config FitWeight_fNIWGRW_veto=''/>
<config FitWeight_fNuwroRW_veto=''/>
<config FitWeight_fNeutRW_veto=''/>
<config FitWeight_fGenieRW_veto=''/>
<!-- # Output Options -->
<!-- # ###################################################### -->
<!-- # Save Nominal prediction with all rw engines at default -->
<config savenominal='0'/>
<!-- # Save prefit with values at starting values -->
<config saveprefit='0'/>
<!-- # Here's the full list of drawing options -->
<!-- DATA/MC/EVT/FINE/RATIO/MODES/SHAPE/WGHT/WEIGHTS/FLUX/XSEC/MASK/COV/INCOV/DECOMP/CANVPDG/CANVMC/SETTINGS'/>
<!-- #config drawopts DATA/MC/EVT/FINE/RATIO/MODES/SHAPE/RESIDUAL/MATRIX/FLUX/MASK/MAP -->
<!-- #config drawopts DATA/MC -->
<config drawopts='DATA/MC/EVT/FINE/RATIO/MODES/SHAPE/FLUX/XSEC/MASK/COV/INCOV/DECOMP/CANVPDG/CANVMC/SETTINGS/PROJ/CANVSLICEMC'/>
+<config nuisflat_SavePreFSI='true' />
+
<config InterpolateSigmaQ0Histogram='1' />
<config InterpolateSigmaQ0HistogramRes='100' />
<config InterpolateSigmaQ0HistogramThrow='1' />
<config InterpolateSigmaQ0HistogramNTHROWS='100000' />
<!-- # Save the shape scaling applied with option SHAPE into the main MC hist -->
<config saveshapescaling='0'/>
<config CorrectGENIEMECNorm='1'/>
<!-- # Set style of 1D output histograms -->
<config linecolour='1'/>
<config linestyle='1'/>
<config linewidth='1'/>
<!-- # For GenericFlux -->
<config isLiteMode='0'/>
<!-- # Statistical Options -->
<!-- # ###################################################### -->
<!-- # Add MC Statistical error to likelihoods -->
<config addmcerror='0'/>
<!-- # NUISMIN Configurations -->
<!-- # ###################################################### -->
<config MAXCALLS='1000000'/>
<config MAXITERATIONS='1000000'/>
<config TOLERANCE='0.001'/>
<!-- # Number of events required in low stats routines -->
<config LOWSTATEVENTS='25000'/>
<!-- # Error band generator configs -->
<!-- # ###################################################### -->
<!-- # For -f ErrorBands creates error bands for given measurements -->
<!-- # How many throws do we want (The higher the better precision) -->
<config error_throws='500'/>
<!-- # Are we throwing uniform or according to Gaussian? -->
<!-- # Only use uniform if wanting to study the limits of a dial. -->
<config error_uniform='0'/>
<config WriteSeperateStacks='1'/>
<!-- # Other Individual Case Configs -->
<!-- # ###################################################### -->
<!-- # Covariance throw options for fake data studies with MiniBooNE data. -->
<config thrown_covariance='FULL'/>
<config throw_mc_stat='0.0'/>
<config throw_diag_syst='0'/>
<config throw_corr_syst='0'/>
<config throw_mc_stat='0'/>
<!-- # Apply a shift to the muon momentum before calculation of Q2 -->
<config muon_momentum_shift='0.0'/>
<config muon_momentum_throw='0'/>
<!-- # MINERvA Specific Configs -->
<config MINERvA_CCinc_XSec_2DEavq3_nu_hadron_cut='0'/>
<config MINERvA_CCinc_XSec_2DEavq3_nu_useq3true='0'/>
<config Modes_split_PN_NN='0'/>
<!-- Use only signal events when reconfiguring -->
<config SignalReconfigures='false'/>
<config FullEventOnSignalReconfigure="true"/>
<!-- # SciBooNE specific -->
<config SciBarDensity='1.04'/>
<config SciBarRecoDist='12.0'/>
<config PenetratingMuonEnergy='1.4'/>
<config FlatEfficiency='1.0'/>
<config NumRangeSteps='50'/>
<config UseProton='true'/>
<config UseZackEff='false'/>
<config MINERvADensity='1.04'/>
<config MINERvARecoDist='10.0'/>
<!-- Different way of reweighting in GENIE -->
<config GENIEWeightEngine_CCRESMode="kModeMaMv"/>
<!--
<config GENIEWeightEngine_CCRESMode="kModeMa"/>
+<<<<<<< HEAD
-->
-<config GENIEWeightEngine_CCQEMode="kModeZExp"/> Using z-expansion
+<config GENIEWeightEngine_CCQEMode="kModeMa"/>
<!--
-<config GENIEWeightEngine_CCQEMode="kModeMa"/> Normal MAQE Llewelyn-Smith
<config GENIEWeightEngine_CCQEMode="kModeNormAndMaShape"/> Taking shape into account (don't really use this...)
<config GENIEWeightEngine_CCQEMode="kModeZExp"/> Using z-expansion
-->
<!-- CCQE/2p2h/1pi Gaussian enhancement method -->
<!-- Apply tilt-shift weights or normal Gaussian parameters-->
<config Gaussian_Enhancement="Normal" />
<!--
<config Gaussian_Enhancement="Tilt-Shift" />
-->
<config NToyThrows='100000' />
<!-- Use NOvA Weights or not -->
<config NOvA_Weights="false" />
<!-- Tune name, for GENIE v3+ -->
<config GENIETune="G18_02a_00_000" />
<config GENIEEventGeneratorList="Default" />
<!--
<config GENIEXSecModelCCRES="genie::ReinSehgalRESPXSec" />
<config GENIEXSecModelCOH="genie::ReinSehgalCOHPiPXSec" />
<config GENIEXSecModelCCQE="genie::LwlynSmithQELCCPXSec" />
-->
<!--
<config GENIEXSecModelCCQE="genie::NievesQELCCPXSec" />
<config GENIEXSecModelCCRES="genie::BergerSehgalRESPXSec2014" />
<config GENIEXSecModelCOH="genie::BergerSehgalCOHPiPXSec2015" />
-->
<config UseShapeCovar="0" />
<config CC0piNBINS="156" />
</nuisance>
diff --git a/src/InputHandler/FitEvent.cxx b/src/InputHandler/FitEvent.cxx
index 50ef1bc..d075570 100644
--- a/src/InputHandler/FitEvent.cxx
+++ b/src/InputHandler/FitEvent.cxx
@@ -1,429 +1,437 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is pddrt 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 "FitEvent.h"
#include <iostream>
#include "TObjArray.h"
FitEvent::FitEvent() {
fGenInfo = NULL;
kRemoveFSIParticles = true;
kRemoveUndefParticles = true;
AllocateParticleStack(400);
};
void FitEvent::AddGeneratorInfo(GeneratorInfoBase* gen) {
fGenInfo = gen;
gen->AllocateParticleStack(kMaxParticles);
}
void FitEvent::AllocateParticleStack(int stacksize) {
LOG(DEB) << "Allocating particle stack of size: " << stacksize << std::endl;
kMaxParticles = stacksize;
fParticleList = new FitParticle*[kMaxParticles];
fParticleMom = new double*[kMaxParticles];
fParticleState = new UInt_t[kMaxParticles];
fParticlePDG = new int[kMaxParticles];
+ fPrimaryVertex = new bool[kMaxParticles];
fOrigParticleMom = new double*[kMaxParticles];
fOrigParticleState = new UInt_t[kMaxParticles];
fOrigParticlePDG = new int[kMaxParticles];
+ fOrigPrimaryVertex = new bool[kMaxParticles];
for (size_t i = 0; i < kMaxParticles; i++) {
fParticleList[i] = NULL;
fParticleMom[i] = new double[4];
fOrigParticleMom[i] = new double[4];
}
if (fGenInfo) fGenInfo->AllocateParticleStack(kMaxParticles);
}
void FitEvent::ExpandParticleStack(int stacksize) {
DeallocateParticleStack();
AllocateParticleStack(stacksize);
}
void FitEvent::DeallocateParticleStack() {
for (size_t i = 0; i < kMaxParticles; i++) {
if (fParticleList[i]) delete fParticleList[i];
delete fParticleMom[i];
delete fOrigParticleMom[i];
}
delete fParticleMom;
delete fOrigParticleMom;
delete fParticleList;
delete fParticleState;
delete fParticlePDG;
+ delete fPrimaryVertex;
delete fOrigParticleState;
delete fOrigParticlePDG;
+ delete fOrigPrimaryVertex;
if (fGenInfo) fGenInfo->DeallocateParticleStack();
kMaxParticles = 0;
}
void FitEvent::ClearFitParticles() {
for (size_t i = 0; i < kMaxParticles; i++) {
fParticleList[i] = NULL;
}
}
void FitEvent::FreeFitParticles() {
for (size_t i = 0; i < kMaxParticles; i++) {
FitParticle* fp = fParticleList[i];
if (fp) delete fp;
fParticleList[i] = NULL;
}
}
void FitEvent::ResetParticleList() {
for (unsigned int i = 0; i < kMaxParticles; i++) {
FitParticle* fp = fParticleList[i];
if (fp) delete fp;
fParticleList[i] = NULL;
}
}
void FitEvent::HardReset() {
for (unsigned int i = 0; i < kMaxParticles; i++) {
fParticleList[i] = NULL;
}
}
void FitEvent::ResetEvent() {
Mode = 9999;
fEventNo = -1;
fTotCrs = -1.0;
fTargetA = -1;
fTargetZ = -1;
fTargetH = -1;
fBound = false;
fNParticles = 0;
if (fGenInfo) fGenInfo->Reset();
for (unsigned int i = 0; i < kMaxParticles; i++) {
if (fParticleList[i]) delete fParticleList[i];
fParticleList[i] = NULL;
continue;
fParticlePDG[i] = 0;
fParticleState[i] = kUndefinedState;
fParticleMom[i][0] = 0.0;
fParticleMom[i][1] = 0.0;
fParticleMom[i][2] = 0.0;
fParticleMom[i][3] = 0.0;
+ fPrimaryVertex[i] = false;
fOrigParticlePDG[i] = 0;
fOrigParticleState[i] = kUndefinedState;
fOrigParticleMom[i][0] = 0.0;
fOrigParticleMom[i][1] = 0.0;
fOrigParticleMom[i][2] = 0.0;
fOrigParticleMom[i][3] = 0.0;
+ fOrigPrimaryVertex[i] = false;
}
}
void FitEvent::OrderStack() {
// Copy current stack
int npart = fNParticles;
for (int i = 0; i < npart; i++) {
fOrigParticlePDG[i] = fParticlePDG[i];
fOrigParticleState[i] = fParticleState[i];
fOrigParticleMom[i][0] = fParticleMom[i][0];
fOrigParticleMom[i][1] = fParticleMom[i][1];
fOrigParticleMom[i][2] = fParticleMom[i][2];
fOrigParticleMom[i][3] = fParticleMom[i][3];
+ fOrigPrimaryVertex[i] = fPrimaryVertex[i];
}
// Now run loops for each particle
fNParticles = 0;
int stateorder[6] = {kInitialState, kFinalState, kFSIState,
kNuclearInitial, kNuclearRemnant, kUndefinedState};
for (int s = 0; s < 6; s++) {
for (int i = 0; i < npart; i++) {
if ((UInt_t)fOrigParticleState[i] != (UInt_t)stateorder[s]) continue;
fParticlePDG[fNParticles] = fOrigParticlePDG[i];
fParticleState[fNParticles] = fOrigParticleState[i];
fParticleMom[fNParticles][0] = fOrigParticleMom[i][0];
fParticleMom[fNParticles][1] = fOrigParticleMom[i][1];
fParticleMom[fNParticles][2] = fOrigParticleMom[i][2];
fParticleMom[fNParticles][3] = fOrigParticleMom[i][3];
+ fPrimaryVertex[fNParticles] = fOrigPrimaryVertex[i];
fNParticles++;
}
}
if (LOG_LEVEL(DEB)) {
LOG(DEB) << "Ordered stack" << std::endl;
for (int i = 0; i < fNParticles; i++) {
LOG(DEB) << "Particle " << i << ". " << fParticlePDG[i] << " "
<< fParticleMom[i][0] << " " << fParticleMom[i][1] << " "
<< fParticleMom[i][2] << " " << fParticleMom[i][3] << " "
<< fParticleState[i] << std::endl;
}
}
if (fNParticles != npart) {
ERR(FTL) << "Dropped some particles when ordering the stack!" << std::endl;
}
return;
}
void FitEvent::Print() {
if (LOG_LEVEL(FIT)) {
LOG(FIT) << "FITEvent print" << std::endl;
LOG(FIT) << "Mode: " << Mode << ", Weight: " << InputWeight << std::endl;
LOG(FIT) << "Particles: " << fNParticles << std::endl;
LOG(FIT) << " -> Particle Stack " << std::endl;
for (int i = 0; i < fNParticles; i++) {
LOG(FIT) << " -> -> " << i << ". " << fParticlePDG[i] << " "
<< fParticleState[i] << " "
<< " Mom(" << fParticleMom[i][0] << ", " << fParticleMom[i][1]
<< ", " << fParticleMom[i][2] << ", " << fParticleMom[i][3]
<< ")." << std::endl;
}
}
return;
}
/* Read/Write own event class */
void FitEvent::SetBranchAddress(TChain* tn) {
tn->SetBranchAddress("Mode", &Mode);
tn->SetBranchAddress("EventNo", &fEventNo);
tn->SetBranchAddress("TotCrs", &fTotCrs);
tn->SetBranchAddress("TargetA", &fTargetA);
tn->SetBranchAddress("TargetH", &fTargetH);
tn->SetBranchAddress("Bound", &fBound);
tn->SetBranchAddress("RWWeight", &SavedRWWeight);
tn->SetBranchAddress("InputWeight", &InputWeight);
}
void FitEvent::AddBranchesToTree(TTree* tn) {
tn->Branch("Mode", &Mode, "Mode/I");
tn->Branch("EventNo", &fEventNo, "EventNo/i");
tn->Branch("TotCrs", &fTotCrs, "TotCrs/D");
tn->Branch("TargetA", &fTargetA, "TargetA/I");
tn->Branch("TargetH", &fTargetH, "TargetH/I");
tn->Branch("Bound", &fBound, "Bound/O");
tn->Branch("RWWeight", &RWWeight, "RWWeight/D");
tn->Branch("InputWeight", &InputWeight, "InputWeight/D");
tn->Branch("NParticles", &fNParticles, "NParticles/I");
tn->Branch("ParticleState", fOrigParticleState,
"ParticleState[NParticles]/i");
tn->Branch("ParticlePDG", fOrigParticlePDG, "ParticlePDG[NParticles]/I");
tn->Branch("ParticleMom", fOrigParticleMom, "ParticleMom[NParticles][4]/D");
}
// ------- EVENT ACCESS FUNCTION --------- //
TLorentzVector FitEvent::GetParticleP4(int index) const {
if (index == -1 or index >= fNParticles) return TLorentzVector();
return TLorentzVector(fParticleMom[index][0], fParticleMom[index][1],
fParticleMom[index][2], fParticleMom[index][3]);
}
TVector3 FitEvent::GetParticleP3(int index) const {
if (index == -1 or index >= fNParticles) return TVector3();
return TVector3(fParticleMom[index][0], fParticleMom[index][1],
fParticleMom[index][2]);
}
double FitEvent::GetParticleMom(int index) const {
if (index == -1 or index >= fNParticles) return 0.0;
return sqrt(fParticleMom[index][0] * fParticleMom[index][0] +
fParticleMom[index][1] * fParticleMom[index][1] +
fParticleMom[index][2] * fParticleMom[index][2]);
}
double FitEvent::GetParticleMom2(int index) const {
if (index == -1 or index >= fNParticles) return 0.0;
return fabs((fParticleMom[index][0] * fParticleMom[index][0] +
fParticleMom[index][1] * fParticleMom[index][1] +
fParticleMom[index][2] * fParticleMom[index][2]));
}
double FitEvent::GetParticleE(int index) const {
if (index == -1 or index >= fNParticles) return 0.0;
return fParticleMom[index][3];
}
int FitEvent::GetParticleState(int index) const {
if (index == -1 or index >= fNParticles) return kUndefinedState;
return (fParticleState[index]);
}
int FitEvent::GetParticlePDG(int index) const {
if (index == -1 or index >= fNParticles) return 0;
return (fParticlePDG[index]);
}
FitParticle* FitEvent::GetParticle(int const i) {
// Check Valid Index
if (i == -1) {
return NULL;
}
// Check Valid
if (i > fNParticles) {
ERR(FTL) << "Requesting particle beyond stack!" << std::endl
<< "i = " << i << " N = " << fNParticles << std::endl
<< "Mode = " << Mode << std::endl;
throw;
}
if (!fParticleList[i]) {
/*
std::cout << "Creating particle with values i " << i << " ";
std::cout << fParticleMom[i][0] << " " << fParticleMom[i][1] << " " <<
fParticleMom[i][2] << " " << fParticleMom[i][3] << " ";
std::cout << fParticlePDG[i] << " " << fParticleState[i] << std::endl;
*/
fParticleList[i] = new FitParticle(fParticleMom[i][0], fParticleMom[i][1],
fParticleMom[i][2], fParticleMom[i][3],
fParticlePDG[i], fParticleState[i]);
} else {
/*
std::cout << "Filling particle with values i " << i << " ";
std::cout << fParticleMom[i][0] << " " << fParticleMom[i][1] << " " <<
fParticleMom[i][2] << " " << fParticleMom[i][3] << " ";
std::cout << fParticlePDG[i] << " "<< fParticleState[i] <<std::endl;
*/
fParticleList[i]->SetValues(fParticleMom[i][0], fParticleMom[i][1],
fParticleMom[i][2], fParticleMom[i][3],
fParticlePDG[i], fParticleState[i]);
}
return fParticleList[i];
}
bool FitEvent::HasParticle(int const pdg, int const state) const {
bool found = false;
for (int i = 0; i < fNParticles; i++) {
if (state != -1 && fParticleState[i] != (uint)state) continue;
if (fParticlePDG[i] == pdg) found = true;
}
return found;
}
int FitEvent::NumParticle(int const pdg, int const state) const {
int nfound = 0;
for (int i = 0; i < fNParticles; i++) {
if (state != -1 and fParticleState[i] != (uint)state) continue;
if (pdg == 0 or fParticlePDG[i] == pdg) nfound += 1;
}
return nfound;
}
std::vector<int> FitEvent::GetAllParticleIndices(int const pdg,
int const state) const {
std::vector<int> indexlist;
for (int i = 0; i < fNParticles; i++) {
if (state != -1 and fParticleState[i] != (uint)state) continue;
if (pdg == 0 or fParticlePDG[i] == pdg) {
indexlist.push_back(i);
}
}
return indexlist;
}
std::vector<FitParticle*> FitEvent::GetAllParticle(int const pdg,
int const state) {
std::vector<int> indexlist = GetAllParticleIndices(pdg, state);
std::vector<FitParticle*> plist;
for (std::vector<int>::iterator iter = indexlist.begin();
iter != indexlist.end(); iter++) {
plist.push_back(GetParticle((*iter)));
}
return plist;
}
int FitEvent::GetHMParticleIndex(int const pdg, int const state) const {
double maxmom2 = -9999999.9;
int maxind = -1;
for (int i = 0; i < fNParticles; i++) {
if (state != -1 and fParticleState[i] != (uint)state) continue;
if (pdg == 0 or fParticlePDG[i] == pdg) {
double newmom2 = GetParticleMom2(i);
if (newmom2 > maxmom2) {
maxind = i;
maxmom2 = newmom2;
}
}
}
return maxind;
}
int FitEvent::GetBeamNeutrinoIndex(void) const {
for (int i = 0; i < fNParticles; i++) {
if (fParticleState[i] != kInitialState) continue;
int pdg = abs(fParticlePDG[i]);
if (pdg == 12 or pdg == 14 or pdg == 16) {
return i;
}
}
return 0;
}
int FitEvent::GetBeamElectronIndex(void) const {
return GetHMISParticleIndex(11);
}
int FitEvent::GetBeamPionIndex(void) const {
return GetHMISParticleIndex(PhysConst::pdg_pions);
}
int FitEvent::NumFSMesons() {
int nMesons = 0;
for (int i = 0; i < fNParticles; i++) {
if (fParticleState[i] != kFinalState) continue;
if (abs(fParticlePDG[i]) >= 111 && abs(fParticlePDG[i]) <= 557)
nMesons += 1;
}
return nMesons;
}
int FitEvent::NumFSLeptons(void) const {
int nLeptons = 0;
for (int i = 0; i < fNParticles; i++) {
if (fParticleState[i] != kFinalState) continue;
if (abs(fParticlePDG[i]) == 11 || abs(fParticlePDG[i]) == 13 ||
abs(fParticlePDG[i]) == 15)
nLeptons += 1;
}
return nLeptons;
}
diff --git a/src/InputHandler/FitEvent.h b/src/InputHandler/FitEvent.h
index f0b3db6..a1970f3 100644
--- a/src/InputHandler/FitEvent.h
+++ b/src/InputHandler/FitEvent.h
@@ -1,637 +1,639 @@
// 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 FITEVENT2_H_SEEN
#define FITEVENT2_H_SEEN
/*!
* \addtogroup InputHandler
* @{
*/
#include <algorithm>
#include <iterator>
#include <vector>
#include "FitParticle.h"
#include "TLorentzVector.h"
#include "TSpline.h"
#include "BaseFitEvt.h"
#include "FitLogger.h"
#include "TArrayD.h"
#include "TTree.h"
#include "TChain.h"
#include "TargetUtils.h"
#include "PhysConst.h"
/// Common container for event particles
class FitEvent : public BaseFitEvt {
public:
FitEvent();
~FitEvent() {};
void FreeFitParticles();
void ClearFitParticles();
void ResetEvent(void);
void ResetParticleList(void);
void HardReset();
void OrderStack();
void SetBranchAddress(TChain* tn);
void AddBranchesToTree(TTree* tn);
void Print();
void PrintChris();
void DeallocateParticleStack();
void AllocateParticleStack(int stacksize);
void ExpandParticleStack(int stacksize);
void AddGeneratorInfo(GeneratorInfoBase* gen);
// ---- HELPER/ACCESS FUNCTIONS ---- //
/// Return True Interaction ID
inline int GetMode (void) const { return Mode; };
/// Return Target Atomic Number
inline int GetTargetA (void) const { return fTargetA; };
/// Return Target Nuclear Charge
inline int GetTargetZ (void) const { return fTargetZ; };
/// Get Event Total Cross-section
inline int GetTotCrs (void) const { return fTotCrs; };
/// Is Event Charged Current?
inline bool IsCC() const { return (abs(Mode) <= 30); };
/// Is Event Neutral Current?
inline bool IsNC() const { return (abs(Mode) > 30); };
/// Return Particle 4-momentum for given index in particle stack
TLorentzVector GetParticleP4 (int index) const;
/// Return Particle 3-momentum for given index in particle stack
TVector3 GetParticleP3 (int index) const;
/// Return Particle absolute momentum for given index in particle stack
double GetParticleMom (int index) const;
/// Return Particle absolute momentum-squared for given index in particle stack
double GetParticleMom2 (int index) const;
/// Return Particle energy for given index in particle stack
double GetParticleE (int index) const;
/// Return Particle State for given index in particle stack
int GetParticleState (int index) const;
/// Return Particle PDG for given index in particle stack
int GetParticlePDG (int index) const;
/// Allows the removal of KE up to total KE.
inline void RemoveKE(int index, double KE){
FitParticle *fp = GetParticle(index);
double mass = fp->M();
double oKE = fp->KE();
double nE = mass + (oKE - KE);
if(nE < mass){ // Can't take more KE than it has
nE = mass;
}
double n3Mom = sqrt(nE*nE - mass*mass);
TVector3 np3 = fp->P3().Unit()*n3Mom;
fParticleMom[index][0] = np3[0];
fParticleMom[index][1] = np3[1];
fParticleMom[index][2] = np3[2];
fParticleMom[index][3] = nE;
}
/// Allows the removal of KE up to total KE.
inline void GiveKE(int index, double KE){
RemoveKE(index,-KE);
}
/// Return Particle for given index in particle stack
FitParticle* GetParticle(int const index);
/// Get Total Number of Particles in stack
inline uint NParticles (void) const { return fNParticles; };
/// Check if event contains a particle given a pdg and state.
/// If no state is passed all states are considered.
bool HasParticle (int const pdg = 0, int const state = -1) const ;
template <size_t N>
inline bool HasParticle(int const (&pdgs)[N], int const state = -1) const {
for (size_t i = 0; i < N; i++) {
if (HasParticle( pdgs[i], state )) {
return false;
}
}
return false;
}
/// Get total number of particles given a pdg and state.
/// If no state is passed all states are considered.
int NumParticle (int const pdg = 0, int const state = -1) const;
template <size_t N>
inline int NumParticle(int const (&pdgs)[N], int const state = -1) const {
int ncount = 0;
for (size_t i = 0; i < N; i++) {
ncount += NumParticle( pdgs[i], state );
}
return ncount;
}
/// Return a vector of particle indices that can be used with 'GetParticle'
/// functions given a particle pdg and state.
/// If no state is passed all states are considered.
std::vector<int> GetAllParticleIndices (int const pdg = -1, int const state = -1) const;
template <size_t N>
inline std::vector<int> GetAllParticleIndices(int const (&pdgs)[N], const int state = -1) const {
std::vector<int> plist;
for (size_t i = 0; i < N; i++) {
std::vector<int> plisttemp = GetAllParticleIndices(pdgs[i], state);
plist.insert( plist.end(), plisttemp.begin(), plisttemp.end() );
}
return plist;
}
/// Return a vector of FitParticles given a particle pdg and state.
/// This is memory intensive and slow than GetAllParticleIndices,
/// but is slightly easier to use.
std::vector<FitParticle*> GetAllParticle (int const pdg = -1, int const state = -1);
template <size_t N>
inline std::vector<FitParticle*> GetAllParticle(int const (&pdgs)[N], int const state = -1) {
std::vector<FitParticle*> plist;
for (size_t i = 0; i < N; i++) {
std::vector<FitParticle*> plisttemp = GetAllParticle(pdgs[i], state);
plist.insert( plist.end(), plisttemp.begin(), plisttemp.end() );
}
return plist;
}
inline std::vector<int> GetAllNuElectronIndices (void) { return GetAllParticleIndices(12); };
inline std::vector<int> GetAllNuMuonIndices (void) { return GetAllParticleIndices(14); };
inline std::vector<int> GetAllNuTauIndices (void) { return GetAllParticleIndices(16); };
inline std::vector<int> GetAllElectronIndices (void) { return GetAllParticleIndices(11); };
inline std::vector<int> GetAllMuonIndices (void) { return GetAllParticleIndices(13); };
inline std::vector<int> GetAllTauIndices (void) { return GetAllParticleIndices(15); };
inline std::vector<int> GetAllProtonIndices (void) { return GetAllParticleIndices(2212); };
inline std::vector<int> GetAllNeutronIndices (void) { return GetAllParticleIndices(2112); };
inline std::vector<int> GetAllPiZeroIndices (void) { return GetAllParticleIndices(111); };
inline std::vector<int> GetAllPiPlusIndices (void) { return GetAllParticleIndices(211); };
inline std::vector<int> GetAllPiMinusIndices (void) { return GetAllParticleIndices(-211); };
inline std::vector<int> GetAllPhotonIndices (void) { return GetAllParticleIndices(22); };
inline std::vector<FitParticle*> GetAllNuElectron (void) { return GetAllParticle(12); };
inline std::vector<FitParticle*> GetAllNuMuon (void) { return GetAllParticle(14); };
inline std::vector<FitParticle*> GetAllNuTau (void) { return GetAllParticle(16); };
inline std::vector<FitParticle*> GetAllElectron (void) { return GetAllParticle(11); };
inline std::vector<FitParticle*> GetAllMuon (void) { return GetAllParticle(13); };
inline std::vector<FitParticle*> GetAllTau (void) { return GetAllParticle(15); };
inline std::vector<FitParticle*> GetAllProton (void) { return GetAllParticle(2212); };
inline std::vector<FitParticle*> GetAllNeutron (void) { return GetAllParticle(2112); };
inline std::vector<FitParticle*> GetAllPiZero (void) { return GetAllParticle(111); };
inline std::vector<FitParticle*> GetAllPiPlus (void) { return GetAllParticle(211); };
inline std::vector<FitParticle*> GetAllPiMinus (void) { return GetAllParticle(-211); };
inline std::vector<FitParticle*> GetAllPhoton (void) { return GetAllParticle(22); };
// --- Highest Momentum Search Functions --- //
/// Returns the Index of the highest momentum particle given a pdg and state.
/// If no state is given all states are considered, but that will just return the
/// momentum of the beam in most cases so is not advised.
int GetHMParticleIndex (int const pdg = 0, int const state = -1) const;
template <size_t N>
inline int GetHMParticleIndex (int const (&pdgs)[N], int const state = -1) const {
double mom = -999.9;
int rtnindex = -1;
for (size_t i = 0; i < N; ++i) {
// Use ParticleMom as doesn't require Particle Mem alloc
int pindex = GetHMParticleIndex(pdgs[i], state);
if (pindex != -1){
double momnew = GetParticleMom2(pindex);
if (momnew > mom) {
rtnindex = pindex;
mom = momnew;
}
}
}
return rtnindex;
};
/// Returns the highest momentum particle given a pdg and state.
/// If no state is given all states are considered, but that will just return the
/// momentum of the beam in most cases so is not advised.
inline FitParticle* GetHMParticle(int const pdg = 0, int const state = -1) {
return GetParticle( GetHMParticleIndex(pdg, state) );
}
template <size_t N>
inline FitParticle* GetHMParticle(int const (&pdgs)[N], int const state) {
return GetParticle(GetHMParticleIndex(pdgs, state));
};
// ---- Initial State Helpers --- //
/// Checks the event has a particle of a given pdg in the initial state.
inline bool HasISParticle(int const pdg) const {
return HasParticle(pdg, kInitialState);
};
template <size_t N>
inline bool HasISParticle(int const (&pdgs)[N]) const {
return HasParticle(pdgs, kInitialState);
};
/// Returns the number of particles with a given pdg in the initial state.
inline int NumISParticle(int const pdg = 0) const {
return NumParticle(pdg, kInitialState);
};
template <size_t N>
inline int NumISParticle(int const (&pdgs)[N]) const {
return NumParticle(pdgs, kInitialState);
};
/// Returns a list of indices for all particles with a given pdg
/// in the initial state. These can be used with the 'GetParticle' functions.
inline std::vector<int> GetAllISParticleIndices(int const pdg = -1) const {
return GetAllParticleIndices(pdg, kInitialState);
};
template <size_t N>
inline std::vector<int> GetAllISParticleIndices(int const (&pdgs)[N]) const {
return GetAllParticleIndices(pdgs, kInitialState);
};
/// Returns a list of particles with a given pdg in the initial state.
/// This function is more memory intensive and slower than the Indices function.
inline std::vector<FitParticle*> GetAllISParticle(int const pdg = -1) {
return GetAllParticle(pdg, kInitialState);
};
template <size_t N>
inline std::vector<FitParticle*> GetAllISParticle(int const (&pdgs)[N]) {
return GetAllParticle(pdgs, kInitialState);
};
/// Returns the highest momentum particle with a given pdg in the initial state.
inline FitParticle* GetHMISParticle(int const pdg) {
return GetHMParticle(pdg, kInitialState);
};
template <size_t N>
inline FitParticle* GetHMISParticle(int const (&pdgs)[N]) {
return GetHMParticle(pdgs, kInitialState);
};
/// Returns the highest momentum particle index with a given pdg in the initial state.
inline int GetHMISParticleIndex(int const pdg) const {
return GetHMParticleIndex(pdg, kInitialState);
};
template <size_t N>
inline int GetHMISParticleIndex(int const (&pdgs)[N]) const {
return GetHMParticleIndex(pdgs, kInitialState);
};
inline bool HasISNuElectron (void) const { return HasISParticle(12); };
inline bool HasISNuMuon (void) const { return HasISParticle(14); };
inline bool HasISNuTau (void) const { return HasISParticle(16); };
inline bool HasISElectron (void) const { return HasISParticle(11); };
inline bool HasISMuon (void) const { return HasISParticle(13); };
inline bool HasISTau (void) const { return HasISParticle(15); };
inline bool HasISProton (void) const { return HasISParticle(2212); };
inline bool HasISNeutron (void) const { return HasISParticle(2112); };
inline bool HasISPiZero (void) const { return HasISParticle(111); };
inline bool HasISPiPlus (void) const { return HasISParticle(211); };
inline bool HasISPiMinus (void) const { return HasISParticle(-211); };
inline bool HasISPhoton (void) const { return HasISParticle(22); };
inline bool HasISLeptons (void) const { return HasISParticle(PhysConst::pdg_leptons); };
inline bool HasISPions (void) const { return HasISParticle(PhysConst::pdg_pions); };
inline bool HasISChargePions (void) const { return HasISParticle(PhysConst::pdg_charged_pions); };
inline int NumISNuElectron (void) const { return NumISParticle(12); };
inline int NumISNuMuon (void) const { return NumISParticle(14); };
inline int NumISNuTau (void) const { return NumISParticle(16); };
inline int NumISElectron (void) const { return NumISParticle(11); };
inline int NumISMuon (void) const { return NumISParticle(13); };
inline int NumISTau (void) const { return NumISParticle(15); };
inline int NumISProton (void) const { return NumISParticle(2212); };
inline int NumISNeutron (void) const { return NumISParticle(2112); };
inline int NumISPiZero (void) const { return NumISParticle(111); };
inline int NumISPiPlus (void) const { return NumISParticle(211); };
inline int NumISPiMinus (void) const { return NumISParticle(-211); };
inline int NumISPhoton (void) const { return NumISParticle(22); };
inline int NumISLeptons (void) const { return NumISParticle(PhysConst::pdg_leptons); };
inline int NumISPions (void) const { return NumISParticle(PhysConst::pdg_pions); };
inline int NumISChargePions (void) const { return NumISParticle(PhysConst::pdg_charged_pions); };
inline std::vector<int> GetAllISNuElectronIndices (void) const { return GetAllISParticleIndices(12); };
inline std::vector<int> GetAllISNuMuonIndices (void) const { return GetAllISParticleIndices(14); };
inline std::vector<int> GetAllISNuTauIndices (void) const { return GetAllISParticleIndices(16); };
inline std::vector<int> GetAllISElectronIndices (void) const { return GetAllISParticleIndices(11); };
inline std::vector<int> GetAllISMuonIndices (void) const { return GetAllISParticleIndices(13); };
inline std::vector<int> GetAllISTauIndices (void) const { return GetAllISParticleIndices(15); };
inline std::vector<int> GetAllISProtonIndices (void) const { return GetAllISParticleIndices(2212); };
inline std::vector<int> GetAllISNeutronIndices (void) const { return GetAllISParticleIndices(2112); };
inline std::vector<int> GetAllISPiZeroIndices (void) const { return GetAllISParticleIndices(111); };
inline std::vector<int> GetAllISPiPlusIndices (void) const { return GetAllISParticleIndices(211); };
inline std::vector<int> GetAllISPiMinusIndices (void) const { return GetAllISParticleIndices(-211); };
inline std::vector<int> GetAllISPhotonIndices (void) const { return GetAllISParticleIndices(22); };
inline std::vector<int> GetAllISLeptonsIndices (void) const { return GetAllISParticleIndices(PhysConst::pdg_leptons); };
inline std::vector<int> GetAllISPionsIndices (void) const { return GetAllISParticleIndices(PhysConst::pdg_pions); };
inline std::vector<int> GetAllISChargePionsIndices(void) const { return GetAllISParticleIndices(PhysConst::pdg_charged_pions); };
inline std::vector<FitParticle*> GetAllISNuElectron (void) { return GetAllISParticle(12); };
inline std::vector<FitParticle*> GetAllISNuMuon (void) { return GetAllISParticle(14); };
inline std::vector<FitParticle*> GetAllISNuTau (void) { return GetAllISParticle(16); };
inline std::vector<FitParticle*> GetAllISElectron (void) { return GetAllISParticle(11); };
inline std::vector<FitParticle*> GetAllISMuon (void) { return GetAllISParticle(13); };
inline std::vector<FitParticle*> GetAllISTau (void) { return GetAllISParticle(15); };
inline std::vector<FitParticle*> GetAllISProton (void) { return GetAllISParticle(2212); };
inline std::vector<FitParticle*> GetAllISNeutron (void) { return GetAllISParticle(2112); };
inline std::vector<FitParticle*> GetAllISPiZero (void) { return GetAllISParticle(111); };
inline std::vector<FitParticle*> GetAllISPiPlus (void) { return GetAllISParticle(211); };
inline std::vector<FitParticle*> GetAllISPiMinus (void) { return GetAllISParticle(-211); };
inline std::vector<FitParticle*> GetAllISPhoton (void) { return GetAllISParticle(22); };
inline std::vector<FitParticle*> GetAllISLeptons (void) { return GetAllISParticle(PhysConst::pdg_leptons); };
inline std::vector<FitParticle*> GetAllISPions (void) { return GetAllISParticle(PhysConst::pdg_pions); };
inline std::vector<FitParticle*> GetAllISChargePions(void) { return GetAllISParticle(PhysConst::pdg_charged_pions); };
inline FitParticle* GetHMISNuElectron (void) { return GetHMISParticle(12); };
inline FitParticle* GetHMISNuMuon (void) { return GetHMISParticle(14); };
inline FitParticle* GetHMISNuTau (void) { return GetHMISParticle(16); };
inline FitParticle* GetHMISElectron (void) { return GetHMISParticle(11); };
inline FitParticle* GetHMISMuon (void) { return GetHMISParticle(13); };
inline FitParticle* GetHMISTau (void) { return GetHMISParticle(15); };
inline FitParticle* GetHMISAnyLeptons (void) { return GetHMISParticle(PhysConst::pdg_all_leptons); };
inline FitParticle* GetHMISProton (void) { return GetHMISParticle(2212); };
inline FitParticle* GetHMISNeutron (void) { return GetHMISParticle(2112); };
inline FitParticle* GetHMISPiZero (void) { return GetHMISParticle(111); };
inline FitParticle* GetHMISPiPlus (void) { return GetHMISParticle(211); };
inline FitParticle* GetHMISPiMinus (void) { return GetHMISParticle(-211); };
inline FitParticle* GetHMISPhoton (void) { return GetHMISParticle(22); };
inline FitParticle* GetHMISLepton (void) { return GetHMISParticle(PhysConst::pdg_leptons); };
inline FitParticle* GetHMISPions (void) { return GetHMISParticle(PhysConst::pdg_pions); };
inline FitParticle* GetHMISChargePions(void) { return GetHMISParticle(PhysConst::pdg_charged_pions); };
inline int GetHMISNuElectronIndex (void) { return GetHMISParticleIndex(12); };
inline int GetHMISNuMuonIndex (void) { return GetHMISParticleIndex(14); };
inline int GetHMISNuTauIndex (void) { return GetHMISParticleIndex(16); };
inline int GetHMISElectronIndex (void) { return GetHMISParticleIndex(11); };
inline int GetHMISMuonIndex (void) { return GetHMISParticleIndex(13); };
inline int GetHMISTauIndex (void) { return GetHMISParticleIndex(15); };
inline int GetHMISProtonIndex (void) { return GetHMISParticleIndex(2212); };
inline int GetHMISNeutronIndex (void) { return GetHMISParticleIndex(2112); };
inline int GetHMISPiZeroIndex (void) { return GetHMISParticleIndex(111); };
inline int GetHMISPiPlusIndex (void) { return GetHMISParticleIndex(211); };
inline int GetHMISPiMinusIndex (void) { return GetHMISParticleIndex(-211); };
inline int GetHMISPhotonIndex (void) { return GetHMISParticleIndex(22); };
inline int GetHMISLeptonsIndex (void) { return GetHMISParticleIndex(PhysConst::pdg_leptons); };
inline int GetHMISPionsIndex (void) { return GetHMISParticleIndex(PhysConst::pdg_pions); };
inline int GetHMISChargePionsIndex(void) { return GetHMISParticleIndex(PhysConst::pdg_charged_pions); };
// ---- Final State Helpers --- //
inline bool HasFSParticle(int const pdg) const {
return HasParticle(pdg, kFinalState);
};
template <size_t N>
inline bool HasFSParticle(int const (&pdgs)[N]) const {
return HasParticle(pdgs, kFinalState);
};
inline int NumFSParticle(int const pdg = 0) const {
return NumParticle(pdg, kFinalState);
};
template <size_t N>
inline int NumFSParticle(int const (&pdgs)[N]) const {
return NumParticle(pdgs, kFinalState);
};
inline std::vector<int> GetAllFSParticleIndices(int const pdg = -1) const {
return GetAllParticleIndices(pdg, kFinalState);
};
template <size_t N>
inline std::vector<int> GetAllFSParticleIndices(int const (&pdgs)[N]) const {
return GetAllParticleIndices(pdgs, kFinalState);
};
inline std::vector<FitParticle*> GetAllFSParticle(int const pdg = -1) {
return GetAllParticle(pdg, kFinalState);
};
template <size_t N>
inline std::vector<FitParticle*> GetAllFSParticle(int const (&pdgs)[N]) {
return GetAllParticle(pdgs, kFinalState);
};
inline FitParticle* GetHMFSParticle(int const pdg) {
return GetHMParticle(pdg, kFinalState);
};
template <size_t N>
inline FitParticle* GetHMFSParticle(int const (&pdgs)[N]) {
return GetHMParticle(pdgs, kFinalState);
};
inline int GetHMFSParticleIndex(int const pdg) const {
return GetHMParticleIndex(pdg, kFinalState);
};
template <size_t N>
inline int GetHMFSParticleIndex(int const (&pdgs)[N]) const {
return GetHMParticleIndex(pdgs, kFinalState);
};
inline bool HasFSNuElectron (void) const { return HasFSParticle(12); };
inline bool HasFSNuMuon (void) const { return HasFSParticle(14); };
inline bool HasFSNuTau (void) const { return HasFSParticle(16); };
inline bool HasFSElectron (void) const { return HasFSParticle(11); };
inline bool HasFSMuon (void) const { return HasFSParticle(13); };
inline bool HasFSTau (void) const { return HasFSParticle(15); };
inline bool HasFSProton (void) const { return HasFSParticle(2212); };
inline bool HasFSNeutron (void) const { return HasFSParticle(2112); };
inline bool HasFSPiZero (void) const { return HasFSParticle(111); };
inline bool HasFSPiPlus (void) const { return HasFSParticle(211); };
inline bool HasFSPiMinus (void) const { return HasFSParticle(-211); };
inline bool HasFSPhoton (void) const { return HasFSParticle(22); };
inline bool HasFSLeptons (void) const { return HasFSParticle(PhysConst::pdg_leptons); };
inline bool HasFSPions (void) const { return HasFSParticle(PhysConst::pdg_pions); };
inline bool HasFSChargePions (void) const { return HasFSParticle(PhysConst::pdg_charged_pions); };
inline int NumFSNuElectron (void) const { return NumFSParticle(12); };
inline int NumFSNuMuon (void) const { return NumFSParticle(14); };
inline int NumFSNuTau (void) const { return NumFSParticle(16); };
inline int NumFSElectron (void) const { return NumFSParticle(11); };
inline int NumFSMuon (void) const { return NumFSParticle(13); };
inline int NumFSTau (void) const { return NumFSParticle(15); };
inline int NumFSProton (void) const { return NumFSParticle(2212); };
inline int NumFSNeutron (void) const { return NumFSParticle(2112); };
inline int NumFSPiZero (void) const { return NumFSParticle(111); };
inline int NumFSPiPlus (void) const { return NumFSParticle(211); };
inline int NumFSPiMinus (void) const { return NumFSParticle(-211); };
inline int NumFSPhoton (void) const { return NumFSParticle(22); };
int NumFSLeptons (void) const; // { return NumFSParticle(PhysConst::pdg_leptons); };
inline int NumFSPions (void) const { return NumFSParticle(PhysConst::pdg_pions); };
inline int NumFSChargePions (void) const { return NumFSParticle(PhysConst::pdg_charged_pions); };
inline std::vector<int> GetAllFSNuElectronIndices (void) const { return GetAllFSParticleIndices(12); };
inline std::vector<int> GetAllFSNuMuonIndices (void) const { return GetAllFSParticleIndices(14); };
inline std::vector<int> GetAllFSNuTauIndices (void) const { return GetAllFSParticleIndices(16); };
inline std::vector<int> GetAllFSElectronIndices (void) const { return GetAllFSParticleIndices(11); };
inline std::vector<int> GetAllFSMuonIndices (void) const { return GetAllFSParticleIndices(13); };
inline std::vector<int> GetAllFSTauIndices (void) const { return GetAllFSParticleIndices(15); };
inline std::vector<int> GetAllFSProtonIndices (void) const { return GetAllFSParticleIndices(2212); };
inline std::vector<int> GetAllFSNeutronIndices (void) const { return GetAllFSParticleIndices(2112); };
inline std::vector<int> GetAllFSPiZeroIndices (void) const { return GetAllFSParticleIndices(111); };
inline std::vector<int> GetAllFSPiPlusIndices (void) const { return GetAllFSParticleIndices(211); };
inline std::vector<int> GetAllFSPiMinusIndices (void) const { return GetAllFSParticleIndices(-211); };
inline std::vector<int> GetAllFSPhotonIndices (void) const { return GetAllFSParticleIndices(22); };
inline std::vector<int> GetAllFSLeptonsIndices (void) const { return GetAllFSParticleIndices(PhysConst::pdg_leptons); };
inline std::vector<int> GetAllFSPionsIndices (void) const { return GetAllFSParticleIndices(PhysConst::pdg_pions); };
inline std::vector<int> GetAllFSChargePionsIndices(void) const { return GetAllFSParticleIndices(PhysConst::pdg_charged_pions); };
inline std::vector<FitParticle*> GetAllFSNuElectron (void) { return GetAllFSParticle(12); };
inline std::vector<FitParticle*> GetAllFSNuMuon (void) { return GetAllFSParticle(14); };
inline std::vector<FitParticle*> GetAllFSNuTau (void) { return GetAllFSParticle(16); };
inline std::vector<FitParticle*> GetAllFSElectron (void) { return GetAllFSParticle(11); };
inline std::vector<FitParticle*> GetAllFSMuon (void) { return GetAllFSParticle(13); };
inline std::vector<FitParticle*> GetAllFSTau (void) { return GetAllFSParticle(15); };
inline std::vector<FitParticle*> GetAllFSProton (void) { return GetAllFSParticle(2212); };
inline std::vector<FitParticle*> GetAllFSNeutron (void) { return GetAllFSParticle(2112); };
inline std::vector<FitParticle*> GetAllFSPiZero (void) { return GetAllFSParticle(111); };
inline std::vector<FitParticle*> GetAllFSPiPlus (void) { return GetAllFSParticle(211); };
inline std::vector<FitParticle*> GetAllFSPiMinus (void) { return GetAllFSParticle(-211); };
inline std::vector<FitParticle*> GetAllFSPhoton (void) { return GetAllFSParticle(22); };
inline std::vector<FitParticle*> GetAllFSLeptons (void) { return GetAllFSParticle(PhysConst::pdg_leptons); };
inline std::vector<FitParticle*> GetAllFSPions (void) { return GetAllFSParticle(PhysConst::pdg_pions); };
inline std::vector<FitParticle*> GetAllFSChargePions (void) { return GetAllFSParticle(PhysConst::pdg_charged_pions); };
inline FitParticle* GetHMFSNuElectron (void) { return GetHMFSParticle(12); };
inline FitParticle* GetHMFSNuMuon (void) { return GetHMFSParticle(14); };
inline FitParticle* GetHMFSNuTau (void) { return GetHMFSParticle(16); };
inline FitParticle* GetHMFSElectron (void) { return GetHMFSParticle(11); };
inline FitParticle* GetHMFSMuon (void) { return GetHMFSParticle(13); };
inline FitParticle* GetHMFSTau (void) { return GetHMFSParticle(15); };
inline FitParticle* GetHMFSAnyLeptons (void) { return GetHMFSParticle(PhysConst::pdg_all_leptons); };
inline FitParticle* GetHMFSProton (void) { return GetHMFSParticle(2212); };
inline FitParticle* GetHMFSNeutron (void) { return GetHMFSParticle(2112); };
inline FitParticle* GetHMFSPiZero (void) { return GetHMFSParticle(111); };
inline FitParticle* GetHMFSPiPlus (void) { return GetHMFSParticle(211); };
inline FitParticle* GetHMFSPiMinus (void) { return GetHMFSParticle(-211); };
inline FitParticle* GetHMFSPhoton (void) { return GetHMFSParticle(22); };
inline FitParticle* GetHMFSLeptons (void) { return GetHMFSParticle(PhysConst::pdg_leptons); };
inline FitParticle* GetHMFSAnyLepton (void) { return GetHMFSParticle(PhysConst::pdg_all_leptons); };
inline FitParticle* GetHMFSPions (void) { return GetHMFSParticle(PhysConst::pdg_pions); };
inline FitParticle* GetHMFSChargePions(void) { return GetHMFSParticle(PhysConst::pdg_charged_pions); };
inline int GetHMFSNuElectronIndex (void) const { return GetHMFSParticleIndex(12); };
inline int GetHMFSNuMuonIndex (void) const { return GetHMFSParticleIndex(14); };
inline int GetHMFSNuTauIndex (void) const { return GetHMFSParticleIndex(16); };
inline int GetHMFSElectronIndex (void) const { return GetHMFSParticleIndex(11); };
inline int GetHMFSMuonIndex (void) const { return GetHMFSParticleIndex(13); };
inline int GetHMFSTauIndex (void) const { return GetHMFSParticleIndex(15); };
inline int GetHMFSProtonIndex (void) const { return GetHMFSParticleIndex(2212); };
inline int GetHMFSNeutronIndex (void) const { return GetHMFSParticleIndex(2112); };
inline int GetHMFSPiZeroIndex (void) const { return GetHMFSParticleIndex(111); };
inline int GetHMFSPiPlusIndex (void) const { return GetHMFSParticleIndex(211); };
inline int GetHMFSPiMinusIndex (void) const { return GetHMFSParticleIndex(-211); };
inline int GetHMFSPhotonIndex (void) const { return GetHMFSParticleIndex(22); };
inline int GetHMFSLeptonsIndex (void) const { return GetHMFSParticleIndex(PhysConst::pdg_leptons); };
inline int GetHMFSAnyLeptonIndex (void) const { return GetHMFSParticleIndex(PhysConst::pdg_all_leptons); };
inline int GetHMFSPionsIndex (void) const { return GetHMFSParticleIndex(PhysConst::pdg_pions); };
inline int GetHMFSChargePionsIndex(void) const { return GetHMFSParticleIndex(PhysConst::pdg_charged_pions); };
// ---- NEUTRINO INCOMING Related Functions
int GetBeamNeutrinoIndex (void) const;
inline TLorentzVector GetBeamNeutrinoP4 (void) const { return GetParticleP4(GetBeamNeutrinoIndex()); };
inline TVector3 GetBeamNeutrinoP3 (void) const { return GetParticleP3(GetBeamNeutrinoIndex()); };
inline double GetBeamNeutrinoMom (void) const { return GetParticleMom(GetBeamNeutrinoIndex()); };
inline double GetBeamNeutrinoMom2 (void) const { return GetParticleMom2(GetBeamNeutrinoIndex()); };
inline double GetBeamNeutrinoE (void) const { return GetParticleE(GetBeamNeutrinoIndex()); };
inline double Enu (void) const { return GetBeamNeutrinoE(); };
inline int GetBeamNeutrinoPDG (void) const { return GetParticlePDG(GetBeamNeutrinoIndex()); };
inline int PDGnu (void) const { return GetBeamNeutrinoPDG(); };
inline int GetNeutrinoInPos (void) const { return GetBeamNeutrinoIndex(); };
inline FitParticle* GetBeamNeutrino (void) { return GetParticle(GetBeamNeutrinoIndex()); };
inline FitParticle* GetNeutrinoIn (void) { return GetParticle(GetBeamNeutrinoIndex()); };
// ---- Electron INCOMING Related Functions
int GetBeamElectronIndex (void) const;
inline TLorentzVector GetBeamElectronP4 (void) const { return GetParticleP4(GetBeamElectronIndex()); };
inline TVector3 GetBeamElectronP3 (void) const { return GetParticleP3(GetBeamElectronIndex()); };
inline double GetBeamElectronMom (void) const { return GetParticleMom(GetBeamElectronIndex()); };
inline double GetBeamElectronMom2 (void) const { return GetParticleMom2(GetBeamElectronIndex()); };
inline double GetBeamElectronE (void) const { return GetParticleE(GetBeamElectronIndex()); };
inline FitParticle* GetBeamElectron (void) { return GetParticle(GetBeamElectronIndex()); };
// ---- Pion INCOMING Functions
int GetBeamPionIndex (void) const;
inline TLorentzVector GetBeamPionP4 (void) const { return GetParticleP4(GetBeamPionIndex()); };
inline TVector3 GetBeamPionP3 (void) const { return GetParticleP3(GetBeamPionIndex()); };
inline double GetBeamPionMom (void) const { return GetParticleMom(GetBeamPionIndex()); };
inline double GetBeamPionMom2 (void) const { return GetParticleMom2(GetBeamPionIndex()); };
inline double GetBeamPionE (void) const { return GetParticleE(GetBeamPionIndex()); };
inline FitParticle* GetBeamPion (void) { return GetParticle(GetBeamPionIndex()); };
/// Legacy Functions
inline FitParticle* PartInfo(uint i) { return GetParticle(i); };
inline UInt_t Npart (void) const { return NPart(); };
inline UInt_t NPart (void) const { return fNParticles; };
// Other Functions
int NumFSMesons();
inline int GetBeamPartPos(void) const { return 0; };
FitParticle* GetBeamPart(void);
int GetLeptonOutPos(void) const;
FitParticle* GetLeptonOut(void);
// Event Information
UInt_t fEventNo;
double fTotCrs;
int fTargetA;
int fTargetZ;
int fTargetH;
bool fBound;
int fDistance;
int fTargetPDG;
// Reduced Particle Stack
UInt_t kMaxParticles;
int fNParticles;
double** fParticleMom;
UInt_t* fParticleState;
int* fParticlePDG;
FitParticle** fParticleList;
+ bool *fPrimaryVertex;
double** fOrigParticleMom;
UInt_t* fOrigParticleState;
int* fOrigParticlePDG;
+ bool* fOrigPrimaryVertex;
double* fNEUT_ParticleStatusCode;
double* fNEUT_ParticleAliveCode;
GeneratorInfoBase* fGenInfo;
// Config Options for this class
bool kRemoveFSIParticles;
bool kRemoveUndefParticles;
};
/*! @} */
#endif
diff --git a/src/InputHandler/FitParticle.h b/src/InputHandler/FitParticle.h
index fa9c263..5141e81 100644
--- a/src/InputHandler/FitParticle.h
+++ b/src/InputHandler/FitParticle.h
@@ -1,109 +1,109 @@
// 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 FITPARTICLE_H_SEEN
#define FITPARTICLE_H_SEEN
/*!
* \addtogroup InputHandler
* @{
*/
#include "TLorentzVector.h"
/// Partle state flags for its position in the event
enum particle_state{
kUndefinedState = 5,
kInitialState = 0,
kFSIState = 1,
kFinalState = 2,
kNuclearInitial = 3,
- kNuclearRemnant = 4
+ kNuclearRemnant = 4,
};
/// Condensed FitParticle class which acts a common format between the generators
class FitParticle {
public:
/// Create particle of given pdg from momentum variables and state
FitParticle(double x, double y, double z, double t, int pdg, Int_t state);
/// Create empty particle (zero momentum)
FitParticle(){};
~FitParticle(){};
/// Used to change values after creation
void SetValues(double x, double y, double z, double t, int pdg, Int_t state);
/// Return Status Code according to particle_state enum
inline int Status (void) const { return fStatus; };
/// Get Particle PDF
inline int PDG (void) const { return fPID; };
/// Check if particle makes it to final state
inline bool IsFinalState (void) const { return (fStatus == kFinalState); };
/// Was particle involved in intermediate state
inline bool IsFSIState (void) const { return (fStatus == kFSIState); };
/// Was particle incoming
inline bool IsInitialState (void) const { return (fStatus == kInitialState); };
/// Get Mass
inline double M (void){ return fP.Mag(); };
/// Get Kinetic Energy
inline double KE (void){ return fP.E() - fP.Mag(); };
/// Get Total Energy
inline double E (void){ return fP.E(); };
/// Get 4 Momentum
inline TLorentzVector P4(void) {return fP;};
/// Get 3 Momentum
inline TVector3 P3(void) {return fP.Vect();};
/// Get 3 momentum magnitude
inline double p(void) { return fP.Vect().Mag(); };
/// Get 3 momentum magnitude squared
inline double p2(void) { return fP.Vect().Mag2(); };
/// Data Members
TLorentzVector fP; ///< Particle 4 Momentum
int fPID; ///< Particle PDG Code
int fIsAlive; ///< Whether the particle is alive at the end of the event (Yes 1, No 0, Other? -1)
int fNEUTStatusCode; ///< Particle Status (Incoming 1, FSI 2, Outgoing 0, Other 3)
double fMass; ///< Particle Mass
int fStatus; ///< State corresponding to particle_state enum
-
+ bool fIsPrimary; ///< Primary target
};
inline std::ostream& operator<<(std::ostream& os, FitParticle const& p){
return os << " Particle[pdgc:" << p.fPID
<< ", stat:"<<p.fStatus
<< ", 4mom:("<< p.fP.X() << "," << p.fP.Y() << "," << p.fP.Z() << "," << p.fP.T() << ")]";
}
/*! @} */
#endif
diff --git a/src/InputHandler/GENIEInputHandler.cxx b/src/InputHandler/GENIEInputHandler.cxx
index 240611c..b3b5420 100644
--- a/src/InputHandler/GENIEInputHandler.cxx
+++ b/src/InputHandler/GENIEInputHandler.cxx
@@ -1,580 +1,583 @@
// 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/>.
*******************************************************************************/
#ifdef __GENIE_ENABLED__
#include "GENIEInputHandler.h"
#pragma push_macro("LOG")
#undef LOG
#pragma push_macro("ERROR")
#undef ERROR
#ifdef GENIE_PRE_R3
#include "Messenger/Messenger.h"
#else
#include "Framework/Messenger/Messenger.h"
#endif
#pragma pop_macro("LOG")
#pragma pop_macro("ERROR")
#include "InputUtils.h"
GENIEGeneratorInfo::~GENIEGeneratorInfo() { DeallocateParticleStack(); }
void GENIEGeneratorInfo::AddBranchesToTree(TTree* tn) {
tn->Branch("GenieParticlePDGs", &fGenieParticlePDGs, "GenieParticlePDGs/I");
}
void GENIEGeneratorInfo::SetBranchesFromTree(TTree* tn) {
tn->SetBranchAddress("GenieParticlePDGs", &fGenieParticlePDGs);
}
void GENIEGeneratorInfo::AllocateParticleStack(int stacksize) {
fGenieParticlePDGs = new int[stacksize];
}
void GENIEGeneratorInfo::DeallocateParticleStack() {
delete fGenieParticlePDGs;
}
void GENIEGeneratorInfo::FillGeneratorInfo(NtpMCEventRecord* ntpl) {
Reset();
// Check for GENIE Event
if (!ntpl) return;
if (!ntpl->event) return;
// Cast Event Record
GHepRecord* ghep = static_cast<GHepRecord*>(ntpl->event);
if (!ghep) return;
// Fill Particle Stack
GHepParticle* p = 0;
TObjArrayIter iter(ghep);
// Loop over all particles
int i = 0;
while ((p = (dynamic_cast<genie::GHepParticle*>((iter).Next())))) {
if (!p) continue;
// Get PDG
fGenieParticlePDGs[i] = p->Pdg();
i++;
}
}
void GENIEGeneratorInfo::Reset() {
for (int i = 0; i < kMaxParticles; i++) {
fGenieParticlePDGs[i] = 0;
}
}
GENIEInputHandler::GENIEInputHandler(std::string const& handle,
std::string const& rawinputs) {
LOG(SAM) << "Creating GENIEInputHandler : " << handle << std::endl;
genie::Messenger::Instance()->SetPriorityLevel("GHepUtils", pFATAL);
// Run a joint input handling
fName = handle;
// Setup the TChain
fGENIETree = new TChain("gtree");
fSaveExtra = FitPar::Config().GetParB("SaveExtraGenie");
fCacheSize = FitPar::Config().GetParI("CacheSize");
fMaxEvents = FitPar::Config().GetParI("MAXEVENTS");
// Are we running with NOvA weights
fNOvAWeights = FitPar::Config().GetParB("NOvA_Weights");
MAQEw = 1.0;
NonResw = 1.0;
RPAQEw = 1.0;
RPARESw = 1.0;
MECw = 1.0;
DISw = 1.0;
NOVAw = 1.0;
// Loop over all inputs and grab flux, eventhist, and nevents
std::vector<std::string> inputs = InputUtils::ParseInputFileList(rawinputs);
for (size_t inp_it = 0; inp_it < inputs.size(); ++inp_it) {
// Open File for histogram access
TFile* inp_file = new TFile(
InputUtils::ExpandInputDirectories(inputs[inp_it]).c_str(), "READ");
if (!inp_file or inp_file->IsZombie()) {
THROW("GENIE File IsZombie() at : '"
<< inputs[inp_it] << "'" << std::endl
<< "Check that your file paths are correct and the file exists!"
<< std::endl
<< "$ ls -lh " << inputs[inp_it]);
}
// Get Flux/Event hist
TH1D* fluxhist = (TH1D*)inp_file->Get("nuisance_flux");
TH1D* eventhist = (TH1D*)inp_file->Get("nuisance_events");
if (!fluxhist or !eventhist) {
ERROR(FTL, "Input File Contents: " << inputs[inp_it]);
inp_file->ls();
THROW("GENIE FILE doesn't contain flux/xsec info."
<< std::endl
<< "Try running the app PrepareGENIE first on :" << inputs[inp_it]
<< std::endl
<< "$ PrepareGENIE -h");
}
// Get N Events
TTree* genietree = (TTree*)inp_file->Get("gtree");
if (!genietree) {
ERROR(FTL, "gtree not located in GENIE file: " << inputs[inp_it]);
THROW("Check your inputs, they may need to be completely regenerated!");
throw;
}
int nevents = genietree->GetEntries();
if (nevents <= 0) {
THROW("Trying to a TTree with "
<< nevents << " to TChain from : " << inputs[inp_it]);
}
// Check for precomputed weights
TTree *weighttree = (TTree*)inp_file->Get("nova_wgts");
if (fNOvAWeights) {
if (!weighttree) {
THROW("Did not find nova_wgts tree in file " << inputs[inp_it] << " but you specified it" << std::endl);
} else {
LOG(FIT) << "Found nova_wgts tree in file " << inputs[inp_it] << std::endl;
}
}
// Register input to form flux/event rate hists
RegisterJointInput(inputs[inp_it], nevents, fluxhist, eventhist);
// Add To TChain
fGENIETree->AddFile(inputs[inp_it].c_str());
if (weighttree != NULL) fGENIETree->AddFriend(weighttree);
}
// Registor all our file inputs
SetupJointInputs();
// Assign to tree
fEventType = kGENIE;
fGenieNtpl = NULL;
fGENIETree->SetBranchAddress("gmcrec", &fGenieNtpl);
// Set up the custom weights
if (fNOvAWeights) {
fGENIETree->SetBranchAddress("MAQEwgt", &MAQEw);
fGENIETree->SetBranchAddress("nonResNormWgt", &NonResw);
fGENIETree->SetBranchAddress("RPAQEWgt", &RPAQEw);
fGENIETree->SetBranchAddress("RPARESWgt", &RPARESw);
fGENIETree->SetBranchAddress("MECWgt", &MECw);
fGENIETree->SetBranchAddress("DISWgt", &DISw);
fGENIETree->SetBranchAddress("nova2018CVWgt", &NOVAw);
}
// Libraries should be seen but not heard...
StopTalking();
fGENIETree->GetEntry(0);
StartTalking();
// Create Fit Event
fNUISANCEEvent = new FitEvent();
fNUISANCEEvent->SetGenieEvent(fGenieNtpl);
if (fSaveExtra) {
fGenieInfo = new GENIEGeneratorInfo();
fNUISANCEEvent->AddGeneratorInfo(fGenieInfo);
}
fNUISANCEEvent->HardReset();
};
GENIEInputHandler::~GENIEInputHandler() {
// if (fGenieGHep) delete fGenieGHep;
// if (fGenieNtpl) delete fGenieNtpl;
// if (fGENIETree) delete fGENIETree;
// if (fGenieInfo) delete fGenieInfo;
}
void GENIEInputHandler::CreateCache() {
if (fCacheSize > 0) {
// fGENIETree->SetCacheEntryRange(0, fNEvents);
fGENIETree->AddBranchToCache("*", 1);
fGENIETree->SetCacheSize(fCacheSize);
}
}
void GENIEInputHandler::RemoveCache() {
// fGENIETree->SetCacheEntryRange(0, fNEvents);
fGENIETree->AddBranchToCache("*", 0);
fGENIETree->SetCacheSize(0);
}
FitEvent* GENIEInputHandler::GetNuisanceEvent(const UInt_t entry, const bool lightweight) {
if (entry >= (UInt_t)fNEvents) return NULL;
// Clear the previous event (See Note 1 in ROOT TClonesArray documentation)
if (fGenieNtpl) {
fGenieNtpl->Clear();
}
// Read Entry from TTree to fill NEUT Vect in BaseFitEvt;
fGENIETree->GetEntry(entry);
// Run NUISANCE Vector Filler
if (!lightweight) {
CalcNUISANCEKinematics();
}
#ifdef __PROB3PP_ENABLED__
else {
// Check for GENIE Event
if (!fGenieNtpl) return NULL;
if (!fGenieNtpl->event) return NULL;
// Cast Event Record
fGenieGHep = static_cast<GHepRecord*>(fGenieNtpl->event);
if (!fGenieGHep) return NULL;
TObjArrayIter iter(fGenieGHep);
genie::GHepParticle* p;
while ((p = (dynamic_cast<genie::GHepParticle*>((iter).Next())))) {
if (!p) {
continue;
}
// Get Status
int state = GetGENIEParticleStatus(p, fNUISANCEEvent->Mode);
if (state != genie::kIStInitialState) {
continue;
}
fNUISANCEEvent->probe_E = p->E() * 1.E3;
fNUISANCEEvent->probe_pdg = p->Pdg();
break;
}
}
#endif
// Setup Input scaling for joint inputs
fNUISANCEEvent->InputWeight = GetInputWeight(entry);
return fNUISANCEEvent;
}
int GENIEInputHandler::GetGENIEParticleStatus(genie::GHepParticle* p, int mode) {
/*
kIStUndefined = -1,
kIStInitialState = 0, / generator-level initial state /
kIStStableFinalState = 1, / generator-level final state:
particles to be tracked by detector-level MC /
kIStIntermediateState = 2,
kIStDecayedState = 3,
kIStCorrelatedNucleon = 10,
kIStNucleonTarget = 11,
kIStDISPreFragmHadronicState = 12,
kIStPreDecayResonantState = 13,
kIStHadronInTheNucleus = 14, / hadrons inside the nucleus: marked
for hadron transport modules to act on /
kIStFinalStateNuclearRemnant = 15, / low energy nuclear fragments
entering the record collectively as a 'hadronic blob' pseudo-particle /
kIStNucleonClusterTarget = 16, // for composite nucleons before
phase space decay
*/
int state = kUndefinedState;
switch (p->Status()) {
case genie::kIStNucleonTarget:
case genie::kIStInitialState:
case genie::kIStCorrelatedNucleon:
case genie::kIStNucleonClusterTarget:
state = kInitialState;
break;
case genie::kIStStableFinalState:
state = kFinalState;
break;
case genie::kIStHadronInTheNucleus:
if (abs(mode) == 2)
state = kInitialState;
else
state = kFSIState;
break;
case genie::kIStPreDecayResonantState:
case genie::kIStDISPreFragmHadronicState:
case genie::kIStIntermediateState:
state = kFSIState;
break;
case genie::kIStFinalStateNuclearRemnant:
case genie::kIStUndefined:
case genie::kIStDecayedState:
default:
break;
}
// Flag to remove nuclear part in genie
if (p->Pdg() > 1000000) {
if (state == kInitialState)
state = kNuclearInitial;
else if (state == kFinalState)
state = kNuclearRemnant;
}
return state;
}
#endif
#ifdef __GENIE_ENABLED__
int GENIEInputHandler::ConvertGENIEReactionCode(GHepRecord* gheprec) {
// Electron Scattering
if (gheprec->Summary()->ProcInfo().IsEM()) {
if (gheprec->Summary()->InitState().ProbePdg() == 11) {
if (gheprec->Summary()->ProcInfo().IsQuasiElastic())
return 1;
else if (gheprec->Summary()->ProcInfo().IsMEC())
return 2;
else if (gheprec->Summary()->ProcInfo().IsResonant())
return 13;
else if (gheprec->Summary()->ProcInfo().IsDeepInelastic())
return 26;
else {
ERROR(WRN,
"Unknown GENIE Electron Scattering Mode!"
<< std::endl
<< "ScatteringTypeId = "
<< gheprec->Summary()->ProcInfo().ScatteringTypeId() << " "
<< "InteractionTypeId = "
<< gheprec->Summary()->ProcInfo().InteractionTypeId()
<< std::endl
<< genie::ScatteringType::AsString(
gheprec->Summary()->ProcInfo().ScatteringTypeId())
<< " "
<< genie::InteractionType::AsString(
gheprec->Summary()->ProcInfo().InteractionTypeId())
<< " " << gheprec->Summary()->ProcInfo().IsMEC());
return 0;
}
}
// Weak CC
} else if (gheprec->Summary()->ProcInfo().IsWeakCC()) {
// CC MEC
if (gheprec->Summary()->ProcInfo().IsMEC()) {
if (pdg::IsNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return 2;
else if (pdg::IsAntiNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return -2;
// CC OTHER
} else {
return utils::ghep::NeutReactionCode(gheprec);
}
// Weak NC
} else if (gheprec->Summary()->ProcInfo().IsWeakNC()) {
// NC MEC
if (gheprec->Summary()->ProcInfo().IsMEC()) {
if (pdg::IsNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return 32;
else if (pdg::IsAntiNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return -32;
// NC OTHER
} else {
return utils::ghep::NeutReactionCode(gheprec);
}
}
return 0;
}
void GENIEInputHandler::CalcNUISANCEKinematics() {
// Reset all variables
fNUISANCEEvent->ResetEvent();
// Check for GENIE Event
if (!fGenieNtpl) return;
if (!fGenieNtpl->event) return;
// Cast Event Record
fGenieGHep = static_cast<GHepRecord*>(fGenieNtpl->event);
if (!fGenieGHep) return;
// Convert GENIE Reaction Code
fNUISANCEEvent->Mode = ConvertGENIEReactionCode(fGenieGHep);
// Set Event Info
fNUISANCEEvent->fEventNo = 0.0;
fNUISANCEEvent->fTotCrs = fGenieGHep->XSec();
// Have a bool storing if interaction happened on free or bound nucleon
bool IsFree = false;
// Set the TargetPDG
if (fGenieGHep->TargetNucleus() != NULL) {
fNUISANCEEvent->fTargetPDG = fGenieGHep->TargetNucleus()->Pdg();
IsFree = false;
// Sometimes GENIE scatters off free nucleons, electrons, photons
// In which TargetNucleus is NULL and we need to find the initial state particle
} else {
// Check the particle is an initial state particle
// Follows GHepRecord::TargetNucleusPosition but doesn't do check on pdg::IsIon
GHepParticle *p = fGenieGHep->Particle(1);
// Check that particle 1 actually exists
if (!p) {
ERR(FTL) << "Can't find particle 1 for GHepRecord" << std::endl;
throw;
}
// If not an ion but is an initial state particle
if (!pdg::IsIon(p->Pdg()) &&
p->Status() == kIStInitialState) {
IsFree = true;
fNUISANCEEvent->fTargetPDG = p->Pdg();
// Catch if something strange happens:
// Here particle 1 is not an initial state particle OR
// particle 1 is an ion OR
// both
} else {
if (pdg::IsIon(p->Pdg())) {
ERR(FTL) << "Particle 1 in GHepRecord stack is an ion but isn't an initial state particle" << std::endl;
throw;
} else {
ERR(FTL) << "Particle 1 in GHepRecord stack is not an ion but is an initial state particle" << std::endl;
throw;
}
}
}
// Set the A and Z and H from the target PDG
// Depends on if we scattered off a free or bound nucleon
if (!IsFree) {
fNUISANCEEvent->fTargetA = TargetUtils::GetTargetAFromPDG(fNUISANCEEvent->fTargetPDG);
fNUISANCEEvent->fTargetZ = TargetUtils::GetTargetZFromPDG(fNUISANCEEvent->fTargetPDG);
fNUISANCEEvent->fTargetH = 0;
} else {
// If free proton scattering
if (fNUISANCEEvent->fTargetPDG == 2212) {
fNUISANCEEvent->fTargetA = 1;
fNUISANCEEvent->fTargetZ = 1;
fNUISANCEEvent->fTargetH = 1;
// If free neutron scattering
} else if (fNUISANCEEvent->fTargetPDG == 2112) {
fNUISANCEEvent->fTargetA = 0;
fNUISANCEEvent->fTargetZ = 1;
fNUISANCEEvent->fTargetH = 0;
// If neither
} else {
fNUISANCEEvent->fTargetA = 0;
fNUISANCEEvent->fTargetZ = 0;
fNUISANCEEvent->fTargetH = 0;
}
}
fNUISANCEEvent->fBound = !IsFree;
fNUISANCEEvent->InputWeight = 1.0; //(1E+38 / genie::units::cm2) * fGenieGHep->XSec();
// And the custom weights
if (fNOvAWeights) {
fNUISANCEEvent->CustomWeight = NOVAw;
fNUISANCEEvent->CustomWeightArray[0] = MAQEw;
fNUISANCEEvent->CustomWeightArray[1] = NonResw;
fNUISANCEEvent->CustomWeightArray[2] = RPAQEw;
fNUISANCEEvent->CustomWeightArray[3] = RPARESw;
fNUISANCEEvent->CustomWeightArray[4] = MECw;
fNUISANCEEvent->CustomWeightArray[5] = NOVAw;
} else {
fNUISANCEEvent->CustomWeight = 1.0;
fNUISANCEEvent->CustomWeightArray[0] = 1.0;
fNUISANCEEvent->CustomWeightArray[1] = 1.0;
fNUISANCEEvent->CustomWeightArray[2] = 1.0;
fNUISANCEEvent->CustomWeightArray[3] = 1.0;
fNUISANCEEvent->CustomWeightArray[4] = 1.0;
fNUISANCEEvent->CustomWeightArray[5] = 1.0;
}
// Get N Particle Stack
unsigned int npart = fGenieGHep->GetEntries();
unsigned int kmax = fNUISANCEEvent->kMaxParticles;
if (npart > kmax) {
ERR(WRN) << "GENIE has too many particles, expanding stack." << std::endl;
fNUISANCEEvent->ExpandParticleStack(npart);
}
// Fill Particle Stack
GHepParticle* p = 0;
TObjArrayIter iter(fGenieGHep);
fNUISANCEEvent->fNParticles = 0;
// Loop over all particles
while ((p = (dynamic_cast<genie::GHepParticle*>((iter).Next())))) {
if (!p) continue;
// Get Status
int state = GetGENIEParticleStatus(p, fNUISANCEEvent->Mode);
// Remove Undefined
if (kRemoveUndefParticles && state == kUndefinedState) continue;
// Remove FSI
if (kRemoveFSIParticles && state == kFSIState) continue;
if (kRemoveNuclearParticles &&
(state == kNuclearInitial || state == kNuclearRemnant))
continue;
// Fill Vectors
int curpart = fNUISANCEEvent->fNParticles;
fNUISANCEEvent->fParticleState[curpart] = state;
// Mom
fNUISANCEEvent->fParticleMom[curpart][0] = p->Px() * 1.E3;
fNUISANCEEvent->fParticleMom[curpart][1] = p->Py() * 1.E3;
fNUISANCEEvent->fParticleMom[curpart][2] = p->Pz() * 1.E3;
fNUISANCEEvent->fParticleMom[curpart][3] = p->E() * 1.E3;
// PDG
fNUISANCEEvent->fParticlePDG[curpart] = p->Pdg();
+ // Set if the particle was on the fundamental vertex
+ fNUISANCEEvent->fPrimaryVertex[curpart] = (p->FirstMother() < 2);
+
// Add to N particle count
fNUISANCEEvent->fNParticles++;
// Extra Check incase GENIE fails.
if ((UInt_t)fNUISANCEEvent->fNParticles == kmax) {
ERR(WRN) << "Number of GENIE Particles exceeds maximum!" << std::endl;
ERR(WRN) << "Extend kMax, or run without including FSI particles!"
<< std::endl;
break;
}
}
// Fill Extra Stack
if (fSaveExtra) fGenieInfo->FillGeneratorInfo(fGenieNtpl);
// Run Initial, FSI, Final, Other ordering.
fNUISANCEEvent->OrderStack();
FitParticle* ISNeutralLepton =
fNUISANCEEvent->GetHMISParticle(PhysConst::pdg_neutrinos);
if (ISNeutralLepton) {
fNUISANCEEvent->probe_E = ISNeutralLepton->E();
fNUISANCEEvent->probe_pdg = ISNeutralLepton->PDG();
}
return;
}
void GENIEInputHandler::Print() {}
#endif
diff --git a/src/InputHandler/NEUTInputHandler.cxx b/src/InputHandler/NEUTInputHandler.cxx
index 7df81d1..b3fec8f 100644
--- a/src/InputHandler/NEUTInputHandler.cxx
+++ b/src/InputHandler/NEUTInputHandler.cxx
@@ -1,503 +1,510 @@
#ifdef __NEUT_ENABLED__
#include "NEUTInputHandler.h"
#include "InputUtils.h"
NEUTGeneratorInfo::~NEUTGeneratorInfo() { DeallocateParticleStack(); }
void NEUTGeneratorInfo::AddBranchesToTree(TTree* tn) {
tn->Branch("NEUTParticleN", fNEUTParticleN, "NEUTParticleN/I");
tn->Branch("NEUTParticleStatusCode", fNEUTParticleStatusCode,
"NEUTParticleStatusCode[NEUTParticleN]/I");
tn->Branch("NEUTParticleAliveCode", fNEUTParticleAliveCode,
"NEUTParticleAliveCode[NEUTParticleN]/I");
}
void NEUTGeneratorInfo::SetBranchesFromTree(TTree* tn) {
tn->SetBranchAddress("NEUTParticleN", &fNEUTParticleN);
tn->SetBranchAddress("NEUTParticleStatusCode", &fNEUTParticleStatusCode);
tn->SetBranchAddress("NEUTParticleAliveCode", &fNEUTParticleAliveCode);
}
void NEUTGeneratorInfo::AllocateParticleStack(int stacksize) {
fNEUTParticleN = 0;
fNEUTParticleStatusCode = new int[stacksize];
fNEUTParticleStatusCode = new int[stacksize];
}
void NEUTGeneratorInfo::DeallocateParticleStack() {
delete fNEUTParticleStatusCode;
delete fNEUTParticleAliveCode;
}
void NEUTGeneratorInfo::FillGeneratorInfo(NeutVect* nevent) {
Reset();
for (int i = 0; i < nevent->Npart(); i++) {
fNEUTParticleStatusCode[i] = nevent->PartInfo(i)->fStatus;
fNEUTParticleAliveCode[i] = nevent->PartInfo(i)->fIsAlive;
fNEUTParticleN++;
}
}
void NEUTGeneratorInfo::Reset() {
for (int i = 0; i < fNEUTParticleN; i++) {
fNEUTParticleStatusCode[i] = -1;
fNEUTParticleAliveCode[i] = 9;
}
fNEUTParticleN = 0;
}
NEUTInputHandler::NEUTInputHandler(std::string const& handle,
std::string const& rawinputs) {
LOG(SAM) << "Creating NEUTInputHandler : " << handle << std::endl;
// Run a joint input handling
fName = handle;
// Setup the TChain
fNEUTTree = new TChain("neuttree");
fSaveExtra = FitPar::Config().GetParB("SaveExtraNEUT");
fCacheSize = FitPar::Config().GetParI("CacheSize");
fMaxEvents = FitPar::Config().GetParI("MAXEVENTS");
// Loop over all inputs and grab flux, eventhist, and nevents
std::vector<std::string> inputs = InputUtils::ParseInputFileList(rawinputs);
for (size_t inp_it = 0; inp_it < inputs.size(); ++inp_it) {
// Open File for histogram access
TFile* inp_file = new TFile(inputs[inp_it].c_str(), "READ");
if (!inp_file or inp_file->IsZombie()) {
THROW("NEUT File IsZombie() at : '"
<< inputs[inp_it] << "'" << std::endl
<< "Check that your file paths are correct and the file exists!"
<< std::endl
<< "$ ls -lh " << inputs[inp_it]);
}
// Get Flux/Event hist
TH1D* fluxhist = (TH1D*)inp_file->Get(
(PlotUtils::GetObjectWithName(inp_file, "flux")).c_str());
TH1D* eventhist = (TH1D*)inp_file->Get(
(PlotUtils::GetObjectWithName(inp_file, "evt")).c_str());
if (!fluxhist or !eventhist) {
ERROR(FTL, "Input File Contents: " << inputs[inp_it]);
inp_file->ls();
THROW(
"NEUT FILE doesn't contain flux/xsec info. You may have to "
"regenerate your MC!");
}
// Get N Events
TTree* neuttree = (TTree*)inp_file->Get("neuttree");
if (!neuttree) {
ERROR(FTL, "neuttree not located in NEUT file: " << inputs[inp_it]);
THROW("Check your inputs, they may need to be completely regenerated!");
throw;
}
int nevents = neuttree->GetEntries();
if (nevents <= 0) {
THROW("Trying to a TTree with "
<< nevents << " to TChain from : " << inputs[inp_it]);
}
// Register input to form flux/event rate hists
RegisterJointInput(inputs[inp_it], nevents, fluxhist, eventhist);
// Add To TChain
fNEUTTree->AddFile(inputs[inp_it].c_str());
}
// Registor all our file inputs
SetupJointInputs();
// Assign to tree
fEventType = kNEUT;
fNeutVect = NULL;
fNEUTTree->SetBranchAddress("vectorbranch", &fNeutVect);
fNEUTTree->GetEntry(0);
// Create Fit Event
fNUISANCEEvent = new FitEvent();
fNUISANCEEvent->SetNeutVect(fNeutVect);
if (fSaveExtra) {
fNeutInfo = new NEUTGeneratorInfo();
fNUISANCEEvent->AddGeneratorInfo(fNeutInfo);
}
fNUISANCEEvent->HardReset();
};
NEUTInputHandler::~NEUTInputHandler(){
// if (fNEUTTree) delete fNEUTTree;
// if (fNeutVect) delete fNeutVect;
// if (fNeutInfo) delete fNeutInfo;
};
void NEUTInputHandler::CreateCache() {
if (fCacheSize > 0) {
// fNEUTTree->SetCacheEntryRange(0, fNEvents);
fNEUTTree->AddBranchToCache("vectorbranch", 1);
fNEUTTree->SetCacheSize(fCacheSize);
}
}
void NEUTInputHandler::RemoveCache() {
// fNEUTTree->SetCacheEntryRange(0, fNEvents);
fNEUTTree->AddBranchToCache("vectorbranch", 0);
fNEUTTree->SetCacheSize(0);
}
FitEvent* NEUTInputHandler::GetNuisanceEvent(const UInt_t entry,
const bool lightweight) {
// Catch too large entries
if (entry >= (UInt_t)fNEvents) return NULL;
// Read Entry from TTree to fill NEUT Vect in BaseFitEvt;
fNEUTTree->GetEntry(entry);
// Run NUISANCE Vector Filler
if (!lightweight) {
CalcNUISANCEKinematics();
}
#ifdef __PROB3PP_ENABLED__
else {
UInt_t npart = fNeutVect->Npart();
for (size_t i = 0; i < npart; i++) {
NeutPart* part = fNUISANCEEvent->fNeutVect->PartInfo(i);
if ((part->fIsAlive == false) && (part->fStatus == -1) &&
std::count(PhysConst::pdg_neutrinos, PhysConst::pdg_neutrinos + 4,
part->fPID)) {
fNUISANCEEvent->probe_E = part->fP.T();
fNUISANCEEvent->probe_pdg = part->fPID;
break;
} else {
continue;
}
}
}
#endif
// Setup Input scaling for joint inputs
fNUISANCEEvent->InputWeight = GetInputWeight(entry);
// Return event pointer
return fNUISANCEEvent;
}
// From NEUT neutclass/neutpart.h
// Bool_t fIsAlive; // Particle should be tracked or not
// ( in the detector simulator )
//
// Int_t fStatus; // Status flag of this particle
// -2: Non existing particle
// -1: Initial state particle
// 0: Normal
// 1: Decayed to the other particle
// 2: Escaped from the detector
// 3: Absorped
// 4: Charge exchanged
// 5: Pauli blocked
// 6: N/A
// 7: Produced child particles
// 8: Inelastically scattered
//
int NEUTInputHandler::GetNeutParticleStatus(NeutPart* part) {
// State
int state = kUndefinedState;
// Remove Pauli blocked events, probably just single pion events
if (part->fStatus == 5) {
state = kFSIState;
// fStatus == -1 means initial state
} else if (part->fIsAlive == false && part->fStatus == -1) {
state = kInitialState;
// NEUT has a bit of a strange convention for fIsAlive and fStatus
// combinations
// for NC and neutrino particle isAlive true/false and status 2 means
// final state particle
// for other particles in NC status 2 means it's an FSI particle
// for CC it means it was an FSI particle
} else if (part->fStatus == 2) {
// NC case is a little strange... The outgoing neutrino might be alive or
// not alive. Remaining particles with status 2 are FSI particles that
// reinteracted
if (abs(fNeutVect->Mode) > 30 &&
(abs(part->fPID) == 16 || abs(part->fPID) == 14 || abs(part->fPID) == 12)) {
state = kFinalState;
// The usual CC case
} else if (part->fIsAlive == true) {
state = kFSIState;
}
} else if (part->fIsAlive == true && part->fStatus == 2 &&
(abs(part->fPID) == 16 || abs(part->fPID) == 14 || abs(part->fPID) == 12)) {
state = kFinalState;
} else if (part->fIsAlive == true && part->fStatus == 0) {
state = kFinalState;
} else if (!part->fIsAlive && (part->fStatus == 1 || part->fStatus == 3 || part->fStatus == 4 || part->fStatus == 7 || part->fStatus == 8)) {
state = kFSIState;
// There's one hyper weird case where fStatus = -3. This apparently corresponds to a nucleon being ejected via pion FSI when there is "data available"
} else if (!part->fIsAlive && (part->fStatus == -3)) {
state = kUndefinedState;
// NC neutrino outgoing
} else if (!part->fIsAlive && part->fStatus == 0 && (abs(part->fPID) == 16 || abs(part->fPID) == 14 || abs(part->fPID) == 12)) {
state = kFinalState;
// Warn if we still find alive particles without classifying them
} else if (part->fIsAlive == true) {
ERR(WRN) << "Undefined NEUT state "
<< " Alive: " << part->fIsAlive << " Status: " << part->fStatus
<< " PDG: " << part->fPID << std::endl;
throw;
// Warn if we find dead particles that we haven't classified
} else {
ERR(WRN) << "Undefined NEUT state "
<< " Alive: " << part->fIsAlive << " Status: " << part->fStatus
<< " PDG: " << part->fPID << std::endl;
throw;
}
return state;
}
void NEUTInputHandler::CalcNUISANCEKinematics() {
// Reset all variables
fNUISANCEEvent->ResetEvent();
// Fill Globals
fNUISANCEEvent->Mode = fNeutVect->Mode;
fNUISANCEEvent->fEventNo = fNeutVect->EventNo;
fNUISANCEEvent->fTargetA = fNeutVect->TargetA;
fNUISANCEEvent->fTargetZ = fNeutVect->TargetZ;
fNUISANCEEvent->fTargetH = fNeutVect->TargetH;
fNUISANCEEvent->fBound = bool(fNeutVect->Ibound);
if (fNUISANCEEvent->fBound) {
fNUISANCEEvent->fTargetPDG = TargetUtils::GetTargetPDGFromZA(
fNUISANCEEvent->fTargetZ, fNUISANCEEvent->fTargetA);
} else {
fNUISANCEEvent->fTargetPDG = 1000010010;
}
// Check Particle Stack
UInt_t npart = fNeutVect->Npart();
UInt_t kmax = fNUISANCEEvent->kMaxParticles;
if (npart > kmax) {
ERR(FTL) << "NEUT has too many particles. Expanding stack." << std::endl;
fNUISANCEEvent->ExpandParticleStack(npart);
throw;
}
+ int nprimary = fNeutVect->Nprimary();
// Fill Particle Stack
for (size_t i = 0; i < npart; i++) {
// Get Current Count
int curpart = fNUISANCEEvent->fNParticles;
// Get NEUT Particle
NeutPart* part = fNeutVect->PartInfo(i);
// State
int state = GetNeutParticleStatus(part);
// Remove Undefined
if (kRemoveUndefParticles && state == kUndefinedState) continue;
// Remove FSI
if (kRemoveFSIParticles && state == kFSIState) continue;
// Remove Nuclear
if (kRemoveNuclearParticles &&
(state == kNuclearInitial || state == kNuclearRemnant))
continue;
// State
fNUISANCEEvent->fParticleState[curpart] = state;
+ // Is the paricle associated with the primary vertex?
+ bool primary = false;
+ // NEUT events are just popped onto the stack as primary, then continues to be non-primary
+ if (i < nprimary) primary = true;
+ fNUISANCEEvent->fPrimaryVertex[curpart] = primary;
+
// Mom
fNUISANCEEvent->fParticleMom[curpart][0] = part->fP.X();
fNUISANCEEvent->fParticleMom[curpart][1] = part->fP.Y();
fNUISANCEEvent->fParticleMom[curpart][2] = part->fP.Z();
fNUISANCEEvent->fParticleMom[curpart][3] = part->fP.T();
// PDG
fNUISANCEEvent->fParticlePDG[curpart] = part->fPID;
// Add up particle count
fNUISANCEEvent->fNParticles++;
}
// Save Extra Generator Info
if (fSaveExtra) {
fNeutInfo->FillGeneratorInfo(fNeutVect);
}
// Run Initial, FSI, Final, Other ordering.
fNUISANCEEvent->OrderStack();
FitParticle* ISNeutralLepton =
fNUISANCEEvent->GetHMISParticle(PhysConst::pdg_neutrinos);
if (ISNeutralLepton) {
fNUISANCEEvent->probe_E = ISNeutralLepton->E();
fNUISANCEEvent->probe_pdg = ISNeutralLepton->PDG();
}
return;
}
void NEUTUtils::FillNeutCommons(NeutVect* nvect) {
// WARNING: This has only been implemented for a neuttree and not GENIE
// This should be kept in sync with T2KNIWGUtils::GetNIWGEvent(TTree)
// NEUT version info. Can't get it to compile properly with this yet
// neutversion_.corev = nvect->COREVer;
// neutversion_.nucev = nvect->NUCEVer;
// neutversion_.nuccv = nvect->NUCCVer;
// Documentation: See nework.h
nework_.modene = nvect->Mode;
nework_.numne = nvect->Npart();
#ifdef NEUT_COMMON_QEAV
nemdls_.mdlqeaf = nvect->QEAVForm;
#else
nemdls_.mdlqeaf = nvect->QEVForm;
#endif
nemdls_.mdlqe = nvect->QEModel;
nemdls_.mdlspi = nvect->SPIModel;
nemdls_.mdldis = nvect->DISModel;
nemdls_.mdlcoh = nvect->COHModel;
neutcoh_.necohepi = nvect->COHModel;
nemdls_.xmaqe = nvect->QEMA;
nemdls_.xmvqe = nvect->QEMV;
nemdls_.kapp = nvect->KAPPA;
// nemdls_.sccfv = SCCFVdef;
// nemdls_.sccfa = SCCFAdef;
// nemdls_.fpqe = FPQEdef;
nemdls_.xmaspi = nvect->SPIMA;
nemdls_.xmvspi = nvect->SPIMV;
nemdls_.xmares = nvect->RESMA;
nemdls_.xmvres = nvect->RESMV;
neut1pi_.xmanffres = nvect->SPIMA;
neut1pi_.xmvnffres = nvect->SPIMV;
neut1pi_.xmarsres = nvect->RESMA;
neut1pi_.xmvrsres = nvect->RESMV;
neut1pi_.neiff = nvect->SPIForm;
neut1pi_.nenrtype = nvect->SPINRType;
neut1pi_.rneca5i = nvect->SPICA5I;
neut1pi_.rnebgscl = nvect->SPIBGScale;
nemdls_.xmacoh = nvect->COHMA;
nemdls_.rad0nu = nvect->COHR0;
// nemdls_.fa1coh = nvect->COHA1err;
// nemdls_.fb1coh = nvect->COHb1err;
// neutdis_.nepdf = NEPDFdef;
// neutdis_.nebodek = NEBODEKdef;
neutcard_.nefrmflg = nvect->FrmFlg;
neutcard_.nepauflg = nvect->PauFlg;
neutcard_.nenefo16 = nvect->NefO16;
neutcard_.nemodflg = nvect->ModFlg;
// neutcard_.nenefmodl = 1;
// neutcard_.nenefmodh = 1;
// neutcard_.nenefkinh = 1;
// neutpiabs_.neabspiemit = 1;
nenupr_.iformlen = nvect->FormLen;
neutpiless_.ipilessdcy = nvect->IPilessDcy;
neutpiless_.rpilessdcy = nvect->RPilessDcy;
neutpiless_.ipilessdcy = nvect->IPilessDcy;
neutpiless_.rpilessdcy = nvect->RPilessDcy;
neffpr_.fefqe = nvect->NuceffFactorPIQE;
neffpr_.fefqeh = nvect->NuceffFactorPIQEH;
neffpr_.fefinel = nvect->NuceffFactorPIInel;
neffpr_.fefabs = nvect->NuceffFactorPIAbs;
neffpr_.fefcx = nvect->NuceffFactorPICX;
neffpr_.fefcxh = nvect->NuceffFactorPICXH;
neffpr_.fefcoh = nvect->NuceffFactorPICoh;
neffpr_.fefqehf = nvect->NuceffFactorPIQEHKin;
neffpr_.fefcxhf = nvect->NuceffFactorPICXKin;
neffpr_.fefcohf = nvect->NuceffFactorPIQELKin;
for (int i = 0; i < nework_.numne; i++) {
nework_.ipne[i] = nvect->PartInfo(i)->fPID;
nework_.pne[i][0] =
(float)nvect->PartInfo(i)->fP.X() / 1000; // VC(NE)WORK in M(G)eV
nework_.pne[i][1] =
(float)nvect->PartInfo(i)->fP.Y() / 1000; // VC(NE)WORK in M(G)eV
nework_.pne[i][2] =
(float)nvect->PartInfo(i)->fP.Z() / 1000; // VC(NE)WORK in M(G)eV
}
// fsihist.h
// neutroot fills a dummy object for events with no FSI to prevent memory leak
// when
// reading the TTree, so check for it here
if ((int)nvect->NfsiVert() ==
1) { // An event with FSI must have at least two vertices
// if (nvect->NfsiPart()!=1 || nvect->Fsiprob!=-1)
// ERR(WRN) << "T2KNeutUtils::fill_neut_commons(TTree) NfsiPart!=1 or
// Fsiprob!=-1 when NfsiVert==1" << std::endl;
fsihist_.nvert = 0;
fsihist_.nvcvert = 0;
fsihist_.fsiprob = 1;
} else { // Real FSI event
fsihist_.nvert = (int)nvect->NfsiVert();
for (int ivert = 0; ivert < fsihist_.nvert; ivert++) {
fsihist_.iflgvert[ivert] = nvect->FsiVertInfo(ivert)->fVertID;
fsihist_.posvert[ivert][0] = (float)nvect->FsiVertInfo(ivert)->fPos.X();
fsihist_.posvert[ivert][1] = (float)nvect->FsiVertInfo(ivert)->fPos.Y();
fsihist_.posvert[ivert][2] = (float)nvect->FsiVertInfo(ivert)->fPos.Z();
}
fsihist_.nvcvert = nvect->NfsiPart();
for (int ip = 0; ip < fsihist_.nvcvert; ip++) {
fsihist_.abspvert[ip] = (float)nvect->FsiPartInfo(ip)->fMomLab;
fsihist_.abstpvert[ip] = (float)nvect->FsiPartInfo(ip)->fMomNuc;
fsihist_.ipvert[ip] = nvect->FsiPartInfo(ip)->fPID;
fsihist_.iverti[ip] = nvect->FsiPartInfo(ip)->fVertStart;
fsihist_.ivertf[ip] = nvect->FsiPartInfo(ip)->fVertEnd;
fsihist_.dirvert[ip][0] = (float)nvect->FsiPartInfo(ip)->fDir.X();
fsihist_.dirvert[ip][1] = (float)nvect->FsiPartInfo(ip)->fDir.Y();
fsihist_.dirvert[ip][2] = (float)nvect->FsiPartInfo(ip)->fDir.Z();
}
fsihist_.fsiprob = nvect->Fsiprob;
}
neutcrscom_.crsx = nvect->Crsx;
neutcrscom_.crsy = nvect->Crsy;
neutcrscom_.crsz = nvect->Crsz;
neutcrscom_.crsphi = nvect->Crsphi;
neutcrscom_.crsq2 = nvect->Crsq2;
neuttarget_.numbndn = nvect->TargetA - nvect->TargetZ;
neuttarget_.numbndp = nvect->TargetZ;
neuttarget_.numfrep = nvect->TargetH;
neuttarget_.numatom = nvect->TargetA;
posinnuc_.ibound = nvect->Ibound;
// put empty nucleon FSI history (since it is not saved in the NeutVect
// format)
// Comment out as NEUT does not have the necessary proton FSI information yet
// nucleonfsihist_.nfnvert = 0;
// nucleonfsihist_.nfnstep = 0;
}
#endif
diff --git a/src/InputHandler/NuWroInputHandler.cxx b/src/InputHandler/NuWroInputHandler.cxx
index 9c7bf05..1853d95 100644
--- a/src/InputHandler/NuWroInputHandler.cxx
+++ b/src/InputHandler/NuWroInputHandler.cxx
@@ -1,480 +1,507 @@
#ifdef __NUWRO_ENABLED__
#include "NuWroInputHandler.h"
#include "InputUtils.h"
NuWroGeneratorInfo::~NuWroGeneratorInfo() { delete fNuWroParticlePDGs; }
void NuWroGeneratorInfo::AddBranchesToTree(TTree* tn) {
tn->Branch("NuWroParticlePDGs", &fNuWroParticlePDGs, "NuWroParticlePDGs/I");
}
void NuWroGeneratorInfo::SetBranchesFromTree(TTree* tn) {
tn->SetBranchAddress("NuWroParticlePDGs", &fNuWroParticlePDGs);
}
void NuWroGeneratorInfo::AllocateParticleStack(int stacksize) {
fNuWroParticlePDGs = new int[stacksize];
}
void NuWroGeneratorInfo::DeallocateParticleStack() {
delete fNuWroParticlePDGs;
}
void NuWroGeneratorInfo::FillGeneratorInfo(event* e) { Reset(); }
void NuWroGeneratorInfo::Reset() {
for (int i = 0; i < kMaxParticles; i++) {
fNuWroParticlePDGs[i] = 0;
}
}
int event1_nof(event* e, int pdg) {
int c = 0;
for (size_t i = 0; i < e->out.size(); i++)
if (e->out[i].pdg == pdg) c++;
return c;
}
NuWroInputHandler::NuWroInputHandler(std::string const& handle,
std::string const& rawinputs) {
LOG(SAM) << "Creating NuWroInputHandler : " << handle << std::endl;
// Run a joint input handling
fName = handle;
fMaxEvents = FitPar::Config().GetParI("MAXEVENTS");
fSaveExtra = false; // FitPar::Config().GetParB("NuWroSaveExtra");
// Setup the TChain
fNuWroTree = new TChain("treeout");
// Loop over all inputs and grab flux, eventhist, and nevents
std::vector<std::string> inputs = InputUtils::ParseInputFileList(rawinputs);
for (size_t inp_it = 0; inp_it < inputs.size(); ++inp_it) {
// Open File for histogram access
TFile* inp_file = new TFile(inputs[inp_it].c_str(), "READ");
if (!inp_file or inp_file->IsZombie()) {
ERR(FTL) << "nuwro File IsZombie() at " << inputs[inp_it] << std::endl;
throw;
}
// Get Flux/Event hist
TH1D* fluxhist = (TH1D*)inp_file->Get(
(PlotUtils::GetObjectWithName(inp_file, "FluxHist")).c_str());
TH1D* eventhist = (TH1D*)inp_file->Get(
(PlotUtils::GetObjectWithName(inp_file, "EvtHist")).c_str());
if (!fluxhist or !eventhist) {
ERR(FTL) << "nuwro FILE doesn't contain flux/xsec info" << std::endl;
if (FitPar::Config().GetParB("regennuwro")) {
ERR(FTL) << "Regen NuWro has not been added yet. Email the developers!"
<< std::endl;
// ProcessNuWroInputFlux(inputs[inp_it]);
throw;
} else {
ERR(FTL) << "If you would like NUISANCE to generate these for you "
<< "please set parameter regennuwro=1 and re-run."
<< std::endl;
throw;
}
}
// Get N Events
TTree* nuwrotree = (TTree*)inp_file->Get("treeout");
if (!nuwrotree) {
ERR(FTL) << "treeout not located in nuwro file! " << inputs[inp_it]
<< std::endl;
throw;
}
int nevents = nuwrotree->GetEntries();
// Register input to form flux/event rate hists
RegisterJointInput(inputs[inp_it], nevents, fluxhist, eventhist);
// Add to TChain
fNuWroTree->Add(inputs[inp_it].c_str());
}
// Registor all our file inputs
SetupJointInputs();
// Setup Events
fNuWroEvent = NULL;
fNuWroTree->SetBranchAddress("e", &fNuWroEvent);
fNuWroTree->GetEntry(0);
fNUISANCEEvent = new FitEvent();
fNUISANCEEvent->fType = kNUWRO;
fNUISANCEEvent->fNuwroEvent = fNuWroEvent;
fNUISANCEEvent->HardReset();
if (fSaveExtra) {
fNuWroInfo = new NuWroGeneratorInfo();
fNUISANCEEvent->AddGeneratorInfo(fNuWroInfo);
}
};
NuWroInputHandler::~NuWroInputHandler() {
if (fNuWroTree) delete fNuWroTree;
}
void NuWroInputHandler::CreateCache() {
// fNuWroTree->SetCacheEntryRange(0, fNEvents);
// fNuWroTree->AddBranchToCache("*", 1);
// fNuWroTree->SetCacheSize(fCacheSize);
}
void NuWroInputHandler::RemoveCache() {
// fNuWroTree->SetCacheEntryRange(0, fNEvents);
// fNuWroTree->AddBranchToCache("*", 0);
// fNuWroTree->SetCacheSize(0);
}
void NuWroInputHandler::ProcessNuWroInputFlux(const std::string file) {}
FitEvent* NuWroInputHandler::GetNuisanceEvent(const UInt_t entry,
const bool lightweight) {
// Catch too large entries
if (entry >= (UInt_t)fNEvents) return NULL;
// Read Entry from TTree to fill NEUT Vect in BaseFitEvt;
fNuWroTree->GetEntry(entry);
// Run NUISANCE Vector Filler
if (!lightweight) {
CalcNUISANCEKinematics();
}
#ifdef __PROB3PP_ENABLED__
for (size_t i = 0; i < fNUISANCEEvent->fNuwroEvent->in.size(); i++) {
if (std::count(PhysConst::pdg_neutrinos, PhysConst::pdg_neutrinos + 4,
fNUISANCEEvent->fNuwroEvent->in[i].pdg)) {
fNUISANCEEvent->probe_E = fNUISANCEEvent->fNuwroEvent->in[i].t;
fNUISANCEEvent->probe_pdg = fNUISANCEEvent->fNuwroEvent->in[i].pdg;
break;
}
}
#endif
// Setup Input scaling for joint inputs
fNUISANCEEvent->InputWeight = GetInputWeight(entry);
#ifdef __USE_NUWRO_SRW_EVENTS__
if (!rwEvs.size()) {
fNuwroParams = fNuWroEvent->par;
}
if (entry >= rwEvs.size()) {
rwEvs.push_back(BaseFitEvt());
rwEvs.back().fType = kNUWRO;
rwEvs.back().Mode = fNUISANCEEvent->Mode;
rwEvs.back().fNuwroSRWEvent = SRW::SRWEvent(*fNuWroEvent);
rwEvs.back().fNuwroEvent = NULL;
rwEvs.back().fNuwroParams = &fNuwroParams;
rwEvs.back().probe_E = rwEvs.back().fNuwroSRWEvent.NeutrinoEnergy;
rwEvs.back().probe_pdg = rwEvs.back().fNuwroSRWEvent.NeutrinoPDG;
}
fNUISANCEEvent->fNuwroSRWEvent = SRW::SRWEvent(*fNuWroEvent);
fNUISANCEEvent->fNuwroParams = &fNuwroParams;
fNUISANCEEvent->probe_E = fNUISANCEEvent->fNuwroSRWEvent.NeutrinoEnergy;
fNUISANCEEvent->probe_pdg = fNUISANCEEvent->fNuwroSRWEvent.NeutrinoPDG;
#endif
return fNUISANCEEvent;
}
int NuWroInputHandler::ConvertNuwroMode(event* e) {
Int_t proton_pdg, neutron_pdg, pion_pdg, pion_plus_pdg, pion_minus_pdg,
lambda_pdg, eta_pdg, kaon_pdg, kaon_plus_pdg;
proton_pdg = 2212;
eta_pdg = 221;
neutron_pdg = 2112;
pion_pdg = 111;
pion_plus_pdg = 211;
pion_minus_pdg = -211;
// O_16_pdg = 100069; // oznacznie z Neuta
lambda_pdg = 3122;
kaon_pdg = 311;
kaon_plus_pdg = 321;
if (e->flag.qel) // kwiazielastyczne oddziaływanie
{
if (e->flag.anty) // jeśli jest to oddziaływanie z antyneutrinem
{
if (e->flag.cc)
return -1;
else {
if (event1_nof(e, proton_pdg))
return -51;
else if (event1_nof(e, neutron_pdg))
return -52; // sprawdzam dodatkowo ?
}
} else // oddziaływanie z neutrinem
{
if (e->flag.cc)
return 1;
else {
if (event1_nof(e, proton_pdg))
return 51;
else if (event1_nof(e, neutron_pdg))
return 52;
}
}
}
if (e->flag.mec) {
if (e->flag.anty)
return -2;
else
return 2;
}
if (e->flag.res) // rezonansowa produkcja: pojedynczy pion, pojed.eta, kaon,
// multipiony
{
Int_t liczba_pionow, liczba_kaonow;
liczba_pionow = event1_nof(e, pion_pdg) + event1_nof(e, pion_plus_pdg) +
event1_nof(e, pion_minus_pdg);
liczba_kaonow = event1_nof(e, kaon_pdg) + event1_nof(e, kaon_pdg);
if (liczba_pionow > 1 || liczba_pionow == 0) // multipiony
{
if (e->flag.anty) {
if (e->flag.cc)
return -21;
else
return -41;
} else {
if (e->flag.cc)
return 21;
else
return 41;
}
}
if (liczba_pionow == 1) {
if (e->flag.anty) // jeśli jest to oddziaływanie z antyneutrinem
{
if (e->flag.cc) {
if (event1_nof(e, neutron_pdg) && event1_nof(e, pion_minus_pdg))
return -11;
if (event1_nof(e, neutron_pdg) && event1_nof(e, pion_pdg)) return -12;
if (event1_nof(e, proton_pdg) && event1_nof(e, pion_minus_pdg))
return -13;
} else {
if (event1_nof(e, proton_pdg)) {
if (event1_nof(e, pion_minus_pdg))
return -33;
else if (event1_nof(e, pion_pdg))
return -32;
} else if (event1_nof(e, neutron_pdg)) {
if (event1_nof(e, pion_plus_pdg))
return -34;
else if (event1_nof(e, pion_pdg))
return -31;
}
}
} else // oddziaływanie z neutrinem
{
if (e->flag.cc) {
if (event1_nof(e, proton_pdg) && event1_nof(e, pion_plus_pdg))
return 11;
if (event1_nof(e, proton_pdg) && event1_nof(e, pion_pdg)) return 12;
if (event1_nof(e, neutron_pdg) && event1_nof(e, pion_plus_pdg))
return 13;
} else {
if (event1_nof(e, proton_pdg)) {
if (event1_nof(e, pion_minus_pdg))
return 33;
else if (event1_nof(e, pion_pdg))
return 32;
} else if (event1_nof(e, neutron_pdg)) {
if (event1_nof(e, pion_plus_pdg))
return 34;
else if (event1_nof(e, pion_pdg))
return 31;
}
}
}
}
if (event1_nof(e, eta_pdg)) // produkcja rezonansowa ety
{
if (e->flag.anty) // jeśli jest to oddziaływanie z antyneutrinem
{
if (e->flag.cc)
return -22;
else {
if (event1_nof(e, neutron_pdg))
return -42;
else if (event1_nof(e, proton_pdg))
return -43; // sprawdzam dodatkowo ?
}
} else // oddziaływanie z neutrinem
{
if (e->flag.cc)
return 22;
else {
if (event1_nof(e, neutron_pdg))
return 42;
else if (event1_nof(e, proton_pdg))
return 43;
}
}
}
if (event1_nof(e, lambda_pdg) == 1 &&
liczba_kaonow == 1) // produkcja rezonansowa kaonu
{
if (e->flag.anty) // jeśli jest to oddziaływanie z antyneutrinem
{
if (e->flag.cc && event1_nof(e, kaon_pdg))
return -23;
else {
if (event1_nof(e, kaon_pdg))
return -44;
else if (event1_nof(e, kaon_plus_pdg))
return -45;
}
} else // oddziaływanie z neutrinem
{
if (e->flag.cc && event1_nof(e, kaon_plus_pdg))
return 23;
else {
if (event1_nof(e, kaon_pdg))
return 44;
else if (event1_nof(e, kaon_plus_pdg))
return 45;
}
}
}
}
if (e->flag.coh) // koherentne oddziaływanie tylko na O(16)
{
Int_t _target;
_target = e->par.nucleus_p + e->par.nucleus_n; // liczba masowa O(16)
if (_target == 16) {
if (e->flag.anty) // jeśli jest to oddziaływanie z antyneutrinem
{
if (e->flag.cc && event1_nof(e, pion_minus_pdg))
return -16;
else if (event1_nof(e, pion_pdg))
return -36;
} else // oddziaływanie z neutrinem
{
if (e->flag.cc && event1_nof(e, pion_plus_pdg))
return 16;
else if (event1_nof(e, pion_pdg))
return 36;
}
}
}
// gleboko nieelastyczne rozpraszanie
if (e->flag.dis) {
if (e->flag.anty) {
if (e->flag.cc)
return -26;
else
return -46;
} else {
if (e->flag.cc)
return 26;
else
return 46;
}
}
return 9999;
}
void NuWroInputHandler::CalcNUISANCEKinematics() {
// std::cout << "NuWro Event Address " << fNuWroEvent << std::endl;
// Reset all variables
fNUISANCEEvent->ResetEvent();
FitEvent* evt = fNUISANCEEvent;
// Sort Event Info
evt->Mode = ConvertNuwroMode(fNuWroEvent);
if (abs(evt->Mode) > 60) {
evt->Mode = 0;
}
evt->fEventNo = 0.0;
evt->fTotCrs = 0.0;
evt->fTargetA = fNuWroEvent->par.nucleus_p + fNuWroEvent->par.nucleus_n;
evt->fTargetZ = fNuWroEvent->par.nucleus_p;
evt->fTargetPDG = TargetUtils::GetTargetPDGFromZA(evt->fTargetZ, evt->fTargetA);
evt->fTargetH = 0;
evt->fBound = (evt->fTargetA != 1);
// Check Particle Stack
UInt_t npart_in = fNuWroEvent->in.size();
UInt_t npart_out = fNuWroEvent->out.size();
UInt_t npart_post = fNuWroEvent->post.size();
UInt_t npart = npart_in + npart_out + npart_post;
UInt_t kmax = evt->kMaxParticles;
if (npart > kmax) {
ERR(WRN) << "NUWRO has too many particles. Expanding stack." << std::endl;
fNUISANCEEvent->ExpandParticleStack(npart);
}
- // Sort Particles
evt->fNParticles = 0;
std::vector<particle>::iterator p_iter;
- // Initial State
- for (p_iter = fNuWroEvent->in.begin(); p_iter != fNuWroEvent->in.end();
- p_iter++) {
- AddNuWroParticle(fNUISANCEEvent, (*p_iter), kInitialState);
+ // Get the Initial State
+ for (p_iter = fNuWroEvent->in.begin(); p_iter != fNuWroEvent->in.end(); p_iter++) {
+ AddNuWroParticle(fNUISANCEEvent, (*p_iter), kInitialState, true);
}
- // FSI State
- // for (size_t i = 0; i < npart_in; i++ ) {
- // AddNuWroParticle(fNUISANCEEvent, (*p_iter), kFSIState);
- // }
+ // Try to find the FSI state particles
+ // Loop over the primary vertex particles
+ // If they match the post-FSI they haven't undergone FSI.
+ // If they don't match post-FSI they have undergone FSI.
+ for (p_iter = fNuWroEvent->out.begin(); p_iter != fNuWroEvent->out.end(); p_iter++) {
+ // Get the particle
+ particle p = (*p_iter);
+ // Check against all the post particles, match them
+ std::vector<particle>::iterator p2_iter;
+ bool match = false;
+ for (p2_iter = fNuWroEvent->post.begin(); p2_iter != fNuWroEvent->post.end(); p2_iter++) {
+ particle p2 = (*p2_iter);
+ // Check energy and pdg
+ // A very small cascade which changes the energy by 1E-5 MeV should be matched
+ match = (fabs(p2.E()-p.E()) < 1E-5 && p2.pdg == p.pdg);
+ // If we match p to p2 break the loop
+ if (match) break;
+ }
+ // If we've looped through the whole particle stack of post-FSI and haven't found a match it's a primary particle that has been FSIed
+ if (!match) AddNuWroParticle(fNUISANCEEvent, (*p_iter), kFSIState, true);
+ }
- // Final State
- for (p_iter = fNuWroEvent->post.begin(); p_iter != fNuWroEvent->post.end();
- p_iter++) {
- AddNuWroParticle(fNUISANCEEvent, (*p_iter), kFinalState);
+ // Loop over the final state particles
+ for (p_iter = fNuWroEvent->post.begin(); p_iter != fNuWroEvent->post.end(); p_iter++) {
+ particle p = (*p_iter);
+ // To find if it's primary or not we have to loop through the primary ones and match, just like above
+ bool match = false;
+ std::vector<particle>::iterator p2_iter;
+ for (p2_iter = fNuWroEvent->out.begin(); p2_iter != fNuWroEvent->out.end(); p2_iter++) {
+ particle p2 = (*p2_iter);
+ match = (fabs(p2.E()-p.E()) < 1E-5 && p2.pdg == p.pdg);
+ if (match) break;
+ }
+ AddNuWroParticle(fNUISANCEEvent, (*p_iter), kFinalState, match);
}
// Fill Generator Info
if (fSaveExtra) fNuWroInfo->FillGeneratorInfo(fNuWroEvent);
// Run Initial, FSI, Final, Other ordering.
fNUISANCEEvent->OrderStack();
FitParticle* ISNeutralLepton =
fNUISANCEEvent->GetHMISParticle(PhysConst::pdg_neutrinos);
if (ISNeutralLepton) {
fNUISANCEEvent->probe_E = ISNeutralLepton->E();
fNUISANCEEvent->probe_pdg = ISNeutralLepton->PDG();
}
return;
}
void NuWroInputHandler::AddNuWroParticle(FitEvent* evt, particle& p,
- int state) {
+ int state, bool primary = false) {
// Add Mom
evt->fParticleMom[evt->fNParticles][0] = static_cast<vect&>(p).x;
evt->fParticleMom[evt->fNParticles][1] = static_cast<vect&>(p).y;
evt->fParticleMom[evt->fNParticles][2] = static_cast<vect&>(p).z;
evt->fParticleMom[evt->fNParticles][3] = static_cast<vect&>(p).t;
+ // For NuWro a particle that we've given a FSI state is a pre-FSI particle
+ // An initial state particle is also a primary vertex praticle
+ evt->fPrimaryVertex[evt->fNParticles] = primary;
+
// Status/PDG
evt->fParticleState[evt->fNParticles] = state;
evt->fParticlePDG[evt->fNParticles] = p.pdg;
// Add to particle count
evt->fNParticles++;
}
void NuWroInputHandler::Print() {}
#endif
diff --git a/src/InputHandler/NuWroInputHandler.h b/src/InputHandler/NuWroInputHandler.h
index 5620a0f..43d7443 100644
--- a/src/InputHandler/NuWroInputHandler.h
+++ b/src/InputHandler/NuWroInputHandler.h
@@ -1,95 +1,95 @@
#ifndef NuWroINPUTHANDLER_H
#define NuWroINPUTHANDLER_H
#ifdef __NUWRO_ENABLED__
#include "GeneratorUtils.h"
#include "InputHandler.h"
#include "PlotUtils.h"
/// NuWro Generator Container to save extra particle status codes.
class NuWroGeneratorInfo : public GeneratorInfoBase {
public:
NuWroGeneratorInfo(){};
virtual ~NuWroGeneratorInfo();
/// Assigns information to branches
void AddBranchesToTree(TTree* tn);
/// Setup reading information from branches
void SetBranchesFromTree(TTree* tn);
/// Allocate any dynamic arrays for a new particle stack size
void AllocateParticleStack(int stacksize);
/// Clear any dynamic arrays
void DeallocateParticleStack();
/// Read extra genie information from the event
void FillGeneratorInfo(event* e);
/// Reset extra information to default/empty values
void Reset();
int kMaxParticles; ///< Number of particles in stack
int* fNuWroParticlePDGs; ///< NuWro Particle PDGs (example)
};
/// Main NuWro Input Reader. Requires events have flux and xsec TH1Ds saved into
/// them.
class NuWroInputHandler : public InputHandlerBase {
#ifdef __USE_NUWRO_SRW_EVENTS__
params fNuwroParams;
std::vector<BaseFitEvt> rwEvs;
#endif
public:
/// Constructor. Can handle single and joint inputs.
NuWroInputHandler(std::string const& handle, std::string const& rawinputs);
~NuWroInputHandler();
/// Create a TTree Cache to speed up file read
void CreateCache();
/// Remove TTree Cache to save memory
void RemoveCache();
/// Returns filled NUISANCEEvent for given entry.
FitEvent* GetNuisanceEvent(const UInt_t entry,
const bool lightweight = false);
#ifdef __USE_NUWRO_SRW_EVENTS__
// Returns filled BaseFitEvent for a given entry;
BaseFitEvt* GetBaseEvent(const UInt_t entry) {
if (rwEvs.size() <= entry) {
THROW("Tried to get cached BaseFitEv[" << entry << "], but only have "
<< rwEvs.size()
<< " in the cache.");
}
return &rwEvs[entry];
}
#endif
/// Fills fNUISANCEEvent from fNuWroEvent
void CalcNUISANCEKinematics();
/// (LEGACY) Automatically creates nuwro flux/event histograms that
/// nuisance needs to normalise events.
void ProcessNuWroInputFlux(const std::string file);
/// Calculates a True Interaction code for NuWro events
int ConvertNuwroMode(event* e);
/// Adds a new particle to NUISANCE stack for given NuWro particle
- void AddNuWroParticle(FitEvent* evt, particle& p, int state);
+ void AddNuWroParticle(FitEvent* evt, particle& p, int state, bool primary);
event* fNuWroEvent; ///< Pointer to NuWro Format Events
/// Print Event Information
void Print();
TChain* fNuWroTree; ///< TTree for reading NuWro event vectors
bool fSaveExtra; ///< Save Extra NuWro info into Nuisance Event
NuWroGeneratorInfo* fNuWroInfo; ///< Extra NuWro Generator Info
};
/*! @} */
#endif
#endif
diff --git a/src/MCStudies/GenericFlux_Vectors.cxx b/src/MCStudies/GenericFlux_Vectors.cxx
index 810d261..bb2e1e8 100644
--- a/src/MCStudies/GenericFlux_Vectors.cxx
+++ b/src/MCStudies/GenericFlux_Vectors.cxx
@@ -1,383 +1,401 @@
// 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_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 = 1E10; // Arbritrarily high energy limit
if (Config::HasPar("EnuMin")) {
EnuMin = Config::GetParD("EnuMin");
}
if (Config::HasPar("EnuMax")) {
EnuMax = Config::GetParD("EnuMax");
}
+ SavePreFSI = Config::Get().GetParB("nuisflat_SavePreFSI");
+ LOG(SAM) << "Running GenericFlux_Vectors saving pre-FSI particles? " << SavePreFSI << std::endl;
+
// 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, EnuMax) / double(fNEvents)) /
this->TotalIntegratedFlux("width");
LOG(SAM) << " Generic Flux Scaling Factor = " << fScaleFactor
<< " [= " << (GetEventHistogram()->Integral("width") * 1E-38) << "/("
<< (fNEvents + 0.) << "*" << TotalIntegratedFlux("width") << ")]"
<< 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("tgta", &tgta, "tgta/I");
eventVariables->Branch("tgtz", &tgtz, "tgtz/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("W_genie", &W_genie, "W_genie/F");
eventVariables->Branch("x", &x, "x/F");
eventVariables->Branch("y", &y, "y/F");
eventVariables->Branch("Eav", &Eav, "Eav/F");
eventVariables->Branch("EavAlt", &EavAlt, "EavAlt/F");
eventVariables->Branch("dalphat", &dalphat, "dalphat/F");
eventVariables->Branch("dpt", &dpt, "dpt/F");
eventVariables->Branch("dphit", &dphit, "dphit/F");
eventVariables->Branch("pnreco_C", &pnreco_C, "pnreco_C/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");
eventVariables->Branch("pdg_rank", pdg_rank, "pdg_rank[nfsp]/I");
+ eventVariables->Branch("status", status, "status[nfsp]/I");
+ eventVariables->Branch("isalive", isalive, "isalive[nfsp]/O");
+ eventVariables->Branch("isprimary", isprimary, "isprimary[nfsp]/O");
// Event Scaling Information
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");
// The customs
eventVariables->Branch("CustomWeight", &CustomWeight, "CustomWeight/F");
eventVariables->Branch("CustomWeightArray", CustomWeightArray, "CustomWeightArray[6]/F");
return;
}
void GenericFlux_Vectors::FillEventVariables(FitEvent *event) {
ResetVariables();
// Fill Signal Variables
FillSignalFlags(event);
LOG(DEB) << "Filling signal" << std::endl;
// 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;
tgta = event->fTargetA;
tgtz = event->fTargetZ;
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);
Eav = FitUtils::GetErecoil_MINERvA_LowRecoil(event)/1.E3;
EavAlt = FitUtils::Eavailable(event)/1.E3;
// Get W_true with assumption of initial state nucleon at rest
float m_n = (float)PhysConst::mass_proton;
// Q2 assuming nucleon at rest
W_nuc_rest = sqrt(-Q2 + 2 * m_n * q0 + m_n * m_n);
// True Q2
W = sqrt(-Q2 + 2 * m_n * q0 + m_n * m_n);
x = Q2 / (2 * m_n * q0);
y = 1 - ELep / Enu_true;
dalphat = FitUtils::Get_STV_dalphat(event, PDGnu, true);
dpt = FitUtils::Get_STV_dpt(event, PDGnu, true);
dphit = FitUtils::Get_STV_dphit(event, PDGnu, true);
pnreco_C = FitUtils::Get_pn_reco_C(event, PDGnu, 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 &&
event->PartInfo(i)->Status() == kFinalState;
- if (!part_alive) continue;
+ if (!SavePreFSI) {
+ if (!part_alive) continue;
+ }
partList.push_back(event->PartInfo(i));
}
// Save outgoing particle vectors
nfsp = (int)partList.size();
std::map<int, std::vector<std::pair<double, int> > > pdgMap;
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;
+ status[i] = partList[i]->Status();
+ isalive[i] = partList[i]->fIsAlive;
+ isprimary[i] = event->fPrimaryVertex[i];
pdgMap[pdg[i]].push_back(std::make_pair(partList[i]->fP.Vect().Mag(), i));
}
for(std::map<int, std::vector<std::pair<double, int> > >::iterator iter = pdgMap.begin();
iter != pdgMap.end(); ++iter){
std::vector<std::pair<double, int> > thisVect = iter->second;
std::sort(thisVect.begin(), thisVect.end());
// Now save the order... a bit funky to avoid inverting
int nPart = (int)thisVect.size()-1;
for (int i=nPart; i >= 0; --i) {
pdg_rank[thisVect[i].second] = nPart-i;
}
}
#ifdef __GENIE_ENABLED__
if (event->fType == kGENIE) {
EventRecord * gevent = static_cast<EventRecord*>(event->genie_event->event);
const Interaction * interaction = gevent->Summary();
const Kinematics & kine = interaction->Kine();
double W_genie = kine.W();
}
#endif
// Fill event weights
Weight = event->RWWeight * event->InputWeight;
RWWeight = event->RWWeight;
InputWeight = event->InputWeight;
// And the Customs
CustomWeight = event->CustomWeight;
for (int i = 0; i < 6; ++i) {
CustomWeightArray[i] = event->CustomWeightArray[i];
}
// Fill the eventVariables Tree
eventVariables->Fill();
return;
};
//********************************************************************
void GenericFlux_Vectors::ResetVariables() {
//********************************************************************
cc = false;
// Reset all Function used to extract any variables of interest to the event
Mode = PDGnu = tgt = tgta = tgtz = PDGLep = 0;
Enu_true = ELep = CosLep = Q2 = q0 = q3 = Enu_QE = Q2_QE = W_nuc_rest = W = x = y = Eav = EavAlt = -999.9;
W_genie = -999;
// Other fun variables
// MINERvA-like ones
dalphat = dpt = dphit = pnreco_C = -999.99;
nfsp = 0;
for (int i = 0; i < kMAX; ++i){
px[i] = py[i] = pz[i] = E[i] = -999;
+<<<<<<< HEAD
pdg[i] = pdg_rank[i] = 0;
+=======
+ pdg[i] = 0;
+ status[i] = -999;
+ isalive[i] = false;
+ isprimary[i] = false;
+>>>>>>> 8cb8d610b53c4930caea73bde1186b7dbb04981c
}
Weight = InputWeight = RWWeight = 0.0;
CustomWeight = 0.0;
for (int i = 0; i < 6; ++i) CustomWeightArray[i] = 0.0;
partList.clear();
flagCCINC = flagNCINC = flagCCQE = flagCC0pi = flagCCQELike = flagNCEL = flagNC0pi = flagCCcoh = flagNCcoh = flagCC1pip = flagNC1pip = flagCC1pim = flagNC1pim = flagCC1pi0 = flagNC1pi0 = false;
}
//********************************************************************
void GenericFlux_Vectors::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);
}
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 9c9a9b7..5a091ed 100644
--- a/src/MCStudies/GenericFlux_Vectors.h
+++ b/src/MCStudies/GenericFlux_Vectors.h
@@ -1,136 +1,141 @@
// 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 GenericFlux_Vectors_H_SEEN
#define GenericFlux_Vectors_H_SEEN
#include "Measurement1D.h"
#include "FitEvent.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 signal flags
void FillSignalFlags(FitEvent *event);
void ResetVariables();
//! 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<FitParticle*> partList;
+ bool SavePreFSI;
+
int Mode;
bool cc;
int PDGnu;
int tgt;
int tgta;
int tgtz;
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;
float Eav;
float EavAlt;
float dalphat;
float W_genie;
float dpt;
float dphit;
float pnreco_C;
// 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];
int pdg_rank[kMAX];
+ int status[kMAX];
+ bool isalive[kMAX];
+ bool isprimary[kMAX];
// Basic event info
float Weight;
float InputWeight;
float RWWeight;
double fScaleFactor;
// Custom weights
float CustomWeight;
float CustomWeightArray[6];
// 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