Public Member Functions | Public Attributes | Private Types | Private Member Functions | Private Attributes | List of all members
CAF Class Reference

#include <CAF.h>

Public Member Functions

 CAF ()
 
 CAF (std::string infile, std::string filename, int correct4origin, double *originTPC)
 
 CAF (const CAF &)=default
 
 ~CAF ()
 
bool BookTFile ()
 
void FillTTree ()
 
void WriteTTree ()
 
void CloseTFile ()
 
void ClearVectors ()
 
void loop ()
 
bool CheckVectorSize ()
 
 CAF (std::string filename, bool isGas=false)
 
 ~CAF ()
 
void fill ()
 
void fillPOT ()
 
void write ()
 
void addRWbranch (int parId, std::string name, std::string wgt_var, std::vector< double > &vars)
 
void Print ()
 
void setToBS ()
 
 CAF (std::string filename, bool isGas=false)
 
 ~CAF ()
 
void fill ()
 
void fillPOT ()
 
void write ()
 
void addRWbranch (int parId, std::string name, std::string wgt_var, std::vector< double > &vars)
 
void Print ()
 
void setToBS ()
 

Public Attributes

int isFD
 
int isFHC
 
int run
 
int subrun
 
int event
 
int ievt
 
int nwgt [100]
 
double cvwgt [100]
 
double wgt [100][100]
 
bool iswgt [100]
 
genie::NtpMCEventRecord * mcrec
 
double pot
 
int meta_run
 
int meta_subrun
 
int version
 
TTree * cafPOT
 
TTree * genie
 
int isCC
 
int neutrinoPDG
 
int neutrinoPDGunosc
 
int mode
 
int LepPDG
 
double Ev
 
double Q2
 
double W
 
double X
 
double Y
 
double NuMomX
 
double NuMomY
 
double NuMomZ
 
double LepMomX
 
double LepMomY
 
double LepMomZ
 
double LepE
 
double LepNuAngle
 
int nP
 
int nN
 
int nipip
 
int nipim
 
int nipi0
 
int nikp
 
int nikm
 
int nik0
 
int niem
 
int niother
 
int nNucleus
 
int nUNKNOWN
 
double eP
 
double eN
 
double ePip
 
double ePim
 
double ePi0
 
double eOther
 
double eRecoP
 
double eRecoN
 
double eRecoPip
 
double eRecoPim
 
double eRecoPi0
 
double eRecoOther
 
double vtx_x
 
double vtx_y
 
double vtx_z
 
double det_x
 
double Ev_reco
 
double Elep_reco
 
double theta_reco
 
int reco_numu
 
int reco_nue
 
int reco_nc
 
int reco_q
 
int muon_contained
 
int muon_tracker
 
int muon_ecal
 
int muon_exit
 
int reco_lepton_pdg
 
float muon_endpoint [3]
 
std::string * muon_endVolName
 
double Ehad_veto
 
double pileup_energy
 
int gastpc_pi_min_mult
 
int gastpc_pi_pl_mult
 
int nFSP
 
int pdg [100]
 
double trkLen [100]
 
double trkLenPerp [100]
 
double ptrue [100]
 
double partEvReco [100]
 
std::vector< std::vector< std::vector< uint64_t > > > * geoEffThrowResults
 

Private Types

typedef std::pair< float, std::string > P
 

Private Member Functions

float calcGluck (double sigmaX, double B, double X0, float nHits, double mom, double length, double &ratio)
 

Private Attributes

TFile * cafFile
 The output TFile pointer */. More...
 
TTree * cafMVA
 The output TTree pointer */. More...
 
std::unordered_map< int, TH2F * > m_pidinterp
 
TFile * _intfile
 
TTree * _inttree
 
Utils * _util
 
std::string _inputfile
 
std::string _outputFile
 The output TFile name */. More...
 
unsigned int _correct4origin
 sets the string for the coordinates origins (World or TPC) More...
 
unsigned int _Run
 
unsigned int _Event
 
unsigned int _SubRun
 
std::vector< int > mode
 
std::vector< int > ccnc
 
std::vector< int > ntype
 
std::vector< int > gint
 
std::vector< int > weight
 
std::vector< int > tgtpdg
 
std::vector< int > gt_t
 
std::vector< int > intert
 
std::vector< int > detected
 
std::vector< int > nGPart
 
std::vector< int > GPartPdg
 
std::vector< int > GPartStatus
 
std::vector< int > GPartFirstMom
 
std::vector< int > GPartLastMom
 
std::vector< int > GPartFirstDaugh
 
std::vector< int > GPartLastDaugh
 
std::vector< float > GPartPx
 
std::vector< float > GPartPy
 
std::vector< float > GPartPz
 
std::vector< float > GPartE
 
std::vector< float > GPartMass
 
std::vector< std::string > GPartName
 
std::vector< double > q2
 
std::vector< double > w
 
std::vector< double > y
 
std::vector< double > x
 
std::vector< double > theta
 
std::vector< double > t
 
std::vector< double > mctime
 
std::vector< double > mcnupx
 
std::vector< double > mcnupy
 
std::vector< double > mcnupz
 
std::vector< double > vertx
 
std::vector< double > verty
 
std::vector< double > vertz
 
std::vector< unsigned int > _nFSP
 
std::vector< int > mctrkid
 
std::vector< int > motherid
 
std::vector< int > pdgmother
 
std::vector< int > truepdg
 
std::vector< int > _MCPStartX
 
std::vector< int > _MCPStartY
 
std::vector< int > _MCPStartZ
 
std::vector< int > _MCPEndX
 
std::vector< int > _MCPEndY
 
std::vector< int > _MCPEndZ
 
std::vector< std::string > _MCProc
 
std::vector< std::string > _MCEndProc
 
std::vector< double > trkLen
 
std::vector< double > trkLenPerp
 
std::vector< double > truep
 
std::vector< double > truepx
 
std::vector< double > truepy
 
std::vector< double > truepz
 
std::vector< double > _angle
 
std::vector< int > recopid
 
std::vector< int > recopidecal
 
std::vector< double > prob_arr
 
std::vector< double > anglereco
 
std::vector< double > _preco
 
std::vector< double > erecon
 
std::vector< double > etime
 
std::vector< unsigned int > isFidStart
 
std::vector< unsigned int > isTPCStart
 
std::vector< unsigned int > isCaloStart
 
std::vector< unsigned int > isInBetweenStart
 
std::vector< unsigned int > isThroughCaloStart
 
std::vector< unsigned int > isFidEnd
 
std::vector< unsigned int > isTPCEnd
 
std::vector< unsigned int > isCaloEnd
 
std::vector< unsigned int > isInBetweenEnd
 
std::vector< unsigned int > isThroughCaloEnd
 
std::vector< unsigned int > isBarrelStart
 
std::vector< unsigned int > isEndcapStart
 
std::vector< unsigned int > isBarrelEnd
 
std::vector< unsigned int > isEndcapEnd
 

Detailed Description

Definition at line 12 of file CAF.h.

Member Typedef Documentation

typedef std::pair<float, std::string> CAF::P
private

Definition at line 14 of file CAF.h.

Constructor & Destructor Documentation

CAF::CAF ( )

Definition at line 22 of file CAF.cpp.

23 : cafFile(nullptr), _intfile(nullptr), _inttree(nullptr), _util(new Utils()), _inputfile(""), _outputFile(""), _correct4origin(0)
24 {
25 
26 }
CAF::_outputFile
std::string _outputFile
The output TFile name */.
Definition: CAF.h:63
CAF::cafFile
TFile * cafFile
The output TFile pointer */.
Definition: CAF.h:53
CAF::_inputfile
std::string _inputfile
Definition: CAF.h:62
CAF::_intfile
TFile * _intfile
Definition: CAF.h:57
CAF::_util
Utils * _util
Definition: CAF.h:60
CAF::_correct4origin
unsigned int _correct4origin
sets the string for the coordinates origins (World or TPC)
Definition: CAF.h:64
CAF::_inttree
TTree * _inttree
Definition: CAF.h:58
Utils
Definition: Utils.h:7
CAF::CAF ( std::string  infile,
std::string  filename,
int  correct4origin,
double *  originTPC 
)

Definition at line 28 of file CAF.cpp.

29 : cafFile(nullptr), _intfile(nullptr), _inttree(nullptr), _util(new Utils()), _inputfile(infile), _outputFile(filename), _correct4origin(correct4origin)
30 {
31  _util->SetOrigin(originTPC);
32 }
CAF::_outputFile
std::string _outputFile
The output TFile name */.
Definition: CAF.h:63
train.filename
string filename
Definition: train.py:213
CAF::cafFile
TFile * cafFile
The output TFile pointer */.
Definition: CAF.h:53
CAF::_inputfile
std::string _inputfile
Definition: CAF.h:62
CAF::_intfile
TFile * _intfile
Definition: CAF.h:57
infile
string infile
Definition: INukeNucleonCorr.cxx:62
Utils::SetOrigin
void SetOrigin(double *origin)
Definition: Utils.cpp:22
CAF::_util
Utils * _util
Definition: CAF.h:60
CAF::_correct4origin
unsigned int _correct4origin
sets the string for the coordinates origins (World or TPC)
Definition: CAF.h:64
CAF::_inttree
TTree * _inttree
Definition: CAF.h:58
Utils
Definition: Utils.h:7
CAF::CAF ( const CAF &  )
default
CAF::~CAF ( )

Definition at line 34 of file CAF.cpp.

35 {
36  delete _util;
37  delete cafMVA;
38  delete cafFile;
39  delete _inttree;
40  delete _intfile;
41 }
CAF::cafFile
TFile * cafFile
The output TFile pointer */.
Definition: CAF.h:53
CAF::_intfile
TFile * _intfile
Definition: CAF.h:57
CAF::_util
Utils * _util
Definition: CAF.h:60
CAF::_inttree
TTree * _inttree
Definition: CAF.h:58
CAF::cafMVA
TTree * cafMVA
The output TTree pointer */.
Definition: CAF.h:54
CAF::CAF ( std::string  filename,
bool  isGas = false 
)
CAF::~CAF ( )
CAF::CAF ( std::string  filename,
bool  isGas = false 
)
CAF::~CAF ( )

Member Function Documentation

void CAF::addRWbranch ( int  parId,
std::string  name,
std::string  wgt_var,
std::vector< double > &  vars 
)
void CAF::addRWbranch ( int  parId,
std::string  name,
std::string  wgt_var,
std::vector< double > &  vars 
)
bool CAF::BookTFile ( )

Definition at line 43 of file CAF.cpp.

44 {
45  _intfile = new TFile( _inputfile.c_str() );
46  if(nullptr == _intfile){
47  std::cout << "Cannot open file " << _inputfile.c_str() << std::endl;
48  return false;
49  }
50 
51  _inttree = (TTree*) _intfile->Get( "anatree/GArAnaTree" );
52  if(nullptr == _inttree){
53  std::cout << "Cannot find tree anatree/GArAnaTree" << std::endl;
54  return false;
55  }
56 
57  cafFile = new TFile(_outputFile.c_str(), "RECREATE");
58  if(nullptr == cafFile){
59  return false;
60  }
61  else
62  {
63  cafMVA = new TTree( "caf", "caf tree" );
64  cafMVA->SetDirectory(0);
65 
66  //Event number
67  cafMVA->Branch("Event", &_Event);
68  //Run number
69  cafMVA->Branch("Run", &_Run);
70  //Sub Run number
71  cafMVA->Branch("SubRun", &_SubRun);
72 
73  //MC Truth info (Generator)
74  cafMVA->Branch("mode", &mode);
75  cafMVA->Branch("q2", &q2);
76  cafMVA->Branch("w", &w);
77  cafMVA->Branch("y", &y);
78  cafMVA->Branch("x", &x);
79  cafMVA->Branch("theta", &theta);
80  cafMVA->Branch("t", &t);
81  cafMVA->Branch("ntype", &ntype);
82  cafMVA->Branch("ccnc", &ccnc);
83  cafMVA->Branch("gint", &gint);
84  cafMVA->Branch("tgtpdg", &tgtpdg);
85  cafMVA->Branch("weight", &weight);
86  cafMVA->Branch("gt_t", &gt_t);
87  cafMVA->Branch("intert", &intert);
88  cafMVA->Branch("mcnupx", &mcnupx);
89  cafMVA->Branch("mcnupy", &mcnupy);
90  cafMVA->Branch("mcnupz", &mcnupz);
91  cafMVA->Branch("vertx", &vertx);
92  cafMVA->Branch("verty", &verty);
93  cafMVA->Branch("vertz", &vertz);
94 
95  //List of particles from GENIE from the interaction and their status flag
96  cafMVA->Branch("nGPart", &nGPart);
97  cafMVA->Branch("GPartName", &GPartName);
98  cafMVA->Branch("GPartPdg", &GPartPdg);
99  cafMVA->Branch("GPartStatus", &GPartStatus);
100  cafMVA->Branch("GPartFirstMom", &GPartFirstMom);
101  cafMVA->Branch("GPartLastMom", &GPartLastMom);
102  cafMVA->Branch("GPartFirstDaugh", &GPartFirstDaugh);
103  cafMVA->Branch("GPartLastDaugh", &GPartLastDaugh);
104  cafMVA->Branch("GPartPx", &GPartPx);
105  cafMVA->Branch("GPartPy", &GPartPy);
106  cafMVA->Branch("GPartPz", &GPartPz);
107  cafMVA->Branch("GPartE", &GPartE);
108  cafMVA->Branch("GPartMass", &GPartMass);
109 
110  //Number of final state particle (primaries)
111  cafMVA->Branch("nFSP", &_nFSP);
112  //MC Particle info
113  cafMVA->Branch("detected", &detected);
114  cafMVA->Branch("pdgmother", &pdgmother);
115  cafMVA->Branch("MCPTime", &mctime);
116  cafMVA->Branch("MCPStartX", &_MCPStartX);
117  cafMVA->Branch("MCPStartY", &_MCPStartY);
118  cafMVA->Branch("MCPStartZ", &_MCPStartZ);
119  cafMVA->Branch("motherid", &motherid);
120  cafMVA->Branch("mctrkid", &mctrkid);
121  cafMVA->Branch("truepx", &truepx);
122  cafMVA->Branch("truepy", &truepy);
123  cafMVA->Branch("truepz", &truepz);
124  cafMVA->Branch("MCPEndX", &_MCPEndX);
125  cafMVA->Branch("MCPEndY", &_MCPEndY);
126  cafMVA->Branch("MCPEndZ", &_MCPEndZ);
127  cafMVA->Branch("MCProc", &_MCProc);
128  cafMVA->Branch("MCEndProc", &_MCEndProc);
129  cafMVA->Branch("angle", &_angle);
130  cafMVA->Branch("truep", &truep);
131  cafMVA->Branch("truepdg", &truepdg);
132  //Reco info
133  cafMVA->Branch("recopid", &recopid);
134  cafMVA->Branch("recopidecal", &recopidecal);
135  cafMVA->Branch("trkLen", &trkLen);
136  cafMVA->Branch("trkLenPerp", &trkLenPerp);
137  cafMVA->Branch("preco", &_preco);
138  cafMVA->Branch("anglereco", &anglereco);
139  cafMVA->Branch("erecon", &erecon);
140  cafMVA->Branch("etime", &etime);
141  cafMVA->Branch("prob_arr", &prob_arr);
142  //Geometry
143  cafMVA->Branch("isFidStart", &isFidStart);
144  cafMVA->Branch("isTPCStart", &isTPCStart);
145  cafMVA->Branch("isCaloStart", &isCaloStart);
146  cafMVA->Branch("isInBetweenStart", &isInBetweenStart);
147  cafMVA->Branch("isThroughCaloStart", &isThroughCaloStart);
148  cafMVA->Branch("isBarrelStart", &isBarrelStart);
149  cafMVA->Branch("isEndcapStart", &isEndcapStart);
150 
151  cafMVA->Branch("isFidEnd", &isFidEnd);
152  cafMVA->Branch("isTPCEnd", &isTPCEnd);
153  cafMVA->Branch("isCaloEnd", &isCaloEnd);
154  cafMVA->Branch("isThroughCaloEnd", &isThroughCaloEnd);
155  cafMVA->Branch("isInBetweenEnd", &isInBetweenEnd);
156  cafMVA->Branch("isBarrelEnd", &isBarrelEnd);
157  cafMVA->Branch("isEndcapEnd", &isEndcapEnd);
158 
159  return true;
160  }
161 }
CAF::_angle
std::vector< double > _angle
Definition: CAF.h:78
CAF::tgtpdg
std::vector< int > tgtpdg
Definition: CAF.h:69
CAF::GPartFirstDaugh
std::vector< int > GPartFirstDaugh
Definition: CAF.h:70
CAF::w
std::vector< double > w
Definition: CAF.h:73
CAF::isInBetweenStart
std::vector< unsigned int > isInBetweenStart
Definition: CAF.h:83
CAF::GPartStatus
std::vector< int > GPartStatus
Definition: CAF.h:70
CAF::_outputFile
std::string _outputFile
The output TFile name */.
Definition: CAF.h:63
CAF::recopid
std::vector< int > recopid
Definition: CAF.h:80
CAF::GPartPz
std::vector< float > GPartPz
Definition: CAF.h:71
CAF::x
std::vector< double > x
Definition: CAF.h:73
CAF::_MCEndProc
std::vector< std::string > _MCEndProc
Definition: CAF.h:77
CAF::isBarrelStart
std::vector< unsigned int > isBarrelStart
Definition: CAF.h:85
CAF::ntype
std::vector< int > ntype
Definition: CAF.h:69
CAF::truep
std::vector< double > truep
Definition: CAF.h:78
CAF::etime
std::vector< double > etime
Definition: CAF.h:81
CAF::_SubRun
unsigned int _SubRun
Definition: CAF.h:67
CAF::GPartPy
std::vector< float > GPartPy
Definition: CAF.h:71
CAF::trkLen
std::vector< double > trkLen
Definition: CAF.h:78
CAF::trkLenPerp
std::vector< double > trkLenPerp
Definition: CAF.h:78
CAF::_MCPStartY
std::vector< int > _MCPStartY
Definition: CAF.h:76
CAF::_nFSP
std::vector< unsigned int > _nFSP
Definition: CAF.h:75
CAF::GPartName
std::vector< std::string > GPartName
Definition: CAF.h:72
CAF::_preco
std::vector< double > _preco
Definition: CAF.h:81
CAF::isEndcapEnd
std::vector< unsigned int > isEndcapEnd
Definition: CAF.h:85
CAF::motherid
std::vector< int > motherid
Definition: CAF.h:76
CAF::isFidStart
std::vector< unsigned int > isFidStart
Definition: CAF.h:83
CAF::GPartMass
std::vector< float > GPartMass
Definition: CAF.h:71
CAF::mctime
std::vector< double > mctime
Definition: CAF.h:73
CAF::isBarrelEnd
std::vector< unsigned int > isBarrelEnd
Definition: CAF.h:85
CAF::gint
std::vector< int > gint
Definition: CAF.h:69
CAF::mcnupy
std::vector< double > mcnupy
Definition: CAF.h:73
CAF::t
std::vector< double > t
Definition: CAF.h:73
CAF::recopidecal
std::vector< int > recopidecal
Definition: CAF.h:80
CAF::isThroughCaloStart
std::vector< unsigned int > isThroughCaloStart
Definition: CAF.h:83
CAF::erecon
std::vector< double > erecon
Definition: CAF.h:81
CAF::GPartFirstMom
std::vector< int > GPartFirstMom
Definition: CAF.h:70
CAF::mcnupz
std::vector< double > mcnupz
Definition: CAF.h:73
CAF::theta
std::vector< double > theta
Definition: CAF.h:73
CAF::cafFile
TFile * cafFile
The output TFile pointer */.
Definition: CAF.h:53
CAF::_MCPEndZ
std::vector< int > _MCPEndZ
Definition: CAF.h:76
CAF::_MCPStartX
std::vector< int > _MCPStartX
Definition: CAF.h:76
CAF::_inputfile
std::string _inputfile
Definition: CAF.h:62
CAF::_intfile
TFile * _intfile
Definition: CAF.h:57
CAF::detected
std::vector< int > detected
Definition: CAF.h:69
CAF::isThroughCaloEnd
std::vector< unsigned int > isThroughCaloEnd
Definition: CAF.h:84
CAF::_Event
unsigned int _Event
Definition: CAF.h:67
CAF::intert
std::vector< int > intert
Definition: CAF.h:69
CAF::anglereco
std::vector< double > anglereco
Definition: CAF.h:81
CAF::_MCPStartZ
std::vector< int > _MCPStartZ
Definition: CAF.h:76
CAF::truepdg
std::vector< int > truepdg
Definition: CAF.h:76
CAF::verty
std::vector< double > verty
Definition: CAF.h:73
CAF::GPartE
std::vector< float > GPartE
Definition: CAF.h:71
CAF::ccnc
std::vector< int > ccnc
Definition: CAF.h:69
CAF::_Run
unsigned int _Run
Definition: CAF.h:67
CAF::q2
std::vector< double > q2
Definition: CAF.h:73
CAF::isInBetweenEnd
std::vector< unsigned int > isInBetweenEnd
Definition: CAF.h:84
CAF::nGPart
std::vector< int > nGPart
Definition: CAF.h:70
CAF::pdgmother
std::vector< int > pdgmother
Definition: CAF.h:76
CAF::_MCPEndY
std::vector< int > _MCPEndY
Definition: CAF.h:76
CAF::isTPCStart
std::vector< unsigned int > isTPCStart
Definition: CAF.h:83
CAF::vertz
std::vector< double > vertz
Definition: CAF.h:73
CAF::truepz
std::vector< double > truepz
Definition: CAF.h:78
CAF::vertx
std::vector< double > vertx
Definition: CAF.h:73
CAF::truepx
std::vector< double > truepx
Definition: CAF.h:78
CAF::weight
std::vector< int > weight
Definition: CAF.h:69
CAF::GPartLastMom
std::vector< int > GPartLastMom
Definition: CAF.h:70
CAF::GPartPx
std::vector< float > GPartPx
Definition: CAF.h:71
CAF::_inttree
TTree * _inttree
Definition: CAF.h:58
CAF::isCaloEnd
std::vector< unsigned int > isCaloEnd
Definition: CAF.h:84
CAF::cafMVA
TTree * cafMVA
The output TTree pointer */.
Definition: CAF.h:54
CAF::isFidEnd
std::vector< unsigned int > isFidEnd
Definition: CAF.h:84
CAF::GPartPdg
std::vector< int > GPartPdg
Definition: CAF.h:70
CAF::gt_t
std::vector< int > gt_t
Definition: CAF.h:69
CAF::mode
std::vector< int > mode
Definition: CAF.h:69
CAF::isTPCEnd
std::vector< unsigned int > isTPCEnd
Definition: CAF.h:84
CAF::y
std::vector< double > y
Definition: CAF.h:73
CAF::isCaloStart
std::vector< unsigned int > isCaloStart
Definition: CAF.h:83
CAF::truepy
std::vector< double > truepy
Definition: CAF.h:78
CAF::mcnupx
std::vector< double > mcnupx
Definition: CAF.h:73
CAF::isEndcapStart
std::vector< unsigned int > isEndcapStart
Definition: CAF.h:85
CAF::mctrkid
std::vector< int > mctrkid
Definition: CAF.h:76
CAF::GPartLastDaugh
std::vector< int > GPartLastDaugh
Definition: CAF.h:70
CAF::prob_arr
std::vector< double > prob_arr
Definition: CAF.h:81
CAF::_MCPEndX
std::vector< int > _MCPEndX
Definition: CAF.h:76
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
CAF::_MCProc
std::vector< std::string > _MCProc
Definition: CAF.h:77
float CAF::calcGluck ( double  sigmaX,
double  B,
double  X0,
float  nHits,
double  mom,
double  length,
double &  ratio 
)
private

