NueAna_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: NueAna
3 // Module Type: analyzer
4 // File: NueAna_module.cc
5 //
6 // dorota.stefan@cern.ch, robert.sulej@cern.ch, tjyang@fnal.gov
7 //
8 ////////////////////////////////////////////////////////////////////////
9 
16 #include "art_root_io/TFileService.h"
18 #include "fhiclcpp/ParameterSet.h"
20 
21 // LArSoft includes
31 #include "lardata/ArtDataHelper/TrackUtils.h" // lar::util::TrackPitchInView()
43 
44 // ROOT includes
45 #include "TTree.h"
46 #include "TTimeStamp.h"
47 
48 // C/C++ libraries
49 #include <memory>
50 #include <utility>
51 
52 constexpr int kMaxTrack = 1000; //maximum number of tracks
53 constexpr int kMaxShower = 1000; //maximum number of showers
54 constexpr int kMaxHits = 40000; //maximum number of hits;
55 constexpr int kMaxVertices = 100; //max number of 3D vertices
56 constexpr int kMaxPrimaries = 20000; //maximum number of primary particles
57 constexpr int kMaxFlash = 1000; //maximum number of flashes
58 //constexpr int kMaxTrackHits = 1000; //maximum number of track trajectory points
59 namespace dunefd {
60  class Hit2D;
61  class NueAna;
62  class IniSegAlg;
63 }
64 
66 public:
67  explicit NueAna(fhicl::ParameterSet const & p);
68  // The destructor generated by the compiler is fine for classes
69  // without bare pointers or other resource use.
70 
71  // Plugins should not be copied or assigned.
72  NueAna(NueAna const &) = delete;
73  NueAna(NueAna &&) = delete;
74  NueAna & operator = (NueAna const &) = delete;
75  NueAna & operator = (NueAna &&) = delete;
76 
77  // Required functions.
78  void analyze(art::Event const & e) override;
79 
80  // Selected optional functions.
81  void beginJob() override;
82  void endJob() override;
83 
84  //void beginSubRun(const art::SubRun& sr);
85  void endSubRun(const art::SubRun& sr) override;
86 
87  void reconfigure(fhicl::ParameterSet const& p) ;
88 
89 private:
90 
91  void ResetVars();
92  bool insideFidVol(const double posX, const double posY, const double posZ);
93 
94  // Declare member data here.
95  TTree *fTree;
96  TTree* fPOT;
97 
98  // TTree variables
99 
100  // Run information
101  int run;
102  int subrun;
103  int event;
104  float evttime;
105  float taulife;
106  short isdata;
107  double pot;
108 
109  // Track information
110  int ntracks_reco; //number of reconstructed tracks
111  int trkid[kMaxTrack]; //track id from recob::Track::ID(), this does not have to be the index of track
112  float trkstartx[kMaxTrack]; //track start position (cm)
115  float trkendx[kMaxTrack]; //track end position (cm)
118  float trkstartdcosx[kMaxTrack]; //track start direction cosine
121  float trkenddcosx[kMaxTrack]; //track end direction cosine
124  float trklen[kMaxTrack]; //track length (cm)
125  float trkke[kMaxTrack][3]; //track kinetic energy (in 3 planes)
126  float trkpida[kMaxTrack][3]; //track PIDA (in 3 planes)
127  int trkbestplane[kMaxTrack]; //best plane for trkke and trkpida
128  float trkmomrange[kMaxTrack]; //Track momentum using range
129 
130  //geant information for the track
131  int trkg4id[kMaxTrack]; //geant track id for the track
132  int trkg4pdg[kMaxTrack]; //pdg of geant particle
133  float trkg4startx[kMaxTrack]; //start position of geant particle
136  float trkg4initdedx[kMaxTrack]; //initial dE/dx of the track using true energy (MeV/cm)
137 
138  int nhits;
140  Short_t hit_plane[kMaxHits]; //plane number
141  Short_t hit_wire[kMaxHits]; //wire number
142  Int_t hit_channel[kMaxHits]; //channel ID
143  Short_t hit_tpc[kMaxHits]; //tpc
144  Float_t hit_peakT[kMaxHits]; //peak time
145  Float_t hit_charge[kMaxHits]; //charge (area)
146  Float_t hit_summedADC[kMaxHits]; //summed ADC
147  Float_t hit_startT[kMaxHits]; //hit start time
148  Float_t hit_endT[kMaxHits]; //hit end time
149  Int_t hit_trkkey[kMaxHits]; //track index if hit is associated with a track
150  Float_t hit_dQds[kMaxHits]; //hit dQ/ds
151  Float_t hit_dEds[kMaxHits]; //hit dE/ds
152  Float_t hit_resrange[kMaxHits]; //hit residual range
153  Int_t hit_shwkey[kMaxHits]; //shower index if hit is associated with a shower
154 
155  // vertex information
156  int infidvol;
157  Short_t nvtx; //number of vertices
158  Float_t vtx[kMaxVertices][3]; //vtx[3]
159 
160  Float_t vtxrecomc; // distance between mc and reco vtx
161  Float_t vtxrecomcx;
162  Float_t vtxrecomcy;
163  Float_t vtxrecomcz;
164 
165  // shower information
166  int nshws; //number of showers
167  int shwid[kMaxShower]; //recob::Shower::ID()
168  Float_t shwdcosx[kMaxShower]; //shower direction cosine
171  Float_t shwstartx[kMaxShower]; //shower start position (cm)
174  Float_t shwenergy[kMaxShower][3]; //shower energy measured on the 3 planes (GeV)
175  Float_t shwdedx[kMaxShower][3]; //shower dE/dx of the initial track measured on the 3 plane (MeV/cm)
176  int shwbestplane[kMaxShower]; //recommended plane for energy and dE/dx information
177  int shwg4id[kMaxShower]; //geant track id for the shower
178 
179  // flash information
180  int flash_total; //total number of flashes
181  Float_t flash_time[kMaxFlash]; //flash time
182  Float_t flash_width[kMaxFlash]; //flash width
183  Float_t flash_abstime[kMaxFlash]; //flash absolute time
184  Float_t flash_YCenter[kMaxFlash]; //flash y center (cm)
185  Float_t flash_YWidth[kMaxFlash]; //flash y width (cm)
186  Float_t flash_ZCenter[kMaxFlash]; //flash z center (cm)
187  Float_t flash_ZWidth[kMaxFlash]; //flash z width (cm)
188  Float_t flash_TotalPE[kMaxFlash]; //flash total pe
189 
190  //mctruth information
191  Int_t mcevts_truth; //number of neutrino Int_teractions in the spill
192  Int_t nuPDG_truth; //neutrino PDG code
193  Int_t ccnc_truth; //0=CC 1=NC
194  Int_t mode_truth; //0=QE/El, 1=RES, 2=DIS, 3=Coherent production
195  Float_t enu_truth; //true neutrino energy
196  Float_t Q2_truth; //Momentum transfer squared
197  Float_t W_truth; //hadronic invariant mass
198  Float_t X_truth;
199  Float_t Y_truth;
200  Int_t hitnuc_truth; //hit nucleon
201  Int_t target_truth; //hit nucleus
202  Float_t nuvtxx_truth; //neutrino vertex x
203  Float_t nuvtxy_truth; //neutrino vertex y
204  Float_t nuvtxz_truth; //neutrino vertex z
205  Float_t nu_dcosx_truth; //neutrino dcos x
206  Float_t nu_dcosy_truth; //neutrino dcos y
207  Float_t nu_dcosz_truth; //neutrino dcos z
208  Float_t lep_mom_truth; //lepton momentum
209  Float_t lep_dcosx_truth; //lepton dcos x
210  Float_t lep_dcosy_truth; //lepton dcos y
211  Float_t lep_dcosz_truth; //lepton dcos z
212  Float_t t0_truth; // t0
213 
214  // === Storing Geant4 MC Truth Information ===
215  int no_primaries; //<---Number of primary Geant4 particles in the event
216  int geant_list_size; //<---Number of Geant4 particles tracked
217  int pdg[kMaxPrimaries]; //<---PDG Code number of this particle
218  Float_t Eng[kMaxPrimaries]; //<---Energy of the particle
219  Float_t Px[kMaxPrimaries]; //<---Px momentum of the particle
220  Float_t Py[kMaxPrimaries]; //<---Py momentum of the particle
221  Float_t Pz[kMaxPrimaries]; //<---Pz momentum of the particle
222  Float_t StartPointx[kMaxPrimaries]; //<---X position that this Geant4 particle started at
223  Float_t StartPointy[kMaxPrimaries]; //<---Y position that this Geant4 particle started at
224  Float_t StartPointz[kMaxPrimaries]; //<---Z position that this Geant4 particle started at
225  Float_t EndPointx[kMaxPrimaries]; //<---X position that this Geant4 particle ended at
226  Float_t EndPointy[kMaxPrimaries]; //<---Y position that this Geant4 particle ended at
227  Float_t EndPointz[kMaxPrimaries]; //<---Z position that this Geant4 particle ended at
228  Float_t Startdcosx[kMaxPrimaries]; //<---X direction cosine that Geant4 particle started at
229  Float_t Startdcosy[kMaxPrimaries]; //<---Y direction cosine that Geant4 particle started at
230  Float_t Startdcosz[kMaxPrimaries]; //<---Z direction cosine that Geant4 particle started at
231  int NumberDaughters[kMaxPrimaries]; //<---Number of Daughters this particle has
232  int TrackId[kMaxPrimaries]; //<---Geant4 TrackID number
233  int Mother[kMaxPrimaries]; //<---TrackID of the mother of this particle
234  int process_primary[kMaxPrimaries]; //<---Is this particle primary (primary = 1, non-primary = 1)
235  std::vector<std::string> G4Process; //<---The process which created this particle
236  std::vector<std::string> G4FinalProcess; //<---The last process which this particle went under
237 
238  //flux information
239  Int_t ptype_flux; //Parent GEANT code particle ID
240  Float_t pdpx_flux; //Parent X momentum at decay point (GeV)
241  Float_t pdpy_flux; //Parent Y momentum at decay point (GeV)
242  Float_t pdpz_flux; //Parent Z momentum at decay point (GeV)
243  Int_t pntype_flux; //oscillated neutrino type
244  Float_t vx_flux; //X position of hadron/muon decay (cm)
245  Float_t vy_flux; //Y position of hadron/muon decay (cm)
246  Float_t vz_flux; //Z position of hadron/muon decay (cm)
247 
248  //end of ttree variables
249 
250  //Module labels to get data products
260 
261  double fFidVolCut;
262 
264 
265 };
266 
267 
269  : EDAnalyzer(pset)
270  , fCalorimetryAlg(pset.get<fhicl::ParameterSet>("CalorimetryAlg"))
271 {
272  reconfigure(pset);
273 }
274 
275 //void dunefd::NueAna::beginSubRun(const art::SubRun& sr){
276 //
277 // auto potListHandle = sr.getHandle< sumdata::POTSummary >(fPOTModuleLabel);
278 // if (potListHandle)
279 // pot = potListHandle->totpot;
280 // else
281 // pot = 0.;
282 // if (fPOT) fPOT->Fill();
283 //}
284 
286 
287  auto potListHandle = sr.getHandle< sumdata::POTSummary >(fPOTModuleLabel);
288  if (potListHandle)
289  pot = potListHandle->totpot;
290  else
291  pot = 0.;
292  if (fPOT) fPOT->Fill();
293 }
294 
296 {
297 
298  std::cout << " **************************************** analyze ************ " << std::endl;
299  // Implementation of required member function here.
300  ResetVars();
304  const sim::ParticleList& plist = pi_serv->ParticleList();
305 
306  run = evt.run();
307  subrun = evt.subRun();
308  event = evt.id().event();
309  art::Timestamp ts = evt.time();
310  TTimeStamp tts(ts.timeHigh(), ts.timeLow());
311  evttime = tts.AsDouble();
312  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
313  auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(evt, clockData);
314  taulife = detProp.ElectronLifetime();
315  isdata = evt.isRealData();
316 
317  // * hits
318  std::vector<art::Ptr<recob::Hit> > hitlist;
319  auto hitListHandle = evt.getHandle< std::vector<recob::Hit> >(fHitsModuleLabel);
320  if (hitListHandle)
321  art::fill_ptr_vector(hitlist, hitListHandle);
322 
323  // * tracks
324  std::vector<art::Ptr<recob::Track> > tracklist;
325  auto trackListHandle = evt.getHandle< std::vector<recob::Track> >(fTrackModuleLabel);
326  if (trackListHandle)
327  art::fill_ptr_vector(tracklist, trackListHandle);
328 
329  // * vertices
330  std::vector<art::Ptr<recob::Vertex> > vtxlist;
331  auto vtxListHandle = evt.getHandle< std::vector<recob::Vertex> >(fVertexModuleLabel);
332  if (vtxListHandle)
333  art::fill_ptr_vector(vtxlist, vtxListHandle);
334 
335  // * showers
336  std::vector<art::Ptr<recob::Shower>> shwlist;
337  auto shwListHandle = evt.getHandle<std::vector<recob::Shower> >(fShowerModuleLabel);
338  if (shwListHandle)
339  art::fill_ptr_vector(shwlist, shwListHandle);
340 
341  // * flashes
342  std::vector<art::Ptr<recob::OpFlash> > flashlist;
343  auto flashListHandle = evt.getHandle< std::vector<recob::OpFlash> >(fFlashModuleLabel);
344  if (flashListHandle)
345  art::fill_ptr_vector(flashlist, flashListHandle);
346 
347  // * associations
348  art::FindManyP<recob::Hit> fmth(trackListHandle, evt, fTrackModuleLabel);
349  art::FindManyP<recob::Hit, recob::TrackHitMeta> fmthm(trackListHandle, evt, fTrackModuleLabel);
350  art::FindManyP<recob::SpacePoint> fmhs(hitListHandle, evt, fTrackModuleLabel);
351  art::FindMany<anab::Calorimetry> fmcal(trackListHandle, evt, fCalorimetryModuleLabel);
352 
353  //hit information
354  nhits = hitlist.size();
356  for (int i = 0; i < nhits && i < kMaxHits ; ++i){//loop over hits
357  hit_channel[i] = hitlist[i]->Channel();
358  hit_plane[i] = hitlist[i]->WireID().Plane;
359  hit_wire[i] = hitlist[i]->WireID().Wire;
360  hit_tpc[i] = hitlist[i]->WireID().TPC;
361  hit_peakT[i] = hitlist[i]->PeakTime();
362  hit_charge[i] = hitlist[i]->Integral();
363  hit_summedADC[i] = hitlist[i]->SummedADC();
364  hit_startT[i] = hitlist[i]->PeakTimeMinusRMS();
365  hit_endT[i] = hitlist[i]->PeakTimePlusRMS();
366 
367  }
368 
369  //track information
370  ntracks_reco=tracklist.size();
371 
372  recob::Track::Vector_t larStart;
373  recob::Track::Vector_t larEnd;
375  for(int i=0; i<std::min(int(tracklist.size()),kMaxTrack);++i){
376  recob::Track::Point_t trackStart, trackEnd;
377  std::tie(trackStart, trackEnd) = tracklist[i]->Extent();
378  larStart = tracklist[i]->VertexDirection();
379  larEnd = tracklist[i]->EndDirection();
380 
381  trkid[i] = tracklist[i]->ID();
382  trkstartx[i] = trackStart.X();
383  trkstarty[i] = trackStart.Y();
384  trkstartz[i] = trackStart.Z();
385  trkendx[i] = trackEnd.X();
386  trkendy[i] = trackEnd.Y();
387  trkendz[i] = trackEnd.Z();
388  trkstartdcosx[i] = larStart.X();
389  trkstartdcosy[i] = larStart.Y();
390  trkstartdcosz[i] = larStart.Z();
391  trkenddcosx[i] = larEnd.X();
392  trkenddcosy[i] = larEnd.Y();
393  trkenddcosz[i] = larEnd.Z();
394  trklen[i] = tracklist[i]->Length();
395  if (!std::isnan(trklen[i])) trkmomrange[i] = trkm.GetTrackMomentum(trklen[i],13);
396  if (fmthm.isValid()){
397  auto vhit = fmthm.at(i);
398  auto vmeta = fmthm.data(i);
399  for (size_t h = 0; h < vhit.size(); ++h){
400  if (vhit[h].key()<kMaxHits){
401  hit_trkkey[vhit[h].key()] = tracklist[i].key();
402  if (vmeta[h]->Dx()){
403  hit_dQds[vhit[h].key()] = vhit[h]->Integral()*fCalorimetryAlg.LifetimeCorrection(clockData, detProp, vhit[h]->PeakTime())/vmeta[h]->Dx();
404  hit_dEds[vhit[h].key()] = fCalorimetryAlg.dEdx_AREA(clockData, detProp, *vhit[h], vmeta[h]->Dx());
405  }
406  hit_resrange[vhit[h].key()] = tracklist[i]->Length(vmeta[h]->Index());
407  }
408  }//loop over all hits
409  }//fmthm is valid
410  else if (fmth.isValid()){
411  std::vector< art::Ptr<recob::Hit> > vhit = fmth.at(i);
412  for (size_t h = 0; h < vhit.size(); ++h){
413  if (vhit[h].key()<kMaxHits){
414  hit_trkkey[vhit[h].key()] = tracklist[i].key();
415  }
416  }
417  }
418  if (fmcal.isValid()){
419  unsigned maxnumhits = 0;
420  std::vector<const anab::Calorimetry*> calos = fmcal.at(i);
421  for (auto const& calo : calos){
422  if (calo->PlaneID().isValid){
423  trkke[i][calo->PlaneID().Plane] = calo->KineticEnergy();
424  if (calo->dEdx().size()>maxnumhits){
425  maxnumhits = calo->dEdx().size();
426  trkbestplane[i] = calo->PlaneID().Plane;
427  }
428  double pida = 0;
429  int used_trkres = 0;
430  for (size_t ip = 0; ip<calo->dEdx().size(); ++ip){
431  if (calo->ResidualRange()[ip]<30){
432  pida += calo->dEdx()[ip]*pow(calo->ResidualRange()[ip],0.42);
433  ++used_trkres;
434  }
435  }
436  if (used_trkres) pida/=used_trkres;
437  trkpida[i][calo->PlaneID().Plane] = pida;
438  }
439  }
440  }
441  if (!isdata&&fmth.isValid()){
442  // Find true track for each reconstructed track
443  int TrackID = 0;
444  std::vector< art::Ptr<recob::Hit> > allHits = fmth.at(i);
445 
446  std::map<int,double> trkide;
447  for(size_t h = 0; h < allHits.size(); ++h){
448  art::Ptr<recob::Hit> hit = allHits[h];
449  std::vector<sim::TrackIDE> TrackIDs = bt_serv->HitToTrackIDEs(clockData, hit);
450  for(size_t e = 0; e < TrackIDs.size(); ++e){
451  trkide[TrackIDs[e].trackID] += TrackIDs[e].energy;
452  }
453  }
454  // Work out which IDE despoited the most charge in the hit if there was more than one.
455  double maxe = -1;
456  double tote = 0;
457  for (std::map<int,double>::iterator ii = trkide.begin(); ii!=trkide.end(); ++ii){
458  tote += ii->second;
459  if ((ii->second)>maxe){
460  maxe = ii->second;
461  TrackID = ii->first;
462  }
463  }
464  // Now have trackID, so get PdG code and T0 etc.
465  const simb::MCParticle *particle = pi_serv->TrackIdToParticle_P(TrackID);
466  if (particle){
467  trkg4id[i] = TrackID;
468  trkg4pdg[i] = particle->PdgCode();
469  trkg4startx[i] = particle->Vx();
470  trkg4starty[i] = particle->Vy();
471  trkg4startz[i] = particle->Vz();
472  float sum_energy = 0;
473  int numhits = 0;
474  //std::map<float,float> hite;
475  double x = 0;
476  double y = 0;
477  double z = 0;
478  double mindis = 1e10;
479  //find the closest point to the neutrino vertex
480  for(size_t h = 0; h < allHits.size(); ++h){
481  art::Ptr<recob::Hit> hit = allHits[h];
482  if (hit->WireID().Plane==2){
483  std::vector<art::Ptr<recob::SpacePoint> > spts = fmhs.at(hit.key());
484  if (spts.size()){
485  double dis = sqrt(pow(spts[0]->XYZ()[0]-trkg4startx[i],2)+
486  pow(spts[0]->XYZ()[1]-trkg4starty[i],2)+
487  pow(spts[0]->XYZ()[2]-trkg4startz[i],2));
488  if (dis<mindis){
489  mindis = dis;
490  x = spts[0]->XYZ()[0];
491  y = spts[0]->XYZ()[1];
492  z = spts[0]->XYZ()[2];
493  }
494  }
495  }
496  }
497  for(size_t h = 0; h < allHits.size(); ++h){
498  art::Ptr<recob::Hit> hit = allHits[h];
499  if (hit->WireID().Plane==2){
500  std::vector<art::Ptr<recob::SpacePoint> > spts = fmhs.at(hit.key());
501  if (spts.size()){
502  if (sqrt(pow(spts[0]->XYZ()[0]-x,2)+
503  pow(spts[0]->XYZ()[1]-y,2)+
504  pow(spts[0]->XYZ()[2]-z,2))<3){
505  std::vector<sim::TrackIDE> TrackIDs = bt_serv->HitToTrackIDEs(clockData, hit);
506  float toten = 0;
507  for(size_t e = 0; e < TrackIDs.size(); ++e){
508  //sum_energy += TrackIDs[e].energy;
509  toten+=TrackIDs[e].energy;
510  }
511  if (toten){
512  sum_energy += toten;
513  ++numhits;
514  }
515  }
516  }
517  }
518  }
519 
520  float pitch = 0;
521  float dis1 = sqrt(pow(trkstartx[i]-trkg4startx[i],2)+
522  pow(trkstarty[i]-trkg4starty[i],2)+
523  pow(trkstartz[i]-trkg4startz[i],2));
524  float dis2 = sqrt(pow(trkendx[i]-trkg4startx[i],2)+
525  pow(trkendy[i]-trkg4starty[i],2)+
526  pow(trkendz[i]-trkg4startz[i],2));
527  if (dis1<dis2){
528  try{
529  pitch = lar::util::TrackPitchInView(*(tracklist[i]),geo::kZ,0);
530  }
531  catch(...){
532  pitch = 0;
533  }
534  }
535  else{
536  try{
537  pitch = lar::util::TrackPitchInView(*(tracklist[i]), geo::kZ,tracklist[i]->NumberTrajectoryPoints()-1);
538  }
539  catch(...){
540  pitch = 0;
541  }
542  }
543  if ( pitch && numhits ) {
544  trkg4initdedx[i] = sum_energy/(numhits*pitch);
545  }
546  else{
547  trkg4initdedx[i] = 0;
548  }
549  }//if (particle)
550  }//MC
551  }
552 
553  //vertex information
554  nvtx = vtxlist.size();
555  for (int i = 0; i < nvtx && i < kMaxVertices ; ++i){//loop over hits
556  Double_t xyz[3] = {};
557  vtxlist[i]->XYZ(xyz);
558  for (size_t j = 0; j<3; ++j) vtx[i][j] = xyz[j];
559  }
560 
561  //shower information
562  if (shwListHandle.isValid()){
563  art::FindManyP<recob::Hit> fmsh(shwListHandle, evt, fShowerModuleLabel);
564 
565  nshws = shwlist.size();
566 
567  for (int i = 0; i<std::min(int(shwlist.size()),kMaxShower); ++i){
568  shwid[i] = shwlist[i]->ID();
569  shwdcosx[i] = shwlist[i]->Direction().X();
570  shwdcosy[i] = shwlist[i]->Direction().Y();
571  shwdcosz[i] = shwlist[i]->Direction().Z();
572  shwstartx[i] = shwlist[i]->ShowerStart().X();
573  shwstarty[i] = shwlist[i]->ShowerStart().Y();
574  shwstartz[i] = shwlist[i]->ShowerStart().Z();
575  for (size_t j = 0; j<(shwlist[i]->Energy()).size(); ++j){
576  shwenergy[i][j] = shwlist[i]->Energy()[j];
577  }
578  for (size_t j = 0; j<(shwlist[i]->dEdx()).size(); ++j){
579  shwdedx[i][j] = shwlist[i]->dEdx()[j];
580  }
581  shwbestplane[i] = shwlist[i]->best_plane();
582  if (fmsh.isValid()){
583  auto vhit = fmsh.at(i);
584  for (size_t h = 0; h < vhit.size(); ++h){
585  if (vhit[h].key()<kMaxHits){
586  hit_shwkey[vhit[h].key()] = shwlist[i].key();
587  }
588  }
589  }
590  if (!isdata&&fmsh.isValid()){
591  // Find true track for each reconstructed track
592  int TrackID = 0;
593  std::vector< art::Ptr<recob::Hit> > allHits = fmsh.at(i);
594  std::map<int,double> trkide;
595  for(size_t h = 0; h < allHits.size(); ++h){
596  art::Ptr<recob::Hit> hit = allHits[h];
597  std::vector<sim::TrackIDE> TrackIDs = bt_serv->HitToTrackIDEs(clockData, hit);
598  for(size_t e = 0; e < TrackIDs.size(); ++e){
599  trkide[TrackIDs[e].trackID] += TrackIDs[e].energy;
600  }
601  }
602  // Work out which IDE despoited the most charge in the hit if there was more than one.
603  double maxe = -1;
604  double tote = 0;
605  for (std::map<int,double>::iterator ii = trkide.begin(); ii!=trkide.end(); ++ii){
606  tote += ii->second;
607  if ((ii->second)>maxe){
608  maxe = ii->second;
609  TrackID = ii->first;
610  }
611  }
612  // Now have trackID, so get PdG code and T0 etc.
613  const simb::MCParticle *particle = pi_serv->TrackIdToParticle_P(TrackID);
614  if (particle){
615  shwg4id[i] = TrackID;
616  }
617  }
618  }
619  }
620  // flash information
621  flash_total = flashlist.size();
622  for ( int f = 0; f < std::min(flash_total,kMaxHits); ++f ) {
623  flash_time[f] = flashlist[f]->Time();
624  flash_width[f] = flashlist[f]->TimeWidth();
625  flash_abstime[f] = flashlist[f]->AbsTime();
626  flash_YCenter[f] = flashlist[f]->YCenter();
627  flash_YWidth[f] = flashlist[f]->YWidth();
628  flash_ZCenter[f] = flashlist[f]->ZCenter();
629  flash_ZWidth[f] = flashlist[f]->ZWidth();
630  flash_TotalPE[f] = flashlist[f]->TotalPE();
631  }
632 
633 
634 
635  if (!isdata){
636 
637  // * MC truth information
638  std::vector<art::Ptr<simb::MCTruth> > mclist;
639  auto mctruthListHandle = evt.getHandle< std::vector<simb::MCTruth> >(fGenieGenModuleLabel);
640  if (mctruthListHandle)
641  art::fill_ptr_vector(mclist, mctruthListHandle);
642 
643  std::vector<art::Ptr<simb::MCFlux> > fluxlist;
644  auto mcfluxListHandle = evt.getHandle< std::vector<simb::MCFlux> >(fGenieGenModuleLabel);
645  if (mcfluxListHandle)
646  art::fill_ptr_vector(fluxlist, mcfluxListHandle);
647 
648 
649  mcevts_truth=mclist.size();
650  if (mcevts_truth){
651  art::Ptr<simb::MCTruth> mctruth = mclist[0];
652  if (mctruth->Origin() == simb::kBeamNeutrino){
653  nuPDG_truth = mctruth->GetNeutrino().Nu().PdgCode();
654  ccnc_truth = mctruth->GetNeutrino().CCNC();
655  mode_truth = mctruth->GetNeutrino().Mode();
656  Q2_truth = mctruth->GetNeutrino().QSqr();
657  W_truth = mctruth->GetNeutrino().W();
658  X_truth = mctruth->GetNeutrino().X();
659  Y_truth = mctruth->GetNeutrino().Y();
660  hitnuc_truth = mctruth->GetNeutrino().HitNuc();
661  target_truth = mctruth->GetNeutrino().Target();
662  enu_truth = mctruth->GetNeutrino().Nu().E();
663  nuvtxx_truth = mctruth->GetNeutrino().Nu().Vx();
664  nuvtxy_truth = mctruth->GetNeutrino().Nu().Vy();
665  nuvtxz_truth = mctruth->GetNeutrino().Nu().Vz();
666  if (mctruth->GetNeutrino().Nu().P()){
667  nu_dcosx_truth = mctruth->GetNeutrino().Nu().Px()/mctruth->GetNeutrino().Nu().P();
668  nu_dcosy_truth = mctruth->GetNeutrino().Nu().Py()/mctruth->GetNeutrino().Nu().P();
669  nu_dcosz_truth = mctruth->GetNeutrino().Nu().Pz()/mctruth->GetNeutrino().Nu().P();
670  }
671  lep_mom_truth = mctruth->GetNeutrino().Lepton().P();
672  if (mctruth->GetNeutrino().Lepton().P()){
673  lep_dcosx_truth = mctruth->GetNeutrino().Lepton().Px()/mctruth->GetNeutrino().Lepton().P();
674  lep_dcosy_truth = mctruth->GetNeutrino().Lepton().Py()/mctruth->GetNeutrino().Lepton().P();
675  lep_dcosz_truth = mctruth->GetNeutrino().Lepton().Pz()/mctruth->GetNeutrino().Lepton().P();
676  }
677 
678  if (mctruth->NParticles()){
679  simb::MCParticle particle = mctruth->GetParticle(0);
680  t0_truth = particle.T();
681  }
682 
683 
684  float mindist2 = 9999; // cm;
685  TVector3 nuvtx(nuvtxx_truth, nuvtxy_truth, nuvtxz_truth);
687  //find the closest reco vertex to the neutrino mc truth
688  if (infidvol)
689  {
690  // vertex is when at least two tracks meet
691  for(size_t i = 0; i < vtxlist.size(); ++i){ // loop over vertices
692  Double_t xyz[3] = {};
693  vtxlist[i]->XYZ(xyz);
694  TVector3 vtxreco(xyz);
695  float dist2 = pma::Dist2(vtxreco, nuvtx);
696  if (dist2 < mindist2)
697  {
698  mindist2 = dist2;
699  vtxrecomc = std::sqrt(dist2);
700  vtxrecomcx = vtxreco.X() - nuvtxx_truth;
701  vtxrecomcy = vtxreco.Y() - nuvtxy_truth;
702  vtxrecomcz = vtxreco.Z() - nuvtxz_truth;
703  }
704  }
705 
706  // two endpoints of tracks are somehow also vertices...
707  for (size_t i = 0; i < tracklist.size(); ++i){ // loop over tracks
708  float dist2 = pma::Dist2(tracklist[i]->Vertex(), nuvtx);
709  if (dist2 < mindist2)
710  {
711  mindist2 = dist2;
712  vtxrecomc = std::sqrt(dist2);
713  vtxrecomcx = tracklist[i]->Vertex().X() - nuvtxx_truth;
714  vtxrecomcy = tracklist[i]->Vertex().Y() - nuvtxy_truth;
715  vtxrecomcz = tracklist[i]->Vertex().Z() - nuvtxz_truth;
716 
717  }
718  dist2 = pma::Dist2(tracklist[i]->End(), nuvtx);
719  if (dist2 < mindist2)
720  {
721  mindist2 = dist2;
722  vtxrecomc = std::sqrt(dist2);
723  vtxrecomcx = tracklist[i]->End().X() - nuvtxx_truth;
724  vtxrecomcy = tracklist[i]->End().Y() - nuvtxy_truth;
725  vtxrecomcz = tracklist[i]->End().Z() - nuvtxz_truth;
726 
727  }
728  }
729  }
730  }//is neutrino
731  }
732 
733  if (fluxlist.size()){
734  ptype_flux = fluxlist[0]->fptype;
735  pdpx_flux = fluxlist[0]->fpdpx;
736  pdpy_flux = fluxlist[0]->fpdpy;
737  pdpz_flux = fluxlist[0]->fpdpz;
738  pntype_flux = fluxlist[0]->fntype;
739  vx_flux = fluxlist[0]->fvx;
740  vy_flux = fluxlist[0]->fvy;
741  vz_flux = fluxlist[0]->fvz;
742  }
743 
744  //save g4 particle information
745  std::vector<const simb::MCParticle* > geant_part;
746 
747  // ### Looping over all the Geant4 particles from the BackTrackerService ###
748  for(size_t p = 0; p < plist.size(); ++p)
749  {
750  // ### Filling the vector with MC Particles ###
751  geant_part.push_back(plist.Particle(p));
752  }
753 
754  //std::cout<<"No of geant part= "<<geant_part.size()<<std::endl;
755 
756  // ### Setting a string for primary ###
757  std::string pri("primary");
758 
759  int primary=0;
760  int geant_particle=0;
761 
762  // ############################################################
763  // ### Determine the number of primary particles from geant ###
764  // ############################################################
765  for( unsigned int i = 0; i < geant_part.size(); ++i ){
766  geant_particle++;
767  // ### Counting the number of primary particles ###
768  if(geant_part[i]->Process()==pri)
769  { primary++;}
770  }//<---End i loop
771 
772 
773  // ### Saving the number of primary particles ###
774  no_primaries=primary;
775  // ### Saving the number of Geant4 particles ###
776  geant_list_size=geant_particle;
777 
778  // ### Looping over all the Geant4 particles ###
779  for( unsigned int i = 0; i < geant_part.size(); ++i ){
780 
781  // ### If this particle is primary, set = 1 ###
782  if(geant_part[i]->Process()==pri)
783  {process_primary[i]=1;}
784  // ### If this particle is not-primary, set = 0 ###
785  else
786  {process_primary[i]=0;}
787 
788  // ### Saving the particles mother TrackID ###
789  Mother[i]=geant_part[i]->Mother();
790  // ### Saving the particles TrackID ###
791  TrackId[i]=geant_part[i]->TrackId();
792  // ### Saving the PDG Code ###
793  pdg[i]=geant_part[i]->PdgCode();
794  // ### Saving the particles Energy ###
795  Eng[i]=geant_part[i]->E();
796 
797  // ### Saving the Px, Py, Pz info ###
798  Px[i]=geant_part[i]->Px();
799  Py[i]=geant_part[i]->Py();
800  Pz[i]=geant_part[i]->Pz();
801 
802  // ### Saving the Start and End Point for this particle ###
803  StartPointx[i]=geant_part[i]->Vx();
804  StartPointy[i]=geant_part[i]->Vy();
805  StartPointz[i]=geant_part[i]->Vz();
806  EndPointx[i]=geant_part[i]->EndPosition()[0];
807  EndPointy[i]=geant_part[i]->EndPosition()[1];
808  EndPointz[i]=geant_part[i]->EndPosition()[2];
809 
810  // ### Saving the processes for this particle ###
811  //std::cout<<"finding proc"<<std::endl;
812  G4Process.push_back( geant_part[i]->Process() );
813  G4FinalProcess.push_back( geant_part[i]->EndProcess() );
814  //std::cout<<"found proc"<<std::endl;
815 // std::cout << "ID " << TrackId[i] << ", pdg " << pdg[i] << ", Start X,Y,Z " << StartPointx[i] << ", " << StartPointy[i] << ", " << StartPointz[i]
816 // << ", End XYZ " << EndPointx[i] << ", " << EndPointy[i] << ", " << EndPointz[i] << ", Start Proc " << G4Process[i] << ", End Proc " << G4FinalProcess[i]
817 // << std::endl;
818 
819  // ### Saving the Start direction cosines for this particle ###
820  Startdcosx[i] = geant_part[i]->Momentum(0).Px() / geant_part[i]->Momentum(0).P();
821  Startdcosy[i] = geant_part[i]->Momentum(0).Py() / geant_part[i]->Momentum(0).P();
822  Startdcosz[i] = geant_part[i]->Momentum(0).Pz() / geant_part[i]->Momentum(0).P();
823  // ### Saving the number of Daughters for this particle ###
824  NumberDaughters[i]=geant_part[i]->NumberDaughters();
825 
826  } //geant particles
827 
828 
829  }//is neutrino
830  fTree->Fill();
831 }
832 
834 {
835  // Implementation of optional member function here.
837  fTree = tfs->make<TTree>("nueana","analysis tree");
838  fTree->Branch("run",&run,"run/I");
839  fTree->Branch("subrun",&subrun,"subrun/I");
840  fTree->Branch("event",&event,"event/I");
841  fTree->Branch("evttime",&evttime,"evttime/F");
842  fTree->Branch("taulife",&taulife,"taulife/F");
843  fTree->Branch("isdata",&isdata,"isdata/S");
844  fTree->Branch("ntracks_reco",&ntracks_reco,"ntracks_reco/I");
845  fTree->Branch("trkid",trkid,"trkid[ntracks_reco]/I");
846  fTree->Branch("trkstartx",trkstartx,"trkstartx[ntracks_reco]/F");
847  fTree->Branch("trkstarty",trkstarty,"trkstarty[ntracks_reco]/F");
848  fTree->Branch("trkstartz",trkstartz,"trkstartz[ntracks_reco]/F");
849  fTree->Branch("trkendx",trkendx,"trkendx[ntracks_reco]/F");
850  fTree->Branch("trkendy",trkendy,"trkendy[ntracks_reco]/F");
851  fTree->Branch("trkendz",trkendz,"trkendz[ntracks_reco]/F");
852  fTree->Branch("trkstartdcosx",trkstartdcosx,"trkstartdcosx[ntracks_reco]/F");
853  fTree->Branch("trkstartdcosy",trkstartdcosy,"trkstartdcosy[ntracks_reco]/F");
854  fTree->Branch("trkstartdcosz",trkstartdcosz,"trkstartdcosz[ntracks_reco]/F");
855  fTree->Branch("trkenddcosx",trkenddcosx,"trkenddcosx[ntracks_reco]/F");
856  fTree->Branch("trkenddcosy",trkenddcosy,"trkenddcosy[ntracks_reco]/F");
857  fTree->Branch("trkenddcosz",trkenddcosz,"trkenddcosz[ntracks_reco]/F");
858  fTree->Branch("trklen",trklen,"trklen[ntracks_reco]/F");
859  fTree->Branch("trkbestplane",trkbestplane,"trkbestplane[ntracks_reco]/I");
860  fTree->Branch("trkmomrange",trkmomrange,"trkmomrange[ntracks_reco]/F");
861  fTree->Branch("trkke",trkke,"trkke[ntracks_reco][3]/F");
862  fTree->Branch("trkpida",trkpida,"trkpida[ntracks_reco][3]/F");
863  fTree->Branch("trkg4id",trkg4id,"trkg4id[ntracks_reco]/I");
864  fTree->Branch("trkg4pdg",trkg4pdg,"trkg4pdg[ntracks_reco]/I");
865  fTree->Branch("trkg4startx",trkg4startx,"trkg4startx[ntracks_reco]/F");
866  fTree->Branch("trkg4starty",trkg4starty,"trkg4starty[ntracks_reco]/F");
867  fTree->Branch("trkg4startz",trkg4startz,"trkg4startz[ntracks_reco]/F");
868  fTree->Branch("trkg4initdedx",trkg4initdedx,"trkg4initdedx[ntracks_reco]/F");
869  fTree->Branch("nshws",&nshws,"nshws/I");
870  fTree->Branch("shwid",shwid,"shwid[nshws]/I");
871  fTree->Branch("shwdcosx",shwdcosx,"shwdcosx[nshws]/F");
872  fTree->Branch("shwdcosy",shwdcosy,"shwdcosy[nshws]/F");
873  fTree->Branch("shwdcosz",shwdcosz,"shwdcosz[nshws]/F");
874  fTree->Branch("shwstartx",shwstartx,"shwstartx[nshws]/F");
875  fTree->Branch("shwstarty",shwstarty,"shwstarty[nshws]/F");
876  fTree->Branch("shwstartz",shwstartz,"shwstartz[nshws]/F");
877  fTree->Branch("shwenergy",shwenergy,"shwenergy[nshws][3]/F");
878  fTree->Branch("shwdedx",shwdedx,"shwdedx[nshws][3]/F");
879  fTree->Branch("shwbestplane",shwbestplane,"shwbestplane[nshws]/I");
880  fTree->Branch("shwg4id",shwg4id,"shwg4id[nshws]/I");
881  fTree->Branch("flash_total" ,&flash_total ,"flash_total/I");
882  fTree->Branch("flash_time" ,flash_time ,"flash_time[flash_total]/F");
883  fTree->Branch("flash_width" ,flash_width ,"flash_width[flash_total]/F");
884  fTree->Branch("flash_abstime",flash_abstime,"flash_abstime[flash_total]/F");
885  fTree->Branch("flash_YCenter",flash_YCenter,"flash_YCenter[flash_total]/F");
886  fTree->Branch("flash_YWidth" ,flash_YWidth ,"flash_YWidth[flash_total]/F");
887  fTree->Branch("flash_ZCenter",flash_ZCenter,"flash_ZCenter[flash_total]/F");
888  fTree->Branch("flash_ZWidth" ,flash_ZWidth ,"flash_ZWidth[flash_total]/F");
889  fTree->Branch("flash_TotalPE",flash_TotalPE,"flash_TotalPE[flash_total]/F");
890  fTree->Branch("nhits",&nhits,"nhits/I");
891  fTree->Branch("nhits_stored",&nhits_stored,"nhits_stored/I");
892  fTree->Branch("hit_plane",hit_plane,"hit_plane[nhits_stored]/S");
893  fTree->Branch("hit_tpc",hit_tpc,"hit_tpc[nhits_stored]/S");
894  fTree->Branch("hit_wire",hit_wire,"hit_wire[nhits_stored]/S");
895  fTree->Branch("hit_channel",hit_channel,"hit_channel[nhits_stored]/I");
896  fTree->Branch("hit_peakT",hit_peakT,"hit_peakT[nhits_stored]/F");
897  fTree->Branch("hit_charge",hit_charge,"hit_charge[nhits_stored]/F");
898  fTree->Branch("hit_summedADC",hit_summedADC,"hit_summedADC[nhits_stored]/F");
899  fTree->Branch("hit_startT",hit_startT,"hit_startT[nhits_stored]/F");
900  fTree->Branch("hit_endT",hit_endT,"hit_endT[nhits_stored]/F");
901  fTree->Branch("hit_trkkey",hit_trkkey,"hit_trkkey[nhits_stored]/I");
902  fTree->Branch("hit_dQds",hit_dQds,"hit_dQds[nhits_stored]/F");
903  fTree->Branch("hit_dEds",hit_dEds,"hit_dEds[nhits_stored]/F");
904  fTree->Branch("hit_resrange",hit_resrange,"hit_resrange[nhits_stored]/F");
905  fTree->Branch("hit_shwkey",hit_shwkey,"hit_shwkey[nhits_stored]/I");
906  fTree->Branch("infidvol",&infidvol,"infidvol/I");
907  fTree->Branch("nvtx",&nvtx,"nvtx/S");
908  fTree->Branch("vtx",vtx,"vtx[nvtx][3]/F");
909  fTree->Branch("vtxrecomc",&vtxrecomc,"vtxrecomc/F");
910  fTree->Branch("vtxrecomcx",&vtxrecomcx,"vtxrecomcx/F");
911  fTree->Branch("vtxrecomcy",&vtxrecomcy,"vtxrecomcy/F");
912  fTree->Branch("vtxrecomcz",&vtxrecomcz,"vtxrecomcz/F");
913  fTree->Branch("mcevts_truth",&mcevts_truth,"mcevts_truth/I");
914  fTree->Branch("nuPDG_truth",&nuPDG_truth,"nuPDG_truth/I");
915  fTree->Branch("ccnc_truth",&ccnc_truth,"ccnc_truth/I");
916  fTree->Branch("mode_truth",&mode_truth,"mode_truth/I");
917  fTree->Branch("enu_truth",&enu_truth,"enu_truth/F");
918  fTree->Branch("Q2_truth",&Q2_truth,"Q2_truth/F");
919  fTree->Branch("W_truth",&W_truth,"W_truth/F");
920  fTree->Branch("X_truth",&X_truth,"X_truth/F");
921  fTree->Branch("Y_truth",&Y_truth,"Y_truth/F");
922  fTree->Branch("hitnuc_truth",&hitnuc_truth,"hitnuc_truth/I");
923  fTree->Branch("target_truth",&target_truth,"target_truth/I");
924  fTree->Branch("nuvtxx_truth",&nuvtxx_truth,"nuvtxx_truth/F");
925  fTree->Branch("nuvtxy_truth",&nuvtxy_truth,"nuvtxy_truth/F");
926  fTree->Branch("nuvtxz_truth",&nuvtxz_truth,"nuvtxz_truth/F");
927  fTree->Branch("nu_dcosx_truth",&nu_dcosx_truth,"nu_dcosx_truth/F");
928  fTree->Branch("nu_dcosy_truth",&nu_dcosy_truth,"nu_dcosy_truth/F");
929  fTree->Branch("nu_dcosz_truth",&nu_dcosz_truth,"nu_dcosz_truth/F");
930  fTree->Branch("lep_mom_truth",&lep_mom_truth,"lep_mom_truth/F");
931  fTree->Branch("lep_dcosx_truth",&lep_dcosx_truth,"lep_dcosx_truth/F");
932  fTree->Branch("lep_dcosy_truth",&lep_dcosy_truth,"lep_dcosy_truth/F");
933  fTree->Branch("lep_dcosz_truth",&lep_dcosz_truth,"lep_dcosz_truth/F");
934  fTree->Branch("t0_truth",&t0_truth,"t0_truth/F");
935  fTree->Branch("no_primaries",&no_primaries,"no_primaries/I");
936  fTree->Branch("geant_list_size",&geant_list_size,"geant_list_size/I");
937  fTree->Branch("pdg",pdg,"pdg[geant_list_size]/I");
938  fTree->Branch("Eng",Eng,"Eng[geant_list_size]/F");
939  fTree->Branch("Px",Px,"Px[geant_list_size]/F");
940  fTree->Branch("Py",Py,"Py[geant_list_size]/F");
941  fTree->Branch("Pz",Pz,"Pz[geant_list_size]/F");
942  fTree->Branch("StartPointx",StartPointx,"StartPointx[geant_list_size]/F");
943  fTree->Branch("StartPointy",StartPointy,"StartPointy[geant_list_size]/F");
944  fTree->Branch("StartPointz",StartPointz,"StartPointz[geant_list_size]/F");
945  fTree->Branch("EndPointx",EndPointx,"EndPointx[geant_list_size]/F");
946  fTree->Branch("EndPointy",EndPointy,"EndPointy[geant_list_size]/F");
947  fTree->Branch("EndPointz",EndPointz,"EndPointz[geant_list_size]/F");
948  fTree->Branch("Startdcosx",Startdcosx,"Startdcosx[geant_list_size]/F");
949  fTree->Branch("Startdcosy",Startdcosy,"Startdcosy[geant_list_size]/F");
950  fTree->Branch("Startdcosz",Startdcosz,"Startdcosz[geant_list_size]/F");
951  fTree->Branch("NumberDaughters",NumberDaughters,"NumberDaughters[geant_list_size]/I");
952  fTree->Branch("Mother",Mother,"Mother[geant_list_size]/I");
953  fTree->Branch("TrackId",TrackId,"TrackId[geant_list_size]/I");
954  fTree->Branch("process_primary",process_primary,"process_primary[geant_list_size]/I");
955  fTree->Branch("G4Process",&G4Process);//,"G4Process[geant_list_size]");
956  fTree->Branch("G4FinalProcess",&G4FinalProcess);//,"G4FinalProcess[geant_list_size]");
957  fTree->Branch("ptype_flux",&ptype_flux,"ptype_flux/I");
958  fTree->Branch("pdpx_flux",&pdpx_flux,"pdpx_flux/F");
959  fTree->Branch("pdpy_flux",&pdpy_flux,"pdpy_flux/F");
960  fTree->Branch("pdpz_flux",&pdpz_flux,"pdpz_flux/F");
961  fTree->Branch("pntype_flux",&pntype_flux,"pntype_flux/I");
962  fTree->Branch("vx_flux",&vx_flux,"vx_flux/F");
963  fTree->Branch("vy_flux",&vy_flux,"vy_flux/F");
964  fTree->Branch("vz_flux",&vz_flux,"vz_flux/F");
965 
966  fPOT = tfs->make<TTree>("pottree","pot tree");
967  fPOT->Branch("pot",&pot,"pot/D");
968  fPOT->Branch("run",&run,"run/I");
969  fPOT->Branch("subrun",&subrun,"subrun/I");
970 }
971 
973 
974  G4Process.clear();
975  G4FinalProcess.clear();
976 
977  run = -9999;
978  subrun = -9999;
979  event = -9999;
980  evttime = -9999;
981  taulife = 0;
982  isdata = -9999;
983 
984  ntracks_reco = 0;
985  for (int i = 0; i < kMaxTrack; ++i){
986  trkid[i] = -9999;
987  trkstartx[i] = -9999;
988  trkstarty[i] = -9999;
989  trkstartz[i] = -9999;
990  trkendx[i] = -9999;
991  trkendy[i] = -9999;
992  trkendz[i] = -9999;
993  trkstartdcosx[i] = -9999;
994  trkstartdcosy[i] = -9999;
995  trkstartdcosz[i] = -9999;
996  trkenddcosx[i] = -9999;
997  trkenddcosy[i] = -9999;
998  trkenddcosz[i] = -9999;
999  trklen[i] = -9999;
1000  trkbestplane[i] = -9999;
1001  trkmomrange[i] = -9999;
1002  for (int j = 0; j<3; ++j){
1003  trkke[i][j] = -9999;
1004  trkpida[i][j] = -9999;
1005  }
1006  trkg4id[i] = -9999;
1007  trkg4pdg[i] = -9999;
1008  trkg4startx[i] = -9999;
1009  trkg4starty[i] = -9999;
1010  trkg4startz[i] = -9999;
1011  trkg4initdedx[i] = -9999;
1012  }
1013 
1014  nshws = 0;
1015  for (int i = 0; i<kMaxShower; ++i){
1016  shwid[i] = -9999;
1017  shwdcosx[i] = -9999;
1018  shwdcosy[i] = -9999;
1019  shwdcosz[i] = -9999;
1020  shwstartx[i] = -9999;
1021  shwstarty[i] = -9999;
1022  shwstartz[i] = -9999;
1023  for (int j = 0; j<3; ++j){
1024  shwenergy[i][j] = -9999;
1025  shwdedx[i][j] = -9999;
1026  }
1027  shwbestplane[i] = -9999;
1028  trkg4id[i] = -9999;
1029  }
1030 
1031  flash_total = 0;
1032  for (int f = 0; f < kMaxFlash; ++f) {
1033  flash_time[f] = -9999;
1034  flash_width[f] = -9999;
1035  flash_abstime[f] = -9999;
1036  flash_YCenter[f] = -9999;
1037  flash_YWidth[f] = -9999;
1038  flash_ZCenter[f] = -9999;
1039  flash_ZWidth[f] = -9999;
1040  flash_TotalPE[f] = -9999;
1041  }
1042 
1043  nhits = 0;
1044  nhits_stored = 0;
1045  for (int i = 0; i<kMaxHits; ++i){
1046  hit_plane[i] = -9999;
1047  hit_wire[i] = -9999;
1048  hit_tpc[i] = -9999;
1049  hit_channel[i] = -9999;
1050  hit_peakT[i] = -9999;
1051  hit_charge[i] = -9999;
1052  hit_summedADC[i] = -9999;
1053  hit_startT[i] = -9999;
1054  hit_endT[i] = -9999;
1055  hit_trkkey[i] = -9999;
1056  hit_dQds[i] = -9999;
1057  hit_dEds[i] = -9999;
1058  hit_resrange[i] = -9999;
1059  hit_shwkey[i] = -9999;
1060  }
1061 
1062  infidvol = 0;
1063  nvtx = 0;
1064  for (int i = 0; i<kMaxVertices; ++i){
1065  vtx[i][0] = -9999;
1066  vtx[i][1] = -9999;
1067  vtx[i][2] = -9999;
1068  }
1069  vtxrecomc = 9999;
1070  vtxrecomcx = 9999;
1071  vtxrecomcy = 9999;
1072  vtxrecomcz = 9999;
1073 
1074  mcevts_truth = -9999;
1075  nuPDG_truth = -9999;
1076  ccnc_truth = -9999;
1077  mode_truth = -9999;
1078  enu_truth = -9999;
1079  Q2_truth = -9999;
1080  W_truth = -9999;
1081  X_truth = -9999;
1082  Y_truth = -9999;
1083  hitnuc_truth = -9999;
1084  target_truth = -9999;
1085  nuvtxx_truth = -9999;
1086  nuvtxy_truth = -9999;
1087  nuvtxz_truth = -9999;
1088  nu_dcosx_truth = -9999;
1089  nu_dcosy_truth = -9999;
1090  nu_dcosz_truth = -9999;
1091  lep_mom_truth = -9999;
1092  lep_dcosx_truth = -9999;
1093  lep_dcosy_truth = -9999;
1094  lep_dcosz_truth = -9999;
1095  t0_truth = -9999;
1096 
1097  no_primaries = -99999;
1098  geant_list_size=-9999;
1099  for (int i = 0; i<kMaxPrimaries; ++i){
1100  pdg[i] = -99999;
1101  Eng[i] = -99999;
1102  Px[i] = -99999;
1103  Py[i] = -99999;
1104  Pz[i] = -99999;
1105  StartPointx[i] = -99999;
1106  StartPointy[i] = -99999;
1107  StartPointz[i] = -99999;
1108  EndPointx[i] = -99999;
1109  EndPointy[i] = -99999;
1110  EndPointz[i] = -99999;
1111  Startdcosx[i]= -99999;
1112  Startdcosy[i]= -99999;
1113  Startdcosz[i]= -99999;
1114  NumberDaughters[i] = -99999;
1115  Mother[i] = -99999;
1116  TrackId[i] = -99999;
1117  process_primary[i] = -99999;
1118  }
1119 
1120  ptype_flux = -99999;
1121  pdpx_flux = -99999;
1122  pdpy_flux = -99999;
1123  pdpz_flux = -99999;
1124  pntype_flux = -99999;
1125  vx_flux = -99999;
1126  vy_flux = -99999;
1127  vz_flux = -99999;
1128 }
1129 
1131 {
1132  // Implementation of optional member function here.
1133 }
1134 
1136 {
1137  fHitsModuleLabel = p.get< std::string >("HitsModuleLabel");
1138  fTrackModuleLabel = p.get< std::string >("TrackModuleLabel");
1139  fShowerModuleLabel = p.get< std::string >("ShowerModuleLabel");
1140  fClusterModuleLabel = p.get< std::string >("ClusterModuleLabel");
1141  fVertexModuleLabel = p.get< std::string >("VertexModuleLabel");
1142  fGenieGenModuleLabel = p.get< std::string >("GenieGenModuleLabel");
1143  fFidVolCut = p.get< double >("FidVolCut");
1144  fPOTModuleLabel = p.get< std::string >("POTModuleLabel");
1145  fFlashModuleLabel = p.get< std::string >("FlashModuleLabel");
1146  fCalorimetryModuleLabel = p.get< std::string >("CalorimetryModuleLabel");
1147  return;
1148 }
1149 
1150 /***********************************************************************/
1151 
1152 bool dunefd::NueAna::insideFidVol(const double posX, const double posY, const double posZ)
1153 {
1154 
1156  double vtx[3] = {posX, posY, posZ};
1157  bool inside = false;
1158 
1159  geo::TPCID idtpc = geom->FindTPCAtPosition(vtx);
1160 
1161  if (geom->HasTPC(idtpc))
1162  {
1163  const geo::TPCGeo& tpcgeo = geom->GetElement(idtpc);
1164  double minx = tpcgeo.MinX(); double maxx = tpcgeo.MaxX();
1165  double miny = tpcgeo.MinY(); double maxy = tpcgeo.MaxY();
1166  double minz = tpcgeo.MinZ(); double maxz = tpcgeo.MaxZ();
1167 
1168  for (size_t c = 0; c < geom->Ncryostats(); c++)
1169  {
1170  const geo::CryostatGeo& cryostat = geom->Cryostat(c);
1171  for (size_t t = 0; t < cryostat.NTPC(); t++)
1172  {
1173  const geo::TPCGeo& tpcg = cryostat.TPC(t);
1174  if (tpcg.MinX() < minx) minx = tpcg.MinX();
1175  if (tpcg.MaxX() > maxx) maxx = tpcg.MaxX();
1176  if (tpcg.MinY() < miny) miny = tpcg.MinY();
1177  if (tpcg.MaxY() > maxy) maxy = tpcg.MaxY();
1178  if (tpcg.MinZ() < minz) minz = tpcg.MinZ();
1179  if (tpcg.MaxZ() > maxz) maxz = tpcg.MaxZ();
1180  }
1181  }
1182 
1183 
1184  //x
1185  double dista = fabs(minx - posX);
1186  double distb = fabs(posX - maxx);
1187  if ((posX > minx) && (posX < maxx) &&
1188  (dista > fFidVolCut) && (distb > fFidVolCut)) inside = true;
1189  //y
1190  dista = fabs(maxy - posY);
1191  distb = fabs(posY - miny);
1192  if (inside && (posY > miny) && (posY < maxy) &&
1193  (dista > fFidVolCut) && (distb > fFidVolCut)) inside = true;
1194  else inside = false;
1195 
1196  //z
1197  dista = fabs(maxz - posZ);
1198  distb = fabs(posZ - minz);
1199  if (inside && (posZ > minz) && (posZ < maxz) &&
1200  (dista > fFidVolCut) && (distb > fFidVolCut)) inside = true;
1201  else inside = false;
1202  }
1203 
1204  return inside;
1205 }
1206 
1207 /***********************************************************************/
1208 
1209 
double E(const int i=0) const
Definition: MCParticle.h:233
Float_t nuvtxx_truth
bool insideFidVol(const double posX, const double posY, const double posZ)
intermediate_table::iterator iterator
std::string fCalorimetryModuleLabel
Float_t hit_startT[kMaxHits]
int TrackId[kMaxPrimaries]
void beginJob() override
constexpr int kMaxFlash
int PdgCode() const
Definition: MCParticle.h:212
int CCNC() const
Definition: MCNeutrino.h:148
constexpr std::uint32_t timeLow() const
Definition: Timestamp.h:29
double QSqr() const
Definition: MCNeutrino.h:157
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:77
double Py(const int i=0) const
Definition: MCParticle.h:231
CryostatGeo const & GetElement(geo::CryostatID const &cryoid) const
float trkendy[kMaxTrack]
Float_t flash_YWidth[kMaxFlash]
float trkenddcosy[kMaxTrack]
float trkstartdcosz[kMaxTrack]
std::vector< sim::TrackIDE > HitToTrackIDEs(detinfo::DetectorClocksData const &clockData, recob::Hit const &hit) const
std::string fClusterModuleLabel
Short_t hit_wire[kMaxHits]
int Mother[kMaxPrimaries]
double Dist2(const TVector2 &v1, const TVector2 &v2)
Definition: Utilities.cxx:37
float trkendz[kMaxTrack]
Float_t hit_dQds[kMaxHits]
float trkg4startx[kMaxTrack]
Float_t shwdedx[kMaxShower][3]
Float_t flash_ZWidth[kMaxFlash]
Handle< PROD > getHandle(SelectorBase const &) const
Definition: DataViewImpl.h:382
Float_t shwstartx[kMaxShower]
std::string string
Definition: nybbler.cc:12
void analyze(art::Event const &e) override
std::vector< TrackID > TrackIDs
const simb::MCParticle * TrackIdToParticle_P(int id) const
float trkg4startz[kMaxTrack]
geo::WireID WireID() const
Definition: Hit.h:233
Float_t shwdcosz[kMaxShower]
Float_t EndPointy[kMaxPrimaries]
Float_t hit_summedADC[kMaxHits]
Float_t flash_TotalPE[kMaxFlash]
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
constexpr std::uint32_t timeHigh() const
Definition: Timestamp.h:34
Float_t shwstarty[kMaxShower]
simb::Origin_t Origin() const
Definition: MCTruth.h:74
constexpr T pow(T x)
Definition: pow.h:72
void endJob() override
double MinX() const
Returns the world x coordinate of the start of the box.
Definition: BoxBoundedGeo.h:88
Geometry information for a single TPC.
Definition: TPCGeo.h:38
double Px(const int i=0) const
Definition: MCParticle.h:230
Short_t hit_plane[kMaxHits]
Float_t hit_dEds[kMaxHits]
int HitNuc() const
Definition: MCNeutrino.h:152
int trkg4id[kMaxTrack]
Float_t Pz[kMaxPrimaries]
Class to keep data related to recob::Hit associated with recob::Track.
Float_t lep_dcosy_truth
Float_t shwdcosx[kMaxShower]
float trkmomrange[kMaxTrack]
int shwbestplane[kMaxShower]
Int_t hit_trkkey[kMaxHits]
int pdg[kMaxPrimaries]
Planes which measure Z direction.
Definition: geo_types.h:132
Float_t flash_abstime[kMaxFlash]
int process_primary[kMaxPrimaries]
double MaxX() const
Returns the world x coordinate of the end of the box.
Definition: BoxBoundedGeo.h:91
int NParticles() const
Definition: MCTruth.h:75
float trkstartx[kMaxTrack]
float trkstartz[kMaxTrack]
unsigned int Index
Float_t hit_endT[kMaxHits]
Geometry information for a single cryostat.
Definition: CryostatGeo.h:43
std::string fVertexModuleLabel
EDAnalyzer(fhicl::ParameterSet const &pset)
Definition: EDAnalyzer.h:25
NueAna & operator=(NueAna const &)=delete
unsigned int Ncryostats() const
Returns the number of cryostats in the detector.
int NumberDaughters[kMaxPrimaries]
float trklen[kMaxTrack]
void reconfigure(fhicl::ParameterSet const &p)
int trkg4pdg[kMaxTrack]
constexpr int kMaxShower
void endSubRun(const art::SubRun &sr) override
object containing MC flux information
art framework interface to geometry description
int shwg4id[kMaxShower]
Float_t hit_charge[kMaxHits]
Float_t nu_dcosy_truth
tracking::Vector_t Vector_t
Definition: Track.h:54
Float_t nuvtxy_truth
geo::TPCID FindTPCAtPosition(double const worldLoc[3]) const
Returns the ID of the TPC at specified location.
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:92
float trkpida[kMaxTrack][3]
float trkg4starty[kMaxTrack]
bool isRealData() const
int trkid[kMaxTrack]
const double e
Float_t shwdcosy[kMaxShower]
const simb::MCParticle & Lepton() const
Definition: MCNeutrino.h:147
int shwid[kMaxShower]
float trkendx[kMaxTrack]
double W() const
Definition: MCNeutrino.h:154
Timestamp time() const
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:67
std::string fHitsModuleLabel
Float_t Eng[kMaxPrimaries]
constexpr int kMaxPrimaries
double Y() const
Definition: MCNeutrino.h:156
def key(type, name=None)
Definition: graph.py:13
std::string fPOTModuleLabel
Float_t nu_dcosx_truth
Float_t flash_time[kMaxFlash]
double P(const int i=0) const
Definition: MCParticle.h:234
key_type key() const noexcept
Definition: Ptr.h:216
T get(std::string const &key) const
Definition: ParameterSet.h:271
double T(const int i=0) const
Definition: MCParticle.h:224
double MinZ() const
Returns the world z coordinate of the start of the box.
p
Definition: test.py:223
unsigned int NTPC() const
Number of TPCs in this cryostat.
Definition: CryostatGeo.h:181
CryostatGeo const & Cryostat(geo::CryostatID const &cryoid) const
Returns the specified cryostat.
SubRunNumber_t subRun() const
Definition: DataViewImpl.cc:78
double X() const
Definition: MCNeutrino.h:155
Float_t lep_mom_truth
RunNumber_t run() const
Definition: DataViewImpl.cc:71
The data type to uniquely identify a TPC.
Definition: geo_types.h:386
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:493
Float_t hit_peakT[kMaxHits]
std::string fGenieGenModuleLabel
float trkke[kMaxTrack][3]
Float_t lep_dcosx_truth
Float_t shwenergy[kMaxShower][3]
Float_t StartPointx[kMaxPrimaries]
Detector simulation of raw signals on wires.
double MaxY() const
Returns the world y coordinate of the end of the box.
const sim::ParticleList & ParticleList() const
const simb::MCParticle & GetParticle(int i) const
Definition: MCTruth.h:76
Float_t flash_ZCenter[kMaxFlash]
std::vector< std::string > G4FinalProcess
bool HasTPC(geo::TPCID const &tpcid) const
Returns whether we have the specified TPC.
float trkenddcosz[kMaxTrack]
std::string fTrackModuleLabel
Float_t hit_resrange[kMaxHits]
double Vx(const int i=0) const
Definition: MCParticle.h:221
Float_t vtx[kMaxVertices][3]
std::string fShowerModuleLabel
Implementation of the Projection Matching Algorithm.
std::string fFlashModuleLabel
Declaration of signal hit object.
calo::CalorimetryAlg fCalorimetryAlg
int Target() const
Definition: MCNeutrino.h:151
float trkstartdcosx[kMaxTrack]
Float_t Startdcosx[kMaxPrimaries]
Float_t Py[kMaxPrimaries]
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:55
constexpr int kMaxVertices
Float_t Startdcosy[kMaxPrimaries]
const TPCGeo & TPC(unsigned int itpc) const
Return the itpc&#39;th TPC in the cryostat.
Definition: CryostatGeo.cxx:93
void End(void)
Definition: gXSecComp.cxx:210
constexpr int kMaxTrack
int trkbestplane[kMaxTrack]
NueAna(fhicl::ParameterSet const &p)
constexpr int kMaxHits
Float_t shwstartz[kMaxShower]
Float_t flash_width[kMaxFlash]
double MaxZ() const
Returns the world z coordinate of the end of the box.
Float_t nuvtxz_truth
float trkg4initdedx[kMaxTrack]
tracking::Point_t Point_t
Definition: Track.h:53
double TrackPitchInView(recob::Track const &track, geo::View_t view, size_t trajectory_point=0U)
Returns the projected length of track on a wire pitch step [cm].
Definition: TrackUtils.cxx:78
double Pz(const int i=0) const
Definition: MCParticle.h:232
Provides recob::Track data product.
double dEdx_AREA(detinfo::DetectorClocksData const &clock_data, detinfo::DetectorPropertiesData const &det_prop, recob::Hit const &hit, double pitch, double T0=0) const
double Vz(const int i=0) const
Definition: MCParticle.h:223
Int_t hit_shwkey[kMaxHits]
Float_t EndPointx[kMaxPrimaries]
EventNumber_t event() const
Definition: EventID.h:116
Float_t StartPointz[kMaxPrimaries]
list x
Definition: train.py:276
boost::graph_traits< ModuleGraph >::vertex_descriptor Vertex
Definition: ModuleGraph.h:25
Float_t Px[kMaxPrimaries]
TCEvent evt
Definition: DataStructs.cxx:7
Float_t Startdcosz[kMaxPrimaries]
auto const & get(AssnsNode< L, R, D > const &r)
Definition: AssnsNode.h:115
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:297
Float_t flash_YCenter[kMaxFlash]
double LifetimeCorrection(detinfo::DetectorClocksData const &clock_data, detinfo::DetectorPropertiesData const &det_prop, double time, double T0=0) const
Float_t StartPointy[kMaxPrimaries]
static constexpr double sr
Definition: Units.h:166
Float_t lep_dcosz_truth
float trkstarty[kMaxTrack]
Float_t nu_dcosz_truth
double MinY() const
Returns the world y coordinate of the start of the box.
int Mode() const
Definition: MCNeutrino.h:149
float trkenddcosx[kMaxTrack]
Utility functions to extract information from recob::Track
Float_t EndPointz[kMaxPrimaries]
Int_t hit_channel[kMaxHits]
EventID id() const
Definition: Event.cc:34
double Vy(const int i=0) const
Definition: MCParticle.h:222
calorimetry
double GetTrackMomentum(double trkrange, int pdg) const
float trkstartdcosy[kMaxTrack]
Short_t hit_tpc[kMaxHits]
QTextStream & endl(QTextStream &s)
Event finding and building.
std::vector< std::string > G4Process
Beam neutrinos.
Definition: MCTruth.h:23