Macros | Enumerations | Functions | Variables
makeDST.cxx File Reference
#include <list>
#include "globbing.h"
#include "treeReader.h"
#include "includeROOT.h"
#include "nusimdata/SimulationBase/MCNeutrino.h"
#include "ReconstructionDataProducts/Track.h"
#include "ReconstructionDataProducts/IDNumberGen.h"

Go to the source code of this file.

Macros

#define writeECAL   true
 
#define killNoRecoVerts   true
 
#define OUTDEV   TextFileOut
 

Enumerations

enum  CounterTag {
  nEvent, nDatainFid, nCCQEinFid, nNCQEinFid,
  nRESinFid, nCCDISinFid, nNCDISinFid, nCCCOHinFid,
  nNCCOHinFid, nNuEinFid, nCounterTag
}
 

Functions

void TransMogrifyTree (int fileNumber, bool firstCall)
 
bool inFiducial (double x, double y, double z)
 
int main (int argc, const char *argv[])
 

Variables

bool viaXrootd = true
 
float PmagCut = 0.150
 
TFile * outFile
 
TFile * inFile
 
TTree * outTree
 
TTree * inTree
 
int Counter [nCounterTag]
 

Macro Definition Documentation

#define killNoRecoVerts   true

Definition at line 30 of file makeDST.cxx.

#define OUTDEV   TextFileOut
#define writeECAL   true

Definition at line 29 of file makeDST.cxx.

Enumeration Type Documentation

enum CounterTag
Enumerator
nEvent 
nDatainFid 
nCCQEinFid 
nNCQEinFid 
nRESinFid 
nCCDISinFid 
nNCDISinFid 
nCCCOHinFid 
nNCCOHinFid 
nNuEinFid 
nCounterTag 

Definition at line 41 of file makeDST.cxx.

41  {nEvent, nDatainFid,
42  nCCQEinFid, nNCQEinFid,
43  nRESinFid, nCCDISinFid,
44  nNCDISinFid, nCCCOHinFid,
45  nNCCOHinFid, nNuEinFid,
46 nCounterTag}; // Must be last enum

Function Documentation

bool inFiducial ( double  x,
double  y,
double  z 
)

Definition at line 17 of file makeDST.cxx.

17  {
18  double const Rcut = 445.0 /2.0; double const Xcut = 430.0 /2.0;
19  double r = sqrt( z*z + y*y );
20  if ( r > Rcut ) return false;
21  if (fabs(x) > Xcut ) return false;
22  return true;
23 }
train.x
list x
Definition: train.py:276
plot_model.r
r
Definition: plot_model.py:57
int main ( int  argc,
const char *  argv[] 
)

Definition at line 53 of file makeDST.cxx.

