ProtoDUNEAnalysisTree_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: ProtoDUNEAnalysisTree
3 // File: ProtoDUNEAnalysisTree_module.cc
4 //
5 // Extract protoDUNE useful information, do a first pre-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 //
11 ////////////////////////////////////////////////////////////////////////
12 
13 #include "art/Framework/Core/EDAnalyzer.h"
14 #include "art/Framework/Core/ModuleMacros.h"
15 #include "art/Framework/Principal/Event.h"
16 #include "art/Framework/Principal/Handle.h"
17 #include "art/Framework/Principal/Run.h"
18 #include "art/Framework/Principal/SubRun.h"
19 #include "art_root_io/TFileService.h"
20 #include "canvas/Utilities/InputTag.h"
21 #include "canvas/Persistency/Common/FindManyP.h"
22 #include "fhiclcpp/ParameterSet.h"
23 #include "messagefacility/MessageLogger/MessageLogger.h"
24 
25 #include "larsim/MCCheater/BackTrackerService.h"
26 #include "larsim/MCCheater/ParticleInventoryService.h"
27 #include "larsim/Simulation/LArG4Parameters.h"
28 #include "lardataobj/RecoBase/Hit.h"
29 #include "lardataobj/RecoBase/Track.h"
30 #include "lardataobj/RecoBase/Shower.h"
31 #include "lardataobj/RecoBase/PFParticle.h"
32 #include "lardataobj/AnalysisBase/CosmicTag.h"
33 #include "lardataobj/AnalysisBase/T0.h"
34 #include "lardataobj/AnalysisBase/Calorimetry.h"
35 #include "larreco/RecoAlg/TrackMomentumCalculator.h"
36 #include "larcore/Geometry/Geometry.h"
37 #include "larcorealg/Geometry/GeometryCore.h"
38 #include "larcoreobj/SimpleTypesAndConstants/geo_types.h"
39 #include "lardata/DetectorInfoServices/DetectorClocksService.h"
40 #include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
41 
42 #include "nusimdata/SimulationBase/MCParticle.h"
43 #include "nusimdata/SimulationBase/MCTruth.h"
44 
45 #include "dune/DuneObj/ProtoDUNEBeamEvent.h"
46 #include "protoduneana/Utilities/ProtoDUNETrackUtils.h"
47 #include "protoduneana/Utilities/ProtoDUNEShowerUtils.h"
48 #include "protoduneana/Utilities/ProtoDUNETruthUtils.h"
49 #include "protoduneana/Utilities/ProtoDUNEPFParticleUtils.h"
50 #include "dune/Protodune/singlephase/DataUtils/ProtoDUNEDataUtils.h"
51 #include "protoduneana/Utilities/ProtoDUNEBeamlineUtils.h"
52 
53 // ROOT includes
54 #include "TTree.h"
55 #include "TFile.h"
56 #include "TString.h"
57 #include "TTimeStamp.h"
58 
59 // C++ Includes
60 #include <stdio.h>
61 #include <stdlib.h>
62 #include <fstream>
63 #include <string>
64 #include <sstream>
65 #include <cmath>
66 #include <algorithm>
67 #include <iostream>
68 #include <vector>
69 
70 // Maximum number of beam particles to save
71 const int NMAXDAUGTHERS = 15;
72 const int NMAXTRUTHDAUGTHERS = 30;
73 
74 namespace protoana {
75  class ProtoDUNEAnalysisTree;
76 }
77 
78 // -----------------------------------------------------------------------------
79 class protoana::ProtoDUNEAnalysisTree : public art::EDAnalyzer {
80 public:
81 
82  explicit ProtoDUNEAnalysisTree(fhicl::ParameterSet const & p);
83 
84  ProtoDUNEAnalysisTree(ProtoDUNEAnalysisTree const &) = delete;
85  ProtoDUNEAnalysisTree(ProtoDUNEAnalysisTree &&) = delete;
86  ProtoDUNEAnalysisTree & operator = (ProtoDUNEAnalysisTree const &) = delete;
87  ProtoDUNEAnalysisTree & operator = (ProtoDUNEAnalysisTree &&) = delete;
88 
89  virtual void beginJob() override;
90  virtual void endJob() override;
91 
92  void beginRun(const art::Run& run) override;
93 
94  // Required functions.
95  void analyze(art::Event const & evt) override;
96 
97 private:
98 
99  // Helper utility functions
100  protoana::ProtoDUNEDataUtils dataUtil;
101  protoana::ProtoDUNEPFParticleUtils pfpUtil;
102  protoana::ProtoDUNETrackUtils trackUtil;
103  protoana::ProtoDUNETruthUtils truthUtil;
104  protoana::ProtoDUNEBeamlineUtils beamlineUtil;
105  art::ServiceHandle<cheat::ParticleInventoryService> pi_service;
106 
107  // Track momentum algorithm calculates momentum based on track range
108  trkf::TrackMomentumCalculator trmom;
109 
110  // Geometry
111  geo::GeometryCore const * fGeometry = &*(art::ServiceHandle<geo::Geometry>());
112 
113  // Initialise tree variables
114  void Initialise();
115 
116  // Fill configuration tree
117  void FillConfigTree();
118 
119  bool FillPrimaryBeamParticle(art::Event const & evt);
120  void FillPrimaryPFParticle(art::Event const & evt,
121  detinfo::DetectorClocksData const& clockData,
122  detinfo::DetectorPropertiesData const& detProp,
123  const recob::PFParticle* particle);
124  void FillPrimaryDaughterPFParticle(art::Event const & evt,
125  detinfo::DetectorClocksData const& clockData,
126  detinfo::DetectorPropertiesData const& detProp,
127  const recob::PFParticle* daughterParticle,
128  int daughterID);
129  void FillPrimaryGrandDaughterPFParticle(art::Event const & evt,
130  detinfo::DetectorClocksData const& clockData,
131  detinfo::DetectorPropertiesData const& detProp,
132  const recob::PFParticle* gdaughterParticle,
133  int daughterID, int gdaughterID);
134 
135  // fcl parameters
136  const art::InputTag fBeamModuleLabel;
137  std::string fCalorimetryTag;
138  std::string fParticleIDTag;
139  std::string fTrackerTag;
140  std::string fShowerTag;
141  std::string fPFParticleTag;
142  std::string fGeneratorTag;
143  std::string fSimulationTag;
144  int fVerbose;
145 
146  // Define configuration tree
147  TTree *fConfigTree;
148 
149  int fNAPAs;
150  int fNChansPerAPA;
151  int fNCryostats;
152  int fNTPCs;
153  int fNChannels;
154  int fNPlanes;
155  double fActiveTPCBoundsX[2];
156  double fActiveTPCBoundsY[2];
157  double fActiveTPCBoundsZ[2];
158 
159  // Beam track tree
160  TTree *fBeamTrack;
161 
162  // Primary reco track/showers tree
163  TTree *fPandoraBeam;
164 
165  // Reco daughters tree
166  TTree* fDaughterTree;
167 
168  // Reco granddaughters tree
169  TTree* fGrandDaughterTree;
170 
171  int fRun;
172  int fSubRun;
173  int fevent;
174  double fTimeStamp;
175  int fNactivefembs[6];
176 
177  // Beam track
178  int fbeamNTracks;
179  int fbeamtrigger;
180  double ftof;
181  double ftof_expElec;
182  double ftof_expMuon;
183  double ftof_expPion;
184  double ftof_expKaon;
185  double ftof_expProt;
186  double ftof_expDeut;
187  int fcerenkovStatus[2];
188  double fcerenkovTime[2];
189  double fcerenkovPressure[2];
190  double fbeamtrackMomentum;
191  double fbeamtrackEnergy;
192  double fbeamtrackPos[4];
193  double fbeamtrackDir[3];
194  int fbeamtrackPdg;
195  int fbeamtrackID;
196  int fBIAndTimingMatched;
197 
198  // Check what the beam track will do in the detector. The beam track couldn't be the same with the primary beam PFParticle.
199  int fbeamtrack_truth_Origin;
200  double fbeamtrack_truth_EndPos[4];
201  double fbeamtrack_truth_Momentum[4];
202  double fbeamtrack_truth_KinEnergy;
203  double fbeamtrack_truth_Pt;
204  double fbeamtrack_truth_Mass;
205  double fbeamtrack_truth_Theta;
206  double fbeamtrack_truth_Phi;
207  double fbeamtrack_truth_TotalLength;
208  int fbeamtrack_truth_Process;
209  int fbeamtrack_truth_EndProcess;
210  int fbeamtrack_truth_Pdg;
211 
212  double fbeamtrack_truth_Pos_InTPCActive[4];
213  double fbeamtrack_truth_Momentum_InTPCActive[4];
214  double fbeamtrack_truth_P_InTPCActive;
215  double fbeamtrack_truth_Pt_InTPCActive;
216  double fbeamtrack_truth_Theta_InTPCActive;
217  double fbeamtrack_truth_Phi_InTPCActive;
218  double fbeamtrack_truth_TotalLength_InTPCActive;
219  double fbeamtrack_truth_KinEnergy_InTPCActive;
220 
221  int fbeamtrack_truth_NDAUGTHERS;
222  int fbeamtrack_truthdaughter_TrackId[NMAXTRUTHDAUGTHERS];
223  int fbeamtrack_truthdaughter_Pdg[NMAXTRUTHDAUGTHERS];
224  int fbeamtrack_truthdaughter_Mother[NMAXTRUTHDAUGTHERS];
225  double fbeamtrack_truthdaughter_StartPosition[NMAXTRUTHDAUGTHERS][4];
226  double fbeamtrack_truthdaughter_EndPosition[NMAXTRUTHDAUGTHERS][4];
227  double fbeamtrack_truthdaughter_P[NMAXTRUTHDAUGTHERS];
228  double fbeamtrack_truthdaughter_Momentum[NMAXTRUTHDAUGTHERS][4];
229  double fbeamtrack_truthdaughter_EndMomentum[NMAXTRUTHDAUGTHERS][4];
230  double fbeamtrack_truthdaughter_Pt[NMAXTRUTHDAUGTHERS];
231  double fbeamtrack_truthdaughter_Mass[NMAXTRUTHDAUGTHERS];
232  double fbeamtrack_truthdaughter_Theta[NMAXTRUTHDAUGTHERS];
233  double fbeamtrack_truthdaughter_Phi[NMAXTRUTHDAUGTHERS];
234  double fbeamtrack_truthdaughter_TotalLength[NMAXTRUTHDAUGTHERS];
235  int fbeamtrack_truthdaughter_Process[NMAXTRUTHDAUGTHERS];
236  int fbeamtrack_truthdaughter_EndProcess[NMAXTRUTHDAUGTHERS];
237 
238  int fbeamtrack_truth_NDECAYDAUGTHERS;
239  int fbeamtrack_truthdecaydaughter_TrackId[NMAXTRUTHDAUGTHERS];
240  int fbeamtrack_truthdecaydaughter_Pdg[NMAXTRUTHDAUGTHERS];
241  int fbeamtrack_truthdecaydaughter_Mother[NMAXTRUTHDAUGTHERS];
242  double fbeamtrack_truthdecaydaughter_StartPosition[NMAXTRUTHDAUGTHERS][4];
243  double fbeamtrack_truthdecaydaughter_EndPosition[NMAXTRUTHDAUGTHERS][4];
244  double fbeamtrack_truthdecaydaughter_P[NMAXTRUTHDAUGTHERS];
245  double fbeamtrack_truthdecaydaughter_Momentum[NMAXTRUTHDAUGTHERS][4];
246  double fbeamtrack_truthdecaydaughter_EndMomentum[NMAXTRUTHDAUGTHERS][4];
247  double fbeamtrack_truthdecaydaughter_Pt[NMAXTRUTHDAUGTHERS];
248  double fbeamtrack_truthdecaydaughter_Mass[NMAXTRUTHDAUGTHERS];
249  double fbeamtrack_truthdecaydaughter_Theta[NMAXTRUTHDAUGTHERS];
250  double fbeamtrack_truthdecaydaughter_Phi[NMAXTRUTHDAUGTHERS];
251  double fbeamtrack_truthdecaydaughter_TotalLength[NMAXTRUTHDAUGTHERS];
252  int fbeamtrack_truthdecaydaughter_Process[NMAXTRUTHDAUGTHERS];
253  int fbeamtrack_truthdecaydaughter_EndProcess[NMAXTRUTHDAUGTHERS];
254 
255  // Reconstructed tracks/showers
256  int fNpfParticles;
257  double fprimaryVertex[3];
258  int fprimaryIstrack;
259  int fprimaryIsshower;
260  double fprimaryBDTScore;
261  int fprimaryNHits;
262  double fprimaryTheta;
263  double fprimaryPhi;
264  double fprimaryLength;
265  double fprimaryMomentum;
266  double fprimaryEndMomentum;
267  double fprimaryEndPosition[3];
268  double fprimaryStartPosition[3];
269  double fprimaryEndDirection[3];
270  double fprimaryStartDirection[3];
271  double fprimaryOpeningAngle;
272  int fprimaryShowerBestPlane;
273  double fprimaryShowerEnergy[3];
274  double fprimaryShowerMIPEnergy[3];
275  double fprimaryShowerdEdx[3];
276  double fprimaryMomentumByRangeProton;
277  double fprimaryMomentumByRangeMuon;
278  double fprimaryKineticEnergy[3];
279  double fprimaryRange[3];
280  double fprimaryTrkPitchC[3];
281  int fprimaryID;
282  double fprimaryT0;
283  double fprimaryEarliestHitPeakTime;
284 
285  int fprimaryPID_Pdg[3];
286  int fprimaryPID_Ndf[3];
287  double fprimaryPID_MinChi2[3];
288  double fprimaryPID_DeltaChi2[3];
289  double fprimaryPID_Chi2Proton[3];
290  double fprimaryPID_Chi2Kaon[3];
291  double fprimaryPID_Chi2Pion[3];
292  double fprimaryPID_Chi2Muon[3];
293  double fprimaryPID_MissingE[3];
294  double fprimaryPID_MissingEavg[3];
295  double fprimaryPID_PIDA[3];
296 
297  int fprimary_truth_Origin;
298  int fprimary_truth_TrackId;
299  int fprimary_truth_Pdg;
300  int fprimary_truth_Mother;
301  double fprimary_truth_StartPosition[4];
302  double fprimary_truth_EndPosition[4];
303  double fprimary_truth_P;
304  double fprimary_truth_Momentum[4];
305  double fprimary_truth_EndMomentum[4];
306  double fprimary_truth_Pt;
307  double fprimary_truth_Mass;
308  double fprimary_truth_Theta;
309  double fprimary_truth_Phi;
310  double fprimary_truth_TotalLength;
311  int fprimary_truth_Process;
312  int fprimary_truth_EndProcess;
313  int fprimary_truth_Isbeammatched;
314  double fprimary_truth_EkinAtVertex;
315  double fprimary_truth_EkinAtVertex_notcorrected;
316  double fprimary_truth_Pos_InTPCActive[4];
317  double fprimary_truth_Momentum_InTPCActive[4];
318  double fprimary_truth_P_InTPCActive;
319  double fprimary_truth_Pt_InTPCActive;
320  double fprimary_truth_Theta_InTPCActive;
321  double fprimary_truth_Phi_InTPCActive;
322  double fprimary_truth_TotalLength_InTPCActive;
323  double fprimary_truth_KinEnergy_InTPCActive;
324 
325  int fprimary_truth_NDAUGTHERS;
326  int fprimary_truthdaughter_TrackId[NMAXTRUTHDAUGTHERS];
327  int fprimary_truthdaughter_Pdg[NMAXTRUTHDAUGTHERS];
328  int fprimary_truthdaughter_Mother[NMAXTRUTHDAUGTHERS];
329  double fprimary_truthdaughter_StartPosition[NMAXTRUTHDAUGTHERS][4];
330  double fprimary_truthdaughter_EndPosition[NMAXTRUTHDAUGTHERS][4];
331  double fprimary_truthdaughter_P[NMAXTRUTHDAUGTHERS];
332  double fprimary_truthdaughter_Momentum[NMAXTRUTHDAUGTHERS][4];
333  double fprimary_truthdaughter_EndMomentum[NMAXTRUTHDAUGTHERS][4];
334  double fprimary_truthdaughter_Pt[NMAXTRUTHDAUGTHERS];
335  double fprimary_truthdaughter_Mass[NMAXTRUTHDAUGTHERS];
336  double fprimary_truthdaughter_Theta[NMAXTRUTHDAUGTHERS];
337  double fprimary_truthdaughter_Phi[NMAXTRUTHDAUGTHERS];
338  double fprimary_truthdaughter_TotalLength[NMAXTRUTHDAUGTHERS];
339  int fprimary_truthdaughter_Process[NMAXTRUTHDAUGTHERS];
340  int fprimary_truthdaughter_EndProcess[NMAXTRUTHDAUGTHERS];
341 
342  int fprimary_truth_NDECAYDAUGTHERS;
343  int fprimary_truthdecaydaughter_TrackId[NMAXTRUTHDAUGTHERS];
344  int fprimary_truthdecaydaughter_Pdg[NMAXTRUTHDAUGTHERS];
345  int fprimary_truthdecaydaughter_Mother[NMAXTRUTHDAUGTHERS];
346  double fprimary_truthdecaydaughter_StartPosition[NMAXTRUTHDAUGTHERS][4];
347  double fprimary_truthdecaydaughter_EndPosition[NMAXTRUTHDAUGTHERS][4];
348  double fprimary_truthdecaydaughter_P[NMAXTRUTHDAUGTHERS];
349  double fprimary_truthdecaydaughter_Momentum[NMAXTRUTHDAUGTHERS][4];
350  double fprimary_truthdecaydaughter_EndMomentum[NMAXTRUTHDAUGTHERS][4];
351  double fprimary_truthdecaydaughter_Pt[NMAXTRUTHDAUGTHERS];
352  double fprimary_truthdecaydaughter_Mass[NMAXTRUTHDAUGTHERS];
353  double fprimary_truthdecaydaughter_Theta[NMAXTRUTHDAUGTHERS];
354  double fprimary_truthdecaydaughter_Phi[NMAXTRUTHDAUGTHERS];
355  double fprimary_truthdecaydaughter_TotalLength[NMAXTRUTHDAUGTHERS];
356  int fprimary_truthdecaydaughter_Process[NMAXTRUTHDAUGTHERS];
357  int fprimary_truthdecaydaughter_EndProcess[NMAXTRUTHDAUGTHERS];
358 
359  // Daughters from primary
360  int fNDAUGHTERS;
361  double fdaughterVertex[3];
362  int fdaughterIstrack[NMAXDAUGTHERS];
363  int fdaughterIsshower[NMAXDAUGTHERS];
364  int fdaughterNHits[NMAXDAUGTHERS];
365  double fdaughterTheta[NMAXDAUGTHERS];
366  double fdaughterPhi[NMAXDAUGTHERS];
367  double fdaughterLength[NMAXDAUGTHERS];
368  double fdaughterMomentum[NMAXDAUGTHERS];
369  double fdaughterEndMomentum[NMAXDAUGTHERS];
370  double fdaughterEndPosition[NMAXDAUGTHERS][3];
371  double fdaughterStartPosition[NMAXDAUGTHERS][3];
372  double fdaughterEndDirection[NMAXDAUGTHERS][3];
373  double fdaughterStartDirection[NMAXDAUGTHERS][3];
374  double fdaughterOpeningAngle[NMAXDAUGTHERS];
375  double fdaughterShowerEnergy[NMAXDAUGTHERS][3];
376  double fdaughterShowerMIPEnergy[NMAXDAUGTHERS][3];
377  double fdaughterShowerdEdx[NMAXDAUGTHERS][3];
378  int fdaughterShowerBestPlane[NMAXDAUGTHERS];
379  double fdaughterMomentumByRangeProton[NMAXDAUGTHERS];
380  double fdaughterMomentumByRangeMuon[NMAXDAUGTHERS];
381  double fdaughterKineticEnergy[NMAXDAUGTHERS][3];
382  double fdaughterRange[NMAXDAUGTHERS][3];
383  double fdaughterTrkPitchC[NMAXDAUGTHERS][3];
384  int fdaughterID[NMAXDAUGTHERS];
385  double fdaughterT0[NMAXDAUGTHERS];
386 
387  int fdaughterPID_Pdg[NMAXDAUGTHERS][3];
388  int fdaughterPID_Ndf[NMAXDAUGTHERS][3];
389  double fdaughterPID_MinChi2[NMAXDAUGTHERS][3];
390  double fdaughterPID_DeltaChi2[NMAXDAUGTHERS][3];
391  double fdaughterPID_Chi2Proton[NMAXDAUGTHERS][3];
392  double fdaughterPID_Chi2Kaon[NMAXDAUGTHERS][3];
393  double fdaughterPID_Chi2Pion[NMAXDAUGTHERS][3];
394  double fdaughterPID_Chi2Muon[NMAXDAUGTHERS][3];
395  double fdaughterPID_MissingE[NMAXDAUGTHERS][3];
396  double fdaughterPID_MissingEavg[NMAXDAUGTHERS][3];
397  double fdaughterPID_PIDA[NMAXDAUGTHERS][3];
398 
399  int fdaughter_truth_TrackId[NMAXDAUGTHERS];
400  int fdaughter_truth_Pdg[NMAXDAUGTHERS];
401  int fdaughter_truth_Mother[NMAXDAUGTHERS];
402  double fdaughter_truth_StartPosition[NMAXDAUGTHERS][4];
403  double fdaughter_truth_EndPosition[NMAXDAUGTHERS][4];
404  double fdaughter_truth_P[NMAXDAUGTHERS];
405  double fdaughter_truth_Momentum[NMAXDAUGTHERS][4];
406  double fdaughter_truth_EndMomentum[NMAXDAUGTHERS][4];
407  double fdaughter_truth_Pt[NMAXDAUGTHERS];
408  double fdaughter_truth_Mass[NMAXDAUGTHERS];
409  double fdaughter_truth_Theta[NMAXDAUGTHERS];
410  double fdaughter_truth_Phi[NMAXDAUGTHERS];
411  double fdaughter_truth_TotalLength[NMAXDAUGTHERS];
412  int fdaughter_truth_Process[NMAXDAUGTHERS];
413  int fdaughter_truth_EndProcess[NMAXDAUGTHERS];
414 
415  int fNGRANDDAUGHTERS;
416  double fgranddaughterVertex[NMAXDAUGTHERS][3];
417  int fgranddaughterIstrack[NMAXDAUGTHERS];
418  int fgranddaughterIsshower[NMAXDAUGTHERS];
419  int fgranddaughterNHits[NMAXDAUGTHERS];
420  double fgranddaughterTheta[NMAXDAUGTHERS];
421  double fgranddaughterPhi[NMAXDAUGTHERS];
422  double fgranddaughterLength[NMAXDAUGTHERS];
423  double fgranddaughterMomentum[NMAXDAUGTHERS];
424  double fgranddaughterEndMomentum[NMAXDAUGTHERS];
425  double fgranddaughterEndPosition[NMAXDAUGTHERS][3];
426  double fgranddaughterStartPosition[NMAXDAUGTHERS][3];
427  double fgranddaughterEndDirection[NMAXDAUGTHERS][3];
428  double fgranddaughterStartDirection[NMAXDAUGTHERS][3];
429  double fgranddaughterOpeningAngle[NMAXDAUGTHERS];
430  double fgranddaughterShowerEnergy[NMAXDAUGTHERS][3];
431  double fgranddaughterShowerMIPEnergy[NMAXDAUGTHERS][3];
432  double fgranddaughterShowerdEdx[NMAXDAUGTHERS][3];
433  int fgranddaughterShowerBestPlane[NMAXDAUGTHERS];
434  double fgranddaughterMomentumByRangeProton[NMAXDAUGTHERS];
435  double fgranddaughterMomentumByRangeMuon[NMAXDAUGTHERS];
436  double fgranddaughterKineticEnergy[NMAXDAUGTHERS][3];
437  double fgranddaughterRange[NMAXDAUGTHERS][3];
438  double fgranddaughterTrkPitchC[NMAXDAUGTHERS][3];
439  int fgranddaughterID[NMAXDAUGTHERS];
440  int fgranddaughterMotherID[NMAXDAUGTHERS];
441  double fgranddaughterT0[NMAXDAUGTHERS];
442 
443  int fgranddaughterPID_Pdg[NMAXDAUGTHERS][3];
444  int fgranddaughterPID_Ndf[NMAXDAUGTHERS][3];
445  double fgranddaughterPID_MinChi2[NMAXDAUGTHERS][3];
446  double fgranddaughterPID_DeltaChi2[NMAXDAUGTHERS][3];
447  double fgranddaughterPID_Chi2Proton[NMAXDAUGTHERS][3];
448  double fgranddaughterPID_Chi2Kaon[NMAXDAUGTHERS][3];
449  double fgranddaughterPID_Chi2Pion[NMAXDAUGTHERS][3];
450  double fgranddaughterPID_Chi2Muon[NMAXDAUGTHERS][3];
451  double fgranddaughterPID_MissingE[NMAXDAUGTHERS][3];
452  double fgranddaughterPID_MissingEavg[NMAXDAUGTHERS][3];
453  double fgranddaughterPID_PIDA[NMAXDAUGTHERS][3];
454 
455  int fgranddaughter_truth_TrackId[NMAXDAUGTHERS];
456  int fgranddaughter_truth_Pdg[NMAXDAUGTHERS];
457  int fgranddaughter_truth_Mother[NMAXDAUGTHERS];
458  double fgranddaughter_truth_StartPosition[NMAXDAUGTHERS][4];
459  double fgranddaughter_truth_EndPosition[NMAXDAUGTHERS][4];
460  double fgranddaughter_truth_P[NMAXDAUGTHERS];
461  double fgranddaughter_truth_Momentum[NMAXDAUGTHERS][4];
462  double fgranddaughter_truth_EndMomentum[NMAXDAUGTHERS][4];
463  double fgranddaughter_truth_Pt[NMAXDAUGTHERS];
464  double fgranddaughter_truth_Mass[NMAXDAUGTHERS];
465  double fgranddaughter_truth_Theta[NMAXDAUGTHERS];
466  double fgranddaughter_truth_Phi[NMAXDAUGTHERS];
467  double fgranddaughter_truth_TotalLength[NMAXDAUGTHERS];
468  int fgranddaughter_truth_Process[NMAXDAUGTHERS];
469  int fgranddaughter_truth_EndProcess[NMAXDAUGTHERS];
470 
471 };
472 
473 // -----------------------------------------------------------------------------
474 protoana::ProtoDUNEAnalysisTree::ProtoDUNEAnalysisTree(fhicl::ParameterSet const & p)
475  :
476  EDAnalyzer(p),
477  dataUtil(p.get<fhicl::ParameterSet>("DataUtils")),
478  beamlineUtil(p.get<fhicl::ParameterSet>("BeamLineUtils")),
479  fBeamModuleLabel(p.get< art::InputTag >("BeamModuleLabel")),
480  fCalorimetryTag(p.get<std::string>("CalorimetryTag")),
481  fParticleIDTag(p.get<std::string>("ParticleIDTag")),
482  fTrackerTag(p.get<std::string>("TrackerTag")),
483  fShowerTag(p.get<std::string>("ShowerTag")),
484  fPFParticleTag(p.get<std::string>("PFParticleTag")),
485  fGeneratorTag(p.get<std::string>("GeneratorTag")),
486  fSimulationTag(p.get<std::string>("SimulationTag")),
487  fVerbose(p.get<int>("Verbose"))
488 {
489 
490 }
491 
492 // -----------------------------------------------------------------------------
493 void protoana::ProtoDUNEAnalysisTree::beginRun(const art::Run&){
494 
495 }
496 
497 // -----------------------------------------------------------------------------
498 void protoana::ProtoDUNEAnalysisTree::analyze(art::Event const & evt){
499 
500  // Initialise tree parameters
501  Initialise();
502 
503  fRun = evt.run();
504  fSubRun = evt.subRun();
505  fevent = evt.id().event();
506  art::Timestamp ts = evt.time();
507  if (ts.timeHigh() == 0){
508  TTimeStamp ts2(ts.timeLow());
509  fTimeStamp = ts2.AsDouble();
510  }
511  else{
512  TTimeStamp ts2(ts.timeHigh(), ts.timeLow());
513  fTimeStamp = ts2.AsDouble();
514  }
515 
516  // Get number of active fembs
517  if(!evt.isRealData()){
518  for(int k=0; k < 6; k++)
519  fNactivefembs[k] = 20;
520  }
521  else{
522  for(int k=0; k < 6; k++)
523  fNactivefembs[k] = dataUtil.GetNActiveFembsForAPA(evt, k);
524  }
525 
526  // Fill beam track info
527  bool isbeamTriggerEvent = FillPrimaryBeamParticle(evt);
528  // If this is not a beam trigger then do not save
529  if(!isbeamTriggerEvent) return;
530 
531  /*
532  // Now we want to access the output from Pandora. This comes in the form of particle flow objects (recob::PFParticle).
533  // The primary PFParticles are those we want to consider and these PFParticles then have a hierarchy of daughters that
534  // describe the whole interaction of a given primary particle
535  //
536  // / daughter track / grand-daughter track
537  // / /
538  // primary track / /
539  // ---------------- ---- daughter track ----- grand-daughter track
540  // \
541  // /\-
542  // /\\-- daughter shower
543  //
544  // The above primary PFParticle will have links to three daughter particles, two track-like and one shower-like
545  */
546 
547  //trmom.SetMinLength(100);
548 
549  // Get all of the PFParticles, by default from the "pandora" product
550  auto recoParticles = evt.getValidHandle<std::vector<recob::PFParticle>>(fPFParticleTag);
551  //std::cout << "All primary pfParticles = " << pfpUtil.GetNumberPrimaryPFParticle(evt,fPFParticleTag) << std::endl;
552 
553  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService>()->DataFor(evt);
554  auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService>()->DataFor(evt, clockData);
555 
556  // We'd like to find the beam particle. Pandora tries to do this for us, so let's use the PFParticle utility
557  // to look for it. Pandora reconstructs slices containing one (or sometimes more) primary PFParticles. These
558  // are tagged as either beam or cosmic for ProtoDUNE. This function automatically considers only those
559  // PFParticles considered as primary
560  std::vector<const recob::PFParticle*> pfParticles = pfpUtil.GetPFParticlesFromBeamSlice(evt,fPFParticleTag);
561  fNpfParticles = pfParticles.size();
562 
563  // We can now look at these particles
564  for(const recob::PFParticle* particle : pfParticles){
565 
566  FillPrimaryPFParticle(evt, clockData, detProp, particle);
567 
568  // Find the particle vertex. We need the tracker tag here because we need to do a bit of
569  // additional work if the PFParticle is track-like to find the vertex.
570  const TVector3 vtx = pfpUtil.GetPFParticleVertex(*particle,evt,fPFParticleTag,fTrackerTag);
571  fprimaryVertex[0] = vtx.X(); fprimaryVertex[1] = vtx.Y(); fprimaryVertex[2] = vtx.Z();
572 
573  // Now we can look for the interaction point of the particle if one exists, i.e where the particle
574  // scatters off an argon nucleus. Shower-like objects won't have an interaction point, so we can
575  // check this by making sure we get a sensible position
576  const TVector3 interactionVtx = pfpUtil.GetPFParticleSecondaryVertex(*particle,evt,fPFParticleTag,fTrackerTag);
577  fdaughterVertex[0] = interactionVtx.X(); fdaughterVertex[1] = interactionVtx.Y(); fdaughterVertex[2] = interactionVtx.Z();
578 
579  // Maximum number of daugthers to be processed
580  if(particle->NumDaughters() > NMAXDAUGTHERS)
581  std::cout << "INFO::Number of daughters is " << particle->NumDaughters() << ". Only the first NMAXDAUGTHERS are processed." << std::endl;
582 
583  // Let's get the daughter PFParticles... we can do this simply without the utility
584  for(const int daughterID : particle->Daughters()){
585  // Daughter ID is the element of the original recoParticle vector
586  const recob::PFParticle *daughterParticle = &(recoParticles->at(daughterID));
587  //std::cout << "Daughter " << daughterID << " has " << daughterParticle->NumDaughters() << " daughters" << std::endl;
588 
589  // Fill tree with daughter info
590  FillPrimaryDaughterPFParticle(evt, clockData, detProp, daughterParticle, daughterID);
591 
592  // Get the secondary vertex from daughter interactions
593  const TVector3 secinteractionVtx = pfpUtil.GetPFParticleSecondaryVertex(*daughterParticle,evt,fPFParticleTag,fTrackerTag);
594  fgranddaughterVertex[fNDAUGHTERS][0] = secinteractionVtx.X(); fgranddaughterVertex[fNDAUGHTERS][1] = secinteractionVtx.Y(); fgranddaughterVertex[fNDAUGHTERS][2] = secinteractionVtx.Z();
595 
596  for(const int gdaughterID : daughterParticle->Daughters()){
597  FillPrimaryGrandDaughterPFParticle(evt, clockData, detProp, daughterParticle, gdaughterID, daughterID);
598 
599  // Only process NMAXDAUGTHERS
600  if(fNGRANDDAUGHTERS > NMAXDAUGTHERS) break;
601  }
602 
603  // Only process NMAXDAUGTHERS
604  if(fNDAUGHTERS > NMAXDAUGTHERS) break;
605 
606  }
607 
608  // For now only consider the first primary track. Need a proper treatment if more than one primary particles are found.
609  break;
610  }
611 
612  // Fill trees
613  fPandoraBeam->Fill();
614  fBeamTrack->Fill();
615  fDaughterTree->Fill();
616  fGrandDaughterTree->Fill();
617 
618 }
619 
620 // -----------------------------------------------------------------------------
621 void protoana::ProtoDUNEAnalysisTree::endJob(){
622 
623 }
624 
625 // -----------------------------------------------------------------------------
626 bool protoana::ProtoDUNEAnalysisTree::FillPrimaryBeamParticle(art::Event const & evt){
627 
628  bool beamTriggerEvent = false;
629  // If this event is MC then we can check what the true beam particle is
630  if(!evt.isRealData()){
631  // Firstly we need to get the list of MCTruth objects from the generator. The standard protoDUNE
632  // simulation has fGeneratorTag = "generator"
633  auto mcTruths = evt.getValidHandle<std::vector<simb::MCTruth>>(fGeneratorTag);
634  fbeamtrack_truth_Origin = (*mcTruths)[0].Origin();
635  fbeamNTracks = 1;
636 
637  // mcTruths is basically a pointer to an std::vector of simb::MCTruth objects. There should only be one
638  // of these, so we pass the first element into the function to get the good particle
639  const simb::MCParticle* geantGoodParticle = truthUtil.GetGeantGoodParticle((*mcTruths)[0],evt);
640 
641  if(geantGoodParticle != 0x0){
642 
643  beamTriggerEvent = true;
644  fbeamtrigger = 12;
645  fbeamtrackPos[0] = geantGoodParticle->Vx();
646  fbeamtrackPos[1] = geantGoodParticle->Vy();
647  fbeamtrackPos[2] = geantGoodParticle->Vz();
648  fbeamtrackPos[3] = geantGoodParticle->T();
649  fbeamtrackMomentum = geantGoodParticle->P();
650  fbeamtrackEnergy = geantGoodParticle->E();
651  fbeamtrackPdg = geantGoodParticle->PdgCode();
652  fbeamtrackID = geantGoodParticle->TrackId();
653 
654  fbeamtrack_truth_EndPos[0] = geantGoodParticle->EndX();
655  fbeamtrack_truth_EndPos[1] = geantGoodParticle->EndY();
656  fbeamtrack_truth_EndPos[2] = geantGoodParticle->EndZ();
657  fbeamtrack_truth_EndPos[3] = geantGoodParticle->EndT();
658  fbeamtrack_truth_Momentum[0] = geantGoodParticle->Px();
659  fbeamtrack_truth_Momentum[1] = geantGoodParticle->Py();
660  fbeamtrack_truth_Momentum[2] = geantGoodParticle->Pz();
661  fbeamtrack_truth_Momentum[3] = geantGoodParticle->E();
662  fbeamtrack_truth_Pt = geantGoodParticle->Pt();
663  fbeamtrack_truth_Mass = geantGoodParticle->Mass();
664  fbeamtrack_truth_Theta = geantGoodParticle->Momentum().Theta();
665  fbeamtrack_truth_Phi = geantGoodParticle->Momentum().Phi();
666  fbeamtrack_truth_TotalLength = geantGoodParticle->Trajectory().TotalLength();
667  fbeamtrack_truth_KinEnergy = std::sqrt(geantGoodParticle->P()*geantGoodParticle->P() + geantGoodParticle->Mass()*geantGoodParticle->Mass()) - geantGoodParticle->Mass();
668  //fbeamtrack_truth_Process = int(geantGoodParticle->Trajectory().ProcessToKey(geantGoodParticle->Process()));
669  //fbeamtrack_truth_EndProcess = int(geantGoodParticle->Trajectory().ProcessToKey(geantGoodParticle->EndProcess()));
670  fbeamtrack_truth_Process = truthUtil.GetProcessKey(geantGoodParticle->Process());
671  fbeamtrack_truth_EndProcess = truthUtil.GetProcessKey(geantGoodParticle->EndProcess());
672  // Truth Pdg may not be the same as the beam PID.
673  fbeamtrack_truth_Pdg = geantGoodParticle->PdgCode();
674  // Get the position and kinematics in the TPC active volume
675  int trajintpcactivevol = truthUtil.GetFirstTrajectoryPointInTPCActiveVolume(*geantGoodParticle, fActiveTPCBoundsX[0], fActiveTPCBoundsX[1], fActiveTPCBoundsY[0], fActiveTPCBoundsY[1], fActiveTPCBoundsZ[0], fActiveTPCBoundsZ[1]);
676  if(trajintpcactivevol > 0){
677  fbeamtrack_truth_Pos_InTPCActive[0] = geantGoodParticle->Vx(trajintpcactivevol);
678  fbeamtrack_truth_Pos_InTPCActive[1] = geantGoodParticle->Vy(trajintpcactivevol);
679  fbeamtrack_truth_Pos_InTPCActive[2] = geantGoodParticle->Vz(trajintpcactivevol);
680  fbeamtrack_truth_Pos_InTPCActive[3] = geantGoodParticle->T(trajintpcactivevol);
681  fbeamtrack_truth_Momentum_InTPCActive[0] = geantGoodParticle->Px(trajintpcactivevol);
682  fbeamtrack_truth_Momentum_InTPCActive[1] = geantGoodParticle->Py(trajintpcactivevol);
683  fbeamtrack_truth_Momentum_InTPCActive[2] = geantGoodParticle->Pz(trajintpcactivevol);
684  fbeamtrack_truth_Momentum_InTPCActive[3] = geantGoodParticle->E(trajintpcactivevol);
685  fbeamtrack_truth_P_InTPCActive = geantGoodParticle->P(trajintpcactivevol);
686  fbeamtrack_truth_Pt_InTPCActive = geantGoodParticle->Pt(trajintpcactivevol);
687  fbeamtrack_truth_Theta_InTPCActive = geantGoodParticle->Momentum(trajintpcactivevol).Theta();
688  fbeamtrack_truth_Phi_InTPCActive = geantGoodParticle->Momentum(trajintpcactivevol).Phi();
689  fbeamtrack_truth_TotalLength_InTPCActive = truthUtil.GetMCParticleLengthInTPCActiveVolume(*geantGoodParticle, fActiveTPCBoundsX[0], fActiveTPCBoundsX[1], fActiveTPCBoundsY[0], fActiveTPCBoundsY[1], fActiveTPCBoundsZ[0], fActiveTPCBoundsZ[1]);
690  fbeamtrack_truth_KinEnergy_InTPCActive = std::sqrt(geantGoodParticle->P(trajintpcactivevol)*geantGoodParticle->P(trajintpcactivevol) + geantGoodParticle->Mass()*geantGoodParticle->Mass()) - geantGoodParticle->Mass();
691  //unsigned int ntrajpoints = geantGoodParticle->NumberTrajectoryPoints();
692  }
693 
694  // If this beam particle has interacted in the detector, find its daugthers.
695  std::vector<art::Ptr<simb::MCParticle> > mcbeampartilelist;
696  auto mcbeamHandle = evt.getHandle<std::vector<simb::MCParticle> >(fSimulationTag);
697  if (mcbeamHandle) art::fill_ptr_vector(mcbeampartilelist, mcbeamHandle);
698 
699  for(size_t p=0; p < mcbeampartilelist.size(); p++){
700  auto & beamprim_mcparticle = mcbeampartilelist[p];
701  if(beamprim_mcparticle->Mother() == geantGoodParticle->TrackId()){
702  fbeamtrack_truthdaughter_TrackId[fbeamtrack_truth_NDAUGTHERS] = beamprim_mcparticle->TrackId();
703  fbeamtrack_truthdaughter_Pdg[fbeamtrack_truth_NDAUGTHERS] = beamprim_mcparticle->PdgCode();
704  fbeamtrack_truthdaughter_Mother[fbeamtrack_truth_NDAUGTHERS] = beamprim_mcparticle->Mother();
705  fbeamtrack_truthdaughter_StartPosition[fbeamtrack_truth_NDAUGTHERS][0] = beamprim_mcparticle->Vx();
706  fbeamtrack_truthdaughter_StartPosition[fbeamtrack_truth_NDAUGTHERS][1] = beamprim_mcparticle->Vy();
707  fbeamtrack_truthdaughter_StartPosition[fbeamtrack_truth_NDAUGTHERS][2] = beamprim_mcparticle->Vz();
708  fbeamtrack_truthdaughter_StartPosition[fbeamtrack_truth_NDAUGTHERS][3] = beamprim_mcparticle->T();
709  fbeamtrack_truthdaughter_EndPosition[fbeamtrack_truth_NDAUGTHERS][0] = beamprim_mcparticle->EndX();
710  fbeamtrack_truthdaughter_EndPosition[fbeamtrack_truth_NDAUGTHERS][1] = beamprim_mcparticle->EndY();
711  fbeamtrack_truthdaughter_EndPosition[fbeamtrack_truth_NDAUGTHERS][2] = beamprim_mcparticle->EndZ();
712  fbeamtrack_truthdaughter_EndPosition[fbeamtrack_truth_NDAUGTHERS][3] = beamprim_mcparticle->EndT();
713  fbeamtrack_truthdaughter_P[fbeamtrack_truth_NDAUGTHERS] = beamprim_mcparticle->P();
714  fbeamtrack_truthdaughter_Momentum[fbeamtrack_truth_NDAUGTHERS][0] = beamprim_mcparticle->Px();
715  fbeamtrack_truthdaughter_Momentum[fbeamtrack_truth_NDAUGTHERS][1] = beamprim_mcparticle->Py();
716  fbeamtrack_truthdaughter_Momentum[fbeamtrack_truth_NDAUGTHERS][2] = beamprim_mcparticle->Pz();
717  fbeamtrack_truthdaughter_Momentum[fbeamtrack_truth_NDAUGTHERS][3] = beamprim_mcparticle->E();
718  fbeamtrack_truthdaughter_Pt[fbeamtrack_truth_NDAUGTHERS] = beamprim_mcparticle->Pt();
719  fbeamtrack_truthdaughter_Mass[fbeamtrack_truth_NDAUGTHERS] = beamprim_mcparticle->Mass();
720  fbeamtrack_truthdaughter_EndMomentum[fbeamtrack_truth_NDAUGTHERS][0] = beamprim_mcparticle->EndPx();
721  fbeamtrack_truthdaughter_EndMomentum[fbeamtrack_truth_NDAUGTHERS][1] = beamprim_mcparticle->EndPy();
722  fbeamtrack_truthdaughter_EndMomentum[fbeamtrack_truth_NDAUGTHERS][2] = beamprim_mcparticle->EndPz();
723  fbeamtrack_truthdaughter_EndMomentum[fbeamtrack_truth_NDAUGTHERS][3] = beamprim_mcparticle->EndE();
724  fbeamtrack_truthdaughter_Theta[fbeamtrack_truth_NDAUGTHERS] = beamprim_mcparticle->Momentum().Theta();
725  fbeamtrack_truthdaughter_Phi[fbeamtrack_truth_NDAUGTHERS] = beamprim_mcparticle->Momentum().Phi();
726  fbeamtrack_truthdaughter_TotalLength[fbeamtrack_truth_NDAUGTHERS] = beamprim_mcparticle->Trajectory().TotalLength();
727  fbeamtrack_truthdaughter_Process[fbeamtrack_truth_NDAUGTHERS] = truthUtil.GetProcessKey(beamprim_mcparticle->Process());
728  fbeamtrack_truthdaughter_EndProcess[fbeamtrack_truth_NDAUGTHERS] = truthUtil.GetProcessKey(beamprim_mcparticle->EndProcess());
729 
730  fbeamtrack_truth_NDAUGTHERS++;
731  if(fbeamtrack_truth_NDAUGTHERS > NMAXTRUTHDAUGTHERS) break;
732 
733  if(truthUtil.GetProcessKey(beamprim_mcparticle->EndProcess()) == 30){ // Decay particle
734  for(size_t pp=0; pp < mcbeampartilelist.size(); pp++){
735  auto & decay_beamprim_mcparticle = mcbeampartilelist[pp];
736  if(decay_beamprim_mcparticle->Mother() == beamprim_mcparticle->TrackId()){
737  fbeamtrack_truthdecaydaughter_TrackId[fbeamtrack_truth_NDECAYDAUGTHERS] = decay_beamprim_mcparticle->TrackId();
738  fbeamtrack_truthdecaydaughter_Pdg[fbeamtrack_truth_NDECAYDAUGTHERS] = decay_beamprim_mcparticle->PdgCode();
739  fbeamtrack_truthdecaydaughter_Mother[fbeamtrack_truth_NDECAYDAUGTHERS] = decay_beamprim_mcparticle->Mother();
740  fbeamtrack_truthdecaydaughter_StartPosition[fbeamtrack_truth_NDECAYDAUGTHERS][0] = decay_beamprim_mcparticle->Vx();
741  fbeamtrack_truthdecaydaughter_StartPosition[fbeamtrack_truth_NDECAYDAUGTHERS][1] = decay_beamprim_mcparticle->Vy();
742  fbeamtrack_truthdecaydaughter_StartPosition[fbeamtrack_truth_NDECAYDAUGTHERS][2] = decay_beamprim_mcparticle->Vz();
743  fbeamtrack_truthdecaydaughter_StartPosition[fbeamtrack_truth_NDECAYDAUGTHERS][3] = decay_beamprim_mcparticle->T();
744  fbeamtrack_truthdecaydaughter_EndPosition[fbeamtrack_truth_NDECAYDAUGTHERS][0] = decay_beamprim_mcparticle->EndX();
745  fbeamtrack_truthdecaydaughter_EndPosition[fbeamtrack_truth_NDECAYDAUGTHERS][1] = decay_beamprim_mcparticle->EndY();
746  fbeamtrack_truthdecaydaughter_EndPosition[fbeamtrack_truth_NDECAYDAUGTHERS][2] = decay_beamprim_mcparticle->EndZ();
747  fbeamtrack_truthdecaydaughter_EndPosition[fbeamtrack_truth_NDECAYDAUGTHERS][3] = decay_beamprim_mcparticle->EndT();
748  fbeamtrack_truthdecaydaughter_P[fbeamtrack_truth_NDECAYDAUGTHERS] = decay_beamprim_mcparticle->P();
749  fbeamtrack_truthdecaydaughter_Momentum[fbeamtrack_truth_NDECAYDAUGTHERS][0] = decay_beamprim_mcparticle->Px();
750  fbeamtrack_truthdecaydaughter_Momentum[fbeamtrack_truth_NDECAYDAUGTHERS][1] = decay_beamprim_mcparticle->Py();
751  fbeamtrack_truthdecaydaughter_Momentum[fbeamtrack_truth_NDECAYDAUGTHERS][2] = decay_beamprim_mcparticle->Pz();
752  fbeamtrack_truthdecaydaughter_Momentum[fbeamtrack_truth_NDECAYDAUGTHERS][3] = decay_beamprim_mcparticle->E();
753  fbeamtrack_truthdecaydaughter_Pt[fbeamtrack_truth_NDECAYDAUGTHERS] = decay_beamprim_mcparticle->Pt();
754  fbeamtrack_truthdecaydaughter_Mass[fbeamtrack_truth_NDECAYDAUGTHERS] = decay_beamprim_mcparticle->Mass();
755  fbeamtrack_truthdecaydaughter_EndMomentum[fbeamtrack_truth_NDECAYDAUGTHERS][0] = decay_beamprim_mcparticle->EndPx();
756  fbeamtrack_truthdecaydaughter_EndMomentum[fbeamtrack_truth_NDECAYDAUGTHERS][1] = decay_beamprim_mcparticle->EndPy();
757  fbeamtrack_truthdecaydaughter_EndMomentum[fbeamtrack_truth_NDECAYDAUGTHERS][2] = decay_beamprim_mcparticle->EndPz();
758  fbeamtrack_truthdecaydaughter_EndMomentum[fbeamtrack_truth_NDECAYDAUGTHERS][3] = decay_beamprim_mcparticle->EndE();
759  fbeamtrack_truthdecaydaughter_Theta[fbeamtrack_truth_NDECAYDAUGTHERS] = decay_beamprim_mcparticle->Momentum().Theta();
760  fbeamtrack_truthdecaydaughter_Phi[fbeamtrack_truth_NDECAYDAUGTHERS] = decay_beamprim_mcparticle->Momentum().Phi();
761  fbeamtrack_truthdecaydaughter_TotalLength[fbeamtrack_truth_NDECAYDAUGTHERS] = decay_beamprim_mcparticle->Trajectory().TotalLength();
762  fbeamtrack_truthdecaydaughter_Process[fbeamtrack_truth_NDECAYDAUGTHERS] = truthUtil.GetProcessKey(decay_beamprim_mcparticle->Process());
763  fbeamtrack_truthdecaydaughter_EndProcess[fbeamtrack_truth_NDECAYDAUGTHERS] = truthUtil.GetProcessKey(decay_beamprim_mcparticle->EndProcess());
764 
765  fbeamtrack_truth_NDECAYDAUGTHERS++;
766  if(fbeamtrack_truth_NDECAYDAUGTHERS > NMAXTRUTHDAUGTHERS) break;
767  }
768  }
769  } // decay particle
770 
771  }
772 
773  } // mc particle vector
774 
775  } // beam geant particle found
776  } // is mc
777  else{ // data
778  // For data we can see if this event comes from a beam trigger
779  beamTriggerEvent = dataUtil.IsBeamTrigger(evt);
780 
781  std::vector< art::Ptr<beam::ProtoDUNEBeamEvent> > beaminfo;
782  auto pdbeamHandle = evt.getHandle< std::vector<beam::ProtoDUNEBeamEvent> >(fBeamModuleLabel);
783  if (pdbeamHandle) art::fill_ptr_vector(beaminfo, pdbeamHandle);
784 
785  fbeamNTracks = (int)beaminfo.size();
786  // For now only consider events with one beam track only
787  if(fbeamNTracks != 1) return false;
788 
789  for(unsigned int i = 0; i < beaminfo.size(); ++i){
790  fBIAndTimingMatched = beaminfo[i]->CheckIsMatched();
791  fbeamtrigger = beaminfo[i]->GetTimingTrigger();
792  fbeamtrackPos[3] = (double)beaminfo[i]->GetRDTimestamp();
793 
794  // If ToF is 0-3 there was a good match corresponding to the different pair-wise combinations of the upstream and downstream channels
795  if(beaminfo[i]->GetTOFChan() >= 0)
796  ftof = beaminfo[i]->GetTOF();
797 
798  // Get Cerenkov
799  if(beaminfo[i]->GetBITrigger() == 1){
800  fcerenkovStatus[0] = beaminfo[i]->GetCKov0Status();
801  fcerenkovStatus[1] = beaminfo[i]->GetCKov1Status();
802  fcerenkovTime[0] = beaminfo[i]->GetCKov0Time();
803  fcerenkovTime[1] = beaminfo[i]->GetCKov1Time();
804  fcerenkovPressure[0] = beaminfo[i]->GetCKov0Pressure();
805  fcerenkovPressure[1] = beaminfo[i]->GetCKov1Pressure();
806  }
807 
808  // Beam particle could have more than one tracks - for now take the first one, need to do this properly
809  auto & tracks = beaminfo[i]->GetBeamTracks();
810  if(!tracks.empty()){
811  fbeamtrackPos[0] = tracks[0].End().X();
812  fbeamtrackPos[1] = tracks[0].End().Y();
813  fbeamtrackPos[2] = tracks[0].End().Z();
814  fbeamtrackDir[0] = tracks[0].StartDirection().X();
815  fbeamtrackDir[1] = tracks[0].StartDirection().Y();
816  fbeamtrackDir[2] = tracks[0].StartDirection().Z();
817  }
818 
819  // Beam momentum
820  auto & beammom = beaminfo[i]->GetRecoBeamMomenta();
821  if(!beammom.empty()){
822  fbeamtrackMomentum = beammom[0];
823  ftof_expElec = beamlineUtil.ComputeTOF(11, beammom[0]);
824  ftof_expMuon = beamlineUtil.ComputeTOF(13, beammom[0]);
825  ftof_expPion = beamlineUtil.ComputeTOF(211, beammom[0]);
826  ftof_expKaon = beamlineUtil.ComputeTOF(321, beammom[0]);
827  ftof_expProt = beamlineUtil.ComputeTOF(2212, beammom[0]);
828  ftof_expDeut = beamlineUtil.ComputeTOF(1000010020, beammom[0]);
829  }
830 
831  // For now only take the first beam particle - need to add some criteria if more than one are found
832  break;
833 
834  }
835  }
836 
837  return beamTriggerEvent;
838 
839 }
840 
841 // -----------------------------------------------------------------------------
842 void protoana::ProtoDUNEAnalysisTree::FillPrimaryPFParticle(art::Event const & evt,
843  detinfo::DetectorClocksData const& clockData,
844  detinfo::DetectorPropertiesData const& detProp,
845  const recob::PFParticle* particle){
846 
847  // Pandora's BDT beam-cosmic score
848  fprimaryBDTScore = (double)pfpUtil.GetBeamCosmicScore(*particle,evt,fPFParticleTag);
849 
850  // NHits associated with this pfParticle
851  fprimaryNHits = (pfpUtil.GetPFParticleHits(*particle,evt,fPFParticleTag)).size();
852 
853  // Earliest hit peak time
854  fprimaryEarliestHitPeakTime = pfpUtil.GetPFParticleEarliestHitPeakTime(*particle,evt,fPFParticleTag);
855 
856  // Get the T0 for this pfParticle
857  std::vector<anab::T0> pfT0vec = pfpUtil.GetPFParticleT0(*particle,evt,fPFParticleTag);
858  if(!pfT0vec.empty())
859  fprimaryT0 = pfT0vec[0].Time();
860 
861  //std::vector<art::Ptr<simb::MCParticle> > mcpartilelist;
862  //mf::LogWarning("ProtoDUNEAnalysisTree") << " Getting MC Particle failed. Ignore if you run on data. " << std::endl;
863  //auto mcHandle = evt.getHandle<std::vector<simb::MCParticle> >(fSimulationTag);
864  //if (mcHandle)
865  //art::fill_ptr_vector(mcpartilelist, mcHandle);
866 
867  // "particle" is the pointer to our beam particle. The recob::Track or recob::Shower object
868  // of this particle might be more helpful. These return null pointers if not track-like / shower-like
869  const recob::Track* thisTrack = pfpUtil.GetPFParticleTrack(*particle, evt,fPFParticleTag,fTrackerTag);
870  const recob::Shower* thisShower = pfpUtil.GetPFParticleShower(*particle,evt,fPFParticleTag,fShowerTag);
871 
872  const simb::MCParticle* mcparticle = NULL;
873 
874  if(thisTrack != 0x0){
875  fprimaryIstrack = 1;
876  fprimaryIsshower = 0;
877  fprimaryID = thisTrack->ParticleId();
878  fprimaryTheta = thisTrack->Theta();
879  fprimaryPhi = thisTrack->Phi();
880  fprimaryLength = thisTrack->Length();
881  fprimaryMomentum = thisTrack->StartMomentum();
882  fprimaryEndMomentum = thisTrack->EndMomentum();
883  fprimaryEndPosition[0] = thisTrack->Trajectory().End().X();
884  fprimaryEndPosition[1] = thisTrack->Trajectory().End().Y();
885  fprimaryEndPosition[2] = thisTrack->Trajectory().End().Z();
886  fprimaryStartPosition[0] = thisTrack->Trajectory().Start().X();
887  fprimaryStartPosition[1] = thisTrack->Trajectory().Start().Y();
888  fprimaryStartPosition[2] = thisTrack->Trajectory().Start().Z();
889  fprimaryEndDirection[0] = thisTrack->Trajectory().EndDirection().X();
890  fprimaryEndDirection[1] = thisTrack->Trajectory().EndDirection().Y();
891  fprimaryEndDirection[2] = thisTrack->Trajectory().EndDirection().Z();
892  fprimaryStartDirection[0] = thisTrack->Trajectory().StartDirection().X();
893  fprimaryStartDirection[1] = thisTrack->Trajectory().StartDirection().Y();
894  fprimaryStartDirection[2] = thisTrack->Trajectory().StartDirection().Z();
895 
896  fprimaryMomentumByRangeMuon = trmom.GetTrackMomentum(thisTrack->Length(),13);
897  fprimaryMomentumByRangeProton = trmom.GetTrackMomentum(thisTrack->Length(),2212);
898 
899  // Calorimetry
900  std::vector<anab::Calorimetry> calovector = trackUtil.GetRecoTrackCalorimetry(*thisTrack, evt, fTrackerTag, fCalorimetryTag);
901  if(calovector.size() != 3 && fVerbose > 0)
902  std::cerr << "WARNING::Calorimetry vector size for primary is = " << calovector.size() << std::endl;
903 
904  for(size_t k = 0; k < calovector.size() && k<3; k++){
905  int plane = calovector[k].PlaneID().Plane;
906  if(plane < 0) continue;
907  if(plane > 2) continue;
908  fprimaryKineticEnergy[plane] = calovector[k].KineticEnergy();
909  fprimaryRange[plane] = calovector[k].Range();
910  fprimaryTrkPitchC[plane] = calovector[k].TrkPitchC();
911  }
912 
913  // PID
914  std::vector<anab::ParticleID> pids = trackUtil.GetRecoTrackPID(*thisTrack, evt, fTrackerTag, fParticleIDTag);
915  if(pids.size() != 3 && fVerbose > 0)
916  std::cerr << "WARNING::PID vector size for primary is = " << pids.size() << std::endl;
917 
918  for(size_t k = 0; k < pids.size() && k<3; k++){
919  int plane = pids[k].PlaneID().Plane;
920  if(plane < 0) continue;
921  if(plane > 2) continue;
922  fprimaryPID_Pdg[plane] = pids[plane].Pdg();
923  fprimaryPID_Ndf[plane] = pids[plane].Ndf();
924  fprimaryPID_MinChi2[plane] = pids[plane].MinChi2();
925  fprimaryPID_DeltaChi2[plane] = pids[plane].DeltaChi2();
926  fprimaryPID_Chi2Proton[plane] = pids[plane].Chi2Proton();
927  fprimaryPID_Chi2Kaon[plane] = pids[plane].Chi2Kaon();
928  fprimaryPID_Chi2Pion[plane] = pids[plane].Chi2Pion();
929  fprimaryPID_Chi2Muon[plane] = pids[plane].Chi2Muon();
930  fprimaryPID_MissingE[plane] = pids[plane].MissingE();
931  fprimaryPID_MissingEavg[plane] = pids[plane].MissingEavg();
932  fprimaryPID_PIDA[plane] = pids[plane].PIDA();
933  }
934 
935  // Get the true mc particle
936  //const simb::MCParticle* mcparticle = truthUtil.GetMCParticleFromRecoTrack(*thisTrack, evt, fTrackerTag);
937  mcparticle = truthUtil.GetMCParticleFromRecoTrack(clockData, *thisTrack, evt, fTrackerTag);
938  } // end is track
939  else if(thisShower != 0x0){
940  fprimaryIstrack = 0;
941  fprimaryIsshower = 1;
942 
943  fprimaryID = thisShower->ID();
944  fprimaryLength = thisShower->Length();
945  fprimaryShowerBestPlane = thisShower->best_plane();
946  fprimaryOpeningAngle = thisShower->OpenAngle();
947  fprimaryStartPosition[0] = thisShower->ShowerStart().X();
948  fprimaryStartPosition[1] = thisShower->ShowerStart().Y();
949  fprimaryStartPosition[2] = thisShower->ShowerStart().Z();
950  fprimaryStartDirection[0] = thisShower->Direction().X();
951  fprimaryStartDirection[1] = thisShower->Direction().Y();
952  fprimaryStartDirection[2] = thisShower->Direction().Z();
953  for(size_t k = 0; k < (thisShower->Energy()).size() && k<3; k++)
954  fprimaryShowerEnergy[k] = thisShower->Energy()[k];
955  for(size_t k = 0; k < (thisShower->MIPEnergy()).size() && k<3; k++)
956  fprimaryShowerMIPEnergy[k] = thisShower->MIPEnergy()[k];
957  for(size_t k = 0; k < (thisShower->dEdx()).size() && k<3; k++)
958  fprimaryShowerdEdx[k] = thisShower->dEdx()[k];
959 
960  //const simb::MCParticle* mcparticle = truthUtil.GetMCParticleFromRecoShower(*thisShower, evt, fShowerTag);
961  mcparticle = truthUtil.GetMCParticleFromRecoShower(clockData, *thisShower, evt, fShowerTag);
962  } // end is shower
963  else{
964  if(fVerbose > 0){
965  std::cout << "INFO::Primary pfParticle is not track or shower!" << std::endl;
966  }
967  return;
968  }
969 
970  if(mcparticle != 0x0){
971  fprimary_truth_Origin = pi_service->TrackIdToMCTruth_P(mcparticle->TrackId())->Origin();
972  fprimary_truth_TrackId = mcparticle->TrackId();
973  fprimary_truth_Pdg = mcparticle->PdgCode();
974  fprimary_truth_Mother = mcparticle->Mother();
975  fprimary_truth_StartPosition[0] = mcparticle->Vx();
976  fprimary_truth_StartPosition[1] = mcparticle->Vy();
977  fprimary_truth_StartPosition[2] = mcparticle->Vz();
978  fprimary_truth_StartPosition[3] = mcparticle->T();
979  fprimary_truth_EndPosition[0] = mcparticle->EndX();
980  fprimary_truth_EndPosition[1] = mcparticle->EndY();
981  fprimary_truth_EndPosition[2] = mcparticle->EndZ();
982  fprimary_truth_EndPosition[3] = mcparticle->EndT();
983  fprimary_truth_P = mcparticle->P();
984  fprimary_truth_Momentum[0] = mcparticle->Px();
985  fprimary_truth_Momentum[1] = mcparticle->Py();
986  fprimary_truth_Momentum[2] = mcparticle->Pz();
987  fprimary_truth_Momentum[3] = mcparticle->E();
988  fprimary_truth_Pt = mcparticle->Pt();
989  fprimary_truth_Mass = mcparticle->Mass();
990  fprimary_truth_EndMomentum[0] = mcparticle->EndPx();
991  fprimary_truth_EndMomentum[1] = mcparticle->EndPy();
992  fprimary_truth_EndMomentum[2] = mcparticle->EndPz();
993  fprimary_truth_EndMomentum[3] = mcparticle->EndE();
994  fprimary_truth_Theta = mcparticle->Momentum().Theta();
995  fprimary_truth_Phi = mcparticle->Momentum().Phi();
996  //fprimary_truth_NDaughters = mcparticle->NumberDaughters();
997  fprimary_truth_TotalLength = mcparticle->Trajectory().TotalLength();
998  fprimary_truth_Process = truthUtil.GetProcessKey(mcparticle->Process());
999  fprimary_truth_EndProcess = truthUtil.GetProcessKey(mcparticle->EndProcess());
1000  fprimary_truth_EkinAtVertex_notcorrected = truthUtil.GetKinEnergyAtVertex(*mcparticle);
1001  double lastpointEloss = truthUtil.GetDepEnergyAtLastTrajPoint(*mcparticle,fActiveTPCBoundsX[0], fActiveTPCBoundsX[1], fActiveTPCBoundsY[0], fActiveTPCBoundsY[1], fActiveTPCBoundsZ[0], fActiveTPCBoundsZ[1]);
1002  fprimary_truth_EkinAtVertex = truthUtil.GetKinEnergyAtVertex(*mcparticle,lastpointEloss);
1003 
1004  // Get the position and kinematics in the TPC active volume
1005  int trajintpcactivevol = truthUtil.GetFirstTrajectoryPointInTPCActiveVolume(*mcparticle, fActiveTPCBoundsX[0], fActiveTPCBoundsX[1], fActiveTPCBoundsY[0], fActiveTPCBoundsY[1], fActiveTPCBoundsZ[0], fActiveTPCBoundsZ[1]);
1006  if(trajintpcactivevol > 0){
1007  fprimary_truth_Pos_InTPCActive[0] = mcparticle->Vx(trajintpcactivevol);
1008  fprimary_truth_Pos_InTPCActive[1] = mcparticle->Vy(trajintpcactivevol);
1009  fprimary_truth_Pos_InTPCActive[2] = mcparticle->Vz(trajintpcactivevol);
1010  fprimary_truth_Pos_InTPCActive[3] = mcparticle->T(trajintpcactivevol);
1011  fprimary_truth_Momentum_InTPCActive[0] = mcparticle->Px(trajintpcactivevol);
1012  fprimary_truth_Momentum_InTPCActive[1] = mcparticle->Py(trajintpcactivevol);
1013  fprimary_truth_Momentum_InTPCActive[2] = mcparticle->Pz(trajintpcactivevol);
1014  fprimary_truth_Momentum_InTPCActive[3] = mcparticle->E(trajintpcactivevol);
1015  fprimary_truth_P_InTPCActive = mcparticle->P(trajintpcactivevol);
1016  fprimary_truth_Pt_InTPCActive = mcparticle->Pt(trajintpcactivevol);
1017  fprimary_truth_Theta_InTPCActive = mcparticle->Momentum(trajintpcactivevol).Theta();
1018  fprimary_truth_Phi_InTPCActive = mcparticle->Momentum(trajintpcactivevol).Phi();
1019  fprimary_truth_TotalLength_InTPCActive = truthUtil.GetMCParticleLengthInTPCActiveVolume(*mcparticle, fActiveTPCBoundsX[0], fActiveTPCBoundsX[1], fActiveTPCBoundsY[0], fActiveTPCBoundsY[1], fActiveTPCBoundsZ[0], fActiveTPCBoundsZ[1]);
1020  fprimary_truth_KinEnergy_InTPCActive = std::sqrt(mcparticle->P(trajintpcactivevol)*mcparticle->P(trajintpcactivevol) + mcparticle->Mass()*mcparticle->Mass()) - mcparticle->Mass();
1021  }
1022 
1023  if(fbeamtrackID != -999 && fprimary_truth_TrackId == fbeamtrackID)
1024  fprimary_truth_Isbeammatched = 1;
1025  else
1026  fprimary_truth_Isbeammatched = 0;
1027 
1028  // If this beam particle has interacted in the detector, find its daugthers.
1029  std::vector<art::Ptr<simb::MCParticle> > mcprimarypartilelist;
1030  auto mcprimaryHandle = evt.getHandle<std::vector<simb::MCParticle> >(fSimulationTag);
1031  if (mcprimaryHandle) art::fill_ptr_vector(mcprimarypartilelist, mcprimaryHandle);
1032 
1033  for(size_t p=0; p < mcprimarypartilelist.size(); p++){
1034  auto & primary_mcparticle = mcprimarypartilelist[p];
1035 
1036  if(primary_mcparticle->Mother() == mcparticle->TrackId()){
1037  fprimary_truthdaughter_TrackId[fprimary_truth_NDAUGTHERS] = primary_mcparticle->TrackId();
1038  fprimary_truthdaughter_Pdg[fprimary_truth_NDAUGTHERS] = primary_mcparticle->PdgCode();
1039  fprimary_truthdaughter_Mother[fprimary_truth_NDAUGTHERS] = primary_mcparticle->Mother();
1040  fprimary_truthdaughter_StartPosition[fprimary_truth_NDAUGTHERS][0] = primary_mcparticle->Vx();
1041  fprimary_truthdaughter_StartPosition[fprimary_truth_NDAUGTHERS][1] = primary_mcparticle->Vy();
1042  fprimary_truthdaughter_StartPosition[fprimary_truth_NDAUGTHERS][2] = primary_mcparticle->Vz();
1043  fprimary_truthdaughter_StartPosition[fprimary_truth_NDAUGTHERS][3] = primary_mcparticle->T();
1044  fprimary_truthdaughter_EndPosition[fprimary_truth_NDAUGTHERS][0] = primary_mcparticle->EndX();
1045  fprimary_truthdaughter_EndPosition[fprimary_truth_NDAUGTHERS][1] = primary_mcparticle->EndY();
1046  fprimary_truthdaughter_EndPosition[fprimary_truth_NDAUGTHERS][2] = primary_mcparticle->EndZ();
1047  fprimary_truthdaughter_EndPosition[fprimary_truth_NDAUGTHERS][3] = primary_mcparticle->EndT();
1048  fprimary_truthdaughter_P[fprimary_truth_NDAUGTHERS] = primary_mcparticle->P();
1049  fprimary_truthdaughter_Momentum[fprimary_truth_NDAUGTHERS][0] = primary_mcparticle->Px();
1050  fprimary_truthdaughter_Momentum[fprimary_truth_NDAUGTHERS][1] = primary_mcparticle->Py();
1051  fprimary_truthdaughter_Momentum[fprimary_truth_NDAUGTHERS][2] = primary_mcparticle->Pz();
1052  fprimary_truthdaughter_Momentum[fprimary_truth_NDAUGTHERS][3] = primary_mcparticle->E();
1053  fprimary_truthdaughter_Pt[fprimary_truth_NDAUGTHERS] = primary_mcparticle->Pt();
1054  fprimary_truthdaughter_Mass[fprimary_truth_NDAUGTHERS] = primary_mcparticle->Mass();
1055  fprimary_truthdaughter_EndMomentum[fprimary_truth_NDAUGTHERS][0] = primary_mcparticle->EndPx();
1056  fprimary_truthdaughter_EndMomentum[fprimary_truth_NDAUGTHERS][1] = primary_mcparticle->EndPy();
1057  fprimary_truthdaughter_EndMomentum[fprimary_truth_NDAUGTHERS][2] = primary_mcparticle->EndPz();
1058  fprimary_truthdaughter_EndMomentum[fprimary_truth_NDAUGTHERS][3] = primary_mcparticle->EndE();
1059  fprimary_truthdaughter_Theta[fprimary_truth_NDAUGTHERS] = primary_mcparticle->Momentum().Theta();
1060  fprimary_truthdaughter_Phi[fprimary_truth_NDAUGTHERS] = primary_mcparticle->Momentum().Phi();
1061  fprimary_truthdaughter_TotalLength[fprimary_truth_NDAUGTHERS] = primary_mcparticle->Trajectory().TotalLength();
1062  fprimary_truthdaughter_Process[fprimary_truth_NDAUGTHERS] = truthUtil.GetProcessKey(primary_mcparticle->Process());
1063  fprimary_truthdaughter_EndProcess[fprimary_truth_NDAUGTHERS] = truthUtil.GetProcessKey(primary_mcparticle->EndProcess());
1064 
1065  fprimary_truth_NDAUGTHERS++;
1066  if(fprimary_truth_NDAUGTHERS > NMAXTRUTHDAUGTHERS) break;
1067 
1068  if(truthUtil.GetProcessKey(primary_mcparticle->EndProcess()) == 30){ // Decay particle
1069  for(size_t pp=0; pp < mcprimarypartilelist.size(); pp++){
1070  auto & decay_primary_mcparticle = mcprimarypartilelist[pp];
1071 
1072  if(decay_primary_mcparticle->Mother() == primary_mcparticle->TrackId()){
1073  fprimary_truthdecaydaughter_TrackId[fprimary_truth_NDECAYDAUGTHERS] = decay_primary_mcparticle->TrackId();
1074  fprimary_truthdecaydaughter_Pdg[fprimary_truth_NDECAYDAUGTHERS] = decay_primary_mcparticle->PdgCode();
1075  fprimary_truthdecaydaughter_Mother[fprimary_truth_NDECAYDAUGTHERS] = decay_primary_mcparticle->Mother();
1076  fprimary_truthdecaydaughter_StartPosition[fprimary_truth_NDECAYDAUGTHERS][0] = decay_primary_mcparticle->Vx();
1077  fprimary_truthdecaydaughter_StartPosition[fprimary_truth_NDECAYDAUGTHERS][1] = decay_primary_mcparticle->Vy();
1078  fprimary_truthdecaydaughter_StartPosition[fprimary_truth_NDECAYDAUGTHERS][2] = decay_primary_mcparticle->Vz();
1079  fprimary_truthdecaydaughter_StartPosition[fprimary_truth_NDECAYDAUGTHERS][3] = decay_primary_mcparticle->T();
1080  fprimary_truthdecaydaughter_EndPosition[fprimary_truth_NDECAYDAUGTHERS][0] = decay_primary_mcparticle->EndX();
1081  fprimary_truthdecaydaughter_EndPosition[fprimary_truth_NDECAYDAUGTHERS][1] = decay_primary_mcparticle->EndY();
1082  fprimary_truthdecaydaughter_EndPosition[fprimary_truth_NDECAYDAUGTHERS][2] = decay_primary_mcparticle->EndZ();
1083  fprimary_truthdecaydaughter_EndPosition[fprimary_truth_NDECAYDAUGTHERS][3] = decay_primary_mcparticle->EndT();
1084  fprimary_truthdecaydaughter_P[fprimary_truth_NDECAYDAUGTHERS] = decay_primary_mcparticle->P();
1085  fprimary_truthdecaydaughter_Momentum[fprimary_truth_NDECAYDAUGTHERS][0] = decay_primary_mcparticle->Px();
1086  fprimary_truthdecaydaughter_Momentum[fprimary_truth_NDECAYDAUGTHERS][1] = decay_primary_mcparticle->Py();
1087  fprimary_truthdecaydaughter_Momentum[fprimary_truth_NDECAYDAUGTHERS][2] = decay_primary_mcparticle->Pz();
1088  fprimary_truthdecaydaughter_Momentum[fprimary_truth_NDECAYDAUGTHERS][3] = decay_primary_mcparticle->E();
1089  fprimary_truthdecaydaughter_Pt[fprimary_truth_NDECAYDAUGTHERS] = decay_primary_mcparticle->Pt();
1090  fprimary_truthdecaydaughter_Mass[fprimary_truth_NDECAYDAUGTHERS] = decay_primary_mcparticle->Mass();
1091  fprimary_truthdecaydaughter_EndMomentum[fprimary_truth_NDECAYDAUGTHERS][0] = decay_primary_mcparticle->EndPx();
1092  fprimary_truthdecaydaughter_EndMomentum[fprimary_truth_NDECAYDAUGTHERS][1] = decay_primary_mcparticle->EndPy();
1093  fprimary_truthdecaydaughter_EndMomentum[fprimary_truth_NDECAYDAUGTHERS][2] = decay_primary_mcparticle->EndPz();
1094  fprimary_truthdecaydaughter_EndMomentum[fprimary_truth_NDECAYDAUGTHERS][3] = decay_primary_mcparticle->EndE();
1095  fprimary_truthdecaydaughter_Theta[fprimary_truth_NDECAYDAUGTHERS] = decay_primary_mcparticle->Momentum().Theta();
1096  fprimary_truthdecaydaughter_Phi[fprimary_truth_NDECAYDAUGTHERS] = decay_primary_mcparticle->Momentum().Phi();
1097  fprimary_truthdecaydaughter_TotalLength[fprimary_truth_NDECAYDAUGTHERS] = decay_primary_mcparticle->Trajectory().TotalLength();
1098  fprimary_truthdecaydaughter_Process[fprimary_truth_NDECAYDAUGTHERS] = truthUtil.GetProcessKey(decay_primary_mcparticle->Process());
1099  fprimary_truthdecaydaughter_EndProcess[fprimary_truth_NDECAYDAUGTHERS] = truthUtil.GetProcessKey(decay_primary_mcparticle->EndProcess());
1100 
1101  fprimary_truth_NDECAYDAUGTHERS++;
1102  if(fprimary_truth_NDECAYDAUGTHERS > NMAXTRUTHDAUGTHERS) break;
1103  }
1104  }
1105  } // decay particle
1106  }
1107  }
1108 
1109  if(fVerbose > 1){
1110  std::cout << "INFO::Process = " << (mcparticle->Process()).c_str() << " , End process = " << (mcparticle->EndProcess()).c_str()
1111  << " , track ID = " << mcparticle->TrackId()
1112  << std::endl;
1113  }
1114  }
1115  else{
1116  if(fVerbose > 0){
1117  std::cout << "INFO::No MCParticle for primary found!" << std::endl;
1118  }
1119  }
1120 
1121 }
1122 
1123 // -----------------------------------------------------------------------------
1124 void protoana::ProtoDUNEAnalysisTree::FillPrimaryDaughterPFParticle(art::Event const & evt,
1125  detinfo::DetectorClocksData const& clockData,
1126  detinfo::DetectorPropertiesData const& detProp,
1127  const recob::PFParticle* daughterParticle,
1128  int daughterID){
1129 
1130  const recob::Track* daughterTrack = pfpUtil.GetPFParticleTrack(*daughterParticle,evt, fPFParticleTag,fTrackerTag);
1131  const recob::Shower* daughterShower = pfpUtil.GetPFParticleShower(*daughterParticle,evt,fPFParticleTag,fShowerTag);
1132  const simb::MCParticle* mcdaughterparticle = NULL;
1133 
1134  if(daughterTrack != 0x0){
1135  fdaughterIstrack[fNDAUGHTERS] = 1;
1136  fdaughterIsshower[fNDAUGHTERS] = 0;
1137  fdaughterTheta[fNDAUGHTERS] = daughterTrack->Theta();
1138  fdaughterPhi[fNDAUGHTERS] = daughterTrack->Phi();
1139  fdaughterLength[fNDAUGHTERS] = daughterTrack->Length();
1140  fdaughterMomentum[fNDAUGHTERS] = daughterTrack->StartMomentum();
1141  fdaughterEndMomentum[fNDAUGHTERS] = daughterTrack->EndMomentum();
1142  fdaughterStartPosition[fNDAUGHTERS][0] = daughterTrack->Trajectory().Start().X();
1143  fdaughterStartPosition[fNDAUGHTERS][1] = daughterTrack->Trajectory().Start().Y();
1144  fdaughterStartPosition[fNDAUGHTERS][2] = daughterTrack->Trajectory().Start().Z();
1145  fdaughterEndPosition[fNDAUGHTERS][0] = daughterTrack->Trajectory().End().X();
1146  fdaughterEndPosition[fNDAUGHTERS][1] = daughterTrack->Trajectory().End().Y();
1147  fdaughterEndPosition[fNDAUGHTERS][2] = daughterTrack->Trajectory().End().Z();
1148  fdaughterStartDirection[fNDAUGHTERS][0] = daughterTrack->Trajectory().StartDirection().X();
1149  fdaughterStartDirection[fNDAUGHTERS][1] = daughterTrack->Trajectory().StartDirection().Y();
1150  fdaughterStartDirection[fNDAUGHTERS][2] = daughterTrack->Trajectory().StartDirection().Z();
1151  fdaughterEndDirection[fNDAUGHTERS][0] = daughterTrack->Trajectory().EndDirection().X();
1152  fdaughterEndDirection[fNDAUGHTERS][1] = daughterTrack->Trajectory().EndDirection().Y();
1153  fdaughterEndDirection[fNDAUGHTERS][2] = daughterTrack->Trajectory().EndDirection().Z();
1154 
1155  fdaughterMomentumByRangeMuon[fNDAUGHTERS] = trmom.GetTrackMomentum(daughterTrack->Length(),13);
1156  fdaughterMomentumByRangeProton[fNDAUGHTERS] = trmom.GetTrackMomentum(daughterTrack->Length(),2212);
1157 
1158  // Calorimetry
1159  std::vector<anab::Calorimetry> daughtercalovector = trackUtil.GetRecoTrackCalorimetry(*daughterTrack, evt, fTrackerTag, fCalorimetryTag);
1160  if(daughtercalovector.size() != 3 && fVerbose > 0)
1161  std::cerr << "WARNING::Calorimetry vector size for daughter is = " << daughtercalovector.size() << std::endl;
1162 
1163  for(size_t k = 0; k < daughtercalovector.size() && k<3; k++){
1164  int plane = daughtercalovector[k].PlaneID().Plane;
1165  if(plane < 0) continue;
1166  if(plane > 2) continue;
1167  fdaughterKineticEnergy[fNDAUGHTERS][plane] = daughtercalovector[k].KineticEnergy();
1168  fdaughterRange[fNDAUGHTERS][plane] = daughtercalovector[k].Range();
1169  fdaughterTrkPitchC[fNDAUGHTERS][plane] = daughtercalovector[k].TrkPitchC();
1170  }
1171 
1172  // PID
1173  std::vector<anab::ParticleID> daughterpids = trackUtil.GetRecoTrackPID(*daughterTrack, evt, fTrackerTag, fParticleIDTag);
1174  if(daughterpids.size() != 3 && fVerbose > 0)
1175  std::cerr << "WARNING::PID vector size for daughter is = " << daughterpids.size() << std::endl;
1176 
1177  for(size_t k = 0; k < daughterpids.size() && k<3; k++){
1178  int plane = daughterpids[k].PlaneID().Plane;
1179  if(plane < 0) continue;
1180  if(plane > 2) continue;
1181  fdaughterPID_Pdg[fNDAUGHTERS][plane] = daughterpids[plane].Pdg();
1182  fdaughterPID_Ndf[fNDAUGHTERS][plane] = daughterpids[plane].Ndf();
1183  fdaughterPID_MinChi2[fNDAUGHTERS][plane] = daughterpids[plane].MinChi2();
1184  fdaughterPID_DeltaChi2[fNDAUGHTERS][plane] = daughterpids[plane].DeltaChi2();
1185  fdaughterPID_Chi2Proton[fNDAUGHTERS][plane] = daughterpids[plane].Chi2Proton();
1186  fdaughterPID_Chi2Kaon[fNDAUGHTERS][plane] = daughterpids[plane].Chi2Kaon();
1187  fdaughterPID_Chi2Pion[fNDAUGHTERS][plane] = daughterpids[plane].Chi2Pion();
1188  fdaughterPID_Chi2Muon[fNDAUGHTERS][plane] = daughterpids[plane].Chi2Muon();
1189  fdaughterPID_MissingE[fNDAUGHTERS][plane] = daughterpids[plane].MissingE();
1190  fdaughterPID_MissingEavg[fNDAUGHTERS][plane] = daughterpids[plane].MissingEavg();
1191  fdaughterPID_PIDA[fNDAUGHTERS][plane] = daughterpids[plane].PIDA();
1192  }
1193 
1194  // Get the true mc particle
1195  //const simb::MCParticle* mcdaughterparticle = truthUtil.GetMCParticleFromRecoTrack(*daughterTrack, evt, fTrackerTag);
1196  mcdaughterparticle = truthUtil.GetMCParticleFromRecoTrack(clockData, *daughterTrack, evt, fTrackerTag);
1197  }
1198  else if(daughterShower != 0x0){
1199  fdaughterIstrack[fNDAUGHTERS] = 0;
1200  fdaughterIsshower[fNDAUGHTERS] = 1;
1201  fdaughterLength[fNDAUGHTERS] = daughterShower->Length();
1202  fdaughterShowerBestPlane[fNDAUGHTERS] = daughterShower->best_plane();
1203  fdaughterOpeningAngle[fNDAUGHTERS] = daughterShower->OpenAngle();
1204  fdaughterStartPosition[fNDAUGHTERS][0] = daughterShower->ShowerStart().X();
1205  fdaughterStartPosition[fNDAUGHTERS][1] = daughterShower->ShowerStart().Y();
1206  fdaughterStartPosition[fNDAUGHTERS][2] = daughterShower->ShowerStart().Z();
1207  fdaughterStartDirection[fNDAUGHTERS][0] = daughterShower->Direction().X();
1208  fdaughterStartDirection[fNDAUGHTERS][1] = daughterShower->Direction().Y();
1209  fdaughterStartDirection[fNDAUGHTERS][2] = daughterShower->Direction().Z();
1210  for(size_t k = 0; k < (daughterShower->Energy()).size() && k<3; k++)
1211  fdaughterShowerEnergy[fNDAUGHTERS][k] = daughterShower->Energy()[k];
1212  for(size_t k = 0; k < (daughterShower->MIPEnergy()).size() && k<3; k++)
1213  fdaughterShowerMIPEnergy[fNDAUGHTERS][k] = daughterShower->MIPEnergy()[k];
1214  for(size_t k = 0; k < (daughterShower->dEdx()).size() && k<3; k++)
1215  fdaughterShowerdEdx[fNDAUGHTERS][k] = daughterShower->dEdx()[k];
1216 
1217  // Get the true mc particle
1218  //const simb::MCParticle* mcdaughterparticle = truthUtil.GetMCParticleFromRecoShower(*daughterShower, evt, fShowerTag);
1219  mcdaughterparticle = truthUtil.GetMCParticleFromRecoShower(clockData, *daughterShower, evt, fShowerTag);
1220  }
1221  else{
1222  if(fVerbose > 0){
1223  std::cout << "INFO::Daughter pfParticle is not track or shower!" << std::endl;
1224  }
1225  //return;
1226  }
1227 
1228  fdaughterID[fNDAUGHTERS] = daughterID;
1229  // NHits associated with this pfParticle
1230  fdaughterNHits[fNDAUGHTERS] = (pfpUtil.GetPFParticleHits(*daughterParticle,evt,fPFParticleTag)).size();
1231  // T0
1232  std::vector<anab::T0> pfdaughterT0vec = pfpUtil.GetPFParticleT0(*daughterParticle,evt,fPFParticleTag);
1233  if(!pfdaughterT0vec.empty())
1234  fdaughterT0[fNDAUGHTERS] = pfdaughterT0vec[0].Time();
1235 
1236  if(mcdaughterparticle != 0x0){
1237  fdaughter_truth_TrackId[fNDAUGHTERS] = mcdaughterparticle->TrackId();
1238  fdaughter_truth_Pdg[fNDAUGHTERS] = mcdaughterparticle->PdgCode();
1239  fdaughter_truth_Mother[fNDAUGHTERS] = mcdaughterparticle->Mother();
1240  fdaughter_truth_StartPosition[fNDAUGHTERS][0] = mcdaughterparticle->Vx();
1241  fdaughter_truth_StartPosition[fNDAUGHTERS][1] = mcdaughterparticle->Vy();
1242  fdaughter_truth_StartPosition[fNDAUGHTERS][2] = mcdaughterparticle->Vz();
1243  fdaughter_truth_StartPosition[fNDAUGHTERS][3] = mcdaughterparticle->T();
1244  fdaughter_truth_EndPosition[fNDAUGHTERS][0] = mcdaughterparticle->EndX();
1245  fdaughter_truth_EndPosition[fNDAUGHTERS][1] = mcdaughterparticle->EndY();
1246  fdaughter_truth_EndPosition[fNDAUGHTERS][2] = mcdaughterparticle->EndZ();
1247  fdaughter_truth_EndPosition[fNDAUGHTERS][3] = mcdaughterparticle->EndT();
1248  fdaughter_truth_P[fNDAUGHTERS] = mcdaughterparticle->P();
1249  fdaughter_truth_Momentum[fNDAUGHTERS][0] = mcdaughterparticle->Px();
1250  fdaughter_truth_Momentum[fNDAUGHTERS][1] = mcdaughterparticle->Py();
1251  fdaughter_truth_Momentum[fNDAUGHTERS][2] = mcdaughterparticle->Pz();
1252  fdaughter_truth_Momentum[fNDAUGHTERS][3] = mcdaughterparticle->E();
1253  fdaughter_truth_Pt[fNDAUGHTERS] = mcdaughterparticle->Pt();
1254  fdaughter_truth_Mass[fNDAUGHTERS] = mcdaughterparticle->Mass();
1255  fdaughter_truth_EndMomentum[fNDAUGHTERS][0] = mcdaughterparticle->EndPx();
1256  fdaughter_truth_EndMomentum[fNDAUGHTERS][1] = mcdaughterparticle->EndPy();
1257  fdaughter_truth_EndMomentum[fNDAUGHTERS][2] = mcdaughterparticle->EndPz();
1258  fdaughter_truth_EndMomentum[fNDAUGHTERS][3] = mcdaughterparticle->EndE();
1259  fdaughter_truth_Theta[fNDAUGHTERS] = mcdaughterparticle->Momentum().Theta();
1260  fdaughter_truth_Phi[fNDAUGHTERS] = mcdaughterparticle->Momentum().Phi();
1261  fdaughter_truth_TotalLength[fNDAUGHTERS] = mcdaughterparticle->Trajectory().TotalLength();
1262  fdaughter_truth_Process[fNDAUGHTERS] = truthUtil.GetProcessKey(mcdaughterparticle->Process());
1263  fdaughter_truth_EndProcess[fNDAUGHTERS] = truthUtil.GetProcessKey(mcdaughterparticle->EndProcess());
1264  }
1265  else{
1266  if(fVerbose > 0){
1267  std::cout << "INFO::No MCParticle for daughter found!" << std::endl;
1268  }
1269  }
1270 
1271  // Increment counter
1272  fNDAUGHTERS++;
1273 
1274 }
1275 
1276 // -----------------------------------------------------------------------------
1277 void protoana::ProtoDUNEAnalysisTree::FillPrimaryGrandDaughterPFParticle(art::Event const & evt,
1278  detinfo::DetectorClocksData const& clockData,
1279  detinfo::DetectorPropertiesData const& detProp,
1280  const recob::PFParticle* gdaughterParticle,
1281  int daughterID, int gdaughterID){
1282 
1283  const recob::Track* gdaughterTrack = pfpUtil.GetPFParticleTrack(*gdaughterParticle,evt, fPFParticleTag,fTrackerTag);
1284  const recob::Shower* gdaughterShower = pfpUtil.GetPFParticleShower(*gdaughterParticle,evt,fPFParticleTag,fShowerTag);
1285  const simb::MCParticle* mcdaughterparticle = NULL;
1286 
1287  if(gdaughterTrack != 0x0){
1288  fgranddaughterIstrack[fNGRANDDAUGHTERS] = 1;
1289  fgranddaughterIsshower[fNGRANDDAUGHTERS] = 0;
1290  fgranddaughterTheta[fNGRANDDAUGHTERS] = gdaughterTrack->Theta();
1291  fgranddaughterPhi[fNGRANDDAUGHTERS] = gdaughterTrack->Phi();
1292  fgranddaughterLength[fNGRANDDAUGHTERS] = gdaughterTrack->Length();
1293  fgranddaughterMomentum[fNGRANDDAUGHTERS] = gdaughterTrack->StartMomentum();
1294  fgranddaughterEndMomentum[fNGRANDDAUGHTERS] = gdaughterTrack->EndMomentum();
1295  fgranddaughterStartPosition[fNGRANDDAUGHTERS][0] = gdaughterTrack->Trajectory().Start().X();
1296  fgranddaughterStartPosition[fNGRANDDAUGHTERS][1] = gdaughterTrack->Trajectory().Start().Y();
1297  fgranddaughterStartPosition[fNGRANDDAUGHTERS][2] = gdaughterTrack->Trajectory().Start().Z();
1298  fgranddaughterEndPosition[fNGRANDDAUGHTERS][0] = gdaughterTrack->Trajectory().End().X();
1299  fgranddaughterEndPosition[fNGRANDDAUGHTERS][1] = gdaughterTrack->Trajectory().End().Y();
1300  fgranddaughterEndPosition[fNGRANDDAUGHTERS][2] = gdaughterTrack->Trajectory().End().Z();
1301  fgranddaughterStartDirection[fNGRANDDAUGHTERS][0] = gdaughterTrack->Trajectory().StartDirection().X();
1302  fgranddaughterStartDirection[fNGRANDDAUGHTERS][1] = gdaughterTrack->Trajectory().StartDirection().Y();
1303  fgranddaughterStartDirection[fNGRANDDAUGHTERS][2] = gdaughterTrack->Trajectory().StartDirection().Z();
1304  fgranddaughterEndDirection[fNGRANDDAUGHTERS][0] = gdaughterTrack->Trajectory().EndDirection().X();
1305  fgranddaughterEndDirection[fNGRANDDAUGHTERS][1] = gdaughterTrack->Trajectory().EndDirection().Y();
1306  fgranddaughterEndDirection[fNGRANDDAUGHTERS][2] = gdaughterTrack->Trajectory().EndDirection().Z();
1307 
1308  fgranddaughterMomentumByRangeMuon[fNGRANDDAUGHTERS] = trmom.GetTrackMomentum(gdaughterTrack->Length(),13);
1309  fgranddaughterMomentumByRangeProton[fNGRANDDAUGHTERS] = trmom.GetTrackMomentum(gdaughterTrack->Length(),2212);
1310 
1311  // Calorimetry
1312  std::vector<anab::Calorimetry> daughtercalovector = trackUtil.GetRecoTrackCalorimetry(*gdaughterTrack, evt, fTrackerTag, fCalorimetryTag);
1313  if(daughtercalovector.size() != 3 && fVerbose > 0)
1314  std::cerr << "WARNING::Calorimetry vector size for grand-daughter is = " << daughtercalovector.size() << std::endl;
1315 
1316  for(size_t k = 0; k < daughtercalovector.size() && k<3; k++){
1317  int plane = daughtercalovector[k].PlaneID().Plane;
1318  if(plane < 0) continue;
1319  if(plane > 2) continue;
1320  fgranddaughterKineticEnergy[fNGRANDDAUGHTERS][plane] = daughtercalovector[k].KineticEnergy();
1321  fgranddaughterRange[fNGRANDDAUGHTERS][plane] = daughtercalovector[k].Range();
1322  fgranddaughterTrkPitchC[fNGRANDDAUGHTERS][plane] = daughtercalovector[k].TrkPitchC();
1323  }
1324 
1325  // PID
1326  std::vector<anab::ParticleID> daughterpids = trackUtil.GetRecoTrackPID(*gdaughterTrack, evt, fTrackerTag, fParticleIDTag);
1327  if(daughterpids.size() != 3 && fVerbose > 0)
1328  std::cerr << "WARNING::PID vector size for grand-daughter is = " << daughterpids.size() << std::endl;
1329 
1330  for(size_t k = 0; k < daughterpids.size() && k<3; k++){
1331  int plane = daughterpids[k].PlaneID().Plane;
1332  if(plane < 0) continue;
1333  if(plane > 2) continue;
1334  fgranddaughterPID_Pdg[fNGRANDDAUGHTERS][plane] = daughterpids[plane].Pdg();
1335  fgranddaughterPID_Ndf[fNGRANDDAUGHTERS][plane] = daughterpids[plane].Ndf();
1336  fgranddaughterPID_MinChi2[fNGRANDDAUGHTERS][plane] = daughterpids[plane].MinChi2();
1337  fgranddaughterPID_DeltaChi2[fNGRANDDAUGHTERS][plane] = daughterpids[plane].DeltaChi2();
1338  fgranddaughterPID_Chi2Proton[fNGRANDDAUGHTERS][plane] = daughterpids[plane].Chi2Proton();
1339  fgranddaughterPID_Chi2Kaon[fNGRANDDAUGHTERS][plane] = daughterpids[plane].Chi2Kaon();
1340  fgranddaughterPID_Chi2Pion[fNGRANDDAUGHTERS][plane] = daughterpids[plane].Chi2Pion();
1341  fgranddaughterPID_Chi2Muon[fNGRANDDAUGHTERS][plane] = daughterpids[plane].Chi2Muon();
1342  fgranddaughterPID_MissingE[fNGRANDDAUGHTERS][plane] = daughterpids[plane].MissingE();
1343  fgranddaughterPID_MissingEavg[fNGRANDDAUGHTERS][plane] = daughterpids[plane].MissingEavg();
1344  fgranddaughterPID_PIDA[fNGRANDDAUGHTERS][plane] = daughterpids[plane].PIDA();
1345  }
1346 
1347  // Get the true mc particle
1348  //const simb::MCParticle* mcdaughterparticle = truthUtil.GetMCParticleFromRecoTrack(*gdaughterTrack, evt, fTrackerTag);
1349  mcdaughterparticle = truthUtil.GetMCParticleFromRecoTrack(clockData, *gdaughterTrack, evt, fTrackerTag);
1350  }
1351  else if(gdaughterShower != 0x0){
1352  fgranddaughterIstrack[fNGRANDDAUGHTERS] = 0;
1353  fgranddaughterIsshower[fNGRANDDAUGHTERS] = 1;
1354  fgranddaughterLength[fNGRANDDAUGHTERS] = gdaughterShower->Length();
1355  fgranddaughterShowerBestPlane[fNGRANDDAUGHTERS] = gdaughterShower->best_plane();
1356  fgranddaughterOpeningAngle[fNGRANDDAUGHTERS] = gdaughterShower->OpenAngle();
1357  fgranddaughterStartPosition[fNGRANDDAUGHTERS][0] = gdaughterShower->ShowerStart().X();
1358  fgranddaughterStartPosition[fNGRANDDAUGHTERS][1] = gdaughterShower->ShowerStart().Y();
1359  fgranddaughterStartPosition[fNGRANDDAUGHTERS][2] = gdaughterShower->ShowerStart().Z();
1360  fgranddaughterStartDirection[fNGRANDDAUGHTERS][0] = gdaughterShower->Direction().X();
1361  fgranddaughterStartDirection[fNGRANDDAUGHTERS][1] = gdaughterShower->Direction().Y();
1362  fgranddaughterStartDirection[fNGRANDDAUGHTERS][2] = gdaughterShower->Direction().Z();
1363  for(size_t k = 0; k < (gdaughterShower->Energy()).size() && k<3; k++)
1364  fgranddaughterShowerEnergy[fNGRANDDAUGHTERS][k] = gdaughterShower->Energy()[k];
1365  for(size_t k = 0; k < (gdaughterShower->MIPEnergy()).size() && k<3; k++)
1366  fgranddaughterShowerMIPEnergy[fNGRANDDAUGHTERS][k] = gdaughterShower->MIPEnergy()[k];
1367  for(size_t k = 0; k < (gdaughterShower->dEdx()).size() && k<3; k++)
1368  fgranddaughterShowerdEdx[fNGRANDDAUGHTERS][k] = gdaughterShower->dEdx()[k];
1369 
1370  // Get the true mc particle
1371  //const simb::MCParticle* mcdaughterparticle = truthUtil.GetMCParticleFromRecoShower(*gdaughterShower, evt, fShowerTag);
1372  mcdaughterparticle = truthUtil.GetMCParticleFromRecoShower(clockData, *gdaughterShower, evt, fShowerTag);
1373  }
1374  else{
1375  if(fVerbose > 0){
1376  std::cout << "INFO::GrandDaughter pfParticle is not track or shower!" << std::endl;
1377  }
1378  //return;
1379  }
1380 
1381  fgranddaughterID[fNGRANDDAUGHTERS] = gdaughterID;
1382  fgranddaughterMotherID[fNGRANDDAUGHTERS] = daughterID;
1383  // NHits associated with this pfParticle
1384  fgranddaughterNHits[fNGRANDDAUGHTERS] = (pfpUtil.GetPFParticleHits(*gdaughterParticle,evt,fPFParticleTag)).size();
1385  // T0
1386  std::vector<anab::T0> pfdaughterT0vec = pfpUtil.GetPFParticleT0(*gdaughterParticle,evt,fPFParticleTag);
1387  if(!pfdaughterT0vec.empty())
1388  fgranddaughterT0[fNGRANDDAUGHTERS] = pfdaughterT0vec[0].Time();
1389 
1390  if(mcdaughterparticle != 0x0){
1391  fgranddaughter_truth_TrackId[fNGRANDDAUGHTERS] = mcdaughterparticle->TrackId();
1392  fgranddaughter_truth_Pdg[fNGRANDDAUGHTERS] = mcdaughterparticle->PdgCode();
1393  fgranddaughter_truth_Mother[fNGRANDDAUGHTERS] = mcdaughterparticle->Mother();
1394  fgranddaughter_truth_StartPosition[fNGRANDDAUGHTERS][0] = mcdaughterparticle->Vx();
1395  fgranddaughter_truth_StartPosition[fNGRANDDAUGHTERS][1] = mcdaughterparticle->Vy();
1396  fgranddaughter_truth_StartPosition[fNGRANDDAUGHTERS][2] = mcdaughterparticle->Vz();
1397  fgranddaughter_truth_StartPosition[fNGRANDDAUGHTERS][3] = mcdaughterparticle->T();
1398  fgranddaughter_truth_EndPosition[fNGRANDDAUGHTERS][0] = mcdaughterparticle->EndX();
1399  fgranddaughter_truth_EndPosition[fNGRANDDAUGHTERS][1] = mcdaughterparticle->EndY();
1400  fgranddaughter_truth_EndPosition[fNGRANDDAUGHTERS][2] = mcdaughterparticle->EndZ();
1401  fgranddaughter_truth_EndPosition[fNGRANDDAUGHTERS][3] = mcdaughterparticle->EndT();
1402  fgranddaughter_truth_P[fNGRANDDAUGHTERS] = mcdaughterparticle->P();
1403  fgranddaughter_truth_Momentum[fNGRANDDAUGHTERS][0] = mcdaughterparticle->Px();
1404  fgranddaughter_truth_Momentum[fNGRANDDAUGHTERS][1] = mcdaughterparticle->Py();
1405  fgranddaughter_truth_Momentum[fNGRANDDAUGHTERS][2] = mcdaughterparticle->Pz();
1406  fgranddaughter_truth_Momentum[fNGRANDDAUGHTERS][3] = mcdaughterparticle->E();
1407  fgranddaughter_truth_Pt[fNGRANDDAUGHTERS] = mcdaughterparticle->Pt();
1408  fgranddaughter_truth_Mass[fNGRANDDAUGHTERS] = mcdaughterparticle->Mass();
1409  fgranddaughter_truth_EndMomentum[fNGRANDDAUGHTERS][0] = mcdaughterparticle->EndPx();
1410  fgranddaughter_truth_EndMomentum[fNGRANDDAUGHTERS][1] = mcdaughterparticle->EndPy();
1411  fgranddaughter_truth_EndMomentum[fNGRANDDAUGHTERS][2] = mcdaughterparticle->EndPz();
1412  fgranddaughter_truth_EndMomentum[fNGRANDDAUGHTERS][3] = mcdaughterparticle->EndE();
1413  fgranddaughter_truth_Theta[fNGRANDDAUGHTERS] = mcdaughterparticle->Momentum().Theta();
1414  fgranddaughter_truth_Phi[fNGRANDDAUGHTERS] = mcdaughterparticle->Momentum().Phi();
1415  fgranddaughter_truth_TotalLength[fNGRANDDAUGHTERS] = mcdaughterparticle->Trajectory().TotalLength();
1416  fgranddaughter_truth_Process[fNGRANDDAUGHTERS] = truthUtil.GetProcessKey(mcdaughterparticle->Process());
1417  fgranddaughter_truth_EndProcess[fNGRANDDAUGHTERS] = truthUtil.GetProcessKey(mcdaughterparticle->EndProcess());
1418  }
1419  else{
1420  if(fVerbose > 0){
1421  std::cout << "INFO::No MCParticle for granddaughter found!" << std::endl;
1422  }
1423  }
1424 
1425  // Increment counter
1426  fNGRANDDAUGHTERS++;
1427 
1428 }
1429 
1430 // -----------------------------------------------------------------------------
1431 void protoana::ProtoDUNEAnalysisTree::FillConfigTree(){
1432 
1433  fNCryostats = fGeometry->Ncryostats();
1434  fNTPCs = fGeometry->NTPC();
1435  fNChannels = fGeometry->Nchannels();
1436  fNPlanes = fGeometry->Nplanes();
1437  fNAPAs = fGeometry->NTPC()*fGeometry->Ncryostats()/2;
1438  if(fNAPAs > 0)
1439  fNChansPerAPA = fGeometry->Nchannels()/fNAPAs;
1440 
1441  fActiveTPCBoundsX[0] = fActiveTPCBoundsY[0] = fActiveTPCBoundsZ[0] = 1000000.0;
1442  fActiveTPCBoundsX[1] = fActiveTPCBoundsY[1] = fActiveTPCBoundsZ[1] = -1000000.0;
1443 
1444  if(fVerbose > 0){
1445  std::cout << "Detector Name: " << fGeometry->DetectorName() << std::endl;
1446  std::cout << "GDML file: " << fGeometry->GDMLFile() << std::endl;
1447  }
1448 
1449  for (geo::TPCGeo const& TPC: fGeometry->IterateTPCs()) {
1450  // get center in world coordinates
1451  double origin[3] = {0.};
1452  double center[3] = {0.};
1453  TPC.LocalToWorld(origin, center);
1454  double tpc[3] = {TPC.HalfWidth(), TPC.HalfHeight(), 0.5*TPC.Length() };
1455  if(center[0] - tpc[0] < fActiveTPCBoundsX[0]) fActiveTPCBoundsX[0] = center[0] - tpc[0];
1456  if(center[0] + tpc[0] > fActiveTPCBoundsX[1]) fActiveTPCBoundsX[1] = center[0] + tpc[0];
1457  if(center[1] - tpc[1] < fActiveTPCBoundsY[0]) fActiveTPCBoundsY[0] = center[1] - tpc[1];
1458  if(center[1] + tpc[1] > fActiveTPCBoundsY[1]) fActiveTPCBoundsY[1] = center[1] + tpc[1];
1459  if(center[2] - tpc[2] < fActiveTPCBoundsZ[0]) fActiveTPCBoundsZ[0] = center[2] - tpc[2];
1460  if(center[2] + tpc[2] > fActiveTPCBoundsZ[1]) fActiveTPCBoundsZ[1] = center[2] + tpc[2];
1461  }
1462 
1463  fConfigTree->Fill();
1464 
1465 }
1466 
1467 // -----------------------------------------------------------------------------
1468 void protoana::ProtoDUNEAnalysisTree::beginJob(){
1469 
1470  art::ServiceHandle<art::TFileService> tfs;
1471 
1472  // Configuration Tree
1473  fConfigTree = tfs->make<TTree>("Config", "Configuration Tree");
1474  fConfigTree->Branch("NAPAs", &fNAPAs, "NAPAs/I");
1475  fConfigTree->Branch("NChansPerAPA", &fNChansPerAPA, "NChansPerAPA/I");
1476  fConfigTree->Branch("NCryostats", &fNCryostats, "NCryostats/I");
1477  fConfigTree->Branch("NTPCs", &fNTPCs, "NTPCs/I");
1478  fConfigTree->Branch("NChannels", &fNChannels, "NChannels/I");
1479  fConfigTree->Branch("NPlanes", &fNPlanes, "NPlanes/I");
1480  fConfigTree->Branch("ActiveTPCBoundsX", &fActiveTPCBoundsX, "ActiveTPCBoundsX[2]/D");
1481  fConfigTree->Branch("ActiveTPCBoundsY", &fActiveTPCBoundsY, "ActiveTPCBoundsY[2]/D");
1482  fConfigTree->Branch("ActiveTPCBoundsZ", &fActiveTPCBoundsZ, "ActiveTPCBoundsZ[2]/D");
1483 
1484  // Fill configuration tree
1485  FillConfigTree();
1486 
1487  fBeamTrack = tfs->make<TTree>("BeamTrack", "Beam track tree");
1488  fBeamTrack->Branch("run", &fRun, "run/I");
1489  fBeamTrack->Branch("subrun", &fSubRun, "subrun/I");
1490  fBeamTrack->Branch("event", &fevent, "event/I");
1491  fBeamTrack->Branch("beamNTracks", &fbeamNTracks, "beamNTracks/I");
1492  fBeamTrack->Branch("beamtrigger", &fbeamtrigger, "beamtrigger/I");
1493  fBeamTrack->Branch("tof", &ftof, "tof/D");
1494  fBeamTrack->Branch("tof_expElec", &ftof_expElec, "tof_expElec/D");
1495  fBeamTrack->Branch("tof_expMuon", &ftof_expMuon, "tof_expMuon/D");
1496  fBeamTrack->Branch("tof_expPion", &ftof_expPion, "tof_expPion/D");
1497  fBeamTrack->Branch("tof_expKaon", &ftof_expKaon, "tof_expKaon/D");
1498  fBeamTrack->Branch("tof_expProt", &ftof_expProt, "tof_expProt/D");
1499  fBeamTrack->Branch("tof_expDeut", &ftof_expDeut, "tof_expDeut/D");
1500  fBeamTrack->Branch("cerenkovStatus", &fcerenkovStatus, "cerenkovStatus[2]/I");
1501  fBeamTrack->Branch("cerenkovTime", &fcerenkovTime, "cerenkovTime[2]/D");
1502  fBeamTrack->Branch("cerenkovPressure", &fcerenkovPressure, "cerenkovPressure[2]/D");
1503  fBeamTrack->Branch("beamtrackMomentum", &fbeamtrackMomentum, "beamtrackMomentum/D");
1504  fBeamTrack->Branch("beamtrackEnergy", &fbeamtrackEnergy, "beamtrackEnergy/D");
1505  fBeamTrack->Branch("beamtrackPos", &fbeamtrackPos, "beamtrackPos[4]/D");
1506  fBeamTrack->Branch("beamtrackDir", &fbeamtrackDir, "beamtrackDir[3]/D");
1507  fBeamTrack->Branch("BIAndTimingMatched", &fBIAndTimingMatched, "BIAndTimingMatched/I");
1508 
1509  fBeamTrack->Branch("beamtrack_truth_Origin", &fbeamtrack_truth_Origin, "beamtrack_truth_Origin/I");
1510  fBeamTrack->Branch("beamtrackPdg", &fbeamtrackPdg, "beamtrackPdg/I");
1511  fBeamTrack->Branch("beamtrackID", &fbeamtrackID, "beamtrackID/I");
1512  fBeamTrack->Branch("beamtrack_truth_EndPos", &fbeamtrack_truth_EndPos, "beamtrack_truth_EndPos[4]/D");
1513  fBeamTrack->Branch("beamtrack_truth_Momentum", &fbeamtrack_truth_Momentum, "beamtrack_truth_Momentum[4]/D");
1514  fBeamTrack->Branch("beamtrack_truth_KinEnergy", &fbeamtrack_truth_KinEnergy, "beamtrack_truth_KinEnergy/D");
1515  fBeamTrack->Branch("beamtrack_truth_Pt", &fbeamtrack_truth_Pt, "beamtrack_truth_Pt/D");
1516  fBeamTrack->Branch("beamtrack_truth_Mass", &fbeamtrack_truth_Mass, "beamtrack_truth_Mass/D");
1517  fBeamTrack->Branch("beamtrack_truth_Theta", &fbeamtrack_truth_Theta, "beamtrack_truth_Theta/D");
1518  fBeamTrack->Branch("beamtrack_truth_Phi", &fbeamtrack_truth_Phi, "beamtrack_truth_Phi/D");
1519  fBeamTrack->Branch("beamtrack_truth_TotalLength", &fbeamtrack_truth_TotalLength, "beamtrack_truth_TotalLength/D");
1520  fBeamTrack->Branch("beamtrack_truth_Process", &fbeamtrack_truth_Process, "beamtrack_truth_Process/I");
1521  fBeamTrack->Branch("beamtrack_truth_EndProcess", &fbeamtrack_truth_EndProcess, "beamtrack_truth_EndProcess/I");
1522  fBeamTrack->Branch("beamtrack_truth_Pdg", &fbeamtrack_truth_Pdg, "beamtrack_truth_Pdg/I");
1523  fBeamTrack->Branch("beamtrack_truth_NDAUGTHERS", &fbeamtrack_truth_NDAUGTHERS, "beamtrack_truth_NDAUGTHERS/I");
1524  fBeamTrack->Branch("beamtrack_truth_NDECAYDAUGTHERS", &fbeamtrack_truth_NDECAYDAUGTHERS, "beamtrack_truth_NDECAYDAUGTHERS/I");
1525  fBeamTrack->Branch("beamtrack_truth_Pos_InTPCActive", &fbeamtrack_truth_Pos_InTPCActive, "beamtrack_truth_Pos_InTPCActive[4]/D");
1526  fBeamTrack->Branch("beamtrack_truth_Momentum_InTPCActive", &fbeamtrack_truth_Momentum_InTPCActive, "beamtrack_truth_Momentum_InTPCActive[4]/D");
1527  fBeamTrack->Branch("beamtrack_truth_P_InTPCActive", &fbeamtrack_truth_P_InTPCActive, "beamtrack_truth_P_InTPCActive/D");
1528  fBeamTrack->Branch("beamtrack_truth_Pt_InTPCActive", &fbeamtrack_truth_Pt_InTPCActive, "beamtrack_truth_Pt_InTPCActive/D");
1529  fBeamTrack->Branch("beamtrack_truth_Theta_InTPCActive", &fbeamtrack_truth_Theta_InTPCActive, "beamtrack_truth_Theta_InTPCActive/D");
1530  fBeamTrack->Branch("beamtrack_truth_Phi_InTPCActive", &fbeamtrack_truth_Phi_InTPCActive, "beamtrack_truth_Phi_InTPCActive/D");
1531  fBeamTrack->Branch("beamtrack_truth_TotalLength_InTPCActive", &fbeamtrack_truth_TotalLength_InTPCActive, "beamtrack_truth_TotalLength_InTPCActive/D");
1532  fBeamTrack->Branch("beamtrack_truth_KinEnergy_InTPCActive", &fbeamtrack_truth_KinEnergy_InTPCActive, "beamtrack_truth_KinEnergy_InTPCActive/D");
1533 
1534  fBeamTrack->Branch("beamtrack_truthdaughter_TrackId", &fbeamtrack_truthdaughter_TrackId, "beamtrack_truthdaughter_TrackId[beamtrack_truth_NDAUGTHERS]/I");
1535  fBeamTrack->Branch("beamtrack_truthdaughter_Pdg", &fbeamtrack_truthdaughter_Pdg, "beamtrack_truthdaughter_Pdg[beamtrack_truth_NDAUGTHERS]/I");
1536  fBeamTrack->Branch("beamtrack_truthdaughter_Mother", &fbeamtrack_truthdaughter_Mother, "beamtrack_truthdaughter_Mother[beamtrack_truth_NDAUGTHERS]/I");
1537  fBeamTrack->Branch("beamtrack_truthdaughter_StartPosition", &fbeamtrack_truthdaughter_StartPosition, "beamtrack_truthdaughter_StartPosition[beamtrack_truth_NDAUGTHERS][4]/D");
1538  fBeamTrack->Branch("beamtrack_truthdaughter_EndPosition", &fbeamtrack_truthdaughter_EndPosition, "beamtrack_truthdaughter_EndPosition[beamtrack_truth_NDAUGTHERS][4]/D");
1539  fBeamTrack->Branch("beamtrack_truthdaughter_P", &fbeamtrack_truthdaughter_P, "beamtrack_truthdaughter_P[beamtrack_truth_NDAUGTHERS]/D");
1540  fBeamTrack->Branch("beamtrack_truthdaughter_Momentum", &fbeamtrack_truthdaughter_Momentum, "beamtrack_truthdaughter_Momentum[beamtrack_truth_NDAUGTHERS][4]/D");
1541  fBeamTrack->Branch("beamtrack_truthdaughter_EndMomentum", &fbeamtrack_truthdaughter_EndMomentum, "beamtrack_truthdaughter_EndMomentum[beamtrack_truth_NDAUGTHERS][4]/D");
1542  fBeamTrack->Branch("beamtrack_truthdaughter_Pt", &fbeamtrack_truthdaughter_Pt, "beamtrack_truthdaughter_Pt[beamtrack_truth_NDAUGTHERS]/D");
1543  fBeamTrack->Branch("beamtrack_truthdaughter_Mass", &fbeamtrack_truthdaughter_Mass, "beamtrack_truthdaughter_Mass[beamtrack_truth_NDAUGTHERS]/D");
1544  fBeamTrack->Branch("beamtrack_truthdaughter_Theta", &fbeamtrack_truthdaughter_Theta, "beamtrack_truthdaughter_Theta[beamtrack_truth_NDAUGTHERS]/D");
1545  fBeamTrack->Branch("beamtrack_truthdaughter_Phi", &fbeamtrack_truthdaughter_Phi, "beamtrack_truthdaughter_Phi[beamtrack_truth_NDAUGTHERS]/D");
1546  fBeamTrack->Branch("beamtrack_truthdaughter_TotalLength", &fbeamtrack_truthdaughter_TotalLength, "beamtrack_truthdaughter_TotalLength[beamtrack_truth_NDAUGTHERS]/D");
1547  fBeamTrack->Branch("beamtrack_truthdaughter_Process", &fbeamtrack_truthdaughter_Process, "beamtrack_truthdaughter_Process[beamtrack_truth_NDAUGTHERS]/I");
1548  fBeamTrack->Branch("beamtrack_truthdaughter_EndProcess", &fbeamtrack_truthdaughter_EndProcess, "beamtrack_truthdaughter_EndProcess[beamtrack_truth_NDAUGTHERS]/I");
1549 
1550  fBeamTrack->Branch("beamtrack_truthdecaydaughter_TrackId", &fbeamtrack_truthdecaydaughter_TrackId, "beamtrack_truthdecaydaughter_TrackId[beamtrack_truth_NDECAYDAUGTHERS]/I");
1551  fBeamTrack->Branch("beamtrack_truthdecaydaughter_Pdg", &fbeamtrack_truthdecaydaughter_Pdg, "beamtrack_truthdecaydaughter_Pdg[beamtrack_truth_NDECAYDAUGTHERS]/I");
1552  fBeamTrack->Branch("beamtrack_truthdecaydaughter_Mother", &fbeamtrack_truthdecaydaughter_Mother, "beamtrack_truthdecaydaughter_Mother[beamtrack_truth_NDECAYDAUGTHERS]/I");
1553  fBeamTrack->Branch("beamtrack_truthdecaydaughter_StartPosition", &fbeamtrack_truthdecaydaughter_StartPosition, "beamtrack_truthdecaydaughter_StartPosition[beamtrack_truth_NDECAYDAUGTHERS][4]/D");
1554  fBeamTrack->Branch("beamtrack_truthdecaydaughter_EndPosition", &fbeamtrack_truthdecaydaughter_EndPosition, "beamtrack_truthdecaydaughter_EndPosition[beamtrack_truth_NDECAYDAUGTHERS][4]/D");
1555  fBeamTrack->Branch("beamtrack_truthdecaydaughter_P", &fbeamtrack_truthdecaydaughter_P, "beamtrack_truthdecaydaughter_P[beamtrack_truth_NDECAYDAUGTHERS]/D");
1556  fBeamTrack->Branch("beamtrack_truthdecaydaughter_Momentum", &fbeamtrack_truthdecaydaughter_Momentum, "beamtrack_truthdecaydaughter_Momentum[beamtrack_truth_NDECAYDAUGTHERS][4]/D");
1557  fBeamTrack->Branch("beamtrack_truthdecaydaughter_EndMomentum", &fbeamtrack_truthdecaydaughter_EndMomentum, "beamtrack_truthdecaydaughter_EndMomentum[beamtrack_truth_NDECAYDAUGTHERS][4]/D");
1558  fBeamTrack->Branch("beamtrack_truthdecaydaughter_Pt", &fbeamtrack_truthdecaydaughter_Pt, "beamtrack_truthdecaydaughter_Pt[beamtrack_truth_NDECAYDAUGTHERS]/D");
1559  fBeamTrack->Branch("beamtrack_truthdecaydaughter_Mass", &fbeamtrack_truthdecaydaughter_Mass, "beamtrack_truthdecaydaughter_Mass[beamtrack_truth_NDECAYDAUGTHERS]/D");
1560  fBeamTrack->Branch("beamtrack_truthdecaydaughter_Theta", &fbeamtrack_truthdecaydaughter_Theta, "beamtrack_truthdecaydaughter_Theta[beamtrack_truth_NDECAYDAUGTHERS]/D");
1561  fBeamTrack->Branch("beamtrack_truthdecaydaughter_Phi", &fbeamtrack_truthdecaydaughter_Phi, "beamtrack_truthdecaydaughter_Phi[beamtrack_truth_NDECAYDAUGTHERS]/D");
1562  fBeamTrack->Branch("beamtrack_truthdecaydaughter_TotalLength", &fbeamtrack_truthdecaydaughter_TotalLength, "beamtrack_truthdecaydaughter_TotalLength[beamtrack_truth_NDECAYDAUGTHERS]/D");
1563  fBeamTrack->Branch("beamtrack_truthdecaydaughter_Process", &fbeamtrack_truthdecaydaughter_Process, "beamtrack_truthdecaydaughter_Process[beamtrack_truth_NDECAYDAUGTHERS]/I");
1564  fBeamTrack->Branch("beamtrack_truthdecaydaughter_EndProcess", &fbeamtrack_truthdecaydaughter_EndProcess, "beamtrack_truthdecaydaughter_EndProcess[beamtrack_truth_NDECAYDAUGTHERS]/I");
1565 
1566  // Analysis tree
1567  fPandoraBeam = tfs->make<TTree>("ReconstructedBeamEvents", "Reconstructed Beam events");
1568  fPandoraBeam->Branch("run", &fRun, "run/I");
1569  fPandoraBeam->Branch("subrun", &fSubRun, "subrun/I");
1570  fPandoraBeam->Branch("event", &fevent, "event/I");
1571  fPandoraBeam->Branch("timestamp", &fTimeStamp, "timestamp/D");
1572  fPandoraBeam->Branch("Nactivefembs", &fNactivefembs, "Nactivefembs[5]/I");
1573  fPandoraBeam->Branch("NpfParticles", &fNpfParticles, "NpfParticles/I");
1574  fPandoraBeam->Branch("primaryEarliestHitPeakTim", &fprimaryEarliestHitPeakTime, "primaryEarliestHitPeakTim/D");
1575  fPandoraBeam->Branch("primaryVertex", &fprimaryVertex, "primaryVertex[3]/D");
1576  fPandoraBeam->Branch("primaryIstrack", &fprimaryIstrack, "primaryIstrack/I");
1577  fPandoraBeam->Branch("primaryIsshower", &fprimaryIsshower, "primaryIsshower/I");
1578  fPandoraBeam->Branch("primaryBDTScore", &fprimaryBDTScore, "primaryBDTScore/D");
1579  fPandoraBeam->Branch("primaryNHits", &fprimaryNHits, "primaryNHits/I");
1580  fPandoraBeam->Branch("primaryTheta", &fprimaryTheta, "primaryTheta/D");
1581  fPandoraBeam->Branch("primaryPhi", &fprimaryPhi, "primaryPhi/D");
1582  fPandoraBeam->Branch("primaryLength", &fprimaryLength, "primaryLength/D");
1583  fPandoraBeam->Branch("primaryMomentum", &fprimaryMomentum, "primaryMomentum/D");
1584  fPandoraBeam->Branch("primaryEndMomentum", &fprimaryEndMomentum, "primaryEndMomentum/D");
1585  fPandoraBeam->Branch("primaryEndPosition", &fprimaryEndPosition, "primaryEndPosition[3]/D");
1586  fPandoraBeam->Branch("primaryStartPosition", &fprimaryStartPosition, "primaryStartPosition[3]/D");
1587  fPandoraBeam->Branch("primaryEndDirection", &fprimaryEndDirection, "primaryEndDirection[3]/D");
1588  fPandoraBeam->Branch("primaryStartDirection", &fprimaryStartDirection, "primaryStartDirection[3]/D");
1589  fPandoraBeam->Branch("primaryOpeningAngle", &fprimaryOpeningAngle, "primaryOpeningAngle/D");
1590  fPandoraBeam->Branch("primaryID", &fprimaryID, "primaryID/I");
1591  fPandoraBeam->Branch("primaryShowerBestPlane", &fprimaryShowerBestPlane, "primaryShowerBestPlane/I");
1592  fPandoraBeam->Branch("primaryShowerEnergy", &fprimaryShowerEnergy, "primaryShowerEnergy[3]/I");
1593  fPandoraBeam->Branch("primaryShowerMIPEnergy", &fprimaryShowerMIPEnergy, "primaryShowerMIPEnergy[3]/I");
1594  fPandoraBeam->Branch("primaryShowerdEdx", &fprimaryShowerdEdx, "primaryShowerdEdx[3]/I");
1595  fPandoraBeam->Branch("primaryMomentumByRangeProton", &fprimaryMomentumByRangeProton, "primaryMomentumByRangeProton/D");
1596  fPandoraBeam->Branch("primaryMomentumByRangeMuon", &fprimaryMomentumByRangeMuon, "primaryMomentumByRangeMuon/D");
1597  fPandoraBeam->Branch("primaryKineticEnergy", &fprimaryKineticEnergy, "primaryKineticEnergy[3]/D");
1598  fPandoraBeam->Branch("primaryRange", &fprimaryRange, "primaryRange[3]/D");
1599  fPandoraBeam->Branch("primaryTrkPitchC", &fprimaryTrkPitchC, "primaryTrkPitchC[3]/D");
1600  fPandoraBeam->Branch("primaryT0", &fprimaryT0, "primaryT0/D");
1601  fPandoraBeam->Branch("primaryPID_Pdg", &fprimaryPID_Pdg, "primaryPID_Pdg[3]/I");
1602  fPandoraBeam->Branch("primaryPID_Ndf", &fprimaryPID_Ndf, "primaryPID_Ndf[3]/I");
1603  fPandoraBeam->Branch("primaryPID_MinChi2", &fprimaryPID_MinChi2, "primaryPID_MinChi2[3]/D");
1604  fPandoraBeam->Branch("primaryPID_DeltaChi2", &fprimaryPID_DeltaChi2, "primaryPID_DeltaChi2[3]/D");
1605  fPandoraBeam->Branch("primaryPID_Chi2Proton", &fprimaryPID_Chi2Proton, "primaryPID_Chi2Proton[3]/D");
1606  fPandoraBeam->Branch("primaryPID_Chi2Kaon", &fprimaryPID_Chi2Kaon, "primaryPID_Chi2Kaon[3]/D");
1607  fPandoraBeam->Branch("primaryPID_Chi2Pion", &fprimaryPID_Chi2Pion, "primaryPID_Chi2Pion[3]/D");
1608  fPandoraBeam->Branch("primaryPID_Chi2Muon", &fprimaryPID_Chi2Muon, "primaryPID_Chi2Muon[3]/D");
1609  fPandoraBeam->Branch("primaryPID_MissingE", &fprimaryPID_MissingE, "primaryPID_MissingE[3]/D");
1610  fPandoraBeam->Branch("primaryPID_MissingEavg", &fprimaryPID_MissingEavg, "primaryPID_MissingEavg[3]/D");
1611  fPandoraBeam->Branch("primaryPID_PIDA", &fprimaryPID_PIDA, "primaryPID_PIDA[3]/D");
1612 
1613  fPandoraBeam->Branch("primary_truth_Origin", &fprimary_truth_Origin, "primary_truth_Origin/I");
1614  fPandoraBeam->Branch("primary_truth_TrackId", &fprimary_truth_TrackId, "primary_truth_TrackId/I");
1615  fPandoraBeam->Branch("primary_truth_Pdg", &fprimary_truth_Pdg, "primary_truth_Pdg/I");
1616  fPandoraBeam->Branch("primary_truth_Mother", &fprimary_truth_Mother, "primary_truth_Mother/I");
1617  fPandoraBeam->Branch("primary_truth_StartPosition", &fprimary_truth_StartPosition, "primary_truth_StartPosition[4]/D");
1618  fPandoraBeam->Branch("primary_truth_EndPosition", &fprimary_truth_EndPosition, "primary_truth_EndPosition[4]/D");
1619  fPandoraBeam->Branch("primary_truth_Momentum", &fprimary_truth_Momentum, "primary_truth_Momentum[4]/D");
1620  fPandoraBeam->Branch("primary_truth_EndMomentum", &fprimary_truth_EndMomentum, "primary_truth_EndMomentum[4]/D");
1621  fPandoraBeam->Branch("primary_truth_P", &fprimary_truth_P, "primary_truth_P/D");
1622  fPandoraBeam->Branch("primary_truth_Pt", &fprimary_truth_Pt, "primary_truth_Pt/D");
1623  fPandoraBeam->Branch("primary_truth_Mass", &fprimary_truth_Mass, "primary_truth_Mass/D");
1624  fPandoraBeam->Branch("primary_truth_Theta", &fprimary_truth_Theta, "primary_truth_Theta/D");
1625  fPandoraBeam->Branch("primary_truth_Phi", &fprimary_truth_Phi, "primary_truth_Phi/D");
1626  fPandoraBeam->Branch("primary_truth_TotalLength", &fprimary_truth_TotalLength, "primary_truth_TotalLength/D");
1627  fPandoraBeam->Branch("primary_truth_Process", &fprimary_truth_Process, "primary_truth_Process/I");
1628  fPandoraBeam->Branch("primary_truth_EndProcess", &fprimary_truth_EndProcess, "primary_truth_EndProcess/I");
1629  fPandoraBeam->Branch("primary_truth_Isbeammatched", &fprimary_truth_Isbeammatched, "primary_truth_Isbeammatched/I");
1630 
1631  fPandoraBeam->Branch("primary_truth_EkinAtVertex_notcorrected", &fprimary_truth_EkinAtVertex_notcorrected, "primary_truth_EkinAtVertex_notcorrected/D");
1632  fPandoraBeam->Branch("primary_truth_EkinAtVertex", &fprimary_truth_EkinAtVertex, "primary_truth_EkinAtVertex/D");
1633  fPandoraBeam->Branch("primary_truth_Pos_InTPCActive", &fprimary_truth_Pos_InTPCActive, "primary_truth_Pos_InTPCActive[4]/D");
1634  fPandoraBeam->Branch("primary_truth_Momentum_InTPCActive", &fprimary_truth_Momentum_InTPCActive, "primary_truth_Momentum_InTPCActive[4]/D");
1635  fPandoraBeam->Branch("primary_truth_P_InTPCActive", &fprimary_truth_P_InTPCActive, "primary_truth_P_InTPCActive/D");
1636  fPandoraBeam->Branch("primary_truth_Pt_InTPCActive", &fprimary_truth_Pt_InTPCActive, "primary_truth_Pt_InTPCActive/D");
1637  fPandoraBeam->Branch("primary_truth_Theta_InTPCActive", &fprimary_truth_Theta_InTPCActive, "primary_truth_Theta_InTPCActive/D");
1638  fPandoraBeam->Branch("primary_truth_Phi_InTPCActive", &fprimary_truth_Phi_InTPCActive, "primary_truth_Phi_InTPCActive/D");
1639  fPandoraBeam->Branch("primary_truth_TotalLength_InTPCActive", &fprimary_truth_TotalLength_InTPCActive, "primary_truth_TotalLength_InTPCActive/D");
1640  fPandoraBeam->Branch("primary_truth_KinEnergy_InTPCActive", &fprimary_truth_KinEnergy_InTPCActive, "primary_truth_KinEnergy_InTPCActive/D");
1641 
1642  fPandoraBeam->Branch("primary_truth_NDAUGTHERS", &fprimary_truth_NDAUGTHERS, "primary_truth_NDAUGTHERS/I");
1643  fPandoraBeam->Branch("primary_truth_NDECAYDAUGTHERS", &fprimary_truth_NDECAYDAUGTHERS, "primary_truth_NDECAYDAUGTHERS/I");
1644  fPandoraBeam->Branch("primary_truthdaughter_TrackId", &fprimary_truthdaughter_TrackId, "primary_truthdaughter_TrackId[primary_truth_NDAUGTHERS]/I");
1645  fPandoraBeam->Branch("primary_truthdaughter_Pdg", &fprimary_truthdaughter_Pdg, "primary_truthdaughter_Pdg[primary_truth_NDAUGTHERS]/I");
1646  fPandoraBeam->Branch("primary_truthdaughter_Mother", &fprimary_truthdaughter_Mother, "primary_truthdaughter_Mother[primary_truth_NDAUGTHERS]/I");
1647  fPandoraBeam->Branch("primary_truthdaughter_StartPosition", &fprimary_truthdaughter_StartPosition, "primary_truthdaughter_StartPosition[primary_truth_NDAUGTHERS][4]/D");
1648  fPandoraBeam->Branch("primary_truthdaughter_EndPosition", &fprimary_truthdaughter_EndPosition, "primary_truthdaughter_EndPosition[primary_truth_NDAUGTHERS][4]/D");
1649  fPandoraBeam->Branch("primary_truthdaughter_P", &fprimary_truthdaughter_P, "primary_truthdaughter_P[primary_truth_NDAUGTHERS]/D");
1650  fPandoraBeam->Branch("primary_truthdaughter_Momentum", &fprimary_truthdaughter_Momentum, "primary_truthdaughter_Momentum[primary_truth_NDAUGTHERS][4]/D");
1651  fPandoraBeam->Branch("primary_truthdaughter_EndMomentum", &fprimary_truthdaughter_EndMomentum, "primary_truthdaughter_EndMomentum[primary_truth_NDAUGTHERS][4]/D");
1652  fPandoraBeam->Branch("primary_truthdaughter_Pt", &fprimary_truthdaughter_Pt, "primary_truthdaughter_Pt[primary_truth_NDAUGTHERS]/D");
1653  fPandoraBeam->Branch("primary_truthdaughter_Mass", &fprimary_truthdaughter_Mass, "primary_truthdaughter_Mass[primary_truth_NDAUGTHERS]/D");
1654  fPandoraBeam->Branch("primary_truthdaughter_Theta", &fprimary_truthdaughter_Theta, "primary_truthdaughter_Theta[primary_truth_NDAUGTHERS]/D");
1655  fPandoraBeam->Branch("primary_truthdaughter_Phi", &fprimary_truthdaughter_Phi, "primary_truthdaughter_Phi[primary_truth_NDAUGTHERS]/D");
1656  fPandoraBeam->Branch("primary_truthdaughter_TotalLength", &fprimary_truthdaughter_TotalLength, "primary_truthdaughter_TotalLength[primary_truth_NDAUGTHERS]/D");
1657  fPandoraBeam->Branch("primary_truthdaughter_Process", &fprimary_truthdaughter_Process, "primary_truthdaughter_Process[primary_truth_NDAUGTHERS]/I");
1658  fPandoraBeam->Branch("primary_truthdaughter_EndProcess", &fprimary_truthdaughter_EndProcess, "primary_truthdaughter_EndProcess[primary_truth_NDAUGTHERS]/I");
1659 
1660  fPandoraBeam->Branch("primary_truthdecaydaughter_TrackId", &fprimary_truthdecaydaughter_TrackId, "primary_truthdecaydaughter_TrackId[primary_truth_NDAUGTHERS]/I");
1661  fPandoraBeam->Branch("primary_truthdecaydaughter_Pdg", &fprimary_truthdecaydaughter_Pdg, "primary_truthdecaydaughter_Pdg[primary_truth_NDAUGTHERS]/I");
1662  fPandoraBeam->Branch("primary_truthdecaydaughter_Mother", &fprimary_truthdecaydaughter_Mother, "primary_truthdecaydaughter_Mother[primary_truth_NDAUGTHERS]/I");
1663  fPandoraBeam->Branch("primary_truthdecaydaughter_StartPosition", &fprimary_truthdecaydaughter_StartPosition, "primary_truthdecaydaughter_StartPosition[primary_truth_NDAUGTHERS][4]/D");
1664  fPandoraBeam->Branch("primary_truthdecaydaughter_EndPosition", &fprimary_truthdecaydaughter_EndPosition, "primary_truthdecaydaughter_EndPosition[primary_truth_NDAUGTHERS][4]/D");
1665  fPandoraBeam->Branch("primary_truthdecaydaughter_P", &fprimary_truthdecaydaughter_P, "primary_truthdecaydaughter_P[primary_truth_NDAUGTHERS]/D");
1666  fPandoraBeam->Branch("primary_truthdecaydaughter_Momentum", &fprimary_truthdecaydaughter_Momentum, "primary_truthdecaydaughter_Momentum[primary_truth_NDAUGTHERS][4]/D");
1667  fPandoraBeam->Branch("primary_truthdecaydaughter_EndMomentum", &fprimary_truthdecaydaughter_EndMomentum, "primary_truthdecaydaughter_EndMomentum[primary_truth_NDAUGTHERS][4]/D");
1668  fPandoraBeam->Branch("primary_truthdecaydaughter_Pt", &fprimary_truthdecaydaughter_Pt, "primary_truthdecaydaughter_Pt[primary_truth_NDAUGTHERS]/D");
1669  fPandoraBeam->Branch("primary_truthdecaydaughter_Mass", &fprimary_truthdecaydaughter_Mass, "primary_truthdecaydaughter_Mass[primary_truth_NDAUGTHERS]/D");
1670  fPandoraBeam->Branch("primary_truthdecaydaughter_Theta", &fprimary_truthdecaydaughter_Theta, "primary_truthdecaydaughter_Theta[primary_truth_NDAUGTHERS]/D");
1671  fPandoraBeam->Branch("primary_truthdecaydaughter_Phi", &fprimary_truthdecaydaughter_Phi, "primary_truthdecaydaughter_Phi[primary_truth_NDAUGTHERS]/D");
1672  fPandoraBeam->Branch("primary_truthdecaydaughter_TotalLength", &fprimary_truthdecaydaughter_TotalLength, "primary_truthdecaydaughter_TotalLength[primary_truth_NDAUGTHERS]/D");
1673  fPandoraBeam->Branch("primary_truthdecaydaughter_Process", &fprimary_truthdecaydaughter_Process, "primary_truthdecaydaughter_Process[primary_truth_NDAUGTHERS]/I");
1674  fPandoraBeam->Branch("primary_truthdecaydaughter_EndProcess", &fprimary_truthdecaydaughter_EndProcess, "primary_truthdecaydaughter_EndProcess[primary_truth_NDAUGTHERS]/I");
1675 
1676  fDaughterTree = tfs->make<TTree>("ReconstructedDaughters", "Reconstructed Daughters");
1677  fDaughterTree->Branch("NDAUGHTERS", &fNDAUGHTERS, "NDAUGHTERS/I");
1678  fDaughterTree->Branch("daughterVertex", &fdaughterVertex, "daughterVertex[3]/D");
1679  fDaughterTree->Branch("daughterIstrack", &fdaughterIstrack, "daughterIstrack[NDAUGHTERS]/I");
1680  fDaughterTree->Branch("daughterIsshower", &fdaughterIsshower, "daughterIsshower[NDAUGHTERS]/I");
1681  fDaughterTree->Branch("daughterNHits", &fdaughterNHits, "daughterNHits[NDAUGHTERS]/I");
1682  fDaughterTree->Branch("daughterTheta", &fdaughterTheta, "daughterTheta[NDAUGHTERS]/D");
1683  fDaughterTree->Branch("daughterPhi", &fdaughterPhi, "daughterPhi[NDAUGHTERS]/D");
1684  fDaughterTree->Branch("daughterLength", &fdaughterLength, "daughterLength[NDAUGHTERS]/D");
1685  fDaughterTree->Branch("daughterMomentum", &fdaughterMomentum, "daughterMomentum[NDAUGHTERS]/D");
1686  fDaughterTree->Branch("daughterEndMomentum", &fdaughterEndMomentum, "daughterEndMomentum[NDAUGHTERS]/D");
1687  fDaughterTree->Branch("daughterEndPosition", &fdaughterEndPosition, "daughterEndPosition[NDAUGHTERS][3]/D");
1688  fDaughterTree->Branch("daughterStartPosition", &fdaughterStartPosition, "daughterStartPosition[NDAUGHTERS][3]/D");
1689  fDaughterTree->Branch("daughterStartDirection", &fdaughterStartDirection, "daughterStartDirection[NDAUGHTERS][3]/D");
1690  fDaughterTree->Branch("daughterEndDirection", &fdaughterEndDirection, "daughterEndDirection[NDAUGHTERS][3]/D");
1691  fDaughterTree->Branch("daughterOpeningAngle", &fdaughterOpeningAngle, "daughterOpeningAngle[NDAUGHTERS]/D");
1692  fDaughterTree->Branch("daughterShowerBestPlane", &fdaughterShowerBestPlane, "daughterShowerBestPlane[NDAUGHTERS]/I");
1693  fDaughterTree->Branch("daughterShowerEnergy", &fdaughterShowerEnergy, "daughterShowerEnergy[NDAUGHTERS][3]/D");
1694  fDaughterTree->Branch("daughterShowerMIPEnergy", &fdaughterShowerMIPEnergy, "daughterShowerMIPEnergy[NDAUGHTERS][3]/D");
1695  fDaughterTree->Branch("daughterShowerdEdx", &fdaughterShowerdEdx, "daughterShowerdEdx[NDAUGHTERS][3]/D");
1696  fDaughterTree->Branch("daughterMomentumByRangeProton", &fdaughterMomentumByRangeProton, "daughterMomentumByRangeProton[NDAUGHTERS]/D");
1697  fDaughterTree->Branch("daughterMomentumByRangeMuon", &fdaughterMomentumByRangeMuon, "daughterMomentumByRangeMuon[NDAUGHTERS]/D");
1698  fDaughterTree->Branch("daughterKineticEnergy", &fdaughterKineticEnergy, "daughterKineticEnergy[NDAUGHTERS][3]/D");
1699  fDaughterTree->Branch("daughterRange", &fdaughterRange, "daughterRange[NDAUGHTERS][3]/D");
1700  fDaughterTree->Branch("daughterTrkPitchC", &fdaughterTrkPitchC, "daughterTrkPitchC[NDAUGHTERS][3]/D");
1701  fDaughterTree->Branch("daughterID", &fdaughterID, "daughterID[NDAUGHTERS]/I");
1702  fDaughterTree->Branch("daughterT0", &fdaughterT0, "daughterT0[NDAUGHTERS]/D");
1703 
1704  fDaughterTree->Branch("daughterPID_Pdg", &fdaughterPID_Pdg, "daughterPID_Pdg[NDAUGHTERS][3]/I");
1705  fDaughterTree->Branch("daughterPID_Ndf", &fdaughterPID_Ndf, "daughterPID_Ndf[NDAUGHTERS][3]/I");
1706  fDaughterTree->Branch("daughterPID_MinChi2", &fdaughterPID_MinChi2, "daughterPID_MinChi2[NDAUGHTERS][3]/D");
1707  fDaughterTree->Branch("daughterPID_DeltaChi2", &fdaughterPID_DeltaChi2, "daughterPID_DeltaChi2[NDAUGHTERS][3]/D");
1708  fDaughterTree->Branch("daughterPID_Chi2Proton", &fdaughterPID_Chi2Proton, "daughterPID_Chi2Proton[NDAUGHTERS][3]/D");
1709  fDaughterTree->Branch("daughterPID_Chi2Kaon", &fdaughterPID_Chi2Kaon, "daughterPID_Chi2Kaon[NDAUGHTERS][3]/D");
1710  fDaughterTree->Branch("daughterPID_Chi2Pion", &fdaughterPID_Chi2Pion, "daughterPID_Chi2Pion[NDAUGHTERS][3]/D");
1711  fDaughterTree->Branch("daughterPID_Chi2Muon", &fdaughterPID_Chi2Muon, "daughterPID_Chi2Muon[NDAUGHTERS][3]/D");
1712  fDaughterTree->Branch("daughterPID_MissingE", &fdaughterPID_MissingE, "daughterPID_MissingE[NDAUGHTERS][3]/D");
1713  fDaughterTree->Branch("daughterPID_MissingEavg", &fdaughterPID_MissingEavg, "daughterPID_MissingEavg[NDAUGHTERS][3]/D");
1714  fDaughterTree->Branch("daughterPID_PIDA", &fdaughterPID_PIDA, "daughterPID_PIDA[NDAUGHTERS][3]/D");
1715 
1716  fDaughterTree->Branch("daughter_truth_TrackId", &fdaughter_truth_TrackId, "daughter_truth_TrackId[NDAUGHTERS]/I");
1717  fDaughterTree->Branch("daughter_truth_Pdg", &fdaughter_truth_Pdg, "daughter_truth_Pdg[NDAUGHTERS]/I");
1718  fDaughterTree->Branch("daughter_truth_Mother", &fdaughter_truth_Mother, "daughter_truth_Mother[NDAUGHTERS]/I");
1719  fDaughterTree->Branch("daughter_truth_StartPosition", &fdaughter_truth_StartPosition, "daughter_truth_StartPosition[NDAUGHTERS][4]/D");
1720  fDaughterTree->Branch("daughter_truth_EndPosition", &fdaughter_truth_EndPosition, "daughter_truth_EndPosition[NDAUGHTERS][4]/D");
1721  fDaughterTree->Branch("daughter_truth_Momentum", &fdaughter_truth_Momentum, "daughter_truth_Momentum[NDAUGHTERS][4]/D");
1722  fDaughterTree->Branch("daughter_truth_EndMomentum", &fdaughter_truth_EndMomentum, "daughter_truth_EndMomentum[NDAUGHTERS][4]/D");
1723  fDaughterTree->Branch("daughter_truth_P", &fdaughter_truth_P, "daughter_truth_P[NDAUGHTERS]/D");
1724  fDaughterTree->Branch("daughter_truth_Pt", &fdaughter_truth_Pt, "daughter_truth_Pt[NDAUGHTERS]/D");
1725  fDaughterTree->Branch("daughter_truth_Mass", &fdaughter_truth_Mass, "daughter_truth_Mass[NDAUGHTERS]/D");
1726  fDaughterTree->Branch("daughter_truth_Theta", &fdaughter_truth_Theta, "daughter_truth_Theta[NDAUGHTERS]/D");
1727  fDaughterTree->Branch("daughter_truth_Phi", &fdaughter_truth_Phi, "daughter_truth_Phi[NDAUGHTERS]/D");
1728  fDaughterTree->Branch("daughter_truth_TotalLength", &fdaughter_truth_TotalLength, "daughter_truth_TotalLength[NDAUGHTERS]/D");
1729  fDaughterTree->Branch("daughter_truth_Process", &fdaughter_truth_Process, "daughter_truth_Process[NDAUGHTERS]/I");
1730  fDaughterTree->Branch("daughter_truth_EndProcess", &fdaughter_truth_EndProcess, "daughter_truth_EndProcess[NDAUGHTERS]/I");
1731 
1732  fGrandDaughterTree = tfs->make<TTree>("ReconstructedGrandDaughters", "Reconstructed GrandDaughters");
1733  fGrandDaughterTree->Branch("NDAUGHTERS", &fNDAUGHTERS, "NDAUGHTERS/I");
1734  fGrandDaughterTree->Branch("NGRANDDAUGHTERS", &fNGRANDDAUGHTERS, "NGRANDDAUGHTERS/I");
1735  fGrandDaughterTree->Branch("granddaughterVertex", &fgranddaughterVertex, "granddaughterVertex[NDAUGHTERS][3]/D");
1736  fGrandDaughterTree->Branch("granddaughterIstrack", &fgranddaughterIstrack, "granddaughterIstrack[NGRANDDAUGHTERS]/I");
1737  fGrandDaughterTree->Branch("granddaughterIsshower", &fgranddaughterIsshower, "granddaughterIsshower[NGRANDDAUGHTERS]/I");
1738  fGrandDaughterTree->Branch("granddaughterNHits", &fgranddaughterNHits, "granddaughterNHits[NGRANDDAUGHTERS]/I");
1739  fGrandDaughterTree->Branch("granddaughterTheta", &fgranddaughterTheta, "granddaughterTheta[NGRANDDAUGHTERS]/D");
1740  fGrandDaughterTree->Branch("granddaughterPhi", &fgranddaughterPhi, "granddaughterPhi[NGRANDDAUGHTERS]/D");
1741  fGrandDaughterTree->Branch("granddaughterLength", &fgranddaughterLength, "granddaughterLength[NGRANDDAUGHTERS]/D");
1742  fGrandDaughterTree->Branch("granddaughterMomentum", &fgranddaughterMomentum, "granddaughterMomentum[NGRANDDAUGHTERS]/D");
1743  fGrandDaughterTree->Branch("granddaughterEndMomentum", &fgranddaughterEndMomentum, "granddaughterEndMomentum[NGRANDDAUGHTERS]/D");
1744  fGrandDaughterTree->Branch("granddaughterEndPosition", &fgranddaughterEndPosition, "granddaughterEndPosition[NGRANDDAUGHTERS][3]/D");
1745  fGrandDaughterTree->Branch("granddaughterStartPosition", &fgranddaughterStartPosition, "granddaughterStartPosition[NGRANDDAUGHTERS][3]/D");
1746  fGrandDaughterTree->Branch("granddaughterStartDirection", &fgranddaughterStartDirection, "granddaughterStartDirection[NGRANDDAUGHTERS][3]/D");
1747  fGrandDaughterTree->Branch("granddaughterEndDirection", &fgranddaughterEndDirection, "granddaughterEndDirection[NGRANDDAUGHTERS][3]/D");
1748  fGrandDaughterTree->Branch("granddaughterOpeningAngle", &fgranddaughterOpeningAngle, "granddaughterOpeningAngle[NGRANDDAUGHTERS]/D");
1749  fGrandDaughterTree->Branch("granddaughterShowerBestPlane", &fgranddaughterShowerBestPlane, "granddaughterShowerBestPlane[NGRANDDAUGHTERS]/I");
1750  fGrandDaughterTree->Branch("granddaughterShowerEnergy", &fgranddaughterShowerEnergy, "granddaughterShowerEnergy[NGRANDDAUGHTERS][3]/D");
1751  fGrandDaughterTree->Branch("granddaughterShowerMIPEnergy", &fgranddaughterShowerMIPEnergy, "granddaughterShowerMIPEnergy[NGRANDDAUGHTERS][3]/D");
1752  fGrandDaughterTree->Branch("granddaughterShowerdEdx", &fgranddaughterShowerdEdx, "granddaughterShowerdEdx[NGRANDDAUGHTERS][3]/D");
1753  fGrandDaughterTree->Branch("granddaughterMomentumByRangeProton", &fgranddaughterMomentumByRangeProton, "granddaughterMomentumByRangeProton[NGRANDDAUGHTERS]/D");
1754  fGrandDaughterTree->Branch("granddaughterMomentumByRangeMuon", &fgranddaughterMomentumByRangeMuon, "granddaughterMomentumByRangeMuon[NGRANDDAUGHTERS]/D");
1755  fGrandDaughterTree->Branch("granddaughterKineticEnergy", &fgranddaughterKineticEnergy, "granddaughterKineticEnergy[NGRANDDAUGHTERS][3]/D");
1756  fGrandDaughterTree->Branch("granddaughterRange", &fgranddaughterRange, "granddaughterRange[NGRANDDAUGHTERS][3]/D");
1757  fGrandDaughterTree->Branch("granddaughterTrkPitchC", &fgranddaughterTrkPitchC, "granddaughterTrkPitchC[NGRANDDAUGHTERS][3]/D");
1758  fGrandDaughterTree->Branch("granddaughterID", &fgranddaughterID, "granddaughterID[NGRANDDAUGHTERS]/I");
1759  fGrandDaughterTree->Branch("granddaughterMotherID", &fgranddaughterMotherID, "granddaughterMotherID[NGRANDDAUGHTERS]/I");
1760  fGrandDaughterTree->Branch("granddaughterT0", &fgranddaughterT0, "granddaughterT0[NGRANDDAUGHTERS]/D");
1761 
1762  fGrandDaughterTree->Branch("granddaughterPID_Pdg", &fgranddaughterPID_Pdg, "granddaughterPID_Pdg[NGRANDDAUGHTERS][3]/I");
1763  fGrandDaughterTree->Branch("granddaughterPID_Ndf", &fgranddaughterPID_Ndf, "granddaughterPID_Ndf[NGRANDDAUGHTERS][3]/I");
1764  fGrandDaughterTree->Branch("granddaughterPID_MinChi2", &fgranddaughterPID_MinChi2, "granddaughterPID_MinChi2[NGRANDDAUGHTERS][3]/D");
1765  fGrandDaughterTree->Branch("granddaughterPID_DeltaChi2", &fgranddaughterPID_DeltaChi2, "granddaughterPID_DeltaChi2[NGRANDDAUGHTERS][3]/D");
1766  fGrandDaughterTree->Branch("granddaughterPID_Chi2Proton", &fgranddaughterPID_Chi2Proton, "granddaughterPID_Chi2Proton[NGRANDDAUGHTERS][3]/D");
1767  fGrandDaughterTree->Branch("granddaughterPID_Chi2Kaon", &fgranddaughterPID_Chi2Kaon, "granddaughterPID_Chi2Kaon[NGRANDDAUGHTERS][3]/D");
1768  fGrandDaughterTree->Branch("granddaughterPID_Chi2Pion", &fgranddaughterPID_Chi2Pion, "granddaughterPID_Chi2Pion[NGRANDDAUGHTERS][3]/D");
1769  fGrandDaughterTree->Branch("granddaughterPID_Chi2Muon", &fgranddaughterPID_Chi2Muon, "granddaughterPID_Chi2Muon[NGRANDDAUGHTERS][3]/D");
1770  fGrandDaughterTree->Branch("granddaughterPID_MissingE", &fgranddaughterPID_MissingE, "granddaughterPID_MissingE[NGRANDDAUGHTERS][3]/D");
1771  fGrandDaughterTree->Branch("granddaughterPID_MissingEavg", &fgranddaughterPID_MissingEavg, "granddaughterPID_MissingEavg[NGRANDDAUGHTERS][3]/D");
1772  fGrandDaughterTree->Branch("granddaughterPID_PIDA", &fgranddaughterPID_PIDA, "granddaughterPID_PIDA[NGRANDDAUGHTERS][3]/D");
1773 
1774  fGrandDaughterTree->Branch("granddaughter_truth_TrackId", &fgranddaughter_truth_TrackId, "granddaughter_truth_TrackId[NGRANDDAUGHTERS]/I");
1775  fGrandDaughterTree->Branch("granddaughter_truth_Pdg", &fgranddaughter_truth_Pdg, "granddaughter_truth_Pdg[NGRANDDAUGHTERS]/I");
1776  fGrandDaughterTree->Branch("granddaughter_truth_Mother", &fgranddaughter_truth_Mother, "granddaughter_truth_Mother[NGRANDDAUGHTERS]/I");
1777  fGrandDaughterTree->Branch("granddaughter_truth_StartPosition", &fgranddaughter_truth_StartPosition, "granddaughter_truth_StartPosition[NGRANDDAUGHTERS][4]/D");
1778  fGrandDaughterTree->Branch("granddaughter_truth_EndPosition", &fgranddaughter_truth_EndPosition, "granddaughter_truth_EndPosition[NGRANDDAUGHTERS][4]/D");
1779  fGrandDaughterTree->Branch("granddaughter_truth_Momentum", &fgranddaughter_truth_Momentum, "granddaughter_truth_Momentum[NGRANDDAUGHTERS][4]/D");
1780  fGrandDaughterTree->Branch("granddaughter_truth_EndMomentum", &fgranddaughter_truth_EndMomentum, "granddaughter_truth_EndMomentum[NGRANDDAUGHTERS][4]/D");
1781  fGrandDaughterTree->Branch("granddaughter_truth_P", &fgranddaughter_truth_P, "granddaughter_truth_P[NGRANDDAUGHTERS]/D");
1782  fGrandDaughterTree->Branch("granddaughter_truth_Pt", &fgranddaughter_truth_Pt, "granddaughter_truth_Pt[NGRANDDAUGHTERS]/D");
1783  fGrandDaughterTree->Branch("granddaughter_truth_Mass", &fgranddaughter_truth_Mass, "granddaughter_truth_Mass[NGRANDDAUGHTERS]/D");
1784  fGrandDaughterTree->Branch("granddaughter_truth_Theta", &fgranddaughter_truth_Theta, "granddaughter_truth_Theta[NGRANDDAUGHTERS]/D");
1785  fGrandDaughterTree->Branch("granddaughter_truth_Phi", &fgranddaughter_truth_Phi, "granddaughter_truth_Phi[NGRANDDAUGHTERS]/D");
1786  fGrandDaughterTree->Branch("granddaughter_truth_TotalLength", &fgranddaughter_truth_TotalLength, "granddaughter_truth_TotalLength[NGRANDDAUGHTERS]/D");
1787  fGrandDaughterTree->Branch("granddaughter_truth_Process", &fgranddaughter_truth_Process, "granddaughter_truth_Process[NGRANDDAUGHTERS]/I");
1788  fGrandDaughterTree->Branch("granddaughter_truth_EndProcess", &fgranddaughter_truth_EndProcess, "granddaughter_truth_EndProcess[NGRANDDAUGHTERS]/I");
1789 
1790 }
1791 
1792 // -----------------------------------------------------------------------------
1793 void protoana::ProtoDUNEAnalysisTree::Initialise(){
1794 
1795  fRun = -999;
1796  fSubRun = -999;
1797  fevent = -999;
1798  fTimeStamp = -999.0;
1799  fNpfParticles = 0;
1800  for(int k=0; k < 5; k++)
1801  fNactivefembs[k] = -999;
1802 
1803  for(int k=0; k < 3; k++){
1804  fprimaryVertex[k] = -999.0;
1805  fdaughterVertex[k] = -999.0;
1806  fprimaryEndPosition[k] = -999.0;
1807  fprimaryStartPosition[k] = -999.0;
1808  fprimaryEndDirection[k] = -999.0;
1809  fprimaryStartDirection[k] = -999.0;
1810  fprimaryKineticEnergy[k] = -999.0;
1811  fprimaryRange[k] = -999.0;
1812  fprimaryTrkPitchC[k] = -999.0;
1813 
1814  fprimaryPID_Pdg[k] = -999.0;
1815  fprimaryPID_Ndf[k] = -999.0;
1816  fprimaryPID_MinChi2[k] = -999.0;
1817  fprimaryPID_DeltaChi2[k] = -999.0;
1818  fprimaryPID_Chi2Proton[k] = -999.0;
1819  fprimaryPID_Chi2Kaon[k] = -999.0;
1820  fprimaryPID_Chi2Pion[k] = -999.0;
1821  fprimaryPID_Chi2Muon[k] = -999.0;
1822  fprimaryPID_MissingE[k] = -999.0;
1823  fprimaryPID_MissingEavg[k] = -999.0;
1824  fprimaryPID_PIDA[k] = -999.0;
1825  }
1826 
1827  fbeamtrack_truth_Origin = -999;
1828  fbeamtrigger = -999;
1829  ftof = -999.0;
1830  fbeamNTracks = -999;
1831  ftof_expElec = -999.0;
1832  ftof_expMuon = -999.0;
1833  ftof_expPion = -999.0;
1834  ftof_expKaon = -999.0;
1835  ftof_expProt = -999.0;
1836  ftof_expDeut = -999.0;
1837  for(int k=0; k < 2; k++){
1838  fcerenkovPressure[k] = -999.0;
1839  fcerenkovTime[k] = -999.0;
1840  fcerenkovStatus[k] = -999;
1841  }
1842  fbeamtrackMomentum = -999.0;
1843  fbeamtrackEnergy = 999.0;
1844  fbeamtrackPdg = -999;
1845  fbeamtrackID = -999;
1846  fbeamtrack_truth_Pt = -999.0;
1847  fbeamtrack_truth_Mass = -999.0;
1848  fbeamtrack_truth_Theta = -999.0;
1849  fbeamtrack_truth_Phi = -999.0;
1850  fbeamtrack_truth_TotalLength = -999.0;
1851  fbeamtrack_truth_KinEnergy = -999.0;
1852  fbeamtrack_truth_P_InTPCActive = -999.0;
1853  fbeamtrack_truth_Pt_InTPCActive = -999.0;
1854  fbeamtrack_truth_Theta_InTPCActive = -999.0;
1855  fbeamtrack_truth_Phi_InTPCActive = -999.0;
1856  fbeamtrack_truth_TotalLength_InTPCActive = -999.0;
1857  fbeamtrack_truth_KinEnergy_InTPCActive = -999.0;
1858 
1859  for(int l=0; l < 3; l++){
1860  fbeamtrackPos[l] = -999.0;
1861  fbeamtrackDir[l] = -999.0;
1862  }
1863  for(int l=0; l < 4; l++){
1864  fbeamtrack_truth_EndPos[l] = -999.0;
1865  fbeamtrack_truth_Momentum[l] = -999.0;
1866  fbeamtrack_truth_Pos_InTPCActive[l] = -999.0;
1867  fbeamtrack_truth_Momentum_InTPCActive[l] = -999.0;
1868  }
1869  fbeamtrack_truth_NDAUGTHERS = 0;
1870  fbeamtrack_truth_NDECAYDAUGTHERS = 0;
1871  for(int k=0; k < NMAXTRUTHDAUGTHERS; k++){
1872  fbeamtrack_truthdaughter_TrackId[k] = -999;
1873  fbeamtrack_truthdaughter_Pdg[k] = -999;
1874  fbeamtrack_truthdaughter_Mother[k] = -999;
1875  fbeamtrack_truthdaughter_P[k] = -999.0;
1876  fbeamtrack_truthdaughter_Pt[k] = -999.0;
1877  fbeamtrack_truthdaughter_Mass[k] = -999.0;
1878  fbeamtrack_truthdaughter_Theta[k] = -999.0;
1879  fbeamtrack_truthdaughter_Phi[k] = -999.0;
1880  fbeamtrack_truthdaughter_TotalLength[k] = -999.0;
1881  fbeamtrack_truthdaughter_Process[k] = -999;
1882  fbeamtrack_truthdaughter_EndProcess[k] = -999;
1883 
1884  fbeamtrack_truthdecaydaughter_TrackId[k] = -999;
1885  fbeamtrack_truthdecaydaughter_Pdg[k] = -999;
1886  fbeamtrack_truthdecaydaughter_Mother[k] = -999;
1887  fbeamtrack_truthdecaydaughter_P[k] = -999.0;
1888  fbeamtrack_truthdecaydaughter_Pt[k] = -999.0;
1889  fbeamtrack_truthdecaydaughter_Mass[k] = -999.0;
1890  fbeamtrack_truthdecaydaughter_Theta[k] = -999.0;
1891  fbeamtrack_truthdecaydaughter_Phi[k] = -999.0;
1892  fbeamtrack_truthdecaydaughter_TotalLength[k] = -999.0;
1893  fbeamtrack_truthdecaydaughter_Process[k] = -999;
1894  fbeamtrack_truthdecaydaughter_EndProcess[k] = -999;
1895  for(int l=0; l < 4; l++){
1896  fbeamtrack_truthdaughter_StartPosition[k][l] = -999.0;
1897  fbeamtrack_truthdaughter_EndPosition[k][l] = -999.0;
1898  fbeamtrack_truthdaughter_Momentum[k][l] = -999.0;
1899  fbeamtrack_truthdaughter_EndMomentum[k][l] = -999.0;
1900 
1901  fbeamtrack_truthdecaydaughter_StartPosition[k][l] = -999.0;
1902  fbeamtrack_truthdecaydaughter_EndPosition[k][l] = -999.0;
1903  fbeamtrack_truthdecaydaughter_Momentum[k][l] = -999.0;
1904  fbeamtrack_truthdecaydaughter_EndMomentum[k][l] = -999.0;
1905  }
1906  }
1907 
1908  fprimaryIstrack = 0;
1909  fprimaryIsshower = 0;
1910 
1911  fprimaryBDTScore = -999.0;
1912  fprimaryNHits = -999;
1913  fprimaryTheta = -999.0;
1914  fprimaryPhi = -999.0;
1915  fprimaryLength = -999.0;
1916  fprimaryMomentum = -999.0;
1917  fprimaryEndMomentum = -999.0;
1918  fprimaryOpeningAngle = -999.0;
1919  fprimaryShowerBestPlane = -999;
1920  fprimaryID = -999;
1921  fprimaryMomentumByRangeProton = -999.0;
1922  fprimaryMomentumByRangeMuon = -999.0;
1923  fprimaryT0 = -999.0;
1924  fprimaryEarliestHitPeakTime = -999.0;
1925 
1926  fprimary_truth_Origin = -999;
1927  fprimary_truth_TrackId = -999;
1928  fprimary_truth_Pdg = -999;
1929  fprimary_truth_Mother = -999;
1930  fprimary_truth_P = -999.0;
1931  fprimary_truth_Pt = -999.0;
1932  fprimary_truth_Mass = -999.0;
1933  fprimary_truth_Theta = -999.0;
1934  fprimary_truth_Phi = -999.0;
1935  fprimary_truth_Process = -999;
1936  fprimary_truth_Isbeammatched = -999;
1937  fprimary_truth_TotalLength = -999;
1938 
1939  for(int l=0; l < 4; l++){
1940  fprimary_truth_StartPosition[l] = -999.0;
1941  fprimary_truth_EndPosition[l] = -999.0;
1942  fprimary_truth_Momentum[l] = -999.0;
1943  fprimary_truth_EndMomentum[l] = -999.0;
1944  fprimary_truth_Pos_InTPCActive[l] = -999.0;
1945  fprimary_truth_Momentum_InTPCActive[l] = -999.0;
1946  }
1947  for(int l=0; l < 3; l++){
1948  fprimaryShowerEnergy[l] = -999;
1949  fprimaryShowerMIPEnergy[l] = -999;
1950  fprimaryShowerdEdx[l] = -999;
1951  }
1952 
1953  fprimary_truth_NDAUGTHERS = 0;
1954  fprimary_truth_NDECAYDAUGTHERS = 0;
1955  for(int k=0; k < NMAXTRUTHDAUGTHERS; k++){
1956  fprimary_truthdaughter_TrackId[k] = -999;
1957  fprimary_truthdaughter_Pdg[k] = -999;
1958  fprimary_truthdaughter_Mother[k] = -999;
1959  fprimary_truthdaughter_P[k] = -999.0;
1960  fprimary_truthdaughter_Pt[k] = -999.0;
1961  fprimary_truthdaughter_Mass[k] = -999.0;
1962  fprimary_truthdaughter_Theta[k] = -999.0;
1963  fprimary_truthdaughter_Phi[k] = -999.0;
1964  fprimary_truthdaughter_TotalLength[k] = -999.0;
1965  fprimary_truthdaughter_Process[k] = -999;
1966  fprimary_truthdaughter_EndProcess[k] = -999;
1967 
1968  fprimary_truthdecaydaughter_TrackId[k] = -999;
1969  fprimary_truthdecaydaughter_Pdg[k] = -999;
1970  fprimary_truthdecaydaughter_Mother[k] = -999;
1971  fprimary_truthdecaydaughter_P[k] = -999.0;
1972  fprimary_truthdecaydaughter_Pt[k] = -999.0;
1973  fprimary_truthdecaydaughter_Mass[k] = -999.0;
1974  fprimary_truthdecaydaughter_Theta[k] = -999.0;
1975  fprimary_truthdecaydaughter_Phi[k] = -999.0;
1976  fprimary_truthdecaydaughter_TotalLength[k] = -999.0;
1977  fprimary_truthdecaydaughter_Process[k] = -999;
1978  fprimary_truthdecaydaughter_EndProcess[k] = -999;
1979  for(int l=0; l < 4; l++){
1980  fprimary_truthdaughter_StartPosition[k][l] = -999.0;
1981  fprimary_truthdaughter_EndPosition[k][l] = -999.0;
1982  fprimary_truthdaughter_Momentum[k][l] = -999.0;
1983  fprimary_truthdaughter_EndMomentum[k][l] = -999.0;
1984 
1985  fprimary_truthdecaydaughter_StartPosition[k][l] = -999.0;
1986  fprimary_truthdecaydaughter_EndPosition[k][l] = -999.0;
1987  fprimary_truthdecaydaughter_Momentum[k][l] = -999.0;
1988  fprimary_truthdecaydaughter_EndMomentum[k][l] = -999.0;
1989  }
1990  }
1991 
1992  fNDAUGHTERS = 0;
1993  fNGRANDDAUGHTERS = 0;
1994  for(int k=0; k < NMAXDAUGTHERS; k++){
1995  fdaughterIstrack[k] = -999;
1996  fdaughterIsshower[k] = -999;
1997  fdaughterNHits[k] = -999;
1998  fdaughterTheta[k] = -999.0;
1999  fdaughterPhi[k] = -999.0;
2000  fdaughterLength[k] = -999.0;
2001  fdaughterMomentum[k] = -999.0;
2002  fdaughterEndMomentum[k] = -999.0;
2003  for(int l=0; l < 3; l++){
2004  fdaughterEndPosition[k][l] = -999.0;
2005  fdaughterStartPosition[k][l] = -999.0;
2006  fdaughterEndDirection[k][l] = -999.0;
2007  fdaughterStartDirection[k][l] = -999.0;
2008  fdaughterKineticEnergy[k][l] = -999.0;
2009  fdaughterRange[k][l] = -999.0;
2010  fdaughterTrkPitchC[k][l] = -999.0;
2011 
2012  fdaughterPID_Pdg[k][l] = -999.0;
2013  fdaughterPID_Ndf[k][l] = -999.0;
2014  fdaughterPID_MinChi2[k][l] = -999.0;
2015  fdaughterPID_DeltaChi2[k][l] = -999.0;
2016  fdaughterPID_Chi2Proton[k][l] = -999.0;
2017  fdaughterPID_Chi2Kaon[k][l] = -999.0;
2018  fdaughterPID_Chi2Pion[k][l] = -999.0;
2019  fdaughterPID_Chi2Muon[k][l] = -999.0;
2020  fdaughterPID_MissingE[k][l] = -999.0;
2021  fdaughterPID_MissingEavg[k][l] = -999.0;
2022  fdaughterPID_PIDA[k][l] = -999.0;
2023  fgranddaughterVertex[k][l] = -999.0;
2024 
2025  fdaughterShowerEnergy[k][l] = -999;
2026  fdaughterShowerMIPEnergy[k][l] = -999;
2027  fdaughterShowerdEdx[k][l] = -999;
2028  }
2029  fdaughterOpeningAngle[k] = -999.0;
2030  fdaughterShowerBestPlane[k] = -999;
2031  fdaughterMomentumByRangeProton[k] = -999.0;
2032  fdaughterMomentumByRangeMuon[k] = -999.0;
2033  fdaughterID[k] = -999;
2034  fdaughterT0[k] = -999;
2035 
2036  fdaughter_truth_TrackId[k] = -999;
2037  fdaughter_truth_Pdg[k] = -999;
2038  fdaughter_truth_Mother[k] = -999;
2039  fdaughter_truth_P[k] = -999.0;
2040  fdaughter_truth_Pt[k] = -999.0;
2041  fdaughter_truth_Mass[k] = -999.0;
2042  fdaughter_truth_Theta[k] = -999.0;
2043  fdaughter_truth_Phi[k] = -999.0;
2044  fdaughter_truth_Process[k] = -999;
2045  fdaughter_truth_TotalLength[k] = -999.0;
2046  for(int l=0; l < 4; l++){
2047  fdaughter_truth_StartPosition[k][l] = -999.0;
2048  fdaughter_truth_EndPosition[k][l] = -999.0;
2049  fdaughter_truth_Momentum[k][l] = -999.0;
2050  fdaughter_truth_EndMomentum[k][l] = -999.0;
2051  }
2052 
2053  // Grand-daugthers
2054  fgranddaughterIstrack[k] = -999;
2055  fgranddaughterIsshower[k] = -999;
2056  fgranddaughterNHits[k] = -999;
2057  fgranddaughterTheta[k] = -999.0;
2058  fgranddaughterPhi[k] = -999.0;
2059  fgranddaughterLength[k] = -999.0;
2060  fgranddaughterMomentum[k] = -999.0;
2061  fgranddaughterEndMomentum[k] = -999.0;
2062  for(int l=0; l < 3; l++){
2063  fgranddaughterEndPosition[k][l] = -999.0;
2064  fgranddaughterStartPosition[k][l] = -999.0;
2065  fgranddaughterEndDirection[k][l] = -999.0;
2066  fgranddaughterStartDirection[k][l] = -999.0;
2067  fgranddaughterKineticEnergy[k][l] = -999.0;
2068  fgranddaughterRange[k][l] = -999.0;
2069  fgranddaughterTrkPitchC[k][l] = -999.0;
2070 
2071  fgranddaughterPID_Pdg[k][l] = -999.0;
2072  fgranddaughterPID_Ndf[k][l] = -999.0;
2073  fgranddaughterPID_MinChi2[k][l] = -999.0;
2074  fgranddaughterPID_DeltaChi2[k][l] = -999.0;
2075  fgranddaughterPID_Chi2Proton[k][l] = -999.0;
2076  fgranddaughterPID_Chi2Kaon[k][l] = -999.0;
2077  fgranddaughterPID_Chi2Pion[k][l] = -999.0;
2078  fgranddaughterPID_Chi2Muon[k][l] = -999.0;
2079  fgranddaughterPID_MissingE[k][l] = -999.0;
2080  fgranddaughterPID_MissingEavg[k][l] = -999.0;
2081  fgranddaughterPID_PIDA[k][l] = -999.0;
2082 
2083  fgranddaughterShowerEnergy[k][l] = -999;
2084  fgranddaughterShowerMIPEnergy[k][l] = -999;
2085  fgranddaughterShowerdEdx[k][l] = -999;
2086  }
2087  fgranddaughterOpeningAngle[k] = -999.0;
2088  fgranddaughterShowerBestPlane[k] = -999;
2089  fgranddaughterMomentumByRangeProton[k] = -999.0;
2090  fgranddaughterMomentumByRangeMuon[k] = -999.0;
2091  fgranddaughterID[k] = -999;
2092  fgranddaughterMotherID[k] = -999;
2093  fgranddaughterT0[k] = -999;
2094 
2095  fgranddaughter_truth_TrackId[k] = -999;
2096  fgranddaughter_truth_Pdg[k] = -999;
2097  fgranddaughter_truth_Mother[k] = -999;
2098  fgranddaughter_truth_P[k] = -999.0;
2099  fgranddaughter_truth_Pt[k] = -999.0;
2100  fgranddaughter_truth_Mass[k] = -999.0;
2101  fgranddaughter_truth_Theta[k] = -999.0;
2102  fgranddaughter_truth_Phi[k] = -999.0;
2103  fgranddaughter_truth_Process[k] = -999;
2104  fgranddaughter_truth_TotalLength[k] = -999.0;
2105  for(int l=0; l < 4; l++){
2106  fgranddaughter_truth_StartPosition[k][l] = -999.0;
2107  fgranddaughter_truth_EndPosition[k][l] = -999.0;
2108  fgranddaughter_truth_Momentum[k][l] = -999.0;
2109  fgranddaughter_truth_EndMomentum[k][l] = -999.0;
2110  }
2111 
2112  }
2113 
2114 
2115 }
2116 
2117 DEFINE_ART_MODULE(protoana::ProtoDUNEAnalysisTree)
simb::MCParticle::E
double E(const int i=0) const
Definition: MCParticle.h:233
recob::Shower::best_plane
int best_plane() const
Definition: Shower.h:200
LArG4Parameters.h
Store parameters for running LArG4.
recob::PFParticle::Daughters
const std::vector< size_t > & Daughters() const
Returns the collection of daughter particles.
Definition: PFParticle.h:114
recob::Shower::ShowerStart
const TVector3 & ShowerStart() const
Definition: Shower.h:192
protoana::ProtoDUNETrackUtils::GetRecoTrackPID
std::vector< anab::ParticleID > GetRecoTrackPID(const recob::Track &track, art::Event const &evt, const std::string trackModule, const std::string pidModule) const
Get the PID from a given track.
Definition: ProtoDUNETrackUtils.cxx:223
recob::Track::EndMomentum
double EndMomentum() const
Definition: Track.h:144
simb::MCParticle::PdgCode
int PdgCode() const
Definition: MCParticle.h:212
art::Timestamp::timeLow
constexpr std::uint32_t timeLow() const
Definition: Timestamp.h:29
simb::MCParticle::Py
double Py(const int i=0) const
Definition: MCParticle.h:231
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
protoana::ProtoDUNEAnalysisTree::fprimary_truthdecaydaughter_Pt
double fprimary_truthdecaydaughter_Pt[NMAXTRUTHDAUGTHERS]
Definition: ProtoDUNEAnalysisTree_module.cc:351
protoana::ProtoDUNEAnalysisTree::fbeamtrack_truthdaughter_Process
int fbeamtrack_truthdaughter_Process[NMAXTRUTHDAUGTHERS]
Definition: ProtoDUNEAnalysisTree_module.cc:235
protoana::ProtoDUNEAnalysisTree::fbeamtrack_truthdecaydaughter_Momentum
double fbeamtrack_truthdecaydaughter_Momentum[NMAXTRUTHDAUGTHERS][4]
Definition: ProtoDUNEAnalysisTree_module.cc:245
protoana::ProtoDUNETruthUtils::GetGeantGoodParticle
const simb::MCParticle * GetGeantGoodParticle(const simb::MCTruth &genTruth, const art::Event &evt) const
Definition: ProtoDUNETruthUtils.cxx:860
protoana::ProtoDUNEAnalysisTree::fbeamtrack_truthdaughter_Momentum
double fbeamtrack_truthdaughter_Momentum[NMAXTRUTHDAUGTHERS][4]
Definition: ProtoDUNEAnalysisTree_module.cc:228
simb::MCParticle::EndE
double EndE() const
Definition: MCParticle.h:244
simb::MCParticle::EndZ
double EndZ() const
Definition: MCParticle.h:228
protoana::ProtoDUNEAnalysisTree::fprimary_truthdecaydaughter_TotalLength
double fprimary_truthdecaydaughter_TotalLength[NMAXTRUTHDAUGTHERS]
Definition: ProtoDUNEAnalysisTree_module.cc:355
protoana::ProtoDUNEAnalysisTree::fprimary_truthdaughter_EndProcess
int fprimary_truthdaughter_EndProcess[NMAXTRUTHDAUGTHERS]
Definition: ProtoDUNEAnalysisTree_module.cc:340
protoana::ProtoDUNEAnalysisTree::fprimary_truthdecaydaughter_StartPosition
double fprimary_truthdecaydaughter_StartPosition[NMAXTRUTHDAUGTHERS][4]
Definition: ProtoDUNEAnalysisTree_module.cc:346
protoana::ProtoDUNEAnalysisTree::ProtoDUNEAnalysisTree
ProtoDUNEAnalysisTree(fhicl::ParameterSet const &p)
Definition: ProtoDUNEAnalysisTree_module.cc:474
recob::Shower::Length
double Length() const
Definition: Shower.h:201
recob::Track::ParticleId
int ParticleId() const
Definition: Track.h:171
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
protoana::ProtoDUNEAnalysisTree::fprimary_truthdaughter_TotalLength
double fprimary_truthdaughter_TotalLength[NMAXTRUTHDAUGTHERS]
Definition: ProtoDUNEAnalysisTree_module.cc:338
art::DataViewImpl::getHandle
Handle< PROD > getHandle(SelectorBase const &) const
Definition: DataViewImpl.h:382
protoana::ProtoDUNEAnalysisTree::fgranddaughter_truth_StartPosition
double fgranddaughter_truth_StartPosition[NMAXDAUGTHERS][4]
Definition: ProtoDUNEAnalysisTree_module.cc:458
string
std::string string
Definition: nybbler.cc:12
protoana::ProtoDUNEAnalysisTree::fgranddaughter_truth_EndMomentum
double fgranddaughter_truth_EndMomentum[NMAXDAUGTHERS][4]
Definition: ProtoDUNEAnalysisTree_module.cc:462
protoana::ProtoDUNEAnalysisTree::fprimary_truthdecaydaughter_P
double fprimary_truthdecaydaughter_P[NMAXTRUTHDAUGTHERS]
Definition: ProtoDUNEAnalysisTree_module.cc:348
simb::MCParticle::Trajectory
const simb::MCTrajectory & Trajectory() const
Definition: MCParticle.h:253
simb::MCParticle::Mother
int Mother() const
Definition: MCParticle.h:213
recob::Track::Trajectory
const recob::TrackTrajectory & Trajectory() const
Access to the stored recob::TrackTrajectory.
Definition: Track.h:98
protoana::ProtoDUNEAnalysisTree::fbeamtrack_truthdaughter_TrackId
int fbeamtrack_truthdaughter_TrackId[NMAXTRUTHDAUGTHERS]
Definition: ProtoDUNEAnalysisTree_module.cc:222
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::ProtoDUNEAnalysisTree::fbeamtrack_truthdaughter_StartPosition
double fbeamtrack_truthdaughter_StartPosition[NMAXTRUTHDAUGTHERS][4]
Definition: ProtoDUNEAnalysisTree_module.cc:225
protoana::ProtoDUNEAnalysisTree::fbeamtrack_truthdecaydaughter_Process
int fbeamtrack_truthdecaydaughter_Process[NMAXTRUTHDAUGTHERS]
Definition: ProtoDUNEAnalysisTree_module.cc:252
art::Timestamp::timeHigh
constexpr std::uint32_t timeHigh() const
Definition: Timestamp.h:34
simb::MCParticle::Mass
double Mass() const
Definition: MCParticle.h:239
protoana::ProtoDUNEAnalysisTree::fbeamtrack_truthdaughter_Theta
double fbeamtrack_truthdaughter_Theta[NMAXTRUTHDAUGTHERS]
Definition: ProtoDUNEAnalysisTree_module.cc:232
geo::TPCGeo
Geometry information for a single TPC.
Definition: TPCGeo.h:38
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::ProtoDUNEAnalysisTree::fbeamtrack_truthdecaydaughter_Mother
int fbeamtrack_truthdecaydaughter_Mother[NMAXTRUTHDAUGTHERS]
Definition: ProtoDUNEAnalysisTree_module.cc:241
protoana::ProtoDUNEAnalysisTree::FillPrimaryDaughterPFParticle
void FillPrimaryDaughterPFParticle(art::Event const &evt, detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const recob::PFParticle *daughterParticle, int daughterID)
Definition: ProtoDUNEAnalysisTree_module.cc:1124
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
protoana::ProtoDUNEAnalysisTree::fgranddaughterPID_Chi2Pion
double fgranddaughterPID_Chi2Pion[NMAXDAUGTHERS][3]
Definition: ProtoDUNEAnalysisTree_module.cc:449
protoana::ProtoDUNEAnalysisTree::fprimary_truthdecaydaughter_Theta
double fprimary_truthdecaydaughter_Theta[NMAXTRUTHDAUGTHERS]
Definition: ProtoDUNEAnalysisTree_module.cc:353
protoana::ProtoDUNEAnalysisTree::fbeamtrack_truthdecaydaughter_Phi
double fbeamtrack_truthdecaydaughter_Phi[NMAXTRUTHDAUGTHERS]
Definition: ProtoDUNEAnalysisTree_module.cc:250
protoana::ProtoDUNEAnalysisTree::fprimary_truthdaughter_Phi
double fprimary_truthdaughter_Phi[NMAXTRUTHDAUGTHERS]
Definition: ProtoDUNEAnalysisTree_module.cc:337
std
STL namespace.
protoana::ProtoDUNEAnalysisTree::fprimary_truthdecaydaughter_EndMomentum
double fprimary_truthdecaydaughter_EndMomentum[NMAXTRUTHDAUGTHERS][4]
Definition: ProtoDUNEAnalysisTree_module.cc:350
protoana::ProtoDUNEAnalysisTree::fprimary_truthdaughter_Process
int fprimary_truthdaughter_Process[NMAXTRUTHDAUGTHERS]
Definition: ProtoDUNEAnalysisTree_module.cc:339
protoana::ProtoDUNEAnalysisTree::fprimary_truthdaughter_StartPosition
double fprimary_truthdaughter_StartPosition[NMAXTRUTHDAUGTHERS][4]
Definition: ProtoDUNEAnalysisTree_module.cc:329
protoana::ProtoDUNEAnalysisTree::fprimary_truthdecaydaughter_Phi
double fprimary_truthdecaydaughter_Phi[NMAXTRUTHDAUGTHERS]
Definition: ProtoDUNEAnalysisTree_module.cc:354
protoana::ProtoDUNEAnalysisTree::fprimary_truthdaughter_Pt
double fprimary_truthdaughter_Pt[NMAXTRUTHDAUGTHERS]
Definition: ProtoDUNEAnalysisTree_module.cc:334
art::EDAnalyzer::EDAnalyzer
EDAnalyzer(fhicl::ParameterSet const &pset)
Definition: EDAnalyzer.h:25
protoana::ProtoDUNEPFParticleUtils::GetPFParticleEarliestHitPeakTime
const double GetPFParticleEarliestHitPeakTime(const recob::PFParticle &particle, art::Event const &evt, const std::string particleLabel) const
Get the earliest hit peak time.
Definition: ProtoDUNEPFParticleUtils.cxx:681
geo::GeometryCore::Ncryostats
unsigned int Ncryostats() const
Returns the number of cryostats in the detector.
Definition: GeometryCore.h:1994
simb::MCParticle::Process
std::string Process() const
Definition: MCParticle.h:215
protoana::ProtoDUNEAnalysisTree::fprimary_truthdaughter_Momentum
double fprimary_truthdaughter_Momentum[NMAXTRUTHDAUGTHERS][4]
Definition: ProtoDUNEAnalysisTree_module.cc:332
protoana::ProtoDUNEAnalysisTree::fprimary_truthdaughter_EndMomentum
double fprimary_truthdaughter_EndMomentum[NMAXTRUTHDAUGTHERS][4]
Definition: ProtoDUNEAnalysisTree_module.cc:333
MCParticle.h
Particle class.
simb::MCParticle::EndY
double EndY() const
Definition: MCParticle.h:227
protoana::ProtoDUNEAnalysisTree::FillPrimaryGrandDaughterPFParticle
void FillPrimaryGrandDaughterPFParticle(art::Event const &evt, detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const recob::PFParticle *gdaughterParticle, int daughterID, int gdaughterID)
Definition: ProtoDUNEAnalysisTree_module.cc:1277
protoana::ProtoDUNEAnalysisTree::fbeamtrack_truthdecaydaughter_EndProcess
int fbeamtrack_truthdecaydaughter_EndProcess[NMAXTRUTHDAUGTHERS]
Definition: ProtoDUNEAnalysisTree_module.cc:253
l
static QStrList * l
Definition: config.cpp:1044
protoana::ProtoDUNEAnalysisTree::fgranddaughterStartDirection
double fgranddaughterStartDirection[NMAXDAUGTHERS][3]
Definition: ProtoDUNEAnalysisTree_module.cc:428
protoana::ProtoDUNEAnalysisTree::pi_service
art::ServiceHandle< cheat::ParticleInventoryService > pi_service
Definition: ProtoDUNEAnalysisTree_module.cc:105
protoana::ProtoDUNEAnalysisTree::fgranddaughterPID_MissingEavg
double fgranddaughterPID_MissingEavg[NMAXDAUGTHERS][3]
Definition: ProtoDUNEAnalysisTree_module.cc:452
protoana::ProtoDUNEAnalysisTree::fbeamtrack_truthdecaydaughter_Pdg
int fbeamtrack_truthdecaydaughter_Pdg[NMAXTRUTHDAUGTHERS]
Definition: ProtoDUNEAnalysisTree_module.cc:240
protoana::ProtoDUNEAnalysisTree::fbeamtrack_truthdaughter_TotalLength
double fbeamtrack_truthdaughter_TotalLength[NMAXTRUTHDAUGTHERS]
Definition: ProtoDUNEAnalysisTree_module.cc:234
Geometry.h
art framework interface to geometry description
protoana::ProtoDUNEAnalysisTree::fbeamtrack_truthdecaydaughter_P
double fbeamtrack_truthdecaydaughter_P[NMAXTRUTHDAUGTHERS]
Definition: ProtoDUNEAnalysisTree_module.cc:244
art::Run
Definition: Run.h:17
simb::MCParticle::TrackId
int TrackId() const
Definition: MCParticle.h:210
protoana::ProtoDUNEAnalysisTree::fgranddaughterShowerEnergy
double fgranddaughterShowerEnergy[NMAXDAUGTHERS][3]
Definition: ProtoDUNEAnalysisTree_module.cc:430
protoana::ProtoDUNETruthUtils::GetKinEnergyAtVertex
double GetKinEnergyAtVertex(const simb::MCParticle &mcpart, double kinene_lastpoint=0.0)
Definition: ProtoDUNETruthUtils.cxx:1069
geo::GeometryCore::Nchannels
unsigned int Nchannels() const
Returns the number of TPC readout channels in the detector.
Definition: GeometryCore.cxx:208
protoana::ProtoDUNETruthUtils::GetDepEnergyAtLastTrajPoint
double GetDepEnergyAtLastTrajPoint(const simb::MCParticle &mcpart, double tpcactiveXlow, double tpcactiveXhigh, double tpcactiveYlow, double tpcactiveYhigh, double tpcactiveZlow, double tpcactiveZhigh)
Definition: ProtoDUNETruthUtils.cxx:1029
protoana::ProtoDUNETruthUtils::GetMCParticleLengthInTPCActiveVolume
double GetMCParticleLengthInTPCActiveVolume(const simb::MCParticle &mcpart, double tpcactiveXlow, double tpcactiveXhigh, double tpcactiveYlow, double tpcactiveYhigh, double tpcactiveZlow, double tpcactiveZhigh)
Definition: ProtoDUNETruthUtils.cxx:1008
protoana::ProtoDUNEAnalysisTree::fprimary_truthdecaydaughter_Pdg
int fprimary_truthdecaydaughter_Pdg[NMAXTRUTHDAUGTHERS]
Definition: ProtoDUNEAnalysisTree_module.cc:344
util::size
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:92
geo::GeometryCore::GDMLFile
std::string GDMLFile() const
Returns the full directory path to the GDML file source.
Definition: GeometryCore.h:1749
protoana::ProtoDUNEAnalysisTree::FillPrimaryPFParticle
void FillPrimaryPFParticle(art::Event const &evt, detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const recob::PFParticle *particle)
Definition: ProtoDUNEAnalysisTree_module.cc:842
art::DataViewImpl::isRealData
bool isRealData() const
Definition: DataViewImpl.cc:113
protoana::ProtoDUNEAnalysisTree::fgranddaughterStartPosition
double fgranddaughterStartPosition[NMAXDAUGTHERS][3]
Definition: ProtoDUNEAnalysisTree_module.cc:426
simb::MCParticle::Pt
double Pt(const int i=0) const
Definition: MCParticle.h:236
NMAXTRUTHDAUGTHERS
const int NMAXTRUTHDAUGTHERS
Definition: ProtoDUNEAnalysisTree_module.cc:72
protoana::ProtoDUNEAnalysisTree::fgranddaughterMomentumByRangeProton
double fgranddaughterMomentumByRangeProton[NMAXDAUGTHERS]
Definition: ProtoDUNEAnalysisTree_module.cc:434
geo::GeometryCore::Nplanes
unsigned int Nplanes(unsigned int tpc=0, unsigned int cstat=0) const
Returns the total number of wire planes in the specified TPC.
Definition: GeometryCore.h:2848
recob::Track::Phi
double Phi() const
Definition: Track.h:178
recob::Track::Length
double Length(size_t p=0) const
Access to various track properties.
Definition: Track.h:167
protoana::ProtoDUNEAnalysisTree::fprimary_truthdecaydaughter_Process
int fprimary_truthdecaydaughter_Process[NMAXTRUTHDAUGTHERS]
Definition: ProtoDUNEAnalysisTree_module.cc:356
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
protoana::ProtoDUNEAnalysisTree::fprimary_truthdaughter_TrackId
int fprimary_truthdaughter_TrackId[NMAXTRUTHDAUGTHERS]
Definition: ProtoDUNEAnalysisTree_module.cc:326
recob::Shower::MIPEnergy
const std::vector< double > & MIPEnergy() const
Definition: Shower.h:198
simb::MCParticle::EndProcess
std::string EndProcess() const
Definition: MCParticle.h:216
protoana::ProtoDUNEAnalysisTree::beamlineUtil
protoana::ProtoDUNEBeamlineUtils beamlineUtil
Definition: ProtoDUNEAnalysisTree_module.cc:104
geo::GeometryCore::DetectorName
std::string DetectorName() const
Returns a string with the name of the detector, as configured.
Definition: GeometryCore.h:1761
NMAXDAUGTHERS
const int NMAXDAUGTHERS
Definition: ProtoDUNEAnalysisTree_module.cc:71
simb::MCParticle::EndPz
double EndPz() const
Definition: MCParticle.h:243
protoana::ProtoDUNEAnalysisTree::fgranddaughterKineticEnergy
double fgranddaughterKineticEnergy[NMAXDAUGTHERS][3]
Definition: ProtoDUNEAnalysisTree_module.cc:436
geo::GeometryCore::IterateTPCs
IteratorBox< TPC_iterator,&GeometryCore::begin_TPC,&GeometryCore::end_TPC > IterateTPCs() const
Enables ranged-for loops on all TPCs of the detector.
Definition: GeometryCore.h:2655
cheat::ParticleInventoryService::TrackIdToMCTruth_P
const art::Ptr< simb::MCTruth > & TrackIdToMCTruth_P(int id) const
Definition: ParticleInventoryService_service.cc:142
protoana::ProtoDUNEAnalysisTree::fprimary_truthdaughter_P
double fprimary_truthdaughter_P[NMAXTRUTHDAUGTHERS]
Definition: ProtoDUNEAnalysisTree_module.cc:331
recob::Track::Theta
double Theta() const
Access to spherical or geographical angles at vertex or at any point.
Definition: Track.h:176
protoana::ProtoDUNEAnalysisTree::fgranddaughterShowerMIPEnergy
double fgranddaughterShowerMIPEnergy[NMAXDAUGTHERS][3]
Definition: ProtoDUNEAnalysisTree_module.cc:431
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
protoana::ProtoDUNEAnalysisTree::fprimary_truthdaughter_Theta
double fprimary_truthdaughter_Theta[NMAXTRUTHDAUGTHERS]
Definition: ProtoDUNEAnalysisTree_module.cc:336
simb::MCParticle::T
double T(const int i=0) const
Definition: MCParticle.h:224
protoana::ProtoDUNEAnalysisTree::fbeamtrack_truthdecaydaughter_TrackId
int fbeamtrack_truthdecaydaughter_TrackId[NMAXTRUTHDAUGTHERS]
Definition: ProtoDUNEAnalysisTree_module.cc:239
protoana::ProtoDUNEAnalysisTree::fgranddaughterPID_MissingE
double fgranddaughterPID_MissingE[NMAXDAUGTHERS][3]
Definition: ProtoDUNEAnalysisTree_module.cc:451
protoana::ProtoDUNEAnalysisTree::fprimary_truthdaughter_Mother
int fprimary_truthdaughter_Mother[NMAXTRUTHDAUGTHERS]
Definition: ProtoDUNEAnalysisTree_module.cc:328
art::DataViewImpl::getValidHandle
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
Definition: DataViewImpl.h:441
protoana::ProtoDUNEAnalysisTree::fbeamtrack_truthdaughter_Mass
double fbeamtrack_truthdaughter_Mass[NMAXTRUTHDAUGTHERS]
Definition: ProtoDUNEAnalysisTree_module.cc:231
test.p
p
Definition: test.py:223
protoana::ProtoDUNEAnalysisTree::fprimary_truthdecaydaughter_Momentum
double fprimary_truthdecaydaughter_Momentum[NMAXTRUTHDAUGTHERS][4]
Definition: ProtoDUNEAnalysisTree_module.cc:349
protoana::ProtoDUNEAnalysisTree::fgranddaughter_truth_EndPosition
double fgranddaughter_truth_EndPosition[NMAXDAUGTHERS][4]
Definition: ProtoDUNEAnalysisTree_module.cc:459
art::DataViewImpl::subRun
SubRunNumber_t subRun() const
Definition: DataViewImpl.cc:78
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
recob::Shower::Direction
const TVector3 & Direction() const
Definition: Shower.h:189
protoana::ProtoDUNEAnalysisTree::fbeamtrack_truthdaughter_Mother
int fbeamtrack_truthdaughter_Mother[NMAXTRUTHDAUGTHERS]
Definition: ProtoDUNEAnalysisTree_module.cc:224
protoana::ProtoDUNEAnalysisTree::fbeamtrack_truthdaughter_Pt
double fbeamtrack_truthdaughter_Pt[NMAXTRUTHDAUGTHERS]
Definition: ProtoDUNEAnalysisTree_module.cc:230
art::DataViewImpl::run
RunNumber_t run() const
Definition: DataViewImpl.cc:71
simb::MCParticle::EndT
double EndT() const
Definition: MCParticle.h:229
protoana::ProtoDUNEAnalysisTree::fbeamtrack_truthdecaydaughter_Pt
double fbeamtrack_truthdecaydaughter_Pt[NMAXTRUTHDAUGTHERS]
Definition: ProtoDUNEAnalysisTree_module.cc:247
geo::GeometryCore
Description of geometry of one entire detector.
Definition: GeometryCore.h:1468
protoana::ProtoDUNEAnalysisTree::fgranddaughterPID_MinChi2
double fgranddaughterPID_MinChi2[NMAXDAUGTHERS][3]
Definition: ProtoDUNEAnalysisTree_module.cc:445
protoana::ProtoDUNEAnalysisTree::fprimary_truthdecaydaughter_Mass
double fprimary_truthdecaydaughter_Mass[NMAXTRUTHDAUGTHERS]
Definition: ProtoDUNEAnalysisTree_module.cc:352
geo_types.h
Definition of data types for geometry description.
protoana::ProtoDUNEAnalysisTree::FillPrimaryBeamParticle
bool FillPrimaryBeamParticle(art::Event const &evt)
Definition: ProtoDUNEAnalysisTree_module.cc:626
protoana::ProtoDUNEAnalysisTree::fprimary_truthdaughter_Mass
double fprimary_truthdaughter_Mass[NMAXTRUTHDAUGTHERS]
Definition: ProtoDUNEAnalysisTree_module.cc:335
keras_to_tensorflow.int
int
Definition: keras_to_tensorflow.py:69
geo::GeometryCore::NTPC
unsigned int NTPC(unsigned int cstat=0) const
Returns the total number of TPCs in the specified cryostat.
Definition: GeometryCore.h:2303
protoana::ProtoDUNEAnalysisTree::fbeamtrack_truthdecaydaughter_TotalLength
double fbeamtrack_truthdecaydaughter_TotalLength[NMAXTRUTHDAUGTHERS]
Definition: ProtoDUNEAnalysisTree_module.cc:251
protoana::ProtoDUNEAnalysisTree::fprimary_truthdaughter_EndPosition
double fprimary_truthdaughter_EndPosition[NMAXTRUTHDAUGTHERS][4]
Definition: ProtoDUNEAnalysisTree_module.cc:330
recob::TrackTrajectory::EndDirection
Vector_t EndDirection() const
Returns the direction of the trajectory at the last point.
Definition: TrackTrajectory.h:321
tracks
Definition: tracks.py:1
protoana::ProtoDUNEAnalysisTree::fbeamtrack_truthdecaydaughter_EndMomentum
double fbeamtrack_truthdecaydaughter_EndMomentum[NMAXTRUTHDAUGTHERS][4]
Definition: ProtoDUNEAnalysisTree_module.cc:246
simb::MCParticle::Vx
double Vx(const int i=0) const
Definition: MCParticle.h:221
protoana::ProtoDUNEAnalysisTree::fbeamtrack_truthdaughter_EndPosition
double fbeamtrack_truthdaughter_EndPosition[NMAXTRUTHDAUGTHERS][4]
Definition: ProtoDUNEAnalysisTree_module.cc:226
protoana::ProtoDUNEAnalysisTree::fbeamtrack_truthdecaydaughter_Mass
double fbeamtrack_truthdecaydaughter_Mass[NMAXTRUTHDAUGTHERS]
Definition: ProtoDUNEAnalysisTree_module.cc:248
Hit.h
Declaration of signal hit object.
recob::Track::StartMomentum
double StartMomentum() const
Definition: Track.h:143
protoana::ProtoDUNEAnalysisTree::fgranddaughterPID_Chi2Kaon
double fgranddaughterPID_Chi2Kaon[NMAXDAUGTHERS][3]
Definition: ProtoDUNEAnalysisTree_module.cc:448
recob::PFParticle
Hierarchical representation of particle flow.
Definition: PFParticle.h:44
protoana::ProtoDUNEAnalysisTree::fbeamtrack_truthdaughter_EndMomentum
double fbeamtrack_truthdaughter_EndMomentum[NMAXTRUTHDAUGTHERS][4]
Definition: ProtoDUNEAnalysisTree_module.cc:229
protoana::ProtoDUNEAnalysisTree::fdaughter_truth_EndPosition
double fdaughter_truth_EndPosition[NMAXDAUGTHERS][4]
Definition: ProtoDUNEAnalysisTree_module.cc:403
protoana::ProtoDUNEAnalysisTree::fdaughter_truth_StartPosition
double fdaughter_truth_StartPosition[NMAXDAUGTHERS][4]
Definition: ProtoDUNEAnalysisTree_module.cc:402
detinfo::DetectorClocksData
Contains all timing reference information for the detector.
Definition: DetectorClocksData.h:283
simb::MCParticle::EndPy
double EndPy() const
Definition: MCParticle.h:242
protoana::ProtoDUNEAnalysisTree::fdaughter_truth_EndMomentum
double fdaughter_truth_EndMomentum[NMAXDAUGTHERS][4]
Definition: ProtoDUNEAnalysisTree_module.cc:406
protoana::ProtoDUNEDataUtils::GetNActiveFembsForAPA
int GetNActiveFembsForAPA(art::Event const &evt, int apa) const
Get number of active fembs in an APA.
Definition: ProtoDUNEDataUtils.cxx:49
simb::MCTrajectory::TotalLength
double TotalLength() const
Definition: MCTrajectory.cxx:52
protoana::ProtoDUNEAnalysisTree::fgranddaughter_truth_Momentum
double fgranddaughter_truth_Momentum[NMAXDAUGTHERS][4]
Definition: ProtoDUNEAnalysisTree_module.cc:461
protoana::ProtoDUNEBeamlineUtils::ComputeTOF
double ComputeTOF(int pdg, double momentum)
Definition: ProtoDUNEBeamlineUtils.cxx:614
protoana::ProtoDUNEAnalysisTree::fbeamtrack_truthdecaydaughter_EndPosition
double fbeamtrack_truthdecaydaughter_EndPosition[NMAXTRUTHDAUGTHERS][4]
Definition: ProtoDUNEAnalysisTree_module.cc:243
protoana::ProtoDUNEAnalysisTree::pfpUtil
protoana::ProtoDUNEPFParticleUtils pfpUtil
Definition: ProtoDUNEAnalysisTree_module.cc:101
recob::TrackTrajectory::End
Point_t const & End() const
Returns the position of the last valid point of the trajectory [cm].
Definition: TrackTrajectory.h:257
simb::MCParticle::Momentum
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:220
wirecell.gen.depos.center
def center(depos, point)
Definition: depos.py:117
simb::MCParticle::Pz
double Pz(const int i=0) const
Definition: MCParticle.h:232
Track.h
Provides recob::Track data product.
simb::MCParticle::Vz
double Vz(const int i=0) const
Definition: MCParticle.h:223
protoana::ProtoDUNEAnalysisTree::analyze
void analyze(art::Event const &evt) override
Definition: ProtoDUNEAnalysisTree_module.cc:498
protoana::ProtoDUNETruthUtils::GetFirstTrajectoryPointInTPCActiveVolume
int GetFirstTrajectoryPointInTPCActiveVolume(const simb::MCParticle &mcpart, double tpcactiveXlow, double tpcactiveXhigh, double tpcactiveYlow, double tpcactiveYhigh, double tpcactiveZlow, double tpcactiveZhigh)
Definition: ProtoDUNETruthUtils.cxx:993
protoana::ProtoDUNEAnalysisTree::fgranddaughterPID_DeltaChi2
double fgranddaughterPID_DeltaChi2[NMAXDAUGTHERS][3]
Definition: ProtoDUNEAnalysisTree_module.cc:446
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::ProtoDUNEAnalysisTree::fgranddaughterPID_Chi2Muon
double fgranddaughterPID_Chi2Muon[NMAXDAUGTHERS][3]
Definition: ProtoDUNEAnalysisTree_module.cc:450
protoana::ProtoDUNEAnalysisTree::fbeamtrack_truthdaughter_Pdg
int fbeamtrack_truthdaughter_Pdg[NMAXTRUTHDAUGTHERS]
Definition: ProtoDUNEAnalysisTree_module.cc:223
protoana::ProtoDUNEAnalysisTree::fprimary_truthdecaydaughter_EndPosition
double fprimary_truthdecaydaughter_EndPosition[NMAXTRUTHDAUGTHERS][4]
Definition: ProtoDUNEAnalysisTree_module.cc:347
GeometryCore.h
Access the description of detector geometry.
simb::MCParticle::EndPx
double EndPx() const
Definition: MCParticle.h:241
protoana::ProtoDUNEAnalysisTree::fgranddaughterEndDirection
double fgranddaughterEndDirection[NMAXDAUGTHERS][3]
Definition: ProtoDUNEAnalysisTree_module.cc:427
protoana::ProtoDUNEAnalysisTree::fbeamtrack_truthdecaydaughter_Theta
double fbeamtrack_truthdecaydaughter_Theta[NMAXTRUTHDAUGTHERS]
Definition: ProtoDUNEAnalysisTree_module.cc:249
protoana::ProtoDUNEAnalysisTree::operator=
ProtoDUNEAnalysisTree & operator=(ProtoDUNEAnalysisTree const &)=delete
tca::evt
TCEvent evt
Definition: DataStructs.cxx:7
protoana::ProtoDUNEAnalysisTree::fbeamtrack_truthdaughter_P
double fbeamtrack_truthdaughter_P[NMAXTRUTHDAUGTHERS]
Definition: ProtoDUNEAnalysisTree_module.cc:227
protoana::ProtoDUNETruthUtils::GetMCParticleFromRecoShower
const simb::MCParticle * GetMCParticleFromRecoShower(detinfo::DetectorClocksData const &clockData, const recob::Shower &shower, art::Event const &evt, std::string showerModule) const
Definition: ProtoDUNETruthUtils.cxx:653
art::get
auto const & get(AssnsNode< L, R, D > const &r)
Definition: AssnsNode.h:115
protoana::ProtoDUNEAnalysisTree::fprimary_truthdecaydaughter_Mother
int fprimary_truthdecaydaughter_Mother[NMAXTRUTHDAUGTHERS]
Definition: ProtoDUNEAnalysisTree_module.cc:345
art::fill_ptr_vector
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:297
protoana::ProtoDUNEAnalysisTree::fprimary_truthdecaydaughter_TrackId
int fprimary_truthdecaydaughter_TrackId[NMAXTRUTHDAUGTHERS]
Definition: ProtoDUNEAnalysisTree_module.cc:343
protoana::ProtoDUNEAnalysisTree::beginRun
void beginRun(const art::Run &run) override
Definition: ProtoDUNEAnalysisTree_module.cc:493
simb::MCParticle::EndX
double EndX() const
Definition: MCParticle.h:226
protoana::ProtoDUNEAnalysisTree::fbeamtrack_truthdaughter_Phi
double fbeamtrack_truthdaughter_Phi[NMAXTRUTHDAUGTHERS]
Definition: ProtoDUNEAnalysisTree_module.cc:233
recob::TrackTrajectory::StartDirection
Vector_t StartDirection() const
Returns the direction of the trajectory at the first point.
Definition: TrackTrajectory.h:317
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::ProtoDUNEAnalysisTree::fbeamtrack_truthdaughter_EndProcess
int fbeamtrack_truthdaughter_EndProcess[NMAXTRUTHDAUGTHERS]
Definition: ProtoDUNEAnalysisTree_module.cc:236
protoana::ProtoDUNEAnalysisTree::fbeamtrack_truthdecaydaughter_StartPosition
double fbeamtrack_truthdecaydaughter_StartPosition[NMAXTRUTHDAUGTHERS][4]
Definition: ProtoDUNEAnalysisTree_module.cc:242
recob::Shower::ID
int ID() const
Definition: Shower.h:187
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
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::ProtoDUNEAnalysisTree::fprimary_truthdecaydaughter_EndProcess
int fprimary_truthdecaydaughter_EndProcess[NMAXTRUTHDAUGTHERS]
Definition: ProtoDUNEAnalysisTree_module.cc:357
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
geo::origin
constexpr Point origin()
Returns a origin position with a point of the specified type.
Definition: geo_vectors.h:227
trkf::TrackMomentumCalculator::GetTrackMomentum
double GetTrackMomentum(double trkrange, int pdg) const
Definition: TrackMomentumCalculator.cxx:144
protoana::ProtoDUNETruthUtils::GetProcessKey
int GetProcessKey(std::string process)
Definition: ProtoDUNETruthUtils.cxx:921
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
protoana::ProtoDUNEDataUtils::IsBeamTrigger
bool IsBeamTrigger(art::Event const &evt) const
Definition: ProtoDUNEDataUtils.cxx:27
protoana::ProtoDUNEAnalysisTree::fgranddaughterEndPosition
double fgranddaughterEndPosition[NMAXDAUGTHERS][3]
Definition: ProtoDUNEAnalysisTree_module.cc:425
protoana::ProtoDUNEAnalysisTree::fgranddaughterPID_Chi2Proton
double fgranddaughterPID_Chi2Proton[NMAXDAUGTHERS][3]
Definition: ProtoDUNEAnalysisTree_module.cc:447
Generated by   doxygen 1.8.11