Page MenuHomeHEPForge
Files F9501507

No OneTemporary
Actions

View Options

diff --git a/src/InputHandler/FitEvent.h b/src/InputHandler/FitEvent.h
index 55a15bb..b467835 100644
--- a/src/InputHandler/FitEvent.h
+++ b/src/InputHandler/FitEvent.h
@@ -1,647 +1,653 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************/
#ifndef FITEVENT2_H_SEEN
#define FITEVENT2_H_SEEN
/*!
* \addtogroup InputHandler
* @{
*/
#include <algorithm>
#include <iterator>
#include <vector>
#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); };
/// 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 <size_t N>
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 <size_t N>
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<int> GetAllParticleIndices (int const pdg = -1, int const state = -1) const;
template <size_t N>
inline std::vector<int> GetAllParticleIndices(int const (&pdgs)[N], const int state = -1) const {
std::vector<int> plist;
for (size_t i = 0; i < N; i++) {
std::vector<int> 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<FitParticle*> GetAllParticle (int const pdg = -1, int const state = -1);
template <size_t N>
inline std::vector<FitParticle*> GetAllParticle(int const (&pdgs)[N], int const state = -1) {
std::vector<FitParticle*> plist;
for (size_t i = 0; i < N; i++) {
std::vector<FitParticle*> plisttemp = GetAllParticle(pdgs[i], state);
plist.insert( plist.end(), plisttemp.begin(), plisttemp.end() );
}
return plist;
}
inline std::vector<int> GetAllNuElectronIndices (void) { return GetAllParticleIndices(12); };
inline std::vector<int> GetAllNuMuonIndices (void) { return GetAllParticleIndices(14); };
inline std::vector<int> GetAllNuTauIndices (void) { return GetAllParticleIndices(16); };
inline std::vector<int> GetAllElectronIndices (void) { return GetAllParticleIndices(11); };
inline std::vector<int> GetAllMuonIndices (void) { return GetAllParticleIndices(13); };
inline std::vector<int> GetAllTauIndices (void) { return GetAllParticleIndices(15); };
inline std::vector<int> GetAllProtonIndices (void) { return GetAllParticleIndices(2212); };
inline std::vector<int> GetAllNeutronIndices (void) { return GetAllParticleIndices(2112); };
inline std::vector<int> GetAllPiZeroIndices (void) { return GetAllParticleIndices(111); };
inline std::vector<int> GetAllPiPlusIndices (void) { return GetAllParticleIndices(211); };
inline std::vector<int> GetAllPiMinusIndices (void) { return GetAllParticleIndices(-211); };
inline std::vector<int> GetAllPhotonIndices (void) { return GetAllParticleIndices(22); };
inline std::vector<FitParticle*> GetAllNuElectron (void) { return GetAllParticle(12); };
inline std::vector<FitParticle*> GetAllNuMuon (void) { return GetAllParticle(14); };
inline std::vector<FitParticle*> GetAllNuTau (void) { return GetAllParticle(16); };
inline std::vector<FitParticle*> GetAllElectron (void) { return GetAllParticle(11); };
inline std::vector<FitParticle*> GetAllMuon (void) { return GetAllParticle(13); };
inline std::vector<FitParticle*> GetAllTau (void) { return GetAllParticle(15); };
inline std::vector<FitParticle*> GetAllProton (void) { return GetAllParticle(2212); };
inline std::vector<FitParticle*> GetAllNeutron (void) { return GetAllParticle(2112); };
inline std::vector<FitParticle*> GetAllPiZero (void) { return GetAllParticle(111); };
inline std::vector<FitParticle*> GetAllPiPlus (void) { return GetAllParticle(211); };
inline std::vector<FitParticle*> GetAllPiMinus (void) { return GetAllParticle(-211); };
inline std::vector<FitParticle*> 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 <size_t N>
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 <size_t N>
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 <size_t N>
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 <size_t N>
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<int> GetAllISParticleIndices(int const pdg = -1) const {
return GetAllParticleIndices(pdg, kInitialState);
};
template <size_t N>
inline std::vector<int> 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<FitParticle*> GetAllISParticle(int const pdg = -1) {
return GetAllParticle(pdg, kInitialState);
};
template <size_t N>
inline std::vector<FitParticle*> 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 <size_t N>
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 <size_t N>
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<int> GetAllISNuElectronIndices (void) const { return GetAllISParticleIndices(12); };
inline std::vector<int> GetAllISNuMuonIndices (void) const { return GetAllISParticleIndices(14); };
inline std::vector<int> GetAllISNuTauIndices (void) const { return GetAllISParticleIndices(16); };
inline std::vector<int> GetAllISElectronIndices (void) const { return GetAllISParticleIndices(11); };
inline std::vector<int> GetAllISMuonIndices (void) const { return GetAllISParticleIndices(13); };
inline std::vector<int> GetAllISTauIndices (void) const { return GetAllISParticleIndices(15); };
inline std::vector<int> GetAllISProtonIndices (void) const { return GetAllISParticleIndices(2212); };
inline std::vector<int> GetAllISNeutronIndices (void) const { return GetAllISParticleIndices(2112); };
inline std::vector<int> GetAllISPiZeroIndices (void) const { return GetAllISParticleIndices(111); };
inline std::vector<int> GetAllISPiPlusIndices (void) const { return GetAllISParticleIndices(211); };
inline std::vector<int> GetAllISPiMinusIndices (void) const { return GetAllISParticleIndices(-211); };
inline std::vector<int> GetAllISPhotonIndices (void) const { return GetAllISParticleIndices(22); };
inline std::vector<int> GetAllISLeptonsIndices (void) const { return GetAllISParticleIndices(PhysConst::pdg_leptons); };
inline std::vector<int> GetAllISPionsIndices (void) const { return GetAllISParticleIndices(PhysConst::pdg_pions); };
inline std::vector<int> GetAllISChargePionsIndices(void) const { return GetAllISParticleIndices(PhysConst::pdg_charged_pions); };
inline std::vector<FitParticle*> GetAllISNuElectron (void) { return GetAllISParticle(12); };
inline std::vector<FitParticle*> GetAllISNuMuon (void) { return GetAllISParticle(14); };
inline std::vector<FitParticle*> GetAllISNuTau (void) { return GetAllISParticle(16); };
inline std::vector<FitParticle*> GetAllISElectron (void) { return GetAllISParticle(11); };
inline std::vector<FitParticle*> GetAllISMuon (void) { return GetAllISParticle(13); };
inline std::vector<FitParticle*> GetAllISTau (void) { return GetAllISParticle(15); };
inline std::vector<FitParticle*> GetAllISProton (void) { return GetAllISParticle(2212); };
inline std::vector<FitParticle*> GetAllISNeutron (void) { return GetAllISParticle(2112); };
inline std::vector<FitParticle*> GetAllISPiZero (void) { return GetAllISParticle(111); };
inline std::vector<FitParticle*> GetAllISPiPlus (void) { return GetAllISParticle(211); };
inline std::vector<FitParticle*> GetAllISPiMinus (void) { return GetAllISParticle(-211); };
inline std::vector<FitParticle*> GetAllISPhoton (void) { return GetAllISParticle(22); };
inline std::vector<FitParticle*> GetAllISLeptons (void) { return GetAllISParticle(PhysConst::pdg_leptons); };
inline std::vector<FitParticle*> GetAllISPions (void) { return GetAllISParticle(PhysConst::pdg_pions); };
inline std::vector<FitParticle*> 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 <size_t N>
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 <size_t N>
inline int NumFSParticle(int const (&pdgs)[N]) const {
return NumParticle(pdgs, kFinalState);
};
inline std::vector<int> GetAllFSParticleIndices(int const pdg = -1) const {
return GetAllParticleIndices(pdg, kFinalState);
};
template <size_t N>
inline std::vector<int> GetAllFSParticleIndices(int const (&pdgs)[N]) const {
return GetAllParticleIndices(pdgs, kFinalState);
};
inline std::vector<FitParticle*> GetAllFSParticle(int const pdg = -1) {
return GetAllParticle(pdg, kFinalState);
};
template <size_t N>
inline std::vector<FitParticle*> GetAllFSParticle(int const (&pdgs)[N]) {
return GetAllParticle(pdgs, kFinalState);
};
inline FitParticle* GetHMFSParticle(int const pdg) {
return GetHMParticle(pdg, kFinalState);
};
template <size_t N>
inline FitParticle* GetHMFSParticle(int const (&pdgs)[N]) {
return GetHMParticle(pdgs, kFinalState);
};
inline int GetHMFSParticleIndex(int const pdg) const {
return GetHMParticleIndex(pdg, kFinalState);
};
template <size_t N>
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<int> GetAllFSNuElectronIndices (void) const { return GetAllFSParticleIndices(12); };
inline std::vector<int> GetAllFSNuMuonIndices (void) const { return GetAllFSParticleIndices(14); };
inline std::vector<int> GetAllFSNuTauIndices (void) const { return GetAllFSParticleIndices(16); };
inline std::vector<int> GetAllFSElectronIndices (void) const { return GetAllFSParticleIndices(11); };
inline std::vector<int> GetAllFSMuonIndices (void) const { return GetAllFSParticleIndices(13); };
inline std::vector<int> GetAllFSTauIndices (void) const { return GetAllFSParticleIndices(15); };
inline std::vector<int> GetAllFSProtonIndices (void) const { return GetAllFSParticleIndices(2212); };
inline std::vector<int> GetAllFSNeutronIndices (void) const { return GetAllFSParticleIndices(2112); };
inline std::vector<int> GetAllFSPiZeroIndices (void) const { return GetAllFSParticleIndices(111); };
inline std::vector<int> GetAllFSPiPlusIndices (void) const { return GetAllFSParticleIndices(211); };
inline std::vector<int> GetAllFSPiMinusIndices (void) const { return GetAllFSParticleIndices(-211); };
inline std::vector<int> GetAllFSPhotonIndices (void) const { return GetAllFSParticleIndices(22); };
inline std::vector<int> GetAllFSLeptonsIndices (void) const { return GetAllFSParticleIndices(PhysConst::pdg_leptons); };
inline std::vector<int> GetAllFSPionsIndices (void) const { return GetAllFSParticleIndices(PhysConst::pdg_pions); };
inline std::vector<int> GetAllFSChargePionsIndices(void) const { return GetAllFSParticleIndices(PhysConst::pdg_charged_pions); };
+ inline std::vector<int> GetAllFSNucleonIndices(void) const { return GetAllFSParticleIndices(PhysConst::pdg_nucleons); };
inline std::vector<FitParticle*> GetAllFSNuElectron (void) { return GetAllFSParticle(12); };
inline std::vector<FitParticle*> GetAllFSNuMuon (void) { return GetAllFSParticle(14); };
inline std::vector<FitParticle*> GetAllFSNuTau (void) { return GetAllFSParticle(16); };
inline std::vector<FitParticle*> GetAllFSElectron (void) { return GetAllFSParticle(11); };
inline std::vector<FitParticle*> GetAllFSMuon (void) { return GetAllFSParticle(13); };
inline std::vector<FitParticle*> GetAllFSTau (void) { return GetAllFSParticle(15); };
inline std::vector<FitParticle*> GetAllFSProton (void) { return GetAllFSParticle(2212); };
inline std::vector<FitParticle*> GetAllFSNeutron (void) { return GetAllFSParticle(2112); };
inline std::vector<FitParticle*> GetAllFSPiZero (void) { return GetAllFSParticle(111); };
inline std::vector<FitParticle*> GetAllFSPiPlus (void) { return GetAllFSParticle(211); };
inline std::vector<FitParticle*> GetAllFSPiMinus (void) { return GetAllFSParticle(-211); };
inline std::vector<FitParticle*> GetAllFSPhoton (void) { return GetAllFSParticle(22); };
inline std::vector<FitParticle*> GetAllFSLeptons (void) { return GetAllFSParticle(PhysConst::pdg_leptons); };
inline std::vector<FitParticle*> GetAllFSPions (void) { return GetAllFSParticle(PhysConst::pdg_pions); };
inline std::vector<FitParticle*> GetAllFSChargePions (void) { return GetAllFSParticle(PhysConst::pdg_charged_pions); };
+ inline std::vector<FitParticle*> 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);
// 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/MCStudies/GenericFlux_Vectors.cxx b/src/MCStudies/GenericFlux_Vectors.cxx
index 5f05293..69e07ed 100644
--- a/src/MCStudies/GenericFlux_Vectors.cxx
+++ b/src/MCStudies/GenericFlux_Vectors.cxx
@@ -1,445 +1,458 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "GenericFlux_Vectors.h"
GenericFlux_Vectors::GenericFlux_Vectors(std::string name,
std::string inputfile, FitWeight *rw,
std::string type,
std::string fakeDataFile) {
// Measurement Details
fName = name;
eventVariables = NULL;
// Define our energy range for flux calcs
EnuMin = 0.;
EnuMax = 1E10; // Arbritrarily high energy limit
if (Config::HasPar("EnuMin")) {
EnuMin = Config::GetParD("EnuMin");
}
if (Config::HasPar("EnuMax")) {
EnuMax = Config::GetParD("EnuMax");
}
SavePreFSI = Config::Get().GetParB("nuisflat_SavePreFSI");
NUIS_LOG(SAM,
"Running GenericFlux_Vectors saving pre-FSI particles? " << SavePreFSI);
// Set default fitter flags
fIsDiag = true;
fIsShape = false;
fIsRawEvents = false;
// This function will sort out the input files automatically and parse all the
// inputs,flags,etc.
// There may be complex cases where you have to do this by hand, but usually
// this will do.
Measurement1D::SetupMeasurement(inputfile, type, rw, fakeDataFile);
eventVariables = NULL;
// Setup fDataHist as a placeholder
this->fDataHist = new TH1D(("empty_data"), ("empty-data"), 1, 0, 1);
this->SetupDefaultHist();
fFullCovar = StatUtils::MakeDiagonalCovarMatrix(fDataHist);
covar = StatUtils::GetInvert(fFullCovar);
// 1. The generator is organised in SetupMeasurement so it gives the
// cross-section in "per nucleon" units.
// So some extra scaling for a specific measurement may be required. For
// Example to get a "per neutron" measurement on carbon
// which we do here, we have to multiple by the number of nucleons 12 and
// divide by the number of neutrons 6.
// N.B. MeasurementBase::PredictedEventRate includes the 1E-38 factor that is
// often included here in other classes that directly integrate the event
// histogram. This method is used here as it now respects EnuMin and EnuMax
// correctly.
this->fScaleFactor =
(this->PredictedEventRate("width", 0, EnuMax) / double(fNEvents)) /
this->TotalIntegratedFlux("width");
NUIS_LOG(SAM, " Generic Flux Scaling Factor = "
<< fScaleFactor
<< " [= " << (GetEventHistogram()->Integral("width") * 1E-38)
<< "/(" << (fNEvents + 0.) << "*"
<< TotalIntegratedFlux("width") << ")]");
if (fScaleFactor <= 0.0) {
NUIS_ABORT("SCALE FACTOR TOO LOW");
}
// Setup our TTrees
this->AddEventVariablesToTree();
this->AddSignalFlagsToTree();
}
void GenericFlux_Vectors::AddEventVariablesToTree() {
// Setup the TTree to save everything
if (!eventVariables) {
Config::Get().out->cd();
eventVariables = new TTree((this->fName + "_VARS").c_str(),
(this->fName + "_VARS").c_str());
}
NUIS_LOG(SAM, "Adding Event Variables");
eventVariables->Branch("Mode", &Mode, "Mode/I");
eventVariables->Branch("cc", &cc, "cc/B");
eventVariables->Branch("PDGnu", &PDGnu, "PDGnu/I");
eventVariables->Branch("Enu_true", &Enu_true, "Enu_true/F");
eventVariables->Branch("tgt", &tgt, "tgt/I");
eventVariables->Branch("tgta", &tgta, "tgta/I");
eventVariables->Branch("tgtz", &tgtz, "tgtz/I");
eventVariables->Branch("PDGLep", &PDGLep, "PDGLep/I");
eventVariables->Branch("ELep", &ELep, "ELep/F");
eventVariables->Branch("CosLep", &CosLep, "CosLep/F");
// Basic interaction kinematics
eventVariables->Branch("Q2", &Q2, "Q2/F");
eventVariables->Branch("q0", &q0, "q0/F");
eventVariables->Branch("q3", &q3, "q3/F");
eventVariables->Branch("Enu_QE", &Enu_QE, "Enu_QE/F");
eventVariables->Branch("Q2_QE", &Q2_QE, "Q2_QE/F");
eventVariables->Branch("W_nuc_rest", &W_nuc_rest, "W_nuc_rest/F");
eventVariables->Branch("W", &W, "W/F");
eventVariables->Branch("W_genie", &W_genie, "W_genie/F");
eventVariables->Branch("x", &x, "x/F");
eventVariables->Branch("y", &y, "y/F");
eventVariables->Branch("Eav", &Eav, "Eav/F");
eventVariables->Branch("EavAlt", &EavAlt, "EavAlt/F");
+ eventVariables->Branch("CosThetaAdler", &CosThetaAdler, "CosThetaAdler/F");
+ eventVariables->Branch("PhiAdler", &PhiAdler, "PhiAdler/F");
+
eventVariables->Branch("dalphat", &dalphat, "dalphat/F");
eventVariables->Branch("dpt", &dpt, "dpt/F");
eventVariables->Branch("dphit", &dphit, "dphit/F");
eventVariables->Branch("pnreco_C", &pnreco_C, "pnreco_C/F");
// Save outgoing particle vectors
eventVariables->Branch("nfsp", &nfsp, "nfsp/I");
eventVariables->Branch("px", px, "px[nfsp]/F");
eventVariables->Branch("py", py, "py[nfsp]/F");
eventVariables->Branch("pz", pz, "pz[nfsp]/F");
eventVariables->Branch("E", E, "E[nfsp]/F");
eventVariables->Branch("pdg", pdg, "pdg[nfsp]/I");
eventVariables->Branch("pdg_rank", pdg_rank, "pdg_rank[nfsp]/I");
// Save init particle vectors
eventVariables->Branch("ninitp", &ninitp, "ninitp/I");
eventVariables->Branch("px_init", px_init, "px_init[ninitp]/F");
eventVariables->Branch("py_init", py_init, "py_init[ninitp]/F");
eventVariables->Branch("pz_init", pz_init, "pz_init[ninitp]/F");
eventVariables->Branch("E_init", E_init, "E_init[ninitp]/F");
eventVariables->Branch("pdg_init", pdg_init, "pdg_init[ninitp]/I");
// Save pre-FSI vectors
eventVariables->Branch("nvertp", &nvertp, "nvertp/I");
eventVariables->Branch("px_vert", px_vert, "px_vert[nvertp]/F");
eventVariables->Branch("py_vert", py_vert, "py_vert[nvertp]/F");
eventVariables->Branch("pz_vert", pz_vert, "pz_vert[nvertp]/F");
eventVariables->Branch("E_vert", E_vert, "E_vert[nvertp]/F");
eventVariables->Branch("pdg_vert", pdg_vert, "pdg_vert[nvertp]/I");
// Event Scaling Information
eventVariables->Branch("Weight", &Weight, "Weight/F");
eventVariables->Branch("InputWeight", &InputWeight, "InputWeight/F");
eventVariables->Branch("RWWeight", &RWWeight, "RWWeight/F");
// Should be a double because may be 1E-39 and less
eventVariables->Branch("fScaleFactor", &fScaleFactor, "fScaleFactor/D");
// The customs
eventVariables->Branch("CustomWeight", &CustomWeight, "CustomWeight/F");
eventVariables->Branch("CustomWeightArray", CustomWeightArray,
"CustomWeightArray[6]/F");
return;
}
void GenericFlux_Vectors::FillEventVariables(FitEvent *event) {
ResetVariables();
// Fill Signal Variables
FillSignalFlags(event);
NUIS_LOG(DEB, "Filling signal");
// Now fill the information
Mode = event->Mode;
cc = event->IsCC();
// Get the incoming neutrino and outgoing lepton
FitParticle *nu = event->GetBeamPart();
FitParticle *lep = event->GetHMFSAnyLepton();
PDGnu = nu->fPID;
Enu_true = nu->fP.E() / 1E3;
tgt = event->fTargetPDG;
tgta = event->fTargetA;
tgtz = event->fTargetZ;
if (lep != NULL) {
PDGLep = lep->fPID;
ELep = lep->fP.E() / 1E3;
CosLep = cos(nu->fP.Vect().Angle(lep->fP.Vect()));
// Basic interaction kinematics
Q2 = -1 * (nu->fP - lep->fP).Mag2() / 1E6;
q0 = (nu->fP - lep->fP).E() / 1E3;
q3 = (nu->fP - lep->fP).Vect().Mag() / 1E3;
// These assume C12 binding from MINERvA... not ideal
Enu_QE = FitUtils::EnuQErec(lep->fP, CosLep, 34., true);
Q2_QE = FitUtils::Q2QErec(lep->fP, CosLep, 34., true);
Eav = FitUtils::GetErecoil_MINERvA_LowRecoil(event) / 1.E3;
EavAlt = FitUtils::Eavailable(event) / 1.E3;
+ // Check if this is a 1pi+ or 1pi0 event
+ if ((SignalDef::isCC1pi(event, PDGnu, 211) || SignalDef::isCC1pi(event, PDGnu, -211) || SignalDef::isCC1pi(event, PDGnu, 111)) && event->NumFSNucleons() == 1) {
+ TLorentzVector Pnu = nu->fP;
+ TLorentzVector Pmu = lep->fP;
+ TLorentzVector Ppi = event->GetHMFSPions()->fP;
+ TLorentzVector Pprot = event->GetHMFSNucleons()->fP;
+ CosThetaAdler = FitUtils::CosThAdler(Pnu, Pmu, Ppi, Pprot);
+ PhiAdler = FitUtils::PhiAdler(Pnu, Pmu, Ppi, Pprot);
+ }
+
// Get W_true with assumption of initial state nucleon at rest
float m_n = (float)PhysConst::mass_proton;
// Q2 assuming nucleon at rest
W_nuc_rest = sqrt(-Q2 + 2 * m_n * q0 + m_n * m_n);
// True Q2
W = sqrt(-Q2 + 2 * m_n * q0 + m_n * m_n);
x = Q2 / (2 * m_n * q0);
y = 1 - ELep / Enu_true;
dalphat = FitUtils::Get_STV_dalphat(event, PDGnu, true);
dpt = FitUtils::Get_STV_dpt(event, PDGnu, true);
dphit = FitUtils::Get_STV_dphit(event, PDGnu, true);
pnreco_C = FitUtils::Get_pn_reco_C(event, PDGnu, true);
}
// Loop over the particles and store all the final state particles in a vector
for (UInt_t i = 0; i < event->Npart(); ++i) {
if (event->PartInfo(i)->fIsAlive &&
event->PartInfo(i)->Status() == kFinalState)
partList.push_back(event->PartInfo(i));
if (SavePreFSI && event->fPrimaryVertex[i])
vertList.push_back(event->PartInfo(i));
if (SavePreFSI && event->PartInfo(i)->IsInitialState())
initList.push_back(event->PartInfo(i));
}
// Save outgoing particle vectors
nfsp = (int)partList.size();
std::map<int, std::vector<std::pair<double, int> > > pdgMap;
for (int i = 0; i < nfsp; ++i) {
px[i] = partList[i]->fP.X() / 1E3;
py[i] = partList[i]->fP.Y() / 1E3;
pz[i] = partList[i]->fP.Z() / 1E3;
E[i] = partList[i]->fP.E() / 1E3;
pdg[i] = partList[i]->fPID;
pdgMap[pdg[i]].push_back(std::make_pair(partList[i]->fP.Vect().Mag(), i));
}
for (std::map<int, std::vector<std::pair<double, int> > >::iterator iter =
pdgMap.begin();
iter != pdgMap.end(); ++iter) {
std::vector<std::pair<double, int> > thisVect = iter->second;
std::sort(thisVect.begin(), thisVect.end());
// Now save the order... a bit funky to avoid inverting
int nPart = (int)thisVect.size() - 1;
for (int i = nPart; i >= 0; --i) {
pdg_rank[thisVect[i].second] = nPart - i;
}
}
// Save pre-FSI particles
nvertp = (int)vertList.size();
for (int i = 0; i < nvertp; ++i) {
px_vert[i] = vertList[i]->fP.X() / 1E3;
py_vert[i] = vertList[i]->fP.Y() / 1E3;
pz_vert[i] = vertList[i]->fP.Z() / 1E3;
E_vert[i] = vertList[i]->fP.E() / 1E3;
pdg_vert[i] = vertList[i]->fPID;
}
// Save init particles
ninitp = (int)initList.size();
for (int i = 0; i < ninitp; ++i) {
px_init[i] = initList[i]->fP.X() / 1E3;
py_init[i] = initList[i]->fP.Y() / 1E3;
pz_init[i] = initList[i]->fP.Z() / 1E3;
E_init[i] = initList[i]->fP.E() / 1E3;
pdg_init[i] = initList[i]->fPID;
}
#ifdef __GENIE_ENABLED__
if (event->fType == kGENIE) {
EventRecord *gevent = static_cast<EventRecord *>(event->genie_event->event);
const Interaction *interaction = gevent->Summary();
const Kinematics &kine = interaction->Kine();
double W_genie = kine.W();
}
#endif
// Fill event weights
Weight = event->RWWeight * event->InputWeight;
RWWeight = event->RWWeight;
InputWeight = event->InputWeight;
// And the Customs
CustomWeight = event->CustomWeight;
for (int i = 0; i < 6; ++i) {
CustomWeightArray[i] = event->CustomWeightArray[i];
}
// Fill the eventVariables Tree
eventVariables->Fill();
return;
};
//********************************************************************
void GenericFlux_Vectors::ResetVariables() {
//********************************************************************
cc = false;
// Reset all Function used to extract any variables of interest to the event
Mode = PDGnu = tgt = tgta = tgtz = PDGLep = 0;
Enu_true = ELep = CosLep = Q2 = q0 = q3 = Enu_QE = Q2_QE = W_nuc_rest = W =
- x = y = Eav = EavAlt = -999.9;
+ x = y = Eav = EavAlt = CosThetaAdler = PhiAdler = -999.9;
W_genie = -999;
// Other fun variables
// MINERvA-like ones
dalphat = dpt = dphit = pnreco_C = -999.99;
nfsp = ninitp = nvertp = 0;
for (int i = 0; i < kMAX; ++i) {
px[i] = py[i] = pz[i] = E[i] = -999;
pdg[i] = pdg_rank[i] = 0;
px_init[i] = py_init[i] = pz_init[i] = E_init[i] = -999;
pdg_init[i] = 0;
px_vert[i] = py_vert[i] = pz_vert[i] = E_vert[i] = -999;
pdg_vert[i] = 0;
}
Weight = InputWeight = RWWeight = 0.0;
CustomWeight = 0.0;
for (int i = 0; i < 6; ++i)
CustomWeightArray[i] = 0.0;
partList.clear();
initList.clear();
vertList.clear();
flagCCINC = flagNCINC = flagCCQE = flagCC0pi = flagCCQELike = flagNCEL =
flagNC0pi = flagCCcoh = flagNCcoh = flagCC1pip = flagNC1pip = flagCC1pim =
flagNC1pim = flagCC1pi0 = flagNC1pi0 = flagCC0piMINERvA = false;
}
//********************************************************************
void GenericFlux_Vectors::FillSignalFlags(FitEvent *event) {
//********************************************************************
// Some example flags are given from SignalDef.
// See src/Utils/SignalDef.cxx for more.
int nuPDG = event->PartInfo(0)->fPID;
// Generic signal flags
flagCCINC = SignalDef::isCCINC(event, nuPDG);
flagNCINC = SignalDef::isNCINC(event, nuPDG);
flagCCQE = SignalDef::isCCQE(event, nuPDG);
flagCCQELike = SignalDef::isCCQELike(event, nuPDG);
flagCC0pi = SignalDef::isCC0pi(event, nuPDG);
flagNCEL = SignalDef::isNCEL(event, nuPDG);
flagNC0pi = SignalDef::isNC0pi(event, nuPDG);
flagCCcoh = SignalDef::isCCCOH(event, nuPDG, 211);
flagNCcoh = SignalDef::isNCCOH(event, nuPDG, 111);
flagCC1pip = SignalDef::isCC1pi(event, nuPDG, 211);
flagNC1pip = SignalDef::isNC1pi(event, nuPDG, 211);
flagCC1pim = SignalDef::isCC1pi(event, nuPDG, -211);
flagNC1pim = SignalDef::isNC1pi(event, nuPDG, -211);
flagCC1pi0 = SignalDef::isCC1pi(event, nuPDG, 111);
flagNC1pi0 = SignalDef::isNC1pi(event, nuPDG, 111);
flagCC0piMINERvA = SignalDef::isCC0pi_MINERvAPTPZ(event, 14);
}
void GenericFlux_Vectors::AddSignalFlagsToTree() {
if (!eventVariables) {
Config::Get().out->cd();
eventVariables = new TTree((this->fName + "_VARS").c_str(),
(this->fName + "_VARS").c_str());
}
NUIS_LOG(SAM, "Adding signal flags");
// Signal Definitions from SignalDef.cxx
eventVariables->Branch("flagCCINC", &flagCCINC, "flagCCINC/O");
eventVariables->Branch("flagNCINC", &flagNCINC, "flagNCINC/O");
eventVariables->Branch("flagCCQE", &flagCCQE, "flagCCQE/O");
eventVariables->Branch("flagCC0pi", &flagCC0pi, "flagCC0pi/O");
eventVariables->Branch("flagCCQELike", &flagCCQELike, "flagCCQELike/O");
eventVariables->Branch("flagNCEL", &flagNCEL, "flagNCEL/O");
eventVariables->Branch("flagNC0pi", &flagNC0pi, "flagNC0pi/O");
eventVariables->Branch("flagCCcoh", &flagCCcoh, "flagCCcoh/O");
eventVariables->Branch("flagNCcoh", &flagNCcoh, "flagNCcoh/O");
eventVariables->Branch("flagCC1pip", &flagCC1pip, "flagCC1pip/O");
eventVariables->Branch("flagNC1pip", &flagNC1pip, "flagNC1pip/O");
eventVariables->Branch("flagCC1pim", &flagCC1pim, "flagCC1pim/O");
eventVariables->Branch("flagNC1pim", &flagNC1pim, "flagNC1pim/O");
eventVariables->Branch("flagCC1pi0", &flagCC1pi0, "flagCC1pi0/O");
eventVariables->Branch("flagNC1pi0", &flagNC1pi0, "flagNC1pi0/O");
eventVariables->Branch("flagCC0piMINERvA", &flagCC0piMINERvA, "flagCC0piMINERvA/O");
};
void GenericFlux_Vectors::Write(std::string drawOpt) {
// First save the TTree
eventVariables->Write();
// Save Flux and Event Histograms too
GetInput()->GetFluxHistogram()->Write();
GetInput()->GetEventHistogram()->Write();
return;
}
// Override functions which aren't really necessary
bool GenericFlux_Vectors::isSignal(FitEvent *event) {
(void)event;
return true;
};
void GenericFlux_Vectors::ScaleEvents() { return; }
void GenericFlux_Vectors::ApplyNormScale(float norm) {
this->fCurrentNorm = norm;
return;
}
void GenericFlux_Vectors::FillHistograms() { return; }
void GenericFlux_Vectors::ResetAll() {
eventVariables->Reset();
return;
}
float GenericFlux_Vectors::GetChi2() { return 0.0; }
diff --git a/src/MCStudies/GenericFlux_Vectors.h b/src/MCStudies/GenericFlux_Vectors.h
index aeae5c0..6dfe824 100644
--- a/src/MCStudies/GenericFlux_Vectors.h
+++ b/src/MCStudies/GenericFlux_Vectors.h
@@ -1,160 +1,163 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************/
#ifndef GenericFlux_Vectors_H_SEEN
#define GenericFlux_Vectors_H_SEEN
#include "Measurement1D.h"
#include "FitEvent.h"
#include "MINERvA_SignalDef.h"
class GenericFlux_Vectors : public Measurement1D {
public:
GenericFlux_Vectors(std::string name, std::string inputfile, FitWeight *rw, std::string type, std::string fakeDataFile);
virtual ~GenericFlux_Vectors() {};
//! Grab info from event
void FillEventVariables(FitEvent *event);
//! Fill signal flags
void FillSignalFlags(FitEvent *event);
void ResetVariables();
//! Fill Custom Histograms
void FillHistograms();
//! ResetAll
void ResetAll();
//! Scale
void ScaleEvents();
//! Norm
void ApplyNormScale(float norm);
//! Define this samples signal
bool isSignal(FitEvent *nvect);
//! Write Files
void Write(std::string drawOpt);
//! Get Chi2
float GetChi2();
void AddEventVariablesToTree();
void AddSignalFlagsToTree();
private:
TTree* eventVariables;
std::vector<FitParticle*> partList;
std::vector<FitParticle*> initList;
std::vector<FitParticle*> vertList;
bool SavePreFSI;
int Mode;
bool cc;
int PDGnu;
int tgt;
int tgta;
int tgtz;
int PDGLep;
float ELep;
float CosLep;
// Basic interaction kinematics
float Q2;
float q0;
float q3;
float Enu_QE;
float Enu_true;
float Q2_QE;
float W_nuc_rest;
float W;
float x;
float y;
float Eav;
float EavAlt;
float dalphat;
float W_genie;
float dpt;
float dphit;
float pnreco_C;
+ float CosThetaAdler;
+ float PhiAdler;
+
// Save outgoing particle vectors
int nfsp;
static const int kMAX = 200;
float px[kMAX];
float py[kMAX];
float pz[kMAX];
float E[kMAX];
int pdg[kMAX];
int pdg_rank[kMAX];
// Save incoming particle info
int ninitp;
float px_init[kMAX];
float py_init[kMAX];
float pz_init[kMAX];
float E_init[kMAX];
int pdg_init[kMAX];
// Save pre-FSI particle info
int nvertp;
float px_vert[kMAX];
float py_vert[kMAX];
float pz_vert[kMAX];
float E_vert[kMAX];
int pdg_vert[kMAX];
// Basic event info
float Weight;
float InputWeight;
float RWWeight;
double fScaleFactor;
// Custom weights
float CustomWeight;
float CustomWeightArray[6];
// Generic signal flags
bool flagCCINC;
bool flagNCINC;
bool flagCCQE;
bool flagCC0pi;
bool flagCCQELike;
bool flagNCEL;
bool flagNC0pi;
bool flagCCcoh;
bool flagNCcoh;
bool flagCC1pip;
bool flagNC1pip;
bool flagCC1pim;
bool flagNC1pim;
bool flagCC1pi0;
bool flagNC1pi0;
bool flagCC0piMINERvA;
};
#endif
diff --git a/src/MINERvA/MINERvA_CC0pinp_STV_XSec_1D_nu.cxx b/src/MINERvA/MINERvA_CC0pinp_STV_XSec_1D_nu.cxx
index 58a2413..06646d9 100644
--- a/src/MINERvA/MINERvA_CC0pinp_STV_XSec_1D_nu.cxx
+++ b/src/MINERvA/MINERvA_CC0pinp_STV_XSec_1D_nu.cxx
@@ -1,359 +1,353 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************/
// Implementation of 2018 MINERvA numu CC0pi STV analysis
// arxiv 1805.05486.pdf
// Clarence Wret
// cwret@fnal.gov
// Stephen Dolan
// Stephen.Dolan@llr.in2p3.fr
#include "MINERvA_CC0pinp_STV_XSec_1D_nu.h"
#include "MINERvA_SignalDef.h"
//********************************************************************
void MINERvA_CC0pinp_STV_XSec_1D_nu::SetupDataSettings() {
//********************************************************************
// Set Distribution
// See header file for enum and some descriptions
std::string name = fSettings.GetS("name");
if (!name.compare("MINERvA_CC0pinp_STV_XSec_1Dpmu_nu"))
fDist = kMuonMom;
else if (!name.compare("MINERvA_CC0pinp_STV_XSec_1Dthmu_nu"))
fDist = kMuonTh;
else if (!name.compare("MINERvA_CC0pinp_STV_XSec_1Dpprot_nu"))
fDist = kPrMom;
else if (!name.compare("MINERvA_CC0pinp_STV_XSec_1Dthprot_nu"))
fDist = kPrTh;
else if (!name.compare("MINERvA_CC0pinp_STV_XSec_1Dpnreco_nu"))
fDist = kNeMom;
else if (!name.compare("MINERvA_CC0pinp_STV_XSec_1Ddalphat_nu"))
fDist = kDalphaT;
else if (!name.compare("MINERvA_CC0pinp_STV_XSec_1Ddpt_nu"))
fDist = kDpT;
else if (!name.compare("MINERvA_CC0pinp_STV_XSec_1Ddphit_nu"))
fDist = kDphiT;
// Location of data, correlation matrices and the titles
std::string titles = "MINERvA_CC0pinp_STV_XSec_1D";
std::string foldername;
std::string distdescript;
// Data release is a single file
std::string rootfile = "MINERvA_1805.05486.root";
fMin = -999;
fMax = 999;
switch (fDist) {
case (kMuonMom):
titles += "pmu";
foldername = "muonmomentum";
distdescript = "Muon momentum in lab frame";
/*
fMin = 2.0;
fMax = 6.0;
*/
fMin = 1.5;
fMax = 10.0;
break;
case (kMuonTh):
titles += "thmu";
foldername = "muontheta";
distdescript = "Muon angle relative neutrino in lab frame";
fMin = 0.0;
fMax = 20.0;
break;
case (kPrMom):
titles += "pprot";
foldername = "protonmomentum";
distdescript = "Proton momentum in lab frame";
// fMin = 0.5;
fMin = 0.45;
fMax = 1.2;
break;
case (kPrTh):
titles += "thprot";
foldername = "protontheta";
distdescript = "Proton angle relative neutrino in lab frame";
fMin = 0.0;
fMax = 70.0;
break;
case (kNeMom):
titles += "pnreco";
foldername = "neutronmomentum";
distdescript = "Neutron momentum in lab frame";
fMin = 0.0;
// fMax = 0.9;
fMax = 2.0;
break;
case (kDalphaT):
foldername = "dalphat";
titles += foldername;
distdescript = "Delta Alpha_T";
fMin = 0.0;
// fMax = 170;
fMax = 180;
break;
case (kDpT):
foldername = "dpt";
titles += foldername;
distdescript = "Delta p_T";
fMin = 0.0;
fMax = 2.0;
break;
case (kDphiT):
foldername = "dphit";
titles += foldername;
distdescript = "Delta phi_T";
fMin = 0.0;
// fMax = 60.0;
fMax = 180.0;
break;
default:
NUIS_ERR(FTL,
"Did not find your specified distribution implemented, exiting");
NUIS_ABORT("You gave " << fName);
}
titles += "_nu";
// All have the same name
std::string dataname = foldername;
// Sample overview ---------------------------------------------------
std::string descrip = distdescript +
"Target: CH \n"
"Flux: MINERvA Forward Horn Current numu ONLY \n"
"Signal: Any event with 1 muon, and 0pi0 in FS, no "
"mesons, at least one proton with: \n"
"1.5GeV < p_mu < 10 GeV\n"
"theta_mu < 20 degrees\n"
"0.45GeV < p_prot < 1.2 GeV\n"
"theta_prot < 70 degrees\n"
"arXiv 1805.05486";
fSettings.SetDescription(descrip);
std::string filename =
GeneralUtils::GetTopLevelDir() +
"/data/MINERvA/CC0pi/CC0pi_STV/MINERvA_1805.05486.root";
// Specify the data
fSettings.SetDataInput(filename + ";" + dataname);
// And the correlations
fSettings.SetCovarInput(filename + ";" + dataname);
// Set titles
fSettings.SetTitle(titles);
};
//********************************************************************
MINERvA_CC0pinp_STV_XSec_1D_nu::MINERvA_CC0pinp_STV_XSec_1D_nu(
nuiskey samplekey) {
//********************************************************************
// A few different distributinos
// Muon momentum, muon angle, proton momentum, proton angle, neutron momentum,
// dalphat, dpt, dphit
+
+ // Hard-code the cuts
+ ProtonMinCut = 450; // MeV
+ ProtonMaxCut = 1200; // MeV
+ ProtonThetaCut = 70; // degrees
+
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
// Load up the data paths and sample descriptions
SetupDataSettings();
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG");
// No Enu cut
fSettings.SetEnuRange(0.0, 100.0);
fSettings.DefineAllowedTargets("C,H");
fSettings.DefineAllowedSpecies("numu");
// Finalise the settings
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = GetEventHistogram()->Integral("width") * double(1E-38) /
(double(fNEvents) * TotalIntegratedFlux("width"));
// Set the data and covariance matrix
SetDataFromRootFile(fSettings.GetDataInput());
SetCovarianceFromRootFile(fSettings.GetCovarInput());
fSettings.SetXTitle(fDataHist->GetXaxis()->GetTitle());
fSettings.SetYTitle(fDataHist->GetYaxis()->GetTitle());
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
//********************************************************************
// Data comes in a TList
// Some bins need stripping out because of zero bin content. Why oh why
void MINERvA_CC0pinp_STV_XSec_1D_nu::SetDataFromRootFile(std::string filename) {
//********************************************************************
std::vector<std::string> tempfile = GeneralUtils::ParseToStr(filename, ";");
TFile *File = new TFile(tempfile[0].c_str(), "READ");
// First object is the data
TH1D *temp = (TH1D *)(((TList *)(File->Get(tempfile[1].c_str())))->At(0));
// Garh, some of the data points are zero in the TH1D (WHY?!) so messes with
// the covariance entries to data bins check Skim through the data and check
// for zero bins
std::vector<double> CrossSection;
std::vector<double> Error;
std::vector<double> BinEdges;
int lastbin = 0;
startbin = 0;
for (int i = 0; i < temp->GetXaxis()->GetNbins() + 2; ++i) {
if (temp->GetBinContent(i + 1) > 0 && temp->GetBinLowEdge(i + 1) >= fMin &&
temp->GetBinLowEdge(i + 1) <= fMax) {
if (startbin == 0)
startbin = i;
lastbin = i;
CrossSection.push_back(temp->GetBinContent(i + 1));
BinEdges.push_back(temp->GetXaxis()->GetBinLowEdge(i + 1));
Error.push_back(temp->GetBinError(i + 1));
}
}
BinEdges.push_back(temp->GetXaxis()->GetBinLowEdge(lastbin + 2));
fDataHist = new TH1D((fSettings.GetName() + "_data").c_str(),
(fSettings.GetFullTitles()).c_str(), BinEdges.size() - 1,
&BinEdges[0]);
fDataHist->SetDirectory(0);
for (unsigned int i = 0; i < BinEdges.size() - 1; ++i) {
fDataHist->SetBinContent(i + 1, CrossSection[i]);
fDataHist->SetBinError(i + 1, Error[i]);
}
fDataHist->GetXaxis()->SetTitle(temp->GetXaxis()->GetTitle());
fDataHist->GetYaxis()->SetTitle(temp->GetYaxis()->GetTitle());
fDataHist->SetTitle(temp->GetTitle());
File->Close();
}
//********************************************************************
// Covariance also needs stripping out
// There's padding (two bins...) and overflow (last bin before the two empty
// bins)
void MINERvA_CC0pinp_STV_XSec_1D_nu::SetCovarianceFromRootFile(
std::string filename) {
//********************************************************************
std::vector<std::string> tempfile = GeneralUtils::ParseToStr(filename, ";");
TFile *File = new TFile(tempfile[0].c_str(), "READ");
// First object is the data, second is data with statistical error only, third
// is the covariance matrix
TMatrixDSym *tempcov =
(TMatrixDSym *)((TList *)File->Get(tempfile[1].c_str()))->At(2);
// Count the number of zero entries
int ngood = 0;
int nstart = -1;
int nend = -1;
// Scan through the middle bin and look for entries
int middle = tempcov->GetNrows() / 2;
int nbinsdata = fDataHist->GetXaxis()->GetNbins();
for (int j = 0; j < tempcov->GetNrows(); ++j) {
if ((*tempcov)(middle, j) > 0 && ngood < nbinsdata) {
ngood++;
if (nstart == -1)
nstart = j;
if (j > nend)
nend = j;
}
}
fFullCovar = new TMatrixDSym(ngood);
for (int i = 0; i < fFullCovar->GetNrows(); ++i) {
for (int j = 0; j < fFullCovar->GetNrows(); ++j) {
(*fFullCovar)(i, j) =
(*tempcov)(i + nstart + startbin - 1, j + nstart + startbin - 1);
}
}
(*fFullCovar) *= 1E38 * 1E38;
File->Close();
// Fill other covars.
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
}
void MINERvA_CC0pinp_STV_XSec_1D_nu::FillEventVariables(FitEvent *event) {
fXVar = -999.99;
// Need a proton and a muon
if (event->NumFSParticle(2212) == 0 || event->NumFSParticle(13) == 0) {
return;
}
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
// Find the highest momentum proton in the event between 450 and 1200 MeV with
// theta_p < 70
- int HMFSProton = 0;
- double HighestMomentum = 0.0;
- // Get the stack of protons
- std::vector<FitParticle *> Protons = event->GetAllFSProton();
- for (size_t i = 0; i < Protons.size(); ++i) {
- if (Protons[i]->p() > 450 && Protons[i]->p() < 1200 &&
- Protons[i]->P3().Angle(Pnu.Vect()) < (M_PI / 180.0) * 70.0 &&
- Protons[i]->p() > HighestMomentum) {
- HMFSProton = i;
- }
- }
- // Now get the proton
- TLorentzVector Pprot = Protons[HMFSProton]->fP;
+ TLorentzVector Pprot = FitUtils::GetProtonInRange(event, ProtonMinCut, ProtonMaxCut, cos(ProtonThetaCut/180.0*M_PI));
switch (fDist) {
case (kMuonMom):
fXVar = Pmu.Vect().Mag() / 1000.0;
break;
case (kMuonTh):
fXVar = Pmu.Vect().Angle(Pnu.Vect()) * (180.0 / M_PI);
break;
case (kPrMom):
fXVar = Pprot.Vect().Mag() / 1000.0;
break;
case (kPrTh):
fXVar = Pprot.Vect().Angle(Pnu.Vect()) * (180.0 / M_PI);
break;
// Use Stephen's validated functions
case (kNeMom):
- fXVar = FitUtils::Get_pn_reco_C(event, 14, true);
+ fXVar = FitUtils::Get_pn_reco_C(event, ProtonMinCut, ProtonMaxCut, cos(ProtonThetaCut/180.0*M_PI), 14, true);
break;
case (kDalphaT):
- fXVar = FitUtils::Get_STV_dalphat(event, 14, true) * (180.0 / M_PI);
+ fXVar = FitUtils::Get_STV_dalphat(event, ProtonMinCut, ProtonMaxCut, cos(ProtonThetaCut/180.0*M_PI), 14, true) * (180.0 / M_PI);
break;
case (kDpT):
- fXVar = FitUtils::Get_STV_dpt(event, 14, true) / 1000.0;
+ fXVar = FitUtils::Get_STV_dpt(event, ProtonMinCut, ProtonMaxCut, cos(ProtonThetaCut/180.0*M_PI), 14, true) / 1000.0;
break;
case (kDphiT):
- fXVar = FitUtils::Get_STV_dphit(event, 14, true) * (180.0 / M_PI);
+ fXVar = FitUtils::Get_STV_dphit(event, ProtonMinCut, ProtonMaxCut, cos(ProtonThetaCut/180.0*M_PI), 14, true) * (180.0 / M_PI);
break;
}
return;
};
//********************************************************************
bool MINERvA_CC0pinp_STV_XSec_1D_nu::isSignal(FitEvent *event)
//********************************************************************
{
return SignalDef::isCC0piNp_MINERvA_STV(event, EnuMin, EnuMax);
}
diff --git a/src/MINERvA/MINERvA_CC0pinp_STV_XSec_1D_nu.h b/src/MINERvA/MINERvA_CC0pinp_STV_XSec_1D_nu.h
index 5b88e40..e769b0b 100644
--- a/src/MINERvA/MINERvA_CC0pinp_STV_XSec_1D_nu.h
+++ b/src/MINERvA/MINERvA_CC0pinp_STV_XSec_1D_nu.h
@@ -1,58 +1,63 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************/
#ifndef MINERvA_CC0PINP_STV_XSEC_1D_NU_H_SEEN
#define MINERvA_CC0PINP_STV_XSEC_1D_NU_H_SEEN
#include "Measurement1D.h"
class MINERvA_CC0pinp_STV_XSec_1D_nu : public Measurement1D {
public:
MINERvA_CC0pinp_STV_XSec_1D_nu(nuiskey samplekey);
virtual ~MINERvA_CC0pinp_STV_XSec_1D_nu() {};
void FillEventVariables(FitEvent *event);
bool isSignal(FitEvent *event);
private:
+
+ double ProtonMinCut;
+ double ProtonMaxCut;
+ double ProtonThetaCut;
+
void SetupDataSettings();
void SetDataFromRootFile(std::string);
void SetCovarianceFromRootFile(std::string);
enum DataDistribution {
kMuonMom,
kMuonTh,
kPrMom,
kPrTh,
kNeMom,
kDalphaT,
kDpT,
kDphiT
} MINERvA_CC0pi_STV_DataDistributions;
DataDistribution fDist;
// For truncating
double fMin;
double fMax;
int startbin;
};
#endif
diff --git a/src/MINERvA/MINERvA_SignalDef.cxx b/src/MINERvA/MINERvA_SignalDef.cxx
index bd762df..59f385f 100644
--- a/src/MINERvA/MINERvA_SignalDef.cxx
+++ b/src/MINERvA/MINERvA_SignalDef.cxx
@@ -1,476 +1,471 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "SignalDef.h"
#include "FitUtils.h"
#include "MINERvA_SignalDef.h"
namespace SignalDef {
// *********************************
// MINERvA CC1pi+/- signal definition (2015 release)
// Note: There is a 2016 release which is different to this (listed below), but
// it is CCNpi+ and has a different W cut
// Note2: The W cut is implemented in the class implementation in MINERvA/
// rather than here so we can draw events that don't pass the W cut (W cut is
// 1.4 GeV)
// Could possibly be changed for slight speed increase since less events
// would be used
//
// MINERvA signal is slightly different to MiniBooNE
//
// Exactly one negative muon
// Exactly one charged pion (both + and -); however, there is a Michel e-
// requirement but UNCLEAR IF UNFOLDED OR NOT (so don't know if should be
// applied)
// Exactly 1 charged pion exits (so + and - charge), however, has Michel
// electron requirement, so look for + only here?
// No restriction on neutral pions or other mesons
// MINERvA has unfolded and not unfolded muon phase space for 2015
//
// Possible issues with the data:
// 1) pi- is allowed in signal even when Michel cut included; most pi- is efficiency corrected in GENIE
// 2) There is a T_pi < 350 MeV cut coming from requiring a stopping pion; this is efficiency corrected in GENIE
// 3) There is a 1.5 < Enu < 10.0 cut in signal definition
// 4) There is an angular muon cut which is sometimes efficiency corrected (why we have bool isRestricted below)
//
// Nice things:
// Much data given: with and without muon angle cuts and with and without shape
// only data + covariance
//
bool isCC1pip_MINERvA(FitEvent *event, double EnuMin, double EnuMax,
bool isRestricted) {
// *********************************
// Signal is both pi+ and pi-
// WARNING: PI- CONTAMINATION IS FULLY GENIE BECAUSE THE MICHEL TAG
// First, make sure it's CCINC
if (!isCCINC(event, 14, EnuMin, EnuMax)) return false;
// Allow pi+/pi-
int piPDG[] = {211, -211};
int nLeptons = event->NumFSLeptons();
int nPion = event->NumFSParticle(piPDG);
// Check that the desired pion exists and is the only meson
if (nPion != 1) return false;
// Check that there is only one final state lepton
if (nLeptons != 1) return false;
// MINERvA released another CC1pi+ xsec without muon unfolding!
// here the muon angle is < 20 degrees (seen in MINOS)
TLorentzVector pnu = event->GetHMISParticle(14)->fP;
TLorentzVector pmu = event->GetHMFSParticle(13)->fP;
if (isRestricted) {
double th_nu_mu = FitUtils::th(pmu, pnu) * 180. / M_PI;
if (th_nu_mu >= 20) return false;
}
// Extract Hadronic Mass
double hadMass = FitUtils::Wrec(pnu, pmu);
// Actual cut is True GENIE Ws! Arg.! Use gNtpcConv definition.
#ifdef __GENIE_ENABLED__
if (event->fType == kGENIE){
EventRecord * gevent = static_cast<EventRecord*>(event->genie_event->event);
const Interaction * interaction = gevent->Summary();
const Kinematics & kine = interaction->Kine();
double Ws = kine.W (true);
// std::cout << "Ws versus WRec = " << Ws << " vs " << hadMass << " " << kine.W(false) << std::endl;
hadMass = Ws * 1000.0;
}
#endif
if (hadMass > 1400.0) return false;
return true;
};
// Updated MINERvA 2017 Signal using Wexp and no restriction on angle
bool isCC1pip_MINERvA_2017(FitEvent *event, double EnuMin, double EnuMax){
// Signal is both pi+ and pi-
// WARNING: PI- CONTAMINATION IS FULLY GENIE BECAUSE THE MICHEL TAG
// First, make sure it's CCINC
if (!isCCINC(event, 14, EnuMin, EnuMax)) return false;
// Allow pi+/pi-
int piPDG[] = {211, -211};
int nLeptons = event->NumFSLeptons();
int nPion = event->NumFSParticle(piPDG);
// Check that the desired pion exists and is the only meson
if (nPion != 1) return false;
// Check that there is only one final state lepton
if (nLeptons != 1) return false;
// Get Muon and Lepton Kinematics
TLorentzVector pnu = event->GetHMISParticle(14)->fP;
TLorentzVector pmu = event->GetHMFSParticle(13)->fP;
// Extract Hadronic Mass
double hadMass = FitUtils::Wrec(pnu, pmu);
// Cut on 2017 data is still 1.4 GeV
if (hadMass > 1400.0) return false;
return true;
};
// *********************************
// MINERvA CCNpi+/- signal definition from 2016 publication
// Different to CC1pi+/- listed above; additional has W < 1.8 GeV
//
// For notes on strangeness of signal definition, see CC1pip_MINERvA
//
// One negative muon
// At least one charged pion
// 1.5 < Enu < 10
// No restrictions on pi0 or other mesons or baryons
// W_reconstructed (ignoring initial state motion) cut at 1.8 GeV
//
// Also writes number of pions (nPions) if studies on this want to be done...
bool isCCNpip_MINERvA(FitEvent *event, double EnuMin,
double EnuMax, bool isRestricted, bool isWtrue) {
// *********************************
// First, make sure it's CCINC
if (!isCCINC(event, 14, EnuMin, EnuMax)) return false;
int nLeptons = event->NumFSLeptons();
// Write the number of pions to the measurement class...
// Maybe better to just do that inside the class?
int nPions = event->NumFSParticle(PhysConst::pdg_charged_pions);
// Check that there is a pion!
if (nPions == 0) return false;
// Check that there is only one final state lepton
if (nLeptons != 1) return false;
// Need the muon and the neutrino to check angles and W
TLorentzVector pnu = event->GetNeutrinoIn()->fP;
TLorentzVector pmu = event->GetHMFSParticle(13)->fP;
// MINERvA released some data with restricted muon angle
// Here the muon angle is < 20 degrees (seen in MINOS)
if (isRestricted) {
double th_nu_mu = FitUtils::th(pmu, pnu) * 180. / M_PI;
if (th_nu_mu >= 20.) return false;
}
// Lastly check the W cut (always at 1.8 GeV)
double Wrec = FitUtils::Wrec(pnu, pmu) + 0.;
// Actual cut is True GENIE Ws! Arg.! Use gNtpcConv definition.
if (isWtrue){
#ifdef __GENIE_ENABLED__
if (event->fType == kGENIE){
GHepRecord* ghep = static_cast<GHepRecord*>(event->genie_event->event);
const Interaction * interaction = ghep->Summary();
const Kinematics & kine = interaction->Kine();
double Ws = kine.W (true);
Wrec = Ws * 1000.0; // Say Wrec is Ws
}
#endif
}
if (Wrec > 1800. || Wrec < 0.0) return false;
return true;
};
//********************************************************************
bool isCCQEnumu_MINERvA(FitEvent *event, double EnuMin, double EnuMax,
bool fullphasespace) {
//********************************************************************
if (!isCCQELike(event, 14, EnuMin, EnuMax)) return false;
TLorentzVector pnu = event->GetHMISParticle(14)->fP;
TLorentzVector pmu = event->GetHMFSParticle(13)->fP;
double ThetaMu = pnu.Vect().Angle(pmu.Vect());
double Enu_rec = FitUtils::EnuQErec(pmu, cos(ThetaMu), 34., true);
// If Restricted phase space
if (!fullphasespace && ThetaMu > 0.34906585) return false;
// restrict energy range
if (Enu_rec < EnuMin || Enu_rec > EnuMax) return false;
return true;
};
//********************************************************************
bool isCCQEnumubar_MINERvA(FitEvent *event, double EnuMin, double EnuMax,
bool fullphasespace) {
//********************************************************************
if (!isCCQELike(event, -14, EnuMin, EnuMax)) return false;
TLorentzVector pnu = event->GetHMISParticle(-14)->fP;
TLorentzVector pmu = event->GetHMFSParticle(-13)->fP;
double ThetaMu = pnu.Vect().Angle(pmu.Vect());
double Enu_rec = FitUtils::EnuQErec(pmu, cos(ThetaMu), 30., true);
// If Restricted phase space
if (!fullphasespace && ThetaMu > 0.34906585) return false;
// restrict energy range
if (Enu_rec < EnuMin || Enu_rec > EnuMax) return false;
return true;
}
//********************************************************************
bool isCCincLowRecoil_MINERvA(FitEvent *event, double EnuMin, double EnuMax) {
//********************************************************************
if (!isCCINC(event, 14, EnuMin, EnuMax)) return false;
// Need at least one muon
if (event->NumFSParticle(13) < 1) return false;
TLorentzVector pmu = event->GetHMFSParticle(13)->fP;
TLorentzVector pnu = event->GetHMISParticle(14)->fP;
// Cut on muon angle greated than 20deg
if (cos(pnu.Vect().Angle(pmu.Vect())) < 0.93969262078) return false;
// Cut on muon energy < 1.5 GeV
if (pmu.E()/1000.0 < 1.5) return false;
return true;
}
// Used in 2014 muon+proton analysis
// Events with muon angles up to 70 degrees
// One right sign muon, at least one proton, no pions
// proton kinetic energies greater than 100 MeV
bool isCC0pi1p_MINERvA(FitEvent *event, double enumin, double enumax) {
bool signal = (
isCC0pi(event, 14, enumin, enumax) && // Require numu CC0pi event
HasProtonKEAboveThreshold(event, 110.0) && // With proton above threshold
(event->GetHMFSMuon())->P3().Angle(
(event->GetNeutrinoIn())->P3())*180./M_PI < 70 // And muon within production angle
);
return signal;
}
// 2015 analysis just asks for 1pi0 and no pi+/pi-
bool isCC1pi0_MINERvA_2015(FitEvent *event, double EnuMin, double EnuMax) {
bool CC1pi0_anu = SignalDef::isCC1pi(event, -14, 111, EnuMin, EnuMax);
return CC1pi0_anu;
}
// 2016 analysis just asks for 1pi0 and no other charged tracks. Half-open to interpretation: we go with "charged tracks" meaning pions. You'll be forgiven for thinking proton tracks should be included here too but we checked with MINERvA
bool isCC1pi0_MINERvA_2016(FitEvent *event, double EnuMin, double EnuMax) {
bool CC1pi0_anu = SignalDef::isCC1pi(event, -14, 111, EnuMin, EnuMax);
/*
// Additionally look for charged proton track
// Comment: This is _NOT_ in the signal definition but was tested
bool HasProton = event->HasFSParticle(2212);
if (CC1pi0_anu) {
if (!HasProton) {
return true;
} else {
return false;
}
} else {
return false;
}
*/
return CC1pi0_anu;
}
// 2016 analysis just asks for 1pi0 and no other charged tracks
bool isCC1pi0_MINERvA_nu(FitEvent *event, double EnuMin, double EnuMax) {
bool CC1pi0_nu = SignalDef::isCC1pi(event, 14, 111, EnuMin, EnuMax);
return CC1pi0_nu;
}
//********************************************************************
bool isCC0pi_MINERvAPTPZ(FitEvent* event, int nuPDG, double emin, double emax){
//********************************************************************
// Check it's CCINC
if (!SignalDef::isCCINC(event, nuPDG, emin, emax)) return false;
// Make Angle Cut > 20.0
TLorentzVector pnu = event->GetHMISParticle(14)->fP;
TLorentzVector pmu = event->GetHMFSParticle(13)->fP;
double th_nu_mu = FitUtils::th(pmu, pnu) * 180. / M_PI;
if (th_nu_mu >= 20.0) return false;
int genie_n_muons = 0;
int genie_n_mesons = 0;
int genie_n_heavy_baryons_plus_pi0s = 0;
int genie_n_photons = 0;
for(unsigned int i = 0; i < event->NParticles(); ++i) {
FitParticle* p = event->GetParticle(i);
if (p->Status() != kFinalState) continue;
int pdg = p->fPID;
double energy = p->fP.E();
if( pdg == 13 ) {
genie_n_muons++;
} else if( pdg == 22 && energy > 10.0 ) {
genie_n_photons++;
} else if( abs(pdg) == 211 || abs(pdg) == 321 || abs(pdg) == 323 || pdg == 111 || pdg == 130 || pdg == 310 || pdg == 311 || pdg == 313 || abs(pdg) == 221 || abs(pdg) == 331 ) {
genie_n_mesons++;
} else if( pdg == 3112 || pdg == 3122 || pdg == 3212 || pdg == 3222 || pdg == 4112 ||
pdg == 4122 || pdg == 4212 || pdg == 4222 || pdg == 411 || pdg == 421 || pdg == 111 ) {
genie_n_heavy_baryons_plus_pi0s++;
}
}
if( genie_n_muons == 1 &&
genie_n_mesons == 0 &&
genie_n_heavy_baryons_plus_pi0s == 0 &&
genie_n_photons == 0 ) return true;
return false;
}
// **************************************************
// Section VI Event Selection of https://journals.aps.org/prd/pdf/10.1103/PhysRevD.97.052002
// Anti-neutrino charged-current
// Post-FSI final states without
// mesons,
// prompt photons above nuclear deexcitation energies
// heavy baryons
// protons above kinetic energy of 120 MeV
// Muon-neutrino angle of 20 degrees
// Parallel muon momentum: 1.5 GeV < P|| < 15 GeV --- N.B. APPARENTLY NOT INCLUDED, see below
// Transverse muon momentum: pT < 1.5 GeV --- N.B. APPARENTLY NOT INCLUDED, see below
bool isCC0pi_anti_MINERvAPTPZ(FitEvent* event, int nuPDG, double emin, double emax){
// **************************************************
// Check it's CCINC
if (!SignalDef::isCCINC(event, nuPDG, emin, emax)) return false;
TLorentzVector pnu = event->GetNeutrinoIn()->fP;
TLorentzVector pmu = event->GetHMFSParticle(-13)->fP;
// Make Angle Cut > 20.0
double th_nu_mu = FitUtils::th(pmu, pnu) * 180. / M_PI;
if (th_nu_mu >= 20.0) return false;
// Heidi Schellman (schellmh@science.oregonstate.edu) assured me that the p_t and p_z (or p_||) cuts aren't actually implemented as a signal definition: they're only implemented in the binning for p_t and p_z (but not Q2QE and EnuQE)
/*
// Cut on pT and pZ
Double_t px = pmu.X()/1.E3;
Double_t py = pmu.Y()/1.E3;
Double_t pt = sqrt(px*px+py*py);
// Don't want to assume the event generators all have neutrino coming along z
// pz is muon momentum projected onto the neutrino direction
Double_t pz = pmu.Vect().Dot(pnu.Vect()*(1.0/pnu.Vect().Mag()))/1.E3;
if (pz > 15 || pz < 1.5) return false;
if (pt > 1.5) return false;
*/
// Find if there are any protons above 120 MeV kinetic energy
if (HasProtonKEAboveThreshold(event, 120.0)) return false;
// Particle counters
int genie_n_muons = 0;
int genie_n_mesons = 0;
int genie_n_heavy_baryons_plus_pi0s = 0;
int genie_n_photons = 0;
// Loop over the particles in the event and count them up
for(unsigned int i = 0; i < event->NParticles(); ++i) {
FitParticle* p = event->GetParticle(i);
if (p->Status() != kFinalState) continue;
int pdg = p->fPID;
double energy = p->fP.E();
// Any charged muons
if( abs(pdg) == 13 ) {
genie_n_muons++;
// De-excitation photons
} else if( pdg == 22 && energy > 10.0 ) {
genie_n_photons++;
// Mesons
} else if( abs(pdg) == 211 || abs(pdg) == 321 || abs(pdg) == 323 || pdg == 111 || pdg == 130 || pdg == 310 || pdg == 311 || pdg == 313 || abs(pdg) == 221 || abs(pdg) == 331 ) {
genie_n_mesons++;
// Heavy baryons and pi0s
} else if( abs(pdg) == 3112 || abs(pdg) == 3122 || abs(pdg) == 3212 || abs(pdg) == 3222 || abs(pdg) == 4112 ||
abs(pdg) == 4122 || abs(pdg) == 4212 || abs(pdg) == 4222 || abs(pdg) == 411 || abs(pdg) == 421 ||
abs(pdg) == 111 ) {
genie_n_heavy_baryons_plus_pi0s++;
}
}
// Look for one muon with no mesons, heavy baryons or deexcitation photons
if( genie_n_muons == 1 && genie_n_mesons == 0 && genie_n_heavy_baryons_plus_pi0s == 0 && genie_n_photons == 0 ) return true;
return false;
}
// MINERvA CC0pi transverse variables signal defintion
bool isCC0piNp_MINERvA_STV(FitEvent *event, double EnuMin, double EnuMax) {
// Require a numu CC0pi event
if (!isCC0pi(event, 14, EnuMin, EnuMax)) return false;
// Require at least one FS proton
if (event->NumFSParticle(2212) == 0) return false;
TLorentzVector pnu = event->GetHMISParticle(14)->fP;
TLorentzVector pmu = event->GetHMFSParticle(13)->fP;
// Muon momentum cuts
if (pmu.Vect().Mag() < 1500 || pmu.Vect().Mag() > 10000) return false;
// Muon angle cuts
if (pmu.Vect().Angle(pnu.Vect()) > (M_PI/180.0)*20.0) return false;
// Since there are emany protons allowed we can't just use the highest proton momentum candidate and place cuts on it
// Get the stack of protons
std::vector<FitParticle*> Protons = event->GetAllFSProton();
// Count how many protons pass the threshold
int nProtonsAboveThreshold = 0;
for (size_t i = 0; i < Protons.size(); ++i) {
if (Protons[i]->p() > 450 &&
Protons[i]->p() < 1200 &&
Protons[i]->P3().Angle(pnu.Vect()) < (M_PI/180.0)*70.0) {
nProtonsAboveThreshold++;
}
}
- // Proton momentum cuts
- //if (pp.Vect().Mag() < 450 || pp.Vect().Mag() > 1200) return false;
- // Proton angle cuts
- //if (pp.Vect().Angle(pnu.Vect()) > (M_PI/180.0)*70) return false;
if (nProtonsAboveThreshold == 0) return false;
-
return true;
};
-}
+} // end namespace
diff --git a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddat_nu.cxx b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddat_nu.cxx
index c9d1b62..fbd7423 100644
--- a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddat_nu.cxx
+++ b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddat_nu.cxx
@@ -1,90 +1,94 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "T2K_SignalDef.h"
#include "T2K_CC0pinp_STV_XSec_1Ddat_nu.h"
//********************************************************************
T2K_CC0pinp_STV_XSec_1Ddat_nu::T2K_CC0pinp_STV_XSec_1Ddat_nu(
nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "T2K_CC0pinp_STV_XSec_1Ddpt_nu sample. \n"
"Target: CH \n"
"Flux: T2K 2.5 degree off-axis (ND280) \n"
"Signal: CC0piNp (N>=1) with 450 MeV < p_p < 1 GeV \n"
" p_mu > 250 MeV \n"
" cth_p > 0.4 \n"
" cth_mu > -0.6 \n"
"https://doi.org/10.1103/PhysRevD.98.032003 \n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("#delta#it{#alpha}_{T} (GeV c^{-1})");
fSettings.SetYTitle(
"#frac{d#sigma}{d#delta#it{#alpha}_{T}} (cm^{2} nucleon^{-1} rads^{-1})");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG");
fSettings.SetEnuRange(0.0, 50.0);
fSettings.DefineAllowedTargets("C,H");
// CCQELike plot information
fSettings.SetTitle("T2K_CC0pinp_STV_XSec_1Ddat_nu");
// fSettings.SetDataInput( GeneralUtils::GetTopLevelDir() +
// "/data/T2K/T2K_CC0pinp_STV_XSec_1Ddat_nu.dat");
fSettings.SetDataInput(FitPar::GetDataBase() +
"/T2K/CC0pi/STV/datResults.root;Result");
fSettings.SetCovarInput(FitPar::GetDataBase() +
"/T2K/CC0pi/STV/datResults.root;Correlation_Matrix");
fSettings.DefineAllowedSpecies("numu");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = GetEventHistogram()->Integral("width") * double(1E-38) /
(double(fNEvents) * TotalIntegratedFlux("width"));
// Plot Setup -------------------------------------------------------
// SetDataFromTextFile( fSettings.GetDataInput() );
// ScaleData(1E-38);
// SetCovarFromDiagonal();
SetDataFromRootFile(fSettings.GetDataInput());
SetCorrelationFromRootFile(fSettings.GetCovarInput());
+ // Proton cuts: 450 MeV to 1 GeV, cos(theta_proton_neutrino) < 0.4
+ ProtonMinCut = 450;
+ ProtonMaxCut = 1000;
+ ProtonCosThetaCut = 0.4;
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
void T2K_CC0pinp_STV_XSec_1Ddat_nu::FillEventVariables(FitEvent *event) {
- fXVar = FitUtils::Get_STV_dalphat(event, 14, true);
+ fXVar = FitUtils::Get_STV_dalphat(event, ProtonMinCut, ProtonMaxCut, ProtonCosThetaCut, 14, true);
return;
};
//********************************************************************
bool T2K_CC0pinp_STV_XSec_1Ddat_nu::isSignal(FitEvent *event)
//********************************************************************
{
return SignalDef::isT2K_CC0pi_STV(event, EnuMin, EnuMax);
}
diff --git a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddat_nu.h b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddat_nu.h
index 8432261..f226217 100644
--- a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddat_nu.h
+++ b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddat_nu.h
@@ -1,36 +1,40 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************/
#ifndef T2K_CC0pinp_STV_XSec_1Ddat_nu_H_SEEN
#define T2K_CC0pinp_STV_XSec_1Ddat_nu_H_SEEN
#include "Measurement1D.h"
class T2K_CC0pinp_STV_XSec_1Ddat_nu : public Measurement1D {
public:
T2K_CC0pinp_STV_XSec_1Ddat_nu(nuiskey samplekey);
virtual ~T2K_CC0pinp_STV_XSec_1Ddat_nu() {};
void FillEventVariables(FitEvent *event);
bool isSignal(FitEvent *event);
private:
+ // Proton cuts
+ double ProtonMinCut;
+ double ProtonMaxCut;
+ double ProtonCosThetaCut;
};
#endif
diff --git a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddphit_nu.cxx b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddphit_nu.cxx
index 9dda5ff..6f7e260 100644
--- a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddphit_nu.cxx
+++ b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddphit_nu.cxx
@@ -1,90 +1,93 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "T2K_SignalDef.h"
#include "T2K_CC0pinp_STV_XSec_1Ddphit_nu.h"
//********************************************************************
T2K_CC0pinp_STV_XSec_1Ddphit_nu::T2K_CC0pinp_STV_XSec_1Ddphit_nu(
nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "T2K_CC0pinp_STV_XSec_1Ddpt_nu sample. \n"
"Target: CH \n"
"Flux: T2K 2.5 degree off-axis (ND280) \n"
"Signal: CC0piNp (N>=1) with 450 MeV < p_p < 1 GeV \n"
" p_mu > 250 MeV \n"
" cth_p > 0.4 \n"
" cth_mu > -0.6 \n"
"https://doi.org/10.1103/PhysRevD.98.032003 \n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("#delta#it{#phi}_{T} (GeV c^{-1})");
fSettings.SetYTitle(
"#frac{d#sigma}{d#delta#it{#phi}_{T}} (cm^{2} nucleon^{-1} rads^{-1})");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG");
fSettings.SetEnuRange(0.0, 50.0);
fSettings.DefineAllowedTargets("C,H");
// CCQELike plot information
fSettings.SetTitle("T2K_CC0pinp_STV_XSec_1Ddphit_nu");
// fSettings.SetDataInput( GeneralUtils::GetTopLevelDir() +
// "/data/T2K/T2K_CC0pinp_STV_XSec_1Ddphit_nu.dat");
fSettings.SetDataInput(FitPar::GetDataBase() +
"/T2K/CC0pi/STV/dphitResults.root;Result");
fSettings.SetCovarInput(
FitPar::GetDataBase() +
"/T2K/CC0pi/STV/dphitResults.root;Correlation_Matrix");
fSettings.DefineAllowedSpecies("numu");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = GetEventHistogram()->Integral("width") * double(1E-38) /
(double(fNEvents) * TotalIntegratedFlux("width"));
// Plot Setup -------------------------------------------------------
// SetDataFromTextFile( fSettings.GetDataInput() );
// ScaleData(1E-38);
// SetCovarFromDiagonal();
SetDataFromRootFile(fSettings.GetDataInput());
SetCorrelationFromRootFile(fSettings.GetCovarInput());
+ ProtonMinCut = 450;
+ ProtonMaxCut = 1000;
+ ProtonCosThetaCut = 0.4;
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
void T2K_CC0pinp_STV_XSec_1Ddphit_nu::FillEventVariables(FitEvent *event) {
- fXVar = FitUtils::Get_STV_dphit(event, 14, true);
+ fXVar = FitUtils::Get_STV_dphit(event, ProtonMinCut, ProtonMaxCut, ProtonCosThetaCut, 14, true);
return;
};
//********************************************************************
bool T2K_CC0pinp_STV_XSec_1Ddphit_nu::isSignal(FitEvent *event)
//********************************************************************
{
return SignalDef::isT2K_CC0pi_STV(event, EnuMin, EnuMax);
}
diff --git a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddphit_nu.h b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddphit_nu.h
index ee17568..8f21503 100644
--- a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddphit_nu.h
+++ b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddphit_nu.h
@@ -1,36 +1,40 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************/
#ifndef T2K_CC0pinp_STV_XSec_1Ddphit_nu_H_SEEN
#define T2K_CC0pinp_STV_XSec_1Ddphit_nu_H_SEEN
#include "Measurement1D.h"
class T2K_CC0pinp_STV_XSec_1Ddphit_nu : public Measurement1D {
public:
T2K_CC0pinp_STV_XSec_1Ddphit_nu(nuiskey samplekey);
virtual ~T2K_CC0pinp_STV_XSec_1Ddphit_nu() {};
void FillEventVariables(FitEvent *event);
bool isSignal(FitEvent *event);
private:
+ // Proton cuts
+ double ProtonMinCut;
+ double ProtonMaxCut;
+ double ProtonCosThetaCut;
};
#endif
diff --git a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.cxx b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.cxx
index 65d6559..1f463ab 100644
--- a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.cxx
+++ b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.cxx
@@ -1,90 +1,93 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "T2K_SignalDef.h"
#include "T2K_CC0pinp_STV_XSec_1Ddpt_nu.h"
//********************************************************************
T2K_CC0pinp_STV_XSec_1Ddpt_nu::T2K_CC0pinp_STV_XSec_1Ddpt_nu(
nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "T2K_CC0pinp_STV_XSec_1Ddpt_nu sample. \n"
"Target: CH \n"
"Flux: T2K 2.5 degree off-axis (ND280) \n"
"Signal: CC0piNp (N>=1) with 450 MeV < p_p < 1 GeV \n"
" p_mu > 250 MeV \n"
" cth_p > 0.4 \n"
" cth_mu > -0.6 \n"
"https://doi.org/10.1103/PhysRevD.98.032003 \n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("#delta#it{p}_{T} (GeV c^{-1})");
fSettings.SetYTitle(
"#frac{d#sigma}{d#delta#it{p}_{T}} (cm^{2} nucleon^{-1} GeV^{-1} c)");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG");
fSettings.SetEnuRange(0.0, 50.0);
fSettings.DefineAllowedTargets("C,H");
// CCQELike plot information
fSettings.SetTitle("T2K_CC0pinp_STV_XSec_1Ddpt_nu");
// fSettings.SetDataInput( GeneralUtils::GetTopLevelDir() +
// "/data/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.dat");
fSettings.SetDataInput(FitPar::GetDataBase() +
"/T2K/CC0pi/STV/dptResults.root;Result");
fSettings.SetCovarInput(FitPar::GetDataBase() +
"/T2K/CC0pi/STV/dptResults.root;Correlation_Matrix");
fSettings.DefineAllowedSpecies("numu");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = GetEventHistogram()->Integral("width") * double(1E-38) /
(double(fNEvents) * TotalIntegratedFlux("width"));
// Plot Setup -------------------------------------------------------
// SetDataFromTextFile( fSettings.GetDataInput() );
// SetCovarFromDiagonal();
// ScaleData(1E-38);
SetDataFromRootFile(fSettings.GetDataInput());
SetCorrelationFromRootFile(fSettings.GetCovarInput());
// SetCovarianceFromRootFile(fSettings.GetCovarInput() );
+ ProtonMinCut = 450;
+ ProtonMaxCut = 1000;
+ ProtonCosThetaCut = 0.4;
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
void T2K_CC0pinp_STV_XSec_1Ddpt_nu::FillEventVariables(FitEvent *event) {
- fXVar = FitUtils::Get_STV_dpt(event, 14, true) / 1000.0;
+ fXVar = FitUtils::Get_STV_dpt(event, ProtonMinCut, ProtonMaxCut, ProtonCosThetaCut, 14, true) / 1000.0;
return;
};
//********************************************************************
bool T2K_CC0pinp_STV_XSec_1Ddpt_nu::isSignal(FitEvent *event)
//********************************************************************
{
return SignalDef::isT2K_CC0pi_STV(event, EnuMin, EnuMax);
}
diff --git a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.h b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.h
index 51e61dd..d828684 100644
--- a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.h
+++ b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.h
@@ -1,36 +1,40 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************/
#ifndef T2K_CC0pinp_STV_XSec_1Ddpt_nu_H_SEEN
#define T2K_CC0pinp_STV_XSec_1Ddpt_nu_H_SEEN
#include "Measurement1D.h"
class T2K_CC0pinp_STV_XSec_1Ddpt_nu : public Measurement1D {
public:
T2K_CC0pinp_STV_XSec_1Ddpt_nu(nuiskey samplekey);
virtual ~T2K_CC0pinp_STV_XSec_1Ddpt_nu() {};
void FillEventVariables(FitEvent *event);
bool isSignal(FitEvent *event);
private:
+ // Proton cuts
+ double ProtonMinCut;
+ double ProtonMaxCut;
+ double ProtonCosThetaCut;
};
#endif
diff --git a/src/Utils/FitUtils.cxx b/src/Utils/FitUtils.cxx
index 332cf4a..e658bae 100644
--- a/src/Utils/FitUtils.cxx
+++ b/src/Utils/FitUtils.cxx
@@ -1,1159 +1,1098 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "FitUtils.h"
/*
MISC Functions
*/
//********************************************************************
double *FitUtils::GetArrayFromMap(std::vector<std::string> invals,
std::map<std::string, double> inmap) {
//********************************************************************
double *outarr = new double[invals.size()];
int count = 0;
for (size_t i = 0; i < invals.size(); i++) {
outarr[count++] = inmap[invals.at(i)];
}
return outarr;
}
/*
MISC Event
*/
//********************************************************************
// Returns the kinetic energy of a particle in GeV
double FitUtils::T(TLorentzVector part) {
//********************************************************************
double E_part = part.E() / 1000.;
double p_part = part.Vect().Mag() / 1000.;
double m_part = sqrt(E_part * E_part - p_part * p_part);
double KE_part = E_part - m_part;
return KE_part;
};
//********************************************************************
// Returns the momentum of a particle in GeV
double FitUtils::p(TLorentzVector part) {
//********************************************************************
double p_part = part.Vect().Mag() / 1000.;
return p_part;
};
double FitUtils::p(FitParticle *part) { return FitUtils::p(part->fP); };
//********************************************************************
// Returns the angle between two particles in radians
double FitUtils::th(TLorentzVector part1, TLorentzVector part2) {
//********************************************************************
double th = part1.Vect().Angle(part2.Vect());
return th;
};
double FitUtils::th(FitParticle *part1, FitParticle *part2) {
return FitUtils::th(part1->fP, part2->fP);
};
// T2K CC1pi+ helper functions
//
//********************************************************************
// Returns the angle between q3 and the pion defined in Raquel's CC1pi+ on CH
// paper
// Uses "MiniBooNE formula" for Enu, here called EnuCC1pip_T2K_MB
//********************************************************************
double FitUtils::thq3pi_CC1pip_T2K(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi) {
// Want this in GeV
TVector3 p_mu = pmu.Vect() * (1. / 1000.);
// Get the reconstructed Enu
// We are not using Michel e sample, so we have pion kinematic information
double Enu = EnuCC1piprec(pnu, pmu, ppi, true);
// Get neutrino unit direction, multiply by reconstructed Enu
TVector3 p_nu = pnu.Vect() * (1. / (pnu.Vect().Mag())) * Enu;
TVector3 p_pi = ppi.Vect() * (1. / 1000.);
// This is now in GeV
TVector3 q3 = (p_nu - p_mu);
// Want this in GeV
double th_q3_pi = q3.Angle(p_pi);
return th_q3_pi;
}
//********************************************************************
// Returns the q3 defined in Raquel's CC1pi+ on CH paper
// Uses "MiniBooNE formula" for Enu
//********************************************************************
double FitUtils::q3_CC1pip_T2K(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi) {
// Can't use the true Enu here; need to reconstruct it
// We do have Michel e- here so reconstruct Enu by "MiniBooNE formula" without
// pion kinematics
// The bool false refers to that we don't have pion kinematics
// Last bool refers to if we have pion kinematic information or not
double Enu = EnuCC1piprec(pnu, pmu, ppi, false);
TVector3 p_mu = pmu.Vect() * (1. / 1000.);
TVector3 p_nu = pnu.Vect() * (1. / pnu.Vect().Mag()) * Enu;
double q3 = (p_nu - p_mu).Mag();
return q3;
}
//********************************************************************
// Returns the W reconstruction from Raquel CC1pi+ CH thesis
// Uses the MiniBooNE formula Enu
//********************************************************************
double FitUtils::WrecCC1pip_T2K_MB(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi) {
double E_mu = pmu.E() / 1000.;
double p_mu = pmu.Vect().Mag() / 1000.;
double E_nu = EnuCC1piprec(pnu, pmu, ppi, false);
double a1 = (E_nu + PhysConst::mass_neutron) - E_mu;
double a2 = E_nu - p_mu;
double wrec = sqrt(a1 * a1 - a2 * a2);
return wrec;
}
//********************************************************
double FitUtils::ProtonQ2QErec(double pE, double binding) {
//********************************************************
const double V = binding / 1000.; // binding potential
const double mn = PhysConst::mass_neutron; // neutron mass
const double mp = PhysConst::mass_proton; // proton mass
const double mn_eff = mn - V; // effective proton mass
const double pki = (pE / 1000.0) - mp;
double q2qe = mn_eff * mn_eff - mp * mp + 2 * mn_eff * (pki + mp - mn_eff);
return q2qe;
};
//********************************************************************
double FitUtils::EnuQErec(TLorentzVector pmu, double costh, double binding,
bool neutrino) {
//********************************************************************
// Convert all values to GeV
const double V = binding / 1000.; // binding potential
const double mn = PhysConst::mass_neutron; // neutron mass
const double mp = PhysConst::mass_proton; // proton mass
double mN_eff = mn - V;
double mN_oth = mp;
if (!neutrino) {
mN_eff = mp - V;
mN_oth = mn;
}
double el = pmu.E() / 1000.;
double pl = (pmu.Vect().Mag()) / 1000.; // momentum of lepton
double ml = sqrt(el * el - pl * pl); // lepton mass
double rEnu =
(2 * mN_eff * el - ml * ml + mN_oth * mN_oth - mN_eff * mN_eff) /
(2 * (mN_eff - el + pl * costh));
return rEnu;
};
// Another good old helper function
double FitUtils::EnuQErec(TLorentzVector pmu, TLorentzVector pnu, double binding, bool neutrino) {
return EnuQErec(pmu, cos(pnu.Vect().Angle(pmu.Vect())), binding, neutrino);
}
double FitUtils::Q2QErec(TLorentzVector pmu, double costh, double binding, bool neutrino) {
double el = pmu.E() / 1000.;
double pl = (pmu.Vect().Mag()) / 1000.; // momentum of lepton
double ml = sqrt(el * el - pl * pl); // lepton mass
double rEnu = EnuQErec(pmu, costh, binding, neutrino);
double q2 = -ml * ml + 2. * rEnu * (el - pl * costh);
return q2;
};
double FitUtils::Q2QErec(TLorentzVector Pmu, TLorentzVector Pnu, double binding, bool neutrino) {
double q2qe = Q2QErec(Pmu, cos(Pnu.Vect().Angle(Pmu.Vect())), binding, neutrino);
return q2qe;
}
double FitUtils::EnuQErec(double pl, double costh, double binding,
bool neutrino) {
if (pl < 0) return 0.; // Make sure nobody is silly
double mN_eff = PhysConst::mass_neutron - binding / 1000.;
double mN_oth = PhysConst::mass_proton;
if (!neutrino) {
mN_eff = PhysConst::mass_proton - binding / 1000.;
mN_oth = PhysConst::mass_neutron;
}
double ml = PhysConst::mass_muon;
double el = sqrt(pl * pl + ml * ml);
double rEnu =
(2 * mN_eff * el - ml * ml + mN_oth * mN_oth - mN_eff * mN_eff) /
(2 * (mN_eff - el + pl * costh));
return rEnu;
};
double FitUtils::Q2QErec(double pl, double costh, double binding,
bool neutrino) {
if (pl < 0) return 0.; // Make sure nobody is silly
double ml = PhysConst::mass_muon;
double el = sqrt(pl * pl + ml * ml);
double rEnu = EnuQErec(pl, costh, binding, neutrino);
double q2 = -ml * ml + 2. * rEnu * (el - pl * costh);
return q2;
};
//********************************************************************
// Reconstructs Enu for CC1pi0
// Very similar for CC1pi+ reconstruction
double FitUtils::EnuCC1pi0rec(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi0) {
//********************************************************************
double E_mu = pmu.E() / 1000;
double p_mu = pmu.Vect().Mag() / 1000;
double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu);
double th_nu_mu = pnu.Vect().Angle(pmu.Vect());
double E_pi0 = ppi0.E() / 1000;
double p_pi0 = ppi0.Vect().Mag() / 1000;
double m_pi0 = sqrt(E_pi0 * E_pi0 - p_pi0 * p_pi0);
double th_nu_pi0 = pnu.Vect().Angle(ppi0.Vect());
const double m_n = PhysConst::mass_neutron; // neutron mass
const double m_p = PhysConst::mass_proton; // proton mass
double th_pi0_mu = ppi0.Vect().Angle(pmu.Vect());
double rEnu = (m_mu * m_mu + m_pi0 * m_pi0 + m_n * m_n - m_p * m_p -
2 * m_n * (E_pi0 + E_mu) + 2 * E_pi0 * E_mu -
2 * p_pi0 * p_mu * cos(th_pi0_mu)) /
(2 * (E_pi0 + E_mu - p_pi0 * cos(th_nu_pi0) -
p_mu * cos(th_nu_mu) - m_n));
return rEnu;
};
//********************************************************************
// Reconstruct Q2 for CC1pi0
// Beware: uses true Enu, not reconstructed Enu
double FitUtils::Q2CC1pi0rec(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi0) {
//********************************************************************
double E_mu = pmu.E() / 1000.; // energy of lepton in GeV
double p_mu = pmu.Vect().Mag() / 1000.; // momentum of lepton
double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu); // lepton mass
double th_nu_mu = pnu.Vect().Angle(pmu.Vect());
// double rEnu = EnuCC1pi0rec(pnu, pmu, ppi0); //reconstructed neutrino energy
// Use true neutrino energy
double rEnu = pnu.E() / 1000.;
double q2 = -m_mu * m_mu + 2. * rEnu * (E_mu - p_mu * cos(th_nu_mu));
return q2;
};
//********************************************************************
// Reconstruct Enu for CC1pi+
// pionInfo reflects if we're using pion kinematics or not
// In T2K CC1pi+ CH the Michel tag is used for pion in which pion kinematic info
// is lost and Enu is reconstructed without pion kinematics
double FitUtils::EnuCC1piprec(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi, bool pionInfo) {
//********************************************************************
double E_mu = pmu.E() / 1000.;
double p_mu = pmu.Vect().Mag() / 1000.;
double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu);
double E_pi = ppi.E() / 1000.;
double p_pi = ppi.Vect().Mag() / 1000.;
double m_pi = sqrt(E_pi * E_pi - p_pi * p_pi);
const double m_n = PhysConst::mass_neutron; // neutron/proton mass
// should really take proton mass for proton interaction, neutron for neutron
// interaction. However, difference is pretty much negligable here!
// need this angle too
double th_nu_pi = pnu.Vect().Angle(ppi.Vect());
double th_nu_mu = pnu.Vect().Angle(pmu.Vect());
double th_pi_mu = ppi.Vect().Angle(pmu.Vect());
double rEnu = -999;
// If experiment doesn't have pion kinematic information (T2K CC1pi+ CH
// example of this; Michel e sample has no directional information on pion)
// We'll need to modify the reconstruction variables
if (!pionInfo) {
// Assumes that pion angle contribution contributes net zero
// Also assumes the momentum of reconstructed pions when using Michel
// electrons is 130 MeV
// So need to redefine E_pi and p_pi
// It's a little unclear what happens to the pmu . ppi term (4-vector); do
// we include the angular contribution?
// This below doesn't
th_nu_pi = M_PI / 2.;
p_pi = 0.130;
E_pi = sqrt(p_pi * p_pi + m_pi * m_pi);
// rEnu = (m_mu*m_mu + m_pi*m_pi - 2*m_n*(E_mu + E_pi) + 2*E_mu*E_pi -
// 2*p_mu*p_pi) / (2*(E_mu + E_pi - p_mu*cos(th_nu_mu) - m_n));
}
rEnu =
(m_mu * m_mu + m_pi * m_pi - 2 * m_n * (E_pi + E_mu) + 2 * E_pi * E_mu -
2 * p_pi * p_mu * cos(th_pi_mu)) /
(2 * (E_pi + E_mu - p_pi * cos(th_nu_pi) - p_mu * cos(th_nu_mu) - m_n));
return rEnu;
};
//********************************************************************
// Reconstruct neutrino energy from outgoing particles; will differ from the
// actual neutrino energy. Here we use assumption of a Delta resonance
double FitUtils::EnuCC1piprecDelta(TLorentzVector pnu, TLorentzVector pmu) {
//********************************************************************
const double m_Delta =
PhysConst::mass_delta; // PDG value for Delta mass in GeV
const double m_n = PhysConst::mass_neutron; // neutron/proton mass
// should really take proton mass for proton interaction, neutron for neutron
// interaction. However, difference is pretty much negligable here!
double E_mu = pmu.E() / 1000;
double p_mu = pmu.Vect().Mag() / 1000;
double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu);
double th_nu_mu = pnu.Vect().Angle(pmu.Vect());
double rEnu = (m_Delta * m_Delta - m_n * m_n - m_mu * m_mu + 2 * m_n * E_mu) /
(2 * (m_n - E_mu + p_mu * cos(th_nu_mu)));
return rEnu;
};
// MOVE TO T2K UTILS!
//********************************************************************
// Reconstruct Enu using "extended MiniBooNE" as defined in Raquel's T2K TN
//
// Supposedly includes pion direction and binding energy of target nucleon
// I'm not convinced (yet), maybe
double FitUtils::EnuCC1piprec_T2K_eMB(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi) {
//********************************************************************
// Unit vector for neutrino momentum
TVector3 p_nu_vect_unit = pnu.Vect() * (1. / pnu.E());
double E_mu = pmu.E() / 1000.;
TVector3 p_mu_vect = pmu.Vect() * (1. / 1000.);
double E_pi = ppi.E() / 1000.;
TVector3 p_pi_vect = ppi.Vect() * (1. / 1000.);
double E_bind = 25. / 1000.;
double m_p = PhysConst::mass_proton;
// Makes life a little easier, gonna square this one
double a1 = m_p - E_bind - E_mu - E_pi;
// Just gets the magnitude square of the muon and pion momentum vectors
double a2 = (p_mu_vect + p_pi_vect).Mag2();
// Gets the somewhat complicated scalar product between neutrino and
// (p_mu+p_pi), scaled to Enu
double a3 = p_nu_vect_unit * (p_mu_vect + p_pi_vect);
double rEnu =
(m_p * m_p - a1 * a1 + a2) / (2 * (m_p - E_bind - E_mu - E_pi + a3));
return rEnu;
}
//********************************************************************
// Reconstructed Q2 for CC1pi+
//
// enuType describes how the neutrino energy is reconstructed
// 0 uses true neutrino energy; case for MINERvA and MiniBooNE
// 1 uses "extended MiniBooNE" reconstructed (T2K CH)
// 2 uses "MiniBooNE" reconstructed (EnuCC1piprec with pionInfo = true) (T2K CH)
// "MiniBooNE" reconstructed (EnuCC1piprec with pionInfo = false, the
// case for T2K when using Michel tag) (T2K CH)
// 3 uses Delta for reconstruction (T2K CH)
double FitUtils::Q2CC1piprec(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi, int enuType, bool pionInfo) {
//********************************************************************
double E_mu = pmu.E() / 1000.; // energy of lepton in GeV
double p_mu = pmu.Vect().Mag() / 1000.; // momentum of lepton
double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu); // lepton mass
double th_nu_mu = pnu.Vect().Angle(pmu.Vect());
// Different ways of reconstructing the neutrino energy; default is just to
// use the truth (case 0)
double rEnu = -999;
switch (enuType) {
case 0: // True neutrino energy, default; this is the case for when Q2
// defined as Q2 true (MiniBooNE, MINERvA)
rEnu = pnu.E() / 1000.;
break;
case 1: // Extended MiniBooNE reconstructed, as defined by Raquel's CC1pi+
// CH T2K analysis
rEnu = EnuCC1piprec_T2K_eMB(pnu, pmu, ppi);
break;
case 2: // MiniBooNE reconstructed, as defined by MiniBooNE and Raquel's
// CC1pi+ CH T2K analysis and Linda's CC1pi+ H2O
// Can have this with and without pion info, depending on if Michel tag
// used (Raquel)
rEnu = EnuCC1piprec(pnu, pmu, ppi, pionInfo);
break;
case 3:
rEnu = EnuCC1piprecDelta(pnu, pmu);
break;
} // No need for default here since default value of enuType = 0, defined in
// header
double q2 = -m_mu * m_mu + 2. * rEnu * (E_mu - p_mu * cos(th_nu_mu));
return q2;
};
//********************************************************************
// Returns the reconstructed W from a nucleon and an outgoing pion
//
// Could do this a lot more clever (pp + ppi).Mag() would do the job, but this
// would be less instructive
//********************************************************************
double FitUtils::MpPi(TLorentzVector pp, TLorentzVector ppi) {
double E_p = pp.E();
double p_p = pp.Vect().Mag();
double m_p = sqrt(E_p * E_p - p_p * p_p);
double E_pi = ppi.E();
double p_pi = ppi.Vect().Mag();
double m_pi = sqrt(E_pi * E_pi - p_pi * p_pi);
double th_p_pi = pp.Vect().Angle(ppi.Vect());
// fairly easy thing to derive since bubble chambers measure the proton!
double invMass = sqrt(m_p * m_p + m_pi * m_pi + 2 * E_p * E_pi -
2 * p_pi * p_p * cos(th_p_pi));
return invMass;
};
//********************************************************
// Reconstruct the hadronic mass using neutrino and muon
// Could technically do E_nu = EnuCC1pipRec(pnu,pmu,ppi) too, but this wwill be
// reconstructed Enu; so gives reconstructed Wrec which most of the time isn't
// used!
//
// Only MINERvA uses this so far; and the Enu is Enu_true
// If we want W_true need to take initial state nucleon motion into account
// Return value is in MeV!!!
double FitUtils::Wrec(TLorentzVector pnu, TLorentzVector pmu) {
//********************************************************
double E_mu = pmu.E();
double p_mu = pmu.Vect().Mag();
double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu);
double th_nu_mu = pnu.Vect().Angle(pmu.Vect());
// The factor of 1000 is necessary for downstream functions
const double m_p = PhysConst::mass_proton * 1000;
// MINERvA cut on W_exp which is tuned to W_true; so use true Enu from
// generators
double E_nu = pnu.E();
double w_rec = 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)));
return w_rec;
};
//********************************************************
// Reconstruct the true hadronic mass using the initial state and muon
// Could technically do E_nu = EnuCC1pipRec(pnu,pmu,ppi) too, but this wwill be
// reconstructed Enu; so gives reconstructed Wrec which most of the time isn't
// used!
//
// No one seems to use this because it's fairly MC dependent!
// Return value is in MeV!!!
double FitUtils::Wtrue(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector pnuc) {
//********************************************************
// Could simply do the TLorentzVector operators here but this is more
// instructive?
// ... and prone to errors ...
double E_mu = pmu.E();
double p_mu = pmu.Vect().Mag();
double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu);
double th_nu_mu = pnu.Vect().Angle(pmu.Vect());
double E_nuc = pnuc.E();
double p_nuc = pnuc.Vect().Mag();
double m_nuc = sqrt(E_nuc * E_nuc - p_nuc * p_nuc);
double th_nuc_mu = pmu.Vect().Angle(pnuc.Vect());
double th_nu_nuc = pnu.Vect().Angle(pnuc.Vect());
double E_nu = pnu.E();
double w_rec = sqrt(m_nuc * m_nuc + m_mu * m_mu - 2 * E_nu * E_mu +
2 * E_nu * p_mu * cos(th_nu_mu) - 2 * E_nuc * E_mu +
2 * p_nuc * p_mu * cos(th_nuc_mu) + 2 * E_nu * E_nuc -
2 * E_nu * p_nuc * cos(th_nu_nuc));
return w_rec;
};
double FitUtils::SumKE_PartVect(std::vector<FitParticle *> const fps) {
double sum = 0.0;
for (size_t p_it = 0; p_it < fps.size(); ++p_it) {
sum += fps[p_it]->KE();
}
return sum;
}
double FitUtils::SumTE_PartVect(std::vector<FitParticle *> const fps) {
double sum = 0.0;
for (size_t p_it = 0; p_it < fps.size(); ++p_it) {
sum += fps[p_it]->E();
}
return sum;
}
/*
E Recoil
*/
double FitUtils::GetErecoil_TRUE(FitEvent *event) {
// Get total energy of hadronic system.
double Erecoil = 0.0;
for (unsigned int i = 2; i < event->Npart(); i++) {
// Only final state
if (!event->PartInfo(i)->fIsAlive) continue;
if (event->PartInfo(i)->fNEUTStatusCode != 0) continue;
// Skip Lepton
if (abs(event->PartInfo(i)->fPID) == abs(event->PartInfo(0)->fPID) - 1)
continue;
// Add Up KE of protons and TE of everything else
if (event->PartInfo(i)->fPID == 2212 || event->PartInfo(i)->fPID == 2112) {
Erecoil +=
fabs(event->PartInfo(i)->fP.E()) - fabs(event->PartInfo(i)->fP.Mag());
} else {
Erecoil += event->PartInfo(i)->fP.E();
}
}
return Erecoil;
}
double FitUtils::GetErecoil_CHARGED(FitEvent *event) {
// Get total energy of hadronic system.
double Erecoil = 0.0;
for (unsigned int i = 2; i < event->Npart(); i++) {
// Only final state
if (!event->PartInfo(i)->fIsAlive) continue;
if (event->PartInfo(i)->fNEUTStatusCode != 0) continue;
// Skip Lepton
if (abs(event->PartInfo(i)->fPID) == abs(event->PartInfo(0)->fPID) - 1)
continue;
// Skip Neutral particles
if (event->PartInfo(i)->fPID == 2112 || event->PartInfo(i)->fPID == 111 ||
event->PartInfo(i)->fPID == 22)
continue;
// Add Up KE of protons and TE of everything else
if (event->PartInfo(i)->fPID == 2212) {
Erecoil +=
fabs(event->PartInfo(i)->fP.E()) - fabs(event->PartInfo(i)->fP.Mag());
} else {
Erecoil += event->PartInfo(i)->fP.E();
}
}
return Erecoil;
}
// MOVE TO MINERVA Utils!
double FitUtils::GetErecoil_MINERvA_LowRecoil(FitEvent *event) {
// Get total energy of hadronic system.
double Erecoil = 0.0;
for (unsigned int i = 2; i < event->Npart(); i++) {
// Only final state
if (!event->PartInfo(i)->fIsAlive) continue;
if (event->PartInfo(i)->Status() != kFinalState) continue;
// Skip Lepton
if (abs(event->PartInfo(i)->fPID) == 13) continue;
// Skip Neutrons particles
if (event->PartInfo(i)->fPID == 2112) continue;
int PID = event->PartInfo(i)->fPID;
// KE of Protons and charged pions
if (PID == 2212 or PID == 211 or PID == -211) {
// Erecoil += FitUtils::T(event->PartInfo(i)->fP);
Erecoil +=
fabs(event->PartInfo(i)->fP.E()) - fabs(event->PartInfo(i)->fP.Mag());
// Total Energy of non-neutrons
// } else if (PID != 2112 and PID < 999 and PID != 22 and abs(PID) !=
// 14) {
} else if (PID == 111 || PID == 11 || PID == -11 || PID == 22) {
Erecoil += (event->PartInfo(i)->fP.E());
}
}
return Erecoil;
}
// MOVE TO MINERVA Utils!
// The alternative Eavailble definition takes true q0 and subtracts the kinetic energy of neutrons and pion masses
// returns in MeV
double FitUtils::Eavailable(FitEvent *event) {
double Eav = 0.0;
// Now take q0 and subtract Eav
double q0 = event->GetBeamPart()->fP.E();
// Get the pdg of incoming neutrino
int ISPDG = event->GetBeamPartPDG();
// For CC
if (event->IsCC()) q0 -= event->GetHMFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1))->fP.E();
else q0 -= event->GetHMFSParticle(ISPDG)->fP.E();
for (unsigned int i = 2; i < event->Npart(); i++) {
// Only final state
if (!event->PartInfo(i)->fIsAlive) continue;
if (event->PartInfo(i)->Status() != kFinalState) continue;
int PID = event->PartInfo(i)->fPID;
// Neutrons
if (PID == 2112) {
// Adding kinetic energy of neutron
Eav += FitUtils::T(event->PartInfo(i)->fP)*1000.;
// All pion masses
} else if (abs(PID) == 211 || PID == 111) {
Eav += event->PartInfo(i)->fP.M();
}
}
return q0-Eav;
}
TVector3 GetVectorInTPlane(const TVector3 &inp, const TVector3 &planarNormal) {
TVector3 pnUnit = planarNormal.Unit();
double inpProjectPN = inp.Dot(pnUnit);
TVector3 InPlanarInput = inp - (inpProjectPN * pnUnit);
return InPlanarInput;
}
TVector3 GetUnitVectorInTPlane(const TVector3 &inp,
const TVector3 &planarNormal) {
return GetVectorInTPlane(inp, planarNormal).Unit();
}
Double_t GetDeltaPhiT(TVector3 const &V_lepton, TVector3 const &V_other,
TVector3 const &Normal, bool PiMinus = false) {
TVector3 V_lepton_T = GetUnitVectorInTPlane(V_lepton, Normal);
TVector3 V_other_T = GetUnitVectorInTPlane(V_other, Normal);
return PiMinus ? acos(V_lepton_T.Dot(V_other_T))
: (M_PI - acos(V_lepton_T.Dot(V_other_T)));
}
TVector3 GetDeltaPT(TVector3 const &V_lepton, TVector3 const &V_other,
TVector3 const &Normal) {
TVector3 V_lepton_T = GetVectorInTPlane(V_lepton, Normal);
TVector3 V_other_T = GetVectorInTPlane(V_other, Normal);
return V_lepton_T + V_other_T;
}
Double_t GetDeltaAlphaT(TVector3 const &V_lepton, TVector3 const &V_other,
TVector3 const &Normal, bool PiMinus = false) {
TVector3 DeltaPT = GetDeltaPT(V_lepton, V_other, Normal);
return GetDeltaPhiT(V_lepton, DeltaPT, Normal, PiMinus);
}
-double FitUtils::Get_STV_dpt(FitEvent *event, int ISPDG, bool Is0pi) {
+double FitUtils::Get_STV_dpt(FitEvent *event, double ProtonMinCut, double ProtonMaxCut, double ProtonCosThetaCut, int ISPDG, bool Is0pi) {
// Check that the neutrino exists
if (event->NumISParticle(ISPDG) == 0) {
return -9999;
}
// Return 0 if the proton or muon are missing
if (event->NumFSParticle(2212) == 0 ||
event->NumFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1)) == 0) {
return -9999;
}
// Now get the TVector3s for each particle
TVector3 const &NuP = event->GetHMISParticle(ISPDG)->fP.Vect();
TVector3 const &LeptonP =
event->GetHMFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1))->fP.Vect();
- // Find the highest momentum proton in the event between 450 and 1200 MeV with theta_p < 70
- TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
- int HMFSProton = 0;
- double HighestMomentum = 0.0;
- // Get the stack of protons
- std::vector<FitParticle*> Protons = event->GetAllFSProton();
- for (size_t i = 0; i < Protons.size(); ++i) {
- if (Protons[i]->p() > 450 &&
- Protons[i]->p() < 1200 &&
- Protons[i]->P3().Angle(Pnu.Vect()) < (M_PI/180.0)*70.0 &&
- Protons[i]->p() > HighestMomentum) {
- HighestMomentum = Protons[i]->p();
- HMFSProton = i;
- }
- }
- // Now get the proton
- TLorentzVector Pprot = Protons[HMFSProton]->fP;
+ // Find the highest momentum proton in the event between ProtonMinCut and ProtonMaxCut MeV with cos(theta_p) > ProtonCosThetaCut
+ TLorentzVector Pprot = FitUtils::GetProtonInRange(event, ProtonMinCut, ProtonMaxCut, ProtonCosThetaCut);
// Get highest momentum proton in allowed proton range
TVector3 HadronP = Pprot.Vect();
// If we don't have a CC0pi signal definition we also add in pion momentum
if (!Is0pi) {
if (event->NumFSParticle(PhysConst::pdg_pions) == 0) {
return -9999;
}
// Count up pion momentum
TLorentzVector pp = event->GetHMFSParticle(PhysConst::pdg_pions)->fP;
HadronP += pp.Vect();
}
return GetDeltaPT(LeptonP, HadronP, NuP).Mag();
}
-double FitUtils::Get_STV_dphit(FitEvent *event, int ISPDG, bool Is0pi) {
+double FitUtils::Get_STV_dphit(FitEvent *event, double ProtonMinCut, double ProtonMaxCut, double ProtonCosThetaCut, int ISPDG, bool Is0pi) {
// Check that the neutrino exists
if (event->NumISParticle(ISPDG) == 0) {
return -9999;
}
// Return 0 if the proton or muon are missing
if (event->NumFSParticle(2212) == 0 ||
event->NumFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1)) == 0) {
return -9999;
}
// Now get the TVector3s for each particle
TVector3 const &NuP = event->GetHMISParticle(ISPDG)->fP.Vect();
TVector3 const &LeptonP =
event->GetHMFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1))->fP.Vect();
- // Find the highest momentum proton in the event between 450 and 1200 MeV with theta_p < 70
- TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
- int HMFSProton = 0;
- double HighestMomentum = 0.0;
- // Get the stack of protons
- std::vector<FitParticle*> Protons = event->GetAllFSProton();
- for (size_t i = 0; i < Protons.size(); ++i) {
- if (Protons[i]->p() > 450 &&
- Protons[i]->p() < 1200 &&
- Protons[i]->P3().Angle(Pnu.Vect()) < (M_PI/180.0)*70.0 &&
- Protons[i]->p() > HighestMomentum) {
- HighestMomentum = Protons[i]->p();
- HMFSProton = i;
- }
- }
- // Now get the proton
- TLorentzVector Pprot = Protons[HMFSProton]->fP;
- // Get highest momentum proton in allowed proton range
+ // Find the highest momentum proton in the event between ProtonMinCut and ProtonMaxCut MeV with cos(theta_p) > ProtonCosThetaCut
+ TLorentzVector Pprot = FitUtils::GetProtonInRange(event, ProtonMinCut, ProtonMaxCut, ProtonCosThetaCut);
TVector3 HadronP = Pprot.Vect();
if (!Is0pi) {
if (event->NumFSParticle(PhysConst::pdg_pions) == 0) {
return -9999;
}
TLorentzVector pp = event->GetHMFSParticle(PhysConst::pdg_pions)->fP;
HadronP += pp.Vect();
}
return GetDeltaPhiT(LeptonP, HadronP, NuP);
}
-double FitUtils::Get_STV_dalphat(FitEvent *event, int ISPDG, bool Is0pi) {
+double FitUtils::Get_STV_dalphat(FitEvent *event, double ProtonMinCut, double ProtonMaxCut, double ProtonCosThetaCut, int ISPDG, bool Is0pi) {
// Check that the neutrino exists
if (event->NumISParticle(ISPDG) == 0) {
return -9999;
}
// Return 0 if the proton or muon are missing
if (event->NumFSParticle(2212) == 0 ||
event->NumFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1)) == 0) {
return -9999;
}
// Now get the TVector3s for each particle
TVector3 const &NuP = event->GetHMISParticle(ISPDG)->fP.Vect();
TVector3 const &LeptonP =
event->GetHMFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1))->fP.Vect();
- // Find the highest momentum proton in the event between 450 and 1200 MeV with theta_p < 70
- TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
- int HMFSProton = 0;
- double HighestMomentum = 0.0;
- // Get the stack of protons
- std::vector<FitParticle*> Protons = event->GetAllFSProton();
- for (size_t i = 0; i < Protons.size(); ++i) {
- if (Protons[i]->p() > 450 &&
- Protons[i]->p() < 1200 &&
- Protons[i]->P3().Angle(Pnu.Vect()) < (M_PI/180.0)*70.0 &&
- Protons[i]->p() > HighestMomentum) {
- HighestMomentum = Protons[i]->p();
- HMFSProton = i;
- }
- }
- // Now get the proton
- TLorentzVector Pprot = Protons[HMFSProton]->fP;
- // Get highest momentum proton in allowed proton range
+ // Find the highest momentum proton in the event between ProtonMinCut and ProtonMaxCut MeV with cos(theta_p) > ProtonCosThetaCut
+ TLorentzVector Pprot = FitUtils::GetProtonInRange(event, ProtonMinCut, ProtonMaxCut, ProtonCosThetaCut);
TVector3 HadronP = Pprot.Vect();
if (!Is0pi) {
if (event->NumFSParticle(PhysConst::pdg_pions) == 0) {
return -9999;
}
TLorentzVector pp = event->GetHMFSParticle(PhysConst::pdg_pions)->fP;
HadronP += pp.Vect();
}
return GetDeltaAlphaT(LeptonP, HadronP, NuP);
}
// As defined in PhysRevC.95.065501
// Using prescription from arXiv 1805.05486
// Returns in GeV
-double FitUtils::Get_pn_reco_C(FitEvent *event, int ISPDG, bool Is0pi) {
+double FitUtils::Get_pn_reco_C(FitEvent *event, double ProtonMinCut, double ProtonMaxCut, double ProtonCosThetaCut, int ISPDG, bool Is0pi) {
const double mn = PhysConst::mass_neutron; // neutron mass
const double mp = PhysConst::mass_proton; // proton mass
// Check that the neutrino exists
if (event->NumISParticle(ISPDG) == 0) {
return -9999;
}
// Return 0 if the proton or muon are missing
if (event->NumFSParticle(2212) == 0 ||
event->NumFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1)) == 0) {
return -9999;
}
// Now get the TVector3s for each particle
TVector3 const &NuP = event->GetHMISParticle(ISPDG)->fP.Vect();
TVector3 const &LeptonP =
event->GetHMFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1))->fP.Vect();
- // Find the highest momentum proton in the event between 450 and 1200 MeV with theta_p < 70
- TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
- int HMFSProton = 0;
- double HighestMomentum = 0.0;
- // Get the stack of protons
- std::vector<FitParticle*> Protons = event->GetAllFSProton();
- for (size_t i = 0; i < Protons.size(); ++i) {
- // Update the highest momentum particle
- if (Protons[i]->p() > 450 &&
- Protons[i]->p() < 1200 &&
- Protons[i]->P3().Angle(Pnu.Vect()) < (M_PI/180.0)*70.0 &&
- Protons[i]->p() > HighestMomentum) {
- HighestMomentum = Protons[i]->p();
- HMFSProton = i;
- }
- }
- // Now get the proton
- TLorentzVector Pprot = Protons[HMFSProton]->fP;
- // Get highest momentum proton in allowed proton range
+ // Find the highest momentum proton in the event between ProtonMinCut and ProtonMaxCut MeV with cos(theta_p) < ProtonCosThetaCut
+ TLorentzVector Pprot = FitUtils::GetProtonInRange(event, ProtonMinCut, ProtonMaxCut, ProtonCosThetaCut);
TVector3 HadronP = Pprot.Vect();
- //TVector3 HadronP = event->GetHMFSParticle(2212)->fP.Vect();
double const el = event->GetHMFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1))->E()/1000.;
double const eh = Pprot.E()/1000.;
if (!Is0pi) {
if (event->NumFSParticle(PhysConst::pdg_pions) == 0) {
return -9999;
}
TLorentzVector pp = event->GetHMFSParticle(PhysConst::pdg_pions)->fP;
HadronP += pp.Vect();
}
TVector3 dpt = GetDeltaPT(LeptonP, HadronP, NuP);
double dptMag = dpt.Mag()/1000.;
double ma = 6*mn + 6*mp - 0.09216; // target mass (E is from PhysRevC.95.065501)
double map = ma - mn + 0.02713; // remnant mass
double pmul = LeptonP.Dot(NuP.Unit())/1000.;
double phl = HadronP.Dot(NuP.Unit())/1000.;
//double pmul = GetVectorInTPlane(LeptonP, dpt).Mag()/1000.;
//double phl = GetVectorInTPlane(HadronP, dpt).Mag()/1000.;
double R = ma + pmul + phl - el - eh;
double dpl = 0.5*R - (map*map + dptMag*dptMag)/(2*R);
//double dpl = ((R*R)-(dptMag*dptMag)-(map*map))/(2*R); // as in in PhysRevC.95.065501 - gives same result
double pn_reco = sqrt((dptMag*dptMag) + (dpl*dpl));
//std::cout << "Diagnostics: " << std::endl;
//std::cout << "mn: " << mn << std::endl;
//std::cout << "ma: " << ma << std::endl;
//std::cout << "map: " << map << std::endl;
//std::cout << "pmu: " << LeptonP.Mag()/1000. << std::endl;
//std::cout << "ph: " << HadronP.Mag()/1000. << std::endl;
//std::cout << "pmul: " << pmul << std::endl;
//std::cout << "phl: " << phl << std::endl;
//std::cout << "el: " << el << std::endl;
//std::cout << "eh: " << eh << std::endl;
//std::cout << "R: " << R << std::endl;
//std::cout << "dptMag: " << dptMag << std::endl;
//std::cout << "dpl: " << dpl << std::endl;
//std::cout << "pn_reco: " << pn_reco << std::endl;
return pn_reco;
}
-double FitUtils::Get_pn_reco_Ar(FitEvent *event, int ISPDG, bool Is0pi) {
+double FitUtils::Get_pn_reco_Ar(FitEvent *event, double ProtonMinCut, double ProtonMaxCut, double ProtonCosThetaCut, int ISPDG, bool Is0pi) {
const double mn = PhysConst::mass_neutron; // neutron mass
const double mp = PhysConst::mass_proton; // proton mass
// Check that the neutrino exists
if (event->NumISParticle(ISPDG) == 0) {
return -9999;
}
// Return 0 if the proton or muon are missing
if (event->NumFSParticle(2212) == 0 ||
event->NumFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1)) == 0) {
return -9999;
}
// Now get the TVector3s for each particle
TVector3 const &NuP = event->GetHMISParticle(ISPDG)->fP.Vect();
TVector3 const &LeptonP =
event->GetHMFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1))->fP.Vect();
- // Find the highest momentum proton in the event between 450 and 1200 MeV with theta_p < 70
- TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
- int HMFSProton = 0;
- double HighestMomentum = 0.0;
- // Get the stack of protons
- std::vector<FitParticle*> Protons = event->GetAllFSProton();
- for (size_t i = 0; i < Protons.size(); ++i) {
- // Update the highest momentum particle
- if (Protons[i]->p() > 450 &&
- Protons[i]->p() < 1200 &&
- Protons[i]->P3().Angle(Pnu.Vect()) < (M_PI/180.0)*70.0 &&
- Protons[i]->p() > HighestMomentum) {
- HighestMomentum = Protons[i]->p();
- HMFSProton = i;
- }
- }
- // Now get the proton
- TLorentzVector Pprot = Protons[HMFSProton]->fP;
- // Get highest momentum proton in allowed proton range
+ // Find the highest momentum proton in the event between ProtonMinCut and ProtonMaxCut MeV with cos(theta_p) > ProtonCosThetaCut
+ TLorentzVector Pprot = FitUtils::GetProtonInRange(event, ProtonMinCut, ProtonMaxCut, ProtonCosThetaCut);
TVector3 HadronP = Pprot.Vect();
- //TVector3 HadronP = event->GetHMFSParticle(2212)->fP.Vect();
double const el = event->GetHMFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1))->E()/1000.;
double const eh = Pprot.E()/1000.;
if (!Is0pi) {
if (event->NumFSParticle(PhysConst::pdg_pions) == 0) {
return -9999;
}
TLorentzVector pp = event->GetHMFSParticle(PhysConst::pdg_pions)->fP;
HadronP += pp.Vect();
}
TVector3 dpt = GetDeltaPT(LeptonP, HadronP, NuP);
double dptMag = dpt.Mag()/1000.;
//double ma = 6*mn + 6*mp - 0.09216; // target mass (E is from PhysRevC.95.065501)
//double map = ma - mn + 0.02713; // remnant mass
double ma = 6*mn + 6*mp - 0.34381; // target mass (E is from PhysRevC.95.065501)
double map = ma - mn + 0.02713; // remnant mass
double pmul = LeptonP.Dot(NuP.Unit())/1000.;
double phl = HadronP.Dot(NuP.Unit())/1000.;
//double pmul = GetVectorInTPlane(LeptonP, dpt).Mag()/1000.;
//double phl = GetVectorInTPlane(HadronP, dpt).Mag()/1000.;
double R = ma + pmul + phl - el - eh;
double dpl = 0.5*R - (map*map + dptMag*dptMag)/(2*R);
//double dpl = ((R*R)-(dptMag*dptMag)-(map*map))/(2*R); // as in in PhysRevC.95.065501 - gives same result
double pn_reco = sqrt((dptMag*dptMag) + (dpl*dpl));
//std::cout << "Diagnostics: " << std::endl;
//std::cout << "mn: " << mn << std::endl;
//std::cout << "ma: " << ma << std::endl;
//std::cout << "map: " << map << std::endl;
//std::cout << "pmu: " << LeptonP.Mag()/1000. << std::endl;
//std::cout << "ph: " << HadronP.Mag()/1000. << std::endl;
//std::cout << "pmul: " << pmul << std::endl;
//std::cout << "phl: " << phl << std::endl;
//std::cout << "el: " << el << std::endl;
//std::cout << "eh: " << eh << std::endl;
//std::cout << "R: " << R << std::endl;
//std::cout << "dptMag: " << dptMag << std::endl;
//std::cout << "dpl: " << dpl << std::endl;
//std::cout << "pn_reco: " << pn_reco << std::endl;
return pn_reco;
}
+// Find the highest momentum proton in the event between ProtonMinCut and ProtonMaxCut MeV with cos(theta_p) > ProtonCosThetaCut
+TLorentzVector FitUtils::GetProtonInRange(FitEvent *event, double ProtonMinCut, double ProtonMaxCut, double ProtonCosThetaCut) {
+ // Get the neutrino
+ TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
+ int HMFSProton = 0;
+ double HighestMomentum = 0.0;
+ // Get the stack of protons
+ std::vector<FitParticle*> Protons = event->GetAllFSProton();
+ for (size_t i = 0; i < Protons.size(); ++i) {
+ if (Protons[i]->p() > ProtonMinCut &&
+ Protons[i]->p() < ProtonMaxCut &&
+ cos(Protons[i]->P3().Angle(Pnu.Vect())) > ProtonCosThetaCut &&
+ Protons[i]->p() > HighestMomentum) {
+ HighestMomentum = Protons[i]->p();
+ HMFSProton = i;
+ }
+ }
+ // Now get the proton
+ TLorentzVector Pprot = Protons[HMFSProton]->fP;
+ return Pprot;
+}
+
// Get Cos theta with Adler angles
double FitUtils::CosThAdler(TLorentzVector Pnu, TLorentzVector Pmu, TLorentzVector Ppi, TLorentzVector Pprot) {
// Get the "resonance" lorentz vector (pion proton system)
TLorentzVector Pres = Pprot + Ppi;
// Boost the particles into the resonance rest frame so we can define the x,y,z axis
Pnu.Boost(-Pres.BoostVector());
Pmu.Boost(-Pres.BoostVector());
Ppi.Boost(-Pres.BoostVector());
// The z-axis is defined as the axis of three-momentum transfer, \vec{k}
// Also unit normalise the axis
TVector3 zAxis = (Pnu.Vect()-Pmu.Vect())*(1.0/((Pnu.Vect()-Pmu.Vect()).Mag()));
// Then the angle between the pion in the "resonance" rest-frame and the z-axis is the theta Adler angle
double costhAdler = cos(Ppi.Vect().Angle(zAxis));
return costhAdler;
}
// Get phi with Adler angles, a bit more complicated...
double FitUtils::PhiAdler(TLorentzVector Pnu, TLorentzVector Pmu, TLorentzVector Ppi, TLorentzVector Pprot) {
// Get the "resonance" lorentz vector (pion proton system)
TLorentzVector Pres = Pprot + Ppi;
// Boost the particles into the resonance rest frame so we can define the x,y,z axis
Pnu.Boost(-Pres.BoostVector());
Pmu.Boost(-Pres.BoostVector());
Ppi.Boost(-Pres.BoostVector());
// The z-axis is defined as the axis of three-momentum transfer, \vec{k}
// Also unit normalise the axis
TVector3 zAxis = (Pnu.Vect()-Pmu.Vect())*(1.0/((Pnu.Vect()-Pmu.Vect()).Mag()));
// The y-axis is then defined perpendicular to z and muon momentum in the resonance frame
TVector3 yAxis = Pnu.Vect().Cross(Pmu.Vect());
yAxis *= 1.0/double(yAxis.Mag());
// And the x-axis is then simply perpendicular to z and x
TVector3 xAxis = yAxis.Cross(zAxis);
xAxis *= 1.0/double(xAxis.Mag());
// Project the pion on to x and y axes
double x = Ppi.Vect().Dot(xAxis);
double y = Ppi.Vect().Dot(yAxis);
double newphi = atan2(y, x)*(180./M_PI);
// Convert negative angles to positive
if (newphi < 0.0) newphi += 360.0;
return newphi;
}
//********************************************************************
double FitUtils::ppInfK(TLorentzVector pmu, double costh, double binding,
bool neutrino) {
//********************************************************************
// Convert all values to GeV
//const double V = binding / 1000.; // binding potential
//const double mn = PhysConst::mass_neutron; // neutron mass
const double mp = PhysConst::mass_proton; // proton mass
double el = pmu.E() / 1000.;
//double pl = (pmu.Vect().Mag()) / 1000.; // momentum of lepton
double enu = EnuQErec(pmu, costh, binding, neutrino);
double ep_inf = enu - el + mp;
double pp_inf = sqrt(ep_inf * ep_inf - mp * mp);
return pp_inf;
};
//********************************************************************
TVector3 FitUtils::tppInfK(TLorentzVector pmu, double costh, double binding,
bool neutrino) {
//********************************************************************
// Convert all values to GeV
//const double V = binding / 1000.; // binding potential
//const double mn = PhysConst::mass_neutron; // neutron mass
//const double mp = PhysConst::mass_proton; // proton mass
double pl_x = pmu.X() / 1000.;
double pl_y = pmu.Y() / 1000.;
double pl_z= pmu.Z() / 1000.;
double enu = EnuQErec(pmu, costh, binding, neutrino);
TVector3 tpp_inf(-pl_x, -pl_y, -pl_z+enu);
return tpp_inf;
};
//********************************************************************
double FitUtils::cthpInfK(TLorentzVector pmu, double costh, double binding,
bool neutrino) {
//********************************************************************
// Convert all values to GeV
//const double V = binding / 1000.; // binding potential
//const double mn = PhysConst::mass_neutron; // neutron mass
const double mp = PhysConst::mass_proton; // proton mass
double el = pmu.E() / 1000.;
double pl = (pmu.Vect().Mag()) / 1000.; // momentum of lepton
double enu = EnuQErec(pmu, costh, binding, neutrino);
double ep_inf = enu - el + mp;
double pp_inf = sqrt(ep_inf * ep_inf - mp * mp);
double cth_inf = (enu*enu + pp_inf*pp_inf - pl*pl)/(2*enu*pp_inf);
return cth_inf;
};
diff --git a/src/Utils/FitUtils.h b/src/Utils/FitUtils.h
index 6c3cb7f..cc9326c 100644
--- a/src/Utils/FitUtils.h
+++ b/src/Utils/FitUtils.h
@@ -1,190 +1,193 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************/
#ifndef FITUTILS_H_SEEN
#define FITUTILS_H_SEEN
#include <math.h>
#include <stdlib.h>
#include <unistd.h>
#include <ctime>
#include <iostream>
#include <numeric>
#include <TChain.h>
#include <TFile.h>
#include <TH1D.h>
#include <TH2D.h>
#include <THStack.h>
#include <TKey.h>
#include <TLegend.h>
#include <TList.h>
#include <TLorentzVector.h>
#include <TObjArray.h>
#include <TROOT.h>
#include <TRandom3.h>
#include <TTree.h>
#include "FitEvent.h"
#include "TGraph.h"
#include "TH2Poly.h"
#include "FitEvent.h"
#include "FitLogger.h"
/*!
* \addtogroup Utils
* @{
*/
/// Functions needed by individual samples for calculating kinematic quantities.
namespace FitUtils {
/// Return a vector of all values saved in map
double *GetArrayFromMap(std::vector<std::string> invals,
std::map<std::string, double> inmap);
/// Returns kinetic energy of particle
double T(TLorentzVector part);
/// Returns momentum of particle
double p(TLorentzVector part);
double p(FitParticle* part);
/// Returns angle between particles (_NOT_ cosine!)
double th(TLorentzVector part, TLorentzVector part2);
double th(FitParticle* part1, FitParticle* part2);
/// Hadronic mass reconstruction
double Wrec(TLorentzVector pnu, TLorentzVector pmu);
/// Hadronic mass true from initial state particles and muon; useful if the full
/// FSI vectors aren't not saved and we for some reasons need W_true
double Wtrue(TLorentzVector pnu, TLorentzVector pmu, TLorentzVector pnuc);
double SumKE_PartVect(std::vector<FitParticle *> const fps);
double SumTE_PartVect(std::vector<FitParticle *> const fps);
/// Return E Hadronic for all FS Particles in Hadronic System
double GetErecoil_TRUE(FitEvent *event);
/// Return E Hadronic for all Charged FS Particles in Hadronic System
double GetErecoil_CHARGED(FitEvent *event);
double Eavailable(FitEvent *event);
/*
CCQE MiniBooNE/MINERvA
*/
/// Function to calculate the reconstructed Q^{2}_{QE}
double Q2QErec(TLorentzVector pmu, double costh, double binding, bool neutrino = true);
/// Function returns the reconstructed E_{nu} values
double EnuQErec(TLorentzVector pmu, double costh, double binding, bool neutrino = true);
/// Function returns the reconstructed E_{nu} values
double EnuQErec(TLorentzVector pmu, TLorentzVector pnu, double binding, bool neutrino = true);
//! Function to calculate the reconstructed Q^{2}_{QE}
double Q2QErec(double pl, double costh, double binding, bool neutrino = true);
//! Function to calculate the reconstructed Q^{2}_{QE}
double Q2QErec(TLorentzVector Pmu, TLorentzVector Pnu, double binding, bool neutrino = true);
//! Function returns the reconstructed E_{nu} values
double EnuQErec(double pl, double costh, double binding,
bool neutrino = true);
/*
CCQE1p MINERvA
*/
/// Reconstruct Q2QE given just the maximum energy proton.
double ProtonQ2QErec(double pE, double binding);
/*
E Recoil MINERvA
*/
double GetErecoil_MINERvA_LowRecoil(FitEvent *event);
/*
CC1pi0 MiniBooNE
*/
/// Reconstruct Enu from CCpi0 vectors and binding energy
double EnuCC1pi0rec(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi0 = TLorentzVector(0, 0, 0, 0));
/// Reconstruct Q2 from CCpi0 vectors and binding energy
double Q2CC1pi0rec(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi0 = TLorentzVector(0, 0, 0, 0));
/*
CC1pi+ MiniBooNE
*/
/// returns reconstructed Enu a la MiniBooNE CCpi+
/// returns reconstructed Enu a la MiniBooNE CCpi+
// Also for when not having pion info (so when we have a Michel tag in T2K)
double EnuCC1piprec(TLorentzVector pnu, TLorentzVector pmu, TLorentzVector ppip,
bool pionInfo = true);
/// returns reconstructed Enu assumming resonance interaction where intermediate
/// resonance was a Delta
double EnuCC1piprecDelta(TLorentzVector pnu, TLorentzVector pmu);
/// returns reconstructed in a variety of flavours
double Q2CC1piprec(TLorentzVector pnu, TLorentzVector pmu, TLorentzVector ppip,
int enuType = 0, bool pionInfo = true);
/*
T2K CC1pi+ on CH
*/
double thq3pi_CC1pip_T2K(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi);
double q3_CC1pip_T2K(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi);
double WrecCC1pip_T2K_MB(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppip);
double EnuCC1piprec_T2K_eMB(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi);
/*
nucleon single pion
*/
double MpPi(TLorentzVector pp, TLorentzVector ppi);
/// Gets delta p T as defined in Phys.Rev. C94 (2016) no.1, 015503
-double Get_STV_dpt(FitEvent *event, int ISPDG, bool Is0pi);
+double Get_STV_dpt(FitEvent *event, double ProtonMinCut=0, double ProtonMaxCut=100E3, double ProtonCosThetaCut=-1, int ISPDG=14, bool Is0pi=true);
/// Gets delta phi T as defined in Phys.Rev. C94 (2016) no.1, 015503
-double Get_STV_dphit(FitEvent *event, int ISPDG, bool Is0pi);
+double Get_STV_dphit(FitEvent *event, double ProtonMinCut=0, double ProtonMaxCut=100E3, double ProtonCosThetaCut=-1, int ISPDG=14, bool Is0pi=true);
/// Gets delta alpha T as defined in Phys.Rev. C94 (2016) no.1, 015503
-double Get_STV_dalphat(FitEvent *event, int ISPDG, bool Is0pi);
+double Get_STV_dalphat(FitEvent *event, double ProtonMinCut=0, double ProtonMaxCut=100E3, double ProtonCosThetaCut=-1, int ISPDG=14, bool Is0pi=true);
// As defined in PhysRevC.95.065501
// Using prescription from arXiv 1805.05486
-double Get_pn_reco_C(FitEvent *event, int ISPDG, bool Is0pi);
-double Get_pn_reco_Ar(FitEvent *event, int ISPDG, bool Is0pi);
+double Get_pn_reco_C(FitEvent *event, double ProtonMinCut=0, double ProtonMaxCut=100E3, double ProtonCosThetaCut=-1, int ISPDG=14, bool Is0pi=true);
+double Get_pn_reco_Ar(FitEvent *event, double ProtonMinCut=0, double ProtonMaxCut=100E3, double ProtonCosThetaCut=-1, int ISPDG=14, bool Is0pi=true);
+
+// Helper for STV; get a proton from an event within some cut range (different for T2K and MINERvA)
+TLorentzVector GetProtonInRange(FitEvent *event, double ProtonMinCut, double ProtonMaxCut, double ProtonCosThetaCut);
//For T2K inferred kinematics analyis - variables defined as on page 7 of T2K TN287v11 (and now arXiv 1802.05078)
double ppInfK(TLorentzVector pmu, double costh, double binding, bool neutrino);
TVector3 tppInfK(TLorentzVector pmu, double costh, double binding, bool neutrino);
double cthpInfK(TLorentzVector pmu, double costh, double binding, bool neutrino);
double CosThAdler(TLorentzVector Pnu, TLorentzVector Pmu, TLorentzVector Ppi, TLorentzVector Pprot);
double PhiAdler(TLorentzVector Pnu, TLorentzVector Pmu, TLorentzVector Ppi, TLorentzVector Pprot);
}
/*! @} */
#endif
diff --git a/src/Utils/PhysConst.h b/src/Utils/PhysConst.h
index f1841b4..e093ceb 100644
--- a/src/Utils/PhysConst.h
+++ b/src/Utils/PhysConst.h
@@ -1,112 +1,114 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************/
#ifndef PHYSCONST_H_SEEN
#define PHYSCONST_H_SEEN
#include <math.h>
#include <stdlib.h>
#include <unistd.h>
#include <cstring>
#include <fstream>
#include <iostream>
#include <iostream>
#include <limits>
#include <numeric>
#include <sstream>
#include <string>
#include <vector>
#include "FitLogger.h"
/// namespace to contain all physical constants used by NUISANCE
namespace PhysConst {
const double mass_proton = 0.93827203; // Proton mass in GeV
const double mass_neutron = 0.93956536; // Neutron mass in GeV
const double mass_nucleon = (mass_proton + mass_neutron) / 2.;
const double mass_proton_kg = 1.6727E-27; // Proton mass in kg
const double mass_neutron_kg = 1.6750E-27; // Neutron mass in kg
const double mass_nucleon_kg = (mass_proton_kg + mass_neutron_kg) / 2.;
const double mass_delta = 1.232; // Delta mass in GeV
const double mass_muon = 0.10565837; // Muon mass in GeV
const double mass_electron = 0.000510998928; // Electron mass in GeV
const double mass_cpi = 0.13957; // charged pion mass in GeV
const double mass_pi0 = 0.13498; // neutral pion mass in GeV
const double mass_cK = 0.493677; // charged kaon mass in GeV
const double mass_K0 = 0.497611; // neutral kaon mass in GeV
inline double GetMass(int pdg) {
switch (abs(pdg)) {
case 11:
return mass_electron;
case 13:
return mass_muon;
case 111:
return mass_pi0;
case 211:
return mass_cpi;
case 321:
return mass_cK;
case 311:
case 310:
case 130:
return mass_K0;
case 2112:
return mass_neutron;
case 2212:
return mass_proton;
default: {
NUIS_ERR(WRN, "Attempted to get mass for PDG: "
<< pdg << ", but it is not catered for. Please add it to "
"src/Utils/PhysConst.h");
return -1;
}
}
}
const double mass_MeV = 1000; // MeV/GeV
const int pdg_neutrinos[] = {12, -12, 14, -14 /*, 16, -16*/};
const int pdg_muons[] = {13, -13};
const int pdg_leptons[] = {11, -11, 13, -13, 15, -15};
const int pdg_all_leptons[] = {11, -11, 13, -13, 15, -15,
12, -12, 14, -14, 16, -16};
const int pdg_pions[] = {211, -211, 111};
const int pdg_charged_pions[] = {211, -211};
const int pdg_strangemesons[] = {
130, 310, 311, 321, 9000311, 9000321, 10311,
10321, 100311, 100321, 9010311, 9010321, 9020311, 9020321,
313, 323, 10313, 10323, 20313, 20323, 100313,
100323, 9000313, 9000323, 30313, 30323, 315, 325,
9000315, 9000325, 10315, 10325, 20315, 20325, 9010315,
9010325, 9020315, 9020325, 317, 327, 9010317, 9010327};
+const int pdg_nucleons[] = {2112, 2212};
+
// Just *-1 to cover possibility
const int pdg_kplus = 321;
const int pdg_antistrangemesons[] = {
-130, -310, -311, -321, -9000311, -9000321, -10311,
-10321, -100311, -100321, -9010311, -9010321, -9020311, -9020321,
-313, -323, -10313, -10323, -20313, -20323, -100313,
-100323, -9000313, -9000323, -30313, -30323, -315, -325,
-9000315, -9000325, -10315, -10325, -20315, -20325, -9010315,
-9010325, -9020315, -9020325, -317, -327, -9010317, -9010327};
}
/*! @} */
#endif

File Metadata

Mime Type
text/x-diff
Expires
Sun, Feb 23, 2:33 PM (22 h, 28 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4485428
Default Alt Text
(166 KB)

Event Timeline