Definition at line 337 of file CAF.cpp.

338 {
339  double sig_meas = sqrt(720./(nHits+4))*( (0.01*sigmaX*mom)/(0.0001*0.3*B*length*length) );
340  double sig_mcs = (0.052/B)*sqrt( 1.43/(0.0001*X0*length) );
341  double sigma = sqrt(sig_meas*sig_meas + sig_mcs*sig_mcs);
342  ratio = sig_meas/sig_mcs;
343 
344  return sigma;
345 }
bool CAF::CheckVectorSize ( )

Definition at line 282 of file CAF.cpp.

283 {
284  bool isOK = true;
285 
286  if(_nFSP.at(0) != pdgmother.size()) isOK = false;
287  if(_nFSP.at(0) != truepdg.size()) isOK = false;
288  if(_nFSP.at(0) != mctime.size()) isOK = false;
289  if(_nFSP.at(0) != mctrkid.size()) isOK = false;
290  if(_nFSP.at(0) != motherid.size()) isOK = false;
291  if(_nFSP.at(0) != _MCPStartX.size()) isOK = false;
292  if(_nFSP.at(0) != _MCPStartY.size()) isOK = false;
293  if(_nFSP.at(0) != _MCPStartZ.size()) isOK = false;
294  if(_nFSP.at(0) != _MCPEndX.size()) isOK = false;
295  if(_nFSP.at(0) != _MCPEndY.size()) isOK = false;
296  if(_nFSP.at(0) != _MCPEndZ.size()) isOK = false;
297  if(_nFSP.at(0) != _MCProc.size()) isOK = false;
298  if(_nFSP.at(0) != _MCEndProc.size()) isOK = false;
299  if(_nFSP.at(0) != trkLen.size()) isOK = false;
300  if(_nFSP.at(0) != trkLenPerp.size()) isOK = false;
301  if(_nFSP.at(0) != truep.size()) isOK = false;
302  if(_nFSP.at(0) != truepx.size()) isOK = false;
303  if(_nFSP.at(0) != truepy.size()) isOK = false;
304  if(_nFSP.at(0) != truepz.size()) isOK = false;
305  if(_nFSP.at(0) != _angle.size()) isOK = false;
306 
307  //Reco values
308  if(_nFSP.at(0) != recopid.size()) isOK = false;
309  if(_nFSP.at(0) != detected.size()) isOK = false;
310  if(_nFSP.at(0) != recopidecal.size()) isOK = false;
311  // if(_nFSP.at(0) != prob_arr.size()) isOK = false;
312  if(_nFSP.at(0) != anglereco.size()) isOK = false;
313  if(_nFSP.at(0) != _preco.size()) isOK = false;
314  if(_nFSP.at(0) != erecon.size()) isOK = false;
315  if(_nFSP.at(0) != etime.size()) isOK = false;
316 
317  //Geometry
318  if(_nFSP.at(0) != isFidStart.size()) isOK = false;
319  if(_nFSP.at(0) != isTPCStart.size()) isOK = false;
320  if(_nFSP.at(0) != isCaloStart.size()) isOK = false;
321  if(_nFSP.at(0) != isThroughCaloStart.size()) isOK = false;
322  if(_nFSP.at(0) != isInBetweenStart.size()) isOK = false;
323  if(_nFSP.at(0) != isBarrelStart.size()) isOK = false;
324  if(_nFSP.at(0) != isEndcapStart.size()) isOK = false;
325 
326  if(_nFSP.at(0) != isFidEnd.size()) isOK = false;
327  if(_nFSP.at(0) != isTPCEnd.size()) isOK = false;
328  if(_nFSP.at(0) != isCaloEnd.size()) isOK = false;
329  if(_nFSP.at(0) != isThroughCaloEnd.size()) isOK = false;
330  if(_nFSP.at(0) != isInBetweenEnd.size()) isOK = false;
331  if(_nFSP.at(0) != isBarrelEnd.size()) isOK = false;
332  if(_nFSP.at(0) != isEndcapEnd.size()) isOK = false;
333 
334  return isOK;
335 }
CAF::_angle
std::vector< double > _angle
Definition: CAF.h:78
CAF::isInBetweenStart
std::vector< unsigned int > isInBetweenStart
Definition: CAF.h:83
CAF::recopid
std::vector< int > recopid
Definition: CAF.h:80
CAF::_MCEndProc
std::vector< std::string > _MCEndProc
Definition: CAF.h:77
CAF::isBarrelStart
std::vector< unsigned int > isBarrelStart
Definition: CAF.h:85
CAF::truep
std::vector< double > truep
Definition: CAF.h:78
CAF::etime
std::vector< double > etime
Definition: CAF.h:81
CAF::trkLen
std::vector< double > trkLen
Definition: CAF.h:78
CAF::trkLenPerp
std::vector< double > trkLenPerp
Definition: CAF.h:78
CAF::_MCPStartY
std::vector< int > _MCPStartY
Definition: CAF.h:76
CAF::_nFSP
std::vector< unsigned int > _nFSP
Definition: CAF.h:75
CAF::_preco
std::vector< double > _preco
Definition: CAF.h:81
CAF::isEndcapEnd
std::vector< unsigned int > isEndcapEnd
Definition: CAF.h:85
CAF::motherid
std::vector< int > motherid
Definition: CAF.h:76
CAF::isFidStart
std::vector< unsigned int > isFidStart
Definition: CAF.h:83
CAF::mctime
std::vector< double > mctime
Definition: CAF.h:73
CAF::isBarrelEnd
std::vector< unsigned int > isBarrelEnd
Definition: CAF.h:85
CAF::recopidecal
std::vector< int > recopidecal
Definition: CAF.h:80
CAF::isThroughCaloStart
std::vector< unsigned int > isThroughCaloStart
Definition: CAF.h:83
CAF::erecon
std::vector< double > erecon
Definition: CAF.h:81
CAF::_MCPEndZ
std::vector< int > _MCPEndZ
Definition: CAF.h:76
CAF::_MCPStartX
std::vector< int > _MCPStartX
Definition: CAF.h:76
CAF::detected
std::vector< int > detected
Definition: CAF.h:69
CAF::isThroughCaloEnd
std::vector< unsigned int > isThroughCaloEnd
Definition: CAF.h:84
CAF::anglereco
std::vector< double > anglereco
Definition: CAF.h:81
CAF::_MCPStartZ
std::vector< int > _MCPStartZ
Definition: CAF.h:76
CAF::truepdg
std::vector< int > truepdg
Definition: CAF.h:76
CAF::isInBetweenEnd
std::vector< unsigned int > isInBetweenEnd
Definition: CAF.h:84
CAF::pdgmother
std::vector< int > pdgmother
Definition: CAF.h:76
CAF::_MCPEndY
std::vector< int > _MCPEndY
Definition: CAF.h:76
CAF::isTPCStart
std::vector< unsigned int > isTPCStart
Definition: CAF.h:83
CAF::truepz
std::vector< double > truepz
Definition: CAF.h:78
CAF::truepx
std::vector< double > truepx
Definition: CAF.h:78
CAF::isCaloEnd
std::vector< unsigned int > isCaloEnd
Definition: CAF.h:84
CAF::isFidEnd
std::vector< unsigned int > isFidEnd
Definition: CAF.h:84
CAF::isTPCEnd
std::vector< unsigned int > isTPCEnd
Definition: CAF.h:84
CAF::isCaloStart
std::vector< unsigned int > isCaloStart
Definition: CAF.h:83
CAF::truepy
std::vector< double > truepy
Definition: CAF.h:78
CAF::isEndcapStart
std::vector< unsigned int > isEndcapStart
Definition: CAF.h:85
CAF::mctrkid
std::vector< int > mctrkid
Definition: CAF.h:76
CAF::_MCPEndX
std::vector< int > _MCPEndX
Definition: CAF.h:76
CAF::_MCProc
std::vector< std::string > _MCProc
Definition: CAF.h:77
void CAF::ClearVectors ( )

Definition at line 193 of file CAF.cpp.

194 {
195  //Generator values
196  mode.clear();
197  ccnc.clear();
198  ntype.clear();
199  gint.clear();
200  weight.clear();
201  tgtpdg.clear();
202  gt_t.clear();
203  intert.clear();
204  q2.clear();
205  w.clear();
206  y.clear();
207  x.clear();
208  theta.clear();
209  t.clear();
210  mcnupx.clear();
211  mcnupy.clear();
212  mcnupz.clear();
213  vertx.clear();
214  verty.clear();
215  vertz.clear();
216 
217  nGPart.clear();
218  GPartPdg.clear();
219  GPartStatus.clear();
220  GPartName.clear();
221  GPartFirstMom.clear();
222  GPartLastMom.clear();
223  GPartFirstDaugh.clear();
224  GPartLastDaugh.clear();
225  GPartPx.clear();
226  GPartPy.clear();
227  GPartPz.clear();
228  GPartE.clear();
229  GPartMass.clear();
230 
231  //MC Particle values
232  _nFSP.clear();
233  pdgmother.clear();
234  truepdg.clear();
235  mctime.clear();
236  mctrkid.clear();
237  motherid.clear();
238  _MCPStartX.clear();
239  _MCPStartY.clear();
240  _MCPStartZ.clear();
241  _MCPEndX.clear();
242  _MCPEndY.clear();
243  _MCPEndZ.clear();
244  _MCProc.clear();
245  _MCEndProc.clear();
246  trkLen.clear();
247  trkLenPerp.clear();
248  truep.clear();
249  truepx.clear();
250  truepy.clear();
251  truepz.clear();
252  _angle.clear();
253 
254  //Reco values
255  detected.clear();
256  recopid.clear();
257  recopidecal.clear();
258  prob_arr.clear();
259  anglereco.clear();
260  _preco.clear();
261  erecon.clear();
262  etime.clear();
263 
264  //Geometry
265  isFidStart.clear();
266  isTPCStart.clear();
267  isCaloStart.clear();
268  isThroughCaloStart.clear();
269  isInBetweenStart.clear();
270  isBarrelStart.clear();
271  isEndcapStart.clear();
272 
273  isFidEnd.clear();
274  isTPCEnd.clear();
275  isCaloEnd.clear();
276  isThroughCaloEnd.clear();
277  isInBetweenEnd.clear();
278  isBarrelEnd.clear();
279  isEndcapEnd.clear();
280 }
CAF::_angle
std::vector< double > _angle
Definition: CAF.h:78
CAF::tgtpdg
std::vector< int > tgtpdg
Definition: CAF.h:69
CAF::GPartFirstDaugh
std::vector< int > GPartFirstDaugh
Definition: CAF.h:70
CAF::w
std::vector< double > w
Definition: CAF.h:73
CAF::isInBetweenStart
std::vector< unsigned int > isInBetweenStart
Definition: CAF.h:83
CAF::GPartStatus
std::vector< int > GPartStatus
Definition: CAF.h:70
CAF::recopid
std::vector< int > recopid
Definition: CAF.h:80
CAF::GPartPz
std::vector< float > GPartPz
Definition: CAF.h:71
CAF::x
std::vector< double > x
Definition: CAF.h:73
CAF::_MCEndProc
std::vector< std::string > _MCEndProc
Definition: CAF.h:77
CAF::isBarrelStart
std::vector< unsigned int > isBarrelStart
Definition: CAF.h:85
CAF::ntype
std::vector< int > ntype
Definition: CAF.h:69
CAF::truep
std::vector< double > truep
Definition: CAF.h:78
CAF::etime
std::vector< double > etime
Definition: CAF.h:81
CAF::GPartPy
std::vector< float > GPartPy
Definition: CAF.h:71
CAF::trkLen
std::vector< double > trkLen
Definition: CAF.h:78
CAF::trkLenPerp
std::vector< double > trkLenPerp
Definition: CAF.h:78
CAF::_MCPStartY
std::vector< int > _MCPStartY
Definition: CAF.h:76
CAF::_nFSP
std::vector< unsigned int > _nFSP
Definition: CAF.h:75
CAF::GPartName
std::vector< std::string > GPartName
Definition: CAF.h:72
CAF::_preco
std::vector< double > _preco
Definition: CAF.h:81
CAF::isEndcapEnd
std::vector< unsigned int > isEndcapEnd
Definition: CAF.h:85
CAF::motherid
std::vector< int > motherid
Definition: CAF.h:76
CAF::isFidStart
std::vector< unsigned int > isFidStart
Definition: CAF.h:83
CAF::GPartMass
std::vector< float > GPartMass
Definition: CAF.h:71
CAF::mctime
std::vector< double > mctime
Definition: CAF.h:73
CAF::isBarrelEnd
std::vector< unsigned int > isBarrelEnd
Definition: CAF.h:85
CAF::gint
std::vector< int > gint
Definition: CAF.h:69
CAF::mcnupy
std::vector< double > mcnupy
Definition: CAF.h:73
CAF::t
std::vector< double > t
Definition: CAF.h:73
CAF::recopidecal
std::vector< int > recopidecal
Definition: CAF.h:80
CAF::isThroughCaloStart
std::vector< unsigned int > isThroughCaloStart
Definition: CAF.h:83
CAF::erecon
std::vector< double > erecon
Definition: CAF.h:81
CAF::GPartFirstMom
std::vector< int > GPartFirstMom
Definition: CAF.h:70
CAF::mcnupz
std::vector< double > mcnupz
Definition: CAF.h:73
CAF::theta
std::vector< double > theta
Definition: CAF.h:73
CAF::_MCPEndZ
std::vector< int > _MCPEndZ
Definition: CAF.h:76
CAF::_MCPStartX
std::vector< int > _MCPStartX
Definition: CAF.h:76
CAF::detected
std::vector< int > detected
Definition: CAF.h:69
CAF::isThroughCaloEnd
std::vector< unsigned int > isThroughCaloEnd
Definition: CAF.h:84
CAF::intert
std::vector< int > intert
Definition: CAF.h:69
CAF::anglereco
std::vector< double > anglereco
Definition: CAF.h:81
CAF::_MCPStartZ
std::vector< int > _MCPStartZ
Definition: CAF.h:76
CAF::truepdg
std::vector< int > truepdg
Definition: CAF.h:76
CAF::verty
std::vector< double > verty
Definition: CAF.h:73
CAF::GPartE
std::vector< float > GPartE
Definition: CAF.h:71
CAF::ccnc
std::vector< int > ccnc
Definition: CAF.h:69
CAF::q2
std::vector< double > q2
Definition: CAF.h:73
CAF::isInBetweenEnd
std::vector< unsigned int > isInBetweenEnd
Definition: CAF.h:84
CAF::nGPart
std::vector< int > nGPart
Definition: CAF.h:70
CAF::pdgmother
std::vector< int > pdgmother
Definition: CAF.h:76
CAF::_MCPEndY
std::vector< int > _MCPEndY
Definition: CAF.h:76
CAF::isTPCStart
std::vector< unsigned int > isTPCStart
Definition: CAF.h:83
CAF::vertz
std::vector< double > vertz
Definition: CAF.h:73
CAF::truepz
std::vector< double > truepz
Definition: CAF.h:78
CAF::vertx
std::vector< double > vertx
Definition: CAF.h:73
CAF::truepx
std::vector< double > truepx
Definition: CAF.h:78
CAF::weight
std::vector< int > weight
Definition: CAF.h:69
CAF::GPartLastMom
std::vector< int > GPartLastMom
Definition: CAF.h:70
CAF::GPartPx
std::vector< float > GPartPx
Definition: CAF.h:71
CAF::isCaloEnd
std::vector< unsigned int > isCaloEnd
Definition: CAF.h:84
CAF::isFidEnd
std::vector< unsigned int > isFidEnd
Definition: CAF.h:84
CAF::GPartPdg
std::vector< int > GPartPdg
Definition: CAF.h:70
CAF::gt_t
std::vector< int > gt_t
Definition: CAF.h:69
CAF::mode
std::vector< int > mode
Definition: CAF.h:69
CAF::isTPCEnd
std::vector< unsigned int > isTPCEnd
Definition: CAF.h:84
CAF::y
std::vector< double > y
Definition: CAF.h:73
CAF::isCaloStart
std::vector< unsigned int > isCaloStart
Definition: CAF.h:83
CAF::truepy
std::vector< double > truepy
Definition: CAF.h:78
CAF::mcnupx
std::vector< double > mcnupx
Definition: CAF.h:73
CAF::isEndcapStart
std::vector< unsigned int > isEndcapStart
Definition: CAF.h:85
CAF::mctrkid
std::vector< int > mctrkid
Definition: CAF.h:76
CAF::GPartLastDaugh
std::vector< int > GPartLastDaugh
Definition: CAF.h:70
CAF::prob_arr
std::vector< double > prob_arr
Definition: CAF.h:81
CAF::_MCPEndX
std::vector< int > _MCPEndX
Definition: CAF.h:76
CAF::_MCProc
std::vector< std::string > _MCProc
Definition: CAF.h:77
void CAF::CloseTFile ( )

