diff --git a/src/InputHandler/FitEvent.cxx b/src/InputHandler/FitEvent.cxx index 0150633..eabd34b 100644 --- a/src/InputHandler/FitEvent.cxx +++ b/src/InputHandler/FitEvent.cxx @@ -1,465 +1,493 @@ // 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 "TObjArray.h" #include 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) { NUIS_LOG(DEB, "Allocating particle stack of size: " << stacksize); 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)) { NUIS_LOG(DEB, "Ordered stack"); for (int i = 0; i < fNParticles; i++) { NUIS_LOG(DEB, "Particle " << i << ". " << fParticlePDG[i] << " " << fParticleMom[i][0] << " " << fParticleMom[i][1] << " " << fParticleMom[i][2] << " " << fParticleMom[i][3] << " " << fParticleState[i]); } } if (fNParticles != npart) { NUIS_ABORT("Dropped some particles when ordering the stack!"); } return; } void FitEvent::Print() { if (LOG_LEVEL(FIT)) { NUIS_LOG(FIT, "FITEvent print"); NUIS_LOG(FIT, "Mode: " << Mode << ", Weight: " << InputWeight); NUIS_LOG(FIT, "Particles: " << fNParticles); NUIS_LOG(FIT, " -> Particle Stack "); for (int i = 0; i < fNParticles; i++) { NUIS_LOG(FIT, " -> -> " << i << ". " << fParticlePDG[i] << " " << fParticleState[i] << " " << " Mom(" << fParticleMom[i][0] << ", " << fParticleMom[i][1] << ", " << fParticleMom[i][2] << ", " << fParticleMom[i][3] << ")."); } } 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) { NUIS_ABORT("Requesting particle beyond stack!" << std::endl << "i = " << i << " N = " << fNParticles << std::endl << "Mode = " << Mode); } 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; } + +// Get the outgoing lepton PDG depending on if it's a CC or NC event +int FitEvent::GetLeptonOutPDG() { + // Make sure the outgoing lepton has the correct PDG + int pdgnu = PDGnu(); + int PDGout; + if (IsCC()) { + if (pdgnu > 0) PDGout = pdgnu-1; + else PDGout = pdgnu+1; + } else { + PDGout = pdgnu; + } + return PDGout; +} + +//******************************************************************** +// Returns the true Q2 of an event +double FitEvent::GetQ2() { + FitParticle *neutrino = GetNeutrinoIn(); + FitParticle *lepton = GetLeptonOut(); + // Figure out why this happens + if (!neutrino || !lepton) return 0.0; + double Q2 = + -1.0 * (lepton->P4() - neutrino->P4()) * + (lepton->P4() - neutrino->P4()) / 1.E6; + return Q2; +} +//******************************************************************** diff --git a/src/InputHandler/FitEvent.h b/src/InputHandler/FitEvent.h index b467835..73b876a 100644 --- a/src/InputHandler/FitEvent.h +++ b/src/InputHandler/FitEvent.h @@ -1,653 +1,660 @@ // 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 { if (abs(this->probe_pdg) == 11) return false; return (abs(Mode) <= 30); }; /// Is Event Neutral Current? inline bool IsNC() const { if (abs(this->probe_pdg) == 11) return true; return (abs(Mode) > 30); }; + // Is Event resonant? + inline bool IsResonant() const { if (Mode != 11 && Mode != 12 && Mode != 13 && Mode != 31 && Mode != 32 && Mode != 33 && Mode != 34) return false; return true; }; /// 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 bool HasFSNucleons (void) const { return HasFSParticle(PhysConst::pdg_nucleons); }; 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 int NumFSNucleons (void) const { return NumFSParticle(PhysConst::pdg_nucleons); }; 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 GetAllFSNucleonIndices(void) const { return GetAllFSParticleIndices(PhysConst::pdg_nucleons); }; 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 std::vector GetAllFSNucleons (void) { return GetAllFSParticle(PhysConst::pdg_nucleons); }; 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 FitParticle* GetHMFSNucleons(void) { return GetHMFSParticle(PhysConst::pdg_nucleons); }; 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); }; inline int GetHMFSChargeNucleonIndex(void) const { return GetHMFSParticleIndex(PhysConst::pdg_nucleons); }; // ---- 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(); - int GetLeptonOutPos(void) const; - FitParticle* GetLeptonOut(void); + // Get outgoing lepton matching PDG of neutrino + int GetLeptonOutPDG(); + FitParticle* GetLeptonOut() {return GetHMFSParticle(GetLeptonOutPDG()); }; + // Get the outgoing lepton index + int GetLeptonIndex() { return GetHMFSParticleIndex(GetLeptonOutPDG()); }; + + double GetQ2(); // 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/Reweight/MINERvAWeightCalcs.cxx b/src/Reweight/MINERvAWeightCalcs.cxx index b9228d6..0f555b4 100644 --- a/src/Reweight/MINERvAWeightCalcs.cxx +++ b/src/Reweight/MINERvAWeightCalcs.cxx @@ -1,1109 +1,871 @@ #ifdef __MINERVA_RW_ENABLED__ #ifdef __GENIE_ENABLED__ #include "MINERvAWeightCalcs.h" #include "BaseFitEvt.h" namespace nuisance { namespace reweight { //******************************************************* MINERvAReWeight_QE::MINERvAReWeight_QE() { //******************************************************* fTwk_NormCCQE = 0.0; fDef_NormCCQE = 1.0; fCur_NormCCQE = fDef_NormCCQE; } //******************************************************* MINERvAReWeight_QE::~MINERvAReWeight_QE(){}; //******************************************************* //******************************************************* double MINERvAReWeight_QE::CalcWeight(BaseFitEvt *evt) { //******************************************************* // Check GENIE if (evt->fType != kGENIE) return 1.0; // Extract the GENIE Record GHepRecord *ghep = static_cast(evt->genie_event->event); const Interaction *interaction = ghep->Summary(); // const InitialState& init_state = interaction->InitState(); const ProcessInfo &proc_info = interaction->ProcInfo(); // const Target& tgt = init_state.Tgt(); // If the event is not QE this Calc doesn't handle it if (!proc_info.IsQuasiElastic()) return 1.0; // WEIGHT CALCULATIONS ------------- double w = 1.0; // CCQE Dial if (!proc_info.IsWeakCC()) w *= fCur_NormCCQE; // Return Combined Weight return w; } //******************************************************* void MINERvAReWeight_QE::SetDialValue(std::string name, double val) { //******************************************************* SetDialValue(Reweight::ConvDial(name, kCUSTOM), val); } //******************************************************* void MINERvAReWeight_QE::SetDialValue(int rwenum, double val) { //******************************************************* // Check Handled int curenum = rwenum % 1000; if (!IsHandled(curenum)) return; // Set Values if (curenum == Reweight::kMINERvARW_NormCCQE) { fTwk_NormCCQE = val; fCur_NormCCQE = fDef_NormCCQE + fTwk_NormCCQE; } // Define Tweaked fTweaked = ((fTwk_NormCCQE != 0.0)); } //******************************************************* bool MINERvAReWeight_QE::IsHandled(int rwenum) { //******************************************************* int curenum = rwenum % 1000; switch (curenum) { case Reweight::kMINERvARW_NormCCQE: return true; default: return false; } } //******************************************************* MINERvAReWeight_MEC::MINERvAReWeight_MEC() { //******************************************************* fTwk_NormCCMEC = 0.0; fDef_NormCCMEC = 1.0; fCur_NormCCMEC = fDef_NormCCMEC; } //******************************************************* MINERvAReWeight_MEC::~MINERvAReWeight_MEC(){}; //******************************************************* //******************************************************* double MINERvAReWeight_MEC::CalcWeight(BaseFitEvt *evt) { //******************************************************* // Check GENIE if (evt->fType != kGENIE) return 1.0; // Extract the GENIE Record GHepRecord *ghep = static_cast(evt->genie_event->event); const Interaction *interaction = ghep->Summary(); // const InitialState& init_state = interaction->InitState(); const ProcessInfo &proc_info = interaction->ProcInfo(); // const Target& tgt = init_state.Tgt(); // If the event is not MEC this Calc doesn't handle it if (!proc_info.IsMEC()) return 1.0; // WEIGHT CALCULATIONS ------------- double w = 1.0; // CCMEC Dial if (!proc_info.IsWeakCC()) w *= fCur_NormCCMEC; // Return Combined Weight return w; } //******************************************************* void MINERvAReWeight_MEC::SetDialValue(std::string name, double val) { //******************************************************* SetDialValue(Reweight::ConvDial(name, kCUSTOM), val); } //******************************************************* void MINERvAReWeight_MEC::SetDialValue(int rwenum, double val) { //******************************************************* // Check Handled int curenum = rwenum % 1000; if (!IsHandled(curenum)) return; // Set Values if (curenum == Reweight::kMINERvARW_NormCCMEC) { fTwk_NormCCMEC = val; fCur_NormCCMEC = fDef_NormCCMEC + fTwk_NormCCMEC; } // Define Tweaked fTweaked = ((fTwk_NormCCMEC != 0.0)); } //******************************************************* bool MINERvAReWeight_MEC::IsHandled(int rwenum) { //******************************************************* int curenum = rwenum % 1000; switch (curenum) { case Reweight::kMINERvARW_NormCCMEC: return true; default: return false; } } //******************************************************* MINERvAReWeight_RES::MINERvAReWeight_RES() { //******************************************************* fTwk_NormCCRES = 0.0; fDef_NormCCRES = 1.0; fCur_NormCCRES = fDef_NormCCRES; } //******************************************************* MINERvAReWeight_RES::~MINERvAReWeight_RES(){}; //******************************************************* //******************************************************* double MINERvAReWeight_RES::CalcWeight(BaseFitEvt *evt) { //******************************************************* // std::cout << "Caculating RES" << std::endl; // Check GENIE if (evt->fType != kGENIE) return 1.0; // Extract the GENIE Record GHepRecord *ghep = static_cast(evt->genie_event->event); const Interaction *interaction = ghep->Summary(); // const InitialState& init_state = interaction->InitState(); const ProcessInfo &proc_info = interaction->ProcInfo(); // const Target& tgt = init_state.Tgt(); // If the event is not RES this Calc doesn't handle it if (!proc_info.IsResonant()) return 1.0; // WEIGHT CALCULATIONS ------------- double w = 1.0; // CCRES Dial if (proc_info.IsWeakCC()) w *= fCur_NormCCRES; // Return Combined Weight return w; } //******************************************************* void MINERvAReWeight_RES::SetDialValue(std::string name, double val) { //******************************************************* SetDialValue(Reweight::ConvDial(name, kCUSTOM), val); } //******************************************************* void MINERvAReWeight_RES::SetDialValue(int rwenum, double val) { //******************************************************* // Check Handled int curenum = rwenum % 1000; if (!IsHandled(curenum)) return; // Set Values if (curenum == Reweight::kMINERvARW_NormCCRES) { fTwk_NormCCRES = val; fCur_NormCCRES = fDef_NormCCRES + fTwk_NormCCRES; } // Define Tweaked fTweaked = ((fTwk_NormCCRES != 0.0)); } //******************************************************* bool MINERvAReWeight_RES::IsHandled(int rwenum) { //******************************************************* int curenum = rwenum % 1000; switch (curenum) { case Reweight::kMINERvARW_NormCCRES: return true; default: return false; } } -//***************************************************************************** -MINOSRPA::MINOSRPA() { - - fApply_MINOSRPA = false; - - fDef_MINOSRPA_A = 1.010; - fTwk_MINOSRPA_A = 0.000; - fDef_MINOSRPA_B = 0.156; - fTwk_MINOSRPA_B = 0.000; - -} -// -double MINOSRPA::CalcWeight(BaseFitEvt* evt) { - if (!fTweaked) return 1.0; - double w = 1.0; - - // Extract the GENIE Record - GHepRecord* ghep = static_cast(evt->genie_event->event); - const Interaction* interaction = ghep->Summary(); - const InitialState& init_state = interaction->InitState(); - const ProcessInfo& proc_info = interaction->ProcInfo(); - const Target& tgt = init_state.Tgt(); - - // If not QE return 1.0 - if (!tgt.IsNucleus()) return 1.0; - if (!proc_info.IsQuasiElastic() && !proc_info.IsResonant()) return 1.0; - - // Extract Beam and Target PDG - GHepParticle* neutrino = ghep->Probe(); - // int bpdg = neutrino->Pdg(); - - GHepParticle* target = ghep->Particle(1); - assert(target); - // int tpdg = target->Pdg(); - - // Extract Q0-Q3 - GHepParticle* fsl = ghep->FinalStatePrimaryLepton(); - const TLorentzVector& k1 = *(neutrino->P4()); - const TLorentzVector& k2 = *(fsl->P4()); - //double q0 = fabs((k1 - k2).E()); - //double q3 = fabs((k1 - k2).Vect().Mag()); - double Q2 = fabs((k1 - k2).Mag2()); - - // Resonant Events - if (proc_info.IsResonant() && fApply_MINOSRPA) { - w *= GetRPAWeight( - Q2, - fCur_MINOSRPA_A, - fCur_MINOSRPA_B - ); - } - - return w; -} -// -double MINOSRPA::GetRPAWeight(double Q2, double A, double B) { - if (Q2 > 0.7) return 1.0; - double w = A / (1.0 + TMath::Exp(1.0 - TMath::Sqrt(Q2) / B)); - return w; -} -// -bool MINOSRPA::IsHandled(int rwenum) { - int curenum = rwenum % 1000; - switch (curenum) { - case Reweight::kMINERvARW_MINOSRPA_Apply: - case Reweight::kMINERvARW_MINOSRPA_A: - case Reweight::kMINERvARW_MINOSRPA_B: - return true; - default: - return false; - } -} -// -void MINOSRPA::SetDialValue(std::string name, double val) { - SetDialValue(Reweight::ConvDial(name, kCUSTOM), val); -} -// -void MINOSRPA::SetDialValue(int rwenum, double val) { - int curenum = rwenum % 1000; - - // Check Handled - if (!IsHandled(curenum)) return; - if (curenum == Reweight::kMINERvARW_MINOSRPA_Apply) fApply_MINOSRPA = (val > 0.5); - if (curenum == Reweight::kMINERvARW_MINOSRPA_A) fTwk_MINOSRPA_A = val; - if (curenum == Reweight::kMINERvARW_MINOSRPA_B) fTwk_MINOSRPA_B = val; - - // Check for changes - fTweaked = (fApply_MINOSRPA || - fabs(fTwk_MINOSRPA_A) > 0.0 || - fabs(fTwk_MINOSRPA_B) > 0.0); - - // Update Values - fCur_MINOSRPA_A = fDef_MINOSRPA_A * (1.0 + 0.1 * fTwk_MINOSRPA_A); - fCur_MINOSRPA_B = fDef_MINOSRPA_B * (1.0 + 0.1 * fTwk_MINOSRPA_B); - -} -// -// -//***************************************************************************** -LagrangeRPA::LagrangeRPA() { - - fApplyRPA = false; - - /* - fI1_Def = 4.18962e-01; - fI1 = fI1_Def; - - fI2_Def = 7.39927e-01; - fI2 = fI2_Def; - */ - - // Table VIII https://arxiv.org/pdf/1903.01558.pdf - fR1_Def = 0.37; - fR1 = fR1_Def; - - fR2_Def = 0.60; - fR2 = fR2_Def; - -} -// -double LagrangeRPA::CalcWeight(BaseFitEvt* evt) { - - if (!fTweaked) return 1.0; - double w = 1.0; - - // Extract the GENIE Record - GHepRecord* ghep = static_cast(evt->genie_event->event); - const Interaction* interaction = ghep->Summary(); - const InitialState& init_state = interaction->InitState(); - const ProcessInfo& proc_info = interaction->ProcInfo(); - const Target& tgt = init_state.Tgt(); - - // If not QE return 1.0 - if (!tgt.IsNucleus()) return 1.0; - if (!proc_info.IsQuasiElastic() && !proc_info.IsResonant()) return 1.0; - - // Extract Beam and Target PDG - GHepParticle* neutrino = ghep->Probe(); - // int bpdg = neutrino->Pdg(); - - GHepParticle* target = ghep->Particle(1); - assert(target); - // int tpdg = target->Pdg(); - - // Extract Q0-Q3 - GHepParticle* fsl = ghep->FinalStatePrimaryLepton(); - const TLorentzVector& k1 = *(neutrino->P4()); - const TLorentzVector& k2 = *(fsl->P4()); - //double q0 = fabs((k1 - k2).E()); - //double q3 = fabs((k1 - k2).Vect().Mag()); - double Q2 = fabs((k1 - k2).Mag2()); - - if (proc_info.IsResonant() && fApplyRPA) { - w *= GetRPAWeight(Q2); - } - - return w; -} -// -// -double LagrangeRPA::GetRPAWeight(double Q2) { - //std::cout << "Getting RPA Weight : " << Q2 << std::endl; - if (Q2 > 0.7) return 1.0; - - // Keep original Lagrange RPA for documentation - /* - double x1 = 0.00; - double x2 = 0.30; - double x3 = 0.70; - - double y1 = 0.00; - double y2 = fI2; - double y3 = 1.00; - - double xv = Q2; - - // Automatically 0 because y1 choice - double W1 = y1 * (xv-x2)*(xv-x3)/((x1-x2)*(x1-x3)); - double W2 = y2 * (xv-x1)*(xv-x3)/((x2-x1)*(x2-x3)); - double W3 = y3 * (xv-x1)*(xv-x2)/((x3-x1)*(x3-x2)); - - double P = W1 + W2 + W3; - double A1 = (1.0 - sqrt(1.0 - fI1)); - double R = P * (1.0 - A1) + A1; - - return 1.0 - (1.0-R)*(1.0-R); - */ - - // Equation 7 https://arxiv.org/pdf/1903.01558.pdf - const double x1 = 0.00; - const double x2 = 0.35; - const double x3 = 0.70; - - // Equation 6 https://arxiv.org/pdf/1903.01558.pdf - double RQ2 = fR2 *( (Q2-x1)*(Q2-x3)/((x2-x1)*(x2-x3)) ) - + (Q2-x1)*(Q2-x2)/((x3-x1)*(x3-x2)); - double weight = 1-(1-fR1)*(1-RQ2)*(1-RQ2); - - // Check range of R1 and R2 - // Commented out because this is implementation dependent: user may want strange double peaks - /* - if (fR1 > 1) return 1; - if (fR2 > 1 || fR2 < 0.5) return 1; - */ - - return weight; -} -// -bool LagrangeRPA::IsHandled(int rwenum) { - int curenum = rwenum % 1000; - switch (curenum) { - case Reweight::kMINERvARW_LagrangeRPA_Apply: - case Reweight::kMINERvARW_LagrangeRPA_R1: - case Reweight::kMINERvARW_LagrangeRPA_R2: - return true; - default: - return false; - } -} -// -void LagrangeRPA::SetDialValue(std::string name, double val) { - SetDialValue(Reweight::ConvDial(name, kCUSTOM), val); -} -// -void LagrangeRPA::SetDialValue(int rwenum, double val) { - int curenum = rwenum % 1000; - - // Check Handled - if (!IsHandled(curenum)) return; - if (curenum == Reweight::kMINERvARW_LagrangeRPA_Apply) fApplyRPA = (val > 0.5); - if (curenum == Reweight::kMINERvARW_LagrangeRPA_R1) fR1 = val; - if (curenum == Reweight::kMINERvARW_LagrangeRPA_R2) fR2 = val; - - // Check for changes - fTweaked = (fApplyRPA); -} -// - //******************************************************* RikRPA::RikRPA() { //******************************************************* // - Syst : kMINERvA_RikRPA_ApplyRPA // - Type : Binary // - Limits : 0.0 (false) -> 1.0 (true) // - Default : 0.0 fApplyDial_RPACorrection = false; // - Syst : kMINERvA_RikRPA_LowQ2 // - Type : Absolute // - Limits : 1.0 -> 1.0 // - Default : 0.0 // - Frac Error : 100% fDefDial_RPALowQ2 = 0.0; fCurDial_RPALowQ2 = fDefDial_RPALowQ2; fErrDial_RPALowQ2 = 0.0; // - Syst : kMINERvA_RikRPA_HighQ2 // - Type : Absolute // - Limits : 1.0 -> 1.0 // - Default : 0.0 // - Frac Error : 100% fDefDial_RPAHighQ2 = 0.0; fCurDial_RPAHighQ2 = fDefDial_RPAHighQ2; fErrDial_RPAHighQ2 = 1.0; // - Syst : kMINERvA_RikRESRPA_ApplyRPA // - Type : Binary // - Limits : 0.0 (false) -> 1.0 (true) // - Default : 0.0 fApplyDial_RESRPACorrection = false; // - Syst : kMINERvA_RikRESRPA_LowQ2 // - Type : Absolute // - Limits : 1.0 -> 1.0 // - Default : 0.0 // - Frac Error : 100% fDefDial_RESRPALowQ2 = 0.0; fCurDial_RESRPALowQ2 = fDefDial_RESRPALowQ2; fErrDial_RESRPALowQ2 = 0.0; // - Syst : kMINERvA_RikRESRPA_HighQ2 // - Type : Absolute // - Limits : 1.0 -> 1.0 // - Default : 0.0 // - Frac Error : 100% fDefDial_RESRPAHighQ2 = 0.0; fCurDial_RESRPAHighQ2 = fDefDial_RESRPAHighQ2; fErrDial_RESRPAHighQ2 = 1.0; // Setup Calculators fEventWeights = new double[5]; for (int i = 0; i < kMaxCalculators; i++) { fRPACalculators[i] = NULL; } fTweaked = false; } //******************************************************* RikRPA::~RikRPA() { //******************************************************* // delete fEventWeights; // for (size_t i = 0; i < kMaxCalculators; i++) { // if (fRPACalculators[i]) delete fRPACalculators[i]; // fRPACalculators[i] = NULL; // } } //******************************************************* double RikRPA::CalcWeight(BaseFitEvt *evt) { //******************************************************* // LOG(FIT) << "Calculating RikRPA" << std::endl; // Return 1.0 if not tweaked if (!fTweaked) return 1.0; double w = 1.0; // Extract the GENIE Record GHepRecord *ghep = static_cast(evt->genie_event->event); const Interaction *interaction = ghep->Summary(); const InitialState &init_state = interaction->InitState(); const ProcessInfo &proc_info = interaction->ProcInfo(); // const Kinematics & kine = interaction->Kine(); // const XclsTag & xcls = interaction->ExclTag(); const Target &tgt = init_state.Tgt(); // If not QE return 1.0 // LOG(FIT) << "RikRPA : Event QE = " << proc_info.IsQuasiElastic() << // std::endl; if (!tgt.IsNucleus()) { return 1.0; } if (!proc_info.IsQuasiElastic() && !proc_info.IsResonant()) return 1.0; // Extract Beam and Target PDG GHepParticle *neutrino = ghep->Probe(); int bpdg = neutrino->Pdg(); GHepParticle *target = ghep->TargetNucleus(); assert(target); int tpdg = target->Pdg(); // Find the enum we need int calcenum = GetRPACalcEnum(bpdg, tpdg); if (calcenum == -1) return 1.0; // Check we have the RPA Calc setup for this enum // if not, set it up at that point if (!fRPACalculators[calcenum]) SetupRPACalculator(calcenum); weightRPA *rpacalc = fRPACalculators[calcenum]; if (!rpacalc) { NUIS_ABORT("Failed to grab the RPA Calculator : " << calcenum); } // Extract Q0-Q3 GHepParticle *fsl = ghep->FinalStatePrimaryLepton(); const TLorentzVector &k1 = *(neutrino->P4()); const TLorentzVector &k2 = *(fsl->P4()); double q0 = fabs((k1 - k2).E()); double q3 = fabs((k1 - k2).Vect().Mag()); double Q2 = fabs((k1 - k2).Mag2()); // Quasielastic if (proc_info.IsQuasiElastic()) { // Now use q0-q3 and RPA Calculator to fill fWeights rpacalc->getWeight(q0, q3, fEventWeights); if (fApplyDial_RPACorrection) { w *= fEventWeights[0]; // CV } // Syst Application : kMINERvA_RikRPA_LowQ2 if (fabs(fCurDial_RPALowQ2) > 0.0) { double interpw = fEventWeights[0]; if (fCurDial_RPALowQ2 > 0.0 && Q2 < 2.0) { interpw = fEventWeights[0] - (fEventWeights[0] - fEventWeights[1]) * fCurDial_RPALowQ2; // WLow+ } else if (fCurDial_RPALowQ2 < 0.0 && Q2 < 2.0) { interpw = fEventWeights[0] - (fEventWeights[2] - fEventWeights[0]) * fCurDial_RPALowQ2; // WLow- } w *= interpw / fEventWeights[0]; // Div by CV again } // Syst Application : kMINERvA_RikRPA_HighQ2 if (fabs(fCurDial_RPAHighQ2) > 0.0) { double interpw = fEventWeights[0]; if (fCurDial_RPAHighQ2 > 0.0) { interpw = fEventWeights[0] - (fEventWeights[0] - fEventWeights[3]) * fCurDial_RPAHighQ2; // WHigh+ } else if (fCurDial_RPAHighQ2 < 0.0) { interpw = fEventWeights[0] - (fEventWeights[4] - fEventWeights[0]) * fCurDial_RPAHighQ2; // WHigh- } w *= interpw / fEventWeights[0]; // Div by CV again } } // Resonant Events if (proc_info.IsResonant()) { // Now use Q2 and RESRPA Calculator to fill fWeights double CV = rpacalc->getWeight(Q2); if (fApplyDial_RESRPACorrection) { w *= CV; // fEventWeights[0]; // CVa } /* // Syst Application : kMINERvA_RikRESRPA_LowQ2 if (fabs(fCurDial_RESRPAHighQ2) > 0.0) { double interpw = fEventWeights[0]; if (fCurDial_RESRPAHighQ2 > 0.0) { interpw = fEventWeights[0] - (fEventWeights[0] - fEventWeights[3]) * fCurDial_RESRPAHighQ2; // WHigh+ } else if (fCurDial_RESRPAHighQ2 < 0.0) { interpw = fEventWeights[0] - (fEventWeights[4] - fEventWeights[0]) * fCurDial_RESRPAHighQ2; // WHigh- } w *= interpw / fEventWeights[0]; // Div by CV again } // Syst Application : kMINERvA_RikRESRPA_HighQ2 if (fabs(fCurDial_RESRPAHighQ2) > 0.0) { double interpw = fEventWeights[0]; if (fCurDial_RESRPAHighQ2 > 0.0) { interpw = fEventWeights[0] - (fEventWeights[0] - fEventWeights[3]) * fCurDial_RESRPAHighQ2; // WHigh+ } else if (fCurDial_RESRPAHighQ2 < 0.0) { interpw = fEventWeights[0] - (fEventWeights[4] - fEventWeights[0]) * fCurDial_RESRPAHighQ2; // WHigh- } w *= interpw / fEventWeights[0]; // Div by CV again } */ } // LOG(FIT) << "RPA Weight = " << w << std::endl; return w; } // namespace reweight //******************************************************* void RikRPA::SetDialValue(std::string name, double val) { //******************************************************* SetDialValue(Reweight::ConvDial(name, kCUSTOM), val); } //******************************************************* void RikRPA::SetDialValue(int rwenum, double val) { //******************************************************* int curenum = rwenum % 1000; // Check Handled if (!IsHandled(curenum)) return; if (curenum == Reweight::kMINERvARW_RikRPA_ApplyRPA) fApplyDial_RPACorrection = (val > 0.5); if (curenum == Reweight::kMINERvARW_RikRPA_LowQ2) fCurDial_RPALowQ2 = val; if (curenum == Reweight::kMINERvARW_RikRPA_HighQ2) fCurDial_RPAHighQ2 = val; if (curenum == Reweight::kMINERvARW_RikRESRPA_ApplyRPA) fApplyDial_RESRPACorrection = (val > 0.5); if (curenum == Reweight::kMINERvARW_RikRESRPA_LowQ2) fCurDial_RESRPALowQ2 = val; if (curenum == Reweight::kMINERvARW_RikRESRPA_HighQ2) fCurDial_RESRPAHighQ2 = val; // Assign flag to say stuff has changed fTweaked = (fApplyDial_RPACorrection || fabs(fCurDial_RPAHighQ2 - fDefDial_RPAHighQ2) > 0.0 || fabs(fCurDial_RPALowQ2 - fDefDial_RPALowQ2) > 0.0 || fApplyDial_RESRPACorrection || fabs(fCurDial_RESRPAHighQ2 - fDefDial_RESRPAHighQ2) > 0.0 || fabs(fCurDial_RESRPALowQ2 - fDefDial_RESRPALowQ2) > 0.0); } //******************************************************* bool RikRPA::IsHandled(int rwenum) { //******************************************************* int curenum = rwenum % 1000; switch (curenum) { case Reweight::kMINERvARW_RikRESRPA_ApplyRPA: return true; case Reweight::kMINERvARW_RikRESRPA_LowQ2: return true; case Reweight::kMINERvARW_RikRESRPA_HighQ2: return true; case Reweight::kMINERvARW_RikRPA_ApplyRPA: return true; case Reweight::kMINERvARW_RikRPA_LowQ2: return true; case Reweight::kMINERvARW_RikRPA_HighQ2: return true; default: return false; } } //******************************************************* void RikRPA::SetupRPACalculator(int calcenum) { //******************************************************* std::string rwdir = FitPar::GetDataBase() + "reweight/MINERvA/RikRPA/"; std::string fidir = ""; switch (calcenum) { case kNuMuC12: fidir = "outNievesRPAratio-nu12C-20GeV-20170202.root"; break; case kNuMuO16: fidir = "outNievesRPAratio-nu16O-20GeV-20170202.root"; break; case kNuMuAr40: fidir = "outNievesRPAratio-nu40Ar-20GeV-20170202.root"; break; case kNuMuCa40: fidir = "outNievesRPAratio-nu40Ca-20GeV-20170202.root"; break; case kNuMuFe56: fidir = "outNievesRPAratio-nu56Fe-20GeV-20170202.root"; break; case kNuMuBarC12: fidir = "outNievesRPAratio-anu12C-20GeV-20170202.root"; break; case kNuMuBarO16: fidir = "outNievesRPAratio-anu16O-20GeV-20170202.root"; break; case kNuMuBarAr40: fidir = "outNievesRPAratio-anu40Ar-20GeV-20170202.root"; break; case kNuMuBarCa40: fidir = "outNievesRPAratio-anu40Ca-20GeV-20170202.root"; break; case kNuMuBarFe56: fidir = "outNievesRPAratio-anu56Fe-20GeV-20170202.root"; break; } NUIS_LOG(FIT, "Loading RPA CALC : " << fidir); TDirectory *olddir = gDirectory; NUIS_LOG(FIT, "***********************************************"); NUIS_LOG(FIT, "Loading a new weightRPA calculator"); NUIS_LOG(FIT, "Authors: Rik Gran, Heidi Schellman"); NUIS_LOG(FIT, "Citation: arXiv:1705.02932 [hep-ex]"); NUIS_LOG(FIT, "***********************************************"); // Test the file exists std::ifstream infile((rwdir + fidir).c_str()); if (!infile.good()) { NUIS_ERR(FTL, "*** ERROR ***"); NUIS_ERR(FTL, "RikRPA file " << rwdir + fidir << " does not exist!"); NUIS_ERR(FTL, "These can be found at https://nuisance.hepforge.org/files/RikRPA/"); NUIS_ERR(FTL, "Please run: wget -r -nH --cut-dirs=2 -np -e robots=off -R " "\"index.html*\" https://nuisance.hepforge.org/files/RikRPA/ -P " << rwdir); NUIS_ERR(FTL, "And try again"); NUIS_ABORT("*************"); } fRPACalculators[calcenum] = new weightRPA(rwdir + fidir); olddir->cd(); return; } //******************************************************* int RikRPA::GetRPACalcEnum(int bpdg, int tpdg) { //******************************************************* if (bpdg == 14 && tpdg == 1000060120) return kNuMuC12; else if (bpdg == 14 && tpdg == 1000080160) return kNuMuO16; else if (bpdg == 14 && tpdg == 1000180400) return kNuMuAr40; else if (bpdg == 14 && tpdg == 1000200400) return kNuMuCa40; else if (bpdg == 14 && tpdg == 1000280560) return kNuMuFe56; else if (bpdg == -14 && tpdg == 1000060120) return kNuMuBarC12; else if (bpdg == -14 && tpdg == 1000080160) return kNuMuBarO16; else if (bpdg == -14 && tpdg == 1000180400) return kNuMuBarAr40; else if (bpdg == -14 && tpdg == 1000200400) return kNuMuBarCa40; else if (bpdg == -14 && tpdg == 1000280560) return kNuMuBarFe56; else { // NUIS_ERR(WRN, "Unknown beam and target combination for RPA Calcs! " //<< bpdg << " " << tpdg); } return -1; } //***************************************************************************** COHBrandon::COHBrandon() { fApply_COHNorm = false; fDef_COHNorm = 0.5; fCur_COHNorm = fDef_COHNorm; fTwk_COHNorm = 0.0; fDef_COHCut = 0.450; fCur_COHCut = fDef_COHCut; fTwk_COHNorm = 0.0; } COHBrandon::~COHBrandon() {}; double COHBrandon::CalcWeight(BaseFitEvt* evt) { // Check GENIE if (evt->fType != kGENIE) return 1.0; // Extract the GENIE Record GHepRecord* ghep = static_cast(evt->genie_event->event); const Interaction* interaction = ghep->Summary(); //const InitialState& init_state = interaction->InitState(); const ProcessInfo& proc_info = interaction->ProcInfo(); //const Target& tgt = init_state.Tgt(); // If the event is not QE this Calc doesn't handle it if (!proc_info.IsCoherent()) return 1.0; // WEIGHT CALCULATIONS ------------- double w = 1.0; double pionE = -999.9; TObjArrayIter piter(ghep); GHepParticle * p = 0; //int ip = -1; while ( (p = (GHepParticle *) piter.Next()) ) { int pdgc = p->Pdg(); int ist = p->Status(); // if (ghep->Particle(ip)->FirstMother()==0) continue; if (pdg::IsPseudoParticle(p->Pdg())) continue; if (ist == kIStStableFinalState) { if (pdgc == 211 || pdgc == -211 || pdgc == 111) { pionE = p->Energy(); break; } } } // std::cout << "Coherent Pion Energy : " << pionE << std::endl; // Apply Weight if (fApply_COHNorm && pionE > 0.0 && pionE < fCur_COHCut) { w *= fCur_COHNorm; } // Return Combined Weight return w; } void COHBrandon::SetDialValue(std::string name, double val) { SetDialValue(Reweight::ConvDial(name, kCUSTOM), val); } void COHBrandon::SetDialValue(int rwenum, double val) { // Check Handled int curenum = rwenum % 1000; if (!IsHandled(curenum)) return; // Set Values if (curenum == Reweight::kMINERvARW_NormCOH) { fTwk_COHNorm = val; fCur_COHNorm = fDef_COHNorm * (1.0 + 0.1 * fTwk_COHNorm); } if (curenum == Reweight::kMINERvARW_CutCOH) { fTwk_COHCut = val; fCur_COHCut = fDef_COHCut * (1.0 + 0.1 * fTwk_COHCut); } if (curenum == Reweight::kMINERvARW_ApplyCOH) { fApply_COHNorm = (val > 0.5); } // Define Tweaked fTweaked = (fApply_COHNorm); } bool COHBrandon::IsHandled(int rwenum) { int curenum = rwenum % 1000; switch (curenum) { case Reweight::kMINERvARW_NormCOH: case Reweight::kMINERvARW_CutCOH: case Reweight::kMINERvARW_ApplyCOH: return true; default: return false; } } //***************************************************************************** WEnhancement::WEnhancement() { fApply_Enhancement = false; fDef_WNorm = 1.0; fCur_WNorm = fDef_WNorm; fTwk_WNorm = 0.0; fDef_WMean = 0.95; fCur_WMean = fDef_WMean; fTwk_WMean = 0.0; fDef_WSigma = 0.225; fCur_WSigma = fDef_WSigma; fTwk_WSigma = 0.0; } WEnhancement::~WEnhancement() {}; double WEnhancement::CalcWeight(BaseFitEvt* evt) { // Check GENIE if (evt->fType != kGENIE) return 1.0; // Extract the GENIE Record GHepRecord* ghep = static_cast(evt->genie_event->event); const Interaction* interaction = ghep->Summary(); //const Kinematics & kine = interaction->Kine(); //const InitialState& init_state = interaction->InitState(); const ProcessInfo& proc_info = interaction->ProcInfo(); //const Target& tgt = init_state.Tgt(); // If the event is not QE this Calc doesn't handle it if (!proc_info.IsWeakCC()) return 1.0; if (!proc_info.IsResonant()) return 1.0; if (!fApply_Enhancement) return 1.0; // WEIGHT CALCULATIONS ------------- double w = 1.0; bool isCC1pi0AtVertex = false; // Get W GHepParticle* neutrino = ghep->Probe(); GHepParticle* fsl = ghep->FinalStatePrimaryLepton(); const TLorentzVector& k1 = *(neutrino->P4()); const TLorentzVector& k2 = *(fsl->P4()); double E_mu = k2.E(); double p_mu = k2.Vect().Mag(); double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu); double th_nu_mu = k1.Vect().Angle(k2.Vect()); // The factor of 1000 is necessary for downstream functions const double m_p = PhysConst::mass_proton; double E_nu = k1.E(); // double hadMass = kine.W (true); double hadMass = sqrt(m_p * m_p + m_mu * m_mu - 2 * m_p * E_mu + \ 2 * E_nu * (m_p - E_mu + p_mu * cos(th_nu_mu))); // Determine if event is CC1pi0 at vertex TObjArrayIter piter(ghep); GHepParticle * p = 0; int pi0 = 0; int piother = 0; while ( (p = (GHepParticle *) piter.Next()) ) { int pdgc = p->Pdg(); int ist = p->Status(); if (pdg::IsPseudoParticle(p->Pdg())) continue; if (ist == genie::kIStStableFinalState) { if (pdgc == 111) { pi0++; } else if (pdgc == 211 || pdgc == -211) { piother++; } } } // std::cout << "Npi0 = " << pi0 << std::endl; isCC1pi0AtVertex = (pi0 == 1 && piother == 0); // Apply Weight if (fApply_Enhancement && isCC1pi0AtVertex) { double enhancement = 1.0 + fCur_WNorm * (exp( -0.5 * pow((fCur_WMean - hadMass) / (fCur_WSigma), 2) )); w *= enhancement; } // Return Combined Weight if (isnan(w) || w < 0.0 || w > 400.0) { w = 1.0; } return w; } void WEnhancement::SetDialValue(std::string name, double val) { SetDialValue(Reweight::ConvDial(name, kCUSTOM), val); } void WEnhancement::SetDialValue(int rwenum, double val) { // Check Handled int curenum = rwenum % 1000; if (!IsHandled(curenum)) return; // Set Values if (curenum == Reweight::kMINERvARW_ApplyWTune) { fApply_Enhancement = (val > 0.5); } if (curenum == Reweight::kMINERvARW_NormWTune) { fTwk_WNorm = val; fCur_WNorm = fDef_WNorm * (1.0 + 0.1 * fTwk_WNorm); } if (curenum == Reweight::kMINERvARW_MeanWTune) { fTwk_WMean = val; fCur_WMean = fDef_WMean * (1.0 + 0.1 * fTwk_WMean); } if (curenum == Reweight::kMINERvARW_SigmaWTune) { fTwk_WSigma = val; fCur_WSigma = fDef_WSigma * (1.0 + 0.1 * fTwk_WSigma); } // Define Tweaked fTweaked = (fApply_Enhancement); } bool WEnhancement::IsHandled(int rwenum) { int curenum = rwenum % 1000; switch (curenum) { case Reweight::kMINERvARW_ApplyWTune: case Reweight::kMINERvARW_NormWTune: case Reweight::kMINERvARW_MeanWTune: case Reweight::kMINERvARW_SigmaWTune: return true; default: return false; } } } // namespace reweight } // namespace nuisance #endif #endif diff --git a/src/Reweight/MINERvAWeightCalcs.h b/src/Reweight/MINERvAWeightCalcs.h index 0daef1e..9c6efe0 100644 --- a/src/Reweight/MINERvAWeightCalcs.h +++ b/src/Reweight/MINERvAWeightCalcs.h @@ -1,259 +1,209 @@ #ifndef MINERVA_WEIGHT_CALCS #define MINERVA_WEIGHT_CALCS #include #ifdef __MINERVA_RW_ENABLED__ #ifdef __GENIE_ENABLED__ #ifdef GENIE_PRE_R3 #include "Conventions/Units.h" #include "EVGCore/EventRecord.h" #include "GHEP/GHepParticle.h" #include "GHEP/GHepRecord.h" #include "GHEP/GHepUtils.h" #include "Ntuple/NtpMCEventRecord.h" #include "PDG/PDGUtils.h" #else #include "Framework/Conventions/Units.h" #include "Framework/EventGen/EventRecord.h" #include "Framework/GHEP/GHepParticle.h" #include "Framework/GHEP/GHepRecord.h" #include "Framework/GHEP/GHepUtils.h" #include "Framework/Ntuple/NtpMCEventRecord.h" #include "Framework/ParticleData/PDGUtils.h" #endif #include "NUISANCEWeightCalcs.h" #include "GeneralUtils.h" #include "NUISANCESyst.h" #include "FitEvent.h" #include "WeightUtils.h" #include "weightRPA.h" using namespace genie; class BaseFitEvt; namespace nuisance { namespace reweight { // QE Dials class MINERvAReWeight_QE : public NUISANCEWeightCalc { public: MINERvAReWeight_QE(); virtual ~MINERvAReWeight_QE(); double CalcWeight(BaseFitEvt* evt); void SetDialValue(std::string name, double val); void SetDialValue(int rwenum, double val); bool IsHandled(int rwenum); double fTwk_NormCCQE; double fCur_NormCCQE; double fDef_NormCCQE; bool fTweaked; }; // MEC Dials class MINERvAReWeight_MEC : public NUISANCEWeightCalc { public: MINERvAReWeight_MEC(); virtual ~MINERvAReWeight_MEC(); double CalcWeight(BaseFitEvt* evt); void SetDialValue(std::string name, double val); void SetDialValue(int rwenum, double val); bool IsHandled(int rwenum); double fTwk_NormCCMEC; double fCur_NormCCMEC; double fDef_NormCCMEC; bool fTweaked; }; // RES Dials class MINERvAReWeight_RES : public NUISANCEWeightCalc { public: MINERvAReWeight_RES(); virtual ~MINERvAReWeight_RES(); double CalcWeight(BaseFitEvt* evt); void SetDialValue(std::string name, double val); void SetDialValue(int rwenum, double val); bool IsHandled(int rwenum); double fTwk_NormCCRES; double fCur_NormCCRES; double fDef_NormCCRES; bool fTweaked; }; - // MINOS pion tuning, https://arxiv.org/pdf/1903.01558.pdf - class MINOSRPA : public NUISANCEWeightCalc { - public: - MINOSRPA(); - ~MINOSRPA(){}; - - double CalcWeight(BaseFitEvt* evt); - void SetDialValue(std::string name, double val); - void SetDialValue(int rwenum, double val); - bool IsHandled(int rwenum); - - double GetRPAWeight(double Q2, double A, double B); - - bool fTweaked; - - bool fApply_MINOSRPA; - - double fTwk_MINOSRPA_A; - double fDef_MINOSRPA_A; - double fCur_MINOSRPA_A; - - double fTwk_MINOSRPA_B; - double fDef_MINOSRPA_B; - double fCur_MINOSRPA_B; - - }; - - // MINERvA pion tuning, https://arxiv.org/pdf/1903.01558.pdf - class LagrangeRPA : public NUISANCEWeightCalc { - public: - LagrangeRPA(); - ~LagrangeRPA(){}; - - double CalcWeight(BaseFitEvt* evt); - void SetDialValue(std::string name, double val); - void SetDialValue(int rwenum, double val); - bool IsHandled(int rwenum); - - double GetRPAWeight(double Q2); - - bool fTweaked; - - bool fApplyRPA; - - double fR1; - double fR2; - double fR1_Def; - double fR2_Def; - }; - /// RPA Weight Calculator that applies RPA systematics /// to GENIE events. GENIE EVENTS ONLY! class RikRPA : public NUISANCEWeightCalc { public: RikRPA(); ~RikRPA(); double CalcWeight(BaseFitEvt* evt); void SetDialValue(std::string name, double val); void SetDialValue(int rwenum, double val); bool IsHandled(int rwenum); void SetupRPACalculator(int calcenum); int GetRPACalcEnum(int bpdg, int tpdg); bool fApplyDial_RPACorrection; double fTwkDial_RPALowQ2; double fDefDial_RPALowQ2; double fCurDial_RPALowQ2; double fErrDial_RPALowQ2; double fTwkDial_RPAHighQ2; double fDefDial_RPAHighQ2; double fCurDial_RPAHighQ2; double fErrDial_RPAHighQ2; bool fApplyDial_RESRPACorrection; double fTwkDial_RESRPALowQ2; double fDefDial_RESRPALowQ2; double fCurDial_RESRPALowQ2; double fErrDial_RESRPALowQ2; double fTwkDial_RESRPAHighQ2; double fDefDial_RESRPAHighQ2; double fCurDial_RESRPAHighQ2; double fErrDial_RESRPAHighQ2; double* fEventWeights; bool fTweaked; const static int kMaxCalculators = 10; enum rpacalcenums { kNuMuC12, kNuMuO16, kNuMuAr40, kNuMuCa40, kNuMuFe56, kNuMuBarC12, kNuMuBarO16, kNuMuBarAr40, kNuMuBarCa40, kNuMuBarFe56 }; weightRPA* fRPACalculators[kMaxCalculators]; }; // Custom coherent tune from MINERvA class COHBrandon : public NUISANCEWeightCalc { public: COHBrandon(); ~COHBrandon(); double CalcWeight(BaseFitEvt* evt); void SetDialValue(std::string name, double val); void SetDialValue(int rwenum, double val); bool IsHandled(int rwenum); bool fApply_COHNorm; double fDef_COHNorm; double fCur_COHNorm; double fTwk_COHNorm; double fDef_COHCut; double fCur_COHCut; double fTwk_COHCut; bool fTweaked; }; class WEnhancement : public NUISANCEWeightCalc { public: WEnhancement(); ~WEnhancement(); double CalcWeight(BaseFitEvt* evt); void SetDialValue(std::string name, double val); void SetDialValue(int rwenum, double val); bool IsHandled(int rwenum); bool fTweaked; bool fApply_Enhancement; double fDef_WNorm; double fCur_WNorm; double fTwk_WNorm; double fDef_WMean; double fCur_WMean; double fTwk_WMean; double fDef_WSigma; double fCur_WSigma; double fTwk_WSigma; }; }; // namespace reweight }; // namespace nuisance #endif // __GENIE_ENABLED__ #endif //__MINERVA_RW_ENABLED__ #endif diff --git a/src/Reweight/NUISANCEWeightCalcs.cxx b/src/Reweight/NUISANCEWeightCalcs.cxx index 9ae3551..e523b9b 100644 --- a/src/Reweight/NUISANCEWeightCalcs.cxx +++ b/src/Reweight/NUISANCEWeightCalcs.cxx @@ -1,545 +1,794 @@ #include "NUISANCEWeightCalcs.h" #include "FitEvent.h" #include "GeneralUtils.h" #include "NUISANCESyst.h" #include "WeightUtils.h" using namespace Reweight; ModeNormCalc::ModeNormCalc() { fNormRES = 1.0; } double ModeNormCalc::CalcWeight(BaseFitEvt *evt) { int mode = abs(evt->Mode); double w = 1.0; if (mode == 11 or mode == 12 or mode == 13) { w *= fNormRES; } return w; } void ModeNormCalc::SetDialValue(std::string name, double val) { SetDialValue(Reweight::ConvDial(name, kCUSTOM), val); } void ModeNormCalc::SetDialValue(int rwenum, double val) { int curenum = rwenum % 1000; // Check Handled if (!IsHandled(curenum)) return; if (curenum == kModeNorm_NormRES) fNormRES = val; } bool ModeNormCalc::IsHandled(int rwenum) { int curenum = rwenum % 1000; switch (curenum) { case kModeNorm_NormRES: return true; default: return false; } } +//***************************************************************************** +MINOSRPA::MINOSRPA() { + + fApply_MINOSRPA = false; + + fDef_MINOSRPA_A = 1.010; + fTwk_MINOSRPA_A = 0.000; + fDef_MINOSRPA_B = 0.156; + fTwk_MINOSRPA_B = 0.000; + +} +// +double MINOSRPA::CalcWeight(BaseFitEvt* evt) { + if (!fTweaked || !fApply_MINOSRPA) return 1.0; + + double w = 1.0; + + // If GENIE is enabled, use old code +#ifdef __GENIE_ENABLED__ + // Extract the GENIE Record + GHepRecord* ghep = static_cast(evt->genie_event->event); + const Interaction* interaction = ghep->Summary(); + const InitialState& init_state = interaction->InitState(); + const ProcessInfo& proc_info = interaction->ProcInfo(); + const Target& tgt = init_state.Tgt(); + + // If not QE return 1.0 + if (!tgt.IsNucleus()) return 1.0; + if (!proc_info.IsQuasiElastic() && !proc_info.IsResonant()) return 1.0; + + // Extract Beam and Target PDG + GHepParticle* neutrino = ghep->Probe(); + // int bpdg = neutrino->Pdg(); + + GHepParticle* target = ghep->Particle(1); + assert(target); + // int tpdg = target->Pdg(); + + // Extract Q0-Q3 + GHepParticle* fsl = ghep->FinalStatePrimaryLepton(); + const TLorentzVector& k1 = *(neutrino->P4()); + const TLorentzVector& k2 = *(fsl->P4()); + //double q0 = fabs((k1 - k2).E()); + //double q3 = fabs((k1 - k2).Vect().Mag()); + double Q2 = fabs((k1 - k2).Mag2()); + + w *= GetRPAWeight(Q2); +#else + // Get the Q2 from NUISANCE if not GENIE + FitEvent *fevt = static_cast(evt); + // Check the event is resonant + if (!fevt->IsResonant()) return 1.0; + int targeta = fevt->GetTargetA(); + int targetz = fevt->GetTargetZ(); + // Apply only to nuclear targets, ignore free protons + if (targeta == 1 || targetz == 1) return 1.0; + // Q2 in GeV2 + double Q2 = fevt->GetQ2(); + w *= GetRPAWeight(Q2); +#endif + + return w; +} + +// Do the actual weight calculation +double MINOSRPA::GetRPAWeight(double Q2) { + if (Q2 > 0.7) return 1.0; + double w = fCur_MINOSRPA_A / (1.0 + TMath::Exp(1.0 - TMath::Sqrt(Q2) / fCur_MINOSRPA_B)); + return w; +} + +bool MINOSRPA::IsHandled(int rwenum) { + int curenum = rwenum % 1000; + switch (curenum) { + case Reweight::kMINERvARW_MINOSRPA_Apply: + case Reweight::kMINERvARW_MINOSRPA_A: + case Reweight::kMINERvARW_MINOSRPA_B: + return true; + default: + return false; + } +} +// +void MINOSRPA::SetDialValue(std::string name, double val) { + SetDialValue(Reweight::ConvDial(name, kCUSTOM), val); +} +// +void MINOSRPA::SetDialValue(int rwenum, double val) { + int curenum = rwenum % 1000; + + // Check Handled + if (!IsHandled(curenum)) return; + if (curenum == Reweight::kMINERvARW_MINOSRPA_Apply) fApply_MINOSRPA = (val > 0.5); + if (curenum == Reweight::kMINERvARW_MINOSRPA_A) fTwk_MINOSRPA_A = val; + if (curenum == Reweight::kMINERvARW_MINOSRPA_B) fTwk_MINOSRPA_B = val; + + // Check for changes + fTweaked = (fApply_MINOSRPA || + fabs(fTwk_MINOSRPA_A) > 0.0 || + fabs(fTwk_MINOSRPA_B) > 0.0); + + // Update Values + fCur_MINOSRPA_A = fDef_MINOSRPA_A * (1.0 + 0.1 * fTwk_MINOSRPA_A); + fCur_MINOSRPA_B = fDef_MINOSRPA_B * (1.0 + 0.1 * fTwk_MINOSRPA_B); +} + +// +// +//***************************************************************************** +LagrangeRPA::LagrangeRPA() { + + fApplyRPA = false; + + /* + fI1_Def = 4.18962e-01; + fI1 = fI1_Def; + + fI2_Def = 7.39927e-01; + fI2 = fI2_Def; + */ + + // Table VIII https://arxiv.org/pdf/1903.01558.pdf + fR1_Def = 0.37; + fR1 = fR1_Def; + + fR2_Def = 0.60; + fR2 = fR2_Def; + +} +// +double LagrangeRPA::CalcWeight(BaseFitEvt* evt) { + + if (!fTweaked || !fApplyRPA) return 1.0; + double w = 1.0; + + // If GENIE is enabled, use old code +#ifdef __GENIE_ENABLED__ + // Extract the GENIE Record + GHepRecord* ghep = static_cast(evt->genie_event->event); + const Interaction* interaction = ghep->Summary(); + const InitialState& init_state = interaction->InitState(); + const ProcessInfo& proc_info = interaction->ProcInfo(); + const Target& tgt = init_state.Tgt(); + + // If not QE return 1.0 + if (!tgt.IsNucleus()) return 1.0; + if (!proc_info.IsQuasiElastic() && !proc_info.IsResonant()) return 1.0; + + // Extract Beam and Target PDG + GHepParticle* neutrino = ghep->Probe(); + // int bpdg = neutrino->Pdg(); + + GHepParticle* target = ghep->Particle(1); + assert(target); + // int tpdg = target->Pdg(); + + // Extract Q0-Q3 + GHepParticle* fsl = ghep->FinalStatePrimaryLepton(); + const TLorentzVector& k1 = *(neutrino->P4()); + const TLorentzVector& k2 = *(fsl->P4()); + //double q0 = fabs((k1 - k2).E()); + //double q3 = fabs((k1 - k2).Vect().Mag()); + double Q2 = fabs((k1 - k2).Mag2()); + w *= GetRPAWeight(Q2); +#else + // Get the Q2 from NUISANCE if not GENIE + FitEvent *fevt = static_cast(evt); + // Check the event is resonant + if (!fevt->IsResonant()) return 1.0; + int targeta = fevt->GetTargetA(); + int targetz = fevt->GetTargetZ(); + // Apply only to nuclear targets, ignore free protons + if (targeta == 1 || targetz == 1) return 1.0; + // Q2 in GeV2 + double Q2 = fevt->GetQ2(); + w *= GetRPAWeight(Q2); +#endif + + return w; +} +// +// +double LagrangeRPA::GetRPAWeight(double Q2) { + //std::cout << "Getting RPA Weight : " << Q2 << std::endl; + if (Q2 > 0.7) return 1.0; + + // Keep original Lagrange RPA for documentation + /* + double x1 = 0.00; + double x2 = 0.30; + double x3 = 0.70; + + double y1 = 0.00; + double y2 = fI2; + double y3 = 1.00; + + double xv = Q2; + + // Automatically 0 because y1 choice + double W1 = y1 * (xv-x2)*(xv-x3)/((x1-x2)*(x1-x3)); + double W2 = y2 * (xv-x1)*(xv-x3)/((x2-x1)*(x2-x3)); + double W3 = y3 * (xv-x1)*(xv-x2)/((x3-x1)*(x3-x2)); + + double P = W1 + W2 + W3; + double A1 = (1.0 - sqrt(1.0 - fI1)); + double R = P * (1.0 - A1) + A1; + + return 1.0 - (1.0-R)*(1.0-R); + */ + + // Equation 7 https://arxiv.org/pdf/1903.01558.pdf + const double x1 = 0.00; + const double x2 = 0.35; + const double x3 = 0.70; + + // Equation 6 https://arxiv.org/pdf/1903.01558.pdf + double RQ2 = fR2 *( (Q2-x1)*(Q2-x3)/((x2-x1)*(x2-x3)) ) + + (Q2-x1)*(Q2-x2)/((x3-x1)*(x3-x2)); + double weight = 1-(1-fR1)*(1-RQ2)*(1-RQ2); + + // Check range of R1 and R2 + // Commented out because this is implementation dependent: user may want strange double peaks + /* + if (fR1 > 1) return 1; + if (fR2 > 1 || fR2 < 0.5) return 1; + */ + + return weight; +} +// +bool LagrangeRPA::IsHandled(int rwenum) { + int curenum = rwenum % 1000; + switch (curenum) { + case Reweight::kMINERvARW_LagrangeRPA_Apply: + case Reweight::kMINERvARW_LagrangeRPA_R1: + case Reweight::kMINERvARW_LagrangeRPA_R2: + return true; + default: + return false; + } +} +// +void LagrangeRPA::SetDialValue(std::string name, double val) { + SetDialValue(Reweight::ConvDial(name, kCUSTOM), val); +} +// +void LagrangeRPA::SetDialValue(int rwenum, double val) { + int curenum = rwenum % 1000; + + // Check Handled + if (!IsHandled(curenum)) return; + if (curenum == Reweight::kMINERvARW_LagrangeRPA_Apply) fApplyRPA = (val > 0.5); + if (curenum == Reweight::kMINERvARW_LagrangeRPA_R1) fR1 = val; + if (curenum == Reweight::kMINERvARW_LagrangeRPA_R2) fR2 = val; + + // Check for changes + fTweaked = (fApplyRPA); +} +// + + BeRPACalc::BeRPACalc() - : fBeRPA_A(0.59), fBeRPA_B(1.05), fBeRPA_D(1.13), fBeRPA_E(0.88), - fBeRPA_U(1.2), nParams(0) { - // A = 0.59 +/- 20% - // B = 1.05 +/- 20% - // D = 1.13 +/- 15% - // E = 0.88 +/- 40% - // U = 1.2 + : fBeRPA_A(0.59), fBeRPA_B(1.05), fBeRPA_D(1.13), fBeRPA_E(0.88), + fBeRPA_U(1.2), nParams(0) { + // A = 0.59 +/- 20% + // B = 1.05 +/- 20% + // D = 1.13 +/- 15% + // E = 0.88 +/- 40% + // U = 1.2 } double BeRPACalc::CalcWeight(BaseFitEvt *evt) { - FitEvent *fevt = static_cast(evt); int mode = abs(evt->Mode); double w = 1.0; if (nParams == 0) { return w; } // Get Q2 // Get final state lepton if (mode == 1) { - double Q2 = - -1.0 * (fevt->GetHMFSAnyLeptons()->P4() - fevt->GetNeutrinoIn()->P4()) * - (fevt->GetHMFSAnyLeptons()->P4() - fevt->GetNeutrinoIn()->P4()) / 1.E6; + FitEvent *fevt = static_cast(evt); + double Q2 = fevt->GetQ2(); // Only CCQE events w *= calcRPA(Q2, fBeRPA_A, fBeRPA_B, fBeRPA_D, fBeRPA_E, fBeRPA_U); } return w; } void BeRPACalc::SetDialValue(std::string name, double val) { SetDialValue(Reweight::ConvDial(name, kCUSTOM), val); } void BeRPACalc::SetDialValue(int rwenum, double val) { int curenum = rwenum % 1000; // Check Handled if (!IsHandled(curenum)) return; // Need 4 or 5 reconfigures if (curenum == kBeRPA_A) fBeRPA_A = val; else if (curenum == kBeRPA_B) fBeRPA_B = val; else if (curenum == kBeRPA_D) fBeRPA_D = val; else if (curenum == kBeRPA_E) fBeRPA_E = val; else if (curenum == kBeRPA_U) fBeRPA_U = val; nParams++; } bool BeRPACalc::IsHandled(int rwenum) { int curenum = rwenum % 1000; switch (curenum) { - case kBeRPA_A: - case kBeRPA_B: - case kBeRPA_D: - case kBeRPA_E: - case kBeRPA_U: - return true; - default: - return false; + case kBeRPA_A: + case kBeRPA_B: + case kBeRPA_D: + case kBeRPA_E: + case kBeRPA_U: + return true; + default: + return false; } } SBLOscWeightCalc::SBLOscWeightCalc() { fDistance = 0.0; fMassSplitting = 0.0; fSin2Theta = 0.0; } double SBLOscWeightCalc::CalcWeight(BaseFitEvt *evt) { FitEvent *fevt = static_cast(evt); FitParticle *pnu = fevt->PartInfo(0); double E = pnu->fP.E() / 1.E3; // Extract energy return GetSBLOscWeight(E); } void SBLOscWeightCalc::SetDialValue(std::string name, double val) { SetDialValue(Reweight::ConvDial(name, kCUSTOM), val); } void SBLOscWeightCalc::SetDialValue(int rwenum, double val) { int curenum = rwenum % 1000; if (!IsHandled(curenum)) return; if (curenum == kSBLOsc_Distance) fDistance = val; if (curenum == kSBLOsc_MassSplitting) fMassSplitting = val; if (curenum == kSBLOsc_Sin2Theta) fSin2Theta = val; } bool SBLOscWeightCalc::IsHandled(int rwenum) { int curenum = rwenum % 1000; switch (curenum) { - case kSBLOsc_Distance: - return true; - case kSBLOsc_MassSplitting: - return true; - case kSBLOsc_Sin2Theta: - return true; - default: - return false; + case kSBLOsc_Distance: + return true; + case kSBLOsc_MassSplitting: + return true; + case kSBLOsc_Sin2Theta: + return true; + default: + return false; } } double SBLOscWeightCalc::GetSBLOscWeight(double E) { if (E <= 0.0) return 1.0 - 0.5 * fSin2Theta; return 1.0 - fSin2Theta * pow(sin(1.267 * fMassSplitting * fDistance / E), 2); } GaussianModeCorr::GaussianModeCorr() { // Apply the tilt-shift Gauss by Patrick // Alternatively set in config fMethod = true; // Init fApply_CCQE = false; fGausVal_CCQE[kPosNorm] = 0.0; fGausVal_CCQE[kPosTilt] = 0.0; fGausVal_CCQE[kPosPq0] = 1.0; fGausVal_CCQE[kPosWq0] = 1.0; fGausVal_CCQE[kPosPq3] = 1.0; fGausVal_CCQE[kPosWq3] = 1.0; fApply_2p2h = false; fGausVal_2p2h[kPosNorm] = 0.0; fGausVal_2p2h[kPosTilt] = 0.0; fGausVal_2p2h[kPosPq0] = 1.0; fGausVal_2p2h[kPosWq0] = 1.0; fGausVal_2p2h[kPosPq3] = 1.0; fGausVal_2p2h[kPosWq3] = 1.0; fApply_2p2h_PPandNN = false; fGausVal_2p2h_PPandNN[kPosNorm] = 0.0; fGausVal_2p2h_PPandNN[kPosTilt] = 0.0; fGausVal_2p2h_PPandNN[kPosPq0] = 1.0; fGausVal_2p2h_PPandNN[kPosWq0] = 1.0; fGausVal_2p2h_PPandNN[kPosPq3] = 1.0; fGausVal_2p2h_PPandNN[kPosWq3] = 1.0; fApply_2p2h_NP = false; fGausVal_2p2h_NP[kPosNorm] = 0.0; fGausVal_2p2h_NP[kPosTilt] = 0.0; fGausVal_2p2h_NP[kPosPq0] = 1.0; fGausVal_2p2h_NP[kPosWq0] = 1.0; fGausVal_2p2h_NP[kPosPq3] = 1.0; fGausVal_2p2h_NP[kPosWq3] = 1.0; fApply_CC1pi = false; fGausVal_CC1pi[kPosNorm] = 0.0; fGausVal_CC1pi[kPosTilt] = 0.0; fGausVal_CC1pi[kPosPq0] = 1.0; fGausVal_CC1pi[kPosWq0] = 1.0; fGausVal_CC1pi[kPosPq3] = 1.0; fGausVal_CC1pi[kPosWq3] = 1.0; fAllowSuppression = false; fDebugStatements = FitPar::Config().GetParB("GaussianModeCorr_DEBUG"); // fDebugStatements = true; } double GaussianModeCorr::CalcWeight(BaseFitEvt *evt) { FitEvent *fevt = static_cast(evt); double rw_weight = 1.0; // Get Neutrino if (!fevt->Npart()) { NUIS_ABORT("NO particles found in stack!"); } - FitParticle *pnu = fevt->GetHMISAnyLeptons(); + FitParticle *pnu = fevt->GetNeutrinoIn(); if (!pnu) { NUIS_ABORT("NO Starting particle found in stack!"); } - int pdgnu = pnu->fPID; - - int expect_fsleppdg = 0; - - if (pdgnu & 1) { - expect_fsleppdg = pdgnu; - } else { - expect_fsleppdg = abs(pdgnu) - 1; - } - FitParticle *plep = fevt->GetHMFSParticle(expect_fsleppdg); - if (!plep) - return 1.0; + FitParticle *plep = fevt->GetLeptonOut(); + if (!plep) return 1.0; TLorentzVector q = pnu->fP - plep->fP; // Extra q0,q3 double q0 = fabs(q.E()) / 1.E3; double q3 = fabs(q.Vect().Mag()) / 1.E3; int initialstate = -1; // Undef if (abs(fevt->Mode) == 2) { int npr = 0; int nne = 0; for (UInt_t j = 0; j < fevt->Npart(); j++) { if ((fevt->PartInfo(j))->fIsAlive) continue; if (fevt->PartInfo(j)->fPID == 2212) npr++; else if (fevt->PartInfo(j)->fPID == 2112) nne++; } if (fevt->Mode == 2 && npr == 1 && nne == 1) { initialstate = 2; } else if (fevt->Mode == 2 && ((npr == 0 && nne == 2) || (npr == 2 && nne == 0))) { initialstate = 1; } } // Apply weighting if (fApply_CCQE && abs(fevt->Mode) == 1) { if (fDebugStatements) std::cout << "Getting CCQE Weight" << std::endl; double g = GetGausWeight(q0, q3, fGausVal_CCQE); if (g < 1.0) g = 1.0; rw_weight *= g; } if (fApply_2p2h && abs(fevt->Mode) == 2) { if (fDebugStatements) std::cout << "Getting 2p2h Weight" << std::endl; if (fDebugStatements) std::cout << "Got q0 q3 = " << q0 << " " << q3 << " mode = " << fevt->Mode << std::endl; rw_weight *= GetGausWeight(q0, q3, fGausVal_2p2h); if (fDebugStatements) std::cout << "Returning Weight " << rw_weight << std::endl; } if (fApply_2p2h_PPandNN && abs(fevt->Mode) == 2 && initialstate == 1) { if (fDebugStatements) std::cout << "Getting 2p2h PPandNN Weight" << std::endl; rw_weight *= GetGausWeight(q0, q3, fGausVal_2p2h_PPandNN); } if (fApply_2p2h_NP && abs(fevt->Mode) == 2 && initialstate == 2) { if (fDebugStatements) std::cout << "Getting 2p2h NP Weight" << std::endl; rw_weight *= GetGausWeight(q0, q3, fGausVal_2p2h_NP); } if (fApply_CC1pi && abs(fevt->Mode) >= 11 && abs(fevt->Mode) <= 13) { if (fDebugStatements) std::cout << "Getting CC1pi Weight" << std::endl; rw_weight *= GetGausWeight(q0, q3, fGausVal_CC1pi); } return rw_weight; } void GaussianModeCorr::SetMethod(bool method) { fMethod = method; if (fMethod == true) { NUIS_LOG(FIT, - " Using tilt-shift Gaussian parameters for Gaussian enhancement..."); + " Using tilt-shift Gaussian parameters for Gaussian enhancement..."); } else { NUIS_LOG(FIT, " Using Normal Gaussian parameters for Gaussian enhancement..."); } }; double GaussianModeCorr::GetGausWeight(double q0, double q3, double vals[]) { // The weight double w = 1.0; // Use tilt-shift method by Patrick if (fMethod) { if (fDebugStatements) { std::cout << "Using Patrick gaussian" << std::endl; } // // CCQE Without Suppression // double Norm = 4.82788679036; // double Tilt = 2.3501416116; // double Pq0 = 0.363964889702; // double Wq0 = 0.133976806938; // double Pq3 = 0.431769740224; // double Wq3 = 0.207666663434; // // Also add for CCQE at the end // return (w > 1.0) ? w : 1.0; // // 2p2h with suppression // double Norm = 15.967; // double Tilt = -0.455655; // double Pq0 = 0.214598; // double Wq0 = 0.0291061; // double Pq3 = 0.480194; // double Wq3 = 0.134588; double Norm = vals[kPosNorm]; double Tilt = vals[kPosTilt]; double Pq0 = vals[kPosPq0]; double Wq0 = vals[kPosWq0]; double Pq3 = vals[kPosPq3]; double Wq3 = vals[kPosWq3]; double a = cos(Tilt) * cos(Tilt) / (2 * Wq0 * Wq0); a += sin(Tilt) * sin(Tilt) / (2 * Wq3 * Wq3); double b = -sin(2 * Tilt) / (4 * Wq0 * Wq0); b += sin(2 * Tilt) / (4 * Wq3 * Wq3); double c = sin(Tilt) * sin(Tilt) / (2 * Wq0 * Wq0); c += cos(Tilt) * cos(Tilt) / (2 * Wq3 * Wq3); w = Norm; w *= exp(-a * (q0 - Pq0) * (q0 - Pq0)); w *= exp(+2.0 * b * (q0 - Pq0) * (q3 - Pq3)); w *= exp(-c * (q3 - Pq3) * (q3 - Pq3)); if (fDebugStatements) { std::cout << "Applied Tilt " << Tilt << " " << cos(Tilt) << " " - << sin(Tilt) << std::endl; + << sin(Tilt) << std::endl; std::cout << "abc = " << a << " " << b << " " << c << std::endl; std::cout << "Returning " << Norm << " " << Pq0 << " " << Wq0 << " " - << Pq3 << " " << Wq3 << " " << w << std::endl; + << Pq3 << " " << Wq3 << " " << w << std::endl; } if (w != w || std::isnan(w) || w < 0.0) { w = 0.0; } if (w < 1.0 and !fAllowSuppression) { w = 1.0; } return w; // Use the MINERvA Gaussian method } else { /* * From MINERvA and Daniel Ruterbories: * Old notes here: * * http://cdcvs.fnal.gov/cgi-bin/public-cvs/cvsweb-public.cgi/AnalysisFramework/Ana/CCQENu/ana_common/data/?cvsroot=mnvsoft * These parameters are slightly altered * * FRESH: * 10.5798 * 0.254032 * 0.50834 * 0.0571035 * 0.129051 * 0.875287 */ if (fDebugStatements) { std::cout << "Using MINERvA Gaussian" << std::endl; } double norm = vals[kPosNorm]; double meanq0 = vals[kPosTilt]; double meanq3 = vals[kPosPq0]; double sigmaq0 = vals[kPosWq0]; double sigmaq3 = vals[kPosPq3]; double corr = vals[kPosWq3]; double z = (q0 - meanq0) * (q0 - meanq0) / (sigmaq0 * sigmaq0) + (q3 - meanq3) * (q3 - meanq3) / (sigmaq3 * sigmaq3) - 2 * corr * (q0 - meanq0) * (q3 - meanq3) / (sigmaq0 * sigmaq3); double ret = 1.0; if ( fabs(1 - corr*corr) < 1.E-5 ) { return 1.0; } if ( (-0.5 * z / (1 - corr*corr)) > 200 or (-0.5 * z / (1 - corr*corr)) < -200 ) { return 1.0; } else { ret = norm * exp( -0.5 * z / (1 - corr*corr) ); } if (ret != ret or ret < 0.0 or isnan(ret)) { return 1.0; } if (fAllowSuppression) return ret; return ret + 1.0; // Need to add 1 to the results } return w; } void GaussianModeCorr::SetDialValue(std::string name, double val) { SetDialValue(Reweight::ConvDial(name, kCUSTOM), val); } void GaussianModeCorr::SetDialValue(int rwenum, double val) { int curenum = rwenum % 1000; // Check Handled if (!IsHandled(curenum)) return; // CCQE Setting for (int i = kGaussianCorr_CCQE_norm; i <= kGaussianCorr_CCQE_Wq3; i++) { if (i == curenum) { int index = i - kGaussianCorr_CCQE_norm; fGausVal_CCQE[index] = val; fApply_CCQE = true; } } // 2p2h Setting for (int i = kGaussianCorr_2p2h_norm; i <= kGaussianCorr_2p2h_Wq3; i++) { if (i == curenum) { int index = i - kGaussianCorr_2p2h_norm; fGausVal_2p2h[index] = val; fApply_2p2h = true; } } // 2p2h_PPandNN Setting for (int i = kGaussianCorr_2p2h_PPandNN_norm; i <= kGaussianCorr_2p2h_PPandNN_Wq3; i++) { if (i == curenum) { int index = i - kGaussianCorr_2p2h_PPandNN_norm; fGausVal_2p2h_PPandNN[index] = val; fApply_2p2h_PPandNN = true; } } // 2p2h_NP Setting for (int i = kGaussianCorr_2p2h_NP_norm; i <= kGaussianCorr_2p2h_NP_Wq3; i++) { if (i == curenum) { int index = i - kGaussianCorr_2p2h_NP_norm; fGausVal_2p2h_NP[index] = val; fApply_2p2h_NP = true; } } // CC1pi Setting for (int i = kGaussianCorr_CC1pi_norm; i <= kGaussianCorr_CC1pi_Wq3; i++) { if (i == curenum) { int index = i - kGaussianCorr_CC1pi_norm; fGausVal_CC1pi[index] = val; fApply_CC1pi = true; } } if (curenum == kGaussianCorr_AllowSuppression) { fAllowSuppression = (val > 0.5); } } bool GaussianModeCorr::IsHandled(int rwenum) { int curenum = rwenum % 1000; switch (curenum) { case kGaussianCorr_CCQE_norm: case kGaussianCorr_CCQE_tilt: case kGaussianCorr_CCQE_Pq0: case kGaussianCorr_CCQE_Wq0: case kGaussianCorr_CCQE_Pq3: case kGaussianCorr_CCQE_Wq3: case kGaussianCorr_2p2h_norm: case kGaussianCorr_2p2h_tilt: case kGaussianCorr_2p2h_Pq0: case kGaussianCorr_2p2h_Wq0: case kGaussianCorr_2p2h_Pq3: case kGaussianCorr_2p2h_Wq3: case kGaussianCorr_2p2h_PPandNN_norm: case kGaussianCorr_2p2h_PPandNN_tilt: case kGaussianCorr_2p2h_PPandNN_Pq0: case kGaussianCorr_2p2h_PPandNN_Wq0: case kGaussianCorr_2p2h_PPandNN_Pq3: case kGaussianCorr_2p2h_PPandNN_Wq3: case kGaussianCorr_2p2h_NP_norm: case kGaussianCorr_2p2h_NP_tilt: case kGaussianCorr_2p2h_NP_Pq0: case kGaussianCorr_2p2h_NP_Wq0: case kGaussianCorr_2p2h_NP_Pq3: case kGaussianCorr_2p2h_NP_Wq3: case kGaussianCorr_CC1pi_norm: case kGaussianCorr_CC1pi_tilt: case kGaussianCorr_CC1pi_Pq0: case kGaussianCorr_CC1pi_Wq0: case kGaussianCorr_CC1pi_Pq3: case kGaussianCorr_CC1pi_Wq3: case kGaussianCorr_AllowSuppression: return true; default: return false; } } diff --git a/src/Reweight/NUISANCEWeightCalcs.h b/src/Reweight/NUISANCEWeightCalcs.h index a1d2b9e..520acfa 100644 --- a/src/Reweight/NUISANCEWeightCalcs.h +++ b/src/Reweight/NUISANCEWeightCalcs.h @@ -1,129 +1,179 @@ #ifndef NUISANCE_WEIGHT_CALCS #define NUISANCE_WEIGHT_CALCS #include "BaseFitEvt.h" #include "BeRPA.h" class NUISANCEWeightCalc { public: NUISANCEWeightCalc() {}; virtual ~NUISANCEWeightCalc() {}; virtual double CalcWeight(BaseFitEvt* evt){return 1.0;}; virtual void SetDialValue(std::string name, double val){}; virtual void SetDialValue(int rwenum, double val){}; virtual bool IsHandled(int rwenum){return false;}; virtual void Print(){}; std::map fDialNameIndex; std::map fDialEnumIndex; std::vector fDialValues; std::string fName; }; class ModeNormCalc : public NUISANCEWeightCalc { public: ModeNormCalc(); ~ModeNormCalc(){}; double CalcWeight(BaseFitEvt* evt); void SetDialValue(std::string name, double val); void SetDialValue(int rwenum, double val); bool IsHandled(int rwenum); double fNormRES; }; +// MINOS pion tuning, https://arxiv.org/pdf/1903.01558.pdf +class MINOSRPA : public NUISANCEWeightCalc { + public: + MINOSRPA(); + ~MINOSRPA(){}; + + double CalcWeight(BaseFitEvt* evt); + void SetDialValue(std::string name, double val); + void SetDialValue(int rwenum, double val); + bool IsHandled(int rwenum); + + double GetRPAWeight(double Q2); + + bool fTweaked; + + bool fApply_MINOSRPA; + + double fTwk_MINOSRPA_A; + double fDef_MINOSRPA_A; + double fCur_MINOSRPA_A; + + double fTwk_MINOSRPA_B; + double fDef_MINOSRPA_B; + double fCur_MINOSRPA_B; + +}; + +// MINERvA pion tuning, https://arxiv.org/pdf/1903.01558.pdf +class LagrangeRPA : public NUISANCEWeightCalc { + public: + LagrangeRPA(); + ~LagrangeRPA(){}; + + double CalcWeight(BaseFitEvt* evt); + void SetDialValue(std::string name, double val); + void SetDialValue(int rwenum, double val); + bool IsHandled(int rwenum); + + double GetRPAWeight(double Q2); + + bool fTweaked; + + bool fApplyRPA; + + double fR1; + double fR2; + double fR1_Def; + double fR2_Def; +}; + class BeRPACalc : public NUISANCEWeightCalc { public: BeRPACalc(); ~BeRPACalc(){}; double CalcWeight(BaseFitEvt* evt); void SetDialValue(std::string name, double val); void SetDialValue(int rwenum, double val); bool IsHandled(int rwenum); private: // Parameter values double fBeRPA_A; double fBeRPA_B; double fBeRPA_D; double fBeRPA_E; double fBeRPA_U; // Counts of enabled parameters int nParams; }; class SBLOscWeightCalc : public NUISANCEWeightCalc { public: SBLOscWeightCalc(); ~SBLOscWeightCalc(){}; double CalcWeight(BaseFitEvt* evt); void SetDialValue(std::string name, double val); void SetDialValue(int rwenum, double val); bool IsHandled(int rwenum); double GetSBLOscWeight(double E); double fDistance; double fMassSplitting; double fSin2Theta; }; class GaussianModeCorr : public NUISANCEWeightCalc { public: GaussianModeCorr(); ~GaussianModeCorr(){}; double CalcWeight(BaseFitEvt* evt); void SetDialValue(std::string name, double val); void SetDialValue(int rwenum, double val); bool IsHandled(int rwenum); double GetGausWeight(double q0, double q3, double vals[]); // Set the Gaussian method (tilt-shift or normal Gaussian parameters) void SetMethod(bool method); // 5 pars describe the Gaussain // 0 norm. // 1 q0 pos // 2 q0 width // 3 q3 pos // 4 q3 width static const int kPosNorm = 0; static const int kPosTilt = 1; static const int kPosPq0 = 2; static const int kPosWq0 = 3; static const int kPosPq3 = 4; static const int kPosWq3 = 5; bool fApply_CCQE; double fGausVal_CCQE[6]; bool fApply_2p2h; double fGausVal_2p2h[6]; bool fApply_2p2h_PPandNN; double fGausVal_2p2h_PPandNN[6]; bool fApply_2p2h_NP; double fGausVal_2p2h_NP[6]; bool fApply_CC1pi; double fGausVal_CC1pi[6]; bool fAllowSuppression; bool fDebugStatements; bool fMethod; }; #endif diff --git a/src/Reweight/NUISANCEWeightEngine.cxx b/src/Reweight/NUISANCEWeightEngine.cxx index 033ed0c..f74ccd7 100644 --- a/src/Reweight/NUISANCEWeightEngine.cxx +++ b/src/Reweight/NUISANCEWeightEngine.cxx @@ -1,137 +1,139 @@ #include "NUISANCEWeightEngine.h" #include "NUISANCEWeightCalcs.h" #ifdef __MINERVA_RW_ENABLED__ #ifdef __GENIE_ENABLED__ #include "MINERvAWeightCalcs.h" #endif #endif NUISANCEWeightEngine::NUISANCEWeightEngine(std::string name) { // Setup the NUISANCE Reweight engine fCalcName = name; NUIS_LOG(FIT, "Setting up NUISANCE Custom RW : " << fCalcName); // Load in all Weight Calculations GaussianModeCorr *GaussianMode = new GaussianModeCorr(); std::string Gaussian_Method = FitPar::Config().GetParS("Gaussian_Enhancement"); if (Gaussian_Method == "Tilt-Shift") { GaussianMode->SetMethod(true); } else if (Gaussian_Method == "Normal") { GaussianMode->SetMethod(false); } else { NUIS_ABORT("I do not recognise method " << Gaussian_Method << " for the Gaussian enhancement, so will die now..."); } + // The NUISANCE calculators fWeightCalculators.push_back(GaussianMode); fWeightCalculators.push_back(new ModeNormCalc()); fWeightCalculators.push_back(new SBLOscWeightCalc()); fWeightCalculators.push_back(new BeRPACalc()); + fWeightCalculators.push_back(new MINOSRPA()); + fWeightCalculators.push_back(new LagrangeRPA()); + // The MINERvA calculators that rely on GENIE variables (so need GENIE support) #ifdef __MINERVA_RW_ENABLED__ #ifdef __GENIE_ENABLED__ fWeightCalculators.push_back(new nuisance::reweight::MINERvAReWeight_QE()); fWeightCalculators.push_back(new nuisance::reweight::MINERvAReWeight_MEC()); fWeightCalculators.push_back(new nuisance::reweight::MINERvAReWeight_RES()); - fWeightCalculators.push_back(new nuisance::reweight::MINOSRPA()); - fWeightCalculators.push_back(new nuisance::reweight::LagrangeRPA()); fWeightCalculators.push_back(new nuisance::reweight::RikRPA()); fWeightCalculators.push_back(new nuisance::reweight::COHBrandon()); fWeightCalculators.push_back(new nuisance::reweight::WEnhancement()); #endif #endif // Set Abs Twk Config fIsAbsTwk = true; }; void NUISANCEWeightEngine::IncludeDial(std::string name, double startval) { // Get First enum int nuisenum = Reweight::ConvDial(name, kCUSTOM); // Setup Maps fEnumIndex[nuisenum]; // = std::vector(0); fNameIndex[name]; // = std::vector(0); // Split by commas std::vector allnames = GeneralUtils::ParseToStr(name, ","); for (uint i = 0; i < allnames.size(); i++) { std::string singlename = allnames[i]; // Get RW int singleenum = Reweight::ConvDial(singlename, kCUSTOM); // Fill Maps int index = fValues.size(); fValues.push_back(0.0); fNUISANCEEnums.push_back(singleenum); // Initialize dial NUIS_LOG(FIT, "Registering " << singlename << " from " << name); // Setup index fEnumIndex[nuisenum].push_back(index); fNameIndex[name].push_back(index); } // Set Value if given if (startval != _UNDEF_DIAL_VALUE_) { SetDialValue(nuisenum, startval); } }; void NUISANCEWeightEngine::SetDialValue(int nuisenum, double val) { std::vector indices = fEnumIndex[nuisenum]; for (uint i = 0; i < indices.size(); i++) { fValues[indices[i]] = val; } } void NUISANCEWeightEngine::SetDialValue(std::string name, double val) { std::vector indices = fNameIndex[name]; for (uint i = 0; i < indices.size(); i++) { fValues[indices[i]] = val; } } void NUISANCEWeightEngine::Reconfigure(bool silent) { for (size_t i = 0; i < fNUISANCEEnums.size(); i++) { // Is this parameter handled bool IsHandledSomewhere = false; for (std::vector::iterator calciter = fWeightCalculators.begin(); calciter != fWeightCalculators.end(); calciter++) { NUISANCEWeightCalc *nuiscalc = static_cast(*calciter); if (nuiscalc->IsHandled(fNUISANCEEnums[i])) { nuiscalc->SetDialValue(fNUISANCEEnums[i], fValues[i]); IsHandledSomewhere = true; } } // Check if this parameter is actually handled if (!IsHandledSomewhere) { std::string name = GetNameFromEnum(i); NUIS_ABORT("NUISANCE parameter " << name << " (enum " << i << ") not enabled, can not continue!"); } // End check of it's handled } // End loop over enums } // Return double NUISANCEWeightEngine::CalcWeight(BaseFitEvt *evt) { double rw_weight = 1.0; // Cast as usable class for (std::vector::iterator iter = fWeightCalculators.begin(); iter != fWeightCalculators.end(); iter++) { NUISANCEWeightCalc *nuiscalc = static_cast(*iter); rw_weight *= nuiscalc->CalcWeight(evt); } // Return rw_weight return rw_weight; }