53  {
54  int maxNfiles = -1; string outName, SearchDir; bool gotOutDir = false;
55  int iArg = 1;
56  while ((iArg < argc) && (argv[iArg][0] == '-')){
57  switch (argv[iArg][1]) {
58 
59  case 'f':
60  maxNfiles = atoi(argv[iArg+1]);
61  cout << "Only processing " << maxNfiles << " files." << endl;
62  if (maxNfiles < 1) {
63  cout << "Wait! " << maxNfiles << "files! That's not right!" << endl;
64  exit(1);
65  }
66  ++iArg;
67  break;
68 
69  case 'c':
70  PmagCut = std::stof(argv[iArg+1]);
71  cout << "Removing tracks and vertices with tracks with P below " <<
72  PmagCut << " MeV." << endl;
73  PmagCut /= 1000.0;
74  ++iArg;
75  break;
76 
77  case 'L':
78  SearchDir = argv[iArg+1];
79  viaXrootd = false;
80  ++iArg;
81  break;
82 
83  case 'o':
84  // Cleverly, argv comes to us with environmental variables expanded
85  // and we don't need wordexp.h.
86  outName = argv[iArg+1];
87  gotOutDir = true;
88  ++iArg;
89  break;
90 
91  case 'h':
92  default:
93  cout << "[-h] This message" << endl;
94  cout << "[-o name] Name for output files; TTree will be in <name>.root," << endl;
95  cout << " statistics will be in <name>.txt. Environmental" << endl;
96  cout << " variables in <name> will be expanded." << endl;
97  cout << "[-f n] Analyze n background files, max" << endl;
98  cout << "[-c m] Set Ptrk cut to Ptrk > m MeV for output." << endl;
99  cout << " The default cut is 150 MeV." << endl;
100  cout << "[-L name] Input files from directory <name> without using" << endl;
101  cout << " xrootd rather than from files in $PWD, with" << endl;
102  cout << " xrootd, which is the default." << endl;
103  exit(0);
104  }
105  ++iArg;
106  }
107  if (!gotOutDir) {
108  cout << "Did not see -o <dir> so no output directory selected." << endl;
109  exit(1);
110  }
111 
112 
113 
114  // Get list of input files
115  string LookFor;
116  vector<string> fileList;
117  int nFiles;
118 
119  if (viaXrootd){
120  SearchDir = "";
121  LookFor = "anatree_*.root";
122  fileList = globbing(SearchDir +LookFor);
123  nFiles = (int)fileList.size();
124  string newAccess = "root://fndca1.fnal.gov:1094/pnfs/fnal.gov/usr";
125  string PWD = getenv("PWD"); PWD.erase(0,5); PWD += "/";
126  newAccess += PWD;
127 
128  for (int iFile=0; iFile<nFiles; ++iFile) {
129  fileList[iFile].insert(0,newAccess);
130  }
131  } else {
132  LookFor = "anatree_*.root";
133  fileList = globbing(SearchDir +LookFor);
134  nFiles = (int)fileList.size();
135  }
136  cout << ( nFiles = (int)fileList.size() ) << " total files found.\n";
137  if (nFiles==0) exit(1);
138 
139  if ( maxNfiles>0 && maxNfiles<nFiles) {
140  nFiles = maxNfiles;
141  }
142  cout << endl;
143 
144 
145 
146  // Output file made and opened here
147  string RootOutName = outName + ".root";
148  if (!viaXrootd && !filehamna(RootOutName)) {
149  cout << "Already existing " << RootOutName << endl;
150  exit(2);
151  }
152  outFile = TFile::Open(RootOutName.c_str(),"NEW");
153  if (outFile->IsZombie()) {cout << "Eat you brains!" << endl; return 1;};
154  outFile->mkdir("anatree");
155  outFile->cd("anatree");
156  outTree = new TTree("GArAnaTree","GArAnaTree");
157 
158 
159 
160  time_t unixT; struct tm* localT;
161  time(&unixT); localT = localtime(&unixT);
162  cout << "Starting at " << asctime(localT) << endl;
163 
164  for (int iFile=0; iFile<nFiles; ++iFile) {
165  cout << "Reading file " << iFile+1 << ": " << fileList[iFile].c_str() << " ";
166  inFile = TFile::Open(fileList[iFile].c_str(),"READ");
167  if (!inFile || !inFile->IsOpen()) {
168  cout << "...but that file can not be opened!" << endl;
169  exit(3);
170  }
171 
172  size_t strtLoc = fileList[iFile].find("_") +1;
173  size_t lenLoc = fileList[iFile].find(".root") -strtLoc;
174  int fileNumber;
175  string numberAsString = fileList[iFile].substr(strtLoc,lenLoc);
176  sscanf(numberAsString.c_str(), "%d", &fileNumber);
177 
178  TDirectory* dir = (TDirectory*)(inFile)->Get("anatree");
179  dir->GetObject("GArAnaTree",inTree);
180  TransMogrifyTree(fileNumber, iFile==0);
181  inFile->Close();
182  }
183  outFile->Write();
184  outFile->Close();
185 
186  time(&unixT); localT = localtime(&unixT);
187  cout << "Stopping at " << asctime(localT) << endl;
188 
189 
190 
191  string TextOutName = outName + ".txt";
192  #define OUTDEV TextFileOut
193  std::ofstream TextFileOut(TextOutName.c_str());
194  // #define OUTDEV cout
195  cout << endl << endl;
196 
197  OUTDEV << Counter[nEvent] << "\tevents" << endl;
198 
199  OUTDEV << Counter[nDatainFid] << "\tevents with MC vertex in fiducial"
200  << endl;
201 
202  OUTDEV << Counter[nCCQEinFid] << "\tCCQE events with MC vertex in fiducial"
203  << endl;
204 
205  OUTDEV << Counter[nNCQEinFid] << "\tNCQE events with MC vertex in fiducial"
206  << endl;
207 
208  OUTDEV << Counter[nRESinFid] << "\tRES events with MC vertex in fiducial"
209  << endl;
210 
211  OUTDEV << Counter[nCCDISinFid] << "\tCC DIS events with MC vertex in fiducial"
212  << endl;
213 
214  OUTDEV << Counter[nNCDISinFid] << "\tNC DIS events with MC vertex in fiducial"
215  << endl;
216 
217  OUTDEV << Counter[nCCCOHinFid] << "\tCC coherent pi events with MC vertex in fiducial"
218  << endl;
219 
220  OUTDEV << Counter[nNCCOHinFid] << "\tNC coherent pi events with MC vertex in fiducial"
221  << endl;
222 
223  OUTDEV << Counter[nNuEinFid] << "\tElastic nu-e events with MC vertex in fiducial"
224  << endl;
225 
226  TextFileOut.close();
227 
228  exit(0);
229 }
PmagCut
float PmagCut
Definition: makeDST.cxx:32
outTree
TTree * outTree
Definition: makeDST.cxx:37
Counter
int Counter[nCounterTag]
Definition: makeDST.cxx:47
viaXrootd
bool viaXrootd
Definition: makeDST.cxx:31
outFile
TFile * outFile
Definition: makeDST.cxx:36
dir
string dir
Definition: INukeNucleonCorr.cxx:61
demo.tm
tm
Definition: demo.py:21
filehamna
bool filehamna(const std::string &filename)
Definition: includeROOT.h:136
TransMogrifyTree
void TransMogrifyTree(int fileNumber, bool firstCall)
Definition: makeDST.cxx:237
freeze_graph.argv
argv
Definition: freeze_graph.py:218
inTree
TTree * inTree
Definition: makeDST.cxx:37
OUTDEV
#define OUTDEV
globbing
std::vector< std::string > globbing(const std::string &pat)
Definition: globbing.h:16
cet::getenv
std::string getenv(std::string const &name)
Definition: getenv.cc:15
keras_to_tensorflow.int
int
Definition: keras_to_tensorflow.py:69
inFile
TFile * inFile
Definition: makeDST.cxx:36
reco_momentum_tuples.nFiles
int nFiles
Definition: reco_momentum_tuples.py:30
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
void TransMogrifyTree ( int  fileNumber,
bool  firstCall 
)

Definition at line 237 of file makeDST.cxx.

