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