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