ProtoDUNEelectronAnaTree_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: ProtoDUNEelectronAnaTree
3 // File: ProtoDUNEelectronAnaTree_module.cc
4 //
5 // Extract protoDUNE useful information, do a firs tpre-selection and save output to a flat tree
6 //
7 //
8 // Georgios Christodoulou - georgios.christodoulou at cern.ch
9 // modified by Aaron Higuera ahiguera@central.uh.edu
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 "lardataobj/RecoBase/PFParticleMetadata.h"
20 #include "art_root_io/TFileService.h"
21 #include "canvas/Utilities/InputTag.h"
22 #include "fhiclcpp/ParameterSet.h"
23 #include "messagefacility/MessageLogger/MessageLogger.h"
24 
25 #include "larreco/Calorimetry/CalorimetryAlg.h"
26 #include "larsim/MCCheater/BackTrackerService.h"
27 #include "larsim/MCCheater/ParticleInventoryService.h"
28 #include "lardataobj/RecoBase/Wire.h"
29 #include "lardataobj/RecoBase/Hit.h"
30 #include "lardataobj/RecoBase/Track.h"
31 #include "lardataobj/RecoBase/Shower.h"
32 #include "lardataobj/RecoBase/PFParticle.h"
33 #include "nusimdata/SimulationBase/MCParticle.h"
34 #include "nusimdata/SimulationBase/MCTruth.h"
35 #include "lardataobj/AnalysisBase/CosmicTag.h"
36 #include "lardataobj/AnalysisBase/T0.h"
37 #include "larreco/RecoAlg/TrackMomentumCalculator.h"
38 #include "lardataobj/AnalysisBase/Calorimetry.h"
39 #include "larcore/Geometry/Geometry.h"
40 #include "dune/DuneObj/ProtoDUNEBeamEvent.h"
41 #include "lardata/Utilities/GeometryUtilities.h"
42 #include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
43 #include "dune/CalibServices/XYZCalibServiceProtoDUNE.h"
44 #include "lardata/ArtDataHelper/MVAReader.h"
45 #include "canvas/Persistency/Common/FindManyP.h"
46 
47 #include "protoduneana/Utilities/ProtoDUNETrackUtils.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/ProtoDUNEShowerUtils.h"
52 #include "protoduneana/Utilities/ProtoDUNEBeamlineUtils.h"
53 
54 // ROOT includes
55 #include "TTree.h"
56 #include "TFile.h"
57 #include "TString.h"
58 #include "TTimeStamp.h"
59 
60 // C++ Includes
61 #include <stdio.h>
62 #include <stdlib.h>
63 #include <fstream>
64 #include <string>
65 #include <sstream>
66 #include <cmath>
67 #include <algorithm>
68 #include <iostream>
69 #include <vector>
70 
71 // Maximum number of beam particles to save
72 const int NMAXDAUGTHERS = 25;
73 const int MAXHits = 6000;
74 namespace protoana {
75  class ProtoDUNEelectronAnaTree;
76 }
77 
78 
79 class protoana::ProtoDUNEelectronAnaTree : public art::EDAnalyzer {
80 public:
81 
82  explicit ProtoDUNEelectronAnaTree(fhicl::ParameterSet const & p);
83 
84  ProtoDUNEelectronAnaTree(ProtoDUNEelectronAnaTree const &) = delete;
85  ProtoDUNEelectronAnaTree(ProtoDUNEelectronAnaTree &&) = delete;
86  ProtoDUNEelectronAnaTree & operator = (ProtoDUNEelectronAnaTree const &) = delete;
87  ProtoDUNEelectronAnaTree & operator = (ProtoDUNEelectronAnaTree &&) = delete;
88 
89  virtual void beginJob() override;
90  virtual void endJob() override;
91 
92  // Required functions.
93  void analyze(art::Event const & evt) override;
94 
95 private:
96 
97  // Helper utility functions
98  protoana::ProtoDUNEDataUtils dataUtil;
99  protoana::ProtoDUNEPFParticleUtils pfpUtil;
100  protoana::ProtoDUNETrackUtils trackUtil;
101  protoana::ProtoDUNETruthUtils truthUtil;
102  protoana::ProtoDUNEShowerUtils showerUtil;
103  protoana::ProtoDUNEBeamlineUtils beamlineUtil;
104 
105 
106  // Track momentum algorithm calculates momentum based on track range
107  trkf::TrackMomentumCalculator trmom;
108  geo::GeometryCore const * fGeometry = &*(art::ServiceHandle<geo::Geometry>());
109 
110  // Initialize tree variables
111  void Initialize();
112  void FillPrimaryPFParticle(art::Event const & evt,
113  detinfo::DetectorClocksData const& clockData,
114  detinfo::DetectorPropertiesData const& detProp,
115  const recob::PFParticle* particle);
116  void FillPrimaryDaughterPFParticle(art::Event const & evt, const recob::PFParticle* daughterParticle, int daughterID);
117 
118  // Fill cosmics tree
119 
120  // fcl parameters
121  const art::InputTag fBeamModuleLabel;
122  std::string fCalorimetryTag;
123  std::string fParticleIDTag;
124  std::string fTrackerTag;
125  std::string fShowerTag;
126  std::string fHitTag;
127  std::string fShowerCaloTag;
128  std::string fPFParticleTag;
129  std::string fGeneratorTag;
130  int fVerbose;
131 
132  TTree *fPandoraBeam;
133 
134  // Tree variables
135  int fRun;
136  int fSubRun;
137  int fevent;
138  double fTimeStamp;
139  int fNactivefembs[6];
140 
141  // Beam track
142  int fbeamtrigger;
143  int fbeamCheckIsMatched;
144  double ftof;
145  int fcerenkovStatus[2];
146  double fcerenkovTime[2];
147  double fcerenkovPressure[2];
148  double fbeamtrackMomentum;
149  double fbeamtrackEnergy;
150  double fbeamtrackPos[3];
151  double fbeamtrackEndPos[3];
152  double fbeamtrackDir[3];
153  double fbeamtrackTime;
154  int fbeam_ntrjPoints;
155  int fbeamtrackPdg;
156  int fbeamtrackID;
157  double fbeamtrackPos_at[1000][3];
158  double fbeamtrackMom_at[1000][4];
159  std::vector<std::string> fbeamtrackEndProcess;
160  int fbeamtrackNDaughters;
161 
162  // Reconstructed tracks/showers
163  double fprimaryVertex[3];
164  int fprimaryIstrack;
165  int fprimaryIsshower;
166  double fprimaryBDTScore;
167  int fprimaryNHits;
168  double fprimaryTheta;
169  double fprimaryPhi;
170  double fprimaryLength;
171  double fprimaryMomentum;
172  double fprimaryEndMomentum;
173  double fprimaryEndPosition[3];
174  double fprimaryStartPosition[3];
175  double fprimaryEndDirection[3];
176  double fprimaryStartDirection[3];
177  double fprimaryOpeningAngle;
178  int fprimaryShowerBestPlane;
179  double fprimaryShowerEnergy;
180  double fprimaryShowerCharge;
181  double fprimaryShowerMIPEnergy;
182  double fprimaryMomentumByRangeProton;
183  double fprimaryMomentumByRangeMuon;
184  int fprimaryIsBeamparticle;
185  int fprimaryTruth_trkID;
186  int fprimaryTruth_pdg;
187  double fprimaryTruth_Edepo;
188  double fprimaryTruth_purity;
189  double fprimaryTruth_E;
190  double fprimaryTruth_vtx[3];
191  double fprimaryKineticEnergy[3];
192  double fprimaryRange[3];
193  int fprimarynCal;
194  double fprimarydEdx[MAXHits];
195  double fprimarydQdx[MAXHits];
196  double fprimary_calX[MAXHits];
197  double fprimary_calY[MAXHits];
198  double fprimary_calZ[MAXHits];
199  double fprimary_cal_pitch[MAXHits];
200  double fprimaryResidualRange[MAXHits];
201  int fprimaryTruthShower_nHits;
202  double fprimaryShowerTruth_Charge;
203  int fprimaryShower_nHits; //collection only
204  int fprimaryShower_hit_w[MAXHits];
205  double fprimaryShower_hit_q[MAXHits];
206  double fprimaryShower_hit_t[MAXHits];
207  double fprimaryShower_hit_X[MAXHits];
208  int fprimaryShower_hit_ch[MAXHits];
209  double fprimaryShower_hit_Y[MAXHits];
210  double fprimaryShower_hit_Z[MAXHits];
211  double fprimaryShower_hit_pitch[MAXHits];
212  double fprimaryShower_hit_cnn[MAXHits];
213  int fprimaryTruthShower_hit_w[MAXHits];
214  double fprimaryTruthShower_hit_q[MAXHits];
215  double fprimaryTruthShower_hit_t[MAXHits];
216  double fprimaryTruthShower_hit_X[MAXHits];
217  double fprimaryTruthShower_hit_Y[MAXHits];
218  double fprimaryTruthShower_hit_Z[MAXHits];
219 
220 
221  int fprimaryID;
222  double fprimaryT0;
223 
224  int fNDAUGHTERS =0;
225  double fdaughterVertex[3];
226  int fdaughterIstrack[NMAXDAUGTHERS];
227  int fdaughterIsshower[NMAXDAUGTHERS];
228  int fdaughterNHits[NMAXDAUGTHERS];
229  double fdaughterTheta[NMAXDAUGTHERS];
230  double fdaughterPhi[NMAXDAUGTHERS];
231  double fdaughterLength[NMAXDAUGTHERS];
232  double fdaughterEndPosition[NMAXDAUGTHERS][3];
233  double fdaughterStartPosition[NMAXDAUGTHERS][3];
234  double fdaughterEndDirection[NMAXDAUGTHERS][3];
235  double fdaughterStartDirection[NMAXDAUGTHERS][3];
236  double fdaughterOpeningAngle[NMAXDAUGTHERS];
237  int fdaughterShowerBestPlane[NMAXDAUGTHERS];
238  int fdaughterID[NMAXDAUGTHERS];
239  double fdaughterT0[NMAXDAUGTHERS];
240 
241 };
242 
243 
244 protoana::ProtoDUNEelectronAnaTree::ProtoDUNEelectronAnaTree(fhicl::ParameterSet const & p)
245  :
246  EDAnalyzer(p),
247  dataUtil(p.get<fhicl::ParameterSet>("DataUtils")),
248  beamlineUtil(p.get<fhicl::ParameterSet>("BeamLineUtils")),
249  fBeamModuleLabel(p.get< art::InputTag >("BeamModuleLabel")),
250  fCalorimetryTag(p.get<std::string>("CalorimetryTag")),
251  fParticleIDTag(p.get<std::string>("ParticleIDTag")),
252  fTrackerTag(p.get<std::string>("TrackerTag")),
253  fShowerTag(p.get<std::string>("ShowerTag")),
254  fHitTag(p.get<std::string>("HitTag")),
255  fShowerCaloTag(p.get<std::string>("ShowerCalorimetryTag")),
256  fPFParticleTag(p.get<std::string>("PFParticleTag")),
257  fGeneratorTag(p.get<std::string>("GeneratorTag")),
258 
259  fVerbose(p.get<int>("Verbose"))
260 {
261 
262 }
263 
264 void protoana::ProtoDUNEelectronAnaTree::beginJob(){
265 
266  art::ServiceHandle<art::TFileService> tfs;
267 
268  fPandoraBeam = tfs->make<TTree>("PandoraBeam", "Beam events reconstructed with Pandora");
269  fPandoraBeam->Branch("run", &fRun, "run/I");
270  fPandoraBeam->Branch("subrun", &fSubRun, "subrun/I");
271  fPandoraBeam->Branch("event", &fevent, "event/I");
272  fPandoraBeam->Branch("timestamp", &fTimeStamp, "timestamp/D");
273  fPandoraBeam->Branch("Nactivefembs", &fNactivefembs, "Nactivefembs[5]/I");
274 
275  fPandoraBeam->Branch("beamtrigger", &fbeamtrigger, "beamtrigger/I");
276  fPandoraBeam->Branch("beamCheckIsMatched", &fbeamCheckIsMatched, "beamCheckIsMatched/I");
277  fPandoraBeam->Branch("tof", &ftof, "tof/D");
278  fPandoraBeam->Branch("cerenkovStatus", &fcerenkovStatus, "cerenkovStatus[2]/I");
279  fPandoraBeam->Branch("cerenkovTime", &fcerenkovTime, "cerenkovTime[2]/D");
280  fPandoraBeam->Branch("cerenkovPressure", &fcerenkovPressure, "cerenkovPressure[2]/D");
281  fPandoraBeam->Branch("beamtrackMomentum", &fbeamtrackMomentum, "beamtrackMomentum/D");
282  fPandoraBeam->Branch("beamtrackEnergy", &fbeamtrackEnergy, "beamtrackEnergy/D");
283  fPandoraBeam->Branch("beamtrackPos", &fbeamtrackPos, "beamtrackPos[3]/D");
284  fPandoraBeam->Branch("beamtrackEndPos", &fbeamtrackEndPos, "beamtrackEndPos[3]/D");
285  fPandoraBeam->Branch("beamtrackDir", &fbeamtrackDir, "beamtrackDir[3]/D");
286  fPandoraBeam->Branch("beamtrackTime", &fbeamtrackTime, "beamtrackTime/D");
287  fPandoraBeam->Branch("beamtrackPdg", &fbeamtrackPdg, "beamtrackPdg/I");
288  fPandoraBeam->Branch("beamtrackID", &fbeamtrackID, "beamtrackID/I");
289  fPandoraBeam->Branch("beam_ntrjPoints", &fbeam_ntrjPoints, "beam_ntrjPoints/I");
290  fPandoraBeam->Branch("beamtrackNDaughters", &fbeamtrackNDaughters, "beamtrackNDaughters/I");
291  fPandoraBeam->Branch("beamtrackPos_at", &fbeamtrackPos_at, "beamtrackPos_at[beam_ntrjPoints][3]/D");
292  fPandoraBeam->Branch("beamtrackMom_at", &fbeamtrackMom_at, "beamtrackMom_at[beam_ntrjPoints][4]/D");
293  fPandoraBeam->Branch("beamtrackEndProcess", &fbeamtrackEndProcess);
294 
295  fPandoraBeam->Branch("primaryVertex", &fprimaryVertex, "primaryVertex[3]/D");
296  fPandoraBeam->Branch("primaryIsBeamparticle", &fprimaryIsBeamparticle, "primaryIsBeamparticle/I");
297  fPandoraBeam->Branch("primaryIstrack", &fprimaryIstrack, "primaryIstrack/I");
298  fPandoraBeam->Branch("primaryIsshower", &fprimaryIsshower, "primaryIsshower/I");
299  fPandoraBeam->Branch("primaryBDTScore", &fprimaryBDTScore, "primaryBDTScore/D");
300  fPandoraBeam->Branch("primaryNHits", &fprimaryNHits, "primaryNHits/I");
301  fPandoraBeam->Branch("primaryTheta", &fprimaryTheta, "primaryTheta/D");
302  fPandoraBeam->Branch("primaryPhi", &fprimaryPhi, "primaryPhi/D");
303  fPandoraBeam->Branch("primaryLength", &fprimaryLength, "primaryLength/D");
304  fPandoraBeam->Branch("primaryMomentum", &fprimaryMomentum, "primaryMomentum/D");
305  fPandoraBeam->Branch("primaryEndMomentum", &fprimaryEndMomentum, "primaryEndMomentum/D");
306  fPandoraBeam->Branch("primaryEndPosition", &fprimaryEndPosition, "primaryEndPosition[3]/D");
307  fPandoraBeam->Branch("primaryStartPosition", &fprimaryStartPosition, "primaryStartPosition[3]/D");
308  fPandoraBeam->Branch("primaryEndDirection", &fprimaryEndDirection, "primaryEndDirection[3]/D");
309  fPandoraBeam->Branch("primaryStartDirection", &fprimaryStartDirection, "primaryStartDirection[3]/D");
310  fPandoraBeam->Branch("primaryOpeningAngle", &fprimaryOpeningAngle, "primaryOpeningAngle/D");
311  fPandoraBeam->Branch("primaryID", &fprimaryID, "primaryID/I");
312  fPandoraBeam->Branch("primaryTruth_Edepo", &fprimaryTruth_Edepo, "primaryTruth_Edepo/D");
313  fPandoraBeam->Branch("primaryTruth_purity", &fprimaryTruth_purity, "primaryTruth_purity/D");
314  fPandoraBeam->Branch("primaryTruth_E", &fprimaryTruth_E, "primaryTruth_E/D");
315  fPandoraBeam->Branch("primaryTruth_vtx", &fprimaryTruth_vtx, "primaryTruth_vtx[3]/D");
316  fPandoraBeam->Branch("primaryTruth_pdg", &fprimaryTruth_pdg, "primaryTruth_pdg/I");
317  fPandoraBeam->Branch("primaryTruth_trkID", &fprimaryTruth_trkID, "primaryTruth_trkID/I");
318  fPandoraBeam->Branch("primaryShowerBestPlane", &fprimaryShowerBestPlane, "primaryShowerBestPlane/I");
319  fPandoraBeam->Branch("primaryShowerCharge", &fprimaryShowerCharge, "primaryShowerCharge/D");
320  fPandoraBeam->Branch("primaryShowerEnergy", &fprimaryShowerEnergy, "primaryShowerEnergy/D");
321  fPandoraBeam->Branch("primaryShowerMIPEnergy", &fprimaryShowerMIPEnergy, "primaryShowerMIPEnergy/D");
322 
323  fPandoraBeam->Branch("primaryShower_nHits", &fprimaryShower_nHits, "primaryShower_nHits/I");
324  fPandoraBeam->Branch("primaryTruthShower_nHits", &fprimaryTruthShower_nHits, "primaryTruthShower_nHits/I");
325  fPandoraBeam->Branch("primaryShowerTruth_Charge", &fprimaryShowerTruth_Charge, "primaryShowerTruth_Charge/D");
326  fPandoraBeam->Branch("primaryShower_hit_q", &fprimaryShower_hit_q, "primaryShower_hit_q[primaryShower_nHits]/D");
327  fPandoraBeam->Branch("primaryShower_hit_w", &fprimaryShower_hit_w, "primaryShower_hit_w[primaryShower_nHits]/I");
328  fPandoraBeam->Branch("primaryShower_hit_t", &fprimaryShower_hit_t, "primaryShower_hit_t[primaryShower_nHits]/D");
329  fPandoraBeam->Branch("primaryShower_hit_X", &fprimaryShower_hit_X, "primaryShower_hit_X[primaryShower_nHits]/D");
330  fPandoraBeam->Branch("primaryShower_hit_Y", &fprimaryShower_hit_Y, "primaryShower_hit_Y[primaryShower_nHits]/D");
331  fPandoraBeam->Branch("primaryShower_hit_Z", &fprimaryShower_hit_Z, "primaryShower_hit_Z[primaryShower_nHits]/D");
332  fPandoraBeam->Branch("primaryShower_hit_ch", &fprimaryShower_hit_ch, "primaryShower_hit_ch[primaryShower_nHits]/I");
333 
334  fPandoraBeam->Branch("primaryTruthShower_hit_q", &fprimaryTruthShower_hit_q, "primaryTruthShower_hit_q[primaryTruthShower_nHits]/D");
335  fPandoraBeam->Branch("primaryTruthShower_hit_w", &fprimaryTruthShower_hit_w, "primaryTruthShower_hit_w[primaryTruthShower_nHits]/I");
336  fPandoraBeam->Branch("primaryTruthShower_hit_t", &fprimaryTruthShower_hit_t, "primaryTruthShower_hit_t[primaryTruthShower_nHits]/D");
337  fPandoraBeam->Branch("primaryTruthShower_hit_X", &fprimaryTruthShower_hit_X, "primaryTruthShower_hit_X[primaryTruthShower_nHits]/D");
338  fPandoraBeam->Branch("primaryTruthShower_hit_Y", &fprimaryTruthShower_hit_Y, "primaryTruthShower_hit_Y[primaryTruthShower_nHits]/D");
339  fPandoraBeam->Branch("primaryTruthShower_hit_Z", &fprimaryTruthShower_hit_Z, "primaryTruthShower_hit_Z[primaryTruthShower_nHits]/D");
340 
341  fPandoraBeam->Branch("primaryShower_hit_pitch", &fprimaryShower_hit_pitch, "primaryShower_hit_pitch[primaryShower_nHits]/D");
342  fPandoraBeam->Branch("primaryShower_hit_cnn", &fprimaryShower_hit_cnn, "primaryShower_hit_cnn[primaryShower_nHits]/D");
343 
344  fPandoraBeam->Branch("primaryMomentumByRangeProton", &fprimaryMomentumByRangeProton, "primaryMomentumByRangeProton/D");
345  fPandoraBeam->Branch("primaryMomentumByRangeMuon", &fprimaryMomentumByRangeMuon, "primaryMomentumByRangeMuon/D");
346  fPandoraBeam->Branch("primaryKineticEnergy", &fprimaryKineticEnergy, "primaryKineticEnergy[3]/D");
347  fPandoraBeam->Branch("primaryRange", &fprimaryRange, "primaryRange[3]/D");
348  fPandoraBeam->Branch("primarynCal", &fprimarynCal, "primarynCal/I");
349  fPandoraBeam->Branch("primarydQdx", &fprimarydQdx, "primarydQdx[primarynCal]/D");
350  fPandoraBeam->Branch("primary_calX", &fprimary_calX, "primary_calX[primarynCal]/D");
351  fPandoraBeam->Branch("primary_calY", &fprimary_calY, "primary_calY[primarynCal]/D");
352  fPandoraBeam->Branch("primary_calZ", &fprimary_calZ, "primary_calZ[primarynCal]/D");
353  fPandoraBeam->Branch("primary_cal_pitch", &fprimary_cal_pitch, "primary_cal_pitch[primarynCal]/D");
354  fPandoraBeam->Branch("primarydEdx", &fprimarydEdx, "primarydEdx[primarynCal]/D");
355  fPandoraBeam->Branch("primaryResidualRange", &fprimaryResidualRange, "primaryResidualRange[primarynCal]/D");
356  fPandoraBeam->Branch("primaryT0", &fprimaryT0, "primaryT0/D");
357 
358  fPandoraBeam->Branch("NDAUGHTERS", &fNDAUGHTERS, "NDAUGHTERS/I");
359  fPandoraBeam->Branch("daughterVertex", &fdaughterVertex, "daughterVertex[3]/D");
360  fPandoraBeam->Branch("daughterIstrack", &fdaughterIstrack, "daughterIstrack[NDAUGHTERS]/I");
361  fPandoraBeam->Branch("daughterIsshower", &fdaughterIsshower, "daughterIsshower[NDAUGHTERS]/I");
362  fPandoraBeam->Branch("daughterNHits", &fdaughterNHits, "daughterNHits[NDAUGHTERS]/I");
363  fPandoraBeam->Branch("daughterTheta", &fdaughterTheta, "daughterTheta[NDAUGHTERS]/D");
364  fPandoraBeam->Branch("daughterPhi", &fdaughterPhi, "daughterPhi[NDAUGHTERS]/D");
365  fPandoraBeam->Branch("daughterLength", &fdaughterLength, "daughterLength[NDAUGHTERS]/D");
366  fPandoraBeam->Branch("daughterEndPosition", &fdaughterEndPosition, "daughterEndPosition[NDAUGHTERS][3]/D");
367  fPandoraBeam->Branch("daughterStartPosition", &fdaughterStartPosition, "daughterStartPosition[NDAUGHTERS][3]/D");
368  fPandoraBeam->Branch("daughterStartDirection", &fdaughterStartDirection, "daughterStartDirection[NDAUGHTERS][3]/D");
369  fPandoraBeam->Branch("daughterEndDirection", &fdaughterEndDirection, "daughterEndDirection[NDAUGHTERS][3]/D");
370  fPandoraBeam->Branch("daughterOpeningAngle", &fdaughterOpeningAngle, "daughterOpeningAngle[NDAUGHTERS]/D");
371  fPandoraBeam->Branch("daughterShowerBestPlane", &fdaughterShowerBestPlane, "daughterShowerBestPlane[NDAUGHTERS]/D");
372  fPandoraBeam->Branch("daughterID", &fdaughterID, "daughterID[NDAUGHTERS]/I");
373  fPandoraBeam->Branch("daughterT0", &fdaughterT0, "daughterT0[NDAUGHTERS]/D");
374 
375 }
376 
377 void protoana::ProtoDUNEelectronAnaTree::analyze(art::Event const & evt){
378 
379  // Initialize tree parameters
380  Initialize();
381  fRun = evt.run();
382  fSubRun = evt.subRun();
383  fevent = evt.id().event();
384  art::Timestamp ts = evt.time();
385  if (ts.timeHigh() == 0){
386  TTimeStamp ts2(ts.timeLow());
387  fTimeStamp = ts2.AsDouble();
388  }
389  else{
390  TTimeStamp ts2(ts.timeHigh(), ts.timeLow());
391  fTimeStamp = ts2.AsDouble();
392  }
393 
394  // Get number of active fembs
395  if(!evt.isRealData()){
396  for(int k=0; k < 6; k++)
397  fNactivefembs[k] = 20;
398  }
399  else{
400  for(int k=0; k < 6; k++)
401  fNactivefembs[k] = dataUtil.GetNActiveFembsForAPA(evt, k);
402  }
403 
404  bool beamTriggerEvent = false;
405  if(!evt.isRealData()){
406  // Firstly we need to get the list of MCTruth objects from the generator. The standard protoDUNE
407  auto mcTruths = evt.getValidHandle<std::vector<simb::MCTruth>>(fGeneratorTag);
408 
409  // Also get the reconstructed beam information in the MC - TO DO
410  const simb::MCParticle* geantGoodParticle = truthUtil.GetGeantGoodParticle((*mcTruths)[0],evt);
411  art::ServiceHandle<cheat::ParticleInventoryService> pi_serv;
412  if(geantGoodParticle != 0x0){
413  beamTriggerEvent = true;
414  fbeamtrigger = 12;
415  fbeamtrackPos[0] = geantGoodParticle->Vx();
416  fbeamtrackPos[1] = geantGoodParticle->Vy();
417  fbeamtrackPos[2] = geantGoodParticle->Vz();
418  fbeamtrackEndPos[0]= geantGoodParticle->EndX();
419  fbeamtrackEndPos[1]= geantGoodParticle->EndY();
420  fbeamtrackEndPos[2]= geantGoodParticle->EndZ();
421  fbeamtrackMomentum = geantGoodParticle->P();
422  fbeamtrackEnergy = geantGoodParticle->E();
423  fbeamtrackPdg = geantGoodParticle->PdgCode();
424  fbeamtrackTime = geantGoodParticle->T();
425  fbeamtrackID = geantGoodParticle->TrackId();
426  fbeamtrackEndProcess.push_back(geantGoodParticle->EndProcess());
427  fbeam_ntrjPoints = geantGoodParticle->NumberTrajectoryPoints();
428  fbeamtrackNDaughters = geantGoodParticle->NumberDaughters();
429  for( size_t i=0; i<geantGoodParticle->NumberTrajectoryPoints() && i<1000; ++i){
430  fbeamtrackPos_at[i][0] = geantGoodParticle->Position(i).X();
431  fbeamtrackPos_at[i][1] = geantGoodParticle->Position(i).Y();
432  fbeamtrackPos_at[i][2] = geantGoodParticle->Position(i).Z();
433  fbeamtrackMom_at[i][0] = geantGoodParticle->Momentum(i).Px();
434  fbeamtrackMom_at[i][1] = geantGoodParticle->Momentum(i).Py();
435  fbeamtrackMom_at[i][2] = geantGoodParticle->Momentum(i).Pz();
436  fbeamtrackMom_at[i][3] = geantGoodParticle->Momentum(i).E();
437  }
438  }
439 
440  if( abs(geantGoodParticle->PdgCode()) != 11 ) return;
441  } // MC
442  else{ //data
443  // For data we can see if this event comes from a beam trigger
444  beamTriggerEvent = dataUtil.IsBeamTrigger(evt);
445  if( !beamTriggerEvent ) return;
446 
447  std::vector< art::Ptr<beam::ProtoDUNEBeamEvent> > beaminfo;
448  auto pdbeamHandle = evt.getHandle< std::vector<beam::ProtoDUNEBeamEvent> >(fBeamModuleLabel);
449  if (pdbeamHandle)
450  art::fill_ptr_vector(beaminfo, pdbeamHandle);
451  else{
452  std::cout<<"No beam information from "<<fBeamModuleLabel<<std::endl;
453  }
454  if(beaminfo.size()){
455  if( beamlineUtil.IsGoodBeamlineTrigger( evt ) ){
456  fbeamCheckIsMatched = beaminfo[0]->CheckIsMatched();
457  fbeamtrigger = beaminfo[0]->GetTimingTrigger();
458  if(beaminfo[0]->GetTOFChan() != -1){//if TOFChan == -1, then there was not a successful match, if it's 0, 1, 2, or 3, then there was a good match corresponding to the different pair-wise combinations of the upstream and downstream channels
459  ftof = beaminfo[0]->GetTOF();
460  }
461  // Get Cerenkov
462  fcerenkovStatus[0] = beaminfo[0]->GetCKov0Status();
463  fcerenkovStatus[1] = beaminfo[0]->GetCKov1Status();
464  fcerenkovTime[0] = beaminfo[0]->GetCKov0Time();
465  fcerenkovTime[1] = beaminfo[0]->GetCKov1Time();
466  fcerenkovPressure[0] = beaminfo[0]->GetCKov0Pressure();
467  fcerenkovPressure[1] = beaminfo[0]->GetCKov1Pressure();
468  auto & tracks = beaminfo[0]->GetBeamTracks();
469  if(!tracks.empty()){
470  fbeamtrackPos[0] = tracks[0].Start().X();
471  fbeamtrackPos[1] = tracks[0].Start().Y();
472  fbeamtrackPos[2] = tracks[0].Start().Z();
473  fbeamtrackEndPos[0] = tracks[0].End().X();
474  fbeamtrackEndPos[1] = tracks[0].End().Y();
475  fbeamtrackEndPos[2] = tracks[0].End().Z();
476  fbeamtrackDir[0] = tracks[0].StartDirection().X();
477  fbeamtrackDir[1] = tracks[0].StartDirection().Y();
478  fbeamtrackDir[2] = tracks[0].StartDirection().Z();
479  }
480  auto & beammom = beaminfo[0]->GetRecoBeamMomenta();
481  if(!beammom.empty()) fbeamtrackMomentum = beammom[0];
482  } //good beam trigger
483  }
484  }//for data
485 
486  //check for reco pandora stuff
487  auto recoParticleHandle = evt.getHandle<std::vector<recob::PFParticle>>(fPFParticleTag);
488  if( !recoParticleHandle ) return;
489 
490  // Get all of the PFParticles, by default from the "pandora" product
491  auto recoParticles = evt.getValidHandle<std::vector<recob::PFParticle>>(fPFParticleTag);
492 
493  // We'd like to find the beam particle. Pandora tries to do this for us, so let's use the PFParticle utility
494  // to look for it. Pandora reconstructs slices containing one (or sometimes more) primary PFParticles. These
495  // are tagged as either beam or cosmic for ProtoDUNE. This function automatically considers only those
496  // PFParticles considered as primary
497  std::vector<const recob::PFParticle*> pfParticles = pfpUtil.GetPFParticlesFromBeamSlice(evt,fPFParticleTag);
498  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService>()->DataFor(evt);
499  auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService>()->DataFor(evt, clockData);
500  for(const recob::PFParticle* particle : pfParticles){
501 
502  FillPrimaryPFParticle(evt, clockData, detProp, particle);
503  const TVector3 vtx = pfpUtil.GetPFParticleVertex(*particle,evt,fPFParticleTag,fTrackerTag);
504 
505  fprimaryVertex[0] = vtx.X(); fprimaryVertex[1] = vtx.Y(); fprimaryVertex[2] = vtx.Z();
506  const TVector3 interactionVtx = pfpUtil.GetPFParticleSecondaryVertex(*particle,evt,fPFParticleTag,fTrackerTag);
507  fdaughterVertex[0] = interactionVtx.X(); fdaughterVertex[1] = interactionVtx.Y(); fdaughterVertex[2] = interactionVtx.Z();
508 
509  fNDAUGHTERS =0;
510  for( const int didx : particle->Daughters() ){
511  const recob::PFParticle *daughterParticle =&(recoParticles->at(didx));
512  FillPrimaryDaughterPFParticle(evt, daughterParticle, didx);
513  // Only process NMAXDAUGTHERS
514  if(fNDAUGHTERS > NMAXDAUGTHERS) break;
515  fNDAUGHTERS ++;
516  }
517 
518  // For now only consider the first primary track. Need a proper treatment if more than one primary particles are found.
519  break;
520  }
521  if(beamTriggerEvent) fPandoraBeam->Fill();
522 }
523 
524 // -----------------------------------------------------------------------------
525 void protoana::ProtoDUNEelectronAnaTree::endJob(){
526 
527 }
528 
529 // -----------------------------------------------------------------------------
530 void protoana::ProtoDUNEelectronAnaTree::FillPrimaryPFParticle(art::Event const & evt,
531  detinfo::DetectorClocksData const& clockData,
532  detinfo::DetectorPropertiesData const& detProp,
533  const recob::PFParticle* particle){
534 
535  // Pandora's BDT beam-cosmic score
536  fprimaryBDTScore = (double)pfpUtil.GetBeamCosmicScore(*particle,evt,fPFParticleTag);
537  // NHits associated with this pfParticle
538  fprimaryNHits = (pfpUtil.GetPFParticleHits(*particle,evt,fPFParticleTag)).size();
539 
540  // Get the T0 for this pfParticle
541  std::vector<anab::T0> pfT0vec = pfpUtil.GetPFParticleT0(*particle,evt,fPFParticleTag);
542  if(!pfT0vec.empty())
543  fprimaryT0 = pfT0vec[0].Time();
544 
545  // "particle" is the pointer to our beam particle. The recob::Track or recob::Shower object
546  // of this particle might be more helpful. These return null pointers if not track-like / shower-like
547  const recob::Track* thisTrack = pfpUtil.GetPFParticleTrack(*particle, evt,fPFParticleTag,fTrackerTag);
548  const recob::Shower* thisShower = pfpUtil.GetPFParticleShower(*particle,evt,fPFParticleTag,fShowerTag);
549 
550  const simb::MCParticle* mcparticle = NULL;
551  if(thisTrack != 0x0){
552 
553  mcparticle = truthUtil.GetMCParticleFromRecoTrack(clockData, *thisTrack, evt, fTrackerTag);
554  if(mcparticle){
555  fprimaryTruth_pdg = mcparticle->PdgCode();
556  fprimaryTruth_trkID = mcparticle->TrackId();
557  fprimaryTruth_vtx[0] =mcparticle->Vx();
558  fprimaryTruth_vtx[1] =mcparticle->Vy();
559  fprimaryTruth_vtx[2] =mcparticle->Vx();
560  fprimaryTruth_E = mcparticle->E();
561  if( fbeamtrackID != -999 && fbeamtrackID == fprimaryTruth_trkID )
562  fprimaryIsBeamparticle = 1;
563  }
564  fprimaryIstrack = 1;
565  fprimaryIsshower = 0;
566  fprimaryID = thisTrack->ParticleId();
567  fprimaryTheta = thisTrack->Theta();
568  fprimaryPhi = thisTrack->Phi();
569  fprimaryLength = thisTrack->Length();
570  fprimaryMomentum = thisTrack->StartMomentum();
571  fprimaryEndMomentum = thisTrack->EndMomentum();
572  fprimaryEndPosition[0] = thisTrack->Trajectory().End().X();
573  fprimaryEndPosition[1] = thisTrack->Trajectory().End().Y();
574  fprimaryEndPosition[2] = thisTrack->Trajectory().End().Z();
575  fprimaryStartPosition[0] = thisTrack->Trajectory().Start().X();
576  fprimaryStartPosition[1] = thisTrack->Trajectory().Start().Y();
577  fprimaryStartPosition[2] = thisTrack->Trajectory().Start().Z();
578  fprimaryEndDirection[0] = thisTrack->Trajectory().EndDirection().X();
579  fprimaryEndDirection[1] = thisTrack->Trajectory().EndDirection().Y();
580  fprimaryEndDirection[2] = thisTrack->Trajectory().EndDirection().Z();
581  fprimaryStartDirection[0] = thisTrack->Trajectory().StartDirection().X();
582  fprimaryStartDirection[1] = thisTrack->Trajectory().StartDirection().Y();
583  fprimaryStartDirection[2] = thisTrack->Trajectory().StartDirection().Z();
584 
585  fprimaryMomentumByRangeMuon = trmom.GetTrackMomentum(thisTrack->Length(),13);
586  fprimaryMomentumByRangeProton = trmom.GetTrackMomentum(thisTrack->Length(),2212);
587 
588  // Calorimetry only colleciton plane
589  std::vector<anab::Calorimetry> calovector = trackUtil.GetRecoTrackCalorimetry(*thisTrack, evt, fTrackerTag, fCalorimetryTag);
590  for(size_t k = 0; k < calovector.size() && k<3; k++){
591  int plane = calovector[k].PlaneID().Plane;
592  if(plane !=2 ) continue;
593  fprimaryKineticEnergy[plane] = calovector[k].KineticEnergy();
594  fprimaryRange[plane] = calovector[k].Range();
595  fprimarynCal = calovector[k].dEdx().size();
596  for(size_t l=0; l<calovector[k].dEdx().size() && l<MAXHits; ++l){
597  fprimarydEdx[l]= calovector[k].dEdx()[l];
598  fprimarydQdx[l]= calovector[k].dQdx()[l];
599  fprimaryResidualRange[l]= calovector[k].ResidualRange()[l];
600  const auto &pos=(calovector[k].XYZ())[l];
601  fprimary_calX[l] = pos.X();
602  fprimary_calY[l] = pos.Y();
603  fprimary_calZ[l] = pos.Z();
604  fprimary_cal_pitch[l] =calovector[k].TrkPitchVec()[l];
605  }
606  }
607  } // end is track
608  else if(thisShower != 0x0){
609 
610  mcparticle = truthUtil.GetMCParticleFromRecoShower(clockData, *thisShower, evt, fShowerTag);
611  if(mcparticle){
612  fprimaryTruth_pdg = mcparticle->PdgCode();
613  fprimaryTruth_trkID = mcparticle->TrackId();
614  fprimaryTruth_vtx[0] =mcparticle->Vx();
615  fprimaryTruth_vtx[1] =mcparticle->Vy();
616  fprimaryTruth_vtx[2] =mcparticle->Vx();
617  fprimaryTruth_E = mcparticle->E();
618  if( fbeamtrackID != -999 && fbeamtrackID == fprimaryTruth_trkID )
619  fprimaryIsBeamparticle = 1;
620 
621  double Edepo = truthUtil.GetDepEnergyMC(evt, fGeometry, mcparticle->TrackId(), 2 );
622  double purity = truthUtil.GetPurity(clockData, *thisShower, evt, fShowerTag);
623  fprimaryTruth_Edepo = Edepo;
624  fprimaryTruth_purity = purity; //not really useful in this case
625 
626  double tot_ch =0;
627  std::vector<const recob::Hit*> hitsFromMCPart = truthUtil.GetMCParticleHits(clockData, *mcparticle, evt, fHitTag);
628  art::FindManyP<recob::SpacePoint> spFromMCPartHits(hitsFromMCPart,evt,fPFParticleTag);
629 
630  if( hitsFromMCPart.size()){
631  int n_hits =0;
632  for( size_t i=0; i<hitsFromMCPart.size(); ++i){
633  if( hitsFromMCPart[i]->WireID().Plane != 2 ) continue;
634  const geo::WireGeo* pwire = fGeometry->WirePtr(hitsFromMCPart[i]->WireID());
635  TVector3 xyzWire = pwire->GetCenter<TVector3>();
636  tot_ch += hitsFromMCPart[i]->Integral();
637  fprimaryTruthShower_hit_w[n_hits]=hitsFromMCPart[i]->WireID().Wire;
638  fprimaryTruthShower_hit_t[n_hits]=hitsFromMCPart[i]->PeakTime();
639  fprimaryTruthShower_hit_q[n_hits]=hitsFromMCPart[i]->Integral();
640  fprimaryTruthShower_hit_X[n_hits]=detProp.ConvertTicksToX(hitsFromMCPart[i]->PeakTime(),hitsFromMCPart[i]->WireID().Plane,hitsFromMCPart[i]->WireID().TPC,0);
641  fprimaryTruthShower_hit_Z[n_hits] = xyzWire.Z();
642  std::vector<art::Ptr<recob::SpacePoint>> sp = spFromMCPartHits.at(i);
643  if(!sp.empty() ){
644  //fprimaryShower_hit_X[idx]= sp[0]->XYZ()[0];
645  fprimaryTruthShower_hit_Y[n_hits]= sp[0]->XYZ()[1];
646  //fprimaryTruthShower_hit_Z[n_hits]= sp[0]->XYZ()[2];
647  }
648  n_hits ++;
649  }
650  fprimaryShowerTruth_Charge = tot_ch;
651  fprimaryTruthShower_nHits = n_hits; //collection only
652  }
653  }
654  fprimaryIstrack = 0;
655  fprimaryIsshower = 1;
656  fprimaryID = thisShower->ID();
657  fprimaryLength = thisShower->Length();
658  fprimaryShowerBestPlane = thisShower->best_plane();
659  fprimaryOpeningAngle = thisShower->OpenAngle();
660  fprimaryStartPosition[0] = thisShower->ShowerStart().X();
661  fprimaryStartPosition[1] = thisShower->ShowerStart().Y();
662  fprimaryStartPosition[2] = thisShower->ShowerStart().Z();
663  fprimaryStartDirection[0] = thisShower->Direction().X();
664  fprimaryStartDirection[1] = thisShower->Direction().Y();
665  fprimaryStartDirection[2] = thisShower->Direction().Z();
666 
667  const std::vector<const recob::Hit*> sh_hits = showerUtil.GetRecoShowerHits(*thisShower, evt, fShowerTag);
668  art::FindManyP<recob::SpacePoint> spFromShowerHits(sh_hits,evt,fPFParticleTag);
669 
670  //work around to save the CNN score
671  auto recoShowers = evt.getValidHandle< std::vector< recob::Shower > >(fShowerTag);
672  art::FindManyP<recob::Hit> findHitsFromShowers(recoShowers,evt,fShowerTag);
673  anab::MVAReader<recob::Hit,4> hitResults(evt, "emtrkmichelid:emtrkmichel" );
674  //int trk_idx = hitResults.getIndex("track");
675  //int em_idx = hitResults.getIndex("em");
676  std::vector< art::Ptr< recob::Hit > > tmp_sh_hits = findHitsFromShowers.at( thisShower->ID() );
677 
678  int idx =0;
679  fprimaryShowerCharge =0.0;
680  for( size_t j=0; j<sh_hits.size() && j<MAXHits; ++j){
681  if( sh_hits[j]->WireID().Plane != 2 ) continue;
682  art::FindManyP<recob::Wire> wFromHits(sh_hits,evt,"hitpdune");
683  std::vector<art::Ptr<recob::Wire>> wires = wFromHits.at(j);
684  const geo::WireGeo* pwire = fGeometry->WirePtr(sh_hits[j]->WireID());
685  TVector3 xyzWire = pwire->GetCenter<TVector3>();
686  std::array<float,4> cnn_out = hitResults.getOutput( tmp_sh_hits[j] );
687  double p_trk_or_sh = cnn_out[ hitResults.getIndex("track") ]+ cnn_out[ hitResults.getIndex("em") ];
688  double cnn_score = cnn_out[ hitResults.getIndex("em") ]/p_trk_or_sh;
689  fprimaryShower_hit_cnn[idx] = cnn_score;
690  fprimaryShowerCharge += sh_hits[j]->Integral();
691  fprimaryShower_hit_w[idx]=sh_hits[j]->WireID().Wire;
692  fprimaryShower_hit_t[idx]=sh_hits[j]->PeakTime();
693  fprimaryShower_hit_q[idx]=sh_hits[j]->Integral();
694  fprimaryShower_hit_X[idx]=detProp.ConvertTicksToX(sh_hits[j]->PeakTime(),sh_hits[j]->WireID().Plane,sh_hits[j]->WireID().TPC,0);
695  fprimaryShower_hit_Z[idx]= xyzWire.Z();
696  fprimaryShower_hit_ch[idx]= wires[0]->Channel(); //only one channel per wire at collection plane
697  std::vector<art::Ptr<recob::SpacePoint>> sp = spFromShowerHits.at(j);
698 
699  if(!sp.empty() ){
700  //fprimaryShower_hit_X[idx]= sp[0]->XYZ()[0];
701  fprimaryShower_hit_Y[idx]= sp[0]->XYZ()[1];
702  //fprimaryShower_hit_Z[idx]= sp[0]->XYZ()[2];
703  }
704 
705  idx ++;
706  }
707  fprimaryShower_nHits = idx; //only collection hits
708  // Calorimetry only colleciton plane
709  std::vector<anab::Calorimetry> calovector = showerUtil.GetRecoShowerCalorimetry(*thisShower, evt, fShowerTag, fShowerCaloTag);
710  if(calovector.size() != 3 && fVerbose > 0)
711  std::cerr << "WARNING::Calorimetry vector size for primary is = " << calovector.size() << std::endl;
712 
713  for(size_t k = 0; k < calovector.size() && k<3; k++){
714  int plane = calovector[k].PlaneID().Plane;
715  if(plane !=2 ) continue;
716  fprimaryKineticEnergy[plane] = calovector[k].KineticEnergy();
717  fprimaryRange[plane] = calovector[k].Range();
718  fprimarynCal = calovector[k].dEdx().size();
719  for(size_t l=0; l<calovector[k].dEdx().size() && l<MAXHits; ++l){
720  fprimarydEdx[l]= calovector[k].dEdx()[l];
721  fprimarydQdx[l]= calovector[k].dQdx()[l];
722  fprimaryResidualRange[l]= calovector[k].ResidualRange()[l];
723  const auto &pos=(calovector[k].XYZ())[l];
724  fprimary_calX[l] = pos.X();
725  fprimary_calY[l] = pos.Y();
726  fprimary_calZ[l] = pos.Z();
727  fprimary_cal_pitch[l] =calovector[k].TrkPitchVec()[l];
728  }
729  }
730 
731  } // end is shower
732 
733 }
734 
735 // -----------------------------------------------------------------------------
736 void protoana::ProtoDUNEelectronAnaTree::FillPrimaryDaughterPFParticle(art::Event const & evt, const recob::PFParticle* daughterParticle, int daughterID){
737  art::Handle<std::vector<recob::Hit>> HitHandle;
738 
739  const recob::Track* daughterTrack = pfpUtil.GetPFParticleTrack(*daughterParticle,evt, fPFParticleTag,fTrackerTag);
740  const recob::Shower* daughterShower = pfpUtil.GetPFParticleShower(*daughterParticle,evt,fPFParticleTag,fShowerTag);
741  if(daughterTrack != 0x0){
742  fdaughterIstrack[fNDAUGHTERS] = 1;
743  fdaughterIsshower[fNDAUGHTERS] = 0;
744  fdaughterTheta[fNDAUGHTERS] = daughterTrack->Theta();
745  fdaughterPhi[fNDAUGHTERS] = daughterTrack->Phi();
746  fdaughterLength[fNDAUGHTERS] = daughterTrack->Length();
747  fdaughterStartPosition[fNDAUGHTERS][0] = daughterTrack->Trajectory().Start().X();
748  fdaughterStartPosition[fNDAUGHTERS][1] = daughterTrack->Trajectory().Start().Y();
749  fdaughterStartPosition[fNDAUGHTERS][2] = daughterTrack->Trajectory().Start().Z();
750  fdaughterEndPosition[fNDAUGHTERS][0] = daughterTrack->Trajectory().End().X();
751  fdaughterEndPosition[fNDAUGHTERS][1] = daughterTrack->Trajectory().End().Y();
752  fdaughterEndPosition[fNDAUGHTERS][2] = daughterTrack->Trajectory().End().Z();
753  fdaughterStartDirection[fNDAUGHTERS][0] = daughterTrack->Trajectory().StartDirection().X();
754  fdaughterStartDirection[fNDAUGHTERS][1] = daughterTrack->Trajectory().StartDirection().Y();
755  fdaughterStartDirection[fNDAUGHTERS][2] = daughterTrack->Trajectory().StartDirection().Z();
756  fdaughterEndDirection[fNDAUGHTERS][0] = daughterTrack->Trajectory().EndDirection().X();
757  fdaughterEndDirection[fNDAUGHTERS][1] = daughterTrack->Trajectory().EndDirection().Y();
758  fdaughterEndDirection[fNDAUGHTERS][2] = daughterTrack->Trajectory().EndDirection().Z();
759 
760  }
761  else if(daughterShower != 0x0){
762  fdaughterIstrack[fNDAUGHTERS] = 0;
763  fdaughterIsshower[fNDAUGHTERS] = 1;
764  fdaughterLength[fNDAUGHTERS] = daughterShower->Length();
765  fdaughterShowerBestPlane[fNDAUGHTERS] = daughterShower->best_plane();
766  fdaughterOpeningAngle[fNDAUGHTERS] = daughterShower->OpenAngle();
767  fdaughterStartPosition[fNDAUGHTERS][0] = daughterShower->ShowerStart().X();
768  fdaughterStartPosition[fNDAUGHTERS][1] = daughterShower->ShowerStart().Y();
769  fdaughterStartPosition[fNDAUGHTERS][2] = daughterShower->ShowerStart().Z();
770  fdaughterStartDirection[fNDAUGHTERS][0] = daughterShower->Direction().X();
771  fdaughterStartDirection[fNDAUGHTERS][1] = daughterShower->Direction().Y();
772  fdaughterStartDirection[fNDAUGHTERS][2] = daughterShower->Direction().Z();
773  }
774  fdaughterID[fNDAUGHTERS] = daughterID;
775  // NHits associated with this pfParticle
776  fdaughterNHits[fNDAUGHTERS] = (pfpUtil.GetPFParticleHits(*daughterParticle,evt,fPFParticleTag)).size();
777  // T0
778  std::vector<anab::T0> pfdaughterT0vec = pfpUtil.GetPFParticleT0(*daughterParticle,evt,fPFParticleTag);
779  if(!pfdaughterT0vec.empty())
780  fdaughterT0[fNDAUGHTERS] = pfdaughterT0vec[0].Time();
781 
782  // Increment counter
783  fNDAUGHTERS++;
784 
785 }
786 // -----------------------------------------------------------------------------
787 // -----------------------------------------------------------------------------
788 void protoana::ProtoDUNEelectronAnaTree::Initialize(){
789  fRun = -999;
790  fSubRun = -999;
791  fevent = -999;
792  fTimeStamp = -999.0;
793  for(int k=0; k < 5; k++)
794  fNactivefembs[k] = -999;
795 
796  for(int k=0; k < 3; k++){
797  fbeamtrackPos[k] = -999.0;
798  fbeamtrackEndPos[k] = -999.0;
799  fbeamtrackDir[k] = -999.0;
800  fprimaryTruth_vtx[k]= -999.0;
801  fprimaryVertex[k] = -999.0;
802  fdaughterVertex[k] = -999.0;
803  fprimaryEndPosition[k] = -999.0;
804  fprimaryStartPosition[k] = -999.0;
805  fprimaryEndDirection[k] = -999.0;
806  fprimaryStartDirection[k] = -999.0;
807  fprimaryKineticEnergy[k] = -999.0;
808  fprimaryRange[k] = -999.0;
809  }
810  for( int l=0; l<1000; l ++) {
811  fbeamtrackPos_at[l][0] = -999.;
812  fbeamtrackPos_at[l][1] = -999.;
813  fbeamtrackPos_at[l][2] = -999.;
814  fbeamtrackMom_at[l][0] = -999.;
815  fbeamtrackMom_at[l][1] = -999.;
816  fbeamtrackMom_at[l][2] = -999.;
817  fbeamtrackMom_at[l][3] = -999.;
818  }
819  fprimaryShower_nHits =0;
820  fprimaryTruthShower_nHits =0;
821  for( int m=0; m<MAXHits; m ++){
822  fprimarydEdx[m]= -999.0;
823  fprimarydQdx[m]= -999.0;
824  fprimary_calX[m] =-999.0;
825  fprimary_calY[m] =-999.0;
826  fprimary_calZ[m] =-999.0;
827  fprimary_cal_pitch[m]=-999.0;
828  fprimaryResidualRange[m] = -999.0;
829 
830  fprimaryShower_hit_w[m] =-999.0;
831  fprimaryShower_hit_q[m] =-999.0;
832  fprimaryShower_hit_t[m] =-999.0;
833  fprimaryShower_hit_X[m] =-999.0;
834  fprimaryShower_hit_Y[m] =-999.0;
835  fprimaryShower_hit_Z[m] =-999.0;
836  fprimaryShower_hit_ch[m] =-999.0;
837  fprimaryTruthShower_hit_w[m] =-999.0;
838  fprimaryTruthShower_hit_q[m] =-999.0;
839  fprimaryTruthShower_hit_t[m] =-999.0;
840  fprimaryTruthShower_hit_X[m] =-999.0;
841  fprimaryTruthShower_hit_Y[m] =-999.0;
842  fprimaryTruthShower_hit_Z[m] =-999.0;
843  fprimaryShower_hit_cnn[m] = -999.0;
844  }
845  fprimaryTruth_trkID =-999;
846  fprimaryTruth_pdg = -999;
847  fprimaryTruth_E = -999;
848  fprimaryTruth_Edepo = -999;
849  fprimaryTruth_purity = -999;
850  fprimarynCal = 0;
851  fbeamCheckIsMatched = -999;
852  fbeamtrigger = -999;
853  ftof = -999.0;
854  for(int k=0; k < 2; k++){
855  fcerenkovPressure[k] = -999.0;
856  fcerenkovTime[k] = -999.0;
857  fcerenkovStatus[k] = -999;
858  }
859  fbeamtrackMomentum = -999.0;
860  fbeamtrackEnergy = 999.0;
861  fbeamtrackPdg = -999;
862  fbeamtrackTime = -999.0;
863  fbeamtrackID = -999;
864  fbeam_ntrjPoints =0;
865  fbeamtrackNDaughters= 0;
866  fprimaryIstrack = -999;
867  fprimaryIsshower = -999;
868 
869  fprimaryBDTScore = -999.0;
870  fprimaryNHits = 0;
871  fprimaryTheta = -999.0;
872  fprimaryPhi = -999.0;
873  fprimaryLength = -999.0;
874  fprimaryMomentum = -999.0;
875  fprimaryEndMomentum = -999.0;
876  fprimaryOpeningAngle = -999.0;
877  fprimaryShowerBestPlane = -999;
878  fprimaryShowerEnergy = -999.0;
879  fprimaryShowerCharge = -999.0;
880  fprimaryShowerTruth_Charge = -999.0;
881  fprimaryShowerMIPEnergy = -999.0;
882  fprimaryID = -999;
883  fprimaryMomentumByRangeProton = -999.0;
884  fprimaryMomentumByRangeMuon = -999.0;
885  fprimaryT0 = -999.0;
886 
887  fNDAUGHTERS = 0;
888 
889  for(int k=0; k < NMAXDAUGTHERS; k++){
890  fdaughterIstrack[k] = -999;
891  fdaughterIsshower[k] = -999;
892  fdaughterNHits[k] = -999;
893  fdaughterTheta[k] = -999.0;
894  fdaughterPhi[k] = -999.0;
895  fdaughterLength[k] = -999.0;
896  for(int l=0; l < 3; l++){
897  fdaughterEndPosition[k][l] = -999.0;
898  fdaughterStartPosition[k][l] = -999.0;
899  fdaughterEndDirection[k][l] = -999.0;
900  fdaughterStartDirection[k][l] = -999.0;
901 
902  }
903  fdaughterOpeningAngle[k] = -999.0;
904  fdaughterShowerBestPlane[k] = -999;
905  fdaughterID[k] = -999;
906  fdaughterT0[k] = -999;
907 
908  }
909 
910 }
911 
912 DEFINE_ART_MODULE(protoana::ProtoDUNEelectronAnaTree)
simb::MCParticle::E
double E(const int i=0) const
Definition: MCParticle.h:233
protoana::ProtoDUNEelectronAnaTree::ProtoDUNEelectronAnaTree
ProtoDUNEelectronAnaTree(fhicl::ParameterSet const &p)
Definition: ProtoDUNEelectronAnaTree_module.cc:244
recob::Shower::best_plane
int best_plane() const
Definition: Shower.h:200
geo::WireGeo
Geometry description of a TPC wireThe wire is a single straight segment on a wire plane...
Definition: WireGeo.h:65
recob::Shower::ShowerStart
const TVector3 & ShowerStart() const
Definition: Shower.h:192
MAXHits
const int MAXHits
Definition: ProtoDUNEelectronAnaTree_module.cc:73
protoana::ProtoDUNEShowerUtils::GetRecoShowerHits
const std::vector< const recob::Hit * > GetRecoShowerHits(const recob::Shower &shower, art::Event const &evt, const std::string showerModule) const
Get the hits from a given reco shower.
Definition: ProtoDUNEShowerUtils.cxx:17
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
protoana::ProtoDUNETruthUtils::GetDepEnergyMC
double GetDepEnergyMC(const art::Event &evt, geo::GeometryCore const *fGeom, int trackid, int whichview) const
Definition: ProtoDUNETruthUtils.cxx:964
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::ProtoDUNETruthUtils::GetGeantGoodParticle
const simb::MCParticle * GetGeantGoodParticle(const simb::MCTruth &genTruth, const art::Event &evt) const
Definition: ProtoDUNETruthUtils.cxx:860
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
art::DataViewImpl::getHandle
Handle< PROD > getHandle(SelectorBase const &) const
Definition: DataViewImpl.h:382
string
std::string string
Definition: nybbler.cc:12
recob::Track::Trajectory
const recob::TrackTrajectory & Trajectory() const
Access to the stored recob::TrackTrajectory.
Definition: Track.h:98
protoana::ProtoDUNETrackUtils::GetRecoTrackCalorimetry
std::vector< anab::Calorimetry > GetRecoTrackCalorimetry(const recob::Track &track, art::Event const &evt, const std::string trackModule, const std::string caloModule) const
Get the Calorimetry(s) from a given reco track.
Definition: ProtoDUNETrackUtils.cxx:72
protoana::ProtoDUNEShowerUtils::GetRecoShowerCalorimetry
std::vector< anab::Calorimetry > GetRecoShowerCalorimetry(const recob::Shower &shower, art::Event const &evt, const std::string showerModule, const std::string caloModule) const
Get shower calo info.
Definition: ProtoDUNEShowerUtils.cxx:130
art::Timestamp::timeHigh
constexpr std::uint32_t timeHigh() const
Definition: Timestamp.h:34
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
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
anab::FVectorReader::getIndex
int getIndex(const std::string &name) const
Index of column with given name, or -1 if name not found.
Definition: MVAReader.h:82
std
STL namespace.
art::EDAnalyzer::EDAnalyzer
EDAnalyzer(fhicl::ParameterSet const &pset)
Definition: EDAnalyzer.h:25
MCParticle.h
Particle class.
protoana::ProtoDUNEelectronAnaTree::FillPrimaryPFParticle
void FillPrimaryPFParticle(art::Event const &evt, detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const recob::PFParticle *particle)
Definition: ProtoDUNEelectronAnaTree_module.cc:530
l
static QStrList * l
Definition: config.cpp:1044
Geometry.h
art framework interface to geometry description
NMAXDAUGTHERS
const int NMAXDAUGTHERS
Definition: ProtoDUNEelectronAnaTree_module.cc:72
simb::MCParticle::TrackId
int TrackId() const
Definition: MCParticle.h:210
protoana::ProtoDUNETruthUtils::GetMCParticleHits
std::vector< const recob::Hit * > GetMCParticleHits(detinfo::DetectorClocksData const &clockData, const simb::MCParticle &mcpart, const art::Event &evt, std::string hitModule, bool use_eve=true) const
Definition: ProtoDUNETruthUtils.cxx:103
util::size
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:92
art::DataViewImpl::isRealData
bool isRealData() const
Definition: DataViewImpl.cc:113
protoana::ProtoDUNEelectronAnaTree::analyze
void analyze(art::Event const &evt) override
Definition: ProtoDUNEelectronAnaTree_module.cc:377
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
recob::Track::Phi
double Phi() const
Definition: Track.h:178
protoana::ProtoDUNEelectronAnaTree::pfpUtil
protoana::ProtoDUNEPFParticleUtils pfpUtil
Definition: ProtoDUNEelectronAnaTree_module.cc:99
recob::Track::Length
double Length(size_t p=0) const
Access to various track properties.
Definition: Track.h:167
art::DataViewImpl::time
Timestamp time() const
Definition: DataViewImpl.cc:106
DEFINE_ART_MODULE
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:67
geo::fhicl::WireID
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
Definition: geo_types_fhicl.h:407
protoana::ProtoDUNETruthUtils::GetPurity
double GetPurity(detinfo::DetectorClocksData const &clockData, const recob::PFParticle &pfpart, const art::Event &evt, std::string pfparticleModule) const
Definition: ProtoDUNETruthUtils.cxx:177
recob::Track::Theta
double Theta() const
Access to spherical or geographical angles at vertex or at any point.
Definition: Track.h:176
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
art::DataViewImpl::getValidHandle
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
Definition: DataViewImpl.h:441
test.p
p
Definition: test.py:223
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::ProtoDUNEBeamlineUtils::IsGoodBeamlineTrigger
bool IsGoodBeamlineTrigger(art::Event const &evt) const
Definition: ProtoDUNEBeamlineUtils.cxx:649
art::DataViewImpl::run
RunNumber_t run() const
Definition: DataViewImpl.cc:71
geo::GeometryCore
Description of geometry of one entire detector.
Definition: GeometryCore.h:1468
keras_to_tensorflow.int
int
Definition: keras_to_tensorflow.py:69
detinfo::DetectorPropertiesData::ConvertTicksToX
double ConvertTicksToX(double ticks, int p, int t, int c) const
Definition: DetectorPropertiesData.cc:115
cache_state.m
m
Definition: cache_state.py:415
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
simb::MCParticle::Vx
double Vx(const int i=0) const
Definition: MCParticle.h:221
Hit.h
Declaration of signal hit object.
recob::Track::StartMomentum
double StartMomentum() const
Definition: Track.h:143
recob::PFParticle
Hierarchical representation of particle flow.
Definition: PFParticle.h:44
detinfo::DetectorClocksData
Contains all timing reference information for the detector.
Definition: DetectorClocksData.h:283
protoana::ProtoDUNEDataUtils::GetNActiveFembsForAPA
int GetNActiveFembsForAPA(art::Event const &evt, int apa) const
Get number of active fembs in an APA.
Definition: ProtoDUNEDataUtils.cxx:49
recob::TrackTrajectory::End
Point_t const & End() const
Returns the position of the last valid point of the trajectory [cm].
Definition: TrackTrajectory.h:257
Track.h
Provides recob::Track data product.
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
Wire.h
Declaration of basic channel signal object.
geo::WireGeo::GetCenter
void GetCenter(double *xyz, double localz=0.0) const
Fills the world coordinate of a point on the wire.
Definition: WireGeo.cxx:73
protoana::ProtoDUNEelectronAnaTree::FillPrimaryDaughterPFParticle
void FillPrimaryDaughterPFParticle(art::Event const &evt, const recob::PFParticle *daughterParticle, int daughterID)
Definition: ProtoDUNEelectronAnaTree_module.cc:736
tca::evt
TCEvent evt
Definition: DataStructs.cxx:7
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
art::fill_ptr_vector
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:297
recob::TrackTrajectory::StartDirection
Vector_t StartDirection() const
Returns the direction of the trajectory at the first point.
Definition: TrackTrajectory.h:317
trkf::Plane
recob::tracking::Plane Plane
Definition: TrackState.h:17
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::ProtoDUNEelectronAnaTree::operator=
ProtoDUNEelectronAnaTree & operator=(ProtoDUNEelectronAnaTree const &)=delete
anab::MVAReader::getOutput
std::array< float, N > getOutput(size_t key) const
Get copy of the MVA output vector at index "key".
Definition: MVAReader.h:129
recob::Shower::ID
int ID() const
Definition: Shower.h:187
protoana::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
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
trkf::TrackMomentumCalculator::GetTrackMomentum
double GetTrackMomentum(double trkrange, int pdg) const
Definition: TrackMomentumCalculator.cxx:144
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
protoana::ProtoDUNEDataUtils::IsBeamTrigger
bool IsBeamTrigger(art::Event const &evt) const
Definition: ProtoDUNEDataUtils.cxx:27
geo::GeometryCore::WirePtr
WireGeo const * WirePtr(geo::WireID const &wireid) const
Returns the specified wire.
Definition: GeometryCore.h:3400
Generated by   doxygen 1.8.11