ProtoDUNEPizeroAnaTree_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: ProtoDUNEPizeroAnaTree
3 // File: ProtoDUNEPizeroAnaTree_module.cc
4 //
5 // Extract protoDUNE useful information, do a first pre-selection and
6 // save output to a flat tree
7 //
8 //
9 // Georgios Christodoulou - georgios.christodoulou at cern.ch
10 // modified by Aaron Higuera ahiguera@central.uh.edu
11 // and Milo Vermeulen milov@nikhef.nl
12 //
13 ////////////////////////////////////////////////////////////////////////
14 
15 #include "art/Framework/Core/EDAnalyzer.h"
16 #include "art/Framework/Core/ModuleMacros.h"
17 #include "art/Framework/Principal/Event.h"
18 #include "art/Framework/Principal/Handle.h"
19 #include "art/Framework/Principal/Run.h"
20 #include "art/Framework/Principal/SubRun.h"
21 #include "art/Framework/Services/Registry/ServiceHandle.h"
22 #include "lardataobj/RecoBase/PFParticleMetadata.h"
23 #include "art_root_io/TFileService.h"
24 #include "canvas/Utilities/InputTag.h"
25 #include "canvas/Persistency/Common/FindManyP.h"
26 #include "fhiclcpp/ParameterSet.h"
27 #include "messagefacility/MessageLogger/MessageLogger.h"
28 
29 #include "larreco/Calorimetry/CalorimetryAlg.h"
30 #include "larsim/MCCheater/BackTrackerService.h"
31 #include "larsim/MCCheater/ParticleInventoryService.h"
32 #include "lardataobj/RecoBase/Hit.h"
33 #include "lardataobj/RecoBase/Track.h"
34 #include "lardataobj/RecoBase/Shower.h"
35 #include "lardataobj/RecoBase/PFParticle.h"
36 #include "nusimdata/SimulationBase/MCParticle.h"
37 #include "nusimdata/SimulationBase/MCTruth.h"
38 #include "lardataobj/AnalysisBase/CosmicTag.h"
39 #include "lardataobj/AnalysisBase/T0.h"
40 #include "lardataobj/AnalysisBase/Calorimetry.h"
41 #include "larcore/Geometry/Geometry.h"
42 #include "dune/DuneObj/ProtoDUNEBeamEvent.h"
43 #include "lardata/Utilities/GeometryUtilities.h"
44 #include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
45 #include "dune/Calib/XYZCalib.h"
46 #include "dune/CalibServices/XYZCalibService.h"
47 #include "dune/CalibServices/XYZCalibServiceProtoDUNE.h"
48 #include "lardata/ArtDataHelper/MVAReader.h"
49 #include "larevt/SpaceCharge/SpaceCharge.h"
50 #include "larevt/SpaceChargeServices/SpaceChargeService.h"
51 
52 
53 #include "protoduneana/Utilities/ProtoDUNETrackUtils.h"
54 #include "protoduneana/Utilities/ProtoDUNETruthUtils.h"
55 #include "protoduneana/Utilities/ProtoDUNEPFParticleUtils.h"
56 #include "dune/Protodune/singlephase/DataUtils/ProtoDUNEDataUtils.h"
57 #include "protoduneana/Utilities/ProtoDUNEShowerUtils.h"
58 #include "protoduneana/Utilities/ProtoDUNEBeamlineUtils.h"
59 #include "PiZeroProcess.h"
60 
61 // ROOT includes
62 #include "TTree.h"
63 #include "TFile.h"
64 #include "TString.h"
65 #include "TTimeStamp.h"
66 
67 // C++ Includes
68 #include <stdio.h>
69 #include <stdlib.h>
70 #include <fstream>
71 #include <string>
72 #include <sstream>
73 #include <cmath>
74 #include <algorithm>
75 #include <iostream>
76 #include <vector>
77 
78 // Maximum number of beam particle daughters to save
79 const int NMAXDAUGHTERS = 25;
80 // // Maximum number of objects in APA3 to save
81 // const int NMAXOBJECTS = 100;
82 // Maximum number of primary MC pi0s to save
83 const int NMAXPIZEROS = 25;
84 // Maximum number of hits to save
85 const int NMAXHITS = 5000;
86 
87 namespace protoana {
88  class ProtoDUNEPizeroAnaTree;
89 }
90 
91 
92 class protoana::ProtoDUNEPizeroAnaTree : public art::EDAnalyzer {
93 public:
94 
95  explicit ProtoDUNEPizeroAnaTree(fhicl::ParameterSet const & p);
96 
97  ProtoDUNEPizeroAnaTree(ProtoDUNEPizeroAnaTree const &) = delete;
98  ProtoDUNEPizeroAnaTree(ProtoDUNEPizeroAnaTree &&) = delete;
99  ProtoDUNEPizeroAnaTree & operator = (ProtoDUNEPizeroAnaTree const &) = delete;
100  ProtoDUNEPizeroAnaTree & operator = (ProtoDUNEPizeroAnaTree &&) = delete;
101 
102  virtual void beginJob() override;
103  virtual void endJob() override;
104 
105  // Required functions.
106  void analyze(art::Event const & evt) override;
107 
108 private:
109 
110  // Helper utility functions
111  protoana::ProtoDUNEDataUtils dataUtil;
112  protoana::ProtoDUNEPFParticleUtils pfpUtil;
113  protoana::ProtoDUNETrackUtils trackUtil;
114  protoana::ProtoDUNETruthUtils truthUtil;
115  protoana::ProtoDUNEShowerUtils showerUtil;
116  protoana::ProtoDUNEBeamlineUtils beamlineUtil;
117  art::ServiceHandle<cheat::BackTrackerService> bt_serv;
118 
119 
120  // Initialise tree variables
121  void Initialise();
122  void ResetPi0Vars();
123  void FillPrimaryPFParticle(art::Event const & evt,
124  detinfo::DetectorClocksData const& clockData,
125  detinfo::DetectorPropertiesData const& detProp,
126  const recob::PFParticle* particle);
127  void FillAPA3Object(art::Event const & evt, art::Handle<std::vector<recob::PFParticle>> pfpHandle);
128  void setPiZeroInfo(art::Event const & evt,
129  detinfo::DetectorClocksData const& clockData,
130  detinfo::DetectorPropertiesData const& detProp,
131  const std::vector<const simb::MCParticle*>& pi0s);
132 
133  // fcl parameters
134  const art::InputTag fBeamModuleLabel;
135  std::string fSimulationTag;
136  std::string fCalorimetryTag;
137  calo::CalorimetryAlg fCalorimetryAlg;
138  std::string fParticleIDTag;
139  std::string fTrackTag;
140  std::string fShowerTag;
141  std::string fShowerCaloTag;
142  std::string fPFParticleTag;
143  std::string fGeneratorTag;
144  std::string fHitTag;
145  int fVerbose;
146 
147  TTree *fPandoraBeam;
148 
149  // Tree variables
150  int fRun;
151  int fSubRun;
152  int fevent;
153  double fTimeStamp;
154  int fNactivefembs[6];
155 
156  // Beam track
157  int fbeamtrigger;
158  int fbeamCheckIsMatched;
159  double ftof;
160  int fcerenkovStatus[2];
161  double fcerenkovTime[2];
162  double fcerenkovPressure[2];
163  double fbeamtrackMomentum;
164  double fbeamtrackEnergy;
165  double fbeamtrackPos[3];
166  double fbeamtrackEndPos[3];
167  double fbeamtrackDir[3];
168  double fbeamtrackTime;
169  int fbeam_ntrjPoints;
170  int fbeamtrackPdg;
171  int fbeamtrackID;
172  double fbeamtrackPos_at[1000][3];
173  double fbeamtrackMom_at[1000][4];
174  std::vector<std::string> fbeamtrackEndProcess;
175  int fbeamtrackNDaughters;
176 
177  // Reconstructed primary particle
178  double fprimaryVertex[3];
179  int fprimaryIstrack;
180  int fprimaryIsshower;
181  double fprimaryBDTScore;
182  double fprimaryCNNScore;
183  int fprimaryNHits;
184  double fprimaryTheta;
185  double fprimaryPhi;
186  double fprimaryLength;
187  double fprimaryMomentum;
188  double fprimaryEndMomentum;
189  double fprimaryEndPosition[3];
190  double fprimaryStartPosition[3];
191  double fprimaryEndDirection[3];
192  double fprimaryStartDirection[3];
193  double fprimaryOpeningAngle;
194  int fprimaryShowerBestPlane;
195  double fprimaryShowerEnergy;
196  double fprimaryShowerCharge;
197  double fprimaryShowerMIPEnergy;
198  double fprimaryMomentumByRangeProton;
199  int fprimaryIsBeamparticle;
200  int fprimaryTruth_trkID;
201  int fprimaryTruth_pdg;
202  double fprimaryTruth_E;
203  double fprimaryTruth_vtx[3];
204  double fprimaryKineticEnergy[3];
205  double fprimaryRange[3];
206  int fprimarynCal;
207  double fprimarydEdx[NMAXHITS];
208  double fprimarydQdx[NMAXHITS];
209  double fprimary_cal_pos[NMAXHITS][3];
210  double fprimary_cal_pitch[NMAXHITS];
211  double fprimaryResidualRange[NMAXHITS];
212  int fprimaryShower_nHits; //collection only
213  int fprimaryShower_hit_w[NMAXHITS];
214  double fprimaryShower_hit_q[NMAXHITS];
215  double fprimaryShower_hit_t[NMAXHITS];
216  double fprimaryShower_hit_pos[NMAXHITS][3];
217  int fprimaryID;
218  double fprimaryT0;
219  int fprimaryNDaughters;
220 
221  // Primary particle daughters
222  int fprimaryDaughterID[NMAXDAUGHTERS];
223  int fprimaryDaughterIstrack[NMAXDAUGHTERS];
224  int fprimaryDaughterIsshower[NMAXDAUGHTERS];
225  double fprimaryDaughterBDTScore[NMAXDAUGHTERS];
226  double fprimaryDaughterCNNScore[NMAXDAUGHTERS];
227  int fprimaryDaughterNCollHits[NMAXDAUGHTERS];
228  int fprimaryDaughterNHits[NMAXDAUGHTERS];
229  double fprimaryDaughterEnergy[NMAXDAUGHTERS];
230  double fprimaryDaughterEnergyFromHits[NMAXDAUGHTERS];
231  double fprimaryDaughterMomentum[NMAXDAUGHTERS];
232  double fprimaryDaughterEndMomentum[NMAXDAUGHTERS];
233  double fprimaryDaughterLength[NMAXDAUGHTERS];
234  double fprimaryDaughterStartPosition[NMAXDAUGHTERS][3];
235  double fprimaryDaughterEndPosition[NMAXDAUGHTERS][3];
236  double fprimaryDaughterStartDirection[NMAXDAUGHTERS][3];
237  double fprimaryDaughterEndDirection[NMAXDAUGHTERS][3];
238  int fprimaryDaughterParentPdg[NMAXDAUGHTERS];
239  int fprimaryDaughterGrandparentPdg[NMAXDAUGHTERS];
240  int fprimaryDaughterParentID[NMAXDAUGHTERS];
241  int fprimaryDaughterGrandparentID[NMAXDAUGHTERS];
242  double fprimaryDaughterParentE[NMAXDAUGHTERS];
243  double fprimaryDaughterGrandparentE[NMAXDAUGHTERS];
244  double fprimaryDaughterParentStart[NMAXDAUGHTERS][3];
245  double fprimaryDaughterGrandparentStart[NMAXDAUGHTERS][3];
246  double fprimaryDaughterParentEnd[NMAXDAUGHTERS][3];
247  double fprimaryDaughterGrandparentEnd[NMAXDAUGHTERS][3];
248  int fprimaryDaughterParentEndProcess[NMAXDAUGHTERS];
249  // double fprimaryDaughterHitdQdx[NMAXDAUGHTERS][NMAXHITS];
250  double fprimaryDaughterCaloHitdQdx[NMAXDAUGHTERS][NMAXHITS];
251  double fprimaryDaughterCaloHitPitch[NMAXDAUGHTERS][NMAXHITS];
252  double fprimaryDaughterOwnHitPitch[NMAXDAUGHTERS][NMAXHITS];
253  double fprimaryDaughterHitCharge[NMAXDAUGHTERS][NMAXHITS];
254  // double fprimaryDaughterHitdEdxArea[NMAXDAUGHTERS][NMAXHITS];
255  int fprimaryDaughterHitWire[NMAXDAUGHTERS][NMAXHITS];
256  double fprimaryDaughterHitTime[NMAXDAUGHTERS][NMAXHITS];
257  double fprimaryDaughterHitPosition[NMAXDAUGHTERS][NMAXHITS][3];
258  // double fprimaryDaughterCaloHitPosition[NMAXDAUGHTERS][NMAXHITS][3];
259 
260  // // Instead of primary daughters, look at all objects in APA3
261  // int fObjID[NMAXOBJECTS];
262  // int fObjIsTrack[NMAXOBJECTS];
263  // int fObjIsShower[NMAXOBJECTS];
264  // int fObjIsPrimaryDaughter[NMAXOBJECTS];
265  // double fObjBDTscore[NMAXOBJECTS];
266  // double fObjCNNscore[NMAXOBJECTS];
267  // double fObjMomentum[NMAXOBJECTS];
268  // double fObjEndMomentum[NMAXOBJECTS];
269  // double fObjLength[NMAXOBJECTS];
270  // double fObjStartPosition[NMAXOBJECTS][3];
271  // double fObjEndPosition[NMAXOBJECTS][3];
272  // double fObjStartDirection[NMAXOBJECTS][3];
273  // double fObjEndDirection[NMAXOBJECTS][3];
274  // int fObjMCParentPdg[NMAXOBJECTS];
275  // int fObjMCGrandparentPdg[NMAXOBJECTS];
276  // int fObjMCParentID[NMAXOBJECTS];
277  // int fObjMCGrandparentID[NMAXOBJECTS];
278  // // Hit information
279  // int fObjNumHits[NMAXOBJECTS];
280  // double fObjHitPosition[NMAXOBJECTS][NMAXHITS][3];
281  // double fObjHitTime[NMAXOBJECTS][NMAXHITS];
282  // int fObjHitWire[NMAXOBJECTS][NMAXHITS];
283  // double fObjHitCharge[NMAXOBJECTS][NMAXHITS];
284  // double fObjHitPitch[NMAXOBJECTS][NMAXHITS];
285  // double fObjHitdEdx[NMAXOBJECTS][NMAXHITS];
286  // double fObjHitdQdx[NMAXOBJECTS][NMAXHITS];
287 
288  // MC pi0 variables
289  int fMCPi0ID[NMAXPIZEROS];
290  double fMCPi0Energy[NMAXPIZEROS];
291  double fMCPi0StartPosition[NMAXPIZEROS][3];
292  double fMCPi0EndPosition[NMAXPIZEROS][3];
293  int fMCPhoton1ID[NMAXPIZEROS];
294  int fMCPhoton2ID[NMAXPIZEROS];
295  double fMCPhoton1Energy[NMAXPIZEROS];
296  double fMCPhoton2Energy[NMAXPIZEROS];
297  double fMCPhoton1StartPosition[NMAXPIZEROS][3];
298  double fMCPhoton2StartPosition[NMAXPIZEROS][3];
299  double fMCPhoton1EndPosition[NMAXPIZEROS][3];
300  double fMCPhoton2EndPosition[NMAXPIZEROS][3];
301  // Reco hits related to photons
302  int fMCPhoton1NumHits[NMAXPIZEROS];
303  int fMCPhoton2NumHits[NMAXPIZEROS];
304  int fMCPhoton1_hit_w[NMAXPIZEROS][NMAXHITS];
305  int fMCPhoton2_hit_w[NMAXPIZEROS][NMAXHITS];
306  double fMCPhoton1_hit_t[NMAXPIZEROS][NMAXHITS];
307  double fMCPhoton2_hit_t[NMAXPIZEROS][NMAXHITS];
308  double fMCPhoton1_hit_q[NMAXPIZEROS][NMAXHITS];
309  double fMCPhoton2_hit_q[NMAXPIZEROS][NMAXHITS];
310  double fMCPhoton1_hit_pos[NMAXPIZEROS][NMAXHITS][3];
311  double fMCPhoton2_hit_pos[NMAXPIZEROS][NMAXHITS][3];
312 
313  // Reco pi0 variables
314  // Showers
315  int fShower1ID[NMAXPIZEROS];
316  int fShower2ID[NMAXPIZEROS];
317  double fShower1StartPosition[NMAXPIZEROS][3];
318  double fShower2StartPosition[NMAXPIZEROS][3];
319  double fShower1Direction[NMAXPIZEROS][3];
320  double fShower2Direction[NMAXPIZEROS][3];
321  double fShower1Length[NMAXPIZEROS];
322  double fShower2Length[NMAXPIZEROS];
323  double fShower1Completeness[NMAXPIZEROS];
324  double fShower2Completeness[NMAXPIZEROS];
325  double fShower1Purity[NMAXPIZEROS];
326  double fShower2Purity[NMAXPIZEROS];
327  // Reco hits related to showers
328  int fShower1NumHits[NMAXPIZEROS];
329  double fShower1_cnn_sc[NMAXPIZEROS][NMAXHITS];
330  double fShower2_cnn_sc[NMAXPIZEROS][NMAXHITS];
331  int fShower2NumHits[NMAXPIZEROS];
332  double fShower1_cal_pos[NMAXPIZEROS][NMAXHITS][3];
333  double fShower2_cal_pos[NMAXPIZEROS][NMAXHITS][3];
334  double fShower1_cal_E[NMAXPIZEROS];
335  double fShower2_cal_E[NMAXPIZEROS];
336  double fShower1Energy[NMAXPIZEROS];
337  double fShower2Energy[NMAXPIZEROS];
338  double fShower1EnergyFromHits[NMAXPIZEROS];
339  double fShower2EnergyFromHits[NMAXPIZEROS];
340  int fShower1HasBeamParent[NMAXPIZEROS];
341  int fShower2HasBeamParent[NMAXPIZEROS];
342  int fShower1HasBeamGrandparent[NMAXPIZEROS];
343  int fShower2HasBeamGrandparent[NMAXPIZEROS];
344  double fShower1_cal_pitch[NMAXPIZEROS][NMAXHITS];
345  double fShower2_cal_pitch[NMAXPIZEROS][NMAXHITS];
346  double fShower1_cal_dEdx[NMAXPIZEROS][NMAXHITS];
347  double fShower2_cal_dEdx[NMAXPIZEROS][NMAXHITS];
348  double fShower1_cal_dQdx[NMAXPIZEROS][NMAXHITS];
349  double fShower2_cal_dQdx[NMAXPIZEROS][NMAXHITS];
350 
351  // Other
352  geo::GeometryCore const* fGeometryService;
353  spacecharge::SpaceChargeService::provider_type const* fSCEService;
354 };
355 
356 
357 protoana::ProtoDUNEPizeroAnaTree::ProtoDUNEPizeroAnaTree(fhicl::ParameterSet const & p)
358  :
359  EDAnalyzer(p),
360  dataUtil(p.get<fhicl::ParameterSet>("DataUtils")),
361  beamlineUtil(p.get<fhicl::ParameterSet>("BeamLineUtils")),
362  fBeamModuleLabel(p.get< art::InputTag >("BeamModuleLabel")),
363  fSimulationTag(p.get<std::string>("SimulationTag")),
364  fCalorimetryTag(p.get<std::string>("CalorimetryTag")),
365  fCalorimetryAlg(p.get<fhicl::ParameterSet>("CalorimetryAlg")),
366  fParticleIDTag(p.get<std::string>("ParticleIDTag")),
367  fTrackTag(p.get<std::string>("TrackTag")),
368  fShowerTag(p.get<std::string>("ShowerTag")),
369  fShowerCaloTag(p.get<std::string>("ShowerCalorimetryTag")),
370  fPFParticleTag(p.get<std::string>("PFParticleTag")),
371  fGeneratorTag(p.get<std::string>("GeneratorTag")),
372  fHitTag(p.get<std::string>("HitTag")),
373  fVerbose(p.get<int>("Verbose")) {
374  // Get a pointer to the geometry service provider.
375  fGeometryService = lar::providerFrom<geo::Geometry>();
376  fSCEService = lar::providerFrom<spacecharge::SpaceChargeService>();
377 }
378 
379 void protoana::ProtoDUNEPizeroAnaTree::beginJob(){
380 
381  art::ServiceHandle<art::TFileService> tfs;
382 
383  fPandoraBeam = tfs->make<TTree>("PandoraBeam", "Beam events reconstructed with Pandora");
384  fPandoraBeam->Branch("run", &fRun, "run/I");
385  fPandoraBeam->Branch("subrun", &fSubRun, "subrun/I");
386  fPandoraBeam->Branch("event", &fevent, "event/I");
387  fPandoraBeam->Branch("timestamp", &fTimeStamp, "timestamp/D");
388  fPandoraBeam->Branch("Nactivefembs", &fNactivefembs, "Nactivefembs[5]/I");
389 
390  fPandoraBeam->Branch("beamtrigger", &fbeamtrigger, "beamtrigger/I");
391  fPandoraBeam->Branch("beamCheckIsMatched", &fbeamCheckIsMatched, "beamCheckIsMatched/I");
392  fPandoraBeam->Branch("tof", &ftof, "tof/D");
393  fPandoraBeam->Branch("cerenkovStatus", &fcerenkovStatus, "cerenkovStatus[2]/I");
394  fPandoraBeam->Branch("cerenkovTime", &fcerenkovTime, "cerenkovTime[2]/D");
395  fPandoraBeam->Branch("cerenkovPressure", &fcerenkovPressure, "cerenkovPressure[2]/D");
396  fPandoraBeam->Branch("beamtrackMomentum", &fbeamtrackMomentum, "beamtrackMomentum/D");
397  fPandoraBeam->Branch("beamtrackEnergy", &fbeamtrackEnergy, "beamtrackEnergy/D");
398  fPandoraBeam->Branch("beamtrackPos", &fbeamtrackPos, "beamtrackPos[3]/D");
399  fPandoraBeam->Branch("beamtrackEndPos", &fbeamtrackEndPos, "beamtrackEndPos[3]/D");
400  fPandoraBeam->Branch("beamtrackDir", &fbeamtrackDir, "beamtrackDir[3]/D");
401  fPandoraBeam->Branch("beamtrackTime", &fbeamtrackTime, "beamtrackTime/D");
402  fPandoraBeam->Branch("beamtrackPdg", &fbeamtrackPdg, "beamtrackPdg/I");
403  fPandoraBeam->Branch("beamtrackID", &fbeamtrackID, "beamtrackID/I");
404  fPandoraBeam->Branch("beam_ntrjPoints", &fbeam_ntrjPoints, "beam_ntrjPoints/I");
405  fPandoraBeam->Branch("beamtrackNDaughters", &fbeamtrackNDaughters, "beamtrackNDaughters/I");
406  fPandoraBeam->Branch("beamtrackPos_at", &fbeamtrackPos_at, "beamtrackPos_at[beam_ntrjPoints][3]/D");
407  fPandoraBeam->Branch("beamtrackMom_at", &fbeamtrackMom_at, "beamtrackMom_at[beam_ntrjPoints][4]/D");
408  fPandoraBeam->Branch("beamtrackEndProcess", &fbeamtrackEndProcess);
409 
410  fPandoraBeam->Branch("primaryVertex", &fprimaryVertex, "primaryVertex[3]/D");
411  fPandoraBeam->Branch("primaryIsBeamparticle", &fprimaryIsBeamparticle, "primaryIsBeamparticle/I");
412  fPandoraBeam->Branch("primaryIstrack", &fprimaryIstrack, "primaryIstrack/I");
413  fPandoraBeam->Branch("primaryIsshower", &fprimaryIsshower, "primaryIsshower/I");
414  fPandoraBeam->Branch("primaryBDTScore", &fprimaryBDTScore, "primaryBDTScore/D");
415  fPandoraBeam->Branch("primaryCNNScore", &fprimaryCNNScore, "primaryCNNScore/D");
416  fPandoraBeam->Branch("primaryNHits", &fprimaryNHits, "primaryNHits/I");
417  fPandoraBeam->Branch("primaryTheta", &fprimaryTheta, "primaryTheta/D");
418  fPandoraBeam->Branch("primaryPhi", &fprimaryPhi, "primaryPhi/D");
419  fPandoraBeam->Branch("primaryLength", &fprimaryLength, "primaryLength/D");
420  fPandoraBeam->Branch("primaryMomentum", &fprimaryMomentum, "primaryMomentum/D");
421  fPandoraBeam->Branch("primaryEndMomentum", &fprimaryEndMomentum, "primaryEndMomentum/D");
422  fPandoraBeam->Branch("primaryEndPosition", &fprimaryEndPosition, "primaryEndPosition[3]/D");
423  fPandoraBeam->Branch("primaryStartPosition", &fprimaryStartPosition, "primaryStartPosition[3]/D");
424  fPandoraBeam->Branch("primaryEndDirection", &fprimaryEndDirection, "primaryEndDirection[3]/D");
425  fPandoraBeam->Branch("primaryStartDirection", &fprimaryStartDirection, "primaryStartDirection[3]/D");
426  fPandoraBeam->Branch("primaryOpeningAngle", &fprimaryOpeningAngle, "primaryOpeningAngle/D");
427  fPandoraBeam->Branch("primaryID", &fprimaryID, "primaryID/I");
428  fPandoraBeam->Branch("primaryTruth_E", &fprimaryTruth_E, "primaryTruth_E/D");
429  fPandoraBeam->Branch("primaryTruth_vtx", &fprimaryTruth_vtx, "primaryTruth_vtx[3]/D");
430  fPandoraBeam->Branch("primaryTruth_pdg", &fprimaryTruth_pdg, "primaryTruth_pdg/I");
431  fPandoraBeam->Branch("primaryTruth_trkID", &fprimaryTruth_trkID, "primaryTruth_trkID/I");
432  fPandoraBeam->Branch("primaryShowerBestPlane", &fprimaryShowerBestPlane, "primaryShowerBestPlane/I");
433  fPandoraBeam->Branch("primaryShowerCharge", &fprimaryShowerCharge, "primaryShowerCharge/D");
434  fPandoraBeam->Branch("primaryShowerEnergy", &fprimaryShowerEnergy, "primaryShowerEnergy/D");
435  fPandoraBeam->Branch("primaryShowerMIPEnergy", &fprimaryShowerMIPEnergy, "primaryShowerMIPEnergy/D");
436 
437  fPandoraBeam->Branch("primaryKineticEnergy", &fprimaryKineticEnergy, "primaryKineticEnergy[3]/D");
438  fPandoraBeam->Branch("primaryRange", &fprimaryRange, "primaryRange[3]/D");
439  fPandoraBeam->Branch("primarynCal", &fprimarynCal, "primarynCal/I");
440  fPandoraBeam->Branch("primarydQdx", &fprimarydQdx, "primarydQdx[primarynCal]/D");
441  fPandoraBeam->Branch("primary_cal_pos", &fprimary_cal_pos, "primary_cal_pos[primarynCal][3]/D");
442  fPandoraBeam->Branch("primary_cal_pitch", &fprimary_cal_pitch, "primary_cal_pitch[primarynCal]/D");
443  fPandoraBeam->Branch("primarydEdx", &fprimarydEdx, "primarydEdx[primarynCal]/D");
444  fPandoraBeam->Branch("primaryResidualRange", &fprimaryResidualRange, "primaryResidualRange[primarynCal]/D");
445  fPandoraBeam->Branch("primaryT0", &fprimaryT0, "primaryT0/D");
446  fPandoraBeam->Branch("primaryNDaughters", &fprimaryNDaughters, "primaryNDaughters/I");
447 
448  fPandoraBeam->Branch("primaryDaughterID", &fprimaryDaughterID, "primaryDaughterID[primaryNDaughters]/I");
449  fPandoraBeam->Branch("primaryDaughterIstrack", &fprimaryDaughterIstrack, "primaryDaughterIstrack[primaryNDaughters]/I");
450  fPandoraBeam->Branch("primaryDaughterIsshower", &fprimaryDaughterIsshower, "primaryDaughterIsshower[primaryNDaughters]/I");
451  fPandoraBeam->Branch("primaryDaughterBDTScore", &fprimaryDaughterBDTScore, "primaryDaughterBDTScore[primaryNDaughters]/D");
452  fPandoraBeam->Branch("primaryDaughterCNNScore", &fprimaryDaughterCNNScore, "primaryDaughterCNNScore[primaryNDaughters]/D");
453  fPandoraBeam->Branch("primaryDaughterNCollHits", &fprimaryDaughterNCollHits, "primaryDaughterNCollHits[primaryNDaughters]/I");
454  fPandoraBeam->Branch("primaryDaughterNHits", &fprimaryDaughterNHits, "primaryDaughterNHits[primaryNDaughters]/I");
455  fPandoraBeam->Branch("primaryDaughterEnergy", &fprimaryDaughterEnergy, "primaryDaughterEnergy[primaryNDaughters]/D");
456  fPandoraBeam->Branch("primaryDaughterEnergyFromHits", &fprimaryDaughterEnergyFromHits, "primaryDaughterEnergyFromHits[primaryNDaughters]/D");
457  fPandoraBeam->Branch("primaryDaughterMomentum", &fprimaryDaughterMomentum, "primaryDaughterMomentum[primaryNDaughters]/D");
458  fPandoraBeam->Branch("primaryDaughterEndMomentum", &fprimaryDaughterEndMomentum, "primaryDaughterEndMomentum[primaryNDaughters]/D");
459  fPandoraBeam->Branch("primaryDaughterLength", &fprimaryDaughterLength, "primaryDaughterLength[primaryNDaughters]/D");
460  fPandoraBeam->Branch("primaryDaughterStartPosition", &fprimaryDaughterStartPosition, "primaryDaughterStartPosition[primaryNDaughters][3]/D");
461  fPandoraBeam->Branch("primaryDaughterEndPosition", &fprimaryDaughterEndPosition, "primaryDaughterEndPosition[primaryNDaughters][3]/D");
462  fPandoraBeam->Branch("primaryDaughterStartDirection", &fprimaryDaughterStartDirection, "primaryDaughterStartDirection[primaryNDaughters][3]/D");
463  fPandoraBeam->Branch("primaryDaughterEndDirection", &fprimaryDaughterEndDirection, "primaryDaughterEndDirection[primaryNDaughters][3]/D");
464  fPandoraBeam->Branch("primaryDaughterParentPdg", &fprimaryDaughterParentPdg, "primaryDaughterParentPdg[primaryNDaughters]/I");
465  fPandoraBeam->Branch("primaryDaughterGrandparentPdg", &fprimaryDaughterGrandparentPdg, "primaryDaughterGrandparentPdg[primaryNDaughters]/I");
466  fPandoraBeam->Branch("primaryDaughterParentID", &fprimaryDaughterParentID, "primaryDaughterParentID[primaryNDaughters]/I");
467  fPandoraBeam->Branch("primaryDaughterGrandparentID", &fprimaryDaughterGrandparentID, "primaryDaughterGrandparentID[primaryNDaughters]/I");
468  fPandoraBeam->Branch("primaryDaughterParentE", &fprimaryDaughterParentE, "primaryDaughterParentE[primaryNDaughters]/D");
469  fPandoraBeam->Branch("primaryDaughterGrandparentE", &fprimaryDaughterGrandparentE, "primaryDaughterGrandparentE[primaryNDaughters]/D");
470  fPandoraBeam->Branch("primaryDaughterParentStart", &fprimaryDaughterParentStart, "primaryDaughterParentStart[primaryNDaughters][3]/D");
471  fPandoraBeam->Branch("primaryDaughterGrandparentStart", &fprimaryDaughterGrandparentStart, "primaryDaughterGrandparentStart[primaryNDaughters][3]/D");
472  fPandoraBeam->Branch("primaryDaughterParentEnd", &fprimaryDaughterParentEnd, "primaryDaughterParentEnd[primaryNDaughters][3]/D");
473  fPandoraBeam->Branch("primaryDaughterGrandparentEnd", &fprimaryDaughterGrandparentEnd, "primaryDaughterGrandparentEnd[primaryNDaughters][3]/D");
474  fPandoraBeam->Branch("primaryDaughterParentEndProcess", &fprimaryDaughterParentEndProcess, "primaryDaughterParentEndProcess[primaryNDaughters]/I");
475  fPandoraBeam->Branch("primaryDaughterHitPosition", &fprimaryDaughterHitPosition, ("primaryDaughterHitPosition[primaryNDaughters][" + std::to_string(NMAXHITS) + "][3]/D").c_str());
476  // fPandoraBeam->Branch("primaryDaughterCaloHitPosition", &fprimaryDaughterCaloHitPosition, ("primaryDaughterCaloHitPosition[primaryNDaughters][" + std::to_string(NMAXHITS) + "][3]/D").c_str());
477  fPandoraBeam->Branch("primaryDaughterHitCharge", &fprimaryDaughterHitCharge, ("primaryDaughterHitCharge[primaryNDaughters][" + std::to_string(NMAXHITS) + "]/D").c_str());
478  // fPandoraBeam->Branch("primaryDaughterHitdEdxArea", &fprimaryDaughterHitdEdxArea, ("primaryDaughterHitdEdxArea[primaryNDaughters][" + std::to_string(NMAXHITS) + "]/D").c_str());
479  // fPandoraBeam->Branch("primaryDaughterHitdQdx", &fprimaryDaughterHitdQdx, ("primaryDaughterHitdQdx[primaryNDaughters][" + std::to_string(NMAXHITS) + "]/D").c_str());
480  fPandoraBeam->Branch("primaryDaughterCaloHitdQdx", &fprimaryDaughterCaloHitdQdx, ("primaryDaughterCaloHitdQdx[primaryNDaughters][" + std::to_string(NMAXHITS) + "]/D").c_str());
481  fPandoraBeam->Branch("primaryDaughterCaloHitPitch", &fprimaryDaughterCaloHitPitch, ("primaryDaughterCaloHitPitch[primaryNDaughters][" + std::to_string(NMAXHITS) + "]/D").c_str());
482  fPandoraBeam->Branch("primaryDaughterOwnHitPitch", &fprimaryDaughterOwnHitPitch, ("primaryDaughterOwnHitPitch[primaryNDaughters][" + std::to_string(NMAXHITS) + "]/D").c_str());
483  fPandoraBeam->Branch("primaryDaughterHitWire", &fprimaryDaughterHitWire, ("primaryDaughterHitWire[primaryNDaughters][" + std::to_string(NMAXHITS) + "]/I").c_str());
484  fPandoraBeam->Branch("primaryDaughterHitTime", &fprimaryDaughterHitTime, ("primaryDaughterHitTime[primaryNDaughters][" + std::to_string(NMAXHITS) + "]/D").c_str());
485 
486  // fPandoraBeam->Branch("ObjID", &fObjID, ("ObjID[" + std::to_string(NMAXOBJECTS) + "]/I").c_str());
487  // fPandoraBeam->Branch("ObjIsTrack", &fObjIsTrack, ("ObjIsTrack[" + std::to_string(NMAXOBJECTS) + "]/I").c_str());
488  // fPandoraBeam->Branch("ObjIsShower", &fObjIsShower, ("ObjIsShower[" + std::to_string(NMAXOBJECTS) + "]/I").c_str());
489  // fPandoraBeam->Branch("ObjIsPrimaryDaughter", &fObjIsPrimaryDaughter, ("ObjIsPrimaryDaughter[" + std::to_string(NMAXOBJECTS) + "]/I").c_str());
490  // fPandoraBeam->Branch("ObjBDTscore", &fObjBDTscore, ("ObjBDTscore[" + std::to_string(NMAXOBJECTS) + "]/D").c_str());
491  // fPandoraBeam->Branch("ObjCNNscore", &fObjCNNscore, ("ObjCNNscore[" + std::to_string(NMAXOBJECTS) + "]/D").c_str());
492  // fPandoraBeam->Branch("ObjMomentum", &fObjMomentum, ("ObjMomentum[" + std::to_string(NMAXOBJECTS) + "]/D").c_str());
493  // fPandoraBeam->Branch("ObjEndMomentum", &fObjEndMomentum, ("ObjEndMomentum[" + std::to_string(NMAXOBJECTS) + "]/D").c_str());
494  // fPandoraBeam->Branch("ObjLength", &fObjLength, ("ObjLength[" + std::to_string(NMAXOBJECTS) + "]/D").c_str());
495  // fPandoraBeam->Branch("ObjStartPosition", &fObjStartPosition, ("ObjStartPosition[" + std::to_string(NMAXOBJECTS) + "][3]/D").c_str());
496  // fPandoraBeam->Branch("ObjEndPosition", &fObjEndPosition, ("ObjEndPosition[" + std::to_string(NMAXOBJECTS) + "][3]/D").c_str());
497  // fPandoraBeam->Branch("ObjStartDirection", &fObjStartDirection, ("ObjStartDirection[" + std::to_string(NMAXOBJECTS) + "][3]/D").c_str());
498  // fPandoraBeam->Branch("ObjEndDirection", &fObjEndDirection, ("ObjEndDirection[" + std::to_string(NMAXOBJECTS) + "][3]/D").c_str());
499  // fPandoraBeam->Branch("ObjParentPdg", &fObjMCParentPdg, ("ObjParentPdg[" + std::to_string(NMAXOBJECTS) + "]/I").c_str());
500  // fPandoraBeam->Branch("ObjGrandparentPdg", &fObjMCGrandparentPdg, ("ObjGrandparentPdg[" + std::to_string(NMAXOBJECTS) + "]/I").c_str());
501  // fPandoraBeam->Branch("ObjParentID", &fObjMCParentID, ("ObjParentID[" + std::to_string(NMAXOBJECTS) + "]/I").c_str());
502  // fPandoraBeam->Branch("ObjGrandparentID", &fObjMCGrandparentID, ("ObjGrandparentID[" + std::to_string(NMAXOBJECTS) + "]/I").c_str());
503  // // Object hit information.
504  // fPandoraBeam->Branch("ObjNumHits", &fObjNumHits, ("ObjNumHits[" + std::to_string(NMAXOBJECTS) + "]/I").c_str());
505  // fPandoraBeam->Branch("ObjHitPosition", &fObjHitPosition, ("ObjHitPosition[" + std::to_string(NMAXOBJECTS) + "][" + std::to_string(NMAXHITS) + "][3]/D").c_str());
506  // fPandoraBeam->Branch("ObjHitTime", &fObjHitTime, ("ObjHitTime[" + std::to_string(NMAXOBJECTS) + "][" + std::to_string(NMAXHITS) + "]/D").c_str());
507  // fPandoraBeam->Branch("ObjHitWire", &fObjHitWire, ("ObjHitWire[" + std::to_string(NMAXOBJECTS) + "][" + std::to_string(NMAXHITS) + "]/I").c_str());
508  // fPandoraBeam->Branch("ObjHitCharge", &fObjHitCharge, ("ObjHitCharge[" + std::to_string(NMAXOBJECTS) + "][" + std::to_string(NMAXHITS) + "]/D").c_str());
509  // fPandoraBeam->Branch("ObjHitPitch", &fObjHitPitch, ("ObjHitPitch[" + std::to_string(NMAXOBJECTS) + "][" + std::to_string(NMAXHITS) + "]/D").c_str());
510  // fPandoraBeam->Branch("ObjHitdEdx", &fObjHitdEdx, ("ObjHitdEdx[" + std::to_string(NMAXOBJECTS) + "][" + std::to_string(NMAXHITS) + "]/D").c_str());
511  // fPandoraBeam->Branch("ObjHitdQdx", &fObjHitdQdx, ("ObjHitdQdx[" + std::to_string(NMAXOBJECTS) + "][" + std::to_string(NMAXHITS) + "]/D").c_str());
512 
513  // MC pi0 variables
514  // fPandoraBeam->Branch("MCPi0ID", &fMCPi0ID, ("MCPi0ID[" + std::to_string(NMAXPIZEROS) + "]/I").c_str());
515  // fPandoraBeam->Branch("MCPi0Energy", &fMCPi0Energy, ("MCPi0Energy[" + std::to_string(NMAXPIZEROS) + "]/D").c_str());
516  // fPandoraBeam->Branch("MCPi0StartPosition", &fMCPi0StartPosition, ("MCPi0StartPosition[" + std::to_string(NMAXPIZEROS) + "][3]/D").c_str());
517  // fPandoraBeam->Branch("MCPi0EndPosition", &fMCPi0EndPosition, ("MCPi0EndPosition[" + std::to_string(NMAXPIZEROS) + "][3]/D").c_str());
518  // fPandoraBeam->Branch("MCPhoton1ID", &fMCPhoton1ID, ("MCPhoton1ID[" + std::to_string(NMAXPIZEROS) + "]/I").c_str());
519  // fPandoraBeam->Branch("MCPhoton2ID", &fMCPhoton2ID, ("MCPhoton2ID[" + std::to_string(NMAXPIZEROS) + "]/I").c_str());
520  // fPandoraBeam->Branch("MCPhoton1Energy", &fMCPhoton1Energy, ("MCPhoton1Energy[" + std::to_string(NMAXPIZEROS) + "]/D").c_str());
521  // fPandoraBeam->Branch("MCPhoton2Energy", &fMCPhoton2Energy, ("MCPhoton2Energy[" + std::to_string(NMAXPIZEROS) + "]/D").c_str());
522  // fPandoraBeam->Branch("MCPhoton1StartPosition", &fMCPhoton1StartPosition, ("MCPhoton1StartPosition[" + std::to_string(NMAXPIZEROS) + "][3]/D").c_str());
523  // fPandoraBeam->Branch("MCPhoton2StartPosition", &fMCPhoton2StartPosition, ("MCPhoton2StartPosition[" + std::to_string(NMAXPIZEROS) + "][3]/D").c_str());
524  // fPandoraBeam->Branch("MCPhoton1EndPosition", &fMCPhoton1EndPosition, ("MCPhoton1EndPosition[" + std::to_string(NMAXPIZEROS) + "][3]/D").c_str());
525  // fPandoraBeam->Branch("MCPhoton2EndPosition", &fMCPhoton2EndPosition, ("MCPhoton2EndPosition[" + std::to_string(NMAXPIZEROS) + "][3]/D").c_str());
526  // // Reco hits related to photons
527  // fPandoraBeam->Branch("MCPhoton1NumHits", &fMCPhoton1NumHits, ("MCPhoton1NumHits[" + std::to_string(NMAXPIZEROS) + "]/I").c_str());
528  // fPandoraBeam->Branch("MCPhoton2NumHits", &fMCPhoton2NumHits, ("MCPhoton2NumHits[" + std::to_string(NMAXPIZEROS) + "]/I").c_str());
529  // fPandoraBeam->Branch("MCPhoton1_hit_w", &fMCPhoton1_hit_w, ("MCPhoton1_hit_w[" + std::to_string(NMAXPIZEROS) + "][" + std::to_string(NMAXHITS) + "]/I").c_str());
530  // fPandoraBeam->Branch("MCPhoton2_hit_w", &fMCPhoton2_hit_w, ("MCPhoton2_hit_w[" + std::to_string(NMAXPIZEROS) + "][" + std::to_string(NMAXHITS) + "]/I").c_str());
531  // fPandoraBeam->Branch("MCPhoton1_hit_q", &fMCPhoton1_hit_q, ("MCPhoton1_hit_q[" + std::to_string(NMAXPIZEROS) + "][" + std::to_string(NMAXHITS) + "]/D").c_str());
532  // fPandoraBeam->Branch("MCPhoton2_hit_q", &fMCPhoton2_hit_q, ("MCPhoton2_hit_q[" + std::to_string(NMAXPIZEROS) + "][" + std::to_string(NMAXHITS) + "]/D").c_str());
533  // fPandoraBeam->Branch("MCPhoton1_hit_t", &fMCPhoton1_hit_t, ("MCPhoton1_hit_t[" + std::to_string(NMAXPIZEROS) + "][" + std::to_string(NMAXHITS) + "]/D").c_str());
534  // fPandoraBeam->Branch("MCPhoton2_hit_t", &fMCPhoton2_hit_t, ("MCPhoton2_hit_t[" + std::to_string(NMAXPIZEROS) + "][" + std::to_string(NMAXHITS) + "]/D").c_str());
535  // fPandoraBeam->Branch("MCPhoton1_hit_pos", &fMCPhoton1_hit_pos, ("MCPhoton1_hit_pos[" + std::to_string(NMAXPIZEROS) + "][" + std::to_string(NMAXHITS) + "][3]/D").c_str());
536  // fPandoraBeam->Branch("MCPhoton2_hit_pos", &fMCPhoton2_hit_pos, ("MCPhoton2_hit_pos[" + std::to_string(NMAXPIZEROS) + "][" + std::to_string(NMAXHITS) + "][3]/D").c_str());
537 
538  // // Shower variables
539  // fPandoraBeam->Branch("Shower1ID", &fShower1ID, ("Shower1ID[" + std::to_string(NMAXPIZEROS) + "]/I").c_str());
540  // fPandoraBeam->Branch("Shower2ID", &fShower2ID, ("Shower2ID[" + std::to_string(NMAXPIZEROS) + "]/I").c_str());
541  // fPandoraBeam->Branch("Shower1StartPosition", &fShower1StartPosition, ("Shower1StartPosition[" + std::to_string(NMAXPIZEROS) + "][3]/D").c_str());
542  // fPandoraBeam->Branch("Shower2StartPosition", &fShower2StartPosition, ("Shower2StartPosition[" + std::to_string(NMAXPIZEROS) + "][3]/D").c_str());
543  // fPandoraBeam->Branch("Shower1Direction", &fShower1Direction, ("Shower1Direction[" + std::to_string(NMAXPIZEROS) + "][3]/D").c_str());
544  // fPandoraBeam->Branch("Shower2Direction", &fShower2Direction, ("Shower2Direction[" + std::to_string(NMAXPIZEROS) + "][3]/D").c_str());
545  // fPandoraBeam->Branch("Shower1Length", &fShower1Length, ("Shower1Length[" + std::to_string(NMAXPIZEROS) + "]/D").c_str());
546  // fPandoraBeam->Branch("Shower2Length", &fShower2Length, ("Shower2Length[" + std::to_string(NMAXPIZEROS) + "]/D").c_str());
547  // fPandoraBeam->Branch("Shower1Completeness", &fShower1Completeness, ("Shower1Completeness[" + std::to_string(NMAXPIZEROS) + "]/D").c_str());
548  // fPandoraBeam->Branch("Shower2Completeness", &fShower2Completeness, ("Shower2Completeness[" + std::to_string(NMAXPIZEROS) + "]/D").c_str());
549  // fPandoraBeam->Branch("Shower1Purity", &fShower1Purity, ("Shower1Purity[" + std::to_string(NMAXPIZEROS) + "]/D").c_str());
550  // fPandoraBeam->Branch("Shower2Purity", &fShower2Purity, ("Shower2Purity[" + std::to_string(NMAXPIZEROS) + "]/D").c_str());
551  // // Reco hits related to showers
552  // fPandoraBeam->Branch("Shower1NumHits", &fShower1NumHits, ("Shower1NumHits[" + std::to_string(NMAXPIZEROS) + "]/I").c_str());
553  // fPandoraBeam->Branch("Shower2NumHits", &fShower2NumHits, ("Shower2NumHits[" + std::to_string(NMAXPIZEROS) + "]/I").c_str());
554  // fPandoraBeam->Branch("Shower1_cnn_sc", &fShower1_cnn_sc, ("Shower1_cnn_sc[" + std::to_string(NMAXPIZEROS) + "][" + std::to_string(NMAXHITS) + "]/D").c_str());
555  // fPandoraBeam->Branch("Shower2_cnn_sc", &fShower2_cnn_sc, ("Shower2_cnn_sc[" + std::to_string(NMAXPIZEROS) + "][" + std::to_string(NMAXHITS) + "]/D").c_str());
556  // fPandoraBeam->Branch("Shower1_cal_pos", &fShower1_cal_pos, ("Shower1_cal_pos[" + std::to_string(NMAXPIZEROS) + "][" + std::to_string(NMAXHITS) + "][3]/D").c_str());
557  // fPandoraBeam->Branch("Shower2_cal_pos", &fShower2_cal_pos, ("Shower2_cal_pos[" + std::to_string(NMAXPIZEROS) + "][" + std::to_string(NMAXHITS) + "][3]/D").c_str());
558  // fPandoraBeam->Branch("Shower1_cal_E", &fShower1_cal_E, ("Shower1_cal_E[" + std::to_string(NMAXPIZEROS) + "]/D").c_str());
559  // fPandoraBeam->Branch("Shower2_cal_E", &fShower2_cal_E, ("Shower2_cal_E[" + std::to_string(NMAXPIZEROS) + "]/D").c_str());
560  // fPandoraBeam->Branch("Shower1Energy", &fShower1Energy, ("Shower1Energy[" + std::to_string(NMAXPIZEROS) + "]/D").c_str());
561  // fPandoraBeam->Branch("Shower2Energy", &fShower2Energy, ("Shower2Energy[" + std::to_string(NMAXPIZEROS) + "]/D").c_str());
562  // fPandoraBeam->Branch("Shower1EnergyFromHits", &fShower1EnergyFromHits, ("Shower1EnergyFromHits[" + std::to_string(NMAXPIZEROS) + "]/D").c_str());
563  // fPandoraBeam->Branch("Shower2EnergyFromHits", &fShower2EnergyFromHits, ("Shower2EnergyFromHits[" + std::to_string(NMAXPIZEROS) + "]/D").c_str());
564  // fPandoraBeam->Branch("Shower1HasBeamParent", &fShower1HasBeamParent, ("Shower1HasBeamParent[" + std::to_string(NMAXPIZEROS) + "]/I").c_str());
565  // fPandoraBeam->Branch("Shower2HasBeamParent", &fShower2HasBeamParent, ("Shower2HasBeamParent[" + std::to_string(NMAXPIZEROS) + "]/I").c_str());
566  // fPandoraBeam->Branch("Shower1HasBeamGrandparent", &fShower1HasBeamGrandparent, ("Shower1HasBeamGrandparent[" + std::to_string(NMAXPIZEROS) + "]/I").c_str());
567  // fPandoraBeam->Branch("Shower2HasBeamGrandparent", &fShower2HasBeamGrandparent, ("Shower2HasBeamGrandparent[" + std::to_string(NMAXPIZEROS) + "]/I").c_str());
568  // fPandoraBeam->Branch("Shower1_cal_pitch", &fShower1_cal_pitch, ("Shower1_cal_pitch[" + std::to_string(NMAXPIZEROS) + "][" + std::to_string(NMAXHITS) + "]/D").c_str());
569  // fPandoraBeam->Branch("Shower2_cal_pitch", &fShower2_cal_pitch, ("Shower2_cal_pitch[" + std::to_string(NMAXPIZEROS) + "][" + std::to_string(NMAXHITS) + "]/D").c_str());
570  // fPandoraBeam->Branch("Shower1_cal_dEdx", &fShower1_cal_dEdx, ("Shower1_cal_dEdx[" + std::to_string(NMAXPIZEROS) + "][" + std::to_string(NMAXHITS) + "]/D").c_str());
571  // fPandoraBeam->Branch("Shower2_cal_dEdx", &fShower2_cal_dEdx, ("Shower2_cal_dEdx[" + std::to_string(NMAXPIZEROS) + "][" + std::to_string(NMAXHITS) + "]/D").c_str());
572  // fPandoraBeam->Branch("Shower1_cal_dQdx", &fShower1_cal_dQdx, ("Shower1_cal_dQdx[" + std::to_string(NMAXPIZEROS) + "][" + std::to_string(NMAXHITS) + "]/D").c_str());
573  // fPandoraBeam->Branch("Shower2_cal_dQdx", &fShower2_cal_dQdx, ("Shower2_cal_dQdx[" + std::to_string(NMAXPIZEROS) + "][" + std::to_string(NMAXHITS) + "]/D").c_str());
574 }
575 
576 void protoana::ProtoDUNEPizeroAnaTree::analyze(art::Event const & evt){
577  // Initialise tree parameters
578  Initialise();
579  fRun = evt.run();
580  fSubRun = evt.subRun();
581  fevent = evt.id().event();
582  art::Timestamp ts = evt.time();
583  if (ts.timeHigh() == 0){
584  TTimeStamp ts2(ts.timeLow());
585  fTimeStamp = ts2.AsDouble();
586  }
587  else{
588  TTimeStamp ts2(ts.timeHigh(), ts.timeLow());
589  fTimeStamp = ts2.AsDouble();
590  }
591 
592  // Get number of active fembs
593  if(!evt.isRealData()){
594  for(int k=0; k < 6; k++)
595  fNactivefembs[k] = 20;
596  }
597  else{
598  for(int k=0; k < 6; k++)
599  fNactivefembs[k] = dataUtil.GetNActiveFembsForAPA(evt, k);
600  }
601 
602  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService>()->DataFor(evt);
603  auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService>()->DataFor(evt, clockData);
604 
605  //check for reco pandora stuff
606  auto recoParticleHandle = evt.getHandle<std::vector<recob::PFParticle>>(fPFParticleTag);
607  if( !recoParticleHandle ) return;
608 
609  // Get all of the PFParticles, by default from the "pandora" product
610  // auto recoParticles = evt.getValidHandle<std::vector<recob::PFParticle>>(fPFParticleTag);
611 
612  // We'd like to find the beam particle. Pandora tries to do this for us, so let's use the PFParticle utility
613  // to look for it. Pandora reconstructs slices containing one (or sometimes more) primary PFParticles. These
614  // are tagged as either beam or cosmic for ProtoDUNE. This function automatically considers only those
615  // PFParticles considered as primary
616  std::vector<const recob::PFParticle*> pfParticles = pfpUtil.GetPFParticlesFromBeamSlice(evt,fPFParticleTag);
617  for(const recob::PFParticle* particle : pfParticles){
618 
619  FillPrimaryPFParticle(evt, clockData, detProp, particle);
620  const TVector3 vtx = pfpUtil.GetPFParticleVertex(*particle,evt,fPFParticleTag,fTrackTag);
621  // std::cout << "Primary PFParticle ID: " << particle->Self() << '\n';
622 
623  fprimaryVertex[0] = vtx.X(); fprimaryVertex[1] = vtx.Y(); fprimaryVertex[2] = vtx.Z();
624  const TVector3 interactionVtx = pfpUtil.GetPFParticleSecondaryVertex(*particle,evt,fPFParticleTag,fTrackTag);
625 
626  // For now only consider the first primary track. Need a proper treatment if more than one primary particles are found.
627  break;
628  }
629 
630  // // Look at all the reconstructed particles that are fully contained in APA3.
631  // FillAPA3Object(evt, recoParticleHandle);
632 
633  bool beamTriggerEvent = false;
634  if(!evt.isRealData()){
635  // Firstly we need to get the list of MCTruth objects from the generator. The standard protoDUNE
636  auto mcTruths = evt.getValidHandle<std::vector<simb::MCTruth>>(fGeneratorTag);
637  // Define and fill a handle to point to a vector of the MCParticles
638  auto MCParticleHandle = evt.getHandle<std::vector<simb::MCParticle>>(fSimulationTag);
639  if (!MCParticleHandle) {
640  // Handle no simb::MCParticles.
641  throw cet::exception("ProtoDUNEPizeroAnaTree")
642  << " No simb::MCParticle objects in this event - "
643  << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
644  }
645  std::map<int, const simb::MCParticle*> MCPmap;
646  for(const simb::MCParticle& mcp : *MCParticleHandle) {
647  MCPmap[mcp.TrackId()] = &mcp;
648  }
649 
650  // Also get the reconstructed beam information in the MC - TO DO
651  const simb::MCParticle* geantGoodParticle = truthUtil.GetGeantGoodParticle((*mcTruths)[0],evt);
652  art::ServiceHandle<cheat::ParticleInventoryService> pi_serv;
653  if(geantGoodParticle != 0x0){
654  beamTriggerEvent = true;
655  fbeamtrigger = 12;
656  fbeamtrackPos[0] = geantGoodParticle->Vx();
657  fbeamtrackPos[1] = geantGoodParticle->Vy();
658  fbeamtrackPos[2] = geantGoodParticle->Vz();
659  fbeamtrackEndPos[0]= geantGoodParticle->EndX();
660  fbeamtrackEndPos[1]= geantGoodParticle->EndY();
661  fbeamtrackEndPos[2]= geantGoodParticle->EndZ();
662  fbeamtrackMomentum = geantGoodParticle->P();
663  fbeamtrackEnergy = geantGoodParticle->E();
664  fbeamtrackPdg = geantGoodParticle->PdgCode();
665  fbeamtrackTime = geantGoodParticle->T();
666  fbeamtrackID = geantGoodParticle->TrackId();
667  fbeamtrackEndProcess.push_back(geantGoodParticle->EndProcess());
668  fbeam_ntrjPoints = std::min(geantGoodParticle->NumberTrajectoryPoints(), (unsigned)1000);
669  fbeamtrackNDaughters = geantGoodParticle->NumberDaughters();
670  for( size_t i=0; i<geantGoodParticle->NumberTrajectoryPoints() && i<1000; ++i){
671  fbeamtrackPos_at[i][0] = geantGoodParticle->Position(i).X();
672  fbeamtrackPos_at[i][1] = geantGoodParticle->Position(i).Y();
673  fbeamtrackPos_at[i][2] = geantGoodParticle->Position(i).Z();
674  fbeamtrackMom_at[i][0] = geantGoodParticle->Momentum(i).Px();
675  fbeamtrackMom_at[i][1] = geantGoodParticle->Momentum(i).Py();
676  fbeamtrackMom_at[i][2] = geantGoodParticle->Momentum(i).Pz();
677  fbeamtrackMom_at[i][3] = geantGoodParticle->Momentum(i).E();
678  }
679 
680  std::cout << "Primary beam particle PDG code: " << fbeamtrackPdg << " and E: "
681  << fbeamtrackEnergy << '\n';
682 
683  std::vector<const simb::MCParticle*> pi0s;
684 
685  // Search for daughter pi0s if the primary particle is a pion.
686  if(abs(fbeamtrackPdg) == 211) {
687  std::cout << "Found a pion! (ID " << fbeamtrackID << ", PDG "
688  << fbeamtrackPdg << ")\n";
689 
690  std::cout << "Pion has pi0 daughters:\n";
691  std::vector<const simb::MCParticle*> daughters;
692  daughters.push_back(geantGoodParticle);
693  while(daughters.size()) {
694  const simb::MCParticle* parent = daughters[0];
695  for(int di = 0; di < parent->NumberDaughters(); ++di) {
696  const int daughter_ID = parent->Daughter(di);
697  const simb::MCParticle* daughter = MCPmap[daughter_ID];
698  daughters.push_back(daughter);
699  if(daughter->PdgCode() == 111) {
700  std::cout << " ID " << daughter_ID << ", PDG " << daughter->PdgCode()
701  << ", P " << daughter->P() << '\n';
702  pi0s.push_back(daughter);
703  }
704  }
705  daughters.erase(daughters.begin());
706  }
707  }
708 
709  // Record MC and shower information
710  setPiZeroInfo(evt, clockData, detProp, pi0s);
711 
712  } //geantGoodParticle
713  } // MC
714  else { //data
715  // For data we can see if this event comes from a beam trigger
716  beamTriggerEvent = dataUtil.IsBeamTrigger(evt);
717  if( !beamTriggerEvent ) return;
718 
719  std::vector< art::Ptr<beam::ProtoDUNEBeamEvent> > beaminfo;
720  auto pdbeamHandle = evt.getHandle< std::vector<beam::ProtoDUNEBeamEvent> >(fBeamModuleLabel);
721  if (pdbeamHandle)
722  art::fill_ptr_vector(beaminfo, pdbeamHandle);
723  else{
724  std::cout<<"No beam information from "<<fBeamModuleLabel<<std::endl;
725  }
726 
727  if(beaminfo.size()){
728  if( beamlineUtil.IsGoodBeamlineTrigger( evt ) ){
729  fbeamCheckIsMatched = beaminfo[0]->CheckIsMatched();
730  fbeamtrigger = beaminfo[0]->GetTimingTrigger();
731  // if TOFChan == -1, then there was not a successful match, if it's
732  // 0, 1, 2, or 3, then there was a good match corresponding to the
733  // different pair-wise combinations of the upstream and downstream
734  // channels
735  if(beaminfo[0]->GetTOFChan() != -1){
736  ftof = beaminfo[0]->GetTOF();
737  }
738  // If possible PIDs contain pion, set PDG code to pion.
739  std::vector<int> PIDs = beamlineUtil.GetPID(evt, 6);
740  if(std::find(PIDs.begin(), PIDs.end(), 211) != PIDs.end() && fprimaryIsshower != 1) {
741  fbeamtrackPdg = 211;
742  } else if(std::find(PIDs.begin(), PIDs.end(), 2212) != PIDs.end() && fprimaryIsshower != 1) {
743  fbeamtrackPdg = 2212;
744  } else if(std::find(PIDs.begin(), PIDs.end(), 321) != PIDs.end() && fprimaryIsshower != 1) {
745  fbeamtrackPdg = 321;
746  }
747  std::cout << "Set PDG to " << fbeamtrackPdg << '\n';
748  // Get Cerenkov
749  fcerenkovStatus[0] = beaminfo[0]->GetCKov0Status();
750  fcerenkovStatus[1] = beaminfo[0]->GetCKov1Status();
751  fcerenkovTime[0] = beaminfo[0]->GetCKov0Time();
752  fcerenkovTime[1] = beaminfo[0]->GetCKov1Time();
753  fcerenkovPressure[0] = beaminfo[0]->GetCKov0Pressure();
754  fcerenkovPressure[1] = beaminfo[0]->GetCKov1Pressure();
755  auto & tracks = beaminfo[0]->GetBeamTracks();
756  if(!tracks.empty()){
757  fbeamtrackPos[0] = tracks[0].Start().X();
758  fbeamtrackPos[1] = tracks[0].Start().Y();
759  fbeamtrackPos[2] = tracks[0].Start().Z();
760  fbeamtrackEndPos[0] = tracks[0].End().X();
761  fbeamtrackEndPos[1] = tracks[0].End().Y();
762  fbeamtrackEndPos[2] = tracks[0].End().Z();
763  fbeamtrackDir[0] = tracks[0].StartDirection().X();
764  fbeamtrackDir[1] = tracks[0].StartDirection().Y();
765  fbeamtrackDir[2] = tracks[0].StartDirection().Z();
766  }
767  auto & beammom = beaminfo[0]->GetRecoBeamMomenta();
768  if(!beammom.empty()) fbeamtrackMomentum = beammom[0];
769 
770  // Record pi0s reconstructed from showers
771  // Make recob::Showers accessible
772  auto rShowerHandle = evt.getHandle<std::vector<recob::Shower>>(fShowerTag);
773  if (!rShowerHandle) {
774  // Handle no recob::Showers.
775  throw cet::exception("PiZeroExtract")
776  << " No recob::Shower objects in this event - "
777  << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
778  }
779  // // Loop through showers
780  // for(const recob::Shower& shower : *rShowerHandle) {
781  // pizero::PiZeroProcess recopzproc(shower, evt, fShowerTag);
782  // // Only record if the showers come somewhat close
783  // if(!recopzproc.allRecoSet()) continue;
784  // const double shdist = pizero::ClosestDistance(recopzproc.shower1(), recopzproc.shower2());
785  // if(shdist > 0 && shdist < 100) {
786  // setPiZeroInfo(evt, recopzproc);
787  // break;
788  // }
789  // }
790  } //good beam trigger
791  } //good beaminfo
792  }//for data
793 
794  // Fill tree if event came from beam trigger.
795  if(beamTriggerEvent) fPandoraBeam->Fill();
796 
797  // Put some space between events.
798  std::cout << "\n\n";
799 }
800 
801 // -----------------------------------------------------------------------------
802 void protoana::ProtoDUNEPizeroAnaTree::endJob(){
803 
804 }
805 
806 // -----------------------------------------------------------------------------
807 void
808 protoana::ProtoDUNEPizeroAnaTree::FillPrimaryPFParticle(art::Event const & evt,
809  detinfo::DetectorClocksData const& clockData,
810  detinfo::DetectorPropertiesData const& detProp,
811  const recob::PFParticle* particle){
812 
813  // Pandora's BDT beam-cosmic score
814  fprimaryBDTScore = (double)pfpUtil.GetBeamCosmicScore(*particle,evt,fPFParticleTag);
815 
816  // "particle" is the pointer to our beam particle. The recob::Track or recob::Shower object
817  // of this particle might be more helpful. These return null pointers if not track-like / shower-like
818  const recob::Track* thisTrack = pfpUtil.GetPFParticleTrack(*particle, evt,fPFParticleTag,fTrackTag);
819  const recob::Shower* thisShower = pfpUtil.GetPFParticleShower(*particle,evt,fPFParticleTag,fShowerTag);
820  auto recoTracks = evt.getValidHandle<std::vector<recob::Track>>(fTrackTag);
821  auto recoShowers = evt.getValidHandle<std::vector<recob::Shower>>(fShowerTag);
822  art::FindManyP<recob::Hit> findHitsFromShowers(recoShowers,evt,fShowerTag);
823  art::FindManyP<recob::Hit> findHitsFromTracks(recoTracks,evt,fTrackTag);
824  // Set some common variables that have different access methods.
825  std::vector<art::Ptr<recob::Hit>> pfpHits;
826  std::vector<anab::Calorimetry> calovector;
827  const simb::MCParticle* mcparticle = NULL;
828  if(thisTrack != 0x0) {
829  pfpHits = findHitsFromTracks.at(thisTrack->ID());
830  calovector = trackUtil.GetRecoTrackCalorimetry(*thisTrack, evt, fTrackTag, fCalorimetryTag);
831  mcparticle = truthUtil.GetMCParticleFromRecoTrack(clockData, *thisTrack, evt, fTrackTag);
832  } else if(thisShower != 0x0) {
833  pfpHits = findHitsFromShowers.at(thisShower->ID());
834  calovector = showerUtil.GetRecoShowerCalorimetry(*thisShower, evt, fShowerTag, fShowerCaloTag);
835  mcparticle = truthUtil.GetMCParticleFromRecoShower(clockData, *thisShower, evt, fShowerTag);
836  }
837  // NHits associated with this pfParticle.
838  fprimaryNHits = pfpHits.size();
839 
840  if(mcparticle){
841  fprimaryTruth_pdg = mcparticle->PdgCode();
842  fprimaryTruth_trkID = mcparticle->TrackId();
843  fprimaryTruth_vtx[0] = mcparticle->Vx();
844  fprimaryTruth_vtx[1] = mcparticle->Vy();
845  fprimaryTruth_vtx[2] = mcparticle->Vx();
846  fprimaryTruth_E = mcparticle->E();
847  if( fbeamtrackID != -999 && fbeamtrackID == fprimaryTruth_trkID ) {
848  fprimaryIsBeamparticle = 1;
849  }
850  }
851 
852  // Calorimetry only collection plane
853  for(size_t k = 0; k < calovector.size() && k<3; k++){
854  int plane = calovector[k].PlaneID().Plane;
855  if(plane !=2 ) continue;
856  fprimaryKineticEnergy[plane] = calovector[k].KineticEnergy();
857  fprimaryRange[plane] = calovector[k].Range();
858  fprimarynCal = std::min(calovector[k].dEdx().size(), (size_t)NMAXHITS);
859  for(size_t l=0; l<calovector[k].dEdx().size() && l<NMAXHITS; ++l){
860  fprimarydEdx[l]= calovector[k].dEdx()[l];
861  fprimarydQdx[l]= calovector[k].dQdx()[l];
862  fprimaryResidualRange[l]= calovector[k].ResidualRange()[l];
863  const auto &pos=(calovector[k].XYZ())[l];
864  fprimary_cal_pos[l][0] = pos.X();
865  fprimary_cal_pos[l][1] = pos.Y();
866  fprimary_cal_pos[l][2] = pos.Z();
867  fprimary_cal_pitch[l] =calovector[k].TrkPitchVec()[l];
868  }
869  }
870 
871  // Get the CNN score for the primary particle.
872  anab::MVAReader<recob::Hit,4> hitResults(evt, "emtrkmichelid:emtrkmichel" );
873  fprimaryCNNScore = 0;
874  for(const art::Ptr<recob::Hit>& hit : pfpHits) {
875  std::array<float,4> cnn_out = hitResults.getOutput( hit );
876  const double p_trk_or_sh = cnn_out[ hitResults.getIndex("track") ]+ cnn_out[ hitResults.getIndex("em") ];
877  const double cnn_score = cnn_out[ hitResults.getIndex("em") ]/p_trk_or_sh;
878  fprimaryCNNScore += cnn_score;
879  }
880  fprimaryCNNScore /= fprimaryNHits;
881  std::cout << "Got primary CNN score " << fprimaryCNNScore << ", true PDG code: " << fprimaryTruth_pdg << '\n';
882  std::cout << "Pandora track: " << thisTrack << ", shower: " << thisShower << '\n';
883 
884  // Get the T0 for this pfParticle
885  std::vector<anab::T0> pfT0vec = pfpUtil.GetPFParticleT0(*particle,evt,fPFParticleTag);
886  if(!pfT0vec.empty())
887  fprimaryT0 = pfT0vec[0].Time();
888 
889  if(thisTrack != 0x0){
890 
891  fprimaryIstrack = 1;
892  fprimaryIsshower = 0;
893  fprimaryID = thisTrack->ParticleId();
894  fprimaryTheta = thisTrack->Theta();
895  fprimaryPhi = thisTrack->Phi();
896  fprimaryLength = thisTrack->Length();
897  fprimaryMomentum = thisTrack->StartMomentum();
898  fprimaryEndMomentum = thisTrack->EndMomentum();
899  fprimaryEndPosition[0] = thisTrack->Trajectory().End().X();
900  fprimaryEndPosition[1] = thisTrack->Trajectory().End().Y();
901  fprimaryEndPosition[2] = thisTrack->Trajectory().End().Z();
902  fprimaryStartPosition[0] = thisTrack->Trajectory().Start().X();
903  fprimaryStartPosition[1] = thisTrack->Trajectory().Start().Y();
904  fprimaryStartPosition[2] = thisTrack->Trajectory().Start().Z();
905  fprimaryEndDirection[0] = thisTrack->Trajectory().EndDirection().X();
906  fprimaryEndDirection[1] = thisTrack->Trajectory().EndDirection().Y();
907  fprimaryEndDirection[2] = thisTrack->Trajectory().EndDirection().Z();
908  fprimaryStartDirection[0] = thisTrack->Trajectory().StartDirection().X();
909  fprimaryStartDirection[1] = thisTrack->Trajectory().StartDirection().Y();
910  fprimaryStartDirection[2] = thisTrack->Trajectory().StartDirection().Z();
911 
912  } // end is track
913  else if(thisShower != 0x0){
914 
915  fprimaryIstrack = 0;
916  fprimaryIsshower = 1;
917  fprimaryID = thisShower->ID();
918  fprimaryLength = thisShower->Length();
919  fprimaryShowerBestPlane = thisShower->best_plane();
920  fprimaryOpeningAngle = thisShower->OpenAngle();
921  fprimaryStartPosition[0] = thisShower->ShowerStart().X();
922  fprimaryStartPosition[1] = thisShower->ShowerStart().Y();
923  fprimaryStartPosition[2] = thisShower->ShowerStart().Z();
924  fprimaryStartDirection[0] = thisShower->Direction().X();
925  fprimaryStartDirection[1] = thisShower->Direction().Y();
926  fprimaryStartDirection[2] = thisShower->Direction().Z();
927 
928  } // end is shower
929 
930  // Record children of primary track
931  auto recoParticleHandle = evt.getHandle<std::vector<recob::PFParticle>>(fPFParticleTag);
932  if( !recoParticleHandle ) return;
933 
934  fprimaryNDaughters = std::min(particle->NumDaughters(), NMAXDAUGHTERS);
935  for(int di = 0; di < fprimaryNDaughters; ++di) {
936  fprimaryDaughterID[di] = particle->Daughter(di);
937  const recob::PFParticle* daughter = &recoParticleHandle->at(fprimaryDaughterID[di]);
938  if(daughter == 0x0) continue;
939 
940  const recob::Track* daughterTrack = pfpUtil.GetPFParticleTrack(*daughter, evt,fPFParticleTag,fTrackTag);
941  const recob::Shower* daughterShower = pfpUtil.GetPFParticleShower(*daughter,evt,fPFParticleTag,fShowerTag);
942 
943  // Pandora's BDT beam-cosmic score
944  fprimaryDaughterBDTScore[di] = (double)pfpUtil.GetBeamCosmicScore(*daughter,evt,fPFParticleTag);
945  // Hits and calorimetry associated with this pfParticle
946  std::vector<art::Ptr<recob::Hit>> pfpDHits;
947  std::vector<anab::Calorimetry> dcalovector;
948  if(daughterTrack != 0x0) {
949  pfpDHits = findHitsFromTracks.at( daughterTrack->ID() );
950  dcalovector = trackUtil.GetRecoTrackCalorimetry(*daughterTrack, evt, fTrackTag, fCalorimetryTag);
951  } else if(daughterShower != 0x0) {
952  pfpDHits = findHitsFromShowers.at( daughterShower->ID() );
953  dcalovector = showerUtil.GetRecoShowerCalorimetry(*daughterShower, evt, fShowerTag, fShowerCaloTag);
954  }
955  if(dcalovector.size() != 3 && fVerbose > 0) {
956  std::cerr << "WARNING::Calorimetry vector size for track is = " << dcalovector.size() << std::endl;
957  }
958 
959  // Loop over calorimetry.
960  std::map<int, double> calo_hit_dQdx;
961  std::map<int, double> calo_hit_pitch;
962  for(size_t k = 0; k < dcalovector.size() && k<3; k++){
963  int plane = dcalovector[k].PlaneID().Plane;
964  if(plane !=2 ) continue;
965  for(size_t l=0; l<dcalovector[k].dEdx().size() && l<NMAXHITS; ++l){
966  // fprimaryDaughterHitdEdx[di][l]= dcalovector[k].dEdx()[l];
967  // fprimaryDaughterCaloHitdQdx[di][l]= dcalovector[k].dQdx()[l];
968  calo_hit_dQdx[dcalovector[k].TpIndices()[l]] = dcalovector[k].dQdx()[l];
969  calo_hit_pitch[dcalovector[k].TpIndices()[l]] = dcalovector[k].TrkPitchVec()[l];
970  // tmp_dQdx.push_back(fprimaryDaughterCaloHitdQdx[di][l]);
971  // fprimaryDaughterResidualRange[l]= dcalovector[k].ResidualRange()[l];
972  const geo::Point_t &pos=(dcalovector[k].XYZ())[l];
973  if(daughterShower != 0) {
974  const TVector3 hitvec(pos.X(), pos.Y(), pos.Z());
975  }
976  // fprimaryDaughterCaloHitPosition[di][l][0] = pos.X();
977  // fprimaryDaughterCaloHitPosition[di][l][1] = pos.Y();
978  // fprimaryDaughterCaloHitPosition[di][l][2] = pos.Z();
979  // fprimaryDaughterHitPitch[di][l] = dcalovector[k].TrkPitchVec()[l];
980  } // Loop over calorimetry hits.
981  } // Loop over calorimetry.
982 
983  // Loop over hits, record Aidan's CNN shower-like score.
984  art::FindManyP<recob::SpacePoint> spFromHits(pfpDHits, evt, fHitTag);
985  fprimaryDaughterCNNScore[di] = 0;
986  unsigned rec_hiti = 0;
987  for( unsigned i = 0; i < pfpDHits.size() && rec_hiti < NMAXHITS; ++i){
988  if(pfpDHits[i]->WireID().Plane != 2) continue;
989 
990  const geo::WireGeo* pwire = fGeometryService->WirePtr(pfpDHits[i]->WireID());
991  TVector3 xyzWire = pwire->GetCenter<TVector3>();
992  fprimaryDaughterHitWire[di][rec_hiti] = pfpDHits[i]->WireID().Wire;
993  fprimaryDaughterHitTime[di][rec_hiti] = pfpDHits[i]->PeakTime();
994  fprimaryDaughterHitPosition[di][rec_hiti][0] = detProp.ConvertTicksToX(pfpDHits[i]->PeakTime(),pfpDHits[i]->WireID().Plane,pfpDHits[i]->WireID().TPC,0);
995  fprimaryDaughterHitPosition[di][rec_hiti][2] = xyzWire.Z();
996  std::vector<art::Ptr<recob::SpacePoint>> sps = spFromHits.at(i);
997  if(!sps.empty()) {
998  fprimaryDaughterHitPosition[di][rec_hiti][1] = sps[0]->XYZ()[1];
999  }
1000  TVector3 hitpos(fprimaryDaughterHitPosition[di][rec_hiti][0],
1001  fprimaryDaughterHitPosition[di][rec_hiti][1],
1002  fprimaryDaughterHitPosition[di][rec_hiti][2]);
1003 
1004  std::array<float,4> cnn_out = hitResults.getOutput( pfpDHits[i] );
1005  const double p_trk_or_sh = cnn_out[ hitResults.getIndex("track") ]+ cnn_out[ hitResults.getIndex("em") ];
1006  const double cnn_score = cnn_out[ hitResults.getIndex("em") ]/p_trk_or_sh;
1007  fprimaryDaughterCNNScore[di] += cnn_score;
1008  // Record general hit information.
1009  fprimaryDaughterHitCharge[di][rec_hiti] = pfpDHits[i]->Integral();
1010  // Determine own hit pitch with SCE corrections.
1011  const double wirePitch = fGeometryService->WirePitch(pfpDHits[i]->WireID().planeID());
1012  TVector3 objdir(0,0,0);
1013  if(daughterTrack != 0x0) {
1014  objdir = {daughterTrack->Trajectory().StartDirection().X(),
1015  daughterTrack->Trajectory().StartDirection().Y(),
1016  daughterTrack->Trajectory().StartDirection().Z()};
1017  } else if(daughterShower != 0x0) {
1018  objdir = daughterShower->Direction();
1019  }
1020  TVector3 hitdir = objdir.Unit();
1021  if(objdir.Mag2() != 0 && hitpos[1] != -999.0) {
1022  // Here we have the object's direction and full hit position.
1023  // Set hit end to point on next wire in object's direction.
1024  objdir = objdir.Unit();
1025  TVector3 hitend = hitpos + objdir * (wirePitch/objdir.Z());
1026  if(fSCEService->EnableCalSpatialSCE()) {
1027  // Do SCE corrections on hitpos and hitend to get proper hit direction.
1028  // std::cout << "Doing SCE corrections.\n";
1029  const geo::Vector_t poscorr = fSCEService->GetCalPosOffsets(geo::Point_t(hitpos), pfpDHits[i]->WireID().TPC);
1030  const geo::Vector_t endcorr = fSCEService->GetCalPosOffsets(geo::Point_t(hitend), pfpDHits[i]->WireID().TPC);
1031  hitpos += TVector3(poscorr.X(), poscorr.Y(), poscorr.Z());
1032  hitend += TVector3(endcorr.X(), endcorr.Y(), endcorr.Z());
1033  }
1034  hitdir = (hitend-hitpos).Unit();
1035  } //else if(objdir.Mag2() != 0) {
1036  // std::cout << "Hit position incomplete: " << hitpos[0] << ',' << hitpos[1] << ',' << hitpos[2] << '\n';
1037  // } else if(hitpos[1] != -999) {
1038  // std::cout << "No object direction.\n";
1039  // } else { std::cout << "??????\n"; }
1040  fprimaryDaughterOwnHitPitch[di][rec_hiti] = wirePitch / abs(hitdir.Z());
1041  fprimaryDaughterCaloHitdQdx[di][rec_hiti] = calo_hit_dQdx[pfpDHits[i].key()];
1042  fprimaryDaughterCaloHitPitch[di][rec_hiti] = calo_hit_pitch[pfpDHits[i].key()];
1043  // std::cout << "Bare dx: " << wirePitch / abs(objdir.Z())
1044  // << ", own dx: " << fprimaryDaughterOwnHitPitch[di][rec_hiti]
1045  // << ", calo dx: " << fprimaryDaughterCaloHitPitch[di][rec_hiti] << '\n';
1046  ++rec_hiti;
1047  } // Loop over hits.
1048  fprimaryDaughterNCollHits[di] = rec_hiti;
1049  fprimaryDaughterNHits[di] = pfpDHits.size();
1050  fprimaryDaughterCNNScore[di] /= rec_hiti;
1051 
1052  if(daughterTrack != 0x0) {
1053  fprimaryDaughterIstrack[di] = 1;
1054  fprimaryDaughterIsshower[di] = 0;
1055 
1056  fprimaryDaughterMomentum[di] = daughterTrack->StartMomentum();
1057  fprimaryDaughterEndMomentum[di] = daughterTrack->EndMomentum();
1058  fprimaryDaughterLength[di] = daughterTrack->Length();
1059  fprimaryDaughterStartPosition[di][0] = daughterTrack->Trajectory().Start().X();
1060  fprimaryDaughterStartPosition[di][1] = daughterTrack->Trajectory().Start().Y();
1061  fprimaryDaughterStartPosition[di][2] = daughterTrack->Trajectory().Start().Z();
1062  fprimaryDaughterEndPosition[di][0] = daughterTrack->Trajectory().End().X();
1063  fprimaryDaughterEndPosition[di][1] = daughterTrack->Trajectory().End().Y();
1064  fprimaryDaughterEndPosition[di][2] = daughterTrack->Trajectory().End().Z();
1065  fprimaryDaughterStartDirection[di][0] = daughterTrack->Trajectory().StartDirection().X();
1066  fprimaryDaughterStartDirection[di][1] = daughterTrack->Trajectory().StartDirection().Y();
1067  fprimaryDaughterStartDirection[di][2] = daughterTrack->Trajectory().StartDirection().Z();
1068  fprimaryDaughterEndDirection[di][0] = daughterTrack->Trajectory().EndDirection().X();
1069  fprimaryDaughterEndDirection[di][1] = daughterTrack->Trajectory().EndDirection().Y();
1070  fprimaryDaughterEndDirection[di][2] = daughterTrack->Trajectory().EndDirection().Z();
1071  } else if(daughterShower != 0x0) {
1072  fprimaryDaughterIstrack[di] = 0;
1073  fprimaryDaughterIsshower[di] = 1;
1074 
1075  std::vector<const recob::Hit*> shHits = showerUtil.GetRecoShowerHits(*daughterShower, evt, fShowerTag);
1076  fprimaryDaughterEnergyFromHits[di] = showerUtil.EstimateEnergyFromHitCharge(clockData,
1077  detProp,
1078  shHits,
1079  fCalorimetryAlg)[2];
1080  fprimaryDaughterLength[di] = daughterShower->Length();
1081  fprimaryDaughterStartDirection[di][0] = daughterShower->Direction().X();
1082  fprimaryDaughterStartDirection[di][1] = daughterShower->Direction().Y();
1083  fprimaryDaughterStartDirection[di][2] = daughterShower->Direction().Z();
1084  fprimaryDaughterStartPosition[di][0] = daughterShower->ShowerStart().X();
1085  fprimaryDaughterStartPosition[di][1] = daughterShower->ShowerStart().Y();
1086  fprimaryDaughterStartPosition[di][2] = daughterShower->ShowerStart().Z();
1087  }
1088 
1089  // Attempt to get the (grand)parent truth information of this daughter.
1090  art::ServiceHandle<cheat::ParticleInventoryService> pi_serv;
1091  const simb::MCParticle* parent = truthUtil.GetMCParticleFromReco(clockData, *daughter, evt, fPFParticleTag);
1092  if(parent != 0x0) {
1093  fprimaryDaughterParentPdg[di] = parent->PdgCode();
1094  // std::cout << "True parent PDG: " << fprimaryDaughterParentPdg[di] << '\n';
1095  fprimaryDaughterParentID[di] = parent->TrackId();
1096  // std::cout << "Daughter parent: " << fprimaryDaughterParentID[di] << '\n';
1097  fprimaryDaughterParentE[di] = parent->E();
1098  // std::cout << "True parent E: " << fprimaryDaughterParentE[di] << '\n';
1099  if(parent->EndProcess() == "conv") {
1100  fprimaryDaughterParentEndProcess[di] = 0;
1101  } else if(parent->EndProcess() == "phot") {
1102  fprimaryDaughterParentEndProcess[di] = 1;
1103  }
1104  parent->Position().Vect().GetXYZ(fprimaryDaughterParentStart[di]);
1105  parent->EndPosition().Vect().GetXYZ(fprimaryDaughterParentEnd[di]);
1106  if(parent->Mother() != 0) {
1107  const simb::MCParticle* gparent = pi_serv->TrackIdToParticle_P(parent->Mother());
1108  fprimaryDaughterGrandparentPdg[di] = gparent->PdgCode();
1109  fprimaryDaughterGrandparentID[di] = gparent->TrackId();
1110  fprimaryDaughterGrandparentE[di] = gparent->E();
1111  // std::cout << "True grandparent PDG: " << fprimaryDaughterGrandparentPdg[di] << '\n';
1112  gparent->Position().Vect().GetXYZ(fprimaryDaughterGrandparentStart[di]);
1113  gparent->EndPosition().Vect().GetXYZ(fprimaryDaughterGrandparentEnd[di]);
1114  }
1115  } // If MCParticle parent.
1116  } // for PFParticle daughters
1117 
1118 } // FillPrimaryPFParticle
1119 
1120 // -----------------------------------------------------------------------------
1121 // -----------------------------------------------------------------------------
1122 // void protoana::ProtoDUNEPizeroAnaTree::FillAPA3Object(art::Event const & evt, art::Handle<std::vector<recob::PFParticle>> pfpHandle) {
1123 // // Put all PFParticles in a map by ID.
1124 // std::unordered_map<size_t, const recob::PFParticle*> PFPmap;
1125 // for(const recob::PFParticle& pfp : *pfpHandle) {
1126 // PFPmap[pfp.Self()] = &pfp;
1127 // }
1128 // // Primary particle ID.
1129 // const size_t beamID = pfpUtil.GetPFParticlesFromBeamSlice(evt,fPFParticleTag).size() > 0? pfpUtil.GetPFParticlesFromBeamSlice(evt,fPFParticleTag)[0]->Self(): -1;
1130 // // Keep track of the number of saved objects
1131 // unsigned obji = 0;
1132 // for(const recob::PFParticle& pfp : *pfpHandle) {
1133 // if(obji >= NMAXOBJECTS) break;
1134 //
1135 // // Determine whether particle is classified as track or shower.
1136 // const recob::Track* track = pfpUtil.GetPFParticleTrack(pfp, evt,fPFParticleTag,fTrackTag);
1137 // const recob::Shower* shower = pfpUtil.GetPFParticleShower(pfp,evt,fPFParticleTag,fShowerTag);
1138 //
1139 // // Do initial check whether object is fully contained in APA3.
1140 // if(track != 0x0) {
1141 // const int TPCstart = fGeometryService->FindTPCAtPosition(track->Start()).deepestIndex();
1142 // const int TPCend = fGeometryService->FindTPCAtPosition(track->End()).deepestIndex();
1143 // // APA3 == TPCID 1??
1144 // if(TPCstart != 1 || TPCend != 1) continue;
1145 // } else if(shower != 0x0) {
1146 // const int TPCstart = fGeometryService->FindTPCAtPosition(shower->ShowerStart()).deepestIndex();
1147 // const int TPCend = fGeometryService->FindTPCAtPosition(shower->ShowerStart() + shower->Length() * shower->Direction()).deepestIndex();
1148 // // APA3 == TPCID 1??
1149 // if(TPCstart != 1 || TPCend != 1) continue;
1150 // } else {
1151 // continue;
1152 // }
1153 //
1154 // // General PFParticle attributes
1155 // fObjID[obji] = pfp.Self();
1156 // if(PFPmap[pfp.Parent()] != 0x0 && PFPmap[pfp.Parent()]->Self() == beamID) {
1157 // fObjIsPrimaryDaughter[obji] = 1;
1158 // // std::cout << "Yes, ID " << fObjID[obji] << '\n';
1159 // } else {
1160 // // std::cout << "No, ID " << fObjID[obji] << '\n';
1161 // fObjIsPrimaryDaughter[obji] = 0;
1162 // }
1163 //
1164 // // Pandora's BDT beam-cosmic score
1165 // fObjBDTscore[obji] = (double)pfpUtil.GetBeamCosmicScore(pfp,evt,fPFParticleTag);
1166 //
1167 // // Get associations between PFParticles, clusters and hits.
1168 // const std::vector<const recob::Cluster*> pfpClusters = pfpUtil.GetPFParticleClusters(pfp,evt,fPFParticleTag);
1169 // const auto allClusters = evt.getValidHandle<std::vector<recob::Cluster>>(fPFParticleTag);
1170 // const art::FindManyP<recob::Hit> findHits(allClusters, evt, fPFParticleTag);
1171 // std::vector<art::Ptr<recob::Hit>> pfpHits;
1172 // for(auto cluster : pfpClusters){
1173 // const std::vector<art::Ptr<recob::Hit>> clusterHits = findHits.at(cluster->ID());
1174 // for(art::Ptr<recob::Hit> hit : clusterHits){
1175 // pfpHits.push_back(hit);
1176 // }
1177 // }
1178 // // NHits associated with this pfParticle
1179 // fObjNumHits[obji] = pfpHits.size();
1180 //
1181 // // Aidan's CNN track-like score: determine the average for this object and save.
1182 // anab::MVAReader<recob::Hit,4> hitResults(evt, "emtrkmichelid:emtrkmichel" );
1183 // fObjCNNscore[obji] = 0;
1184 // // Loop over PFP hits.
1185 // for( unsigned i = 0; i < pfpHits.size() && i < NMAXHITS; ++i){
1186 // // Collection plane only.
1187 // if(pfpHits[i]->WireID().Plane != 2) continue;
1188 // // Record CNN score per hit.
1189 // const std::array<float,4> cnn_out = hitResults.getOutput( pfpHits[i] );
1190 // const double p_trk_or_sh = cnn_out[ hitResults.getIndex("track") ]+ cnn_out[ hitResults.getIndex("em") ];
1191 // const double cnn_score = cnn_out[ hitResults.getIndex("em") ]/p_trk_or_sh;
1192 // fObjCNNscore[obji] += cnn_score;
1193 // // Record other hit information.
1194 // fObjHitTime[obji][i] = pfpHits[i]->PeakTime();
1195 // fObjHitWire[obji][i] = pfpHits[i]->WireID().Wire;
1196 // fObjHitCharge[obji][i] = pfpHits[i]->Integral();
1197 // }
1198 // fObjCNNscore[obji] /= fObjNumHits[obji];
1199 // // std::cout << "CNN score: " << fObjCNNscore[obji] << ". Pandora says:\n";
1200 //
1201 // // Distinctive track/shower information.
1202 // const simb::MCParticle* mcparticle;
1203 // if(track != 0) {
1204 // // PFParticle is marked as track by Pandora.
1205 // mcparticle = truthUtil.GetMCParticleFromReco(*track, evt, fTrackTag);
1206 //
1207 // fObjIsShower[obji] = 0;
1208 // fObjIsTrack[obji] = 1;
1209 // fObjLength[obji] = track->Length();
1210 // fObjStartPosition[obji][0] = track->Start().X();
1211 // fObjStartPosition[obji][1] = track->Start().Y();
1212 // fObjStartPosition[obji][2] = track->Start().Z();
1213 // fObjEndPosition[obji][0] = track->End().X();
1214 // fObjEndPosition[obji][1] = track->End().Y();
1215 // fObjEndPosition[obji][2] = track->End().Z();
1216 // fObjStartDirection[obji][0] = track->StartDirection().X();
1217 // fObjStartDirection[obji][1] = track->StartDirection().Y();
1218 // fObjStartDirection[obji][2] = track->StartDirection().Z();
1219 // fObjEndDirection[obji][0] = track->EndDirection().X();
1220 // fObjEndDirection[obji][1] = track->EndDirection().Y();
1221 // fObjEndDirection[obji][2] = track->EndDirection().Z();
1222 // fObjMomentum[obji] = track->StartMomentum();
1223 // fObjEndMomentum[obji] = track->EndMomentum();
1224 //
1225 // // Calorimetry only collection plane
1226 // std::vector<anab::Calorimetry> calovector = trackUtil.GetRecoTrackCalorimetry(*track, evt, fTrackTag, fCalorimetryTag);
1227 // if(calovector.size() != 3 && fVerbose > 0)
1228 // std::cerr << "WARNING::Calorimetry vector size for track is = " << calovector.size() << std::endl;
1229 //
1230 // for(size_t k = 0; k < calovector.size() && k<3; k++){
1231 // int plane = calovector[k].PlaneID().Plane;
1232 // if(plane !=2 ) continue;
1233 // for(size_t l=0; l<calovector[k].dEdx().size() && l<NMAXHITS; ++l){
1234 // fObjHitdEdx[obji][l]= calovector[k].dEdx()[l];
1235 // fObjHitdQdx[obji][l]= calovector[k].dQdx()[l];
1236 // // fObjResidualRange[l]= calovector[k].ResidualRange()[l];
1237 // const auto &pos=(calovector[k].XYZ())[l];
1238 // fObjHitPosition[obji][l][0] = pos.X();
1239 // fObjHitPosition[obji][l][1] = pos.Y();
1240 // fObjHitPosition[obji][l][2] = pos.Z();
1241 // fObjHitPitch[obji][l] = calovector[k].TrkPitchVec()[l];
1242 // }
1243 // }
1244 // } else if(shower != 0) {
1245 // // PFParticle is marked as shower by Pandora.
1246 // mcparticle = truthUtil.GetMCParticleFromReco(*shower, evt, fShowerTag);
1247 //
1248 // fObjIsTrack[obji] = 0;
1249 // fObjIsShower[obji] = 1;
1250 // fObjLength[obji] = shower->Length();
1251 // // fObjShowerBestPlane[obji] = shower->best_plane();
1252 // // fObjOpeningAngle[obji] = shower->OpenAngle();
1253 // fObjStartPosition[obji][0] = shower->ShowerStart().X();
1254 // fObjStartPosition[obji][1] = shower->ShowerStart().Y();
1255 // fObjStartPosition[obji][2] = shower->ShowerStart().Z();
1256 // const TVector3 shEnd = shower->ShowerStart() + shower->Length()*shower->Direction();
1257 // fObjEndPosition[obji][0] = shEnd.X();
1258 // fObjEndPosition[obji][1] = shEnd.Y();
1259 // fObjEndPosition[obji][2] = shEnd.Z();
1260 // fObjStartDirection[obji][0] = shower->Direction().X();
1261 // fObjStartDirection[obji][1] = shower->Direction().Y();
1262 // fObjStartDirection[obji][2] = shower->Direction().Z();
1263 // fObjEndDirection[obji][0] = shower->Direction().X();
1264 // fObjEndDirection[obji][1] = shower->Direction().Y();
1265 // fObjEndDirection[obji][2] = shower->Direction().Z();
1266 //
1267 // // Calorimetry only collection plane
1268 // std::vector<anab::Calorimetry> calovector = showerUtil.GetRecoShowerCalorimetry(*shower, evt, fShowerTag, fShowerCaloTag);
1269 // if(calovector.size() != 3 && fVerbose > 0)
1270 // std::cerr << "WARNING::Calorimetry vector size for shower is = " << calovector.size() << std::endl;
1271 //
1272 // for(size_t k = 0; k < calovector.size() && k<3; k++){
1273 // int plane = calovector[k].PlaneID().Plane;
1274 // if(plane !=2 ) continue;
1275 // for(size_t l=0; l<calovector[k].dEdx().size() && l<NMAXHITS; ++l){
1276 // fObjHitdEdx[obji][l]= calovector[k].dEdx()[l];
1277 // fObjHitdQdx[obji][l]= calovector[k].dQdx()[l];
1278 // // fObjResidualRange[l]= calovector[k].ResidualRange()[l];
1279 // const auto &pos=(calovector[k].XYZ())[l];
1280 // fObjHitPosition[obji][l][0] = pos.X();
1281 // fObjHitPosition[obji][l][1] = pos.Y();
1282 // fObjHitPosition[obji][l][2] = pos.Z();
1283 // fObjHitPitch[obji][l] = calovector[k].TrkPitchVec()[l];
1284 // }
1285 // }
1286 // } else {
1287 // continue;
1288 // }
1289 // // MCParticle matching
1290 // if(mcparticle){
1291 // art::ServiceHandle<cheat::ParticleInventoryService> pi_serv;
1292 // fObjMCParentPdg[obji] = mcparticle->PdgCode();
1293 // fObjMCParentID[obji] = mcparticle->TrackId();
1294 // if(mcparticle->Mother() > 0 && pi_serv->TrackIdToParticle_P(mcparticle->Mother()) != 0x0) {
1295 // fObjMCGrandparentPdg[obji] = pi_serv->TrackIdToParticle_P(mcparticle->Mother())->PdgCode();
1296 // fObjMCGrandparentID[obji] = pi_serv->TrackIdToParticle_P(mcparticle->Mother())->TrackId();
1297 // }
1298 // }
1299 //
1300 // // int fObjID[NMAXOBJECTS];
1301 // // int fObjIsTrack[NMAXOBJECTS];
1302 // // int fObjIsShower[NMAXOBJECTS];
1303 // // int fObjIsPrimaryDaughter[NMAXOBJECTS];
1304 // // double fObjBDTscore[NMAXOBJECTS];
1305 // // double fObjCNNscore[NMAXOBJECTS];
1306 // // double fObjMomentum[NMAXOBJECTS];
1307 // // double fObjEndMomentum[NMAXOBJECTS];
1308 // // double fObjLength[NMAXOBJECTS];
1309 // // double fObjStartPosition[NMAXOBJECTS][3];
1310 // // double fObjEndPosition[NMAXOBJECTS][3];
1311 // // double fObjStartDirection[NMAXOBJECTS][3];
1312 // // double fObjEndDirection[NMAXOBJECTS][3];
1313 // // int fObjMCParentPdg[NMAXOBJECTS];
1314 // // int fObjMCGrandparentPdg[NMAXOBJECTS];
1315 // // int fObjMCParentID[NMAXOBJECTS];
1316 // // int fObjMCGrandparentID[NMAXOBJECTS];
1317 // // // Hit information
1318 // // int fObjNumHits[NMAXOBJECTS];
1319 // // double fObjHitPosition[NMAXOBJECTS][NMAXHITS][3];
1320 // // double fObjHitTime[NMAXOBJECTS][NMAXHITS];
1321 // // int fObjHitWire[NMAXOBJECTS][NMAXHITS];
1322 // // double fObjHitCharge[NMAXOBJECTS][NMAXHITS];
1323 // // double fObjHitPitch[NMAXOBJECTS][NMAXHITS];
1324 // // double fObjHitdEdx[NMAXOBJECTS][NMAXHITS];
1325 // // double fObjHitdQdx[NMAXOBJECTS][NMAXHITS];
1326 //
1327 // // Only increment the counter if the object was contained in APA3.
1328 // ++obji;
1329 // }
1330 // }
1331 
1332 // -----------------------------------------------------------------------------
1333 // -----------------------------------------------------------------------------
1334 void protoana::ProtoDUNEPizeroAnaTree::setPiZeroInfo(art::Event const & evt,
1335  detinfo::DetectorClocksData const& clockData,
1336  detinfo::DetectorPropertiesData const& detProp,
1337  const std::vector<const simb::MCParticle*>& pi0s) {
1338  // Clean slate
1339  ResetPi0Vars();
1340 
1341  // Loop over all pi0s in the array by index.
1342  for(unsigned pi0i = 0; pi0i < NMAXPIZEROS && pi0i < pi0s.size(); ++pi0i) {
1343  pizero::PiZeroProcess pzproc(*(pi0s[pi0i]), evt, clockData, detProp, fShowerTag);
1344 
1345  // Set MC object variables
1346  if(pzproc.allMCSet()) {
1347 
1348  const simb::MCParticle* pi0 = pzproc.pi0();
1349  const simb::MCParticle* photon1 = pzproc.photon1();
1350  const simb::MCParticle* photon2 = pzproc.photon2();
1351 
1352  // MC pi0 variables
1353  fMCPi0ID[pi0i] = pi0->TrackId();
1354  fMCPi0Energy[pi0i] = pi0->E();
1355  pi0->Position().Vect().GetXYZ(fMCPi0StartPosition[pi0i]);
1356  pi0->EndPosition().Vect().GetXYZ(fMCPi0EndPosition[pi0i]);
1357  fMCPhoton1ID[pi0i] = photon1->TrackId();
1358  fMCPhoton2ID[pi0i] = photon2->TrackId();
1359  fMCPhoton1Energy[pi0i] = photon1->E();
1360  fMCPhoton2Energy[pi0i] = photon2->E();
1361  photon1->Position().Vect().GetXYZ(fMCPhoton1StartPosition[pi0i]);
1362  photon2->Position().Vect().GetXYZ(fMCPhoton2StartPosition[pi0i]);
1363  photon1->EndPosition().Vect().GetXYZ(fMCPhoton1EndPosition[pi0i]);
1364  photon2->EndPosition().Vect().GetXYZ(fMCPhoton2EndPosition[pi0i]);
1365  // Reco hits related to photons
1366  std::vector<const recob::Hit*> ph1hits = truthUtil.GetMCParticleHits(clockData, *photon1, evt, fHitTag);
1367  std::vector<const recob::Hit*> ph2hits = truthUtil.GetMCParticleHits(clockData, *photon2, evt, fHitTag);
1368  // Photon 1 hits
1369  art::FindManyP<recob::SpacePoint> ph1sps(ph1hits, evt, fPFParticleTag);
1370  unsigned phi1 = 0;
1371  for(unsigned i=0; i<ph1hits.size(); ++i) {
1372  if(ph1hits[i]->WireID().Plane != 2) continue;
1373  fMCPhoton1_hit_w[pi0i][phi1] = ph1hits[phi1]->WireID().Wire;
1374  fMCPhoton1_hit_t[pi0i][phi1] = ph1hits[phi1]->PeakTime();
1375  fMCPhoton1_hit_q[pi0i][phi1] = ph1hits[phi1]->Integral();
1376  std::vector<art::Ptr<recob::SpacePoint>> sps = ph1sps.at(i);
1377  if(!sps.empty()) {
1378  fMCPhoton1_hit_pos[pi0i][phi1][0] = sps[0]->XYZ()[0];
1379  fMCPhoton1_hit_pos[pi0i][phi1][1] = sps[0]->XYZ()[1];
1380  fMCPhoton1_hit_pos[pi0i][phi1][2] = sps[0]->XYZ()[2];
1381  }
1382  ++phi1;
1383  }
1384  // Photon 2 hits
1385  art::FindManyP<recob::SpacePoint> ph2sps(ph2hits, evt, fPFParticleTag);
1386  unsigned phi2 = 0;
1387  for(unsigned i=0; i<ph2hits.size(); ++i) {
1388  if(ph2hits[i]->WireID().Plane != 2) continue;
1389  fMCPhoton2_hit_w[pi0i][phi2] = ph2hits[phi2]->WireID().Wire;
1390  fMCPhoton2_hit_t[pi0i][phi2] = ph2hits[phi2]->PeakTime();
1391  fMCPhoton2_hit_q[pi0i][phi2] = ph2hits[phi2]->Integral();
1392  std::vector<art::Ptr<recob::SpacePoint>> sps = ph2sps.at(i);
1393  if(!sps.empty()) {
1394  fMCPhoton2_hit_pos[pi0i][phi2][0] = sps[0]->XYZ()[0];
1395  fMCPhoton2_hit_pos[pi0i][phi2][1] = sps[0]->XYZ()[1];
1396  fMCPhoton2_hit_pos[pi0i][phi2][2] = sps[0]->XYZ()[2];
1397  }
1398  ++phi2;
1399  }
1400  fMCPhoton1NumHits[pi0i] = phi1;
1401  fMCPhoton2NumHits[pi0i] = phi2;
1402 
1403  } // if MC set
1404 
1405  //This is how we get cnn score for now
1406  auto recoShowers = evt.getValidHandle< std::vector< recob::Shower > >(fShowerTag);
1407  auto recoTracks = evt.getValidHandle< std::vector< recob::Track > >(fTrackTag);
1408  art::FindManyP<recob::Hit> findHitsFromShowers(recoShowers,evt,fShowerTag);
1409  art::FindManyP<recob::Hit> findHitsFromTracks(recoTracks,evt,fTrackTag);
1410  anab::MVAReader<recob::Hit,4> hitResults(evt, "emtrkmichelid:emtrkmichel" );
1411 
1412  // First shower
1413  const recob::Shower* shower1 = pzproc.shower1();
1414  if(shower1 != 0x0) {
1415  fShower1ID[pi0i] = showerUtil.GetShowerIndex(*shower1, evt, fShowerTag);
1416  std::vector< art::Ptr< recob::Hit > > sh1_hits = findHitsFromShowers.at(fShower1ID[pi0i]);
1417  for( size_t i=0; i<sh1_hits.size(); ++i){
1418  std::array<float,4> cnn_out = hitResults.getOutput( sh1_hits[i] );
1419  double p_trk_or_sh = cnn_out[ hitResults.getIndex("track") ]+ cnn_out[ hitResults.getIndex("em") ];
1420  double cnn_score = cnn_out[ hitResults.getIndex("em") ]/p_trk_or_sh;
1421  fShower1_cnn_sc[pi0i][i] = cnn_score;
1422  }
1423  shower1->ShowerStart().GetXYZ(fShower1StartPosition[pi0i]);
1424  shower1->Direction().GetXYZ(fShower1Direction[pi0i]);
1425  fShower1Length[pi0i] = shower1->Length();
1426  unsigned numMChits = truthUtil.GetMCParticleHits(clockData, *pzproc.photon1(), evt, fHitTag).size();
1427  unsigned showerHits = showerUtil.GetRecoShowerHits(*shower1, evt, fShowerTag).size();
1428  unsigned sharedHits = truthUtil.GetSharedHits(clockData, *pzproc.photon1(), *shower1, evt, fShowerTag, true).size();
1429  fShower1Completeness[pi0i] = (double)sharedHits/numMChits;
1430  fShower1Purity[pi0i] = (double)sharedHits/showerHits;
1431  // Find out whether (grand)parent of shower was primary.
1432  // Go through all PFParticles until the right one is found.
1433  const recob::PFParticle* shpfp = 0x0;
1434  std::vector<const recob::PFParticle*> beamPFParticles = pfpUtil.GetPFParticlesFromBeamSlice(evt,fPFParticleTag);
1435  auto PFParticleHandle = evt.getHandle<std::vector<recob::PFParticle> >(fPFParticleTag);
1436  if (PFParticleHandle && beamPFParticles.size() != 0) {
1437  const recob::PFParticle* beampfp = beamPFParticles[0];
1438  for(const recob::PFParticle& pfp : *PFParticleHandle) {
1439  const recob::Shower* thish = pfpUtil.GetPFParticleShower(pfp,evt,fPFParticleTag,fShowerTag);
1440  if(thish == 0x0) continue;
1441  // Test whether it's the right shower by comparing lengths.
1442  if(thish->Length() - shower1->Length() < 1e-5) {
1443  shpfp = &pfp;
1444  break;
1445  }
1446  }
1447  // Is PFP's parent or PFP's grandparent the beam particle?
1448  const recob::PFParticle* parent = 0x0;
1449  const recob::PFParticle* gparent = 0x0;
1450  if(shpfp && shpfp->Parent() != 0 && shpfp->Parent() < PFParticleHandle->size())
1451  parent = &PFParticleHandle->at(shpfp->Parent());
1452  if(parent && parent->Parent() != 0 && parent->Parent() < PFParticleHandle->size()) gparent = &PFParticleHandle->at(parent->Parent());
1453 
1454  fShower1HasBeamParent[pi0i] = parent && parent->Self() == beampfp->Self()? 1: 0;
1455  fShower1HasBeamGrandparent[pi0i] = gparent && gparent->Self() == beampfp->Self()? 1: 0;
1456  } // if get pfparticles and beampfparticles
1457 
1458  // Calorimetry related to showers
1459  std::vector<anab::Calorimetry> sh1calo = showerUtil.GetRecoShowerCalorimetry(*shower1, evt, fShowerTag, fShowerCaloTag);
1460  for(unsigned k=0; k<sh1calo.size(); ++k) {
1461  if(sh1calo[k].PlaneID().Plane != 2) continue;
1462 
1463  fShower1NumHits[pi0i] = std::min((int)sh1calo[k].dEdx().size(), NMAXHITS);
1464  fShower1_cal_E[pi0i] = sh1calo[k].KineticEnergy();
1465  std::vector<const recob::Hit*> shHitPtrs = showerUtil.GetRecoShowerHits(
1466  *shower1, evt, fShowerTag);
1467  fShower1EnergyFromHits[pi0i] = showerUtil.EstimateEnergyFromHitCharge(clockData,
1468  detProp,
1469  shHitPtrs,
1470  fCalorimetryAlg)[2];
1471 
1472  for(int i=0; i<fShower1NumHits[pi0i]; ++i) {
1473  const geo::Point_t& pos = sh1calo[k].XYZ()[i];
1474  fShower1_cal_pos[pi0i][i][0] = pos.X();
1475  fShower1_cal_pos[pi0i][i][1] = pos.Y();
1476  fShower1_cal_pos[pi0i][i][2] = pos.Z();
1477  fShower1_cal_pitch[pi0i][i] = sh1calo[k].TrkPitchVec()[i];
1478  fShower1_cal_dEdx[pi0i][i] = sh1calo[k].dEdx()[i];
1479  fShower1_cal_dQdx[pi0i][i] = sh1calo[k].dQdx()[i];
1480  }
1481  }
1482  } // if shower1 != 0x0
1483  // Second shower
1484  const recob::Shower* shower2 = pzproc.shower2();
1485  if(shower2 != 0x0) {
1486  // Reco pi0 variables
1487  fShower2ID[pi0i] = showerUtil.GetShowerIndex(*shower2, evt, fShowerTag);
1488  std::vector< art::Ptr< recob::Hit > > sh2_hits = findHitsFromShowers.at(fShower2ID[pi0i]);
1489  for( size_t i=0; i<sh2_hits.size(); ++i){
1490  std::array<float,4> cnn_out = hitResults.getOutput( sh2_hits[i] );
1491  double p_trk_or_sh = cnn_out[ hitResults.getIndex("track") ]+ cnn_out[ hitResults.getIndex("em") ];
1492  double cnn_score = cnn_out[ hitResults.getIndex("em") ]/p_trk_or_sh;
1493  fShower2_cnn_sc[pi0i][i] = cnn_score;
1494  }
1495  shower2->ShowerStart().GetXYZ(fShower2StartPosition[pi0i]);
1496  shower2->Direction().GetXYZ(fShower2Direction[pi0i]);
1497  fShower2Length[pi0i] = shower2->Length();
1498  unsigned numMChits = truthUtil.GetMCParticleHits(clockData, *pzproc.photon2(), evt, fHitTag).size();
1499  unsigned showerHits = showerUtil.GetRecoShowerHits(*shower2, evt, fShowerTag).size();
1500  unsigned sharedHits = truthUtil.GetSharedHits(clockData, *pzproc.photon2(), *shower2, evt, fShowerTag, true).size();
1501  fShower2Completeness[pi0i] = (double)sharedHits/numMChits;
1502  fShower2Purity[pi0i] = (double)sharedHits/showerHits;
1503  // Find out whether (grand)parent of shower was primary.
1504  // Go through all PFParticles until the right one is found.
1505  const recob::PFParticle* shpfp = 0x0;
1506  std::vector<const recob::PFParticle*> beamPFParticles = pfpUtil.GetPFParticlesFromBeamSlice(evt,fPFParticleTag);
1507  auto PFParticleHandle = evt.getHandle<std::vector<recob::PFParticle>>(fPFParticleTag);
1508  if (PFParticleHandle && beamPFParticles.size() != 0) {
1509  const recob::PFParticle* beampfp = beamPFParticles[0];
1510  for(const recob::PFParticle& pfp : *PFParticleHandle) {
1511  const recob::Shower* thish = pfpUtil.GetPFParticleShower(pfp,evt,fPFParticleTag,fShowerTag);
1512  if(thish == 0x0) continue;
1513  // Test whether it's the right shower by comparing lengths.
1514  if(thish->Length() - shower2->Length() < 1e-5) {
1515  shpfp = &pfp;
1516  break;
1517  }
1518  }
1519  // Is PFP's parent or PFP's grandparent the beam particle?
1520  const recob::PFParticle* parent = 0x0;
1521  const recob::PFParticle* gparent = 0x0;
1522  if(shpfp && shpfp->Parent() != 0 && shpfp->Parent() < PFParticleHandle->size()) parent = &PFParticleHandle->at(shpfp->Parent());
1523  if(parent && parent->Parent() != 0 && parent->Parent() < PFParticleHandle->size()) gparent = &PFParticleHandle->at(parent->Parent());
1524 
1525  fShower2HasBeamParent[pi0i] = parent && parent->Self() == beampfp->Self()? 1: 0;
1526  fShower2HasBeamGrandparent[pi0i] = gparent && gparent->Self() == beampfp->Self()? 1: 0;
1527  } // if get pfparticles and beampfparticles
1528 
1529  // Calorimetry related to showers
1530  std::vector<anab::Calorimetry> sh2calo = showerUtil.GetRecoShowerCalorimetry(*shower2, evt, fShowerTag, fShowerCaloTag);
1531  for(unsigned k=0; k<sh2calo.size(); ++k) {
1532  if(sh2calo[k].PlaneID().Plane != 2) continue;
1533 
1534  fShower2NumHits[pi0i] = std::min((int)sh2calo[k].dEdx().size(), NMAXHITS);
1535  fShower2_cal_E[pi0i] = sh2calo[k].KineticEnergy();
1536  std::vector<const recob::Hit*> shHitPtrs = showerUtil.GetRecoShowerHits(
1537  *shower2, evt, fShowerTag);
1538  fShower2EnergyFromHits[pi0i] = showerUtil.EstimateEnergyFromHitCharge(clockData, detProp,
1539  shHitPtrs, fCalorimetryAlg)[2];
1540  for(int i=0; i<fShower2NumHits[pi0i]; ++i) {
1541  const geo::Point_t& pos = sh2calo[k].XYZ()[i];
1542  fShower2_cal_pos[pi0i][i][0] = pos.X();
1543  fShower2_cal_pos[pi0i][i][1] = pos.Y();
1544  fShower2_cal_pos[pi0i][i][2] = pos.Z();
1545  fShower2_cal_pitch[pi0i][i] = sh2calo[k].TrkPitchVec()[i];
1546  fShower2_cal_dEdx[pi0i][i] = sh2calo[k].dEdx()[i];
1547  fShower2_cal_dQdx[pi0i][i] = sh2calo[k].dQdx()[i];
1548  }
1549  }
1550  } // if shower2 != 0x0
1551  } // Loop over pi0s
1552 
1553 } // setPiZeroInfo
1554 
1555 
1556 // -----------------------------------------------------------------------------
1557 // -----------------------------------------------------------------------------
1558 void protoana::ProtoDUNEPizeroAnaTree::Initialise(){
1559  fRun = -999;
1560  fSubRun = -999;
1561  fevent = -999;
1562  fTimeStamp = -999.0;
1563  for(int k=0; k < 5; k++)
1564  fNactivefembs[k] = -999;
1565 
1566  for(int k=0; k < 3; k++){
1567  fbeamtrackPos[k] = -999.0;
1568  fbeamtrackEndPos[k] = -999.0;
1569  fbeamtrackDir[k] = -999.0;
1570  fprimaryTruth_vtx[k]= -999.0;
1571  fprimaryVertex[k] = -999.0;
1572  fprimaryEndPosition[k] = -999.0;
1573  fprimaryStartPosition[k] = -999.0;
1574  fprimaryEndDirection[k] = -999.0;
1575  fprimaryStartDirection[k] = -999.0;
1576  fprimaryKineticEnergy[k] = -999.0;
1577  fprimaryRange[k] = -999.0;
1578  }
1579  for( int l=0; l<1000; l ++) {
1580  fbeamtrackPos_at[l][0] = -999.;
1581  fbeamtrackPos_at[l][1] = -999.;
1582  fbeamtrackPos_at[l][2] = -999.;
1583  fbeamtrackMom_at[l][0] = -999.;
1584  fbeamtrackMom_at[l][1] = -999.;
1585  fbeamtrackMom_at[l][2] = -999.;
1586  fbeamtrackMom_at[l][3] = -999.;
1587  }
1588  fprimaryShower_nHits =0;
1589  for( int m=0; m<NMAXHITS; m ++){
1590  fprimarydEdx[m]= -999.0;
1591  fprimarydQdx[m]= -999.0;
1592  fprimary_cal_pos[m][0] = -999.0;
1593  fprimary_cal_pos[m][1] = -999.0;
1594  fprimary_cal_pos[m][2] = -999.0;
1595  fprimary_cal_pitch[m]=-999.0;
1596  fprimaryResidualRange[m] = -999.0;
1597 
1598  fprimaryShower_hit_w[m] =-999.0;
1599  fprimaryShower_hit_q[m] =-999.0;
1600  fprimaryShower_hit_t[m] =-999.0;
1601  fprimaryShower_hit_pos[m][0] =-999.0;
1602  fprimaryShower_hit_pos[m][1] =-999.0;
1603  fprimaryShower_hit_pos[m][2] =-999.0;
1604  }
1605  fprimaryTruth_trkID =-999;
1606  fprimaryTruth_pdg = -999;
1607  fprimaryTruth_E = -999;
1608  fprimarynCal = 0;
1609  fbeamCheckIsMatched = -999;
1610  fbeamtrigger = -999;
1611  ftof = -999.0;
1612  for(int k=0; k < 2; k++){
1613  fcerenkovPressure[k] = -999.0;
1614  fcerenkovTime[k] = -999.0;
1615  fcerenkovStatus[k] = -999;
1616  }
1617  fbeamtrackMomentum = -999.0;
1618  fbeamtrackEnergy = 999.0;
1619  fbeamtrackPdg = -999;
1620  fbeamtrackTime = -999.0;
1621  fbeamtrackID = -999;
1622  fbeam_ntrjPoints =0;
1623  fbeamtrackNDaughters= 0;
1624  fprimaryIstrack = -999;
1625  fprimaryIsshower = -999;
1626 
1627  fprimaryBDTScore = -999.0;
1628  fprimaryCNNScore = -999.0;
1629  fprimaryNHits = 0;
1630  fprimaryTheta = -999.0;
1631  fprimaryPhi = -999.0;
1632  fprimaryLength = -999.0;
1633  fprimaryMomentum = -999.0;
1634  fprimaryEndMomentum = -999.0;
1635  fprimaryOpeningAngle = -999.0;
1636  fprimaryShowerBestPlane = -999;
1637  fprimaryShowerEnergy = -999.0;
1638  fprimaryShowerCharge = -999.0;
1639  fprimaryShowerMIPEnergy = -999.0;
1640  fprimaryID = -999;
1641  fprimaryMomentumByRangeProton = -999.0;
1642  fprimaryT0 = -999.0;
1643  fprimaryNDaughters = 0;
1644 
1645  for(unsigned k = 0; k < NMAXDAUGHTERS; ++k) {
1646  fprimaryDaughterID[k] = -999;
1647  fprimaryDaughterIstrack[k] = -999;
1648  fprimaryDaughterIsshower[k] = -999;
1649  fprimaryDaughterBDTScore[k] = -999.0;
1650  fprimaryDaughterCNNScore[k] = -999.0;
1651  fprimaryDaughterNCollHits[k] = 0;
1652  fprimaryDaughterNHits[k] = 0;
1653  fprimaryDaughterEnergy[k] = -999.0;
1654  fprimaryDaughterEnergyFromHits[k] = -999.0;
1655  fprimaryDaughterMomentum[k] = -999.0;
1656  fprimaryDaughterEndMomentum[k] = -999.0;
1657  fprimaryDaughterLength[k] = -999.0;
1658  for(unsigned l = 0; l < 3; ++l) {
1659  fprimaryDaughterStartPosition[k][l] = -999.0;
1660  fprimaryDaughterEndPosition[k][l] = -999.0;
1661  fprimaryDaughterStartDirection[k][l] = -999.0;
1662  fprimaryDaughterEndDirection[k][l] = -999.0;
1663  }
1664  fprimaryDaughterParentPdg[k] = -999;
1665  fprimaryDaughterGrandparentPdg[k] = -999;
1666  fprimaryDaughterParentID[k] = -999;
1667  fprimaryDaughterGrandparentID[k] = -999;
1668  fprimaryDaughterParentE[k] = -999.0;
1669  fprimaryDaughterGrandparentE[k] = -999.0;
1670  fprimaryDaughterParentStart[k][0] = -999.0;
1671  fprimaryDaughterParentStart[k][1] = -999.0;
1672  fprimaryDaughterParentStart[k][2] = -999.0;
1673  fprimaryDaughterGrandparentStart[k][0] = -999.0;
1674  fprimaryDaughterGrandparentStart[k][1] = -999.0;
1675  fprimaryDaughterGrandparentStart[k][2] = -999.0;
1676  fprimaryDaughterParentEnd[k][0] = -999.0;
1677  fprimaryDaughterParentEnd[k][1] = -999.0;
1678  fprimaryDaughterParentEnd[k][2] = -999.0;
1679  fprimaryDaughterGrandparentEnd[k][0] = -999.0;
1680  fprimaryDaughterGrandparentEnd[k][1] = -999.0;
1681  fprimaryDaughterGrandparentEnd[k][2] = -999.0;
1682  fprimaryDaughterParentEndProcess[k] = -999.0;
1683  for(unsigned l = 0; l < NMAXHITS; ++l) {
1684  // fprimaryDaughterHitdQdx[k][l] = -999.0;
1685  fprimaryDaughterCaloHitdQdx[k][l] = -999.0;
1686  fprimaryDaughterCaloHitPitch[k][l] = -999.0;
1687  fprimaryDaughterOwnHitPitch[k][l] = -999.0;
1688  fprimaryDaughterHitCharge[k][l] = -999.0;
1689  // fprimaryDaughterHitdEdxArea[k][l] = -999.0;
1690  fprimaryDaughterHitWire[k][l] = -999;
1691  fprimaryDaughterHitTime[k][l] = -999.0;
1692  fprimaryDaughterHitPosition[k][l][0] = -999.0;
1693  fprimaryDaughterHitPosition[k][l][1] = -999.0;
1694  fprimaryDaughterHitPosition[k][l][2] = -999.0;
1695  // fprimaryDaughterCaloHitPosition[k][l][0] = -999.0;
1696  // fprimaryDaughterCaloHitPosition[k][l][1] = -999.0;
1697  // fprimaryDaughterCaloHitPosition[k][l][2] = -999.0;
1698  }
1699  }
1700 
1701  // // Instead of primary daughters, look at all objects in APA3
1702  // for(int obji = 0; obji < NMAXOBJECTS; ++obji) {
1703  // fObjID[obji] = -999;
1704  // fObjIsTrack[obji] = -999;
1705  // fObjIsShower[obji] = -999;
1706  // fObjIsPrimaryDaughter[obji] = -999;
1707  // fObjBDTscore[obji] = -999.0;
1708  // fObjCNNscore[obji] = -999.0;
1709  // fObjMomentum[obji] = -999.0;
1710  // fObjEndMomentum[obji] = -999.0;
1711  // fObjLength[obji] = -999.0;
1712  // for(int pi = 0; pi < 3; ++pi) {
1713  // fObjStartPosition[obji][pi] = -999.0;
1714  // fObjEndPosition[obji][pi] = -999.0;
1715  // fObjStartDirection[obji][pi] = -999.0;
1716  // fObjEndDirection[obji][pi] = -999.0;
1717  // }
1718  // fObjMCParentPdg[obji] = -999;
1719  // fObjMCGrandparentPdg[obji] = -999;
1720  // fObjMCParentID[obji] = -999;
1721  // fObjMCGrandparentID[obji] = -999;
1722  // // Hit information
1723  // fObjNumHits[obji] = -999;
1724  // for(int hiti = 0; hiti < NMAXHITS; ++hiti) {
1725  // fObjHitPosition[obji][hiti][0] = -999.0;
1726  // fObjHitPosition[obji][hiti][1] = -999.0;
1727  // fObjHitPosition[obji][hiti][2] = -999.0;
1728  // fObjHitTime[obji][hiti] = -999.0;
1729  // fObjHitWire[obji][hiti] = -999;
1730  // fObjHitCharge[obji][hiti] = -999.0;
1731  // fObjHitPitch[obji][hiti] = -999.0;
1732  // fObjHitdEdx[obji][hiti] = -999.0;
1733  // fObjHitdQdx[obji][hiti] = -999.0;
1734  // }
1735  // }
1736 
1737  ResetPi0Vars();
1738 } // Initialise
1739 
1740 void protoana::ProtoDUNEPizeroAnaTree::ResetPi0Vars() {
1741  for(int pi0i = 0; pi0i < NMAXPIZEROS; ++pi0i) {
1742  // MC pi0 variables
1743  fMCPi0ID[pi0i] = -999;
1744  fMCPi0Energy[pi0i] = -999.0;
1745  for(int i=0; i<3; ++i) {
1746  fMCPi0StartPosition[pi0i][i] = -999.0;
1747  fMCPi0EndPosition[pi0i][i] = -999.0;
1748  }
1749  fMCPhoton1ID[pi0i] = -999;
1750  fMCPhoton2ID[pi0i] = -999;
1751  fMCPhoton1Energy[pi0i] = -999.0;
1752  fMCPhoton2Energy[pi0i] = -999.0;
1753  for(int i=0; i<3; ++i) {
1754  fMCPhoton1StartPosition[pi0i][i] = -999.0;
1755  fMCPhoton2StartPosition[pi0i][i] = -999.0;
1756  fMCPhoton1EndPosition[pi0i][i] = -999.0;
1757  fMCPhoton2EndPosition[pi0i][i] = -999.0;
1758  }
1759  // Reco hits related to photons
1760  fMCPhoton1NumHits[pi0i] = -999;
1761  fMCPhoton2NumHits[pi0i] = -999;
1762  for(int i=0; i < NMAXHITS; ++i) {
1763  fMCPhoton1_hit_w[pi0i][i] = -999;
1764  fMCPhoton2_hit_w[pi0i][i] = -999;
1765  fMCPhoton1_hit_t[pi0i][i] = -999.0;
1766  fMCPhoton2_hit_t[pi0i][i] = -999.0;
1767  fMCPhoton1_hit_q[pi0i][i] = -999.0;
1768  fMCPhoton2_hit_q[pi0i][i] = -999.0;
1769  fMCPhoton1_hit_pos[pi0i][i][0] = -999.0;
1770  fMCPhoton2_hit_pos[pi0i][i][0] = -999.0;
1771  fMCPhoton1_hit_pos[pi0i][i][1] = -999.0;
1772  fMCPhoton2_hit_pos[pi0i][i][1] = -999.0;
1773  fMCPhoton1_hit_pos[pi0i][i][2] = -999.0;
1774  fMCPhoton2_hit_pos[pi0i][i][2] = -999.0;
1775  }
1776 
1777  // Reco pi0 variables
1778  fShower1ID[pi0i] = -999;
1779  fShower2ID[pi0i] = -999;
1780  for(int i=0; i<3; ++i) {
1781  fShower1StartPosition[pi0i][i] = -999.0;
1782  fShower2StartPosition[pi0i][i] = -999.0;
1783  fShower1Direction[pi0i][i] = -999.0;
1784  fShower2Direction[pi0i][i] = -999.0;
1785  }
1786  fShower1Length[pi0i] = -999.0;
1787  fShower2Length[pi0i] = -999.0;
1788  fShower1Completeness[pi0i] = -999.0;
1789  fShower2Completeness[pi0i] = -999.0;
1790  fShower1Purity[pi0i] = -999.0;
1791  fShower2Purity[pi0i] = -999.0;
1792  // Reco hits related to showers
1793  fShower1NumHits[pi0i] = -999;
1794  fShower2NumHits[pi0i] = -999;
1795  fShower1_cal_E[pi0i] = -999.0;
1796  fShower2_cal_E[pi0i] = -999.0;
1797  fShower1Energy[pi0i] = -999.0;
1798  fShower2Energy[pi0i] = -999.0;
1799  fShower1EnergyFromHits[pi0i] = -999.0;
1800  fShower2EnergyFromHits[pi0i] = -999.0;
1801  fShower1HasBeamParent[pi0i] = -999;
1802  fShower2HasBeamParent[pi0i] = -999;
1803  fShower1HasBeamGrandparent[pi0i] = -999;
1804  fShower2HasBeamGrandparent[pi0i] = -999;
1805  for(int i=0; i<NMAXHITS; ++i) {
1806  fShower1_cnn_sc[pi0i][i] = -999.0;
1807  fShower2_cnn_sc[pi0i][i] = -999.0;
1808  fShower1_cal_pos[pi0i][i][0] = -999.0;
1809  fShower2_cal_pos[pi0i][i][0] = -999.0;
1810  fShower1_cal_pos[pi0i][i][1] = -999.0;
1811  fShower2_cal_pos[pi0i][i][1] = -999.0;
1812  fShower1_cal_pos[pi0i][i][2] = -999.0;
1813  fShower2_cal_pos[pi0i][i][2] = -999.0;
1814  fShower1_cal_pitch[pi0i][i] = -999.0;
1815  fShower2_cal_pitch[pi0i][i] = -999.0;
1816  fShower1_cal_dEdx[pi0i][i] = -999.0;
1817  fShower2_cal_dEdx[pi0i][i] = -999.0;
1818  fShower1_cal_dQdx[pi0i][i] = -999.0;
1819  fShower2_cal_dQdx[pi0i][i] = -999.0;
1820  }
1821  } // Loop over NMAXPIZEROS
1822 } // ResetPi0Vars
1823 
1824 DEFINE_ART_MODULE(protoana::ProtoDUNEPizeroAnaTree)
simb::MCParticle::E
double E(const int i=0) const
Definition: MCParticle.h:233
protoana::ProtoDUNEPizeroAnaTree::fShower2_cal_dEdx
double fShower2_cal_dEdx[NMAXPIZEROS][NMAXHITS]
Definition: ProtoDUNEPizeroAnaTree_module.cc:347
geo::WireGeo
Geometry description of a TPC wireThe wire is a single straight segment on a wire plane...
Definition: WireGeo.h:65
simb::MCParticle::NumberTrajectoryPoints
unsigned int NumberTrajectoryPoints() const
Definition: MCParticle.h:218
recob::Shower::ShowerStart
const TVector3 & ShowerStart() const
Definition: Shower.h:192
simb::MCParticle::Position
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:219
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::ProtoDUNEPizeroAnaTree::fShower1_cal_pos
double fShower1_cal_pos[NMAXPIZEROS][NMAXHITS][3]
Definition: ProtoDUNEPizeroAnaTree_module.cc:332
protoana::ProtoDUNEShowerUtils::EstimateEnergyFromHitCharge
std::vector< double > EstimateEnergyFromHitCharge(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const std::vector< const recob::Hit * > &hits, calo::CalorimetryAlg caloAlg)
Definition: ProtoDUNEShowerUtils.cxx:78
protoana::ProtoDUNEPizeroAnaTree::fprimaryDaughterStartDirection
double fprimaryDaughterStartDirection[NMAXDAUGHTERS][3]
Definition: ProtoDUNEPizeroAnaTree_module.cc:236
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::ProtoDUNEPizeroAnaTree::fShower1_cal_pitch
double fShower1_cal_pitch[NMAXPIZEROS][NMAXHITS]
Definition: ProtoDUNEPizeroAnaTree_module.cc:344
recob::PFParticle::NumDaughters
int NumDaughters() const
Returns the number of daughter particles flowing from this one.
Definition: PFParticle.h:89
protoana::ProtoDUNETruthUtils::GetGeantGoodParticle
const simb::MCParticle * GetGeantGoodParticle(const simb::MCTruth &genTruth, const art::Event &evt) const
Definition: ProtoDUNETruthUtils.cxx:860
protoana::ProtoDUNEPizeroAnaTree::fShower2_cal_pitch
double fShower2_cal_pitch[NMAXPIZEROS][NMAXHITS]
Definition: ProtoDUNEPizeroAnaTree_module.cc:345
simb::MCParticle::EndPosition
const TLorentzVector & EndPosition() const
Definition: MCParticle.h:225
simb::MCParticle::EndZ
double EndZ() const
Definition: MCParticle.h:228
recob::PFParticle::Self
size_t Self() const
Returns the index of this particle.
Definition: PFParticle.h:92
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
protoana::ProtoDUNEPizeroAnaTree::setPiZeroInfo
void setPiZeroInfo(art::Event const &evt, detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const std::vector< const simb::MCParticle * > &pi0s)
Definition: ProtoDUNEPizeroAnaTree_module.cc:1334
string
std::string string
Definition: nybbler.cc:12
protoana::ProtoDUNEPizeroAnaTree::fSCEService
spacecharge::SpaceChargeService::provider_type const * fSCEService
Definition: ProtoDUNEPizeroAnaTree_module.cc:353
protoana::ProtoDUNEPizeroAnaTree::fprimaryDaughterCaloHitdQdx
double fprimaryDaughterCaloHitdQdx[NMAXDAUGHTERS][NMAXHITS]
Definition: ProtoDUNEPizeroAnaTree_module.cc:250
spacecharge::SpaceCharge::GetCalPosOffsets
virtual geo::Vector_t GetCalPosOffsets(geo::Point_t const &point, int const &TPCid) const =0
cheat::ParticleInventoryService::TrackIdToParticle_P
const simb::MCParticle * TrackIdToParticle_P(int id) const
Definition: ParticleInventoryService_service.cc:128
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::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::ProtoDUNEPizeroAnaTree::FillAPA3Object
void FillAPA3Object(art::Event const &evt, art::Handle< std::vector< recob::PFParticle >> pfpHandle)
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::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::ProtoDUNEPizeroAnaTree::fprimaryDaughterCaloHitPitch
double fprimaryDaughterCaloHitPitch[NMAXDAUGHTERS][NMAXHITS]
Definition: ProtoDUNEPizeroAnaTree_module.cc:251
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.
geo::Vector_t
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space.
Definition: geo_vectors.h:164
art::EDAnalyzer::EDAnalyzer
EDAnalyzer(fhicl::ParameterSet const &pset)
Definition: EDAnalyzer.h:25
MCParticle.h
Particle class.
simb::MCParticle::EndY
double EndY() const
Definition: MCParticle.h:227
l
static QStrList * l
Definition: config.cpp:1044
protoana::ProtoDUNEPizeroAnaTree::fShower2_cal_dQdx
double fShower2_cal_dQdx[NMAXPIZEROS][NMAXHITS]
Definition: ProtoDUNEPizeroAnaTree_module.cc:349
geo::GeometryCore::WirePitch
geo::Length_t WirePitch(geo::PlaneID const &planeid) const
Returns the distance between two consecutive wires.
Definition: GeometryCore.cxx:824
protoana::ProtoDUNEPizeroAnaTree::fprimaryDaughterHitTime
double fprimaryDaughterHitTime[NMAXDAUGHTERS][NMAXHITS]
Definition: ProtoDUNEPizeroAnaTree_module.cc:256
simb::MCParticle::NumberDaughters
int NumberDaughters() const
Definition: MCParticle.h:217
Geometry.h
art framework interface to geometry description
simb::MCParticle::TrackId
int TrackId() const
Definition: MCParticle.h:210
simb::MCParticle::Daughter
int Daughter(const int i) const
Definition: MCParticle.cxx:112
protoana::ProtoDUNEPizeroAnaTree::fMCPhoton2_hit_pos
double fMCPhoton2_hit_pos[NMAXPIZEROS][NMAXHITS][3]
Definition: ProtoDUNEPizeroAnaTree_module.cc:311
NMAXHITS
const int NMAXHITS
Definition: ProtoDUNEPizeroAnaTree_module.cc:85
pizero::PiZeroProcess::shower2
const recob::Shower * shower2() const
Definition: PiZeroProcess.h:153
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
protoana::ProtoDUNETruthUtils::GetMCParticleFromReco
const simb::MCParticle * GetMCParticleFromReco(detinfo::DetectorClocksData const &clockData, const recob::PFParticle &pfpart, art::Event const &evt, std::string pfparticleModule) const
Definition: ProtoDUNETruthUtils.cxx:696
art::DataViewImpl::isRealData
bool isRealData() const
Definition: DataViewImpl.cc:113
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
protoana::ProtoDUNEPizeroAnaTree::fMCPhoton2_hit_q
double fMCPhoton2_hit_q[NMAXPIZEROS][NMAXHITS]
Definition: ProtoDUNEPizeroAnaTree_module.cc:309
dEdx
double dEdx(float dqdx, float Efield)
Definition: doAna.cpp:21
pizero::PiZeroProcess::photon1
const simb::MCParticle * photon1() const
Definition: PiZeroProcess.h:144
recob::Track::Phi
double Phi() const
Definition: Track.h:178
pizero::PiZeroProcess::allMCSet
bool allMCSet() const
Definition: PiZeroProcess.h:281
NMAXPIZEROS
const int NMAXPIZEROS
Definition: ProtoDUNEPizeroAnaTree_module.cc:83
e
const double e
Definition: gUpMuFluxGen.cxx:165
recob::Track::Length
double Length(size_t p=0) const
Access to various track properties.
Definition: Track.h:167
protoana::ProtoDUNEPizeroAnaTree::fShower1_cal_dQdx
double fShower1_cal_dQdx[NMAXPIZEROS][NMAXHITS]
Definition: ProtoDUNEPizeroAnaTree_module.cc:348
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
simb::MCParticle::EndProcess
std::string EndProcess() const
Definition: MCParticle.h:216
pizero::PiZeroProcess::pi0
const simb::MCParticle * pi0() const
Definition: PiZeroProcess.h:137
protoana::ProtoDUNEShowerUtils::GetShowerIndex
int GetShowerIndex(const recob::Shower &shower, art::Event const &evt, const std::string showerModule) const
If the shower.ID() isn&#39;t filled we must find the actual shower index ourselves.
Definition: ProtoDUNEShowerUtils.cxx:106
protoana::ProtoDUNEPizeroAnaTree::fShower1_cal_dEdx
double fShower1_cal_dEdx[NMAXPIZEROS][NMAXHITS]
Definition: ProtoDUNEPizeroAnaTree_module.cc:346
recob::PFParticle::Parent
size_t Parent() const
Definition: PFParticle.h:96
recob::Track::Theta
double Theta() const
Access to spherical or geographical angles at vertex or at any point.
Definition: Track.h:176
protoana::ProtoDUNEPizeroAnaTree::operator=
ProtoDUNEPizeroAnaTree & operator=(ProtoDUNEPizeroAnaTree const &)=delete
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
protoana::ProtoDUNEPizeroAnaTree::fprimaryDaughterHitPosition
double fprimaryDaughterHitPosition[NMAXDAUGHTERS][NMAXHITS][3]
Definition: ProtoDUNEPizeroAnaTree_module.cc:257
simb::MCParticle::T
double T(const int i=0) const
Definition: MCParticle.h:224
art::DataViewImpl::getValidHandle
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
Definition: DataViewImpl.h:441
test.p
p
Definition: test.py:223
art::DataViewImpl::subRun
SubRunNumber_t subRun() const
Definition: DataViewImpl.cc:78
protoana::ProtoDUNEPizeroAnaTree::fMCPhoton1_hit_t
double fMCPhoton1_hit_t[NMAXPIZEROS][NMAXHITS]
Definition: ProtoDUNEPizeroAnaTree_module.cc:306
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
NMAXDAUGHTERS
const int NMAXDAUGHTERS
Definition: ProtoDUNEPizeroAnaTree_module.cc:79
protoana::ProtoDUNEPizeroAnaTree::fMCPhoton2_hit_t
double fMCPhoton2_hit_t[NMAXPIZEROS][NMAXHITS]
Definition: ProtoDUNEPizeroAnaTree_module.cc:307
art::DataViewImpl::run
RunNumber_t run() const
Definition: DataViewImpl.cc:71
protoana::ProtoDUNEPizeroAnaTree::fprimaryDaughterHitWire
int fprimaryDaughterHitWire[NMAXDAUGHTERS][NMAXHITS]
Definition: ProtoDUNEPizeroAnaTree_module.cc:255
geo::GeometryCore
Description of geometry of one entire detector.
Definition: GeometryCore.h:1468
protoana::ProtoDUNEPizeroAnaTree::fprimaryDaughterOwnHitPitch
double fprimaryDaughterOwnHitPitch[NMAXDAUGHTERS][NMAXHITS]
Definition: ProtoDUNEPizeroAnaTree_module.cc:252
protoana::ProtoDUNEPizeroAnaTree::beamlineUtil
protoana::ProtoDUNEBeamlineUtils beamlineUtil
Definition: ProtoDUNEPizeroAnaTree_module.cc:116
hit
Detector simulation of raw signals on wires.
Definition: NumberOfHitsFilter_module.cc:25
protoana::ProtoDUNEPizeroAnaTree::fShower2_cal_pos
double fShower2_cal_pos[NMAXPIZEROS][NMAXHITS][3]
Definition: ProtoDUNEPizeroAnaTree_module.cc:333
protoana::ProtoDUNEPizeroAnaTree::fprimaryDaughterHitCharge
double fprimaryDaughterHitCharge[NMAXDAUGHTERS][NMAXHITS]
Definition: ProtoDUNEPizeroAnaTree_module.cc:253
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
protoana::ProtoDUNEPizeroAnaTree::bt_serv
art::ServiceHandle< cheat::BackTrackerService > bt_serv
Definition: ProtoDUNEPizeroAnaTree_module.cc:117
protoana::ProtoDUNEPizeroAnaTree::FillPrimaryPFParticle
void FillPrimaryPFParticle(art::Event const &evt, detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const recob::PFParticle *particle)
Definition: ProtoDUNEPizeroAnaTree_module.cc:808
cache_state.m
m
Definition: cache_state.py:415
geo::Point_t
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
Definition: geo_vectors.h:184
recob::TrackTrajectory::EndDirection
Vector_t EndDirection() const
Returns the direction of the trajectory at the last point.
Definition: TrackTrajectory.h:321
pizero::PiZeroProcess::shower1
const recob::Shower * shower1() const
Definition: PiZeroProcess.h:150
tracks
Definition: tracks.py:1
simb::MCParticle::Vx
double Vx(const int i=0) const
Definition: MCParticle.h:221
recob::Track::ID
int ID() const
Definition: Track.h:198
spacecharge::SpaceCharge::EnableCalSpatialSCE
virtual bool EnableCalSpatialSCE() const =0
Hit.h
Declaration of signal hit object.
protoana::ProtoDUNEBeamlineUtils::GetPID
std::vector< int > GetPID(beam::ProtoDUNEBeamEvent const &beamevt, double nominal_momentum)
Definition: ProtoDUNEBeamlineUtils.cxx:459
protoana::ProtoDUNEPizeroAnaTree::pfpUtil
protoana::ProtoDUNEPFParticleUtils pfpUtil
Definition: ProtoDUNEPizeroAnaTree_module.cc:112
recob::Track::StartMomentum
double StartMomentum() const
Definition: Track.h:143
recob::PFParticle
Hierarchical representation of particle flow.
Definition: PFParticle.h:44
cet::sqlite::min
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:55
pizero::PiZeroProcess::photon2
const simb::MCParticle * photon2() const
Definition: PiZeroProcess.h:147
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
protoana::ProtoDUNEPizeroAnaTree::fShower2_cnn_sc
double fShower2_cnn_sc[NMAXPIZEROS][NMAXHITS]
Definition: ProtoDUNEPizeroAnaTree_module.cc:330
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
protoana::ProtoDUNEPizeroAnaTree::fprimaryDaughterGrandparentEnd
double fprimaryDaughterGrandparentEnd[NMAXDAUGHTERS][3]
Definition: ProtoDUNEPizeroAnaTree_module.cc:247
protoana::ProtoDUNEPizeroAnaTree::fprimaryDaughterStartPosition
double fprimaryDaughterStartPosition[NMAXDAUGHTERS][3]
Definition: ProtoDUNEPizeroAnaTree_module.cc:234
Track.h
Provides recob::Track data product.
simb::MCParticle::Vz
double Vz(const int i=0) const
Definition: MCParticle.h:223
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::ProtoDUNEPizeroAnaTree::fShower1_cnn_sc
double fShower1_cnn_sc[NMAXPIZEROS][NMAXHITS]
Definition: ProtoDUNEPizeroAnaTree_module.cc:329
CRT::PlaneID
detail::Node< FrameID, bool > PlaneID
Definition: CRTID.h:125
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
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
simb::MCParticle::EndX
double EndX() const
Definition: MCParticle.h:226
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::ProtoDUNEPizeroAnaTree::ProtoDUNEPizeroAnaTree
ProtoDUNEPizeroAnaTree(fhicl::ParameterSet const &p)
Definition: ProtoDUNEPizeroAnaTree_module.cc:357
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
protoana::ProtoDUNEPizeroAnaTree::fMCPhoton1_hit_pos
double fMCPhoton1_hit_pos[NMAXPIZEROS][NMAXHITS][3]
Definition: ProtoDUNEPizeroAnaTree_module.cc:310
recob::Shower::ID
int ID() const
Definition: Shower.h:187
art::to_string
std::string to_string(ModuleType const mt)
Definition: ModuleType.h:34
protoana::ProtoDUNEPizeroAnaTree::analyze
void analyze(art::Event const &evt) override
Definition: ProtoDUNEPizeroAnaTree_module.cc:576
protoana::ProtoDUNETruthUtils::GetSharedHits
std::vector< const recob::Hit * > GetSharedHits(detinfo::DetectorClocksData const &clockData, const simb::MCParticle &mcpart, const recob::PFParticle &pfpart, const art::Event &evt, std::string pfparticleModule, bool delta_ray=false) const
Definition: ProtoDUNETruthUtils.cxx:67
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
wirecell.util.wires.graph.parent
def parent(G, child, parent_type)
Definition: graph.py:67
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
fhicl::exception
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
protoana::ProtoDUNEDataUtils::IsBeamTrigger
bool IsBeamTrigger(art::Event const &evt) const
Definition: ProtoDUNEDataUtils.cxx:27
protoana::ProtoDUNEPizeroAnaTree::fMCPhoton1_hit_q
double fMCPhoton1_hit_q[NMAXPIZEROS][NMAXHITS]
Definition: ProtoDUNEPizeroAnaTree_module.cc:308
protoana::ProtoDUNEPizeroAnaTree::fprimaryDaughterGrandparentStart
double fprimaryDaughterGrandparentStart[NMAXDAUGHTERS][3]
Definition: ProtoDUNEPizeroAnaTree_module.cc:245
recob::PFParticle::Daughter
size_t Daughter(size_t idx) const
Returns the ID of the specified daughter.
Definition: PFParticle.h:111
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