Definition at line 163 of file CAF.cpp.

164 {
165  if(cafFile != nullptr) {
166  cafFile->Close();
167  return;
168  }
169  else{ return; }
170 }
CAF::cafFile
TFile * cafFile
The output TFile pointer */.
Definition: CAF.h:53
void CAF::fill ( )
void CAF::fill ( )
void CAF::fillPOT ( )
void CAF::fillPOT ( )
void CAF::FillTTree ( )

Definition at line 183 of file CAF.cpp.

184 {
185  if(cafMVA != nullptr) {
186  cafFile->cd();
187  cafMVA->Fill();
188  }
189 
190  return;
191 }
CAF::cafFile
TFile * cafFile
The output TFile pointer */.
Definition: CAF.h:53
CAF::cafMVA
TTree * cafMVA
The output TTree pointer */.
Definition: CAF.h:54
void CAF::loop ( )

Definition at line 348 of file CAF.cpp.

349 {
350  //double gastpc_len = 5.; // track length cut in cm
351  const float gastpc_len = 2.; // new track length cut in cm based on Thomas' study of low energy protons
352  // dont care about electrons -- check momentum and see if hit ECAL
353  const float gastpc_B = 0.5; // B field strength in Tesla
354  const float gastpc_padPitch = 1.0; // 1 mm. Actual pad pitch varies, which is going to be impossible to implement
355  const float gastpc_X0 = 1300.; // cm = 13m radiation length
356  //Resolution for short tracks //TODO check this numbers!
357  const float sigmaP_short = 0.1; //in GeV
358  // point resolution
359  const float sigma_x = 0.1;
360 
361  std::vector<float> v;
362  for (float pit = 0.040; pit < 20.0; pit += 0.001)
363  {
364  v.push_back(pit);
365  }
366 
367  //as function of KE
368  //0 -> 50 MeV ~ 20%
369  //> 50 MeV ~ 40%
370  const float NeutronECAL_detEff[2] = {0.2, 0.4};
371  const float sigmaNeutronECAL_first = 0.11;
372  //TODO fraction of rescatters
373  // float sigmaNeutronECAL_rescatter = 0.26;
374  //ECAL energy resolution sigmaE/E
375  const float ECAL_stock = 0.06; //in %
376  const float ECAL_const = 0.02;
377  TF1 *fRes = new TF1("fRes", "TMath::Sqrt ( [0]*[0]/x + [1]*[1] )");
378  fRes->FixParameter(0, ECAL_stock);
379  fRes->FixParameter(1, ECAL_const);
380  //ECAL sampling fraction
381  // double sampling_frac = 4.32;
382  //ECAL nlayers
383  const int nLayers = 60;
384  //ECAL MIP resolution (based on AHCAL numbers)
385  const double ECAL_MIP_Res = 0.23;
386  //MIP2GeV conversion factor
387  const double MIP2GeV_factor = 0.814 / 1000;
388  //float ECAL_pi0_resolution = 0.13; //sigmaE/E in between at rest (17%) and high energy (~few %)
389  const float ECAL_time_resolution = 1.; // 1 ns time resolution
390  TParticlePDG *neutron = TDatabasePDG::Instance()->GetParticle(2112);
391  const float neutron_mass = neutron->Mass(); //in GeV
392  //TParticlePDG *pi0 = TDatabasePDG::Instance()->GetParticle(111);
393  //float pi0_mass = pi0->Mass(); //in GeV
394 
395  TString filename="${DUNE_PARDATA_DIR}/MPD/dedxPID/dedxpidmatrices8kevcm.root";
396  TFile infile(filename,"READ");
397 
398  m_pidinterp.clear();
399  char str[11];
400  for (int q = 0; q < 501; ++q)
401  {
402  sprintf(str, "%d", q);
403  std::string s = "pidmatrix";
404  s.append(str);
405  // read the 500 histograms one by one; each histogram is a
406  // 6 by 6 matrix of probabilities for a given momentum value
407  m_pidinterp.insert( std::make_pair(q, (TH2F*) infile.Get(s.c_str())->Clone("pidinterp")) );
408  }
409 
410  // infile.Close();
411 
412  //------------------------------------------------------------------------
413 
414  int Event = 0;
415  int SubRun = 0;
416  int Run = 0;
417 
418  std::vector<float> *MC_Q2 = 0;
419  TBranch *b_MC_Q2 = 0;
420  std::vector<float> *MC_W = 0;
421  TBranch *b_MC_W = 0;
422  std::vector<float> *MC_Y = 0;
423  TBranch *b_MC_Y = 0;
424  std::vector<float> *MC_X = 0;
425  TBranch *b_MC_X = 0;
426  std::vector<float> *MC_Theta = 0;
427  TBranch *b_MC_Theta = 0;
428  std::vector<float> *MCVertX = 0;
429  TBranch *b_MCVertX = 0;
430  std::vector<float> *MCVertY = 0;
431  TBranch *b_MCVertY = 0;
432  std::vector<float> *MCVertZ = 0;
433  TBranch *b_MCVertZ = 0;
434  std::vector<float> *MCNuPx = 0;
435  TBranch *b_MCNuPx = 0;
436  std::vector<float> *MCNuPy = 0;
437  TBranch *b_MCNuPy = 0;
438  std::vector<float> *MCNuPz = 0;
439  TBranch *b_MCNuPz = 0;
440  std::vector<int> *NType = 0;
441  TBranch *b_NType = 0;
442  std::vector<int> *CCNC = 0;
443  TBranch *b_CCNC = 0;
444  std::vector<int> *Mode = 0;
445  TBranch *b_Mode = 0;
446  std::vector<int> *Gint=0;
447  TBranch *b_Gint = 0;
448  std::vector<int> *TgtPDG = 0;
449  TBranch *b_TgtPDG = 0;
450  std::vector<int> *GT_T = 0;
451  TBranch *b_GT_T = 0;
452  std::vector<int> *InterT=0;
453  TBranch *b_InterT = 0;
454  std::vector<float> *Weight = 0;
455  TBranch *b_Weight = 0;
456 
457  std::vector<int> *_nGPart = 0;
458  TBranch *b_nGPart = 0;
459  std::vector<int> *_GPartPdg = 0;
460  TBranch *b_GPartPdg = 0;
461  std::vector<int> *_GPartStatus = 0;
462  TBranch *b_GPartStatus = 0;
463  std::vector<std::string> *_GPartName = 0;
464  TBranch *b_GPartName = 0;
465  std::vector<int> *_GPartFirstMom = 0;
466  TBranch *b_GPartFirstMom = 0;
467  std::vector<int> *_GPartLastMom = 0;
468  TBranch *b_GPartLastMom = 0;
469  std::vector<int> *_GPartFirstDaugh = 0;
470  TBranch *b_GPartFirstDaugh = 0;
471  std::vector<int> *_GPartLastDaugh = 0;
472  TBranch *b_GPartLastDaugh = 0;
473  std::vector<float> *_GPartPx = 0;
474  TBranch *b_GPartPx = 0;
475  std::vector<float> *_GPartPy = 0;
476  TBranch *b_GPartPy = 0;
477  std::vector<float> *_GPartPz = 0;
478  TBranch *b_GPartPz = 0;
479  std::vector<float> *_GPartE = 0;
480  TBranch *b_GPartE = 0;
481  std::vector<float> *_GPartMass = 0;
482  TBranch *b_GPartMass = 0;
483 
484  std::vector<int> *PDG = 0;
485  TBranch *b_PDG = 0;
486  std::vector<int> *MCPTrkID = 0;
487  TBranch *b_MCPTrkID = 0;
488  std::vector<int> *MCMotherTrkID = 0;
489  TBranch *b_MCMotherTrkID = 0;
490  std::vector<int> *PDGMother = 0;
491  TBranch *b_PDGMother = 0;
492  std::vector<float> *MCPTime = 0;
493  TBranch *b_MCPTime = 0;
494  std::vector<float> *MCPStartX = 0;
495  TBranch *b_MCPStartX = 0;
496  std::vector<float> *MCPStartY = 0;
497  TBranch *b_MCPStartY = 0;
498  std::vector<float> *MCPStartZ = 0;
499  TBranch *b_MCPStartZ = 0;
500  std::vector<float> *MCPStartPX = 0;
501  TBranch *b_MCPStartPX = 0;
502  std::vector<float> *MCPStartPY = 0;
503  TBranch *b_MCPStartPY = 0;
504  std::vector<float> *MCPStartPZ = 0;
505  TBranch *b_MCPStartPZ = 0;
506  std::vector<std::string> *MCPProc = 0;
507  TBranch *b_MCPProc = 0;
508  std::vector<std::string> *MCPEndProc = 0;
509  TBranch *b_MCPEndProc = 0;
510  std::vector<float> *MCPEndX = 0;
511  TBranch *b_MCPEndX = 0;
512  std::vector<float> *MCPEndY = 0;
513  TBranch *b_MCPEndY = 0;
514  std::vector<float> *MCPEndZ = 0;
515  TBranch *b_MCPEndZ = 0;
516  std::vector<float> *TrajMCPX = 0;
517  TBranch *b_TrajMCPX = 0;
518  std::vector<float> *TrajMCPY = 0;
519  TBranch *b_TrajMCPY = 0;
520  std::vector<float> *TrajMCPZ = 0;
521  TBranch *b_TrajMCPZ = 0;
522  std::vector<int> *TrajMCPTrajIndex = 0;
523  TBranch *b_TrajMCPTrajIndex = 0;
524 
525  _inttree->SetBranchAddress("Event", &Event);
526  _inttree->SetBranchAddress("SubRun", &SubRun);
527  _inttree->SetBranchAddress("Run", &Run);
528 
529  //Generator info
530  _inttree->SetBranchAddress("NType", &NType, &b_NType);
531  _inttree->SetBranchAddress("CCNC", &CCNC, &b_CCNC);
532  _inttree->SetBranchAddress("MC_Q2", &MC_Q2, &b_MC_Q2);
533  _inttree->SetBranchAddress("MC_W", &MC_W, &b_MC_W);
534  _inttree->SetBranchAddress("MC_Y", &MC_Y, &b_MC_Y);
535  _inttree->SetBranchAddress("MC_X", &MC_X, &b_MC_X);
536  _inttree->SetBranchAddress("MC_Theta", &MC_Theta, &b_MC_Theta);
537  _inttree->SetBranchAddress("Mode", &Mode, &b_Mode);
538  _inttree->SetBranchAddress("Gint", &Gint, &b_Gint);
539  _inttree->SetBranchAddress("TgtPDG", &TgtPDG, &b_TgtPDG);
540  _inttree->SetBranchAddress("GT_T", &GT_T, &b_GT_T);
541  _inttree->SetBranchAddress("MCVertX", &MCVertX, &b_MCVertX);
542  _inttree->SetBranchAddress("MCVertY", &MCVertY, &b_MCVertY);
543  _inttree->SetBranchAddress("MCVertZ", &MCVertZ, &b_MCVertZ);
544  _inttree->SetBranchAddress("MCNuPx", &MCNuPx, &b_MCNuPx);
545  _inttree->SetBranchAddress("MCNuPy", &MCNuPy, &b_MCNuPy);
546  _inttree->SetBranchAddress("MCNuPz", &MCNuPz, &b_MCNuPz);
547  _inttree->SetBranchAddress("InterT", &InterT, &b_InterT);
548  _inttree->SetBranchAddress("Weight", &Weight, &b_Weight);
549 
550  _inttree->SetBranchAddress("nGPart", &_nGPart, &b_nGPart);
551  _inttree->SetBranchAddress("GPartPdg", &_GPartPdg, &b_GPartPdg);
552  _inttree->SetBranchAddress("GPartStatus", &_GPartStatus, &b_GPartStatus);
553  _inttree->SetBranchAddress("GPartName", &_GPartName, &b_GPartName);
554  _inttree->SetBranchAddress("GPartFirstMom", &_GPartFirstMom, &b_GPartFirstMom);
555  _inttree->SetBranchAddress("GPartLastMom", &_GPartLastMom, &b_GPartLastMom);
556  _inttree->SetBranchAddress("GPartFirstDaugh", &_GPartFirstDaugh, &b_GPartFirstDaugh);
557  _inttree->SetBranchAddress("GPartLastDaugh", &_GPartLastDaugh, &b_GPartLastDaugh);
558  _inttree->SetBranchAddress("GPartPx", &_GPartPx, &b_GPartPx);
559  _inttree->SetBranchAddress("GPartPy", &_GPartPy, &b_GPartPy);
560  _inttree->SetBranchAddress("GPartPz", &_GPartPz, &b_GPartPz);
561  _inttree->SetBranchAddress("GPartE", &_GPartE, &b_GPartE);
562  _inttree->SetBranchAddress("GPartMass", &_GPartMass, &b_GPartMass);
563 
564  //MC info
565  _inttree->SetBranchAddress("PDG", &PDG, &b_PDG);
566  _inttree->SetBranchAddress("MCTrkID", &MCPTrkID, &b_MCPTrkID);
567  _inttree->SetBranchAddress("MotherTrkID", &MCMotherTrkID, &b_MCMotherTrkID);
568  _inttree->SetBranchAddress("MCPTime", &MCPTime, &b_MCPTime);
569  _inttree->SetBranchAddress("MCPStartX", &MCPStartX, &b_MCPStartX);
570  _inttree->SetBranchAddress("MCPStartY", &MCPStartY, &b_MCPStartY);
571  _inttree->SetBranchAddress("MCPStartZ", &MCPStartZ, &b_MCPStartZ);
572  _inttree->SetBranchAddress("MCPEndX", &MCPEndX, &b_MCPEndX);
573  _inttree->SetBranchAddress("MCPEndY", &MCPEndY, &b_MCPEndY);
574  _inttree->SetBranchAddress("MCPEndZ", &MCPEndZ, &b_MCPEndZ);
575  _inttree->SetBranchAddress("PDGMother", &PDGMother, &b_PDGMother);
576  _inttree->SetBranchAddress("MCPStartPX", &MCPStartPX, &b_MCPStartPX);
577  _inttree->SetBranchAddress("MCPStartPY", &MCPStartPY, &b_MCPStartPY);
578  _inttree->SetBranchAddress("MCPStartPZ", &MCPStartPZ, &b_MCPStartPZ);
579  _inttree->SetBranchAddress("MCPProc", &MCPProc, &b_MCPProc);
580  _inttree->SetBranchAddress("MCPEndProc", &MCPEndProc, &b_MCPEndProc);
581  _inttree->SetBranchAddress("TrajMCPX", &TrajMCPX, &b_TrajMCPX);
582  _inttree->SetBranchAddress("TrajMCPY", &TrajMCPY, &b_TrajMCPY);
583  _inttree->SetBranchAddress("TrajMCPZ", &TrajMCPZ, &b_TrajMCPZ);
584  // _inttree->SetBranchAddress("TrajMCPTrajIndex", &TrajMCPTrajIndex);
585  _inttree->SetBranchAddress("TrajMCPTrackID", &TrajMCPTrajIndex, &b_TrajMCPTrajIndex);
586 
587  //gamma, neutron, pi0, k0L, k0S, k0, delta0
588  std::vector<int> neutrinos = {12, 14, 16};
589  std::vector<int> pdg_neutral = {22, 2112, 111, 130, 310, 311, 2114};
590  //pion, muon, proton, kaon, deuteron, electron
591  std::vector<int> pdg_charged = {211, 13, 2212, 321, 1000010020, 11};
592  //-------------------------------------------------------------------
593 
594  // Main event loop
595  for( int entry = 0; entry < _inttree->GetEntries(); entry++ )
596  {
597  _inttree->GetEntry(entry);
598  this->ClearVectors();
599 
600  //Filling MCTruth values
601  // std::cout << "Event " << Event << " Run " << Run << std::endl;
602  // if(Event < 259) continue;
603 
604  _Event = Event;
605  _Run = Run;
606  _SubRun = SubRun;
607 
608  for(size_t i = 0; i < NType->size(); i++)
609  {
610  ccnc.push_back(CCNC->at(i));
611  ntype.push_back(NType->at(i));
612  q2.push_back(MC_Q2->at(i));
613  w.push_back(MC_W->at(i));
614  y.push_back(MC_Y->at(i));
615  x.push_back(MC_X->at(i));
616  theta.push_back(MC_Theta->at(i));
617  mode.push_back(Mode->at(i));
618  intert.push_back(InterT->at(i));
619  if(_correct4origin){
620  vertx.push_back(MCVertX->at(i) - _util->GetOrigin()[0]);
621  verty.push_back(MCVertY->at(i) - _util->GetOrigin()[1]);
622  vertz.push_back(MCVertZ->at(i) - _util->GetOrigin()[2]);
623  } else {
624  vertx.push_back(MCVertX->at(i));
625  verty.push_back(MCVertY->at(i));
626  vertz.push_back(MCVertZ->at(i));
627  }
628  mcnupx.push_back(MCNuPx->at(i));
629  mcnupy.push_back(MCNuPy->at(i));
630  mcnupz.push_back(MCNuPz->at(i));
631  }
632 
633  for(size_t i = 0; i < Gint->size(); i++)
634  {
635  gint.push_back(Gint->at(i));
636  tgtpdg.push_back(TgtPDG->at(i));//target pdg
637  gt_t.push_back(GT_T->at(i));
638  weight.push_back(Weight->at(i));
639  }
640 
641  for(size_t i = 0; i < _nGPart->size(); i++)
642  {
643  nGPart.push_back(_nGPart->at(i));
644  for(int j = 0; j < _nGPart->at(i); j++) {
645  GPartPdg.push_back(_GPartPdg->at(j));
646  GPartStatus.push_back(_GPartStatus->at(j));
647  GPartName.push_back(_GPartName->at(j));
648  GPartFirstMom.push_back(_GPartFirstMom->at(j));
649  GPartLastMom.push_back(_GPartLastMom->at(j));
650  GPartFirstDaugh.push_back(_GPartFirstDaugh->at(j));
651  GPartLastDaugh.push_back(_GPartLastDaugh->at(j));
652  GPartPx.push_back(_GPartPx->at(j));
653  GPartPy.push_back(_GPartPy->at(j));
654  GPartPz.push_back(_GPartPz->at(j));
655  GPartE.push_back(_GPartE->at(j));
656  GPartMass.push_back(_GPartMass->at(j));
657  }
658  }
659 
660  //--------------------------------------------------------------------------
661  // Start of Parameterized Reconstruction
662  //--------------------------------------------------------------------------
663  unsigned int nFSP = 0;
664 
665  //TODO we should skip particles that we don't see or cannot reconstruct! change filling caf with i to index counting the particle
666  //have to be careful with indexes and continue
667  //---------------------------------------------------------------
668  // all Gluckstern calculations happen in the following loop
669  for(size_t i = 0; i < MCPStartPX->size(); i++ )
670  {
671  //Get the creating process
672  std::string mcp_process = MCPProc->at(i);
673  //Get ending process
674  std::string mcp_endprocess = MCPEndProc->at(i);
675  int mctrackid = MCPTrkID->at(i);
676  int pdg = PDG->at(i);
677  const TVector3 spoint(MCPStartX->at(i) - _util->GetOrigin()[0], MCPStartY->at(i) - _util->GetOrigin()[1], MCPStartZ->at(i) - _util->GetOrigin()[2]);
678  const TVector3 epoint(MCPEndX->at(i) - _util->GetOrigin()[0], MCPEndY->at(i) - _util->GetOrigin()[1], MCPEndZ->at(i) - _util->GetOrigin()[2]);
679 
680  TVector3 mcp(MCPStartPX->at(i), MCPStartPY->at(i), MCPStartPZ->at(i));
681  float ptrue = (mcp).Mag();
682  float pypz = std::sqrt( MCPStartPY->at(i)*MCPStartPY->at(i) + MCPStartPZ->at(i)*MCPStartPZ->at(i));
683 
684  //need to ignore neutrals for this - put the value to 0
685  auto result = std::find(pdg_neutral.begin(), pdg_neutral.end(), abs(pdg));
686  bool isNeutral = (result != pdg_neutral.end()) ? true : false;
687 
688  //start track length
689  //***************************************************************************************************************/
690 
691  if( isNeutral )
692  {
693  trkLen.push_back(-1);
694  trkLenPerp.push_back(-1);
695  }
696  else
697  {
698  // calculate the total and the transverse track lengths and restrict the
699  // tracklength to be above the gas TPC track length threshold
700  double tracklen = 0.;
701  double tracklen_perp = 0.;
702 
703  //CAREFUL No offset for the trajectory points (origin for them is the TPC?)??????
704  //TODO check if the mcp point is within the TPC volume! Skip for mcp in the ECAL (showers)
705  //TODO Link showers to original mcp?
706  for(size_t itraj = 1; itraj < TrajMCPX->size(); itraj++)
707  {
708  //check that it is the correct mcp
709  if(TrajMCPTrajIndex->at(itraj) == mctrackid)
710  {
711  //Traj point+1
712  TVector3 point(TrajMCPX->at(itraj) - _util->GetOrigin()[0], TrajMCPY->at(itraj) - _util->GetOrigin()[1], TrajMCPZ->at(itraj) - _util->GetOrigin()[2]);
713 
714  //point is not in the TPC anymore - stop traj loop
715  if(not _util->PointInTPC(point))
716  {
717  // // std::cout << "Point not within the TPC: " << point.X() << " r " << std::sqrt(point.Y()*point.Y() + point.Z()*point.Z()) << std::endl;
718  continue;
719  }
720 
721  // find the length of the track by getting the distance between each hit
722  TVector3 diff(TrajMCPX->at(itraj) - TrajMCPX->at(itraj-1), TrajMCPY->at(itraj) - TrajMCPY->at(itraj-1), TrajMCPZ->at(itraj) - TrajMCPZ->at(itraj-1));
723  // perp length
724  TVector2 tracklen_perp_vec(TrajMCPZ->at(itraj) - TrajMCPZ->at(itraj-1), TrajMCPY->at(itraj) - TrajMCPY->at(itraj-1));
725  // Summing up
726  tracklen += diff.Mag();
727  tracklen_perp += tracklen_perp_vec.Mod();
728  }
729  }
730 
731  trkLen.push_back(tracklen);
732  trkLenPerp.push_back(tracklen_perp);
733  }
734 
735  //end track length
736  //***************************************************************************************************************/
737 
738  TVector3 xhat(1, 0, 0);
739  // float pz = mcp.Z();
740  // float pt = (mcp.Cross(xhat)).Mag();
741  // float px = mcp.X();
742  // float py = mcp.Y();
743  //float mctrackid = MCPTrkID->at(i);
744  // angle with respect to the incoming neutrino
745  float angle = atan(mcp.X() / mcp.Z());
746  float ecaltime = _util->GaussianSmearing(MCPTime->at(i), ECAL_time_resolution);
747  float time = MCPTime->at(i);
748 
749  //Check where start point is for the mcp
750  isFidStart.push_back(_util->PointInFiducial(spoint));
751  isTPCStart.push_back(_util->PointInTPC(spoint));
752  isCaloStart.push_back(_util->PointInCalo(spoint));
753  isThroughCaloStart.push_back(_util->isThroughCalo(spoint));
754  isInBetweenStart.push_back(_util->PointStopBetween(spoint));
755  isBarrelStart.push_back(_util->isBarrel(spoint));
756  isEndcapStart.push_back(_util->isEndcap(spoint));
757  //Check where endpoint of mcp is
758  isFidEnd.push_back(_util->PointInFiducial(epoint));
759  isTPCEnd.push_back(_util->PointInTPC(epoint));
760  isCaloEnd.push_back(_util->PointInCalo(epoint));
761  isThroughCaloEnd.push_back(_util->isThroughCalo(epoint));
762  isInBetweenEnd.push_back(_util->PointStopBetween(epoint));
763  isBarrelEnd.push_back(_util->isBarrel(epoint));
764  isEndcapEnd.push_back(_util->isEndcap(epoint));
765 
766  //start tpc
767  //***************************************************************************************************************/
768 
769  //Visible in the TPC
770  if( trkLen.at(nFSP) > gastpc_len )
771  {
772  for (int pidm = 0; pidm < 6; ++pidm)
773  {
774  if ( abs(pdg) == pdg_charged.at(pidm) )
775  {
776  // // std::cout << "Entered reco TPC" << std::endl;
777 
778  //Use range instead of Gluckstern for stopping tracks
779  //TODO is that correct? What if it is a scatter in the TPC? Need to check if daughter is same particle
780  float preco = 0;
781 
782  // save the true PDG, parametrized PID comes later
783  truepdg.push_back(pdg);
784  truepx.push_back(MCPStartPX->at(i));
785  truepy.push_back(MCPStartPY->at(i));
786  truepz.push_back(MCPStartPZ->at(i));
787  if(_correct4origin){
788  _MCPStartX.push_back(MCPStartX->at(i) - _util->GetOrigin()[0]);
789  _MCPStartY.push_back(MCPStartY->at(i) - _util->GetOrigin()[1]);
790  _MCPStartZ.push_back(MCPStartZ->at(i) - _util->GetOrigin()[2]);
791  _MCPEndX.push_back(MCPEndX->at(i) - _util->GetOrigin()[0]);
792  _MCPEndY.push_back(MCPEndY->at(i) - _util->GetOrigin()[1]);
793  _MCPEndZ.push_back(MCPEndZ->at(i) - _util->GetOrigin()[2]);
794  } else {
795  _MCPStartX.push_back(MCPStartX->at(i));
796  _MCPStartY.push_back(MCPStartY->at(i));
797  _MCPStartZ.push_back(MCPStartZ->at(i));
798  _MCPEndX.push_back(MCPEndX->at(i));
799  _MCPEndY.push_back(MCPEndY->at(i));
800  _MCPEndZ.push_back(MCPEndZ->at(i));
801  }
802  pdgmother.push_back(PDGMother->at(i));
803  // save the true momentum
804  truep.push_back(ptrue);
805  // save the true angle
806  _angle.push_back(angle);
807  //Save MC process
808  _MCProc.push_back(mcp_process);
809  _MCEndProc.push_back(mcp_endprocess);
810  mctime.push_back(time);
811  mctrkid.push_back(MCPTrkID->at(i));
812  motherid.push_back(MCMotherTrkID->at(i));
813 
814  //Case for range, the end point of the mcp is in the TPC, does not reach the ecal
815  if( _util->PointInTPC(epoint) )
816  {
817  // calculate number of trackpoints
818  float nHits = round (trkLen.at(nFSP) / gastpc_padPitch);
819  // angular resolution first term
820  float sigma_angle_1 = ((sigma_x * sigma_x * 0.0001) / trkLen.at(nFSP)*trkLen.at(nFSP)*0.0001) * (12*(nHits-1))/(nHits*(nHits+1));
821  // scattering term in Gluckstern formula
822  float sigma_angle_2 = (0.015*0.015 / (3. * ptrue * ptrue)) * (trkLen.at(nFSP)/gastpc_X0);
823  // angular resolution from the two terms above
824  float sigma_angle_short = sqrt(sigma_angle_1 + sigma_angle_2);
825  //reconstructed angle
826  float angle_reco = _util->GaussianSmearing(angle, sigma_angle_short);
827  //reconstructed momentum
828  preco = _util->GaussianSmearing( ptrue, sigmaP_short );
829 
830  if(preco > 0)
831  _preco.push_back(preco);
832  else
833  _preco.push_back(-1);
834 
835  anglereco.push_back(angle_reco);
836 
837  erecon.push_back(-1);
838  recopidecal.push_back(-1);
839  etime.push_back(-1);
840  detected.push_back(-1);
841  }
842  else
843  {
844  //Case where the endpoint is not in the TPC, should be able to use the Gluckstern formula
845  // calculate number of trackpoints
846  // float nHits = round (trkLen.at(nFSP) / gastpc_padPitch);
847  // // measurement term in Gluckstern formula
848  // float fracSig_meas = sqrt(720./(nHits+4)) * ((0.01*gastpc_padPitch*ptrue) / (0.3 * gastpc_B * 0.0001 *trkLenPerp.at(nFSP)*trkLenPerp.at(nFSP)));
849  // // multiple Coulomb scattering term in Gluckstern formula
850  // float fracSig_MCS = (0.052*sqrt(1.43)) / (gastpc_B * sqrt(gastpc_X0*trkLenPerp.at(nFSP)*0.0001));
851  // // momentum resoltion from the two terms above
852  // float sigmaP = ptrue * sqrt( fracSig_meas*fracSig_meas + fracSig_MCS*fracSig_MCS );
853  // // now Gaussian smear the true momentum using the momentum resolution
854  // preco = _util->GaussianSmearing( ptrue, sigmaP );
855 
856  //Vivek Gluckstern
857  int nHits = round (trkLen.at(nFSP) / gastpc_padPitch);
858  double ratio = 0.;
859  float sigmaPt = pypz*calcGluck(sigma_x, gastpc_B, gastpc_X0, nHits, pypz, trkLenPerp.at(nFSP), ratio);
860  float tan_lambda = MCPStartPX->at(i)/pypz;
861  float lambda = atan2(MCPStartPX->at(i), pypz); // between -pi and +pi
862  float del_lambda = 0.0062;
863  float temp = pypz*tan_lambda*del_lambda;
864  float sigmaP = (1./cos(lambda))*sqrt(sigmaPt*sigmaPt + temp*temp);
865  // now Gaussian smear the true momentum using the momentum resolution
866  preco = _util->GaussianSmearing( ptrue, sigmaP );
867 
868  // measurement term in the Gluckstern formula for calculating the
869  // angular resolution
870  float sigma_angle_1 = ((sigma_x * sigma_x * 0.0001) / trkLen.at(nFSP)*trkLen.at(nFSP)*0.0001) * (12*(nHits-1))/(nHits*(nHits+1));
871  // scattering term in Gluckstern formula
872  float sigma_angle_2 = (0.015*0.015 / (3. * ptrue * ptrue)) * (trkLen.at(nFSP)/gastpc_X0);
873  // angular resolution from the two terms above
874  float sigma_angle = sqrt(sigma_angle_1 + sigma_angle_2);
875  // now Gaussian smear the true angle using the angular resolution
876  float angle_reco = _util->GaussianSmearing(angle, sigma_angle);
877 
878  // save reconstructed momentum and angle to cafanatree
879  if(preco > 0)
880  _preco.push_back(preco);
881  else
882  _preco.push_back(-1);
883 
884  anglereco.push_back(angle_reco);
885 
886  //Reaches the ECAL and stops there
887  if( _util->PointInCalo(epoint) )
888  {
889  //Need energy measurement in ecal
890  TParticlePDG *part = TDatabasePDG::Instance()->GetParticle(abs(pdg));
891  if(nullptr == part)
892  {
893  // std::cout << "Could not find particle in root pdg table, pdg " << pdg << std::endl;
894  //deuteron
895  if( pdg == 1000010020 ) {
896  float mass = 1.8756;//in GeV mass deuteron
897  float etrue = std::sqrt(ptrue*ptrue + mass*mass) - mass;//needs to be KE here
898  float ECAL_resolution = fRes->Eval(etrue)*etrue;
899  float ereco = _util->GaussianSmearing(etrue, ECAL_resolution);
900  erecon.push_back(ereco);
901  recopidecal.push_back(-1);
902  detected.push_back(1);
903  etime.push_back(ecaltime);
904  }
905  else
906  {
907  erecon.push_back(-1);
908  recopidecal.push_back(-1);
909  detected.push_back(0);
910  etime.push_back(-1);
911  }
912  }
913  else
914  {
915  //by default should be tagged as an electron as it has a track,
916  //otherwise tag as gamma if not track -> need mis association rate, and use dE/dX in Scintillator?
917  //separation between e and mu/pi should be around 100%
918  //separation mu/pi -> based on Chris study with only the ECAL (no Muon ID detector)
919  //separation with p and mu/pi/e ?? high energy -> confusion with mu/pi, low energy confusion with e
920  //using E/p to ID?
921  float mass = part->Mass();//in GeV
922  float etrue = std::sqrt(ptrue*ptrue + mass*mass);//Just energy here is fine
923  float ECAL_resolution = fRes->Eval(etrue)*etrue;
924  float ereco = _util->GaussianSmearing(etrue, ECAL_resolution);
925  erecon.push_back((ereco > 0) ? ereco : 0.);
926  detected.push_back(1);
927  etime.push_back(ecaltime);
928  // // std::cout << "E/p " << ereco/preco << " true pdg " << pdg << std::endl;
929 
930  //Electron
931  if( abs(pdg) == 11 ){
932  recopidecal.push_back(11);
933  }
934  else if( abs(pdg) == 13 || abs(pdg) == 211 )
935  {
936  //Muons and Pions
937  //ptrue < 480 MeV/c 100% separation
938  //80% from 480 to 750
939  //90% up to 750 to 900
940  //95% over 900
941  float random_number = _util->GetRamdomNumber();
942 
943  if(ptrue < 0.48) {
944  recopidecal.push_back(abs(pdg));//100% efficiency by range
945  }
946  else if(ptrue >= 0.48 && ptrue < 0.75)
947  {
948  //case muon
949  if(abs(pdg) == 13)
950  {
951  if(random_number > (1 - 0.8)) {
952  recopidecal.push_back(13);
953  }
954  else{
955  recopidecal.push_back(211);
956  }
957  }
958 
959  //case pion
960  if(abs(pdg) == 211)
961  {
962  if(random_number > (1 - 0.8)) {
963  recopidecal.push_back(211);
964  }
965  else{
966  recopidecal.push_back(13);
967  }
968  }
969  }
970  else if(ptrue >= 0.75 && ptrue < 0.9)
971  {
972  //case muon
973  if(abs(pdg) == 13){
974  if(random_number > (1 - 0.9)) {
975  recopidecal.push_back(13);
976  }
977  else{
978  recopidecal.push_back(211);
979  }
980  }
981  //case pion
982  if(abs(pdg) == 211) {
983  if(random_number > (1 - 0.9)) {
984  recopidecal.push_back(211);
985  }
986  else{
987  recopidecal.push_back(13);
988  }
989  }
990  }
991  else
992  {
993  //case muon
994  if(abs(pdg) == 13){
995  if(random_number > (1 - 0.95)) {
996  recopidecal.push_back(13);
997  }
998  else{
999  recopidecal.push_back(211);
1000  }
1001  }
1002  //case pion
1003  if(abs(pdg) == 211){
1004  if(random_number > (1 - 0.95)) {
1005  recopidecal.push_back(211);
1006  }
1007  else{
1008  recopidecal.push_back(13);
1009  }
1010  }
1011  }
1012  }
1013  else if( abs(pdg) == 2212 )
1014  {
1015  recopidecal.push_back(2212);//TODO for p/pi separation
1016  }
1017  else {
1018  recopidecal.push_back(-1);
1019  }
1020  }
1021  }
1022  else if( _util->isThroughCalo(epoint) )
1023  {
1024  //Case the endpoint is outside the CALO -> it went through the ECAL (mu/pi/p possible)
1025  //the ECAL will see 60 MIPs on average
1026  double Evis = (double)nLayers; //in MIP
1027  //Smearing to account for Gaussian detector noise (Landau negligible)
1028  Evis = _util->GaussianSmearing(Evis, ECAL_MIP_Res);
1029  //1 MIP = 0.814 MeV
1030  double Erec = Evis * MIP2GeV_factor;
1031  erecon.push_back((Erec > 0) ? Erec : 0.);
1032  etime.push_back(ecaltime);
1033  detected.push_back(1);
1034 
1035  //Muon/Pions/Protons are reco as Muons (without MuID detector)
1036  if( abs(pdg) == 13 || abs(pdg) == 211 || abs(pdg) == 2212 ) {
1037  recopidecal.push_back(13);
1038  }
1039  else{
1040  recopidecal.push_back(-1);
1041  }
1042  }
1043  else
1044  {
1045  //Does not reach the ECAL???
1046  erecon.push_back(-1);
1047  recopidecal.push_back(-1);
1048  etime.push_back(-1);
1049  detected.push_back(0);
1050  }
1051  } //end endpoint is not in TPC
1052 
1053  //end tpc
1054  //***************************************************************************************************************/
1055 
1056  //start pid
1057  //***************************************************************************************************************/
1058 
1059  //--------------------------------------------------------------------------
1060  // Start of PID Parametrization
1061  //--------------------------------------------------------------------------
1062 
1063  float p = preco;
1064  // read the PID parametrization ntuple from T. Junk
1065  // TString filename="pid.root";
1066 
1067  std::vector<double> vec;
1068  std::vector<std::string> pnamelist = {"#pi", "#mu", "p", "K", "d", "e"};
1069  std::vector<std::string> recopnamelist = {"#pi", "#mu", "p", "K", "d", "e"};
1070 
1071  int qclosest = 0;
1072  float dist = 100000000.;
1073 
1074  // // std::cout << "preco " << p << " qclosest " << qclosest << std::endl;
1075 
1076  for (int q = 0; q < 501; ++q)
1077  {
1078  //Check the title and the reco momentum take only the one that fits
1079  std::string fulltitle = m_pidinterp[q]->GetTitle();
1080  unsigned first = fulltitle.find("=");
1081  unsigned last = fulltitle.find("GeV");
1082  std::string substr = fulltitle.substr(first+1, last - first-1);
1083  float pidinterp_mom = std::atof(substr.c_str());
1084  //calculate the distance between the bin and mom, store the q the closest
1085  float disttemp = std::abs(pidinterp_mom - p);
1086  //// std::cout << disttemp << " " << dist << std::endl;
1087 
1088  // // std::cout << "preco " << p << " ptitle " << pidinterp_mom << " dist " << disttemp << " q " << q << std::endl;
1089 
1090  if( disttemp < dist ){
1091  dist = disttemp;
1092  qclosest = q;
1093  // // std::cout << "pid mom " << pidinterp_mom << " reco mom " << p << " dist " << dist << " qclosest " << qclosest << std::endl;
1094  }
1095  } // closes the "pidmatrix" loop
1096 
1097  // // std::cout << "Started pid" << std::endl;
1098  //loop over the columns (true pid)
1099  std::vector< std::pair<float, std::string> > v_prob;
1100 
1101  // // std::cout << "pidm " << pidm << std::endl;
1102  //get true particle name
1103  std::string trueparticlename = m_pidinterp[qclosest]->GetXaxis()->GetBinLabel(pidm+1);
1104  // // std::cout << trueparticlename << std::endl;
1105  if ( trueparticlename == pnamelist[pidm] )
1106  {
1107  //loop over the rows (reco pid)
1108  for (int pidr = 0; pidr < 6; ++pidr)
1109  {
1110  // // std::cout << "pidr " << pidr << std::endl;
1111  std::string recoparticlename = m_pidinterp[qclosest]->GetYaxis()->GetBinLabel(pidr+1);
1112  if (recoparticlename == recopnamelist[pidr])
1113  {
1114  float prob = m_pidinterp[qclosest]->GetBinContent(pidm+1,pidr+1);
1115  prob_arr.push_back(prob);
1116 
1117  // // std::cout << "true part " << trueparticlename << " true pid " << pdg_charged.at(pidm) << " reco name " << recoparticlename << " reco part list "
1118  // << recopnamelist[pidr] << " true mom " << ptrue << " reco mom " << p << " prob " << pidinterp->GetBinContent(pidm+1,pidr+1) << '\n';
1119  //Need to check random number value and prob value then associate the recopdg to the reco prob
1120  v_prob.push_back( std::make_pair(prob, recoparticlename) );
1121  }
1122  }
1123 
1124  int pid = -1;
1125  if(v_prob.size() > 1)
1126  {
1127  //Order the vector of prob
1128  std::sort(v_prob.begin(), v_prob.end());
1129  //Throw a random number between 0 and 1
1130  float random_number = _util->GetRamdomNumber();
1131  //Make cumulative sum to get the range
1132  std::partial_sum(v_prob.begin(), v_prob.end(), v_prob.begin(), [](const P& _x, const P& _y){return P(_x.first + _y.first, _y.second);});
1133 
1134  for(size_t ivec = 0; ivec < v_prob.size()-1; ivec++)
1135  {
1136  if( random_number < v_prob.at(ivec+1).first && random_number >= v_prob.at(ivec).first )
1137  {
1138  pid = pdg_charged.at( std::distance( recopnamelist.begin(), std::find(recopnamelist.begin(), recopnamelist.end(), v_prob.at(ivec+1).second) ) );
1139  }
1140  }
1141  }
1142  else
1143  {
1144  pid = pdg_charged.at( std::distance( recopnamelist.begin(), std::find(recopnamelist.begin(), recopnamelist.end(), v_prob.at(0).second) ) );
1145  }
1146 
1147  recopid.push_back( pid );
1148  } // closes the if statement
1149 
1150  //end pid
1151  //***************************************************************************************************************/
1152 
1153  } // closes the conditional statement of trueparticlename == MC true pdg
1154  else
1155  {
1156  //not in the pdglist of particles but visible in TPC?
1157  auto found = std::find(pdg_charged.begin(), pdg_charged.end(), abs(pdg));
1158  if(found == pdg_charged.end())
1159  {
1160  // // std::cout << "Maybe visible but not {#pi, #mu, p, K, d, e};" << std::endl;
1161  // // std::cout << "pdg " << pdg << std::endl;
1162 
1163  truepdg.push_back(pdg);
1164  detected.push_back(0);
1165  truepx.push_back(MCPStartPX->at(i));
1166  truepy.push_back(MCPStartPY->at(i));
1167  truepz.push_back(MCPStartPZ->at(i));
1168  if(_correct4origin){
1169  _MCPStartX.push_back(MCPStartX->at(i) - _util->GetOrigin()[0]);
1170  _MCPStartY.push_back(MCPStartY->at(i) - _util->GetOrigin()[1]);
1171  _MCPStartZ.push_back(MCPStartZ->at(i) - _util->GetOrigin()[2]);
1172  _MCPEndX.push_back(MCPEndX->at(i) - _util->GetOrigin()[0]);
1173  _MCPEndY.push_back(MCPEndY->at(i) - _util->GetOrigin()[1]);
1174  _MCPEndZ.push_back(MCPEndZ->at(i) - _util->GetOrigin()[2]);
1175  } else {
1176  _MCPStartX.push_back(MCPStartX->at(i));
1177  _MCPStartY.push_back(MCPStartY->at(i));
1178  _MCPStartZ.push_back(MCPStartZ->at(i));
1179  _MCPEndX.push_back(MCPEndX->at(i));
1180  _MCPEndY.push_back(MCPEndY->at(i));
1181  _MCPEndZ.push_back(MCPEndZ->at(i));
1182  }
1183  pdgmother.push_back(PDGMother->at(i));
1184  // save the true momentum
1185  truep.push_back(ptrue);
1186  // save the true angle
1187  _angle.push_back(angle);
1188  //Save MC process
1189  _MCProc.push_back(mcp_process);
1190  _MCEndProc.push_back(mcp_endprocess);
1191  mctime.push_back(time);
1192  etime.push_back(-1);
1193  erecon.push_back(-1);
1194  _preco.push_back(-1);
1195  anglereco.push_back(-1);
1196  recopid.push_back(-1);
1197  recopidecal.push_back(-1);
1198  for (int pidr = 0; pidr < 6; ++pidr) prob_arr.push_back(-1);
1199  mctrkid.push_back(MCPTrkID->at(i));
1200  motherid.push_back(MCMotherTrkID->at(i));
1201  break; //break out of the loop over the vertical bining!
1202  }
1203  }
1204  } // closes the vertical bining loop of the pid matrix
1205  }//close if track_length > tpc_min_length
1206  else
1207  {
1208  //Not visible in the TPC
1209  //for neutrons
1210  if(std::abs(pdg) == 2112)
1211  {
1212  //start neutrons
1213  //***************************************************************************************************************/
1214 
1215  if(_util->PointInCalo(epoint)) //needs to stop in the ECAL
1216  {
1217  //check if it can be detected by the ECAL
1218  //Assumes 40% efficiency to detect
1219  float random_number = _util->GetRamdomNumber();
1220  float true_KE = std::sqrt(ptrue*ptrue + neutron_mass*neutron_mass) - neutron_mass;//KE here
1221  // float true_KE = ptrue*ptrue / (2*neutron_mass); // in GeV
1222  int index = (true_KE >= 0.05) ? 1 : 0;
1223 
1224  // // std::cout << "KE " << true_KE << " index " << index << " 1 - eff " << 1-NeutronECAL_detEff[index] << " rdnm " << random_number << std::endl;
1225 
1226  if(random_number > (1 - NeutronECAL_detEff[index]) && true_KE > 0.003)//Threshold of 3 MeV
1227  {
1228  //TODO random is first interaction or rescatter and smear accordingly to Chris's study
1229  //Detected in the ECAL
1230  // recopid.push_back(2112);
1231  recopid.push_back(-1); //reco pid set to 0?
1232  detected.push_back(1);
1233  float eres = sigmaNeutronECAL_first * true_KE;
1234  float ereco_KE = _util->GaussianSmearing( true_KE, eres );
1235  float ereco = ereco_KE + neutron_mass;
1236  erecon.push_back(ereco > 0 ? ereco : 0.);
1237  // // std::cout << "true part n true energy " << std::sqrt(ptrue*ptrue + neutron_mass*neutron_mass) << " ereco " << erecon[i] << std::endl;
1238  truepdg.push_back(pdg);
1239  truep.push_back(ptrue);
1240  truepx.push_back(MCPStartPX->at(i));
1241  truepy.push_back(MCPStartPY->at(i));
1242  truepz.push_back(MCPStartPZ->at(i));
1243  if(_correct4origin){
1244  _MCPStartX.push_back(MCPStartX->at(i) - _util->GetOrigin()[0]);
1245  _MCPStartY.push_back(MCPStartY->at(i) - _util->GetOrigin()[1]);
1246  _MCPStartZ.push_back(MCPStartZ->at(i) - _util->GetOrigin()[2]);
1247  _MCPEndX.push_back(MCPEndX->at(i) - _util->GetOrigin()[0]);
1248  _MCPEndY.push_back(MCPEndY->at(i) - _util->GetOrigin()[1]);
1249  _MCPEndZ.push_back(MCPEndZ->at(i) - _util->GetOrigin()[2]);
1250  } else {
1251  _MCPStartX.push_back(MCPStartX->at(i));
1252  _MCPStartY.push_back(MCPStartY->at(i));
1253  _MCPStartZ.push_back(MCPStartZ->at(i));
1254  _MCPEndX.push_back(MCPEndX->at(i));
1255  _MCPEndY.push_back(MCPEndY->at(i));
1256  _MCPEndZ.push_back(MCPEndZ->at(i));
1257  }
1258  pdgmother.push_back(PDGMother->at(i));
1259  //Save MC process
1260  _MCProc.push_back(mcp_process);
1261  _MCEndProc.push_back(mcp_endprocess);
1262  mctime.push_back(time);
1263  etime.push_back(ecaltime);
1264  _angle.push_back(angle);
1265  _preco.push_back(-1);
1266  anglereco.push_back(-1);
1267  recopidecal.push_back(2112);
1268  for (int pidr = 0; pidr < 6; ++pidr) prob_arr.push_back(-1);
1269  mctrkid.push_back(MCPTrkID->at(i));
1270  motherid.push_back(MCMotherTrkID->at(i));
1271  }
1272  else
1273  {
1274  //neutron not detected
1275  detected.push_back(0);
1276  truep.push_back(ptrue);
1277  recopid.push_back(-1);
1278  erecon.push_back(-1);
1279  truepdg.push_back(pdg);
1280  truepx.push_back(MCPStartPX->at(i));
1281  truepy.push_back(MCPStartPY->at(i));
1282  truepz.push_back(MCPStartPZ->at(i));
1283  if(_correct4origin){
1284  _MCPStartX.push_back(MCPStartX->at(i) - _util->GetOrigin()[0]);
1285  _MCPStartY.push_back(MCPStartY->at(i) - _util->GetOrigin()[1]);
1286  _MCPStartZ.push_back(MCPStartZ->at(i) - _util->GetOrigin()[2]);
1287  _MCPEndX.push_back(MCPEndX->at(i) - _util->GetOrigin()[0]);
1288  _MCPEndY.push_back(MCPEndY->at(i) - _util->GetOrigin()[1]);
1289  _MCPEndZ.push_back(MCPEndZ->at(i) - _util->GetOrigin()[2]);
1290  } else {
1291  _MCPStartX.push_back(MCPStartX->at(i));
1292  _MCPStartY.push_back(MCPStartY->at(i));
1293  _MCPStartZ.push_back(MCPStartZ->at(i));
1294  _MCPEndX.push_back(MCPEndX->at(i));
1295  _MCPEndY.push_back(MCPEndY->at(i));
1296  _MCPEndZ.push_back(MCPEndZ->at(i));
1297  }
1298  pdgmother.push_back(PDGMother->at(i));
1299  //Save MC process
1300  _MCProc.push_back(mcp_process);
1301  _MCEndProc.push_back(mcp_endprocess);
1302  mctime.push_back(time);
1303  etime.push_back(-1);
1304  _angle.push_back(angle);
1305  _preco.push_back(-1);
1306  anglereco.push_back(-1);
1307  recopidecal.push_back(-1);
1308  for (int pidr = 0; pidr < 6; ++pidr) prob_arr.push_back(-1);
1309  mctrkid.push_back(MCPTrkID->at(i));
1310  motherid.push_back(MCMotherTrkID->at(i));
1311  }
1312  } //endpoint is in ECAL
1313  else
1314  {
1315  //Endpoint is not in calo (TPC/isInBetween or outside Calo)
1316  detected.push_back(0);
1317  truep.push_back(ptrue);
1318  recopid.push_back(-1);
1319  erecon.push_back(-1);
1320  truepdg.push_back(pdg);
1321  truepx.push_back(MCPStartPX->at(i));
1322  truepy.push_back(MCPStartPY->at(i));
1323  truepz.push_back(MCPStartPZ->at(i));
1324  if(_correct4origin){
1325  _MCPStartX.push_back(MCPStartX->at(i) - _util->GetOrigin()[0]);
1326  _MCPStartY.push_back(MCPStartY->at(i) - _util->GetOrigin()[1]);
1327  _MCPStartZ.push_back(MCPStartZ->at(i) - _util->GetOrigin()[2]);
1328  _MCPEndX.push_back(MCPEndX->at(i) - _util->GetOrigin()[0]);
1329  _MCPEndY.push_back(MCPEndY->at(i) - _util->GetOrigin()[1]);
1330  _MCPEndZ.push_back(MCPEndZ->at(i) - _util->GetOrigin()[2]);
1331  } else {
1332  _MCPStartX.push_back(MCPStartX->at(i));
1333  _MCPStartY.push_back(MCPStartY->at(i));
1334  _MCPStartZ.push_back(MCPStartZ->at(i));
1335  _MCPEndX.push_back(MCPEndX->at(i));
1336  _MCPEndY.push_back(MCPEndY->at(i));
1337  _MCPEndZ.push_back(MCPEndZ->at(i));
1338  }
1339  pdgmother.push_back(PDGMother->at(i));
1340  //Save MC process
1341  _MCProc.push_back(mcp_process);
1342  _MCEndProc.push_back(mcp_endprocess);
1343  mctime.push_back(time);
1344  etime.push_back(-1);
1345  _angle.push_back(angle);
1346  _preco.push_back(-1);
1347  anglereco.push_back(-1);
1348  recopidecal.push_back(-1);
1349  for (int pidr = 0; pidr < 6; ++pidr) prob_arr.push_back(-1);
1350  mctrkid.push_back(MCPTrkID->at(i));
1351  motherid.push_back(MCMotherTrkID->at(i));
1352  }
1353  //End neutrons
1354  //***************************************************************************************************************/
1355  }
1356  else if(std::abs(pdg) == 111) //for pi0s
1357  {
1358  //start pi0s
1359  //***************************************************************************************************************/
1360  erecon.push_back(-1);
1361  recopid.push_back(-1);
1362  detected.push_back(0);
1363  recopidecal.push_back(-1);
1364 
1365  truep.push_back(ptrue);
1366  truepdg.push_back(pdg);
1367  truepx.push_back(MCPStartPX->at(i));
1368  truepy.push_back(MCPStartPY->at(i));
1369  truepz.push_back(MCPStartPZ->at(i));
1370  if(_correct4origin){
1371  _MCPStartX.push_back(MCPStartX->at(i) - _util->GetOrigin()[0]);
1372  _MCPStartY.push_back(MCPStartY->at(i) - _util->GetOrigin()[1]);
1373  _MCPStartZ.push_back(MCPStartZ->at(i) - _util->GetOrigin()[2]);
1374  _MCPEndX.push_back(MCPEndX->at(i) - _util->GetOrigin()[0]);
1375  _MCPEndY.push_back(MCPEndY->at(i) - _util->GetOrigin()[1]);
1376  _MCPEndZ.push_back(MCPEndZ->at(i) - _util->GetOrigin()[2]);
1377  } else {
1378  _MCPStartX.push_back(MCPStartX->at(i));
1379  _MCPStartY.push_back(MCPStartY->at(i));
1380  _MCPStartZ.push_back(MCPStartZ->at(i));
1381  _MCPEndX.push_back(MCPEndX->at(i));
1382  _MCPEndY.push_back(MCPEndY->at(i));
1383  _MCPEndZ.push_back(MCPEndZ->at(i));
1384  }
1385  pdgmother.push_back(PDGMother->at(i));
1386  _angle.push_back(angle);
1387  //Save MC process
1388  _MCProc.push_back(mcp_process);
1389  _MCEndProc.push_back(mcp_endprocess);
1390  mctime.push_back(time);
1391  etime.push_back(-1);
1392  _preco.push_back(-1);
1393  anglereco.push_back(-1);
1394 
1395  for (int pidr = 0; pidr < 6; ++pidr) prob_arr.push_back(-1);
1396 
1397  mctrkid.push_back(MCPTrkID->at(i));
1398  motherid.push_back(MCMotherTrkID->at(i));
1399 
1400  //end pi0s
1401  //***************************************************************************************************************/
1402  }
1403  else if(std::abs(pdg) == 22)//for gammas
1404  {
1405  //start gammas
1406  //***************************************************************************************************************/
1407 
1408  if( PDGMother->at(i) != 111 )
1409  {
1410  //Endpoint is in the ECAL
1411  if(_util->PointInCalo(epoint))
1412  {
1413  //if they hit the ECAL and smear their energy
1414  float ECAL_resolution = fRes->Eval(ptrue)*ptrue;
1415  float ereco = _util->GaussianSmearing(ptrue, ECAL_resolution);
1416  erecon.push_back( (ereco > 0) ? ereco : 0. );
1417  recopid.push_back(-1);
1418  detected.push_back(1);
1419 
1420  truep.push_back(ptrue);
1421  truepdg.push_back(pdg);
1422  truepx.push_back(MCPStartPX->at(i));
1423  truepy.push_back(MCPStartPY->at(i));
1424  truepz.push_back(MCPStartPZ->at(i));
1425  if(_correct4origin){
1426  _MCPStartX.push_back(MCPStartX->at(i) - _util->GetOrigin()[0]);
1427  _MCPStartY.push_back(MCPStartY->at(i) - _util->GetOrigin()[1]);
1428  _MCPStartZ.push_back(MCPStartZ->at(i) - _util->GetOrigin()[2]);
1429  _MCPEndX.push_back(MCPEndX->at(i) - _util->GetOrigin()[0]);
1430  _MCPEndY.push_back(MCPEndY->at(i) - _util->GetOrigin()[1]);
1431  _MCPEndZ.push_back(MCPEndZ->at(i) - _util->GetOrigin()[2]);
1432  } else {
1433  _MCPStartX.push_back(MCPStartX->at(i));
1434  _MCPStartY.push_back(MCPStartY->at(i));
1435  _MCPStartZ.push_back(MCPStartZ->at(i));
1436  _MCPEndX.push_back(MCPEndX->at(i));
1437  _MCPEndY.push_back(MCPEndY->at(i));
1438  _MCPEndZ.push_back(MCPEndZ->at(i));
1439  }
1440  pdgmother.push_back(PDGMother->at(i));
1441  _angle.push_back(angle);
1442  //Save MC process
1443  _MCProc.push_back(mcp_process);
1444  _MCEndProc.push_back(mcp_endprocess);
1445  mctime.push_back(time);
1446  etime.push_back(ecaltime);
1447  _preco.push_back(-1);
1448  anglereco.push_back(-1);
1449 
1450  //reach the ECAL, should be tagged as gamma
1451  recopidecal.push_back(22);
1452 
1453  for (int pidr = 0; pidr < 6; ++pidr) prob_arr.push_back(-1);
1454  mctrkid.push_back(MCPTrkID->at(i));
1455  motherid.push_back(MCMotherTrkID->at(i));
1456  }
1457  else if(_util->PointInTPC(epoint) || _util->PointStopBetween(epoint) || _util->isThroughCalo(epoint))
1458  {
1459  //case endpoint is in the TPC (Converted!) or in between the TPC/ECAL
1460  erecon.push_back(-1);
1461  recopid.push_back(-1);
1462  detected.push_back(0);
1463 
1464  truep.push_back(ptrue);
1465  truepdg.push_back(pdg);
1466  truepx.push_back(MCPStartPX->at(i));
1467  truepy.push_back(MCPStartPY->at(i));
1468  truepz.push_back(MCPStartPZ->at(i));
1469  if(_correct4origin){
1470  _MCPStartX.push_back(MCPStartX->at(i) - _util->GetOrigin()[0]);
1471  _MCPStartY.push_back(MCPStartY->at(i) - _util->GetOrigin()[1]);
1472  _MCPStartZ.push_back(MCPStartZ->at(i) - _util->GetOrigin()[2]);
1473  _MCPEndX.push_back(MCPEndX->at(i) - _util->GetOrigin()[0]);
1474  _MCPEndY.push_back(MCPEndY->at(i) - _util->GetOrigin()[1]);
1475  _MCPEndZ.push_back(MCPEndZ->at(i) - _util->GetOrigin()[2]);
1476  } else {
1477  _MCPStartX.push_back(MCPStartX->at(i));
1478  _MCPStartY.push_back(MCPStartY->at(i));
1479  _MCPStartZ.push_back(MCPStartZ->at(i));
1480  _MCPEndX.push_back(MCPEndX->at(i));
1481  _MCPEndY.push_back(MCPEndY->at(i));
1482  _MCPEndZ.push_back(MCPEndZ->at(i));
1483  }
1484  pdgmother.push_back(PDGMother->at(i));
1485  _angle.push_back(angle);
1486  //Save MC process
1487  _MCProc.push_back(mcp_process);
1488  _MCEndProc.push_back(mcp_endprocess);
1489  mctime.push_back(time);
1490  etime.push_back(-1);
1491  _preco.push_back(-1);
1492  anglereco.push_back(-1);
1493  //converted so not seen in ECAL
1494  recopidecal.push_back(-1);
1495  for (int pidr = 0; pidr < 6; ++pidr) prob_arr.push_back(-1);
1496  mctrkid.push_back(MCPTrkID->at(i));
1497  motherid.push_back(MCMotherTrkID->at(i));
1498  }
1499  else{
1500  }
1501  }
1502  else
1503  {
1504  //case they are from pi0
1505  //Endpoint is not in the tracker, reaches the ecal
1506  if(_util->PointInCalo(epoint))
1507  {
1508  //if they hit the ECAL and smear their energy
1509  float ECAL_resolution = fRes->Eval(ptrue)*ptrue;
1510  float ereco = _util->GaussianSmearing(ptrue, ECAL_resolution);
1511  erecon.push_back((ereco > 0) ? ereco : 0.);
1512  recopid.push_back(-1);
1513  detected.push_back(1);
1514 
1515  truep.push_back(ptrue);
1516  truepdg.push_back(pdg);
1517  truepx.push_back(MCPStartPX->at(i));
1518  truepy.push_back(MCPStartPY->at(i));
1519  truepz.push_back(MCPStartPZ->at(i));
1520  if(_correct4origin){
1521  _MCPStartX.push_back(MCPStartX->at(i) - _util->GetOrigin()[0]);
1522  _MCPStartY.push_back(MCPStartY->at(i) - _util->GetOrigin()[1]);
1523  _MCPStartZ.push_back(MCPStartZ->at(i) - _util->GetOrigin()[2]);
1524  _MCPEndX.push_back(MCPEndX->at(i) - _util->GetOrigin()[0]);
1525  _MCPEndY.push_back(MCPEndY->at(i) - _util->GetOrigin()[1]);
1526  _MCPEndZ.push_back(MCPEndZ->at(i) - _util->GetOrigin()[2]);
1527  } else {
1528  _MCPStartX.push_back(MCPStartX->at(i));
1529  _MCPStartY.push_back(MCPStartY->at(i));
1530  _MCPStartZ.push_back(MCPStartZ->at(i));
1531  _MCPEndX.push_back(MCPEndX->at(i));
1532  _MCPEndY.push_back(MCPEndY->at(i));
1533  _MCPEndZ.push_back(MCPEndZ->at(i));
1534  }
1535  pdgmother.push_back(PDGMother->at(i));
1536  _angle.push_back(angle);
1537  //Save MC process
1538  _MCProc.push_back(mcp_process);
1539  _MCEndProc.push_back(mcp_endprocess);
1540  mctime.push_back(time);
1541  etime.push_back(ecaltime);
1542  _preco.push_back(-1);
1543  anglereco.push_back(-1);
1544 
1545  //reaches the ecal
1546  recopidecal.push_back(22);
1547 
1548  for (int pidr = 0; pidr < 6; ++pidr) prob_arr.push_back(-1);
1549  mctrkid.push_back(MCPTrkID->at(i));
1550  motherid.push_back(MCMotherTrkID->at(i));
1551  }
1552  else if(_util->PointInTPC(epoint) || _util->PointStopBetween(epoint) || _util->isThroughCalo(epoint))
1553  {
1554  //from pi0 and converted in TPC or stopped between TPC/ECAL
1555  erecon.push_back(-1);
1556  recopid.push_back(-1);
1557  detected.push_back(0);
1558 
1559  truep.push_back(ptrue);
1560  truepdg.push_back(pdg);
1561  truepx.push_back(MCPStartPX->at(i));
1562  truepy.push_back(MCPStartPY->at(i));
1563  truepz.push_back(MCPStartPZ->at(i));
1564  if(_correct4origin){
1565  _MCPStartX.push_back(MCPStartX->at(i) - _util->GetOrigin()[0]);
1566  _MCPStartY.push_back(MCPStartY->at(i) - _util->GetOrigin()[1]);
1567  _MCPStartZ.push_back(MCPStartZ->at(i) - _util->GetOrigin()[2]);
1568  _MCPEndX.push_back(MCPEndX->at(i) - _util->GetOrigin()[0]);
1569  _MCPEndY.push_back(MCPEndY->at(i) - _util->GetOrigin()[1]);
1570  _MCPEndZ.push_back(MCPEndZ->at(i) - _util->GetOrigin()[2]);
1571  } else {
1572  _MCPStartX.push_back(MCPStartX->at(i));
1573  _MCPStartY.push_back(MCPStartY->at(i));
1574  _MCPStartZ.push_back(MCPStartZ->at(i));
1575  _MCPEndX.push_back(MCPEndX->at(i));
1576  _MCPEndY.push_back(MCPEndY->at(i));
1577  _MCPEndZ.push_back(MCPEndZ->at(i));
1578  }
1579  pdgmother.push_back(PDGMother->at(i));
1580  _angle.push_back(angle);
1581  //Save MC process
1582  _MCProc.push_back(mcp_process);
1583  _MCEndProc.push_back(mcp_endprocess);
1584  mctime.push_back(time);
1585  etime.push_back(-1);
1586  _preco.push_back(-1);
1587  anglereco.push_back(-1);
1588  //converted not seen by ecal
1589  recopidecal.push_back(-1);
1590  for (int pidr = 0; pidr < 6; ++pidr) prob_arr.push_back(-1);
1591  mctrkid.push_back(MCPTrkID->at(i));
1592  motherid.push_back(MCMotherTrkID->at(i));
1593  }
1594  else{
1595  }
1596  }
1597  //end gammas
1598  //***************************************************************************************************************/
1599  }
1600  else if(std::find(neutrinos.begin(), neutrinos.end(), abs(pdg)) != neutrinos.end())
1601  {
1602  //Case for neutrinos from NC interactions
1603  truepdg.push_back(pdg);
1604  detected.push_back(0);
1605  truepx.push_back(MCPStartPX->at(i));
1606  truepy.push_back(MCPStartPY->at(i));
1607  truepz.push_back(MCPStartPZ->at(i));
1608  if(_correct4origin){
1609  _MCPStartX.push_back(MCPStartX->at(i) - _util->GetOrigin()[0]);
1610  _MCPStartY.push_back(MCPStartY->at(i) - _util->GetOrigin()[1]);
1611  _MCPStartZ.push_back(MCPStartZ->at(i) - _util->GetOrigin()[2]);
1612  _MCPEndX.push_back(MCPEndX->at(i) - _util->GetOrigin()[0]);
1613  _MCPEndY.push_back(MCPEndY->at(i) - _util->GetOrigin()[1]);
1614  _MCPEndZ.push_back(MCPEndZ->at(i) - _util->GetOrigin()[2]);
1615  } else {
1616  _MCPStartX.push_back(MCPStartX->at(i));
1617  _MCPStartY.push_back(MCPStartY->at(i));
1618  _MCPStartZ.push_back(MCPStartZ->at(i));
1619  _MCPEndX.push_back(MCPEndX->at(i));
1620  _MCPEndY.push_back(MCPEndY->at(i));
1621  _MCPEndZ.push_back(MCPEndZ->at(i));
1622  }
1623  pdgmother.push_back(PDGMother->at(i));
1624  // save the true momentum
1625  truep.push_back(ptrue);
1626  // save the true angle
1627  _angle.push_back(angle);
1628  //Save MC process
1629  _MCProc.push_back(mcp_process);
1630  _MCEndProc.push_back(mcp_endprocess);
1631  mctime.push_back(time);
1632 
1633  etime.push_back(-1);
1634  erecon.push_back(0.);
1635  recopidecal.push_back(-1);
1636  _preco.push_back(-1);
1637  anglereco.push_back(-1);
1638  recopid.push_back(-1);
1639  for (int pidr = 0; pidr < 6; ++pidr) prob_arr.push_back(-1);
1640 
1641  mctrkid.push_back(MCPTrkID->at(i));
1642  motherid.push_back(MCMotherTrkID->at(i));
1643  }
1644  else
1645  {
1646  //Case for particles that stop or go through ECAL (problematic particles with no track length????)
1647  //Not visible in the TPC and not neutron or gamma or pi0 (otherwise it has been already done above)
1648 
1649  if(_util->PointInCalo(epoint))
1650  {
1651  truepdg.push_back(pdg);
1652  detected.push_back(1);
1653  truepx.push_back(MCPStartPX->at(i));
1654  truepy.push_back(MCPStartPY->at(i));
1655  truepz.push_back(MCPStartPZ->at(i));
1656  if(_correct4origin){
1657  _MCPStartX.push_back(MCPStartX->at(i) - _util->GetOrigin()[0]);
1658  _MCPStartY.push_back(MCPStartY->at(i) - _util->GetOrigin()[1]);
1659  _MCPStartZ.push_back(MCPStartZ->at(i) - _util->GetOrigin()[2]);
1660  _MCPEndX.push_back(MCPEndX->at(i) - _util->GetOrigin()[0]);
1661  _MCPEndY.push_back(MCPEndY->at(i) - _util->GetOrigin()[1]);
1662  _MCPEndZ.push_back(MCPEndZ->at(i) - _util->GetOrigin()[2]);
1663  } else {
1664  _MCPStartX.push_back(MCPStartX->at(i));
1665  _MCPStartY.push_back(MCPStartY->at(i));
1666  _MCPStartZ.push_back(MCPStartZ->at(i));
1667  _MCPEndX.push_back(MCPEndX->at(i));
1668  _MCPEndY.push_back(MCPEndY->at(i));
1669  _MCPEndZ.push_back(MCPEndZ->at(i));
1670  }
1671  pdgmother.push_back(PDGMother->at(i));
1672  // save the true momentum
1673  truep.push_back(ptrue);
1674  // save the true angle
1675  _angle.push_back(angle);
1676  //Save MC process
1677  _MCProc.push_back(mcp_process);
1678  _MCEndProc.push_back(mcp_endprocess);
1679  mctime.push_back(time);
1680 
1681  etime.push_back(ecaltime);
1682 
1683  TParticlePDG *part = TDatabasePDG::Instance()->GetParticle(abs(pdg));
1684  float mass = 0.;
1685  if(nullptr != part)
1686  mass = part->Mass();//in GeV
1687 
1688  float etrue = std::sqrt(ptrue*ptrue + mass*mass);
1689  float ECAL_resolution = fRes->Eval(etrue)*etrue;
1690  float ereco = _util->GaussianSmearing(etrue, ECAL_resolution);
1691  erecon.push_back((ereco > 0) ? ereco : 0.);
1692 
1693  //Electron
1694  if( abs(pdg) == 11 ){
1695  recopidecal.push_back(11);
1696  }
1697  else if( abs(pdg) == 13 || abs(pdg) == 211 )
1698  {
1699  //Muons and Pions
1700  //ptrue < 480 MeV/c 100% separation
1701  //80% from 480 to 750
1702  //90% up to 750 to 900
1703  //95% over 900
1704  float random_number = _util->GetRamdomNumber();
1705 
1706  if(ptrue < 0.48) {
1707  recopidecal.push_back(abs(pdg));//100% efficiency by range
1708  }
1709  else if(ptrue >= 0.48 && ptrue < 0.75)
1710  {
1711  //case muon
1712  if(abs(pdg) == 13)
1713  {
1714  if(random_number > (1 - 0.8)) {
1715  recopidecal.push_back(13);
1716  }
1717  else{
1718  recopidecal.push_back(211);
1719  }
1720  }
1721 
1722  //case pion
1723  if(abs(pdg) == 211)
1724  {
1725  if(random_number > (1 - 0.8)) {
1726  recopidecal.push_back(211);
1727  }
1728  else{
1729  recopidecal.push_back(13);
1730  }
1731  }
1732  }
1733  else if(ptrue >= 0.75 && ptrue < 0.9)
1734  {
1735  //case muon
1736  if(abs(pdg) == 13){
1737  if(random_number > (1 - 0.9)) {
1738  recopidecal.push_back(13);
1739  }
1740  else{
1741  recopidecal.push_back(211);
1742  }
1743  }
1744  //case pion
1745  if(abs(pdg) == 211) {
1746  if(random_number > (1 - 0.9)) {
1747  recopidecal.push_back(211);
1748  }
1749  else{
1750  recopidecal.push_back(13);
1751  }
1752  }
1753  }
1754  else
1755  {
1756  //case muon
1757  if(abs(pdg) == 13){
1758  if(random_number > (1 - 0.95)) {
1759  recopidecal.push_back(13);
1760  }
1761  else{
1762  recopidecal.push_back(211);
1763  }
1764  }
1765  //case pion
1766  if(abs(pdg) == 211){
1767  if(random_number > (1 - 0.95)) {
1768  recopidecal.push_back(211);
1769  }
1770  else{
1771  recopidecal.push_back(13);
1772  }
1773  }
1774  }
1775  }
1776  else if( abs(pdg) == 2212 )
1777  {
1778  recopidecal.push_back(2212);//TODO for p/pi separation
1779  }
1780  else {
1781  recopidecal.push_back(-1);
1782  }
1783 
1784  _preco.push_back(-1);
1785  anglereco.push_back(-1);
1786  recopid.push_back(-1);
1787  for (int pidr = 0; pidr < 6; ++pidr) prob_arr.push_back(-1);
1788 
1789  mctrkid.push_back(MCPTrkID->at(i));
1790  motherid.push_back(MCMotherTrkID->at(i));
1791  }
1792  else if (_util->isThroughCalo(epoint))
1793  {
1794  truepdg.push_back(pdg);
1795  truepx.push_back(MCPStartPX->at(i));
1796  truepy.push_back(MCPStartPY->at(i));
1797  truepz.push_back(MCPStartPZ->at(i));
1798  if(_correct4origin){
1799  _MCPStartX.push_back(MCPStartX->at(i) - _util->GetOrigin()[0]);
1800  _MCPStartY.push_back(MCPStartY->at(i) - _util->GetOrigin()[1]);
1801  _MCPStartZ.push_back(MCPStartZ->at(i) - _util->GetOrigin()[2]);
1802  _MCPEndX.push_back(MCPEndX->at(i) - _util->GetOrigin()[0]);
1803  _MCPEndY.push_back(MCPEndY->at(i) - _util->GetOrigin()[1]);
1804  _MCPEndZ.push_back(MCPEndZ->at(i) - _util->GetOrigin()[2]);
1805  } else {
1806  _MCPStartX.push_back(MCPStartX->at(i));
1807  _MCPStartY.push_back(MCPStartY->at(i));
1808  _MCPStartZ.push_back(MCPStartZ->at(i));
1809  _MCPEndX.push_back(MCPEndX->at(i));
1810  _MCPEndY.push_back(MCPEndY->at(i));
1811  _MCPEndZ.push_back(MCPEndZ->at(i));
1812  }
1813  pdgmother.push_back(PDGMother->at(i));
1814  // save the true momentum
1815  truep.push_back(ptrue);
1816  // save the true angle
1817  _angle.push_back(angle);
1818  //Save MC process
1819  _MCProc.push_back(mcp_process);
1820  _MCEndProc.push_back(mcp_endprocess);
1821  mctime.push_back(time);
1822  //Case the endpoint is outside the CALO -> it went through the ECAL (mu/pi/p possible)
1823  //the ECAL will see 60 MIPs on average
1824  double Evis = (double)nLayers; //in MIP
1825  //Smearing to account for Gaussian detector noise (Landau negligible)
1826  Evis = _util->GaussianSmearing(Evis, ECAL_MIP_Res);
1827  //1 MIP = 0.814 MeV
1828  double Erec = Evis * MIP2GeV_factor;
1829  erecon.push_back((Erec > 0) ? Erec : 0.);
1830  etime.push_back(ecaltime);
1831  detected.push_back(1);
1832 
1833  //Muon/Pions/Protons are reco as Muons (without MuID detector)
1834  if( abs(pdg) == 13 || abs(pdg) == 211 || abs(pdg) == 2212 ) {
1835  recopidecal.push_back(13);
1836  }
1837  else {
1838  recopidecal.push_back(-1);
1839  }
1840 
1841  _preco.push_back(-1);
1842  anglereco.push_back(-1);
1843  recopid.push_back(-1);
1844  for (int pidr = 0; pidr < 6; ++pidr) prob_arr.push_back(-1);
1845 
1846  mctrkid.push_back(MCPTrkID->at(i));
1847  motherid.push_back(MCMotherTrkID->at(i));
1848  }
1849  else if(_util->PointInTPC(epoint) || _util->PointStopBetween(epoint))
1850  {
1851  truepdg.push_back(pdg);
1852  detected.push_back(0);
1853  truepx.push_back(MCPStartPX->at(i));
1854  truepy.push_back(MCPStartPY->at(i));
1855  truepz.push_back(MCPStartPZ->at(i));
1856  if(_correct4origin){
1857  _MCPStartX.push_back(MCPStartX->at(i) - _util->GetOrigin()[0]);
1858  _MCPStartY.push_back(MCPStartY->at(i) - _util->GetOrigin()[1]);
1859  _MCPStartZ.push_back(MCPStartZ->at(i) - _util->GetOrigin()[2]);
1860  _MCPEndX.push_back(MCPEndX->at(i) - _util->GetOrigin()[0]);
1861  _MCPEndY.push_back(MCPEndY->at(i) - _util->GetOrigin()[1]);
1862  _MCPEndZ.push_back(MCPEndZ->at(i) - _util->GetOrigin()[2]);
1863  } else {
1864  _MCPStartX.push_back(MCPStartX->at(i));
1865  _MCPStartY.push_back(MCPStartY->at(i));
1866  _MCPStartZ.push_back(MCPStartZ->at(i));
1867  _MCPEndX.push_back(MCPEndX->at(i));
1868  _MCPEndY.push_back(MCPEndY->at(i));
1869  _MCPEndZ.push_back(MCPEndZ->at(i));
1870  }
1871  pdgmother.push_back(PDGMother->at(i));
1872  // save the true momentum
1873  truep.push_back(ptrue);
1874  // save the true angle
1875  _angle.push_back(angle);
1876  //Save MC process
1877  _MCProc.push_back(mcp_process);
1878  _MCEndProc.push_back(mcp_endprocess);
1879  mctime.push_back(time);
1880  etime.push_back(-1);
1881  erecon.push_back(-1);
1882  _preco.push_back(-1);
1883  anglereco.push_back(-1);
1884  recopid.push_back(-1);
1885  recopidecal.push_back(-1);
1886  for (int pidr = 0; pidr < 6; ++pidr) prob_arr.push_back(-1);
1887  mctrkid.push_back(MCPTrkID->at(i));
1888  motherid.push_back(MCMotherTrkID->at(i));
1889  }
1890  else {
1891 
1892  }
1893  }// end is not neutron, pi0 or gamma
1894  }// end not visible in TPC
1895  nFSP++;
1896  } // closes the MC truth loop
1897 
1898  _nFSP.push_back(nFSP);
1899 
1900  //Check if vectors have good size
1901  if(this->CheckVectorSize()) {
1902  this->FillTTree();
1903  } else {
1904  std::cerr << "Event " << _Event << std::endl;
1905  std::cerr << "Number of FSP " << nFSP << std::endl;
1906 
1907  std::cerr << "Size of pdgmother " << pdgmother.size() << std::endl;
1908  std::cerr << "Size of truepdg " << truepdg.size() << std::endl;
1909  std::cerr << "Size of mctime " << mctime.size() << std::endl;
1910  std::cerr << "Size of mctrkid " << mctrkid.size() << std::endl;
1911  std::cerr << "Size of motherid " << motherid.size() << std::endl;
1912  std::cerr << "Size of _MCPStartX " << _MCPStartX.size() << std::endl;
1913  std::cerr << "Size of _MCPStartY " << _MCPStartY.size() << std::endl;
1914  std::cerr << "Size of _MCPStartZ " << _MCPStartZ.size() << std::endl;
1915  std::cerr << "Size of _MCPEndX " << _MCPEndX.size() << std::endl;
1916  std::cerr << "Size of _MCPEndY " << _MCPEndY.size() << std::endl;
1917  std::cerr << "Size of _MCPEndZ " << _MCPEndZ.size() << std::endl;
1918  std::cerr << "Size of _MCProc " << _MCProc.size() << std::endl;
1919  std::cerr << "Size of _MCEndProc " << _MCEndProc.size() << std::endl;
1920  std::cerr << "Size of trkLen " << trkLen.size() << std::endl;
1921  std::cerr << "Size of trkLenPerp " << trkLenPerp.size() << std::endl;
1922  std::cerr << "Size of truep " << truep.size() << std::endl;
1923  std::cerr << "Size of truepx " << truepx.size() << std::endl;
1924  std::cerr << "Size of truepy " << truepy.size() << std::endl;
1925  std::cerr << "Size of truepz " << truepz.size() << std::endl;
1926  std::cerr << "Size of _angle " << _angle.size() << std::endl;
1927 
1928  //Reco values
1929  std::cerr << "Size of detected " << detected.size() << std::endl;
1930  std::cerr << "Size of recopid " << recopid.size() << std::endl;
1931  std::cerr << "Size of recopidecal " << recopidecal.size() << std::endl;
1932  // std::cerr << "Size of prob_arr " << prob_arr.size() << std::endl;
1933  std::cerr << "Size of anglereco " << anglereco.size() << std::endl;
1934  std::cerr << "Size of _preco " << _preco.size() << std::endl;
1935  std::cerr << "Size of erecon " << erecon.size() << std::endl;
1936  std::cerr << "Size of etime " << etime.size() << std::endl;
1937 
1938  //Geometry
1939  std::cerr << "Size of isFidStart " << isFidStart.size() << std::endl;
1940  std::cerr << "Size of isTPCStart " << isTPCStart.size() << std::endl;
1941  std::cerr << "Size of isCaloStart " << isCaloStart.size() << std::endl;
1942  std::cerr << "Size of isThroughCaloStart " << isThroughCaloStart.size() << std::endl;
1943  std::cerr << "Size of isInBetweenStart " << isInBetweenStart.size() << std::endl;
1944  std::cerr << "Size of isBarrelStart " << isBarrelStart.size() << std::endl;
1945  std::cerr << "Size of isEndcapStart " << isEndcapStart.size() << std::endl;
1946 
1947  std::cerr << "Size of isFidEnd " << isFidEnd.size() << std::endl;
1948  std::cerr << "Size of isTPCEnd " << isTPCEnd.size() << std::endl;
1949  std::cerr << "Size of isCaloEnd " << isCaloEnd.size() << std::endl;
1950  std::cerr << "Size of isThroughCaloEnd " << isThroughCaloEnd.size() << std::endl;
1951  std::cerr << "Size of isInBetweenEnd " << isInBetweenEnd.size() << std::endl;
1952  std::cerr << "Size of isBarrelEnd " << isBarrelEnd.size() << std::endl;
1953  std::cerr << "Size of isEndcapEnd " << isEndcapEnd.size() << std::endl;
1954 
1955  std::cerr << "Event with wrong vector sizes... skipped" << std::endl;
1956  }
1957  } // closes the event loop
1958 } // closes the main loop function
Utils::isBarrel
bool isBarrel(const TVector3 &point)
Definition: Utils.cpp:96
CAF::_angle
std::vector< double > _angle
Definition: CAF.h:78
CAF::tgtpdg
std::vector< int > tgtpdg
Definition: CAF.h:69
CAF::GPartFirstDaugh
std::vector< int > GPartFirstDaugh
Definition: CAF.h:70
CAF::w
std::vector< double > w
Definition: CAF.h:73
CAF::isInBetweenStart
std::vector< unsigned int > isInBetweenStart
Definition: CAF.h:83
entry
QList< Entry > entry
Definition: tclscanner.cpp:1150
CAF::GPartStatus
std::vector< int > GPartStatus
Definition: CAF.h:70
CAF::recopid
std::vector< int > recopid
Definition: CAF.h:80
CAF::GPartPz
std::vector< float > GPartPz
Definition: CAF.h:71
result
static QCString result
Definition: fortranscanner.cpp:56614
CAF::x
std::vector< double > x
Definition: CAF.h:73
CAF::_MCEndProc
std::vector< std::string > _MCEndProc
Definition: CAF.h:77
CAF::isBarrelStart
std::vector< unsigned int > isBarrelStart
Definition: CAF.h:85
CAF::pdg
int pdg[100]
Definition: CAF.h:50
CAF::ntype
std::vector< int > ntype
Definition: CAF.h:69
string
std::string string
Definition: nybbler.cc:12
CAF::truep
std::vector< double > truep
Definition: CAF.h:78
CAF::etime
std::vector< double > etime
Definition: CAF.h:81
CAF::_SubRun
unsigned int _SubRun
Definition: CAF.h:67
CAF::GPartPy
std::vector< float > GPartPy
Definition: CAF.h:71
CAF::trkLen
std::vector< double > trkLen
Definition: CAF.h:78
CAF::trkLenPerp
std::vector< double > trkLenPerp
Definition: CAF.h:78
CAF::_MCPStartY
std::vector< int > _MCPStartY
Definition: CAF.h:76
CAF::nFSP
int nFSP
Definition: CAF.h:49
CAF::_nFSP
std::vector< unsigned int > _nFSP
Definition: CAF.h:75
CAF::GPartName
std::vector< std::string > GPartName
Definition: CAF.h:72
Utils::GaussianSmearing
float GaussianSmearing(const float &mean, const float &sigma)
Definition: Utils.h:58
CAF::_preco
std::vector< double > _preco
Definition: CAF.h:81
CAF::calcGluck
float calcGluck(double sigmaX, double B, double X0, float nHits, double mom, double length, double &ratio)
Definition: CAF.cpp:337
CAF::isEndcapEnd
std::vector< unsigned int > isEndcapEnd
Definition: CAF.h:85
CAF::motherid
std::vector< int > motherid
Definition: CAF.h:76
CAF::isFidStart
std::vector< unsigned int > isFidStart
Definition: CAF.h:83
CAF::GPartMass
std::vector< float > GPartMass
Definition: CAF.h:71
train.filename
string filename
Definition: train.py:213
CAF::mctime
std::vector< double > mctime
Definition: CAF.h:73
CAF::isBarrelEnd
std::vector< unsigned int > isBarrelEnd
Definition: CAF.h:85
CAF::gint
std::vector< int > gint
Definition: CAF.h:69
CAF::mcnupy
std::vector< double > mcnupy
Definition: CAF.h:73
CAF::recopidecal
std::vector< int > recopidecal
Definition: CAF.h:80
CAF::isThroughCaloStart
std::vector< unsigned int > isThroughCaloStart
Definition: CAF.h:83
CAF::erecon
std::vector< double > erecon
Definition: CAF.h:81
PDG
const uint PDG
Definition: qregexp.cpp:140
CAF::GPartFirstMom
std::vector< int > GPartFirstMom
Definition: CAF.h:70
CAF::mcnupz
std::vector< double > mcnupz
Definition: CAF.h:73
CAF::theta
std::vector< double > theta
Definition: CAF.h:73
Utils::GetOrigin
double * GetOrigin()
Definition: Utils.h:60
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
CAF::_MCPEndZ
std::vector< int > _MCPEndZ
Definition: CAF.h:76
CAF::_MCPStartX
std::vector< int > _MCPStartX
Definition: CAF.h:76
infile
string infile
Definition: INukeNucleonCorr.cxx:62
CAF::detected
std::vector< int > detected
Definition: CAF.h:69
CAF::P
std::pair< float, std::string > P
Definition: CAF.h:14
Utils::PointInTPC
bool PointInTPC(const TVector3 &point)
Definition: Utils.cpp:43
CAF::isThroughCaloEnd
std::vector< unsigned int > isThroughCaloEnd
Definition: CAF.h:84
submit_ppfxjobs.temp
temp
Definition: submit_ppfxjobs.py:110
Utils::isThroughCalo
bool isThroughCalo(const TVector3 &point)
Definition: Utils.cpp:86
CAF::_Event
unsigned int _Event
Definition: CAF.h:67
CAF::intert
std::vector< int > intert
Definition: CAF.h:69
CAF::anglereco
std::vector< double > anglereco
Definition: CAF.h:81
CAF::_MCPStartZ
std::vector< int > _MCPStartZ
Definition: CAF.h:76
CAF::truepdg
std::vector< int > truepdg
Definition: CAF.h:76
test.p
p
Definition: test.py:223
Utils::PointInCalo
bool PointInCalo(const TVector3 &point)
Definition: Utils.cpp:58
CAF::verty
std::vector< double > verty
Definition: CAF.h:73
CAF::m_pidinterp
std::unordered_map< int, TH2F * > m_pidinterp
Definition: CAF.h:55
protoana::distance
double distance(double x1, double y1, double z1, double x2, double y2, double z2)
Definition: truepionXsection_module.cc:217
CAF::GPartE
std::vector< float > GPartE
Definition: CAF.h:71
CAF::ccnc
std::vector< int > ccnc
Definition: CAF.h:69
CAF::_Run
unsigned int _Run
Definition: CAF.h:67
CAF::_util
Utils * _util
Definition: CAF.h:60
CAF::q2
std::vector< double > q2
Definition: CAF.h:73
CAF::isInBetweenEnd
std::vector< unsigned int > isInBetweenEnd
Definition: CAF.h:84
CAF::nGPart
std::vector< int > nGPart
Definition: CAF.h:70
CAF::pdgmother
std::vector< int > pdgmother
Definition: CAF.h:76
CAF::_MCPEndY
std::vector< int > _MCPEndY
Definition: CAF.h:76
CAF::ptrue
double ptrue[100]
Definition: CAF.h:51
CAF::isTPCStart
std::vector< unsigned int > isTPCStart
Definition: CAF.h:83
CAF::vertz
std::vector< double > vertz
Definition: CAF.h:73
CAF::truepz
std::vector< double > truepz
Definition: CAF.h:78
CAF::vertx
std::vector< double > vertx
Definition: CAF.h:73
larg4::dist
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
Definition: OpFastScintillation.hh:563
CAF::truepx
std::vector< double > truepx
Definition: CAF.h:78
CAF::weight
std::vector< int > weight
Definition: CAF.h:69
Event
Definition: types.h:32
CAF::GPartLastMom
std::vector< int > GPartLastMom
Definition: CAF.h:70
CAF::_correct4origin
unsigned int _correct4origin
sets the string for the coordinates origins (World or TPC)
Definition: CAF.h:64
CAF::GPartPx
std::vector< float > GPartPx
Definition: CAF.h:71
CAF::_inttree
TTree * _inttree
Definition: CAF.h:58
CAF::CheckVectorSize
bool CheckVectorSize()
Definition: CAF.cpp:282
CAF::isCaloEnd
std::vector< unsigned int > isCaloEnd
Definition: CAF.h:84
CAF::isFidEnd
std::vector< unsigned int > isFidEnd
Definition: CAF.h:84
CAF::ClearVectors
void ClearVectors()
Definition: CAF.cpp:193
CAF::GPartPdg
std::vector< int > GPartPdg
Definition: CAF.h:70
Utils::PointInFiducial
bool PointInFiducial(const TVector3 &point)
Definition: Utils.cpp:29
CAF::gt_t
std::vector< int > gt_t
Definition: CAF.h:69
CAF::mode
std::vector< int > mode
Definition: CAF.h:69
CAF::FillTTree
void FillTTree()
Definition: CAF.cpp:183
CAF::isTPCEnd
std::vector< unsigned int > isTPCEnd
Definition: CAF.h:84
CAF::y
std::vector< double > y
Definition: CAF.h:73
CAF::isCaloStart
std::vector< unsigned int > isCaloStart
Definition: CAF.h:83
CAF::truepy
std::vector< double > truepy
Definition: CAF.h:78
CAF::mcnupx
std::vector< double > mcnupx
Definition: CAF.h:73
Utils::isEndcap
bool isEndcap(const TVector3 &point)
Definition: Utils.cpp:107
CAF::isEndcapStart
std::vector< unsigned int > isEndcapStart
Definition: CAF.h:85
CAF::mctrkid
std::vector< int > mctrkid
Definition: CAF.h:76
CAF::GPartLastDaugh
std::vector< int > GPartLastDaugh
Definition: CAF.h:70
s
static QCString * s
Definition: config.cpp:1042
Utils::PointStopBetween
bool PointStopBetween(const TVector3 &point)
Definition: Utils.cpp:72
str
static QCString str
Definition: fortrancode.cpp:27098
CAF::prob_arr
std::vector< double > prob_arr
Definition: CAF.h:81
CAF::_MCPEndX
std::vector< int > _MCPEndX
Definition: CAF.h:76
Utils::GetRamdomNumber
float GetRamdomNumber()
Definition: Utils.h:56
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
QCString::append
QCString & append(const char *s)
Definition: qcstring.cpp:383
CAF::_MCProc
std::vector< std::string > _MCProc
Definition: CAF.h:77
void CAF::Print ( )
void CAF::Print ( )
void CAF::setToBS ( )
void CAF::setToBS ( )
void CAF::write ( )
void CAF::write ( )
void CAF::WriteTTree ( )

