LArSoft  v10_04_05
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 
93 };
94 
96 {
97  ++numHitsPreClustering[trackID];
98 }
99 
101 {
102  ++numSignalHitsPostClustering[clusID];
103 }
104 
106 {
107  ++numNoiseHitsPostClustering[clusID];
108 }
109 
111  TrackID trackID)
112 {
113  clusterToTrackID[clusID] = trackID;
114  trackToClusterIDs[trackID].push_back(clusID);
115 }
116 
118 {
119  return (double)numSignalHitsPostClustering[clusID] /
120  (double)numHitsPreClustering[clusterToTrackID[clusID]];
121 }
122 
124 {
125  return (double)numSignalHitsPostClustering[clusID] / (double)(GetNumberHitsInCluster(clusID));
126 }
127 
129 {
130  return 1 / (double)trackToClusterIDs.at(trackID).size();
131 }
132 
134 {
135  return numHitsPreClustering[trackID];
136 }
137 
139 {
140  return numSignalHitsPostClustering[clusID] + numNoiseHitsPostClustering[clusID];
141 }
142 
144 {
145  int nHits = 0;
146  for (auto& trackHits : numHitsPreClustering)
147  nHits += trackHits.second;
148  return nHits;
149 }
150 
152 {
153  ClusterIDs v;
154  for (std::map<ClusterID, TrackID>::iterator i = clusterToTrackID.begin();
155  i != clusterToTrackID.end();
156  i++)
157  v.push_back(i->first);
158  return v;
159 }
160 
162 {
163  TrackIDs v;
164  for (std::map<TrackID, ClusterIDs>::iterator i = trackToClusterIDs.begin();
165  i != trackToClusterIDs.end();
166  i++)
167  v.push_back(i->first);
168  return v;
169 }
170 
171 std::vector<std::pair<TrackID, ClusterIDs>> ClusteringValidation::ClusterCounter::GetPhotons()
172 {
173  std::vector<std::pair<TrackID, ClusterIDs>> photonVector;
174  for (unsigned int track = 0; track < GetListOfTrackIDs().size(); ++track)
175  if (!IsNoise(GetListOfTrackIDs().at(track)) &&
176  pi_serv->TrackIdToParticle_P((int)GetListOfTrackIDs().at(track))->PdgCode() == 22)
177  photonVector.push_back(std::pair<TrackID, ClusterIDs>(
178  GetListOfTrackIDs().at(track), trackToClusterIDs.at(GetListOfTrackIDs().at(track))));
179  return photonVector;
180 }
181 
183 {
184  return clusterToTrackID.at(id);
185 }
186 
188 {
189  return IsNoise(clusterToTrackID.at(clusID));
190 }
191 
193 {
194  return (int)trackID == 0 ? true : false;
195 }
196 
198 {
199  if (GetPhotons().size() > 2 || GetPhotons().size() == 0) return false;
200  TrackIDs goodPhotons;
201  for (unsigned int photon = 0; photon < GetPhotons().size(); ++photon)
202  for (unsigned int cluster = 0; cluster < GetPhotons().at(photon).second.size(); ++cluster)
203  if (GetCompleteness(GetPhotons().at(photon).second.at(cluster)) > 0.5)
204  goodPhotons.push_back(GetPhotons().at(photon).first);
205  if ((GetPhotons().size() == 2 && goodPhotons.size() > 2) ||
206  (GetPhotons().size() == 1 && goodPhotons.size() > 1))
207  std::cout << "More than 2 with >50%?!" << std::endl;
208  bool pass = ((GetPhotons().size() == 2 && goodPhotons.size() == 2) ||
209  (GetPhotons().size() == 1 && goodPhotons.size() == 1));
210  return pass;
211 }
212 
213 // --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
215 public:
216  explicit ClusterAnalyser(std::string& label);
217 
218  void Analyse(detinfo::DetectorClocksData const& clockData,
221  const art::FindManyP<recob::Hit>& fmh,
222  int numHits);
223  TrackID FindTrackID(detinfo::DetectorClocksData const& clockData, art::Ptr<recob::Hit>& hit);
224  TrackID FindTrueTrack(detinfo::DetectorClocksData const& clockData,
225  std::vector<art::Ptr<recob::Hit>>& clusterHits);
226  double FindPhotonAngle();
227  double GetEndTrackDistance(TrackID id1, TrackID id2);
228  const simb::MCParticle* GetPi0();
229  TObjArray GetHistograms();
230  void MakeHistograms();
231  void WriteFile();
232 
233 private:
234  // Clustering properties
235  std::string fClusterLabel;
236 
237  // hists
241  TH1* hPi0Angle;
256  TProfile* hCleanlinessAngle;
260  TProfile* hComplCleanlAngle;
263  TEfficiency* hEfficiencyAngle;
264  TEfficiency* hEfficiencyEnergy;
267  TObjArray fHistArray;
268 
269  std::map<unsigned int, std::map<unsigned int, std::unique_ptr<ClusterCounter>>> clusterMap;
270  std::map<TrackID, const simb::MCParticle*> trueParticles;
271 
272  // Services
276 };
277 
279 {
280 
281  fClusterLabel = clusterLabel;
282 
283  // Make the histograms
284  hCompleteness = new TH1D("Completeness", ";Completeness;", 101, 0, 1.01);
285  hCompletenessEnergy =
286  new TProfile("CompletenessEnergy", ";True Energy (GeV);Completeness", 100, 0, 2);
287  hCompletenessAngle =
288  new TProfile("CompletenessAngle", ";True Angle (deg);Completeness;", 100, 0, 180);
289  hCompletenessConversionDistance = new TProfile(
290  "CompletenessConversionDistance", ";True Distance from Vertex (cm);Completeness", 100, 0, 200);
291  hCompletenessConversionSeparation = new TProfile("CompletenessConversionSeparation",
292  ";True Conversion Separation (cm);Completeness",
293  100,
294  0,
295  200);
296  hCleanliness = new TH1D("Cleanliness", ";Cleanliness;", 101, 0, 1.01);
297  hCleanlinessEnergy =
298  new TProfile("CleanlinessEnergy", ";True Energy (GeV);Cleanliness", 100, 0, 2);
299  hCleanlinessAngle =
300  new TProfile("CleanlinessAngle", ";True Angle (deg);Cleanliness;", 100, 0, 180);
301  hCleanlinessConversionDistance = new TProfile(
302  "CleanlinessConversionDistance", ";True Distance from Vertex (cm);Cleanliness", 100, 0, 200);
303  hCleanlinessConversionSeparation = new TProfile(
304  "CleanlinessConversionSeparation", ";True Conversion Separation (cm);Cleanliness", 100, 0, 200);
305  hComplCleanl = new TH1D("CompletenessCleanliness", ";Completeness * Cleanliness;", 101, 0, 1.01);
306  hComplCleanlEnergy = new TProfile(
307  "CompletenessCleanlinessEnergy", ";True Energy (GeV);Completeness * Cleanliness", 100, 0, 2);
308  hComplCleanlAngle = new TProfile(
309  "CompletenessCleanlinessAngle", ";True Angle (deg);Completeness * Cleanliness;", 100, 0, 180);
310  hComplCleanlConversionDistance =
311  new TProfile("CompletenessCleanlinessConversionDistance",
312  ";True Distance from Vertex (cm);Completeness * Cleanliness",
313  100,
314  0,
315  200);
316  hComplCleanlConversionSeparation =
317  new TProfile("CompletenessCleanlinessConversionSeparation",
318  ";True Conversion Separation (cm);Completeness * Cleanliness",
319  100,
320  0,
321  200);
322  hPi0Energy = new TH1D("Pi0EnergyCut", ";True Energy (GeV);", 25, 0, 2);
323  hPi0Energy->Sumw2();
324  hPi0Angle = new TH1D("Pi0AngleCut", ";True Angle (deg);", 25, 0, 180);
325  hPi0Angle->Sumw2();
326  hPi0ConversionDistance =
327  new TH1D("Pi0ConversionDistanceCut", ";True Distance from Vertex (cm);", 25, 0, 200);
328  hPi0ConversionDistance->Sumw2();
329  hPi0ConversionSeparation =
330  new TH1D("Pi0ConversionSeparationCut", ";True Separation from Vertex (cm);", 25, 0, 200);
331  hPi0ConversionSeparation->Sumw2();
332  hPi0EnergyCut = new TH1D("Pi0EnergyCut", ";True Energy (GeV);Efficiency", 25, 0, 2);
333  hPi0EnergyCut->Sumw2();
334  hPi0AngleCut = new TH1D("Pi0AngleCut", ";True Angle (deg);Efficiency", 25, 0, 180);
335  hPi0AngleCut->Sumw2();
336  hPi0ConversionDistanceCut =
337  new TH1D("Pi0ConversionDistanceCut", ";True Distance from Vertex (cm);Efficiency", 25, 0, 200);
338  hPi0ConversionDistanceCut->Sumw2();
339  hPi0ConversionSeparationCut = new TH1D(
340  "Pi0ConversionSeparationCut", ";True Separation from Vertex (cm);Efficiency", 25, 0, 200);
341  hPi0ConversionSeparationCut->Sumw2();
342  hNumHitsCompleteness =
343  new TH2D("NumHitsCompleteness", ";Completeness;Size of Cluster", 101, 0, 1.01, 100, 0, 100);
344  hNumHitsEnergy =
345  new TH2D("NumHitsEnergy", ";True Energy (GeV);Size of Cluster", 100, 0, 2, 100, 0, 100);
346 }
347 
351  const art::FindManyP<recob::Hit>& fmh,
352  int minHits)
353 {
354 
355  // Make a map of cluster counters in TPC/plane space
356  for (auto const& id : wireReadoutGeom.Iterate<geo::PlaneID>(geo::CryostatID{0})) {
357  auto const [tpc, plane] = std::make_tuple(id.TPC, id.Plane);
358  clusterMap[tpc][plane] = std::make_unique<ClusterCounter>();
359  }
360 
361  // Save preclustered hits
362  for (size_t hitIt = 0; hitIt < hits.size(); ++hitIt) {
363  art::Ptr<recob::Hit> hit = hits.at(hitIt);
364  TrackID trackID = FindTrackID(clockData, hit);
365  clusterMap[hit->WireID().TPC % 2][hit->WireID().Plane]->AddHitPreClustering(trackID);
366  }
367 
368  // Save true tracks
369  trueParticles.clear();
370  const sim::ParticleList& particles = pi_serv->ParticleList();
371  for (sim::ParticleList::const_iterator particleIt = particles.begin();
372  particleIt != particles.end();
373  ++particleIt) {
374  const simb::MCParticle* particle = particleIt->second;
375  trueParticles[(TrackID)particle->TrackId()] = particle;
376  }
377 
378  // Save the clustered hits
379  for (size_t clusIt = 0; clusIt < clusters.size(); ++clusIt) {
380 
381  // Get cluster information
382  unsigned int tpc = clusters.at(clusIt)->Plane().TPC % 2;
383  unsigned int plane = clusters.at(clusIt)->Plane().Plane;
384  ClusterID id = (ClusterID)clusters.at(clusIt)->ID();
385 
386  // Only analyse planes with enough hits in to fairly represent the clustering
387  if (clusterMap[tpc][plane]->GetNumberHitsInPlane() < minHits) continue;
388 
389  // Get the hits from the cluster
390  std::vector<art::Ptr<recob::Hit>> clusterHits = fmh.at(clusIt);
391 
392  if (clusterHits.size() < 10) continue;
393 
394  // Find which track this cluster belongs to
395  TrackID trueTrackID = FindTrueTrack(clockData, clusterHits);
396 
397  // Save the info for this cluster
398  clusterMap[tpc][plane]->AssociateClusterAndTrack(id, trueTrackID);
399  for (std::vector<art::Ptr<recob::Hit>>::iterator clusHitIt = clusterHits.begin();
400  clusHitIt != clusterHits.end();
401  ++clusHitIt) {
402  art::Ptr<recob::Hit> hit = *clusHitIt;
403  TrackID trackID = FindTrackID(clockData, hit);
404  if (trackID == trueTrackID)
405  clusterMap[tpc][plane]->AddSignalHitPostClustering(id);
406  else
407  clusterMap[tpc][plane]->AddNoiseHitPostClustering(id);
408  }
409  }
410 
411  this->MakeHistograms();
412 }
413 
415  detinfo::DetectorClocksData const& clockData,
417 {
418  double particleEnergy = 0;
419  TrackID likelyTrackID = (TrackID)0;
420  std::vector<sim::TrackIDE> trackIDs = bt_serv->HitToTrackIDEs(clockData, hit);
421  for (unsigned int idIt = 0; idIt < trackIDs.size(); ++idIt) {
422  if (trackIDs.at(idIt).energy > particleEnergy) {
423  particleEnergy = trackIDs.at(idIt).energy;
424  likelyTrackID = (TrackID)TMath::Abs(trackIDs.at(idIt).trackID);
425  }
426  }
427  return likelyTrackID;
428 }
429 
431  detinfo::DetectorClocksData const& clockData,
432  std::vector<art::Ptr<recob::Hit>>& clusterHits)
433 {
434  std::map<TrackID, double> trackMap;
435  for (std::vector<art::Ptr<recob::Hit>>::iterator clusHitIt = clusterHits.begin();
436  clusHitIt != clusterHits.end();
437  ++clusHitIt) {
438  art::Ptr<recob::Hit> hit = *clusHitIt;
439  TrackID trackID = FindTrackID(clockData, hit);
440  trackMap[trackID] += hit->Integral();
441  }
442  //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;
443  double highestCharge = 0;
444  TrackID clusterTrack = (TrackID)0;
445  for (std::map<TrackID, double>::iterator trackIt = trackMap.begin(); trackIt != trackMap.end();
446  ++trackIt)
447  if (trackIt->second > highestCharge) {
448  highestCharge = trackIt->second;
449  clusterTrack = trackIt->first;
450  }
451  return clusterTrack;
452 }
453 
455 {
456  const simb::MCParticle* pi0 = GetPi0();
457  if (pi0->NumberDaughters() != 2) return -999;
458  double angle = (trueParticles.at((TrackID)pi0->Daughter(0))
459  ->Momentum()
460  .Angle(trueParticles.at((TrackID)pi0->Daughter(1))->Momentum().Vect()) *
461  (180 / TMath::Pi()));
462  return angle;
463 }
464 
466 {
467  const simb::MCParticle* pi0 = nullptr;
468  for (std::map<TrackID, const simb::MCParticle*>::iterator particleIt = trueParticles.begin();
469  particleIt != trueParticles.end();
470  ++particleIt)
471  if (particleIt->second->PdgCode() == 111) pi0 = particleIt->second;
472  return pi0;
473 }
474 
476 {
477  return TMath::Sqrt(
478  TMath::Power(
479  trueParticles.at(id1)->EndPosition().X() - trueParticles.at(id2)->EndPosition().X(), 2) +
480  TMath::Power(
481  trueParticles.at(id1)->EndPosition().Y() - trueParticles.at(id2)->EndPosition().Y(), 2) +
482  TMath::Power(
483  trueParticles.at(id1)->EndPosition().Z() - trueParticles.at(id2)->EndPosition().Z(), 2));
484 }
485 
487 {
488 
489  // Make efficiency histograms
490  hEfficiencyEnergy = new TEfficiency(*hPi0EnergyCut, *hPi0Energy);
491  hEfficiencyAngle = new TEfficiency(*hPi0AngleCut, *hPi0Angle);
492  hEfficiencyConversionDistance =
493  new TEfficiency(*hPi0ConversionDistanceCut, *hPi0ConversionDistance);
494  hEfficiencyConversionSeparation =
495  new TEfficiency(*hPi0ConversionSeparationCut, *hPi0ConversionSeparation);
496 
497  hEfficiencyEnergy->SetName("EfficiencyEnergy");
498  hEfficiencyAngle->SetName("EnergyAngle");
499  hEfficiencyConversionDistance->SetName("EfficiencyConversionDistance");
500  hEfficiencyConversionSeparation->SetName("EfficiencyConversionSeparation");
501 
502  // Add all the histograms to the object array
503  fHistArray.Add(hCompleteness);
504  fHistArray.Add(hCompletenessEnergy);
505  fHistArray.Add(hCompletenessAngle);
506  fHistArray.Add(hCompletenessConversionDistance);
507  fHistArray.Add(hCompletenessConversionSeparation);
508  fHistArray.Add(hCleanliness);
509  fHistArray.Add(hCleanlinessEnergy);
510  fHistArray.Add(hCleanlinessAngle);
511  fHistArray.Add(hCleanlinessConversionDistance);
512  fHistArray.Add(hCleanlinessConversionSeparation);
513  fHistArray.Add(hComplCleanl);
514  fHistArray.Add(hComplCleanlEnergy);
515  fHistArray.Add(hComplCleanlAngle);
516  fHistArray.Add(hComplCleanlConversionDistance);
517  fHistArray.Add(hComplCleanlConversionSeparation);
518  fHistArray.Add(hEfficiencyEnergy);
519  fHistArray.Add(hEfficiencyAngle);
520  fHistArray.Add(hEfficiencyConversionDistance);
521  fHistArray.Add(hEfficiencyConversionSeparation);
522  fHistArray.Add(hNumHitsCompleteness);
523  fHistArray.Add(hNumHitsEnergy);
524 
525  return fHistArray;
526 }
527 
529 {
530 
531  // Loop over the tpcs and planes in the geometry
532  for (auto const& id : wireReadoutGeom.Iterate<geo::PlaneID>(geo::CryostatID{0})) {
533  auto const [tpc, plane] = std::make_tuple(id.TPC, id.Plane);
534  auto& counter = clusterMap[tpc][plane];
535  ClusterIDs clusterIDs = counter->GetListOfClusterIDs();
536 
537  // Fill histograms for the efficiency
538  if (counter->GetPhotons().size() == 2) {
539  hPi0Angle->Fill(FindPhotonAngle());
540  hPi0Energy->Fill(GetPi0()->Momentum().E());
541  hPi0ConversionDistance->Fill(std::min(
542  GetEndTrackDistance(counter->GetPhotons().at(0).first, (TrackID)GetPi0()->TrackId()),
543  GetEndTrackDistance(counter->GetPhotons().at(1).first, (TrackID)GetPi0()->TrackId())));
544  hPi0ConversionSeparation->Fill(
545  GetEndTrackDistance(counter->GetPhotons().at(0).first, counter->GetPhotons().at(1).first));
546  if (counter->PassesCut()) {
547  hPi0AngleCut->Fill(FindPhotonAngle());
548  hPi0EnergyCut->Fill(GetPi0()->Momentum().E());
549  hPi0ConversionDistanceCut->Fill(std::min(
550  GetEndTrackDistance(counter->GetPhotons().at(0).first, (TrackID)GetPi0()->TrackId()),
551  GetEndTrackDistance(counter->GetPhotons().at(1).first, (TrackID)GetPi0()->TrackId())));
552  hPi0ConversionSeparationCut->Fill(GetEndTrackDistance(counter->GetPhotons().at(0).first,
553  counter->GetPhotons().at(1).first));
554  }
555  else
556  std::cout << "TPC " << tpc << ", Plane " << plane << " fails the cut" << std::endl;
557  }
558 
559  // Look at all the clusters
560  for (unsigned int cluster = 0; cluster < clusterIDs.size(); ++cluster) {
561 
562  ClusterID clusID = clusterIDs.at(cluster);
563  double completeness = counter->GetCompleteness(clusID);
564  double cleanliness = counter->GetCleanliness(clusID);
565  int numClusterHits = counter->GetNumberHitsInCluster(clusID);
566 
567  // Fill histograms for this cluster
568  hCompleteness->Fill(completeness, numClusterHits);
569  hCleanliness->Fill(cleanliness, numClusterHits);
570  hComplCleanl->Fill(completeness * cleanliness, numClusterHits);
571  hNumHitsCompleteness->Fill(completeness, numClusterHits);
572 
573  // Is this cluster doesn't correspond to a true particle continue
574  if (counter->IsNoise(clusID)) continue;
575 
576  double pi0Energy = GetPi0()->Momentum().E();
577  double pi0DecayAngle = FindPhotonAngle();
578  double conversionDistance =
579  GetEndTrackDistance(counter->GetTrack(clusID), (TrackID)GetPi0()->TrackId());
580 
581  hCompletenessEnergy->Fill(pi0Energy, completeness, numClusterHits);
582  hCompletenessAngle->Fill(pi0DecayAngle, completeness, numClusterHits);
583  hCompletenessConversionDistance->Fill(conversionDistance, completeness, numClusterHits);
584  hCleanlinessEnergy->Fill(pi0Energy, cleanliness, numClusterHits);
585  hCleanlinessAngle->Fill(pi0DecayAngle, cleanliness, numClusterHits);
586  hCleanlinessConversionDistance->Fill(conversionDistance, cleanliness, numClusterHits);
587  hComplCleanlEnergy->Fill(pi0Energy, cleanliness * completeness, numClusterHits);
588  hComplCleanlAngle->Fill(pi0DecayAngle, cleanliness * completeness, numClusterHits);
589  hComplCleanlConversionDistance->Fill(
590  conversionDistance, cleanliness * completeness, numClusterHits);
591  hNumHitsEnergy->Fill(pi0Energy, numClusterHits);
592 
593  // Continue if there are not two photons in the view
594  if (counter->GetPhotons().size() != 2) continue;
595 
596  double conversionSeparation =
597  GetEndTrackDistance(counter->GetPhotons()[0].first, counter->GetPhotons()[1].first);
598 
599  hCompletenessConversionSeparation->Fill(conversionSeparation, completeness, numClusterHits);
600  hCleanlinessConversionSeparation->Fill(conversionSeparation, cleanliness, numClusterHits);
601  }
602  }
603 }
604 
606 {
607 
608  // Average completeness/cleanliness
609  double avCompleteness = hCompleteness->GetMean();
610  double avCleanliness = hCleanliness->GetMean();
611 
612  // Write file
613  std::ofstream outFile("effpur");
614  outFile << avCompleteness << " " << avCleanliness;
615  outFile.close();
616 }
617 
618 // --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
620 public:
621  explicit ClusteringValidation(fhicl::ParameterSet const& p);
622 
623 private:
624  void analyze(art::Event const& evt);
625  void beginJob();
626  void endJob();
627 
628  // Clusterings to compare and the hits which the clustering was run over
629  std::vector<std::string> fClusterModuleLabels;
630  std::string fHitsModuleLabel;
631 
632  // Minimum hits needed to analyse a plane
634 
635  // Canvas on which to save histograms
636  TCanvas* fCanvas;
637 
638  // The cluster analysers
639  // -- one for each clustering type being compared
640  std::map<std::string, std::unique_ptr<ClusterAnalyser>> clusterAnalysis;
641 };
642 
644  : EDAnalyzer(pset)
645 {
646  fMinHitsInPlane = pset.get<int>("MinHitsInPlane");
647  fClusterModuleLabels = pset.get<std::vector<std::string>>("ClusterModuleLabels");
648  fHitsModuleLabel = pset.get<std::string>("HitsModuleLabel");
649 
650  fCanvas = new TCanvas("fCanvas", "", 800, 600);
651  gStyle->SetOptStat(0);
652 }
653 
655 {
656  // Get the hits from the event
658  std::vector<art::Ptr<recob::Hit>> hits;
659  if (evt.getByLabel(fHitsModuleLabel, hitHandle)) art::fill_ptr_vector(hits, hitHandle);
660 
661  // Get clustering information from event
662  // and give to the ClusterAnalyser to analyse
663  for (auto clustering : fClusterModuleLabels) {
664 
665  // Get the clusters from the event
667  std::vector<art::Ptr<recob::Cluster>> clusters;
668  if (evt.getByLabel(clustering, clusterHandle)) art::fill_ptr_vector(clusters, clusterHandle);
669 
670  // Find the associations (the hits) for the clusters
671  art::FindManyP<recob::Hit> fmh(clusterHandle, evt, clustering);
672 
673  // Analyse this particular clustering
674  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
675  clusterAnalysis.at(clustering)->Analyse(clockData, hits, clusters, fmh, fMinHitsInPlane);
676  }
677 }
678 
680 {
681 
682  // Construct the cluster analysers here
683  // -- one for each of the clustering types to compare
684  for (auto clustering : fClusterModuleLabels)
685  clusterAnalysis[clustering] = (std::unique_ptr<ClusterAnalyser>)new ClusterAnalyser(clustering);
686 }
687 
689 {
690 
691  // Make a map of all the histograms for each clustering
692  std::map<std::string, TObjArray> allHistograms;
693  for (auto clustering : fClusterModuleLabels)
694  allHistograms[clustering] = clusterAnalysis.at(clustering)->GetHistograms();
695 
696  // Write histograms to file
697  TFile* file = TFile::Open("validationHistograms.root", "RECREATE");
698  for (auto clustering : fClusterModuleLabels) {
699  file->cd();
700  file->mkdir(clustering.c_str());
701  file->cd(clustering.c_str());
702  for (int histIt = 0; histIt < allHistograms.begin()->second.GetEntriesFast(); ++histIt)
703  allHistograms.at(clustering).At(histIt)->Write();
704  }
705 
706  // Write images of overlaid histograms
707  for (int histIt = 0; histIt < allHistograms.begin()->second.GetEntriesFast(); ++histIt) {
708  fCanvas->cd();
709  fCanvas->Clear();
710  const char* name = allHistograms.begin()->second.At(histIt)->GetName();
711  TLegend* l = new TLegend(0.6, 0.8, 0.8, 0.9, name, "brNDC");
712  int clusterings = 1;
713  for (std::map<std::string, TObjArray>::iterator clusteringIt = allHistograms.begin();
714  clusteringIt != allHistograms.end();
715  ++clusteringIt, ++clusterings) {
716  TH1* h = (TH1*)allHistograms.at(clusteringIt->first).At(histIt);
717  h->SetLineColor(clusterings);
718  h->SetMarkerColor(clusterings);
719  if (clusterings == 1)
720  h->Draw();
721  else
722  h->Draw("same");
723  l->AddEntry(h, clusteringIt->first.c_str(), "lp");
724  }
725  l->Draw("same");
726  //fCanvas->SaveAs(name+TString(".png"));
727  file->cd();
728  fCanvas->Write(name);
729  }
730 
731  file->Close();
732  delete file;
733 
734  if (clusterAnalysis.find("blurredclusteringdc") != clusterAnalysis.end())
735  clusterAnalysis.at("blurredclusteringdc")->WriteFile();
736 }
737 
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.
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:364
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
float Integral() const
Integral under the calibrated signal waveform of the hit, in tick x ADC units.
Definition: Hit.h:254
geo::WireID const & WireID() const
Initial tdc tick for hit.
Definition: Hit.h:290
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:292
#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()
Interface for a class providing readout channel mapping to geometry.
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:373
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:315
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
The data type to uniquely identify a cryostat.
Definition: geo_types.h:187
std::map< std::string, std::unique_ptr< ClusterAnalyser > > clusterAnalysis
art::ServiceHandle< cheat::ParticleInventoryService > pi_serv