makeDST.cxx
Go to the documentation of this file.
1 #include <list>
2 
3 #include "globbing.h"
4 #include "treeReader.h"
5 #include "includeROOT.h"
6 
7 #include "nusimdata/SimulationBase/MCNeutrino.h"
8 #include "ReconstructionDataProducts/Track.h"
9 #include "ReconstructionDataProducts/IDNumberGen.h"
10 
11 void TransMogrifyTree(int fileNumber, bool firstCall);
12 
13 
14 
15 // In a fiducial corresponding to 1 metric ton of Ar nuclei? (Corresponds to
16 // about 35cm from edge of the TPC).
17 bool inFiducial(double x, double y, double z) {
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 }
24 
25 
26 
27 //==============================================================================
28 //==============================================================================
29 #define writeECAL true
30 #define killNoRecoVerts true
31 bool viaXrootd = true;
32 float PmagCut = 0.150; // In GeV.
33 
34 
35 
36 TFile* outFile; TFile* inFile;
37 TTree* outTree; TTree* inTree;
38 
39 
40 
41 enum CounterTag {nEvent, nDatainFid,
42  nCCQEinFid, nNCQEinFid,
43  nRESinFid, nCCDISinFid,
44  nNCDISinFid, nCCCOHinFid,
45  nNCCOHinFid, nNuEinFid,
46 nCounterTag}; // Must be last enum
47 int Counter[nCounterTag];
48 
49 
50 
51 //==============================================================================
52 //==============================================================================
53 int main(int argc , const char* argv[]){
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 }
230 
231 
232 
233 
234 
235 //==============================================================================
236 //==============================================================================
237 void TransMogrifyTree(int fileNumber, bool firstCall) {
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
main
int main(int argc, const char *argv[])
Definition: makeDST.cxx:53
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
vectorFromTree::getDataVector
vector< T > & getDataVector()
Definition: treeReader.h:140
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
viaXrootd
bool viaXrootd
Definition: makeDST.cxx:31
simb::kCCQE
charged current quasi-elastic
Definition: MCNeutrino.h:96
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
PDG
const uint PDG
Definition: qregexp.cpp:140
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
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
keras_to_tensorflow.int
int
Definition: keras_to_tensorflow.py:69
inFile
TFile * inFile
Definition: makeDST.cxx:36
simb::kCCCOH
charged current coherent pion
Definition: MCNeutrino.h:139
Event
Definition: types.h:32
train.x
list x
Definition: train.py:276
reco_momentum_tuples.nFiles
int nFiles
Definition: reco_momentum_tuples.py:30
vectorFromTree::getData
T getData(int i)
Definition: treeReader.h:161
gar::rec::IDNumber
size_t IDNumber
Definition: IDNumberGen.h:71
CounterTag
CounterTag
Definition: makeDST.cxx:41
plot_model.r
r
Definition: plot_model.py:57
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
Generated by   doxygen 1.8.11