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