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