Definition at line 172 of file CAF.cpp.

173 {
174  if(cafMVA != nullptr) {
175  cafFile->cd();
176  cafMVA->Write();
177  cafFile->Close();
178  }
179 
180  return;
181 }
CAF::cafFile
TFile * cafFile
The output TFile pointer */.
Definition: CAF.h:53
CAF::cafMVA
TTree * cafMVA
The output TTree pointer */.
Definition: CAF.h:54

Member Data Documentation

std::vector<double> CAF::_angle
private

Definition at line 78 of file CAF.h.

unsigned int CAF::_correct4origin
private

sets the string for the coordinates origins (World or TPC)

Definition at line 64 of file CAF.h.

unsigned int CAF::_Event
private

Definition at line 67 of file CAF.h.

std::string CAF::_inputfile
private

Definition at line 62 of file CAF.h.

TFile* CAF::_intfile
private

Definition at line 57 of file CAF.h.

TTree* CAF::_inttree
private

Definition at line 58 of file CAF.h.

std::vector<std::string> CAF::_MCEndProc
private

Definition at line 77 of file CAF.h.

std::vector<int> CAF::_MCPEndX
private

Definition at line 76 of file CAF.h.

std::vector<int> CAF::_MCPEndY
private

Definition at line 76 of file CAF.h.

