37 #include "TEfficiency.h" 64 void AddHitPreClustering(
TrackID id);
65 void AddSignalHitPostClustering(
ClusterID id);
66 void AddNoiseHitPostClustering(
ClusterID id);
70 double GetEfficiency(
TrackID id);
73 int GetNumberHitsFromTrack(
TrackID id);
75 int GetNumberHitsInPlane();
76 std::vector<std::pair<TrackID, ClusterIDs>> GetPhotons();
98 ++numHitsPreClustering[trackID];
103 ++numSignalHitsPostClustering[clusID];
108 ++numNoiseHitsPostClustering[clusID];
114 clusterToTrackID[clusID] = trackID;
115 trackToClusterIDs[trackID].push_back(clusID);
120 return (
double)numSignalHitsPostClustering[clusID] /
121 (double)numHitsPreClustering[clusterToTrackID[clusID]];
126 return (
double)numSignalHitsPostClustering[clusID] / (double)(GetNumberHitsInCluster(clusID));
131 return 1 / (double)trackToClusterIDs.at(trackID).size();
136 return numHitsPreClustering[trackID];
141 return numSignalHitsPostClustering[clusID] + numNoiseHitsPostClustering[clusID];
147 for (
auto& trackHits : numHitsPreClustering)
148 nHits += trackHits.second;
156 i != clusterToTrackID.end();
158 v.push_back(i->first);
166 i != trackToClusterIDs.end();
168 v.push_back(i->first);
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))));
185 return clusterToTrackID.at(
id);
190 return IsNoise(clusterToTrackID.at(clusID));
195 return (
int)trackID == 0 ?
true :
false;
200 if (GetPhotons().
size() > 2 || GetPhotons().
size() == 0)
return false;
202 for (
unsigned int photon = 0; photon < GetPhotons().size(); ++photon)
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));
227 double FindPhotonAngle();
230 TObjArray GetHistograms();
231 void MakeHistograms();
270 std::map<unsigned int, std::map<unsigned int, std::unique_ptr<ClusterCounter>>>
clusterMap;
282 fClusterLabel = clusterLabel;
285 hCompleteness =
new TH1D(
"Completeness",
";Completeness;", 101, 0, 1.01);
286 hCompletenessEnergy =
287 new TProfile(
"CompletenessEnergy",
";True Energy (GeV);Completeness", 100, 0, 2);
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",
297 hCleanliness =
new TH1D(
"Cleanliness",
";Cleanliness;", 101, 0, 1.01);
299 new TProfile(
"CleanlinessEnergy",
";True Energy (GeV);Cleanliness", 100, 0, 2);
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",
317 hComplCleanlConversionSeparation =
318 new TProfile(
"CompletenessCleanlinessConversionSeparation",
319 ";True Conversion Separation (cm);Completeness * Cleanliness",
323 hPi0Energy =
new TH1D(
"Pi0EnergyCut",
";True Energy (GeV);", 25, 0, 2);
325 hPi0Angle =
new TH1D(
"Pi0AngleCut",
";True Angle (deg);", 25, 0, 180);
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);
346 new TH2D(
"NumHitsEnergy",
";True Energy (GeV);Size of Cluster", 100, 0, 2, 100, 0, 100);
358 auto const [tpc, plane] = std::make_tuple(
id.TPC,
id.
Plane);
359 clusterMap[tpc][plane] = std::make_unique<ClusterCounter>();
363 for (
size_t hitIt = 0; hitIt <
hits.size(); ++hitIt) {
365 TrackID trackID = FindTrackID(clockData, hit);
370 trueParticles.clear();
373 particleIt != particles.
end();
380 for (
size_t clusIt = 0; clusIt < clusters.size(); ++clusIt) {
383 unsigned int tpc = clusters.at(clusIt)->Plane().TPC % 2;
384 unsigned int plane = clusters.at(clusIt)->Plane().Plane;
388 if (clusterMap[tpc][plane]->GetNumberHitsInPlane() < minHits)
continue;
391 std::vector<art::Ptr<recob::Hit>> clusterHits = fmh.at(clusIt);
393 if (clusterHits.size() < 10)
continue;
396 TrackID trueTrackID = FindTrueTrack(clockData, clusterHits);
399 clusterMap[tpc][plane]->AssociateClusterAndTrack(
id, trueTrackID);
401 clusHitIt != clusterHits.end();
404 TrackID trackID = FindTrackID(clockData, hit);
405 if (trackID == trueTrackID)
406 clusterMap[tpc][plane]->AddSignalHitPostClustering(
id);
408 clusterMap[tpc][plane]->AddNoiseHitPostClustering(
id);
412 this->MakeHistograms();
419 double particleEnergy = 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);
428 return likelyTrackID;
435 std::map<TrackID, double> trackMap;
437 clusHitIt != clusterHits.end();
440 TrackID trackID = FindTrackID(clockData, hit);
441 trackMap[trackID] += hit->
Integral();
444 double highestCharge = 0;
448 if (trackIt->second > highestCharge) {
449 highestCharge = trackIt->second;
450 clusterTrack = trackIt->first;
461 .Angle(trueParticles.at((
TrackID)pi0->
Daughter(1))->Momentum().Vect()) *
462 (180 / TMath::Pi()));
470 particleIt != trueParticles.end();
472 if (particleIt->second->PdgCode() == 111) pi0 = particleIt->second;
480 trueParticles.at(id1)->EndPosition().X() - trueParticles.at(id2)->EndPosition().X(), 2) +
482 trueParticles.at(id1)->EndPosition().Y() - trueParticles.at(id2)->EndPosition().Y(), 2) +
484 trueParticles.at(id1)->EndPosition().Z() - trueParticles.at(id2)->EndPosition().Z(), 2));
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);
498 hEfficiencyEnergy->SetName(
"EfficiencyEnergy");
499 hEfficiencyAngle->SetName(
"EnergyAngle");
500 hEfficiencyConversionDistance->SetName(
"EfficiencyConversionDistance");
501 hEfficiencyConversionSeparation->SetName(
"EfficiencyConversionSeparation");
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);
534 auto const [tpc, plane] = std::make_tuple(
id.TPC,
id.
Plane);
535 auto&
counter = clusterMap[tpc][plane];
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));
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));
557 std::cout <<
"TPC " << tpc <<
", Plane " << plane <<
" fails the cut" << std::endl;
564 double completeness =
counter->GetCompleteness(clusID);
565 double cleanliness =
counter->GetCleanliness(clusID);
566 int numClusterHits =
counter->GetNumberHitsInCluster(clusID);
569 hCompleteness->Fill(completeness, numClusterHits);
570 hCleanliness->Fill(cleanliness, numClusterHits);
571 hComplCleanl->Fill(completeness * cleanliness, numClusterHits);
572 hNumHitsCompleteness->Fill(completeness, numClusterHits);
575 if (
counter->IsNoise(clusID))
continue;
577 double pi0Energy = GetPi0()->Momentum().E();
578 double pi0DecayAngle = FindPhotonAngle();
579 double conversionDistance =
580 GetEndTrackDistance(
counter->GetTrack(clusID), (
TrackID)GetPi0()->TrackId());
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);
595 if (
counter->GetPhotons().size() != 2)
continue;
597 double conversionSeparation =
598 GetEndTrackDistance(
counter->GetPhotons()[0].first,
counter->GetPhotons()[1].first);
600 hCompletenessConversionSeparation->Fill(conversionSeparation, completeness, numClusterHits);
601 hCleanlinessConversionSeparation->Fill(conversionSeparation, cleanliness, numClusterHits);
610 double avCompleteness = hCompleteness->GetMean();
611 double avCleanliness = hCleanliness->GetMean();
614 std::ofstream outFile(
"effpur");
615 outFile << avCompleteness <<
" " << avCleanliness;
647 fMinHitsInPlane = pset.
get<
int>(
"MinHitsInPlane");
648 fClusterModuleLabels = pset.
get<std::vector<std::string>>(
"ClusterModuleLabels");
649 fHitsModuleLabel = pset.
get<std::string>(
"HitsModuleLabel");
651 fCanvas =
new TCanvas(
"fCanvas",
"", 800, 600);
652 gStyle->SetOptStat(0);
659 std::vector<art::Ptr<recob::Hit>>
hits;
664 for (
auto clustering : fClusterModuleLabels) {
668 std::vector<art::Ptr<recob::Cluster>> clusters;
676 clusterAnalysis.at(clustering)->Analyse(clockData, hits, clusters, fmh, fMinHitsInPlane);
685 for (
auto clustering : fClusterModuleLabels)
686 clusterAnalysis[clustering] = (std::unique_ptr<ClusterAnalyser>)
new ClusterAnalyser(clustering);
693 std::map<std::string, TObjArray> allHistograms;
694 for (
auto clustering : fClusterModuleLabels)
695 allHistograms[clustering] = clusterAnalysis.at(clustering)->GetHistograms();
698 TFile*
file = TFile::Open(
"validationHistograms.root",
"RECREATE");
699 for (
auto clustering : fClusterModuleLabels) {
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();
708 for (
int histIt = 0; histIt < allHistograms.begin()->second.GetEntriesFast(); ++histIt) {
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");
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)
724 l->AddEntry(h, clusteringIt->first.c_str(),
"lp");
729 fCanvas->Write(name);
735 if (clusterAnalysis.find(
"blurredclusteringdc") != clusterAnalysis.end())
736 clusterAnalysis.at(
"blurredclusteringdc")->WriteFile();
void AddHitPreClustering(TrackID id)
art::ServiceHandle< cheat::BackTrackerService const > bt_serv
TProfile * hCleanlinessConversionDistance
int GetNumberHitsFromTrack(TrackID id)
TProfile * hCleanlinessConversionSeparation
bool IsNoise(ClusterID id)
std::map< ClusterID, int > numNoiseHitsPostClustering
TProfile * hComplCleanlConversionDistance
double GetCompleteness(ClusterID id)
TProfile * hComplCleanlEnergy
TProfile * hCompletenessAngle
std::string fClusterLabel
art::ServiceHandle< cheat::BackTrackerService const > bt_serv
ClusteringValidation(fhicl::ParameterSet const &p)
TH1 * hPi0ConversionSeparation
std::vector< TrackID > TrackIDs
Declaration of signal hit object.
art::ServiceHandle< geo::Geometry const > geometry
list_type::const_iterator const_iterator
double GetEndTrackDistance(TrackID id1, TrackID id2)
TProfile * hCompletenessConversionSeparation
The data type to uniquely identify a Plane.
float Integral() const
Integral under the calibrated signal waveform of the hit, in tick x ADC units.
double GetEfficiency(TrackID id)
geo::WireID const & WireID() const
Initial tdc tick for hit.
std::map< TrackID, ClusterIDs > trackToClusterIDs
std::map< TrackID, std::map< std::string, double > > particleProperties
std::string fHitsModuleLabel
TEfficiency * hEfficiencyAngle
Cluster finding and building.
std::map< TrackID, int > numHitsPreClustering
TrackIDs GetListOfTrackIDs()
int NumberDaughters() const
ClusterAnalyser(std::string &label)
int Daughter(const int i) const
std::map< ClusterID, int > numSignalHitsPostClustering
TEfficiency * hEfficiencyConversionSeparation
int GetNumberHitsInPlane()
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
void AddNoiseHitPostClustering(ClusterID id)
TrackID FindTrackID(detinfo::DetectorClocksData const &clockData, art::Ptr< recob::Hit > &hit)
auto counter(T begin, T end)
Returns an object to iterate values from begin to end in a range-for loop.
#define DEFINE_ART_MODULE(klass)
std::vector< std::pair< TrackID, ClusterIDs > > GetPhotons()
T get(std::string const &key) const
TProfile * hCleanlinessAngle
art::ServiceHandle< cheat::ParticleInventoryService const > pi_serv
TrackID FindTrueTrack(detinfo::DetectorClocksData const &clockData, std::vector< art::Ptr< recob::Hit >> &clusterHits)
std::vector< ClusterID > ClusterIDs
void AddSignalHitPostClustering(ClusterID id)
PlaneID_t Plane
Index of the plane within its TPC.
Declaration of cluster object.
std::map< unsigned int, std::map< unsigned int, std::unique_ptr< ClusterCounter > > > clusterMap
Detector simulation of raw signals on wires.
TProfile * hComplCleanlConversionSeparation
TProfile * hCleanlinessEnergy
std::vector< std::string > fClusterModuleLabels
TProfile * hCompletenessEnergy
double GetCleanliness(ClusterID id)
const simb::MCParticle * GetPi0()
TEfficiency * hEfficiencyEnergy
TH1 * hPi0ConversionDistanceCut
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.
TrackID GetTrack(ClusterID id)
std::map< TrackID, simb::MCParticle > trueParticles
void AssociateClusterAndTrack(ClusterID clusID, TrackID trackID)
TH2 * hNumHitsCompleteness
TObjArray GetHistograms()
TProfile * hCompletenessConversionDistance
void analyze(art::Event const &evt)
TH1 * hPi0ConversionSeparationCut
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
TPCID_t TPC
Index of the TPC within its cryostat.
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.
recob::tracking::Plane Plane
TH1 * hPi0ConversionDistance
std::map< ClusterID, TrackID > clusterToTrackID
ClusterIDs GetListOfClusterIDs()
art framework interface to geometry description
TEfficiency * hEfficiencyConversionDistance
The data type to uniquely identify a cryostat.
std::map< std::string, std::unique_ptr< ClusterAnalyser > > clusterAnalysis
int GetNumberHitsInCluster(ClusterID id)
art::ServiceHandle< cheat::ParticleInventoryService > pi_serv
TProfile * hComplCleanlAngle
art::ServiceHandle< geo::Geometry const > geometry