diff --git a/src/MCStudies/GenericFlux_Vectors.cxx b/src/MCStudies/GenericFlux_Vectors.cxx
index ad7d530..308144d 100644
--- a/src/MCStudies/GenericFlux_Vectors.cxx
+++ b/src/MCStudies/GenericFlux_Vectors.cxx
@@ -1,397 +1,436 @@
// 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");
+
+ // Save init particle vectors
+ eventVariables->Branch("ninitp", &ninitp, "ninitp/I");
+ eventVariables->Branch("px_init", px_init, "px_init[ninitp]/F");
+ eventVariables->Branch("py_init", py_init, "py_init[ninitp]/F");
+ eventVariables->Branch("pz_init", pz_init, "pz_init[ninitp]/F");
+ eventVariables->Branch("E_init", E_init, "E_init[ninitp]/F");
+ eventVariables->Branch("pdg_init", pdg_init, "pdg_init[ninitp]/I");
+
+ // Save pre-FSI vectors
+ eventVariables->Branch("nvertp", &nvertp, "nvertp/I");
+ eventVariables->Branch("px_vert", px_vert, "px_vert[nvertp]/F");
+ eventVariables->Branch("py_vert", py_vert, "py_vert[nvertp]/F");
+ eventVariables->Branch("pz_vert", pz_vert, "pz_vert[nvertp]/F");
+ eventVariables->Branch("E_vert", E_vert, "E_vert[nvertp]/F");
+ eventVariables->Branch("pdg_vert", pdg_vert, "pdg_vert[nvertp]/I");
// 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 (!SavePreFSI) {
- if (!part_alive) continue;
- }
+ if (event->PartInfo(i)->fIsAlive &&
+ event->PartInfo(i)->Status() == kFinalState)
+ partList.push_back(event->PartInfo(i));
+
+ if (SavePreFSI && event->fPrimaryVertex[i])
+ vertList.push_back(event->PartInfo(i));
- partList.push_back(event->PartInfo(i));
+ if (SavePreFSI && event->PartInfo(i)->IsInitialState())
+ initList.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;
}
}
+ // Save pre-FSI particles
+ nvertp = (int)vertList.size();
+ for (int i = 0; i < nvertp; ++i) {
+ px_vert[i] = vertList[i]->fP.X() / 1E3;
+ py_vert[i] = vertList[i]->fP.Y() / 1E3;
+ pz_vert[i] = vertList[i]->fP.Z() / 1E3;
+ E_vert[i] = vertList[i]->fP.E() / 1E3;
+ pdg_vert[i] = vertList[i]->fPID;
+ }
+
+ // Save init particles
+ ninitp = (int)initList.size();
+ for (int i = 0; i < ninitp; ++i) {
+ px_init[i] = initList[i]->fP.X() / 1E3;
+ py_init[i] = initList[i]->fP.Y() / 1E3;
+ pz_init[i] = initList[i]->fP.Z() / 1E3;
+ E_init[i] = initList[i]->fP.E() / 1E3;
+ pdg_init[i] = initList[i]->fPID;
+ }
+
+
#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;
+ nfsp = ninitp = nvertp = 0;
for (int i = 0; i < kMAX; ++i){
px[i] = py[i] = pz[i] = E[i] = -999;
pdg[i] = pdg_rank[i] = 0;
- status[i] = -999;
- isalive[i] = false;
- isprimary[i] = false;
+
+ px_init[i] = py_init[i] = pz_init[i] = E_init[i] = -999;
+ pdg_init[i] = 0;
+
+ px_vert[i] = py_vert[i] = pz_vert[i] = E_vert[i] = -999;
+ pdg_vert[i] = 0;
+
}
Weight = InputWeight = RWWeight = 0.0;
CustomWeight = 0.0;
for (int i = 0; i < 6; ++i) CustomWeightArray[i] = 0.0;
partList.clear();
+ initList.clear();
+ vertList.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 5a091ed..d388e95 100644
--- a/src/MCStudies/GenericFlux_Vectors.h
+++ b/src/MCStudies/GenericFlux_Vectors.h
@@ -1,141 +1,156 @@
// 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;
+ std::vector<FitParticle*> initList;
+ std::vector<FitParticle*> vertList;
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];
+
+ // Save incoming particle info
+ int ninitp;
+ float px_init[kMAX];
+ float py_init[kMAX];
+ float pz_init[kMAX];
+ float E_init[kMAX];
+ int pdg_init[kMAX];
+
+ // Save pre-FSI particle info
+ int nvertp;
+ float px_vert[kMAX];
+ float py_vert[kMAX];
+ float pz_vert[kMAX];
+ float E_vert[kMAX];
+ int pdg_vert[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