237  {
238 
239 
240 
241  Long64_t iEntry;
242 
243  Int_t fFileNumber;
244  if (firstCall) outTree->Branch("FileNumber",&fFileNumber,"FileNumber/I");
245 
246  scalarFromTree<Int_t> Event (inTree,"Event",&iEntry);
247  Int_t fEvent;
248  if (firstCall) outTree->Branch("Event",&fEvent,"Event/I");
249 
250 
251 
252  vectorFromTree<Int_t> NType(inTree,"NType",&iEntry);
253  vector<Int_t> fNType;
254  if (firstCall) outTree->Branch("NType",&fNType);
255 
256  vectorFromTree<Int_t> Mode(inTree,"Mode",&iEntry);
257  vector<Int_t> fMode;
258  if (firstCall) outTree->Branch("Mode",&fMode);
259 
260  vectorFromTree<int> InterT(inTree,"InterT",&iEntry);
261  vector<int> fInterT;
262  if (firstCall) outTree->Branch("InterT",&fInterT);
263 
264  vectorFromTree<Float_t> MCVertX(inTree,"MCVertX",&iEntry);
265  vector<Float_t> fMCVertX;
266  if (firstCall) outTree->Branch("MCVertX",&fMCVertX);
267 
268  vectorFromTree<Float_t> MCVertY(inTree,"MCVertY",&iEntry);
269  vector<Float_t> fMCVertY;
270  if (firstCall) outTree->Branch("MCVertY",&fMCVertY);
271 
272  vectorFromTree<Float_t> MCVertZ(inTree,"MCVertZ",&iEntry);
273  vector<Float_t> fMCVertZ;
274  if (firstCall) outTree->Branch("MCVertZ",&fMCVertZ);
275 
276  vectorFromTree<Float_t> MCNuPx(inTree,"MCNuPx",&iEntry);
277  vector<Float_t> fMCNuPx;
278  if (firstCall) outTree->Branch("MCNuPx",&fMCNuPx);
279 
280  vectorFromTree<Float_t> MCNuPy(inTree,"MCNuPy",&iEntry);
281  vector<Float_t> fMCNuPy;
282  if (firstCall) outTree->Branch("MCNuPy",&fMCNuPy);
283 
284  vectorFromTree<Float_t> MCNuPz(inTree,"MCNuPz",&iEntry);
285  vector<Float_t> fMCNuPz;
286  if (firstCall) outTree->Branch("MCNuPz",&fMCNuPz);
287 
288  vectorFromTree<Float_t> MC_T(inTree,"MC_T",&iEntry);
289  vector<Float_t> fMC_T;
290  if (firstCall) outTree->Branch("MC_T",&fMC_T);
291 
292 
293 
294  vectorFromTree<Int_t> PDG(inTree,"PDG",&iEntry);
295  vector<Int_t> fPDG;
296  if (firstCall) outTree->Branch("PDG",&fPDG);
297 
298  vectorFromTree<Int_t> PDGMother(inTree,"PDGMother",&iEntry);
299  vector<Int_t> fPDGMother;
300  if (firstCall) outTree->Branch("PDGMother",&fPDGMother);
301 
302  vectorFromTree<Float_t> MCPStartX(inTree,"MCPStartX",&iEntry);
303  vector<Float_t> fMCPStartX;
304  if (firstCall) outTree->Branch("MCPStartX",&fMCPStartX);
305 
306  vectorFromTree<Float_t> MCPStartY(inTree,"MCPStartY",&iEntry);
307  vector<Float_t> fMCPStartY;
308  if (firstCall) outTree->Branch("MCPStartY",&fMCPStartY);
309 
310  vectorFromTree<Float_t> MCPStartZ(inTree,"MCPStartZ",&iEntry);
311  vector<Float_t> fMCPStartZ;
312  if (firstCall) outTree->Branch("MCPStartZ",&fMCPStartZ);
313 
314  vectorFromTree<Float_t> MCPStartPX(inTree,"MCPStartPX",&iEntry);
315  vector<Float_t> fMCPStartPX;
316  if (firstCall) outTree->Branch("MCPStartPX",&fMCPStartPX);
317 
318  vectorFromTree<Float_t> MCPStartPY(inTree,"MCPStartPY",&iEntry);
319  vector<Float_t> fMCPStartPY;
320  if (firstCall) outTree->Branch("MCPStartPY",&fMCPStartPY);
321 
322  vectorFromTree<Float_t> MCPStartPZ(inTree,"MCPStartPZ",&iEntry);
323  vector<Float_t> fMCPStartPZ;
324  if (firstCall) outTree->Branch("MCPStartPZ",&fMCPStartPZ);
325 
326 
327 
328  vectorFromTree<gar::rec::IDNumber>
329  TrackIDNumber(inTree,"TrackIDNumber",&iEntry);
330  vector<gar::rec::IDNumber> fTrackIDNumber;
331  if (firstCall) outTree->Branch("TrackIDNumber",&fTrackIDNumber);
332 
333  vectorFromTree<Float_t> TrackStartX(inTree,"TrackStartX",&iEntry);
334  vector<Float_t> fTrackStartX;
335  if (firstCall) outTree->Branch("TrackStartX",&fTrackStartX);
336 
337  vectorFromTree<Float_t> TrackStartY(inTree,"TrackStartY",&iEntry);
338  vector<Float_t> fTrackStartY;
339  if (firstCall) outTree->Branch("TrackStartY",&fTrackStartY);
340 
341  vectorFromTree<Float_t> TrackStartZ(inTree,"TrackStartZ",&iEntry);
342  vector<Float_t> fTrackStartZ;
343  if (firstCall) outTree->Branch("TrackStartZ",&fTrackStartZ);
344 
345  vectorFromTree<Float_t> TrackStartPX(inTree,"TrackStartPX",&iEntry);
346  vector<Float_t> fTrackStartPX;
347  if (firstCall) outTree->Branch("TrackStartPX",&fTrackStartPX);
348 
349  vectorFromTree<Float_t> TrackStartPY(inTree,"TrackStartPY",&iEntry);
350  vector<Float_t> fTrackStartPY;
351  if (firstCall) outTree->Branch("TrackStartPY",&fTrackStartPY);
352 
353  vectorFromTree<Float_t> TrackStartPZ(inTree,"TrackStartPZ",&iEntry);
354  vector<Float_t> fTrackStartPZ;
355  if (firstCall) outTree->Branch("TrackStartPZ",&fTrackStartPZ);
356 
357  vectorFromTree<Int_t> TrackStartQ(inTree,"TrackStartQ",&iEntry);
358  vector<Int_t> fTrackStartQ;
359  if (firstCall) outTree->Branch("TrackStartQ",&fTrackStartQ);
360 
361  vectorFromTree<Float_t> TrackEndX(inTree,"TrackEndX",&iEntry);
362  vector<Float_t> fTrackEndX;
363  if (firstCall) outTree->Branch("TrackEndX",&fTrackEndX);
364 
365  vectorFromTree<Float_t> TrackEndY(inTree,"TrackEndY",&iEntry);
366  vector<Float_t> fTrackEndY;
367  if (firstCall) outTree->Branch("TrackEndY",&fTrackEndY);
368 
369  vectorFromTree<Float_t> TrackEndZ(inTree,"TrackEndZ",&iEntry);
370  vector<Float_t> fTrackEndZ;
371  if (firstCall) outTree->Branch("TrackEndZ",&fTrackEndZ);
372 
373  vectorFromTree<Float_t> TrackEndPX(inTree,"TrackEndPX",&iEntry);
374  vector<Float_t> fTrackEndPX;
375  if (firstCall) outTree->Branch("TrackEndPX",&fTrackEndPX);
376 
377  vectorFromTree<Float_t> TrackEndPY(inTree,"TrackEndPY",&iEntry);
378  vector<Float_t> fTrackEndPY;
379  if (firstCall) outTree->Branch("TrackEndPY",&fTrackEndPY);
380 
381  vectorFromTree<Float_t> TrackEndPZ(inTree,"TrackEndPZ",&iEntry);
382  vector<Float_t> fTrackEndPZ;
383  if (firstCall) outTree->Branch("TrackEndPZ",&fTrackEndPZ);
384 
385  vectorFromTree<Int_t> TrackEndQ(inTree,"TrackEndQ",&iEntry);
386  vector<Int_t> fTrackEndQ;
387  if (firstCall) outTree->Branch("TrackEndQ",&fTrackEndQ);
388 
389 
390 
391  vectorFromTree<Float_t> TrackLenF(inTree,"TrackLenF",&iEntry);
392  vector<Float_t> fTrackLenF;
393  if (firstCall) outTree->Branch("TrackLenF",&fTrackLenF);
394 
395  vectorFromTree<Float_t> TrackLenB(inTree,"TrackLenB",&iEntry);
396  vector<Float_t> fTrackLenB;
397  if (firstCall) outTree->Branch("TrackLenB",&fTrackLenB);
398 
399  vectorFromTree<Int_t> NTPCClustersOnTrack(inTree,"NTPCClustersOnTrack",&iEntry);
400  vector<Int_t> fNTPCClustersOnTrack;
401  if (firstCall) outTree->Branch("NTPCClustersOnTrack",&fNTPCClustersOnTrack);
402 
403  vectorFromTree<Float_t> TrackAvgIonF(inTree,"TrackAvgIonF",&iEntry);
404  vector<Float_t> fTrackAvgIonF;
405  if (firstCall) outTree->Branch("TrackAvgIonF",&fTrackAvgIonF);
406 
407  vectorFromTree<Float_t> TrackAvgIonB(inTree,"TrackAvgIonB",&iEntry);
408  vector<Float_t> fTrackAvgIonB;
409  if (firstCall) outTree->Branch("TrackAvgIonB",&fTrackAvgIonB);
410 
411 
412 
413  vectorFromTree<gar::rec::IDNumber>
414  VertIDNumber(inTree,"VertIDNumber",&iEntry);
415  vector<gar::rec::IDNumber> fVertIDNumber;
416  if (firstCall) outTree->Branch("VertIDNumber",&fVertIDNumber);
417 
418  vectorFromTree<Float_t> VertX(inTree,"VertX",&iEntry);
419  vector<Float_t> fVertX;
420  if (firstCall) outTree->Branch("VertX",&fVertX);
421 
422  vectorFromTree<Float_t> VertY(inTree,"VertY",&iEntry);
423  vector<Float_t> fVertY;
424  if (firstCall) outTree->Branch("VertY",&fVertY);
425 
426  vectorFromTree<Float_t> VertZ(inTree,"VertZ",&iEntry);
427  vector<Float_t> fVertZ;
428  if (firstCall) outTree->Branch("VertZ",&fVertZ);
429 
430  vectorFromTree<Int_t> VertN(inTree,"VertN",&iEntry);
431  vector<Int_t> fVertN;
432  if (firstCall) outTree->Branch("VertN",&fVertN);
433 
434  vectorFromTree<Int_t> VertQ(inTree,"VertQ",&iEntry);
435  vector<Int_t> fVertQ;
436  if (firstCall) outTree->Branch("VertQ",&fVertQ);
437 
438 
439 
440  vectorFromTree<gar::rec::IDNumber>
441  VT_VertIDNumber(inTree,"VT_VertIDNumber",&iEntry);
442  vector<gar::rec::IDNumber> fVT_VertIDNumber;
443  if (firstCall) outTree->Branch("VT_VertIDNumber",&fVT_VertIDNumber);
444 
445  vectorFromTree<gar::rec::IDNumber>
446  VT_TrackIDNumber(inTree,"VT_TrackIDNumber",&iEntry);
447  vector<gar::rec::IDNumber> fVT_TrackIDNumber;
448  if (firstCall) outTree->Branch("VT_TrackIDNumber",&fVT_TrackIDNumber);
449 
450  vectorFromTree<gar::rec::TrackEnd>
451  VT_TrackEnd(inTree,"VT_TrackEnd",&iEntry);
452  vector<gar::rec::TrackEnd> fVT_TrackEnd;
453  if (firstCall) outTree->Branch("VT_TrackEnd",&fVT_TrackEnd);
454 
455 
456 
457  #ifdef writeECAL
458  vectorFromTree<gar::rec::IDNumber>
459  ClusterIDNumber(inTree,"ClusterIDNumber",&iEntry);
460  vector<gar::rec::IDNumber> fClusterIDNumber;
461  if (firstCall) outTree->Branch("ClusterIDNumber",&fClusterIDNumber);
462 
463  vectorFromTree<UInt_t> ClusterNhits(inTree,"ClusterNhits",&iEntry);
464  vector<UInt_t> fClusterNhits;
465  if (firstCall) outTree->Branch("ClusterNhits",&fClusterNhits);
466 
467  vectorFromTree<Float_t> ClusterEnergy(inTree,"ClusterEnergy",&iEntry);
468  vector<Float_t> fClusterEnergy;
469  if (firstCall) outTree->Branch("ClusterEnergy",&fClusterEnergy);
470 
471  vectorFromTree<Float_t> ClusterX(inTree,"ClusterX",&iEntry);
472  vector<Float_t> fClusterX;
473  if (firstCall) outTree->Branch("ClusterX",&fClusterX);
474 
475  vectorFromTree<Float_t> ClusterY(inTree,"ClusterY",&iEntry);
476  vector<Float_t> fClusterY;
477  if (firstCall) outTree->Branch("ClusterY",&fClusterY);
478 
479  vectorFromTree<Float_t> ClusterZ(inTree,"ClusterZ",&iEntry);
480  vector<Float_t> fClusterZ;
481  if (firstCall) outTree->Branch("ClusterZ",&fClusterZ);
482 
483 
484 
485  vectorFromTree<gar::rec::IDNumber>
486  ECALAssn_ClusIDNumber(inTree,"ECALAssn_ClusIDNumber",&iEntry);
487  vector<gar::rec::IDNumber> fECALAssn_ClusIDNumber;
488  if (firstCall) outTree->Branch("ECALAssn_ClusIDNumber",&fECALAssn_ClusIDNumber);
489 
490  vectorFromTree<gar::rec::IDNumber>
491  ECALAssn_TrackIDNumber(inTree,"ECALAssn_TrackIDNumber",&iEntry);
492  vector<gar::rec::IDNumber> fECALAssn_TrackIDNumber;
493  if (firstCall) outTree->Branch("ECALAssn_TrackIDNumber",&fECALAssn_TrackIDNumber);
494 
495  vectorFromTree<gar::rec::TrackEnd>
496  ECALAssn_TrackEnd(inTree,"ECALAssn_TrackEnd",&iEntry);
497  vector<gar::rec::TrackEnd> fECALAssn_TrackEnd;
498  if (firstCall) outTree->Branch("ECALAssn_TrackEnd",&fECALAssn_TrackEnd);
499  #endif
500 
501 
502 
503 
504  Long64_t nEntry = inTree->GetEntries();
505  cout << "nEntry: " << nEntry << "\t Run Number: " << fileNumber << endl;
506  for (iEntry=0; iEntry<nEntry; ++iEntry) {
507 
508  ++Counter[nEvent];
509 
510 
511 
512  inTree->GetEntry(iEntry);
513 
514  fFileNumber = fileNumber;
515  fEvent = Event.getData();
516 
517  fNType = NType.getDataVector();
518  fMode = Mode.getDataVector();
519  fInterT = InterT.getDataVector();
520  fMCVertX = MCVertX.getDataVector();
521  fMCVertY = MCVertY.getDataVector();
522  fMCVertZ = MCVertZ.getDataVector();
523  fMCNuPx = MCNuPx.getDataVector();
524  fMCNuPy = MCNuPy.getDataVector();
525  fMCNuPz = MCNuPz.getDataVector();
526  fMC_T = MC_T.getDataVector();
527 
528  fPDG = PDG.getDataVector();
529  fPDGMother = PDGMother.getDataVector();
530  fMCPStartX = MCPStartX.getDataVector();
531  fMCPStartY = MCPStartY.getDataVector();
532  fMCPStartZ = MCPStartZ.getDataVector();
533  fMCPStartPX = MCPStartPX.getDataVector();
534  fMCPStartPY = MCPStartPY.getDataVector();
535  fMCPStartPZ = MCPStartPZ.getDataVector();
536 
537  fTrackIDNumber = TrackIDNumber.getDataVector();
538  fTrackStartX = TrackStartX.getDataVector();
539  fTrackStartY = TrackStartY.getDataVector();
540  fTrackStartZ = TrackStartZ.getDataVector();
541  fTrackStartPX = TrackStartPX.getDataVector();
542  fTrackStartPY = TrackStartPY.getDataVector();
543  fTrackStartPZ = TrackStartPZ.getDataVector();
544  fTrackStartQ = TrackStartQ.getDataVector();
545  fTrackEndX = TrackEndX.getDataVector();
546  fTrackEndY = TrackEndY.getDataVector();
547  fTrackEndZ = TrackEndZ.getDataVector();
548  fTrackEndPX = TrackEndPX.getDataVector();
549  fTrackEndPY = TrackEndPY.getDataVector();
550  fTrackEndPZ = TrackEndPZ.getDataVector();
551  fTrackEndQ = TrackEndQ.getDataVector();
552 
553  fTrackLenF = TrackLenF.getDataVector();
554  fTrackLenB = TrackLenB.getDataVector();
555  fNTPCClustersOnTrack = NTPCClustersOnTrack.getDataVector();
556  fTrackAvgIonF = TrackAvgIonF.getDataVector();
557  fTrackAvgIonB = TrackAvgIonB.getDataVector();
558 
559  fVertIDNumber = VertIDNumber.getDataVector();
560  fVertX = VertX.getDataVector();
561  fVertY = VertY.getDataVector();
562  fVertZ = VertZ.getDataVector();
563  fVertN = VertN.getDataVector();
564  fVertQ = VertQ.getDataVector();
565 
566  fVT_VertIDNumber = VT_VertIDNumber.getDataVector();
567  fVT_TrackIDNumber = VT_TrackIDNumber.getDataVector();
568  fVT_TrackEnd = VT_TrackEnd.getDataVector();
569 
570  #ifdef writeECAL
571  fClusterIDNumber = ClusterIDNumber.getDataVector();
572  fClusterNhits = ClusterNhits.getDataVector();
573  fClusterEnergy = ClusterEnergy.getDataVector();
574  fClusterX = ClusterX.getDataVector();
575  fClusterY = ClusterY.getDataVector();
576  fClusterZ = ClusterZ.getDataVector();
577 
578  fECALAssn_ClusIDNumber = ECALAssn_ClusIDNumber.getDataVector();
579  fECALAssn_TrackIDNumber = ECALAssn_TrackIDNumber.getDataVector();
580  fECALAssn_TrackEnd = ECALAssn_TrackEnd.getDataVector();
581  #endif
582 
583 
584 
585 
586  Float_t XvertMC, YvertMC, ZvertMC;
587  int nMCVerts = MCVertX.size();
588  for (int iMCVert=0; iMCVert<nMCVerts; ++iMCVert) {
589  XvertMC = MCVertX.getData(iMCVert);
590  YvertMC = MCVertY.getData(iMCVert);
591  ZvertMC = MCVertZ.getData(iMCVert);
592  if ( inFiducial(XvertMC,YvertMC,ZvertMC) ) {
593  ++Counter[nDatainFid];
594  int interType = InterT.getData(iMCVert);
595  if (interType==simb::kCCQE) ++Counter[nCCQEinFid];
596  if (interType==simb::kNCQE) ++Counter[nNCQEinFid];
597  if (interType==simb::kNuanceOffset) ++Counter[nRESinFid];
598  if (simb::kResCCNuProtonPiPlus<=interType &&
599  interType<=simb::kResCCNuBarProtonPi0Pi0) ++Counter[nRESinFid];
600  if (interType==simb::kCCDIS) ++Counter[nCCDISinFid];
601  if (interType==simb::kNCDIS) ++Counter[nNCDISinFid];
602  if (interType==simb::kCCCOH) ++Counter[nCCCOHinFid];
603  if (interType==simb::kNCCOH) ++Counter[nNCCOHinFid];
604  if (interType==simb::kNuElectronElastic) ++Counter[nNuEinFid];
605  }
606  }
607 
608 
609 
610 
611  #ifdef killNoRecoVerts
612  if (fVertIDNumber.size()==0) continue;
613  #endif
614 
615  // Setup iterators for the vector.erase function on the tracks
616  vector<gar::rec::IDNumber>::iterator itTrackIDNumber = fTrackIDNumber.begin();
617  vector<Float_t>::iterator itTrackStartX = fTrackStartX.begin();
618  vector<Float_t>::iterator itTrackStartY = fTrackStartY.begin();
619  vector<Float_t>::iterator itTrackStartZ = fTrackStartZ.begin();
620  vector<Float_t>::iterator itTrackStartPX = fTrackStartPX.begin();
621  vector<Float_t>::iterator itTrackStartPY = fTrackStartPY.begin();
622  vector<Float_t>::iterator itTrackStartPZ = fTrackStartPZ.begin();
623  vector<Int_t>::iterator itTrackStartQ = fTrackStartQ.begin();
624  vector<Float_t>::iterator itTrackEndX = fTrackEndX.begin();
625  vector<Float_t>::iterator itTrackEndY = fTrackEndY.begin();
626  vector<Float_t>::iterator itTrackEndZ = fTrackEndZ.begin();
627  vector<Float_t>::iterator itTrackEndPX = fTrackEndPX.begin();
628  vector<Float_t>::iterator itTrackEndPY = fTrackEndPY.begin();
629  vector<Float_t>::iterator itTrackEndPZ = fTrackEndPZ.begin();
630  vector<Int_t>::iterator itTrackEndQ = fTrackEndQ.begin();
631  vector<Float_t>::iterator itTrackLenF = fTrackLenF.begin();
632  vector<Float_t>::iterator itTrackLenB = fTrackLenB.begin();
633  vector<Int_t>::iterator itNTPCClustersOnTrack = fNTPCClustersOnTrack.begin();
634  vector<Float_t>::iterator itTrackAvgIonF = fTrackAvgIonF.begin();
635  vector<Float_t>::iterator itTrackAvgIonB = fTrackAvgIonB.begin();
636 
637  // Step through the vectors, deleting as you go (this is slow, ok) be sure to keep
638  // track of the reduced lengths
639  int nTracks = fTrackIDNumber.size();
640  for (int iTrack=0; iTrack<nTracks; ++iTrack) {
641  Float_t Pstart = Qadd(fTrackStartPX[iTrack],fTrackStartPY[iTrack],fTrackStartPZ[iTrack]);
642  Float_t Pend = Qadd(fTrackEndPX[iTrack], fTrackEndPY[iTrack], fTrackEndPZ[iTrack]);
643  if ( Pstart<PmagCut && Pend<PmagCut ) {
644  itTrackIDNumber = fTrackIDNumber.erase(itTrackIDNumber);
645  itTrackStartX = fTrackStartX.erase(itTrackStartX);
646  itTrackStartY = fTrackStartY.erase(itTrackStartY);
647  itTrackStartZ = fTrackStartZ.erase(itTrackStartZ);
648  itTrackStartPX = fTrackStartPX.erase(itTrackStartPX);
649  itTrackStartPY = fTrackStartPY.erase(itTrackStartPY);
650  itTrackStartPZ = fTrackStartPZ.erase(itTrackStartPZ);
651  itTrackStartQ = fTrackStartQ.erase(itTrackStartQ);
652  itTrackEndX = fTrackEndX.erase(itTrackEndX);
653  itTrackEndY = fTrackEndY.erase(itTrackEndY);
654  itTrackEndZ = fTrackEndZ.erase(itTrackEndZ);
655  itTrackEndPX = fTrackEndPX.erase(itTrackEndPX);
656  itTrackEndPY = fTrackEndPY.erase(itTrackEndPY);
657  itTrackEndPZ = fTrackEndPZ.erase(itTrackEndPZ);
658  itTrackEndQ = fTrackEndQ.erase(itTrackEndQ);
659  itTrackLenF = fTrackLenF.erase(itTrackLenF);
660  itTrackLenB = fTrackLenB.erase(itTrackLenB);
661  itNTPCClustersOnTrack = fNTPCClustersOnTrack.erase(itNTPCClustersOnTrack);
662  itTrackAvgIonF = fTrackAvgIonF.erase(itTrackAvgIonF);
663  itTrackAvgIonB = fTrackAvgIonB.erase(itTrackAvgIonB);
664  --nTracks; --iTrack; // for (... iTrack<nTrack...) will break when iTrack==-1, nTrack==0
665  } else {
666  ++itTrackIDNumber;
667  ++itTrackStartX; ++itTrackStartY; ++itTrackStartZ;
668  ++itTrackStartPX; ++itTrackStartPY; ++itTrackStartPZ; ++itTrackStartQ;
669  ++itTrackEndX; ++itTrackEndY; ++itTrackEndZ;
670  ++itTrackEndPX; ++itTrackEndPY; ++itTrackEndPZ; ++itTrackEndQ;
671  ++itTrackLenF; ++itTrackLenB; ++itNTPCClustersOnTrack;
672  ++itTrackAvgIonF; ++itTrackAvgIonB;
673  }
674  }
675  if (fTrackIDNumber.size()==0) continue;
676 
677 
678 
679  // Find the vertices with vert-track associations that are missing one or more tracks
680  std::list<gar::rec::IDNumber> deadVertex;
681  int nAssns = fVT_VertIDNumber.size();
682  for (int iAssn=0; iAssn<nAssns; ++iAssn) {
683  gar::rec::IDNumber thisIDNumber = fVT_TrackIDNumber[iAssn];
684  vector<gar::rec::IDNumber>::iterator found;
685  found = std::find(fTrackIDNumber.begin(),fTrackIDNumber.end(), thisIDNumber);
686  if ( found == fTrackIDNumber.end() ) {
687  // This vertex-track association is dead. Save the dead vertex number
688  deadVertex.push_back(fVT_VertIDNumber[iAssn]);
689  }
690  }
691  if (deadVertex.size()>0) {
692  // Get the list cleaned up
693  std::list<gar::rec::IDNumber>::iterator deadStart = deadVertex.begin();
694  std::list<gar::rec::IDNumber>::iterator deadEnd;
695  deadVertex.sort();
696  deadEnd = std::unique(deadVertex.begin(),deadVertex.end());
697 
698 
699 
700  // Delete the dead vertices
701  vector<gar::rec::IDNumber>::iterator itVertIDNumber = fVertIDNumber.begin();
702  vector<Float_t>::iterator itVertX = fVertX.begin();
703  vector<Float_t>::iterator itVertY = fVertY.begin();
704  vector<Float_t>::iterator itVertZ = fVertZ.begin();
705  vector<Int_t>::iterator itVertN = fVertN.begin();
706  vector<Int_t>::iterator itVertQ = fVertQ.begin();
707  int nVerts = fVertIDNumber.size();
708  for (int iVert=0; iVert<nVerts; ++iVert) {
709  gar::rec::IDNumber thisIDNumber = fVertIDNumber[iVert];
710  std::list<gar::rec::IDNumber>::iterator found;
711  found = std::find(deadStart,deadEnd, thisIDNumber);
712  if ( found != deadEnd ) {
713  // Indeed, this is a dead vertex. Delete it
714  itVertIDNumber = fVertIDNumber.erase(itVertIDNumber);
715  itVertX = fVertX.erase(itVertX);
716  itVertY = fVertY.erase(itVertY);
717  itVertZ = fVertZ.erase(itVertZ);
718  itVertN = fVertN.erase(itVertN);
719  itVertQ = fVertQ.erase(itVertQ);
720  --nVerts; --iVert;
721  } else {
722  ++itVertIDNumber;
723  ++itVertX; ++itVertY; ++itVertZ; ++itVertN; ++itVertQ;
724  }
725  }
726  #ifdef killNoRecoVerts
727  if (fVertIDNumber.size()==0) continue;
728  #endif
729 
730 
731 
732  // Delete the dead vert-track associations
733  vector<gar::rec::IDNumber>::iterator itVT_VertIDNumber = fVT_VertIDNumber.begin();
734  vector<gar::rec::IDNumber>::iterator itVT_TrackIDNumber = fVT_TrackIDNumber.begin();
735  vector<gar::rec::TrackEnd>::iterator itVT_TrackEnd = fVT_TrackEnd.begin();
736  for (int iAssn=0; iAssn<nAssns; ++iAssn) {
737  gar::rec::IDNumber thisIDNumber = fVT_VertIDNumber[iAssn];
738  std::list<gar::rec::IDNumber>::iterator found;
739  found = std::find(deadStart,deadEnd, thisIDNumber);
740  if ( found != deadEnd ) {
741  // Indeed, this is a dead track-vertex association
742  itVT_VertIDNumber = fVT_VertIDNumber.erase(itVT_VertIDNumber);
743  itVT_TrackIDNumber = fVT_TrackIDNumber.erase(itVT_TrackIDNumber);
744  itVT_TrackEnd = fVT_TrackEnd.erase(itVT_TrackEnd);
745  --nAssns; --iAssn;
746  } else {
747  ++itVT_VertIDNumber; ++itVT_TrackIDNumber; ++itVT_TrackEnd;
748  }
749  }
750 
751 
752 
753  #ifdef writeECAL
754  // Delete the dead ECAL-track associations, but not the ECAL cluster info
755  vector<gar::rec::IDNumber>::iterator itECALAssn_ClusIDNumber = fECALAssn_ClusIDNumber.begin();
756  vector<gar::rec::IDNumber>::iterator itECALAssn_TrackIDNumber = fECALAssn_TrackIDNumber.begin();
757  vector<gar::rec::TrackEnd>::iterator itECALAssn_TrackEnd = fECALAssn_TrackEnd.begin();
758  nAssns = fECALAssn_ClusIDNumber.size();
759  for (int iAssn=0; iAssn<nAssns; ++iAssn) {
760  gar::rec::IDNumber thisIDNumber = fECALAssn_TrackIDNumber[iAssn];
761  vector<gar::rec::IDNumber>::iterator found;
762  found = std::find(fTrackIDNumber.begin(),fTrackIDNumber.end(), thisIDNumber);
763  if ( found == fTrackIDNumber.end() ) {
764  // Indeed, this is a dead track-vertex association
765  itECALAssn_ClusIDNumber = fECALAssn_ClusIDNumber.erase(itECALAssn_ClusIDNumber);
766  itECALAssn_TrackIDNumber = fECALAssn_TrackIDNumber.erase(itECALAssn_TrackIDNumber);
767  itECALAssn_TrackEnd = fECALAssn_TrackEnd.erase(itECALAssn_TrackEnd);
768  --nAssns; --iAssn;
769  } else {
770  ++itECALAssn_ClusIDNumber; ++itECALAssn_TrackIDNumber; ++itECALAssn_TrackEnd;
771  }
772  }
773  #endif
774  }
775 
776 
777 
778 
779  outTree->Fill();
780  }
781  return;
782 }
simb::kNCQE
neutral current quasi-elastic
Definition: MCNeutrino.h:97
iterator
intermediate_table::iterator iterator
Definition: intermediate_table.cc:27
simb::kNCDIS
charged current deep inelastic scatter
Definition: MCNeutrino.h:134
PmagCut
float PmagCut
Definition: makeDST.cxx:32
simb::kNuElectronElastic
neutrino electron elastic scatter
Definition: MCNeutrino.h:140
inFiducial
bool inFiducial(double x, double y, double z)
Definition: makeDST.cxx:17
Qadd
double Qadd(double const x, double const y)
Definition: includeROOT.h:77
outTree
TTree * outTree
Definition: makeDST.cxx:37
Counter
int Counter[nCounterTag]
Definition: makeDST.cxx:47
simb::kNuanceOffset
offset to account for adding in Nuance codes to this enum
Definition: MCNeutrino.h:95
simb::kCCQE
charged current quasi-elastic
Definition: MCNeutrino.h:96
PDG
const uint PDG
Definition: qregexp.cpp:140
inTree
TTree * inTree
Definition: makeDST.cxx:37
simb::kCCDIS
charged current deep inelastic scatter
Definition: MCNeutrino.h:133
simb::kResCCNuProtonPiPlus
resonant charged current, nu p -> l- p pi+
Definition: MCNeutrino.h:98
simb::kCCCOH
charged current coherent pion
Definition: MCNeutrino.h:139
Event
Definition: types.h:32
gar::rec::IDNumber
size_t IDNumber
Definition: IDNumberGen.h:71
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030

Variable Documentation

int Counter[nCounterTag]

Definition at line 47 of file makeDST.cxx.

TFile* inFile

Definition at line 36 of file makeDST.cxx.

TTree* inTree

Definition at line 37 of file makeDST.cxx.

TFile* outFile

Definition at line 36 of file makeDST.cxx.

TTree* outTree

Definition at line 37 of file makeDST.cxx.

float PmagCut = 0.150

Definition at line 32 of file makeDST.cxx.

bool viaXrootd = true

Definition at line 31 of file makeDST.cxx.

Generated by   doxygen 1.8.11