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();
97 ++numHitsPreClustering[trackID];
102 ++numSignalHitsPostClustering[clusID];
107 ++numNoiseHitsPostClustering[clusID];
113 clusterToTrackID[clusID] = trackID;
114 trackToClusterIDs[trackID].push_back(clusID);
119 return (
double)numSignalHitsPostClustering[clusID] /
120 (double)numHitsPreClustering[clusterToTrackID[clusID]];
125 return (
double)numSignalHitsPostClustering[clusID] / (double)(GetNumberHitsInCluster(clusID));
130 return 1 / (double)trackToClusterIDs.at(trackID).size();
135 return numHitsPreClustering[trackID];
140 return numSignalHitsPostClustering[clusID] + numNoiseHitsPostClustering[clusID];
146 for (
auto& trackHits : numHitsPreClustering)
147 nHits += trackHits.second;
155 i != clusterToTrackID.end();
157 v.push_back(i->first);
165 i != trackToClusterIDs.end();
167 v.push_back(i->first);
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))));
184 return clusterToTrackID.at(
id);
189 return IsNoise(clusterToTrackID.at(clusID));
194 return (
int)trackID == 0 ?
true :
false;
199 if (GetPhotons().
size() > 2 || GetPhotons().
size() == 0)
return false;
201 for (
unsigned int photon = 0; photon < GetPhotons().size(); ++photon)
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));
226 double FindPhotonAngle();
229 TObjArray GetHistograms();
230 void MakeHistograms();
269 std::map<unsigned int, std::map<unsigned int, std::unique_ptr<ClusterCounter>>>
clusterMap;
281 fClusterLabel = clusterLabel;
284 hCompleteness =
new TH1D(
"Completeness",
";Completeness;", 101, 0, 1.01);
285 hCompletenessEnergy =
286 new TProfile(
"CompletenessEnergy",
";True Energy (GeV);Completeness", 100, 0, 2);
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",
296 hCleanliness =
new TH1D(
"Cleanliness",
";Cleanliness;", 101, 0, 1.01);
298 new TProfile(
"CleanlinessEnergy",
";True Energy (GeV);Cleanliness", 100, 0, 2);
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",
316 hComplCleanlConversionSeparation =
317 new TProfile(
"CompletenessCleanlinessConversionSeparation",
318 ";True Conversion Separation (cm);Completeness * Cleanliness",
322 hPi0Energy =
new TH1D(
"Pi0EnergyCut",
";True Energy (GeV);", 25, 0, 2);
324 hPi0Angle =
new TH1D(
"Pi0AngleCut",
";True Angle (deg);", 25, 0, 180);
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);
345 new TH2D(
"NumHitsEnergy",
";True Energy (GeV);Size of Cluster", 100, 0, 2, 100, 0, 100);
357 auto const [tpc, plane] = std::make_tuple(
id.TPC,
id.
Plane);
358 clusterMap[tpc][plane] = std::make_unique<ClusterCounter>();
362 for (
size_t hitIt = 0; hitIt <
hits.size(); ++hitIt) {
364 TrackID trackID = FindTrackID(clockData, hit);
369 trueParticles.clear();
372 particleIt != particles.
end();
379 for (
size_t clusIt = 0; clusIt < clusters.size(); ++clusIt) {
382 unsigned int tpc = clusters.at(clusIt)->Plane().TPC % 2;
383 unsigned int plane = clusters.at(clusIt)->Plane().Plane;
387 if (clusterMap[tpc][plane]->GetNumberHitsInPlane() < minHits)
continue;
390 std::vector<art::Ptr<recob::Hit>> clusterHits = fmh.at(clusIt);
392 if (clusterHits.size() < 10)
continue;
395 TrackID trueTrackID = FindTrueTrack(clockData, clusterHits);
398 clusterMap[tpc][plane]->AssociateClusterAndTrack(
id, trueTrackID);
400 clusHitIt != clusterHits.end();
403 TrackID trackID = FindTrackID(clockData, hit);
404 if (trackID == trueTrackID)
405 clusterMap[tpc][plane]->AddSignalHitPostClustering(
id);
407 clusterMap[tpc][plane]->AddNoiseHitPostClustering(
id);
411 this->MakeHistograms();
418 double particleEnergy = 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);
427 return likelyTrackID;
434 std::map<TrackID, double> trackMap;
436 clusHitIt != clusterHits.end();
439 TrackID trackID = FindTrackID(clockData, hit);
440 trackMap[trackID] += hit->
Integral();
443 double highestCharge = 0;
447 if (trackIt->second > highestCharge) {
448 highestCharge = trackIt->second;
449 clusterTrack = trackIt->first;
460 .Angle(trueParticles.at((
TrackID)pi0->
Daughter(1))->Momentum().Vect()) *
461 (180 / TMath::Pi()));
469 particleIt != trueParticles.end();
471 if (particleIt->second->PdgCode() == 111) pi0 = particleIt->second;
479 trueParticles.at(id1)->EndPosition().X() - trueParticles.at(id2)->EndPosition().X(), 2) +
481 trueParticles.at(id1)->EndPosition().Y() - trueParticles.at(id2)->EndPosition().Y(), 2) +
483 trueParticles.at(id1)->EndPosition().Z() - trueParticles.at(id2)->EndPosition().Z(), 2));
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);
497 hEfficiencyEnergy->SetName(
"EfficiencyEnergy");
498 hEfficiencyAngle->SetName(
"EnergyAngle");
499 hEfficiencyConversionDistance->SetName(
"EfficiencyConversionDistance");
500 hEfficiencyConversionSeparation->SetName(
"EfficiencyConversionSeparation");
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);
533 auto const [tpc, plane] = std::make_tuple(
id.TPC,
id.
Plane);
534 auto&
counter = clusterMap[tpc][plane];
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));
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));
556 std::cout <<
"TPC " << tpc <<
", Plane " << plane <<
" fails the cut" << std::endl;
563 double completeness =
counter->GetCompleteness(clusID);
564 double cleanliness =
counter->GetCleanliness(clusID);
565 int numClusterHits =
counter->GetNumberHitsInCluster(clusID);
568 hCompleteness->Fill(completeness, numClusterHits);
569 hCleanliness->Fill(cleanliness, numClusterHits);
570 hComplCleanl->Fill(completeness * cleanliness, numClusterHits);
571 hNumHitsCompleteness->Fill(completeness, numClusterHits);
574 if (
counter->IsNoise(clusID))
continue;
576 double pi0Energy = GetPi0()->Momentum().E();
577 double pi0DecayAngle = FindPhotonAngle();
578 double conversionDistance =
579 GetEndTrackDistance(
counter->GetTrack(clusID), (
TrackID)GetPi0()->TrackId());
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);
594 if (
counter->GetPhotons().size() != 2)
continue;
596 double conversionSeparation =
597 GetEndTrackDistance(
counter->GetPhotons()[0].first,
counter->GetPhotons()[1].first);
599 hCompletenessConversionSeparation->Fill(conversionSeparation, completeness, numClusterHits);
600 hCleanlinessConversionSeparation->Fill(conversionSeparation, cleanliness, numClusterHits);
609 double avCompleteness = hCompleteness->GetMean();
610 double avCleanliness = hCleanliness->GetMean();
613 std::ofstream outFile(
"effpur");
614 outFile << avCompleteness <<
" " << avCleanliness;
646 fMinHitsInPlane = pset.
get<
int>(
"MinHitsInPlane");
647 fClusterModuleLabels = pset.
get<std::vector<std::string>>(
"ClusterModuleLabels");
648 fHitsModuleLabel = pset.
get<std::string>(
"HitsModuleLabel");
650 fCanvas =
new TCanvas(
"fCanvas",
"", 800, 600);
651 gStyle->SetOptStat(0);
658 std::vector<art::Ptr<recob::Hit>>
hits;
663 for (
auto clustering : fClusterModuleLabels) {
667 std::vector<art::Ptr<recob::Cluster>> clusters;
675 clusterAnalysis.at(clustering)->Analyse(clockData, hits, clusters, fmh, fMinHitsInPlane);
684 for (
auto clustering : fClusterModuleLabels)
685 clusterAnalysis[clustering] = (std::unique_ptr<ClusterAnalyser>)
new ClusterAnalyser(clustering);
692 std::map<std::string, TObjArray> allHistograms;
693 for (
auto clustering : fClusterModuleLabels)
694 allHistograms[clustering] = clusterAnalysis.at(clustering)->GetHistograms();
697 TFile*
file = TFile::Open(
"validationHistograms.root",
"RECREATE");
698 for (
auto clustering : fClusterModuleLabels) {
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();
707 for (
int histIt = 0; histIt < allHistograms.begin()->second.GetEntriesFast(); ++histIt) {
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");
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)
723 l->AddEntry(h, clusteringIt->first.c_str(),
"lp");
728 fCanvas->Write(name);
734 if (clusterAnalysis.find(
"blurredclusteringdc") != clusterAnalysis.end())
735 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.
list_type::const_iterator const_iterator
double GetEndTrackDistance(TrackID id1, TrackID id2)
TProfile * hCompletenessConversionSeparation
The data type to uniquely identify a Plane.
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
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()
Interface for a class providing readout channel mapping to geometry.
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()
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