LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
PointIdEffTest_module.cc
Go to the documentation of this file.
1 // Class: PointIdEffTest
3 // Module Type: analyzer
4 // File: PointIdEffTest_module.cc
5 //
6 // Author: dorota.stefan@cern.ch
7 //
8 // Generated at Fri Apr 29 06:42:27 2016 by Dorota Stefan using artmod
9 // from cetpkgsupport v1_10_01.
11 
27 
35 #include "art_root_io/TFileService.h"
39 #include "cetlib_except/exception.h"
40 #include "fhiclcpp/types/Atom.h"
41 #include "fhiclcpp/types/Comment.h"
42 #include "fhiclcpp/types/Name.h"
43 #include "fhiclcpp/types/Table.h"
45 
46 #include "TH1.h"
47 #include "TTree.h"
48 
49 #include <array>
50 #include <cmath>
51 #include <fstream>
52 #include <iostream>
53 #include <unordered_map>
54 #include <vector>
55 
56 #define MVA_LENGTH 4
57 
58 namespace nnet {
59 
61  public:
62  enum EId { kShower = 0, kTrack = 1, kMichel = 2 };
63 
64  struct Config {
65  using Name = fhicl::Name;
67 
69  Name("CalorimetryAlg"),
70  Comment("Used to calculate electron lifetime correction.")};
71 
73  Comment("Simulation producer")};
74 
76  Name("PfpModuleLabel"),
77  Comment("PFP producer tag, to compare with NNet results")};
78 
80  Comment("NNet outputs tag")};
81 
82  fhicl::Atom<bool> SaveHitsFile{Name("SaveHitsFile"), Comment("Dump hits info to text file")};
83 
84  fhicl::Atom<unsigned int> View{Name("View"), Comment("Which view is evaluated")};
85  };
87 
88  explicit PointIdEffTest(Parameters const& config);
89 
90  PointIdEffTest(PointIdEffTest const&) = delete;
91  PointIdEffTest(PointIdEffTest&&) = delete;
92  PointIdEffTest& operator=(PointIdEffTest const&) = delete;
94 
95  private:
96  void beginRun(const art::Run& run) override;
97  void beginJob() override;
98  void endJob() override;
99  void analyze(art::Event const& e) override;
100 
101  void cleanup();
102 
103  void countTruthDep(const std::vector<sim::SimChannel>& channels,
104  float& emLike,
105  float& trackLike) const;
106 
107  void countPfpDep(detinfo::DetectorClocksData const& clockData,
108  detinfo::DetectorPropertiesData const& detProp,
109  const std::vector<recob::PFParticle>& pfparticles,
110  const art::FindManyP<recob::Cluster>& pfpclus,
111  const art::FindManyP<recob::Hit>& cluhits,
112  float& emLike,
113  float& trackLike) const;
114 
115  bool isMuonDecaying(const simb::MCParticle& particle,
116  const std::unordered_map<int, const simb::MCParticle*>& particleMap) const;
117 
118  int testCNN(detinfo::DetectorClocksData const& clockData,
119  detinfo::DetectorPropertiesData const& detProp,
120  const std::vector<sim::SimChannel>& channels,
122  const std::array<float, MVA_LENGTH>& cnn_out,
124  size_t cidx);
125 
126  int fRun, fEvent;
130  float fHitEM_mc, fpEM;
135 
136  //const int kEffSize = 100;
140  float fHitRecoEM[100], fHitRecoFractionEM[100];
141 
142  int fMcPid;
143  int fClSize;
144  float fPidValue; // P(track-like)
145 
146  int fTrkOk[100], fTrkBad[100];
147  int fShOk[100], fShBad[100];
148  int fNone, fTotal;
149 
151 
153 
155 
156  std::ofstream fHitsOutFile;
157 
158  unsigned int fView;
159 
160  std::unordered_map<int, const simb::MCParticle*> fParticleMap;
161 
167  };
168 
169 } // namespace nnet
170 
172  : art::EDAnalyzer(config)
173  , fMcPid(-1)
174  , fClSize(0)
175  , fTrkLikeIdx(-1)
176  , fEmLikeIdx(-1)
177  , fNoneIdx(-1)
178  , fMichelLikeIdx(-1)
179  , fView(config().View())
180  , fCalorimetryAlg(config().CalorimetryAlg())
182  , fPfpModuleLabel(config().PfpModuleLabel())
183  , fNNetModuleLabel(config().NNetModuleLabel())
184  , fSaveHitsFile(config().SaveHitsFile())
185 {}
186 
188 {
190  fElectronsToGeV = 1. / larParameters->GeVToElectrons();
191 }
192 
194 {
196 
197  fEventTree = tfs->make<TTree>("event", "event info");
198  fEventTree->Branch("fRun", &fRun, "fRun/I");
199  fEventTree->Branch("fEvent", &fEvent, "fEvent/I");
200  fEventTree->Branch("fMcDepEM", &fMcDepEM, "fMcDepEM/F");
201  fEventTree->Branch("fMcDepTrack", &fMcDepTrack, "fMcDepTrack/F");
202  fEventTree->Branch("fMcFractionEM", &fMcFractionEM, "fMcFractionEM/F");
203  fEventTree->Branch("fPfpDepEM", &fPfpDepEM, "fPfpDepEM/F");
204  fEventTree->Branch("fPfpDepTrack", &fPfpDepTrack, "fPfpDepTrack/F");
205  fEventTree->Branch("fHitEM_0p5", &fHitEM_0p5, "fHitEM_0p5/F");
206  fEventTree->Branch("fHitMichel_0p5", &fHitMichel_0p5, "fHitMichel_0p5/F");
207  fEventTree->Branch("fHitTrack_0p5", &fHitTrack_0p5, "fHitTrack_0p5/F");
208  fEventTree->Branch("fHitMcFractionEM", &fHitMcFractionEM, "fHitMcFractionEM/F");
209  fEventTree->Branch("fHitEM_0p85", &fHitEM_0p85, "fHitEM_0p85/F");
210  fEventTree->Branch("fHitTrack_0p85", &fHitTrack_0p85, "fHitTrack_0p85/F");
211  fEventTree->Branch("fCleanHit", &fCleanHit, "fCleanHit/F");
212 
213  fEventTree->Branch("fHitsEM_OK_0p5", fHitsEM_OK_0p5, "fHitsEM_OK_0p5[100]/F");
214  fEventTree->Branch("fHitsTrack_OK_0p5", fHitsTrack_OK_0p5, "fHitsTrack_OK_0p5[100]/F");
215  fEventTree->Branch("fHitsMichel_OK_0p5", fHitsMichel_OK_0p5, "fHitsMichel_OK_0p5[100]/F");
216  fEventTree->Branch(
217  "fHitsMichel_False_0p5", fHitsMichel_False_0p5, "fHitsMichel_False_0p5[100]/F");
218 
219  fEventTree->Branch("fHitsEM_OK_0p85", fHitsEM_OK_0p85, "fHitsEM_OK_0p85[100]/F");
220  fEventTree->Branch("fHitsTrack_OK_0p85", fHitsTrack_OK_0p85, "fHitsTrack_OK_0p85[100]/F");
221 
222  fEventTree->Branch("fHitRecoEM", fHitRecoEM, "fHitRecoEM[100]/F");
223  fEventTree->Branch("fHitRecoFractionEM", fHitRecoFractionEM, "fHitRecoFractionEM[100]/F");
224 
225  fClusterTree = tfs->make<TTree>("cluster", "clusters info");
226  fClusterTree->Branch("fMcPid", &fMcPid, "fMcPid/I");
227  fClusterTree->Branch("fClSize", &fClSize, "fClSize/I");
228  fClusterTree->Branch("fPidValue", &fPidValue, "fPidValue/F");
229  fClusterTree->Branch("fpMichel_cl", &fpMichel_cl, "fpMichel_cl/F");
230 
231  fHitTree = tfs->make<TTree>("hits", "hits info");
232  fHitTree->Branch("fRun", &fRun, "fRun/I");
233  fHitTree->Branch("fEvent", &fEvent, "fEvent/I");
234  fHitTree->Branch("fHitEM_mc", &fHitEM_mc, "fHitEM_mc/F");
235  fHitTree->Branch("fpEM", &fpEM, "fpEM/F");
236  fHitTree->Branch("fPidValue", &fPidValue, "fPidValue/F");
237  fHitTree->Branch("fHitMichel_mc", &fHitMichel_mc, "fHitMichel_mc/F");
238  fHitTree->Branch("fpMichel_hit", &fpMichel_hit, "fpMichel_hit/F");
239  fHitTree->Branch("fOutTrk", &fOutTrk, "fOutTrk/F");
240  fHitTree->Branch("fOutEM", &fOutEM, "fOutEM/F");
241  fHitTree->Branch("fOutNone", &fOutNone, "fOutNone/F");
242 
243  if (fSaveHitsFile) fHitsOutFile.open("hits_pid.prn");
244 
245  fNone = 0;
246  fTotal = 0;
247  for (size_t i = 0; i < 100; ++i) {
248  fShOk[i] = 0;
249  fShBad[i] = 0;
250  fTrkOk[i] = 0;
251  fTrkBad[i] = 0;
252  }
253 }
254 
256 {
257  if (fSaveHitsFile) fHitsOutFile.close();
258 
260 
261  TTree* thrTree = tfs->make<TTree>("threshold", "error rate vs threshold");
262 
263  float thr, shErr, trkErr;
264  thrTree->Branch("thr", &thr, "thr/F");
265  thrTree->Branch("shErr", &shErr, "shErr/F");
266  thrTree->Branch("trkErr", &trkErr, "trkErr/F");
267 
268  for (size_t i = 0; i < 100; ++i) {
269  thr = 0.01 * i;
270  shErr = fShBad[i] / float(fShBad[i] + fShOk[i]);
271  trkErr = fTrkBad[i] / float(fTrkBad[i] + fTrkOk[i]);
272  thrTree->Fill();
273  }
274 }
275 
277 {
278  fParticleMap.clear();
279 
280  fMcDepEM = 0;
281  fMcDepTrack = 0;
282  fMcFractionEM = 0;
283  fPfpDepEM = 0;
284  fPfpDepTrack = 0;
285  fHitEM_0p5 = 0;
286  fHitTrack_0p5 = 0;
287  fHitEM_0p85 = 0;
288  fHitTrack_0p85 = 0;
289  fHitMichel_0p5 = 0;
290  fHitMcFractionEM = 0;
291  fTotHit = 0;
292  fCleanHit = 0;
293 
294  for (size_t i = 0; i < 100; ++i) {
295  fHitsEM_OK_0p5[i] = 0;
296  fHitsTrack_OK_0p5[i] = 0;
297  fHitsEM_OK_0p85[i] = 0;
298  fHitsTrack_OK_0p85[i] = 0;
299  fHitsMichel_OK_0p5[i] = 0;
300  fHitsMichel_False_0p5[i] = 0;
301  fHitRecoEM[i] = 0;
302  fHitRecoFractionEM[i] = 0;
303  }
304 
305  fTrkLikeIdx = -1;
306  fEmLikeIdx = -1;
307  fNoneIdx = -1;
308  fMichelLikeIdx = -1;
309 }
310 
312 {
313  cleanup(); // remove everything from members
314 
315  fRun = e.run();
316  fEvent = e.id().event();
317 
318  // access to MC information
319 
320  // MC particles list
321  auto particleHandle = e.getValidHandle<std::vector<simb::MCParticle>>(fSimulationProducerLabel);
322  for (auto const& particle : *particleHandle) {
323  fParticleMap[particle.TrackId()] = &particle;
324  }
325 
326  // SimChannels
327  auto simChannelHandle = e.getValidHandle<std::vector<sim::SimChannel>>(fSimulationProducerLabel);
328  countTruthDep(*simChannelHandle, fMcDepEM, fMcDepTrack);
329 
330  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(e);
331  auto const detProp =
333 
334  // PFParticle selection results
336  if (e.getByLabel(fPfpModuleLabel, pfpHandle)) {
337  auto cluHandle = e.getValidHandle<std::vector<recob::Cluster>>(fPfpModuleLabel);
338  const art::FindManyP<recob::Cluster> clusFromPfps(pfpHandle, e, fPfpModuleLabel);
339  const art::FindManyP<recob::Hit> hitsFromClus(cluHandle, e, fPfpModuleLabel);
340  countPfpDep(
341  clockData, detProp, *pfpHandle, clusFromPfps, hitsFromClus, fPfpDepEM, fPfpDepTrack);
342  }
343 
344  // output from cnn's
345 
347  e, fNNetModuleLabel); // hit-by-hit outpus just to be dumped to file for debugging
348  fTrkLikeIdx = hitResults.getIndex("track");
349  fEmLikeIdx = hitResults.getIndex("em");
350  fNoneIdx = hitResults.getIndex("none");
351  fMichelLikeIdx = hitResults.getIndex("michel");
352  if ((fTrkLikeIdx < 0) || (fEmLikeIdx < 0)) {
353  throw cet::exception("PointIdEffTest")
354  << "No em/track labeled columns in MVA data products." << std::endl;
355  }
356 
358  e, fNNetModuleLabel); // outputs for clusters recovered in not-throwing way
359  if (cluResults) {
360  const art::FindManyP<recob::Hit> hitsFromClusters(
361  cluResults->dataHandle(), e, cluResults->dataTag());
362 
363  for (size_t c = 0; c < cluResults->size(); ++c) {
364  const recob::Cluster& clu = cluResults->item(c);
365 
366  if (clu.Plane().Plane != fView) continue;
367 
368  const std::vector<art::Ptr<recob::Hit>>& hits = hitsFromClusters.at(c);
369  std::array<float, MVA_LENGTH> cnn_out = cluResults->getOutput(c);
370 
371  testCNN(clockData,
372  detProp,
373  *simChannelHandle,
374  hits,
375  cnn_out,
376  hitResults.outputs(),
377  c); // test hits in the cluster
378  }
379 
380  if (fTotHit > 0)
382  else
383  fCleanHit = 0;
384 
385  double totMcDep = fMcDepEM + fMcDepTrack;
386  if (totMcDep)
387  fMcFractionEM = fMcDepEM / totMcDep;
388  else
389  fMcFractionEM = 0;
390 
391  double totEmTrk0p5 = fHitEM_0p5 + fHitTrack_0p5;
392  if (totEmTrk0p5 > 0)
393  fHitMcFractionEM = fHitEM_0p5 / totEmTrk0p5;
394  else
395  fHitMcFractionEM = 0;
396 
397  for (size_t i = 0; i < 100; ++i) {
398  if (fHitEM_0p5 > 0)
400  else
401  fHitsEM_OK_0p5[i] = 0;
402 
403  if (fHitTrack_0p5 > 0)
405  else
406  fHitsTrack_OK_0p5[i] = 0;
407 
408  if (fHitEM_0p85 > 0)
410  else
411  fHitsEM_OK_0p85[i] = 0;
412 
413  if (fHitTrack_0p85 > 0)
415  else
416  fHitsTrack_OK_0p85[i] = 0;
417 
418  if (totEmTrk0p5 > 0)
419  fHitRecoFractionEM[i] = fHitRecoEM[i] / totEmTrk0p5;
420  else
421  fHitRecoFractionEM[i] = 0;
422  }
423 
424  fEventTree->Fill();
425  }
426  else {
427  mf::LogWarning("TrainingDataAlg") << "MVA FOR CLUSTERS NOT FOUND";
428  }
429 
430  cleanup(); // remove everything from members
431 }
432 /******************************************/
433 
434 void nnet::PointIdEffTest::countTruthDep(const std::vector<sim::SimChannel>& channels,
435  float& emLike,
436  float& trackLike) const
437 {
438  emLike = 0;
439  trackLike = 0;
440  for (auto const& channel : channels) {
441  // for every time slice in this channel:
442  auto const& timeSlices = channel.TDCIDEMap();
443  for (auto const& timeSlice : timeSlices) {
444  // loop over the energy deposits.
445  auto const& energyDeposits = timeSlice.second;
446  for (auto const& energyDeposit : energyDeposits) {
447  int trackID = energyDeposit.trackID;
448 
449  double energy = energyDeposit.numElectrons * fElectronsToGeV * 1000;
450 
451  if (trackID < 0) { emLike += energy; }
452  else if (trackID > 0) {
453  auto search = fParticleMap.find(trackID);
454  bool found = true;
455  if (search == fParticleMap.end()) {
456  mf::LogWarning("TrainingDataAlg") << "PARTICLE NOT FOUND";
457  found = false;
458  }
459 
460  int pdg = 0;
461  if (found) {
462  const simb::MCParticle& particle = *((*search).second);
463  if (!pdg) pdg = particle.PdgCode(); // not EM activity so read what PDG it is
464  }
465 
466  if ((pdg == 11) || (pdg == -11) || (pdg == 22))
467  emLike += energy;
468  else
469  trackLike += energy;
470  }
471  }
472  }
473  }
474 }
475 /******************************************/
476 
478  detinfo::DetectorPropertiesData const& detProp,
479  const std::vector<recob::PFParticle>& pfparticles,
480  const art::FindManyP<recob::Cluster>& pfpclus,
481  const art::FindManyP<recob::Hit>& cluhits,
482  float& emLike,
483  float& trackLike) const
484 {
485  emLike = 0;
486  trackLike = 0;
487  for (size_t i = 0; i < pfparticles.size(); ++i) {
488  const auto& pfp = pfparticles[i];
489  const auto& clus = pfpclus.at(i);
490 
491  float hitdep = 0;
492  for (const auto& c : clus) {
493  const auto& hits = cluhits.at(c.key());
494  for (const auto& h : hits) {
495  if (h->View() == fView) {
496  hitdep +=
497  h->SummedADC() * fCalorimetryAlg.LifetimeCorrection(clockData, detProp, h->PeakTime());
498  }
499  }
500  }
501 
502  if ((pfp.PdgCode() == 11) || pfp.PdgCode() == -11) { emLike += hitdep; }
503  else {
504  trackLike += hitdep;
505  }
506  }
507 }
508 /******************************************/
509 
511  const simb::MCParticle& particle,
512  const std::unordered_map<int, const simb::MCParticle*>& particleMap) const
513 {
514  bool hasElectron = false, hasNuMu = false, hasNuE = false;
515 
516  int pdg = abs(particle.PdgCode());
517  if ((pdg == 13) && (particle.EndProcess() == "FastScintillation")) // potential muon decay at rest
518  {
519  unsigned int nSec = particle.NumberDaughters();
520  for (size_t d = 0; d < nSec; ++d) {
521  auto d_search = particleMap.find(particle.Daughter(d));
522  if (d_search != particleMap.end()) {
523  auto const& daughter = *((*d_search).second);
524  int d_pdg = abs(daughter.PdgCode());
525  if (d_pdg == 11)
526  hasElectron = true;
527  else if (d_pdg == 14)
528  hasNuMu = true;
529  else if (d_pdg == 12)
530  hasNuE = true;
531  }
532  }
533  }
534 
535  return (hasElectron && hasNuMu && hasNuE);
536 }
537 /******************************************/
538 
540  detinfo::DetectorPropertiesData const& detProp,
541  const std::vector<sim::SimChannel>& channels,
543  const std::array<float, MVA_LENGTH>& cnn_out,
545  size_t cidx)
546 {
547  fClSize = hits.size();
548 
549  std::unordered_map<int, int> mcHitPid;
550 
551  fPidValue = 0;
552  double p_trk_or_sh = cnn_out[fTrkLikeIdx] + cnn_out[fEmLikeIdx];
553  if (p_trk_or_sh > 0) { fPidValue = cnn_out[fTrkLikeIdx] / p_trk_or_sh; }
554 
555  double p_michel = 0;
556  if (fMichelLikeIdx >= 0) { fpMichel_cl = cnn_out[fMichelLikeIdx]; }
557 
558  double totEnSh = 0, totEnTrk = 0;
559  for (auto const& hit : hits) {
560  // the channel associated with this hit.
561  auto hitChannelNumber = hit->Channel();
562 
563  double hitEn = 0, hitEnSh = 0, hitEnTrk = 0, hitEnMichel = 0;
564 
565  auto const& vout = hit_outs[hit.key()];
566  fOutTrk = vout[fTrkLikeIdx];
567  fOutEM = vout[fEmLikeIdx];
568  if (fNoneIdx >= 0) { fOutNone = vout[fNoneIdx]; }
569  if (fMichelLikeIdx >= 0) { p_michel = vout[fMichelLikeIdx]; }
570 
571  for (auto const& channel : channels) {
572  if (channel.Channel() != hitChannelNumber) continue;
573 
574  // for every time slice in this channel:
575  auto const& timeSlices = channel.TDCIDEMap();
576  for (auto const& timeSlice : timeSlices) {
577  int time = timeSlice.first;
578  if (std::abs(hit->TimeDistanceAsRMS(time)) < 1.0) {
579  // loop over the energy deposits.
580  auto const& energyDeposits = timeSlice.second;
581 
582  for (auto const& energyDeposit : energyDeposits) {
583  int trackID = energyDeposit.trackID;
584 
585  double energy = energyDeposit.numElectrons * fElectronsToGeV * 1000;
586  hitEn += energy;
587 
588  if (trackID < 0) { hitEnSh += energy; } // EM activity
589  else if (trackID > 0) {
590  auto search = fParticleMap.find(trackID);
591  if (search != fParticleMap.end()) {
592  const simb::MCParticle& particle = *((*search).second);
593  int pdg = particle.PdgCode(); // not EM activity so read what PDG it is
594 
595  if ((pdg == 11) || (pdg == -11) || (pdg == 22))
596  hitEnSh += energy;
597  else
598  hitEnTrk += energy;
599 
600  if (pdg == 11) // electron, check if it is Michel
601  {
602  auto msearch = fParticleMap.find(particle.Mother());
603  if (msearch != fParticleMap.end()) {
604  auto const& mother = *((*msearch).second);
605  if (isMuonDecaying(mother, fParticleMap)) { hitEnMichel += energy; }
606  }
607  }
608  }
609  else {
610  mf::LogWarning("TrainingDataAlg") << "PARTICLE NOT FOUND";
611  }
612  }
613  }
614  }
615  }
616  }
617  totEnSh += hitEnSh;
618  totEnTrk += hitEnTrk;
619 
620  double hitAdc =
621  hit->SummedADC() * fCalorimetryAlg.LifetimeCorrection(clockData, detProp, hit->PeakTime());
622  fTotHit += hitAdc;
623 
624  int hitPidMc_0p5 = -1;
625  if (hitEnSh > hitEnTrk) {
626  fHitEM_0p5 += hitAdc;
627  hitPidMc_0p5 = nnet::PointIdEffTest::kShower;
628  }
629  else {
630  fHitTrack_0p5 += hitAdc;
631  hitPidMc_0p5 = nnet::PointIdEffTest::kTrack;
632  }
633  mcHitPid[hit.key()] = hitPidMc_0p5;
634  auto const& hout = hit_outs[hit.key()];
635  fpEM = 0;
636  float hit_trk_or_sh = hout[fTrkLikeIdx] + hout[fEmLikeIdx];
637  if (hit_trk_or_sh > 0) fpEM = hout[fEmLikeIdx] / hit_trk_or_sh;
638  fHitEM_mc = hitEnSh / (hitEnSh + hitEnTrk);
639 
640  int hitPidMc_0p85 = -1;
641  double hitDep = hitEnSh + hitEnTrk;
642  if (hitEnSh > 0.85 * hitDep) {
643  fHitEM_0p85 += hitAdc;
644  fCleanHit += hitAdc;
645  hitPidMc_0p85 = nnet::PointIdEffTest::kShower;
646  }
647  else if (hitEnTrk > 0.85 * hitDep) {
648  fHitTrack_0p85 += hitAdc;
649  fCleanHit += hitAdc;
650  hitPidMc_0p85 = nnet::PointIdEffTest::kTrack;
651  }
652 
653  bool mcMichel = false;
654  fpMichel_hit = p_michel;
655  fHitMichel_mc = hitEnMichel / hitEn;
656  if (fHitMichel_mc > 0.5) {
657  fHitMichel_0p5 += hitAdc;
658  mcMichel = true;
659  }
660 
661  fHitTree->Fill();
662 
663  for (size_t i = 0; i < 100; ++i) {
664  double thr = 0.01 * i;
665 
666  if (p_michel > thr) {
667  if (mcMichel) { fHitsMichel_OK_0p5[i] += hitAdc; }
668  else {
669  fHitsMichel_False_0p5[i] += hitAdc;
670  }
671  }
672 
673  int recoPid = -1;
674  if (fPidValue < thr) {
676  fHitRecoEM[i] += hitAdc;
677  }
678  else
680 
681  if ((recoPid == nnet::PointIdEffTest::kShower) &&
682  (hitPidMc_0p5 == nnet::PointIdEffTest::kShower)) {
683  fHitsEM_OK_0p5[i] += hitAdc;
684  }
685  else if ((recoPid == nnet::PointIdEffTest::kTrack) &&
686  (hitPidMc_0p5 == nnet::PointIdEffTest::kTrack)) {
687  fHitsTrack_OK_0p5[i] += hitAdc;
688  }
689 
690  if ((recoPid == nnet::PointIdEffTest::kShower) &&
691  (hitPidMc_0p85 == nnet::PointIdEffTest::kShower)) {
692  fHitsEM_OK_0p85[i] += hitAdc;
693  }
694  else if ((recoPid == nnet::PointIdEffTest::kTrack) &&
695  (hitPidMc_0p85 == nnet::PointIdEffTest::kTrack)) {
696  fHitsTrack_OK_0p85[i] += hitAdc;
697  }
698  }
699  }
700 
701  // ************ count clusters *************
702 
703  fMcPid = -1;
704  if (totEnSh > 1.5 * totEnTrk) // major energy deposit from EM activity
705  {
707  }
708  else if (totEnTrk > 1.5 * totEnSh) {
710  }
711 
712  for (size_t i = 0; i < 100; ++i) {
713  double thr = 0.01 * i;
714 
715  int recoPid = -1;
716  if (fPidValue < thr)
718  else
720 
722  fShOk[i] += fClSize;
723  }
724  else if ((recoPid == nnet::PointIdEffTest::kTrack) &&
726  fTrkOk[i] += fClSize;
727  }
728  else if ((recoPid == nnet::PointIdEffTest::kShower) &&
730  fTrkBad[i] += fClSize;
731  }
732  else if ((recoPid == nnet::PointIdEffTest::kTrack) &&
734  fShBad[i] += fClSize;
735  }
736  else {
737  fNone++;
738  }
739  }
740  fTotal++;
741 
742  if (fSaveHitsFile) {
743  for (auto const& h : hits) {
744  auto const& vout = hit_outs[h.key()];
745  double hitPidValue = 0;
746  double h_trk_or_sh = vout[fTrkLikeIdx] + vout[fEmLikeIdx];
747  if (h_trk_or_sh > 0) hitPidValue = vout[fTrkLikeIdx] / h_trk_or_sh;
748 
749  fHitsOutFile << fRun << " " << fEvent << " " << h->WireID().TPC << " " << h->WireID().Wire
750  << " " << h->PeakTime() << " "
751  << h->SummedADC() *
752  fCalorimetryAlg.LifetimeCorrection(clockData, detProp, h->PeakTime())
753  << " " << mcHitPid[h.key()] << " " << fPidValue << " " << hitPidValue;
754 
755  if (fMichelLikeIdx >= 0) {
756  fHitsOutFile << " " << vout[fMichelLikeIdx]; // is michel?
757  }
758 
759  fHitsOutFile << " " << cidx << std::endl;
760  }
761  }
762 
763  fClusterTree->Fill();
764  return fMcPid;
765 }
766 /******************************************/
767 
art::InputTag fSimulationProducerLabel
Store parameters for running LArG4.
int PdgCode() const
Definition: MCParticle.h:213
int Mother() const
Definition: MCParticle.h:214
Declaration of signal hit object.
constexpr auto abs(T v)
Returns the absolute value of the argument.
int getIndex(const std::string &name) const
Index of column with given name, or -1 if name not found.
Definition: MVAReader.h:82
void countPfpDep(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const std::vector< recob::PFParticle > &pfparticles, const art::FindManyP< recob::Cluster > &pfpclus, const art::FindManyP< recob::Hit > &cluhits, float &emLike, float &trackLike) const
Set of hits with a 2D structure.
Definition: Cluster.h:69
void countTruthDep(const std::vector< sim::SimChannel > &channels, float &emLike, float &trackLike) const
geo::PlaneID Plane() const
Returns the plane ID this cluster lies on.
Definition: Cluster.h:717
PointIdEffTest & operator=(PointIdEffTest const &)=delete
EDAnalyzer(fhicl::ParameterSet const &pset)
Definition: EDAnalyzer.cc:6
Particle class.
int NumberDaughters() const
Definition: MCParticle.h:218
Definition: Run.h:37
int Daughter(const int i) const
Definition: MCParticle.cxx:118
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:289
void analyze(art::Event const &e) override
std::unordered_map< int, const simb::MCParticle * > fParticleMap
PointIdEffTest(Parameters const &config)
void hits()
Definition: readHits.C:15
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
void beginRun(const art::Run &run) override
std::string EndProcess() const
Definition: MCParticle.h:217
int testCNN(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const std::vector< sim::SimChannel > &channels, const std::vector< art::Ptr< recob::Hit >> &hits, const std::array< float, MVA_LENGTH > &cnn_out, const std::vector< anab::FeatureVector< MVA_LENGTH >> &hit_outs, size_t cidx)
bool isMuonDecaying(const simb::MCParticle &particle, const std::unordered_map< int, const simb::MCParticle * > &particleMap) const
Definition: EmTrack.h:40
double energy
Definition: plottest35.C:25
Float_t d
Definition: plot.C:235
std::vector< FeatureVector< N > > const & outputs() const
Access the vector of the feature vectors.
Definition: MVAReader.h:124
Provides recob::Track data product.
fhicl::Atom< art::InputTag > PfpModuleLabel
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:481
fhicl::Table< calo::CalorimetryAlg::Config > CalorimetryAlg
Declaration of cluster object.
Detector simulation of raw signals on wires.
fhicl::Atom< unsigned int > View
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
fhicl::Atom< art::InputTag > SimModuleLabel
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Contains all timing reference information for the detector.
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
object containing MC truth information necessary for making RawDigits and doing back tracking ...
Definition: MVAAlg.h:12
fhicl::Atom< art::InputTag > NNetModuleLabel
EventNumber_t event() const
Definition: EventID.h:116
Declaration of basic channel signal object.
double LifetimeCorrection(detinfo::DetectorClocksData const &clock_data, detinfo::DetectorPropertiesData const &det_prop, double time, double T0=0) const
Collection of Physical constants used in LArSoft.
Float_t e
Definition: plot.C:35
RunNumber_t run() const
Definition: Event.cc:29
calo::CalorimetryAlg fCalorimetryAlg
static std::unique_ptr< MVAReader > create(const art::Event &evt, const art::InputTag &tag)
Definition: MVAReader.h:108
EventID id() const
Definition: Event.cc:23
double GeVToElectrons() const
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33