diff --git a/src/InputHandler/FitEvent.cxx b/src/InputHandler/FitEvent.cxx
index fc302aa..13a4e14 100644
--- a/src/InputHandler/FitEvent.cxx
+++ b/src/InputHandler/FitEvent.cxx
@@ -1,459 +1,459 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is pddrt of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "FitEvent.h"
#include <iostream>
#include "TObjArray.h"
FitEvent::FitEvent() {
fGenInfo = NULL;
kRemoveFSIParticles = true;
kRemoveUndefParticles = true;
AllocateParticleStack(400);
};
void FitEvent::AddGeneratorInfo(GeneratorInfoBase* gen) {
fGenInfo = gen;
gen->AllocateParticleStack(kMaxParticles);
}
void FitEvent::AllocateParticleStack(int stacksize) {
LOG(DEB) << "Allocating particle stack of size: " << stacksize << std::endl;
kMaxParticles = stacksize;
fParticleList = new FitParticle*[kMaxParticles];
fParticleMom = new double*[kMaxParticles];
fParticleState = new UInt_t[kMaxParticles];
fParticlePDG = new int[kMaxParticles];
fOrigParticleMom = new double*[kMaxParticles];
fOrigParticleState = new UInt_t[kMaxParticles];
fOrigParticlePDG = new int[kMaxParticles];
for (size_t i = 0; i < kMaxParticles; i++) {
fParticleList[i] = NULL;
fParticleMom[i] = new double[4];
fOrigParticleMom[i] = new double[4];
}
if (fGenInfo) fGenInfo->AllocateParticleStack(kMaxParticles);
}
void FitEvent::ExpandParticleStack(int stacksize) {
DeallocateParticleStack();
AllocateParticleStack(stacksize);
}
void FitEvent::DeallocateParticleStack() {
for (size_t i = 0; i < kMaxParticles; i++) {
if (fParticleList[i]) delete fParticleList[i];
delete fParticleMom[i];
delete fOrigParticleMom[i];
}
delete fParticleMom;
delete fOrigParticleMom;
delete fParticleList;
delete fParticleState;
delete fParticlePDG;
delete fOrigParticleState;
delete fOrigParticlePDG;
if (fGenInfo) fGenInfo -> DeallocateParticleStack();
kMaxParticles = 0;
}
void FitEvent::ClearFitParticles() {
for (size_t i = 0; i < kMaxParticles; i++) {
fParticleList[i] = NULL;
}
}
void FitEvent::FreeFitParticles() {
for (size_t i = 0; i < kMaxParticles; i++) {
FitParticle* fp = fParticleList[i];
if (fp) delete fp;
fParticleList[i] = NULL;
}
}
void FitEvent::ResetParticleList() {
for (unsigned int i = 0; i < kMaxParticles; i++) {
FitParticle* fp = fParticleList[i];
if (fp) delete fp;
fParticleList[i] = NULL;
}
}
void FitEvent::HardReset() {
for (unsigned int i = 0; i < kMaxParticles; i++) {
fParticleList[i] = NULL;
}
}
void FitEvent::ResetEvent() {
fMode = 9999;
Mode = 9999;
fEventNo = -1;
fTotCrs = -1.0;
fTargetA = -1;
fTargetZ = -1;
fTargetH = -1;
fBound = false;
fNParticles = 0;
if (fGenInfo) fGenInfo->Reset();
for (unsigned int i = 0; i < kMaxParticles; i++) {
if (fParticleList[i]) delete fParticleList[i];
fParticleList[i] = NULL;
continue;
fParticlePDG[i] = 0;
fParticleState[i] = kUndefinedState;
fParticleMom[i][0] = 0.0;
fParticleMom[i][1] = 0.0;
fParticleMom[i][2] = 0.0;
fParticleMom[i][3] = 0.0;
fOrigParticlePDG[i] = 0;
fOrigParticleState[i] = kUndefinedState;
fOrigParticleMom[i][0] = 0.0;
fOrigParticleMom[i][1] = 0.0;
fOrigParticleMom[i][2] = 0.0;
fOrigParticleMom[i][3] = 0.0;
}
}
void FitEvent::OrderStack() {
// Copy current stack
int npart = fNParticles;
for (int i = 0; i < npart; i++) {
fOrigParticlePDG[i] = fParticlePDG[i];
fOrigParticleState[i] = fParticleState[i];
fOrigParticleMom[i][0] = fParticleMom[i][0];
fOrigParticleMom[i][1] = fParticleMom[i][1];
fOrigParticleMom[i][2] = fParticleMom[i][2];
fOrigParticleMom[i][3] = fParticleMom[i][3];
}
// Now run loops for each particle
fNParticles = 0;
int stateorder[6] = {kInitialState, kFinalState, kFSIState,
kNuclearInitial, kNuclearRemnant, kUndefinedState
};
for (int s = 0; s < 6; s++) {
for (int i = 0; i < npart; i++) {
if ((UInt_t)fOrigParticleState[i] != (UInt_t)stateorder[s]) continue;
fParticlePDG[fNParticles] = fOrigParticlePDG[i];
fParticleState[fNParticles] = fOrigParticleState[i];
fParticleMom[fNParticles][0] = fOrigParticleMom[i][0];
fParticleMom[fNParticles][1] = fOrigParticleMom[i][1];
fParticleMom[fNParticles][2] = fOrigParticleMom[i][2];
fParticleMom[fNParticles][3] = fOrigParticleMom[i][3];
fNParticles++;
}
}
if (LOG_LEVEL(DEB)) {
LOG(DEB) << "Ordered stack" << std::endl;
for (int i = 0; i < fNParticles; i++) {
LOG(DEB) << "Particle " << i << ". " << fParticlePDG[i] << " "
<< fParticleMom[i][0] << " " << fParticleMom[i][1] << " "
<< fParticleMom[i][2] << " " << fParticleMom[i][3] << " "
<< fParticleState[i] << std::endl;
}
}
if (fNParticles != npart) {
ERR(FTL) << "Dropped some particles when ordering the stack!" << std::endl;
}
return;
}
void FitEvent::Print() {
if (LOG_LEVEL(FIT)) {
LOG(FIT) << "FITEvent print" << std::endl;
LOG(FIT) << "Mode: " << fMode << std::endl;
LOG(FIT) << "Particles: " << fNParticles << std::endl;
LOG(FIT) << " -> Particle Stack " << std::endl;
for (int i = 0; i < fNParticles; i++) {
LOG(FIT) << " -> -> " << i << ". " << fParticlePDG[i] << " "
<< fParticleState[i] << " "
<< " Mom(" << fParticleMom[i][0] << ", " << fParticleMom[i][1]
<< ", " << fParticleMom[i][2] << ", " << fParticleMom[i][3] << ")."
<< std::endl;
}
}
return;
}
/* Read/Write own event class */
void FitEvent::SetBranchAddress(TChain* tn) {
tn->SetBranchAddress("Mode", &fMode);
tn->SetBranchAddress("Mode", &Mode);
tn->SetBranchAddress("EventNo", &fEventNo);
tn->SetBranchAddress("TotCrs", &fTotCrs);
tn->SetBranchAddress("TargetA", &fTargetA);
tn->SetBranchAddress("TargetH", &fTargetH);
tn->SetBranchAddress("Bound", &fBound);
tn->SetBranchAddress("RWWeight", &SavedRWWeight);
tn->SetBranchAddress("InputWeight", &InputWeight);
// This has to be setup by the handler now :(
// tn->SetBranchAddress("NParticles", &fNParticles);
// tn->SetBranchAddress("ParticleState", &fParticleState);
// tn->SetBranchAddress("ParticlePDG", &fParticlePDG);
// tn->SetBranchAddress("ParticleMom", &fParticleMom);
}
void FitEvent::AddBranchesToTree(TTree* tn) {
tn->Branch("Mode", &Mode, "Mode/I");
tn->Branch("EventNo", &fEventNo, "EventNo/i");
tn->Branch("TotCrs", &fTotCrs, "TotCrs/D");
tn->Branch("TargetA", &fTargetA, "TargetA/I");
tn->Branch("TargetH", &fTargetH, "TargetH/I");
tn->Branch("Bound", &fBound, "Bound/O");
tn->Branch("RWWeight", &RWWeight, "RWWeight/D");
tn->Branch("InputWeight", &InputWeight, "InputWeight/D");
tn->Branch("NParticles", &fNParticles, "NParticles/I");
tn->Branch("ParticleState", fOrigParticleState, "ParticleState[NParticles]/i");
tn->Branch("ParticlePDG", fOrigParticlePDG, "ParticlePDG[NParticles]/I");
tn->Branch("ParticleMom", fOrigParticleMom, "ParticleMom[NParticles][4]/D");
}
// ------- EVENT ACCESS FUNCTION --------- //
TLorentzVector FitEvent::GetParticleP4 (int index) const {
if (index == -1 or index >= fNParticles) return TLorentzVector();
return TLorentzVector( fParticleMom[index][0],
fParticleMom[index][1],
fParticleMom[index][2],
fParticleMom[index][3] );
}
TVector3 FitEvent::GetParticleP3 (int index) const {
if (index == -1 or index >= fNParticles) return TVector3();
return TVector3( fParticleMom[index][0],
fParticleMom[index][1],
fParticleMom[index][2] );
}
double FitEvent::GetParticleMom (int index) const {
if (index == -1 or index >= fNParticles) return 0.0;
return sqrt(fParticleMom[index][0] * fParticleMom[index][0] +
fParticleMom[index][1] * fParticleMom[index][1] +
fParticleMom[index][2] * fParticleMom[index][2]);
}
double FitEvent::GetParticleMom2 (int index) const {
if (index == -1 or index >= fNParticles) return 0.0;
return fabs((fParticleMom[index][0] * fParticleMom[index][0] +
fParticleMom[index][1] * fParticleMom[index][1] +
fParticleMom[index][2] * fParticleMom[index][2]));
}
double FitEvent::GetParticleE (int index) const {
if (index == -1 or index >= fNParticles) return 0.0;
return fParticleMom[index][3];
}
int FitEvent::GetParticleState (int index) const {
if (index == -1 or index >= fNParticles) return kUndefinedState;
return (fParticleState[index]);
}
int FitEvent::GetParticlePDG (int index) const {
if (index == -1 or index >= fNParticles) return 0;
return (fParticlePDG[index]);
}
FitParticle* FitEvent::GetParticle (int const i) {
// Check Valid Index
if (i == -1){
return NULL;
}
// Check Valid
if (i > fNParticles) {
ERR(FTL) << "Requesting particle beyond stack!" << std::endl
<< "i = " << i << " N = " << fNParticles << std::endl
<< "Mode = " << fMode << std::endl;
throw;
}
if (!fParticleList[i]) {
/*
std::cout << "Creating particle with values i " << i << " ";
std::cout << fParticleMom[i][0] << " " << fParticleMom[i][1] << " " << fParticleMom[i][2] << " " << fParticleMom[i][3] << " ";
std::cout << fParticlePDG[i] << " " << fParticleState[i] << std::endl;
*/
fParticleList[i] = new FitParticle(fParticleMom[i][0], fParticleMom[i][1],
fParticleMom[i][2], fParticleMom[i][3],
fParticlePDG[i], fParticleState[i]);
} else {
/*
std::cout << "Filling particle with values i " << i << " ";
std::cout << fParticleMom[i][0] << " " << fParticleMom[i][1] << " " << fParticleMom[i][2] << " " << fParticleMom[i][3] << " ";
std::cout << fParticlePDG[i] << " "<< fParticleState[i] <<std::endl;
*/
- // fParticleList[i]->SetValues(fParticleMom[i][0], fParticleMom[i][1],
- // fParticleMom[i][2], fParticleMom[i][3],
- // fParticlePDG[i], fParticleState[i]);
+ fParticleList[i]->SetValues(fParticleMom[i][0], fParticleMom[i][1],
+ fParticleMom[i][2], fParticleMom[i][3],
+ fParticlePDG[i], fParticleState[i]);
}
return fParticleList[i];
}
bool FitEvent::HasParticle(int const pdg, int const state) const {
bool found = false;
for (int i = 0; i < fNParticles; i++) {
if (state != -1 && fParticleState[i] != (uint)state) continue;
if (fParticlePDG[i] == pdg) found = true;
}
return found;
}
int FitEvent::NumParticle(int const pdg, int const state) const {
int nfound = 0;
for (int i = 0; i < fNParticles; i++) {
if (state != -1 and fParticleState[i] != (uint)state) continue;
if (pdg == 0 or fParticlePDG[i] == pdg) nfound += 1;
}
return nfound;
}
std::vector<int> FitEvent::GetAllParticleIndices (int const pdg, int const state) const {
std::vector<int> indexlist;
for (int i = 0; i < fNParticles; i++) {
if (state != -1 and fParticleState[i] != (uint)state) continue;
if (pdg == 0 or fParticlePDG[i] == pdg) {
indexlist.push_back(i);
}
}
return indexlist;
}
std::vector<FitParticle*> FitEvent::GetAllParticle(int const pdg, int const state) {
std::vector<int> indexlist = GetAllParticleIndices(pdg, state);
std::vector<FitParticle*> plist;
for (std::vector<int>::iterator iter = indexlist.begin();
iter != indexlist.end(); iter++) {
plist.push_back( GetParticle((*iter)) );
}
return plist;
}
int FitEvent::GetHMParticleIndex (int const pdg, int const state) const {
double maxmom2 = -9999999.9;
int maxind = -1;
for (int i = 0; i < fNParticles; i++) {
if (state != -1 and fParticleState[i] != (uint)state) continue;
if (pdg == 0 or fParticlePDG[i] == pdg) {
double newmom2 = GetParticleMom2(i);
if (newmom2 > maxmom2) {
maxind = i;
maxmom2 = newmom2;
}
}
}
return maxind;
}
int FitEvent::GetBeamNeutrinoIndex (void) const {
for (int i = 0; i < fNParticles; i++){
if (fParticleState[i] != kInitialState) continue;
int pdg = abs(fParticlePDG[i]);
if (pdg == 12 or pdg == 14 or pdg == 16){
return i;
}
}
return 0;
}
int FitEvent::GetBeamElectronIndex (void) const {
return GetHMISParticleIndex( 11 );
}
int FitEvent::GetBeamPionIndex (void) const {
return GetHMISParticleIndex( PhysConst::pdg_pions );
}
int FitEvent::NumFSMesons() {
int nMesons = 0;
for (int i = 0; i < fNParticles; i++) {
if (fParticleState[i] != kFinalState) continue;
if (abs(fParticlePDG[i]) >= 111 && abs(fParticlePDG[i]) <= 557)
nMesons += 1;
}
return nMesons;
}
int FitEvent::NumFSLeptons(void) const{
int nLeptons = 0;
for (int i = 0; i < fNParticles; i++) {
if (fParticleState[i] != kFinalState) continue;
if (abs(fParticlePDG[i]) == 11 || abs(fParticlePDG[i]) == 13 ||
abs(fParticlePDG[i]) == 15)
nLeptons += 1;
}
return nLeptons;
}
diff --git a/src/InputHandler/FitEvent.h b/src/InputHandler/FitEvent.h
index 50bd5be..dc20304 100644
--- a/src/InputHandler/FitEvent.h
+++ b/src/InputHandler/FitEvent.h
@@ -1,615 +1,640 @@
// 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 "FitParameters.h"
#include "BaseFitEvt.h"
#include "FitLogger.h"
#include "TArrayD.h"
#include "TTree.h"
#include "TChain.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 fMode; };
/// Return Target Atomic Number
inline int GetTargetA (void) const { return fTargetA; };
/// Return Target Nuclear Charge
inline int GetTargetZ (void) const { return fTargetZ; };
/// Get Event Total Cross-section
inline int GetTotCrs (void) const { return fTotCrs; };
/// Is Event Charged Current?
inline bool IsCC() const { return (abs(fMode) <= 30); };
/// Is Event Neutral Current?
inline bool IsNC() const { return (abs(fMode) > 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 int NumFSNuElectron (void) const { return NumFSParticle(12); };
inline int NumFSNuMuon (void) const { return NumFSParticle(14); };
inline int NumFSNuTau (void) const { return NumFSParticle(16); };
inline int NumFSElectron (void) const { return NumFSParticle(11); };
inline int NumFSMuon (void) const { return NumFSParticle(13); };
inline int NumFSTau (void) const { return NumFSParticle(15); };
inline int NumFSProton (void) const { return NumFSParticle(2212); };
inline int NumFSNeutron (void) const { return NumFSParticle(2112); };
inline int NumFSPiZero (void) const { return NumFSParticle(111); };
inline int NumFSPiPlus (void) const { return NumFSParticle(211); };
inline int NumFSPiMinus (void) const { return NumFSParticle(-211); };
inline int NumFSPhoton (void) const { return NumFSParticle(22); };
int NumFSLeptons (void) const; // { return NumFSParticle(PhysConst::pdg_leptons); };
inline int NumFSPions (void) const { return NumFSParticle(PhysConst::pdg_pions); };
inline int NumFSChargePions (void) const { return NumFSParticle(PhysConst::pdg_charged_pions); };
inline std::vector<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<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 FitParticle* GetHMFSNuElectron (void) { return GetHMFSParticle(12); };
inline FitParticle* GetHMFSNuMuon (void) { return GetHMFSParticle(14); };
inline FitParticle* GetHMFSNuTau (void) { return GetHMFSParticle(16); };
inline FitParticle* GetHMFSElectron (void) { return GetHMFSParticle(11); };
inline FitParticle* GetHMFSMuon (void) { return GetHMFSParticle(13); };
inline FitParticle* GetHMFSTau (void) { return GetHMFSParticle(15); };
inline FitParticle* GetHMFSAnyLeptons (void) { return GetHMFSParticle(PhysConst::pdg_all_leptons); };
inline FitParticle* GetHMFSProton (void) { return GetHMFSParticle(2212); };
inline FitParticle* GetHMFSNeutron (void) { return GetHMFSParticle(2112); };
inline FitParticle* GetHMFSPiZero (void) { return GetHMFSParticle(111); };
inline FitParticle* GetHMFSPiPlus (void) { return GetHMFSParticle(211); };
inline FitParticle* GetHMFSPiMinus (void) { return GetHMFSParticle(-211); };
inline FitParticle* GetHMFSPhoton (void) { return GetHMFSParticle(22); };
inline FitParticle* GetHMFSLeptons (void) { return GetHMFSParticle(PhysConst::pdg_leptons); };
inline FitParticle* GetHMFSAnyLepton (void) { return GetHMFSParticle(PhysConst::pdg_all_leptons); };
inline FitParticle* GetHMFSPions (void) { return GetHMFSParticle(PhysConst::pdg_pions); };
inline FitParticle* GetHMFSChargePions(void) { return GetHMFSParticle(PhysConst::pdg_charged_pions); };
inline int GetHMFSNuElectronIndex (void) const { return GetHMFSParticleIndex(12); };
inline int GetHMFSNuMuonIndex (void) const { return GetHMFSParticleIndex(14); };
inline int GetHMFSNuTauIndex (void) const { return GetHMFSParticleIndex(16); };
inline int GetHMFSElectronIndex (void) const { return GetHMFSParticleIndex(11); };
inline int GetHMFSMuonIndex (void) const { return GetHMFSParticleIndex(13); };
inline int GetHMFSTauIndex (void) const { return GetHMFSParticleIndex(15); };
inline int GetHMFSProtonIndex (void) const { return GetHMFSParticleIndex(2212); };
inline int GetHMFSNeutronIndex (void) const { return GetHMFSParticleIndex(2112); };
inline int GetHMFSPiZeroIndex (void) const { return GetHMFSParticleIndex(111); };
inline int GetHMFSPiPlusIndex (void) const { return GetHMFSParticleIndex(211); };
inline int GetHMFSPiMinusIndex (void) const { return GetHMFSParticleIndex(-211); };
inline int GetHMFSPhotonIndex (void) const { return GetHMFSParticleIndex(22); };
inline int GetHMFSLeptonsIndex (void) const { return GetHMFSParticleIndex(PhysConst::pdg_leptons); };
inline int GetHMFSAnyLeptonIndex (void) const { return GetHMFSParticleIndex(PhysConst::pdg_all_leptons); };
inline int GetHMFSPionsIndex (void) const { return GetHMFSParticleIndex(PhysConst::pdg_pions); };
inline int GetHMFSChargePionsIndex(void) const { return GetHMFSParticleIndex(PhysConst::pdg_charged_pions); };
// ---- NEUTRINO INCOMING Related Functions
int GetBeamNeutrinoIndex (void) const;
inline TLorentzVector GetBeamNeutrinoP4 (void) const { return GetParticleP4(GetBeamNeutrinoIndex()); };
inline TVector3 GetBeamNeutrinoP3 (void) const { return GetParticleP3(GetBeamNeutrinoIndex()); };
inline double GetBeamNeutrinoMom (void) const { return GetParticleMom(GetBeamNeutrinoIndex()); };
inline double GetBeamNeutrinoMom2 (void) const { return GetParticleMom2(GetBeamNeutrinoIndex()); };
inline double GetBeamNeutrinoE (void) const { return GetParticleE(GetBeamNeutrinoIndex()); };
inline double Enu (void) const { return GetBeamNeutrinoE(); };
inline int GetBeamNeutrinoPDG (void) const { return GetParticlePDG(GetBeamNeutrinoIndex()); };
inline int PDGnu (void) const { return GetBeamNeutrinoPDG(); };
inline int GetNeutrinoInPos (void) const { return GetBeamNeutrinoIndex(); };
inline FitParticle* GetBeamNeutrino (void) { return GetParticle(GetBeamNeutrinoIndex()); };
inline FitParticle* GetNeutrinoIn (void) { return GetParticle(GetBeamNeutrinoIndex()); };
// ---- Electron INCOMING Related Functions
int GetBeamElectronIndex (void) const;
inline TLorentzVector GetBeamElectronP4 (void) const { return GetParticleP4(GetBeamElectronIndex()); };
inline TVector3 GetBeamElectronP3 (void) const { return GetParticleP3(GetBeamElectronIndex()); };
inline double GetBeamElectronMom (void) const { return GetParticleMom(GetBeamElectronIndex()); };
inline double GetBeamElectronMom2 (void) const { return GetParticleMom2(GetBeamElectronIndex()); };
inline double GetBeamElectronE (void) const { return GetParticleE(GetBeamElectronIndex()); };
inline FitParticle* GetBeamElectron (void) { return GetParticle(GetBeamElectronIndex()); };
// ---- Pion INCOMING Functions
int GetBeamPionIndex (void) const;
inline TLorentzVector GetBeamPionP4 (void) const { return GetParticleP4(GetBeamPionIndex()); };
inline TVector3 GetBeamPionP3 (void) const { return GetParticleP3(GetBeamPionIndex()); };
inline double GetBeamPionMom (void) const { return GetParticleMom(GetBeamPionIndex()); };
inline double GetBeamPionMom2 (void) const { return GetParticleMom2(GetBeamPionIndex()); };
inline double GetBeamPionE (void) const { return GetParticleE(GetBeamPionIndex()); };
inline FitParticle* GetBeamPion (void) { return GetParticle(GetBeamPionIndex()); };
/// Legacy Functions
inline FitParticle* PartInfo(uint i) { return GetParticle(i); };
inline UInt_t Npart (void) const { return NPart(); };
inline UInt_t NPart (void) const { return fNParticles; };
// Other Functions
int NumFSMesons();
inline int GetBeamPartPos(void) const { return 0; };
FitParticle* GetBeamPart(void);
int GetLeptonOutPos(void) const;
FitParticle* GetLeptonOut(void);
// Event Information
double weight; // need for silly reason
int Mode; // Public access needed
int fMode;
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;
double** fOrigParticleMom;
UInt_t* fOrigParticleState;
int* fOrigParticlePDG;
double* fNEUT_ParticleStatusCode;
double* fNEUT_ParticleAliveCode;
GeneratorInfoBase* fGenInfo;
// Config Options for this class
bool kRemoveFSIParticles;
bool kRemoveUndefParticles;
};
/*! @} */
#endif
diff --git a/src/InputHandler/FitParticle.h b/src/InputHandler/FitParticle.h
index ab1bb4f..fa9c263 100644
--- a/src/InputHandler/FitParticle.h
+++ b/src/InputHandler/FitParticle.h
@@ -1,127 +1,109 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************/
#ifndef FITPARTICLE_H_SEEN
#define FITPARTICLE_H_SEEN
/*!
* \addtogroup InputHandler
* @{
*/
#include "TLorentzVector.h"
/// Partle state flags for its position in the event
enum particle_state{
kUndefinedState = 5,
kInitialState = 0,
kFSIState = 1,
kFinalState = 2,
kNuclearInitial = 3,
kNuclearRemnant = 4
};
/// Condensed FitParticle class which acts a common format between the generators
class FitParticle {
public:
/// Create particle of given pdg from momentum variables and state
FitParticle(double x, double y, double z, double t, int pdg, Int_t state);
/// Create empty particle (zero momentum)
FitParticle(){};
~FitParticle(){};
/// Used to change values after creation
void SetValues(double x, double y, double z, double t, int pdg, Int_t state);
/// Return Status Code according to particle_state enum
inline int Status (void) const { return fStatus; };
/// Get Particle PDF
inline int PDG (void) const { return fPID; };
/// Check if particle makes it to final state
inline bool IsFinalState (void) const { return (fStatus == kFinalState); };
/// Was particle involved in intermediate state
inline bool IsFSIState (void) const { return (fStatus == kFSIState); };
/// Was particle incoming
inline bool IsInitialState (void) const { return (fStatus == kInitialState); };
/// Get Mass
inline double M (void){ return fP.Mag(); };
/// Get Kinetic Energy
inline double KE (void){ return fP.E() - fP.Mag(); };
- /// Allows the removal of KE up to total KE.
- inline void RemoveKE(double KE){
- double mass = M();
- double oKE = this->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 = P3().Unit()*n3Mom;
- fP.SetXYZT(np3[0],np3[1],np3[2],nE);
- }
-
-/// Allows the removal of KE up to total KE.
- inline void GiveKE(double KE){
- RemoveKE(-KE);
- }
-
/// Get Total Energy
inline double E (void){ return fP.E(); };
/// Get 4 Momentum
inline TLorentzVector P4(void) {return fP;};
/// Get 3 Momentum
inline TVector3 P3(void) {return fP.Vect();};
/// Get 3 momentum magnitude
inline double p(void) { return fP.Vect().Mag(); };
/// Get 3 momentum magnitude squared
inline double p2(void) { return fP.Vect().Mag2(); };
/// Data Members
TLorentzVector fP; ///< Particle 4 Momentum
int fPID; ///< Particle PDG Code
int fIsAlive; ///< Whether the particle is alive at the end of the event (Yes 1, No 0, Other? -1)
int fNEUTStatusCode; ///< Particle Status (Incoming 1, FSI 2, Outgoing 0, Other 3)
double fMass; ///< Particle Mass
int fStatus; ///< State corresponding to particle_state enum
};
inline std::ostream& operator<<(std::ostream& os, FitParticle const& p){
return os << " Particle[pdgc:" << p.fPID
<< ", stat:"<<p.fStatus
<< ", 4mom:("<< p.fP.X() << "," << p.fP.Y() << "," << p.fP.Z() << "," << p.fP.T() << ")]";
}
/*! @} */
#endif
diff --git a/src/MCStudies/Smear_SVDUnfold_Propagation_Osc.cxx b/src/MCStudies/Smear_SVDUnfold_Propagation_Osc.cxx
index 8fa021c..88d2530 100644
--- a/src/MCStudies/Smear_SVDUnfold_Propagation_Osc.cxx
+++ b/src/MCStudies/Smear_SVDUnfold_Propagation_Osc.cxx
@@ -1,252 +1,475 @@
// 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 "Smear_SVDUnfold_Propagation_Osc.h"
#include "SmearceptanceUtils.h"
#include "OscWeightEngine.h"
//********************************************************************
Smear_SVDUnfold_Propagation_Osc::Smear_SVDUnfold_Propagation_Osc(
nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip =
"Simple measurement class for doing fake data oscillation studies.\n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetTitle("Osc Studies");
fSettings.SetDescription(descrip);
fSettings.SetXTitle("XXX");
fSettings.SetYTitle("Number of events");
fSettings.SetEnuRange(0.0, 50);
fSettings.SetAllowedTypes("EVT/SHAPE/DIAG", "EVT/SHAPE/DIAG");
fSettings.DefineAllowedTargets("*");
fSettings.DefineAllowedSpecies("numu");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = GetEventHistogram()->Integral("width") / (fNEvents + 0.);
// Plot Setup -------------------------------------------------------
// Check that we have the relevant histograms specified.
if (!samplekey.Has("NDDataHist") || !samplekey.Has("FDDataHist") ||
!samplekey.Has("NDToSpectrumSmearingMatrix")) {
THROW(
"Expected to attributes named: NDDataHist, FDDataHist, "
"NDToSpectrumSmearingMatrix on the Smear_SVDUnfold_Propagation_Osc "
"sample "
"tag.");
}
NDDataHist = NULL;
FDDataHist = NULL;
NDToSpectrumSmearingMatrix = NULL;
fMCHist = NULL;
fDataHist = NULL;
std::vector<std::string> NDHistDescriptor =
GeneralUtils::ParseToStr(samplekey.GetS("NDDataHist"), ",");
if (NDHistDescriptor.size() == 2) {
NDDataHist = PlotUtils::GetTH1DFromRootFile(NDHistDescriptor[0],
NDHistDescriptor[1]);
+ ND_True_Spectrum_Hist =
+ PlotUtils::GetTH1DFromRootFile(NDHistDescriptor[0], "ELep_rate");
}
if (!NDDataHist) {
THROW("Attempted to load NDDataHist from the descriptor: \""
<< samplekey.GetS("NDDataHist")
<< "\", but failed. Does it look correct? \"<ROOTFILE>,<HISTNAME>\"");
}
std::vector<std::string> FDHistDescriptor =
GeneralUtils::ParseToStr(samplekey.GetS("FDDataHist"), ",");
if (FDHistDescriptor.size() == 2) {
FDDataHist = PlotUtils::GetTH1DFromRootFile(FDHistDescriptor[0],
FDHistDescriptor[1]);
+
+ FD_True_Spectrum_Hist =
+ PlotUtils::GetTH1DFromRootFile(FDHistDescriptor[0], "ELep_rate");
}
if (!FDDataHist) {
THROW("Attempted to load FDDataHist from the descriptor: \""
<< samplekey.GetS("FDDataHist")
<< "\", but failed. Does it look correct? \"<ROOTFILE>,<HISTNAME>\"");
}
std::vector<std::string> NDToEvSmearingDescriptor = GeneralUtils::ParseToStr(
samplekey.GetS("NDToSpectrumSmearingMatrix"), ",");
if (NDToEvSmearingDescriptor.size() == 2) {
NDToSpectrumSmearingMatrix = PlotUtils::GetTH2DFromRootFile(
NDToEvSmearingDescriptor[0], NDToEvSmearingDescriptor[1]);
}
if (!NDToSpectrumSmearingMatrix) {
THROW("Attempted to load NDToSpectrumSmearingMatrix from the descriptor: \""
<< samplekey.GetS("NDToSpectrumSmearingMatrix")
<< "\", but failed. Does it look correct? \"<ROOTFILE>,<HISTNAME>\"");
}
TMatrixD NDToSpectrumResponseMatrix_l = SmearceptanceUtils::GetMatrix(
SmearceptanceUtils::SVDGetInverse(NDToSpectrumSmearingMatrix));
NDToSpectrumResponseMatrix.ResizeTo(NDToSpectrumResponseMatrix_l);
NDToSpectrumResponseMatrix = NDToSpectrumResponseMatrix_l;
TMatrixD SpectrumToFDResponseMatrix_l =
SmearceptanceUtils::GetMatrix(NDToSpectrumSmearingMatrix);
SpectrumToFDResponseMatrix.ResizeTo(SpectrumToFDResponseMatrix_l);
SpectrumToFDResponseMatrix = SpectrumToFDResponseMatrix_l;
- fDataHist = FDDataHist;
+ fDataHist = static_cast<TH1D *>(FDDataHist->Clone());
fDataHist->SetNameTitle((fSettings.GetName() + "_data").c_str(),
(fSettings.GetFullTitles()).c_str());
SetCovarFromDiagonal();
TruncateUpTo = 0;
if (samplekey.Has("TruncateUpTo")) {
TruncateUpTo = samplekey.GetI("TruncateUpTo");
QLOG(SAM, "Allowed to truncate unfolding matrix by up to "
<< TruncateUpTo
<< " singular values to limit negative ENu spectra.");
}
+ if (samplekey.Has("FitRegion_Min")) {
+ FitRegion_Min = samplekey.GetD("FitRegion_Min");
+ } else {
+ FitRegion_Min = 0xdeadbeef;
+ }
+ if (samplekey.Has("FitRegion_Max")) {
+ FitRegion_Max = samplekey.GetD("FitRegion_Max");
+ } else {
+ FitRegion_Max = 0xdeadbeef;
+ }
+ TruncateStart = 0;
+ if (Config::Get().GetConfigNode("smear.SVD.truncation")) {
+ TruncateStart = Config::Get().ConfI("smear.SVD.truncation");
+
+ TMatrixD NDToSpectrumResponseMatrix_l =
+ SmearceptanceUtils::GetMatrix(SmearceptanceUtils::SVDGetInverse(
+ NDToSpectrumSmearingMatrix, TruncateStart));
+
+ NDToSpectrumResponseMatrix.ResizeTo(NDToSpectrumResponseMatrix_l);
+ NDToSpectrumResponseMatrix = NDToSpectrumResponseMatrix_l;
+ }
+
+ if (TruncateStart >= TruncateUpTo) {
+ TruncateUpTo = TruncateStart + 1;
+ }
+
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
void Smear_SVDUnfold_Propagation_Osc::FillEventVariables(FitEvent *event){};
bool Smear_SVDUnfold_Propagation_Osc::isSignal(FitEvent *event) {
return false;
}
+void Smear_SVDUnfold_Propagation_Osc::UnfoldToNDETrueSpectrum(void) {
+ ND_Unfolded_Spectrum_Hist = NDToSpectrumSmearingMatrix->ProjectionY();
+ ND_Unfolded_Spectrum_Hist->Clear();
+
+ TRandom3 rnjesus;
+
+ bool HasNegValue = false;
+ int truncations = TruncateStart;
+ do {
+ if (truncations >= TruncateUpTo) {
+ THROW("Unfolded enu spectrum had negative values even after "
+ << truncations << " SVD singular value truncations.");
+ }
+
+ // Unfold ND ERec -> Enu spectrum
+ // ------------------------------
+ SmearceptanceUtils::PushTH1ThroughMatrixWithErrors(
+ NDDataHist, ND_Unfolded_Spectrum_Hist, NDToSpectrumResponseMatrix, 1000,
+ false);
+
+ HasNegValue = false;
+
+ for (Int_t bi_it = 1;
+ bi_it < ND_Unfolded_Spectrum_Hist->GetXaxis()->GetNbins() + 1;
+ ++bi_it) {
+ if ((FitRegion_Min != 0xdeadbeef) &&
+ (ND_Unfolded_Spectrum_Hist->GetXaxis()->GetBinUpEdge(bi_it) <=
+ FitRegion_Min)) {
+ continue;
+ }
+ if ((FitRegion_Max != 0xdeadbeef) &&
+ (ND_Unfolded_Spectrum_Hist->GetXaxis()->GetBinLowEdge(bi_it) >
+ FitRegion_Max)) {
+ continue;
+ }
+
+ if (ND_Unfolded_Spectrum_Hist->GetBinContent(bi_it) < 0) {
+ HasNegValue = true;
+ break;
+ }
+ }
+
+ if (HasNegValue) {
+ TMatrixD NDToSpectrumResponseMatrix_l =
+ SmearceptanceUtils::GetMatrix(SmearceptanceUtils::SVDGetInverse(
+ NDToSpectrumSmearingMatrix, truncations));
+
+ NDToSpectrumResponseMatrix.ResizeTo(NDToSpectrumResponseMatrix_l);
+ NDToSpectrumResponseMatrix = NDToSpectrumResponseMatrix_l;
+ }
+
+ truncations++;
+ } while (HasNegValue);
+}
+
void Smear_SVDUnfold_Propagation_Osc::ConvertEventRates(void) {
+ // Apply FD/ND weights
+ // ------------------------------
+ OscHist->Write("ND_FD_Corrected_ENuTrue", TObject::kOverwrite);
+
+ // Apply Oscillations
+ // ------------------------------
+ FitWeight *fw = FitBase::GetRW();
+ OscWeightEngine *oscWE =
+ dynamic_cast<OscWeightEngine *>(fw->GetRWEngine(kOSCILLATION));
+
+ if (!oscWE) {
+ THROW(
+ "Couldn't load oscillation weight engine for sample: "
+ "Smear_SVDUnfold_Propagation_Osc.");
+ }
+
+ for (Int_t bi_it = 1; bi_it < OscHist->GetXaxis()->GetNbins() + 1; ++bi_it) {
+ double oscWeight =
+ oscWE->CalcWeight(OscHist->GetXaxis()->GetBinCenter(bi_it), 14);
+ OscHist->SetBinContent(bi_it, OscHist->GetBinContent(bi_it) * oscWeight);
+ }
+ OscHist->Write("ENuTrue_FD", TObject::kOverwrite);
+
+ TGraph POsc;
+
+ POsc.Set(1E4 - 1);
+
+ double min = OscHist->GetXaxis()->GetBinLowEdge(1);
+ double step =
+ (OscHist->GetXaxis()->GetBinUpEdge(OscHist->GetXaxis()->GetNbins()) -
+ OscHist->GetXaxis()->GetBinLowEdge(1)) /
+ double(1E4);
+
+ for (size_t i = 1; i < 1E4; ++i) {
+ double enu = min + i * step;
+ double ow = oscWE->CalcWeight(enu, 14);
+ if (ow != ow) {
+ std::cout << "Bad osc weight for ENu: " << enu << std::endl;
+ }
+ POsc.SetPoint(i - 1, enu, ow);
+ }
+
+ POsc.Write("POsc", TObject::kOverwrite);
+
+ // Forward fold Spectrum -> ERec FD
+ // ------------------------------
+ if (fMCHist) {
+ fMCHist->SetDirectory(NULL);
+ delete fMCHist;
+ }
+
+ fMCHist = static_cast<TH1D *>(FDDataHist->Clone());
+
+ SmearceptanceUtils::PushTH1ThroughMatrixWithErrors(
+ OscHist, fMCHist, SpectrumToFDResponseMatrix, 1000, true);
+
+ fMCHist->SetBinContent(0, 0);
+ fMCHist->SetBinError(0, 0);
+ fDataHist->SetBinContent(0, 0);
+ fDataHist->SetBinError(0, 0);
+
+ fMCHist->SetBinContent(fMCHist->GetXaxis()->GetNbins() + 1, 0);
+ fMCHist->SetBinError(fMCHist->GetXaxis()->GetNbins() + 1, 0);
+ fDataHist->SetBinContent(fMCHist->GetXaxis()->GetNbins() + 1, 0);
+ fDataHist->SetBinError(fMCHist->GetXaxis()->GetNbins() + 1, 0);
+
+ for (Int_t bi_it = 1; bi_it < fMCHist->GetXaxis()->GetNbins() + 1; ++bi_it) {
+ if ((FitRegion_Min != 0xdeadbeef) &&
+ (fMCHist->GetXaxis()->GetBinUpEdge(bi_it) <= FitRegion_Min)) {
+ fMCHist->SetBinContent(bi_it, 0);
+ fMCHist->SetBinError(bi_it, 0);
+ fDataHist->SetBinContent(bi_it, 0);
+ fDataHist->SetBinError(bi_it, 0);
+ }
+ if ((FitRegion_Max != 0xdeadbeef) &&
+ (fMCHist->GetXaxis()->GetBinLowEdge(bi_it) > FitRegion_Max)) {
+ fMCHist->SetBinContent(bi_it, 0);
+ fMCHist->SetBinError(bi_it, 0);
+ fDataHist->SetBinContent(bi_it, 0);
+ fDataHist->SetBinError(bi_it, 0);
+ }
+ }
+}
+
+void Smear_SVDUnfold_Propagation_Osc::Write(void) {
+ NDToSpectrumSmearingMatrix->Write("SmearingMatrix_ND", TObject::kOverwrite);
+
TH1D *OscHist = NDToSpectrumSmearingMatrix->ProjectionY();
OscHist->Clear();
TRandom3 rnjesus;
NDDataHist->Write("ERec_ND", TObject::kOverwrite);
bool HasNegValue = false;
- int truncations = 0;
+ int truncations = TruncateStart;
do {
- if (truncations > TruncateUpTo) {
+ if (truncations >= TruncateUpTo) {
THROW("Unfolded enu spectrum had negative values even after "
<< truncations << " SVD singular value truncations.");
}
// Unfold ND ERec -> Enu spectrum
// ------------------------------
SmearceptanceUtils::PushTH1ThroughMatrixWithErrors(
NDDataHist, OscHist, NDToSpectrumResponseMatrix, 1000, false);
std::stringstream ss("");
- ss << "Unfolded_ENuTrue_ND_trunc" << truncations;
+ ss << "Unfolded_ENuTrue_ND_trunc_" << truncations;
OscHist->Write(ss.str().c_str(), TObject::kOverwrite);
+ ss.str("");
+ ss << "Response_ND_trunc_" << truncations;
+ SmearceptanceUtils::SVDGetInverse(NDToSpectrumSmearingMatrix, truncations)
+ ->Write(ss.str().c_str(), TObject::kOverwrite);
HasNegValue = false;
for (Int_t bi_it = 1; bi_it < OscHist->GetXaxis()->GetNbins() + 1;
++bi_it) {
+ if ((FitRegion_Min != 0xdeadbeef) &&
+ (OscHist->GetXaxis()->GetBinUpEdge(bi_it) <= FitRegion_Min)) {
+ continue;
+ }
+ if ((FitRegion_Max != 0xdeadbeef) &&
+ (OscHist->GetXaxis()->GetBinLowEdge(bi_it) > FitRegion_Max)) {
+ continue;
+ }
+
if (OscHist->GetBinContent(bi_it) < 0) {
HasNegValue = true;
break;
}
}
if (HasNegValue) {
TMatrixD NDToSpectrumResponseMatrix_l =
SmearceptanceUtils::GetMatrix(SmearceptanceUtils::SVDGetInverse(
NDToSpectrumSmearingMatrix, truncations));
NDToSpectrumResponseMatrix.ResizeTo(NDToSpectrumResponseMatrix_l);
NDToSpectrumResponseMatrix = NDToSpectrumResponseMatrix_l;
}
truncations++;
} while (HasNegValue);
// Apply FD/ND weights
// ------------------------------
OscHist->Write("ND_FD_Corrected_ENuTrue", TObject::kOverwrite);
// Apply Oscillations
// ------------------------------
FitWeight *fw = FitBase::GetRW();
OscWeightEngine *oscWE =
dynamic_cast<OscWeightEngine *>(fw->GetRWEngine(kOSCILLATION));
if (!oscWE) {
THROW(
"Couldn't load oscillation weight engine for sample: "
"Smear_SVDUnfold_Propagation_Osc.");
}
for (Int_t bi_it = 1; bi_it < OscHist->GetXaxis()->GetNbins() + 1; ++bi_it) {
double oscWeight =
oscWE->CalcWeight(OscHist->GetXaxis()->GetBinCenter(bi_it), 14);
OscHist->SetBinContent(bi_it, OscHist->GetBinContent(bi_it) * oscWeight);
}
OscHist->Write("ENuTrue_FD", TObject::kOverwrite);
TGraph POsc;
POsc.Set(1E4 - 1);
double min = OscHist->GetXaxis()->GetBinLowEdge(1);
double step =
(OscHist->GetXaxis()->GetBinUpEdge(OscHist->GetXaxis()->GetNbins()) -
OscHist->GetXaxis()->GetBinLowEdge(1)) /
double(1E4);
for (size_t i = 1; i < 1E4; ++i) {
double enu = min + i * step;
double ow = oscWE->CalcWeight(enu, 14);
if (ow != ow) {
std::cout << "Bad osc weight for ENu: " << enu << std::endl;
}
POsc.SetPoint(i - 1, enu, ow);
}
POsc.Write("POsc", TObject::kOverwrite);
- return;
-
// Forward fold Spectrum -> ERec FD
// ------------------------------
if (fMCHist) {
fMCHist->SetDirectory(NULL);
delete fMCHist;
}
fMCHist = static_cast<TH1D *>(FDDataHist->Clone());
SmearceptanceUtils::PushTH1ThroughMatrixWithErrors(
- OscHist, fMCHist, SpectrumToFDResponseMatrix, 1000, false);
+ OscHist, fMCHist, SpectrumToFDResponseMatrix, 1000, true);
+
+ NDToSpectrumSmearingMatrix->Write("SmearingMatrix_FD", TObject::kOverwrite);
fMCHist->Write("ERec_FD_pred", TObject::kOverwrite);
- FDDataHist->Write("ERec_FD_obs", TObject::kOverwrite);
+ fDataHist->Write("ERec_FD_obs", TObject::kOverwrite);
+
+ fMCHist->SetBinContent(0, 0);
+ fMCHist->SetBinError(0, 0);
+ fDataHist->SetBinContent(0, 0);
+ fDataHist->SetBinError(0, 0);
+
+ fMCHist->SetBinContent(fMCHist->GetXaxis()->GetNbins() + 1, 0);
+ fMCHist->SetBinError(fMCHist->GetXaxis()->GetNbins() + 1, 0);
+ fDataHist->SetBinContent(fMCHist->GetXaxis()->GetNbins() + 1, 0);
+ fDataHist->SetBinError(fMCHist->GetXaxis()->GetNbins() + 1, 0);
+
+ for (Int_t bi_it = 1; bi_it < fMCHist->GetXaxis()->GetNbins() + 1; ++bi_it) {
+ if ((FitRegion_Min != 0xdeadbeef) &&
+ (fMCHist->GetXaxis()->GetBinUpEdge(bi_it) <= FitRegion_Min)) {
+ fMCHist->SetBinContent(bi_it, 0);
+ fMCHist->SetBinError(bi_it, 0);
+ fDataHist->SetBinContent(bi_it, 0);
+ fDataHist->SetBinError(bi_it, 0);
+ }
+ if ((FitRegion_Max != 0xdeadbeef) &&
+ (fMCHist->GetXaxis()->GetBinLowEdge(bi_it) > FitRegion_Max)) {
+ fMCHist->SetBinContent(bi_it, 0);
+ fMCHist->SetBinError(bi_it, 0);
+ fDataHist->SetBinContent(bi_it, 0);
+ fDataHist->SetBinError(bi_it, 0);
+ }
+ }
+
+ fMCHist->Write("ERec_FD_pred_masked", TObject::kOverwrite);
+ fDataHist->Write("ERec_FD_obs_masked", TObject::kOverwrite);
+
+ if (ND_True_Spectrum_Hist) {
+ ND_True_Spectrum_Hist->Write("ETrue_ND_Spectrum", TObject::kOverwrite);
+ }
+ if (FD_True_Spectrum_Hist) {
+ FD_True_Spectrum_Hist->Write("ETrue_FD_Spectrum", TObject::kOverwrite);
+ }
}
diff --git a/src/MCStudies/Smear_SVDUnfold_Propagation_Osc.h b/src/MCStudies/Smear_SVDUnfold_Propagation_Osc.h
index 0457968..db8f646 100644
--- a/src/MCStudies/Smear_SVDUnfold_Propagation_Osc.h
+++ b/src/MCStudies/Smear_SVDUnfold_Propagation_Osc.h
@@ -1,47 +1,57 @@
// 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 Smear_SVDUnfold_Propagation_Osc_H_SEEN
#define Smear_SVDUnfold_Propagation_Osc_H_SEEN
#include "Measurement1D.h"
class Smear_SVDUnfold_Propagation_Osc : public Measurement1D {
public:
Smear_SVDUnfold_Propagation_Osc(nuiskey samplekey);
virtual ~Smear_SVDUnfold_Propagation_Osc(){};
void FillEventVariables(FitEvent *event);
bool isSignal(FitEvent *event);
void ConvertEventRates(void);
TH1D *NDDataHist;
TH1D *FDDataHist;
+ TH1D *ND_Unfolded_Spectrum_Hist;
+ TH1D *NDFD_Corrected_Spectrum_Hist;
+ TH1D *FD_Propagated_Spectrum_Hist;
+
+ TH1D *ND_True_Spectrum_Hist;
+ TH1D *FD_True_Spectrum_Hist;
+
TH2D *NDToSpectrumSmearingMatrix;
TMatrixD NDToSpectrumResponseMatrix;
TMatrixD SpectrumToFDResponseMatrix;
+ Int_t TruncateStart;
Int_t TruncateUpTo;
+ double FitRegion_Min;
+ double FitRegion_Max;
};
#endif
diff --git a/src/MCStudies/Smearceptance_Tester.cxx b/src/MCStudies/Smearceptance_Tester.cxx
index 6110443..28f498b 100644
--- a/src/MCStudies/Smearceptance_Tester.cxx
+++ b/src/MCStudies/Smearceptance_Tester.cxx
@@ -1,690 +1,690 @@
// 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 "Smearceptance_Tester.h"
#include "SmearceptanceUtils.h"
#include "Smearcepterton.h"
+#define DEBUG_SMEARTESTER
+
//********************************************************************
/// @brief Class to perform smearceptance MC Studies on a custom measurement
Smearceptance_Tester::Smearceptance_Tester(std::string name,
std::string inputfile, FitWeight *rw,
std::string type,
std::string fakeDataFile) {
//********************************************************************
// Measurement Details
std::vector<std::string> splitName = GeneralUtils::ParseToStr(name, "_");
size_t firstUS = name.find_first_of("_");
std::string smearceptorName = name.substr(firstUS + 1);
fName = name.substr(0, firstUS);
eventVariables = NULL;
// Define our energy range for flux calcs
EnuMin = 0.;
EnuMax = 100.; // Arbritrarily high energy limit
// 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);
this->fScaleFactor =
(GetEventHistogram()->Integral("width") * 1E-38 / (fNEvents + 0.)) /
this->TotalIntegratedFlux();
LOG(SAM) << "Smearceptance Flux Scaling Factor = " << fScaleFactor
<< std::endl;
if (fScaleFactor <= 0.0) {
ERR(WRN) << "SCALE FACTOR TOO LOW " << std::endl;
sleep(20);
}
// Setup our TTrees
this->AddEventVariablesToTree();
smearceptor = &Smearcepterton::Get().GetSmearcepter(smearceptorName);
Int_t RecNBins = 20, TrueNBins = 20;
double RecBinL = 0, TrueBinL = 0, RecBinH = 10, TrueBinH = 10;
if (Config::Get().GetConfigNode("smear.reconstructed.binning")) {
std::vector<std::string> args = GeneralUtils::ParseToStr(
Config::Get().ConfS("smear.reconstructed.binning"), ",");
RecNBins = GeneralUtils::StrToInt(args[0]);
RecBinL = GeneralUtils::StrToDbl(args[1]);
RecBinH = GeneralUtils::StrToDbl(args[2]);
TrueNBins = RecNBins;
TrueBinL = RecBinL;
TrueBinH = RecBinH;
}
if (Config::Get().GetConfigNode("smear.true.binning")) {
std::vector<std::string> args = GeneralUtils::ParseToStr(
Config::Get().ConfS("smear.true.binning"), ",");
TrueNBins = GeneralUtils::StrToInt(args[0]);
TrueBinL = GeneralUtils::StrToDbl(args[1]);
TrueBinH = GeneralUtils::StrToDbl(args[2]);
}
SVDTruncation = 0;
if (Config::Get().GetConfigNode("smear.true.binning")) {
SVDTruncation = Config::Get().ConfI("smear.SVD.truncation");
QLOG(SAM, "Applying SVD truncation of: " << SVDTruncation)
}
QLOG(SAM, "Using binning True: " << TrueNBins << ", [" << TrueBinL << " -- "
<< TrueBinH << "], Rec: " << RecNBins
<< ", [" << RecBinL << " -- " << RecBinH
<< "]");
ETrueDistrib = new TH1D("ELep_rate", ";True E_{#nu};Count", TrueNBins,
TrueBinL, TrueBinH);
ERecDistrib = new TH1D("ELepRec_rate", ";Rec E_{#nu};Count", RecNBins,
RecBinL, RecBinH);
ETrueDistrib->Sumw2();
ERecDistrib->Sumw2();
RecoSmear =
- new TH2D("ELepHadVis_Recon", ";Recon. E_{#nu};True E_{#nu}", RecNBins,
+ new TH2D("ELepHadVis_Recon", ";True E_{#nu};Recon. E_{#nu}", RecNBins,
RecBinL, RecBinH, TrueNBins, TrueBinL, TrueBinH);
+ RecoSmear->Sumw2();
}
void Smearceptance_Tester::AddEventVariablesToTree() {
// Setup the TTree to save everything
if (!eventVariables) {
FitPar::Config().out->cd();
eventVariables = new TTree((this->fName + "_VARS").c_str(),
(this->fName + "_VARS").c_str());
}
LOG(SAM) << "Adding Event Variables" << std::endl;
eventVariables->Branch("Omega_true", &Omega_true, "Omega_true/F");
eventVariables->Branch("Q2_true", &Q2_true, "Q2_true/F");
eventVariables->Branch("Mode_true", &Mode_true, "Mode_true/I");
eventVariables->Branch("EISLep_true", &EISLep_true, "EISLep_true/F");
eventVariables->Branch("HMFS_mu_true", &HMFS_mu_true);
eventVariables->Branch("HMFS_pip_true", &HMFS_pip_true);
eventVariables->Branch("HMFS_pim_true", &HMFS_pim_true);
eventVariables->Branch("HMFS_cpi_true", &HMFS_cpi_true);
eventVariables->Branch("HMFS_p_true", &HMFS_p_true);
eventVariables->Branch("KEFSHad_cpip_true", &KEFSHad_cpip_true,
"KEFSHad_cpip_true/F");
eventVariables->Branch("KEFSHad_cpim_true", &KEFSHad_cpim_true,
"KEFSHad_cpim_true/F");
eventVariables->Branch("KEFSHad_cpi_true", &KEFSHad_cpi_true,
"KEFSHad_cpi_true/F");
eventVariables->Branch("TEFSHad_pi0_true", &TEFSHad_pi0_true,
"TEFSHad_pi0_true/F");
eventVariables->Branch("KEFSHad_p_true", &KEFSHad_p_true, "KEFSHad_p_true/F");
eventVariables->Branch("KEFSHad_n_true", &KEFSHad_n_true, "KEFSHad_n_true/F");
eventVariables->Branch("EFSHad_true", &EFSHad_true, "EFSHad_true/F");
eventVariables->Branch("EFSChargedEMHad_true", &EFSChargedEMHad_true,
"EFSChargedEMHad_true/F");
eventVariables->Branch("EFSLep_true", &EFSLep_true, "EFSLep_true/F");
eventVariables->Branch("EFSgamma_true", &EFSgamma_true, "EFSgamma_true/F");
eventVariables->Branch("PDGISLep_true", &PDGISLep_true, "PDGISLep_true/I");
eventVariables->Branch("PDGFSLep_true", &PDGFSLep_true, "PDGFSLep_true/I");
eventVariables->Branch("Nprotons_true", &Nprotons_true, "Nprotons_true/I");
eventVariables->Branch("Nneutrons_true", &Nneutrons_true, "Nneutrons_true/I");
eventVariables->Branch("Ncpiplus_true", &Ncpiplus_true, "Ncpiplus_true/I");
eventVariables->Branch("Ncpiminus_true", &Ncpiminus_true, "Ncpiminus_true/I");
eventVariables->Branch("Ncpi_true", &Ncpi_true, "Ncpi_true/I");
eventVariables->Branch("Npi0_true", &Npi0_true, "Npi0_true/I");
eventVariables->Branch("HMFS_mu_rec", &HMFS_mu_rec);
eventVariables->Branch("HMFS_pip_rec", &HMFS_pip_rec);
eventVariables->Branch("HMFS_pim_rec", &HMFS_pim_rec);
eventVariables->Branch("HMFS_cpi_rec", &HMFS_cpi_rec);
eventVariables->Branch("HMFS_p_rec", &HMFS_p_rec);
eventVariables->Branch("KEFSHad_cpip_rec", &KEFSHad_cpip_rec,
"KEFSHad_cpip_rec/F");
eventVariables->Branch("KEFSHad_cpim_rec", &KEFSHad_cpim_rec,
"KEFSHad_cpim_rec/F");
eventVariables->Branch("KEFSHad_cpi_rec", &KEFSHad_cpi_rec,
"KEFSHad_cpi_rec/F");
eventVariables->Branch("TEFSHad_pi0_rec", &TEFSHad_pi0_rec,
"TEFSHad_pi0_rec/F");
eventVariables->Branch("KEFSHad_p_rec", &KEFSHad_p_rec, "KEFSHad_p_rec/F");
eventVariables->Branch("KEFSHad_n_rec", &KEFSHad_n_rec, "KEFSHad_n_rec/F");
eventVariables->Branch("EFSHad_rec", &EFSHad_rec, "EFSHad_rec/F");
eventVariables->Branch("EFSLep_rec", &EFSLep_rec, "EFSLep_rec/F");
eventVariables->Branch("EFSVis_cpip", &EFSVis_cpip, "EFSVis_cpip/F");
eventVariables->Branch("EFSVis_cpim", &EFSVis_cpim, "EFSVis_cpim/F");
eventVariables->Branch("EFSVis_cpi", &EFSVis_cpi, "EFSVis_cpi/F");
eventVariables->Branch("EFSVis_pi0", &EFSVis_pi0, "EFSVis_pi0/F");
eventVariables->Branch("EFSVis_p", &EFSVis_p, "EFSVis_p/F");
eventVariables->Branch("EFSVis_n", &EFSVis_n, "EFSVis_n/F");
eventVariables->Branch("EFSVis_gamma", &EFSVis_gamma, "EFSVis_gamma/F");
eventVariables->Branch("EFSVis_other", &EFSVis_other, "EFSVis_other/F");
eventVariables->Branch("EFSVis", &EFSVis, "EFSVis/F");
eventVariables->Branch("FSCLep_seen", &FSCLep_seen, "FSCLep_seen/I");
eventVariables->Branch("Nprotons_seen", &Nprotons_seen, "Nprotons_seen/I");
eventVariables->Branch("Nneutrons_seen", &Nneutrons_seen, "Nneutrons_seen/I");
eventVariables->Branch("Ncpip_seen", &Ncpip_seen, "Ncpip_seen/I");
eventVariables->Branch("Ncpim_seen", &Ncpim_seen, "Ncpim_seen/I");
eventVariables->Branch("Ncpi_seen", &Ncpi_seen, "Ncpi_seen/I");
eventVariables->Branch("Npi0_seen", &Npi0_seen, "Npi0_seen/I");
eventVariables->Branch("Nothers_seen", &Nothers_seen, "Nothers_seen/I");
eventVariables->Branch("EISLep_QE_rec", &EISLep_QE_rec, "EISLep_QE_rec/F");
eventVariables->Branch("EISLep_LepHad_rec", &EISLep_LepHad_rec,
"EISLep_LepHad_rec/F");
eventVariables->Branch("EISLep_LepHadVis_rec", &EISLep_LepHadVis_rec,
"EISLep_LepHadVis_rec/F");
eventVariables->Branch("Nprotons_contributed", &Nprotons_contributed,
"Nprotons_contributed/I");
eventVariables->Branch("Nneutrons_contributed", &Nneutrons_contributed,
"Nneutrons_contributed/I");
eventVariables->Branch("Ncpip_contributed", &Ncpip_contributed,
"Ncpip_contributed/I");
eventVariables->Branch("Ncpim_contributed", &Ncpim_contributed,
"Ncpim_contributed/I");
eventVariables->Branch("Ncpi_contributed", &Ncpi_contributed,
"Ncpi_contributed/I");
eventVariables->Branch("Npi0_contributed", &Npi0_contributed,
"Npi0_contributed/I");
eventVariables->Branch("Ngamma_contributed", &Ngamma_contributed,
"Ngamma_contributed/I");
eventVariables->Branch("Nothers_contibuted", &Nothers_contibuted,
"Nothers_contibuted/I");
eventVariables->Branch("Weight", &Weight, "Weight/F");
eventVariables->Branch("RWWeight", &RWWeight, "RWWeight/F");
eventVariables->Branch("InputWeight", &InputWeight, "InputWeight/F");
eventVariables->Branch("FluxWeight", &FluxWeight, "FluxWeight/F");
eventVariables->Branch("EffWeight", &EffWeight, "EffWeight/F");
eventVariables->Branch("xsecScaling", &xsecScaling, "xsecScaling/F");
eventVariables->Branch("flagCCINC_true", &flagCCINC_true, "flagCCINC_true/O");
eventVariables->Branch("flagCC0Pi_true", &flagCC0Pi_true, "flagCC0Pi_true/O");
eventVariables->Branch("flagCCINC_rec", &flagCCINC_rec, "flagCCINC_rec/O");
eventVariables->Branch("flagCC0Pi_rec", &flagCC0Pi_rec, "flagCC0Pi_rec/O");
}
template <size_t N>
int CountNPdgsSeen(RecoInfo ri, int const (&pdgs)[N]) {
int sum = 0;
for (size_t pdg_it = 0; pdg_it < N; ++pdg_it) {
sum +=
std::count(ri.RecObjClass.begin(), ri.RecObjClass.end(), pdgs[pdg_it]);
}
return sum;
}
template <size_t N>
int CountNNotPdgsSeen(RecoInfo ri, int const (&pdgs)[N]) {
int sum = 0;
for (size_t pdg_it = 0; pdg_it < N; ++pdg_it) {
sum +=
(std::count(ri.RecObjClass.begin(), ri.RecObjClass.end(), pdgs[pdg_it])
? 0
: 1);
}
return sum;
}
template <size_t N>
int CountNPdgsContributed(RecoInfo ri, int const (&pdgs)[N]) {
int sum = 0;
for (size_t pdg_it = 0; pdg_it < N; ++pdg_it) {
sum += std::count(ri.TrueContribPDGs.begin(), ri.TrueContribPDGs.end(),
pdgs[pdg_it]);
}
return sum;
}
template <size_t N>
int CountNNotPdgsContributed(RecoInfo ri, int const (&pdgs)[N]) {
int sum = 0;
for (size_t pdg_it = 0; pdg_it < N; ++pdg_it) {
sum += (std::count(ri.TrueContribPDGs.begin(), ri.TrueContribPDGs.end(),
pdgs[pdg_it])
? 0
: 1);
}
return sum;
}
TLorentzVector GetHMFSRecParticles(RecoInfo ri, int pdg) {
TLorentzVector mom(0, 0, 0, 0);
for (size_t p_it = 0; p_it < ri.RecObjMom.size(); ++p_it) {
if ((ri.RecObjClass[p_it] == pdg) &&
(mom.Mag() < ri.RecObjMom[p_it].Mag())) {
mom.SetXYZM(ri.RecObjMom[p_it].X(), ri.RecObjMom[p_it].Y(),
ri.RecObjMom[p_it].Z(),
PhysConst::GetMass(ri.RecObjClass[p_it]) * 1.0E3);
}
}
return mom;
}
template <size_t N>
double SumKE_RecoInfo(RecoInfo ri, int const (&pdgs)[N], double mass) {
double sum = 0;
for (size_t p_it = 0; p_it < ri.RecObjMom.size(); ++p_it) {
if (!std::count(pdgs, pdgs + N,
ri.RecObjClass[p_it])) { // If we don't care about this
// particle type.
continue;
}
sum += sqrt(ri.RecObjMom[p_it].Mag2() + mass * mass) - mass;
}
return sum;
}
template <size_t N>
double SumTE_RecoInfo(RecoInfo ri, int const (&pdgs)[N], double mass) {
double sum = 0;
for (size_t p_it = 0; p_it < ri.RecObjMom.size(); ++p_it) {
if (!std::count(pdgs, pdgs + N,
ri.RecObjClass[p_it])) { // If we don't care about this
// particle type.
continue;
}
sum += sqrt(ri.RecObjMom[p_it].Mag2() + mass * mass);
}
return sum;
}
template <size_t N>
double SumVisE_RecoInfo(RecoInfo ri, int const (&pdgs)[N]) {
double sum = 0;
for (size_t p_it = 0; p_it < ri.RecVisibleEnergy.size(); ++p_it) {
if (!std::count(pdgs, pdgs + N,
ri.TrueContribPDGs[p_it])) { // If we don't care about this
// particle type.
continue;
}
sum += ri.RecVisibleEnergy[p_it];
}
return sum;
}
template <size_t N>
double SumVisE_RecoInfo_NotPdgs(RecoInfo ri, int const (&pdgs)[N]) {
double sum = 0;
for (size_t p_it = 0; p_it < ri.RecVisibleEnergy.size(); ++p_it) {
if (std::count(pdgs, pdgs + N,
ri.TrueContribPDGs[p_it])) { // If we know about this
// particle type.
continue;
}
sum += ri.RecVisibleEnergy[p_it];
}
return sum;
}
//********************************************************************
void Smearceptance_Tester::FillEventVariables(FitEvent *event) {
//********************************************************************
static int const cpipPDG[] = {211};
static int const cpimPDG[] = {-211};
static int const pi0PDG[] = {111};
static int const ProtonPDG[] = {2212};
static int const NeutronPDG[] = {2112};
static int const GammaPDG[] = {22};
static int const CLeptonPDGs[] = {11, 13, 15};
- static int const ExplicitPDGs[] = {211, -211, 11, 2212, 2112, 22,
- 11, 13, 15, 12, 14, 16};
+ static int const ExplicitPDGs[] = {211, -211, 111, 2212, 2112, 22,
+ 11, 13, 15, 12, 14, 16};
RecoInfo *ri = smearceptor->Smearcept(event);
HMFS_mu_true = TLorentzVector(0, 0, 0, 0);
HMFS_mu_rec = TLorentzVector(0, 0, 0, 0);
FitParticle *fsMu = event->GetHMFSMuon();
if (fsMu) {
HMFS_mu_true = fsMu->P4();
HMFS_mu_rec = GetHMFSRecParticles(*ri, 13);
}
HMFS_pip_true = TLorentzVector(0, 0, 0, 0);
HMFS_pip_rec = TLorentzVector(0, 0, 0, 0);
FitParticle *fsPip = event->GetHMFSPiPlus();
if (fsPip) {
HMFS_pip_true = fsPip->P4();
HMFS_pip_rec = GetHMFSRecParticles(*ri, 211);
}
HMFS_pim_true = TLorentzVector(0, 0, 0, 0);
HMFS_pim_rec = TLorentzVector(0, 0, 0, 0);
FitParticle *fsPim = event->GetHMFSPiMinus();
if (fsPim) {
HMFS_pim_true = fsPim->P4();
HMFS_pim_rec = GetHMFSRecParticles(*ri, -211);
}
HMFS_cpi_true = TLorentzVector(0, 0, 0, 0);
HMFS_cpi_rec = TLorentzVector(0, 0, 0, 0);
if (fsPip || fsPim) {
if (!fsPip) {
HMFS_cpi_true = HMFS_pim_true;
HMFS_cpi_rec = HMFS_pim_rec;
} else if (!fsPim) {
HMFS_cpi_true = HMFS_pip_true;
HMFS_cpi_rec = HMFS_pip_rec;
} else {
HMFS_cpi_true =
(fsPip->p2() > fsPim->p2()) ? HMFS_pip_true : HMFS_pim_true;
HMFS_cpi_rec = (fsPip->p2() > fsPim->p2()) ? HMFS_pip_rec : HMFS_pim_rec;
}
}
HMFS_p_true = TLorentzVector(0, 0, 0, 0);
HMFS_p_rec = TLorentzVector(0, 0, 0, 0);
FitParticle *fsP = event->GetHMFSProton();
if (fsP) {
HMFS_p_true = fsP->P4();
HMFS_p_rec = GetHMFSRecParticles(*ri, 2212);
}
TLorentzVector FourMomentumTransfer =
(event->GetHMISAnyLeptons()->P4() - event->GetHMFSAnyLeptons()->P4());
Omega_true = FourMomentumTransfer.E();
Q2_true = -1 * FourMomentumTransfer.Mag2();
Mode_true = event->Mode;
EISLep_true = event->GetHMISAnyLeptons()->E();
KEFSHad_cpip_true = FitUtils::SumTE_PartVect(event->GetAllFSPiPlus());
KEFSHad_cpim_true = FitUtils::SumTE_PartVect(event->GetAllFSPiMinus());
KEFSHad_cpi_true = KEFSHad_cpip_true + KEFSHad_cpim_true;
TEFSHad_pi0_true = FitUtils::SumTE_PartVect(event->GetAllFSPiZero());
KEFSHad_p_true = FitUtils::SumKE_PartVect(event->GetAllFSProton());
KEFSHad_n_true = FitUtils::SumKE_PartVect(event->GetAllFSNeutron());
EFSHad_true =
KEFSHad_cpi_true + TEFSHad_pi0_true + KEFSHad_p_true + KEFSHad_n_true;
EFSChargedEMHad_true = KEFSHad_cpi_true + TEFSHad_pi0_true + KEFSHad_p_true;
EFSLep_true = event->GetHMFSAnyLeptons()->E();
EFSgamma_true = FitUtils::SumTE_PartVect(event->GetAllFSPhoton());
PDGISLep_true = event->GetHMISAnyLeptons()->PDG();
PDGFSLep_true = event->GetHMFSAnyLeptons()->PDG();
Nprotons_true = event->GetAllFSProton().size();
Nneutrons_true = event->GetAllFSNeutron().size();
Ncpiplus_true = event->GetAllFSPiPlus().size();
Ncpiminus_true = event->GetAllFSPiMinus().size();
Ncpi_true = Ncpiplus_true + Ncpiminus_true;
Npi0_true = event->GetAllFSPiZero().size();
KEFSHad_cpip_rec =
SumKE_RecoInfo(*ri, cpipPDG, PhysConst::mass_cpi * PhysConst::mass_MeV);
KEFSHad_cpim_rec =
SumKE_RecoInfo(*ri, cpimPDG, PhysConst::mass_cpi * PhysConst::mass_MeV);
KEFSHad_cpi_rec = KEFSHad_cpip_rec + KEFSHad_cpim_rec;
TEFSHad_pi0_rec =
SumTE_RecoInfo(*ri, pi0PDG, PhysConst::mass_pi0 * PhysConst::mass_MeV);
KEFSHad_p_rec = SumKE_RecoInfo(*ri, ProtonPDG,
PhysConst::mass_proton * PhysConst::mass_MeV);
KEFSHad_n_rec = SumKE_RecoInfo(*ri, NeutronPDG,
PhysConst::mass_neutron * PhysConst::mass_MeV);
EFSHad_rec =
KEFSHad_cpi_rec + TEFSHad_pi0_rec + KEFSHad_p_rec + KEFSHad_n_rec;
TLorentzVector FSLepMom_rec(0, 0, 0, 0);
if (event->GetHMFSAnyLeptons()) {
FSLepMom_rec = GetHMFSRecParticles(*ri, event->GetHMFSAnyLeptons()->PDG());
EFSLep_rec = FSLepMom_rec.E();
} else {
EFSLep_rec = 0;
}
EFSVis_cpip = SumVisE_RecoInfo(*ri, cpipPDG);
EFSVis_cpim = SumVisE_RecoInfo(*ri, cpimPDG);
EFSVis_cpi = EFSVis_cpip + EFSVis_cpim;
EFSVis_pi0 = SumVisE_RecoInfo(*ri, pi0PDG);
EFSVis_p = SumVisE_RecoInfo(*ri, ProtonPDG);
EFSVis_n = SumVisE_RecoInfo(*ri, NeutronPDG);
EFSVis_gamma = SumVisE_RecoInfo(*ri, GammaPDG);
EFSVis_other = SumVisE_RecoInfo_NotPdgs(*ri, ExplicitPDGs);
- EFSVis = EFSVis_cpip + EFSVis_cpim + EFSVis_cpi + EFSVis_pi0 + EFSVis_p +
- EFSVis_n + EFSVis_gamma;
+ EFSVis = EFSVis_cpi + EFSVis_pi0 + EFSVis_p + EFSVis_n + EFSVis_gamma;
FSCLep_seen = CountNPdgsSeen(*ri, CLeptonPDGs);
Nprotons_seen = CountNPdgsSeen(*ri, ProtonPDG);
Nneutrons_seen = CountNPdgsSeen(*ri, NeutronPDG);
Ncpip_seen = CountNPdgsSeen(*ri, cpipPDG);
Ncpim_seen = CountNPdgsSeen(*ri, cpimPDG);
Ncpi_seen = Ncpip_seen + Ncpim_seen;
Npi0_seen = CountNPdgsSeen(*ri, pi0PDG);
Nothers_seen = CountNNotPdgsSeen(*ri, ExplicitPDGs);
if (FSCLep_seen && (FSLepMom_rec.Mag() > 1E-8)) {
EISLep_QE_rec =
FitUtils::EnuQErec(FSLepMom_rec.Mag() / 1000.0, FSLepMom_rec.CosTheta(),
34, PDGFSLep_true > 0) *
1000.0;
} else {
EISLep_QE_rec = 0;
}
EISLep_LepHad_rec = EFSLep_rec + EFSHad_rec;
EISLep_LepHadVis_rec = EFSLep_rec + EFSHad_rec + EFSVis;
Nprotons_contributed = CountNPdgsContributed(*ri, ProtonPDG);
Nneutrons_contributed = CountNPdgsContributed(*ri, NeutronPDG);
Ncpip_contributed = CountNPdgsContributed(*ri, cpipPDG);
Ncpim_contributed = CountNPdgsContributed(*ri, cpimPDG);
Ncpi_contributed = Ncpip_contributed + Ncpim_contributed;
Npi0_contributed = CountNPdgsContributed(*ri, pi0PDG);
Ngamma_contributed = CountNPdgsContributed(*ri, GammaPDG);
Nothers_contibuted = CountNNotPdgsContributed(*ri, ExplicitPDGs);
Weight = event->RWWeight * event->InputWeight;
RWWeight = event->RWWeight;
InputWeight = event->InputWeight;
FluxWeight = GetFluxHistogram()->GetBinContent(
GetFluxHistogram()->FindBin(EISLep_true)) /
GetFluxHistogram()->Integral();
EffWeight = ri->Weight;
xsecScaling = fScaleFactor;
flagCCINC_true = PDGFSLep_true & 1;
flagCC0Pi_true = (Ncpi_true + Npi0_true) == 0;
flagCCINC_rec = FSCLep_seen && PDGFSLep_true & 1;
flagCC0Pi_rec = ((Ncpi_seen + Npi0_seen) == 0) && flagCCINC_rec;
// Fill the eventVariables Tree
eventVariables->Fill();
- RecoSmear->Fill(flagCCINC_rec ? EISLep_LepHadVis_rec / 1000.0 : -1,
- EISLep_true / 1000.0, Weight);
+ RecoSmear->Fill(EISLep_true / 1000.0,
+ flagCCINC_rec ? EISLep_LepHadVis_rec / 1000.0 : -1, Weight);
ETrueDistrib->Fill(EISLep_true / 1000.0, flagCCINC_true ? Weight : 0);
ERecDistrib->Fill(EISLep_LepHadVis_rec / 1000.0, flagCCINC_rec ? Weight : 0);
return;
};
//********************************************************************
void Smearceptance_Tester::Write(std::string drawOpt) {
//********************************************************************
// First save the TTree
eventVariables->Write();
// Save Flux and Event Histograms too
GetInput()->GetFluxHistogram()->Write();
GetInput()->GetEventHistogram()->Write();
TH2D *SmearMatrix_ev =
static_cast<TH2D *>(RecoSmear->Clone("ELepHadVis_Smear_ev"));
for (Int_t trueAxis_it = 1;
- trueAxis_it < RecoSmear->GetYaxis()->GetNbins() + 1; ++trueAxis_it) {
+ trueAxis_it < RecoSmear->GetXaxis()->GetNbins() + 1; ++trueAxis_it) {
double NEISLep = ETrueDistrib->GetBinContent(trueAxis_it);
for (Int_t recoAxis_it = 1;
- recoAxis_it < RecoSmear->GetXaxis()->GetNbins() + 1; ++recoAxis_it) {
+ recoAxis_it < RecoSmear->GetYaxis()->GetNbins() + 1; ++recoAxis_it) {
if (NEISLep > std::numeric_limits<double>::epsilon()) {
SmearMatrix_ev->SetBinContent(
- recoAxis_it, trueAxis_it,
- SmearMatrix_ev->GetBinContent(recoAxis_it, trueAxis_it) / NEISLep);
+ trueAxis_it, recoAxis_it,
+ SmearMatrix_ev->GetBinContent(trueAxis_it, recoAxis_it) / NEISLep);
}
}
}
ETrueDistrib->Write();
ERecDistrib->Write();
RecoSmear->Write();
SmearMatrix_ev->Write();
TH2D *ResponseMatrix_ev =
SmearceptanceUtils::SVDGetInverse(SmearMatrix_ev, SVDTruncation);
- ResponseMatrix_ev = SmearceptanceUtils::SwapXYTH2D(ResponseMatrix_ev);
ResponseMatrix_ev->SetName("ResponseMatrix_ev");
ResponseMatrix_ev->Write();
#ifdef DEBUG_SMEARTESTER
- TMatrixD SmearMatrix_ev_md = SmearceptanceUtils::GetMatrix(
- SmearceptanceUtils::SwapXYTH2D(SmearMatrix_ev));
+ TMatrixD SmearMatrix_ev_md = SmearceptanceUtils::GetMatrix(SmearMatrix_ev);
TH1D *SmearedEvt = static_cast<TH1D *>(ERecDistrib->Clone());
SmearedEvt->SetNameTitle("SmearedEvt", ";Rec E_{#nu}; count");
SmearceptanceUtils::PushTH1ThroughMatrixWithErrors(
ETrueDistrib, SmearedEvt, SmearMatrix_ev_md, 5000, false);
SmearedEvt->Write();
SmearedEvt->Scale(1, "width");
SmearedEvt->SetName("SmearedEvt_bw");
SmearedEvt->Write();
#endif
TMatrixD ResponseMatrix_evt_md =
SmearceptanceUtils::GetMatrix(ResponseMatrix_ev);
TH1D *Unfolded_enu_obs = static_cast<TH1D *>(ETrueDistrib->Clone());
Unfolded_enu_obs->SetNameTitle("UnfoldedENu_evt", ";True E_{#nu};count");
SmearceptanceUtils::PushTH1ThroughMatrixWithErrors(
ERecDistrib, Unfolded_enu_obs, ResponseMatrix_evt_md, 5000, false);
Unfolded_enu_obs->Write();
Unfolded_enu_obs->Scale(1, "width");
Unfolded_enu_obs->SetName("UnfoldedENu_evt_bw");
Unfolded_enu_obs->Write();
ETrueDistrib->Scale(1, "width");
ETrueDistrib->SetName("ELep_rate_bw");
ETrueDistrib->Write();
ERecDistrib->Scale(1, "width");
ERecDistrib->SetName("ELepRec_rate_bw");
ERecDistrib->Write();
}
// -------------------------------------------------------------------
// Purely MC Plot
// Following functions are just overrides to handle this
// -------------------------------------------------------------------
//********************************************************************
/// Everything is classed as signal...
bool Smearceptance_Tester::isSignal(FitEvent *event) {
//********************************************************************
(void)event;
return true;
};
//********************************************************************
void Smearceptance_Tester::ScaleEvents() {
//********************************************************************
// Saving everything to a TTree so no scaling required
return;
}
//********************************************************************
void Smearceptance_Tester::ApplyNormScale(float norm) {
//********************************************************************
// Saving everything to a TTree so no scaling required
this->fCurrentNorm = norm;
return;
}
//********************************************************************
void Smearceptance_Tester::FillHistograms() {
//********************************************************************
// No Histograms need filling........
return;
}
//********************************************************************
void Smearceptance_Tester::ResetAll() {
//********************************************************************
eventVariables->Reset();
return;
}
//********************************************************************
float Smearceptance_Tester::GetChi2() {
//********************************************************************
// No Likelihood to test, purely MC
return 0.0;
}
diff --git a/src/Smearceptance/EfficiencyApplicator.cxx b/src/Smearceptance/EfficiencyApplicator.cxx
index b046cb5..55ab9b3 100644
--- a/src/Smearceptance/EfficiencyApplicator.cxx
+++ b/src/Smearceptance/EfficiencyApplicator.cxx
@@ -1,388 +1,388 @@
// 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 "EfficiencyApplicator.h"
#include "TEfficiency.h"
#include "TH2.h"
#include "TH3.h"
//#define DEBUG_EFFAPP
EfficiencyApplicator::DependVar GetVarType(std::string const &axisvar) {
if (axisvar == "kMomentum") {
return EfficiencyApplicator::kMomentum;
} else if (axisvar == "kKE") {
return EfficiencyApplicator::kKE;
} else if (axisvar == "kTheta") {
return EfficiencyApplicator::kTheta;
} else if (axisvar == "kCosTheta") {
return EfficiencyApplicator::kCosTheta;
} else if (axisvar == "kPhi") {
return EfficiencyApplicator::kPhi;
}
return EfficiencyApplicator::kNoAxis;
}
TH1 *GetEffHist(TFile *inputFile, std::string const &HistName) {
TH1 *hist = dynamic_cast<TH1 *>(inputFile->Get(HistName.c_str()));
if (hist) {
return hist;
}
TEfficiency *effHist =
dynamic_cast<TEfficiency *>(inputFile->Get(HistName.c_str()));
if (effHist) {
TH1D *numer = dynamic_cast<TH1D *>(effHist->GetCopyPassedHisto());
TH1D *denom = dynamic_cast<TH1D *>(effHist->GetCopyTotalHisto());
if (numer && denom) {
numer->Divide(denom);
denom->SetDirectory(NULL);
delete denom;
// Gonna be a memory leak, but I'll get over it
numer->SetDirectory(NULL);
return numer;
}
ERROR(FTL, "TEfficiency internal histograms were not TH1Ds.");
}
THROW("Couldn't get appropriate efficiency object named "
<< HistName << " from input file " << inputFile->GetName());
}
/// Reads particle efficiency nodes
///
/// Nodes look like:
/// <EfficiencyApplicator Name="D00N_ND_LAr">
/// <EfficiencyCurve PDG="211,-211" InputFile="effs.root"
/// HistName="cpion_eff_mom_ctheta" NDims="2" XAxis="kMomentum"
/// YAxis="kCosTheta" Interpolate="false" />
/// <!-- Can also contain ThresholdAccepter elements as below-->
/// <RecoThreshold PDG="2212" RecoThresholdAbsCosTheta_Min="0" />
/// <VisThreshold PDG="2212" VisThresholdKE_MeV="10" Contrib="K" />
/// </EfficiencyApplicator>
void EfficiencyApplicator::SpecifcSetup(nuiskey &nk) {
rand.~TRandom3();
new (&rand) TRandom3();
std::vector<nuiskey> effDescriptors =
nk.GetListOfChildNodes("EfficiencyCurve");
for (size_t t_it = 0; t_it < effDescriptors.size(); ++t_it) {
std::string inputFileName = effDescriptors[t_it].GetS("InputFile");
std::string HistName = effDescriptors[t_it].GetS("HistName");
TFile inputFile(inputFileName.c_str());
if (!inputFile.IsOpen()) {
THROW("Couldn't open specified input root file: " << inputFileName);
}
TH1 *inpHist = GetEffHist(&inputFile, HistName);
if (!inpHist) {
THROW("Couldn't get TH1D named: " << HistName << " from input root file: "
<< inputFileName);
}
int NDims = effDescriptors[t_it].GetI("NDims");
if (NDims < 1 || NDims > 3) {
THROW("Read NDims attribute as: " << NDims << ", efficiency curve can "
"currently have between 1 "
"and 3 dimensions.");
}
EfficiencyApplicator::DependVar XVar =
GetVarType(effDescriptors[t_it].GetS("XAxis"));
double XAxisScale = effDescriptors[t_it].Has("XAxisScaleToInternal")
? effDescriptors[t_it].GetD("XAxisScaleToInternal")
: 1;
EfficiencyApplicator::DependVar YVar =
NDims > 1 ? GetVarType(effDescriptors[t_it].GetS("YAxis"))
: EfficiencyApplicator::kNoAxis;
double YAxisScale = effDescriptors[t_it].Has("YAxisScaleToInternal")
? effDescriptors[t_it].GetD("YAxisScaleToInternal")
: 1;
EfficiencyApplicator::DependVar ZVar =
NDims > 2 ? GetVarType(effDescriptors[t_it].GetS("ZAxis"))
: EfficiencyApplicator::kNoAxis;
double ZAxisScale = effDescriptors[t_it].Has("ZAxisScaleToInternal")
? effDescriptors[t_it].GetD("ZAxisScaleToInternal")
: 1;
bool Interpolate = effDescriptors[t_it].Has("Interpolate") &&
effDescriptors[t_it].GetI("Interpolate");
// bool ApplyAsWeight = effDescriptors[t_it].Has("ApplyAsWeight") &&
// effDescriptors[t_it].GetI("ApplyAsWeight");
std::string pdgs_s = effDescriptors[t_it].GetS("PDG");
std::vector<int> pdgs_i = GeneralUtils::ParseToInt(pdgs_s, ",");
for (size_t pdg_it = 0; pdg_it < pdgs_i.size(); ++pdg_it) {
if (Efficiencies.count(pdgs_i[pdg_it])) {
ERROR(WRN, "Smearceptor " << ElementName << ":" << InstanceName
<< " already has a efficiency for PDG: "
<< pdgs_i[pdg_it]);
}
EffMap em;
em.EffCurve = static_cast<TH1 *>(inpHist->Clone());
em.EffCurve->SetDirectory(NULL);
em.Interpolate = Interpolate;
// em.ApplyAsWeight = ApplyAsWeight;
em.NDims = NDims;
em.DependVars[0] = XVar;
em.DependVars[1] = YVar;
em.DependVars[2] = ZVar;
em.AxisScales[0] = XAxisScale;
em.AxisScales[1] = YAxisScale;
em.AxisScales[2] = ZAxisScale;
Efficiencies[pdgs_i[pdg_it]] = em;
- QLOG(SAM,
+ QLOG(FIT,
"Added reconstruction efficiency curve for PDG: " << pdgs_i[pdg_it]);
}
}
SlaveTA.Setup(nk);
}
RecoInfo *EfficiencyApplicator::Smearcept(FitEvent *fe) {
RecoInfo *ri = new RecoInfo();
for (size_t p_it = 0; p_it < fe->NParticles(); ++p_it) {
FitParticle *fp = fe->GetParticle(p_it);
#ifdef DEBUG_EFFAPP
std::cout << std::endl;
std::cout << "[" << p_it << "]: " << fp->PDG() << ", " << fp->Status()
<< ", " << fp->E() << " -- KE:" << fp->KE()
<< " Mom: " << fp->P3().Mag() << std::flush;
#endif
if (fp->Status() != kFinalState) {
#ifdef DEBUG_EFFAPP
std::cout << " -- Not final state." << std::flush;
#endif
continue;
}
if (!Efficiencies.count(fp->PDG())) {
SlaveTA.SmearceptOneParticle(ri, fp
#ifdef DEBUG_THRESACCEPT
,
p_it
#endif
);
continue;
}
EffMap &em = Efficiencies[fp->PDG()];
double kineProps[3];
for (Int_t dim_it = 0; dim_it < em.NDims; ++dim_it) {
switch (em.DependVars[dim_it]) {
case kMomentum: {
kineProps[dim_it] = fp->P3().Mag();
break;
}
case kKE: {
kineProps[dim_it] = fp->KE();
break;
}
case kTheta: {
kineProps[dim_it] = fp->P3().Theta();
break;
}
case kCosTheta: {
kineProps[dim_it] = fp->P3().CosTheta();
break;
}
case kPhi: {
kineProps[dim_it] = fp->P3().Phi();
break;
}
default: { THROW("Trying to find particle value for a kNoAxis."); }
}
kineProps[dim_it] /= em.AxisScales[dim_it];
}
double effProb = 0;
switch (em.NDims) {
case 1: {
TH1 *hist = em.EffCurve;
if (em.Interpolate &&
(hist->GetXaxis()->GetBinCenter(1) < kineProps[0]) &&
(hist->GetXaxis()->GetBinCenter(hist->GetXaxis()->GetNbins()) >
kineProps[0])) {
effProb = hist->Interpolate(kineProps[0]);
} else {
Int_t xbin = hist->GetXaxis()->FindFixBin(kineProps[0]);
if (!xbin || ((hist->GetXaxis()->GetNbins() + 1) == xbin)) {
ERROR(WRN, "Tried to apply effiency but XBin: "
<< xbin << " is outside range (/"
<< hist->GetXaxis()->GetNbins() << "): Prop "
<< kineProps[0] << ", ["
<< hist->GetXaxis()->GetBinLowEdge(1) << " -- "
<< hist->GetXaxis()->GetBinUpEdge(
hist->GetXaxis()->GetNbins()));
}
effProb = hist->GetBinContent(xbin);
}
break;
}
case 2: {
TH2 *hist = static_cast<TH2 *>(em.EffCurve);
if (em.Interpolate &&
(hist->GetXaxis()->GetBinCenter(1) < kineProps[0]) &&
(hist->GetXaxis()->GetBinCenter(hist->GetXaxis()->GetNbins()) >
kineProps[0]) &&
(hist->GetYaxis()->GetBinCenter(1) < kineProps[1]) &&
(hist->GetYaxis()->GetBinCenter(hist->GetYaxis()->GetNbins()) >
kineProps[1])) {
effProb = hist->Interpolate(kineProps[0], kineProps[1]);
} else {
Int_t xbin = hist->GetXaxis()->FindFixBin(kineProps[0]);
Int_t ybin = hist->GetYaxis()->FindFixBin(kineProps[1]);
if (!xbin || ((hist->GetXaxis()->GetNbins() + 1) == xbin)) {
ERROR(WRN, "Tried to apply effiency but XBin: "
<< xbin << " is outside range (/"
<< hist->GetXaxis()->GetNbins() << "): Prop "
<< kineProps[0] << ", ["
<< hist->GetXaxis()->GetBinLowEdge(1) << " -- "
<< hist->GetXaxis()->GetBinUpEdge(
hist->GetXaxis()->GetNbins()));
}
if (!ybin || ((hist->GetYaxis()->GetNbins() + 1) == ybin)) {
ERROR(WRN, "Tried to apply effiency but XBin: "
<< ybin << " is outside range (/"
<< hist->GetYaxis()->GetNbins() << "): Prop "
<< kineProps[0] << ", ["
<< hist->GetYaxis()->GetBinLowEdge(1) << " -- "
<< hist->GetYaxis()->GetBinUpEdge(
hist->GetYaxis()->GetNbins()));
}
effProb = hist->GetBinContent(xbin, ybin);
#ifdef DEBUG_EFFAPP
std::cout << "\t\t: XProp: " << kineProps[0]
<< ", YProp: " << kineProps[1] << " x/y bins: " << xbin
<< "/" << ybin << ". Prop ? " << effProb << std::endl;
#endif
}
break;
}
case 3: {
TH3 *hist = static_cast<TH3 *>(em.EffCurve);
if (em.Interpolate &&
(hist->GetXaxis()->GetBinCenter(1) < kineProps[0]) &&
(hist->GetXaxis()->GetBinCenter(hist->GetXaxis()->GetNbins()) >
kineProps[0]) &&
(hist->GetYaxis()->GetBinCenter(1) < kineProps[1]) &&
(hist->GetYaxis()->GetBinCenter(hist->GetYaxis()->GetNbins()) >
kineProps[2]) &&
(hist->GetZaxis()->GetBinCenter(hist->GetZaxis()->GetNbins()) >
kineProps[2])) {
effProb = hist->Interpolate(kineProps[0], kineProps[1], kineProps[2]);
} else {
Int_t xbin = hist->GetXaxis()->FindFixBin(kineProps[0]);
Int_t ybin = hist->GetYaxis()->FindFixBin(kineProps[1]);
Int_t zbin = hist->GetZaxis()->FindFixBin(kineProps[2]);
if (!xbin || ((hist->GetXaxis()->GetNbins() + 1) == xbin)) {
ERROR(WRN, "Tried to apply effiency but XBin: "
<< xbin << " is outside range (/"
<< hist->GetXaxis()->GetNbins() << "): Prop "
<< kineProps[0] << ", ["
<< hist->GetXaxis()->GetBinLowEdge(1) << " -- "
<< hist->GetXaxis()->GetBinUpEdge(
hist->GetXaxis()->GetNbins()));
}
if (!ybin || ((hist->GetYaxis()->GetNbins() + 1) == ybin)) {
ERROR(WRN, "Tried to apply effiency but XBin: "
<< ybin << " is outside range (/"
<< hist->GetYaxis()->GetNbins() << "): Prop "
<< kineProps[0] << ", ["
<< hist->GetYaxis()->GetBinLowEdge(1) << " -- "
<< hist->GetYaxis()->GetBinUpEdge(
hist->GetYaxis()->GetNbins()));
}
if (!zbin || ((hist->GetZaxis()->GetNbins() + 1) == zbin)) {
ERROR(WRN, "Tried to apply effiency but ZBin: "
<< zbin << " is outside range (/"
<< hist->GetZaxis()->GetNbins() << "): Prop "
<< kineProps[0] << ", ["
<< hist->GetZaxis()->GetBinLowEdge(1) << " -- "
<< hist->GetZaxis()->GetBinUpEdge(
hist->GetZaxis()->GetNbins()));
}
effProb = hist->GetBinContent(xbin, ybin, zbin);
}
break;
}
}
bool accepted = (rand.Uniform() < effProb);
if (accepted) {
#ifdef DEBUG_EFFAPP
std::cout << " -- Reconstructed with probability: " << effProb
<< std::flush;
#endif
ri->RecObjMom.push_back(fp->P3());
ri->RecObjClass.push_back(fp->PDG());
continue;
}
#ifdef DEBUG_EFFAPP
std::cout << " -- Rejected with probability: " << effProb << std::flush;
#endif
SlaveTA.SmearceptOneParticle(ri, fp
#ifdef DEBUG_THRESACCEPT
,
p_it
#endif
);
}
#ifdef DEBUG_EFFAPP
std::cout << std::endl;
#endif
#ifdef DEBUG_EFFAPP
std::cout << "Reconstructed " << ri->RecObjMom.size() << " particles. "
<< std::endl;
#endif
return ri;
}
diff --git a/src/Smearceptance/EnergyShuffler.cxx b/src/Smearceptance/EnergyShuffler.cxx
index 4fc0977..a4dc6b2 100644
--- a/src/Smearceptance/EnergyShuffler.cxx
+++ b/src/Smearceptance/EnergyShuffler.cxx
@@ -1,155 +1,157 @@
#include "EnergyShuffler.h"
/// Node looks like
/// <EnergyShuffler>
/// <Shuffler From="2212" To="2112" Fraction="0.2" />
/// <Shuffler From="211,-211" To="" Fraction="0.5" />
/// </EnergyShuffler>
void EnergyShuffler::Setup(nuiskey &nk) {
std::vector<nuiskey> shuffleDescriptors = nk.GetListOfChildNodes("Shuffler");
for (size_t t_it = 0; t_it < shuffleDescriptors.size(); ++t_it) {
if (!shuffleDescriptors[t_it].Has("From") ||
!shuffleDescriptors[t_it].Has("Fraction")) {
THROW(
"Shuffler element must have at least the From and Fraction "
"attributes.");
}
std::string from_pdgs_s = shuffleDescriptors[t_it].GetS("From");
std::vector<int> from_pdgs_i = GeneralUtils::ParseToInt(from_pdgs_s, ",");
if (!from_pdgs_i.size()) {
THROW("Shuffler element must have at least one From PDG specified.");
}
std::vector<int> to_pdgs_i;
if (shuffleDescriptors[t_it].Has("To")) {
std::string to_pdgs_s = shuffleDescriptors[t_it].GetS("To");
to_pdgs_i = GeneralUtils::ParseToInt(to_pdgs_s, ",");
}
double Fraction = shuffleDescriptors[t_it].GetD("Fraction");
for (size_t f_it = 0; f_it < from_pdgs_i.size(); ++f_it) {
ShuffleDescriptor sd;
sd.ToPDGs = to_pdgs_i;
sd.EFraction = Fraction;
ShufflersDescriptors.push_back(std::make_pair(from_pdgs_i[f_it], sd));
QLOG(FIT, "\tAdded EnergyShuffler from "
<< from_pdgs_i[f_it] << " to " << to_pdgs_i.size()
<< " particle species at " << sd.EFraction
<< " KE fraction.")
}
}
}
void EnergyShuffler::DoTheShuffle(FitEvent *fe) {
std::map<Int_t, double> TakenEnergy;
std::map<Int_t, Int_t> NumParticlesToShare;
// For each particle.
// If in a from: take energy and add it to a sum.
// Count particles of each species.
for (size_t p_it = 0; p_it < fe->NParticles(); ++p_it) {
FitParticle *fp = fe->GetParticle(p_it);
if (fp->Status() != kFinalState) {
continue;
}
int PDG = fp->PDG();
for (size_t sh_it = 0; sh_it < ShufflersDescriptors.size(); ++sh_it) {
ShuffleDescriptor &sd = ShufflersDescriptors[sh_it].second;
if (std::find(sd.ToPDGs.begin(), sd.ToPDGs.end(), PDG) !=
sd.ToPDGs.end()) { // If this is a particle
// that we're giving energy
// to, take note.
if (!NumParticlesToShare.count(sh_it)) {
NumParticlesToShare[sh_it] = 0;
}
NumParticlesToShare[sh_it]++;
#ifdef DEBUG_ESHUFFLER
std::cout << "Found recieving particle in pool " << sh_it << "."
<< std::endl;
#endif
}
if (ShufflersDescriptors[sh_it].first != PDG) {
continue;
}
double KETaken = sd.EFraction * fp->KE();
if (!TakenEnergy.count(sh_it)) {
TakenEnergy[sh_it] = 0;
}
TakenEnergy[sh_it] += KETaken;
#ifdef DEBUG_ESHUFFLER
std::cout << "Taking " << KETaken << " KE from " << fp->PDG() << " ("
<< fp << ") with " << fp->KE() << " (mom: " << fp->p() << ")."
<< std::endl;
#endif
- fp->RemoveKE(KETaken);
+ fe->RemoveKE(p_it, KETaken);
+ fp = fe->GetParticle(p_it);
#ifdef DEBUG_ESHUFFLER
std::cout << "\t->" << fp->KE() << " (mom: " << fp->p() << ") => "
<< sh_it << "." << std::endl;
#endif
}
}
double LostEnergy = 0;
for (std::map<Int_t, double>::iterator te_it = TakenEnergy.begin();
te_it != TakenEnergy.end(); ++te_it) {
double EToGive = te_it->second;
// Get energy share for each 'to' particle
if (NumParticlesToShare.count(te_it->first)) { // If we have any particles
// to share the energy
// between.
#ifdef DEBUG_ESHUFFLER
std::cout << "Energy from Shuffler " << te_it->first << " " << EToGive
<< " shared between " << NumParticlesToShare[te_it->first]
<< " particle." << std::endl;
#endif
EToGive /= double(NumParticlesToShare[te_it->first]);
} else {
LostEnergy += EToGive;
continue;
}
ShuffleDescriptor &sd = ShufflersDescriptors[te_it->first].second;
for (size_t p_it = 0; p_it < fe->NParticles(); ++p_it) {
FitParticle *fp = fe->GetParticle(p_it);
if (fp->Status() != kFinalState) {
continue;
}
int PDG = fp->PDG();
if (std::find(sd.ToPDGs.begin(), sd.ToPDGs.end(), PDG) ==
sd.ToPDGs.end()) { // This shuffler has no energy for this particle
continue;
}
#ifdef DEBUG_ESHUFFLER
std::cout << "Giving " << EToGive << " KE to " << fp->PDG() << " with "
<< fp->KE() << " (mom: " << fp->p() << ")." << std::endl;
#endif
- fp->GiveKE(EToGive);
+ fe->GiveKE(p_it, EToGive);
+ fp = fe->GetParticle(p_it);
#ifdef DEBUG_ESHUFFLER
std::cout << "\t->" << fp->KE() << " (mom: " << fp->p() << ")."
<< std::endl;
#endif
}
}
#ifdef DEBUG_ESHUFFLER
if (LostEnergy > 0) {
std::cout << "" << LostEnergy << " of KE went nowhere..." << std::endl;
}
#endif
}
diff --git a/src/Smearceptance/MetaSimpleSmearcepter.cxx b/src/Smearceptance/MetaSimpleSmearcepter.cxx
index 32f1e55..eac85e9 100644
--- a/src/Smearceptance/MetaSimpleSmearcepter.cxx
+++ b/src/Smearceptance/MetaSimpleSmearcepter.cxx
@@ -1,73 +1,74 @@
// 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 "EfficiencyApplicator.h"
#include "GaussianSmearer.h"
#include "ThresholdAccepter.h"
#include "TrackedMomentumMatrixSmearer.h"
+#include "VisECoalescer.h"
#include "MetaSimpleSmearcepter.h"
void MetaSimpleSmearcepter::SpecifcSetup(nuiskey &nk) {
std::map<std::string, SmearceptionFactory_fcn> factories;
factories["ThresholdAccepter"] = &BuildSmearcepter<ThresholdAccepter>;
factories["EfficiencyApplicator"] = &BuildSmearcepter<EfficiencyApplicator>;
factories["GaussianSmearer"] = &BuildSmearcepter<GaussianSmearer>;
factories["VisECoalescer"] = &BuildSmearcepter<VisECoalescer>;
factories["TrackedMomentumMatrixSmearer"] =
&BuildSmearcepter<TrackedMomentumMatrixSmearer>;
std::vector<nuiskey> smearcepters = nk.GetListOfChildNodes();
for (size_t smear_it = 0; smear_it < smearcepters.size(); ++smear_it) {
std::string const &smearType = smearcepters[smear_it].GetElementName();
if (smearType == "EnergyShuffler") {
ES = new EnergyShuffler();
ES->Setup(smearcepters[smear_it]);
continue;
}
if (!factories.count(smearType)) {
ERROR(WRN, "No known smearer accepts elements named: \"" << smearType
<< "\"");
continue;
}
Smearcepters.push_back(factories[smearType](smearcepters[smear_it]));
QLOG(FIT, "MetaSimpleSmearcepter adopted child smearcepter: "
<< Smearcepters.back()->GetName()
<< " of type: " << Smearcepters.back()->GetElementName());
}
NSmearcepters = Smearcepters.size();
}
RecoInfo *MetaSimpleSmearcepter::Smearcept(FitEvent *fe) {
if (ES) {
ES->DoTheShuffle(fe);
}
RecoInfo *ri = NULL;
for (size_t sm_it = 0; sm_it < NSmearcepters; ++sm_it) {
if (!sm_it) {
ri = Smearcepters[sm_it]->Smearcept(fe);
} else {
Smearcepters[sm_it]->SmearRecoInfo(ri);
}
}
return ri;
}
diff --git a/src/Smearceptance/SmearceptanceUtils.cxx b/src/Smearceptance/SmearceptanceUtils.cxx
index 6dda8c9..bd00b75 100644
--- a/src/Smearceptance/SmearceptanceUtils.cxx
+++ b/src/Smearceptance/SmearceptanceUtils.cxx
@@ -1,278 +1,284 @@
// 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 "SmearceptanceUtils.h"
#include "TDecompSVD.h"
#include "FitLogger.h"
namespace SmearceptanceUtils {
double Smear1DProp(TH2D *mapping, double TrueProp, TRandom3 *rnjesus) {
bool myrand = false;
if (!rnjesus) {
rnjesus = new TRandom3();
myrand = true;
}
if (myrand) {
delete rnjesus;
}
THROW("NIMPLEMENTED");
return 0;
}
TVectorD SVDInverseSolve(TVectorD *inp, TMatrixD *mapping) {
TDecompSVD svd(*mapping);
bool ok;
TVectorD c_svd = svd.Solve(*inp, ok);
if (!ok) {
THROW("Failed to solve SVD matrix equation.");
}
return c_svd;
}
TVectorD SVDInverseSolve(TH1D *inp, TMatrixD *mapping) {
TVectorD inp_v = GetVector(inp);
return SVDInverseSolve(&inp_v, mapping);
}
TVectorD SVDInverseSolve(TH1D *inp, TH2D *mapping) {
TMatrixD mat(mapping->GetXaxis()->GetNbins(),
mapping->GetYaxis()->GetNbins());
for (Int_t xb_it = 0; xb_it < mapping->GetXaxis()->GetNbins(); ++xb_it) {
for (Int_t yb_it = 0; yb_it < mapping->GetYaxis()->GetNbins(); ++yb_it) {
mat[xb_it][yb_it] = mapping->GetBinContent(xb_it + 1, yb_it + 1);
}
}
return SVDInverseSolve(inp, &mat);
}
TH2D *SVDGetInverse(TH2D *mapping, int NToTruncate) {
TMatrixD mat = GetMatrix(mapping);
if (mat.GetNcols() > mat.GetNrows()) {
THROW("Trying to invert a " << mat.GetNrows() << "x" << mat.GetNcols()
<< " matrix.");
}
- TH2D *inverse = SwapXYTH2D(mapping);
+ TH2D *inverse = dynamic_cast<TH2D *>(mapping->Clone());
inverse->SetName("inverse");
inverse->Reset();
TDecompSVD svd(mat);
svd.Decompose();
if (NToTruncate) {
TVectorD Sig(svd.GetSig());
TMatrixD U(svd.GetU());
TMatrixD V(svd.GetV());
if (svd.GetV().TestBit(TMatrixD::kTransposed)) {
THROW("ARGHH");
}
TMatrixD V_T = V.Transpose(V);
TMatrixD Sig_TruncM(U.GetNrows(), V.GetNrows());
for (Int_t i = 0; i < U.GetNrows(); ++i) {
for (Int_t j = 0; j < V.GetNrows(); ++j) {
Sig_TruncM[i][j] =
((i != j) || (i >= (Sig.GetNrows() - NToTruncate))) ? 0 : Sig[i];
}
}
TMatrixD Trunc = U * Sig_TruncM * V_T;
svd.~TDecompSVD();
new (&svd) TDecompSVD(Trunc);
}
TMatrixD inv = svd.Invert();
if (fabs(inv[mapping->GetXaxis()->GetNbins() / 2]
[mapping->GetXaxis()->GetNbins() / 2] -
mat[mapping->GetXaxis()->GetNbins() / 2]
[mapping->GetXaxis()->GetNbins() / 2]) <
std::numeric_limits<double>::epsilon()) {
THROW("Failed to SVD invert matrix.");
}
for (Int_t xb_it = 0; xb_it < inverse->GetXaxis()->GetNbins(); ++xb_it) {
for (Int_t yb_it = 0; yb_it < inverse->GetYaxis()->GetNbins(); ++yb_it) {
inverse->SetBinContent(xb_it + 1, yb_it + 1, inv[yb_it][xb_it]);
}
}
return inverse;
}
void GetSVDDecomp(TH2D *mapping, TVectorD &Sig, TMatrixD &U, TMatrixD &V) {
TMatrixD mat = GetMatrix(mapping);
TDecompSVD svd(mat);
U.ResizeTo(svd.GetU());
U = svd.GetU();
V.ResizeTo(svd.GetV());
V = svd.GetU();
Sig.ResizeTo(svd.GetSig());
Sig = svd.GetSig();
}
TVectorD GetVector(TH1D *inp) {
TVectorD vec(inp->GetXaxis()->GetNbins());
for (Int_t xb_it = 0; xb_it < inp->GetXaxis()->GetNbins(); ++xb_it) {
vec[xb_it] = inp->GetBinContent(xb_it + 1);
}
return vec;
}
TVectorD GetErrorVector(TH1D *inp) {
TVectorD vec(inp->GetXaxis()->GetNbins());
for (Int_t xb_it = 0; xb_it < inp->GetXaxis()->GetNbins(); ++xb_it) {
vec[xb_it] = inp->GetBinError(xb_it + 1);
}
return vec;
}
TMatrixD GetMatrix(TH2D *inp) {
TMatrixD mat(inp->GetYaxis()->GetNbins(), inp->GetXaxis()->GetNbins());
for (Int_t xb_it = 0; xb_it < inp->GetXaxis()->GetNbins(); ++xb_it) {
for (Int_t yb_it = 0; yb_it < inp->GetYaxis()->GetNbins(); ++yb_it) {
mat[yb_it][xb_it] = inp->GetBinContent(xb_it + 1, yb_it + 1);
}
}
return mat;
}
TH1D *GetTH1FromVector(TVectorD const &inp, TH1D *templ) {
TH1D *hist;
if (templ) {
hist = static_cast<TH1D *>(templ->Clone());
hist->Reset();
hist->SetName("vectHist");
} else {
hist = new TH1D("vectHist", "", inp.GetNrows(), 0, inp.GetNrows());
}
+ hist->SetDirectory(NULL);
for (Int_t xb_it = 0; xb_it < inp.GetNrows(); ++xb_it) {
hist->SetBinContent(xb_it + 1, inp[xb_it]);
}
return hist;
}
TH2D *GetTH2FromMatrix(TMatrixD const &inp, TH2D *templ) {
TH2D *hist;
if (templ) {
hist = static_cast<TH2D *>(templ->Clone());
hist->Reset();
hist->SetName("matHist");
} else {
hist = new TH2D("matHist", "", inp.GetNcols(), 0, inp.GetNcols(),
inp.GetNrows(), 0, inp.GetNrows());
}
+ hist->SetDirectory(NULL);
+
for (Int_t xb_it = 0; xb_it < inp.GetNcols(); ++xb_it) {
for (Int_t yb_it = 0; yb_it < inp.GetNrows(); ++yb_it) {
hist->SetBinContent(xb_it + 1, yb_it + 1, inp[yb_it][xb_it]);
}
}
return hist;
}
TVectorD ThrowVectFromHist(TH1D *inp, TRandom3 *rnjesus, bool allowNeg) {
TVectorD vec(inp->GetXaxis()->GetNbins());
for (Int_t xb_it = 0; xb_it < inp->GetXaxis()->GetNbins(); ++xb_it) {
size_t attempt = 0;
do {
if (attempt > 1000) {
THROW("Looks like we aren't getting anywhere with this bin: "
<< inp->GetBinContent(xb_it + 1) << " +- "
<< inp->GetBinError(xb_it + 1));
}
vec[xb_it] = inp->GetBinContent(xb_it + 1) +
inp->GetBinError(xb_it + 1) * rnjesus->Gaus();
attempt++;
} while ((!allowNeg) && (vec[xb_it] < 0));
}
return vec;
}
void PushTH1ThroughMatrixWithErrors(TH1D *inp, TH1D *oup, TMatrixD &response,
size_t NToys, bool allowNeg) {
TRandom3 rnjesus;
oup->Reset();
TVectorD Mean(oup->GetXaxis()->GetNbins());
TVectorD RMS(oup->GetXaxis()->GetNbins());
std::vector<TVectorD> Toys;
Toys.reserve(NToys);
double NToysFact = 1.0 / double(NToys);
for (size_t t_it = 0; t_it < NToys; ++t_it) {
TVectorD Toy = ThrowVectFromHist(inp, &rnjesus, allowNeg);
TVectorD UnfoldToy = response * Toy;
for (Int_t bi_it = 0; bi_it < oup->GetXaxis()->GetNbins(); ++bi_it) {
Mean[bi_it] += UnfoldToy[bi_it] * NToysFact;
}
Toys.push_back(UnfoldToy);
}
for (size_t t_it = 0; t_it < NToys; ++t_it) {
for (Int_t bi_it = 0; bi_it < oup->GetXaxis()->GetNbins(); ++bi_it) {
RMS[bi_it] += (Mean[bi_it] - Toys[t_it][bi_it]) *
(Mean[bi_it] - Toys[t_it][bi_it]) * NToysFact;
}
}
for (Int_t bi_it = 0; bi_it < oup->GetXaxis()->GetNbins(); ++bi_it) {
oup->SetBinContent(bi_it + 1, Mean[bi_it]);
oup->SetBinError(bi_it + 1, sqrt(RMS[bi_it]));
}
}
-TH2D *SwapXYTH2D(TH2D *templ) {
+TH2D * SwapXYTH2D(TH2D *templ) {
TH2D *Swapped = new TH2D(
(std::string(templ->GetName()) + "_c").c_str(), "",
templ->GetYaxis()->GetNbins(), templ->GetYaxis()->GetXbins()->GetArray(),
templ->GetXaxis()->GetNbins(), templ->GetXaxis()->GetXbins()->GetArray());
Swapped->Reset();
+ Swapped->SetDirectory(NULL);
+
+
std::string title = ";";
title += templ->GetYaxis()->GetTitle();
title += ";";
title += templ->GetXaxis()->GetTitle();
Swapped->SetTitle(title.c_str());
for (Int_t x_it = 0; x_it < templ->GetXaxis()->GetNbins() + 2; ++x_it) {
for (Int_t y_it = 0; y_it < templ->GetYaxis()->GetNbins() + 2; ++y_it) {
Swapped->SetBinContent(y_it, x_it, templ->GetBinContent(x_it, y_it));
Swapped->SetBinError(y_it, x_it, templ->GetBinError(x_it, y_it));
}
}
return Swapped;
}
}
diff --git a/src/Smearceptance/Smearcepterton.cxx b/src/Smearceptance/Smearcepterton.cxx
index 55fdbe7..17ff95d 100644
--- a/src/Smearceptance/Smearcepterton.cxx
+++ b/src/Smearceptance/Smearcepterton.cxx
@@ -1,77 +1,78 @@
// 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 "Smearcepterton.h"
-#include "ThresholdAccepter.h"
#include "EfficiencyApplicator.h"
#include "GaussianSmearer.h"
-#include "TrackedMomentumMatrixSmearer.h"
#include "MetaSimpleSmearcepter.h"
+#include "ThresholdAccepter.h"
+#include "TrackedMomentumMatrixSmearer.h"
+#include "VisECoalescer.h"
#include <vector>
Smearcepterton* Smearcepterton::_inst = NULL;
Smearcepterton& Smearcepterton::Get() {
if (!_inst) {
_inst = new Smearcepterton();
}
return *_inst;
}
Smearcepterton::Smearcepterton() { InitialiserSmearcepters(); }
void Smearcepterton::InitialiserSmearcepters() {
// hard coded list of tag name -> smearcepter factories, add here to add your
// own.
std::map<std::string, SmearceptionFactory_fcn> factories;
factories["ThresholdAccepter"] = &BuildSmearcepter<ThresholdAccepter>;
factories["EfficiencyApplicator"] = &BuildSmearcepter<EfficiencyApplicator>;
factories["GaussianSmearer"] = &BuildSmearcepter<GaussianSmearer>;
- factories["TrackedMomentumMatrixSmearer"] = &BuildSmearcepter<TrackedMomentumMatrixSmearer>;
+ factories["TrackedMomentumMatrixSmearer"] =
+ &BuildSmearcepter<TrackedMomentumMatrixSmearer>;
factories["VisECoalescer"] = &BuildSmearcepter<VisECoalescer>;
factories["MetaSimpleSmearcepter"] = &BuildSmearcepter<MetaSimpleSmearcepter>;
-
std::vector<nuiskey> smearcepterBlocks = Config::QueryKeys("smearcepters");
for (size_t smearB_it = 0; smearB_it < smearcepterBlocks.size();
++smearB_it) {
std::vector<nuiskey> smearcepters =
smearcepterBlocks[smearB_it].GetListOfChildNodes();
for (size_t smear_it = 0; smear_it < smearcepters.size(); ++smear_it) {
std::string const& smearType = smearcepters[smear_it].GetElementName();
if (!factories.count(smearType)) {
ERROR(WRN, "No known smearer accepts elements named: \"" << smearType
<< "\"");
continue;
}
ISmearcepter* smearer = factories[smearType](smearcepters[smear_it]);
Smearcepters[smearer->GetName()] = smearer;
QLOG(FIT, "Configured smearer named: " << smearer->GetName()
<< " of type: "
<< smearer->GetElementName());
}
}
}
diff --git a/src/Smearceptance/ThresholdAccepter.cxx b/src/Smearceptance/ThresholdAccepter.cxx
index 4a4dff3..adc7239 100644
--- a/src/Smearceptance/ThresholdAccepter.cxx
+++ b/src/Smearceptance/ThresholdAccepter.cxx
@@ -1,330 +1,330 @@
// 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 "ThresholdAccepter.h"
namespace {
ThresholdAccepter::KineVar GetKineType(nuiskey &nk) {
if (nk.Has("RecoThresholdMomentum_MeV")) {
return ThresholdAccepter::kMomentum;
} else if (nk.Has("RecoThresholdKE_MeV")) {
return ThresholdAccepter::kKE;
} else if (nk.Has("RecoThresholdCosTheta_Max")) {
return ThresholdAccepter::kCosTheta_Max;
} else if (nk.Has("RecoThresholdCosTheta_Min")) {
return ThresholdAccepter::kCosTheta_Min;
} else if (nk.Has("RecoThresholdAbsCosTheta_Max")) {
return ThresholdAccepter::kAbsCosTheta_Max;
} else if (nk.Has("RecoThresholdAbsCosTheta_Min")) {
return ThresholdAccepter::kAbsCosTheta_Min;
} else {
THROW("Cannot determine the threshold type for Smearcepter element.");
}
return ThresholdAccepter::kNoVar;
}
std::string GetKineTypeName(ThresholdAccepter::KineVar kv) {
switch (kv) {
case ThresholdAccepter::kMomentum:
return "Momentum";
case ThresholdAccepter::kKE:
return "KE";
case ThresholdAccepter::kCosTheta_Max:
return "CosTheta_Max";
case ThresholdAccepter::kCosTheta_Min:
return "CosTheta_Min";
case ThresholdAccepter::kAbsCosTheta_Max:
return "AbsCosTheta_Max";
case ThresholdAccepter::kAbsCosTheta_Min:
return "CosTheta_Min";
default:
return "NoVar";
}
}
double GetKineThreshold(nuiskey &nk, ThresholdAccepter::KineVar kv) {
switch (kv) {
case ThresholdAccepter::kMomentum:
return nk.GetD("RecoThresholdMomentum_MeV");
case ThresholdAccepter::kKE:
return nk.GetD("RecoThresholdKE_MeV");
case ThresholdAccepter::kCosTheta_Max:
return nk.GetD("RecoThresholdCosTheta_Max");
case ThresholdAccepter::kCosTheta_Min:
return nk.GetD("RecoThresholdCosTheta_Min");
case ThresholdAccepter::kAbsCosTheta_Max:
return nk.GetD("RecoThresholdAbsCosTheta_Max");
case ThresholdAccepter::kAbsCosTheta_Min:
return nk.GetD("RecoThresholdAbsCosTheta_Min");
default:
return 0;
}
}
double GetKineVal(FitParticle *fp, ThresholdAccepter::Thresh &rt) {
switch (rt.ThresholdType) {
case ThresholdAccepter::kMomentum:
return fp->P3().Mag();
case ThresholdAccepter::kKE:
return fp->KE();
case ThresholdAccepter::kCosTheta_Max:
return fp->P3().CosTheta();
case ThresholdAccepter::kCosTheta_Min:
return fp->P3().CosTheta();
case ThresholdAccepter::kAbsCosTheta_Max:
return fabs(fp->P3().CosTheta());
case ThresholdAccepter::kAbsCosTheta_Min:
return fabs(fp->P3().CosTheta());
default:
return 0;
}
}
bool PassesThreshold(FitParticle *fp, ThresholdAccepter::Thresh &rt) {
switch (rt.ThresholdType) {
case ThresholdAccepter::kMomentum:
return (fp->P3().Mag() > rt.ThresholdVal);
case ThresholdAccepter::kKE:
return (fp->KE() > rt.ThresholdVal);
case ThresholdAccepter::kCosTheta_Max:
return (fp->P3().CosTheta() < rt.ThresholdVal);
case ThresholdAccepter::kCosTheta_Min:
return (fp->P3().CosTheta() > rt.ThresholdVal);
case ThresholdAccepter::kAbsCosTheta_Max:
return (fabs(fp->P3().CosTheta()) < rt.ThresholdVal);
case ThresholdAccepter::kAbsCosTheta_Min:
return (fabs(fp->P3().CosTheta()) > rt.ThresholdVal);
default:
return 0;
}
}
}
/// Reads particle threshold nodes
///
/// Nodes look like:
/// <ThresholdAccepter Name="D00N_ND_LAr">
/// <RecoThreshold PDG="211,-211" RecoThresholdKE_MeV="30" />
/// <RecoThreshold PDG="211,-211" RecoThresholdCosTheta_Max="1" />
/// <RecoThreshold PDG="2212" RecoThresholdMomentum_MeV="350" />
/// <RecoThreshold PDG="2212" RecoThresholdAbsCosTheta_Min="0" />
/// <VisThreshold PDG="2212" VisThresholdKE_MeV="10" Contrib="K" />
/// <VisThreshold PDG="22" VisThresholdKE_MeV="10" Contrib="T" />
/// </ThresholdAccepter>
void ThresholdAccepter::SpecifcSetup(nuiskey &nk) {
std::vector<nuiskey> recoThresholdDescriptors =
nk.GetListOfChildNodes("RecoThreshold");
for (size_t t_it = 0; t_it < recoThresholdDescriptors.size(); ++t_it) {
std::string pdgs_s = recoThresholdDescriptors[t_it].GetS("PDG");
std::vector<int> pdgs_i = GeneralUtils::ParseToInt(pdgs_s, ",");
for (size_t pdg_it = 0; pdg_it < pdgs_i.size(); ++pdg_it) {
Thresh t;
t.ThresholdType = GetKineType(recoThresholdDescriptors[t_it]);
t.ThresholdVal =
GetKineThreshold(recoThresholdDescriptors[t_it], t.ThresholdType);
ReconThresholds[pdgs_i[pdg_it]].push_back(t);
- QLOG(SAM, "Added reconstruction threshold of type: "
+ QLOG(FIT, "Added reconstruction threshold of type: "
<< ReconThresholds[pdgs_i[pdg_it]].back().ThresholdVal
<< " "
<< GetKineTypeName(
ReconThresholds[pdgs_i[pdg_it]].back().ThresholdType)
<< ", for PDG: " << pdgs_i[pdg_it]);
}
}
std::vector<nuiskey> visThresholdDescriptors =
nk.GetListOfChildNodes("VisThreshold");
for (size_t t_it = 0; t_it < visThresholdDescriptors.size(); ++t_it) {
std::string pdgs_s = visThresholdDescriptors[t_it].GetS("PDG");
std::vector<int> pdgs_i = GeneralUtils::ParseToInt(pdgs_s, ",");
for (size_t pdg_it = 0; pdg_it < pdgs_i.size(); ++pdg_it) {
if (VisThresholds.count(pdgs_i[pdg_it])) {
ERROR(WRN, "Smearceptor " << ElementName << ":" << InstanceName
<< " already has a threshold for PDG: "
<< pdgs_i[pdg_it]);
}
VisThresh vt;
vt.UseKE = visThresholdDescriptors[t_it].Has("Contrib")
? (visThresholdDescriptors[t_it].GetS("Contrib") == "K")
: false;
vt.Fraction = visThresholdDescriptors[t_it].Has("Fraction")
? visThresholdDescriptors[t_it].GetD("Fraction")
: 1;
if (visThresholdDescriptors[t_it].Has("VisThresholdKE_MeV")) {
vt.ThresholdType = ThresholdAccepter::kKE;
vt.ThresholdVal =
visThresholdDescriptors[t_it].GetD("VisThresholdKE_MeV");
} else if (visThresholdDescriptors[t_it].Has(
"VisThresholdMomentum_MeV")) {
vt.ThresholdType = ThresholdAccepter::kMomentum;
vt.ThresholdVal =
visThresholdDescriptors[t_it].GetD("VisThresholdMomentum_MeV");
;
} else {
ERROR(WRN, "Smearceptor "
<< ElementName << ":" << InstanceName
<< " cannot find threshold information for PDG: "
<< pdgs_i[pdg_it]);
continue;
}
VisThresholds[pdgs_i[pdg_it]] = vt;
- QLOG(SAM,
+ QLOG(FIT,
"Added visibility threshold of MeV "
<< VisThresholds[pdgs_i[pdg_it]].ThresholdVal << " "
<< GetKineTypeName(VisThresholds[pdgs_i[pdg_it]].ThresholdType)
<< ", for PDG: " << pdgs_i[pdg_it]
<< ". If visible, particle deposits: "
<< (VisThresholds[pdgs_i[pdg_it]].UseKE ? "KE" : "TE"));
}
}
}
void ThresholdAccepter::SmearceptOneParticle(RecoInfo *ri, FitParticle *fp
#ifdef DEBUG_THRESACCEPT
,
size_t p_it
#endif
) {
#ifdef DEBUG_THRESACCEPT
std::cout << std::endl;
std::cout << "[" << p_it << " = " << fp << "]: " << fp->PDG() << ", "
<< fp->Status() << ", " << fp->E() << " -- KE:" << fp->KE()
<< " Mom: " << fp->P3().Mag() << std::flush;
#endif
if (fp->Status() != kFinalState) {
#ifdef DEBUG_THRESACCEPT
std::cout << " -- Not final state." << std::flush;
#endif
return;
}
if ((ReconThresholds.count(fp->PDG()) + VisThresholds.count(fp->PDG())) ==
0) {
#ifdef DEBUG_THRESACCEPT
std::cout << " -- Undetectable." << std::flush;
#endif
return;
}
// If no reco thresholds it should fall through to EVis
bool Passes = ReconThresholds[fp->PDG()].size();
bool FailEnergyThresh = !ReconThresholds[fp->PDG()].size();
for (size_t rt_it = 0; rt_it < ReconThresholds[fp->PDG()].size(); ++rt_it) {
bool Passed = PassesThreshold(fp, ReconThresholds[fp->PDG()][rt_it]);
if (!Passed) {
#ifdef DEBUG_THRESACCEPT
std::cout << "\n\t -- Rejected. ("
<< GetKineTypeName(
ReconThresholds[fp->PDG()][rt_it].ThresholdType)
<< " Threshold: "
<< ReconThresholds[fp->PDG()][rt_it].ThresholdVal << " | "
<< GetKineVal(fp, ReconThresholds[fp->PDG()][rt_it]) << ")"
<< std::flush;
if ((ReconThresholds[fp->PDG()][rt_it].ThresholdType ==
ThresholdAccepter::kMomentum) ||
(ReconThresholds[fp->PDG()][rt_it].ThresholdType ==
ThresholdAccepter::kKE)) {
FailEnergyThresh = true;
}
#endif
} else {
#ifdef DEBUG_THRESACCEPT
std::cout << "\n\t -- Accepted. ("
<< GetKineTypeName(
ReconThresholds[fp->PDG()][rt_it].ThresholdType)
<< " Threshold: "
<< ReconThresholds[fp->PDG()][rt_it].ThresholdVal << " | "
<< GetKineVal(fp, ReconThresholds[fp->PDG()][rt_it]) << ")"
<< std::flush;
#endif
}
Passes = Passes && Passed;
}
if (Passes) {
#ifdef DEBUG_THRESACCEPT
std::cout << " -- Reconstructed." << std::flush;
#endif
ri->RecObjMom.push_back(fp->P3());
ri->RecObjClass.push_back(fp->PDG());
return;
} else if (!FailEnergyThresh) {
#ifdef DEBUG_THRESACCEPT
std::cout << " -- Failed non-Energy threshold, no chance for EVis."
<< std::flush;
#endif
return;
}
if (((VisThresholds[fp->PDG()].ThresholdType == ThresholdAccepter::kKE) &&
(VisThresholds[fp->PDG()].ThresholdVal <
fp->KE())) // Above KE-style threshold
|| ((VisThresholds[fp->PDG()].ThresholdType ==
ThresholdAccepter::kMomentum) &&
(VisThresholds[fp->PDG()].ThresholdVal <
fp->P3().Mag())) // Above mom-style threshold
) {
#ifdef DEBUG_THRESACCEPT
std::cout << " -- Contributed to VisE. ("
<< GetKineTypeName(VisThresholds[fp->PDG()].ThresholdType) << ": "
<< VisThresholds[fp->PDG()].ThresholdVal << ")" << std::flush;
#endif
ri->RecVisibleEnergy.push_back(
VisThresholds[fp->PDG()].Fraction *
(VisThresholds[fp->PDG()].UseKE ? fp->KE() : fp->E()));
ri->TrueContribPDGs.push_back(fp->PDG());
return;
} else {
#ifdef DEBUG_THRESACCEPT
std::cout << " -- Rejected. "
<< " Vis: ("
<< GetKineTypeName(VisThresholds[fp->PDG()].ThresholdType) << ": "
<< VisThresholds[fp->PDG()].ThresholdVal << ")" << std::flush;
#endif
}
}
RecoInfo *ThresholdAccepter::Smearcept(FitEvent *fe) {
RecoInfo *ri = new RecoInfo();
for (size_t p_it = 0; p_it < fe->NParticles(); ++p_it) {
FitParticle *fp = fe->GetParticle(p_it);
SmearceptOneParticle(ri, fp
#ifdef DEBUG_THRESACCEPT
,
p_it
#endif
);
}
#ifdef DEBUG_THRESACCEPT
std::cout << std::endl;
#endif
return ri;
}
diff --git a/src/Smearceptance/TrackedMomentumMatrixSmearer.cxx b/src/Smearceptance/TrackedMomentumMatrixSmearer.cxx
index 25f1f5f..9dda0cc 100644
--- a/src/Smearceptance/TrackedMomentumMatrixSmearer.cxx
+++ b/src/Smearceptance/TrackedMomentumMatrixSmearer.cxx
@@ -1,431 +1,411 @@
// 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 "TrackedMomentumMatrixSmearer.h"
namespace {
TrackedMomentumMatrixSmearer::DependVar GetVarType(std::string const &axisvar) {
if (axisvar == "Momentum") {
return TrackedMomentumMatrixSmearer::kMomentum;
} else if (axisvar == "KE") {
return TrackedMomentumMatrixSmearer::kKE;
} else if (axisvar == "TE") {
return TrackedMomentumMatrixSmearer::kTE;
}
return TrackedMomentumMatrixSmearer::kNoVar;
}
}
TH1D const *TrackedMomentumMatrixSmearer::SmearMap::GetRecoSlice(double val) {
if ((val < RecoSlices.front().first.first) ||
(val > RecoSlices.back().first.second)) {
ERROR(WRN,
"Kinematic property: " << val << ", not within smearable range: ["
<< RecoSlices.front().first.first << " -- "
<< RecoSlices.back().first.second << "].");
return NULL;
}
int L = 0, U = RecoSlices.size();
-#ifdef DEBUG_MATSMEAR
- std::cout << "\nBin search for " << val
- << " within: " << RecoSlices.front().first.first << " -- "
- << RecoSlices.back().first.second << std::endl;
-#endif
+
while (true) {
-#ifdef DEBUG_MATSMEAR
- std::cout << "\t\t U:" << U << ", L: " << L << " (/" << RecoSlices.size()
- << ")" << std::endl;
-#endif
if (U == L) {
-#ifdef DEBUG_MATSMEAR
- std::cout << "Found: " << L << std::endl;
-#endif
return RecoSlices[L].second;
}
int R = (U - L);
int m = L + (R / 2);
-#ifdef DEBUG_MATSMEAR
- std::cout << "Bin: " << m << "[ " << RecoSlices[m].first.first << " -- "
- << RecoSlices[m].first.second << "] (Search: " << val << ")"
- << std::endl;
-#endif
-
if (val <= RecoSlices[m].first.first) {
U = m - 1;
continue;
}
if (val > RecoSlices[m].first.second) {
L = m + 1;
continue;
}
if ((val > RecoSlices[m].first.first) &&
(val <= RecoSlices[m].first.second)) {
-#ifdef DEBUG_MATSMEAR
- std::cout << "Found: " << m << std::endl;
-#endif
return RecoSlices[m].second;
}
THROW("Binary smearing search failed. Check logic.");
}
}
TH1D *GetMapSlice(TH2D *mp, int SliceBin, bool AlongX) {
int NBins = (AlongX ? mp->GetXaxis() : mp->GetYaxis())->GetNbins();
int NOtherBins = (AlongX ? mp->GetYaxis() : mp->GetXaxis())->GetNbins();
if (SliceBin >= NOtherBins) {
THROW("Asked for slice " << SliceBin << " but the " << (AlongX ? 'Y' : 'X')
<< " axis only has " << NOtherBins);
}
if ((AlongX ? mp->GetXaxis() : mp->GetYaxis())->IsVariableBinSize() &&
((NBins + 1) !=
(AlongX ? mp->GetXaxis() : mp->GetYaxis())->GetXbins()->GetSize())) {
THROW(
"Attemping to take binning slice of variable binning, but NBins+1 != "
"NBinEdges: "
<< (NBins + 1) << " != "
<< (AlongX ? mp->GetXaxis() : mp->GetYaxis())->GetXbins()->GetSize());
}
std::stringstream ss("");
ss << mp->GetName() << (AlongX ? 'Y' : 'X') << "Slice_" << SliceBin;
std::stringstream st("");
st << mp->GetTitle() << ";"
<< (AlongX ? mp->GetXaxis() : mp->GetYaxis())->GetTitle() << ";"
<< "Count";
TH1D *Ret;
if ((AlongX ? mp->GetXaxis() : mp->GetYaxis())->IsVariableBinSize()) {
Ret = new TH1D(
ss.str().c_str(), st.str().c_str(), NBins,
(AlongX ? mp->GetXaxis() : mp->GetYaxis())->GetXbins()->GetArray());
} else {
Ret = new TH1D(ss.str().c_str(), st.str().c_str(), NBins,
(AlongX ? mp->GetXaxis() : mp->GetYaxis())->GetXmin(),
(AlongX ? mp->GetXaxis() : mp->GetYaxis())->GetXmax());
}
for (int bi_it = 0; bi_it < NBins + 2; ++bi_it) {
int X = AlongX ? bi_it : SliceBin + 1;
int Y = AlongX ? SliceBin + 1 : bi_it;
int GBin = mp->GetBin(X, Y);
Ret->SetBinContent(bi_it, mp->GetBinContent(GBin));
Ret->SetBinError(bi_it, mp->GetBinError(GBin));
}
#ifdef DEBUG_MATSMEAR
std::cout << "Took slice: " << SliceBin << " spanning ["
<< Ret->GetXaxis()->GetBinLowEdge(1) << " -- "
<< Ret->GetXaxis()->GetBinUpEdge(Ret->GetXaxis()->GetNbins())
<< "] (Orignal span ["
<< (AlongX ? mp->GetXaxis() : mp->GetYaxis())->GetBinLowEdge(1)
<< " -- "
<< (AlongX ? mp->GetXaxis() : mp->GetYaxis())
->GetBinUpEdge(
(AlongX ? mp->GetXaxis() : mp->GetYaxis())->GetNbins())
<< "]) " << (AlongX ? mp->GetXaxis() : mp->GetYaxis())->GetXmin()
<< " -- " << (AlongX ? mp->GetXaxis() : mp->GetYaxis())->GetXmax()
<< " IsVBin: "
<< (AlongX ? mp->GetXaxis() : mp->GetYaxis())->IsVariableBinSize()
<< std::endl;
#endif
Ret->SetDirectory(NULL);
return Ret;
}
void TrackedMomentumMatrixSmearer::SmearMap::SetSlicesFromMap(TH2D *map,
bool TruthIsY) {
int NSlices = (TruthIsY ? map->GetYaxis() : map->GetXaxis())->GetNbins();
for (Int_t TrueSlice_it = 0; TrueSlice_it < NSlices; ++TrueSlice_it) {
std::pair<double, double> BinEdges;
BinEdges.first = (TruthIsY ? map->GetYaxis() : map->GetXaxis())
->GetBinLowEdge(TrueSlice_it + 1);
BinEdges.second = (TruthIsY ? map->GetYaxis() : map->GetXaxis())
->GetBinUpEdge(TrueSlice_it + 1);
TH1D *slice = GetMapSlice(map, TrueSlice_it, TruthIsY);
RecoSlices.push_back(std::make_pair(BinEdges, slice));
}
QLOG(FIT, "\tAdded " << RecoSlices.size() << " reco slices.");
}
/// Reads particle efficiency nodes
///
/// Nodes look like:
/// <TrackedMomentumMatrixSmearer Name="D00N_ND_LAr">
/// <SmearMatrix PDG="211,-211" InputFile="smear.root"
/// HistName="cpion_mom_smear" Kinematics="[KE|Momentum]"
/// MatrixToInternal="1E3" YIsTrue="true"/>
/// </TrackedMomentumMatrixSmearer>
void TrackedMomentumMatrixSmearer::SpecifcSetup(nuiskey &nk) {
std::vector<nuiskey> effDescriptors = nk.GetListOfChildNodes("SmearMatrix");
for (size_t t_it = 0; t_it < effDescriptors.size(); ++t_it) {
std::string inputFileName = effDescriptors[t_it].GetS("InputFile");
std::string HistName = effDescriptors[t_it].GetS("HistName");
bool YIsTrue = effDescriptors[t_it].Has("YIsTrue")
? effDescriptors[t_it].GetI("YIsTrue")
: true;
double UnitsScale = effDescriptors[t_it].Has("MatrixToInternal")
? effDescriptors[t_it].GetD("MatrixToInternal")
: 1;
TFile inputFile(inputFileName.c_str());
if (!inputFile.IsOpen()) {
THROW("Couldn't open specified input root file: " << inputFileName);
}
TH2D *inpHist = dynamic_cast<TH2D *>(inputFile.Get(HistName.c_str()));
if (!inpHist) {
THROW("Couldn't get TH2D named: " << HistName << " from input root file: "
<< inputFileName);
}
TrackedMomentumMatrixSmearer::DependVar var =
GetVarType(effDescriptors[t_it].GetS("Kinematics"));
std::string pdgs_s = effDescriptors[t_it].GetS("PDG");
std::vector<int> pdgs_i = GeneralUtils::ParseToInt(pdgs_s, ",");
for (size_t pdg_it = 0; pdg_it < pdgs_i.size(); ++pdg_it) {
if (ParticleMappings.count(pdgs_i[pdg_it])) {
ERROR(WRN, "Smearceptor " << ElementName << ":" << InstanceName
<< " already has a smearing for PDG: "
<< pdgs_i[pdg_it]);
}
SmearMap sm;
sm.SetSlicesFromMap(inpHist, YIsTrue);
sm.SmearVar = var;
sm.UnitsScale = UnitsScale;
ParticleMappings[pdgs_i[pdg_it]] = sm;
- QLOG(SAM, "Added smearing map for PDG: " << pdgs_i[pdg_it]);
+ QLOG(FIT, "Added smearing map for PDG: " << pdgs_i[pdg_it]);
}
}
SlaveGS.Setup(nk);
}
RecoInfo *TrackedMomentumMatrixSmearer::Smearcept(FitEvent *fe) {
RecoInfo *ri = new RecoInfo();
for (size_t p_it = 0; p_it < fe->NParticles(); ++p_it) {
FitParticle *fp = fe->GetParticle(p_it);
#ifdef DEBUG_MATSMEAR
std::cout << std::endl;
std::cout << "[" << p_it << "]: " << fp->PDG() << ", " << fp->Status()
<< ", " << fp->E() << " -- KE:" << fp->KE()
<< " Mom: " << fp->P3().Mag() << std::flush;
#endif
if (fp->Status() != kFinalState) {
#ifdef DEBUG_MATSMEAR
std::cout << " -- Not final state." << std::flush;
#endif
continue;
}
if (!ParticleMappings.count(fp->PDG())) {
SlaveGS.SmearceptOneParticle(ri, fp
#ifdef DEBUG_GAUSSSMEAR
,
p_it
#endif
);
continue;
}
SmearMap &sm = ParticleMappings[fp->PDG()];
double kineProp = 0;
switch (sm.SmearVar) {
case kMomentum: {
kineProp = fp->P3().Mag();
break;
}
case kKE: {
kineProp = fp->KE();
break;
}
case kTE: {
kineProp = fp->E();
break;
}
default: { THROW("Trying to find particle value for a kNoAxis."); }
}
TH1 const *recoDistrib = sm.GetRecoSlice(kineProp / sm.UnitsScale);
#ifdef DEBUG_MATSMEAR
std::cout << " -- Got slice spanning ["
<< recoDistrib->GetXaxis()->GetBinLowEdge(1) << " -- "
<< recoDistrib->GetXaxis()->GetBinUpEdge(
recoDistrib->GetXaxis()->GetNbins())
<< "]" << std::endl;
#endif
if (!recoDistrib) {
#ifdef DEBUG_MATSMEAR
std::cout << " -- outside smearable range." << std::flush;
#endif
continue;
}
if (recoDistrib->Integral() == 0) {
ERROR(WRN, "True slice has no reconstructed events. Not smearing.")
continue;
}
double Smeared = recoDistrib->GetRandom() * sm.UnitsScale;
#ifdef DEBUG_MATSMEAR
std::cout << "GotRandom: " << Smeared << ", MPV: "
<< recoDistrib->GetXaxis()->GetBinCenter(
recoDistrib->GetMaximumBin()) *
sm.UnitsScale
<< std::endl;
#endif
switch (sm.SmearVar) {
case kMomentum: {
ri->RecObjMom.push_back(fp->P3().Unit() * Smeared);
#ifdef DEBUG_MATSMEAR
std::cout << " -- Smeared: " << fp->p() << " -> " << Smeared << "."
<< std::flush;
#endif
break;
}
case kKE: {
double mass = fp->P4().M();
double TE = mass + Smeared;
double magP = sqrt(TE * TE - mass * mass);
ri->RecObjMom.push_back(fp->P3().Unit() * magP);
#ifdef DEBUG_MATSMEAR
std::cout << " -- Smeared: " << fp->KE() << " (mass: " << mass
<< ") -> " << Smeared
<< ". Smear Mom: " << ri->RecObjMom.back().Mag() << "."
<< std::flush;
#endif
break;
}
case kTE: {
double mass = fp->P4().M();
double TE = Smeared;
double magP = sqrt(TE * TE - mass * mass);
ri->RecObjMom.push_back(fp->P3().Unit() * magP);
#ifdef DEBUG_MATSMEAR
std::cout << " -- Smeared: " << fp->E() << " (mass: " << mass
<< ") -> " << Smeared
<< ". Smear Mom: " << ri->RecObjMom.back().Mag() << "."
<< std::flush;
#endif
break;
}
default: {}
}
#ifdef DEBUG_MATSMEAR
std::cout << " -- momentum reconstructed as " << ri->RecObjMom.back().Mag()
<< "." << std::endl;
#endif
if (ri->RecObjMom.back().Mag() != ri->RecObjMom.back().Mag()) {
ERROR(WRN, "Invalid particle built.");
ri->RecObjMom.pop_back();
#include "TCanvas.h"
TCanvas *Test = new TCanvas("c1", "");
static_cast<TH1 *>(recoDistrib->Clone())->Draw();
Test->SaveAs("Fail.png");
delete Test;
THROW("ARGH");
} else {
ri->RecObjClass.push_back(fp->PDG());
}
}
#ifdef DEBUG_MATSMEAR
std::cout << std::endl;
#endif
return ri;
}
void TrackedMomentumMatrixSmearer::SmearRecoInfo(RecoInfo *ri) {
for (size_t p_it = 0; p_it < ri->RecObjMom.size(); ++p_it) {
if (!ParticleMappings.count(ri->RecObjClass[p_it])) {
SlaveGS.SmearceptOneParticle(ri->RecObjMom[p_it], ri->RecObjClass[p_it]);
continue;
}
SmearMap &sm = ParticleMappings[ri->RecObjClass[p_it]];
double kineProp = 0;
switch (sm.SmearVar) {
case kMomentum: {
kineProp = ri->RecObjMom[p_it].Mag();
break;
}
case kKE: {
double mass = PhysConst::GetMass(ri->RecObjClass[p_it]) * 1E3;
kineProp = sqrt(ri->RecObjMom[p_it].Mag2() + mass * mass) - mass;
break;
}
case kTE: {
double mass = PhysConst::GetMass(ri->RecObjClass[p_it]) * 1E3;
kineProp = sqrt(ri->RecObjMom[p_it].Mag2() + mass * mass);
break;
}
default: { THROW("Trying to find particle value for a kNoAxis."); }
}
TH1 const *recoDistrib = sm.GetRecoSlice(kineProp / sm.UnitsScale);
if (!recoDistrib) {
continue;
}
double Smeared = recoDistrib->GetRandom() * sm.UnitsScale;
switch (sm.SmearVar) {
case kMomentum: {
ri->RecObjMom[p_it] = ri->RecObjMom[p_it].Unit() * Smeared;
break;
}
case kKE: {
double mass = PhysConst::GetMass(ri->RecObjClass[p_it]) * 1E3;
double TE = mass + Smeared;
double magP = sqrt(TE * TE - mass * mass);
ri->RecObjMom[p_it] = ri->RecObjMom[p_it].Unit() * magP;
break;
}
case kTE: {
double mass = PhysConst::GetMass(ri->RecObjClass[p_it]) * 1E3;
double TE = Smeared;
double magP = sqrt(TE * TE - mass * mass);
ri->RecObjMom[p_it] = ri->RecObjMom[p_it].Unit() * magP;
break;
}
default: {}
}
}
}
diff --git a/src/Utils/FitLogger.cxx b/src/Utils/FitLogger.cxx
index 6819b20..a5d364f 100644
--- a/src/Utils/FitLogger.cxx
+++ b/src/Utils/FitLogger.cxx
@@ -1,310 +1,350 @@
// 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 "FitLogger.h"
#include <fcntl.h>
#include <unistd.h>
-
namespace Logger {
// Logger Variables
-int log_verb = 4;
+int log_verb = 4;
bool use_colors = true;
-bool showtrace = false;
+bool showtrace = false;
std::ostream* __LOG_outstream(&std::cout);
std::ofstream __LOG_nullstream;
// Error Variables
int err_verb = 0;
std::ostream* __ERR_outstream(&std::cerr);
// Extra Variables
bool external_verb = false;
-bool super_rainbow_mode = true; //!< For when fitting gets boring.
+bool super_rainbow_mode = true; //!< For when fitting gets boring.
unsigned int super_rainbow_mode_colour = 0;
-std::streambuf *default_cout = std::cout.rdbuf();
-std::streambuf *default_cerr = std::cerr.rdbuf();
+std::streambuf* default_cout = std::cout.rdbuf();
+std::streambuf* default_cerr = std::cerr.rdbuf();
std::ofstream redirect_stream("/dev/null");
int silentfd = open("/dev/null", O_WRONLY);
int savedstdoutfd = dup(fileno(stdout));
int savedstderrfd = dup(fileno(stderr));
int nloggercalls = 0;
int timelastlog = 0;
-
}
-
-
// -------- Logging Functions --------- //
-
bool LOGGING(int level) {
- // std::cout << "LOGGING : " << __FILENAME__ << " " << __FUNCTION__ << std::endl;
- return (Logger::log_verb >= (uint)__GETLOG_LEVEL(level, __FILENAME__, __FUNCTION__));
+ // std::cout << "LOGGING : " << __FILENAME__ << " " << __FUNCTION__ <<
+ // std::endl;
+ return (Logger::log_verb >=
+ (uint)__GETLOG_LEVEL(level, __FILENAME__, __FUNCTION__));
};
int __GETLOG_LEVEL(int level, const char* filename, const char* funct) {
-
#ifdef __DEBUG__
int logfile = FitPar::Config().GetParI("logging." + std::string(filename));
if (logfile >= DEB and logfile <= EVT) {
level = logfile;
}
int logfunc = FitPar::Config().GetParI("logging." + std::string(funct));
if (logfunc >= DEB and logfunc <= EVT) {
level = logfunc;
}
#endif
return level;
};
-std::ostream& __OUTLOG(int level, const char* filename, const char* funct, int line) {
-
+std::ostream& __OUTLOG(int level, const char* filename, const char* funct,
+ int line) {
if (Logger::log_verb < (unsigned int)level &&
Logger::log_verb != (unsigned int)DEB) {
return (Logger::__LOG_nullstream);
} else {
-
if (Logger::use_colors) {
switch (level) {
- case FIT: std::cout << BOLDGREEN; break;
- case MIN: std::cout << BOLDBLUE; break;
- case SAM: std::cout << MAGENTA; break;
- case REC: std::cout << BLUE; break;
- case SIG: std::cout << GREEN; break;
- case DEB: std::cout << CYAN; break;
- default: break;
+ case FIT:
+ std::cout << BOLDGREEN;
+ break;
+ case MIN:
+ std::cout << BOLDBLUE;
+ break;
+ case SAM:
+ std::cout << MAGENTA;
+ break;
+ case REC:
+ std::cout << BLUE;
+ break;
+ case SIG:
+ std::cout << GREEN;
+ break;
+ case DEB:
+ std::cout << CYAN;
+ break;
+ default:
+ break;
}
}
switch (level) {
- case FIT: std::cout << "[LOG Fitter]"; break;
- case MIN: std::cout << "[LOG Minmzr]"; break;
- case SAM: std::cout << "[LOG Sample]"; break;
- case REC: std::cout << "[LOG Reconf]"; break;
- case SIG: std::cout << "[LOG Signal]"; break;
- case EVT: std::cout << "[LOG Event ]"; break;
- case DEB: std::cout << "[LOG DEBUG ]"; break;
- default: std::cout << "[LOG INFO ]"; break;
+ case FIT:
+ std::cout << "[LOG Fitter]";
+ break;
+ case MIN:
+ std::cout << "[LOG Minmzr]";
+ break;
+ case SAM:
+ std::cout << "[LOG Sample]";
+ break;
+ case REC:
+ std::cout << "[LOG Reconf]";
+ break;
+ case SIG:
+ std::cout << "[LOG Signal]";
+ break;
+ case EVT:
+ std::cout << "[LOG Event ]";
+ break;
+ case DEB:
+ std::cout << "[LOG DEBUG ]";
+ break;
+ default:
+ std::cout << "[LOG INFO ]";
+ break;
}
// Apply indent
if (true) {
switch (level) {
- case FIT: std::cout << ": "; break;
- case MIN: std::cout << ":- "; break;
- case SAM: std::cout << ":-- "; break;
- case REC: std::cout << ":--- "; break;
- case SIG: std::cout << ":---- "; break;
- case EVT: std::cout << ":----- "; break;
- case DEB: std::cout << ":------ "; break;
- default: std::cout << " "; break;
+ case FIT:
+ std::cout << ": ";
+ break;
+ case MIN:
+ std::cout << ":- ";
+ break;
+ case SAM:
+ std::cout << ":-- ";
+ break;
+ case REC:
+ std::cout << ":--- ";
+ break;
+ case SIG:
+ std::cout << ":---- ";
+ break;
+ case EVT:
+ std::cout << ":----- ";
+ break;
+ case DEB:
+ std::cout << ":------ ";
+ break;
+ default:
+ std::cout << " ";
+ break;
}
}
if (Logger::use_colors) std::cout << RESET;
if (Logger::showtrace) {
- std::cout << " : " << filename << "::" << funct << "[l. " << line << "] : ";
+ std::cout << " : " << filename << "::" << funct << "[l. " << line
+ << "] : ";
}
return *(Logger::__LOG_outstream);
}
}
-void SETVERBOSITY(int level) {
- Logger::log_verb = level;
-}
+void SETVERBOSITY(int level) { Logger::log_verb = level; }
void SETVERBOSITY(std::string verb) {
- if (!verb.compare("DEB")) Logger::log_verb = -1;
- else if (!verb.compare("QUIET")) Logger::log_verb = 0;
- else if (!verb.compare("FIT")) Logger::log_verb = 1;
- else if (!verb.compare("MIN")) Logger::log_verb = 2;
- else if (!verb.compare("SAM")) Logger::log_verb = 3;
- else if (!verb.compare("REC")) Logger::log_verb = 4;
- else if (!verb.compare("SIG")) Logger::log_verb = 5;
- else if (!verb.compare("EVT")) Logger::log_verb = 6;
- else Logger::log_verb = std::atoi(verb.c_str());
+ if (!verb.compare("DEB"))
+ Logger::log_verb = -1;
+ else if (!verb.compare("QUIET"))
+ Logger::log_verb = 0;
+ else if (!verb.compare("FIT"))
+ Logger::log_verb = 1;
+ else if (!verb.compare("MIN"))
+ Logger::log_verb = 2;
+ else if (!verb.compare("SAM"))
+ Logger::log_verb = 3;
+ else if (!verb.compare("REC"))
+ Logger::log_verb = 4;
+ else if (!verb.compare("SIG"))
+ Logger::log_verb = 5;
+ else if (!verb.compare("EVT"))
+ Logger::log_verb = 6;
+ else
+ Logger::log_verb = std::atoi(verb.c_str());
}
/// Set Trace Option
-void SETTRACE(bool val) {
- Logger::showtrace = val;
-}
+void SETTRACE(bool val) { Logger::showtrace = val; }
// ------ ERROR FUNCTIONS ---------- //
-std::ostream& __OUTERR(int level, const char* filename, const char* funct, int line) {
+std::ostream& __OUTERR(int level, const char* filename, const char* funct,
+ int line) {
if (Logger::use_colors) std::cerr << RED;
switch (level) {
- case FTL: std::cerr << "[ERR FATAL ]: "; break;
- case WRN: std::cerr << "[ERR WARN ]: "; break;
+ case FTL:
+ std::cerr << "[ERR FATAL ]: ";
+ break;
+ case WRN:
+ std::cerr << "[ERR WARN ]: ";
+ break;
}
if (Logger::use_colors) std::cerr << RESET;
// Allows enable error debugging trace
if (true or Logger::showtrace) {
std::cout << filename << "::" << funct << "[l. " << line << "] : ";
}
return *(Logger::__ERR_outstream);
}
// ----------- External Logging ----------- //
-void SETEXTERNALVERBOSITY(int level) {
- Logger::external_verb = (level > 0);
-}
+void SETEXTERNALVERBOSITY(int level) { Logger::external_verb = (level > 0); }
void StopTalking() {
-
// Check verbosity set correctly
if (!Logger::external_verb) return;
// Only redirect if we're not debugging
if (Logger::log_verb == (unsigned int)DEB) return;
std::cout.rdbuf(Logger::redirect_stream.rdbuf());
std::cerr.rdbuf(Logger::redirect_stream.rdbuf());
shhnuisancepythiaitokay_();
fflush(stdout);
fflush(stderr);
dup2(Logger::silentfd, fileno(stdout));
dup2(Logger::silentfd, fileno(stderr));
}
void StartTalking() {
-
// Check verbosity set correctly
if (!Logger::external_verb) return;
std::cout.rdbuf(Logger::default_cout);
std::cerr.rdbuf(Logger::default_cerr);
canihaznuisancepythia_();
fflush(stdout);
fflush(stderr);
dup2(Logger::savedstdoutfd, fileno(stdout));
dup2(Logger::savedstderrfd, fileno(stderr));
}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
//******************************************
void LOG_VERB(std::string verb) {
-//******************************************
-
- if (!verb.compare("DEB")) Logger::log_verb = -1;
- else if (!verb.compare("QUIET")) Logger::log_verb = 0;
- else if (!verb.compare("FIT")) Logger::log_verb = 1;
- else if (!verb.compare("MIN")) Logger::log_verb = 2;
- else if (!verb.compare("SAM")) Logger::log_verb = 3;
- else if (!verb.compare("REC")) Logger::log_verb = 4;
- else if (!verb.compare("SIG")) Logger::log_verb = 5;
- else if (!verb.compare("EVT")) Logger::log_verb = 6;
+ //******************************************
+
+ if (!verb.compare("DEB"))
+ Logger::log_verb = -1;
+ else if (!verb.compare("QUIET"))
+ Logger::log_verb = 0;
+ else if (!verb.compare("FIT"))
+ Logger::log_verb = 1;
+ else if (!verb.compare("MIN"))
+ Logger::log_verb = 2;
+ else if (!verb.compare("SAM"))
+ Logger::log_verb = 3;
+ else if (!verb.compare("REC"))
+ Logger::log_verb = 4;
+ else if (!verb.compare("SIG"))
+ Logger::log_verb = 5;
+ else if (!verb.compare("EVT"))
+ Logger::log_verb = 6;
// else Logger::log_verb = GeneralUtils::StrToInt(verb);
std::cout << "Set logging verbosity to : " << Logger::log_verb << std::endl;
return;
}
//******************************************
void ERR_VERB(std::string verb) {
-//******************************************
+ //******************************************
std::cout << "Setting ERROR VERB" << std::endl;
- if (!verb.compare("ERRQUIET")) Logger::err_verb = 0;
- else if (!verb.compare("FTL")) Logger::err_verb = 1;
- else if (!verb.compare("WRN")) Logger::err_verb = 2;
+ if (!verb.compare("ERRQUIET"))
+ Logger::err_verb = 0;
+ else if (!verb.compare("FTL"))
+ Logger::err_verb = 1;
+ else if (!verb.compare("WRN"))
+ Logger::err_verb = 2;
// else Logger::err_verb = GeneralUtils::StrToInt(verb);
std::cout << "Set error verbosity to : " << Logger::err_verb << std::endl;
return;
}
//******************************************
bool LOG_LEVEL(int level) {
-//******************************************
+ //******************************************
- if (Logger::log_verb == (unsigned int) DEB) {
+ if (Logger::log_verb == (unsigned int)DEB) {
return true;
}
- if (Logger::log_verb < (unsigned int) level) {
+ if (Logger::log_verb < (unsigned int)level) {
return false;
}
return true;
}
-void SET_TRACE(bool val) {
- Logger::showtrace = val;
-}
-
+void SET_TRACE(bool val) { Logger::showtrace = val; }
//******************************************
std::ostream& _LOG(int level, const char* filename, const char* func, int line)
//******************************************
{
return __OUTLOG(level, filename, func, line);
}
//******************************************
std::ostream& _ERR(int level, const char* filename, const char* func, int line)
//******************************************
{
-
if (Logger::use_colors) std::cerr << RED;
if (Logger::showtrace) {
std::cout << filename << "::" << func << "[l. " << line << "] : ";
}
switch (level) {
- case FTL: std::cerr << "[ERR FATAL ]: "; break;
- case WRN: std::cerr << "[ERR WARN ] : "; break;
+ case FTL:
+ std::cerr << "[ERR FATAL ]: ";
+ break;
+ case WRN:
+ std::cerr << "[ERR WARN ] : ";
+ break;
}
if (Logger::use_colors) std::cerr << RESET;
return *Logger::__ERR_outstream;
}
-
-
diff --git a/src/Utils/FitLogger.h b/src/Utils/FitLogger.h
index 6817bac..d8b26e9 100644
--- a/src/Utils/FitLogger.h
+++ b/src/Utils/FitLogger.h
@@ -1,207 +1,218 @@
// 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 FITLOGGER_HPP
#define FITLOGGER_HPP
/*!
* \addtogroup FitBase
* @{
*/
+#include <fstream>
#include <iosfwd>
#include <iostream>
-#include <fstream>
#include <sstream>
-#include "Initialiser.h"
#include "FitParameters.h"
+#include "Initialiser.h"
#include "TRandom3.h"
-#define RESET "\033[0m"
-#define BLACK "\033[30m" /* Black */
-#define RED "\033[31m" /* Red */
-#define GREEN "\033[32m" /* Green */
-#define YELLOW "\033[33m" /* Yellow */
-#define BLUE "\033[34m" /* Blue */
-#define MAGENTA "\033[35m" /* Magenta */
-#define CYAN "\033[36m" /* Cyan */
-#define WHITE "\033[37m" /* White */
-#define BOLDBLACK "\033[1m\033[30m" /* Bold Black */
-#define BOLDRED "\033[1m\033[31m" /* Bold Red */
-#define BOLDGREEN "\033[1m\033[32m" /* Bold Green */
-#define BOLDYELLOW "\033[1m\033[33m" /* Bold Yellow */
-#define BOLDBLUE "\033[1m\033[34m" /* Bold Blue */
-#define BOLDMAGENTA "\033[1m\033[35m" /* Bold Magenta */
-#define BOLDCYAN "\033[1m\033[36m" /* Bold Cyan */
-#define BOLDWHITE "\033[1m\033[37m" /* Bold White */
-
+#define RESET "\033[0m"
+#define BLACK "\033[30m" /* Black */
+#define RED "\033[31m" /* Red */
+#define GREEN "\033[32m" /* Green */
+#define YELLOW "\033[33m" /* Yellow */
+#define BLUE "\033[34m" /* Blue */
+#define MAGENTA "\033[35m" /* Magenta */
+#define CYAN "\033[36m" /* Cyan */
+#define WHITE "\033[37m" /* White */
+#define BOLDBLACK "\033[1m\033[30m" /* Bold Black */
+#define BOLDRED "\033[1m\033[31m" /* Bold Red */
+#define BOLDGREEN "\033[1m\033[32m" /* Bold Green */
+#define BOLDYELLOW "\033[1m\033[33m" /* Bold Yellow */
+#define BOLDBLUE "\033[1m\033[34m" /* Bold Blue */
+#define BOLDMAGENTA "\033[1m\033[35m" /* Bold Magenta */
+#define BOLDCYAN "\033[1m\033[36m" /* Bold Cyan */
+#define BOLDWHITE "\033[1m\033[37m" /* Bold White */
namespace Logger {
-extern int log_verb; //!< Current VERBOSITY
-extern int err_verb; //!< Current ERROR VERBOSITY
+extern int log_verb; //!< Current VERBOSITY
+extern int err_verb; //!< Current ERROR VERBOSITY
extern bool external_verb;
-extern bool use_colors; //!< Use BASH Terminal Colors Flag
-extern bool super_rainbow_mode; //!< For when fitting gets boring.
+extern bool use_colors; //!< Use BASH Terminal Colors Flag
+extern bool super_rainbow_mode; //!< For when fitting gets boring.
extern unsigned int super_rainbow_mode_colour;
-extern bool showtrace; // Quick Tracing for debugging
+extern bool showtrace; // Quick Tracing for debugging
extern int nloggercalls;
extern int timelastlog;
-extern std::streambuf *default_cout; //!< Where the STDOUT stream is currently directed
-extern std::streambuf *default_cerr; //!< Where the STDERR stream is currently directed
-extern std::ofstream redirect_stream; //!< Where should unwanted messages be thrown
+extern std::streambuf*
+ default_cout; //!< Where the STDOUT stream is currently directed
+extern std::streambuf*
+ default_cerr; //!< Where the STDERR stream is currently directed
+extern std::ofstream
+ redirect_stream; //!< Where should unwanted messages be thrown
}
/// Returns full path to file currently in
-#define __FILENAME__ (strrchr(__FILE__, '/') ? strrchr(__FILE__, '/') + 1 : __FILE__)
-
+#define __FILENAME__ \
+ (strrchr(__FILE__, '/') ? strrchr(__FILE__, '/') + 1 : __FILE__)
// ------ LOGGER FUNCTIONS ------------ //
namespace Logger {
- /// NULL Output Stream
- extern std::ofstream __LOG_nullstream;
+/// NULL Output Stream
+extern std::ofstream __LOG_nullstream;
- /// Logging Stream
- extern std::ostream* __LOG_outstream;
+/// Logging Stream
+extern std::ostream* __LOG_outstream;
}
/// Fitter VERBOSITY Enumerations
/// These go through the different depths of the fitter.
///
/// 0 QUIET - Little output.
/// 1 FIT - Top Level Minimizer Status
/// 2 MIN - Output from the FCN Minimizer Functions
/// 3 SAM - Output from each of the samples during setup etc
/// 4 REC - Output during each reconfigure. Percentage progress etc.
/// 5 SIG - Output during every signal event that is found.
/// 6 EVT - Output during every event.
-/// -1 DEB - Will print only debugging info wherever a LOG(DEB) statement was made
+/// -1 DEB - Will print only debugging info wherever a LOG(DEB) statement was
+/// made
enum __LOG_levels { DEB = -1, QUIET, FIT, MIN, SAM, REC, SIG, EVT };
/// Returns log level for a given file/function
int __GETLOG_LEVEL(int level, const char* filename, const char* funct);
/// Actually runs the logger
-std::ostream& __OUTLOG(int level, const char* filename, const char* funct, int line);
+std::ostream& __OUTLOG(int level, const char* filename, const char* funct,
+ int line);
/// Global Logging Definitions
-#define QLOG(level, stream) \
-{ \
-if (Logger::log_verb >= __GETLOG_LEVEL(level, __FILENAME__, __FUNCTION__)){ \
- __OUTLOG(level, __FILENAME__, __FUNCTION__, __LINE__) << stream << std::endl; \
-} \
-};
-
-#define BREAK(level) \
-{ \ \
- if (Logger::log_verb >= __GETLOG_LEVEL(level, __FILENAME__, __FUNCTION__)){ \
- __OUTLOG(level, __FILENAME__, __FUNCTION__, __LINE__) << std::endl; \
- }\
-};
-
+#define QLOG(level, stream) \
+ { \
+ if (Logger::log_verb >= \
+ __GETLOG_LEVEL(level, __FILENAME__, __FUNCTION__)) { \
+ __OUTLOG(level, __FILENAME__, __FUNCTION__, __LINE__) << stream \
+ << std::endl; \
+ } \
+ };
+
+#define BREAK(level) \
+ { \
+ \ if (Logger::log_verb >= \
+ __GETLOG_LEVEL(level, __FILENAME__, __FUNCTION__)) { \
+ __OUTLOG(level, __FILENAME__, __FUNCTION__, __LINE__) << std::endl; \
+ } \
+ };
/// Return whether logging level is valid
bool LOGGING(int level);
/// Set Global Verbosity
void SETVERBOSITY(int level);
/// Set Global Verbosity from String
void SETVERBOSITY(std::string verb);
/// Set Trace Option
void SETTRACE(bool val);
// ----------- ERROR FUNCTIONS ---------- //
/// Error Stream
extern std::ostream* __ERR_outstream;
/// Fitter ERROR VERBOSITY Enumerations
///
/// 0 QUIET - No Error Output
/// 1 FTL - Show errors only if fatal
/// 2 WRN - Show Warning messages
enum __ERR_levels { ERRQUIET = 0, FTL, WRN };
/// Actually runs the error messager
-std::ostream& __OUTERR(int level, const char* filename, const char* funct, int line);
+std::ostream& __OUTERR(int level, const char* filename, const char* funct,
+ int line);
/// Error Logging Function
-#define ERROR(level, stream) \
-{ \
- __OUTERR(level, __FILENAME__, __FUNCTION__, __LINE__) << stream << std::endl; \
-};
+#define ERROR(level, stream) \
+ { \
+ __OUTERR(level, __FILENAME__, __FUNCTION__, __LINE__) << stream \
+ << std::endl; \
+ };
// ----------- ERROR HANDLING ------------- //
/// Exit the program with given error message stream
-#define THROW(stream) \
-{ \
- __OUTERR(FTL, __FILENAME__, __FUNCTION__, __LINE__) << stream << std::endl; \
- __OUTERR(FTL, __FILENAME__, __FUNCTION__, __LINE__) << "Exiting!" << std::endl; \
- std::abort(); \
-}
-
+#define THROW(stream) \
+ { \
+ __OUTERR(FTL, __FILENAME__, __FUNCTION__, __LINE__) << stream \
+ << std::endl; \
+ __OUTERR(FTL, __FILENAME__, __FUNCTION__, __LINE__) \
+ << "Attempting to save output file." << std::endl; \
+ if (FitPar::Config().out && FitPar::Config().out->IsOpen()) { \
+ FitPar::Config().out->Write(); \
+ FitPar::Config().out->Close(); \
+ __OUTERR(FTL, __FILENAME__, __FUNCTION__, __LINE__) << "Done." \
+ << std::endl; \
+ } else { \
+ __OUTERR(FTL, __FILENAME__, __FUNCTION__, __LINE__) \
+ << "No output file set." << std::endl; \
+ } \
+ __OUTERR(FTL, __FILENAME__, __FUNCTION__, __LINE__) << "Exiting!" \
+ << std::endl; \
+ std::abort(); \
+ }
// ----------- External Logging ----------- //
void SETEXTERNALVERBOSITY(int level);
void StopTalking();
void StartTalking();
extern "C" {
- void shhnuisancepythiaitokay_(void);
- void canihaznuisancepythia_(void);
+void shhnuisancepythiaitokay_(void);
+void canihaznuisancepythia_(void);
}
-
-
-
-
-
-
-
-
-
// ---------- LEGACY FUNCTIONS -------------- //
bool LOG_LEVEL(int level);
//! Set LOG VERBOSITY from a string
void LOG_VERB(std::string verb);
inline void LOG_VERB(int verb) { Logger::log_verb = verb; };
void SET_TRACE(bool val);
//! Set ERROR VERBOSITY from a string
void ERR_VERB(std::string verb);
inline void ERR_VERB(int verb) { Logger::err_verb = verb; };
-
-/// Logging Function. Use as a string stream. e.g. LOG(SAM) << "This sample is dope." << std::endl;
-std::ostream& _LOG(int level, const char* filename, const char* funct, int line);
+/// Logging Function. Use as a string stream. e.g. LOG(SAM) << "This sample is
+/// dope." << std::endl;
+std::ostream& _LOG(int level, const char* filename, const char* funct,
+ int line);
#define LOG(level) _LOG(level, __FILENAME__, __FUNCTION__, __LINE__)
-
-//! Error Function. Use as a string stream. e.g. ERR(FTL) << "The fit is completely buggered." << std::endl;
-std::ostream& _ERR(int level, const char* filename, const char* funct, int line);
+//! Error Function. Use as a string stream. e.g. ERR(FTL) << "The fit is
+//! completely buggered." << std::endl;
+std::ostream& _ERR(int level, const char* filename, const char* funct,
+ int line);
#define ERR(level) _ERR(level, __FILENAME__, __FUNCTION__, __LINE__)
-
/*! @} */
#endif
-
diff --git a/src/Utils/GeneralUtils.cxx b/src/Utils/GeneralUtils.cxx
index c99e8a8..e0b10b7 100644
--- a/src/Utils/GeneralUtils.cxx
+++ b/src/Utils/GeneralUtils.cxx
@@ -1,196 +1,195 @@
// 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 "GeneralUtils.h"
std::string GeneralUtils::BoolToStr(bool val) {
std::ostringstream ss;
ss << val;
return ss.str();
}
std::string GeneralUtils::IntToStr(int val) {
std::ostringstream ss;
ss << val;
return ss.str();
};
std::string GeneralUtils::DblToStr(double val) {
std::ostringstream ss;
ss << val;
return ss.str();
};
-std::vector<std::string> GeneralUtils::LoadCharToVectStr(int argc, char* argv[]){
+std::vector<std::string> GeneralUtils::LoadCharToVectStr(int argc,
+ char* argv[]) {
std::vector<std::string> vect;
- for (int i = 1; i < argc; i++){
- vect.push_back( std::string(argv[i]) );
+ for (int i = 1; i < argc; i++) {
+ vect.push_back(std::string(argv[i]));
}
return vect;
}
-std::vector<std::string> GeneralUtils::ParseToStr(std::string str, const char* del) {
-
+std::vector<std::string> GeneralUtils::ParseToStr(std::string str,
+ const char* del) {
std::istringstream stream(str);
std::string temp_string;
std::vector<std::string> vals;
while (std::getline(stream >> std::ws, temp_string, *del)) {
if (temp_string.empty()) continue;
vals.push_back(temp_string);
}
return vals;
}
std::vector<double> GeneralUtils::ParseToDbl(std::string str, const char* del) {
-
std::istringstream stream(str);
std::string temp_string;
std::vector<double> vals;
while (std::getline(stream >> std::ws, temp_string, *del)) {
if (temp_string.empty()) continue;
std::istringstream stream(temp_string);
double entry;
stream >> entry;
vals.push_back(entry);
-
}
return vals;
}
std::vector<int> GeneralUtils::ParseToInt(std::string str, const char* del) {
-
std::istringstream stream(str);
std::string temp_string;
std::vector<int> vals;
while (std::getline(stream >> std::ws, temp_string, *del)) {
if (temp_string.empty()) continue;
std::istringstream stream(temp_string);
int entry;
stream >> entry;
vals.push_back(entry);
-
}
return vals;
}
double GeneralUtils::StrToDbl(std::string str) {
-
std::istringstream stream(str);
double val;
stream >> val;
return val;
}
int GeneralUtils::StrToInt(std::string str) {
-
std::istringstream stream(str);
int val;
stream >> val;
return val;
}
bool GeneralUtils::StrToBool(std::string str) {
-
// convert result to lower case
for (uint i = 0; i < str.size(); i++) str[i] = std::tolower(str[i]);
// Test for true/false
- if (!str.compare("false")) return false;
- else if (!str.compare("true") ) return true;
+ if (!str.compare("false"))
+ return false;
+ else if (!str.compare("true"))
+ return true;
if (str.empty()) return false;
// Push into bool
std::istringstream stream(str);
bool val;
stream >> val;
return val;
}
-std::vector<std::string> GeneralUtils::ParseFileToStr(std::string str, const char* del) {
-
+std::vector<std::string> GeneralUtils::ParseFileToStr(std::string str,
+ const char* del) {
std::vector<std::string> linevect;
std::string line;
std::ifstream read;
read.open(str.c_str());
if (!read.is_open()) {
- THROW("Cannot open file " << str << " in ParseFileToStr");
+ ERROR(FTL, "Cannot open file " << str << " in ParseFileToStr");
+ throw;
}
- while ( std::getline(read >> std::ws, line, *del) ) {
+ while (std::getline(read >> std::ws, line, *del)) {
linevect.push_back(line);
}
read.close();
return linevect;
}
std::string GeneralUtils::GetTopLevelDir() {
-
static bool first = true;
static std::string topLevelVarVal;
if (first) {
- char * const var = getenv("NUISANCE");
+ char* const var = getenv("NUISANCE");
if (!var) {
- THROW("Cannot find top level directory! Set the NUISANCE environmental variable");
+ ERROR(FTL,
+ "Cannot find top level directory! Set the NUISANCE environmental "
+ "variable");
+ throw;
}
topLevelVarVal = std::string(var);
first = false;
}
return topLevelVarVal;
}
-
-std::string GeneralUtils::ReplaceAll(std::string const &inp, std::string const &from,
- std::string const &to) {
+std::string GeneralUtils::ReplaceAll(std::string const& inp,
+ std::string const& from,
+ std::string const& to) {
std::stringstream ss("");
size_t nextOccurence = 0;
size_t prevOccurence = 0;
bool AtEnd = false;
while (!AtEnd) {
nextOccurence = inp.find(from, prevOccurence);
if (nextOccurence == std::string::npos) {
if (prevOccurence == inp.length()) {
break;
}
AtEnd = true;
}
if ((nextOccurence != prevOccurence) || (nextOccurence == 0)) {
ss << inp.substr(prevOccurence, (nextOccurence - prevOccurence));
if (!AtEnd) {
ss << to;
}
}
prevOccurence = nextOccurence + from.size();
}
return ss.str();
}