std::vector<int> CAF::_MCPEndZ
private

Definition at line 76 of file CAF.h.

std::vector<std::string> CAF::_MCProc
private

Definition at line 77 of file CAF.h.

std::vector<int> CAF::_MCPStartX
private

Definition at line 76 of file CAF.h.

std::vector<int> CAF::_MCPStartY
private

Definition at line 76 of file CAF.h.

std::vector<int> CAF::_MCPStartZ
private

Definition at line 76 of file CAF.h.

std::vector<unsigned int> CAF::_nFSP
private

Definition at line 75 of file CAF.h.

std::string CAF::_outputFile
private

The output TFile name */.

Definition at line 63 of file CAF.h.

std::vector<double> CAF::_preco
private

Definition at line 81 of file CAF.h.

unsigned int CAF::_Run
private

Definition at line 67 of file CAF.h.

unsigned int CAF::_SubRun
private

Definition at line 67 of file CAF.h.

Utils* CAF::_util
private

Definition at line 60 of file CAF.h.

std::vector<double> CAF::anglereco
private

Definition at line 81 of file CAF.h.

TFile * CAF::cafFile
private

The output TFile pointer */.

Definition at line 53 of file CAF.h.

TTree * CAF::cafMVA
private

The output TTree pointer */.

Definition at line 54 of file CAF.h.

