CAF.cpp
Go to the documentation of this file.
1 #include "CAF.h"
2 
3 #include "TRandom3.h"
4 #include "TVector2.h"
5 #include "TVector3.h"
6 #include "TLorentzVector.h"
7 #include "TH1F.h"
8 #include "TRandom.h"
9 #include "TF1.h"
10 #include "TDatabasePDG.h"
11 #include "TParticlePDG.h"
12 
13 #include <stdio.h>
14 #include <vector>
15 #include <iostream>
16 #include <string>
17 #include <stdlib.h>
18 #include <numeric>
19 #include <functional>
20 #include <exception>
21 
22 CAF::CAF()
23 : cafFile(nullptr), _intfile(nullptr), _inttree(nullptr), _util(new Utils()), _inputfile(""), _outputFile(""), _correct4origin(0)
24 {
25 
26 }
27 
28 CAF::CAF( std::string infile, std::string filename, int correct4origin, double *originTPC)
29 : cafFile(nullptr), _intfile(nullptr), _inttree(nullptr), _util(new Utils()), _inputfile(infile), _outputFile(filename), _correct4origin(correct4origin)
30 {
31  _util->SetOrigin(originTPC);
32 }
33 
34 CAF::~CAF()
35 {
36  delete _util;
37  delete cafMVA;
38  delete cafFile;
39  delete _inttree;
40  delete _intfile;
41 }
42 
43 bool CAF::BookTFile()
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 }
162 
163 void CAF::CloseTFile()
164 {
165  if(cafFile != nullptr) {
166  cafFile->Close();
167  return;
168  }
169  else{ return; }
170 }
171 
172 void CAF::WriteTTree()
173 {
174  if(cafMVA != nullptr) {
175  cafFile->cd();
176  cafMVA->Write();
177  cafFile->Close();
178  }
179 
180  return;
181 }
182 
183 void CAF::FillTTree()
184 {
185  if(cafMVA != nullptr) {
186  cafFile->cd();
187  cafMVA->Fill();
188  }
189 
190  return;
191 }
192 
193 void CAF::ClearVectors()
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 }
281 
282 bool CAF::CheckVectorSize()
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 }
336 
337 float CAF::calcGluck(double sigmaX, double B, double X0, float nHits, double mom, double length, double& ratio)
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 }
346 
347 // main loop function
348 void CAF::loop()
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::_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
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::CAF
CAF()
Definition: CAF.cpp:22
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::BookTFile
bool BookTFile()
Definition: CAF.cpp:43
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::WriteTTree
void WriteTTree()
Definition: CAF.cpp:172
CAF::mcnupy
std::vector< double > mcnupy
Definition: CAF.h:73
CAF::t
std::vector< double > t
Definition: CAF.h:73
CAF::CloseTFile
void CloseTFile()
Definition: CAF.cpp:163
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::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
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::loop
void loop()
Definition: CAF.cpp:348
CAF::truepdg
std::vector< int > truepdg
Definition: CAF.h:76
Utils::SetOrigin
void SetOrigin(double *origin)
Definition: Utils.cpp:22
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::~CAF
~CAF()
Definition: CAF.cpp:34
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::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
Utils
Definition: Utils.h:7
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
Generated by   doxygen 1.8.11