mcsXsection_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: mcsXsection //
3 // File: mcsXsection_module.cc //
4 // Ajib Paudel MCS cross section module //
5 // ajib.paudel@phys.ksu.edu //
6 ////////////////////////////////////////////////////////////////////////
7 
8 #include "art/Framework/Core/EDAnalyzer.h"
9 #include "larcorealg/Geometry/Exceptions.h"
10 #include "art/Framework/Core/ModuleMacros.h"
11 #include "art/Framework/Principal/Event.h"
12 #include "art/Framework/Principal/Handle.h"
13 #include "art/Framework/Principal/Run.h"
14 #include "art/Framework/Principal/SubRun.h"
15 #include "art_root_io/TFileService.h"
16 #include "canvas/Utilities/InputTag.h"
17 #include "canvas/Persistency/Common/FindManyP.h"
18 #include "fhiclcpp/ParameterSet.h"
19 #include "messagefacility/MessageLogger/MessageLogger.h"
20 
21 #include "larcore/Geometry/Geometry.h"
22 #include "larcoreobj/SimpleTypesAndConstants/geo_types.h"
23 #include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
24 #include "larsim/MCCheater/BackTrackerService.h"
25 #include "larsim/MCCheater/ParticleInventoryService.h"
26 #include "lardataobj/RecoBase/Hit.h"
27 #include "lardataobj/RecoBase/Track.h"
28 #include "lardataobj/RecoBase/TrackHitMeta.h"
29 #include "lardataobj/RecoBase/SpacePoint.h"
30 #include "lardataobj/RecoBase/Shower.h"
31 #include "lardataobj/RecoBase/PFParticle.h"
32 #include "nusimdata/SimulationBase/MCParticle.h"
33 #include "nusimdata/SimulationBase/MCTruth.h"
34 #include "lardataobj/AnalysisBase/CosmicTag.h"
35 #include "lardataobj/AnalysisBase/T0.h"
36 #include "larreco/RecoAlg/TrackMomentumCalculator.h"
37 #include "lardataobj/AnalysisBase/Calorimetry.h"
38 #include "lardata/ArtDataHelper/MVAReader.h"
39 #include "dune/DuneObj/ProtoDUNEBeamEvent.h"
40 #include "larreco/RecoAlg/TrajectoryMCSFitter.h"
41 //#include "dune/EventGenerator/ProtoDUNEbeamDataProducts/ProtoDUNEbeamsim.h"
42 //#include "dune/EventGenerator/ProtoDUNEbeamDataProducts/ProtoDUNEBeamInstrument.h"
43 //#include "dune/Protodune/singlephase/DataUtils/ProtoDUNETrackUtils.h"
44 //#include "dune/Protodune/singlephase/DataUtils/ProtoDUNEShowerUtils.h"
45 //#include "dune/Protodune/singlephase/DataUtils/ProtoDUNETruthUtils.h"
46 //#include "dune/Protodune/singlephase/DataUtils/ProtoDUNEPFParticleUtils.h"
47 //#include "dune/Protodune/singlephase/DataUtils/ProtoDUNEDataUtils.h"
48 #include "protoduneana/Utilities/ProtoDUNETrackUtils.h"
49 #include "protoduneana/Utilities/ProtoDUNEShowerUtils.h"
50 #include "protoduneana/Utilities/ProtoDUNETruthUtils.h"
51 #include "protoduneana/Utilities/ProtoDUNEPFParticleUtils.h"
52 #include "dune/Protodune/singlephase/DataUtils/ProtoDUNEDataUtils.h"
53 #include "protoduneana/Utilities/ProtoDUNEBeamlineUtils.h"
54 #include "protoduneana/Utilities/ProtoDUNEBeamCuts.h"
55 #include "larsim/MCCheater/ParticleInventoryService.h"
56 
57 
58 //#g4reweight
59 #include "geant4reweight/src/ReweightBase/G4ReweighterFactory.hh"
60 #include "geant4reweight/src/ReweightBase/G4Reweighter.hh"
61 #include "geant4reweight/src/ReweightBase/G4ReweightTraj.hh"
62 #include "geant4reweight/src/PropBase/G4ReweightParameterMaker.hh"
63 #include "geant4reweight/src/ReweightBase/G4MultiReweighter.hh"
64 #include "protoduneana/Utilities/G4ReweightUtils.h"
65 
66 
67 // ROOT includes
68 #include "TTree.h"
69 #include "TFile.h"
70 #include "TString.h"
71 #include "TTimeStamp.h"
72 #include "TH2.h"
73 // C++ Includes
74 #include <stdio.h>
75 #include <stdlib.h>
76 #include <string>
77 #include <sstream>
78 #include <cmath>
79 #include <algorithm>
80 #include <iostream>
81 #include <vector>
82 #include "TComplex.h"
83 #include <fstream>
84 #include "TPaveStats.h"
85 #include <iostream>
86 #include <string>
87 #include "math.h"
88 #include <iterator>
89 #include <map>
90 #include <tuple>
91 
92 // Maximum number of beam particles to save
93 const int NMAXDAUGTHERS = 30;
94 double m_pi=0.1395;//GeV/c^2
95 double prim_energy=0.0;
96 //double geant_energy=0;
97 
98 
99 float soln(float x1, float x2, float x3, float x4, float y1, float y2, float y3, float y4, float result[2]){
100  if(x3==x4){
101  result[0]=x3;
102  result[1]=y3;
103  }
104  else if((x1-x2)*(y3-y4)-(y1-y2)*(x3-x4)==0){
105  result[0]=x3;
106  result[1]=y3;
107  }
108  else{
109  result[0]=((x1*y2-y1*x2)*(x3-x4)-(x1-x2)*(x3*y4-y3*x4))/((x1-x2)*(y3-y4)-(y1-y2)*(x3-x4));
110  result[1]=((x1*y2-y1*x2)*(y3-y4)-(y1-y2)*(x3*y4-y3*x4))/((x1-x2)*(y3-y4)-(y1-y2)*(x3-x4));
111  }
112  return *result;
113 }
114 double angle2d(double x0, double y0, double x1, double y1, double x2, double y2){
115  double theta=((x1-x0)*(x2-x1)+(y1-y0)*(y2-y1))/sqrt(((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0))*((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)));
116  return acos(theta);
117 
118 }
119 
120 
121 double theta12(double x1, double x2, double y1, double y2, double z1, double z2,double x1p, double x2p, double y1p, double y2p, double z1p, double z2p){
122  double numer=(x2-x1)*(x2p-x1p)+(y2-y1)*(y2p-y1p)+(z2-z1)*(z2p-z1p);
123  double den1=(x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)+(z2-z1)*(z2-z1);
124  double den2=(x2p-x1p)*(x2p-x1p)+(y2p-y1p)*(y2p-y1p)+(z2p-z1p)*(z2p-z1p);
125  return 180/3.14*acos(numer/sqrt(den1*den2));
126 }
127 
128 
129 namespace protoana {
130  class mcsXsection;
131 }
132 
133 using protoana::G4ReweightUtils::CreateRWTraj;
134 using protoana::G4ReweightUtils::CreateNRWTrajs;
135 using protoana::G4ReweightUtils::GetNTrajPMSigmaWeights;
136 using protoana::G4ReweightUtils::GetNTrajWeightFromSetPars;
137 
138 class protoana::mcsXsection : public art::EDAnalyzer {
139 public:
140 
141  explicit mcsXsection(fhicl::ParameterSet const & p);
142 
143  mcsXsection(mcsXsection const &) = delete;
144  mcsXsection(mcsXsection &&) = delete;
145  mcsXsection & operator = (mcsXsection const &) = delete;
146  mcsXsection & operator = (mcsXsection &&) = delete;
147 
148  virtual void beginJob() override;
149  virtual void endJob() override;
150 
151  // Required functions.
152  void analyze(art::Event const & evt) override;
153 
154 private:
155  const art::InputTag fTrackModuleLabel;
156  const art::InputTag fBeamModuleLabel;
157  // Helper utility functions
158  protoana::ProtoDUNEBeamCuts beam_cuts;
159  protoana::ProtoDUNEBeamlineUtils fBeamlineUtils;
160  protoana::ProtoDUNEDataUtils dataUtil;
161 
162  protoana::ProtoDUNEPFParticleUtils pfpUtil;
163  protoana::ProtoDUNETrackUtils trackUtil;
164  protoana::ProtoDUNETruthUtils truthUtil;
165 
166 
167  // Initialise tree variables
168  void Initialise();
169 
170  // Fill cosmics tree
171  void FillCosmicsTree(art::Event const & evt, std::string pfParticleTag);
172 
173  // fcl parameters
174  //const art::InputTag fNNetModuleLabel; //label of the module used for CNN tagging
175 
176 
177  std::string fCalorimetryTag;
178  std::string fTrackerTag;
179  std::string fShowerTag;
180  std::string fPFParticleTag;
181  std::string fGeneratorTag;
182  //bool fVerbose;
183 
184  geo::GeometryCore const * fGeometry;
185 
186  TTree *fPandoraBeam;
187  //TTree *fPandoraCosmics;
188  // Tree variables
189  int fRun;
190  int fSubRun;
191  int fevent;
192  double fTimeStamp;
193 
194 
195  double fbeamtrackMomentum;
196  double fbeamtrackP[3]; //Px/Py/Pz
197  double fbeamtrackEnergy;
198  double fbeamtrackPos[3];
199  double fbeamtrackDir[3];
200  double fbeamtrackTime;
201  int fbeamtrackPdg;
202  int fbeamtrackID;
203 
204  //int NumberBeamTrajectoryPoints;
205  std::vector<float> beamtrk_x;
206  std::vector<float> beamtrk_y;
207  std::vector<float> beamtrk_z;
208  std::vector<int> beamtrk_z_wire;
209  std::vector<int> beamtrk_z_tpc;
210 
211  std::vector<float> beamtrk_Px;
212  std::vector<float> beamtrk_Py;
213  std::vector<float> beamtrk_Pz;
214  std::vector<float> beamtrk_Eng;
215  std::vector<float> peakT_2;
216  std::vector<float> hitz_2;
217  // std::vector<float> inelscore;
218  //std::vector<float> elscore;
219  //std::vector<float> nonescore;
220 
221  /////////
222  //beam info from spectrometer
223  std::vector<double> beamMomentum_spec;
224  std::vector<double> beamPosx_spec;
225  std::vector<double> beamPosy_spec;
226  std::vector<double> beamPosz_spec;
227 
228  std::vector<double> beamDirx_spec;
229  std::vector<double> beamDiry_spec;
230  std::vector<double> beamDirz_spec;
231 
232 
233 
234 
235 
236 
237  std::vector<float> int_2;
238  std::vector<float> hitz_1;
239  std::vector<float> int_1;
240  std::vector<float> hitz_0;
241  std::vector<float> int_0;
242 
243  // Reconstructed tracks/showers information
244  double fvertex[3];
245  double fsecvertex[3];
246 
247  int fisprimarytrack;
248  int fisprimaryshower;
249  double fprimaryBDTScore;
250  int fprimaryNHits;
251  double fprimaryTheta;
252  double fprimaryPhi;
253  double fprimaryLength;
254  double fprimaryMomentum;
255  double fprimaryEndMomentum;
256  double fprimaryEndPosition[3];
257  double fprimaryStartPosition[3];
258  double fprimaryEndDirection[3];
259  double fprimaryStartDirection[3];
260  double fprimaryOpeningAngle;
261  int fprimaryShowerBestPlane;
262  double fprimaryShowerEnergy;
263  double fprimaryShowerMIPEnergy;
264  double fprimaryShowerdEdx;
265  double fprimaryMomentumByRangeProton;
266  double fprimaryMomentumByRangeMuon;
267  double fprimaryKineticEnergy[3];
268  double fprimaryRange[3];
269  int fprimaryID;
270 
271 
272 
273 
274  //*********************truth information*******************************************//
275  int fprimary_truth_TrackId;
276  int fprimary_truth_Pdg;
277  int ftruthpdg;
278  double fprimary_truth_StartPosition[4];
279  // double fprimary_truth_Startwiretime[4];
280  double fprimary_truth_EndPosition[4];
281  std::string fprimary_truth_EndProcess;
282  std::string truth_last_process;
283  double fprimary_truth_P;
284  double fprimary_truth_Momentum[4];
285  double fprimary_truth_EndMomentum[4];
286  double fprimary_truth_Pt;
287  double fprimary_truth_Mass;
288  double fprimary_truth_Theta;
289  double fprimary_truth_Phi;
290  int fprimary_truth_Process;
291  int fprimary_truth_Isbeammatched;
292  int fprimary_truth_NDaughters;
293  double fprimary_truth_tracklength;
294 
295 
296  //interaction point details
297  std::vector<double> interactionX;
298  std::vector<double> interactionY;
299  std::vector<double> interactionZ;
300  std::vector<double> interactionU;
301  std::vector<double> interactionV;
302  std::vector<double> interactionW;
303  std::vector<double> interactionT;
304 
305 
306 
307  std::vector<std::string> interactionProcesslist;
308  std::vector<double> interactionAngles;
309  std::vector<double> interactionAnglesUT;
310  std::vector<double> interactionAnglesVT;
311  std::vector<double> interactionAnglesZT;
312  std::vector<double> Zintersection;
313  std::vector<double> Zintersection1;
314  std::vector<double> timeintersection;
315  std::vector<double> timeintersection1;
316  std::vector<double> deltaZint;
317  std::vector<double> deltatimeint;
318  //calo info
319  std::vector< std::vector<double> > primtrk_dqdx;
320  std::vector< std::vector<double> > primtrk_resrange;
321  std::vector< std::vector<double> > primtrk_dedx;
322  std::vector<double> primtrk_range;
323  std::vector< std::vector<double> > primtrk_hitx;
324  std::vector< std::vector<double> > primtrk_hity;
325  std::vector< std::vector<double> > primtrk_hitz;
326  std::vector< std::vector<int> > primtrk_hitz_wire;
327  std::vector< std::vector<int> > primtrk_hitz_tpc;
328 
329  std::vector< std::vector<double> > primtrk_pitch;
330  std::vector<std::vector<double> > primtrk_truth_Z;
331  std::vector<std::vector<int> > primtrk_truth_Z_wire;
332  std::vector<std::vector<int> > primtrk_truth_Z_tpc;
333 
334 
335  std::vector<std::vector<double> > primtrk_truth_Eng;
336  std::vector<std::vector<double> > primtrk_truth_trkide;
337  std::vector<std::vector<float> > wireno_0;
338  std::vector<std::vector<float> > wireno_1;
339  std::vector<std::vector<float> > wireno_2;
340  std::vector<std::vector<float> > peakTime_0;
341  std::vector<std::vector<float> > peakTime_1;
342  std::vector<std::vector<float> > peakTime_2;
343  std::vector<std::vector<float> > dq_0;
344  std::vector<std::vector<float> > dq_1;
345  std::vector<std::vector<float> > dq_2;
346 
347  ////Michel tagging info
348 
349  std::vector< std::vector<int> > endhitssecondary;
350  std::vector< std::vector<double> > secondarystartx;
351  std::vector< std::vector<double> > secondaryendx;
352  std::vector< std::vector<double> > secondarystarty;
353  std::vector< std::vector<double> > secondaryendy;
354  std::vector< std::vector<double> > secondarystartz;
355  std::vector< std::vector<double> > secondaryendz;
356  std::vector< std::vector<double> > dQmichel;
357  std::vector< std::vector<double> > dQtrackend;
358  std::vector< std::vector<double> > dQtrackbegin;
359  std::vector< std::vector<double> > primsectheta;
360  std::vector< std::vector<double> > tracklengthsecondary;
361  std::vector< std::vector<int> > MtrackID;
362  ////CNN for michel tagging
363  std::vector< std::vector<double> > trackscore;
364  std::vector< std::vector<double> > emscore;
365  std::vector< std::vector<double> > michelscore;
366  std::vector< std::vector<double> > nonescore;
367 
368 
369 
370  // Daughters from primary
371  int fNDAUGHTERS;
372  int fisdaughtertrack[NMAXDAUGTHERS];
373  int fisdaughtershower[NMAXDAUGTHERS];
374  int fdaughterNHits[NMAXDAUGTHERS];
375  double fdaughterTheta[NMAXDAUGTHERS];
376  double fdaughterPhi[NMAXDAUGTHERS];
377  double fdaughterLength[NMAXDAUGTHERS];
378  double fdaughterMomentum[NMAXDAUGTHERS];
379  double fdaughterEndMomentum[NMAXDAUGTHERS];
380  double fdaughterEndPosition[NMAXDAUGTHERS][3];
381  double fdaughterStartPosition[NMAXDAUGTHERS][3];
382  double fdaughterEndDirection[NMAXDAUGTHERS][3];
383  double fdaughterStartDirection[NMAXDAUGTHERS][3];
384  double fdaughterOpeningAngle[NMAXDAUGTHERS];
385  double fdaughterShowerEnergy[NMAXDAUGTHERS];
386  double fdaughterShowerMIPEnergy[NMAXDAUGTHERS];
387  double fdaughterShowerdEdx[NMAXDAUGTHERS];
388  int fdaughterShowerBestPlane[NMAXDAUGTHERS];
389  double fdaughterMomentumByRangeProton[NMAXDAUGTHERS];
390  double fdaughterMomentumByRangeMuon[NMAXDAUGTHERS];
391  double fdaughterKineticEnergy[NMAXDAUGTHERS][3];
392  double fdaughterRange[NMAXDAUGTHERS][3];
393  int fdaughterID[NMAXDAUGTHERS];
394  //double fdaughterT0[NMAXDAUGTHERS];
395 
396  int fdaughter_truth_TrackId[NMAXDAUGTHERS];
397  int fdaughter_truth_Pdg[NMAXDAUGTHERS];
398  double fdaughter_truth_StartPosition[NMAXDAUGTHERS][4];
399  double fdaughter_truth_EndPosition[NMAXDAUGTHERS][4];
400  double fdaughter_truth_P[NMAXDAUGTHERS];
401  double fdaughter_truth_Momentum[NMAXDAUGTHERS][4];
402  double fdaughter_truth_EndMomentum[NMAXDAUGTHERS][4];
403  double fdaughter_truth_Pt[NMAXDAUGTHERS];
404  double fdaughter_truth_Mass[NMAXDAUGTHERS];
405  double fdaughter_truth_Theta[NMAXDAUGTHERS];
406  double fdaughter_truth_Phi[NMAXDAUGTHERS];
407  int fdaughter_truth_Process[NMAXDAUGTHERS];
408 
409  double minX = -360.0;
410  double maxX = 360.0;
411  double minY =0.0;
412  double maxY = 600.0;
413  double minZ = 0.0;
414  double maxZ = 695.0;
415 
416 
417  int true_beam_PDG;
418  int true_beam_ID;
419 
420  int n_g4rw;
421  double g4rw_react_st;
422  double g4rw_d_react;
423  double g4rw_elast_st;
424  double g4rw_d_elast;
425  std::vector<double> g4rw_react;
426  std::vector<double> g4rw_elast;
427  std::vector<double> g4rw_set_weights;
428  //Geant4Reweight stuff
429  int RW_PDG;
430  TFile FracsFile, XSecFile;
431  std::vector<fhicl::ParameterSet> ParSet;
432  G4ReweightParameterMaker ParMaker;
433  G4ReweightManager RWManager;
434  G4MultiReweighter MultiRW;
435  //G4ReweighterFactory RWFactory;
436  //G4Reweighter * theRW;
437 
438  //bool fVerbose;
439 
440 
441 };
442 
443 
444 protoana::mcsXsection::mcsXsection(fhicl::ParameterSet const & p)
445  :
446  EDAnalyzer(p),
447 
448  // fNNetModuleLabel(p.get<art::InputTag>("NNetModuleLabel")),
449  fBeamModuleLabel(p.get< art::InputTag >("BeamModuleLabel")),
450  beam_cuts(p.get<fhicl::ParameterSet>("BeamCuts")),
451  fBeamlineUtils(p.get<fhicl::ParameterSet>("BeamlineUtils")),
452  dataUtil(p.get<fhicl::ParameterSet>("DataUtils")),
453  fCalorimetryTag(p.get<std::string>("CalorimetryTag")),
454  fTrackerTag(p.get<std::string>("TrackerTag")),
455  fShowerTag(p.get<std::string>("ShowerTag")),
456  fPFParticleTag(p.get<std::string>("PFParticleTag")),
457  fGeneratorTag(p.get<std::string>("GeneratorTag")),
458  RW_PDG(p.get<int>("RW_PDG")),
459  FracsFile( (p.get< std::string >( "FracsFile" )).c_str(), "OPEN" ),
460  XSecFile( (p.get< std::string >( "XSecFile" )).c_str(), "OPEN"),
461  ParSet(p.get<std::vector<fhicl::ParameterSet>>("ParameterSet")),
462  ParMaker(ParSet, RW_PDG),
463  RWManager({p.get<fhicl::ParameterSet>("Material")}),
464  MultiRW(RW_PDG, FracsFile, ParSet,
465  p.get<fhicl::ParameterSet>("Material"),
466  &RWManager) {
467  std::cout << "done" << std::endl;
468 }
469 
470 /* fVerbose(p.get<bool>("Verbose"))
471  {
472 
473  }
474  */
475  void protoana::mcsXsection::beginJob(){
476 
477  art::ServiceHandle<art::TFileService> tfs;
478 
479  fPandoraBeam = tfs->make<TTree>("PandoraBeam", "Beam events reconstructed with Pandora");
480  fPandoraBeam->Branch("run", &fRun, "run/I");
481  fPandoraBeam->Branch("subrun", &fSubRun, "subrun/I");
482  fPandoraBeam->Branch("event", &fevent, "event/I");
483  fPandoraBeam->Branch("timestamp", &fTimeStamp, "timestamp/D");
484  fPandoraBeam->Branch("beamtrackMomentum", &fbeamtrackMomentum, "beamtrackMomentum/D");
485  fPandoraBeam->Branch("beamtrackP", &fbeamtrackP, "beamtrackP[3]/D");
486  fPandoraBeam->Branch("beamtrackEnergy", &fbeamtrackEnergy, "beamtrackEnergy/D");
487  fPandoraBeam->Branch("beamtrackPos", &fbeamtrackPos, "beamtrackPos[3]/D");
488  fPandoraBeam->Branch("beamtrackDir", &fbeamtrackDir, "beamtrackDir[3]/D");
489  fPandoraBeam->Branch("beamtrackTime", &fbeamtrackTime, "beamtrackTime/D");
490  fPandoraBeam->Branch("beamtrackPdg", &fbeamtrackPdg, "beamtrackPdg/I");
491  fPandoraBeam->Branch("beamtrackID", &fbeamtrackID, "beamtrackID/I");
492 
493  //fPandoraBeam->Branch("NumberBeamTrajectoryPoints", &NumberBeamTrajectoryPoints, "NumberBeamTrajectoryPoints/I");
494  fPandoraBeam->Branch("beamtrk_x",&beamtrk_x);
495  fPandoraBeam->Branch("beamtrk_y",&beamtrk_y);
496  fPandoraBeam->Branch("beamtrk_z",&beamtrk_z);
497  fPandoraBeam->Branch("beamtrk_z_wire",&beamtrk_z_wire);
498  fPandoraBeam->Branch("beamtrk_z_tpc",&beamtrk_z_tpc);
499 
500  fPandoraBeam->Branch("beamtrk_Px",&beamtrk_Px);
501  fPandoraBeam->Branch("beamtrk_Py",&beamtrk_Py);
502  fPandoraBeam->Branch("beamtrk_Pz",&beamtrk_Pz);
503  fPandoraBeam->Branch("beamtrk_Eng",&beamtrk_Eng);
504  ////.............spectrometer parameters
505  fPandoraBeam->Branch("beamMomentum_spec",&beamMomentum_spec);
506  fPandoraBeam->Branch("beamPosx_spec",&beamPosx_spec);
507  fPandoraBeam->Branch("beamPosy_spec",&beamPosy_spec);
508  fPandoraBeam->Branch("beamPosz_spec",&beamPosz_spec);
509  fPandoraBeam->Branch("beamDirx_spec",&beamDirx_spec);
510  fPandoraBeam->Branch("beamDiry_spec",&beamDiry_spec);
511  fPandoraBeam->Branch("beamDirz_spec",&beamDirz_spec);
512 
513 
514 
515 
516  fPandoraBeam->Branch("vertex", &fvertex, "vertex[3]/D");
517  fPandoraBeam->Branch("secvertex", &fsecvertex, "secvertex[3]/D");
518  fPandoraBeam->Branch("isprimarytrack", &fisprimarytrack, "isprimarytrack/I");
519  fPandoraBeam->Branch("isprimaryshower", &fisprimaryshower, "isprimaryshower/I");
520  fPandoraBeam->Branch("primaryBDTScore", &fprimaryBDTScore, "primaryBDTScore/D");
521  fPandoraBeam->Branch("primaryNHits", &fprimaryNHits, "fprimaryNHits/I");
522  fPandoraBeam->Branch("primaryTheta", &fprimaryTheta, "primaryTheta/D");
523  fPandoraBeam->Branch("primaryPhi", &fprimaryPhi, "primaryPhi/D");
524  fPandoraBeam->Branch("primaryLength", &fprimaryLength, "primaryLength/D");
525  fPandoraBeam->Branch("primaryMomentum", &fprimaryMomentum, "primaryMomentum/D");
526  fPandoraBeam->Branch("primaryEndMomentum", &fprimaryEndMomentum, "primaryEndMomentum/D");
527  fPandoraBeam->Branch("primaryEndPosition", &fprimaryEndPosition, "primaryEndPosition[3]/D");
528  fPandoraBeam->Branch("primaryStartPosition", &fprimaryStartPosition, "primaryStartPosition[3]/D");
529  fPandoraBeam->Branch("primaryEndDirection", &fprimaryEndDirection, "primaryEndDirection[3]/D");
530  fPandoraBeam->Branch("primaryStartDirection", &fprimaryStartDirection, "primaryStartDirection[3]/D");
531  fPandoraBeam->Branch("primaryOpeningAngle", &fprimaryOpeningAngle, "primaryOpeningAngle/D");
532  fPandoraBeam->Branch("primaryID", &fprimaryID, "primaryID/I");
533  fPandoraBeam->Branch("primaryShowerBestPlane", &fprimaryShowerBestPlane, "primaryShowerBestPlane/I");
534  fPandoraBeam->Branch("primaryShowerEnergy", &fprimaryShowerEnergy, "primaryShowerEnergy/I");
535  fPandoraBeam->Branch("primaryShowerMIPEnergy", &fprimaryShowerMIPEnergy, "primaryShowerMIPEnergy/I");
536  fPandoraBeam->Branch("primaryShowerdEdx", &fprimaryShowerdEdx, "primaryShowerdEdx/I");
537  fPandoraBeam->Branch("primaryMomentumByRangeProton", &fprimaryMomentumByRangeProton, "primaryMomentumByRangeProton/D");
538  fPandoraBeam->Branch("primaryMomentumByRangeMuon", &fprimaryMomentumByRangeMuon, "primaryMomentumByRangeMuon/D");
539  fPandoraBeam->Branch("primaryKineticEnergy", &fprimaryKineticEnergy, "primaryKineticEnergy[3]/D");
540  fPandoraBeam->Branch("primaryRange", &fprimaryRange, "primaryRange[3]/D");
541 
542 
543  fPandoraBeam->Branch("primary_truth_TrackId", &fprimary_truth_TrackId, "primary_truth_TrackId/I");
544  fPandoraBeam->Branch("primary_truth_Pdg", &fprimary_truth_Pdg, "primary_truth_Pdg/I");
545  fPandoraBeam->Branch("truthpdg", &ftruthpdg, "truthpdg/I");
546 
547 
548  fPandoraBeam->Branch("primary_truth_StartPosition", &fprimary_truth_StartPosition, "primary_truth_StartPosition[4]/D");
549  // fPandoraBeam->Branch("primary_truth_Startwiretime", &fprimary_truth_Startwiretime, "primary_truth_Startwiretime[4]/D");
550  fPandoraBeam->Branch("primary_truth_EndPosition", &fprimary_truth_EndPosition, "primary_truth_EndPosition[4]/D");
551  fPandoraBeam->Branch("primary_truth_Momentum", &fprimary_truth_Momentum, "primary_truth_Momentum[4]/D");
552  fPandoraBeam->Branch("primary_truth_EndMomentum", &fprimary_truth_EndMomentum, "primary_truth_EndMomentum[4]/D");
553  fPandoraBeam->Branch("primary_truth_P", &fprimary_truth_P, "primary_truth_P/D");
554  fPandoraBeam->Branch("primary_truth_Pt", &fprimary_truth_Pt, "primary_truth_Pt/D");
555  fPandoraBeam->Branch("primary_truth_Mass", &fprimary_truth_Mass, "primary_truth_Mass/D");
556  fPandoraBeam->Branch("primary_truth_Theta", &fprimary_truth_Theta, "primary_truth_Theta/D");
557  fPandoraBeam->Branch("primary_truth_Phi", &fprimary_truth_Phi, "primary_truth_Phi/D");
558  fPandoraBeam->Branch("primary_truth_Process", &fprimary_truth_Process, "primary_truth_Process/I");
559  fPandoraBeam->Branch("primary_truth_Isbeammatched", &fprimary_truth_Isbeammatched, "primary_truth_Isbeammatched/I");
560  fPandoraBeam->Branch("primary_truth_NDaughters", &fprimary_truth_NDaughters, "primary_truth_NDaughters/I");
561  fPandoraBeam->Branch("primary_truth_EndProcess", &fprimary_truth_EndProcess);
562  fPandoraBeam->Branch("truth_last_process", &truth_last_process);
563  fPandoraBeam->Branch("primary_truth_tracklength", &fprimary_truth_tracklength, "primary_truth_tracklength/D");
564 
565  fPandoraBeam->Branch("interactionX",&interactionX);
566  fPandoraBeam->Branch("interactionY",&interactionY);
567  fPandoraBeam->Branch("interactionZ",&interactionZ);
568  fPandoraBeam->Branch("interactionU",&interactionU);
569  fPandoraBeam->Branch("interactionV",&interactionV);
570  fPandoraBeam->Branch("interactionW",&interactionW);
571  fPandoraBeam->Branch("interactionT",&interactionT);
572  fPandoraBeam->Branch("interactionProcesslist",&interactionProcesslist);
573  fPandoraBeam->Branch("interactionAngles",&interactionAngles);
574  fPandoraBeam->Branch("interactionAnglesUT",&interactionAnglesUT);
575  fPandoraBeam->Branch("interactionAnglesVT",&interactionAnglesVT);
576  fPandoraBeam->Branch("interactionAnglesZT",&interactionAnglesZT);
577 
578  fPandoraBeam->Branch("Zintersection",&Zintersection);
579  fPandoraBeam->Branch("timeintersection",&timeintersection);
580  fPandoraBeam->Branch("Zintersection1",&Zintersection1);
581  fPandoraBeam->Branch("timeintersection1",&timeintersection1);
582  fPandoraBeam->Branch("deltaZint",&deltaZint);
583  fPandoraBeam->Branch("deltatimeint",&deltatimeint);
584  fPandoraBeam->Branch("peakT_2",&peakT_2);
585  fPandoraBeam->Branch("hitz_2",&hitz_2);
586  // fPandoraBeam->Branch("inelscore",&inelscore);
587  //fPandoraBeam->Branch("elscore",&elscore);
588  //fPandoraBeam->Branch("nonescore",&nonescore);
589 
590  fPandoraBeam->Branch("int_2",&int_2);
591  fPandoraBeam->Branch("hitz_1",&hitz_1);
592  fPandoraBeam->Branch("int_1",&int_1);
593  fPandoraBeam->Branch("hitz_0",&hitz_0);
594  fPandoraBeam->Branch("int_0",&int_0);
595  fPandoraBeam->Branch("wireno_0",&wireno_0);
596  fPandoraBeam->Branch("wireno_1",&wireno_1);
597  fPandoraBeam->Branch("wireno_2",&wireno_2);
598  fPandoraBeam->Branch("peakTime_0",&peakTime_0);
599  fPandoraBeam->Branch("peakTime_1",&peakTime_1);
600  fPandoraBeam->Branch("peakTime_2",&peakTime_2);
601  fPandoraBeam->Branch("dq_0",&dq_0);
602  fPandoraBeam->Branch("dq_1",&dq_1);
603  fPandoraBeam->Branch("dq_2",&dq_2);
604  //Michel tagging
605  fPandoraBeam->Branch("endhitssecondary",&endhitssecondary);
606  fPandoraBeam->Branch("secondarystartx",&secondarystartx);
607  fPandoraBeam->Branch("secondarystarty",&secondarystarty);
608  fPandoraBeam->Branch("secondarystartz",&secondarystartz);
609  fPandoraBeam->Branch("secondaryendx",&secondaryendx);
610  fPandoraBeam->Branch("secondaryendy",&secondaryendy);
611  fPandoraBeam->Branch("secondaryendz",&secondaryendz);
612  fPandoraBeam->Branch("dQmichel",&dQmichel);
613  fPandoraBeam->Branch("primsectheta",&primsectheta);
614  fPandoraBeam->Branch("dQtrackbegin",&dQtrackbegin);
615  fPandoraBeam->Branch("dQtrackend",&dQtrackend);
616  fPandoraBeam->Branch("tracklengthsecondary",&tracklengthsecondary);
617  fPandoraBeam->Branch("MtrackID",&MtrackID);
618  fPandoraBeam->Branch("trackscore",&trackscore);
619  fPandoraBeam->Branch("emscore",&emscore);
620  fPandoraBeam->Branch("michelscore",&michelscore);
621  fPandoraBeam->Branch("nonescore",&nonescore);
622  /////////
623 
624 
625 
626  fPandoraBeam->Branch("NDAUGHTERS", &fNDAUGHTERS, "NDAUGHTERS/I");
627  fPandoraBeam->Branch("isdaughtertrack", &fisdaughtertrack, "isdaughtertrack[NDAUGHTERS]/I");
628  fPandoraBeam->Branch("isdaughtershower", &fisdaughtershower, "isdaughtershower[NDAUGHTERS]/I");
629  fPandoraBeam->Branch("daughterNHits", &fdaughterNHits, "daughterNHits[NDAUGHTERS]/I");
630  fPandoraBeam->Branch("daughterTheta", &fdaughterTheta, "daughterTheta[NDAUGHTERS]/D");
631  fPandoraBeam->Branch("daughterPhi", &fdaughterPhi, "daughterPhi[NDAUGHTERS]/D");
632  fPandoraBeam->Branch("daughterLength", &fdaughterLength, "daughterLength[NDAUGHTERS]/D");
633  fPandoraBeam->Branch("daughterMomentum", &fdaughterMomentum, "daughterMomentum[NDAUGHTERS]/D");
634  fPandoraBeam->Branch("daughterEndMomentum", &fdaughterEndMomentum, "daughterEndMomentum[NDAUGHTERS]/D");
635  fPandoraBeam->Branch("daughterEndPosition", &fdaughterEndPosition, "daughterEndPosition[NDAUGHTERS][3]/D");
636  fPandoraBeam->Branch("daughterStartPosition", &fdaughterStartPosition, "daughterStartPosition[NDAUGHTERS][3]/D");
637  fPandoraBeam->Branch("daughterStartDirection", &fdaughterStartDirection, "daughterStartDirection[NDAUGHTERS][3]/D");
638  fPandoraBeam->Branch("daughterEndDirection", &fdaughterEndDirection, "daughterEndDirection[NDAUGHTERS][3]/D");
639  fPandoraBeam->Branch("daughterOpeningAngle", &fdaughterOpeningAngle, "daughterOpeningAngle[NDAUGHTERS]/D");
640  fPandoraBeam->Branch("daughterShowerBestPlane", &fdaughterShowerBestPlane, "daughterShowerBestPlane[NDAUGHTERS]/D");
641  fPandoraBeam->Branch("daughterShowerEnergy", &fdaughterShowerEnergy, "daughterShowerEnergy[NDAUGHTERS]/D");
642  fPandoraBeam->Branch("daughterShowerMIPEnergy", &fdaughterShowerMIPEnergy, "daughterShowerMIPEnergy[NDAUGHTERS]/D");
643  fPandoraBeam->Branch("daughterShowerdEdx", &fdaughterShowerdEdx, "daughterShowerdEdx[NDAUGHTERS]/D");
644  fPandoraBeam->Branch("daughterMomentumByRangeProton", &fdaughterMomentumByRangeProton,"daughterMomentumByRangeProton[NDAUGHTERS]/D");
645  fPandoraBeam->Branch("daughterMomentumByRangeMuon", &fdaughterMomentumByRangeMuon, "daughterMomentumByRangeMuon[NDAUGHTERS]/D");
646  fPandoraBeam->Branch("daughterKineticEnergy", &fdaughterKineticEnergy, "daughterKineticEnergy[NDAUGHTERS][3]/D");
647  fPandoraBeam->Branch("daughterRange", &fdaughterRange, "daughterRange[NDAUGHTERS][3]/D");
648  fPandoraBeam->Branch("daughterID", &fdaughterID, "daughterID[NDAUGHTERS]/I");
649  // fPandoraBeam->Branch("daughterT0", &fdaughterT0, "daughterT0[NDAUGHTERS]/D");
650 
651  fPandoraBeam->Branch("daughter_truth_TrackId", &fdaughter_truth_TrackId, "daughter_truth_TrackId[NDAUGHTERS]/I");
652  fPandoraBeam->Branch("daughter_truth_Pdg", &fdaughter_truth_Pdg, "daughter_truth_Pdg[NDAUGHTERS]/I");
653  fPandoraBeam->Branch("daughter_truth_StartPosition", &fdaughter_truth_StartPosition, "daughter_truth_StartPosition[NDAUGHTERS][4]/D");
654  fPandoraBeam->Branch("daughter_truth_EndPosition", &fdaughter_truth_EndPosition, "daughter_truth_EndPosition[NDAUGHTERS][4]/D");
655  fPandoraBeam->Branch("daughter_truth_Momentum", &fdaughter_truth_Momentum, "daughter_truth_Momentum[NDAUGHTERS][4]/D");
656  fPandoraBeam->Branch("daughter_truth_EndMomentum", &fdaughter_truth_EndMomentum, "daughter_truth_EndMomentum[NDAUGHTERS][4]/D");
657  fPandoraBeam->Branch("daughter_truth_P", &fdaughter_truth_P, "daughter_truth_P[NDAUGHTERS]/D");
658  fPandoraBeam->Branch("daughter_truth_Pt", &fdaughter_truth_Pt, "daughter_truth_Pt[NDAUGHTERS]/D");
659  fPandoraBeam->Branch("daughter_truth_Mass", &fdaughter_truth_Mass, "daughter_truth_Mass[NDAUGHTERS]/D");
660  fPandoraBeam->Branch("daughter_truth_Theta", &fdaughter_truth_Theta, "daughter_truth_Theta[NDAUGHTERS]/D");
661  fPandoraBeam->Branch("daughter_truth_Phi", &fdaughter_truth_Phi, "daughter_truth_Phi[NDAUGHTERS]/D");
662  fPandoraBeam->Branch("daughter_truth_Process", &fdaughter_truth_Process, "daughter_truth_Process[NDAUGHTERS]/I");
663 
664 
665  // hKE_truth_reco = tfs->make<TH2D>("hKE_truth_reco","KE truth vs reco;true KE in MeV;reconstructed KE in MeV" , 3000, -100, 1400 ,3000, -100, 1400);
666  // hEndZ_truth_reco = tfs->make<TH2D>("hEndZ_truth_reco","EndZ truth vs reco;true EndZ (cm);reconstructed EndZ(cm)" , 695, 0, 695 ,695,0,695);
667 
668  fPandoraBeam->Branch("primtrk_dqdx",&primtrk_dqdx);
669  fPandoraBeam->Branch("primtrk_dedx",&primtrk_dedx);
670  fPandoraBeam->Branch("primtrk_resrange",&primtrk_resrange);
671  fPandoraBeam->Branch("primtrk_range",&primtrk_range);
672  fPandoraBeam->Branch("primtrk_hitx",&primtrk_hitx);
673  fPandoraBeam->Branch("primtrk_hity",&primtrk_hity);
674  fPandoraBeam->Branch("primtrk_hitz",&primtrk_hitz);
675  fPandoraBeam->Branch("primtrk_hitz_wire",&primtrk_hitz_wire);
676  fPandoraBeam->Branch("primtrk_hitz_tpc",&primtrk_hitz_tpc);
677  fPandoraBeam->Branch("primtrk_pitch",&primtrk_pitch);
678 
679  ////g4reweight
680  fPandoraBeam->Branch("n_g4rw", &n_g4rw, "n_g4rw/I");
681  fPandoraBeam->Branch("g4rw_react_st", &g4rw_react_st, "g4rw_react_st/D");
682  fPandoraBeam->Branch("g4rw_d_react", &g4rw_d_react, "g4rw_d_react/D");
683  fPandoraBeam->Branch("g4rw_elast_st", &g4rw_elast_st, "g4rw_elast_st/D");
684  fPandoraBeam->Branch("g4rw_d_elast", &g4rw_d_elast, "g4rw_d_elast/D");
685  fPandoraBeam->Branch("g4rw_react", &g4rw_react);
686  fPandoraBeam->Branch("g4rw_elast", &g4rw_elast);
687  fPandoraBeam->Branch("g4rw_set_weights", &g4rw_set_weights);
688 
689 
690 
691  ///////////////////////////////////////////////////
692  fPandoraBeam->Branch("primtrk_truth_Z",&primtrk_truth_Z);//primary track true Z positions
693  fPandoraBeam->Branch("primtrk_truth_Z_wire",&primtrk_truth_Z_wire);//primary track true Z positions
694  fPandoraBeam->Branch("primtrk_truth_Z_tpc",&primtrk_truth_Z_tpc);//primary track true Z positions
695 
696 
697  fPandoraBeam->Branch("primtrk_truth_Eng",&primtrk_truth_Eng);//primary track true Energy deposited for each Z position
698  fPandoraBeam->Branch("primtrk_truth_trkide",&primtrk_truth_trkide);//primary track true Energy deposited for each Z position
699  }//begin job
700 
701  void protoana::mcsXsection::analyze(art::Event const & evt){
702 
703  // Initialise tree parameters
704  Initialise();
705 
706  art::ServiceHandle<cheat::ParticleInventoryService> pi_serv;
707  fGeometry = &*(art::ServiceHandle<geo::Geometry>());
708  art::ServiceHandle < geo::Geometry > fGeometryService_rw;
709  // const sim::ParticleList& plist=pi_serv->ParticleList();
710 
711  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService>()->DataFor(evt);
712  auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService>()->DataFor(evt, clockData);
713  // anab::MVAReader<recob::Hit,3> hitResults(evt, fNNetModuleLabel);
714  anab::MVAReader<recob::Hit,4> hitResults(evt, "emtrkmichelid:emtrkmichel");
715  art::ServiceHandle<geo::Geometry> geom;
716 
717 
718  int beamid=-9999;
719  int truthid=-999;
720 
721  fRun = evt.run();
722  fSubRun = evt.subRun();
723  fevent = evt.id().event();
724  art::Timestamp ts = evt.time();
725  if (ts.timeHigh() == 0){
726  TTimeStamp ts2(ts.timeLow());
727  fTimeStamp = ts2.AsDouble();
728  }
729  else{
730  TTimeStamp ts2(ts.timeHigh(), ts.timeLow());
731  fTimeStamp = ts2.AsDouble();
732  }
733  std::cout<<"hhecking code "<<std::endl;
734  bool beamTriggerEvent = false;
735  // If this event is MC then we can check what the true beam particle is
736  auto mcTruths = evt.getValidHandle<std::vector<simb::MCTruth>>(fGeneratorTag);
737  if(!evt.isRealData()){
738  //Spectrometer values
739  //for prod. 3, new implementation to access the beam momentum from spectrometer //////////////////////////////
740  auto beamHandle = evt.getValidHandle<std::vector<beam::ProtoDUNEBeamEvent>>("generator");
741  std::vector<art::Ptr<beam::ProtoDUNEBeamEvent>> beamVec;
742  if( beamHandle.isValid()){
743  art::fill_ptr_vector(beamVec, beamHandle);
744  }
745  const beam::ProtoDUNEBeamEvent & beamEvent = *(beamVec.at(0)); //Should just have one
746 
747  //Access momentum
748  const std::vector< double > & momenta = beamEvent.GetRecoBeamMomenta();
749  std::cout << "Number of reconstructed beam momenta from spec: " << momenta.size() << std::endl;
750 
751  if( momenta.size() > 0 ) std::cout << "Measured Beam Momentum from spec: " << momenta.at(0) << std::endl;
752 
753  //std::cout<<"beam mom size:"<<momenta.size()<<std::endl;
754  for (size_t i = 0; i<momenta.size(); ++i){
755  beamMomentum_spec.push_back(momenta[i]);
756  //std::cout<<"beam mom["<<i<<"]:"<<momenta[i]<<" [GeV]"<<std::endl;
757  }
758 
759  auto & btracks = beamEvent.GetBeamTracks();
760  std::cout<<"beam trk size:"<<btracks.size()<<std::endl;
761  for (size_t i = 0; i<btracks.size(); ++i){
762  std::cout<<"beamPosx/beamPosy/beamPosz:"<<btracks[i].End().X()<<"/"<<btracks[i].End().Y()<<"/"<<btracks[i].End().Z()<<std::endl;
763  std::cout<<"beamDirx/beamDiry/beamDirz:"<<btracks[i].StartDirection().X()<<"/"<<btracks[i].StartDirection().Y()<<"/"<<btracks[i].StartDirection().Z()<<std::endl;
764 
765  beamPosx_spec.push_back(btracks[i].End().X());
766  beamPosy_spec.push_back(btracks[i].End().Y());
767  beamPosz_spec.push_back(btracks[i].End().Z());
768  beamDirx_spec.push_back(btracks[i].StartDirection().X());
769  beamDiry_spec.push_back(btracks[i].StartDirection().Y());
770  beamDirz_spec.push_back(btracks[i].StartDirection().Z());
771 
772  }
773 
774  //////Spectrometer values
775  const simb::MCParticle* geantGoodParticle = truthUtil.GetGeantGoodParticle((*mcTruths)[0],evt);
776  //std::cout<<"geantGoodParticle "<<geantGoodParticle.size()<<std::endl;
777  if(geantGoodParticle != 0x0){
778  std::cout<<"geant good particle loop "<<std::endl;
779  std::cout << "Found GEANT particle corresponding to the good particle with pdg = " << geantGoodParticle->PdgCode()
780  << " , track id = " << geantGoodParticle->TrackId()
781  << " , Vx/Vy/Vz = " << geantGoodParticle->Vx() << "/"<< geantGoodParticle->Vy() << "/" << geantGoodParticle->Vz()
782 
783  << std::endl;
784 
785  std::vector<double> tmp_primtrk_truth_Z;
786  std::vector<double> tmp_primtrk_truth_Eng;
787  std::vector<double> tmp_primtrk_truth_trkide;
788  std::vector<int> tmp_wire;
789  std::vector<int> tmp_tpc;
790 
791  beamTriggerEvent = true;
792  fbeamtrackPos[0] = geantGoodParticle->Vx();
793  fbeamtrackPos[1] = geantGoodParticle->Vy();
794  fbeamtrackPos[2] = geantGoodParticle->Vz();
795  fbeamtrackMomentum = geantGoodParticle->P();
796  fbeamtrackP[0] = geantGoodParticle->Px();
797  fbeamtrackP[1] = geantGoodParticle->Py();
798  fbeamtrackP[2] = geantGoodParticle->Pz();
799  // fbeamtrackEnergy = geantGoodParticle->E();
800  fbeamtrackPdg = geantGoodParticle->PdgCode();
801  fbeamtrackTime = geantGoodParticle->T();
802  fbeamtrackID = geantGoodParticle->TrackId();
803  beamid = geantGoodParticle->TrackId();
804 
805  ///////////////////////////////////////
806  //g4 reweight -------------------------------------------------------------------------------------------------//
807  true_beam_PDG = geantGoodParticle->PdgCode();
808  true_beam_ID = geantGoodParticle->TrackId();
809  //true_beam_len = geantGoodParticle->Trajectory().TotalLength();
810 
811  std::cout << "Doing reweight" << std::endl;
812  if (true_beam_PDG == RW_PDG) { //if PDG=RW_PDG
813  std::vector<G4ReweightTraj *> trajs = CreateNRWTrajs(*geantGoodParticle, pi_serv->ParticleList(),fGeometryService_rw, fevent, "LAr", true);
814  //G4ReweightTraj theTraj(true_beam_ID, true_beam_PDG, 0, fevent, {0,0});
815  int n_rw_array=200;
816  double react_st=0.0;
817  double d_react=.01;
818  double elast_st=.9;
819  double d_elast=.05;
820  n_g4rw=n_rw_array;
821  g4rw_react_st=react_st; g4rw_d_react=d_react;
822  g4rw_elast_st=elast_st; g4rw_d_elast=d_elast;
823  std::cout<<"G4reweighting loop step 1"<<std::endl;
824  if (ParSet.size()==2) { //if parset size=2 from fichl setting, i.e. dim(el+inel)=2
825  for (int i = 0; i < n_rw_array; ++i) {
826  for (int j = 0; j < n_rw_array; ++j) {
827  std::cout<<"loop i, j "<<i<<" "<<j<<std::endl;
828  double tmp_p1=(react_st + i*d_react); //reac
829  double tmp_p2=(elast_st + j*d_elast); //elast
830  std::vector<double> input_values = {tmp_p1,tmp_p2}; //HY:the index here based on the sequence of the fichl setting
831  // std::vector<double> input_values = {ParSet.size(), 1.};
832  std::cout<<"loop temp1, 2 "<<tmp_p1<<" "<<tmp_p2<<std::endl;
833  bool set_values = MultiRW.SetAllParameterValues(input_values);
834  std::cout<<"setvalues "<<set_values<<std::endl;
835  if (!set_values) continue;
836  double temp_w = GetNTrajWeightFromSetPars(trajs, MultiRW);
837  //double temp_w = MultiRW.GetWeightFromSetParameters(theTraj);
838  std::cout<<"weight "<<temp_w<<std::endl;
839  g4rw_react.push_back(tmp_p1);
840  g4rw_elast.push_back(tmp_p2);
841  g4rw_set_weights.push_back(temp_w);
842  //std::cout<<"(inel,el,w):("<<tmp_p1<<","<<tmp_p2<<","<<temp_w<<")"<<std::endl;
843  }
844  }
845  } //if parset size=2
846  } //if PDG=RW_PDG
847  //g4 reweight -------------------------------------------------------------------------------------------------//
848 
849 
850  ////////////////////////////////////////
851 
852  ///
853 
854  truth_last_process=geantGoodParticle->EndProcess();
855  ////here is another new parameter added
856 
857  fprimary_truth_Process = int(geantGoodParticle->Trajectory().ProcessToKey(geantGoodParticle->Process()));
858  prim_energy=0;
859  for(size_t i_s=0; i_s < geantGoodParticle->NumberTrajectoryPoints(); i_s++){ //loop over beam tracks
860  // if(geantGoodParticle->Position(i_s).Z()>0) break;
861  beamtrk_x.push_back(geantGoodParticle->Position(i_s).X());
862  beamtrk_y.push_back(geantGoodParticle->Position(i_s).Y());
863  beamtrk_z.push_back(geantGoodParticle->Position(i_s).Z());
864  double pos_true[3]={geantGoodParticle->Position(i_s).X(),geantGoodParticle->Position(i_s).Y(),geantGoodParticle->Position(i_s).Z()};
865 
866 
867  /* if(geantGoodParticle->Position(i_s).Z()>=0 && geantGoodParticle->Position(i_s).Z()<=230){
868  geo::WireID wireID = geom->NearestWireID(pos_true, 2);
869  if (!wireID) wireID = geom->Plane(2).ClosestWireID(wireID);
870  beamtrk_z_wire.push_back(wireID.Wire);
871  beamtrk_z_tpc.push_back(wireID.TPC);
872  }
873  if(geantGoodParticle->Position(i_s).Z()<0 || geantGoodParticle->Position(i_s).Z()>230){
874  beamtrk_z_wire.push_back(-99999);
875  beamtrk_z_tpc.push_back(-99999);
876 
877  }*/
878  geo::TPCID tpc = geom->FindTPCAtPosition(pos_true);
879  if(tpc.isValid){
880  int tpc_no=tpc.TPC;
881  geo::PlaneID planeID = geo::PlaneID(0, tpc_no, 2);
882  geo::WireID wireID;
883  // geo::WireID wireID = geom->NearestWireID(pos_true, planeID);
884  try{
885  wireID = geom->NearestWireID(pos_true, planeID);
886  }
887  catch(geo::InvalidWireError const& e) {
888  wireID = e.suggestedWireID(); // pick the closest valid wire
889  }
890  beamtrk_z_wire.push_back(wireID.Wire);
891  beamtrk_z_tpc.push_back(wireID.TPC);
892  }
893  if(!tpc.isValid){
894  beamtrk_z_wire.push_back(-9999);
895  beamtrk_z_tpc.push_back(-9999);
896 
897  }
898 
899 
900 
901 
902 
903  /* geo::PlaneGeo const& plane = geom->Plane(2);
904  geo::WireID wireID;
905  try {
906  wireID = plane.NearestWireID(pos_true);
907  }
908  catch (geo::InvalidWireIDError const& e) {
909  if (!e.hasSuggestedWire()) throw;
910  wireID = plane.ClosestWireID(e.suggestedWireID());
911  }
912  beamtrk_z_wire.push_back(wireID.Wire);
913  beamtrk_z_tpc.push_back(wireID.TPC);*/
914  ///////trying new stuff
915 
916 
917 
918  beamtrk_Px.push_back(geantGoodParticle->Momentum(i_s).X());
919  beamtrk_Py.push_back(geantGoodParticle->Momentum(i_s).Y());
920  beamtrk_Pz.push_back(geantGoodParticle->Momentum(i_s).Z());
921 
922  beamtrk_Eng.push_back(geantGoodParticle->Momentum(i_s).E()-geantGoodParticle->Mass());
923  if(geantGoodParticle->Position(i_s).Z()<0) prim_energy=1000*(geantGoodParticle->Momentum(i_s).E()-geantGoodParticle->Mass());
924  if(geantGoodParticle->Position(i_s).Z()<0) fbeamtrackEnergy=prim_energy; //correct energy at the beginning of track
925  } //loop over beam trks
926 
927  //new section
928  art::ServiceHandle<cheat::BackTrackerService> bt_serv;
929 
930 
931  // art::ServiceHandle<geo::Geometry> geom;
932  simb::MCTrajectory truetraj=geantGoodParticle->Trajectory();
933  auto thisTrajectoryProcessMap1 = truetraj.TrajectoryProcesses();
934  if (thisTrajectoryProcessMap1.size()){
935  for(auto const& couple: thisTrajectoryProcessMap1){
936  // int_label=truetraj.KeyToProcess(couple.second);
937  fprimary_truth_EndPosition[0]=((truetraj.at(couple.first)).first).X();
938  fprimary_truth_EndPosition[1]=((truetraj.at(couple.first)).first).Y();
939  fprimary_truth_EndPosition[2]=((truetraj.at(couple.first)).first).Z();
940  fprimary_truth_EndProcess=truetraj.KeyToProcess(couple.second);
941  fprimary_truth_Momentum[0]=((truetraj.at(couple.first)).second).X();
942  fprimary_truth_Momentum[1]= ((truetraj.at(couple.first)).second).Y();
943  fprimary_truth_Momentum[2]=((truetraj.at(couple.first)).second).Z();
944  break;
945  }
946  }
947 
948  ////saving the complete information of all the interactions
949  std::cout<<"interaction map size "<<thisTrajectoryProcessMap1.size()<<std::endl;
950  if (thisTrajectoryProcessMap1.size()){
951  for(auto const& couple1: thisTrajectoryProcessMap1){
952 
953  if ((truetraj.KeyToProcess(couple1.second)).find("CoulombScat")!= std::string::npos) continue;
954  // Let's check if the interaction is in the the TPC
955  auto interactionPos4D = (truetraj.at(couple1.first)).first ;
956  if (interactionPos4D.Z() < minZ || interactionPos4D.Z() > maxZ ) continue;
957  else if (interactionPos4D.X() < minX || interactionPos4D.X() > maxX ) continue;
958  else if (interactionPos4D.Y() < minY || interactionPos4D.Y() > maxY ) continue;
959  interactionX.push_back(((truetraj.at(couple1.first)).first).X());
960  interactionY.push_back(((truetraj.at(couple1.first)).first).Y());
961  interactionZ.push_back(((truetraj.at(couple1.first)).first).Z());
962  double xval=((truetraj.at(couple1.first)).first).X();
963  double zval=((truetraj.at(couple1.first)).first).Z();
964  unsigned int tpcno=1;
965  if(xval<=0 && zval<232) tpcno=1;
966  if(xval<=0 && zval>232 && zval<464) tpcno=5;
967  if(xval<=0 && zval>=464) tpcno=9;
968  if(xval>0 && zval<232) tpcno=2;
969  if(xval>0 && zval>232 && zval<464) tpcno=6;
970  if(xval>0 && zval>=464) tpcno=10;
971 
972  interactionT.push_back(detProp.ConvertXToTicks(((truetraj.at(couple1.first)).first).X(), 2, tpcno, 0));
973  interactionU.push_back(fGeometry->WireCoordinate(((truetraj.at(couple1.first)).first).Y(), ((truetraj.at(couple1.first)).first).Z(),0, tpcno, 0));
974  interactionV.push_back(fGeometry->WireCoordinate(((truetraj.at(couple1.first)).first).Y(), ((truetraj.at(couple1.first)).first).Z(),1, tpcno, 0));
975  interactionW.push_back(fGeometry->WireCoordinate(((truetraj.at(couple1.first)).first).Y(), ((truetraj.at(couple1.first)).first).Z(),2, tpcno, 0));
976  interactionProcesslist.push_back(truetraj.KeyToProcess(couple1.second));
977  std::cout<<"number of interactions "<<thisTrajectoryProcessMap1.size()<<std::endl;
978  std::cout<<"int X, Y, Z and process "<<((truetraj.at(couple1.first)).first).X()<<" "<<((truetraj.at(couple1.first)).first).Y()<<" "<<((truetraj.at(couple1.first)).first).Z()<<" "<<truetraj.KeyToProcess(couple1.second)<<std::endl;
979  ///get the interaction angle here
980  double interactionAngle = 999999.; // This needs to be changed
981  double ut=999999;
982  double vt=999999;
983  double zt=999999;
984  //--------------------- Int Angle ---------------------------
985  // Try to retreive the interaction angle
986  auto prevInteractionPos4D = (truetraj.at(couple1.first-1)).first ;
987  auto prevInteractionPos3D = prevInteractionPos4D.Vect() ;
988  auto interactionPos3D = interactionPos4D.Vect() ;
989  auto distanceBtwPoint = interactionPos3D - prevInteractionPos3D;
990  //Let's try to see if the next point exists
991  if (truetraj.size() > couple1.first + 1) {
992  // The particle doesn't die. No need to check for anything else.
993  auto nextInteractionPos4D = (truetraj.at(couple1.first+1)).first ;
994  auto nextInteractionPos3D = nextInteractionPos4D.Vect() ;
995  auto distanceBtwPointNext = nextInteractionPos3D - interactionPos3D;
996  std::cout<<"distance between points values "<<distanceBtwPoint.X()<<" "<<distanceBtwPoint.Z()<<" "<<distanceBtwPointNext.X()<<" "<<distanceBtwPointNext.Z()<<std::endl;
997 
998  interactionAngles.push_back(TMath::ACos(distanceBtwPointNext.Dot(distanceBtwPoint)/(distanceBtwPointNext.Mag()*distanceBtwPoint.Mag() ) ));
999 
1000 
1001  double u0=0.4669*fGeometry->WireCoordinate(prevInteractionPos3D.Y(),prevInteractionPos3D.Z(),0,1,0);
1002  double v0=0.4669*fGeometry->WireCoordinate(prevInteractionPos3D.Y(),prevInteractionPos3D.Z(),1,1,0);
1003  double z0=0.4792*prevInteractionPos3D.Z();
1004  double x0=prevInteractionPos3D.X();
1005 
1006  double u1=0.4669*fGeometry->WireCoordinate(interactionPos3D.Y(),interactionPos3D.Z(),0,1,0);
1007  double v1=0.4669*fGeometry->WireCoordinate(interactionPos3D.Y(),interactionPos3D.Z(),1,1,0);
1008  double z1=0.4792*interactionPos3D.Z();
1009  double x1=interactionPos3D.Z();
1010 
1011  double u2=0.4669*fGeometry->WireCoordinate(nextInteractionPos3D.Y(),nextInteractionPos3D.Z(),0,1,0);
1012  double v2=0.4669*fGeometry->WireCoordinate(nextInteractionPos3D.Y(),nextInteractionPos3D.Z(),1,1,0);
1013  double z2=0.4792*nextInteractionPos3D.Z();
1014  double x2=nextInteractionPos3D.Z();
1015 
1016  interactionAnglesUT.push_back(angle2d(u0,x0,u1,x1,u2,x2));
1017  interactionAnglesVT.push_back(angle2d(v0,x0,v1,x1,v2,x2));
1018  interactionAnglesZT.push_back(angle2d(z0,x0,z1,x1,z2,x2));
1019 
1020 
1021  continue;
1022  }
1023  interactionAnglesUT.push_back(ut);
1024  interactionAnglesVT.push_back(vt);
1025  interactionAnglesZT.push_back(zt);
1026  interactionAngles.push_back(interactionAngle);
1027  }
1028  }
1029 
1030  geo::View_t view = geom->View(2);
1031  auto simIDE_prim=bt_serv->TrackIdToSimIDEs_Ps(geantGoodParticle->TrackId(),view);
1032  std::map<double, sim::IDE> orderedSimIDE;
1033  for (auto& ide : simIDE_prim) orderedSimIDE[ide->z]= *ide;
1034  auto inTPCPoint = truetraj.begin();
1035  auto Momentum0 = inTPCPoint->second;
1036  auto old_iter = orderedSimIDE.begin();
1037  double tlen=0.0;
1038  double xi=0.0;double yi=0.0;double zi=0.0;
1039  int count=0;
1040 
1041 
1042  for ( auto iter= orderedSimIDE.begin(); iter!= orderedSimIDE.end(); iter++,old_iter++){
1043  auto currentIde = iter->second;
1044  if(currentIde.z<minZ) continue;
1045  else if (currentIde.x < minX || currentIde.x > maxX ) continue;
1046  else if (currentIde.y < minY || currentIde.y > maxY ) continue;
1047  tmp_primtrk_truth_Z.push_back(currentIde.z);
1048  tmp_primtrk_truth_Eng.push_back(currentIde.energy);
1049  tmp_primtrk_truth_trkide.push_back(currentIde.trackID);
1050  double pos_true[3]={currentIde.x,currentIde.y,currentIde.z};
1051  /*if(currentIde.z>=0 && currentIde.z<=230){
1052  geo::WireID wireID = geom->NearestWireID(pos_true, 2);
1053  if (!wireID) wireID = geom->Plane(2).ClosestWireID(wireID);
1054  tmp_wire.push_back(wireID.Wire);
1055  tmp_tpc.push_back(wireID.TPC);
1056  }
1057  if(currentIde.z<0 && currentIde.z>230){
1058  tmp_wire.push_back(-99999);
1059  tmp_tpc.push_back(-99999);
1060  }*/
1061  geo::TPCID tpc1 = geom->FindTPCAtPosition(pos_true);
1062  if(tpc1.isValid){
1063  int tpc_no=tpc1.TPC;
1064  geo::PlaneID planeID = geo::PlaneID(0, tpc_no, 2);
1065  // geo::WireID wireID = geom->NearestWireID(pos_true, planeID);
1066  geo::WireID wireID;
1067  try{
1068  wireID = geom->NearestWireID(pos_true, planeID);
1069  }
1070  catch(geo::InvalidWireError const& e) {
1071  wireID = e.suggestedWireID(); // pick the closest valid wire
1072  }
1073  tmp_wire.push_back(wireID.Wire);
1074  tmp_tpc.push_back(wireID.TPC);
1075  }
1076  if(!tpc1.isValid){
1077  tmp_wire.push_back(-9999);
1078  tmp_tpc.push_back(-9999);
1079 
1080  }
1081 
1082 
1083  if(count==0){
1084  fprimary_truth_StartPosition[0] = currentIde.x;
1085  fprimary_truth_StartPosition[1] = currentIde.y;
1086  fprimary_truth_StartPosition[2] = currentIde.z;
1087  }
1088  if(currentIde.trackID>=0){
1089  if(count>0){
1090  tlen=tlen+TMath::Sqrt(std::pow(currentIde.x-xi,2)+std::pow(currentIde.y-yi,2)+std::pow(currentIde.z-zi,2));
1091  }//if count
1092  xi=currentIde.x;yi=currentIde.y;zi=currentIde.z;
1093  count++;
1094  }//trackid>0 loop
1095  }// iter loop
1096  fprimary_truth_tracklength=tlen;
1097  primtrk_truth_Z.push_back(tmp_primtrk_truth_Z);
1098  primtrk_truth_Eng.push_back(tmp_primtrk_truth_Eng);
1099  primtrk_truth_trkide.push_back(tmp_primtrk_truth_trkide);
1100  primtrk_truth_Z_wire.push_back(tmp_wire);
1101  primtrk_truth_Z_tpc.push_back(tmp_tpc);
1102 
1103 
1104  tmp_primtrk_truth_Z.clear();
1105  tmp_primtrk_truth_Eng.clear();
1106  tmp_primtrk_truth_trkide.clear();
1107  tmp_wire.clear();
1108  tmp_tpc.clear();
1109  //new section
1110 
1111  }// geantGoodParticle
1112  }//is not real data loop
1113 
1114  /*
1115  // Now we want to access the output from Pandora. This comes in the form of particle flow objects (recob::PFParticle).
1116  // The primary PFParticles are those we want to consider and these PFParticles then have a hierarchy of daughters that
1117  // describe the whole interaction of a given primary particle
1118  //
1119  // / daughter track
1120  // /
1121  // primary track /
1122  // ---------------- ---- daughter track
1123  // \
1124  // /\-
1125  // /\\-- daughter shower
1126  //
1127  // The above primary PFParticle will have links to three daughter particles, two track-like and one shower-like
1128  */
1129 
1130  // Track momentum algorithm calculates momentum based on track range
1131  trkf::TrackMomentumCalculator trmom;
1132  //trmom.SetMinLength(100);
1133 
1134  // Get all of the PFParticles, by default from the "pandora" product
1135  auto recoParticles = evt.getValidHandle<std::vector<recob::PFParticle>>(fPFParticleTag);
1136  std::cout << "All primary pfParticles = " << pfpUtil.GetNumberPrimaryPFParticle(evt,fPFParticleTag) << std::endl;
1137 
1138  // We'd like to find the beam particle. Pandora tries to do this for us, so let's use the PFParticle utility
1139  // to look for it. Pandora reconstructs slices containing one (or sometimes more) primary PFParticles. These
1140  // are tagged as either beam or cosmic for ProtoDUNE. This function automatically considers only those
1141  // PFParticles considered as primary
1142  // std::vector<recob::PFParticle*> pfParticles = pfpUtil.GetPFParticlesFromBeamSlice(evt,fPFParticleTag);
1143 
1144 
1145 
1146  auto pfParticles = pfpUtil.GetPFParticlesFromBeamSlice(evt,fPFParticleTag);
1147  //cluster information
1148 
1149  std::vector<art::Ptr<recob::Hit> > hitlist;
1150  auto hitListHandle = evt.getHandle< std::vector<recob::Hit> >("hitpdune");
1151  if (hitListHandle) art::fill_ptr_vector(hitlist, hitListHandle);
1152 
1153  // Implementation of required member function here.
1154  std::vector<art::Ptr<recob::Track> > tracklist;
1155  auto trackListHandle = evt.getHandle< std::vector<recob::Track> >("pandoraTrack");
1156  if (trackListHandle) art::fill_ptr_vector(tracklist, trackListHandle);
1157  else return;
1158  art::FindManyP<recob::Track> thass(hitListHandle, evt, "pandoraTrack"); //to associate hit just trying
1159 
1160  std::vector<art::Ptr<recob::PFParticle> > pfplist;
1161  auto PFPListHandle = evt.getHandle< std::vector<recob::PFParticle> >("pandora");
1162  if (PFPListHandle) art::fill_ptr_vector(pfplist, PFPListHandle);
1163 
1164  std::vector<art::Ptr<recob::Cluster>> clusterlist;
1165  auto clusterListHandle = evt.getHandle< std::vector<recob::Cluster> >("pandora"); // to get information about the hits
1166  if (clusterListHandle) art::fill_ptr_vector(clusterlist, clusterListHandle);
1167 
1168  art::FindManyP<recob::Cluster> fmcp(PFPListHandle,evt,"pandora");
1169  art::FindManyP<recob::Track> pftrack(PFPListHandle,evt,"pandoraTrack");
1170  art::FindManyP<recob::Hit> clhit(clusterListHandle,evt,"pandora");
1171 
1172  std::cout<<"number of pfp_particles "<<pfplist.size()<<std::endl;
1173  std::cout<<" size of pfParticles testing size "<<pfParticles.size()<<std::endl;
1174  art::FindManyP<recob::Hit, recob::TrackHitMeta> fmthm(trackListHandle, evt,"pandoraTrack"); // to associate tracks and hits
1175 
1176 
1177  std::cout<<"outside pfparticle loop"<<std::endl;
1178  // We can now look at these particles
1179  for(const recob::PFParticle* particle : pfParticles){
1180  // Pandora's BDT beam-cosmic score
1181  fprimaryBDTScore = (double)pfpUtil.GetBeamCosmicScore(*particle,evt,fPFParticleTag);
1182  // NHits associated with this pfParticle
1183  fprimaryNHits = (pfpUtil.GetPFParticleHits(*particle,evt,fPFParticleTag)).size();
1184  // of this particle might be more helpful. These return null pointers if not track-like / shower-like
1185  const recob::Track* thisTrack = pfpUtil.GetPFParticleTrack(*particle,evt,fPFParticleTag,fTrackerTag);
1186  const recob::Shower* thisShower = pfpUtil.GetPFParticleShower(*particle,evt,fPFParticleTag,fShowerTag);
1187  /////new line added here
1188  // std::vector<art::Ptr<recob::Cluster>> allClusters=fmcp.at(particle);
1189  // std::cout<<allClusters.size();
1190 
1191  /////////////////////Michel tagging stage//////////
1192 
1193  double beamstx=-30;
1194  double beamendx=-30;
1195  double beamsty=420;
1196  double beamendy=420;
1197  double beamstz=30;
1198  double beamendz=100;
1199 
1200  if(thisTrack != 0x0){
1201  if(!beam_cuts.IsBeamlike(*thisTrack, evt, "1")) return;
1202  beamstx=thisTrack->Start().X();
1203  beamsty=thisTrack->Start().Y();
1204  beamstz=thisTrack->Start().Z();
1205  beamendx=thisTrack->End().X();
1206  beamendy=thisTrack->End().Y();
1207  beamendz=thisTrack->End().Z();
1208  std::cout<<"beamstx "<<beamstx<<std::endl;
1209 
1210  //////Michel tagging here
1211  std::vector<double> secondarystartx1;
1212  std::vector<double> secondarystarty1;
1213  std::vector<double> secondarystartz1;
1214  std::vector<double> secondaryendx1;
1215  std::vector<double> secondaryendy1;
1216  std::vector<double> secondaryendz1;
1217  std::vector<double> dQmichel1;
1218  std::vector<double> dQtrackbegin1;
1219  std::vector<double> dQtrackend1;
1220  std::vector<double> tracklengthsecondary1;
1221  std::vector<double> primsectheta1;
1222  std::vector<int> endhitssecondary1;
1223  std::vector<int> MtrackID1;
1224  std::vector<double> trks1;
1225  std::vector<double> ems1;
1226  std::vector<double> michels1;
1227  std::vector<double> nones1;
1228  size_t NTracks = tracklist.size();
1229  std::cout<<"number of tracks "<<NTracks<<std::endl;
1230  for(size_t i=0;i<NTracks;i++){
1231  art::Ptr<recob::Track> ptrack(trackListHandle, i);
1232  const recob::Track& track = *ptrack;
1233  auto pos = track.Vertex();
1234  auto end = track.End();
1235  int counter1=0;
1236  double startx=pos.X();
1237  double starty=pos.Y();
1238  double startz=pos.Z();
1239  double endx=end.X();
1240  double endy=end.Y();
1241  double endz=end.Z();
1242  if(track.Length()<5) continue;
1243  if(TMath::Max(endy,starty)>500 || TMath::Min(endy, starty)<200 || TMath::Max(startx, endx)>0||TMath::Min(startx,endx)<-200||TMath::Max(startz,endz)<30) continue;
1244  //if(TMath::Max(endy,starty)>520 || TMath::Min(endy, starty)<150 || TMath::Max(startx, endx)>20||TMath::Min(startx,endx)<-300||TMath::Max(startz,endz)<230) continue;
1245 
1246  std::vector<int> wirenos;
1247  std::vector<float> peakts,dqbuff1;
1248  std::vector<float> dQstart,dQend;
1249  std::vector<double> micheldq;
1250  wirenos.clear();peakts.clear();dqbuff1.clear();
1251  float peaktime=-1;
1252  int wireno=-99999;
1253  int tpcno=-1;
1254  float zlast0=-99999;
1255  float zlast=-99999;
1256  std::vector<std::tuple<double,double,double,double,int,double>> buff_ZYXTWQ;
1257  buff_ZYXTWQ.clear();
1258  double thetavalue=theta12(beamstx,beamendx,beamsty,beamendy,beamstz,beamendz,startx,endx,starty,endy,startz,endz);
1259  if(fmthm.isValid()){
1260  auto vhit=fmthm.at(i);
1261  auto vmeta=fmthm.data(i);
1262  for (size_t ii = 0; ii<vhit.size(); ++ii){ //loop over all meta data hit
1263  bool fBadhit = false;
1264  if (vmeta[ii]->Index() == static_cast<unsigned int>(std::numeric_limits<int>::max())){
1265  fBadhit = true;
1266  //cout<<"fBadHit"<<fBadhit<<endl;
1267  continue;
1268  }
1269  if (vmeta[ii]->Index()>=tracklist[i]->NumberTrajectoryPoints()){
1270  throw cet::exception("Calorimetry_module.cc") << "Requested track trajectory index "<<vmeta[ii]->Index()<<" exceeds the total number of trajectory points "<<tracklist[i]->NumberTrajectoryPoints()<<" for track index "<<i<<". Something is wrong with the track reconstruction. Please contact tjyang@fnal.gov!!";
1271  }
1272  if (!tracklist[i]->HasValidPoint(vmeta[ii]->Index())){
1273  fBadhit = true;
1274  // cout<<"had valid point "<<fBadhit<<endl;
1275  continue;
1276  }
1277 
1278  auto loc = tracklist[i]->LocationAtPoint(vmeta[ii]->Index());
1279  if (fBadhit) continue; //HY::If BAD hit, skip this hit and go next
1280  if (loc.Z()<-100) continue; //hit not on track
1281  if(vhit[ii]->WireID().Plane==2){
1282  buff_ZYXTWQ.push_back(std::make_tuple(loc.Z(),loc.Y(),loc.X(),vhit[ii]->PeakTime(),vhit[ii]->WireID().Wire,vhit[ii]->Integral()));
1283  wirenos.push_back(vhit[ii]->WireID().Wire);
1284  peakts.push_back(vhit[ii]->PeakTime());
1285  zlast=loc.Z();
1286  if(zlast>zlast0){
1287  zlast0=zlast;
1288  wireno=vhit[ii]->WireID().Wire;
1289  peaktime=vhit[ii]->PeakTime();
1290  tpcno=vhit[ii]->WireID().TPC;
1291  }
1292  }//planenum 2
1293  }//loop over vhit
1294  }//fmthm valid
1295  //save start and end point of each track
1296  //taking care of flipped start and end point
1297  if(endz<startz){
1298  startx=end.X();
1299  starty=end.Y();
1300  startz=end.Z();
1301  endx=pos.X();
1302  endy=pos.Y();
1303  endz=pos.Z();
1304  }
1305  double trk_score=0.0;
1306  double em_score=0;
1307  double michel_score=0;
1308  double none_score=0;
1309  for(size_t hitl=0;hitl<hitlist.size();hitl++){
1310  std::array<float,4> cnn_out=hitResults.getOutput(hitlist[hitl]);
1311  auto & tracks = thass.at(hitlist[hitl].key());
1312  if (!tracks.empty() && tracks[0].key()!=ptrack.key() && tracklist[tracks[0].key()]->Length()>25) continue;
1313  // if (!tracks.empty() && tracks[0].key()!=ptrack.key() && tracklist[tracks[0].key()]->Length()>25) continue;
1314  bool test=true;
1315  float peakth1=hitlist[hitl]->PeakTime();
1316  int wireh1=hitlist[hitl]->WireID().Wire;
1317  for(size_t m=0;m<wirenos.size();m++){
1318  if(wireh1==wirenos[m] && peakth1==peakts[m]){
1319  test=false;
1320  break;
1321  }
1322  }
1323  if(!test) continue;
1324  int planeid=hitlist[hitl]->WireID().Plane;
1325  int tpcid=hitlist[hitl]->WireID().TPC;
1326  if(abs(wireh1-wireno)<15 && abs(peakth1-peaktime)<100 && planeid==2 && tpcid==tpcno){
1327  // if(abs(wireh1-wireno)<20 && abs(peakth1-peaktime)<150 && planeid==2 && tpcid==tpcno){
1328  counter1++;
1329  // std::cout<<"wireno, counter "<<wireno<<" "<<counter1<<std::endl;
1330  micheldq.push_back(hitlist[hitl]->Integral());
1331  trk_score+=cnn_out[hitResults.getIndex("track")];
1332  em_score+=cnn_out[hitResults.getIndex("em")];
1333  michel_score+=cnn_out[hitResults.getIndex("michel")];
1334  none_score+=cnn_out[hitResults.getIndex("none")];
1335  // std::cout<<"track, em, michel none"<<cnn_out[hitResults.getIndex("track")]<<" "<<cnn_out[hitResults.getIndex("em")]<<" "<<cnn_out[hitResults.getIndex("michel")]<<" "<<cnn_out[hitResults.getIndex("none")]<<std::endl;
1336  //std::cout<<"wire, peaktime selected hit, wires "<<wireno<<" "<<peaktime<<" "<<hitlist[hitl]->PeakTime()<<" "<<hitlist[hitl]->WireID().Wire<<std::endl;
1337  }
1338  }//hitlist loop
1339  if(buff_ZYXTWQ.size()<10) continue;
1340  sort(buff_ZYXTWQ.begin(),buff_ZYXTWQ.end());
1341  dQstart.clear(); dQend.clear();
1342  int qi11=buff_ZYXTWQ.size();
1343  for(int qi=5;qi<TMath::Min(15,qi11);qi++){
1344  dQstart.push_back(std::get<5>(buff_ZYXTWQ[qi]));
1345  }
1346 
1347  for(int qi=qi11-5;qi<qi11;qi++){
1348  dQend.push_back(std::get<5>(buff_ZYXTWQ[qi]));
1349  }
1350  secondarystartx1.push_back(startx);
1351  secondarystarty1.push_back(starty);
1352  secondarystartz1.push_back(startz);
1353  secondaryendx1.push_back(endx);
1354  secondaryendy1.push_back(endy);
1355  secondaryendz1.push_back(endz);
1356  endhitssecondary1.push_back(counter1);
1357  tracklengthsecondary1.push_back(track.Length());
1358  dQmichel1.push_back(TMath::Median(micheldq.size(),&micheldq[0]));
1359  dQtrackbegin1.push_back(TMath::Median(dQstart.size(),&dQstart[0]));
1360  dQtrackend1.push_back(TMath::Median(dQend.size(),&dQend[0]));
1361  primsectheta1.push_back(thetavalue);
1362  MtrackID1.push_back(track.ID());
1363  trks1.push_back(trk_score);
1364  ems1.push_back(em_score);
1365  michels1.push_back(michel_score);
1366  nones1.push_back(none_score);
1367  std::cout<<"avg, trackscore, emscore, michelscore, nonescore "<<trk_score<<" "<<em_score<<" "<<michel_score<<" "<<none_score<<std::endl;
1368  }//Ntracks
1369  secondarystartx.push_back(secondarystartx1);
1370  secondarystarty.push_back(secondarystarty1);
1371  secondarystartz.push_back(secondarystartz1);
1372  secondaryendx.push_back(secondaryendx1);
1373  secondaryendy.push_back(secondaryendy1);
1374  secondaryendz.push_back(secondaryendz1);
1375  endhitssecondary.push_back(endhitssecondary1);
1376  tracklengthsecondary.push_back(tracklengthsecondary1);
1377  dQmichel.push_back(dQmichel1);
1378  dQtrackbegin.push_back(dQtrackbegin1);
1379  dQtrackend.push_back(dQtrackend1);
1380  primsectheta.push_back(primsectheta1);
1381  MtrackID.push_back(MtrackID1);
1382 
1383  trackscore.push_back(trks1);
1384  emscore.push_back(ems1);
1385  michelscore.push_back(michels1);
1386  nonescore.push_back(nones1);
1387 
1388  endhitssecondary1.clear();
1389  secondarystartx1.clear();
1390  secondaryendx1.clear();
1391  secondarystarty1.clear();
1392  secondaryendy1.clear();
1393  secondarystartz1.clear();
1394  secondaryendz1.clear();
1395  dQmichel1.clear();
1396  primsectheta1.clear();
1397  dQtrackbegin1.clear();
1398  dQtrackend1.clear();
1399  tracklengthsecondary1.clear();
1400  MtrackID1.clear();
1401  trks1.clear();
1402  ems1.clear();
1403  michels1.clear();
1404  nones1.clear();
1405  }
1406  ///////////////End of Michel checking
1407 
1408  /////Michel tagging ends here
1409  if(thisTrack != 0x0){
1410  if(!beam_cuts.IsBeamlike(*thisTrack, evt, "1")) return;
1411  // Get the true mc particle
1412  const simb::MCParticle* mcparticle = truthUtil.GetMCParticleFromRecoTrack(clockData, *thisTrack, evt, fTrackerTag);
1413  if(mcparticle!=0x0){
1414  std::cout<<"ftruth pdg "<<mcparticle->PdgCode()<<std::endl;
1415  ftruthpdg=mcparticle->PdgCode();
1416  truthid=mcparticle->TrackId();
1417  fprimary_truth_Pdg= mcparticle->PdgCode();
1418  std::cout<<"primary_truth_Pdg "<<fprimary_truth_Pdg<<std::endl;
1419  fprimary_truth_Isbeammatched=0;
1420 
1421  fprimary_truth_TrackId = mcparticle->TrackId();
1422  fprimary_truth_StartPosition[3] = mcparticle->T();
1423  fprimary_truth_EndPosition[3] = mcparticle->EndT();
1424  fprimary_truth_P = mcparticle->P();
1425  fprimary_truth_Momentum[3] = mcparticle->E();
1426  fprimary_truth_Pt = mcparticle->Pt();
1427  fprimary_truth_Mass = mcparticle->Mass();
1428  fprimary_truth_EndMomentum[0] = mcparticle->EndPx();
1429  fprimary_truth_EndMomentum[1] = mcparticle->EndPy();
1430  fprimary_truth_EndMomentum[2] = mcparticle->EndPz();
1431  fprimary_truth_EndMomentum[3] = mcparticle->EndE();
1432  fprimary_truth_Theta = mcparticle->Momentum().Theta();
1433  fprimary_truth_Phi = mcparticle->Momentum().Phi();
1434  fprimary_truth_NDaughters = mcparticle->NumberDaughters();
1435  if(beamid==truthid) fprimary_truth_Isbeammatched=1;
1436  }//mcparticle loop
1437 
1438 
1439  fisprimarytrack = 1;
1440  fisprimaryshower = 0;
1441  fprimaryID = thisTrack->ID();
1442  fprimaryTheta = thisTrack->Theta();
1443  fprimaryPhi = thisTrack->Phi();
1444  fprimaryLength = thisTrack->Length();
1445  fprimaryMomentum = thisTrack->StartMomentum();
1446  fprimaryEndMomentum = thisTrack->EndMomentum();
1447  fprimaryEndPosition[0] = thisTrack->End().X();
1448  fprimaryEndPosition[1] = thisTrack->End().Y();
1449  fprimaryEndPosition[2] = thisTrack->End().Z();
1450  fprimaryStartPosition[0] = thisTrack->Start().X();
1451  fprimaryStartPosition[1] = thisTrack->Start().Y();
1452  fprimaryStartPosition[2] = thisTrack->Start().Z();
1453  fprimaryEndDirection[0] = thisTrack->Trajectory().EndDirection().X();
1454  fprimaryEndDirection[1] = thisTrack->Trajectory().EndDirection().Y();
1455  fprimaryEndDirection[2] = thisTrack->Trajectory().EndDirection().Z();
1456  fprimaryStartDirection[0] = thisTrack->Trajectory().StartDirection().X();
1457  fprimaryStartDirection[1] = thisTrack->Trajectory().StartDirection().Y();
1458  fprimaryStartDirection[2] = thisTrack->Trajectory().StartDirection().Z();
1459  fprimaryMomentumByRangeMuon = trmom.GetTrackMomentum(thisTrack->Length(),13);
1460  fprimaryMomentumByRangeProton = trmom.GetTrackMomentum(thisTrack->Length(),2212);
1461  std::vector<anab::Calorimetry> calovector = trackUtil.GetRecoTrackCalorimetry(*thisTrack, evt, fTrackerTag, fCalorimetryTag);
1462  if(calovector.size() != 3)
1463  std::cerr << "WARNING::Calorimetry vector size for primary is = " << calovector.size() << std::endl;
1464  std::vector<double> tmp_primtrk_dqdx;
1465  std::vector<double> tmp_primtrk_resrange;
1466  std::vector<double> tmp_primtrk_dedx;
1467  std::vector<double> tmp_primtrk_hitx;
1468  std::vector<double> tmp_primtrk_hity;
1469  std::vector<double> tmp_primtrk_hitz;
1470  std::vector<double> tmp_primtrk_pitch;
1471  std::vector<int> tmp_zwire;
1472  std::vector<int> tmp_ztpc;
1473 
1474 
1475  for (auto & calo : calovector) {
1476  if (calo.PlaneID().Plane == 2){ //only collection plane
1477  primtrk_range.push_back(calo.Range());
1478  for (size_t ihit = 0; ihit < calo.dQdx().size(); ++ihit){ //loop over hits
1479  tmp_primtrk_dqdx.push_back(calo.dQdx()[ihit]);
1480  tmp_primtrk_resrange.push_back(calo.ResidualRange()[ihit]);
1481  tmp_primtrk_dedx.push_back(calo.dEdx()[ihit]);
1482  tmp_primtrk_pitch.push_back(calo.TrkPitchVec()[ihit]);
1483 
1484  const auto &primtrk_pos=(calo.XYZ())[ihit];
1485  tmp_primtrk_hitx.push_back(primtrk_pos.X());
1486  tmp_primtrk_hity.push_back(primtrk_pos.Y());
1487  tmp_primtrk_hitz.push_back(primtrk_pos.Z());
1488 
1489  double pos_true[3]={primtrk_pos.X(),primtrk_pos.Y(),primtrk_pos.Z()};
1490  /* if(primtrk_pos.Z()>=0 && primtrk_pos.Z()<=230){
1491  geo::WireID wireID = geom->NearestWireID(pos_true, 2);
1492  if (!wireID) wireID = geom->Plane(2).ClosestWireID(wireID);
1493  tmp_zwire.push_back(wireID.Wire);
1494  tmp_ztpc.push_back(wireID.TPC);
1495  }
1496  if(primtrk_pos.Z()<0 || primtrk_pos.Z()>230){
1497  tmp_zwire.push_back(-99999);
1498  tmp_ztpc.push_back(-99999);
1499  }
1500  */
1501  geo::TPCID tpc2 = geom->FindTPCAtPosition(pos_true);
1502  if(tpc2.isValid){
1503  int tpc_no=tpc2.TPC;
1504  geo::PlaneID planeID = geo::PlaneID(0, tpc_no, 2);
1505  geo::WireID wireID;
1506  //geo::WireID wireID = geom->NearestWireID(pos_true, planeID);
1507  try{
1508  wireID = geom->NearestWireID(pos_true, planeID);
1509  }
1510  catch(geo::InvalidWireError const& e) {
1511  wireID = e.suggestedWireID(); // pick the closest valid wire
1512  }
1513  tmp_zwire.push_back(wireID.Wire);
1514  tmp_ztpc.push_back(wireID.TPC);
1515  }
1516  if(!tpc2.isValid){
1517  tmp_zwire.push_back(-9999);
1518  tmp_ztpc.push_back(-9999);
1519  }
1520  } //loop over hits
1521  } //only collection plane
1522  }//calovector
1523 
1524  if(tmp_primtrk_dqdx.size()!=0){
1525  if(tmp_primtrk_hitz[0]>tmp_primtrk_hitz[tmp_primtrk_hitz.size()-1]){
1526  std::reverse(tmp_primtrk_hitz.begin(),tmp_primtrk_hitz.end());
1527  std::reverse(tmp_primtrk_hity.begin(),tmp_primtrk_hity.end());
1528  std::reverse(tmp_primtrk_hitx.begin(),tmp_primtrk_hitx.end());
1529  std::reverse(tmp_primtrk_pitch.begin(),tmp_primtrk_pitch.end());
1530  std::reverse(tmp_primtrk_dedx.begin(),tmp_primtrk_dedx.end());
1531  std::reverse(tmp_primtrk_dqdx.begin(),tmp_primtrk_dqdx.end());
1532  std::reverse(tmp_primtrk_resrange.begin(),tmp_primtrk_resrange.end());
1533 
1534  std::reverse(tmp_ztpc.begin(),tmp_ztpc.end());
1535  std::reverse(tmp_zwire.begin(),tmp_zwire.end());
1536  }
1537  primtrk_dqdx.push_back(tmp_primtrk_dqdx);
1538  primtrk_resrange.push_back(tmp_primtrk_resrange);
1539  primtrk_dedx.push_back(tmp_primtrk_dedx);
1540  primtrk_hitx.push_back(tmp_primtrk_hitx);
1541  primtrk_hity.push_back(tmp_primtrk_hity);
1542  primtrk_hitz.push_back(tmp_primtrk_hitz);
1543  primtrk_pitch.push_back(tmp_primtrk_pitch);
1544 
1545  primtrk_hitz_wire.push_back(tmp_zwire);
1546  primtrk_hitz_tpc.push_back(tmp_ztpc);
1547 
1548  }
1549  tmp_primtrk_dqdx.clear();
1550  tmp_primtrk_resrange.clear();
1551  tmp_primtrk_dedx.clear();
1552  tmp_primtrk_hitx.clear();
1553  tmp_primtrk_hity.clear();
1554  tmp_primtrk_hitz.clear();
1555  tmp_primtrk_pitch.clear();
1556  tmp_ztpc.clear();
1557  tmp_zwire.clear();
1558 
1559  for(size_t k = 0; k < calovector.size() && k<3; k++){
1560  fprimaryKineticEnergy[k] = calovector[k].KineticEnergy();
1561  fprimaryRange[k] = calovector[k].Range();
1562  //const std::vector< double > & dedxvec = calovector[k].dEdx();
1563  }
1564 
1565 
1566 
1567  int planenum=999;
1568  float zpos=-999;
1569 
1570  //////////Section for looking into integral dQ vs hitz position////////////
1571  //hits and calorimetry loop
1572  if(fmthm.isValid()){
1573  auto vhit=fmthm.at(fprimaryID);
1574  auto vmeta=fmthm.data(fprimaryID);
1575  for (size_t ii = 0; ii<vhit.size(); ++ii){ //loop over all meta data hit
1576  bool fBadhit = false;
1577  if (vmeta[ii]->Index() == static_cast<unsigned int>(std::numeric_limits<int>::max())){
1578  fBadhit = true;
1579  //cout<<"fBadHit"<<fBadhit<<endl;
1580  continue;
1581  }
1582  if (vmeta[ii]->Index()>=tracklist[fprimaryID]->NumberTrajectoryPoints()){
1583  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!!";
1584  }
1585  if (!tracklist[fprimaryID]->HasValidPoint(vmeta[ii]->Index())){
1586  fBadhit = true;
1587  // cout<<"had valid point "<<fBadhit<<endl;
1588  continue;
1589  }
1590 
1591  auto loc = tracklist[fprimaryID]->LocationAtPoint(vmeta[ii]->Index());
1592  // xpos=loc.X();
1593  // ypos=loc.Y();
1594  zpos=loc.Z();
1595  // cout<<"x, y, z "<<xpos<<" "<<ypos<<" "<<zpos<<endl;
1596  // cout<<"BadHit"<<fBadhit<<endl;
1597  if (fBadhit) continue; //HY::If BAD hit, skip this hit and go next
1598  if (zpos<-100) continue; //hit not on track
1599  planenum=vhit[ii]->WireID().Plane;
1600  if(planenum==2){
1601  // std::cout<<"inside loop"<<std::endl;
1602  // std::array<float,3> cnn_out=hitResults.getOutput(vhit[ii]);
1603  peakT_2.push_back(vhit[ii]->PeakTime());
1604  int_2.push_back(vhit[ii]->Integral());
1605  hitz_2.push_back(zpos);
1606  // inelscore.push_back(cnn_out[hitResults.getIndex("inel")]);
1607  // elscore.push_back(cnn_out[hitResults.getIndex("el")]);
1608  // nonescore.push_back(cnn_out[hitResults.getIndex("none")]);
1609 
1610  // std::cout<<"inel, el, none score"<<cnn_out[hitResults.getIndex("inel")]<<" "<<cnn_out[hitResults.getIndex("el")]<<" "<<cnn_out[hitResults.getIndex("none")]<<std::endl;
1611  // std::cout<<"peaktime "<<vhit[ii]->PeakTime()<<std::endl;
1612  }//planenum 2
1613  if(planenum==1){
1614  int_1.push_back(vhit[ii]->Integral());
1615  hitz_1.push_back(zpos);
1616  }//planenum 1
1617  if(planenum==0){
1618  int_0.push_back(vhit[ii]->Integral());
1619  hitz_0.push_back(zpos);
1620  }//planenum 0
1621  }//loop over vhit
1622  }//fmthm valid
1623  //hits and calorimetry loop
1624 
1625 
1626  //**************Section for dE/dx correction*****************************************************************************************************************************//
1627  //**********************************************************************************************************************************************************************//
1628  std::vector<float> Stw, Endw, Stt, Endt, Stwires, Endwires, Stticks, Endticks, TPCb, TPCcl, primwireb, primtickb, wirebuf, tickbuf, primdqb;
1629  std::vector<std::vector<float> > primwire;
1630  std::vector<std::vector<float> > clwire;
1631  std::vector<std::vector<float> > primtick;
1632  std::vector<std::vector<float> > cltick;
1633  Stw.clear(); Endw.clear(); Stt.clear(); Endt.clear(); Stwires.clear(); Endwires.clear(); Stticks.clear(); Endticks.clear(); TPCb.clear(); TPCcl.clear();
1634  float den;
1635  float numw, numt,wire_no,ticks_no;
1636  for(size_t p1=0;p1<pfplist.size();p1++){
1637  std::vector<art::Ptr<recob::Track>> trk=pftrack.at(p1);
1638  std::vector<art::Ptr<recob::Cluster>> allClusters=fmcp.at(p1);
1639  for(size_t c1=0;c1<allClusters.size();c1++){
1640  if(trk.size() && int(trk[0].key())==fprimaryID) std::cout<<"cluster tpc and plane "<<allClusters[c1]->Plane().TPC<<" "<<allClusters[c1]->Plane().Plane<<" size "<<clhit.at(allClusters[c1].key()).size()<<std::endl;
1641  if(allClusters[c1]->Plane().TPC!=1) continue;
1642  if(allClusters[c1]->Plane().Plane==0 && trk.size() && int(trk[0].key())==fprimaryID){
1643  std::vector<art::Ptr<recob::Hit>> allHits=clhit.at(allClusters[c1].key());
1644  for(size_t h1=0;h1<allHits.size();h1++){
1645  primwireb.push_back(allHits[h1]->WireID().Wire);
1646  primtickb.push_back(allHits[h1]->PeakTime());
1647  primdqb.push_back(allHits[h1]->Integral());
1648  }
1649  wireno_0.push_back(primwireb);
1650  peakTime_0.push_back(primtickb);
1651  dq_0.push_back(primdqb);
1652  primwireb.clear();
1653  primtickb.clear();
1654  primdqb.clear();
1655  }
1656  if(allClusters[c1]->Plane().Plane==1 && trk.size() && int(trk[0].key())==fprimaryID){
1657  std::vector<art::Ptr<recob::Hit>> allHits=clhit.at(allClusters[c1].key());
1658  for(size_t h1=0;h1<allHits.size();h1++){
1659  primwireb.push_back(allHits[h1]->WireID().Wire);
1660  primtickb.push_back(allHits[h1]->PeakTime());
1661  primdqb.push_back(allHits[h1]->Integral());
1662  }
1663  wireno_1.push_back(primwireb);
1664  peakTime_1.push_back(primtickb);
1665  dq_1.push_back(primdqb);
1666  primwireb.clear();
1667  primtickb.clear();
1668  primdqb.clear();
1669  }
1670 
1671  if(allClusters[c1]->Plane().Plane!=2) continue;
1672  if(trk.size() && int(trk[0].key())==fprimaryID){
1673  std::vector<art::Ptr<recob::Hit>> allHits=clhit.at(allClusters[c1].key());
1674  Stw.push_back(allClusters[c1]->StartWire());
1675  Endw.push_back(allClusters[c1]->EndWire());
1676  Stt.push_back(allClusters[c1]->StartTick());
1677  Endt.push_back(allClusters[c1]->EndTick());
1678  TPCb.push_back(allClusters[c1]->Plane().TPC);
1679  for(size_t h1=0;h1<allHits.size();h1++){
1680  primwireb.push_back(allHits[h1]->WireID().Wire);
1681  primtickb.push_back(allHits[h1]->PeakTime());
1682  primdqb.push_back(allHits[h1]->Integral());
1683  }
1684  wireno_2.push_back(primwireb);
1685  peakTime_2.push_back(primtickb);
1686  dq_2.push_back(primdqb);
1687  primwire.push_back(primwireb);
1688  primtick.push_back(primtickb);
1689  primwireb.clear();
1690  primtickb.clear();
1691  primdqb.clear();
1692  }
1693  else if(trk.size()){
1694  std::vector<art::Ptr<recob::Hit>> allHits=clhit.at(allClusters[c1].key());
1695  // std::cout<<"size of hits "<<allHits.size()<<std::endl;
1696  Stwires.push_back(allClusters[c1]->StartWire());
1697  Endwires.push_back(allClusters[c1]->EndWire());
1698  Stticks.push_back(allClusters[c1]->StartTick());
1699  Endticks.push_back(allClusters[c1]->EndTick());
1700  TPCcl.push_back(allClusters[c1]->Plane().TPC);
1701  for(size_t h2=0;h2<allHits.size();h2++){
1702  wirebuf.push_back(allHits[h2]->WireID().Wire);
1703  tickbuf.push_back(allHits[h2]->PeakTime());
1704  }
1705  clwire.push_back(wirebuf);
1706  cltick.push_back(tickbuf);
1707  wirebuf.clear();
1708  tickbuf.clear();
1709  }
1710  }
1711  }
1712  //solving the equations here
1713  for(size_t clt=0;clt<Stw.size();clt++){
1714  for(size_t cl1=0;cl1<Stwires.size();cl1++){
1715  if(TPCcl[cl1]!=TPCb[clt]) continue;
1716  den=(Stw[clt]-Endw[clt])*(Stticks[cl1]-Endticks[cl1])-(Stt[clt]-Endt[clt])*(Stwires[cl1]-Endwires[cl1]);
1717  bool testv=false;
1718  if(Stwires[cl1]==Endwires[cl1]){
1719  wire_no=Stwires[cl1];
1720  ticks_no=Stticks[cl1];
1721  testv=true;
1722  }
1723  if(den==0 && !testv) continue;
1724  if(!testv){
1725  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]);
1726  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]);
1727  wire_no=numw/den;
1728  ticks_no=numt/den;
1729  }
1730  if(((Stw[clt]-20<wire_no && Endw[clt]+20>wire_no)||(Stw[clt]+20>wire_no && Endw[clt]<wire_no-20))&&((Stt[clt]-100<ticks_no && Endt[clt]+100>ticks_no)||(Stt[clt]+100>ticks_no && Endt[clt]-100<ticks_no)) && ((Stwires[cl1]-20<wire_no && Endwires[cl1]+20>wire_no)||(Stwires[cl1]+20>wire_no && Endwires[cl1]-20<wire_no)) && ((Stticks[cl1]-100<ticks_no && Endticks[cl1]+100>ticks_no)||(Stticks[cl1]+100>ticks_no && Endticks[cl1]-100<ticks_no))){
1731  double xyzStart[3];
1732  double xyzEnd[3];
1733  unsigned int wireno=std::round(wire_no);
1734  // geo::WireID wireid(0,TPCb[clt],2,wireno);
1735  if(wireno>=0 && wireno<=479){
1736  fGeometry->WireEndPoints(0,TPCb[clt],2,wireno, xyzStart, xyzEnd);
1737  Zintersection.push_back(xyzStart[2]);
1738  timeintersection.push_back(ticks_no);
1739  // std::cout<<"intersecting "<<xyzStart[2]<<std::endl;
1740  }
1741  else{
1742  Zintersection.push_back(wireno*0.47);
1743  timeintersection.push_back(ticks_no);
1744  }
1745  //go for the second solution now
1746  float pws=wireno; //p=primary w=wire s=small b=big t=tick c=other clusters
1747  float pwb=wireno;
1748  float pts=ticks_no;
1749  float ptb=ticks_no;
1750  float smallest=wireno;
1751  float biggest=wireno;
1752  float small, big;
1753  //finding new variables for beam track
1754  for(size_t hit1=0;hit1<primwire[clt].size();hit1++){
1755  if(primwire[clt][hit1]>=wireno-30 && primwire[clt][hit1]<=wireno+30){
1756  if(primwire[clt][hit1]>=wireno-30 && primwire[clt][hit1]<=wireno){
1757  small=primwire[clt][hit1];
1758  if(small<smallest){
1759  pws=primwire[clt][hit1];
1760  pts=primtick[clt][hit1];
1761  }
1762  smallest=primwire[clt][hit1];
1763  }
1764  if(primwire[clt][hit1]<=wireno+30 && primwire[clt][hit1]>=wireno){
1765  big=primwire[clt][hit1];
1766  if(big>biggest){
1767  pwb=primwire[clt][hit1];
1768  ptb=primtick[clt][hit1];
1769  }
1770  biggest=primwire[clt][hit1];
1771  }
1772 
1773  }//if +_30 loop
1774 
1775  }//hit1 loop
1776  //finding new variables for remaining clusters
1777  float cws=wireno;
1778  float cwb=wireno;
1779  float cts=ticks_no;
1780  float ctb=ticks_no;
1781  smallest=wireno;
1782  biggest=wireno;
1783  float csmall, cbig;
1784  for(size_t hit1=0;hit1<clwire[cl1].size();hit1++){
1785  if(clwire[cl1][hit1]>=wireno-30 && clwire[cl1][hit1]<=wireno+30){
1786  if(clwire[cl1][hit1]>=wireno-30 && clwire[cl1][hit1]<=wireno){
1787  csmall=clwire[cl1][hit1];
1788  if(csmall<smallest){
1789  cws=clwire[cl1][hit1];
1790  cts=cltick[cl1][hit1];
1791  }
1792  smallest=clwire[cl1][hit1];
1793  }
1794  if(clwire[cl1][hit1]<=wireno+30 && clwire[cl1][hit1]>=wireno){
1795  cbig=clwire[cl1][hit1];
1796  if(cbig>biggest){
1797  cwb=clwire[cl1][hit1];
1798  ctb=cltick[cl1][hit1];
1799  }
1800  biggest=clwire[cl1][hit1];
1801  }
1802 
1803  }//if +_30 loop
1804  }//hit1 loop
1805  float ans[2];
1806  soln(pws, pwb, cws, cwb, pts, ptb, cts, ctb,ans);
1807  // std::cout<<"answer values "<<ans[0]<<" "<<ans[1]<<std::endl;
1808  // std::cout<<"wire number and new start end wire "<<wireno<<" st w "<<pws<<" end w "<<pwb<<" time and new times "<<ticks_no<<" "<<pts<<" "<<ptb<<" same for intersecting cluster "<<cws<<" "<<cwb<<" "<<cts<<" "<<ctb<<" hits in a cluster "<<clwire[cl1].size()<<"start wire and ticks "<<Stwires[cl1]<<" "<<Endwires[cl1]<<" "<<Stticks[cl1]<<" "<<Endticks[cl1]<<std::endl;
1809  int wir=round(ans[0]);
1810  double xyzStart1[3];
1811  double xyzEnd1[3];
1812 
1813  /* fGeometry->WireEndPoints(0,2,2,10, xyzStart1, xyzEnd1);//will remove this line
1814  stdd::cout<<"Wire start X in the postive direction "<<xyzStart1[0]<<std::endl;
1815  fGeometry->WireEndPoints(0,1,2,10, xyzStart1, xyzEnd1);//will remove this line
1816  std::cout<<"Wire start X in the negative direction "<<xyzStart1[0]<<std::endl;*/
1817 
1818  if(((Stw[clt]-5<ans[0] && Endw[clt]+5>ans[0])||(Stw[clt]+5>ans[0] && Endw[clt]<ans[0]-5))&&((Stt[clt]-50<ans[1] && Endt[clt]+50>ans[1])||(Stt[clt]+50>ans[1] && Endt[clt]-50<ans[1])) && ((Stwires[cl1]-5<ans[0] && Endwires[cl1]+5>ans[0])||(Stwires[cl1]+5>ans[0] && Endwires[cl1]-5<ans[0])) && ((Stticks[cl1]-50<ans[1] && Endticks[cl1]+50>ans[1])||(Stticks[cl1]+50>ans[1] && Endticks[cl1]-50<ans[1]))){
1819  if(wir>=0 && wir<=479){
1820  fGeometry->WireEndPoints(0,TPCb[clt],2,wir, xyzStart1, xyzEnd1);
1821  Zintersection1.push_back(xyzStart1[2]);
1822  timeintersection1.push_back(ans[1]);
1823  std::cout<<"wire X position is here "<<xyzStart1[0]<<std::endl;
1824  std::cout<<"wire X position End is here "<<xyzEnd1[0]<<std::endl;
1825  }
1826  else{
1827  Zintersection1.push_back(wire_no*0.47);
1828  timeintersection1.push_back(ticks_no);
1829  ans[0]=wire_no;
1830  ans[1]=ticks_no;
1831  }
1832  }
1833  else{
1834  Zintersection1.push_back(99999);
1835  timeintersection1.push_back(99999);
1836  }
1837 
1838  // if((pws<=ans[0] && pwb>=ans[0])&&(cws<=ans[0] && cwb>=ans[0])&&(pts<=ans[1] && ptb>=ans[1])&&(cts<=ans[1] && ctb>=ans[1])){
1839  if(((Stw[clt]<=ans[0] && Endw[clt]>=ans[0])||(Stw[clt]>=ans[0] && Endw[clt]<=ans[0]))&&((Stt[clt]<=ans[1] && Endt[clt]>ans[1])||(Stt[clt]>ans[1] && Endt[clt]<ans[1])) && ((Stwires[cl1]<ans[0] && Endwires[cl1]>ans[0])||(Stwires[cl1]>ans[0] && Endwires[cl1]<ans[0])) && ((Stticks[cl1]<ans[1] && Endticks[cl1]>ans[1])||(Stticks[cl1]>ans[1] && Endticks[cl1]<ans[1]))){
1840  deltaZint.push_back(0);
1841  deltatimeint.push_back(0);
1842  }
1843  else if((Stw[clt]<=ans[0] && Endw[clt]>=ans[0])||(Stw[clt]>=ans[0] && Endw[clt]<=ans[0])){
1844  deltaZint.push_back(abs(std::min(abs(ans[0]-Stwires[cl1]),abs(ans[0]-Endwires[cl1]))));
1845  deltatimeint.push_back(abs(std::min(abs(ans[1]-Stticks[cl1]),abs(ans[1]-Endticks[cl1]))));
1846  }
1847  else if((Stwires[cl1]<=ans[0] && Endwires[cl1]>=ans[0])||(Stwires[cl1]>=ans[0] && Endwires[cl1]<=ans[0])){
1848  deltaZint.push_back(abs(std::min(abs(ans[0]-Stw[clt]),abs(ans[0]-Endw[clt]))));
1849  deltatimeint.push_back(abs(std::min(abs(ans[1]-Stt[clt]),abs(ans[1]-Endt[clt]))));
1850  }
1851  else{
1852  deltaZint.push_back(99999);
1853  deltatimeint.push_back(99999);
1854  }
1855 
1856  }//if solution lies within wire coordinates
1857  // }//if wire coordinates
1858  }//cl1 loop
1859  }//clt
1860  primwire.clear();
1861  clwire.clear();
1862  primtick.clear();
1863  cltick.clear();
1864 
1865  ////*****************section for dE/dx correction*******************************////
1866 
1867 
1868  }//this track
1869 
1870 
1871 
1872 
1873  //*********************************************//
1874  //*********************************************//
1875  ///////////Showere case and daughter trajectories
1876  else if(thisShower != 0x0){
1877  fisprimarytrack = 0;
1878  fisprimaryshower = 1;
1879 
1880  fprimaryID = thisShower->ID();
1881  fprimaryLength = thisShower->Length();
1882  fprimaryShowerBestPlane = thisShower->best_plane();
1883  fprimaryOpeningAngle = thisShower->OpenAngle();
1884  fprimaryStartPosition[0] = thisShower->ShowerStart().X();
1885  fprimaryStartPosition[1] = thisShower->ShowerStart().Y();
1886  fprimaryStartPosition[2] = thisShower->ShowerStart().Z();
1887  fprimaryStartDirection[0] = thisShower->Direction().X();
1888  fprimaryStartDirection[1] = thisShower->Direction().Y();
1889  fprimaryStartDirection[2] = thisShower->Direction().Z();
1890  if( (thisShower->Energy()).size() > 0 )
1891  fprimaryShowerEnergy = thisShower->Energy()[0]; // thisShower->best_plane()
1892  if( (thisShower->MIPEnergy()).size() > 0 )
1893  fprimaryShowerMIPEnergy = thisShower->MIPEnergy()[0];
1894  if( (thisShower->dEdx()).size() > 0 )
1895  fprimaryShowerdEdx = thisShower->dEdx()[0];
1896  }
1897  else{
1898  std::cout << "INFO::Primary pfParticle is not track or shower. Skip!" << std::endl;
1899  continue;
1900  }
1901 
1902  // Find the particle vertex. We need the tracker tag here because we need to do a bit of
1903  // additional work if the PFParticle is track-like to find the vertex.
1904  const TVector3 vtx = pfpUtil.GetPFParticleVertex(*particle,evt,fPFParticleTag,fTrackerTag);
1905  fvertex[0] = vtx.X(); fvertex[1] = vtx.Y(); fvertex[2] = vtx.Z();
1906 
1907  // Now we can look for the interaction point of the particle if one exists, i.e where the particle
1908  // scatters off an argon nucleus. Shower-like objects won't have an interaction point, so we can
1909  // check this by making sure we get a sensible position
1910  const TVector3 interactionVtx = pfpUtil.GetPFParticleSecondaryVertex(*particle,evt,fPFParticleTag,fTrackerTag);
1911  fsecvertex[0] = interactionVtx.X(); fsecvertex[1] = interactionVtx.Y(); fsecvertex[2] = interactionVtx.Z();
1912 
1913  // Maximum number of daugthers to be processed
1914  if(particle->NumDaughters() > NMAXDAUGTHERS)
1915  std::cout << "INFO::Number of daughters is " << particle->NumDaughters() << ". Only the first NMAXDAUGTHERS are processed." << std::endl;
1916 
1917  // Let's get the daughter PFParticles... we can do this simply without the utility
1918  for(const int daughterID : particle->Daughters()){
1919  // Daughter ID is the element of the original recoParticle vector
1920  const recob::PFParticle *daughterParticle = &(recoParticles->at(daughterID));
1921 
1922  const recob::Track* daughterTrack = pfpUtil.GetPFParticleTrack(*daughterParticle,evt,fPFParticleTag,fTrackerTag);
1923  const recob::Shower* daughterShower = pfpUtil.GetPFParticleShower(*daughterParticle,evt,fPFParticleTag,fShowerTag);
1924 
1925  if(daughterTrack != 0x0){
1926  fisdaughtertrack[fNDAUGHTERS] = 1;
1927  fisdaughtershower[fNDAUGHTERS] = 0;
1928  fdaughterTheta[fNDAUGHTERS] = daughterTrack->Theta();
1929  fdaughterPhi[fNDAUGHTERS] = daughterTrack->Phi();
1930  fdaughterLength[fNDAUGHTERS] = daughterTrack->Length();
1931  fdaughterMomentum[fNDAUGHTERS] = daughterTrack->StartMomentum();
1932  fdaughterEndMomentum[fNDAUGHTERS] = daughterTrack->EndMomentum();
1933  fdaughterStartPosition[fNDAUGHTERS][0] = daughterTrack->Trajectory().Start().X();
1934  fdaughterStartPosition[fNDAUGHTERS][1] = daughterTrack->Trajectory().Start().Y();
1935  fdaughterStartPosition[fNDAUGHTERS][2] = daughterTrack->Trajectory().Start().Z();
1936  fdaughterEndPosition[fNDAUGHTERS][0] = daughterTrack->Trajectory().End().X();
1937  fdaughterEndPosition[fNDAUGHTERS][1] = daughterTrack->Trajectory().End().Y();
1938  fdaughterEndPosition[fNDAUGHTERS][2] = daughterTrack->Trajectory().End().Z();
1939  fdaughterStartDirection[fNDAUGHTERS][0] = daughterTrack->Trajectory().StartDirection().X();
1940  fdaughterStartDirection[fNDAUGHTERS][1] = daughterTrack->Trajectory().StartDirection().Y();
1941  fdaughterStartDirection[fNDAUGHTERS][2] = daughterTrack->Trajectory().StartDirection().Z();
1942  fdaughterEndDirection[fNDAUGHTERS][0] = daughterTrack->Trajectory().EndDirection().X();
1943  fdaughterEndDirection[fNDAUGHTERS][1] = daughterTrack->Trajectory().EndDirection().Y();
1944  fdaughterEndDirection[fNDAUGHTERS][2] = daughterTrack->Trajectory().EndDirection().Z();
1945 
1946  fdaughterMomentumByRangeMuon[fNDAUGHTERS] = trmom.GetTrackMomentum(daughterTrack->Length(),13);
1947  fdaughterMomentumByRangeProton[fNDAUGHTERS] = trmom.GetTrackMomentum(daughterTrack->Length(),2212);
1948 
1949  std::vector<anab::Calorimetry> daughtercalovector = trackUtil.GetRecoTrackCalorimetry(*daughterTrack, evt, fTrackerTag, fCalorimetryTag);
1950  if(daughtercalovector.size() != 3)
1951  std::cerr << "WARNING::Calorimetry vector size for daughter is = " << daughtercalovector.size() << std::endl;
1952 
1953  for(size_t k = 0; k < daughtercalovector.size() && k<3; k++){
1954  fdaughterKineticEnergy[fNDAUGHTERS][k] = daughtercalovector[k].KineticEnergy();
1955  fdaughterRange[fNDAUGHTERS][k] = daughtercalovector[k].Range();
1956  }
1957 
1958  // Get the true mc particle
1959  const simb::MCParticle* mcdaughterparticle = truthUtil.GetMCParticleFromRecoTrack(clockData, *daughterTrack, evt, fTrackerTag);
1960  if(mcdaughterparticle != 0x0){
1961  fdaughter_truth_TrackId[fNDAUGHTERS] = mcdaughterparticle->TrackId();
1962  fdaughter_truth_Pdg[fNDAUGHTERS] = mcdaughterparticle->PdgCode();
1963  fdaughter_truth_StartPosition[fNDAUGHTERS][0] = mcdaughterparticle->Vx();
1964  fdaughter_truth_StartPosition[fNDAUGHTERS][1] = mcdaughterparticle->Vy();
1965  fdaughter_truth_StartPosition[fNDAUGHTERS][2] = mcdaughterparticle->Vz();
1966  fdaughter_truth_StartPosition[fNDAUGHTERS][3] = mcdaughterparticle->T();
1967  fdaughter_truth_EndPosition[fNDAUGHTERS][0] = mcdaughterparticle->EndX();
1968  fdaughter_truth_EndPosition[fNDAUGHTERS][1] = mcdaughterparticle->EndY();
1969  fdaughter_truth_EndPosition[fNDAUGHTERS][2] = mcdaughterparticle->EndZ();
1970  fdaughter_truth_EndPosition[fNDAUGHTERS][3] = mcdaughterparticle->EndT();
1971  fdaughter_truth_P[fNDAUGHTERS] = mcdaughterparticle->P();
1972  fdaughter_truth_Momentum[fNDAUGHTERS][0] = mcdaughterparticle->Px();
1973  fdaughter_truth_Momentum[fNDAUGHTERS][1] = mcdaughterparticle->Py();
1974  fdaughter_truth_Momentum[fNDAUGHTERS][2] = mcdaughterparticle->Pz();
1975  fdaughter_truth_Momentum[fNDAUGHTERS][3] = mcdaughterparticle->E();
1976  fdaughter_truth_Pt[fNDAUGHTERS] = mcdaughterparticle->Pt();
1977  fdaughter_truth_Mass[fNDAUGHTERS] = mcdaughterparticle->Mass();
1978  fdaughter_truth_EndMomentum[fNDAUGHTERS][0] = mcdaughterparticle->EndPx();
1979  fdaughter_truth_EndMomentum[fNDAUGHTERS][1] = mcdaughterparticle->EndPy();
1980  fdaughter_truth_EndMomentum[fNDAUGHTERS][2] = mcdaughterparticle->EndPz();
1981  fdaughter_truth_EndMomentum[fNDAUGHTERS][3] = mcdaughterparticle->EndE();
1982  fdaughter_truth_Theta[fNDAUGHTERS] = mcdaughterparticle->Momentum().Theta();
1983  fdaughter_truth_Phi[fNDAUGHTERS] = mcdaughterparticle->Momentum().Phi();
1984  fdaughter_truth_Process[fNDAUGHTERS] = int(mcdaughterparticle->Trajectory().ProcessToKey(mcdaughterparticle->Process()));
1985  std::cout << "Daughter Process = " << (mcdaughterparticle->Process()).c_str()
1986  << " , mother = " << mcdaughterparticle->Mother()
1987  << std::endl;
1988  }
1989  }
1990  else if(daughterShower != 0x0){
1991  fisdaughtertrack[fNDAUGHTERS] = 0;
1992  fisdaughtershower[fNDAUGHTERS] = 1;
1993  fdaughterLength[fNDAUGHTERS] = daughterShower->Length();
1994  fdaughterShowerBestPlane[fNDAUGHTERS] = daughterShower->best_plane();
1995  fdaughterOpeningAngle[fNDAUGHTERS] = daughterShower->OpenAngle();
1996  fdaughterStartPosition[fNDAUGHTERS][0] = daughterShower->ShowerStart().X();
1997  fdaughterStartPosition[fNDAUGHTERS][1] = daughterShower->ShowerStart().Y();
1998  fdaughterStartPosition[fNDAUGHTERS][2] = daughterShower->ShowerStart().Z();
1999  fdaughterStartDirection[fNDAUGHTERS][0] = daughterShower->Direction().X();
2000  fdaughterStartDirection[fNDAUGHTERS][1] = daughterShower->Direction().Y();
2001  fdaughterStartDirection[fNDAUGHTERS][2] = daughterShower->Direction().Z();
2002  if( (daughterShower->Energy()).size() > 0 )
2003  fdaughterShowerEnergy[fNDAUGHTERS] = daughterShower->Energy()[0]; // thisShower->best_plane()
2004  if( (daughterShower->MIPEnergy()).size() > 0 )
2005  fdaughterShowerMIPEnergy[fNDAUGHTERS] = daughterShower->MIPEnergy()[0];
2006  if( (daughterShower->dEdx()).size() > 0 )
2007  fdaughterShowerdEdx[fNDAUGHTERS] = daughterShower->dEdx()[0];
2008  }
2009  else{
2010  std::cout << "INFO::Daughter pfParticle is not track or shower. Skip!" << std::endl;
2011  continue;
2012  }
2013 
2014  fdaughterID[fNDAUGHTERS] = daughterID;
2015  // NHits associated with this pfParticle
2016  fdaughterNHits[fNDAUGHTERS] = (pfpUtil.GetPFParticleHits(*daughterParticle,evt,fPFParticleTag)).size();
2017  // T0
2018  // std::vector<anab::T0> pfdaughterT0vec = pfpUtil.GetPFParticleT0(*daughterParticle,evt,fPFParticleTag);
2019  // if(!pfT0vec.empty())
2020  // fdaughterT0[fNDAUGHTERS] = pfdaughterT0vec[0].Time();
2021 
2022  fNDAUGHTERS++;
2023 
2024  // Only process NMAXDAUGTHERS
2025  if(fNDAUGHTERS > NMAXDAUGTHERS) break;
2026 
2027  }
2028 
2029  break;
2030  }//particle loop pfparticle
2031 
2032  // Fill trees
2033  if(beamTriggerEvent)
2034  fPandoraBeam->Fill();
2035 
2036  //fPandoraCosmics->Fill();
2037  }//analyzer
2038 
2039  void protoana::mcsXsection::endJob(){
2040 
2041  }
2042 
2043  void protoana::mcsXsection::FillCosmicsTree(art::Event const & evt, std::string pfParticleTag){
2044 
2045  // To fill
2046 
2047  }
2048 
2049 
2050 
2051  void protoana::mcsXsection::Initialise(){
2052 
2053  fRun = -999;
2054  fSubRun = -999;
2055  fevent = -999;
2056  fTimeStamp = -999.0;
2057 
2058  for(int k=0; k < 3; k++){
2059  fvertex[k] = -999.0;
2060  fsecvertex[k] = -999.0;
2061  fprimaryEndPosition[k] = -999.0;
2062  fprimaryStartPosition[k] = -999.0;
2063  fprimaryEndDirection[k] = -999.0;
2064  fprimaryStartDirection[k] = -999.0;
2065  fprimaryKineticEnergy[k] = -999.0;
2066  fprimaryRange[k] = -999.0;
2067  }
2068  fbeamtrackMomentum = -999.0;
2069  fbeamtrackEnergy = 999.0;
2070  fbeamtrackPdg = -9999;
2071  fbeamtrackTime = -999.0;
2072  fbeamtrackID = -999;
2073  for(int l=0; l < 3; l++){
2074  fbeamtrackP[l] = -999.0;
2075  fbeamtrackPos[l] = -999.0;
2076  fbeamtrackDir[l] = -999.0;
2077  }
2078 
2079  //NumberBeamTrajectoryPoints=0;
2080  beamtrk_x.clear();
2081  beamtrk_y.clear();
2082  beamtrk_z.clear();
2083  beamtrk_Px.clear();
2084  beamtrk_Py.clear();
2085  beamtrk_Pz.clear();
2086  beamtrk_Eng.clear();
2087  beamtrk_z_wire.clear();
2088  beamtrk_z_tpc.clear();
2089 
2090 
2091  beamMomentum_spec.clear();
2092  beamPosx_spec.clear();
2093  beamPosy_spec.clear();
2094  beamPosz_spec.clear();
2095  beamDirx_spec.clear();
2096  beamDiry_spec.clear();
2097  beamDirz_spec.clear();
2098 
2099 
2100 
2101  fisprimarytrack = 0;
2102  fisprimaryshower = 0;
2103  fNDAUGHTERS = 0;
2104 
2105  fprimaryBDTScore = -999.0;
2106  fprimaryNHits = -999;
2107  fprimaryTheta = -999.0;
2108  fprimaryPhi = -999.0;
2109  fprimaryLength = -999.0;
2110  fprimaryMomentum = -999.0;
2111  fprimaryEndMomentum = -999.0;
2112  fprimaryOpeningAngle = -999.0;
2113  fprimaryShowerBestPlane = -999;
2114  fprimaryShowerEnergy = -999;
2115  fprimaryShowerMIPEnergy = -999;
2116  fprimaryShowerdEdx = -999;
2117  fprimaryID = -999;
2118  fprimaryMomentumByRangeProton = -999.0;
2119  fprimaryMomentumByRangeMuon = -999.0;
2120 
2121  fprimary_truth_TrackId = -999;
2122  fprimary_truth_Pdg = -999;
2123  ftruthpdg=-999;
2124  fprimary_truth_P = -999.0;
2125  fprimary_truth_Pt = -999.0;
2126  fprimary_truth_Mass = -999.0;
2127  fprimary_truth_Theta = -999.0;
2128  fprimary_truth_Phi = -999.0;
2129  fprimary_truth_Process = -999;
2130  fprimary_truth_Isbeammatched = -999;
2131  fprimary_truth_NDaughters = -999;
2132  fprimary_truth_tracklength=-9999;
2133  fprimary_truth_EndProcess="";
2134  truth_last_process="";
2135  for(int l=0; l < 4; l++){
2136  fprimary_truth_StartPosition[l] = -999.0;
2137  fprimary_truth_EndPosition[l] = -999.0;
2138  fprimary_truth_Momentum[l] = -999.0;
2139  fprimary_truth_EndMomentum[l] = -999.0;
2140  }
2141 
2142  for(int k=0; k < NMAXDAUGTHERS; k++){
2143  fisdaughtertrack[k] = -999;
2144  fisdaughtershower[k] = -999;
2145  fdaughterNHits[k] = -999;
2146  fdaughterTheta[k] = -999.0;
2147  fdaughterPhi[k] = -999.0;
2148  fdaughterLength[k] = -999.0;
2149  fdaughterMomentum[k] = -999.0;
2150  fdaughterEndMomentum[k] = -999.0;
2151  for(int l=0; l < 3; l++){
2152  fdaughterEndPosition[k][l] = -999.0;
2153  fdaughterStartPosition[k][l] = -999.0;
2154  fdaughterEndDirection[k][l] = -999.0;
2155  fdaughterStartDirection[k][l] = -999.0;
2156  fdaughterKineticEnergy[k][l] = -999.0;
2157  fdaughterRange[k][l] = -999.0;
2158  }
2159  fdaughterOpeningAngle[k] = -999.0;
2160  fdaughterShowerBestPlane[k] = -999;
2161  fdaughterShowerEnergy[k] = -999;
2162  fdaughterShowerMIPEnergy[k] = -999;
2163  fdaughterShowerdEdx[k] = -999;
2164  fdaughterMomentumByRangeProton[k] = -999.0;
2165  fdaughterMomentumByRangeMuon[k] = -999.0;
2166  fdaughterID[k] = -999;
2167  // fdaughterT0[k] = -999;
2168 
2169  fdaughter_truth_TrackId[k] = -999;
2170  fdaughter_truth_Pdg[k] = -999;
2171  fdaughter_truth_P[k] = -999.0;
2172  fdaughter_truth_Pt[k] = -999.0;
2173  fdaughter_truth_Mass[k] = -999.0;
2174  fdaughter_truth_Theta[k] = -999.0;
2175  fdaughter_truth_Phi[k] = -999.0;
2176  fdaughter_truth_Process[k] = -999;
2177  for(int l=0; l < 4; l++){
2178  fdaughter_truth_StartPosition[k][l] = -999.0;
2179  fdaughter_truth_EndPosition[k][l] = -999.0;
2180  fdaughter_truth_Momentum[k][l] = -999.0;
2181  fdaughter_truth_EndMomentum[k][l] = -999.0;
2182  }
2183  }
2184  int_2.clear();
2185  hitz_2.clear();
2186  // inelscore.clear();
2187  // elscore.clear();
2188  // nonescore.clear();
2189  int_1.clear();
2190  hitz_1.clear();
2191  int_0.clear();
2192  hitz_0.clear();
2193  peakT_2.clear();
2194  primtrk_dqdx.clear();
2195  primtrk_resrange.clear();
2196  primtrk_dedx.clear();
2197  primtrk_range.clear();
2198  primtrk_hitx.clear();
2199  primtrk_hity.clear();
2200  primtrk_hitz.clear();
2201  primtrk_hitz_wire.clear();
2202  primtrk_hitz_tpc.clear();
2203 
2204  primtrk_pitch.clear();
2205  primtrk_truth_Z.clear();
2206  primtrk_truth_Z_wire.clear();
2207  primtrk_truth_Z_tpc.clear();
2208 
2209 
2210  primtrk_truth_Eng.clear();
2211  primtrk_truth_trkide.clear();
2212  interactionX.clear();
2213  interactionY.clear();
2214  interactionZ.clear();
2215  interactionU.clear();
2216  interactionV.clear();
2217  interactionW.clear();
2218  interactionT.clear();
2219 
2220 
2221  interactionProcesslist.clear();
2222  interactionAngles.clear();
2223  interactionAnglesUT.clear();
2224  interactionAnglesVT.clear();
2225  interactionAnglesZT.clear();
2226  Zintersection.clear();
2227  timeintersection.clear();
2228  Zintersection1.clear();
2229  timeintersection1.clear();
2230  deltaZint.clear();
2231  deltatimeint.clear();
2232  wireno_0.clear();
2233  wireno_1.clear();
2234  wireno_2.clear();
2235  dq_0.clear();
2236  dq_1.clear();
2237  dq_2.clear();
2238  peakTime_0.clear();
2239  peakTime_1.clear();
2240  peakTime_2.clear();
2241  //Michel tagging
2242  endhitssecondary.clear();
2243  secondarystartx.clear();
2244  secondaryendx.clear();
2245  secondarystarty.clear();
2246  secondaryendy.clear();
2247  secondarystartz.clear();
2248  secondaryendz.clear();
2249  dQmichel.clear();
2250  primsectheta.clear();
2251  dQtrackbegin.clear();
2252  dQtrackend.clear();
2253  tracklengthsecondary.clear();
2254  MtrackID.clear();
2255  trackscore.clear();
2256  emscore.clear();
2257  michelscore.clear();
2258  nonescore.clear();
2259  //geant4 reweight
2260  n_g4rw=0;
2261  g4rw_react_st=1;
2262  g4rw_d_react=0;
2263  g4rw_elast_st=1;
2264  g4rw_d_elast=0;
2265 
2266  g4rw_set_weights.clear();
2267  g4rw_elast.clear();
2268  g4rw_react.clear();
2269 
2270  }
2271 
2272  DEFINE_ART_MODULE(protoana::mcsXsection)
simb::MCParticle::E
double E(const int i=0) const
Definition: MCParticle.h:233
geo::GeometryCore::WireCoordinate
geo::Length_t WireCoordinate(double YPos, double ZPos, geo::PlaneID const &planeid) const
Returns the index of the nearest wire to the specified position.
Definition: GeometryCore.cxx:1128
recob::Shower::best_plane
int best_plane() const
Definition: Shower.h:200
ValidateOpDetReco.end
end
while True: pbar.update(maxval-len(onlies[E][S])) #print iS, "/", len(onlies[E][S]) found = False for...
Definition: ValidateOpDetReco.py:548
protoana::mcsXsection::MtrackID
std::vector< std::vector< int > > MtrackID
Definition: mcsXsection_module.cc:361
protoana::mcsXsection::wireno_1
std::vector< std::vector< float > > wireno_1
Definition: mcsXsection_module.cc:338
protoana::mcsXsection::wireno_0
std::vector< std::vector< float > > wireno_0
Definition: mcsXsection_module.cc:337
protoana::mcsXsection::tracklengthsecondary
std::vector< std::vector< double > > tracklengthsecondary
Definition: mcsXsection_module.cc:360
m_pi
double m_pi
Definition: mcsXsection_module.cc:94
protoana::mcsXsection::fdaughterMomentum
double fdaughterMomentum[NMAXDAUGTHERS]
Definition: mcsXsection_module.cc:378
protoana::mcsXsection::g4rw_react
std::vector< double > g4rw_react
Definition: mcsXsection_module.cc:425
simb::MCParticle::NumberTrajectoryPoints
unsigned int NumberTrajectoryPoints() const
Definition: MCParticle.h:218
protoana::mcsXsection::beamtrk_Pz
std::vector< float > beamtrk_Pz
Definition: mcsXsection_module.cc:213
recob::Shower::ShowerStart
const TVector3 & ShowerStart() const
Definition: Shower.h:192
simb::MCParticle::Position
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:219
protoana::mcsXsection::fdaughter_truth_TrackId
int fdaughter_truth_TrackId[NMAXDAUGTHERS]
Definition: mcsXsection_module.cc:396
recob::Track::EndMomentum
double EndMomentum() const
Definition: Track.h:144
simb::MCParticle::PdgCode
int PdgCode() const
Definition: MCParticle.h:212
art::Timestamp::timeLow
constexpr std::uint32_t timeLow() const
Definition: Timestamp.h:29
protoana::mcsXsection::fdaughter_truth_Process
int fdaughter_truth_Process[NMAXDAUGTHERS]
Definition: mcsXsection_module.cc:407
protoana::mcsXsection::peakTime_2
std::vector< std::vector< float > > peakTime_2
Definition: mcsXsection_module.cc:342
protoana::mcsXsection::operator=
mcsXsection & operator=(mcsXsection const &)=delete
protoana::mcsXsection::interactionAngles
std::vector< double > interactionAngles
Definition: mcsXsection_module.cc:308
simb::MCParticle::Py
double Py(const int i=0) const
Definition: MCParticle.h:231
protoana::mcsXsection::peakT_2
std::vector< float > peakT_2
Definition: mcsXsection_module.cc:215
prim_energy
double prim_energy
Definition: mcsXsection_module.cc:95
protoana::mcsXsection::endhitssecondary
std::vector< std::vector< int > > endhitssecondary
Definition: mcsXsection_module.cc:349
protoana::mcsXsection::peakTime_1
std::vector< std::vector< float > > peakTime_1
Definition: mcsXsection_module.cc:341
protoana::mcsXsection::beamDirx_spec
std::vector< double > beamDirx_spec
Definition: mcsXsection_module.cc:228
protoana::mcsXsection::ParSet
std::vector< fhicl::ParameterSet > ParSet
Definition: mcsXsection_module.cc:431
protoana::ProtoDUNETruthUtils::GetGeantGoodParticle
const simb::MCParticle * GetGeantGoodParticle(const simb::MCTruth &genTruth, const art::Event &evt) const
Definition: ProtoDUNETruthUtils.cxx:860
protoana::mcsXsection::beamDirz_spec
std::vector< double > beamDirz_spec
Definition: mcsXsection_module.cc:230
reco_momentum_tuples.x3
x3
Definition: reco_momentum_tuples.py:84
simb::MCParticle::EndE
double EndE() const
Definition: MCParticle.h:244
result
static QCString result
Definition: fortranscanner.cpp:56614
simb::MCParticle::EndZ
double EndZ() const
Definition: MCParticle.h:228
protoana::mcsXsection::fdaughterID
int fdaughterID[NMAXDAUGTHERS]
Definition: mcsXsection_module.cc:393
protoana::mcsXsection::beamMomentum_spec
std::vector< double > beamMomentum_spec
Definition: mcsXsection_module.cc:223
protoana::mcsXsection::dQmichel
std::vector< std::vector< double > > dQmichel
Definition: mcsXsection_module.cc:356
protoana::mcsXsection::trackscore
std::vector< std::vector< double > > trackscore
Definition: mcsXsection_module.cc:363
protoana::mcsXsection::primtrk_hitz_wire
std::vector< std::vector< int > > primtrk_hitz_wire
Definition: mcsXsection_module.cc:326
recob::Shower::Length
double Length() const
Definition: Shower.h:201
protoana::ProtoDUNEPFParticleUtils::GetPFParticleShower
const recob::Shower * GetPFParticleShower(const recob::PFParticle &particle, art::Event const &evt, const std::string particleLabel, const std::string showerLabel) const
Get the shower associated to this particle. Returns a null pointer if not found.
Definition: ProtoDUNEPFParticleUtils.cxx:496
art::DataViewImpl::getHandle
Handle< PROD > getHandle(SelectorBase const &) const
Definition: DataViewImpl.h:382
geo::View_t
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
protoana::mcsXsection::primtrk_dqdx
std::vector< std::vector< double > > primtrk_dqdx
Definition: mcsXsection_module.cc:319
string
std::string string
Definition: nybbler.cc:12
protoana::mcsXsection::beamPosx_spec
std::vector< double > beamPosx_spec
Definition: mcsXsection_module.cc:224
simb::MCParticle::Trajectory
const simb::MCTrajectory & Trajectory() const
Definition: MCParticle.h:253
protoana::mcsXsection::fdaughterEndMomentum
double fdaughterEndMomentum[NMAXDAUGTHERS]
Definition: mcsXsection_module.cc:379
simb::MCParticle::Mother
int Mother() const
Definition: MCParticle.h:213
truth_ana.theta
theta
Definition: truth_ana.py:138
recob::Track::Trajectory
const recob::TrackTrajectory & Trajectory() const
Access to the stored recob::TrackTrajectory.
Definition: Track.h:98
protoana::ProtoDUNETrackUtils::GetRecoTrackCalorimetry
std::vector< anab::Calorimetry > GetRecoTrackCalorimetry(const recob::Track &track, art::Event const &evt, const std::string trackModule, const std::string caloModule) const
Get the Calorimetry(s) from a given reco track.
Definition: ProtoDUNETrackUtils.cxx:72
angle2d
double angle2d(double x0, double y0, double x1, double y1, double x2, double y2)
Definition: mcsXsection_module.cc:114
protoana::mcsXsection::primtrk_truth_Eng
std::vector< std::vector< double > > primtrk_truth_Eng
Definition: mcsXsection_module.cc:335
simb::MCTrajectory::at
const value_type & at(const size_type i) const
Definition: MCTrajectory.h:175
protoana::mcsXsection::peakTime_0
std::vector< std::vector< float > > peakTime_0
Definition: mcsXsection_module.cc:340
art::Timestamp::timeHigh
constexpr std::uint32_t timeHigh() const
Definition: Timestamp.h:34
simb::MCParticle::Mass
double Mass() const
Definition: MCParticle.h:239
cet::pow
constexpr T pow(T x)
Definition: pow.h:72
geo::PlaneID
The data type to uniquely identify a Plane.
Definition: geo_types.h:472
geo::CryostatID::isValid
bool isValid
Whether this ID points to a valid element.
Definition: geo_types.h:211
simb::MCTrajectory::KeyToProcess
std::string KeyToProcess(unsigned char const &key) const
Definition: MCTrajectory.cxx:117
protoana::ProtoDUNEPFParticleUtils::GetPFParticleHits
const std::vector< const recob::Hit * > GetPFParticleHits(const recob::PFParticle &particle, art::Event const &evt, const std::string particleLabel) const
Get the hits associated to the PFParticle.
Definition: ProtoDUNEPFParticleUtils.cxx:564
recob::Shower::Energy
const std::vector< double > & Energy() const
Definition: Shower.h:195
simb::MCParticle::Px
double Px(const int i=0) const
Definition: MCParticle.h:230
protoana::mcsXsection::beamtrk_x
std::vector< float > beamtrk_x
Definition: mcsXsection_module.cc:205
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::mcsXsection::Zintersection1
std::vector< double > Zintersection1
Definition: mcsXsection_module.cc:313
protoana::mcsXsection::fdaughterEndPosition
double fdaughterEndPosition[NMAXDAUGTHERS][3]
Definition: mcsXsection_module.cc:380
TrackHitMeta.h
Class to keep data related to recob::Hit associated with recob::Track.
protoana::mcsXsection::primtrk_hitx
std::vector< std::vector< double > > primtrk_hitx
Definition: mcsXsection_module.cc:323
protoana::mcsXsection::emscore
std::vector< std::vector< double > > emscore
Definition: mcsXsection_module.cc:364
anab::FVectorReader::getIndex
int getIndex(const std::string &name) const
Index of column with given name, or -1 if name not found.
Definition: MVAReader.h:82
cheat::BackTrackerService::TrackIdToSimIDEs_Ps
std::vector< const sim::IDE * > TrackIdToSimIDEs_Ps(int const &id) const
Definition: BackTrackerService_service.cc:134
protoana::mcsXsection::primtrk_truth_Z_wire
std::vector< std::vector< int > > primtrk_truth_Z_wire
Definition: mcsXsection_module.cc:331
protoana::G4ReweightUtils::GetNTrajWeightFromSetPars
double GetNTrajWeightFromSetPars(const std::vector< G4ReweightTraj * > &trajs, G4MultiReweighter &rw)
Definition: G4ReweightUtils.cxx:286
std
STL namespace.
geo::WireID::Wire
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:580
protoana::mcsXsection::deltaZint
std::vector< double > deltaZint
Definition: mcsXsection_module.cc:316
protoana::mcsXsection::beamtrk_Px
std::vector< float > beamtrk_Px
Definition: mcsXsection_module.cc:211
protoana::mcsXsection::dQtrackbegin
std::vector< std::vector< double > > dQtrackbegin
Definition: mcsXsection_module.cc:358
Index
unsigned int Index
Definition: sspmapmaker_v1.c:14
protoana::mcsXsection::primtrk_truth_trkide
std::vector< std::vector< double > > primtrk_truth_trkide
Definition: mcsXsection_module.cc:336
art::EDAnalyzer::EDAnalyzer
EDAnalyzer(fhicl::ParameterSet const &pset)
Definition: EDAnalyzer.h:25
simb::MCParticle::Process
std::string Process() const
Definition: MCParticle.h:215
reco_momentum_tuples.x2
x2
Definition: reco_momentum_tuples.py:83
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::mcsXsection::fdaughterMomentumByRangeProton
double fdaughterMomentumByRangeProton[NMAXDAUGTHERS]
Definition: mcsXsection_module.cc:389
protoana::mcsXsection::fisdaughtershower
int fisdaughtershower[NMAXDAUGTHERS]
Definition: mcsXsection_module.cc:373
protoana::mcsXsection::dataUtil
protoana::ProtoDUNEDataUtils dataUtil
Definition: mcsXsection_module.cc:160
protoana::mcsXsection::fdaughterRange
double fdaughterRange[NMAXDAUGTHERS][3]
Definition: mcsXsection_module.cc:392
l
static QStrList * l
Definition: config.cpp:1044
protoana::mcsXsection::fisdaughtertrack
int fisdaughtertrack[NMAXDAUGTHERS]
Definition: mcsXsection_module.cc:372
protoana::ProtoDUNEBeamCuts::IsBeamlike
bool IsBeamlike(const recob::Track &, const art::Event &, std::string)
protoana::mcsXsection::wireno_2
std::vector< std::vector< float > > wireno_2
Definition: mcsXsection_module.cc:339
simb::MCParticle::NumberDaughters
int NumberDaughters() const
Definition: MCParticle.h:217
Geometry.h
art framework interface to geometry description
protoana::mcsXsection::interactionV
std::vector< double > interactionV
Definition: mcsXsection_module.cc:301
protoana::mcsXsection::fdaughterOpeningAngle
double fdaughterOpeningAngle[NMAXDAUGTHERS]
Definition: mcsXsection_module.cc:384
simb::MCParticle::TrackId
int TrackId() const
Definition: MCParticle.h:210
protoana::mcsXsection::fdaughter_truth_Pt
double fdaughter_truth_Pt[NMAXDAUGTHERS]
Definition: mcsXsection_module.cc:403
protoana::mcsXsection::timeintersection
std::vector< double > timeintersection
Definition: mcsXsection_module.cc:314
protoana::mcsXsection::fdaughter_truth_StartPosition
double fdaughter_truth_StartPosition[NMAXDAUGTHERS][4]
Definition: mcsXsection_module.cc:398
protoana::mcsXsection::secondaryendx
std::vector< std::vector< double > > secondaryendx
Definition: mcsXsection_module.cc:351
protoana::mcsXsection::fdaughterTheta
double fdaughterTheta[NMAXDAUGTHERS]
Definition: mcsXsection_module.cc:375
protoana::mcsXsection::primtrk_range
std::vector< double > primtrk_range
Definition: mcsXsection_module.cc:322
protoana::mcsXsection::interactionY
std::vector< double > interactionY
Definition: mcsXsection_module.cc:298
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::mcsXsection::beamPosy_spec
std::vector< double > beamPosy_spec
Definition: mcsXsection_module.cc:225
protoana::mcsXsection::secondaryendz
std::vector< std::vector< double > > secondaryendz
Definition: mcsXsection_module.cc:355
protoana::mcsXsection::hitz_1
std::vector< float > hitz_1
Definition: mcsXsection_module.cc:238
protoana::mcsXsection::beamtrk_z_tpc
std::vector< int > beamtrk_z_tpc
Definition: mcsXsection_module.cc:209
art::DataViewImpl::isRealData
bool isRealData() const
Definition: DataViewImpl.cc:113
protoana::mcsXsection::hitz_0
std::vector< float > hitz_0
Definition: mcsXsection_module.cc:240
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
simb::MCParticle::Pt
double Pt(const int i=0) const
Definition: MCParticle.h:236
protoana::mcsXsection::fTrackModuleLabel
const art::InputTag fTrackModuleLabel
Definition: mcsXsection_module.cc:155
zpos
TH3F * zpos
Definition: doAna.cpp:31
recob::Track::Phi
double Phi() const
Definition: Track.h:178
e
const double e
Definition: gUpMuFluxGen.cxx:165
protoana::mcsXsection::interactionZ
std::vector< double > interactionZ
Definition: mcsXsection_module.cc:299
recob::Track::Length
double Length(size_t p=0) const
Access to various track properties.
Definition: Track.h:167
protoana::mcsXsection::michelscore
std::vector< std::vector< double > > michelscore
Definition: mcsXsection_module.cc:365
art::DataViewImpl::time
Timestamp time() const
Definition: DataViewImpl.cc:106
recob::Shower::dEdx
const std::vector< double > & dEdx() const
Definition: Shower.h:203
DEFINE_ART_MODULE
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:67
protoana::mcsXsection::fdaughterStartPosition
double fdaughterStartPosition[NMAXDAUGTHERS][3]
Definition: mcsXsection_module.cc:381
geo::fhicl::WireID
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
Definition: geo_types_fhicl.h:407
protoana::mcsXsection::fBeamModuleLabel
const art::InputTag fBeamModuleLabel
Definition: mcsXsection_module.cc:156
recob::Shower::MIPEnergy
const std::vector< double > & MIPEnergy() const
Definition: Shower.h:198
simb::MCParticle::EndProcess
std::string EndProcess() const
Definition: MCParticle.h:216
protoana::mcsXsection::fdaughterLength
double fdaughterLength[NMAXDAUGTHERS]
Definition: mcsXsection_module.cc:377
protoana::mcsXsection::fdaughterShowerBestPlane
int fdaughterShowerBestPlane[NMAXDAUGTHERS]
Definition: mcsXsection_module.cc:388
protoana::mcsXsection::ParMaker
G4ReweightParameterMaker ParMaker
Definition: mcsXsection_module.cc:432
recob::Track::Start
Point_t const & Start() const
Access to track position at different points.
Definition: Track.h:123
wirecell.dfp.graph.key
def key(type, name=None)
Definition: graph.py:13
protoana::mcsXsection::primtrk_truth_Z_tpc
std::vector< std::vector< int > > primtrk_truth_Z_tpc
Definition: mcsXsection_module.cc:332
simb::MCParticle::EndPz
double EndPz() const
Definition: MCParticle.h:243
protoana::mcsXsection::primtrk_hity
std::vector< std::vector< double > > primtrk_hity
Definition: mcsXsection_module.cc:324
Exceptions.h
Collection of exceptions for Geometry system.
beam::ProtoDUNEBeamEvent::GetRecoBeamMomenta
const std::vector< double > & GetRecoBeamMomenta() const
Definition: ProtoDUNEBeamEvent.h:91
protoana::mcsXsection::endJob
virtual void endJob() override
Definition: mcsXsection_module.cc:2039
protoana::mcsXsection::beamtrk_z
std::vector< float > beamtrk_z
Definition: mcsXsection_module.cc:207
protoana::mcsXsection::secondarystarty
std::vector< std::vector< double > > secondarystarty
Definition: mcsXsection_module.cc:352
recob::Track::Theta
double Theta() const
Access to spherical or geographical angles at vertex or at any point.
Definition: Track.h:176
protoana::mcsXsection::fdaughter_truth_Momentum
double fdaughter_truth_Momentum[NMAXDAUGTHERS][4]
Definition: mcsXsection_module.cc:401
simb::MCParticle::P
double P(const int i=0) const
Definition: MCParticle.h:234
art::Ptr::key
key_type key() const noexcept
Definition: Ptr.h:216
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
recob::Shower::OpenAngle
double OpenAngle() const
Definition: Shower.h:202
recob::Track::Vertex
Point_t const & Vertex() const
Definition: Track.h:124
protoana::mcsXsection::primtrk_dedx
std::vector< std::vector< double > > primtrk_dedx
Definition: mcsXsection_module.cc:321
protoana::mcsXsection::fdaughterShowerMIPEnergy
double fdaughterShowerMIPEnergy[NMAXDAUGTHERS]
Definition: mcsXsection_module.cc:386
protoana::mcsXsection::dq_2
std::vector< std::vector< float > > dq_2
Definition: mcsXsection_module.cc:345
simb::MCParticle::T
double T(const int i=0) const
Definition: MCParticle.h:224
art::DataViewImpl::getValidHandle
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
Definition: DataViewImpl.h:441
test.p
p
Definition: test.py:223
protoana::mcsXsection::primsectheta
std::vector< std::vector< double > > primsectheta
Definition: mcsXsection_module.cc:359
protoana::mcsXsection::Zintersection
std::vector< double > Zintersection
Definition: mcsXsection_module.cc:312
protoana::mcsXsection::fGeometry
geo::GeometryCore const * fGeometry
Definition: mcsXsection_module.cc:184
protoana::mcsXsection::interactionAnglesVT
std::vector< double > interactionAnglesVT
Definition: mcsXsection_module.cc:310
protoana::mcsXsection::int_2
std::vector< float > int_2
Definition: mcsXsection_module.cc:237
art::DataViewImpl::subRun
SubRunNumber_t subRun() const
Definition: DataViewImpl.cc:78
protoana::mcsXsection::fBeamlineUtils
protoana::ProtoDUNEBeamlineUtils fBeamlineUtils
Definition: mcsXsection_module.cc:159
protoana::ProtoDUNETruthUtils::GetMCParticleFromRecoTrack
const simb::MCParticle * GetMCParticleFromRecoTrack(detinfo::DetectorClocksData const &clockData, const recob::Track &track, art::Event const &evt, std::string trackModule) const
Definition: ProtoDUNETruthUtils.cxx:609
protoana::ProtoDUNEPFParticleUtils::GetPFParticlesFromBeamSlice
const std::vector< const recob::PFParticle * > GetPFParticlesFromBeamSlice(art::Event const &evt, const std::string particleLabel) const
Return the pointers for the PFParticles in the beam slice. Returns an empty vector is no beam slice w...
Definition: ProtoDUNEPFParticleUtils.cxx:175
recob::Shower::Direction
const TVector3 & Direction() const
Definition: Shower.h:189
protoana::mcsXsection::primtrk_resrange
std::vector< std::vector< double > > primtrk_resrange
Definition: mcsXsection_module.cc:320
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::mcsXsection::beamtrk_z_wire
std::vector< int > beamtrk_z_wire
Definition: mcsXsection_module.cc:208
art::DataViewImpl::run
RunNumber_t run() const
Definition: DataViewImpl.cc:71
max
static int max(int a, int b)
Definition: fortranscanner.cpp:60366
simb::MCParticle::EndT
double EndT() const
Definition: MCParticle.h:229
protoana::mcsXsection::interactionU
std::vector< double > interactionU
Definition: mcsXsection_module.cc:300
geo::TPCID
The data type to uniquely identify a TPC.
Definition: geo_types.h:386
protoana::mcsXsection::beamtrk_y
std::vector< float > beamtrk_y
Definition: mcsXsection_module.cc:206
geo::GeometryCore
Description of geometry of one entire detector.
Definition: GeometryCore.h:1468
protoana::mcsXsection::fdaughter_truth_Mass
double fdaughter_truth_Mass[NMAXDAUGTHERS]
Definition: mcsXsection_module.cc:404
simb::MCTrajectory::TrajectoryProcesses
ProcessMap const & TrajectoryProcesses() const
Definition: MCTrajectory.h:188
geo::GeometryCore::View
View_t View(geo::PlaneID const &pid) const
Returns the view (wire orientation) on the channels of specified TPC plane.
Definition: GeometryCore.cxx:601
geo_types.h
Definition of data types for geometry description.
protoana::mcsXsection::beamDiry_spec
std::vector< double > beamDiry_spec
Definition: mcsXsection_module.cc:229
protoana::mcsXsection::int_0
std::vector< float > int_0
Definition: mcsXsection_module.cc:241
soln
float soln(float x1, float x2, float x3, float x4, float y1, float y2, float y3, float y4, float result[2])
Definition: mcsXsection_module.cc:99
cheat::ParticleInventoryService::ParticleList
const sim::ParticleList & ParticleList() const
Definition: ParticleInventoryService_service.cc:114
protoana::mcsXsection::fdaughter_truth_EndPosition
double fdaughter_truth_EndPosition[NMAXDAUGTHERS][4]
Definition: mcsXsection_module.cc:399
protoana::mcsXsection::fdaughterKineticEnergy
double fdaughterKineticEnergy[NMAXDAUGTHERS][3]
Definition: mcsXsection_module.cc:391
protoana::mcsXsection::g4rw_elast
std::vector< double > g4rw_elast
Definition: mcsXsection_module.cc:426
keras_to_tensorflow.int
int
Definition: keras_to_tensorflow.py:69
theta12
double theta12(double x1, double x2, double y1, double y2, double z1, double z2, double x1p, double x2p, double y1p, double y2p, double z1p, double z2p)
Definition: mcsXsection_module.cc:121
protoana::mcsXsection::primtrk_truth_Z
std::vector< std::vector< double > > primtrk_truth_Z
Definition: mcsXsection_module.cc:330
cache_state.m
m
Definition: cache_state.py:415
recob::TrackTrajectory::EndDirection
Vector_t EndDirection() const
Returns the direction of the trajectory at the last point.
Definition: TrackTrajectory.h:321
protoana::mcsXsection::fprimary_truth_EndProcess
std::string fprimary_truth_EndProcess
Definition: mcsXsection_module.cc:281
protoana::mcsXsection::primtrk_hitz_tpc
std::vector< std::vector< int > > primtrk_hitz_tpc
Definition: mcsXsection_module.cc:327
tracks
Definition: tracks.py:1
simb::MCParticle::Vx
double Vx(const int i=0) const
Definition: MCParticle.h:221
protoana::mcsXsection::fdaughterPhi
double fdaughterPhi[NMAXDAUGTHERS]
Definition: mcsXsection_module.cc:376
recob::Track::ID
int ID() const
Definition: Track.h:198
protoana::mcsXsection::hitz_2
std::vector< float > hitz_2
Definition: mcsXsection_module.cc:216
Hit.h
Declaration of signal hit object.
protoana::mcsXsection::beam_cuts
protoana::ProtoDUNEBeamCuts beam_cuts
Definition: mcsXsection_module.cc:158
recob::Track::StartMomentum
double StartMomentum() const
Definition: Track.h:143
recob::PFParticle
Hierarchical representation of particle flow.
Definition: PFParticle.h:44
protoana::mcsXsection::deltatimeint
std::vector< double > deltatimeint
Definition: mcsXsection_module.cc:317
protoana::mcsXsection::fdaughter_truth_P
double fdaughter_truth_P[NMAXDAUGTHERS]
Definition: mcsXsection_module.cc:400
cet::sqlite::min
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:55
simb::MCTrajectory::size
size_type size() const
Definition: MCTrajectory.h:166
simb::MCParticle::EndPy
double EndPy() const
Definition: MCParticle.h:242
protoana::mcsXsection::fdaughterNHits
int fdaughterNHits[NMAXDAUGTHERS]
Definition: mcsXsection_module.cc:374
End
void End(void)
Definition: gXSecComp.cxx:210
protoana::mcsXsection::fdaughterShowerdEdx
double fdaughterShowerdEdx[NMAXDAUGTHERS]
Definition: mcsXsection_module.cc:387
reverse
static unsigned int reverse(QString &chars, unsigned char *level, unsigned int a, unsigned int b)
Definition: qstring.cpp:11649
protoana::mcsXsection::interactionProcesslist
std::vector< std::string > interactionProcesslist
Definition: mcsXsection_module.cc:307
protoana::mcsXsection::timeintersection1
std::vector< double > timeintersection1
Definition: mcsXsection_module.cc:315
protoana::mcsXsection::fdaughter_truth_Pdg
int fdaughter_truth_Pdg[NMAXDAUGTHERS]
Definition: mcsXsection_module.cc:397
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::mcsXsection::g4rw_set_weights
std::vector< double > g4rw_set_weights
Definition: mcsXsection_module.cc:427
protoana::mcsXsection::beginJob
virtual void beginJob() override
Definition: mcsXsection_module.cc:475
simb::MCParticle::Momentum
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:220
simb::MCParticle::Pz
double Pz(const int i=0) const
Definition: MCParticle.h:232
Track.h
Provides recob::Track data product.
protoana::mcsXsection::primtrk_hitz
std::vector< std::vector< double > > primtrk_hitz
Definition: mcsXsection_module.cc:325
simb::MCParticle::Vz
double Vz(const int i=0) const
Definition: MCParticle.h:223
protoana::mcsXsection::beamtrk_Eng
std::vector< float > beamtrk_Eng
Definition: mcsXsection_module.cc:214
protoana::mcsXsection::beamtrk_Py
std::vector< float > beamtrk_Py
Definition: mcsXsection_module.cc:212
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::mcsXsection::FillCosmicsTree
void FillCosmicsTree(art::Event const &evt, std::string pfParticleTag)
Definition: mcsXsection_module.cc:2043
protoana::mcsXsection::fdaughterMomentumByRangeMuon
double fdaughterMomentumByRangeMuon[NMAXDAUGTHERS]
Definition: mcsXsection_module.cc:390
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
recob::Track::End
Point_t const & End() const
Definition: Track.h:125
protoana::mcsXsection::dq_1
std::vector< std::vector< float > > dq_1
Definition: mcsXsection_module.cc:344
protoana::mcsXsection::primtrk_pitch
std::vector< std::vector< double > > primtrk_pitch
Definition: mcsXsection_module.cc:329
NMAXDAUGTHERS
const int NMAXDAUGTHERS
Definition: mcsXsection_module.cc:93
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
reco_momentum_tuples.count
int count
Definition: reco_momentum_tuples.py:202
beam::ProtoDUNEBeamEvent::GetBeamTracks
const std::vector< recob::Track > & GetBeamTracks() const
Definition: ProtoDUNEBeamEvent.h:192
protoana::mcsXsection::RWManager
G4ReweightManager RWManager
Definition: mcsXsection_module.cc:433
CRT::PlaneID
detail::Node< FrameID, bool > PlaneID
Definition: CRTID.h:125
protoana::mcsXsection::interactionAnglesZT
std::vector< double > interactionAnglesZT
Definition: mcsXsection_module.cc:311
protoana::mcsXsection::interactionT
std::vector< double > interactionT
Definition: mcsXsection_module.cc:303
protoana::mcsXsection::beamPosz_spec
std::vector< double > beamPosz_spec
Definition: mcsXsection_module.cc:226
simb::MCParticle::EndPx
double EndPx() const
Definition: MCParticle.h:241
protoana::mcsXsection::interactionW
std::vector< double > interactionW
Definition: mcsXsection_module.cc:302
protoana::mcsXsection::MultiRW
G4MultiReweighter MultiRW
Definition: mcsXsection_module.cc:434
protoana::mcsXsection::dQtrackend
std::vector< std::vector< double > > dQtrackend
Definition: mcsXsection_module.cc:357
protoana::mcsXsection::fdaughterStartDirection
double fdaughterStartDirection[NMAXDAUGTHERS][3]
Definition: mcsXsection_module.cc:383
protoana::mcsXsection::fdaughter_truth_Theta
double fdaughter_truth_Theta[NMAXDAUGTHERS]
Definition: mcsXsection_module.cc:405
tca::evt
TCEvent evt
Definition: DataStructs.cxx:7
protoana::mcsXsection::interactionX
std::vector< double > interactionX
Definition: mcsXsection_module.cc:297
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
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::mcsXsection::truthUtil
protoana::ProtoDUNETruthUtils truthUtil
Definition: mcsXsection_module.cc:164
protoana::G4ReweightUtils::CreateNRWTrajs
std::vector< G4ReweightTraj * > CreateNRWTrajs(const simb::MCParticle &part, const sim::ParticleList &plist, art::ServiceHandle< geo::Geometry > geo_serv, int event, std::string material_name="LAr", bool fVerbose=false)
Definition: G4ReweightUtils.cxx:115
recob::TrackTrajectory::StartDirection
Vector_t StartDirection() const
Returns the direction of the trajectory at the first point.
Definition: TrackTrajectory.h:317
protoana::mcsXsection::secondarystartx
std::vector< std::vector< double > > secondarystartx
Definition: mcsXsection_module.cc:350
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
trkf::Plane
recob::tracking::Plane Plane
Definition: TrackState.h:17
protoana::mcsXsection::fdaughterEndDirection
double fdaughterEndDirection[NMAXDAUGTHERS][3]
Definition: mcsXsection_module.cc:382
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::mcsXsection::secondaryendy
std::vector< std::vector< double > > secondaryendy
Definition: mcsXsection_module.cc:353
protoana::mcsXsection::int_1
std::vector< float > int_1
Definition: mcsXsection_module.cc:239
protoana::mcsXsection::nonescore
std::vector< std::vector< double > > nonescore
Definition: mcsXsection_module.cc:366
protoana::mcsXsection::mcsXsection
mcsXsection(fhicl::ParameterSet const &p)
Definition: mcsXsection_module.cc:444
anab::MVAReader::getOutput
std::array< float, N > getOutput(size_t key) const
Get copy of the MVA output vector at index "key".
Definition: MVAReader.h:129
recob::Shower::ID
int ID() const
Definition: Shower.h:187
protoana::mcsXsection::pfpUtil
protoana::ProtoDUNEPFParticleUtils pfpUtil
Definition: mcsXsection_module.cc:162
protoana::mcsXsection::analyze
void analyze(art::Event const &evt) override
Definition: mcsXsection_module.cc:701
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
protoana::mcsXsection::interactionAnglesUT
std::vector< double > interactionAnglesUT
Definition: mcsXsection_module.cc:309
art::Event::id
EventID id() const
Definition: Event.cc:34
simb::MCParticle::Vy
double Vy(const int i=0) const
Definition: MCParticle.h:222
protoana::G4ReweightUtils::GetNTrajPMSigmaWeights
std::pair< double, double > GetNTrajPMSigmaWeights(const std::vector< G4ReweightTraj * > &trajs, G4MultiReweighter &rw, size_t iPar)
Definition: G4ReweightUtils.cxx:296
recob::Track
Track from a non-cascading particle.A recob::Track consists of a recob::TrackTrajectory, plus additional members relevant for a "fitted" track:
Definition: Track.h:49
calo
calorimetry
Definition: CaloChecker_module.cc:43
trkf::TrackMomentumCalculator::GetTrackMomentum
double GetTrackMomentum(double trkrange, int pdg) const
Definition: TrackMomentumCalculator.cxx:144
fhicl::exception
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
protoana::mcsXsection::fdaughterShowerEnergy
double fdaughterShowerEnergy[NMAXDAUGTHERS]
Definition: mcsXsection_module.cc:385
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
protoana::mcsXsection::secondarystartz
std::vector< std::vector< double > > secondarystartz
Definition: mcsXsection_module.cc:354
protoana::mcsXsection::fdaughter_truth_EndMomentum
double fdaughter_truth_EndMomentum[NMAXDAUGTHERS][4]
Definition: mcsXsection_module.cc:402
protoana::mcsXsection::fdaughter_truth_Phi
double fdaughter_truth_Phi[NMAXDAUGTHERS]
Definition: mcsXsection_module.cc:406
reco_momentum_tuples.x1
x1
Definition: reco_momentum_tuples.py:82
protoana::mcsXsection::dq_0
std::vector< std::vector< float > > dq_0
Definition: mcsXsection_module.cc:343
protoana::mcsXsection::trackUtil
protoana::ProtoDUNETrackUtils trackUtil
Definition: mcsXsection_module.cc:163
Generated by   doxygen 1.8.11