TTree * CAF::cafPOT

Definition at line 47 of file NUSYS.h.

std::vector<int> CAF::ccnc
private

Definition at line 69 of file CAF.h.

double CAF::cvwgt

Definition at line 33 of file NUSYS.h.

double CAF::det_x

Definition at line 36 of file CAF.h.

std::vector<int> CAF::detected
private

Definition at line 69 of file CAF.h.

double CAF::Ehad_veto

Definition at line 44 of file CAF.h.

double CAF::Elep_reco

Definition at line 39 of file CAF.h.

double CAF::eN

Definition at line 31 of file CAF.h.

double CAF::eOther

Definition at line 31 of file CAF.h.

double CAF::eP

Definition at line 31 of file CAF.h.

double CAF::ePi0

Definition at line 31 of file CAF.h.

double CAF::ePim

Definition at line 31 of file CAF.h.

double CAF::ePip

Definition at line 31 of file CAF.h.

double CAF::eRecoN

Definition at line 32 of file CAF.h.

std::vector<double> CAF::erecon
private

Definition at line 81 of file CAF.h.

double CAF::eRecoOther

Definition at line 32 of file CAF.h.

double CAF::eRecoP

Definition at line 32 of file CAF.h.

double CAF::eRecoPi0

Definition at line 32 of file CAF.h.

