IniSegReco_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////////////////////////
2 // Class: IniSegAlg
3 // Module Type: producer
4 // File: IniSegReco_module.cc
5 // Authors: dorota.stefan@cern.ch, robert.sulej@cern.ch, tjyang@fnal.gov
6 //
7 ////////////////////////////////////////////////////////////////////////////////////////////////////
8 
9 #include "art/Framework/Core/EDAnalyzer.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 "fhiclcpp/ParameterSet.h"
18 #include "messagefacility/MessageLogger/MessageLogger.h"
19 
20 // LArSoft includes
21 #include "larcore/Geometry/Geometry.h"
22 #include "lardataobj/RecoBase/SpacePoint.h"
23 #include "lardataobj/RecoBase/Track.h"
24 #include "lardataobj/RecoBase/Hit.h"
25 #include "lardataobj/RecoBase/Cluster.h"
26 #include "lardataobj/RecoBase/Vertex.h"
27 #include "lardata/Utilities/AssociationUtil.h"
28 #include "larsim/MCCheater/ParticleInventoryService.h"
29 #include "nusimdata/SimulationBase/MCTruth.h"
30 #include "larreco/RecoAlg/ProjectionMatchingAlg.h"
31 #include "larreco/RecoAlg/PMAlg/PmaTrack3D.h"
32 #include "larreco/RecoAlg/PMAlg/Utilities.h"
33 #include "IniSegAlg/IniSegAlg.h"
34 #include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
35 
36 // ROOT includes
37 #include "TTree.h"
38 #include "TTimeStamp.h"
39 
40 // C/C++ libraries
41 #include <memory>
42 #include <utility>
43 #include <fstream>
44 namespace dunefd {
45  class IniSegReco;
46  class Hit2D;
47  class IniSegAlg;
48 }
49 
50 class dunefd::IniSegReco : public art::EDProducer {
51 public:
52  explicit IniSegReco(fhicl::ParameterSet const & p);
53 
54  IniSegReco(IniSegReco const &) = delete;
55  IniSegReco(IniSegReco &&) = delete;
56  IniSegReco & operator = (IniSegReco const &) = delete;
57  IniSegReco & operator = (IniSegReco &&) = delete;
58 
59  void beginJob() override;
60 
61  void reconfigure(fhicl::ParameterSet const& p) ;
62 
63  void produce(art::Event & e) override;
64 
65 
66 private:
67  void ResetVars(void);
68  void chargeParticlesatVtx(art::Event const & evt);
69 
70  bool insideFidVol(TLorentzVector const & pvtx) const;
71 
72  float t0Corr(art::Event const & evt,
73  detinfo::DetectorPropertiesData const& detProp,
74  TLorentzVector const & pvtx);
75 
76  TVector2 getMCvtx2d(TVector3 const & mcvtx3d,
77  const size_t cryo, const size_t tpc, const size_t plane) const;
78  TVector2 getMCdir2d(TVector3 const & mcvtx3d, TVector3 const & mcdir3d,
79  const size_t cryo, const size_t tpc, const size_t plane) const;
80  TVector3 findElDir(art::Ptr<simb::MCTruth> const mctruth) const;
81  std::vector< TVector3 > findPhDir() const;
82  std::vector< TVector3 > findDirs(art::Ptr<simb::MCTruth> const mctruth, int pdg) const;
83 
84  void collectCls(art::Event const & evt,
85  detinfo::DetectorPropertiesData const& detProp,
86  art::Ptr<simb::MCTruth> const mctruth);
87  std::vector< dunefd::Hit2D > reselectCls(std::map<size_t, std::vector< dunefd::Hit2D > > const & cls,
88  art::Ptr<simb::MCTruth> const mctruth,
89  const size_t cryo, const size_t tpc, const size_t plane) const;
90  void make3dseg(art::Event const & evt,
91  detinfo::DetectorPropertiesData const& detProp,
92  std::vector< std::vector<Hit2D> > const & src, TVector3 const & primary);
93 
94  recob::Track convertFrom(pma::Track3D const & src);
95 
96  void UseClusters(art::Event const & evt,
97  detinfo::DetectorPropertiesData const& detProp);
98 
99  void UseTracks(art::Event const & evt,
100  detinfo::DetectorPropertiesData const& detProp);
101 
102  std::vector< pma::Track3D* > pmatracks;
103 
104  TTree *fTree;
105 
106  int run;
107  int subrun;
108  int event;
109  short isdata;
110  int trkindex; // track index in the event
111 
112  Float_t lep_dedx;
113  Float_t dx;
114  Float_t lep_dist;
115  Float_t lep_distx;
116  Float_t lep_disty;
117  Float_t lep_distz;
118  Float_t lep_distreco;
119  Float_t lep_distrecox;
120  Float_t lep_distrecoy;
121  Float_t lep_distrecoz;
122  Float_t cos;
123  Float_t coslx; Float_t cosly; Float_t coslz;
124  Float_t coslrx; Float_t coslry; Float_t coslrz;
125  Float_t t0;
126  Float_t vtxrecomc;
127  Float_t vtxrecomcx;
128  Float_t vtxrecomcy;
129  Float_t vtxrecomcz;
130  Float_t vtxupstream;
131  Float_t vtxupstreamx;
132  Float_t vtxupstreamy;
133  Float_t vtxupstreamz;
134  int ngamma;
135  Float_t convdist;
136 
137  std::string fHitsModuleLabel;
138  std::string fClusterModuleLabel;
139  std::string fVertexModuleLabel;
140  std::string fTrackModuleLabel;
141  std::string fGenieGenModuleLabel;
142 
143  double fFidVolCut;
144  pma::ProjectionMatchingAlg fProjectionMatchingAlg;
145 
146  std::ofstream file;
147  std::ofstream file1;
148 };
149 // ------------------------------------------------------
150 
151 dunefd::IniSegReco::IniSegReco(fhicl::ParameterSet const & pset)
152 : EDProducer{pset}, fProjectionMatchingAlg(pset.get< fhicl::ParameterSet >("ProjectionMatchingAlg"))
153 {
154  this->reconfigure(pset);
155  produces< std::vector<recob::Track> >();
156  produces< std::vector<recob::SpacePoint> >();
157  produces< art::Assns<recob::Track, recob::Hit> >();
158  produces< art::Assns<recob::Track, recob::SpacePoint> >();
159  produces< art::Assns<recob::SpacePoint, recob::Hit> >();
160 }
161 // ------------------------------------------------------
162 
163 void dunefd::IniSegReco::beginJob()
164 {
165  // Implementation of optional member function here.
166  art::ServiceHandle<art::TFileService> tfs;
167  fTree = tfs->make<TTree>("nueana","analysis tree");
168  fTree->Branch("run",&run,"run/I");
169  fTree->Branch("subrun",&subrun,"subrun/I");
170  fTree->Branch("event",&event,"event/I");
171  fTree->Branch("lep_dedx",&lep_dedx,"lep_dedx/F");
172  fTree->Branch("dx",&dx,"dx/F");
173  fTree->Branch("lep_dist",&lep_dist,"lep_dist/F");
174  fTree->Branch("lep_distx",&lep_distx,"lep_distx/F");
175  fTree->Branch("lep_disty",&lep_disty,"lep_disty/F");
176  fTree->Branch("lep_distz",&lep_distz,"lep_distz/F");
177  fTree->Branch("lep_distreco",&lep_distreco,"lep_distreco/F");
178  fTree->Branch("lep_distrecox",&lep_distrecox,"lep_distrecox/F");
179  fTree->Branch("lep_distrecoy",&lep_distrecoy,"lep_distrecoy/F");
180  fTree->Branch("lep_distrecoz",&lep_distrecoz,"lep_distrecoz/F");
181  fTree->Branch("cos", &cos, "cos/F");
182  fTree->Branch("coslx", &coslx, "coslx/F");
183  fTree->Branch("cosly", &cosly, "cosly/F");
184  fTree->Branch("coslz", &coslz, "coslz/F");
185  fTree->Branch("coslrx", &coslrx, "coslrx/F");
186  fTree->Branch("coslry", &coslry, "coslry/F");
187  fTree->Branch("coslrz", &coslrz, "coslrz/F");
188  fTree->Branch("t0", &t0, "t0/F");
189  fTree->Branch("vtxrecomc", &vtxrecomc, "vtxrecomc/F");
190  fTree->Branch("vtxrecomcx", &vtxrecomcx, "vtxrecomcx/F");
191  fTree->Branch("vtxrecomcy", &vtxrecomcy, "vtxrecomcy/F");
192  fTree->Branch("vtxrecomcz", &vtxrecomcz, "vtxrecomcz/F");
193  fTree->Branch("vtxupstream", &vtxupstream, "vtxupstream/F");
194  fTree->Branch("vtxupstreamx", &vtxupstreamx, "vtxupstreamx/F");
195  fTree->Branch("vtxupstreamy", &vtxupstreamy, "vtxupstreamy/F");
196  fTree->Branch("vtxupstreamz", &vtxupstreamz, "vtxupstreamz/F");
197  fTree->Branch("ngamma", &ngamma, "ngamma/I");
198  fTree->Branch("convdist", &convdist, "convdist/F");
199 
200 
201  file.open("data.dat");
202  file1.open("data1.dat");
203 }
204 
205 void dunefd::IniSegReco::reconfigure(fhicl::ParameterSet const& pset)
206 {
207  fHitsModuleLabel = pset.get< std::string >("HitsModuleLabel");
208  fClusterModuleLabel = pset.get< std::string >("ClusterModuleLabel");
209  fTrackModuleLabel = pset.get< std::string >("TrackModuleLabel");
210  fVertexModuleLabel = pset.get< std::string >("VertexModuleLabel");
211  fGenieGenModuleLabel = pset.get< std::string >("GenieGenModuleLabel");
212 
213  fFidVolCut = pset.get< double >("FidVolCut");
214  return;
215 }
216 
217 void dunefd::IniSegReco::ResetVars()
218 {
219  ngamma = 0;
220  convdist = -10.0F;
221  pmatracks.clear();
222  lep_dist = -9999; //cm
223  lep_distx = -9999;
224  lep_disty = -9999;
225  lep_distz = -9999;
226  lep_distreco = -9999;
227  lep_distrecox = -9999;
228  lep_distrecoy = -9999;
229  lep_distrecoz = -9999;
230  lep_dedx = -9999;
231  t0 = -9999;
232  dx = -9999;
233  cos = -9999;
234  coslx = -9999; cosly = -9999; coslz = -9999;
235  coslrx = -9999; coslry = -9999; coslrz = -9999;
236  vtxrecomc = -9999;
237  vtxrecomcx = -9999;
238  vtxrecomcy = -9999;
239  vtxrecomcz = -9999;
240  vtxupstream = -9999;
241  vtxupstreamx = -9999;
242  vtxupstreamy = -9999;
243  vtxupstreamz = -9999;
244  return;
245 }
246 
247 void dunefd::IniSegReco::produce(art::Event& evt)
248 {
249  ResetVars();
250  const art::ServiceHandle<geo::Geometry> geom;
251 
252  run = evt.run();
253  subrun = evt.subRun();
254  event = evt.id().event();
255  isdata = evt.isRealData();
256 
257  std::unique_ptr< std::vector< recob::Track > > tracks(new std::vector< recob::Track >);
258  std::unique_ptr< std::vector< recob::SpacePoint > > allsp(new std::vector< recob::SpacePoint >);
259 
260  std::unique_ptr< art::Assns< recob::Track, recob::Hit > > trk2hit(new art::Assns< recob::Track, recob::Hit >);
261  std::unique_ptr< art::Assns< recob::Track, recob::SpacePoint > > trk2sp(new art::Assns< recob::Track, recob::SpacePoint >);
262  std::unique_ptr< art::Assns< recob::SpacePoint, recob::Hit > > sp2hit(new art::Assns< recob::SpacePoint, recob::Hit >);
263 
264  TVector3 primary(0, 0, 0);
265  bool isinside = false;
266  if (!isdata)
267  {
268  // * MC truth information
269  std::vector<art::Ptr<simb::MCTruth> > mclist;
270  auto mctruthListHandle = evt.getHandle< std::vector<simb::MCTruth> >(fGenieGenModuleLabel);
271  if (mctruthListHandle)
272  art::fill_ptr_vector(mclist, mctruthListHandle);
273 
274  if (mclist.size())
275  {
276  art::Ptr<simb::MCTruth> mctruth = mclist[0];
277  const TLorentzVector& pvtx = mctruth->GetNeutrino().Nu().Position();
278  primary = TVector3(pvtx.X(), pvtx.Y(), pvtx.Z());
279 
280  if (insideFidVol(pvtx))
281  {
282  isinside = true;
283  chargeParticlesatVtx(evt);
284 
285  art::ServiceHandle<cheat::ParticleInventoryService> pi_serv;
286  const sim::ParticleList& plist = pi_serv->ParticleList();
287 
288  bool photon = false;
289  for (sim::ParticleList::const_iterator ipar = plist.begin(); ipar != plist.end(); ++ipar)
290  {
291  simb::MCParticle* particle = ipar->second;
292  TLorentzVector mom = particle->Momentum();
293  TVector3 momvec(mom.Px(), mom.Py(), mom.Pz());
294 
295  if ((particle->PdgCode() == 22) && (momvec.Mag() > 0.030))
296  {
297  ngamma++; photon = true;
298  TLorentzVector conversion = particle->EndPosition();
299  TVector3 convec(conversion.X(), conversion.Y(), conversion.Z());
300  convdist = std::sqrt(pma::Dist2(primary, convec));
301  }
302  }
303  if (!photon) ngamma = -10;
304  }
305  }
306  }
307 
308  if (isinside)
309  {
310  auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(evt);
311  UseTracks(evt, detProp);
312  fTree->Fill();
313  }
314 
315  evt.put(std::move(tracks));
316  evt.put(std::move(allsp));
317  evt.put(std::move(trk2hit));
318  evt.put(std::move(trk2sp));
319  evt.put(std::move(sp2hit));
320 }
321 
322 /***********************************************************************/
323 
324 void dunefd::IniSegReco::UseClusters(art::Event const & evt,
325  detinfo::DetectorPropertiesData const& detProp)
326 {
327  if (!isdata)
328  {
329  // * MC truth information
330  std::vector<art::Ptr<simb::MCTruth> > mclist;
331  auto mctruthListHandle = evt.getHandle< std::vector<simb::MCTruth> >(fGenieGenModuleLabel);
332  if (mctruthListHandle)
333  art::fill_ptr_vector(mclist, mctruthListHandle);
334 
335  if (mclist.size())
336  {
337  art::Ptr<simb::MCTruth> mctruth = mclist[0];
338  collectCls(evt, detProp, mctruth);
339  }
340  }
341 }
342 
343 /***********************************************************************/
344 
345 void dunefd::IniSegReco::UseTracks(art::Event const & evt,
346  detinfo::DetectorPropertiesData const& detProp)
347 {
348  if (!isdata)
349  {
350  // * hits
351  std::vector<art::Ptr<recob::Hit> > hitlist;
352  auto hitListHandle = evt.getHandle< std::vector<recob::Hit> >(fHitsModuleLabel);
353  if (hitListHandle)
354  art::fill_ptr_vector(hitlist, hitListHandle);
355 
356  // * tracks
357  std::vector<art::Ptr<recob::Track> > tracklist;
358  auto trackListHandle = evt.getHandle< std::vector<recob::Track> >(fTrackModuleLabel);
359  if (trackListHandle)
360  art::fill_ptr_vector(tracklist, trackListHandle);
361 
362  // * vertices
363  std::vector<art::Ptr<recob::Vertex> > vtxlist;
364  auto vtxListHandle = evt.getHandle< std::vector<recob::Vertex> >(fVertexModuleLabel);
365  if (vtxListHandle)
366  art::fill_ptr_vector(vtxlist, vtxListHandle);
367 
368  // * monte carlo
369  std::vector<art::Ptr<simb::MCTruth> > mclist;
370  auto mctruthListHandle = evt.getHandle< std::vector<simb::MCTruth> >(fGenieGenModuleLabel);
371  if (mctruthListHandle)
372  art::fill_ptr_vector(mclist, mctruthListHandle);
373 
374  if (mclist.size())
375  {
376  art::Ptr<simb::MCTruth> mctruth = mclist[0];
377 
378  const simb::MCParticle& particle = mctruth->GetNeutrino().Nu();
379  const TLorentzVector& pvtx = particle.Position();
380  TVector3 primary(pvtx.X(), pvtx.Y(), pvtx.Z());
381 
382  // search for the closest reco vertex to mc primary
383  TVector3 closestvtx; TVector3 minzvtx;
384  if (vtxlist.size())
385  {
386  double xyz[3] = {0.0, 0.0, 0.0};
387  vtxlist[0]->XYZ(xyz);
388  float vxreco = xyz[0];
389  if (vxreco > 0) vxreco -= t0Corr(evt, detProp, pvtx);
390  else vxreco += t0Corr(evt, detProp, pvtx);
391  xyz[0] = vxreco;
392 
393  TVector3 vtxreco(xyz);
394  vtxreco.SetXYZ(xyz[0], xyz[1], xyz[2]);
395  closestvtx.SetXYZ(xyz[0], xyz[1], xyz[2]);
396 
397  double mindist2 = pma::Dist2(primary, vtxreco);
398  // loop over vertices to look for the closest to primary
399  for (size_t v = 1; v < vtxlist.size(); ++v)
400  {
401  vtxlist[v]->XYZ(xyz);
402  float temp = xyz[0];
403  if (temp > 0) temp -= t0Corr(evt, detProp, pvtx);
404  else temp += t0Corr(evt, detProp, pvtx);
405  xyz[0] = temp;
406 
407  vtxreco.SetXYZ(xyz[0], xyz[1], xyz[2]);
408  float dist2 = pma::Dist2(primary, vtxreco);
409  if (dist2 < mindist2)
410  {
411  closestvtx.SetXYZ(xyz[0], xyz[1], xyz[2]);
412  mindist2 = dist2;
413  }
414  }
415 
416  // loop over vertices to look the most upstream
417  double minz_xyz[3] = {0.0, 0.0, 0.0};
418  double minz = 9999;
419  for (size_t v = 0; v < vtxlist.size(); ++v)
420  {
421  vtxlist[v]->XYZ(minz_xyz);
422  float temp = minz_xyz[0];
423  if (temp > 0) temp -= t0Corr(evt, detProp, pvtx);
424  else temp += t0Corr(evt, detProp, pvtx);
425  minz_xyz[0] = temp;
426 
427  if (minz_xyz[2] < minz)
428  {
429  minz = minz_xyz[2];
430  minzvtx.SetXYZ(minz_xyz[0], minz_xyz[1], minz_xyz[2]);
431  }
432  }
433 
434  // loop over tracks
435  for (size_t t = 0; t < tracklist.size(); ++t)
436  {
437  if (!tracklist[t]->NumberTrajectoryPoints()) continue;
438  TVector3 pos_p = tracklist[t]->LocationAtPoint<TVector3>(0);
439 
440  float temp = pos_p.X();
441  if (temp > 0) temp -= t0Corr(evt, detProp, pvtx);
442  else temp += t0Corr(evt, detProp, pvtx);
443  pos_p.SetX(temp);
444 
445  float dist2 = pma::Dist2(primary, pos_p);
446  if (dist2 < mindist2)
447  {
448  closestvtx.SetXYZ(pos_p.X(), pos_p.Y(), pos_p.Z());
449  mindist2 = dist2;
450  }
451 
452  if (pos_p.Z() < minzvtx.Z())
453  {
454  minzvtx.SetXYZ(pos_p.X(), pos_p.Y(), pos_p.Z());
455  }
456 
457  TVector3 pos_end = tracklist[t]->LocationAtPoint<TVector3>(tracklist[t]->NumberTrajectoryPoints()-1);
458  temp = pos_end.X();
459  if (temp > 0) temp -= t0Corr(evt, detProp, pvtx);
460  else temp += t0Corr(evt, detProp, pvtx);
461  pos_end.SetX(temp);
462 
463  dist2 = pma::Dist2(primary, pos_end);
464  if (dist2 < mindist2)
465  {
466  closestvtx.SetXYZ(pos_end.X(), pos_end.Y(), pos_end.Z());
467  mindist2 = dist2;
468  }
469 
470  if (pos_end.Z() < minzvtx.Z())
471  {
472  minzvtx.SetXYZ(pos_end.X(), pos_end.Y(), pos_end.Z());
473  }
474  }
475 
476  vtxrecomc = std::sqrt(pma::Dist2(primary, closestvtx));
477  vtxrecomcx = primary.X() - closestvtx.X();
478  vtxrecomcy = primary.Y() - closestvtx.Y();
479  vtxrecomcz = primary.Z() - closestvtx.Z();
480  vtxupstream = std::sqrt(pma::Dist2(primary, minzvtx));
481  vtxupstreamx = primary.X() - minzvtx.X();
482  vtxupstreamy = primary.Y() - minzvtx.Y();
483  vtxupstreamz = primary.Z() - minzvtx.Z();
484  }
485 
486  // vertex region of event
487  // background events
488  if (insideFidVol(pvtx) && (abs(mctruth->GetNeutrino().Lepton().PdgCode()) != 11) && tracklist.size())
489  {
490  // mc
491  TVector3 mcvtx3d(pvtx.X(), pvtx.Y(), pvtx.Z());
492  std::vector< TVector3 > mcdir3d = findPhDir();
493 
494  float minlep_dist = 9999;
495  for (size_t v = 0; v < mcdir3d.size(); ++v)
496  {
497  // reco
498  IniSegAlg recoini(tracklist, mcvtx3d);
499  recoini.FeedwithMc(mcdir3d[v]);
500 
501  if (recoini.IsFound())
502  {
503  art::Ptr<recob::Track> recotrack = recoini.GetTrk();
504  if (!recotrack->NumberTrajectoryPoints()) continue;
505 
506  art::FindManyP< recob::Hit > fb(trackListHandle, evt, fTrackModuleLabel);
507  std::vector< art::Ptr<recob::Hit> > recoinihit = fb.at(recotrack.key());
508 
509  // use recob::Track functionality as much as possible
510  // const double setlength = 2.5; double length = 0.0; // cm
511 
512  TVector3 pos_p = recotrack->LocationAtPoint<TVector3>(0);
513  float px = pos_p.X();
514  if (px > 0) px -= t0Corr(evt, detProp, pvtx);
515  else px += t0Corr(evt, detProp, pvtx);
516  pos_p.SetX(px);
517 
518  float ldist = std::sqrt(pma::Dist2(primary, pos_p));
519  if (ldist < minlep_dist)
520  {
521  minlep_dist = ldist;
522  lep_dist = ldist;
523  lep_distx = primary.X() - pos_p.X();
524  lep_disty = primary.Y() - pos_p.Y();
525  lep_distz = primary.Z() - pos_p.Z();
526  if (mclist.size())
527  {
528  lep_distreco = std::sqrt(pma::Dist2(closestvtx, pos_p));
529  lep_distrecox = closestvtx.X() - pos_p.X();
530  lep_distrecoy = closestvtx.Y() - pos_p.Y();
531  lep_distrecoz = closestvtx.Z() - pos_p.Z();
532  }
533 
534  cos = recoini.GetCos();
535  coslrx = mcdir3d[v].X(); coslry = mcdir3d[v].Y(); coslrz = mcdir3d[v].Z();
536 
537  /*************************************************************/
538  /* WARNING */
539  /*************************************************************/
540  /* The dQdx information in recob::Track has been deprecated */
541  /* since 2016 and in 11/2018 the recob::Track interface was */
542  /* changed and DQdxAtPoint and NumberdQdx were removed. */
543  /* Therefore the code below is now commented out */
544  /* (note that it was most likely not functional anyways). */
545  /* For any issue please contact: larsoft-team@fnal.gov */
546  /*************************************************************/
547  /*
548  size_t fp = 0; bool hitcoll = false;
549  for (size_t p = 0; p < recotrack->NumberTrajectoryPoints(); ++p)
550  if (recotrack->DQdxAtPoint(p, geo::kZ) > 0)
551  {pos_p = recotrack->LocationAtPoint<TVector3>(p); fp = p; hitcoll = true; break;}
552 
553  // loop over trajectory point to get dQdx.
554  if (hitcoll)
555  for (size_t p = (fp+1); p < recotrack->NumberTrajectoryPoints(); ++p)
556  {
557  TVector3 pos = recotrack->LocationAtPoint<TVector3>(p);
558  length += std::sqrt(pma::Dist2(pos_p, pos));
559  pos_p = recotrack->LocationAtPoint<TVector3>(p);
560 
561  if (length > setlength) break;
562  dx = length;
563  double dqdx_p = recotrack->DQdxAtPoint(p, geo::kZ);
564  if (dqdx_p > 0) lep_dedx += recoinihit[p]->SummedADC();
565  }
566 
567  if (dx > 0.0) lep_dedx /= dx;
568  */
569  /*************************************************************/
570  }
571  }
572  }
573  }
574  else if (insideFidVol(pvtx) && (abs(mctruth->GetNeutrino().Lepton().PdgCode()) == 11) && tracklist.size()) // from nue vertex
575  {
576  // mc
577  TVector3 mcvtx3d(pvtx.X(), pvtx.Y(), pvtx.Z());
578  TVector3 mcdir3d = findElDir(mctruth);
579  coslx = mcdir3d.X(); cosly = mcdir3d.Y(); coslz = mcdir3d.Z();
580 
581  // reco
582  IniSegAlg recoini(tracklist, mcvtx3d);
583  recoini.FeedwithMc(mcdir3d); // mcvtx3d == primary
584 
585  if (recoini.IsFound())
586  {
587  art::Ptr<recob::Track> recotrack = recoini.GetTrk();
588  if (!recotrack->NumberTrajectoryPoints()) return;
589 
590  art::FindManyP< recob::Hit > fb(trackListHandle, evt, fTrackModuleLabel);
591  std::vector< art::Ptr<recob::Hit> > recoinihit = fb.at(recotrack.key());
592 
593  // use recob::Track functionality as much as possible
594  // const double setlength = 2.5; double length = 0.0; // cm
595 
596  TVector3 pos_p = recotrack->LocationAtPoint<TVector3>(0);
597  float px = pos_p.X();
598  if (px > 0) px -= t0Corr(evt,detProp, pvtx);
599  else px += t0Corr(evt, detProp, pvtx);
600  pos_p.SetX(px);
601 
602  lep_dist = std::sqrt(pma::Dist2(primary, pos_p));
603  lep_distx = primary.X() - pos_p.X();
604  lep_disty = primary.Y() - pos_p.Y();
605  lep_distz = primary.Z() - pos_p.Z();
606  if (mclist.size())
607  {
608  lep_distreco = std::sqrt(pma::Dist2(closestvtx, pos_p));
609  lep_distrecox = closestvtx.X() - pos_p.X();
610  lep_distrecoy = closestvtx.Y() - pos_p.Y();
611  lep_distrecoz = closestvtx.Z() - pos_p.Z();
612  }
613 
614  cos = recoini.GetCos();
615  coslrx = mcdir3d.X(); coslry = mcdir3d.Y(); coslrz = mcdir3d.Z();
616 
617  /*************************************************************/
618  /* WARNING */
619  /*************************************************************/
620  /* The dQdx information in recob::Track has been deprecated */
621  /* since 2016 and in 11/2018 the recob::Track interface was */
622  /* changed and DQdxAtPoint and NumberdQdx were removed. */
623  /* Therefore the code below is now commented out */
624  /* (note that it was most likely not functional anyways). */
625  /* For any issue please contact: larsoft-team@fnal.gov */
626  /*************************************************************/
627  /*
628  size_t fp = 0; bool hitcoll = false;
629  for (size_t p = 0; p < recotrack->NumberTrajectoryPoints(); ++p)
630  if (recotrack->DQdxAtPoint(p, geo::kZ) > 0) {pos_p = recotrack->LocationAtPoint<TVector3>(p); fp = p; hitcoll = true; break;}
631 
632  // loop over trajectory point to get dQdx.
633  if (hitcoll)
634  for (size_t p = (fp+1); p < recotrack->NumberTrajectoryPoints(); ++p)
635  {
636  TVector3 pos = recotrack->LocationAtPoint<TVector3>(p);
637  length += std::sqrt(pma::Dist2(pos_p, pos));
638  pos_p = recotrack->LocationAtPoint<TVector3>(p);
639 
640  if (length > setlength) break;
641  dx = length;
642  double dqdx_p = recotrack->DQdxAtPoint(p, geo::kZ);
643  if (dqdx_p > 0) lep_dedx += recoinihit[p]->SummedADC();
644  }
645  if (dx > 0.0) lep_dedx /= dx;
646  */
647  /*************************************************************/
648  }
649 
650  }
651  }
652  }
653 }
654 
655 /***********************************************************************/
656 
657 void dunefd::IniSegReco::collectCls(art::Event const & evt,
658  detinfo::DetectorPropertiesData const& detProp,
659  art::Ptr<simb::MCTruth> const mctruth)
660 {
661  art::ServiceHandle<geo::Geometry> geom;
662 
663  // * clusters
664  std::vector<art::Ptr<recob::Cluster> > clusterlist;
665  auto clusterListHandle = evt.getHandle< std::vector<recob::Cluster> >(fClusterModuleLabel);
666  if (clusterListHandle)
667  {
668  art::fill_ptr_vector(clusterlist, clusterListHandle);
669  art::FindManyP< recob::Hit > hc(clusterListHandle, evt, fClusterModuleLabel);
670 
671  for (size_t c = 0; c < geom->Ncryostats(); ++c)
672  {
673  const geo::CryostatGeo& cryo = geom->Cryostat(c);
674  for (size_t t = 0; t < cryo.NTPC(); ++t)
675  {
676  const geo::TPCGeo& tpc = cryo.TPC(t);
677  std::vector< std::vector< dunefd::Hit2D > > clinput;
678  for (size_t p = 0; p < tpc.Nplanes(); ++p)
679  {
680  std::map<size_t, std::vector< dunefd::Hit2D > > cls;
681  for (size_t i = 0; i < clusterlist.size(); ++i)
682  {
683  std::vector< art::Ptr<recob::Hit> > hitscl;
684  std::vector< Hit2D > hits2dcl;
685  hitscl = hc.at(i);
686 
687  bool found = false;
688  for (size_t h = 0; h < hitscl.size(); ++h)
689  {
690  size_t wire = hitscl[h]->WireID().Wire;
691  size_t plane = hitscl[h]->WireID().Plane;
692  size_t tpc = hitscl[h]->WireID().TPC;
693  size_t cryo = hitscl[h]->WireID().Cryostat;
694 
695  if ((plane == p) && (tpc == t) && (cryo == c))
696  {
697  found = true;
698  TVector2 point = pma::WireDriftToCm(detProp, wire, hitscl[h]->PeakTime(), plane, tpc, cryo);
699  Hit2D hit2d(point, hitscl[h].key());
700  hits2dcl.push_back(hit2d);
701  }
702  }
703 
704  if (found) cls[i] = hits2dcl;
705  }
706 
707  if (cls.size())
708  clinput.push_back(reselectCls(cls, mctruth, c, t, p));
709  }
710 
711  if (clinput.size() > 1)
712  {
713  const TLorentzVector& pvtx = mctruth->GetNeutrino().Nu().Position();
714  TVector3 primary(pvtx.X(), pvtx.Y(), pvtx.Z());
715  make3dseg(evt, detProp, clinput, primary);
716  }
717  }
718  }
719 
720  }
721 }
722 
723 /***********************************************************************/
724 
725 std::vector< dunefd::Hit2D > dunefd::IniSegReco::reselectCls(std::map<size_t, std::vector< dunefd::Hit2D > > const & cls, art::Ptr<simb::MCTruth> const mctruth, const size_t cryo, const size_t tpc, const size_t plane) const
726 {
727  std::vector< dunefd::Hit2D > cluster;
728  if (!cls.size()) return cluster;
729 
730  const simb::MCParticle& particle = mctruth->GetNeutrino().Nu();
731  const TLorentzVector& pvtx = particle.Position();
732 
733  if (insideFidVol(pvtx) && (abs(mctruth->GetNeutrino().Lepton().PdgCode()) == 11))
734  {
735  // mc
736  TVector3 mcvtx3d(pvtx.X(), pvtx.Y(), pvtx.Z());
737  TVector3 mcdir3d = findElDir(mctruth);
738  TVector2 mcvtx2d = getMCvtx2d(mcvtx3d, cryo, tpc, plane);
739  TVector2 mcdir2d = getMCdir2d(mcvtx3d, mcdir3d, cryo, tpc, plane);
740 
741  // reco: find the best clusters to proceed with segment reconstruction
742  IniSegAlg recoini(cls);
743  recoini.FeedwithMc(mcvtx2d, mcdir2d, mcdir3d);
744  cluster = recoini.GetCl();
745  }
746 
747  return cluster;
748 }
749 
750 /***********************************************************************/
751 
752 void dunefd::IniSegReco::make3dseg(
753  art::Event const & evt,
754  detinfo::DetectorPropertiesData const& detProp,
755  std::vector< std::vector< dunefd::Hit2D > > const & src,
756  TVector3 const & primary)
757 {
758  if (src.size() < 2) return;
759  if ((src[0].size() < 2) || (src[1].size() < 2)) return;
760 
761  std::vector<art::Ptr<recob::Hit> > hitlist;
762  auto hitListHandle = evt.getHandle< std::vector<recob::Hit> >(fHitsModuleLabel);
763  if (hitListHandle)
764  art::fill_ptr_vector(hitlist, hitListHandle);
765 
766  std::vector< art::Ptr<recob::Hit> > hitsel;
767 
768  for (size_t i = 0; i < src.size(); ++i)
769  for (size_t h = 0; h < src[i].size(); ++h)
770  hitsel.push_back(hitlist[src[i][h].GetKey()]);
771 
772  if (hitsel.size() > 5)
773  {
774  pma::Track3D* trk = fProjectionMatchingAlg.buildSegment(detProp, hitsel);
775  pmatracks.push_back(trk);
776  }
777 }
778 
779 /***********************************************************************/
780 
781 TVector2 dunefd::IniSegReco::getMCdir2d(TVector3 const & mcvtx3d, TVector3 const & mcdir3d,
782  const size_t cryo, const size_t tpc, const size_t plane) const
783 {
784  TVector3 shift3d = mcvtx3d + mcdir3d;
785  TVector2 shift2d = pma::GetProjectionToPlane(shift3d, plane, tpc, cryo) - getMCvtx2d(mcvtx3d, cryo, tpc, plane);
786  TVector2 mcdir2d = shift2d * (1.0 / shift2d.Mod());
787 
788  return mcdir2d;
789 }
790 
791 /***********************************************************************/
792 
793 TVector2 dunefd::IniSegReco::getMCvtx2d(TVector3 const & mcvtx3d, const size_t cryo, const size_t tpc, const size_t plane) const
794 {
795  TVector2 mcvtx2d = pma::GetProjectionToPlane(mcvtx3d, plane, tpc, cryo);
796 
797  return mcvtx2d;
798 }
799 
800 /***********************************************************************/
801 
802 TVector3 dunefd::IniSegReco::findElDir(art::Ptr<simb::MCTruth> const mctruth) const
803 {
804  TVector3 dir(0, 0, 0);
805  const simb::MCParticle& lepton = mctruth->GetNeutrino().Lepton();
806 
807  if (abs(lepton.PdgCode()) == 11)
808  {
809  TLorentzVector mom = lepton.Momentum();
810  TVector3 momvec3(mom.Px(), mom.Py(), mom.Pz());
811  dir = momvec3 * (1 / momvec3.Mag());
812  }
813 
814  return dir;
815 }
816 
817 /***********************************************************************/
818 
819 void dunefd::IniSegReco::chargeParticlesatVtx(art::Event const & evt)
820 {
821  std::cout << " chargeParticleatVtx " << std::endl;
822  art::ServiceHandle<cheat::ParticleInventoryService> pi_serv;
823  const sim::ParticleList& plist = pi_serv->ParticleList();
824 
825  TVector3 pospri; bool primary = false;
826  for (sim::ParticleList::const_iterator ipar = plist.begin(); ipar != plist.end(); ++ipar)
827  {
828  const simb::MCParticle* particle = ipar->second;
829 
830  if (particle->Process() != "primary") continue;
831 
832  pospri.SetXYZ(particle->Position(0).X(), particle->Position(0).Y(), particle->Position(0).Z());
833  primary = true; break;
834  }
835 
836  if (!primary) return;
837 
838  size_t ninteractions = 0;
839  for (sim::ParticleList::const_iterator ipar = plist.begin(); ipar != plist.end(); ++ipar)
840  {
841  std::vector< double > temp;
842  const simb::MCParticle* particle = ipar->second;
843 
844  std::cout << " pdg " << particle->PdgCode() << " " << particle->P(0) << std::endl;
845  std::cout << " position " << particle->Position(0).X() << ", " << particle->Position(0).Y() << ", " << particle->Position(0).Z() << std::endl;
846 
847  const TLorentzVector& startpos = particle->Position(0);
848  TVector3 start(startpos.X(), startpos.Y(), startpos.Z());
849  const TLorentzVector& endpos = particle->EndPosition();
850  TVector3 stop(endpos.X(), endpos.Y(), endpos.Z());
851 
852  double diststartpri = std::sqrt(pma::Dist2(start, pospri));
853  double diststoppri = std::sqrt(pma::Dist2(stop, pospri));
854  if ((diststartpri > 0.1) && (diststartpri < 1.0))
855  {
856  //if (diststoppri < 3.0)
857  //{
858  double ek = particle->E(0) - particle->Mass();
859  file << evt.run() << " " << evt.id().event() << " " << particle->PdgCode() << " " << diststartpri << " " << diststoppri << " " << particle->P(0) << " " << particle->Mass() << " " << ek << std::endl;
860 
861  if (ek > 0.05) ninteractions++;
862  //}
863  }
864  }
865 
866  file1 << evt.run() << " " << evt.id().event() << " " << ninteractions << std::endl;
867 }
868 
869 /***********************************************************************/
870 
871 std::vector< TVector3 > dunefd::IniSegReco::findPhDir() const
872 {
873  std::vector< TVector3 > phdirs;
874  art::ServiceHandle<cheat::ParticleInventoryService> pi_serv;
875  const sim::ParticleList& plist = pi_serv->ParticleList();
876 
877  for (sim::ParticleList::const_iterator ipar = plist.begin(); ipar != plist.end(); ++ipar)
878  {
879  const simb::MCParticle* particle = ipar->second;
880  if (particle->Process() != "primary") continue;
881 
882  TVector3 dir(0, 0, 0);
883  if (particle->PdgCode() == 111)
884  {
885  if (particle->NumberDaughters() != 2) continue;
886 
887  const simb::MCParticle* daughter1 = pi_serv->TrackIdToParticle_P(particle->Daughter(0));
888  if (daughter1->PdgCode() != 22) continue;
889 
890  const simb::MCParticle* daughter2 = pi_serv->TrackIdToParticle_P(particle->Daughter(1));
891  if (daughter2->PdgCode() != 22) continue;
892 
893  if (daughter1->EndProcess() == "conv")
894  {
895  TLorentzVector mom = pi_serv->TrackIdToParticle_P(particle->Daughter(0))->Momentum();
896  TVector3 momvec3(mom.Px(), mom.Py(), mom.Pz());
897  dir = momvec3 * (1 / momvec3.Mag());
898  phdirs.push_back(dir);
899  }
900 
901  if (daughter2->EndProcess() == "conv")
902  {
903  TLorentzVector mom = pi_serv->TrackIdToParticle_P(particle->Daughter(1))->Momentum();
904  TVector3 momvec3(mom.Px(), mom.Py(), mom.Pz());
905  dir = momvec3 * (1 / momvec3.Mag());
906  phdirs.push_back(dir);
907  }
908  }
909  else if (particle->PdgCode() == 22)
910  {
911  TLorentzVector mom = particle->Momentum();
912  TVector3 momvec3(mom.Px(), mom.Py(), mom.Pz());
913  dir = momvec3 * (1 / momvec3.Mag());
914  phdirs.push_back(dir);
915  }
916  }
917 
918  return phdirs;
919 }
920 
921 /***********************************************************************/
922 
923 std::vector< TVector3 > dunefd::IniSegReco::findDirs(art::Ptr<simb::MCTruth> const mctruth, int pdg) const
924 {
925  std::vector< TVector3 > dirs;
926 
927  art::ServiceHandle<cheat::ParticleInventoryService> pi_serv;
928  const sim::ParticleList& plist = pi_serv->ParticleList();
929 
930  for (sim::ParticleList::const_iterator ipar = plist.begin(); ipar != plist.end(); ++ipar)
931  {
932  simb::MCParticle* particle = ipar->second;
933  TLorentzVector mom = particle->Momentum();
934  TVector3 momvec(mom.Px(), mom.Py(), mom.Pz());
935 
936  if ((particle->PdgCode() == pdg) && (momvec.Mag() > 0.030))
937  {
938  TLorentzVector momconv = particle->EndMomentum();
939  TVector3 momconvec3(momconv.Px(), momconv.Py(), momconv.Pz());
940  TVector3 dir = momconvec3 * (1 / momconvec3.Mag());
941  dirs.push_back(dir);
942  }
943  }
944 
945  return dirs;
946 }
947 
948 /***********************************************************************/
949 
950 float dunefd::IniSegReco::t0Corr(art::Event const & evt,
951  detinfo::DetectorPropertiesData const& detProp,
952  TLorentzVector const & pvtx)
953 {
954  float corrt0x = 0.0F;
955 
956  art::ServiceHandle<geo::Geometry> geom;
957 
958  // * MC truth information
959  std::vector<art::Ptr<simb::MCTruth> > mclist;
960  auto mctruthListHandle = evt.getHandle< std::vector<simb::MCTruth> >(fGenieGenModuleLabel);
961  if (mctruthListHandle)
962  art::fill_ptr_vector(mclist, mctruthListHandle);
963 
964  double vtx[3] = {pvtx.X(), pvtx.Y(), pvtx.Z()};
965 
966  geo::TPCID idtpc = geom->FindTPCAtPosition(vtx);
967 
968  if (geom->HasTPC(idtpc))
969  if (mclist.size())
970  {
971  art::Ptr<simb::MCTruth> mctruth = mclist[0];
972  if (mctruth->NParticles())
973  {
974  simb::MCParticle particle = mctruth->GetParticle(0);
975  t0 = particle.T(); // ns
976  corrt0x = t0 * 1.e-3 * detProp.DriftVelocity();
977  }
978  }
979 
980 
981  return corrt0x;
982 }
983 
984 /***********************************************************************/
985 
986 bool dunefd::IniSegReco::insideFidVol(TLorentzVector const & pvtx) const
987 {
988  art::ServiceHandle<geo::Geometry> geom;
989  double vtx[3] = {pvtx.X(), pvtx.Y(), pvtx.Z()};
990  bool inside = false;
991 
992  geo::TPCID idtpc = geom->FindTPCAtPosition(vtx);
993 
994  if (geom->HasTPC(idtpc))
995  {
996  const geo::TPCGeo& tpcgeo = geom->GetElement(idtpc);
997  double minx = tpcgeo.MinX(); double maxx = tpcgeo.MaxX();
998  double miny = tpcgeo.MinY(); double maxy = tpcgeo.MaxY();
999  double minz = tpcgeo.MinZ(); double maxz = tpcgeo.MaxZ();
1000 
1001  /*for (size_t c = 0; c < geom->Ncryostats(); c++)
1002  {
1003  const geo::CryostatGeo& cryostat = geom->Cryostat(c);
1004  for (size_t t = 0; t < cryostat.NTPC(); t++)
1005  {
1006  const geo::TPCGeo& tpcg = cryostat.TPC(t);
1007  if (tpcg.MinX() < minx) minx = tpcg.MinX();
1008  if (tpcg.MaxX() > maxx) maxx = tpcg.MaxX();
1009  if (tpcg.MinY() < miny) miny = tpcg.MinY();
1010  if (tpcg.MaxY() > maxy) maxy = tpcg.MaxY();
1011  if (tpcg.MinZ() < minz) minz = tpcg.MinZ();
1012  if (tpcg.MaxZ() > maxz) maxz = tpcg.MaxZ();
1013  }
1014  }*/
1015 
1016  //x
1017  double dista = fabs(minx - pvtx.X());
1018  double distb = fabs(pvtx.X() - maxx);
1019 ;
1020  if ((pvtx.X() > minx) && (pvtx.X() < maxx) &&
1021  (dista > fFidVolCut) && (distb > fFidVolCut))
1022  {
1023  inside = true;
1024  }
1025  else { inside = false; }
1026  //y
1027  dista = fabs(maxy - pvtx.Y());
1028  distb = fabs(pvtx.Y() - miny);
1029  if (inside && (pvtx.Y() > miny) && (pvtx.Y() < maxy) &&
1030  (dista > fFidVolCut) && (distb > fFidVolCut)) inside = true;
1031  else inside = false;
1032 
1033  //z
1034  dista = fabs(maxz - pvtx.Z());
1035  distb = fabs(pvtx.Z() - minz);
1036  if (inside && (pvtx.Z() > minz) && (pvtx.Z() < maxz) &&
1037  (dista > fFidVolCut) && (distb > fFidVolCut)) inside = true;
1038  else inside = false;
1039  }
1040 
1041  return inside;
1042 }
1043 
1044 /***********************************************************************/
1045 
1046 recob::Track dunefd::IniSegReco::convertFrom(pma::Track3D const & src)
1047 {
1048  std::vector< TVector3 > xyz, dircos;
1049 
1050  for (size_t i = 0; i < src.size(); ++i)
1051  {
1052  xyz.push_back(src[i]->Point3D());
1053 
1054  if (i < src.size() - 1)
1055  {
1056  TVector3 dc(src[i + 1]->Point3D());
1057  dc -= src[i]->Point3D();
1058  dc *= 1.0 / dc.Mag();
1059  dircos.push_back(dc);
1060  }
1061  else dircos.push_back(dircos.back());
1062  }
1063 
1064  if (xyz.size() != dircos.size())
1065  {
1066  mf::LogError("IniSegReco") << "pma::Track3D to recob::Track conversion problem.";
1067  }
1068  return recob::Track(recob::TrackTrajectory(recob::tracking::convertCollToPoint(xyz),
1069  recob::tracking::convertCollToVector(dircos),
1070  recob::Track::Flags_t(xyz.size()), false),
1071  0, -1., 0, recob::tracking::SMatrixSym55(), recob::tracking::SMatrixSym55(), trkindex);
1072 }
1073 
1074 /***********************************************************************/
1075 
1076 DEFINE_ART_MODULE(dunefd::IniSegReco)
simb::MCParticle::E
double E(const int i=0) const
Definition: MCParticle.h:233
t0
code to link reconstructed objects back to the MC truth information
Definition: DirectHitParticleAssns_tool.cc:12
simb::MCParticle::Position
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:219
simb::MCParticle::PdgCode
int PdgCode() const
Definition: MCParticle.h:212
simb::MCTruth::GetNeutrino
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:77
geo::GeometryCore::GetElement
CryostatGeo const & GetElement(geo::CryostatID const &cryoid) const
Definition: GeometryCore.h:2033
gar::rec::alg::cluster
void cluster(In first, In last, Out result, Pred *pred)
Definition: NNClusters.h:41
simb::MCParticle::EndPosition
const TLorentzVector & EndPosition() const
Definition: MCParticle.h:225
PmaTrack3D.h
Implementation of the Projection Matching Algorithm.
dunefd::IniSegReco::operator=
IniSegReco & operator=(IniSegReco const &)=delete
pma::Dist2
double Dist2(const TVector2 &v1, const TVector2 &v2)
Definition: Utilities.cxx:37
art::DataViewImpl::getHandle
Handle< PROD > getHandle(SelectorBase const &) const
Definition: DataViewImpl.h:382
string
std::string string
Definition: nybbler.cc:12
recob::tracking::SMatrixSym55
ROOT::Math::SMatrix< Double32_t, 5, 5, ROOT::Math::MatRepSym< Double32_t, 5 > > SMatrixSym55
Definition: TrackingTypes.h:85
cheat::ParticleInventoryService::TrackIdToParticle_P
const simb::MCParticle * TrackIdToParticle_P(int id) const
Definition: ParticleInventoryService_service.cc:128
dunefd::IniSegReco::make3dseg
void make3dseg(art::Event const &evt, detinfo::DetectorPropertiesData const &detProp, std::vector< std::vector< Hit2D > > const &src, TVector3 const &primary)
Definition: IniSegReco_module.cc:752
recob::tracking::convertCollToVector
std::vector< Vector_t > convertCollToVector(std::vector< Vector > const &coll)
Definition: TrackingTypes.h:58
geo::TPCGeo::Nplanes
unsigned int Nplanes() const
Number of planes in this tpc.
Definition: TPCGeo.h:165
simb::MCNeutrino::Nu
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
recob::Track::LocationAtPoint
Point_t const & LocationAtPoint(size_t i) const
Definition: Track.h:126
art::EDProducer::EDProducer
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.h:20
lar::fhicl::geo::Point3D
::fhicl::TupleAs< Point(::geo::Length_t,::geo::Length_t,::geo::Length_t)> Point3D
Atom object for reading a 3D point or vector (centimeters).
Definition: geo_vectors_fhicl.h:38
simb::MCParticle::Mass
double Mass() const
Definition: MCParticle.h:239
geo::BoxBoundedGeo::MinX
double MinX() const
Returns the world x coordinate of the start of the box.
Definition: BoxBoundedGeo.h:88
geo::TPCGeo
Geometry information for a single TPC.
Definition: TPCGeo.h:38
dunefd::IniSegReco::convertFrom
recob::Track convertFrom(pma::Track3D const &src)
Definition: IniSegReco_module.cc:1046
vector
struct vector vector
recob::Track::Flags_t
TrackTrajectory::Flags_t Flags_t
Definition: Track.h:68
recob::Track::NumberTrajectoryPoints
size_t NumberTrajectoryPoints() const
Various functions related to the presence and the number of (valid) points.
Definition: Track.h:102
pma::WireDriftToCm
TVector2 WireDriftToCm(detinfo::DetectorPropertiesData const &detProp, unsigned int wire, float drift, unsigned int plane, unsigned int tpc, unsigned int cryo)
Definition: Utilities.cxx:294
geo::BoxBoundedGeo::MaxX
double MaxX() const
Returns the world x coordinate of the end of the box.
Definition: BoxBoundedGeo.h:91
run_cnn_1class.tracks
tracks
Definition: run_cnn_1class.py:54
dunefd::IniSegReco::getMCvtx2d
TVector2 getMCvtx2d(TVector3 const &mcvtx3d, const size_t cryo, const size_t tpc, const size_t plane) const
Definition: IniSegReco_module.cc:793
dunefd::IniSegReco::chargeParticlesatVtx
void chargeParticlesatVtx(art::Event const &evt)
Definition: IniSegReco_module.cc:819
const_iterator
intermediate_table::const_iterator const_iterator
Definition: intermediate_table.cc:28
simb::MCTruth::NParticles
int NParticles() const
Definition: MCTruth.h:75
dunefd::IniSegReco::IniSegReco
IniSegReco(fhicl::ParameterSet const &p)
Definition: IniSegReco_module.cc:151
reco_momentum_tuples.px
px
Definition: reco_momentum_tuples.py:218
dir
string dir
Definition: INukeNucleonCorr.cxx:61
mf::LogError
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
Definition: MessageLogger.h:211
geo::CryostatGeo
Geometry information for a single cryostat.
Definition: CryostatGeo.h:43
geo::GeometryCore::Ncryostats
unsigned int Ncryostats() const
Returns the number of cryostats in the detector.
Definition: GeometryCore.h:1994
simb::MCParticle::Process
std::string Process() const
Definition: MCParticle.h:215
dunefd::IniSegReco::UseClusters
void UseClusters(art::Event const &evt, detinfo::DetectorPropertiesData const &detProp)
Definition: IniSegReco_module.cc:324
simb::MCParticle::NumberDaughters
int NumberDaughters() const
Definition: MCParticle.h:217
Geometry.h
art framework interface to geometry description
simb::MCParticle::Daughter
int Daughter(const int i) const
Definition: MCParticle.cxx:112
dunefd::IniSegReco::pmatracks
std::vector< pma::Track3D * > pmatracks
Definition: IniSegReco_module.cc:102
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
art::DataViewImpl::isRealData
bool isRealData() const
Definition: DataViewImpl.cc:113
reco_momentum_tuples.h
h
Definition: reco_momentum_tuples.py:66
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
dunefd::IniSegReco::findDirs
std::vector< TVector3 > findDirs(art::Ptr< simb::MCTruth > const mctruth, int pdg) const
Definition: IniSegReco_module.cc:923
e
const double e
Definition: gUpMuFluxGen.cxx:165
simb::MCNeutrino::Lepton
const simb::MCParticle & Lepton() const
Definition: MCNeutrino.h:147
dunefd::IniSegReco::produce
void produce(art::Event &e) override
Definition: IniSegReco_module.cc:247
DEFINE_ART_MODULE
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:67
simb::MCParticle::EndProcess
std::string EndProcess() const
Definition: MCParticle.h:216
wirecell.dfp.graph.key
def key(type, name=None)
Definition: graph.py:13
recob::TrackTrajectory
A trajectory in space reconstructed from hits.
Definition: TrackTrajectory.h:62
dunefd::IniSegAlg::FeedwithMc
void FeedwithMc(TVector2 const &vtx, TVector2 const &dir, TVector3 const &dir3d)
Definition: IniSegAlg.cxx:47
dunefd::IniSegReco::fGenieGenModuleLabel
std::string fGenieGenModuleLabel
Definition: IniSegReco_module.cc:141
submit_ppfxjobs.temp
temp
Definition: submit_ppfxjobs.py:110
reco_momentum_tuples.t
t
Definition: reco_momentum_tuples.py:25
wirecell.gen.depos.move
def move(depos, offset)
Definition: depos.py:107
dunefd::IniSegReco::insideFidVol
bool insideFidVol(TLorentzVector const &pvtx) const
Definition: IniSegReco_module.cc:986
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
simb::MCParticle::T
double T(const int i=0) const
Definition: MCParticle.h:224
dunefd::IniSegReco::beginJob
void beginJob() override
Definition: IniSegReco_module.cc:163
geo::BoxBoundedGeo::MinZ
double MinZ() const
Returns the world z coordinate of the start of the box.
Definition: BoxBoundedGeo.h:118
test.p
p
Definition: test.py:223
geo::CryostatGeo::NTPC
unsigned int NTPC() const
Number of TPCs in this cryostat.
Definition: CryostatGeo.h:181
geo::GeometryCore::Cryostat
CryostatGeo const & Cryostat(geo::CryostatID const &cryoid) const
Returns the specified cryostat.
Definition: GeometryCore.cxx:321
dunefd::IniSegReco::collectCls
void collectCls(art::Event const &evt, detinfo::DetectorPropertiesData const &detProp, art::Ptr< simb::MCTruth > const mctruth)
Definition: IniSegReco_module.cc:657
art::DataViewImpl::subRun
SubRunNumber_t subRun() const
Definition: DataViewImpl.cc:78
art::DataViewImpl::put
ProductID put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: DataViewImpl.h:686
recob::tracking::convertCollToPoint
std::vector< Point_t > convertCollToPoint(std::vector< Point > const &coll)
Definition: TrackingTypes.h:55
make_THn_beam_input.Momentum
Momentum
Definition: make_THn_beam_input.py:70
art::DataViewImpl::run
RunNumber_t run() const
Definition: DataViewImpl.cc:71
detinfo::DetectorPropertiesData::DriftVelocity
double DriftVelocity(double efield=0., double temperature=0.) const
cm/us
Definition: DetectorPropertiesData.cc:23
pma::ProjectionMatchingAlg::buildSegment
pma::Track3D * buildSegment(const detinfo::DetectorPropertiesData &clockData, const std::vector< art::Ptr< recob::Hit >> &hits_1, const std::vector< art::Ptr< recob::Hit >> &hits_2={}) const
Definition: ProjectionMatchingAlg.cxx:899
geo::TPCID
The data type to uniquely identify a TPC.
Definition: geo_types.h:386
dunefd::IniSegReco::findElDir
TVector3 findElDir(art::Ptr< simb::MCTruth > const mctruth) const
Definition: IniSegReco_module.cc:802
dunefd::IniSegReco::t0Corr
float t0Corr(art::Event const &evt, detinfo::DetectorPropertiesData const &detProp, TLorentzVector const &pvtx)
Definition: IniSegReco_module.cc:950
dunefd::IniSegReco::reconfigure
void reconfigure(fhicl::ParameterSet const &p)
Definition: IniSegReco_module.cc:205
geo::BoxBoundedGeo::MaxY
double MaxY() const
Returns the world y coordinate of the end of the box.
Definition: BoxBoundedGeo.h:106
cheat::ParticleInventoryService::ParticleList
const sim::ParticleList & ParticleList() const
Definition: ParticleInventoryService_service.cc:114
dunefd::IniSegReco::fHitsModuleLabel
std::string fHitsModuleLabel
Definition: IniSegReco_module.cc:137
dunefd::IniSegReco::fVertexModuleLabel
std::string fVertexModuleLabel
Definition: IniSegReco_module.cc:139
simb::MCTruth::GetParticle
const simb::MCParticle & GetParticle(int i) const
Definition: MCTruth.h:76
dunefd::IniSegAlg::IsFound
const bool IsFound()
Definition: IniSegAlg.h:54
dunefd::IniSegReco::fClusterModuleLabel
std::string fClusterModuleLabel
Definition: IniSegReco_module.cc:138
geo::GeometryCore::HasTPC
bool HasTPC(geo::TPCID const &tpcid) const
Returns whether we have the specified TPC.
Definition: GeometryCore.h:2397
Utilities.h
Implementation of the Projection Matching Algorithm.
Hit.h
Declaration of signal hit object.
dunefd::IniSegReco::getMCdir2d
TVector2 getMCdir2d(TVector3 const &mcvtx3d, TVector3 const &mcdir3d, const size_t cryo, const size_t tpc, const size_t plane) const
Definition: IniSegReco_module.cc:781
geo::CryostatGeo::TPC
const TPCGeo & TPC(unsigned int itpc) const
Return the itpc&#39;th TPC in the cryostat.
Definition: CryostatGeo.cxx:93
geo::BoxBoundedGeo::MaxZ
double MaxZ() const
Returns the world z coordinate of the end of the box.
Definition: BoxBoundedGeo.h:121
dunefd::IniSegAlg::GetCl
std::vector< dunefd::Hit2D > const & GetCl() const
Definition: IniSegAlg.h:52
pma::GetProjectionToPlane
TVector2 GetProjectionToPlane(const TVector3 &p, unsigned int plane, unsigned int tpc, unsigned int cryo)
Definition: Utilities.cxx:270
simb::MCParticle::Momentum
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:220
dunefd::IniSegReco::UseTracks
void UseTracks(art::Event const &evt, detinfo::DetectorPropertiesData const &detProp)
Definition: IniSegReco_module.cc:345
Track.h
Provides recob::Track data product.
run_cnn_1class.pdg
pdg
Definition: run_cnn_1class.py:54
art::EventID::event
EventNumber_t event() const
Definition: EventID.h:116
dunefd::IniSegReco::fTrackModuleLabel
std::string fTrackModuleLabel
Definition: IniSegReco_module.cc:140
proxy::Track
TrackCollectionProxyElement< TrackCollProxy > Track
Proxy to an element of a proxy collection of recob::Track objects.
Definition: Track.h:1036
dunefd::IniSegReco::reselectCls
std::vector< dunefd::Hit2D > reselectCls(std::map< size_t, std::vector< dunefd::Hit2D > > const &cls, art::Ptr< simb::MCTruth > const mctruth, const size_t cryo, const size_t tpc, const size_t plane) const
Definition: IniSegReco_module.cc:725
tca::evt
TCEvent evt
Definition: DataStructs.cxx:7
art::fill_ptr_vector
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:297
pma::Track3D::size
size_t size() const
Definition: PmaTrack3D.h:108
dunefd::IniSegReco::fProjectionMatchingAlg
pma::ProjectionMatchingAlg fProjectionMatchingAlg
Definition: IniSegReco_module.cc:144
dunefd::IniSegAlg::GetCos
float const & GetCos() const
Definition: IniSegAlg.h:58
geo::BoxBoundedGeo::MinY
double MinY() const
Returns the world y coordinate of the start of the box.
Definition: BoxBoundedGeo.h:103
simb::MCParticle::EndMomentum
const TLorentzVector & EndMomentum() const
Definition: MCParticle.h:240
art::Event::id
EventID id() const
Definition: Event.cc:34
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
dunefd::IniSegAlg::GetTrk
art::Ptr< recob::Track > const & GetTrk() const
Definition: IniSegAlg.h:55
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
event
Event finding and building.
Definition: EventCheater_module.cc:32
dunefd::IniSegReco::findPhDir
std::vector< TVector3 > findPhDir() const
Definition: IniSegReco_module.cc:871
Generated by   doxygen 1.8.11