diff --git a/parameters/config.xml b/parameters/config.xml
index fb75bda..5ba8803 100644
--- a/parameters/config.xml
+++ b/parameters/config.xml
@@ -1,229 +1,232 @@
+
+
+
Using z-expansion
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 .
*******************************************************************************/
#include "FitEvent.h"
#include
#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] <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 FitEvent::GetAllParticleIndices(int const pdg,
int const state) const {
std::vector 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 FitEvent::GetAllParticle(int const pdg,
int const state) {
std::vector indexlist = GetAllParticleIndices(pdg, state);
std::vector plist;
for (std::vector::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 .
*******************************************************************************/
#ifndef FITEVENT2_H_SEEN
#define FITEVENT2_H_SEEN
/*!
* \addtogroup InputHandler
* @{
*/
#include
#include
#include
#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
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
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 GetAllParticleIndices (int const pdg = -1, int const state = -1) const;
template
inline std::vector GetAllParticleIndices(int const (&pdgs)[N], const int state = -1) const {
std::vector plist;
for (size_t i = 0; i < N; i++) {
std::vector 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 GetAllParticle (int const pdg = -1, int const state = -1);
template
inline std::vector GetAllParticle(int const (&pdgs)[N], int const state = -1) {
std::vector plist;
for (size_t i = 0; i < N; i++) {
std::vector plisttemp = GetAllParticle(pdgs[i], state);
plist.insert( plist.end(), plisttemp.begin(), plisttemp.end() );
}
return plist;
}
inline std::vector GetAllNuElectronIndices (void) { return GetAllParticleIndices(12); };
inline std::vector GetAllNuMuonIndices (void) { return GetAllParticleIndices(14); };
inline std::vector GetAllNuTauIndices (void) { return GetAllParticleIndices(16); };
inline std::vector GetAllElectronIndices (void) { return GetAllParticleIndices(11); };
inline std::vector GetAllMuonIndices (void) { return GetAllParticleIndices(13); };
inline std::vector GetAllTauIndices (void) { return GetAllParticleIndices(15); };
inline std::vector GetAllProtonIndices (void) { return GetAllParticleIndices(2212); };
inline std::vector GetAllNeutronIndices (void) { return GetAllParticleIndices(2112); };
inline std::vector GetAllPiZeroIndices (void) { return GetAllParticleIndices(111); };
inline std::vector GetAllPiPlusIndices (void) { return GetAllParticleIndices(211); };
inline std::vector GetAllPiMinusIndices (void) { return GetAllParticleIndices(-211); };
inline std::vector GetAllPhotonIndices (void) { return GetAllParticleIndices(22); };
inline std::vector GetAllNuElectron (void) { return GetAllParticle(12); };
inline std::vector GetAllNuMuon (void) { return GetAllParticle(14); };
inline std::vector GetAllNuTau (void) { return GetAllParticle(16); };
inline std::vector GetAllElectron (void) { return GetAllParticle(11); };
inline std::vector GetAllMuon (void) { return GetAllParticle(13); };
inline std::vector GetAllTau (void) { return GetAllParticle(15); };
inline std::vector GetAllProton (void) { return GetAllParticle(2212); };
inline std::vector GetAllNeutron (void) { return GetAllParticle(2112); };
inline std::vector GetAllPiZero (void) { return GetAllParticle(111); };
inline std::vector GetAllPiPlus (void) { return GetAllParticle(211); };
inline std::vector GetAllPiMinus (void) { return GetAllParticle(-211); };
inline std::vector 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
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
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
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
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 GetAllISParticleIndices(int const pdg = -1) const {
return GetAllParticleIndices(pdg, kInitialState);
};
template
inline std::vector 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 GetAllISParticle(int const pdg = -1) {
return GetAllParticle(pdg, kInitialState);
};
template
inline std::vector 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
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
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 GetAllISNuElectronIndices (void) const { return GetAllISParticleIndices(12); };
inline std::vector GetAllISNuMuonIndices (void) const { return GetAllISParticleIndices(14); };
inline std::vector GetAllISNuTauIndices (void) const { return GetAllISParticleIndices(16); };
inline std::vector GetAllISElectronIndices (void) const { return GetAllISParticleIndices(11); };
inline std::vector GetAllISMuonIndices (void) const { return GetAllISParticleIndices(13); };
inline std::vector GetAllISTauIndices (void) const { return GetAllISParticleIndices(15); };
inline std::vector GetAllISProtonIndices (void) const { return GetAllISParticleIndices(2212); };
inline std::vector GetAllISNeutronIndices (void) const { return GetAllISParticleIndices(2112); };
inline std::vector GetAllISPiZeroIndices (void) const { return GetAllISParticleIndices(111); };
inline std::vector GetAllISPiPlusIndices (void) const { return GetAllISParticleIndices(211); };
inline std::vector GetAllISPiMinusIndices (void) const { return GetAllISParticleIndices(-211); };
inline std::vector GetAllISPhotonIndices (void) const { return GetAllISParticleIndices(22); };
inline std::vector GetAllISLeptonsIndices (void) const { return GetAllISParticleIndices(PhysConst::pdg_leptons); };
inline std::vector GetAllISPionsIndices (void) const { return GetAllISParticleIndices(PhysConst::pdg_pions); };
inline std::vector GetAllISChargePionsIndices(void) const { return GetAllISParticleIndices(PhysConst::pdg_charged_pions); };
inline std::vector GetAllISNuElectron (void) { return GetAllISParticle(12); };
inline std::vector GetAllISNuMuon (void) { return GetAllISParticle(14); };
inline std::vector GetAllISNuTau (void) { return GetAllISParticle(16); };
inline std::vector GetAllISElectron (void) { return GetAllISParticle(11); };
inline std::vector GetAllISMuon (void) { return GetAllISParticle(13); };
inline std::vector GetAllISTau (void) { return GetAllISParticle(15); };
inline std::vector GetAllISProton (void) { return GetAllISParticle(2212); };
inline std::vector GetAllISNeutron (void) { return GetAllISParticle(2112); };
inline std::vector GetAllISPiZero (void) { return GetAllISParticle(111); };
inline std::vector GetAllISPiPlus (void) { return GetAllISParticle(211); };
inline std::vector GetAllISPiMinus (void) { return GetAllISParticle(-211); };
inline std::vector GetAllISPhoton (void) { return GetAllISParticle(22); };
inline std::vector GetAllISLeptons (void) { return GetAllISParticle(PhysConst::pdg_leptons); };
inline std::vector GetAllISPions (void) { return GetAllISParticle(PhysConst::pdg_pions); };
inline std::vector 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
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
inline int NumFSParticle(int const (&pdgs)[N]) const {
return NumParticle(pdgs, kFinalState);
};
inline std::vector GetAllFSParticleIndices(int const pdg = -1) const {
return GetAllParticleIndices(pdg, kFinalState);
};
template
inline std::vector GetAllFSParticleIndices(int const (&pdgs)[N]) const {
return GetAllParticleIndices(pdgs, kFinalState);
};
inline std::vector GetAllFSParticle(int const pdg = -1) {
return GetAllParticle(pdg, kFinalState);
};
template
inline std::vector GetAllFSParticle(int const (&pdgs)[N]) {
return GetAllParticle(pdgs, kFinalState);
};
inline FitParticle* GetHMFSParticle(int const pdg) {
return GetHMParticle(pdg, kFinalState);
};
template
inline FitParticle* GetHMFSParticle(int const (&pdgs)[N]) {
return GetHMParticle(pdgs, kFinalState);
};
inline int GetHMFSParticleIndex(int const pdg) const {
return GetHMParticleIndex(pdg, kFinalState);
};
template
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 GetAllFSNuElectronIndices (void) const { return GetAllFSParticleIndices(12); };
inline std::vector GetAllFSNuMuonIndices (void) const { return GetAllFSParticleIndices(14); };
inline std::vector GetAllFSNuTauIndices (void) const { return GetAllFSParticleIndices(16); };
inline std::vector GetAllFSElectronIndices (void) const { return GetAllFSParticleIndices(11); };
inline std::vector GetAllFSMuonIndices (void) const { return GetAllFSParticleIndices(13); };
inline std::vector GetAllFSTauIndices (void) const { return GetAllFSParticleIndices(15); };
inline std::vector GetAllFSProtonIndices (void) const { return GetAllFSParticleIndices(2212); };
inline std::vector GetAllFSNeutronIndices (void) const { return GetAllFSParticleIndices(2112); };
inline std::vector GetAllFSPiZeroIndices (void) const { return GetAllFSParticleIndices(111); };
inline std::vector GetAllFSPiPlusIndices (void) const { return GetAllFSParticleIndices(211); };
inline std::vector GetAllFSPiMinusIndices (void) const { return GetAllFSParticleIndices(-211); };
inline std::vector GetAllFSPhotonIndices (void) const { return GetAllFSParticleIndices(22); };
inline std::vector GetAllFSLeptonsIndices (void) const { return GetAllFSParticleIndices(PhysConst::pdg_leptons); };
inline std::vector GetAllFSPionsIndices (void) const { return GetAllFSParticleIndices(PhysConst::pdg_pions); };
inline std::vector GetAllFSChargePionsIndices(void) const { return GetAllFSParticleIndices(PhysConst::pdg_charged_pions); };
inline std::vector GetAllFSNuElectron (void) { return GetAllFSParticle(12); };
inline std::vector GetAllFSNuMuon (void) { return GetAllFSParticle(14); };
inline std::vector GetAllFSNuTau (void) { return GetAllFSParticle(16); };
inline std::vector GetAllFSElectron (void) { return GetAllFSParticle(11); };
inline std::vector GetAllFSMuon (void) { return GetAllFSParticle(13); };
inline std::vector GetAllFSTau (void) { return GetAllFSParticle(15); };
inline std::vector GetAllFSProton (void) { return GetAllFSParticle(2212); };
inline std::vector GetAllFSNeutron (void) { return GetAllFSParticle(2112); };
inline std::vector GetAllFSPiZero (void) { return GetAllFSParticle(111); };
inline std::vector GetAllFSPiPlus (void) { return GetAllFSParticle(211); };
inline std::vector GetAllFSPiMinus (void) { return GetAllFSParticle(-211); };
inline std::vector GetAllFSPhoton (void) { return GetAllFSParticle(22); };
inline std::vector GetAllFSLeptons (void) { return GetAllFSParticle(PhysConst::pdg_leptons); };
inline std::vector GetAllFSPions (void) { return GetAllFSParticle(PhysConst::pdg_pions); };
inline std::vector 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 .
*******************************************************************************/
#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:"<.
*******************************************************************************/
#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(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((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 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(fGenieNtpl->event);
if (!fGenieGHep) return NULL;
TObjArrayIter iter(fGenieGHep);
genie::GHepParticle* p;
while ((p = (dynamic_cast((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(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((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);
+ //std::cout << curpart << " " << p->Pdg() << " " << p->FirstMother() << std::endl;
+ //std::cout << p->LastMother() << std::endl;
+ //std::cout << "**" << std::endl;
+
// 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/NuWroInputHandler.cxx b/src/InputHandler/NuWroInputHandler.cxx
index 9c7bf05..83dc5f2 100644
--- a/src/InputHandler/NuWroInputHandler.cxx
+++ b/src/InputHandler/NuWroInputHandler.cxx
@@ -1,480 +1,497 @@
#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 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::iterator p_iter;
// Initial State
for (p_iter = fNuWroEvent->in.begin(); p_iter != fNuWroEvent->in.end();
p_iter++) {
AddNuWroParticle(fNUISANCEEvent, (*p_iter), kInitialState);
}
- // FSI State
- // for (size_t i = 0; i < npart_in; i++ ) {
- // AddNuWroParticle(fNUISANCEEvent, (*p_iter), kFSIState);
- // }
+ for (p_iter = fNuWroEvent->all.begin(); p_iter != fNuWroEvent->all.end(); p_iter++) {
+ std::cout << (*p_iter)->pdg << std::endl;
+ std::cout << (*p_iter)->ks << std::endl;
+ std::cout << (*p_iter)->orgig << std::endl;
+ std::cout << (*p_iter)->id << std::endl;
+ std::cout << (*p_iter)->mother << std::endl;
+ std::cout << (*p_iter)->endproc << std::endl;
+ std::cout << (*p_iter)->his_fermi << std::endl;
+ std::cout << (*p_iter)->primary << std::endl;
+ std::cout << (*p_iter)->E() << std::endl;
+ }
+
+ // NuWro saves the incoming, pre-FSI and post-FSI particles
+ // So we fill the incoming and outgoing particles, but no the FSI particles
+ // This is because sometimes the pre-FSI particles are the same as the post-FSI particles, so we would double count our particle stack
+ /*
+ for (p_iter = fNuWroEvent->out.begin(); p_iter != fNuWroEvent->out.end();
+ p_iter++) {
+ AddNuWroParticle(fNUISANCEEvent, (*p_iter), kFSIState);
+ }
+ */
// Final State
for (p_iter = fNuWroEvent->post.begin(); p_iter != fNuWroEvent->post.end();
p_iter++) {
AddNuWroParticle(fNUISANCEEvent, (*p_iter), kFinalState);
}
// 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) {
// Add Mom
evt->fParticleMom[evt->fNParticles][0] = static_cast(p).x;
evt->fParticleMom[evt->fNParticles][1] = static_cast(p).y;
evt->fParticleMom[evt->fNParticles][2] = static_cast(p).z;
evt->fParticleMom[evt->fNParticles][3] = static_cast(p).t;
// 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/MCStudies/GenericFlux_Vectors.cxx b/src/MCStudies/GenericFlux_Vectors.cxx
index 6ab61b7..a8cec51 100644
--- a/src/MCStudies/GenericFlux_Vectors.cxx
+++ b/src/MCStudies/GenericFlux_Vectors.cxx
@@ -1,368 +1,383 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#include "GenericFlux_Vectors.h"
GenericFlux_Vectors::GenericFlux_Vectors(std::string name,
std::string inputfile, FitWeight *rw,
std::string type,
std::string fakeDataFile) {
// Measurement Details
fName = name;
eventVariables = NULL;
// Define our energy range for flux calcs
EnuMin = 0.;
EnuMax = 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("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;
+ }
+ //std::cout << event->PartInfo(i)->Status() << std::endl;
partList.push_back(event->PartInfo(i));
}
// Save outgoing particle vectors
nfsp = (int)partList.size();
for (int i = 0; i < nfsp; ++i) {
px[i] = partList[i]->fP.X() / 1E3;
py[i] = partList[i]->fP.Y() / 1E3;
pz[i] = partList[i]->fP.Z() / 1E3;
E[i] = partList[i]->fP.E() / 1E3;
pdg[i] = partList[i]->fPID;
+ status[i] = partList[i]->Status();
+ isalive[i] = partList[i]->fIsAlive;
+ isprimary[i] = event->fPrimaryVertex[i];
}
#ifdef __GENIE_ENABLED__
if (event->fType == kGENIE) {
EventRecord * gevent = static_cast(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;
pdg[i] = 0;
+ status[i] = -999;
+ isalive[i] = false;
+ isprimary[i] = false;
}
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 181a283..770f32d 100644
--- a/src/MCStudies/GenericFlux_Vectors.h
+++ b/src/MCStudies/GenericFlux_Vectors.h
@@ -1,135 +1,140 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#ifndef GenericFlux_Vectors_H_SEEN
#define GenericFlux_Vectors_H_SEEN
#include "Measurement1D.h"
#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 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 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