double CAF::eRecoPim

Definition at line 32 of file CAF.h.

double CAF::eRecoPip

Definition at line 32 of file CAF.h.

std::vector<double> CAF::etime
private

Definition at line 81 of file CAF.h.

double CAF::Ev

Definition at line 28 of file CAF.h.

double CAF::Ev_reco

Definition at line 39 of file CAF.h.

int CAF::event

Definition at line 25 of file NUSYS.h.

int CAF::gastpc_pi_min_mult

Definition at line 48 of file CAF.h.

int CAF::gastpc_pi_pl_mult

Definition at line 48 of file CAF.h.

TTree * CAF::genie

Definition at line 48 of file NUSYS.h.

std::vector< std::vector < std::vector < uint64_t > > >* CAF::geoEffThrowResults

Definition at line 64 of file CAF.h.

std::vector<int> CAF::gint
private

Definition at line 69 of file CAF.h.

std::vector<float> CAF::GPartE
private

Definition at line 71 of file CAF.h.

std::vector<int> CAF::GPartFirstDaugh
private

Definition at line 70 of file CAF.h.

std::vector<int> CAF::GPartFirstMom
private

Definition at line 70 of file CAF.h.

std::vector<int> CAF::GPartLastDaugh
private

Definition at line 70 of file CAF.h.

