LArSoft  v10_04_05
Liquid Argon Software toolkit - https://larsoft.org/
MultiEMShowers_module.cc
Go to the documentation of this file.
1 // Class: MultiEMShowers
3 // Module Type: analyzer
4 // File: MultiEMShowers_module.cc
5 // Author: dorota.stefan@cern.ch robert.sulej@cern.ch
7 
13 #include "art_root_io/TFileService.h"
16 #include "fhiclcpp/ParameterSet.h"
18 
29 
30 #include <algorithm>
31 #include <limits>
32 #include <memory>
33 
34 // ROOT includes
35 #include "TLorentzVector.h"
36 #include "TTree.h"
37 
38 namespace ems {
39  class MCinfo;
40  class MultiEMShowers;
41 }
42 
43 class ems::MCinfo {
44 public:
45  MCinfo();
46  void Info();
47  void Findtpcborders();
48 
49  int GetNgammas() const { return fNgammas; }
50 
51  double GetMompi0() const { return fMompi0; }
52  double GetMomGamma1() const { return fGammamom1; }
53  double GetMomGamma2() const { return fGammamom2; }
54 
55  double GetCosine() { return fCosine; }
56 
57  TVector3 const& GetPrimary() const { return fPrimary; }
58  TVector3 const& GetPospi0() const { return fPi0pos; }
59  TVector3 const& GetPosgamma1() const { return fConvgamma1; }
60  TVector3 const& GetPosgamma2() const { return fConvgamma2; }
61 
62  TVector3 const& GetDirgamma1() const { return fDirgamma1; }
63  TVector3 const& GetDirgamma2() const { return fDirgamma2; }
64 
65  bool IsInside1() const { return fInside1; }
66  bool IsInside2() const { return fInside2; }
67 
68  bool IsCompton() const { return fCompton; }
69 
70 private:
71  bool insideFidVol(const TLorentzVector& pvtx) const;
72  double fMinx{std::numeric_limits<double>::max()};
73  double fMaxx{std::numeric_limits<double>::min()};
74  double fMiny{std::numeric_limits<double>::max()};
75  double fMaxy{std::numeric_limits<double>::min()};
76  double fMinz{std::numeric_limits<double>::max()};
77  double fMaxz{std::numeric_limits<double>::min()};
78 
79  double fFidVolCut;
80 
81  int fNgammas;
82 
83  double fMompi0;
84  double fGammamom1;
85  bool fInside1;
86  double fGammamom2;
87  bool fInside2;
88 
89  double fCosine;
90 
91  bool fCompton;
92 
93  TVector3 fPrimary;
94  TVector3 fPi0pos;
95  TVector3 fConvgamma1;
96  TVector3 fConvgamma2;
97  TVector3 fDirgamma1;
98  TVector3 fDirgamma2;
99 };
100 
102 {
103  Info();
104  Findtpcborders();
105 }
106 
108 {
110  for (const geo::TPCGeo& tpcg : geom->Iterate<geo::TPCGeo>()) {
111  fMinx = std::min(fMinx, tpcg.MinX());
112  fMaxx = std::max(fMaxx, tpcg.MaxX());
113  fMiny = std::min(fMiny, tpcg.MinY());
114  fMaxy = std::max(fMaxy, tpcg.MaxY());
115  fMinz = std::min(fMinz, tpcg.MinZ());
116  fMaxz = std::max(fMaxz, tpcg.MaxZ());
117  }
118 }
119 
121 {
122  fMompi0 = 0.0;
123  fPi0pos.SetXYZ(0, 0, 0);
124  fNgammas = 0;
125  fCosine = 0.0;
126  fInside1 = false;
127  fInside2 = false;
128  fCompton = false;
129 
130  fGammamom1 = 0.0;
131  fGammamom2 = 0.0;
132  fConvgamma1.SetXYZ(0, 0, 0);
133  fConvgamma2.SetXYZ(0, 0, 0);
134  fDirgamma1.SetXYZ(0, 0, 0);
135  fDirgamma2.SetXYZ(0, 0, 0);
136 
138  const sim::ParticleList& plist = pi_serv->ParticleList();
139  for (sim::ParticleList::const_iterator ipar = plist.begin(); ipar != plist.end(); ++ipar) {
140  const simb::MCParticle* particle = ipar->second;
141 
142  if (particle->Process() != "primary") continue;
143 
144  TLorentzVector posvec = particle->Position();
145  TVector3 pose(posvec.X(), posvec.Y(), posvec.Z());
146  fPrimary = pose;
147 
148  if (particle->PdgCode() == 111) {
149  fMompi0 = particle->P();
150 
151  TLorentzVector posvec3 = particle->Position();
152  TVector3 pospi0(posvec3.X(), posvec3.Y(), posvec3.Z());
153  fPi0pos = pospi0;
154 
155  if (particle->NumberDaughters() != 2) continue;
156 
157  const simb::MCParticle* daughter1 = pi_serv->TrackIdToParticle_P(particle->Daughter(0));
158  if (daughter1->PdgCode() != 22) continue;
159 
160  const simb::MCParticle* daughter2 = pi_serv->TrackIdToParticle_P(particle->Daughter(1));
161  if (daughter2->PdgCode() != 22) continue;
162 
163  fNgammas = particle->NumberDaughters();
164  TLorentzVector mom1 = pi_serv->TrackIdToParticle_P(particle->Daughter(0))->Momentum();
165  TLorentzVector mom2 = pi_serv->TrackIdToParticle_P(particle->Daughter(1))->Momentum();
166 
167  // compton process
168  if (daughter1->EndProcess() == "phot") fCompton = true;
169  if (daughter2->EndProcess() == "phot") fCompton = true;
170 
171  TVector3 mom1vec3(mom1.Px(), mom1.Py(), mom1.Pz());
172  fGammamom1 = pi_serv->TrackIdToParticle_P(particle->Daughter(0))->P();
173  TVector3 mom2vec3(mom2.Px(), mom2.Py(), mom2.Pz());
174  fGammamom2 = pi_serv->TrackIdToParticle_P(particle->Daughter(1))->P();
175 
176  TLorentzVector pos1 = pi_serv->TrackIdToParticle_P(particle->Daughter(0))->EndPosition();
177  TLorentzVector pos2 = pi_serv->TrackIdToParticle_P(particle->Daughter(1))->EndPosition();
178 
179  if (insideFidVol(pos1)) fInside1 = true;
180  if (insideFidVol(pos2)) fInside2 = true;
181 
182  fConvgamma1.SetXYZ(pos1.X(), pos1.Y(), pos1.Z());
183  fConvgamma2.SetXYZ(pos2.X(), pos2.Y(), pos2.Z());
184 
185  TVector3 vecnorm1 = mom1vec3.Unit();
186  fDirgamma1 = vecnorm1;
187  TVector3 vecnorm2 = mom2vec3.Unit();
188  fDirgamma2 = vecnorm2;
189 
191  }
192  else {
193  fNgammas = particle->NumberDaughters();
194  }
195  }
196 }
197 
198 bool ems::MCinfo::insideFidVol(const TLorentzVector& pvtx) const
199 {
200 
201  bool inside = false;
202  //x
203  double dista = fabs(fMinx - pvtx.X());
204  double distb = fabs(pvtx.X() - fMaxx);
205  if ((pvtx.X() > fMinx) && (pvtx.X() < fMaxx) && (dista > fFidVolCut) && (distb > fFidVolCut))
206  inside = true;
207  //y
208  dista = fabs(fMaxy - pvtx.Y());
209  distb = fabs(pvtx.Y() - fMiny);
210  if (inside && (pvtx.Y() > fMiny) && (pvtx.Y() < fMaxy) && (dista > fFidVolCut) &&
211  (distb > fFidVolCut))
212  inside = true;
213  else
214  inside = false;
215 
216  //z
217  dista = fabs(fMaxz - pvtx.Z());
218  distb = fabs(pvtx.Z() - fMinz);
219  if (inside && (pvtx.Z() > fMinz) && (pvtx.Z() < fMaxz) && (dista > fFidVolCut) &&
220  (distb > fFidVolCut))
221  inside = true;
222  else
223  inside = false;
224 
225  return inside;
226 }
227 
229 public:
230  explicit MultiEMShowers(fhicl::ParameterSet const& p);
231 
232  MultiEMShowers(MultiEMShowers const&) = delete;
233  MultiEMShowers(MultiEMShowers&&) = delete;
234  MultiEMShowers& operator=(MultiEMShowers const&) = delete;
235  MultiEMShowers& operator=(MultiEMShowers&&) = delete;
236 
237 private:
238  void beginJob() override;
239  void endJob() override;
240 
241  void analyze(art::Event const& e) override;
242 
243  bool convCluster(art::Event const& evt);
244  double getMinDist(detinfo::DetectorPropertiesData const& detProp,
246  TVector3 const& convmc,
247  size_t view,
248  size_t tpc,
249  size_t cryo);
254 
255  // ROOT
256  TTree* fEvTree;
258  int fNGroups;
259 
260  // mc
261  double fPi0mom;
262  double fGmom1;
263  double fGmom2;
264  double fMcth;
265  int fNgammas;
267  int fEvComp;
269  int fEvInput;
270  TVector3 fGdir1;
271  TVector3 fGdir2;
272  TVector3 fPrimary;
273 
274  //reco
275  int fEvReco;
277  int fEv2Good;
278  int fCountph;
280 
281  TTree* fShTree;
282  TTree* fRecoTree;
283  double fStartX;
284  double fStartY;
285  double fStartZ;
286  double fDedxZ;
287  double fDedxV;
288  double fDedxU;
289  double fMCrecovtx;
290  double fMCrecoTh;
293  double fRecth;
294  double fRecthgood;
297  double fGdirmcreco1;
298  double fGdirmcreco2;
301 
307 };
308 
310 {
311  fConvGood = 0;
312  fConvWrong = 0;
313  fConvBothGood = 0;
314  fEvFidVol = 0;
315  fEvComp = 0;
316  fEvGMomCut = 0;
317  fEvReco = 0;
318  fEvInput = 0;
319  fEv2Groups = 0;
320  fEv2Good = 0;
321  fHitsModuleLabel = p.get<art::InputTag>("HitsModuleLabel");
322  fCluModuleLabel = p.get<art::InputTag>("ClustersModuleLabel");
323  fTrk3DModuleLabel = p.get<art::InputTag>("Trk3DModuleLabel");
324  fVtxModuleLabel = p.get<art::InputTag>("VtxModuleLabel");
325  fShsModuleLabel = p.get<art::InputTag>("ShsModuleLabel");
326 }
327 
329 {
331 
332  fEvTree = tfs->make<TTree>("MultiShowers", "showers3d");
333  fEvTree->Branch("fEvNumber", &fEvNumber, "fEvNumber/I");
334  fEvTree->Branch("fNGroups", &fNGroups, "fNGroups/I");
335  fEvTree->Branch("fPi0mom", &fPi0mom, "fPi0mom/D");
336  fEvTree->Branch("fNgammas", &fNgammas, "fNgammas/I");
337  fEvTree->Branch("fGmom1", &fGmom1, "fGmom1/D");
338  fEvTree->Branch("fGmom2", &fGmom2, "fGmom2/D");
339  fEvTree->Branch("fMcth", &fMcth, "fMcth/D");
340  fEvTree->Branch("fRecth", &fRecth, "fRecth/D");
341  fEvTree->Branch("fMCrecovtx", &fMCrecovtx, "fMCrecovtx/D");
342  fEvTree->Branch("fMCrecoTh", &fMCrecoTh, "fMCrecoTh/D");
343  fEvTree->Branch("fConvGood", &fConvGood, "fConvGood/I");
344  fEvTree->Branch("fConvWrong", &fConvWrong, "fConvWrong/I");
345  fEvTree->Branch("fConvBothGood", &fConvBothGood, "fConvBothGood/I");
346  fEvTree->Branch("fGammasInside", &fGammasInside, "fGammasInside/I");
347  fEvTree->Branch("fCountph", &fCountph, "fCountph/I");
348  fEvTree->Branch("fCountreco", &fCountreco, "fCountreco/I");
349 
350  fRecoTree = tfs->make<TTree>("Cascades", "conv points");
351  fRecoTree->Branch("fRecth", &fRecth, "fRecth/D");
352  fRecoTree->Branch("fMCrecovtx", &fMCrecovtx, "fMCrecovtx/D");
353  fRecoTree->Branch("fMCrecoTh", &fMCrecoTh, "fMCrecoTh/D");
354  fRecoTree->Branch("fRecthgood", &fRecthgood, "fRecthgood/D");
355  fRecoTree->Branch("fMCrecovtxgood", &fMCrecovtxgood, "fMCrecovtxgood/D");
356  fRecoTree->Branch("fMCrecoThgood", &fMCrecoThgood, "fMCrecoThgood/D");
357  fRecoTree->Branch("fGdirmcreco1", &fGdirmcreco1, "fGdirmcreco1/D");
358  fRecoTree->Branch("fGdirmcreco2", &fGdirmcreco2, "fGdirmcreco2/D");
359  fRecoTree->Branch("fGdirmcreco1good", &fGdirmcreco1good, "fGdirmcreco1good/D");
360  fRecoTree->Branch("fGdirmcreco2good", &fGdirmcreco2good, "fGdirmcreco2good/D");
361 
362  fShTree = tfs->make<TTree>("Shower", "conv point");
363 
364  fShTree->Branch("fStartX", &fStartX, "fStartX/D");
365  fShTree->Branch("fStartY", &fStartY, "fStartY/D");
366  fShTree->Branch("fStartZ", &fStartZ, "fStartZ/D");
367  fShTree->Branch("fDedxZ", &fDedxZ, "fDedxZ/D");
368  fShTree->Branch("fDedxV", &fDedxV, "fDedxV/D");
369  fShTree->Branch("fDedxU", &fDedxU, "fDedxU/D");
370 }
371 
373 {
374  mf::LogInfo log("MultiEMShower");
375  log << "******************** fEvFidVol = " << fEvFidVol << "\n";
376  log << "******************** fEvGMomCut = " << fEvGMomCut << "\n";
377  log << "******************** fEvComp = " << fEvComp << "\n";
378  log << "******************** fEvReco = " << fEvReco << "\n";
379  log << "******************** fEvInput = " << fEvInput << "\n";
380  log << "******************** fEv2Groups = " << fEv2Groups << "\n";
381  log << "******************** fEv2Good = " << fEv2Good << "\n";
382  if (fEvInput)
383  log << "******************** reco % = " << double(fEvReco) / double(fEvInput) << "\n";
384 }
385 
387 {
388  fEvNumber = e.id().event();
389  fNGroups = 0;
390  fStartX = 0.0;
391  fStartY = 0.0;
392  fStartZ = 0.0;
393  fPi0mom = 0.0;
394  fNgammas = 0;
395  fDistConvrecomc1 = 0.0;
396  fDistConvrecomc2 = 0.0;
397  fMCrecovtx = -400.0;
398  fMCrecovtxgood = -400.0;
399  fRecth = -400.0;
400  fRecthgood = -400.0;
401  fMCrecoTh = -400.0;
402  fMCrecoThgood = -400.0;
403  fGammasInside = 0;
404  fCountph = 0;
405  fCountreco = 0;
406  fGdirmcreco1 = 0.0;
407  fGdirmcreco2 = 0.0;
408  fGdirmcreco1good = 0.0;
409  fGdirmcreco2good = 0.0;
410  fDedxZ = 0.0;
411  fDedxV = 0.0;
412  fDedxU = 0.0;
413 
414  ems::MCinfo mc;
415  fPrimary = mc.GetPrimary();
416  fPi0mom = mc.GetMompi0();
417  fGmom1 = mc.GetMomGamma1();
418  fGmom2 = mc.GetMomGamma2();
419  fGdir1 = mc.GetDirgamma1();
420  fGdir2 = mc.GetDirgamma2();
421  fNgammas = mc.GetNgammas();
422  TVector3 pospi0 = mc.GetPospi0();
423 
424  double cosinemc = mc.GetCosine();
425  fMcth = 180.0F * (std::acos(cosinemc)) / TMath::Pi();
426  TVector3 convp[2];
427  convp[0] = mc.GetPosgamma1();
428  convp[1] = mc.GetPosgamma2();
429  const double maxdist = 2.0; //cm
430 
431  // check whether two photons are inside fid vol
432  if (mc.IsInside1() && mc.IsInside2()) {
433  fGammasInside = 1;
434  fEvFidVol++;
435  }
436 
437  if ((fGmom1 > 0.1) && (fGmom2 > 0.1)) fEvGMomCut++;
438 
439  if (mc.IsCompton()) fEvComp++;
440 
446  if (e.getByLabel(fShsModuleLabel, shsListHandle) &&
447  e.getByLabel(fTrk3DModuleLabel, trkListHandle) &&
448  e.getByLabel(fVtxModuleLabel, vtxListHandle) &&
449  e.getByLabel(fCluModuleLabel, cluListHandle) &&
450  e.getByLabel(fHitsModuleLabel, hitListHandle)) {
451  art::FindManyP<recob::Cluster> cluFromShs(shsListHandle, e, fShsModuleLabel);
452  art::FindManyP<recob::Cluster> cluFromTrk(trkListHandle, e, fTrk3DModuleLabel);
453  art::FindManyP<recob::Vertex> vtxFromTrk(trkListHandle, e, fVtxModuleLabel);
454  art::FindManyP<recob::Hit> hitFromClu(cluListHandle, e, fCluModuleLabel);
455 
456  fNGroups = shsListHandle->size();
457 
458  fCountph = 0;
459  if (fNgammas == 2) // pi0
460  {
461  int idph = -1;
462  for (size_t s = 0; s < shsListHandle->size(); ++s) {
463  const recob::Shower& sh = (*shsListHandle)[s];
464  double mindist = maxdist;
465  bool found = false;
466 
467  for (int i = 0; i < fNgammas; ++i) {
468  double dist = sqrt(pma::Dist2(sh.ShowerStart(), convp[i]));
469  if ((dist < mindist) && (idph != i)) {
470  mindist = dist;
471  idph = i;
472  found = true;
473  }
474  }
475  if (found) {
476  fConvGood++;
477  fCountph++;
478  }
479  else {
480  fConvWrong++;
481  }
482  }
483  if (fCountph == 2) fConvBothGood++;
484 
485  // plot a few variables if there are 2 showers
486  if (fCountph == 2)
487  for (size_t s = 0; s < shsListHandle->size(); ++s) {
488  const recob::Shower& sh = (*shsListHandle)[s];
489  TVector3 pos = sh.ShowerStart();
490  fStartX = pos.X();
491  fStartY = pos.Y();
492  fStartZ = pos.Z();
493  std::vector<double> const& vecdedx = sh.dEdx();
494 
495  if (vecdedx.size() == 3) {
496  fDedxZ = vecdedx[0];
497  fDedxV = vecdedx[1];
498  fDedxU = vecdedx[2];
499  }
500 
501  fShTree->Fill();
502  }
503  }
504  else // other than pi0
505  {
506  for (size_t s = 0; s < shsListHandle->size(); ++s) {
507  const recob::Shower& sh = (*shsListHandle)[s];
508  double mindist = maxdist;
509 
510  double dist = sqrt(pma::Dist2(sh.ShowerStart(), fPrimary));
511  if (dist < mindist) {
512  TVector3 pos = sh.ShowerStart();
513  fStartX = pos.X();
514  fStartY = pos.Y();
515  fStartZ = pos.Z();
516  std::vector<double> vecdedx = sh.dEdx();
517  if (vecdedx.size() == 3) {
518  fDedxZ = vecdedx[0];
519  fDedxV = vecdedx[1];
520  fDedxU = vecdedx[2];
521  }
522  }
523 
524  fShTree->Fill();
525  }
526  }
527  // compute the crossing point
528 
529  //cut from mc and clusters
530 
531  if (mc.IsInside1() && mc.IsInside2() && (fGmom1 > 0.1) && (fGmom2 > 0.1) && (!mc.IsCompton()) &&
532  convCluster(e)) {
533  fCountreco = 1;
534  if (fNGroups == 2) fEv2Groups++;
535  if ((fNGroups == 2) && (fCountph == 2)) fEv2Good++;
536  // cut from reco
537  //if (countph == 2)
538  if (fNGroups == 2) {
539  std::vector<std::pair<TVector3, TVector3>> lines;
540  const recob::Shower& sh1 = (*shsListHandle)[0];
541  const recob::Shower& sh2 = (*shsListHandle)[1];
542 
543  std::pair<TVector3, TVector3> frontback1(sh1.ShowerStart(),
544  sh1.ShowerStart() + sh1.Direction());
545  std::pair<TVector3, TVector3> frontback2(sh2.ShowerStart(),
546  sh2.ShowerStart() + sh2.Direction());
547  lines.push_back(frontback1);
548  lines.push_back(frontback2);
549 
550  TVector3 result;
551  pma::SolveLeastSquares3D(lines, result); // mse.
552 
553  double dist1_0 = pma::Dist2(result, sh1.ShowerStart());
554  double dist2_0 = pma::Dist2(result, sh1.ShowerStart() + sh1.Direction());
555  double dist1_1 = pma::Dist2(result, sh2.ShowerStart());
556  double dist2_1 = pma::Dist2(result, sh2.ShowerStart() + sh2.Direction());
557  if ((dist1_0 > dist2_0) || (dist1_1 > dist2_1)) {}
558  else {
559  fMCrecovtx = std::sqrt(pma::Dist2(pospi0, result));
560 
561  if (fCountph == 2) fMCrecovtxgood = fMCrecovtx;
562 
563  double cosine_reco = sh1.Direction() * sh2.Direction();
564  fRecth = 180.0F * (std::acos(cosine_reco)) / TMath::Pi();
565 
566  fGdirmcreco1 = fGdir1 * sh1.Direction();
567  fGdirmcreco2 = fGdir2 * sh2.Direction();
568  if (fCountph == 2) {
571  }
572 
573  if (fCountph == 2) fRecthgood = fRecth;
574 
575  fMCrecoTh = fRecth - fMcth;
576 
577  if (fCountph == 2) fMCrecoThgood = fMCrecoTh;
578 
579  fEvReco++;
580  fRecoTree->Fill();
581  }
582  }
583  fEvInput++;
584  //fRecoTree->Fill();
585  }
586  }
587 
588  fEvTree->Fill();
589 }
590 
591 // true if there are clusters corresponding to mc conversion points
593 {
594  ems::MCinfo mc;
595  TVector3 const convp[2]{mc.GetPosgamma1(), mc.GetPosgamma2()};
596 
597  geo::Point_t const vtx{convp[0].X(), convp[0].Y(), convp[0].Z()};
598 
600  geo::TPCID idtpc = geom->FindTPCAtPosition(vtx);
601  geo::CryostatID const cryoid = geom->PositionToCryostatID(vtx);
602 
605 
606  //map: conversion point, vec of id clusters in each view
607  std::map<size_t, std::vector<size_t>> used;
608 
609  art::FindManyP<recob::Hit> fbc(cluListHandle, evt, fCluModuleLabel);
610 
611  auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(evt);
612  auto const& readoutGeom = art::ServiceHandle<geo::WireReadout>()->Get();
613  double maxdist = 1.0; // 1 cm
614  if (geom->HasTPC(idtpc)) {
615  if (evt.getByLabel(fHitsModuleLabel, hitListHandle) &&
616  evt.getByLabel(fCluModuleLabel, cluListHandle)) {
617  size_t conv = 0;
618  while (conv < 2) {
619  for (auto const& plane : readoutGeom.Iterate<geo::PlaneGeo>(cryoid)) {
620 
621  double mindist = maxdist;
622  int clid = 0;
623  for (size_t c = 0; c < cluListHandle->size(); ++c) {
624 
625  bool exist = false;
626  for (auto const& ids : used)
627  for (auto i : ids.second)
628  if (i == c) exist = true;
629  if (exist) continue;
630 
631  std::vector<art::Ptr<recob::Hit>> hits = fbc.at(c);
632  if (hits.size() < 20) continue;
633  if (hits[0]->WireID().Plane != plane.View()) continue;
634 
635  double dist =
636  getMinDist(detProp, hits, convp[conv], plane.View(), idtpc.TPC, cryoid.Cryostat);
637  if (dist < mindist) {
638  mindist = dist;
639  clid = c;
640  }
641  }
642  if (mindist < maxdist) used[conv].push_back(clid);
643  }
644  conv++;
645  }
646  }
647  }
648  bool result = false;
649 
650  if (used.size() > 1)
651  for (auto const& ids : used) {
652  if (ids.second.size() > 1)
653  result = true;
654  else {
655  result = false;
656  break;
657  }
658  }
659 
660  return result;
661 }
662 
665  TVector3 const& convmc,
666  size_t view,
667  size_t tpc,
668  size_t cryo)
669 {
670  double mindist = 9999;
671  // MC vertex projected to view
672  TVector2 proj = pma::GetProjectionToPlane(convmc, view, tpc, cryo);
673 
674  // loop over hits to find the closest to MC 2d vtx
675  for (size_t h = 0; h < v.size(); ++h) {
676  if ((v[h]->WireID().Plane == view) && (v[h]->WireID().TPC == tpc)) {
677  TVector2 hpoint =
678  pma::WireDriftToCm(detProp, v[h]->WireID().Wire, v[h]->PeakTime(), view, tpc, cryo);
679 
680  double dist = pma::Dist2(proj, hpoint);
681  if (dist < mindist) { mindist = dist; }
682  }
683  }
684 
685  return mindist;
686 }
687 
const TVector3 & ShowerStart() const
Definition: Shower.h:197
TVector3 const & GetDirgamma2() const
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:220
int PdgCode() const
Definition: MCParticle.h:213
double SolveLeastSquares3D(const std::vector< std::pair< TVector3, TVector3 >> &lines, TVector3 &result)
Definition: Utilities.cxx:158
double Dist2(const TVector2 &v1, const TVector2 &v2)
Definition: Utilities.cxx:39
const simb::MCParticle * TrackIdToParticle_P(int id) const
MultiEMShowers(fhicl::ParameterSet const &p)
Declaration of signal hit object.
list_type::const_iterator const_iterator
Definition: ParticleList.h:132
Geometry information for a single TPC.
Definition: TPCGeo.h:33
void analyze(art::Event const &e) override
bool convCluster(art::Event const &evt)
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:195
TVector2 WireDriftToCm(detinfo::DetectorPropertiesData const &detProp, unsigned int wire, float drift, unsigned int plane, unsigned int tpc, unsigned int cryo)
Definition: Utilities.cxx:286
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
std::string Process() const
Definition: MCParticle.h:216
double GetMomGamma2() const
int NumberDaughters() const
Definition: MCParticle.h:218
double GetMompi0() const
int Daughter(const int i) const
Definition: MCParticle.cxx:118
double getMinDist(detinfo::DetectorPropertiesData const &detProp, std::vector< art::Ptr< recob::Hit >> const &v, TVector3 const &convmc, size_t view, size_t tpc, size_t cryo)
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:289
void hits()
Definition: readHits.C:15
TPCID FindTPCAtPosition(Point_t const &point) const
Returns the ID of the TPC at specified location.
const std::vector< double > & dEdx() const
Definition: Shower.h:235
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
std::string EndProcess() const
Definition: MCParticle.h:217
void beginJob()
Definition: Breakpoints.cc:14
TVector3 const & GetPosgamma1() const
double P(const int i=0) const
Definition: MCParticle.h:235
T get(std::string const &key) const
Definition: ParameterSet.h:314
int GetNgammas() const
Geometry information for a single wire plane.The plane is represented in the geometry by a solid whic...
Definition: PlaneGeo.h:67
iterator begin()
Definition: ParticleList.h:305
const TVector3 & Direction() const
Definition: Shower.h:188
The data type to uniquely identify a TPC.
Definition: geo_types.h:306
Declaration of cluster object.
TVector3 const & GetPrimary() const
TVector3 const & GetDirgamma1() const
Definition: DirOfGamma.h:22
const sim::ParticleList & ParticleList() const
bool IsInside2() const
CryostatID PositionToCryostatID(Point_t const &point) const
Returns the ID of the cryostat at specified location.
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
Definition: geo_vectors.h:180
bool IsCompton() const
Implementation of the Projection Matching Algorithm.
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
TVector3 const & GetPospi0() const
double GetMomGamma1() const
double dist(const TReal *x, const TReal *y, const unsigned int dimension)
Provides recob::Track data product.
TVector2 GetProjectionToPlane(const TVector3 &p, unsigned int plane, unsigned int tpc, unsigned int cryo)
Definition: Utilities.cxx:263
Float_t proj
Definition: plot.C:35
EventNumber_t event() const
Definition: EventID.h:116
TCEvent evt
Definition: DataStructs.cxx:8
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:315
Float_t e
Definition: plot.C:35
TVector3 const & GetPosgamma2() const
recob::tracking::Plane Plane
Definition: TrackState.h:17
bool IsInside1() const
bool HasTPC(TPCID const &tpcid) const
Returns whether we have the specified TPC.
Definition: GeometryCore.h:429
EventID id() const
Definition: Event.cc:23
art framework interface to geometry description
The data type to uniquely identify a cryostat.
Definition: geo_types.h:187
bool insideFidVol(const TLorentzVector &pvtx) const