protonmccnn_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: protonmccnn
3 // File: protonmccnn_module.cc
4 //
5 // Extract protoDUNE useful information, do a firs tpre-selection and save output to a flat tree
6 //
7 // Some parts are copied from the beam module example
8 //
9 // Georgios Christodoulou - georgios.christodoulou at cern.ch
10 // Heng-Ye Liao modified it for his protonanalysis - liao@phys.ksu.edu
11 //
12 ////////////////////////////////////////////////////////////////////////
13 
14 #include "art/Framework/Core/EDAnalyzer.h"
15 #include "art/Framework/Core/ModuleMacros.h"
16 #include "art/Framework/Principal/Event.h"
17 #include "art/Framework/Principal/Handle.h"
18 #include "art/Framework/Principal/Run.h"
19 #include "art/Framework/Principal/SubRun.h"
20 //#include "art/Framework/Services/Optional/TFileService.h"
21 #include "art_root_io/TFileService.h"
22 //#include "art_root_io/TFileDirectory.h"
23 #include "canvas/Utilities/InputTag.h"
24 #include "canvas/Persistency/Common/FindManyP.h"
25 #include "fhiclcpp/ParameterSet.h"
26 #include "messagefacility/MessageLogger/MessageLogger.h"
27 #include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
28 
29 #include "larsim/MCCheater/BackTrackerService.h"
30 #include "larsim/MCCheater/ParticleInventoryService.h"
31 #include "lardataobj/RecoBase/Hit.h"
32 #include "lardataobj/RecoBase/Track.h"
33 #include "lardataobj/RecoBase/Shower.h"
34 #include "lardataobj/RecoBase/TrackHitMeta.h"
35 #include "lardataobj/RecoBase/PFParticle.h"
36 #include "nusimdata/SimulationBase/MCParticle.h"
37 #include "nusimdata/SimulationBase/MCTruth.h"
38 #include "lardataobj/AnalysisBase/CosmicTag.h"
39 #include "lardataobj/AnalysisBase/T0.h"
40 #include "larreco/RecoAlg/TrackMomentumCalculator.h"
41 #include "lardataobj/AnalysisBase/Calorimetry.h"
42 
43 #include "dune/DuneObj/ProtoDUNEBeamEvent.h"
44 //#include "dune/EventGenerator/ProtoDUNEbeamDataProducts/ProtoDUNEbeamsim.h"
45 //#include "dune/EventGenerator/ProtoDUNEbeamDataProducts/ProtoDUNEBeamInstrument.h"
46 
47 #include "lardata/ArtDataHelper/MVAReader.h"
48 
49 #include "protoduneana/Utilities/ProtoDUNETrackUtils.h"
50 #include "protoduneana/Utilities/ProtoDUNEShowerUtils.h"
51 #include "protoduneana/Utilities/ProtoDUNETruthUtils.h"
52 #include "protoduneana/Utilities/ProtoDUNEPFParticleUtils.h"
53 #include "dune/Protodune/singlephase/DataUtils/ProtoDUNEDataUtils.h"
54 
55 #include "larsim/MCCheater/ParticleInventoryService.h"
56 
57 // ROOT includes
58 #include "TTree.h"
59 #include "TFile.h"
60 #include "TString.h"
61 #include "TTimeStamp.h"
62 
63 // C++ Includes
64 #include <stdio.h>
65 #include <stdlib.h>
66 #include <fstream>
67 #include <string>
68 #include <sstream>
69 #include <cmath>
70 #include <algorithm>
71 #include <iostream>
72 #include <vector>
73 
74 // Maximum number of beam particles to save
75 const int NMAXDAUGTHERS = 15;
76 
77 // const int kMaxTracks = 1000;
78 // const int kMaxHits = 10000;
79 //double prim_energy=0.0;
80 
81 namespace protoana {
82  class protonmccnn;
83 }
84 
85 
86 class protoana::protonmccnn : public art::EDAnalyzer {
87 public:
88 
89  explicit protonmccnn(fhicl::ParameterSet const & p);
90 
91  protonmccnn(protonmccnn const &) = delete;
92  protonmccnn(protonmccnn &&) = delete;
93  protonmccnn & operator = (protonmccnn const &) = delete;
94  protonmccnn & operator = (protonmccnn &&) = delete;
95 
96  virtual void beginJob() override;
97  virtual void endJob() override;
98 
99  // Required functions.
100  void analyze(art::Event const & evt) override;
101 
102 private:
103 
104  // Helper utility functions
105  protoana::ProtoDUNEDataUtils dataUtil;
106  protoana::ProtoDUNEPFParticleUtils pfpUtil;
107  protoana::ProtoDUNETrackUtils trackUtil;
108  protoana::ProtoDUNETruthUtils truthUtil;
109 
110  // Initialise tree variables
111  void Initialise();
112 
113  // Fill cosmics tree
114  void FillCosmicsTree(art::Event const & evt, std::string pfParticleTag);
115 
116  // fcl parameters
117  const art::InputTag fNNetModuleLabel; //label of the module used for CNN tagging
118  const art::InputTag fBeamModuleLabel;
119  std::string fCalorimetryTag;
120  std::string fTrackerTag;
121  std::string fShowerTag;
122  std::string fPFParticleTag;
123  std::string fGeneratorTag;
124  bool fVerbose;
125 
126  //geometry
127  geo::GeometryCore const * fGeometry;
128 
129  double MCTruthT0, TickT0;
130  int nT0s;
131 
132  TTree *fPandoraBeam;
133  //TTree *fPandoraCosmics;
134 
135  // Tree variables
136  int fRun;
137  int fSubRun;
138  int fevent;
139  double fTimeStamp;
140  int fNactivefembs[6];
141 
142  // Beam track
143  int fbeamtrigger;
144  double ftof;
145  int fcerenkov1;
146  int fcerenkov2;
147  double fbeamtrackMomentum;
148  double fbeamtrackP[3]; //Px/Py/Pz
149  double fbeamtrackEnergy;
150  double fbeamtrackPos[3];
151  double fbeamtrackDir[3];
152  double fbeamtrackTime;
153  int fbeamtrackPdg;
154  int fbeamtrackID;
155  //int NumberBeamTrajectoryPoints;
156  std::vector<double> beamtrk_x;
157  std::vector<double> beamtrk_y;
158  std::vector<double> beamtrk_z;
159  std::vector<double> beamtrk_Px;
160  std::vector<double> beamtrk_Py;
161  std::vector<double> beamtrk_Pz;
162  std::vector<double> beamtrk_Eng;
163 
164  //Front face KE & Positions
165  double ke_ff;
166  std::vector<double> pos_ff;
167 
168  //KE_true and KE_reco
169  std::vector<double> primtrk_ke_true;
170  std::vector<double> primtrk_ke_reco;
171  std::string primtrk_end_g4process;
172  std::vector<double> primtrk_range_true;
173 
174  //MCS
175  //std::vector<float> primtrk_mcs_angles_reco;
176  //std::vector<float> primtrk_mcs_radlengths_reco;
177  //std::vector<float> primtrk_mcs_angles_true;
178  //std::vector<float> primtrk_mcs_radlengths_true;
179 
180  //true info
181  double ke_ff_true;
182  //std::vector<unsigned char> primtrk_hit_processkey;
183  std::vector< std::vector<double> > primtrk_hitx_true;
184  std::vector< std::vector<double> > primtrk_hity_true;
185  std::vector< std::vector<double> > primtrk_hitz_true;
186  std::vector< std::vector<double> > primtrk_trkid_true;
187  std::vector< std::vector<double> > primtrk_edept_true;
188 
189  //save all the true info
190  std::vector< std::vector<double> > primtrk_true_x;
191  std::vector< std::vector<double> > primtrk_true_y;
192  std::vector< std::vector<double> > primtrk_true_z;
193  std::vector< std::vector<double> > primtrk_true_trkid;
194  std::vector< std::vector<double> > primtrk_true_edept;
195 
196  //mctraj info
197  std::vector< std::vector<double> > primtrj_true_x;
198  std::vector< std::vector<double> > primtrj_true_y;
199  std::vector< std::vector<double> > primtrj_true_z;
200  std::vector< std::vector<double> > primtrj_true_edept;
201 
202  //hits from pfparticles
203  std::vector<double> pfphit_peaktime_c;
204  std::vector<double> pfphit_peaktime_u;
205  std::vector<double> pfphit_peaktime_v;
206  std::vector<double> pfphit_wireid_c;
207  std::vector<double> pfphit_wireid_u;
208  std::vector<double> pfphit_wireid_v;
209 
210 
211  // Reconstructed tracks/showers
212  double fvertex[3];
213  double fsecvertex[3];
214 
215  int fisprimarytrack;
216  int fisprimaryshower;
217  double fprimaryBDTScore;
218  int fprimaryNHits;
219  double fprimaryTheta;
220  double fprimaryPhi;
221  double fprimaryLength;
222  double fprimaryMomentum;
223  double fprimaryEndMomentum;
224  double fprimaryEndPosition[3];
225  double fprimaryStartPosition[3];
226  double fprimaryEndDirection[3];
227  double fprimaryStartDirection[3];
228  double fprimaryOpeningAngle;
229  int fprimaryShowerBestPlane;
230  double fprimaryShowerEnergy;
231  double fprimaryShowerMIPEnergy;
232  double fprimaryShowerdEdx;
233  double fprimaryMomentumByRangeProton;
234  double fprimaryMomentumByRangeMuon;
235  double fprimaryKineticEnergy[3];
236  double fprimaryRange[3];
237  int fprimaryID;
238  double fprimaryT0;
239 
240  //Truth info
241  int ftruthpdg;
242  int fprimary_truth_TrackId;
243  int fprimary_truth_Pdg;
244  double fprimary_truth_StartPosition[4];
245  double fprimary_truth_StartPosition_MC[4];
246  double fprimary_truth_EndPosition[4];
247  double fprimary_truth_EndPosition_MC[4];
248  double fprimary_truth_P;
249  double fprimary_truth_Momentum[4];
250  double fprimary_truth_EndMomentum[4];
251  //std::string G4Process_Primary_End;
252  double fprimary_truth_Pt;
253  double fprimary_truth_Mass;
254  double fprimary_truth_Theta;
255  double fprimary_truth_Phi;
256  //int fprimary_truth_Process;
257  //std::string fprimary_truth_Process;
258  //std::string fprimary_truth_EndProcess;
259  //std::string fprimary_backtrker_truth_Process;
260  std::string fprimary_truth_EndProcess;
261  int fprimary_truth_Isbeammatched;
262  int fprimary_truth_NDaughters;
263 
264  //carlo info
265  std::vector< std::vector<double> > primtrk_dqdx;
266  std::vector< std::vector<double> > primtrk_resrange;
267  std::vector< std::vector<double> > primtrk_dedx;
268  std::vector<double> primtrk_range;
269  std::vector< std::vector<double> > primtrk_hitx;
270  std::vector< std::vector<double> > primtrk_hity;
271  std::vector< std::vector<double> > primtrk_hitz;
272  std::vector< std::vector<double> > primtrk_pitch;
273 
274  //hit and CNN score
275 
276  //cnn score
277  std::vector<float> inelscore_c;
278  std::vector<float> elscore_c;
279  std::vector<float> nonescore_c;
280  //associtated (x,y,z)
281  std::vector<float> x_c;
282  std::vector<float> y_c;
283  std::vector<float> z_c;
284 
285 
286  //(x,y,z) of highest CNN score
287  //double xyz_inelscore_c[3];
288  //double xyz_elscore_c[3];
289 
290 
291  //traj hit(x,y,z) to hit(wid,pt)
292  std::vector<double> wid_c;
293  std::vector<double> tt_c;
294  std::vector<double> wid_v;
295  std::vector<double> tt_v;
296  std::vector<double> wid_u;
297  std::vector<double> tt_u;
298 
299 
300  //info of interactions
301  std::vector<double> interactionX;
302  std::vector<double> interactionY;
303  std::vector<double> interactionZ;
304  std::vector<double> interactionE;
305  std::vector<std::string> interactionProcesslist;
306  std::vector<double> interactionAngles;
307  std::vector<double> Zintersection;
308  std::vector<double> Yintersection;
309  std::vector<double> Xintersection;
310  std::vector<double> timeintersection;
311 
312  std::vector<double> interaction_wid_c;
313  std::vector<double> interaction_tt_c;
314  std::vector<double> interaction_wid_v;
315  std::vector<double> interaction_tt_v;
316  std::vector<double> interaction_wid_u;
317  std::vector<double> interaction_tt_u;
318 
319  // Daughters from primary
320  int fNDAUGHTERS;
321  int fisdaughtertrack[NMAXDAUGTHERS];
322  int fisdaughtershower[NMAXDAUGTHERS];
323  int fdaughterNHits[NMAXDAUGTHERS];
324  double fdaughterTheta[NMAXDAUGTHERS];
325  double fdaughterPhi[NMAXDAUGTHERS];
326  double fdaughterLength[NMAXDAUGTHERS];
327  double fdaughterMomentum[NMAXDAUGTHERS];
328  double fdaughterEndMomentum[NMAXDAUGTHERS];
329  double fdaughterEndPosition[NMAXDAUGTHERS][3];
330  double fdaughterStartPosition[NMAXDAUGTHERS][3];
331  double fdaughterEndDirection[NMAXDAUGTHERS][3];
332  double fdaughterStartDirection[NMAXDAUGTHERS][3];
333  double fdaughterOpeningAngle[NMAXDAUGTHERS];
334  double fdaughterShowerEnergy[NMAXDAUGTHERS];
335  double fdaughterShowerMIPEnergy[NMAXDAUGTHERS];
336  double fdaughterShowerdEdx[NMAXDAUGTHERS];
337  int fdaughterShowerBestPlane[NMAXDAUGTHERS];
338  double fdaughterMomentumByRangeProton[NMAXDAUGTHERS];
339  double fdaughterMomentumByRangeMuon[NMAXDAUGTHERS];
340  double fdaughterKineticEnergy[NMAXDAUGTHERS][3];
341  double fdaughterRange[NMAXDAUGTHERS][3];
342  int fdaughterID[NMAXDAUGTHERS];
343  //double fdaughterT0[NMAXDAUGTHERS];
344 
345  int fdaughter_truth_TrackId[NMAXDAUGTHERS];
346  int fdaughter_truth_Pdg[NMAXDAUGTHERS];
347  double fdaughter_truth_StartPosition[NMAXDAUGTHERS][4];
348  double fdaughter_truth_EndPosition[NMAXDAUGTHERS][4];
349  double fdaughter_truth_P[NMAXDAUGTHERS];
350  double fdaughter_truth_Momentum[NMAXDAUGTHERS][4];
351  double fdaughter_truth_EndMomentum[NMAXDAUGTHERS][4];
352  double fdaughter_truth_Pt[NMAXDAUGTHERS];
353  double fdaughter_truth_Mass[NMAXDAUGTHERS];
354  double fdaughter_truth_Theta[NMAXDAUGTHERS];
355  double fdaughter_truth_Phi[NMAXDAUGTHERS];
356  int fdaughter_truth_Process[NMAXDAUGTHERS];
357 
358  double minX = -360.0;
359  double maxX = 360.0;
360  double minY =0.0;
361  double maxY = 600.0;
362  double minZ = 0.0;
363  double maxZ = 695.0;
364 
365 };
366 
367 
368 protoana::protonmccnn::protonmccnn(fhicl::ParameterSet const & p)
369  :
370  EDAnalyzer(p),
371  dataUtil(p.get<fhicl::ParameterSet>("DataUtils")),
372  fNNetModuleLabel(p.get<art::InputTag>("NNetModuleLabel")),
373  fBeamModuleLabel(p.get< art::InputTag >("BeamModuleLabel")),
374  fCalorimetryTag(p.get<std::string>("CalorimetryTag")),
375  fTrackerTag(p.get<std::string>("TrackerTag")),
376  fShowerTag(p.get<std::string>("ShowerTag")),
377  fPFParticleTag(p.get<std::string>("PFParticleTag")),
378  fGeneratorTag(p.get<std::string>("GeneratorTag")),
379  fVerbose(p.get<bool>("Verbose"))
380 {
381 
382 }
383 
384 void protoana::protonmccnn::beginJob(){
385 
386  art::ServiceHandle<art::TFileService> tfs;
387 
388  fPandoraBeam = tfs->make<TTree>("PandoraBeam", "Beam events reconstructed with Pandora");
389  fPandoraBeam->Branch("run", &fRun, "run/I");
390  fPandoraBeam->Branch("subrun", &fSubRun, "subrun/I");
391  fPandoraBeam->Branch("event", &fevent, "event/I");
392  fPandoraBeam->Branch("timestamp", &fTimeStamp, "timestamp/D");
393  fPandoraBeam->Branch("Nactivefembs", &fNactivefembs, "Nactivefembs[5]/I");
394 
395  fPandoraBeam->Branch("beamtrigger", &fbeamtrigger, "beamtrigger/I");
396  fPandoraBeam->Branch("tof", &ftof, "tof/D");
397  fPandoraBeam->Branch("cerenkov1", &fcerenkov1, "cerenkov1/I");
398  fPandoraBeam->Branch("cerenkov2", &fcerenkov2, "cerenkov2/I");
399  fPandoraBeam->Branch("beamtrackMomentum", &fbeamtrackMomentum, "beamtrackMomentum/D");
400  fPandoraBeam->Branch("beamtrackP", &fbeamtrackP, "beamtrackP[3]/D");
401  fPandoraBeam->Branch("beamtrackEnergy", &fbeamtrackEnergy, "beamtrackEnergy/D");
402  fPandoraBeam->Branch("beamtrackPos", &fbeamtrackPos, "beamtrackPos[3]/D");
403  fPandoraBeam->Branch("beamtrackDir", &fbeamtrackDir, "beamtrackDir[3]/D");
404  fPandoraBeam->Branch("beamtrackTime", &fbeamtrackTime, "beamtrackTime/D");
405  fPandoraBeam->Branch("beamtrackPdg", &fbeamtrackPdg, "beamtrackPdg/I");
406  fPandoraBeam->Branch("beamtrackID", &fbeamtrackID, "beamtrackID/I");
407  //fPandoraBeam->Branch("NumberBeamTrajectoryPoints", &NumberBeamTrajectoryPoints, "NumberBeamTrajectoryPoints/I");
408  fPandoraBeam->Branch("beamtrk_x",&beamtrk_x);
409  fPandoraBeam->Branch("beamtrk_y",&beamtrk_y);
410  fPandoraBeam->Branch("beamtrk_z",&beamtrk_z);
411  fPandoraBeam->Branch("beamtrk_Px",&beamtrk_Px);
412  fPandoraBeam->Branch("beamtrk_Py",&beamtrk_Py);
413  fPandoraBeam->Branch("beamtrk_Pz",&beamtrk_Pz);
414  fPandoraBeam->Branch("beamtrk_Eng",&beamtrk_Eng);
415 
416  fPandoraBeam->Branch("ke_ff", &ke_ff, "ke_ff/D");
417  fPandoraBeam->Branch("pos_ff",&pos_ff);
418 
419  fPandoraBeam->Branch("primtrk_ke_true",&primtrk_ke_true);
420  fPandoraBeam->Branch("primtrk_ke_reco",&primtrk_ke_reco);
421  fPandoraBeam->Branch("primtrk_end_g4process",&primtrk_end_g4process);
422  fPandoraBeam->Branch("primtrk_range_true",&primtrk_range_true);
423 
424  fPandoraBeam->Branch("ke_ff_true", &ke_ff_true, "ke_ff_true/D");
425  //fPandoraBeam->Branch("primtrk_hit_processkey",&primtrk_hit_processkey);
426  fPandoraBeam->Branch("primtrk_hitx_true",&primtrk_hitx_true);
427  fPandoraBeam->Branch("primtrk_hity_true",&primtrk_hity_true);
428  fPandoraBeam->Branch("primtrk_hitz_true",&primtrk_hitz_true);
429  fPandoraBeam->Branch("primtrk_trkid_true",&primtrk_trkid_true);
430  fPandoraBeam->Branch("primtrk_edept_true",&primtrk_edept_true);
431 
432  fPandoraBeam->Branch("pfphit_peaktime_c",&pfphit_peaktime_c);
433  fPandoraBeam->Branch("pfphit_peaktime_u",&pfphit_peaktime_u);
434  fPandoraBeam->Branch("pfphit_peaktime_v",&pfphit_peaktime_v);
435  fPandoraBeam->Branch("pfphit_wireid_c",&pfphit_wireid_c);
436  fPandoraBeam->Branch("pfphit_wireid_u",&pfphit_wireid_u);
437  fPandoraBeam->Branch("pfphit_wireid_v",&pfphit_wireid_v);
438 
439  fPandoraBeam->Branch("primtrk_true_x",&primtrk_true_x);
440  fPandoraBeam->Branch("primtrk_true_y",&primtrk_true_y);
441  fPandoraBeam->Branch("primtrk_true_z",&primtrk_true_z);
442  fPandoraBeam->Branch("primtrk_true_trkid",&primtrk_true_trkid);
443  fPandoraBeam->Branch("primtrk_true_edept",&primtrk_true_edept);
444 
445  fPandoraBeam->Branch("primtrj_true_x",&primtrj_true_x);
446  fPandoraBeam->Branch("primtrj_true_y",&primtrj_true_y);
447  fPandoraBeam->Branch("primtrj_true_z",&primtrj_true_z);
448  fPandoraBeam->Branch("primtrj_true_edept",&primtrj_true_edept);
449 
450  fPandoraBeam->Branch("vertex", &fvertex, "vertex[3]/D");
451  fPandoraBeam->Branch("secvertex", &fsecvertex, "secvertex[3]/D");
452  fPandoraBeam->Branch("isprimarytrack", &fisprimarytrack, "isprimarytrack/I");
453  fPandoraBeam->Branch("isprimaryshower", &fisprimaryshower, "isprimaryshower/I");
454  fPandoraBeam->Branch("primaryBDTScore", &fprimaryBDTScore, "primaryBDTScore/D");
455  fPandoraBeam->Branch("primaryNHits", &fprimaryNHits, "fprimaryNHits/I");
456  fPandoraBeam->Branch("primaryTheta", &fprimaryTheta, "primaryTheta/D");
457  fPandoraBeam->Branch("primaryPhi", &fprimaryPhi, "primaryPhi/D");
458  fPandoraBeam->Branch("primaryLength", &fprimaryLength, "primaryLength/D");
459  fPandoraBeam->Branch("primaryMomentum", &fprimaryMomentum, "primaryMomentum/D");
460  fPandoraBeam->Branch("primaryEndMomentum", &fprimaryEndMomentum, "primaryEndMomentum/D");
461  fPandoraBeam->Branch("primaryEndPosition", &fprimaryEndPosition, "primaryEndPosition[3]/D");
462  fPandoraBeam->Branch("primaryStartPosition", &fprimaryStartPosition, "primaryStartPosition[3]/D");
463  fPandoraBeam->Branch("primaryEndDirection", &fprimaryEndDirection, "primaryEndDirection[3]/D");
464  fPandoraBeam->Branch("primaryStartDirection", &fprimaryStartDirection, "primaryStartDirection[3]/D");
465  fPandoraBeam->Branch("primaryOpeningAngle", &fprimaryOpeningAngle, "primaryOpeningAngle/D");
466  fPandoraBeam->Branch("primaryID", &fprimaryID, "primaryID/I");
467  fPandoraBeam->Branch("primaryShowerBestPlane", &fprimaryShowerBestPlane, "primaryShowerBestPlane/I");
468  fPandoraBeam->Branch("primaryShowerEnergy", &fprimaryShowerEnergy, "primaryShowerEnergy/I");
469  fPandoraBeam->Branch("primaryShowerMIPEnergy", &fprimaryShowerMIPEnergy, "primaryShowerMIPEnergy/I");
470  fPandoraBeam->Branch("primaryShowerdEdx", &fprimaryShowerdEdx, "primaryShowerdEdx/I");
471  fPandoraBeam->Branch("primaryMomentumByRangeProton", &fprimaryMomentumByRangeProton, "primaryMomentumByRangeProton/D");
472  fPandoraBeam->Branch("primaryMomentumByRangeMuon", &fprimaryMomentumByRangeMuon, "primaryMomentumByRangeMuon/D");
473  fPandoraBeam->Branch("primaryKineticEnergy", &fprimaryKineticEnergy, "primaryKineticEnergy[3]/D");
474  fPandoraBeam->Branch("primaryRange", &fprimaryRange, "primaryRange[3]/D");
475  fPandoraBeam->Branch("primaryT0", &fprimaryT0, "primaryT0/D");
476 
477  fPandoraBeam->Branch("primary_truth_TrackId", &fprimary_truth_TrackId, "primary_truth_TrackId/I");
478  fPandoraBeam->Branch("primary_truth_Pdg", &fprimary_truth_Pdg, "primary_truth_Pdg/I");
479  fPandoraBeam->Branch("truthpdg", &ftruthpdg, "truthpdg/I");
480  fPandoraBeam->Branch("primary_truth_StartPosition", &fprimary_truth_StartPosition, "primary_truth_StartPosition[4]/D");
481  fPandoraBeam->Branch("primary_truth_StartPosition_MC", &fprimary_truth_StartPosition_MC, "primary_truth_StartPosition_MC[4]/D");
482  fPandoraBeam->Branch("primary_truth_EndPosition", &fprimary_truth_EndPosition, "primary_truth_EndPosition[4]/D");
483  fPandoraBeam->Branch("primary_truth_EndPosition_MC", &fprimary_truth_EndPosition_MC, "primary_truth_EndPosition_MC[4]/D");
484  fPandoraBeam->Branch("primary_truth_Momentum", &fprimary_truth_Momentum, "primary_truth_Momentum[4]/D");
485  fPandoraBeam->Branch("primary_truth_EndMomentum", &fprimary_truth_EndMomentum, "primary_truth_EndMomentum[4]/D");
486  fPandoraBeam->Branch("primary_truth_P", &fprimary_truth_P, "primary_truth_P/D");
487  fPandoraBeam->Branch("primary_truth_Pt", &fprimary_truth_Pt, "primary_truth_Pt/D");
488  fPandoraBeam->Branch("primary_truth_Mass", &fprimary_truth_Mass, "primary_truth_Mass/D");
489  fPandoraBeam->Branch("primary_truth_Theta", &fprimary_truth_Theta, "primary_truth_Theta/D");
490  fPandoraBeam->Branch("primary_truth_Phi", &fprimary_truth_Phi, "primary_truth_Phi/D");
491  //fPandoraBeam->Branch("primary_truth_Process", &fprimary_truth_Process, "primary_truth_Process/I");
492  //fPandoraBeam->Branch("primary_truth_Process", &fprimary_truth_Process);
493  fPandoraBeam->Branch("primary_truth_EndProcess", &fprimary_truth_EndProcess);
494  //fPandoraBeam->Branch("primary_backtrker_truth_Process", &fprimary_backtrker_truth_Process);
495  //fPandoraBeam->Branch("primary_backtrker_truth_EndProcess", &fprimary_backtrker_truth_EndProcess);
496  fPandoraBeam->Branch("primary_truth_Isbeammatched", &fprimary_truth_Isbeammatched, "primary_truth_Isbeammatched/I");
497  fPandoraBeam->Branch("primary_truth_NDaughters", &fprimary_truth_NDaughters, "primary_truth_NDaughters/I");
498  //fPandoraBeam->Branch("G4Process_Primary_End", &G4Process_Primary_End);
499 
500  fPandoraBeam->Branch("interactionX",&interactionX);
501  fPandoraBeam->Branch("interactionY",&interactionY);
502  fPandoraBeam->Branch("interactionZ",&interactionZ);
503  fPandoraBeam->Branch("interactionE",&interactionE);
504  fPandoraBeam->Branch("interactionProcesslist",&interactionProcesslist);
505  fPandoraBeam->Branch("interactionAngles",&interactionAngles);
506 
507  fPandoraBeam->Branch("Zintersection",&Zintersection);
508  fPandoraBeam->Branch("Yintersection",&Yintersection);
509  fPandoraBeam->Branch("Xintersection",&Xintersection);
510  fPandoraBeam->Branch("timeintersection",&timeintersection);
511 
512  fPandoraBeam->Branch("interaction_wid_c",&interaction_wid_c);
513  fPandoraBeam->Branch("interaction_tt_c",&interaction_tt_c);
514  fPandoraBeam->Branch("interaction_wid_v",&interaction_wid_v);
515  fPandoraBeam->Branch("interaction_tt_v",&interaction_tt_v);
516  fPandoraBeam->Branch("interaction_wid_u",&interaction_wid_u);
517  fPandoraBeam->Branch("interaction_tt_u",&interaction_tt_u);
518 
519  fPandoraBeam->Branch("inelscore_c",&inelscore_c);
520  fPandoraBeam->Branch("elscore_c",&elscore_c);
521  fPandoraBeam->Branch("nonescore_c",&nonescore_c);
522  fPandoraBeam->Branch("x_c",&x_c);
523  fPandoraBeam->Branch("y_c",&y_c);
524  fPandoraBeam->Branch("z_c",&z_c);
525  //fPandoraBeam->Branch("xyz_inelscore_c",&xyz_inelscore_c,"xyz_inelscore_c[3]/D");
526  //fPandoraBeam->Branch("xyz_elscore_c",&xyz_elscore_c,"xyz_elscore_c[3]/D");
527 
528  fPandoraBeam->Branch("wid_c",&wid_c);
529  fPandoraBeam->Branch("tt_c",&tt_c);
530  fPandoraBeam->Branch("wid_v",&wid_v);
531  fPandoraBeam->Branch("tt_v",&tt_v);
532  fPandoraBeam->Branch("wid_u",&wid_u);
533  fPandoraBeam->Branch("tt_u",&tt_u);
534 
535  fPandoraBeam->Branch("NDAUGHTERS", &fNDAUGHTERS, "NDAUGHTERS/I");
536  fPandoraBeam->Branch("isdaughtertrack", &fisdaughtertrack, "isdaughtertrack[NDAUGHTERS]/I");
537  fPandoraBeam->Branch("isdaughtershower", &fisdaughtershower, "isdaughtershower[NDAUGHTERS]/I");
538  fPandoraBeam->Branch("daughterNHits", &fdaughterNHits, "daughterNHits[NDAUGHTERS]/I");
539  fPandoraBeam->Branch("daughterTheta", &fdaughterTheta, "daughterTheta[NDAUGHTERS]/D");
540  fPandoraBeam->Branch("daughterPhi", &fdaughterPhi, "daughterPhi[NDAUGHTERS]/D");
541  fPandoraBeam->Branch("daughterLength", &fdaughterLength, "daughterLength[NDAUGHTERS]/D");
542  fPandoraBeam->Branch("daughterMomentum", &fdaughterMomentum, "daughterMomentum[NDAUGHTERS]/D");
543  fPandoraBeam->Branch("daughterEndMomentum", &fdaughterEndMomentum, "daughterEndMomentum[NDAUGHTERS]/D");
544  fPandoraBeam->Branch("daughterEndPosition", &fdaughterEndPosition, "daughterEndPosition[NDAUGHTERS][3]/D");
545  fPandoraBeam->Branch("daughterStartPosition", &fdaughterStartPosition, "daughterStartPosition[NDAUGHTERS][3]/D");
546  fPandoraBeam->Branch("daughterStartDirection", &fdaughterStartDirection, "daughterStartDirection[NDAUGHTERS][3]/D");
547  fPandoraBeam->Branch("daughterEndDirection", &fdaughterEndDirection, "daughterEndDirection[NDAUGHTERS][3]/D");
548  fPandoraBeam->Branch("daughterOpeningAngle", &fdaughterOpeningAngle, "daughterOpeningAngle[NDAUGHTERS]/D");
549  fPandoraBeam->Branch("daughterShowerBestPlane", &fdaughterShowerBestPlane, "daughterShowerBestPlane[NDAUGHTERS]/D");
550  fPandoraBeam->Branch("daughterShowerEnergy", &fdaughterShowerEnergy, "daughterShowerEnergy[NDAUGHTERS]/D");
551  fPandoraBeam->Branch("daughterShowerMIPEnergy", &fdaughterShowerMIPEnergy, "daughterShowerMIPEnergy[NDAUGHTERS]/D");
552  fPandoraBeam->Branch("daughterShowerdEdx", &fdaughterShowerdEdx, "daughterShowerdEdx[NDAUGHTERS]/D");
553  fPandoraBeam->Branch("daughterMomentumByRangeProton", &fdaughterMomentumByRangeProton,"daughterMomentumByRangeProton[NDAUGHTERS]/D");
554  fPandoraBeam->Branch("daughterMomentumByRangeMuon", &fdaughterMomentumByRangeMuon, "daughterMomentumByRangeMuon[NDAUGHTERS]/D");
555  fPandoraBeam->Branch("daughterKineticEnergy", &fdaughterKineticEnergy, "daughterKineticEnergy[NDAUGHTERS][3]/D");
556  fPandoraBeam->Branch("daughterRange", &fdaughterRange, "daughterRange[NDAUGHTERS][3]/D");
557  fPandoraBeam->Branch("daughterID", &fdaughterID, "daughterID[NDAUGHTERS]/I");
558  //fPandoraBeam->Branch("daughterT0", &fdaughterT0, "daughterT0[NDAUGHTERS]/D");
559 
560  fPandoraBeam->Branch("daughter_truth_TrackId", &fdaughter_truth_TrackId, "daughter_truth_TrackId[NDAUGHTERS]/I");
561  fPandoraBeam->Branch("daughter_truth_Pdg", &fdaughter_truth_Pdg, "daughter_truth_Pdg[NDAUGHTERS]/I");
562  fPandoraBeam->Branch("daughter_truth_StartPosition", &fdaughter_truth_StartPosition, "daughter_truth_StartPosition[NDAUGHTERS][4]/D");
563  fPandoraBeam->Branch("daughter_truth_EndPosition", &fdaughter_truth_EndPosition, "daughter_truth_EndPosition[NDAUGHTERS][4]/D");
564  fPandoraBeam->Branch("daughter_truth_Momentum", &fdaughter_truth_Momentum, "daughter_truth_Momentum[NDAUGHTERS][4]/D");
565  fPandoraBeam->Branch("daughter_truth_EndMomentum", &fdaughter_truth_EndMomentum, "daughter_truth_EndMomentum[NDAUGHTERS][4]/D");
566  fPandoraBeam->Branch("daughter_truth_P", &fdaughter_truth_P, "daughter_truth_P[NDAUGHTERS]/D");
567  fPandoraBeam->Branch("daughter_truth_Pt", &fdaughter_truth_Pt, "daughter_truth_Pt[NDAUGHTERS]/D");
568  fPandoraBeam->Branch("daughter_truth_Mass", &fdaughter_truth_Mass, "daughter_truth_Mass[NDAUGHTERS]/D");
569  fPandoraBeam->Branch("daughter_truth_Theta", &fdaughter_truth_Theta, "daughter_truth_Theta[NDAUGHTERS]/D");
570  fPandoraBeam->Branch("daughter_truth_Phi", &fdaughter_truth_Phi, "daughter_truth_Phi[NDAUGHTERS]/D");
571  fPandoraBeam->Branch("daughter_truth_Process", &fdaughter_truth_Process, "daughter_truth_Process[NDAUGHTERS]/I");
572 
573  fPandoraBeam->Branch("primtrk_dqdx",&primtrk_dqdx);
574  fPandoraBeam->Branch("primtrk_dedx",&primtrk_dedx);
575  fPandoraBeam->Branch("primtrk_resrange",&primtrk_resrange);
576  fPandoraBeam->Branch("primtrk_range",&primtrk_range);
577  fPandoraBeam->Branch("primtrk_hitx",&primtrk_hitx);
578  fPandoraBeam->Branch("primtrk_hity",&primtrk_hity);
579  fPandoraBeam->Branch("primtrk_hitz",&primtrk_hitz);
580  fPandoraBeam->Branch("primtrk_pitch",&primtrk_pitch);
581 
582  //fPandoraCosmics = tfs->make<TTree>("PandoraCosmics", "Cosmic tracks reconstructed with Pandora");
583  //fPandoraCosmics->Branch("run", &fRun, "run/I");
584  //fPandoraCosmics->Branch("subrun", &fSubRun, "subrun/I");
585  //fPandoraCosmics->Branch("event", &fevent, "event/I");
586  //fPandoraCosmics->Branch("timestamp", &fTimeStamp, "timestamp/D");
587  //fPandoraCosmics->Branch("Nactivefembs", &fNactivefembs, "Nactivefembs[5]/I");
588  //fPandoraCosmics->Branch("beamtrigger", &fbeamtrigger, "beamtrigger/I");
589  //fPandoraCosmics->Branch("tof", &ftof, "tof/D");
590  //fPandoraCosmics->Branch("cerenkov1", &fcerenkov1, "cerenkov1/I");
591  //fPandoraCosmics->Branch("cerenkov2", &fcerenkov2, "cerenkov2/I");
592  //fPandoraCosmics->Branch("beamtrackMomentum", &fbeamtrackMomentum, "beamtrackMomentum/D");
593  //fPandoraCosmics->Branch("beamtrackPos", &fbeamtrackPos, "beamtrackPos[3]/D");
594  //fPandoraCosmics->Branch("beamtrackDir", &fbeamtrackDir, "beamtrackDir[3]/D");
595 
596 }
597 
598 void protoana::protonmccnn::analyze(art::Event const & evt){
599 
600  // Initialise tree parameters
601  Initialise();
602 
603  int beamid=-9999;
604  int truthid=-999;
605 
606  art::ServiceHandle<cheat::ParticleInventoryService> pi_serv;
607  fGeometry = &*(art::ServiceHandle<geo::Geometry>());
608  anab::MVAReader<recob::Hit,3> hitResults(evt, fNNetModuleLabel);
609 
610  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService>()->DataFor(evt);
611  auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService>()->DataFor(evt, clockData);
612  //T0
613  std::vector<const anab::T0*> T0s;
614  nT0s = T0s.size();
615  std::cout << "Got " << nT0s << " T0s" <<std::endl;
616 
617  if(nT0s > 0){
618  std::cout << "T0s size: " << nT0s << std::endl;
619  MCTruthT0 = T0s[0]->Time();
620  }
621  else{
622  std::cout << "No T0s found" << std::endl;
623  MCTruthT0 = 0;
624  }
625  TickT0 = MCTruthT0 / sampling_rate(clockData);
626  std::cout<<"TickT0:"<<TickT0<<std::endl;
627 
628 
629  fRun = evt.run();
630  fSubRun = evt.subRun();
631  fevent = evt.id().event();
632  art::Timestamp ts = evt.time();
633  if (ts.timeHigh() == 0){
634  TTimeStamp ts2(ts.timeLow());
635  fTimeStamp = ts2.AsDouble();
636  }
637  else{
638  TTimeStamp ts2(ts.timeHigh(), ts.timeLow());
639  fTimeStamp = ts2.AsDouble();
640  }
641 
642  // Get number of active fembs
643  if(!evt.isRealData()){
644  for(int k=0; k < 6; k++)
645  fNactivefembs[k] = 20;
646  }
647  else{
648  for(int k=0; k < 6; k++)
649  fNactivefembs[k] = dataUtil.GetNActiveFembsForAPA(evt, k);
650  }
651 
652  bool beamTriggerEvent = false;
653  // If this event is MC then we can check what the true beam particle is
654  auto mcTruths = evt.getValidHandle<std::vector<simb::MCTruth>>(fGeneratorTag);
655  if(!evt.isRealData()){
656  // Firstly we need to get the list of MCTruth objects from the generator. The standard protoDUNE
657  // simulation has fGeneratorTag = "generator"
658  auto mcTruths = evt.getValidHandle<std::vector<simb::MCTruth>>(fGeneratorTag);
659 
660  // Also get the reconstructed beam information in the MC - TO DO
661  //auto beamsim = evt.getValidHandle<std::vector<sim::ProtoDUNEbeamsim> >(fGeneratorTag);
662  //const sim::ProtoDUNEbeamsim beamsimobj = (*beamsim)[0];
663  //std::cout << beamsimobj.NInstruments() << std::endl;
664  //sim::ProtoDUNEBeamInstrument beamsim_tof1 = beamsimobj.GetInstrument("TOF1");
665  //sim::ProtoDUNEBeamInstrument beamsim_trig2 = beamsimobj.GetInstrument("TRIG2");
666  //std::cout << beamsim_trig2.GetT() - beamsim_tof1.GetT() << " , " << beamsim_trig2.GetSmearedVar1() - beamsim_tof1.GetSmearedVar1() << std::endl;
667  //std::cout << beamsimobj.GetInstrument("TRIG2").GetT() - beamsimobj.GetInstrument("TOF1").GetT() << " , " << beamsimobj.GetInstrument("TRIG2").GetSmearedVar1() - beamsimobj.GetInstrument("TOF1").GetSmearedVar1() << std::endl;
668 
669  // mcTruths is basically a pointer to an std::vector of simb::MCTruth objects. There should only be one
670  // of these, so we pass the first element into the function to get the good particle
671  const simb::MCParticle* geantGoodParticle = truthUtil.GetGeantGoodParticle((*mcTruths)[0],evt);
672 
673  if(geantGoodParticle != 0x0){
674  std::cout << "Found GEANT particle corresponding to the good particle with pdg = " << geantGoodParticle->PdgCode()
675  << " , track id = " << geantGoodParticle->TrackId()
676  << " , Vx/Vy/Vz = " << geantGoodParticle->Vx() << "/"<< geantGoodParticle->Vy() << "/" << geantGoodParticle->Vz()
677  << std::endl;
678 
679  //NumberBeamTrajectoryPoints=geantGoodParticle->NumberTrajectoryPoints();
680  //std::cout<<"NumberBeamTrajectoryPoints for beam particles:"<<NumberBeamTrajectoryPoints<<std::endl;
681 
682  beamTriggerEvent = true;
683  fbeamtrigger = 12;
684  fbeamtrackPos[0] = geantGoodParticle->Vx();
685  fbeamtrackPos[1] = geantGoodParticle->Vy();
686  fbeamtrackPos[2] = geantGoodParticle->Vz();
687  fbeamtrackMomentum = geantGoodParticle->P();
688  fbeamtrackP[0] = geantGoodParticle->Px();
689  fbeamtrackP[1] = geantGoodParticle->Py();
690  fbeamtrackP[2] = geantGoodParticle->Pz();
691  fbeamtrackEnergy = geantGoodParticle->E();
692  fbeamtrackPdg = geantGoodParticle->PdgCode();
693  fbeamtrackTime = geantGoodParticle->T();
694  fbeamtrackID = geantGoodParticle->TrackId();
695  //prim_energy=0;
696 
697  for(size_t i_s=0; i_s < geantGoodParticle->NumberTrajectoryPoints(); i_s++){ //loop over beam tracks
698  beamtrk_x.push_back(geantGoodParticle->Position(i_s).X());
699  beamtrk_y.push_back(geantGoodParticle->Position(i_s).Y());
700  beamtrk_z.push_back(geantGoodParticle->Position(i_s).Z());
701 
702  beamtrk_Px.push_back(geantGoodParticle->Momentum(i_s).X());
703  beamtrk_Py.push_back(geantGoodParticle->Momentum(i_s).Y());
704  beamtrk_Pz.push_back(geantGoodParticle->Momentum(i_s).Z());
705 
706  beamtrk_Eng.push_back(geantGoodParticle->Momentum(i_s).E()-geantGoodParticle->Mass());
707  } //loop over beam trks
708 
709  //Get KE at front face of TPC --------------------------------------------//
710  int key_reach_tpc=-99;
711  for (size_t kk=0; kk<beamtrk_z.size(); ++kk) { //loop over all beam hits
712  double zpos_beam=beamtrk_z.at(kk);
713  if ((zpos_beam+0.49375)<0.01&&(zpos_beam+0.49375)>-0.01) { //find key at ff
714  key_reach_tpc=kk;
715  } //kind key at ff
716  } //loop over all beam hits
717 
718  //Define KE and Momentum at the entering point of TPC
719  if (key_reach_tpc!=-99) { //ke and pos at front face
720  ke_ff=1000.*(beamtrk_Eng.at(key_reach_tpc)); //unit: MeV
721  pos_ff.push_back(beamtrk_x.at(key_reach_tpc));
722  pos_ff.push_back(beamtrk_y.at(key_reach_tpc));
723  pos_ff.push_back(beamtrk_z.at(key_reach_tpc));
724 
725  fprimary_truth_Momentum[0]=beamtrk_Px.at(key_reach_tpc);
726  fprimary_truth_Momentum[1]=beamtrk_Py.at(key_reach_tpc);
727  fprimary_truth_Momentum[2]=beamtrk_Pz.at(key_reach_tpc);
728  } //ke and pos at front face
729  else {
730  std::cout<<"This particle doesn't enter TPC!!!"<<std::endl;
731  }
732  //Get KE at front face of TPC --------------------------------------------//
733 
734  //Get Truth info
735  fprimary_truth_TrackId = geantGoodParticle->TrackId();
736  fprimary_truth_Pdg = geantGoodParticle->PdgCode();
737  beamid = geantGoodParticle->TrackId();
738  fprimary_truth_StartPosition[3] = geantGoodParticle->T();
739  fprimary_truth_EndPosition[3] = geantGoodParticle->EndT();
740 
741  fprimary_truth_StartPosition_MC[0] = geantGoodParticle->Vx();
742  fprimary_truth_StartPosition_MC[1] = geantGoodParticle->Vy();
743  fprimary_truth_StartPosition_MC[2] = geantGoodParticle->Vz();
744  fprimary_truth_StartPosition_MC[3] = geantGoodParticle->T();
745 
746  fprimary_truth_EndPosition_MC[0] = geantGoodParticle->EndX();
747  fprimary_truth_EndPosition_MC[1] = geantGoodParticle->EndY();
748  fprimary_truth_EndPosition_MC[2] = geantGoodParticle->EndZ();
749  fprimary_truth_EndPosition_MC[3] = geantGoodParticle->EndT();
750 
751  fprimary_truth_P = geantGoodParticle->P();
752  //fprimary_truth_Momentum[0] = geantGoodParticle->Px();
753  //fprimary_truth_Momentum[1] = geantGoodParticle->Py();
754  //fprimary_truth_Momentum[2] = geantGoodParticle->Pz();
755  fprimary_truth_Momentum[3] = geantGoodParticle->E();
756  fprimary_truth_Pt = geantGoodParticle->Pt();
757  fprimary_truth_Mass = geantGoodParticle->Mass();
758  //fprimary_truth_EndMomentum[0] = geantGoodParticle->EndPx();
759  //fprimary_truth_EndMomentum[1] = geantGoodParticle->EndPy();
760  //fprimary_truth_EndMomentum[2] = geantGoodParticle->EndPz();
761  fprimary_truth_EndMomentum[3] = geantGoodParticle->EndE();
762  fprimary_truth_Theta = geantGoodParticle->Momentum().Theta();
763  fprimary_truth_Phi = geantGoodParticle->Momentum().Phi();
764  fprimary_truth_NDaughters = geantGoodParticle->NumberDaughters();
765  //fprimary_truth_Process = int(geantGoodParticle->Trajectory().ProcessToKey(geantGoodParticle->Process()));
766  //fprimary_truth_Process = geantGoodParticle->Process(); //HY::wired result
767  //fprimary_truth_EndProcess = geantGoodParticle->EndProcess(); //HY:wierd result
768 
769  //fprimary_backtrker_truth_Process = int(geantGoodParticle->Trajectory().ProcessToKey(geantGoodParticle->Process()));
770 
771  }
772  }
773  else{
774  // For data we can see if this event comes from a beam trigger
775  beamTriggerEvent = dataUtil.IsBeamTrigger(evt);
776 
777  std::vector< art::Ptr<beam::ProtoDUNEBeamEvent> > beaminfo;
778  auto pdbeamHandle = evt.getHandle< std::vector<beam::ProtoDUNEBeamEvent> >(fBeamModuleLabel);
779  if (pdbeamHandle) art::fill_ptr_vector(beaminfo, pdbeamHandle);
780 
781  for(unsigned int i = 0; i < beaminfo.size(); ++i){
782  //if(!beaminfo[i]->CheckIsMatched()) continue;
783  fbeamtrigger = beaminfo[i]->GetTimingTrigger();
784  fbeamtrackTime = (double)beaminfo[i]->GetRDTimestamp();
785 
786  // If ToF is 0-3 there was a good match corresponding to the different pair-wise combinations of the upstream and downstream channels
787  if(beaminfo[i]->GetTOFChan() >= 0)
788  ftof = beaminfo[i]->GetTOF();
789 
790  // Get Cerenkov
791  if(beaminfo[i]->GetBITrigger() == 1){
792  fcerenkov1 = beaminfo[i]->GetCKov0Status();
793  fcerenkov2 = beaminfo[i]->GetCKov1Status();
794  }
795 
796  // Beam particle could have more than one tracks - for now take the first one, need to do this properly
797  auto & tracks = beaminfo[i]->GetBeamTracks();
798  if(!tracks.empty()){
799  fbeamtrackPos[0] = tracks[0].End().X();
800  fbeamtrackPos[1] = tracks[0].End().Y();
801  fbeamtrackPos[2] = tracks[0].End().Z();
802  fbeamtrackDir[0] = tracks[0].StartDirection().X();
803  fbeamtrackDir[1] = tracks[0].StartDirection().Y();
804  fbeamtrackDir[2] = tracks[0].StartDirection().Z();
805  }
806 
807  // Beam momentum
808  auto & beammom = beaminfo[i]->GetRecoBeamMomenta();
809  if(!beammom.empty())
810  fbeamtrackMomentum = beammom[0];
811 
812  // For now only take the first beam particle - need to add some criteria if more than one are found
813  break;
814 
815  }
816  }
817 
818  /*
819  // Now we want to access the output from Pandora. This comes in the form of particle flow objects (recob::PFParticle).
820  // The primary PFParticles are those we want to consider and these PFParticles then have a hierarchy of daughters that
821  // describe the whole interaction of a given primary particle
822  //
823  // / daughter track
824  // /
825  // primary track /
826  // ---------------- ---- daughter track
827  // \
828  // /\-
829  // /\\-- daughter shower
830  //
831  // The above primary PFParticle will have links to three daughter particles, two track-like and one shower-like
832  */
833 
834  // Track momentum algorithm calculates momentum based on track range
835  trkf::TrackMomentumCalculator trmom;
836  //trmom.SetMinLength(100);
837 
838  // Get all of the PFParticles, by default from the "pandora" product
839  auto recoParticles = evt.getValidHandle<std::vector<recob::PFParticle>>(fPFParticleTag);
840  std::cout << "All primary pfParticles = " << pfpUtil.GetNumberPrimaryPFParticle(evt,fPFParticleTag) << std::endl;
841 
842  // We'd like to find the beam particle. Pandora tries to do this for us, so let's use the PFParticle utility
843  // to look for it. Pandora reconstructs slices containing one (or sometimes more) primary PFParticles. These
844  // are tagged as either beam or cosmic for ProtoDUNE. This function automatically considers only those
845  // PFParticles considered as primary
846  //std::vector<const recob::PFParticle*> pfParticles = pfpUtil.GetPFParticlesFromBeamSlice(evt,fPFParticleTag);
847  auto pfParticles = pfpUtil.GetPFParticlesFromBeamSlice(evt,fPFParticleTag);
848 
849  //cluster information
850  std::vector<art::Ptr<recob::Track> > tracklist;
851  std::vector<art::Ptr<recob::PFParticle> > pfplist;
852 
853  auto trackListHandle = evt.getHandle< std::vector<recob::Track> >("pandoraTrack");
854  if (trackListHandle) art::fill_ptr_vector(tracklist, trackListHandle);
855 
856  auto PFPListHandle = evt.getHandle< std::vector<recob::PFParticle> > ("pandora");
857  if (PFPListHandle) art::fill_ptr_vector(pfplist, PFPListHandle);
858 
859  std::vector<art::Ptr<recob::Cluster>> clusterlist;
860  auto clusterListHandle = evt.getHandle< std::vector<recob::Cluster> >("pandora"); ; // to get information about the hits
861  if(clusterListHandle) art::fill_ptr_vector(clusterlist, clusterListHandle);
862 
863  art::FindManyP<recob::Cluster> fmcp(PFPListHandle,evt,"pandora");
864  art::FindManyP<recob::Track> pftrack(PFPListHandle,evt,"pandoraTrack");
865 
866  std::cout<<"number of pfp_particles "<<pfplist.size()<<std::endl;
867  std::cout<<" size of pfParticles "<<pfParticles.size()<<std::endl;
868  art::FindManyP<recob::Hit, recob::TrackHitMeta> fmthm(trackListHandle, evt,"pandoraTrack"); // to associate tracks and hits
869 
870  /* for(size_t p1=0;p1<pfplist.size();p1++){
871  std::vector<art::Ptr<recob::Track>> trk=pftrack.at(p1);
872  if(trk.size()) std::cout<<" trk key "<<trk[0].key()<<std::endl;
873 
874  std::vector<art::Ptr<recob::Cluster>> allClusters=fmcp.at(p1);
875  std::cout<<"cluster size for each particle "<<allClusters.size()<<std::endl;
876  for(size_t c1=0;c1<allClusters.size();c1++){
877  std::cout<<" cluster ID "<<allClusters[c1]->ID();
878  std::cout<<" plane number "<<allClusters[c1]->Plane().Plane;
879  std::cout<<" TPC number "<<allClusters[c1]->Plane().TPC;
880  std::cout<<" start wire "<<allClusters[c1]->StartWire();
881  std::cout<<" end wire "<<allClusters[c1]->EndWire();
882  std::cout<<" start tick "<<allClusters[c1]->StartTick();
883  std::cout<<" end tick "<<allClusters[c1]->EndTick();
884  }
885  }*/
886 
887 
888 
889 
890  // We can now look at these particles
891  for(const recob::PFParticle* particle : pfParticles){
892 
893  // Pandora's BDT beam-cosmic score
894  fprimaryBDTScore = (double)pfpUtil.GetBeamCosmicScore(*particle,evt,fPFParticleTag);
895 
896  // NHits associated with this pfParticle
897  fprimaryNHits = (pfpUtil.GetPFParticleHits(*particle,evt,fPFParticleTag)).size();
898 
899  // Get the T0 for this pfParticle
900  std::vector<anab::T0> pfT0vec = pfpUtil.GetPFParticleT0(*particle,evt,fPFParticleTag);
901  if(!pfT0vec.empty())
902  fprimaryT0 = pfT0vec[0].Time();
903 
904  //std::cout << "Pdg Code = " << particle->PdgCode() << std::endl;
905  // "particle" is the pointer to our beam particle. The recob::Track or recob::Shower object
906  // of this particle might be more helpful. These return null pointers if not track-like / shower-like
907  const recob::Track* thisTrack = pfpUtil.GetPFParticleTrack(*particle,evt,fPFParticleTag,fTrackerTag);
908  const recob::Shower* thisShower = pfpUtil.GetPFParticleShower(*particle,evt,fPFParticleTag,fShowerTag);
909  if(thisTrack != 0x0) { //this track
910  // Get the true mc particle
911  const simb::MCParticle* mcparticle0 = truthUtil.GetMCParticleFromRecoTrack(clockData, *thisTrack, evt, fTrackerTag);
912  std::cout<<"inside the this track loop "<<std::endl;
913  if(mcparticle0!=0x0) {
914  std::cout<<"fTruth PDG: "<<mcparticle0->PdgCode()<<std::endl;
915  ftruthpdg=mcparticle0->PdgCode();
916 
917  truthid=mcparticle0->TrackId();
918  fprimary_truth_Isbeammatched=0;
919  if(beamid==truthid) fprimary_truth_Isbeammatched=1;
920  }
921 
922  fisprimarytrack = 1;
923  fisprimaryshower = 0;
924 
925  //fprimaryID = thisTrack->ParticleId(); //Do NOT use ParticleID, u got nothing
926  fprimaryID = thisTrack->ID();
927  fprimaryTheta = thisTrack->Theta();
928  fprimaryPhi = thisTrack->Phi();
929  fprimaryLength = thisTrack->Length();
930  fprimaryMomentum = thisTrack->StartMomentum();
931  fprimaryEndMomentum = thisTrack->EndMomentum();
932 
933  fprimaryEndPosition[0] = thisTrack->Trajectory().End().X();
934  fprimaryEndPosition[1] = thisTrack->Trajectory().End().Y();
935  fprimaryEndPosition[2] = thisTrack->Trajectory().End().Z();
936  fprimaryStartPosition[0] = thisTrack->Trajectory().Start().X();
937  fprimaryStartPosition[1] = thisTrack->Trajectory().Start().Y();
938  fprimaryStartPosition[2] = thisTrack->Trajectory().Start().Z();
939  fprimaryEndDirection[0] = thisTrack->Trajectory().EndDirection().X();
940  fprimaryEndDirection[1] = thisTrack->Trajectory().EndDirection().Y();
941  fprimaryEndDirection[2] = thisTrack->Trajectory().EndDirection().Z();
942  fprimaryStartDirection[0] = thisTrack->Trajectory().StartDirection().X();
943  fprimaryStartDirection[1] = thisTrack->Trajectory().StartDirection().Y();
944  fprimaryStartDirection[2] = thisTrack->Trajectory().StartDirection().Z();
945 
946  fprimaryMomentumByRangeMuon = trmom.GetTrackMomentum(thisTrack->Length(),13);
947  fprimaryMomentumByRangeProton = trmom.GetTrackMomentum(thisTrack->Length(),2212);
948 
949 
950  //Get the clusters/hits associated with the PFParticle -------------------------------------------------------------------------//
951  //const std::vector<const recob::Cluster*> thisCluster = pfpUtil.GetPFParticleClusters(*particle,evt,fPFParticleTag);
952  //std::cout<<"size of cluster:"<<thisCluster.size()<<std::endl; //should be three (u,v,c)
953 
954  //Get the number of clusters associated to the PFParticle
955  //unsigned int num_pfp_clusters = pfpUtil.GetNumberPFParticleClusters(*particle, evt, fPFParticleTag);
956  //std::cout<<"num_pfp_clusters:"<<num_pfp_clusters<<std::endl;
957  //std::cout<<"pfplist.size="<<pfplist.size()<<std::endl; //all the pfp particles in a pool
958 
959  //Get the number of hits
960  //unsigned int num_hits=pfpUtil.GetNumberPFParticleHits(*particle, evt, fPFParticleTag);
961 
962  //Get the hits associated to the PFParticle
963  const std::vector<const recob::Hit*> pfpHits=pfpUtil.GetPFParticleHits(*particle, evt, fPFParticleTag);
964  //std::cout<<"size of hit:"<<pfpHits.size()<<std::endl;
965  //std::cout<<"num_hits:"<<num_hits<<std::endl;
966 
967  for(const recob::Hit* hi : pfpHits){
968  //std::cout<<"(w,t,ID,TPC,Cryostat):("<<hi->WireID().Wire<<","<<hi->PeakTime()<<","<<hi->WireID().Plane<<","<<hi->WireID().TPC<<","<<hi->WireID().Cryostat<<")"<<std::endl;
969  if (hi->WireID().Plane==2) {
970  pfphit_peaktime_c.push_back((double)hi->PeakTime());
971  pfphit_wireid_c.push_back((double)hi->WireID().Wire);
972  }
973  if (hi->WireID().Plane==1) {
974  pfphit_peaktime_v.push_back((double)hi->PeakTime());
975  pfphit_wireid_v.push_back((double)hi->WireID().Wire);
976  }
977  if (hi->WireID().Plane==0) {
978  pfphit_peaktime_u.push_back((double)hi->PeakTime());
979  pfphit_wireid_u.push_back((double)hi->WireID().Wire);
980  }
981  }
982  //-------------------------------------------------------------------------------------------------------------------------------//
983 
984 
985  std::vector<anab::Calorimetry> calovector = trackUtil.GetRecoTrackCalorimetry(*thisTrack, evt, fTrackerTag, fCalorimetryTag);
986  if(calovector.size() != 3)
987  std::cerr << "WARNING::Calorimetry vector size for primary is = " << calovector.size() << std::endl;
988 
989  //HY::Get the Calorimetry(s) from thisTrack
990  std::vector<double> tmp_primtrk_dqdx;
991  std::vector<double> tmp_primtrk_resrange;
992  std::vector<double> tmp_primtrk_dedx;
993  std::vector<double> tmp_primtrk_hitx;
994  std::vector<double> tmp_primtrk_hity;
995  std::vector<double> tmp_primtrk_hitz;
996  std::vector<double> tmp_primtrk_pitch;
997  //std::vector<double> tmp_primtrk_truth_Eng;
998  for (auto & calo : calovector) {
999  if (calo.PlaneID().Plane == 2){ //only collection plane
1000  primtrk_range.push_back(calo.Range());
1001  std::cout<<"primtrk_range:"<<calo.Range()<<std::endl;
1002  for (size_t ihit = 0; ihit < calo.dQdx().size(); ++ihit){ //loop over hits
1003  tmp_primtrk_dqdx.push_back(calo.dQdx()[ihit]);
1004  tmp_primtrk_resrange.push_back(calo.ResidualRange()[ihit]);
1005  tmp_primtrk_dedx.push_back(calo.dEdx()[ihit]);
1006  tmp_primtrk_pitch.push_back(calo.TrkPitchVec()[ihit]);
1007 
1008  const auto &primtrk_pos=(calo.XYZ())[ihit];
1009  tmp_primtrk_hitx.push_back(primtrk_pos.X());
1010  tmp_primtrk_hity.push_back(primtrk_pos.Y());
1011  tmp_primtrk_hitz.push_back(primtrk_pos.Z());
1012  //std::cout<<"dqdx="<<calo.dQdx()[ihit]<<"; resrange="<<calo.ResidualRange()[ihit]<<std::endl;
1013  //std::cout<<"(X,Y,Z)="<<"("<<calo.XYZ()[ihit].X()<<","<<calo.XYZ()[ihit].Y()<<","<<calo.XYZ()[ihit].Z()<<")"<<std::endl;
1014 
1015  } //loop over hits
1016  } //only collection plane
1017  }
1018 
1019  //HY::Associtate the reco info with the truth info using backtracker (only track, no shower) -------------------------------------------------------------------------------//
1020  // Get the truth info
1021  //const simb::MCParticle* mcparticle2 = truthUtil.GetMCParticleFromRecoTrack(*thisTrack, evt, fTrackerTag);
1022  const simb::MCParticle* geantGoodParticle1 = truthUtil.GetGeantGoodParticle((*mcTruths)[0],evt);
1023  //if(mcparticle2 != 0x0) { //mcparticle
1024 
1025  //std::vector<unsigned char> tmp_primtrk_hit_processkey;
1026  std::vector<double> tmp_primtrk_hitx_true;
1027  std::vector<double> tmp_primtrk_hity_true;
1028  std::vector<double> tmp_primtrk_hitz_true;
1029  std::vector<double> tmp_primtrk_trkid_true;
1030  std::vector<double> tmp_primtrk_edept_true;
1031 
1032  std::vector<double> tmp_primtrk_true_x;
1033  std::vector<double> tmp_primtrk_true_y;
1034  std::vector<double> tmp_primtrk_true_z;
1035  std::vector<double> tmp_primtrk_true_trkid;
1036  std::vector<double> tmp_primtrk_true_edept;
1037 
1038  std::vector<double> tmp_primtrj_true_x;
1039  std::vector<double> tmp_primtrj_true_y;
1040  std::vector<double> tmp_primtrj_true_z;
1041  std::vector<double> tmp_primtrj_true_edept;
1042 
1043  if(geantGoodParticle1!= 0x0 && geantGoodParticle1->Process()=="primary") { //sansity check
1044  //const simb::MCParticle *geantGoodParticle=pi_serv->TrackIdToMotherParticle_P(mcparticle2->TrackId());
1045  //if(geantGoodParticle != 0x0 && geantGoodParticle->Process()=="primary") { //geatGoodParticle and primary p Loop
1046  double ke_reco=ke_ff; //ke_reco at ff
1047  double ke_true=ke_ff; //ke_true at ff
1048 
1049  art::ServiceHandle<cheat::BackTrackerService> bt_serv;
1050  art::ServiceHandle<geo::Geometry> geom;
1051  simb::MCTrajectory truetraj=geantGoodParticle1->Trajectory();
1052  auto thisTrajectoryProcessMap1 = truetraj.TrajectoryProcesses();
1053  std::string int_label="";
1054  //double endz_true=-999;
1055  if (thisTrajectoryProcessMap1.size()){
1056  for(auto const& couple: thisTrajectoryProcessMap1){
1057  //endz_true=((truetraj.at(couple.first)).first).Z();
1058  int_label=truetraj.KeyToProcess(couple.second);
1059 
1060  fprimary_truth_EndPosition[0]=((truetraj.at(couple.first)).first).X();
1061  fprimary_truth_EndPosition[1]=((truetraj.at(couple.first)).first).Y();
1062  fprimary_truth_EndPosition[2]=((truetraj.at(couple.first)).first).Z();
1063  fprimary_truth_EndProcess=truetraj.KeyToProcess(couple.second);
1064 
1065  fprimary_truth_EndMomentum[0]=((truetraj.at(couple.first)).second).X();
1066  fprimary_truth_EndMomentum[1]= ((truetraj.at(couple.first)).second).Y();
1067  fprimary_truth_EndMomentum[2]=((truetraj.at(couple.first)).second).Z();
1068 
1069  break;
1070  }
1071  }
1072  primtrk_end_g4process=int_label;
1073 
1074  //study of interaction angles
1075  //if (thisTrajectoryProcessMap1.size()) { //TrajectoryProcessMap1
1076  if (thisTrajectoryProcessMap1.size()) { //TrajectoryProcessMap1
1077  for(auto const& couple1: thisTrajectoryProcessMap1) { //go through this traj with all the interaction vertices
1078  interactionX.push_back(((truetraj.at(couple1.first)).first).X());
1079  interactionY.push_back(((truetraj.at(couple1.first)).first).Y());
1080  interactionZ.push_back(((truetraj.at(couple1.first)).first).Z());
1081  interactionE.push_back(((truetraj.at(couple1.first)).first).E());
1082  interactionProcesslist.push_back(truetraj.KeyToProcess(couple1.second));
1083 
1084  //get the TPC num
1085  double xval=((truetraj.at(couple1.first)).first).X();
1086  double yval=((truetraj.at(couple1.first)).first).Y();
1087  double zval=((truetraj.at(couple1.first)).first).Z();
1088  unsigned int tpcno=1;
1089  if(xval<=0 && zval<232) tpcno=1;
1090  if(xval<=0 && zval>232 && zval<464) tpcno=5;
1091  if(xval<=0 && zval>=464) tpcno=9;
1092  if(xval>0 && zval<232) tpcno=2;
1093  if(xval>0 && zval>232 && zval<464) tpcno=6;
1094  if(xval>0 && zval>=464) tpcno=10;
1095 
1096  //convert the position of the interaction vertex to (wireID, peak time)
1097  interaction_wid_c.push_back(fGeometry->WireCoordinate(yval, zval, 2, tpcno, 0));
1098  interaction_wid_v.push_back(fGeometry->WireCoordinate(yval, zval, 1, tpcno, 0));
1099  interaction_wid_u.push_back(fGeometry->WireCoordinate(yval, zval, 0, tpcno, 0));
1100 
1101  interaction_tt_c.push_back(detProp.ConvertXToTicks(xval, 2, tpcno, 0));
1102  interaction_tt_v.push_back(detProp.ConvertXToTicks(xval, 1, tpcno, 0));
1103  interaction_tt_u.push_back(detProp.ConvertXToTicks(xval, 0, tpcno, 0));
1104 
1105  //interactionT.push_back(detProp.ConvertXToTicks(xval, 2, tpcno, 0));
1106  //interactionU.push_back(fGeometry->WireCoordinate(((truetraj.at(couple1.first)).first).Y(), ((truetraj.at(couple1.first)).first).Z(),0, tpcno, 0));
1107  //interactionV.push_back(fGeometry->WireCoordinate(((truetraj.at(couple1.first)).first).Y(), ((truetraj.at(couple1.first)).first).Z(),1, tpcno, 0));
1108  //interactionW.push_back(fGeometry->WireCoordinate(((truetraj.at(couple1.first)).first).Y(), ((truetraj.at(couple1.first)).first).Z(),2, tpcno, 0));
1109 
1110  //convert(x,y,z) to (wireid ,time ticks)
1111  //geo::PlaneID collection_plane = geom->PlaneIDs(2);
1112  //std::cout<<"fprimaryT0:"<<fprimaryT0<<std::endl;
1113  //std::cout<<"\nint_vtx x/y/z:"<<((truetraj.at(couple1.first)).first).X()<<"/"<<((truetraj.at(couple1.first)).first).Y()<<"/"<<((truetraj.at(couple1.first)).first).Z()<<std::endl;
1114  //std::cout<<"(wid,tt)_c:"<<"("<<fGeometry->WireCoordinate(yval, zval, 2, tpcno, 0)<<","<<detProp.ConvertXToTicks(xval, 2, tpcno, 0)<<")"<<std::endl;
1115  //std::cout<<"(wid,tt)_v:"<<"("<<fGeometry->WireCoordinate(yval, zval, 1, tpcno, 0)<<","<<detProp.ConvertXToTicks(xval, 1, tpcno, 0)<<")"<<std::endl;
1116  //std::cout<<"(wid,tt)_u:"<<"("<<fGeometry->WireCoordinate(yval, zval, 0, tpcno, 0)<<","<<detProp.ConvertXToTicks(xval, 0, tpcno, 0)<<")"<<std::endl;
1117 
1118 
1119  //not interested of CoulombScat
1120  if ((truetraj.KeyToProcess(couple1.second)).find("CoulombScat")!= std::string::npos) continue;
1121 
1122  //check if the interaction is in the TPC
1123  auto interactionPos4D = (truetraj.at(couple1.first)).first ;
1124 
1125  if (interactionPos4D.Z() < minZ || interactionPos4D.Z() > maxZ ) continue;
1126  else if (interactionPos4D.X() < minX || interactionPos4D.X() > maxX ) continue;
1127  else if (interactionPos4D.Y() < minY || interactionPos4D.Y() > maxY ) continue;
1128 
1129  ///get the interaction angle here
1130  double interactionAngle = 999999.; // This needs to be changed
1131  //--------------------- Int Angle ---------------------------
1132  // Try to retreive the interaction angle
1133  auto prevInteractionPos4D = (truetraj.at(couple1.first-1)).first ;
1134  auto prevInteractionPos3D = prevInteractionPos4D.Vect() ;
1135  auto interactionPos3D = interactionPos4D.Vect() ;
1136  auto distanceBtwPoint = interactionPos3D - prevInteractionPos3D;
1137 
1138  //see if the next point exists
1139  if (truetraj.size() > couple1.first + 1) {
1140  // The particle doesn't die. No need to check for anything else.
1141  auto nextInteractionPos4D = (truetraj.at(couple1.first+1)).first ;
1142  auto nextInteractionPos3D = nextInteractionPos4D.Vect() ;
1143  auto distanceBtwPointNext = nextInteractionPos3D - interactionPos3D;
1144  interactionAngles.push_back(TMath::ACos(distanceBtwPointNext.Dot(distanceBtwPoint)/(distanceBtwPointNext.Mag()*distanceBtwPoint.Mag() ) ));
1145  }
1146  else { // The particle has come to an end. Let's check the daugthers.
1147  if (geantGoodParticle1->NumberDaughters() == 0 ){
1148  interactionAngles.push_back(interactionAngle);
1149  break;
1150  }
1151  double maxAngle = 0.;
1152  int numberOfTroubleDaugther = 0;
1153  //Loop on the daughters
1154  const sim::ParticleList& plist=pi_serv->ParticleList();
1155  sim::ParticleList::const_iterator itPart1=plist.begin();
1156  for(size_t dPart1=0;(dPart1<plist.size()) && (plist.begin()!=plist.end());++dPart1) {
1157  const simb::MCParticle* dPart=(itPart1++)->second;
1158  if (dPart->Mother() != 1 ) continue;
1159  auto daugthTraj = dPart->Trajectory();
1160  // if (debug) std::cout<< dPart->PdgCode()
1161  // <<" , length: "<< (dPart->Trajectory()).TotalLength() <<"\n";
1162  //<<"First Point: "<< ((daugthTraj[0].first).Vect()).X() <<","<< ((daugthTraj[0].first).Vect()).Y()<<","<<((daugthTraj[0].first).Vect()).Z() <<"\n";
1163  if ((dPart->NumberTrajectoryPoints () < 2 )||!(TMath::Abs(dPart->PdgCode()) == 13||TMath::Abs(dPart->PdgCode()) == 11 || TMath::Abs(dPart->PdgCode()) == 211 || TMath::Abs(dPart->PdgCode()) == 321 || TMath::Abs(dPart->PdgCode()) == 2212)||(daugthTraj.TotalLength() < 0.5 )) {
1164  interactionAngle=999999;
1165  continue;
1166  }
1167 
1168  numberOfTroubleDaugther++;
1169 
1170  auto daughtFirstPt = ((daugthTraj[0]).first).Vect() ;
1171  auto daughtSecondPt = ((daugthTraj[1]).first).Vect() ;
1172  auto distanceBtwPointNext = daughtSecondPt - daughtFirstPt;
1173  interactionAngle = TMath::ACos(distanceBtwPointNext.Dot(distanceBtwPoint)/(distanceBtwPointNext.Mag()*distanceBtwPoint.Mag()));
1174  if ( maxAngle < interactionAngle ) maxAngle = interactionAngle;
1175  interactionAngle = maxAngle;
1176  // If the track finishes without visible daugthers, we're likely to see it. Give it a big angle!
1177  if (!numberOfTroubleDaugther) interactionAngle = 9999999.;//It's a huge angle: it's greater than the cut, so I don't remove this! But it also doesn't go into the angle plot
1178  if (interactionAngle < 0.0001) std::cout<<"-------------------------------------------------> \n";
1179  interactionAngles.push_back(interactionAngle);
1180  break;
1181  }
1182  } // The particle has come to an end. Let's check the daugthers.
1183  } //go through this traj with all the interaction vertices
1184  } //TrajectoryProcessMap1
1185 
1186  //save all the mctraj info
1187  std::cout<<"\n MCTrajectory of this prim. trk has:"<<truetraj.size()<<" hits!!"<<std::endl;
1188  for (size_t tt=0; tt<truetraj.size(); ++tt) { //loop over all the hits in this mc traj
1189  tmp_primtrj_true_x.push_back(truetraj.X(tt));
1190  tmp_primtrj_true_y.push_back(truetraj.Y(tt));
1191  tmp_primtrj_true_z.push_back(truetraj.Z(tt));
1192  tmp_primtrj_true_edept.push_back(truetraj.E(tt));
1193  //std::cout<<"[mctraj] x/y/z/E:"<<truetraj.X(tt)<<"/"<<truetraj.Y(tt)<<"/"<<truetraj.Z(tt)<<"/"<<truetraj.E(tt)<<std::endl;
1194  //std::cout<<"[int] x/y/z/E:"<<interactionX.at(tt)<<"/"<<interactionY.at(tt)<<"/"<<interactionZ.at(tt)<<"/"<<interactionE.at(tt)<<std::endl;
1195  } //loop over all the hits in this mc traj
1196 
1197  //use backtracker to find the associtated true info
1198  geo::View_t view = geom->View(2);
1199  auto simIDE_prim=bt_serv->TrackIdToSimIDEs_Ps(geantGoodParticle1->TrackId(),view);
1200  std::map<double, sim::IDE> orderedSimIDE; //id & e
1201  for (auto& ide : simIDE_prim) orderedSimIDE[ide->z]= *ide; //order in z-direction
1202  auto inTPCPoint = truetraj.begin();
1203  auto Momentum0 = inTPCPoint->second;
1204  //auto old_iter = orderedSimIDE.begin();
1205  //std::cout<<"True KE at front face: "<<ke_true<<"; Reconstructed KE at ff: "<<ke_reco<<std::endl;
1206  primtrk_ke_reco.push_back(ke_reco);
1207  primtrk_ke_true.push_back(ke_true);
1208 
1209  double trklen=0.0;
1210  bool run_me_onetime=true;
1211  double xi=0.0; double yi=0.0; double zi=0.0;
1212  int cnt=0;
1213  //sanity check on the start/end positions
1214  if(tmp_primtrk_dqdx.size()!=0) { //make sure the dqdx vector is not empty
1215  if(tmp_primtrk_hitz[0]>tmp_primtrk_hitz[tmp_primtrk_hitz.size()-1]) { //flip trk vectors if Pandora messes up the trk direction
1216  std::reverse(tmp_primtrk_hitz.begin(), tmp_primtrk_hitz.end());
1217  std::reverse(tmp_primtrk_hity.begin(), tmp_primtrk_hity.end());
1218  std::reverse(tmp_primtrk_hitx.begin(), tmp_primtrk_hitx.end());
1219  std::reverse(tmp_primtrk_pitch.begin(),tmp_primtrk_pitch.end());
1220  std::reverse(tmp_primtrk_dedx.begin(), tmp_primtrk_dedx.end());
1221  std::reverse(tmp_primtrk_dqdx.begin(), tmp_primtrk_dqdx.end());
1222  std::reverse(tmp_primtrk_resrange.begin(), tmp_primtrk_resrange.end());
1223  } //flip trk vectors if Pandora messes up the trk direction
1224  //p.s. simIDE can only get Edept, needs to calculate true KE by your own
1225 
1226  for(size_t idx1=0; idx1<tmp_primtrk_dqdx.size()-1; idx1++) { //energy deposition: reco hit loop
1227  ke_reco-=tmp_primtrk_pitch[idx1]*tmp_primtrk_dedx[idx1];
1228  double currentDepEng = 0.;
1229  for ( auto iter= orderedSimIDE.begin(); iter!= orderedSimIDE.end(); iter++) { //simIde loop ---------------------------------------------------------------------------//
1230  auto currentIde = iter->second;
1231 
1232  //calculation only inside TPC
1233  if(currentIde.z<minZ || currentIde.z > maxZ ) continue;
1234  else if (currentIde.x < minX || currentIde.x > maxX ) continue;
1235  else if (currentIde.y < minY || currentIde.y > maxY ) continue;
1236 
1237  //if(cnt==0) { //start positions
1238  if(cnt==0&&currentIde.trackID>=0) { //start positions
1239  fprimary_truth_StartPosition[0] = currentIde.x;
1240  fprimary_truth_StartPosition[1] = currentIde.y;
1241  fprimary_truth_StartPosition[2] = currentIde.z;
1242 
1243  //tmp_primtrk_hit_processkey.push_back(currentIde.ProcessToKey);
1244  tmp_primtrk_hitx_true.push_back(currentIde.x);
1245  tmp_primtrk_hity_true.push_back(currentIde.y);
1246  tmp_primtrk_hitz_true.push_back(currentIde.z);
1247  tmp_primtrk_trkid_true.push_back(currentIde.trackID);
1248  tmp_primtrk_edept_true.push_back(currentIde.energy);
1249  } //start positions
1250 
1251  //run the simIde loop only one time to sum over all trk length in each segment
1252  //if (currentIde.trackID!=-1) { //discard true shower info
1253  //if(cnt>0 && run_me_onetime) { //run at once
1254  if (currentIde.trackID>=0) { //discard true shower info (trackID<0)
1255  if(cnt>0 && run_me_onetime) { //run at once
1256  trklen+=TMath::Sqrt(std::pow(currentIde.x-xi,2)+std::pow(currentIde.y-yi,2)+std::pow(currentIde.z-zi,2));
1257  //std::cout<<"trklen: "<<trklen<<std::endl;
1258  //std::cout<<"dx, dy and dz: "<<currentIde.x-xi<<" "<<currentIde.y-yi<<" "<<currentIde.z-zi<<std::endl;
1259  //std::cout<<"x,y,z: "<<currentIde.x<<", "<<currentIde.y<<", "<<currentIde.z<<"; trkId: "<<currentIde.trackID<<"; energy: "<<currentIde.energy<<std::endl;
1260 
1261  //fprimary_truth_EndPosition[0]=currentIde.x;
1262  //fprimary_truth_EndPosition[1]=currentIde.y;
1263  //fprimary_truth_EndPosition[2]=currentIde.z;
1264 
1265  //tmp_primtrk_hit_processkey.push_back(currentIde.ProcessToKey);
1266  tmp_primtrk_hitx_true.push_back(currentIde.x);
1267  tmp_primtrk_hity_true.push_back(currentIde.y);
1268  tmp_primtrk_hitz_true.push_back(currentIde.z);
1269  tmp_primtrk_trkid_true.push_back(currentIde.trackID);
1270  tmp_primtrk_edept_true.push_back(currentIde.energy);
1271  } //run at once
1272  xi=currentIde.x; yi=currentIde.y; zi=currentIde.z;
1273  cnt++;
1274  //std::cout<<"cnt:"<<cnt<<std::endl;
1275  } //discard true shower info
1276 
1277  //true E dept within the thin slice (simIde can only get Edept, not KE)
1278  if ( currentIde.z <= tmp_primtrk_hitz[idx1]) continue;
1279  if ( currentIde.z > tmp_primtrk_hitz[idx1+1]) continue;
1280  currentDepEng += currentIde.energy; //total energy deposition in the current slice
1281  } //simIde loop -------------------------------------------------------------------------------------------------------------------------------------------------------//
1282  ke_true -= currentDepEng; //KE in the current slice
1283  run_me_onetime=false; //finished the summation of the true trk length
1284 
1285  if(currentDepEng>0.001) { //remove the zombie tracks
1286  primtrk_ke_reco.push_back(ke_reco);
1287  primtrk_ke_true.push_back(ke_true);
1288  } //remove the zombie tracks
1289 
1290  //std::cout<<"primtrk_ke_reco["<<idx1<<"]:"<<ke_reco<<" ; primtrk_ke_true:"<<ke_true<<std::endl;
1291 
1292  //if(currentDepEng>0.001) hKE_truth_reco->Fill(kineticEnergy,KE_rec);
1293  }//energy deposition: reco hit loop
1294 
1295 
1296  //sime IDE loop again to save all the true info
1297  for ( auto iter2=orderedSimIDE.begin(); iter2!=orderedSimIDE.end(); iter2++) { //simIde loop2 ---------------------------------------------------------------------------//
1298  auto currentIde2 = iter2->second;
1299 
1300  //calculation only inside TPC
1301  if(currentIde2.z<minZ || currentIde2.z > maxZ ) continue;
1302  else if (currentIde2.x < minX || currentIde2.x > maxX ) continue;
1303  else if (currentIde2.y < minY || currentIde2.y > maxY ) continue;
1304 
1305  //if(currentIde2.trackID>=0) { // no shower
1306  //tmp_primtrk_hit_processkey.push_back(currentIde2.ProcessToKey);
1307  tmp_primtrk_true_x.push_back(currentIde2.x);
1308  tmp_primtrk_true_y.push_back(currentIde2.y);
1309  tmp_primtrk_true_z.push_back(currentIde2.z);
1310  tmp_primtrk_true_trkid.push_back(currentIde2.trackID);
1311  tmp_primtrk_true_edept.push_back(currentIde2.energy);
1312 
1313  //std::cout<<"[simeide2] x,y,z: "<<currentIde2.x<<", "<<currentIde2.y<<", "<<currentIde2.z<<"; trkId: "<<currentIde2.trackID<<"; energy: "<<currentIde2.energy<<std::endl;
1314  //} //no shower
1315  } //simIde loop2 -------------------------------------------------------------------------------------------------------------------------------------------------------//
1316 
1317 
1318  } //make sure that primtrk_tmp size g.t. 0
1319 
1320  //} //geatGoodParticle and primary p Loop
1321 
1322  //std::cout<<"primtrk_range_true:"<<trklen<<std::endl;
1323  primtrk_range_true.push_back(trklen);
1324 
1325  } //mcpartice
1326  //---------------------------------------------------------------------------------------------------------------//
1327 
1328 
1329  primtrk_dqdx.push_back(tmp_primtrk_dqdx);
1330  primtrk_resrange.push_back(tmp_primtrk_resrange);
1331  primtrk_dedx.push_back(tmp_primtrk_dedx);
1332  primtrk_hitx.push_back(tmp_primtrk_hitx);
1333  primtrk_hity.push_back(tmp_primtrk_hity);
1334  primtrk_hitz.push_back(tmp_primtrk_hitz);
1335  primtrk_pitch.push_back(tmp_primtrk_pitch);
1336 
1337  //primtrk_hit_processkey.push_back(tmp_primtrk_hit_processkey);
1338  primtrk_hitx_true.push_back(tmp_primtrk_hitx_true);
1339  primtrk_hity_true.push_back(tmp_primtrk_hity_true);
1340  primtrk_hitz_true.push_back(tmp_primtrk_hitz_true);
1341  primtrk_trkid_true.push_back(tmp_primtrk_trkid_true);
1342  primtrk_edept_true.push_back(tmp_primtrk_edept_true);
1343 
1344  primtrk_true_x.push_back(tmp_primtrk_true_x);
1345  primtrk_true_y.push_back(tmp_primtrk_true_y);
1346  primtrk_true_z.push_back(tmp_primtrk_true_z);
1347  primtrk_true_trkid.push_back(tmp_primtrk_true_trkid);
1348  primtrk_true_edept.push_back(tmp_primtrk_true_edept);
1349 
1350  primtrj_true_x.push_back(tmp_primtrj_true_x);
1351  primtrj_true_y.push_back(tmp_primtrj_true_y);
1352  primtrj_true_z.push_back(tmp_primtrj_true_z);
1353  primtrj_true_edept.push_back(tmp_primtrj_true_edept);
1354 
1355  tmp_primtrk_dqdx.clear();
1356  tmp_primtrk_resrange.clear();
1357  tmp_primtrk_dedx.clear();
1358  tmp_primtrk_hitx.clear();
1359  tmp_primtrk_hity.clear();
1360  tmp_primtrk_hitz.clear();
1361  tmp_primtrk_pitch.clear();
1362 
1363  //tmp_primtrk_hit_processkey.clear();
1364  tmp_primtrk_hitx_true.clear();
1365  tmp_primtrk_hity_true.clear();
1366  tmp_primtrk_hitz_true.clear();
1367  tmp_primtrk_trkid_true.clear();
1368  tmp_primtrk_edept_true.clear();
1369 
1370  tmp_primtrk_true_x.clear();
1371  tmp_primtrk_true_y.clear();
1372  tmp_primtrk_true_z.clear();
1373  tmp_primtrk_true_trkid.clear();
1374  tmp_primtrk_true_edept.clear();
1375 
1376  tmp_primtrj_true_x.clear();
1377  tmp_primtrj_true_y.clear();
1378  tmp_primtrj_true_z.clear();
1379  tmp_primtrj_true_edept.clear();
1380 
1381  for(size_t k = 0; k < calovector.size() && k<3; k++){
1382  fprimaryKineticEnergy[k] = calovector[k].KineticEnergy();
1383  fprimaryRange[k] = calovector[k].Range();
1384  //const std::vector< double > & dedxvec = calovector[k].dEdx();
1385  }
1386 
1387  //Get CNN score of each hit ------------------------------------------------------------------------//
1388  int planenum=999;
1389  float zpos=-999;
1390  float ypos=-999;
1391  float xpos=-999;
1392 
1393  //float max_inel_score_c=-999.;
1394  //float max_el_score_c=-999.;
1395  if(fmthm.isValid()){ //if non-empty fmthm
1396  auto vhit=fmthm.at(fprimaryID);
1397  auto vmeta=fmthm.data(fprimaryID); //indices of meta data are the same as data
1398  for (size_t ii = 0; ii<vhit.size(); ++ii){ //loop over all meta data hit
1399  bool fBadhit = false;
1400  if (vmeta[ii]->Index() == static_cast<unsigned int>(std::numeric_limits<int>::max())){
1401  fBadhit = true;
1402  //cout<<"fBadHit"<<fBadhit<<endl;
1403  continue;
1404  }
1405  if (vmeta[ii]->Index()>=tracklist[fprimaryID]->NumberTrajectoryPoints()){
1406  throw cet::exception("Calorimetry_module.cc") << "Requested track trajectory index "<<vmeta[ii]->Index()<<" exceeds the total number of trajectory points "<<tracklist[fprimaryID]->NumberTrajectoryPoints()<<" for track index "<<fprimaryID<<". Something is wrong with the track reconstruction. Please contact tjyang@fnal.gov!!";
1407  }
1408  if (!tracklist[fprimaryID]->HasValidPoint(vmeta[ii]->Index())){
1409  fBadhit = true;
1410  // cout<<"had valid point "<<fBadhit<<endl;
1411  continue;
1412  }
1413 
1414  //get (x,y,z)
1415  auto loc = tracklist[fprimaryID]->LocationAtPoint(vmeta[ii]->Index());
1416  xpos=loc.X();
1417  ypos=loc.Y();
1418  zpos=loc.Z();
1419  //std::cout<<"x, y, z: "<<xpos<<" "<<ypos<<" "<<zpos<<std::endl;
1420  //std::cout<<"BadHit"<<fBadhit<<std::endl;
1421 
1422  //get the TPC num
1423  unsigned int tpc_no=1;
1424  if(xpos<=0 && zpos<232) tpc_no=1;
1425  if(xpos<=0 && zpos>232 && zpos<464) tpc_no=5;
1426  if(xpos<=0 && zpos>=464) tpc_no=9;
1427  if(xpos>0 && zpos<232) tpc_no=2;
1428  if(xpos>0 && zpos>232 && zpos<464) tpc_no=6;
1429  if(xpos>0 && zpos>=464) tpc_no=10;
1430 
1431  //skip the bad hit if any
1432  if (fBadhit) continue; //HY::If BAD hit, skip this hit and go next
1433  if (zpos<-100) continue; //hit not on track
1434  planenum=vhit[ii]->WireID().Plane;
1435 
1436  if(planenum==2){
1437  //std::cout<<"inside loop"<<std::endl;
1438  std::array<float,3> cnn_out=hitResults.getOutput(vhit[ii]);
1439  //peakT_2.push_back(vhit[ii]->PeakTime());
1440  //int_2.push_back(vhit[ii]->Integral());
1441  //hitz_2.push_back(zpos);
1442  inelscore_c.push_back(cnn_out[hitResults.getIndex("inel")]);
1443  elscore_c.push_back(cnn_out[hitResults.getIndex("el")]);
1444  nonescore_c.push_back(cnn_out[hitResults.getIndex("none")]);
1445 
1446  //save the associtated 3d position
1447  x_c.push_back(xpos);
1448  y_c.push_back(ypos);
1449  z_c.push_back(zpos);
1450  //std::cout<<"(x_c/y_c/z_c): ("<<xpos<<","<<ypos<<","<<zpos<<")"<<std::endl;
1451 
1452  //convert the position of the interaction to (wireID, peak time)
1453  wid_c.push_back(fGeometry->WireCoordinate(ypos, zpos, planenum, tpc_no, 0));
1454  tt_c.push_back(detProp.ConvertXToTicks(xpos, planenum, tpc_no, 0));
1455 
1456  //std::cout<<"(w,t):("<<fGeometry->WireCoordinate(ypos, zpos, planenum, tpc_no, 0)<<","<<detProp.ConvertXToTicks(xpos, planenum, tpc_no, 0)<<")|"<<
1457  //"[inel,el,non]:["<<cnn_out[hitResults.getIndex("inel")]<<","<<cnn_out[hitResults.getIndex("el")]<<","<<cnn_out[hitResults.getIndex("none")]<<"]"<<std::endl;
1458  // std::cout<<"peaktime "<<vhit[ii]->PeakTime()<<std::endl;
1459 
1460  //if (cnn_out[hitResults.getIndex("inel")]>max_inel_score_c){ //select max. inel_score
1461  //max_inel_score_c=cnn_out[hitResults.getIndex("inel")];
1462  //xyz_inelscore_c[0]=xpos;
1463  //xyz_inelscore_c[1]=ypos;
1464  //xyz_inelscore_c[2]=zpos;
1465  //} //select max. inel_score
1466 
1467  //if (cnn_out[hitResults.getIndex("el")]>max_el_score_c){ //select max. el_score
1468  //max_el_score_c=cnn_out[hitResults.getIndex("el")];
1469  //xyz_elscore_c[0]=xpos;
1470  //xyz_elscore_c[1]=ypos;
1471  //xyz_elscore_c[2]=zpos;
1472  //} //select max. el_score
1473  }//planenum 2
1474  if(planenum==1){
1475  //int_1.push_back(vhit[ii]->Integral());
1476  //hitz_1.push_back(zpos);
1477 
1478  //convert the position of the interaction to (wireID, peak time)
1479  wid_v.push_back(fGeometry->WireCoordinate(ypos, zpos, planenum, tpc_no, 0));
1480  tt_v.push_back(detProp.ConvertXToTicks(xpos, planenum, tpc_no, 0));
1481 
1482  }//planenum 1
1483  if(planenum==0){
1484  //int_0.push_back(vhit[ii]->Integral());
1485  //hitz_0.push_back(zpos);
1486 
1487  //convert the position of the interaction to (wireID, peak time)
1488  wid_u.push_back(fGeometry->WireCoordinate(ypos, zpos, planenum, tpc_no, 0));
1489  tt_u.push_back(detProp.ConvertXToTicks(xpos, planenum, tpc_no, 0));
1490  }//planenum 0
1491 
1492 
1493 
1494  } //loop over all meta data hit
1495  } //if non-empty fmthm
1496 
1497 
1498 
1499 
1500  //Get CNN score of each hit ------------------------------------------------------------------------//
1501 
1502 
1503  //Use Ajib's intersection calculation function
1504  //[1]identify the beam track and tag other tracks
1505  std::vector<float> Stw, Endw, Stt, Endt, Stwires, Endwires, Stticks, Endticks, TPCb, TPCcl;
1506  Stw.clear(); Endw.clear(); Stt.clear(); Endt.clear(); Stwires.clear(); Endwires.clear(); Stticks.clear(); Endticks.clear(); TPCb.clear(); TPCcl.clear();
1507  float den;
1508  float numw, numt,wire_no,ticks_no;
1509  for(size_t p1=0;p1<pfplist.size();p1++){
1510  std::vector<art::Ptr<recob::Track>> trk=pftrack.at(p1);
1511  std::vector<art::Ptr<recob::Cluster>> allClusters=fmcp.at(p1);
1512  //std::cout<<"fprimaryID:"<<fprimaryID<<std::endl;
1513  for(size_t c1=0;c1<allClusters.size();c1++){
1514  if(allClusters[c1]->Plane().Plane!=2) continue;
1515  if(trk.size() && int(trk[0].key())==fprimaryID){
1516  Stw.push_back(allClusters[c1]->StartWire());
1517  Endw.push_back(allClusters[c1]->EndWire());
1518  Stt.push_back(allClusters[c1]->StartTick());
1519  Endt.push_back(allClusters[c1]->EndTick());
1520  TPCb.push_back(allClusters[c1]->Plane().TPC);
1521  //std::cout<<"\nst/end wire:"<<allClusters[c1]->StartWire()<<"/"<<allClusters[c1]->EndWire()<<std::endl;
1522  }
1523  else{
1524  Stwires.push_back(allClusters[c1]->StartWire());
1525  Endwires.push_back(allClusters[c1]->EndWire());
1526  Stticks.push_back(allClusters[c1]->StartTick());
1527  Endticks.push_back(allClusters[c1]->EndTick());
1528  TPCcl.push_back(allClusters[c1]->Plane().TPC);
1529  }
1530  }
1531  }
1532  //[2]find interaction points if any (assuming all tracks are straight, find the interaction points)
1533  for(size_t clt=0;clt<Stw.size();clt++){
1534  for(size_t cl1=0;cl1<Stwires.size();cl1++){
1535  if(TPCcl[cl1]!=TPCb[clt]) continue;
1536  std::cout<<"tpc are equal "<<std::endl;
1537  den=(Stw[clt]-Endw[clt])*(Stticks[cl1]-Endticks[cl1])-(Stt[clt]-Endt[clt])*(Stwires[cl1]-Endwires[cl1]);
1538  if(den==0) continue;
1539  numw=(Stw[clt]*Endt[clt]-Stt[clt]*Endw[clt])*(Stwires[cl1]-Endwires[cl1])-(Stw[clt]-Endw[clt])*(Stwires[cl1]*Endticks[cl1]-Stticks[cl1]*Endwires[cl1]);
1540  numt=(Stw[clt]*Endt[clt]-Stt[clt]*Endw[clt])*(Stticks[cl1]-Endticks[cl1])-(Stt[clt]-Endt[clt])*(Stwires[cl1]*Endticks[cl1]-Stticks[cl1]*Endwires[cl1]);
1541  wire_no=numw/den;
1542  ticks_no=numt/den;
1543  // std::cout<<"wireno and ticks not solution "<<wire_no<<" "<<ticks_no<<std::endl;
1544  if(((Stw[clt]<wire_no && Endw[clt]>wire_no)||(Stw[clt]>wire_no && Endw[clt]<wire_no))&&((Stt[clt]<ticks_no && Endt[clt]>ticks_no)||(Stt[clt]>ticks_no && Endt[clt]<ticks_no)) && ((Stwires[cl1]<wire_no && Endwires[cl1]>wire_no)||(Stwires[cl1]>wire_no && Endwires[cl1]<wire_no)) && ((Stticks[cl1]<ticks_no && Endticks[cl1]>ticks_no)||(Stticks[cl1]>ticks_no && Endticks[cl1]<ticks_no))) {
1545  std::cout<<"intersection wire and ticks are "<<std::round(wire_no)<<" "<<ticks_no<<" Stw Endw StT EndT "<<Stwires[cl1]<<" "<<Endwires[cl1]<<" "<<Stticks[cl1]<<" "<<Endticks[cl1]<<std::endl;
1546  double xyzStart[3];
1547  double xyzEnd[3];
1548  unsigned int wireno=std::round(wire_no);
1549  geo::WireID wireid(0,TPCb[clt],2,wireno);
1550  fGeometry->WireEndPoints(0,TPCb[clt],2,wireno, xyzStart, xyzEnd);
1551  std::cout<<"Z position of intersection = "<<xyzStart[2]<<" "<<xyzEnd[2]<<" "<<wireno<<std::endl;
1552  Zintersection.push_back(xyzStart[2]);
1553  Yintersection.push_back(xyzStart[1]);
1554  Xintersection.push_back(xyzStart[0]);
1555  timeintersection.push_back(ticks_no);
1556  }
1557 
1558  }
1559  }
1560 
1561 
1562 
1563 
1564 
1565  } //this track
1566  else if(thisShower != 0x0){
1567  fisprimarytrack = 0;
1568  fisprimaryshower = 1;
1569 
1570  fprimaryID = thisShower->ID();
1571  fprimaryLength = thisShower->Length();
1572  fprimaryShowerBestPlane = thisShower->best_plane();
1573  fprimaryOpeningAngle = thisShower->OpenAngle();
1574  fprimaryStartPosition[0] = thisShower->ShowerStart().X();
1575  fprimaryStartPosition[1] = thisShower->ShowerStart().Y();
1576  fprimaryStartPosition[2] = thisShower->ShowerStart().Z();
1577  fprimaryStartDirection[0] = thisShower->Direction().X();
1578  fprimaryStartDirection[1] = thisShower->Direction().Y();
1579  fprimaryStartDirection[2] = thisShower->Direction().Z();
1580  if( (thisShower->Energy()).size() > 0 )
1581  fprimaryShowerEnergy = thisShower->Energy()[0]; // thisShower->best_plane()
1582  if( (thisShower->MIPEnergy()).size() > 0 )
1583  fprimaryShowerMIPEnergy = thisShower->MIPEnergy()[0];
1584  if( (thisShower->dEdx()).size() > 0 )
1585  fprimaryShowerdEdx = thisShower->dEdx()[0];
1586  }
1587  else{
1588  std::cout << "INFO::Primary pfParticle is not track or shower. Skip!" << std::endl;
1589  continue;
1590  }
1591 
1592  // Find the particle vertex. We need the tracker tag here because we need to do a bit of
1593  // additional work if the PFParticle is track-like to find the vertex.
1594  const TVector3 vtx = pfpUtil.GetPFParticleVertex(*particle,evt,fPFParticleTag,fTrackerTag);
1595  fvertex[0] = vtx.X(); fvertex[1] = vtx.Y(); fvertex[2] = vtx.Z();
1596 
1597  // Now we can look for the interaction point of the particle if one exists, i.e where the particle
1598  // scatters off an argon nucleus. Shower-like objects won't have an interaction point, so we can
1599  // check this by making sure we get a sensible position
1600  const TVector3 interactionVtx = pfpUtil.GetPFParticleSecondaryVertex(*particle,evt,fPFParticleTag,fTrackerTag);
1601  fsecvertex[0] = interactionVtx.X(); fsecvertex[1] = interactionVtx.Y(); fsecvertex[2] = interactionVtx.Z();
1602 
1603  // Maximum number of daugthers to be processed
1604  if(particle->NumDaughters() > NMAXDAUGTHERS)
1605  std::cout << "INFO::Number of daughters is " << particle->NumDaughters() << ". Only the first NMAXDAUGTHERS are processed." << std::endl;
1606 
1607  // Let's get the daughter PFParticles... we can do this simply without the utility
1608  for(const int daughterID : particle->Daughters()){
1609  // Daughter ID is the element of the original recoParticle vector
1610  const recob::PFParticle *daughterParticle = &(recoParticles->at(daughterID));
1611  std::cout << "Daughter " << daughterID << " has " << daughterParticle->NumDaughters() << " daughters" << std::endl;
1612 
1613  const recob::Track* daughterTrack = pfpUtil.GetPFParticleTrack(*daughterParticle,evt,fPFParticleTag,fTrackerTag);
1614  const recob::Shower* daughterShower = pfpUtil.GetPFParticleShower(*daughterParticle,evt,fPFParticleTag,fShowerTag);
1615 
1616  if(daughterTrack != 0x0){
1617  fisdaughtertrack[fNDAUGHTERS] = 1;
1618  fisdaughtershower[fNDAUGHTERS] = 0;
1619  fdaughterTheta[fNDAUGHTERS] = daughterTrack->Theta();
1620  fdaughterPhi[fNDAUGHTERS] = daughterTrack->Phi();
1621  fdaughterLength[fNDAUGHTERS] = daughterTrack->Length();
1622  fdaughterMomentum[fNDAUGHTERS] = daughterTrack->StartMomentum();
1623  fdaughterEndMomentum[fNDAUGHTERS] = daughterTrack->EndMomentum();
1624  fdaughterStartPosition[fNDAUGHTERS][0] = daughterTrack->Trajectory().Start().X();
1625  fdaughterStartPosition[fNDAUGHTERS][1] = daughterTrack->Trajectory().Start().Y();
1626  fdaughterStartPosition[fNDAUGHTERS][2] = daughterTrack->Trajectory().Start().Z();
1627  fdaughterEndPosition[fNDAUGHTERS][0] = daughterTrack->Trajectory().End().X();
1628  fdaughterEndPosition[fNDAUGHTERS][1] = daughterTrack->Trajectory().End().Y();
1629  fdaughterEndPosition[fNDAUGHTERS][2] = daughterTrack->Trajectory().End().Z();
1630  fdaughterStartDirection[fNDAUGHTERS][0] = daughterTrack->Trajectory().StartDirection().X();
1631  fdaughterStartDirection[fNDAUGHTERS][1] = daughterTrack->Trajectory().StartDirection().Y();
1632  fdaughterStartDirection[fNDAUGHTERS][2] = daughterTrack->Trajectory().StartDirection().Z();
1633  fdaughterEndDirection[fNDAUGHTERS][0] = daughterTrack->Trajectory().EndDirection().X();
1634  fdaughterEndDirection[fNDAUGHTERS][1] = daughterTrack->Trajectory().EndDirection().Y();
1635  fdaughterEndDirection[fNDAUGHTERS][2] = daughterTrack->Trajectory().EndDirection().Z();
1636 
1637  fdaughterMomentumByRangeMuon[fNDAUGHTERS] = trmom.GetTrackMomentum(daughterTrack->Length(),13);
1638  fdaughterMomentumByRangeProton[fNDAUGHTERS] = trmom.GetTrackMomentum(daughterTrack->Length(),2212);
1639 
1640  std::vector<anab::Calorimetry> daughtercalovector = trackUtil.GetRecoTrackCalorimetry(*daughterTrack, evt, fTrackerTag, fCalorimetryTag);
1641  if(daughtercalovector.size() != 3)
1642  std::cerr << "WARNING::Calorimetry vector size for daughter is = " << daughtercalovector.size() << std::endl;
1643 
1644  for(size_t k = 0; k < daughtercalovector.size() && k<3; k++){
1645  fdaughterKineticEnergy[fNDAUGHTERS][k] = daughtercalovector[k].KineticEnergy();
1646  fdaughterRange[fNDAUGHTERS][k] = daughtercalovector[k].Range();
1647  }
1648 
1649  // Get the true mc particle
1650  const simb::MCParticle* mcdaughterparticle = truthUtil.GetMCParticleFromRecoTrack(clockData, *daughterTrack, evt, fTrackerTag);
1651  if(mcdaughterparticle != 0x0){
1652  fdaughter_truth_TrackId[fNDAUGHTERS] = mcdaughterparticle->TrackId();
1653  fdaughter_truth_Pdg[fNDAUGHTERS] = mcdaughterparticle->PdgCode();
1654  fdaughter_truth_StartPosition[fNDAUGHTERS][0] = mcdaughterparticle->Vx();
1655  fdaughter_truth_StartPosition[fNDAUGHTERS][1] = mcdaughterparticle->Vy();
1656  fdaughter_truth_StartPosition[fNDAUGHTERS][2] = mcdaughterparticle->Vz();
1657  fdaughter_truth_StartPosition[fNDAUGHTERS][3] = mcdaughterparticle->T();
1658  fdaughter_truth_EndPosition[fNDAUGHTERS][0] = mcdaughterparticle->EndX();
1659  fdaughter_truth_EndPosition[fNDAUGHTERS][1] = mcdaughterparticle->EndY();
1660  fdaughter_truth_EndPosition[fNDAUGHTERS][2] = mcdaughterparticle->EndZ();
1661  fdaughter_truth_EndPosition[fNDAUGHTERS][3] = mcdaughterparticle->EndT();
1662  fdaughter_truth_P[fNDAUGHTERS] = mcdaughterparticle->P();
1663  fdaughter_truth_Momentum[fNDAUGHTERS][0] = mcdaughterparticle->Px();
1664  fdaughter_truth_Momentum[fNDAUGHTERS][1] = mcdaughterparticle->Py();
1665  fdaughter_truth_Momentum[fNDAUGHTERS][2] = mcdaughterparticle->Pz();
1666  fdaughter_truth_Momentum[fNDAUGHTERS][3] = mcdaughterparticle->E();
1667  fdaughter_truth_Pt[fNDAUGHTERS] = mcdaughterparticle->Pt();
1668  fdaughter_truth_Mass[fNDAUGHTERS] = mcdaughterparticle->Mass();
1669  fdaughter_truth_EndMomentum[fNDAUGHTERS][0] = mcdaughterparticle->EndPx();
1670  fdaughter_truth_EndMomentum[fNDAUGHTERS][1] = mcdaughterparticle->EndPy();
1671  fdaughter_truth_EndMomentum[fNDAUGHTERS][2] = mcdaughterparticle->EndPz();
1672  fdaughter_truth_EndMomentum[fNDAUGHTERS][3] = mcdaughterparticle->EndE();
1673  fdaughter_truth_Theta[fNDAUGHTERS] = mcdaughterparticle->Momentum().Theta();
1674  fdaughter_truth_Phi[fNDAUGHTERS] = mcdaughterparticle->Momentum().Phi();
1675  fdaughter_truth_Process[fNDAUGHTERS] = int(mcdaughterparticle->Trajectory().ProcessToKey(mcdaughterparticle->Process()));
1676  std::cout << "Daughter Process = " << (mcdaughterparticle->Process()).c_str()
1677  << " , mother = " << mcdaughterparticle->Mother()
1678  << std::endl;
1679  }
1680  }
1681  else if(daughterShower != 0x0){
1682  fisdaughtertrack[fNDAUGHTERS] = 0;
1683  fisdaughtershower[fNDAUGHTERS] = 1;
1684  fdaughterLength[fNDAUGHTERS] = daughterShower->Length();
1685  fdaughterShowerBestPlane[fNDAUGHTERS] = daughterShower->best_plane();
1686  fdaughterOpeningAngle[fNDAUGHTERS] = daughterShower->OpenAngle();
1687  fdaughterStartPosition[fNDAUGHTERS][0] = daughterShower->ShowerStart().X();
1688  fdaughterStartPosition[fNDAUGHTERS][1] = daughterShower->ShowerStart().Y();
1689  fdaughterStartPosition[fNDAUGHTERS][2] = daughterShower->ShowerStart().Z();
1690  fdaughterStartDirection[fNDAUGHTERS][0] = daughterShower->Direction().X();
1691  fdaughterStartDirection[fNDAUGHTERS][1] = daughterShower->Direction().Y();
1692  fdaughterStartDirection[fNDAUGHTERS][2] = daughterShower->Direction().Z();
1693  if( (daughterShower->Energy()).size() > 0 )
1694  fdaughterShowerEnergy[fNDAUGHTERS] = daughterShower->Energy()[0]; // thisShower->best_plane()
1695  if( (daughterShower->MIPEnergy()).size() > 0 )
1696  fdaughterShowerMIPEnergy[fNDAUGHTERS] = daughterShower->MIPEnergy()[0];
1697  if( (daughterShower->dEdx()).size() > 0 )
1698  fdaughterShowerdEdx[fNDAUGHTERS] = daughterShower->dEdx()[0];
1699  }
1700  else{
1701  std::cout << "INFO::Daughter pfParticle is not track or shower. Skip!" << std::endl;
1702  continue;
1703  }
1704 
1705  fdaughterID[fNDAUGHTERS] = daughterID;
1706  // NHits associated with this pfParticle
1707  fdaughterNHits[fNDAUGHTERS] = (pfpUtil.GetPFParticleHits(*daughterParticle,evt,fPFParticleTag)).size();
1708  // T0
1709  //std::vector<anab::T0> pfdaughterT0vec = pfpUtil.GetPFParticleT0(*daughterParticle,evt,fPFParticleTag);
1710  //if(!pfT0vec.empty())
1711  // fdaughterT0[fNDAUGHTERS] = pfdaughterT0vec[0].Time();
1712 
1713  fNDAUGHTERS++;
1714 
1715  // Only process NMAXDAUGTHERS
1716  if(fNDAUGHTERS > NMAXDAUGTHERS) break;
1717 
1718  }
1719 
1720  // For actually studying the objects, it is easier to have the daughters in their track and shower forms.
1721  // We can use the utility to get a vector of track-like and a vector of shower-like daughters
1722  //const std::vector<const recob::Track*> trackDaughters = pfpUtil.GetPFParticleDaughterTracks(*particle,evt,fPFParticleTag,fTrackerTag);
1723  //const std::vector<const recob::Shower*> showerDaughters = pfpUtil.GetPFParticleDaughterShowers(*particle,evt,fPFParticleTag,fShowerTag);
1724 
1725  // For now only consider the first primary track. Need a proper treatment if more than one primary particles are found.
1726  break;
1727  }
1728 
1729  // Fill trees
1730  if(beamTriggerEvent)
1731  fPandoraBeam->Fill();
1732 
1733  //fPandoraCosmics->Fill();
1734 
1735 }
1736 
1737 void protoana::protonmccnn::endJob(){
1738 
1739 }
1740 
1741 void protoana::protonmccnn::FillCosmicsTree(art::Event const & evt, std::string pfParticleTag){
1742 
1743  // To fill
1744 
1745 }
1746 
1747 void protoana::protonmccnn::Initialise(){
1748 
1749  fRun = -999;
1750  fSubRun = -999;
1751  fevent = -999;
1752  fTimeStamp = -999.0;
1753  for(int k=0; k < 5; k++)
1754  fNactivefembs[k] = -999;
1755 
1756  for(int k=0; k < 3; k++){
1757  fvertex[k] = -999.0;
1758  fsecvertex[k] = -999.0;
1759  fprimaryEndPosition[k] = -999.0;
1760  fprimaryStartPosition[k] = -999.0;
1761  fprimaryEndDirection[k] = -999.0;
1762  fprimaryStartDirection[k] = -999.0;
1763  fprimaryKineticEnergy[k] = -999.0;
1764  fprimaryRange[k] = -999.0;
1765 
1766  //xyz_inelscore_c[k] =-999.;
1767  //xyz_elscore_c[k] =-999.;
1768 
1769  }
1770 
1771  fbeamtrigger = -999;
1772  ftof = -999.0;
1773  fcerenkov1 = -999;
1774  fcerenkov2 = -999;
1775  fbeamtrackMomentum = -999.0;
1776  fbeamtrackEnergy = 999.0;
1777  fbeamtrackPdg = -999;
1778  fbeamtrackTime = -999.0;
1779  fbeamtrackID = -999;
1780  for(int l=0; l < 3; l++){
1781  fbeamtrackP[l] = -999.0;
1782  fbeamtrackPos[l] = -999.0;
1783  fbeamtrackDir[l] = -999.0;
1784  }
1785 
1786  //NumberBeamTrajectoryPoints=0;
1787  beamtrk_x.clear();
1788  beamtrk_y.clear();
1789  beamtrk_z.clear();
1790  beamtrk_Px.clear();
1791  beamtrk_Py.clear();
1792  beamtrk_Pz.clear();
1793  beamtrk_Eng.clear();
1794 
1795  x_c.clear();
1796  y_c.clear();
1797  z_c.clear();
1798 
1799 
1800  ke_ff=0.;
1801  pos_ff.clear();
1802 
1803  primtrk_ke_true.clear();
1804  primtrk_ke_reco.clear();
1805  primtrk_end_g4process="";
1806  primtrk_range_true.clear();
1807 
1808  fisprimarytrack = 0;
1809  fisprimaryshower = 0;
1810  fNDAUGHTERS = 0;
1811 
1812  fprimaryBDTScore = -999.0;
1813  fprimaryNHits = -999;
1814  fprimaryTheta = -999.0;
1815  fprimaryPhi = -999.0;
1816  fprimaryLength = -999.0;
1817  fprimaryMomentum = -999.0;
1818  fprimaryEndMomentum = -999.0;
1819  fprimaryOpeningAngle = -999.0;
1820  fprimaryShowerBestPlane = -999;
1821  fprimaryShowerEnergy = -999;
1822  fprimaryShowerMIPEnergy = -999;
1823  fprimaryShowerdEdx = -999;
1824  fprimaryID = -999;
1825  fprimaryMomentumByRangeProton = -999.0;
1826  fprimaryMomentumByRangeMuon = -999.0;
1827  fprimaryT0 = -999.0;
1828 
1829  fprimary_truth_TrackId = -999;
1830  fprimary_truth_Pdg = -999;
1831  fprimary_truth_P = -999.0;
1832  fprimary_truth_Pt = -999.0;
1833  fprimary_truth_Mass = -999.0;
1834  fprimary_truth_Theta = -999.0;
1835  fprimary_truth_Phi = -999.0;
1836  //fprimary_truth_Process = -999;
1837  fprimary_truth_Isbeammatched = -999;
1838  fprimary_truth_NDaughters = -999;
1839 
1840  //fprimary_truth_Process="";
1841  fprimary_truth_EndProcess="";
1842  //fprimary_backtrker_truth_Process="";
1843  //fprimary_backtrker_truth_EndProcess="";
1844 
1845  for(int l=0; l < 4; l++){
1846  fprimary_truth_StartPosition[l] = -999.0;
1847  fprimary_truth_StartPosition_MC[l] = -999.0;
1848  fprimary_truth_EndPosition[l] = -999.0;
1849  fprimary_truth_EndPosition_MC[l] = -999.0;
1850  fprimary_truth_Momentum[l] = -999.0;
1851  fprimary_truth_EndMomentum[l] = -999.0;
1852  }
1853 
1854  for(int k=0; k < NMAXDAUGTHERS; k++){
1855  fisdaughtertrack[k] = -999;
1856  fisdaughtershower[k] = -999;
1857  fdaughterNHits[k] = -999;
1858  fdaughterTheta[k] = -999.0;
1859  fdaughterPhi[k] = -999.0;
1860  fdaughterLength[k] = -999.0;
1861  fdaughterMomentum[k] = -999.0;
1862  fdaughterEndMomentum[k] = -999.0;
1863  for(int l=0; l < 3; l++){
1864  fdaughterEndPosition[k][l] = -999.0;
1865  fdaughterStartPosition[k][l] = -999.0;
1866  fdaughterEndDirection[k][l] = -999.0;
1867  fdaughterStartDirection[k][l] = -999.0;
1868  fdaughterKineticEnergy[k][l] = -999.0;
1869  fdaughterRange[k][l] = -999.0;
1870  }
1871  fdaughterOpeningAngle[k] = -999.0;
1872  fdaughterShowerBestPlane[k] = -999;
1873  fdaughterShowerEnergy[k] = -999;
1874  fdaughterShowerMIPEnergy[k] = -999;
1875  fdaughterShowerdEdx[k] = -999;
1876  fdaughterMomentumByRangeProton[k] = -999.0;
1877  fdaughterMomentumByRangeMuon[k] = -999.0;
1878  fdaughterID[k] = -999;
1879  //fdaughterT0[k] = -999;
1880 
1881  fdaughter_truth_TrackId[k] = -999;
1882  fdaughter_truth_Pdg[k] = -999;
1883  fdaughter_truth_P[k] = -999.0;
1884  fdaughter_truth_Pt[k] = -999.0;
1885  fdaughter_truth_Mass[k] = -999.0;
1886  fdaughter_truth_Theta[k] = -999.0;
1887  fdaughter_truth_Phi[k] = -999.0;
1888  fdaughter_truth_Process[k] = -999;
1889  for(int l=0; l < 4; l++){
1890  fdaughter_truth_StartPosition[k][l] = -999.0;
1891  fdaughter_truth_EndPosition[k][l] = -999.0;
1892  fdaughter_truth_Momentum[k][l] = -999.0;
1893  fdaughter_truth_EndMomentum[k][l] = -999.0;
1894  }
1895  }
1896 
1897  primtrk_dqdx.clear();
1898  primtrk_resrange.clear();
1899  primtrk_dedx.clear();
1900  primtrk_range.clear();
1901  primtrk_hitx.clear();
1902  primtrk_hity.clear();
1903  primtrk_hitz.clear();
1904  primtrk_pitch.clear();
1905 
1906  inelscore_c.clear();
1907  elscore_c.clear();
1908  nonescore_c.clear();
1909 
1910 
1911  ke_ff_true=0.;
1912  //primtrk_hit_processkey.clear();
1913  primtrk_hitx_true.clear();
1914  primtrk_hity_true.clear();
1915  primtrk_hitz_true.clear();
1916  primtrk_trkid_true.clear();
1917  primtrk_edept_true.clear();
1918 
1919  primtrk_true_x.clear();
1920  primtrk_true_y.clear();
1921  primtrk_true_z.clear();
1922  primtrk_true_trkid.clear();
1923  primtrk_true_edept.clear();
1924 
1925  primtrj_true_x.clear();
1926  primtrj_true_y.clear();
1927  primtrj_true_z.clear();
1928  primtrj_true_edept.clear();
1929 
1930  interactionX.clear();
1931  interactionY.clear();
1932  interactionZ.clear();
1933  interactionE.clear();
1934  interactionProcesslist.clear();
1935  interactionAngles.clear();
1936 
1937  Zintersection.clear();
1938  Yintersection.clear();
1939  Xintersection.clear();
1940  timeintersection.clear();
1941 
1942  pfphit_peaktime_c.clear();
1943  pfphit_peaktime_u.clear();
1944  pfphit_peaktime_v.clear();
1945 
1946  pfphit_wireid_c.clear();
1947  pfphit_wireid_u.clear();
1948  pfphit_wireid_v.clear();
1949 
1950  interaction_wid_c.clear();
1951  interaction_tt_c.clear();
1952  interaction_wid_v.clear();
1953  interaction_tt_v.clear();
1954  interaction_wid_u.clear();
1955  interaction_tt_u.clear();
1956 
1957  wid_c.clear();
1958  tt_c.clear();
1959  wid_v.clear();
1960  tt_v.clear();
1961  wid_u.clear();
1962  tt_u.clear();
1963 
1964 }
1965 
1966 DEFINE_ART_MODULE(protoana::protonmccnn)
simb::MCParticle::E
double E(const int i=0) const
Definition: MCParticle.h:233
geo::GeometryCore::WireCoordinate
geo::Length_t WireCoordinate(double YPos, double ZPos, geo::PlaneID const &planeid) const
Returns the index of the nearest wire to the specified position.
Definition: GeometryCore.cxx:1128
protoana::protonmccnn::wid_v
std::vector< double > wid_v
Definition: protonmccnn_module.cc:294
recob::Shower::best_plane
int best_plane() const
Definition: Shower.h:200
protoana::protonmccnn::inelscore_c
std::vector< float > inelscore_c
Definition: protonmccnn_module.cc:277
protoana::protonmccnn::fdaughterNHits
int fdaughterNHits[NMAXDAUGTHERS]
Definition: protonmccnn_module.cc:323
protoana::protonmccnn::primtrj_true_edept
std::vector< std::vector< double > > primtrj_true_edept
Definition: protonmccnn_module.cc:200
NMAXDAUGTHERS
const int NMAXDAUGTHERS
Definition: protonmccnn_module.cc:75
simb::MCParticle::NumberTrajectoryPoints
unsigned int NumberTrajectoryPoints() const
Definition: MCParticle.h:218
protoana::protonmccnn::primtrk_true_y
std::vector< std::vector< double > > primtrk_true_y
Definition: protonmccnn_module.cc:191
protoana::protonmccnn::pfphit_wireid_u
std::vector< double > pfphit_wireid_u
Definition: protonmccnn_module.cc:207
recob::Shower::ShowerStart
const TVector3 & ShowerStart() const
Definition: Shower.h:192
simb::MCTrajectory::Z
double Z(const size_type i) const
Definition: MCTrajectory.h:151
simb::MCTrajectory::X
double X(const size_type i) const
Definition: MCTrajectory.h:149
recob::Track::EndMomentum
double EndMomentum() const
Definition: Track.h:144
simb::MCParticle::PdgCode
int PdgCode() const
Definition: MCParticle.h:212
protoana::protonmccnn::fdaughter_truth_EndPosition
double fdaughter_truth_EndPosition[NMAXDAUGTHERS][4]
Definition: protonmccnn_module.cc:348
art::Timestamp::timeLow
constexpr std::uint32_t timeLow() const
Definition: Timestamp.h:29
protoana::protonmccnn::primtrj_true_x
std::vector< std::vector< double > > primtrj_true_x
Definition: protonmccnn_module.cc:197
protoana::protonmccnn::fisdaughtertrack
int fisdaughtertrack[NMAXDAUGTHERS]
Definition: protonmccnn_module.cc:321
protoana::protonmccnn::beamtrk_Eng
std::vector< double > beamtrk_Eng
Definition: protonmccnn_module.cc:162
protoana::protonmccnn::Xintersection
std::vector< double > Xintersection
Definition: protonmccnn_module.cc:309
protoana::protonmccnn::pfphit_wireid_v
std::vector< double > pfphit_wireid_v
Definition: protonmccnn_module.cc:208
protoana::protonmccnn::interactionAngles
std::vector< double > interactionAngles
Definition: protonmccnn_module.cc:306
simb::MCParticle::Py
double Py(const int i=0) const
Definition: MCParticle.h:231
protoana::protonmccnn::fdaughter_truth_StartPosition
double fdaughter_truth_StartPosition[NMAXDAUGTHERS][4]
Definition: protonmccnn_module.cc:347
protoana::ProtoDUNEPFParticleUtils::GetPFParticleT0
std::vector< anab::T0 > GetPFParticleT0(const recob::PFParticle &particle, art::Event const &evt, std::string particleLabel) const
Get the T0(s) from a given PFParticle.
Definition: ProtoDUNEPFParticleUtils.cxx:101
recob::PFParticle::NumDaughters
int NumDaughters() const
Returns the number of daughter particles flowing from this one.
Definition: PFParticle.h:89
protoana::protonmccnn::beamtrk_y
std::vector< double > beamtrk_y
Definition: protonmccnn_module.cc:157
protoana::ProtoDUNETruthUtils::GetGeantGoodParticle
const simb::MCParticle * GetGeantGoodParticle(const simb::MCTruth &genTruth, const art::Event &evt) const
Definition: ProtoDUNETruthUtils.cxx:860
simb::MCParticle::EndE
double EndE() const
Definition: MCParticle.h:244
simb::MCParticle::EndZ
double EndZ() const
Definition: MCParticle.h:228
protoana::protonmccnn::fGeometry
geo::GeometryCore const * fGeometry
Definition: protonmccnn_module.cc:127
protoana::protonmccnn::x_c
std::vector< float > x_c
Definition: protonmccnn_module.cc:281
simb::MCTrajectory::E
double E(const size_type i) const
Definition: MCTrajectory.h:156
protoana::protonmccnn::interaction_wid_u
std::vector< double > interaction_wid_u
Definition: protonmccnn_module.cc:316
xpos
TH3F * xpos
Definition: doAna.cpp:29
protoana::protonmccnn::fdaughter_truth_EndMomentum
double fdaughter_truth_EndMomentum[NMAXDAUGTHERS][4]
Definition: protonmccnn_module.cc:351
recob::Shower::Length
double Length() const
Definition: Shower.h:201
protoana::protonmccnn::operator=
protonmccnn & operator=(protonmccnn const &)=delete
protoana::ProtoDUNEPFParticleUtils::GetPFParticleShower
const recob::Shower * GetPFParticleShower(const recob::PFParticle &particle, art::Event const &evt, const std::string particleLabel, const std::string showerLabel) const
Get the shower associated to this particle. Returns a null pointer if not found.
Definition: ProtoDUNEPFParticleUtils.cxx:496
art::DataViewImpl::getHandle
Handle< PROD > getHandle(SelectorBase const &) const
Definition: DataViewImpl.h:382
geo::View_t
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
string
std::string string
Definition: nybbler.cc:12
simb::MCParticle::Trajectory
const simb::MCTrajectory & Trajectory() const
Definition: MCParticle.h:253
protoana::protonmccnn::timeintersection
std::vector< double > timeintersection
Definition: protonmccnn_module.cc:310
simb::MCParticle::Mother
int Mother() const
Definition: MCParticle.h:213
protoana::protonmccnn::beamtrk_z
std::vector< double > beamtrk_z
Definition: protonmccnn_module.cc:158
recob::Track::Trajectory
const recob::TrackTrajectory & Trajectory() const
Access to the stored recob::TrackTrajectory.
Definition: Track.h:98
protoana::ProtoDUNETrackUtils::GetRecoTrackCalorimetry
std::vector< anab::Calorimetry > GetRecoTrackCalorimetry(const recob::Track &track, art::Event const &evt, const std::string trackModule, const std::string caloModule) const
Get the Calorimetry(s) from a given reco track.
Definition: ProtoDUNETrackUtils.cxx:72
protoana::protonmccnn::interactionY
std::vector< double > interactionY
Definition: protonmccnn_module.cc:302
protoana::protonmccnn::fdaughter_truth_Momentum
double fdaughter_truth_Momentum[NMAXDAUGTHERS][4]
Definition: protonmccnn_module.cc:350
simb::MCTrajectory::at
const value_type & at(const size_type i) const
Definition: MCTrajectory.h:175
art::Timestamp::timeHigh
constexpr std::uint32_t timeHigh() const
Definition: Timestamp.h:34
simb::MCParticle::Mass
double Mass() const
Definition: MCParticle.h:239
cet::pow
constexpr T pow(T x)
Definition: pow.h:72
simb::MCTrajectory::KeyToProcess
std::string KeyToProcess(unsigned char const &key) const
Definition: MCTrajectory.cxx:117
protoana::protonmccnn::beginJob
virtual void beginJob() override
Definition: protonmccnn_module.cc:384
protoana::protonmccnn::FillCosmicsTree
void FillCosmicsTree(art::Event const &evt, std::string pfParticleTag)
Definition: protonmccnn_module.cc:1741
protoana::ProtoDUNEPFParticleUtils::GetPFParticleHits
const std::vector< const recob::Hit * > GetPFParticleHits(const recob::PFParticle &particle, art::Event const &evt, const std::string particleLabel) const
Get the hits associated to the PFParticle.
Definition: ProtoDUNEPFParticleUtils.cxx:564
recob::Shower::Energy
const std::vector< double > & Energy() const
Definition: Shower.h:195
simb::MCParticle::Px
double Px(const int i=0) const
Definition: MCParticle.h:230
protoana::ProtoDUNEPFParticleUtils::GetPFParticleVertex
const TVector3 GetPFParticleVertex(const recob::PFParticle &particle, art::Event const &evt, const std::string particleLabel, const std::string trackLabel) const
Function to find the interaction vertex of a primary PFParticle.
Definition: ProtoDUNEPFParticleUtils.cxx:352
TrackHitMeta.h
Class to keep data related to recob::Hit associated with recob::Track.
protoana::protonmccnn::interaction_tt_u
std::vector< double > interaction_tt_u
Definition: protonmccnn_module.cc:317
recob::tracking::Plane::Plane
Plane(const Point_t &planePos, const Vector_t &planeDir)
Constructor from reference position on the plane and direction orthogonal to the plane.
Definition: TrackingPlane.h:61
anab::FVectorReader::getIndex
int getIndex(const std::string &name) const
Index of column with given name, or -1 if name not found.
Definition: MVAReader.h:82
cheat::BackTrackerService::TrackIdToSimIDEs_Ps
std::vector< const sim::IDE * > TrackIdToSimIDEs_Ps(int const &id) const
Definition: BackTrackerService_service.cc:134
std
STL namespace.
protoana::protonmccnn::primtrk_ke_true
std::vector< double > primtrk_ke_true
Definition: protonmccnn_module.cc:169
const_iterator
intermediate_table::const_iterator const_iterator
Definition: intermediate_table.cc:28
protoana::protonmccnn::analyze
void analyze(art::Event const &evt) override
Definition: protonmccnn_module.cc:598
protoana::protonmccnn::beamtrk_Py
std::vector< double > beamtrk_Py
Definition: protonmccnn_module.cc:160
protoana::protonmccnn::fdaughter_truth_Process
int fdaughter_truth_Process[NMAXDAUGTHERS]
Definition: protonmccnn_module.cc:356
Index
unsigned int Index
Definition: sspmapmaker_v1.c:14
protoana::protonmccnn::fdaughterLength
double fdaughterLength[NMAXDAUGTHERS]
Definition: protonmccnn_module.cc:326
art::EDAnalyzer::EDAnalyzer
EDAnalyzer(fhicl::ParameterSet const &pset)
Definition: EDAnalyzer.h:25
simb::MCParticle::Process
std::string Process() const
Definition: MCParticle.h:215
protoana::protonmccnn::fdaughterMomentum
double fdaughterMomentum[NMAXDAUGTHERS]
Definition: protonmccnn_module.cc:327
protoana::protonmccnn::primtrk_true_x
std::vector< std::vector< double > > primtrk_true_x
Definition: protonmccnn_module.cc:190
MCParticle.h
Particle class.
simb::MCTrajectory::ProcessToKey
unsigned char ProcessToKey(std::string const &process) const
Definition: MCTrajectory.cxx:97
simb::MCParticle::EndY
double EndY() const
Definition: MCParticle.h:227
l
static QStrList * l
Definition: config.cpp:1044
protoana::protonmccnn::wid_u
std::vector< double > wid_u
Definition: protonmccnn_module.cc:296
protoana::protonmccnn::fprimary_truth_EndProcess
std::string fprimary_truth_EndProcess
Definition: protonmccnn_module.cc:260
protoana::protonmccnn::primtrk_dqdx
std::vector< std::vector< double > > primtrk_dqdx
Definition: protonmccnn_module.cc:265
simb::MCParticle::NumberDaughters
int NumberDaughters() const
Definition: MCParticle.h:217
protoana::protonmccnn::primtrk_hity_true
std::vector< std::vector< double > > primtrk_hity_true
Definition: protonmccnn_module.cc:184
protoana::protonmccnn::primtrj_true_y
std::vector< std::vector< double > > primtrj_true_y
Definition: protonmccnn_module.cc:198
protoana::protonmccnn::endJob
virtual void endJob() override
Definition: protonmccnn_module.cc:1737
ypos
TH3F * ypos
Definition: doAna.cpp:30
simb::MCParticle::TrackId
int TrackId() const
Definition: MCParticle.h:210
protoana::protonmccnn::truthUtil
protoana::ProtoDUNETruthUtils truthUtil
Definition: protonmccnn_module.cc:108
protoana::protonmccnn::tt_u
std::vector< double > tt_u
Definition: protonmccnn_module.cc:297
protoana::protonmccnn::y_c
std::vector< float > y_c
Definition: protonmccnn_module.cc:282
protoana::protonmccnn::wid_c
std::vector< double > wid_c
Definition: protonmccnn_module.cc:292
protoana::protonmccnn::pos_ff
std::vector< double > pos_ff
Definition: protonmccnn_module.cc:166
util::size
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:92
protoana::protonmccnn::fdaughterMomentumByRangeProton
double fdaughterMomentumByRangeProton[NMAXDAUGTHERS]
Definition: protonmccnn_module.cc:338
art::DataViewImpl::isRealData
bool isRealData() const
Definition: DataViewImpl.cc:113
protoana::protonmccnn::dataUtil
protoana::ProtoDUNEDataUtils dataUtil
Definition: protonmccnn_module.cc:105
tt
Definition: type_traits.h:61
simb::MCParticle::Pt
double Pt(const int i=0) const
Definition: MCParticle.h:236
protoana::protonmccnn::fisdaughtershower
int fisdaughtershower[NMAXDAUGTHERS]
Definition: protonmccnn_module.cc:322
zpos
TH3F * zpos
Definition: doAna.cpp:31
recob::Track::Phi
double Phi() const
Definition: Track.h:178
protoana::protonmccnn::fdaughter_truth_P
double fdaughter_truth_P[NMAXDAUGTHERS]
Definition: protonmccnn_module.cc:349
protoana::protonmccnn::trackUtil
protoana::ProtoDUNETrackUtils trackUtil
Definition: protonmccnn_module.cc:107
recob::Track::Length
double Length(size_t p=0) const
Access to various track properties.
Definition: Track.h:167
simb::MCTrajectory::Y
double Y(const size_type i) const
Definition: MCTrajectory.h:150
art::DataViewImpl::time
Timestamp time() const
Definition: DataViewImpl.cc:106
recob::Shower::dEdx
const std::vector< double > & dEdx() const
Definition: Shower.h:203
DEFINE_ART_MODULE
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:67
recob::Shower::MIPEnergy
const std::vector< double > & MIPEnergy() const
Definition: Shower.h:198
protoana::protonmccnn::fdaughter_truth_Mass
double fdaughter_truth_Mass[NMAXDAUGTHERS]
Definition: protonmccnn_module.cc:353
wirecell.dfp.graph.key
def key(type, name=None)
Definition: graph.py:13
simb::MCParticle::EndPz
double EndPz() const
Definition: MCParticle.h:243
protoana::protonmccnn::pfphit_wireid_c
std::vector< double > pfphit_wireid_c
Definition: protonmccnn_module.cc:206
recob::Track::Theta
double Theta() const
Access to spherical or geographical angles at vertex or at any point.
Definition: Track.h:176
protoana::protonmccnn::fdaughterKineticEnergy
double fdaughterKineticEnergy[NMAXDAUGTHERS][3]
Definition: protonmccnn_module.cc:340
simb::MCParticle::P
double P(const int i=0) const
Definition: MCParticle.h:234
protoana::ProtoDUNEPFParticleUtils::GetPFParticleTrack
const recob::Track * GetPFParticleTrack(const recob::PFParticle &particle, art::Event const &evt, const std::string particleLabel, const std::string trackLabel) const
Get the track associated to this particle. Returns a null pointer if not found.
Definition: ProtoDUNEPFParticleUtils.cxx:476
recob::Shower::OpenAngle
double OpenAngle() const
Definition: Shower.h:202
simb::MCParticle::T
double T(const int i=0) const
Definition: MCParticle.h:224
art::DataViewImpl::getValidHandle
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
Definition: DataViewImpl.h:441
test.p
p
Definition: test.py:223
protoana::protonmccnn::beamtrk_x
std::vector< double > beamtrk_x
Definition: protonmccnn_module.cc:156
protoana::protonmccnn::primtrk_hity
std::vector< std::vector< double > > primtrk_hity
Definition: protonmccnn_module.cc:270
art::DataViewImpl::subRun
SubRunNumber_t subRun() const
Definition: DataViewImpl.cc:78
protoana::protonmccnn::interaction_tt_c
std::vector< double > interaction_tt_c
Definition: protonmccnn_module.cc:313
protoana::protonmccnn::fdaughter_truth_Pt
double fdaughter_truth_Pt[NMAXDAUGTHERS]
Definition: protonmccnn_module.cc:352
protoana::protonmccnn::primtrk_hitz_true
std::vector< std::vector< double > > primtrk_hitz_true
Definition: protonmccnn_module.cc:185
protoana::ProtoDUNETruthUtils::GetMCParticleFromRecoTrack
const simb::MCParticle * GetMCParticleFromRecoTrack(detinfo::DetectorClocksData const &clockData, const recob::Track &track, art::Event const &evt, std::string trackModule) const
Definition: ProtoDUNETruthUtils.cxx:609
protoana::ProtoDUNEPFParticleUtils::GetPFParticlesFromBeamSlice
const std::vector< const recob::PFParticle * > GetPFParticlesFromBeamSlice(art::Event const &evt, const std::string particleLabel) const
Return the pointers for the PFParticles in the beam slice. Returns an empty vector is no beam slice w...
Definition: ProtoDUNEPFParticleUtils.cxx:175
protoana::protonmccnn::fdaughterPhi
double fdaughterPhi[NMAXDAUGTHERS]
Definition: protonmccnn_module.cc:325
recob::Shower::Direction
const TVector3 & Direction() const
Definition: Shower.h:189
protoana::protonmccnn::fdaughter_truth_Phi
double fdaughter_truth_Phi[NMAXDAUGTHERS]
Definition: protonmccnn_module.cc:355
protoana::protonmccnn::interaction_wid_c
std::vector< double > interaction_wid_c
Definition: protonmccnn_module.cc:312
art::DataViewImpl::run
RunNumber_t run() const
Definition: DataViewImpl.cc:71
max
static int max(int a, int b)
Definition: fortranscanner.cpp:60366
protoana::protonmccnn::fdaughterOpeningAngle
double fdaughterOpeningAngle[NMAXDAUGTHERS]
Definition: protonmccnn_module.cc:333
simb::MCParticle::EndT
double EndT() const
Definition: MCParticle.h:229
geo::GeometryCore
Description of geometry of one entire detector.
Definition: GeometryCore.h:1468
simb::MCTrajectory::TrajectoryProcesses
ProcessMap const & TrajectoryProcesses() const
Definition: MCTrajectory.h:188
protoana::protonmccnn::interaction_tt_v
std::vector< double > interaction_tt_v
Definition: protonmccnn_module.cc:315
protoana::protonmccnn::fprimary_truth_StartPosition_MC
double fprimary_truth_StartPosition_MC[4]
Definition: protonmccnn_module.cc:245
geo::GeometryCore::View
View_t View(geo::PlaneID const &pid) const
Returns the view (wire orientation) on the channels of specified TPC plane.
Definition: GeometryCore.cxx:601
protoana::protonmccnn::fdaughterTheta
double fdaughterTheta[NMAXDAUGTHERS]
Definition: protonmccnn_module.cc:324
protoana::protonmccnn::fdaughterShowerdEdx
double fdaughterShowerdEdx[NMAXDAUGTHERS]
Definition: protonmccnn_module.cc:336
protoana::protonmccnn::z_c
std::vector< float > z_c
Definition: protonmccnn_module.cc:283
protoana::protonmccnn::interactionZ
std::vector< double > interactionZ
Definition: protonmccnn_module.cc:303
protoana::protonmccnn::primtrk_range
std::vector< double > primtrk_range
Definition: protonmccnn_module.cc:268
protoana::protonmccnn::fdaughterID
int fdaughterID[NMAXDAUGTHERS]
Definition: protonmccnn_module.cc:342
cheat::ParticleInventoryService::ParticleList
const sim::ParticleList & ParticleList() const
Definition: ParticleInventoryService_service.cc:114
protoana::protonmccnn::fdaughterStartPosition
double fdaughterStartPosition[NMAXDAUGTHERS][3]
Definition: protonmccnn_module.cc:330
keras_to_tensorflow.int
int
Definition: keras_to_tensorflow.py:69
protoana::protonmccnn::pfphit_peaktime_c
std::vector< double > pfphit_peaktime_c
Definition: protonmccnn_module.cc:203
recob::TrackTrajectory::EndDirection
Vector_t EndDirection() const
Returns the direction of the trajectory at the last point.
Definition: TrackTrajectory.h:321
protoana::protonmccnn::beamtrk_Pz
std::vector< double > beamtrk_Pz
Definition: protonmccnn_module.cc:161
tracks
Definition: tracks.py:1
simb::MCParticle::Vx
double Vx(const int i=0) const
Definition: MCParticle.h:221
protoana::protonmccnn::pfpUtil
protoana::ProtoDUNEPFParticleUtils pfpUtil
Definition: protonmccnn_module.cc:106
recob::Track::ID
int ID() const
Definition: Track.h:198
protoana::protonmccnn::primtrk_true_edept
std::vector< std::vector< double > > primtrk_true_edept
Definition: protonmccnn_module.cc:194
protoana::protonmccnn::fBeamModuleLabel
const art::InputTag fBeamModuleLabel
Definition: protonmccnn_module.cc:118
Hit.h
Declaration of signal hit object.
protoana::protonmccnn::primtrk_hitz
std::vector< std::vector< double > > primtrk_hitz
Definition: protonmccnn_module.cc:271
recob::Track::StartMomentum
double StartMomentum() const
Definition: Track.h:143
protoana::protonmccnn::Yintersection
std::vector< double > Yintersection
Definition: protonmccnn_module.cc:308
recob::PFParticle
Hierarchical representation of particle flow.
Definition: PFParticle.h:44
protoana::protonmccnn::primtrk_dedx
std::vector< std::vector< double > > primtrk_dedx
Definition: protonmccnn_module.cc:267
simb::MCTrajectory::size
size_type size() const
Definition: MCTrajectory.h:166
simb::MCParticle::EndPy
double EndPy() const
Definition: MCParticle.h:242
protoana::protonmccnn::fdaughter_truth_TrackId
int fdaughter_truth_TrackId[NMAXDAUGTHERS]
Definition: protonmccnn_module.cc:345
protoana::ProtoDUNEDataUtils::GetNActiveFembsForAPA
int GetNActiveFembsForAPA(art::Event const &evt, int apa) const
Get number of active fembs in an APA.
Definition: ProtoDUNEDataUtils.cxx:49
reverse
static unsigned int reverse(QString &chars, unsigned char *level, unsigned int a, unsigned int b)
Definition: qstring.cpp:11649
protoana::protonmccnn::fdaughter_truth_Theta
double fdaughter_truth_Theta[NMAXDAUGTHERS]
Definition: protonmccnn_module.cc:354
protoana::protonmccnn::fdaughterShowerEnergy
double fdaughterShowerEnergy[NMAXDAUGTHERS]
Definition: protonmccnn_module.cc:334
recob::TrackTrajectory::End
Point_t const & End() const
Returns the position of the last valid point of the trajectory [cm].
Definition: TrackTrajectory.h:257
protoana::protonmccnn::elscore_c
std::vector< float > elscore_c
Definition: protonmccnn_module.cc:278
protoana::protonmccnn::primtrk_true_z
std::vector< std::vector< double > > primtrk_true_z
Definition: protonmccnn_module.cc:192
simb::MCParticle::Momentum
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:220
simb::MCParticle::Pz
double Pz(const int i=0) const
Definition: MCParticle.h:232
Track.h
Provides recob::Track data product.
protoana::protonmccnn::fdaughterEndDirection
double fdaughterEndDirection[NMAXDAUGTHERS][3]
Definition: protonmccnn_module.cc:331
simb::MCParticle::Vz
double Vz(const int i=0) const
Definition: MCParticle.h:223
protoana::protonmccnn::fdaughter_truth_Pdg
int fdaughter_truth_Pdg[NMAXDAUGTHERS]
Definition: protonmccnn_module.cc:346
protoana::protonmccnn::fdaughterShowerMIPEnergy
double fdaughterShowerMIPEnergy[NMAXDAUGTHERS]
Definition: protonmccnn_module.cc:335
geo::GeometryCore::WireEndPoints
void WireEndPoints(geo::WireID const &wireid, double *xyzStart, double *xyzEnd) const
Fills two arrays with the coordinates of the wire end points.
Definition: GeometryCore.cxx:1243
protoana::ProtoDUNEPFParticleUtils::GetPFParticleSecondaryVertex
const TVector3 GetPFParticleSecondaryVertex(const recob::PFParticle &particle, art::Event const &evt, const std::string particleLabel, const std::string trackLabel) const
Function to find the secondary interaction vertex of a primary PFParticle.
Definition: ProtoDUNEPFParticleUtils.cxx:407
art::EventID::event
EventNumber_t event() const
Definition: EventID.h:116
protoana::protonmccnn::beamtrk_Px
std::vector< double > beamtrk_Px
Definition: protonmccnn_module.cc:159
protoana::protonmccnn::interactionE
std::vector< double > interactionE
Definition: protonmccnn_module.cc:304
protoana::protonmccnn::fdaughterRange
double fdaughterRange[NMAXDAUGTHERS][3]
Definition: protonmccnn_module.cc:341
simb::MCParticle::EndPx
double EndPx() const
Definition: MCParticle.h:241
protoana::protonmccnn::primtrk_resrange
std::vector< std::vector< double > > primtrk_resrange
Definition: protonmccnn_module.cc:266
protoana::protonmccnn::pfphit_peaktime_u
std::vector< double > pfphit_peaktime_u
Definition: protonmccnn_module.cc:204
protoana::protonmccnn::primtrk_edept_true
std::vector< std::vector< double > > primtrk_edept_true
Definition: protonmccnn_module.cc:187
recob::Hit
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:48
protoana::protonmccnn::primtrk_trkid_true
std::vector< std::vector< double > > primtrk_trkid_true
Definition: protonmccnn_module.cc:186
protoana::protonmccnn::Zintersection
std::vector< double > Zintersection
Definition: protonmccnn_module.cc:307
protoana::protonmccnn::primtrk_true_trkid
std::vector< std::vector< double > > primtrk_true_trkid
Definition: protonmccnn_module.cc:193
tca::evt
TCEvent evt
Definition: DataStructs.cxx:7
protoana::protonmccnn::primtrj_true_z
std::vector< std::vector< double > > primtrj_true_z
Definition: protonmccnn_module.cc:199
protoana::protonmccnn::fdaughterMomentumByRangeMuon
double fdaughterMomentumByRangeMuon[NMAXDAUGTHERS]
Definition: protonmccnn_module.cc:339
protoana::protonmccnn::pfphit_peaktime_v
std::vector< double > pfphit_peaktime_v
Definition: protonmccnn_module.cc:205
art::get
auto const & get(AssnsNode< L, R, D > const &r)
Definition: AssnsNode.h:115
art::fill_ptr_vector
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:297
protoana::protonmccnn::interactionProcesslist
std::vector< std::string > interactionProcesslist
Definition: protonmccnn_module.cc:305
protoana::protonmccnn::nonescore_c
std::vector< float > nonescore_c
Definition: protonmccnn_module.cc:279
protoana::protonmccnn::fNNetModuleLabel
const art::InputTag fNNetModuleLabel
Definition: protonmccnn_module.cc:117
simb::MCParticle::EndX
double EndX() const
Definition: MCParticle.h:226
protoana::ProtoDUNEPFParticleUtils::GetNumberPrimaryPFParticle
unsigned int GetNumberPrimaryPFParticle(art::Event const &evt, const std::string particleLabel) const
Get the number of primary PFParticles.
Definition: ProtoDUNEPFParticleUtils.cxx:24
recob::TrackTrajectory::StartDirection
Vector_t StartDirection() const
Returns the direction of the trajectory at the first point.
Definition: TrackTrajectory.h:317
util::quantities::second
second_as<> second
Type of time stored in seconds, in double precision.
Definition: spacetime.h:85
protoana::protonmccnn::tt_c
std::vector< double > tt_c
Definition: protonmccnn_module.cc:293
detinfo::sampling_rate
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
Definition: DetectorClocksData.h:643
trkf::Plane
recob::tracking::Plane Plane
Definition: TrackState.h:17
protoana::protonmccnn::fdaughterStartDirection
double fdaughterStartDirection[NMAXDAUGTHERS][3]
Definition: protonmccnn_module.cc:332
protoana::protonmccnn::tt_v
std::vector< double > tt_v
Definition: protonmccnn_module.cc:295
recob::TrackTrajectory::Start
Point_t const & Start() const
Returns the position of the first valid point of the trajectory [cm].
Definition: TrackTrajectory.h:253
protoana::protonmccnn::primtrk_ke_reco
std::vector< double > primtrk_ke_reco
Definition: protonmccnn_module.cc:170
bool
int bool
Definition: qglobal.h:345
protoana::protonmccnn::interactionX
std::vector< double > interactionX
Definition: protonmccnn_module.cc:301
anab::MVAReader::getOutput
std::array< float, N > getOutput(size_t key) const
Get copy of the MVA output vector at index "key".
Definition: MVAReader.h:129
recob::Shower::ID
int ID() const
Definition: Shower.h:187
protoana::protonmccnn::interaction_wid_v
std::vector< double > interaction_wid_v
Definition: protonmccnn_module.cc:314
protoana::protonmccnn::primtrk_range_true
std::vector< double > primtrk_range_true
Definition: protonmccnn_module.cc:172
protoana::protonmccnn::fdaughterShowerBestPlane
int fdaughterShowerBestPlane[NMAXDAUGTHERS]
Definition: protonmccnn_module.cc:337
protoana::ProtoDUNEPFParticleUtils::GetBeamCosmicScore
float GetBeamCosmicScore(const recob::PFParticle &particle, art::Event const &evt, const std::string particleLabel) const
Access the BDT output used to decide if a slice is beam-like or cosmic-like.
Definition: ProtoDUNEPFParticleUtils.cxx:182
protoana::protonmccnn::primtrk_hitx
std::vector< std::vector< double > > primtrk_hitx
Definition: protonmccnn_module.cc:269
art::Event::id
EventID id() const
Definition: Event.cc:34
simb::MCParticle::Vy
double Vy(const int i=0) const
Definition: MCParticle.h:222
protoana::protonmccnn::primtrk_hitx_true
std::vector< std::vector< double > > primtrk_hitx_true
Definition: protonmccnn_module.cc:183
recob::Track
Track from a non-cascading particle.A recob::Track consists of a recob::TrackTrajectory, plus additional members relevant for a "fitted" track:
Definition: Track.h:49
calo
calorimetry
Definition: CaloChecker_module.cc:43
trkf::TrackMomentumCalculator::GetTrackMomentum
double GetTrackMomentum(double trkrange, int pdg) const
Definition: TrackMomentumCalculator.cxx:144
fhicl::exception
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
protoana::ProtoDUNEDataUtils::IsBeamTrigger
bool IsBeamTrigger(art::Event const &evt) const
Definition: ProtoDUNEDataUtils.cxx:27
protoana::protonmccnn::primtrk_pitch
std::vector< std::vector< double > > primtrk_pitch
Definition: protonmccnn_module.cc:272
protoana::protonmccnn::protonmccnn
protonmccnn(fhicl::ParameterSet const &p)
Definition: protonmccnn_module.cc:368
protoana::protonmccnn::fdaughterEndMomentum
double fdaughterEndMomentum[NMAXDAUGTHERS]
Definition: protonmccnn_module.cc:328
protoana::protonmccnn::fdaughterEndPosition
double fdaughterEndPosition[NMAXDAUGTHERS][3]
Definition: protonmccnn_module.cc:329
Generated by   doxygen 1.8.11