std::vector<int> CAF::GPartLastMom
private

Definition at line 70 of file CAF.h.

std::vector<float> CAF::GPartMass
private

Definition at line 71 of file CAF.h.

std::vector<std::string> CAF::GPartName
private

Definition at line 72 of file CAF.h.

std::vector<int> CAF::GPartPdg
private

Definition at line 70 of file CAF.h.

std::vector<float> CAF::GPartPx
private

Definition at line 71 of file CAF.h.

std::vector<float> CAF::GPartPy
private

Definition at line 71 of file CAF.h.

std::vector<float> CAF::GPartPz
private

Definition at line 71 of file CAF.h.

std::vector<int> CAF::GPartStatus
private

Definition at line 70 of file CAF.h.

std::vector<int> CAF::gt_t
private

Definition at line 69 of file CAF.h.

int CAF::ievt

Definition at line 25 of file NUSYS.h.

std::vector<int> CAF::intert
private

Definition at line 69 of file CAF.h.

std::vector<unsigned int> CAF::isBarrelEnd
private

Definition at line 85 of file CAF.h.

std::vector<unsigned int> CAF::isBarrelStart
private

Definition at line 85 of file CAF.h.

std::vector<unsigned int> CAF::isCaloEnd
private

Definition at line 84 of file CAF.h.

std::vector<unsigned int> CAF::isCaloStart
private

Definition at line 83 of file CAF.h.

int CAF::isCC

Definition at line 27 of file CAF.h.

std::vector<unsigned int> CAF::isEndcapEnd
private

Definition at line 85 of file CAF.h.

std::vector<unsigned int> CAF::isEndcapStart
private

Definition at line 85 of file CAF.h.

int CAF::isFD

Definition at line 23 of file NUSYS.h.

int CAF::isFHC

Definition at line 23 of file NUSYS.h.

std::vector<unsigned int> CAF::isFidEnd
private

Definition at line 84 of file CAF.h.

std::vector<unsigned int> CAF::isFidStart
private

Definition at line 83 of file CAF.h.

std::vector<unsigned int> CAF::isInBetweenEnd
private

Definition at line 84 of file CAF.h.

std::vector<unsigned int> CAF::isInBetweenStart
private

Definition at line 83 of file CAF.h.

std::vector<unsigned int> CAF::isThroughCaloEnd
private

Definition at line 84 of file CAF.h.

std::vector<unsigned int> CAF::isThroughCaloStart
private

Definition at line 83 of file CAF.h.

std::vector<unsigned int> CAF::isTPCEnd
private

Definition at line 84 of file CAF.h.

std::vector<unsigned int> CAF::isTPCStart
private

Definition at line 83 of file CAF.h.

bool CAF::iswgt

Definition at line 35 of file NUSYS.h.

double CAF::LepE

Definition at line 28 of file CAF.h.

double CAF::LepMomX

Definition at line 28 of file CAF.h.

double CAF::LepMomY

Definition at line 28 of file CAF.h.

double CAF::LepMomZ

Definition at line 28 of file CAF.h.

double CAF::LepNuAngle

Definition at line 28 of file CAF.h.

int CAF::LepPDG

Definition at line 27 of file CAF.h.

std::unordered_map<int, TH2F*> CAF::m_pidinterp
private

Definition at line 55 of file CAF.h.

std::vector<double> CAF::mcnupx
private

Definition at line 73 of file CAF.h.

std::vector<double> CAF::mcnupy
private

Definition at line 73 of file CAF.h.

std::vector<double> CAF::mcnupz
private

Definition at line 73 of file CAF.h.

genie::NtpMCEventRecord * CAF::mcrec

Definition at line 38 of file NUSYS.h.

std::vector<double> CAF::mctime
private

Definition at line 73 of file CAF.h.

std::vector<int> CAF::mctrkid
private

Definition at line 76 of file CAF.h.

int CAF::meta_run

Definition at line 42 of file NUSYS.h.

int CAF::meta_subrun

Definition at line 42 of file NUSYS.h.

int CAF::mode

Definition at line 27 of file CAF.h.

std::vector<int> CAF::mode
private

Definition at line 69 of file CAF.h.

std::vector<int> CAF::motherid
private

Definition at line 76 of file CAF.h.

int CAF::muon_contained

Definition at line 41 of file CAF.h.

int CAF::muon_ecal

Definition at line 41 of file CAF.h.

float CAF::muon_endpoint[3]

Definition at line 42 of file CAF.h.

std::string* CAF::muon_endVolName

Definition at line 43 of file CAF.h.

int CAF::muon_exit

Definition at line 41 of file CAF.h.

int CAF::muon_tracker

Definition at line 41 of file CAF.h.

int CAF::neutrinoPDG

Definition at line 27 of file CAF.h.

int CAF::neutrinoPDGunosc

Definition at line 27 of file CAF.h.

int CAF::nFSP

Definition at line 49 of file CAF.h.

std::vector<int> CAF::nGPart
private

Definition at line 70 of file CAF.h.

int CAF::niem

Definition at line 30 of file CAF.h.

int CAF::nik0

Definition at line 30 of file CAF.h.

int CAF::nikm

Definition at line 30 of file CAF.h.

int CAF::nikp

Definition at line 30 of file CAF.h.

int CAF::niother

Definition at line 30 of file CAF.h.

int CAF::nipi0

Definition at line 30 of file CAF.h.

int CAF::nipim

Definition at line 30 of file CAF.h.

int CAF::nipip

Definition at line 30 of file CAF.h.

int CAF::nN

Definition at line 30 of file CAF.h.

int CAF::nNucleus

Definition at line 30 of file CAF.h.

int CAF::nP

Definition at line 30 of file CAF.h.

std::vector<int> CAF::ntype
private

Definition at line 69 of file CAF.h.

double CAF::NuMomX

Definition at line 28 of file CAF.h.

double CAF::NuMomY

Definition at line 28 of file CAF.h.

double CAF::NuMomZ

Definition at line 28 of file CAF.h.

int CAF::nUNKNOWN

Definition at line 30 of file CAF.h.

int CAF::nwgt

Definition at line 32 of file NUSYS.h.

double CAF::partEvReco[100]

Definition at line 51 of file CAF.h.

int CAF::pdg[100]

Definition at line 50 of file CAF.h.

std::vector<int> CAF::pdgmother
private

Definition at line 76 of file CAF.h.

double CAF::pileup_energy

Definition at line 45 of file CAF.h.

double CAF::pot

Definition at line 41 of file NUSYS.h.

std::vector<double> CAF::prob_arr
private

Definition at line 81 of file CAF.h.

double CAF::ptrue[100]

Definition at line 51 of file CAF.h.

double CAF::Q2

Definition at line 28 of file CAF.h.

std::vector<double> CAF::q2
private

Definition at line 73 of file CAF.h.

int CAF::reco_lepton_pdg

Definition at line 41 of file CAF.h.

int CAF::reco_nc

Definition at line 40 of file CAF.h.

int CAF::reco_nue

Definition at line 40 of file CAF.h.

int CAF::reco_numu

Definition at line 40 of file CAF.h.

int CAF::reco_q

Definition at line 40 of file CAF.h.

std::vector<int> CAF::recopid
private

Definition at line 80 of file CAF.h.

std::vector<int> CAF::recopidecal
private

Definition at line 80 of file CAF.h.

int CAF::run

Definition at line 25 of file NUSYS.h.

int CAF::subrun

Definition at line 25 of file NUSYS.h.

std::vector<double> CAF::t
private

Definition at line 73 of file CAF.h.

std::vector<int> CAF::tgtpdg
private

Definition at line 69 of file CAF.h.

std::vector<double> CAF::theta
private

Definition at line 73 of file CAF.h.

double CAF::theta_reco

Definition at line 39 of file CAF.h.

double CAF::trkLen[100]

Definition at line 51 of file CAF.h.

std::vector<double> CAF::trkLen
private

Definition at line 78 of file CAF.h.

double CAF::trkLenPerp[100]

Definition at line 51 of file CAF.h.

std::vector<double> CAF::trkLenPerp
private

Definition at line 78 of file CAF.h.

std::vector<double> CAF::truep
private

Definition at line 78 of file CAF.h.

std::vector<int> CAF::truepdg
private

Definition at line 76 of file CAF.h.

std::vector<double> CAF::truepx
private

Definition at line 78 of file CAF.h.

std::vector<double> CAF::truepy
private

Definition at line 78 of file CAF.h.

std::vector<double> CAF::truepz
private

Definition at line 78 of file CAF.h.

int CAF::version

Definition at line 43 of file NUSYS.h.

std::vector<double> CAF::vertx
private

Definition at line 73 of file CAF.h.

std::vector<double> CAF::verty
private

Definition at line 73 of file CAF.h.

std::vector<double> CAF::vertz
private

Definition at line 73 of file CAF.h.

double CAF::vtx_x

Definition at line 35 of file CAF.h.

double CAF::vtx_y

Definition at line 35 of file CAF.h.

double CAF::vtx_z

Definition at line 35 of file CAF.h.

double CAF::W

Definition at line 28 of file CAF.h.

std::vector<double> CAF::w
private

Definition at line 73 of file CAF.h.

std::vector<int> CAF::weight
private

Definition at line 69 of file CAF.h.

double CAF::wgt

Definition at line 34 of file NUSYS.h.

double CAF::X

Definition at line 28 of file CAF.h.

std::vector<double> CAF::x
private

Definition at line 73 of file CAF.h.

double CAF::Y

Definition at line 28 of file CAF.h.

std::vector<double> CAF::y
private

Definition at line 73 of file CAF.h.

The documentation for this class was generated from the following files:
Generated by   doxygen 1.8.11