LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
ClusteringValidation_module.cc
Go to the documentation of this file.
1 // Class: ClusteringValidation
3 // Module type: analyser
4 // File: ClusteringValidation_module.cc
5 // Author: Mike Wallbank (m.wallbank@sheffied.ac.uk), May 2015
6 //
7 // A module to validate clustering algorithms.
8 // Compares the output of different clustering algorithms run over a pi0 sample.
9 //
10 // Usage: Specify the hit finder (HitsModuleLabel) and the clustering outputs
11 // to validate (ClusterModuleLabels) in the fhicl file.
12 // Module will make validation plots for all clusterings specified and also
13 // produce comparison plots. Number of clusterings to analyse can be one or more.
14 // Saves everything in the file validationHistograms.root.
16 
24 #include "fhiclcpp/ParameterSet.h"
25 
26 // LArSoft includes
34 
35 // ROOT & STL includes
36 #include "TCanvas.h"
37 #include "TEfficiency.h"
38 #include "TFile.h"
39 #include "TH1.h"
40 #include "TH2D.h"
41 #include "TLegend.h"
42 #include "TProfile.h"
43 #include "TStyle.h"
44 
45 #include <algorithm>
46 #include <fstream>
47 #include <iostream>
48 #include <map>
49 #include <memory>
50 
53  class ClusterAnalyser;
54  class ClusterCounter;
55 }
56 
57 enum class ClusterID : int {};
58 enum class TrackID : int {};
59 typedef std::vector<ClusterID> ClusterIDs;
60 typedef std::vector<TrackID> TrackIDs;
61 
63 public:
64  void AddHitPreClustering(TrackID id);
65  void AddSignalHitPostClustering(ClusterID id);
66  void AddNoiseHitPostClustering(ClusterID id);
67  void AssociateClusterAndTrack(ClusterID clusID, TrackID trackID);
68  double GetCompleteness(ClusterID id);
69  double GetCleanliness(ClusterID id);
70  double GetEfficiency(TrackID id);
71  ClusterIDs GetListOfClusterIDs();
72  TrackIDs GetListOfTrackIDs();
73  int GetNumberHitsFromTrack(TrackID id);
74  int GetNumberHitsInCluster(ClusterID id);
75  int GetNumberHitsInPlane();
76  std::vector<std::pair<TrackID, ClusterIDs>> GetPhotons();
77  TrackID GetTrack(ClusterID id);
78  bool IsNoise(ClusterID id);
79  bool IsNoise(TrackID id);
80  bool PassesCut();
81 
82 private:
83  std::map<TrackID, int> numHitsPreClustering;
84  std::map<ClusterID, int> numSignalHitsPostClustering;
85  std::map<ClusterID, int> numNoiseHitsPostClustering;
86  std::map<ClusterID, TrackID> clusterToTrackID;
87  std::map<TrackID, ClusterIDs> trackToClusterIDs;
88  std::map<TrackID, std::map<std::string, double>> particleProperties;
89  std::map<TrackID, simb::MCParticle> trueParticles;
90 
94 };
95 
97 {
98  ++numHitsPreClustering[trackID];
99 }
100 
102 {
103  ++numSignalHitsPostClustering[clusID];
104 }
105 
107 {
108  ++numNoiseHitsPostClustering[clusID];
109 }
110 
112  TrackID trackID)
113 {
114  clusterToTrackID[clusID] = trackID;
115  trackToClusterIDs[trackID].push_back(clusID);
116 }
117 
119 {
120  return (double)numSignalHitsPostClustering[clusID] /
121  (double)numHitsPreClustering[clusterToTrackID[clusID]];
122 }
123 
125 {
126  return (double)numSignalHitsPostClustering[clusID] / (double)(GetNumberHitsInCluster(clusID));
127 }
128 
130 {
131  return 1 / (double)trackToClusterIDs.at(trackID).size();
132 }
133 
135 {
136  return numHitsPreClustering[trackID];
137 }
138 
140 {
141  return numSignalHitsPostClustering[clusID] + numNoiseHitsPostClustering[clusID];
142 }
143 
145 {
146  int nHits = 0;
147  for (auto& trackHits : numHitsPreClustering)
148  nHits += trackHits.second;
149  return nHits;
150 }
151 
153 {
154  ClusterIDs v;
155  for (std::map<ClusterID, TrackID>::iterator i = clusterToTrackID.begin();
156  i != clusterToTrackID.end();
157  i++)
158  v.push_back(i->first);
159  return v;
160 }
161 
163 {
164  TrackIDs v;
165  for (std::map<TrackID, ClusterIDs>::iterator i = trackToClusterIDs.begin();
166  i != trackToClusterIDs.end();
167  i++)
168  v.push_back(i->first);
169  return v;
170 }
171 
172 std::vector<std::pair<TrackID, ClusterIDs>> ClusteringValidation::ClusterCounter::GetPhotons()
173 {
174  std::vector<std::pair<TrackID, ClusterIDs>> photonVector;
175  for (unsigned int track = 0; track < GetListOfTrackIDs().size(); ++track)
176  if (!IsNoise(GetListOfTrackIDs().at(track)) &&
177  pi_serv->TrackIdToParticle_P((int)GetListOfTrackIDs().at(track))->PdgCode() == 22)
178  photonVector.push_back(std::pair<TrackID, ClusterIDs>(
179  GetListOfTrackIDs().at(track), trackToClusterIDs.at(GetListOfTrackIDs().at(track))));
180  return photonVector;
181 }
182 
184 {
185  return clusterToTrackID.at(id);
186 }
187 
189 {
190  return IsNoise(clusterToTrackID.at(clusID));
191 }
192 
194 {
195  return (int)trackID == 0 ? true : false;
196 }
197 
199 {
200  if (GetPhotons().size() > 2 || GetPhotons().size() == 0) return false;
201  TrackIDs goodPhotons;
202  for (unsigned int photon = 0; photon < GetPhotons().size(); ++photon)
203  for (unsigned int cluster = 0; cluster < GetPhotons().at(photon).second.size(); ++cluster)
204  if (GetCompleteness(GetPhotons().at(photon).second.at(cluster)) > 0.5)
205  goodPhotons.push_back(GetPhotons().at(photon).first);
206  if ((GetPhotons().size() == 2 && goodPhotons.size() > 2) ||
207  (GetPhotons().size() == 1 && goodPhotons.size() > 1))
208  std::cout << "More than 2 with >50%?!" << std::endl;
209  bool pass = ((GetPhotons().size() == 2 && goodPhotons.size() == 2) ||
210  (GetPhotons().size() == 1 && goodPhotons.size() == 1));
211  return pass;
212 }
213 
214 // --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
216 public:
217  explicit ClusterAnalyser(std::string& label);
218 
219  void Analyse(detinfo::DetectorClocksData const& clockData,
222  const art::FindManyP<recob::Hit>& fmh,
223  int numHits);
224  TrackID FindTrackID(detinfo::DetectorClocksData const& clockData, art::Ptr<recob::Hit>& hit);
225  TrackID FindTrueTrack(detinfo::DetectorClocksData const& clockData,
226  std::vector<art::Ptr<recob::Hit>>& clusterHits);
227  double FindPhotonAngle();
228  double GetEndTrackDistance(TrackID id1, TrackID id2);
229  const simb::MCParticle* GetPi0();
230  TObjArray GetHistograms();
231  void MakeHistograms();
232  void WriteFile();
233 
234 private:
235  // Clustering properties
236  std::string fClusterLabel;
237 
238  // hists
242  TH1* hPi0Angle;
257  TProfile* hCleanlinessAngle;
261  TProfile* hComplCleanlAngle;
264  TEfficiency* hEfficiencyAngle;
265  TEfficiency* hEfficiencyEnergy;
268  TObjArray fHistArray;
269 
270  std::map<unsigned int, std::map<unsigned int, std::unique_ptr<ClusterCounter>>> clusterMap;
271  std::map<TrackID, const simb::MCParticle*> trueParticles;
272 
273  // Services
277 };
278 
280 {
281 
282  fClusterLabel = clusterLabel;
283 
284  // Make the histograms
285  hCompleteness = new TH1D("Completeness", ";Completeness;", 101, 0, 1.01);
286  hCompletenessEnergy =
287  new TProfile("CompletenessEnergy", ";True Energy (GeV);Completeness", 100, 0, 2);
288  hCompletenessAngle =
289  new TProfile("CompletenessAngle", ";True Angle (deg);Completeness;", 100, 0, 180);
290  hCompletenessConversionDistance = new TProfile(
291  "CompletenessConversionDistance", ";True Distance from Vertex (cm);Completeness", 100, 0, 200);
292  hCompletenessConversionSeparation = new TProfile("CompletenessConversionSeparation",
293  ";True Conversion Separation (cm);Completeness",
294  100,
295  0,
296  200);
297  hCleanliness = new TH1D("Cleanliness", ";Cleanliness;", 101, 0, 1.01);
298  hCleanlinessEnergy =
299  new TProfile("CleanlinessEnergy", ";True Energy (GeV);Cleanliness", 100, 0, 2);
300  hCleanlinessAngle =
301  new TProfile("CleanlinessAngle", ";True Angle (deg);Cleanliness;", 100, 0, 180);
302  hCleanlinessConversionDistance = new TProfile(
303  "CleanlinessConversionDistance", ";True Distance from Vertex (cm);Cleanliness", 100, 0, 200);
304  hCleanlinessConversionSeparation = new TProfile(
305  "CleanlinessConversionSeparation", ";True Conversion Separation (cm);Cleanliness", 100, 0, 200);
306  hComplCleanl = new TH1D("CompletenessCleanliness", ";Completeness * Cleanliness;", 101, 0, 1.01);
307  hComplCleanlEnergy = new TProfile(
308  "CompletenessCleanlinessEnergy", ";True Energy (GeV);Completeness * Cleanliness", 100, 0, 2);
309  hComplCleanlAngle = new TProfile(
310  "CompletenessCleanlinessAngle", ";True Angle (deg);Completeness * Cleanliness;", 100, 0, 180);
311  hComplCleanlConversionDistance =
312  new TProfile("CompletenessCleanlinessConversionDistance",
313  ";True Distance from Vertex (cm);Completeness * Cleanliness",
314  100,
315  0,
316  200);
317  hComplCleanlConversionSeparation =
318  new TProfile("CompletenessCleanlinessConversionSeparation",
319  ";True Conversion Separation (cm);Completeness * Cleanliness",
320  100,
321  0,
322  200);
323  hPi0Energy = new TH1D("Pi0EnergyCut", ";True Energy (GeV);", 25, 0, 2);
324  hPi0Energy->Sumw2();
325  hPi0Angle = new TH1D("Pi0AngleCut", ";True Angle (deg);", 25, 0, 180);
326  hPi0Angle->Sumw2();
327  hPi0ConversionDistance =
328  new TH1D("Pi0ConversionDistanceCut", ";True Distance from Vertex (cm);", 25, 0, 200);
329  hPi0ConversionDistance->Sumw2();
330  hPi0ConversionSeparation =
331  new TH1D("Pi0ConversionSeparationCut", ";True Separation from Vertex (cm);", 25, 0, 200);
332  hPi0ConversionSeparation->Sumw2();
333  hPi0EnergyCut = new TH1D("Pi0EnergyCut", ";True Energy (GeV);Efficiency", 25, 0, 2);
334  hPi0EnergyCut->Sumw2();
335  hPi0AngleCut = new TH1D("Pi0AngleCut", ";True Angle (deg);Efficiency", 25, 0, 180);
336  hPi0AngleCut->Sumw2();
337  hPi0ConversionDistanceCut =
338  new TH1D("Pi0ConversionDistanceCut", ";True Distance from Vertex (cm);Efficiency", 25, 0, 200);
339  hPi0ConversionDistanceCut->Sumw2();
340  hPi0ConversionSeparationCut = new TH1D(
341  "Pi0ConversionSeparationCut", ";True Separation from Vertex (cm);Efficiency", 25, 0, 200);
342  hPi0ConversionSeparationCut->Sumw2();
343  hNumHitsCompleteness =
344  new TH2D("NumHitsCompleteness", ";Completeness;Size of Cluster", 101, 0, 1.01, 100, 0, 100);
345  hNumHitsEnergy =
346  new TH2D("NumHitsEnergy", ";True Energy (GeV);Size of Cluster", 100, 0, 2, 100, 0, 100);
347 }
348 
352  const art::FindManyP<recob::Hit>& fmh,
353  int minHits)
354 {
355 
356  // Make a map of cluster counters in TPC/plane space
357  for (auto const& id : geometry->Iterate<geo::PlaneID>(geo::CryostatID{0})) {
358  auto const [tpc, plane] = std::make_tuple(id.TPC, id.Plane);
359  clusterMap[tpc][plane] = std::make_unique<ClusterCounter>();
360  }
361 
362  // Save preclustered hits
363  for (size_t hitIt = 0; hitIt < hits.size(); ++hitIt) {
364  art::Ptr<recob::Hit> hit = hits.at(hitIt);
365  TrackID trackID = FindTrackID(clockData, hit);
366  clusterMap[hit->WireID().TPC % 2][hit->WireID().Plane]->AddHitPreClustering(trackID);
367  }
368 
369  // Save true tracks
370  trueParticles.clear();
371  const sim::ParticleList& particles = pi_serv->ParticleList();
372  for (sim::ParticleList::const_iterator particleIt = particles.begin();
373  particleIt != particles.end();
374  ++particleIt) {
375  const simb::MCParticle* particle = particleIt->second;
376  trueParticles[(TrackID)particle->TrackId()] = particle;
377  }
378 
379  // Save the clustered hits
380  for (size_t clusIt = 0; clusIt < clusters.size(); ++clusIt) {
381 
382  // Get cluster information
383  unsigned int tpc = clusters.at(clusIt)->Plane().TPC % 2;
384  unsigned int plane = clusters.at(clusIt)->Plane().Plane;
385  ClusterID id = (ClusterID)clusters.at(clusIt)->ID();
386 
387  // Only analyse planes with enough hits in to fairly represent the clustering
388  if (clusterMap[tpc][plane]->GetNumberHitsInPlane() < minHits) continue;
389 
390  // Get the hits from the cluster
391  std::vector<art::Ptr<recob::Hit>> clusterHits = fmh.at(clusIt);
392 
393  if (clusterHits.size() < 10) continue;
394 
395  // Find which track this cluster belongs to
396  TrackID trueTrackID = FindTrueTrack(clockData, clusterHits);
397 
398  // Save the info for this cluster
399  clusterMap[tpc][plane]->AssociateClusterAndTrack(id, trueTrackID);
400  for (std::vector<art::Ptr<recob::Hit>>::iterator clusHitIt = clusterHits.begin();
401  clusHitIt != clusterHits.end();
402  ++clusHitIt) {
403  art::Ptr<recob::Hit> hit = *clusHitIt;
404  TrackID trackID = FindTrackID(clockData, hit);
405  if (trackID == trueTrackID)
406  clusterMap[tpc][plane]->AddSignalHitPostClustering(id);
407  else
408  clusterMap[tpc][plane]->AddNoiseHitPostClustering(id);
409  }
410  }
411 
412  this->MakeHistograms();
413 }
414 
416  detinfo::DetectorClocksData const& clockData,
418 {
419  double particleEnergy = 0;
420  TrackID likelyTrackID = (TrackID)0;
421  std::vector<sim::TrackIDE> trackIDs = bt_serv->HitToTrackIDEs(clockData, hit);
422  for (unsigned int idIt = 0; idIt < trackIDs.size(); ++idIt) {
423  if (trackIDs.at(idIt).energy > particleEnergy) {
424  particleEnergy = trackIDs.at(idIt).energy;
425  likelyTrackID = (TrackID)TMath::Abs(trackIDs.at(idIt).trackID);
426  }
427  }
428  return likelyTrackID;
429 }
430 
432  detinfo::DetectorClocksData const& clockData,
433  std::vector<art::Ptr<recob::Hit>>& clusterHits)
434 {
435  std::map<TrackID, double> trackMap;
436  for (std::vector<art::Ptr<recob::Hit>>::iterator clusHitIt = clusterHits.begin();
437  clusHitIt != clusterHits.end();
438  ++clusHitIt) {
439  art::Ptr<recob::Hit> hit = *clusHitIt;
440  TrackID trackID = FindTrackID(clockData, hit);
441  trackMap[trackID] += hit->Integral();
442  }
443  //return std::max_element(trackMap.begin(), trackMap.end(), [](const std::pair<int,double>& p1, const std::pair<int,double>& p2) {return p1.second < p2.second;} )->first;
444  double highestCharge = 0;
445  TrackID clusterTrack = (TrackID)0;
446  for (std::map<TrackID, double>::iterator trackIt = trackMap.begin(); trackIt != trackMap.end();
447  ++trackIt)
448  if (trackIt->second > highestCharge) {
449  highestCharge = trackIt->second;
450  clusterTrack = trackIt->first;
451  }
452  return clusterTrack;
453 }
454 
456 {
457  const simb::MCParticle* pi0 = GetPi0();
458  if (pi0->NumberDaughters() != 2) return -999;
459  double angle = (trueParticles.at((TrackID)pi0->Daughter(0))
460  ->Momentum()
461  .Angle(trueParticles.at((TrackID)pi0->Daughter(1))->Momentum().Vect()) *
462  (180 / TMath::Pi()));
463  return angle;
464 }
465 
467 {
468  const simb::MCParticle* pi0 = nullptr;
469  for (std::map<TrackID, const simb::MCParticle*>::iterator particleIt = trueParticles.begin();
470  particleIt != trueParticles.end();
471  ++particleIt)
472  if (particleIt->second->PdgCode() == 111) pi0 = particleIt->second;
473  return pi0;
474 }
475 
477 {
478  return TMath::Sqrt(
479  TMath::Power(
480  trueParticles.at(id1)->EndPosition().X() - trueParticles.at(id2)->EndPosition().X(), 2) +
481  TMath::Power(
482  trueParticles.at(id1)->EndPosition().Y() - trueParticles.at(id2)->EndPosition().Y(), 2) +
483  TMath::Power(
484  trueParticles.at(id1)->EndPosition().Z() - trueParticles.at(id2)->EndPosition().Z(), 2));
485 }
486 
488 {
489 
490  // Make efficiency histograms
491  hEfficiencyEnergy = new TEfficiency(*hPi0EnergyCut, *hPi0Energy);
492  hEfficiencyAngle = new TEfficiency(*hPi0AngleCut, *hPi0Angle);
493  hEfficiencyConversionDistance =
494  new TEfficiency(*hPi0ConversionDistanceCut, *hPi0ConversionDistance);
495  hEfficiencyConversionSeparation =
496  new TEfficiency(*hPi0ConversionSeparationCut, *hPi0ConversionSeparation);
497 
498  hEfficiencyEnergy->SetName("EfficiencyEnergy");
499  hEfficiencyAngle->SetName("EnergyAngle");
500  hEfficiencyConversionDistance->SetName("EfficiencyConversionDistance");
501  hEfficiencyConversionSeparation->SetName("EfficiencyConversionSeparation");
502 
503  // Add all the histograms to the object array
504  fHistArray.Add(hCompleteness);
505  fHistArray.Add(hCompletenessEnergy);
506  fHistArray.Add(hCompletenessAngle);
507  fHistArray.Add(hCompletenessConversionDistance);
508  fHistArray.Add(hCompletenessConversionSeparation);
509  fHistArray.Add(hCleanliness);
510  fHistArray.Add(hCleanlinessEnergy);
511  fHistArray.Add(hCleanlinessAngle);
512  fHistArray.Add(hCleanlinessConversionDistance);
513  fHistArray.Add(hCleanlinessConversionSeparation);
514  fHistArray.Add(hComplCleanl);
515  fHistArray.Add(hComplCleanlEnergy);
516  fHistArray.Add(hComplCleanlAngle);
517  fHistArray.Add(hComplCleanlConversionDistance);
518  fHistArray.Add(hComplCleanlConversionSeparation);
519  fHistArray.Add(hEfficiencyEnergy);
520  fHistArray.Add(hEfficiencyAngle);
521  fHistArray.Add(hEfficiencyConversionDistance);
522  fHistArray.Add(hEfficiencyConversionSeparation);
523  fHistArray.Add(hNumHitsCompleteness);
524  fHistArray.Add(hNumHitsEnergy);
525 
526  return fHistArray;
527 }
528 
530 {
531 
532  // Loop over the tpcs and planes in the geometry
533  for (auto const& id : geometry->Iterate<geo::PlaneID>(geo::CryostatID{0})) {
534  auto const [tpc, plane] = std::make_tuple(id.TPC, id.Plane);
535  auto& counter = clusterMap[tpc][plane];
536  ClusterIDs clusterIDs = counter->GetListOfClusterIDs();
537 
538  // Fill histograms for the efficiency
539  if (counter->GetPhotons().size() == 2) {
540  hPi0Angle->Fill(FindPhotonAngle());
541  hPi0Energy->Fill(GetPi0()->Momentum().E());
542  hPi0ConversionDistance->Fill(std::min(
543  GetEndTrackDistance(counter->GetPhotons().at(0).first, (TrackID)GetPi0()->TrackId()),
544  GetEndTrackDistance(counter->GetPhotons().at(1).first, (TrackID)GetPi0()->TrackId())));
545  hPi0ConversionSeparation->Fill(
546  GetEndTrackDistance(counter->GetPhotons().at(0).first, counter->GetPhotons().at(1).first));
547  if (counter->PassesCut()) {
548  hPi0AngleCut->Fill(FindPhotonAngle());
549  hPi0EnergyCut->Fill(GetPi0()->Momentum().E());
550  hPi0ConversionDistanceCut->Fill(std::min(
551  GetEndTrackDistance(counter->GetPhotons().at(0).first, (TrackID)GetPi0()->TrackId()),
552  GetEndTrackDistance(counter->GetPhotons().at(1).first, (TrackID)GetPi0()->TrackId())));
553  hPi0ConversionSeparationCut->Fill(GetEndTrackDistance(counter->GetPhotons().at(0).first,
554  counter->GetPhotons().at(1).first));
555  }
556  else
557  std::cout << "TPC " << tpc << ", Plane " << plane << " fails the cut" << std::endl;
558  }
559 
560  // Look at all the clusters
561  for (unsigned int cluster = 0; cluster < clusterIDs.size(); ++cluster) {
562 
563  ClusterID clusID = clusterIDs.at(cluster);
564  double completeness = counter->GetCompleteness(clusID);
565  double cleanliness = counter->GetCleanliness(clusID);
566  int numClusterHits = counter->GetNumberHitsInCluster(clusID);
567 
568  // Fill histograms for this cluster
569  hCompleteness->Fill(completeness, numClusterHits);
570  hCleanliness->Fill(cleanliness, numClusterHits);
571  hComplCleanl->Fill(completeness * cleanliness, numClusterHits);
572  hNumHitsCompleteness->Fill(completeness, numClusterHits);
573 
574  // Is this cluster doesn't correspond to a true particle continue
575  if (counter->IsNoise(clusID)) continue;
576 
577  double pi0Energy = GetPi0()->Momentum().E();
578  double pi0DecayAngle = FindPhotonAngle();
579  double conversionDistance =
580  GetEndTrackDistance(counter->GetTrack(clusID), (TrackID)GetPi0()->TrackId());
581 
582  hCompletenessEnergy->Fill(pi0Energy, completeness, numClusterHits);
583  hCompletenessAngle->Fill(pi0DecayAngle, completeness, numClusterHits);
584  hCompletenessConversionDistance->Fill(conversionDistance, completeness, numClusterHits);
585  hCleanlinessEnergy->Fill(pi0Energy, cleanliness, numClusterHits);
586  hCleanlinessAngle->Fill(pi0DecayAngle, cleanliness, numClusterHits);
587  hCleanlinessConversionDistance->Fill(conversionDistance, cleanliness, numClusterHits);
588  hComplCleanlEnergy->Fill(pi0Energy, cleanliness * completeness, numClusterHits);
589  hComplCleanlAngle->Fill(pi0DecayAngle, cleanliness * completeness, numClusterHits);
590  hComplCleanlConversionDistance->Fill(
591  conversionDistance, cleanliness * completeness, numClusterHits);
592  hNumHitsEnergy->Fill(pi0Energy, numClusterHits);
593 
594  // Continue if there are not two photons in the view
595  if (counter->GetPhotons().size() != 2) continue;
596 
597  double conversionSeparation =
598  GetEndTrackDistance(counter->GetPhotons()[0].first, counter->GetPhotons()[1].first);
599 
600  hCompletenessConversionSeparation->Fill(conversionSeparation, completeness, numClusterHits);
601  hCleanlinessConversionSeparation->Fill(conversionSeparation, cleanliness, numClusterHits);
602  }
603  }
604 }
605 
607 {
608 
609  // Average completeness/cleanliness
610  double avCompleteness = hCompleteness->GetMean();
611  double avCleanliness = hCleanliness->GetMean();
612 
613  // Write file
614  std::ofstream outFile("effpur");
615  outFile << avCompleteness << " " << avCleanliness;
616  outFile.close();
617 }
618 
619 // --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
621 public:
622  explicit ClusteringValidation(fhicl::ParameterSet const& p);
623 
624 private:
625  void analyze(art::Event const& evt);
626  void beginJob();
627  void endJob();
628 
629  // Clusterings to compare and the hits which the clustering was run over
630  std::vector<std::string> fClusterModuleLabels;
631  std::string fHitsModuleLabel;
632 
633  // Minimum hits needed to analyse a plane
635 
636  // Canvas on which to save histograms
637  TCanvas* fCanvas;
638 
639  // The cluster analysers
640  // -- one for each clustering type being compared
641  std::map<std::string, std::unique_ptr<ClusterAnalyser>> clusterAnalysis;
642 };
643 
645  : EDAnalyzer(pset)
646 {
647  fMinHitsInPlane = pset.get<int>("MinHitsInPlane");
648  fClusterModuleLabels = pset.get<std::vector<std::string>>("ClusterModuleLabels");
649  fHitsModuleLabel = pset.get<std::string>("HitsModuleLabel");
650 
651  fCanvas = new TCanvas("fCanvas", "", 800, 600);
652  gStyle->SetOptStat(0);
653 }
654 
656 {
657  // Get the hits from the event
659  std::vector<art::Ptr<recob::Hit>> hits;
660  if (evt.getByLabel(fHitsModuleLabel, hitHandle)) art::fill_ptr_vector(hits, hitHandle);
661 
662  // Get clustering information from event
663  // and give to the ClusterAnalyser to analyse
664  for (auto clustering : fClusterModuleLabels) {
665 
666  // Get the clusters from the event
668  std::vector<art::Ptr<recob::Cluster>> clusters;
669  if (evt.getByLabel(clustering, clusterHandle)) art::fill_ptr_vector(clusters, clusterHandle);
670 
671  // Find the associations (the hits) for the clusters
672  art::FindManyP<recob::Hit> fmh(clusterHandle, evt, clustering);
673 
674  // Analyse this particular clustering
675  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
676  clusterAnalysis.at(clustering)->Analyse(clockData, hits, clusters, fmh, fMinHitsInPlane);
677  }
678 }
679 
681 {
682 
683  // Construct the cluster analysers here
684  // -- one for each of the clustering types to compare
685  for (auto clustering : fClusterModuleLabels)
686  clusterAnalysis[clustering] = (std::unique_ptr<ClusterAnalyser>)new ClusterAnalyser(clustering);
687 }
688 
690 {
691 
692  // Make a map of all the histograms for each clustering
693  std::map<std::string, TObjArray> allHistograms;
694  for (auto clustering : fClusterModuleLabels)
695  allHistograms[clustering] = clusterAnalysis.at(clustering)->GetHistograms();
696 
697  // Write histograms to file
698  TFile* file = TFile::Open("validationHistograms.root", "RECREATE");
699  for (auto clustering : fClusterModuleLabels) {
700  file->cd();
701  file->mkdir(clustering.c_str());
702  file->cd(clustering.c_str());
703  for (int histIt = 0; histIt < allHistograms.begin()->second.GetEntriesFast(); ++histIt)
704  allHistograms.at(clustering).At(histIt)->Write();
705  }
706 
707  // Write images of overlaid histograms
708  for (int histIt = 0; histIt < allHistograms.begin()->second.GetEntriesFast(); ++histIt) {
709  fCanvas->cd();
710  fCanvas->Clear();
711  const char* name = allHistograms.begin()->second.At(histIt)->GetName();
712  TLegend* l = new TLegend(0.6, 0.8, 0.8, 0.9, name, "brNDC");
713  int clusterings = 1;
714  for (std::map<std::string, TObjArray>::iterator clusteringIt = allHistograms.begin();
715  clusteringIt != allHistograms.end();
716  ++clusteringIt, ++clusterings) {
717  TH1* h = (TH1*)allHistograms.at(clusteringIt->first).At(histIt);
718  h->SetLineColor(clusterings);
719  h->SetMarkerColor(clusterings);
720  if (clusterings == 1)
721  h->Draw();
722  else
723  h->Draw("same");
724  l->AddEntry(h, clusteringIt->first.c_str(), "lp");
725  }
726  l->Draw("same");
727  //fCanvas->SaveAs(name+TString(".png"));
728  file->cd();
729  fCanvas->Write(name);
730  }
731 
732  file->Close();
733  delete file;
734 
735  if (clusterAnalysis.find("blurredclusteringdc") != clusterAnalysis.end())
736  clusterAnalysis.at("blurredclusteringdc")->WriteFile();
737 }
738 
art::ServiceHandle< cheat::BackTrackerService const > bt_serv
intermediate_table::iterator iterator
std::map< ClusterID, int > numNoiseHitsPostClustering
art::ServiceHandle< cheat::BackTrackerService const > bt_serv
std::vector< TrackID > TrackIDs
Declaration of signal hit object.
art::ServiceHandle< geo::Geometry const > geometry
list_type::const_iterator const_iterator
Definition: ParticleList.h:132
double GetEndTrackDistance(TrackID id1, TrackID id2)
The data type to uniquely identify a Plane.
Definition: geo_types.h:463
float Integral() const
Integral under the calibrated signal waveform of the hit, in tick x ADC units.
Definition: Hit.h:244
geo::WireID const & WireID() const
Initial tdc tick for hit.
Definition: Hit.h:280
std::map< TrackID, ClusterIDs > trackToClusterIDs
std::map< TrackID, std::map< std::string, double > > particleProperties
Cluster finding and building.
Particle class.
int NumberDaughters() const
Definition: MCParticle.h:218
int TrackId() const
Definition: MCParticle.h:211
int Daughter(const int i) const
Definition: MCParticle.cxx:118
std::map< ClusterID, int > numSignalHitsPostClustering
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:289
TrackID FindTrackID(detinfo::DetectorClocksData const &clockData, art::Ptr< recob::Hit > &hit)
void hits()
Definition: readHits.C:15
auto counter(T begin, T end)
Returns an object to iterate values from begin to end in a range-for loop.
Definition: counter.h:295
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
void beginJob()
Definition: Breakpoints.cc:14
Float_t E
Definition: plot.C:20
std::vector< std::pair< TrackID, ClusterIDs > > GetPhotons()
T get(std::string const &key) const
Definition: ParameterSet.h:314
art::ServiceHandle< cheat::ParticleInventoryService const > pi_serv
iterator begin()
Definition: ParticleList.h:305
TrackID FindTrueTrack(detinfo::DetectorClocksData const &clockData, std::vector< art::Ptr< recob::Hit >> &clusterHits)
std::vector< ClusterID > ClusterIDs
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:481
Declaration of cluster object.
std::map< unsigned int, std::map< unsigned int, std::unique_ptr< ClusterCounter > > > clusterMap
Detector simulation of raw signals on wires.
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
std::map< TrackID, const simb::MCParticle * > trueParticles
Contains all timing reference information for the detector.
std::map< TrackID, simb::MCParticle > trueParticles
TFile * file
void AssociateClusterAndTrack(ClusterID clusID, TrackID trackID)
TCEvent evt
Definition: DataStructs.cxx:8
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:306
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:399
void Analyse(detinfo::DetectorClocksData const &clockData, std::vector< art::Ptr< recob::Hit >> &hits, std::vector< art::Ptr< recob::Cluster >> &clusters, const art::FindManyP< recob::Hit > &fmh, int numHits)
second_as<> second
Type of time stored in seconds, in double precision.
Definition: spacetime.h:82
recob::tracking::Plane Plane
Definition: TrackState.h:17
Float_t track
Definition: plot.C:35
std::map< ClusterID, TrackID > clusterToTrackID
art framework interface to geometry description
The data type to uniquely identify a cryostat.
Definition: geo_types.h:192
std::map< std::string, std::unique_ptr< ClusterAnalyser > > clusterAnalysis
art::ServiceHandle< cheat::ParticleInventoryService > pi_serv
art::ServiceHandle< geo::Geometry const > geometry