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;
}