LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
TCShowerTemplateMaker_module.cc
Go to the documentation of this file.
1 // -------------------------------------------------
2 // makes shower profile templates
3 //
4 // Author: Rory Fitzpatrick (roryfitz@umich.edu)
5 // Created: 7/16/18
6 // -------------------------------------------------
7 
8 // Framework includes
14 #include "art_root_io/TFileService.h"
16 #include "fhiclcpp/ParameterSet.h"
17 
29 
30 #include "TH1.h"
31 #include "TH3.h"
32 #include "TProfile.h"
33 #include "TProfile2D.h"
34 
35 namespace shower {
36 
38  public:
39  explicit TCShowerTemplateMaker(fhicl::ParameterSet const& pset);
40 
41  private:
42  void beginJob() override;
43  void analyze(const art::Event& evt) override;
44 
45  void showerProfile(detinfo::DetectorClocksData const& clockData,
46  detinfo::DetectorPropertiesData const& detProp,
48  TVector3 shwvtx,
49  TVector3 shwdir,
50  double elep);
51  void showerProfileTrue(detinfo::DetectorClocksData const& clockData,
52  detinfo::DetectorPropertiesData const& detProp,
54  double elep);
56  simb::MCParticle electron);
57 
61 
65 
69 
73 
77 
81 
85 
89 
93 
95  TH3F* fTransverse;
101 
102  // electrons
109  // photons
116  // other
123 
124  const int LBINS = 20;
125  const int LMIN = 0;
126  const int LMAX = 5;
127 
128  const int TBINS = 20;
129  const int TMIN = -5;
130  const int TMAX = 5;
131 
132  const int EBINS = 20;
133  const double EMIN = 0.5;
134  const double EMAX = 20.5;
135 
136  const double X0 = 14;
137 
138  std::string fHitModuleLabel;
139  std::string fShowerModuleLabel;
140  std::string fGenieGenModuleLabel;
141  std::string fDigitModuleLabel;
142 
144 
145  }; // class TCShowerAnalysis
146 
147 } // shower
148 
149 // -------------------------------------------------
150 
152  : EDAnalyzer(pset)
153  , fHitModuleLabel(pset.get<std::string>("HitModuleLabel", "trajcluster"))
154  , fShowerModuleLabel(pset.get<std::string>("ShowerModuleLabel", "tcshower"))
155  , fGenieGenModuleLabel(pset.get<std::string>("GenieGenModuleLabel", "generator"))
156  , fDigitModuleLabel(pset.get<std::string>("DigitModuleLabel", "largeant"))
157  , fCalorimetryAlg(pset.get<fhicl::ParameterSet>("CalorimetryAlg"))
158 {} // TCShowerTemplateMaker
159 
160 // -------------------------------------------------
161 
163 {
164 
167  tfs->make<TProfile>("fShowerProfileSimLong",
168  "longitudinal e- profile (true, simchannel);t;E (MeV)",
169  LBINS,
170  LMIN,
171  LMAX);
172  fShowerProfileHitLong = tfs->make<TProfile>(
173  "fShowerProfileHitLong", "longitudinal e- profile (true, hit);t;E (MeV)", LBINS, LMIN, LMAX);
174  fShowerProfileRecoLong = tfs->make<TProfile>(
175  "fShowerProfileRecoLong", "longitudinal e- profile (reco);t;Q", LBINS, LMIN, LMAX);
176 
177  fShowerProfileSimLong2D = tfs->make<TProfile2D>(
178  "fShowerProfileSimLong2D",
179  "longitudinal e- profile (true, simchannel);t;electron energy (MeV);E (MeV)",
180  LBINS,
181  LMIN,
182  LMAX,
183  EBINS,
184  EMIN,
185  EMAX);
187  tfs->make<TProfile2D>("fShowerProfileHitLong2D",
188  "longitudinal e- profile (true, hit);t;electron energy (MeV);E (MeV)",
189  LBINS,
190  LMIN,
191  LMAX,
192  EBINS,
193  EMIN,
194  EMAX);
196  tfs->make<TProfile2D>("fShowerProfileRecoLong2D",
197  "longitudinal e- profile (reco);t;electron energy (MeV);Q",
198  LBINS,
199  LMIN,
200  LMAX,
201  EBINS,
202  EMIN,
203  EMAX);
204 
206  tfs->make<TProfile>("fShowerProfileSimTrans",
207  "transverse e- profile (true, simchannel);dist (cm);E (MeV)",
208  TBINS,
209  TMIN,
210  TMAX);
212  tfs->make<TProfile>("fShowerProfileHitTrans",
213  "transverse e- profile (true, hit);dist (cm);E (MeV)",
214  TBINS,
215  TMIN,
216  TMAX);
217  fShowerProfileRecoTrans = tfs->make<TProfile>(
218  "fShowerProfileRecoTrans", "transverse e- profile (reco);dist (cm);Q", TBINS, TMIN, TMAX);
219 
220  fShowerProfileSimTrans2D = tfs->make<TProfile2D>(
221  "fShowerProfileSimTrans2D",
222  "transverse e- profile (true, simchannel);t;electron energy (MeV);E (MeV)",
223  TBINS,
224  TMIN,
225  TMAX,
226  EBINS,
227  EMIN,
228  EMAX);
230  tfs->make<TProfile2D>("fShowerProfileHitTrans2D",
231  "transverse e- profile (true, hit);t;electron energy (MeV);E (MeV)",
232  TBINS,
233  TMIN,
234  TMAX,
235  EBINS,
236  EMIN,
237  EMAX);
239  tfs->make<TProfile2D>("fShowerProfileRecoTrans2D",
240  "transverse e- profile (reco);t;electron energy (MeV);Q",
241  TBINS,
242  TMIN,
243  TMAX,
244  EBINS,
245  EMIN,
246  EMAX);
247 
248  fShowerProfileSimTrans2D_1 = tfs->make<TProfile2D>(
249  "fShowerProfileSimTrans2D_1",
250  "transverse e- profile [0 <= t < 1] (true, simchannel);t;electron energy (MeV);E (MeV)",
251  TBINS,
252  TMIN,
253  TMAX,
254  EBINS,
255  EMIN,
256  EMAX);
257  fShowerProfileHitTrans2D_1 = tfs->make<TProfile2D>(
258  "fShowerProfileHitTrans2D_1",
259  "transverse e- profile [0 <= t < 1] (true, hit);t;electron energy (MeV);E (MeV)",
260  TBINS,
261  TMIN,
262  TMAX,
263  EBINS,
264  EMIN,
265  EMAX);
267  tfs->make<TProfile2D>("fShowerProfileRecoTrans2D_1",
268  "transverse e- profile [0 <= t < 1] (reco);t;electron energy (MeV);Q",
269  TBINS,
270  TMIN,
271  TMAX,
272  EBINS,
273  EMIN,
274  EMAX);
275 
276  fShowerProfileSimTrans2D_2 = tfs->make<TProfile2D>(
277  "fShowerProfileSimTrans2D_2",
278  "transverse e- profile [1 <= t < 2] (true, simchannel);t;electron energy (MeV);E (MeV)",
279  TBINS,
280  TMIN,
281  TMAX,
282  EBINS,
283  EMIN,
284  EMAX);
285  fShowerProfileHitTrans2D_2 = tfs->make<TProfile2D>(
286  "fShowerProfileHitTrans2D_2",
287  "transverse e- profile [1 <= t < 2] (true, hit);t;electron energy (MeV);E (MeV)",
288  TBINS,
289  TMIN,
290  TMAX,
291  EBINS,
292  EMIN,
293  EMAX);
295  tfs->make<TProfile2D>("fShowerProfileRecoTrans2D_2",
296  "transverse e- profile [1 <= t < 2] (reco);t;electron energy (MeV);Q",
297  TBINS,
298  TMIN,
299  TMAX,
300  EBINS,
301  EMIN,
302  EMAX);
303 
304  fShowerProfileSimTrans2D_3 = tfs->make<TProfile2D>(
305  "fShowerProfileSimTrans2D_3",
306  "transverse e- profile [2 <= t < 3] (true, simchannel);t;electron energy (MeV);E (MeV)",
307  TBINS,
308  TMIN,
309  TMAX,
310  EBINS,
311  EMIN,
312  EMAX);
313  fShowerProfileHitTrans2D_3 = tfs->make<TProfile2D>(
314  "fShowerProfileHitTrans2D_3",
315  "transverse e- profile [2 <= t < 3] (true, hit);t;electron energy (MeV);E (MeV)",
316  TBINS,
317  TMIN,
318  TMAX,
319  EBINS,
320  EMIN,
321  EMAX);
323  tfs->make<TProfile2D>("fShowerProfileRecoTrans2D_3",
324  "transverse e- profile [2 <= t < 3] (reco);t;electron energy (MeV);Q",
325  TBINS,
326  TMIN,
327  TMAX,
328  EBINS,
329  EMIN,
330  EMAX);
331 
332  fShowerProfileSimTrans2D_4 = tfs->make<TProfile2D>(
333  "fShowerProfileSimTrans2D_4",
334  "transverse e- profile [3 <= t < 4] (true, simchannel);t;electron energy (MeV);E (MeV)",
335  TBINS,
336  TMIN,
337  TMAX,
338  EBINS,
339  EMIN,
340  EMAX);
341  fShowerProfileHitTrans2D_4 = tfs->make<TProfile2D>(
342  "fShowerProfileHitTrans2D_4",
343  "transverse e- profile [3 <= t < 4] (true, hit);t;electron energy (MeV);E (MeV)",
344  TBINS,
345  TMIN,
346  TMAX,
347  EBINS,
348  EMIN,
349  EMAX);
351  tfs->make<TProfile2D>("fShowerProfileRecoTrans2D_4",
352  "transverse e- profile [3 <= t < 4] (reco);t;electron energy (MeV);Q",
353  TBINS,
354  TMIN,
355  TMAX,
356  EBINS,
357  EMIN,
358  EMAX);
359 
360  fShowerProfileSimTrans2D_5 = tfs->make<TProfile2D>(
361  "fShowerProfileSimTrans2D_5",
362  "transverse e- profile [4 <= t < 5] (true, simchannel);t;electron energy (MeV);E (MeV)",
363  TBINS,
364  TMIN,
365  TMAX,
366  EBINS,
367  EMIN,
368  EMAX);
369  fShowerProfileHitTrans2D_5 = tfs->make<TProfile2D>(
370  "fShowerProfileHitTrans2D_5",
371  "transverse e- profile [4 <= t < 5] (true, hit);t;electron energy (MeV);E (MeV)",
372  TBINS,
373  TMIN,
374  TMAX,
375  EBINS,
376  EMIN,
377  EMAX);
379  tfs->make<TProfile2D>("fShowerProfileRecoTrans2D_5",
380  "transverse e- profile [4 <= t < 5] (reco);t;electron energy (MeV);Q",
381  TBINS,
382  TMIN,
383  TMAX,
384  EBINS,
385  EMIN,
386  EMAX);
387 
388  fLongitudinal = tfs->make<TH3F>("fLongitudinal",
389  "longitudinal e- profile;t;electron energy (MeV);Q",
390  LBINS,
391  LMIN,
392  LMAX,
393  EBINS,
394  EMIN,
395  EMAX,
396  100,
397  0,
398  150000);
399  fTransverse = tfs->make<TH3F>("fTransverse",
400  "transverse e- profile;dist (cm);electron energy (MeV);Q",
401  TBINS,
402  TMIN,
403  TMAX,
404  EBINS,
405  EMIN,
406  EMAX,
407  100,
408  0,
409  150000);
410  fTransverse_1 =
411  tfs->make<TH3F>("fTransverse_1",
412  "transverse e- profile [0 <= t < 1];dist (cm);electron energy (MeV);Q",
413  TBINS,
414  TMIN,
415  TMAX,
416  EBINS,
417  EMIN,
418  EMAX,
419  100,
420  0,
421  100000);
422  fTransverse_2 =
423  tfs->make<TH3F>("fTransverse_2",
424  "transverse e- profile [1 <= t < 2];dist (cm);electron energy (MeV);Q",
425  TBINS,
426  TMIN,
427  TMAX,
428  EBINS,
429  EMIN,
430  EMAX,
431  100,
432  0,
433  100000);
434  fTransverse_3 =
435  tfs->make<TH3F>("fTransverse_3",
436  "transverse e- profile [2 <= t < 3];dist (cm);electron energy (MeV);Q",
437  TBINS,
438  TMIN,
439  TMAX,
440  EBINS,
441  EMIN,
442  EMAX,
443  100,
444  0,
445  100000);
446  fTransverse_4 =
447  tfs->make<TH3F>("fTransverse_4",
448  "transverse e- profile [3 <= t < 4];dist (cm);electron energy (MeV);Q",
449  TBINS,
450  TMIN,
451  TMAX,
452  EBINS,
453  EMIN,
454  EMAX,
455  100,
456  0,
457  100000);
458  fTransverse_5 =
459  tfs->make<TH3F>("fTransverse_5",
460  "transverse e- profile [4 <= t < 5];dist (cm);electron energy (MeV);Q",
461  TBINS,
462  TMIN,
463  TMAX,
464  EBINS,
465  EMIN,
466  EMAX,
467  100,
468  0,
469  100000);
470 
471  // electrons
472  fLongitudinal_electron = tfs->make<TProfile>(
473  "fLongitudinal_electron", "longitudinal e- profile (reco);t;Q", LBINS, LMIN, LMAX);
475  tfs->make<TProfile>("fTransverse1_electron",
476  "transverse e- profile [0 <= t < 1] (reco);dist (cm);Q",
477  TBINS,
478  TMIN,
479  TMAX);
481  tfs->make<TProfile>("fTransverse2_electron",
482  "transverse e- profile [1 <= t < 2] (reco);dist (cm);Q",
483  TBINS,
484  TMIN,
485  TMAX);
487  tfs->make<TProfile>("fTransverse3_electron",
488  "transverse e- profile [2 <= t < 3] (reco);dist (cm);Q",
489  TBINS,
490  TMIN,
491  TMAX);
493  tfs->make<TProfile>("fTransverse4_electron",
494  "transverse e- profile [3 <= t < 4] (reco);dist (cm);Q",
495  TBINS,
496  TMIN,
497  TMAX);
499  tfs->make<TProfile>("fTransverse5_electron",
500  "transverse e- profile [4 <= t < 5] (reco);dist (cm);Q",
501  TBINS,
502  TMIN,
503  TMAX);
504 
505  // photons
506  fLongitudinal_photon = tfs->make<TProfile>(
507  "fLongitudinal_photon", "longitudinal photon profile (reco);t;Q", LBINS, LMIN, LMAX);
509  tfs->make<TProfile>("fTransverse1_photon",
510  "transverse photon profile [0 <= t < 1] (reco);dist (cm);Q",
511  TBINS,
512  TMIN,
513  TMAX);
515  tfs->make<TProfile>("fTransverse2_photon",
516  "transverse photon profile [1 <= t < 2] (reco);dist (cm);Q",
517  TBINS,
518  TMIN,
519  TMAX);
521  tfs->make<TProfile>("fTransverse3_photon",
522  "transverse photon profile [2 <= t < 3] (reco);dist (cm);Q",
523  TBINS,
524  TMIN,
525  TMAX);
527  tfs->make<TProfile>("fTransverse4_photon",
528  "transverse photon profile [3 <= t < 4] (reco);dist (cm);Q",
529  TBINS,
530  TMIN,
531  TMAX);
533  tfs->make<TProfile>("fTransverse5_photon",
534  "transverse photon profile [4 <= t < 5] (reco);dist (cm);Q",
535  TBINS,
536  TMIN,
537  TMAX);
538 
539  // other
540  fLongitudinal_other = tfs->make<TProfile>(
541  "fLongitudinal_other", "longitudinal other profile (reco);t;Q", LBINS, LMIN, LMAX);
543  tfs->make<TProfile>("fTransverse1_other",
544  "transverse other profile [0 <= t < 1] (reco);dist (cm);Q",
545  TBINS,
546  TMIN,
547  TMAX);
549  tfs->make<TProfile>("fTransverse2_other",
550  "transverse other profile [1 <= t < 2] (reco);dist (cm);Q",
551  TBINS,
552  TMIN,
553  TMAX);
555  tfs->make<TProfile>("fTransverse3_other",
556  "transverse other profile [2 <= t < 3] (reco);dist (cm);Q",
557  TBINS,
558  TMIN,
559  TMAX);
561  tfs->make<TProfile>("fTransverse4_other",
562  "transverse other profile [3 <= t < 4] (reco);dist (cm);Q",
563  TBINS,
564  TMIN,
565  TMAX);
567  tfs->make<TProfile>("fTransverse5_other",
568  "transverse other profile [4 <= t < 5] (reco);dist (cm);Q",
569  TBINS,
570  TMIN,
571  TMAX);
572 
573 } // beginJob
574 
575 // -------------------------------------------------
576 
578 {
579 
581  std::vector<art::Ptr<recob::Hit>> hitlist;
582  if (evt.getByLabel(fHitModuleLabel, hitListHandle)) art::fill_ptr_vector(hitlist, hitListHandle);
583 
585  std::vector<art::Ptr<sim::SimChannel>> simchanlist;
586  if (evt.getByLabel(fDigitModuleLabel, scListHandle))
587  art::fill_ptr_vector(simchanlist, scListHandle);
588 
589  art::Handle<std::vector<recob::Shower>> showerListHandle;
590  std::vector<art::Ptr<recob::Shower>> showerlist;
591  if (evt.getByLabel(fShowerModuleLabel, showerListHandle))
592  art::fill_ptr_vector(showerlist, showerListHandle);
593 
594  art::Handle<std::vector<simb::MCTruth>> mctruthListHandle;
595  std::vector<art::Ptr<simb::MCTruth>> mclist;
596  if (evt.getByLabel(fGenieGenModuleLabel, mctruthListHandle))
597  art::fill_ptr_vector(mclist, mctruthListHandle);
598 
599  art::FindManyP<recob::Hit> shwfm(showerListHandle, evt, fShowerModuleLabel);
600 
601  if (empty(mclist)) return;
602 
603  auto const clock_data = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
604  auto const det_prop =
606 
607  art::Ptr<simb::MCTruth> mctruth = mclist[0];
608  if (mctruth->NeutrinoSet()) {
609  if (std::abs(mctruth->GetNeutrino().Nu().PdgCode()) == 12 &&
610  mctruth->GetNeutrino().CCNC() == 0) {
611  double elep = mctruth->GetNeutrino().Lepton().E();
612  if (showerlist.size()) {
613  std::vector<art::Ptr<recob::Hit>> showerhits = shwfm.at(0);
614  showerProfile(clock_data,
615  det_prop,
616  showerhits,
617  showerlist[0]->ShowerStart(),
618  showerlist[0]->Direction(),
619  elep);
620  }
621  showerProfileTrue(clock_data, det_prop, hitlist, elep);
622  showerProfileTrue(simchanlist, mctruth->GetNeutrino().Lepton());
623  }
624  }
625 } // analyze
626 
627 // -------------------------------------------------
628 
630  detinfo::DetectorPropertiesData const& detProp,
631  std::vector<art::Ptr<recob::Hit>> showerhits,
632  TVector3 shwvtx,
633  TVector3 shwdir,
634  double elep)
635 {
637 
638  auto collectionPlane = geo::PlaneID(0, 0, 1);
639 
640  double shwVtxTime = detProp.ConvertXToTicks(shwvtx[0], collectionPlane);
641  using geo::vect::toPoint;
642  double shwVtxWire = geom->WireCoordinate(toPoint(shwvtx), collectionPlane);
643 
644  double shwTwoTime = detProp.ConvertXToTicks(shwvtx[0] + shwdir[0], collectionPlane);
645  double shwTwoWire = geom->WireCoordinate(toPoint(shwvtx + shwdir), collectionPlane);
646 
647  TH1F* ltemp = new TH1F("ltemp", "ltemp", LBINS, LMIN, LMAX);
648  TH1F* ttemp = new TH1F("ttemp", "ttemp", TBINS, TMIN, TMAX);
649 
650  TH1F* ttemp_1 = new TH1F("ttemp_1", "ttemp_1", TBINS, TMIN, TMAX);
651  TH1F* ttemp_2 = new TH1F("ttemp_2", "ttemp_2", TBINS, TMIN, TMAX);
652  TH1F* ttemp_3 = new TH1F("ttemp_3", "ttemp_3", TBINS, TMIN, TMAX);
653  TH1F* ttemp_4 = new TH1F("ttemp_4", "ttemp_4", TBINS, TMIN, TMAX);
654  TH1F* ttemp_5 = new TH1F("ttemp_5", "ttemp_5", TBINS, TMIN, TMAX);
655 
656  for (size_t i = 0; i < showerhits.size(); ++i) {
657  if (showerhits[i]->WireID().Plane != collectionPlane.Plane) continue;
658 
659  double wirePitch = geom->WirePitch(showerhits[i]->WireID());
660  double tickToDist = detProp.DriftVelocity(detProp.Efield(), detProp.Temperature());
661  tickToDist *= 1.e-3 * sampling_rate(clockData); // 1e-3 is conversion of 1/us to 1/ns
662 
663  double xvtx = shwVtxTime * tickToDist;
664  double yvtx = shwVtxWire * wirePitch;
665 
666  double xtwo = shwTwoTime * tickToDist;
667  double ytwo = shwTwoWire * wirePitch;
668 
669  double xtwoorth = (ytwo - yvtx) + xvtx;
670  double ytwoorth = -(xtwo - xvtx) + yvtx;
671 
672  double xhit = showerhits[i]->PeakTime() * tickToDist;
673  double yhit = showerhits[i]->WireID().Wire * wirePitch;
674 
675  double ldist = std::abs((ytwoorth - yvtx) * xhit - (xtwoorth - xvtx) * yhit + xtwoorth * yvtx -
676  ytwoorth * xvtx) /
677  std::hypot(ytwoorth - yvtx, xtwoorth - xvtx);
678  double tdist = ((ytwo - yvtx) * xhit - (xtwo - xvtx) * yhit + xtwo * yvtx - ytwo * xvtx) /
679  std::hypot(ytwo - yvtx, xtwo - xvtx);
680 
681  double to3D = 1. / std::hypot(xvtx - xtwo,
682  yvtx - ytwo); // distance between two points in 3D space is one
683  ldist *= to3D;
684  tdist *= to3D;
685  double t = ldist / X0;
686 
687  double Q = showerhits[i]->Integral() *
688  fCalorimetryAlg.LifetimeCorrection(clockData, detProp, showerhits[i]->PeakTime());
689 
690  ltemp->Fill(t, Q);
691  ttemp->Fill(tdist, Q);
692 
693  if (t < 1)
694  ttemp_1->Fill(tdist, Q);
695  else if (t < 2)
696  ttemp_2->Fill(tdist, Q);
697  else if (t < 3)
698  ttemp_3->Fill(tdist, Q);
699  else if (t < 4)
700  ttemp_4->Fill(tdist, Q);
701  else if (t < 5)
702  ttemp_5->Fill(tdist, Q);
703 
704  } // loop through showerhits
705 
706  for (int i = 0; i < LBINS; ++i) {
707  if (ltemp->GetBinContent(i + 1) == 0) continue;
708  fShowerProfileRecoLong->Fill(ltemp->GetBinCenter(i + 1), ltemp->GetBinContent(i + 1));
709  fShowerProfileRecoLong2D->Fill(ltemp->GetBinCenter(i + 1), elep, ltemp->GetBinContent(i + 1));
710  fLongitudinal->Fill(ltemp->GetBinCenter(i + 1), elep, ltemp->GetBinContent(i + 1));
711  }
712 
713  for (int i = 0; i < TBINS; ++i) {
714  if (ttemp->GetBinContent(i + 1) == 0) continue;
715  fShowerProfileRecoTrans->Fill(ttemp->GetBinCenter(i + 1), ttemp->GetBinContent(i + 1));
716  fShowerProfileRecoTrans2D->Fill(ttemp->GetBinCenter(i + 1), elep, ttemp->GetBinContent(i + 1));
717  fTransverse->Fill(ttemp->GetBinCenter(i + 1), elep, ttemp->GetBinContent(i + 1));
718 
720  ttemp_1->GetBinCenter(i + 1), elep, ttemp_1->GetBinContent(i + 1));
722  ttemp_2->GetBinCenter(i + 1), elep, ttemp_2->GetBinContent(i + 1));
724  ttemp_3->GetBinCenter(i + 1), elep, ttemp_3->GetBinContent(i + 1));
726  ttemp_4->GetBinCenter(i + 1), elep, ttemp_4->GetBinContent(i + 1));
728  ttemp_5->GetBinCenter(i + 1), elep, ttemp_5->GetBinContent(i + 1));
729 
730  fTransverse_1->Fill(ttemp_1->GetBinCenter(i + 1), elep, ttemp_1->GetBinContent(i + 1));
731  fTransverse_2->Fill(ttemp_2->GetBinCenter(i + 1), elep, ttemp_2->GetBinContent(i + 1));
732  fTransverse_3->Fill(ttemp_3->GetBinCenter(i + 1), elep, ttemp_3->GetBinContent(i + 1));
733  fTransverse_4->Fill(ttemp_4->GetBinCenter(i + 1), elep, ttemp_4->GetBinContent(i + 1));
734  fTransverse_5->Fill(ttemp_5->GetBinCenter(i + 1), elep, ttemp_5->GetBinContent(i + 1));
735  }
736 } // showerProfile
737 
738 // -------------------------------------------------
739 
741  detinfo::DetectorClocksData const& clockData,
742  detinfo::DetectorPropertiesData const& detProp,
744  double elep)
745 {
747  auto collectionPlane = geo::PlaneID(0, 0, 1);
750  std::map<int, double> trkID_E;
751 
752  TH1F* ltemp = new TH1F("ltemp", "ltemp", LBINS, LMIN, LMAX);
753  TH1F* ttemp = new TH1F("ttemp", "ttemp", TBINS, TMIN, TMAX);
754 
755  TH1F* ttemp_1 = new TH1F("ttemp_1", "ttemp_1", TBINS, TMIN, TMAX);
756  TH1F* ttemp_2 = new TH1F("ttemp_2", "ttemp_2", TBINS, TMIN, TMAX);
757  TH1F* ttemp_3 = new TH1F("ttemp_3", "ttemp_3", TBINS, TMIN, TMAX);
758  TH1F* ttemp_4 = new TH1F("ttemp_4", "ttemp_4", TBINS, TMIN, TMAX);
759  TH1F* ttemp_5 = new TH1F("ttemp_5", "ttemp_5", TBINS, TMIN, TMAX);
760 
761  double xvtx = -999;
762  double yvtx = -999;
763  double zvtx = -999;
764  double xtwo = -999;
765  double ytwo = -999;
766  double ztwo = -999;
767  double shwvtxT = -999;
768  double shwvtxW = -999;
769  double shwtwoT = -999;
770  double shwtwoW = -999;
771 
772  double shwvtxx = -999;
773  double shwvtxy = -999;
774  double shwtwox = -999;
775  double shwtwoy = -999;
776  double xtwoorth = -999;
777  double ytwoorth = -999;
778 
779  double wirePitch = -999;
780  double tickToDist = -999;
781 
782  bool foundParent = false;
783 
784  for (auto const& hit : allhits) {
785  if (hit->WireID().Plane != collectionPlane.Plane) continue;
786 
787  // art::Ptr<recob::Hit> hit = allhits[i];
788  std::vector<sim::TrackIDE> trackIDs = btserv->HitToEveTrackIDEs(clockData, hit);
789 
790  for (size_t j = 0; j < trackIDs.size(); ++j) {
791  // only want energy associated with the electron and electron must have neutrino mother
792  if (std::abs((piserv->TrackIdToParticle_P(trackIDs[j].trackID))->PdgCode()) != 11) continue;
793  if (std::abs((piserv->TrackIdToParticle_P(trackIDs[j].trackID))->Mother()) != 0) continue;
794 
795  if (!foundParent) {
796  xvtx = (piserv->TrackIdToParticle_P(trackIDs[j].trackID))->Vx();
797  yvtx = (piserv->TrackIdToParticle_P(trackIDs[j].trackID))->Vy();
798  zvtx = (piserv->TrackIdToParticle_P(trackIDs[j].trackID))->Vz();
799 
800  xtwo = xvtx + (piserv->TrackIdToParticle_P(trackIDs[j].trackID))->Px();
801  ytwo = yvtx + (piserv->TrackIdToParticle_P(trackIDs[j].trackID))->Py();
802  ztwo = zvtx + (piserv->TrackIdToParticle_P(trackIDs[j].trackID))->Pz();
803 
804  shwvtxT = detProp.ConvertXToTicks(xvtx, collectionPlane);
805  geo::Point_t const vtx{xvtx, yvtx, zvtx};
806  shwvtxW = geom->WireCoordinate(vtx, collectionPlane);
807 
808  shwtwoT = detProp.ConvertXToTicks(xtwo, collectionPlane);
809  shwtwoW = geom->WireCoordinate(vtx, collectionPlane);
810 
811  wirePitch = geom->WirePitch(hit->WireID());
812  tickToDist = detProp.DriftVelocity(detProp.Efield(), detProp.Temperature());
813  tickToDist *= 1.e-3 * sampling_rate(clockData); // 1e-3 is conversion of 1/us to 1/ns
814 
815  shwvtxx = shwvtxT * tickToDist;
816  shwvtxy = shwvtxW * wirePitch;
817 
818  shwtwox = shwtwoT * tickToDist;
819  shwtwoy = shwtwoW * wirePitch;
820 
821  xtwoorth = (shwtwoy - shwvtxy) + shwvtxx;
822  ytwoorth = -(shwtwox - shwvtxx) + shwvtxy;
823 
824  foundParent = true;
825  }
826  double xhit = hit->PeakTime() * tickToDist;
827  double yhit = hit->WireID().Wire * wirePitch;
828 
829  double ldist = abs((ytwoorth - shwvtxy) * xhit - (xtwoorth - shwvtxx) * yhit +
830  xtwoorth * shwvtxy - ytwoorth * shwvtxx) /
831  std::hypot(ytwoorth - shwvtxy, xtwoorth - shwvtxx);
832  double tdist = ((shwtwoy - shwvtxy) * xhit - (shwtwox - shwvtxx) * yhit + shwtwox * shwvtxy -
833  shwtwoy * shwvtxx) /
834  std::hypot(shwtwoy - shwvtxy, shwtwox - shwvtxx);
835 
836  double to3D = std::hypot(xvtx - xtwo, yvtx - ytwo, zvtx - ztwo) /
837  std::hypot(shwvtxx - shwtwox,
838  shwvtxy - shwtwoy); // distance between two points in 3D space is one
839  ldist *= to3D;
840  tdist *= to3D;
841  double t = ldist / X0;
842 
843  double energy = trackIDs[j].energy;
844 
845  ltemp->Fill(t, energy);
846  ttemp->Fill(tdist, energy);
847 
848  if (t < 1)
849  ttemp_1->Fill(tdist, energy);
850  else if (t < 2)
851  ttemp_2->Fill(tdist, energy);
852  else if (t < 3)
853  ttemp_3->Fill(tdist, energy);
854  else if (t < 4)
855  ttemp_4->Fill(tdist, energy);
856  else if (t < 5)
857  ttemp_5->Fill(tdist, energy);
858 
859  } // loop through track IDE
860 
861  } // loop through all hits
862 
863  for (int i = 0; i < LBINS; ++i) {
864  if (ltemp->GetBinContent(i + 1) == 0) continue;
865  fShowerProfileHitLong->Fill(ltemp->GetBinCenter(i + 1), ltemp->GetBinContent(i + 1));
866  fShowerProfileHitLong2D->Fill(ltemp->GetBinCenter(i + 1), elep, ltemp->GetBinContent(i + 1));
867  }
868 
869  for (int i = 0; i < TBINS; ++i) {
870  if (ttemp->GetBinContent(i + 1) == 0) continue;
871  fShowerProfileHitTrans->Fill(ttemp->GetBinCenter(i + 1), ttemp->GetBinContent(i + 1));
872  fShowerProfileHitTrans2D->Fill(ttemp->GetBinCenter(i + 1), elep, ttemp->GetBinContent(i + 1));
873 
875  ttemp_1->GetBinCenter(i + 1), elep, ttemp_1->GetBinContent(i + 1));
877  ttemp_2->GetBinCenter(i + 1), elep, ttemp_2->GetBinContent(i + 1));
879  ttemp_3->GetBinCenter(i + 1), elep, ttemp_3->GetBinContent(i + 1));
881  ttemp_4->GetBinCenter(i + 1), elep, ttemp_4->GetBinContent(i + 1));
883  ttemp_5->GetBinCenter(i + 1), elep, ttemp_5->GetBinContent(i + 1));
884  }
885 } // showerProfileTrue
886 
887 // -------------------------------------------------
888 
891  simb::MCParticle electron)
892 {
895 
896  std::vector<sim::MCEnDep> alledep;
897 
898  TH1F* ltemp = new TH1F("ltemp", "ltemp", LBINS, LMIN, LMAX);
899  TH1F* ttemp = new TH1F("ttemp", "ttemp", TBINS, TMIN, TMAX);
900 
901  TH1F* ttemp_1 = new TH1F("ttemp_1", "ttemp_1", TBINS, TMIN, TMAX);
902  TH1F* ttemp_2 = new TH1F("ttemp_2", "ttemp_2", TBINS, TMIN, TMAX);
903  TH1F* ttemp_3 = new TH1F("ttemp_3", "ttemp_3", TBINS, TMIN, TMAX);
904  TH1F* ttemp_4 = new TH1F("ttemp_4", "ttemp_4", TBINS, TMIN, TMAX);
905  TH1F* ttemp_5 = new TH1F("ttemp_5", "ttemp_5", TBINS, TMIN, TMAX);
906 
907  // get electron energy depositions
908  for (size_t i = 0; i < allchan.size(); ++i) {
909  art::Ptr<sim::SimChannel> simchan = allchan[i];
910  if (geom->View(simchan->Channel()) != geo::kV) continue;
911  auto tdc_ide_map = simchan->TDCIDEMap();
912 
913  for (auto const& tdc_ide_pair : tdc_ide_map) {
914  for (auto const& ide : tdc_ide_pair.second) {
915  if (piserv->TrackIdToMotherParticle_P(ide.trackID) == NULL) continue;
916  if (std::abs(piserv->TrackIdToMotherParticle_P(ide.trackID)->PdgCode()) != 11) continue;
917 
919  edep.SetVertex(ide.x, ide.y, ide.z);
920  edep.SetEnergy(ide.energy);
921  edep.SetTrackId(ide.trackID);
922 
923  alledep.push_back(edep);
924 
925  } // loop through ide
926  } // loop through tdc_ide
927 
928  } // loop through channels
929 
930  double x0 = electron.Vx();
931  double y0 = electron.Vy();
932  double z0 = electron.Vz();
933 
934  double x2 = electron.Px();
935  double y2 = electron.Py();
936  double z2 = electron.Pz();
937 
938  TVector3 v0(x2, y2, z2);
939  v0 = v0.Unit();
940 
941  for (size_t i = 0; i < alledep.size(); ++i) {
942  double x = (double)alledep[i].Vertex()[0];
943  double y = (double)alledep[i].Vertex()[1];
944  double z = (double)alledep[i].Vertex()[2];
945 
946  TVector3 v1(x - x0, y - y0, z - z0);
947 
948  double ldist = v0.Dot(v1);
949  double t = ldist / X0;
950  double tdist = (v0.Orthogonal()).Dot(v1);
951 
952  double energy = alledep[i].Energy();
953 
954  ltemp->Fill(t, energy);
955  ttemp->Fill(tdist, energy);
956 
957  if (t < 1)
958  ttemp_1->Fill(tdist, energy);
959  else if (t < 2)
960  ttemp_2->Fill(tdist, energy);
961  else if (t < 3)
962  ttemp_3->Fill(tdist, energy);
963  else if (t < 4)
964  ttemp_4->Fill(tdist, energy);
965  else if (t < 5)
966  ttemp_5->Fill(tdist, energy);
967  }
968 
969  for (int i = 0; i < LBINS; ++i) {
970  if (ltemp->GetBinContent(i + 1) == 0) continue;
971  fShowerProfileSimLong->Fill(ltemp->GetBinCenter(i + 1), ltemp->GetBinContent(i + 1));
973  ltemp->GetBinCenter(i + 1), electron.E(), ltemp->GetBinContent(i + 1));
974  }
975 
976  for (int i = 0; i < TBINS; ++i) {
977  if (ttemp->GetBinContent(i + 1) == 0) continue;
978  fShowerProfileSimTrans->Fill(ttemp->GetBinCenter(i + 1), ttemp->GetBinContent(i + 1));
980  ttemp->GetBinCenter(i + 1), electron.E(), ttemp->GetBinContent(i + 1));
981 
983  ttemp_1->GetBinCenter(i + 1), electron.E(), ttemp_1->GetBinContent(i + 1));
985  ttemp_2->GetBinCenter(i + 1), electron.E(), ttemp_2->GetBinContent(i + 1));
987  ttemp_3->GetBinCenter(i + 1), electron.E(), ttemp_3->GetBinContent(i + 1));
989  ttemp_4->GetBinCenter(i + 1), electron.E(), ttemp_4->GetBinContent(i + 1));
991  ttemp_5->GetBinCenter(i + 1), electron.E(), ttemp_5->GetBinContent(i + 1));
992  }
993 } // showerProfileTrue
994 
995 // -------------------------------------------------
996 
Float_t x
Definition: compare.C:6
double E(const int i=0) const
Definition: MCParticle.h:234
TCShowerTemplateMaker(fhicl::ParameterSet const &pset)
Length_t WireCoordinate(Point_t const &pos, PlaneID const &planeid) const
Returns the index of the nearest wire to the specified position.
int PdgCode() const
Definition: MCParticle.h:213
int CCNC() const
Definition: MCNeutrino.h:148
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:77
double Py(const int i=0) const
Definition: MCParticle.h:232
void SetEnergy(float e)
Definition: MCDataHolder.h:36
void analyze(const art::Event &evt) override
Planes which measure V.
Definition: geo_types.h:136
const simb::MCParticle * TrackIdToParticle_P(int id) const
Declaration of signal hit object.
Float_t y
Definition: compare.C:6
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
Double_t z
Definition: plot.C:276
double Px(const int i=0) const
Definition: MCParticle.h:231
double Temperature() const
In kelvin.
::geo::Point_t toPoint(Point const &p)
Convert the specified point into a geo::Point_t.
constexpr auto abs(T v)
Returns the absolute value of the argument.
STL namespace.
void SetTrackId(unsigned int id)
Definition: MCDataHolder.h:38
void showerProfile(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, std::vector< art::Ptr< recob::Hit >> showerhits, TVector3 shwvtx, TVector3 shwdir, double elep)
EDAnalyzer(fhicl::ParameterSet const &pset)
Definition: EDAnalyzer.cc:6
Float_t y2[n_points_geant4]
Definition: compare.C:26
double Efield(unsigned int planegap=0) const
kV/cm
IDparameter< geo::PlaneID > PlaneID
Member type of validated geo::PlaneID parameter.
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:289
void showerProfileTrue(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, std::vector< art::Ptr< recob::Hit >> allhits, double elep)
const simb::MCParticle & Lepton() const
Definition: MCNeutrino.h:147
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
parameter set interface
double energy
Definition: plottest35.C:25
void SetVertex(float x, float y, float z)
Definition: MCDataHolder.h:29
double ConvertXToTicks(double X, int p, int t, int c) const
double DriftVelocity(double efield=0., double temperature=0.) const
cm/us
Double_t edep
Definition: macro.C:13
Detector simulation of raw signals on wires.
raw::ChannelID_t Channel() const
Returns the readout channel this object describes.
Definition: SimChannel.h:323
double Vx(const int i=0) const
Definition: MCParticle.h:222
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
decltype(auto) get(T &&obj)
ADL-aware version of std::to_string.
Definition: StdUtils.h:120
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Contains all timing reference information for the detector.
View_t View(PlaneID const &pid) const
Returns the view (wire orientation) on the channels of specified TPC plane.
std::vector< sim::TrackIDE > HitToEveTrackIDEs(detinfo::DetectorClocksData const &clockData, recob::Hit const &hit) const
double Pz(const int i=0) const
Definition: MCParticle.h:233
object containing MC truth information necessary for making RawDigits and doing back tracking ...
double Vz(const int i=0) const
Definition: MCParticle.h:224
TDCIDEs_t const & TDCIDEMap() const
Returns all the deposited energy information as stored.
Definition: SimChannel.h:319
boost::graph_traits< ModuleGraph >::vertex_descriptor Vertex
Definition: ModuleGraph.h:25
Direction
Definition: types.h:12
bool NeutrinoSet() const
Definition: MCTruth.h:78
TCEvent evt
Definition: DataStructs.cxx:8
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:306
const simb::MCParticle * TrackIdToMotherParticle_P(int id) const
double LifetimeCorrection(detinfo::DetectorClocksData const &clock_data, detinfo::DetectorPropertiesData const &det_prop, double time, double T0=0) const
Float_t x2[n_points_geant4]
Definition: compare.C:26
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
Length_t WirePitch(PlaneID const &planeid=plane_zero) const
Returns the distance between two consecutive wires.
double Vy(const int i=0) const
Definition: MCParticle.h:223
art framework interface to geometry description
decltype(auto) constexpr empty(T &&obj)
ADL-aware version of std::empty.
Definition: StdUtils.h:109