diff --git a/src/InputHandler/FitEvent.cxx b/src/InputHandler/FitEvent.cxx index d075570..f3d5876 100644 --- a/src/InputHandler/FitEvent.cxx +++ b/src/InputHandler/FitEvent.cxx @@ -1,437 +1,441 @@ // 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::GetBeamPartIndex(void) const { + return GetHMISParticleIndex(this->probe_pdg); +} + 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 a1970f3..55a15bb 100644 --- a/src/InputHandler/FitEvent.h +++ b/src/InputHandler/FitEvent.h @@ -1,639 +1,647 @@ // 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); }; + inline bool IsCC() const { if (abs(this->probe_pdg) == 11) return false; return (abs(Mode) <= 30); }; /// Is Event Neutral Current? - inline bool IsNC() const { return (abs(Mode) > 30); }; + inline bool IsNC() const { if (abs(this->probe_pdg) == 11) return true; 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()); }; + // ---- Generic beam incoming functions + // I'm not 100% sure why these can't replace the above (FitEvent knows the type) + int GetBeamPartIndex (void) const; + inline TLorentzVector GetBeamPartP4 (void) const { return GetParticleP4(GetBeamPartIndex()); }; + inline TVector3 GetBeamPartP3 (void) const { return GetParticleP3(GetBeamPartIndex()); }; + inline double GetBeamPartMom (void) const { return GetParticleMom(GetBeamPartIndex()); }; + inline double GetBeamPartMom2 (void) const { return GetParticleMom2(GetBeamPartIndex()); }; + inline double GetBeamPartE (void) const { return GetParticleE(GetBeamPartIndex()); }; + inline int GetBeamPartPDG (void) const { return GetParticlePDG(GetBeamPartIndex()); }; + inline int GetPartInPos (void) const { return GetBeamPartIndex(); }; + inline FitParticle* GetBeamPart (void) { return GetParticle(GetBeamPartIndex()); }; + /// 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