From c4f3e3fc953417500916b973be27f63f46cadf80 Mon Sep 17 00:00:00 2001 From: Asa Nehm Date: Tue, 30 Jul 2024 10:02:33 -0500 Subject: [PATCH 1/7] Translating Xiaoyan's python code for this into C++/ROOT and calling it in TMS_Reco --- src/TMS_Reco.cpp | 3 +++ src/TMS_Reco.h | 3 +++ 2 files changed, 6 insertions(+) diff --git a/src/TMS_Reco.cpp b/src/TMS_Reco.cpp index d08a4e44..9c8ca351 100644 --- a/src/TMS_Reco.cpp +++ b/src/TMS_Reco.cpp @@ -1264,6 +1264,9 @@ std::vector TMS_TrackFinder::TrackMatching3D() { // Sort track SpatialPrio(aTrack.Hits); + // Charge ID + aTrack.Charge = ChargeID.ID_Track_Charge(aTrack.Hits); + // Track Length aTrack.Length = CalculateTrackLength3D(aTrack); #ifdef DEBUG diff --git a/src/TMS_Reco.h b/src/TMS_Reco.h index f3c93e1a..c413363b 100644 --- a/src/TMS_Reco.h +++ b/src/TMS_Reco.h @@ -24,6 +24,8 @@ #include "TMS_Track.h" +#include "TMS_ChargeID.h" + // Hand over to the Kalman reconstruction once we find tracks #include "TMS_Kalman.h" @@ -312,6 +314,7 @@ class TMS_TrackFinder { TMS_Kalman KalmanFitter; TMS_DBScan DBSCAN; + TMS_ChargeID ChargeID;//ID_Track_Charge(const std::vector &Hits); int FindBin(double Rho); // The candidates for each particle From 234b98dc042097f6aba7181c0c87d6d16595f2cf Mon Sep 17 00:00:00 2001 From: Asa Nehm Date: Tue, 30 Jul 2024 10:02:50 -0500 Subject: [PATCH 2/7] Translating Xiaoyan's python code for this into C++/ROOT and calling it in TMS_Reco --- src/TMS_ChargeID.cpp | 163 +++++++++++++++++++++++++++++++++++++++++++ src/TMS_ChargeID.h | 24 +++++++ 2 files changed, 187 insertions(+) create mode 100644 src/TMS_ChargeID.cpp create mode 100644 src/TMS_ChargeID.h diff --git a/src/TMS_ChargeID.cpp b/src/TMS_ChargeID.cpp new file mode 100644 index 00000000..9d672ee8 --- /dev/null +++ b/src/TMS_ChargeID.cpp @@ -0,0 +1,163 @@ +#include "TMS_ChargeID.h" + +//TODO??? TMS_ChargeID::TMS_ChargeID() : + +// Check if hit is in the far negative region of x (detector) +bool TMS_ChargeID::region1(const TMS_Hit &Hit) { + bool is_region1 = true; + if (Hit.GetRecoX() < -3520.0 || Hit.GetRecoX() > -1750.0) is_region1 = false; //TODO these numbers should be called in from TMS_Constants.h once they exist there + return is_region1; +} + +// Check if hit is in the middle region of x (detector) +bool TMS_ChargeID::region2(const TMS_Hit &Hit) { + bool is_region2 = true; + if (Hit.GetRecoX() < -1750.0 || Hit.GetRecoX() > 1750.0) is_region2 = false; //TODO these numbers should be called in from TMS_Constants.h once they exist there + return is_region2; +} + +// Check if hit is in the far positive region of x (detector) +bool TMS_ChargeID::region3(const TMS_Hit &Hit) { + bool is_region3 = true; + if (Hit.GetRecoX() < 1750.0 || Hit.GetRecoX() > 3520.0) is_region3 = false; //TODO these numbers should be called in from TMS_Constants.h once they exist there + return is_region3; +} + +// This takes the hit positions of a track as output and outputs a number representing the charge +// Output number 13 -> muon, -13 -> antimuon, -999999999 -> particle cannot be identified by this method +int TMS_ChargeID::ID_Track_Charge(const std::vector &Track) { + // Initialize some local variables and empty arrays/vectors to be used + int chargeID = -999999999; + std::vector muon_hits = *Track; + std::vector region1_hits; + std::vector region2_hits; + std::vector region3_hits; + std::vector::iterator n_changeregion = muon_hits.begin(); + int n_plus = 0; + int n_minus = 0; + + // Check if there are any hits in the track + if (muon_hits.size() < 3) return chargeID; + + // Check the muon starting region + // below are the algorithms that cut the muon hits that are in different regions, + // and store them separately to do the calculations afterwards + SpatialPrio(muon_hits); + + // If muon starts in region 1 + if (region1(muon_hits.front())) { + for (std::vector::iterator hit = muon_hits.begin(); hit != muon_hits.end(); ++hit) { + // Check it the muon is still in the same region as the starting region. + // Mark down the point where the muon crosses into a different region, + // or if there is an invalid entry. + if (!region1(muon_hits[*hit]) || muon_hits[*hit] == -999) { + n_changeregion = *hit; + break; + } + // Now we have only the hits in region 1 + region1_hits.push_back(*hit); + } + if (region2(muon_hits[*n_changeregion])) { + // Fill region 2 hits if there are any + for (std::vector::iterator it = muon_hits[*n_changeregion]; it != muon_hits.end(); ++it) { + if (!region2(muon_hits[*it]) || muon_hits[*it] == -999) break; + region2_hits.push_back(*it); + } + } + } + // If muon starts in region 3, similar to the above algorithm + if (region3(muon_hits.front())) { + for (std::vector::iterator hit = muon_hits.begin(); hit != muon_hits.end(); ++hit) { + if (!region3(muon_hits[*hit]) || muon_hits[*hit] == -999) { + n_changeregion = *hit; + break; + } + region3_hits.push_back(*hit); + } + if (region2(muon_hits[*n_changeregion])) { + for (std::vector::iterator it = muon_hits[*n_changeregion]; it != muon_hits.end(); ++it) { + if (!region2(muon_hits[*it]) || muon_hits[*it] == -999) break; + region2_hits.push_back(*it); + } + } + } + // If muon starts in region 2, similar to the above algorithm + if (region2(muon_hits.front())) { + for (std::vector::iterator hit = muon_hits.begin(); hit != muon_hits.end(); ++hit) { + if (!region2(muon_hits[*hit]) || muon_hits[*hit] == -999) { + n_changeregion = *hit; + break; + } + region2_hits.push_back(*hit); + } + if (region1(muon_hits[*n_changeregion])) { + for (std::vector::iterator it = muon_hits[*n_changeregion]; it != muon_hits.end(); ++it) { + if (!region1(muon_hits[*it]) || muon_hits[*it] == -999) break; + region1_hits.push_back(*it); + } + } + if (region3(muon_hits[*n_changeregion])) { + for (std::vector::iterator it = muon_hits[*n_changeregion]; it != muon_hits.end(); ++it) { + if (!region3(muon_hits[*it]) || muon_hits[*it] == -999) break; + region3_hits.push_back(*it); + } + } + } + + // Check how many hits there are, three entries means one hit TODO???? + int total_hit_region1 = region1_hits.size(); + int total_hit_region2 = region2_hits.size(); + int total_hit_region3 = region3_hits.size(); + + // Now the muon hits are collected in three different regions, do the calculation + // The calculations use the following method: draw a line from the first to the last hit + // and then count if there are more hits to the right of the line + // or more to the left of the line, to decide which direction the muon is bending + + // Region 1 calculation + if (total_hit_region1 > 2 && region1_hits.front().GetZ() < region1_hits.back().GetZ()) { + double m = (region1_hits.front().GetRecoX() - region1_hits.back().GetRecoX()) / (region1_hits.front().GetZ() - region1_hits.back().GetZ()); + for (int i = 0; i != total_hit_region1; ++i) { + double x_interpolation = m * (region1_hits[i].GetZ() - region1_hits.front().GetZ()) + region1_hits.front().GetRecoX(); + double signed_dist = region1_hits[i].GetRecoX() - x_interpolation; + if (signed_dist > 0) n_plus += 1; + if (signed_dist < 0) n_minus += 1; + } + } + + // Region 2 calculation + if (total_hit_region2 > 2 && region2_hits.front().GetZ() < region2_hits.back().GetZ()) { + double m = (region2_hits.front().GetRecoX() - region2_hits.back().GetRecoX()) / (region2_hits.front().GetZ() - region2_hits.back().GetZ()); + for (int i = 0; i != total_hit_region2; ++i) { + double x_interpolation = m * (region2_hits[i].GetZ() - region2_hits.front().GetZ()) + region2_hits.front().GetRecoX(); + double signed_dist = region2_hits[i].GetRecoX() - x_interpolation; + if (signed_dist < 0) n_plus += 1; + if (signed_dist > 0) n_minus += 1; + } + } + + // Region 3 calculation + if (total_hit_region3 > 2 && region3_hits.front().GetZ() < region3_hits.back().GetZ()) { + double m = (region3_hits.front().GetRecoX() - region3_hits.back().GetRecoX()) / (region3_hits.front().GetZ() - region3_hits.back().GetZ()); + for (int i = 0; i != total_hit_region3; ++i) { + double x_interpolation = m * (region3_hits[i].GetZ() - region3_hits.front().GetZ()) + region3_hits.front().GetRecoX(); + double signed_dist = region3_hits[i].GetRecoX() - x_interpolation; + if (signed_dist > 0) n_plus += 1; + if (signed_dist < 0) n_minus += 1; + } + } + + // Now the calculation is done, identify whether this particle is a muon or an antimuon. + // After the identificationon you can assign the chargeID to be 13 or -13 + // depending whether it is a muon or antimuon + // if no calculation can be done then the chargeID here is defaulted to be -999999999, + // which means this particle cannot be identified with this method + + // Now judge whether the particle is a muon or an antimuon + if (n_plus < n_minus) chargeID = 13; + if (n_plus > n_minus) chargeID = -13; + + return chargeID; +} + + diff --git a/src/TMS_ChargeID.h b/src/TMS_ChargeID.h new file mode 100644 index 00000000..be441e9d --- /dev/null +++ b/src/TMS_ChargeID.h @@ -0,0 +1,24 @@ +#ifndef _TMS_CHARGEID_H_SEEN_ +#define _TMS_CHARGEID_H_SEEN_ + +#include "TMS_Hit.h" +#include "TMS_Constants.h" + +//Charge identification class +class TMS_ChargeID { + public: + TMS_ChargeID(); + // Functions to determine in which magnetic field orientation the track is + bool region1(const TMS_Hit &Hit); + bool region2(const TMS_Hit &Hit); + bool region3(const TMS_Hit &Hit); + + // Function for the actual charge identification + int ID_Track_Charge(const std::vector &Track); + + private: + TMS_Hit Hit; + std::vector Track; +}; + +#endif From ba09f3f1c570357a7527d4effec8a908fabb8cd3 Mon Sep 17 00:00:00 2001 From: Asa Nehm Date: Tue, 30 Jul 2024 10:30:01 -0500 Subject: [PATCH 3/7] Correcting a small mistake from copying --- src/TMS_ChargeID.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/TMS_ChargeID.cpp b/src/TMS_ChargeID.cpp index 9d672ee8..20207c8f 100644 --- a/src/TMS_ChargeID.cpp +++ b/src/TMS_ChargeID.cpp @@ -117,7 +117,7 @@ int TMS_ChargeID::ID_Track_Charge(const std::vector &Track) { // Region 1 calculation if (total_hit_region1 > 2 && region1_hits.front().GetZ() < region1_hits.back().GetZ()) { double m = (region1_hits.front().GetRecoX() - region1_hits.back().GetRecoX()) / (region1_hits.front().GetZ() - region1_hits.back().GetZ()); - for (int i = 0; i != total_hit_region1; ++i) { + for (int i = 1; i != (total_hit_region1 - 1); ++i) { double x_interpolation = m * (region1_hits[i].GetZ() - region1_hits.front().GetZ()) + region1_hits.front().GetRecoX(); double signed_dist = region1_hits[i].GetRecoX() - x_interpolation; if (signed_dist > 0) n_plus += 1; @@ -128,7 +128,7 @@ int TMS_ChargeID::ID_Track_Charge(const std::vector &Track) { // Region 2 calculation if (total_hit_region2 > 2 && region2_hits.front().GetZ() < region2_hits.back().GetZ()) { double m = (region2_hits.front().GetRecoX() - region2_hits.back().GetRecoX()) / (region2_hits.front().GetZ() - region2_hits.back().GetZ()); - for (int i = 0; i != total_hit_region2; ++i) { + for (int i = 1; i != (total_hit_region2 - 1); ++i) { double x_interpolation = m * (region2_hits[i].GetZ() - region2_hits.front().GetZ()) + region2_hits.front().GetRecoX(); double signed_dist = region2_hits[i].GetRecoX() - x_interpolation; if (signed_dist < 0) n_plus += 1; @@ -139,7 +139,7 @@ int TMS_ChargeID::ID_Track_Charge(const std::vector &Track) { // Region 3 calculation if (total_hit_region3 > 2 && region3_hits.front().GetZ() < region3_hits.back().GetZ()) { double m = (region3_hits.front().GetRecoX() - region3_hits.back().GetRecoX()) / (region3_hits.front().GetZ() - region3_hits.back().GetZ()); - for (int i = 0; i != total_hit_region3; ++i) { + for (int i = 1; i != (total_hit_region3 - 1); ++i) { double x_interpolation = m * (region3_hits[i].GetZ() - region3_hits.front().GetZ()) + region3_hits.front().GetRecoX(); double signed_dist = region3_hits[i].GetRecoX() - x_interpolation; if (signed_dist > 0) n_plus += 1; From 2ee7e762fdd86ef47e808c457993633c0e972a03 Mon Sep 17 00:00:00 2001 From: Asa Nehm Date: Thu, 1 Aug 2024 08:07:46 -0500 Subject: [PATCH 4/7] Fixing vector iterator problem and adding TMS_Charge in Makefile --- src/Makefile | 2 +- src/TMS_ChargeID.cpp | 40 +++++++++++++++++++--------------------- 2 files changed, 20 insertions(+), 22 deletions(-) diff --git a/src/Makefile b/src/Makefile index f6d97aef..43232f01 100644 --- a/src/Makefile +++ b/src/Makefile @@ -55,7 +55,7 @@ MKDIR_P := mkdir -p # Library directory where we build to LIB_DIR=../lib -TMS_OBJ = TMS_Bar.o TMS_Hit.o TMS_TrueHit.o TMS_Event.o TMS_Track.o TMS_TrueParticle.o TMS_EventViewer.o TMS_Reco.o TMS_Kalman.o TMS_TreeWriter.o TMS_ReadoutTreeWriter.o TMS_Manager.o TMS_Readout_Manager.o BField_Handler.o TMS_Utils.o +TMS_OBJ = TMS_Bar.o TMS_Hit.o TMS_TrueHit.o TMS_Event.o TMS_Track.o TMS_TrueParticle.o TMS_EventViewer.o TMS_Reco.o TMS_ChargeID.o TMS_Kalman.o TMS_TreeWriter.o TMS_ReadoutTreeWriter.o TMS_Manager.o TMS_Readout_Manager.o BField_Handler.o TMS_Utils.o # Our lovely collection of libs LIB_OBJ = $(EDEP_LIBS) $(ROOT_LIBS) diff --git a/src/TMS_ChargeID.cpp b/src/TMS_ChargeID.cpp index 20207c8f..0b8eb49c 100644 --- a/src/TMS_ChargeID.cpp +++ b/src/TMS_ChargeID.cpp @@ -28,7 +28,7 @@ bool TMS_ChargeID::region3(const TMS_Hit &Hit) { int TMS_ChargeID::ID_Track_Charge(const std::vector &Track) { // Initialize some local variables and empty arrays/vectors to be used int chargeID = -999999999; - std::vector muon_hits = *Track; + std::vector muon_hits = Track; std::vector region1_hits; std::vector region2_hits; std::vector region3_hits; @@ -42,25 +42,23 @@ int TMS_ChargeID::ID_Track_Charge(const std::vector &Track) { // Check the muon starting region // below are the algorithms that cut the muon hits that are in different regions, // and store them separately to do the calculations afterwards - SpatialPrio(muon_hits); - // If muon starts in region 1 if (region1(muon_hits.front())) { for (std::vector::iterator hit = muon_hits.begin(); hit != muon_hits.end(); ++hit) { // Check it the muon is still in the same region as the starting region. // Mark down the point where the muon crosses into a different region, // or if there is an invalid entry. - if (!region1(muon_hits[*hit]) || muon_hits[*hit] == -999) { - n_changeregion = *hit; + if (!region1(*hit) || (*hit).GetRecoX() == -999) { + n_changeregion = hit; break; } // Now we have only the hits in region 1 region1_hits.push_back(*hit); } - if (region2(muon_hits[*n_changeregion])) { + if (region2(*n_changeregion)) { // Fill region 2 hits if there are any - for (std::vector::iterator it = muon_hits[*n_changeregion]; it != muon_hits.end(); ++it) { - if (!region2(muon_hits[*it]) || muon_hits[*it] == -999) break; + for (std::vector::iterator it = n_changeregion; it != muon_hits.end(); ++it) { + if (!region2(*it) || (*it).GetRecoX() == -999) break; region2_hits.push_back(*it); } } @@ -68,15 +66,15 @@ int TMS_ChargeID::ID_Track_Charge(const std::vector &Track) { // If muon starts in region 3, similar to the above algorithm if (region3(muon_hits.front())) { for (std::vector::iterator hit = muon_hits.begin(); hit != muon_hits.end(); ++hit) { - if (!region3(muon_hits[*hit]) || muon_hits[*hit] == -999) { - n_changeregion = *hit; + if (!region3(*hit) || (*hit).GetRecoX() == -999) { + n_changeregion = hit; break; } region3_hits.push_back(*hit); } - if (region2(muon_hits[*n_changeregion])) { - for (std::vector::iterator it = muon_hits[*n_changeregion]; it != muon_hits.end(); ++it) { - if (!region2(muon_hits[*it]) || muon_hits[*it] == -999) break; + if (region2(*n_changeregion)) { + for (std::vector::iterator it = n_changeregion; it != muon_hits.end(); ++it) { + if (!region2(*it) || (*it).GetRecoX() == -999) break; region2_hits.push_back(*it); } } @@ -84,21 +82,21 @@ int TMS_ChargeID::ID_Track_Charge(const std::vector &Track) { // If muon starts in region 2, similar to the above algorithm if (region2(muon_hits.front())) { for (std::vector::iterator hit = muon_hits.begin(); hit != muon_hits.end(); ++hit) { - if (!region2(muon_hits[*hit]) || muon_hits[*hit] == -999) { - n_changeregion = *hit; + if (!region2(*hit) || (*hit).GetRecoX() == -999) { + n_changeregion = hit; break; } region2_hits.push_back(*hit); } - if (region1(muon_hits[*n_changeregion])) { - for (std::vector::iterator it = muon_hits[*n_changeregion]; it != muon_hits.end(); ++it) { - if (!region1(muon_hits[*it]) || muon_hits[*it] == -999) break; + if (region1(*n_changeregion)) { + for (std::vector::iterator it = n_changeregion; it != muon_hits.end(); ++it) { + if (!region1(*it) || (*it).GetRecoX() == -999) break; region1_hits.push_back(*it); } } - if (region3(muon_hits[*n_changeregion])) { - for (std::vector::iterator it = muon_hits[*n_changeregion]; it != muon_hits.end(); ++it) { - if (!region3(muon_hits[*it]) || muon_hits[*it] == -999) break; + if (region3(*n_changeregion)) { + for (std::vector::iterator it = n_changeregion; it != muon_hits.end(); ++it) { + if (!region3(*it) || (*it).GetRecoX() == -999) break; region3_hits.push_back(*it); } } From af38ede706d09f57f99f40b7ccb20824fe69acdb Mon Sep 17 00:00:00 2001 From: Jeffrey Kleykamp Date: Thu, 1 Aug 2024 09:23:56 -0500 Subject: [PATCH 5/7] Fix undefined reference bug --- src/TMS_ChargeID.cpp | 4 +++- src/TMS_ChargeID.h | 5 +++-- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/src/TMS_ChargeID.cpp b/src/TMS_ChargeID.cpp index 0b8eb49c..e1386738 100644 --- a/src/TMS_ChargeID.cpp +++ b/src/TMS_ChargeID.cpp @@ -1,6 +1,8 @@ #include "TMS_ChargeID.h" -//TODO??? TMS_ChargeID::TMS_ChargeID() : +TMS_ChargeID::TMS_ChargeID() { + +} // Check if hit is in the far negative region of x (detector) bool TMS_ChargeID::region1(const TMS_Hit &Hit) { diff --git a/src/TMS_ChargeID.h b/src/TMS_ChargeID.h index be441e9d..7f4e765a 100644 --- a/src/TMS_ChargeID.h +++ b/src/TMS_ChargeID.h @@ -17,8 +17,9 @@ class TMS_ChargeID { int ID_Track_Charge(const std::vector &Track); private: - TMS_Hit Hit; - std::vector Track; + // These are never used, todo delete? + //TMS_Hit Hit; + //std::vector Track; }; #endif From 8d8359fe0875f82a395f1d052f22440e84866f8b Mon Sep 17 00:00:00 2001 From: Asa Nehm Date: Fri, 2 Aug 2024 08:46:12 -0500 Subject: [PATCH 6/7] Adding the charge to TMS_TreeWriter for them to be saved to the output --- src/TMS_TreeWriter.cpp | 3 +++ src/TMS_TreeWriter.h | 1 + 2 files changed, 4 insertions(+) diff --git a/src/TMS_TreeWriter.cpp b/src/TMS_TreeWriter.cpp index 547e6804..3379e7c0 100644 --- a/src/TMS_TreeWriter.cpp +++ b/src/TMS_TreeWriter.cpp @@ -236,6 +236,7 @@ void TMS_TreeWriter::MakeBranches() { Reco_Tree->Branch("EnergyRange", RecoTrackEnergyRange, "EnergyRange[nTracks]/F"); Reco_Tree->Branch("EnergyDeposit",RecoTrackEnergyDeposit, "EnergyDeposit[nTracks]/F"); Reco_Tree->Branch("Length", RecoTrackLength, "Length[nTracks]/F"); + Reco_Tree->Branch("Charge", RecoTrackCharge, "Charge[nTracks]/I"); MakeTruthBranches(Truth_Info); @@ -1204,6 +1205,7 @@ void TMS_TreeWriter::Fill(TMS_Event &event) { RecoTrackEnergyRange[itTrack] = RecoTrack->EnergyRange; RecoTrackLength[itTrack] = 0.5 * (TrackLengthU[itTrack] + TrackLengthV[itTrack]); // RecoTrack->Length;, 2d is better estimate than 3d because of y jumps RecoTrackEnergyDeposit[itTrack] = RecoTrack->EnergyDeposit; + RecoTrackCharge[itTrack] = RecoTrack->Charge; for (int j = 0; j < 3; j++) { RecoTrackStartPos[itTrack][j] = RecoTrack->Start[j]; RecoTrackEndPos[itTrack][j] = RecoTrack->End[j]; @@ -1651,6 +1653,7 @@ void TMS_TreeWriter::Clear() { RecoTrackEnergyRange[i] = DEFAULT_CLEARING_FLOAT; RecoTrackEnergyDeposit[i] = DEFAULT_CLEARING_FLOAT; RecoTrackLength[i] = DEFAULT_CLEARING_FLOAT; + RecoTrackCharge[i] = DEFAULT_CLEARING_FLOAT; } RecoTrackN = 0; diff --git a/src/TMS_TreeWriter.h b/src/TMS_TreeWriter.h index b926339e..2a523cd8 100644 --- a/src/TMS_TreeWriter.h +++ b/src/TMS_TreeWriter.h @@ -52,6 +52,7 @@ class TMS_TreeWriter { float RecoTrackEnergyRange[__TMS_MAX_TRACKS__]; float RecoTrackEnergyDeposit[__TMS_MAX_TRACKS__]; float RecoTrackLength[__TMS_MAX_TRACKS__]; + int RecoTrackCharge[__TMS_MAX_TRACKS__]; private: TMS_TreeWriter(); From 4a5995a5c63d3d594a3efdd14d0bacbfe1713819 Mon Sep 17 00:00:00 2001 From: Asa Nehm Date: Thu, 22 Aug 2024 06:57:41 -0500 Subject: [PATCH 7/7] Added the performance evaluation of the charge ID --- scripts/Reco/performance_reco.py | 455 +++++++++++++++++++++++++++---- 1 file changed, 398 insertions(+), 57 deletions(-) diff --git a/scripts/Reco/performance_reco.py b/scripts/Reco/performance_reco.py index fb0a3503..37f7dd15 100755 --- a/scripts/Reco/performance_reco.py +++ b/scripts/Reco/performance_reco.py @@ -16,12 +16,23 @@ green_cbf = '#009e73' mp.style.use('seaborn-poster') -mp.rc('axes', labelsize = 12) # fontsize of the x and y labels -mp.rc('xtick', labelsize = 12) # fontsize of the tick labels -mp.rc('ytick', labelsize = 12) # fontsize of the tick labels +mp.rc('axes', labelsize = 15) # fontsize of the x and y labels +mp.rc('xtick', labelsize = 15) # fontsize of the tick labels +mp.rc('ytick', labelsize = 15) # fontsize of the tick labels + +# dyslexia friendly background +mp.rcParams['axes.facecolor'] = '#fafafa' +mp.rcParams["figure.facecolor"] = '#fafafa' +mp.rcParams["savefig.facecolor"] = '#fafafa' + +# off-black text +mp.rcParams['text.color'] = '#424242' +mp.rcParams['axes.labelcolor'] = '#424242' +mp.rcParams['xtick.color'] = '#424242' +mp.rcParams['ytick.color'] = '#424242' ### Actual function that loops through the spills -def draw_performance(out_dir, input_filename): +def draw_performance(out_dir, input_filename, Xlayers): if not os.path.exists(input_filename): raise ValueError(f"Cannor find input_filename {input_filename}") # Make sure we read in the correct file and have the output directory @@ -53,6 +64,9 @@ def draw_performance(out_dir, input_filename): Primary_True_End = np.ones((n_events, 5, 3), dtype = float) * -9999. True_TrackLength = np.ones((n_events, 5), dtype = float) * -9999. Reco_TrackLength = np.ones((n_events, 5), dtype = float) * -9999. + Reco_Charge = np.ones((n_events, 5), dtype = float) * -9999. + True_Charge = np.ones((n_events, 5), dtype = float) * -9999. + True_KE = np.ones((n_events, 5), dtype = float) * -9999. for current_spill_number in range(max_n_spills): for i in range(n_events): @@ -83,6 +97,9 @@ def draw_performance(out_dir, input_filename): RecoTrackPrimaryParticleTruePositionTrackEnd = np.frombuffer(true_event.RecoTrackPrimaryParticleTruePositionTrackEnd, dtype = np.float32) Muon_TrueTrackLength = true_event.Muon_TrueTrackLength Reco_Track_Length = np.frombuffer(event.Length, dtype = np.float32) + Reco_Track_Charge = event.Charge + MomentumTrackStart = np.frombuffer(true_event.RecoTrackPrimaryParticleTrueMomentumTrackStart, dtype = np.float32) + True_PDG = true_event.PDG nTracks = event.nTracks if nTracks <= 0: continue @@ -101,8 +118,14 @@ def draw_performance(out_dir, input_filename): Primary_True_Start[i, j, 0] = RecoTrackPrimaryParticleTruePositionTrackStart[j*4 + 0] Primary_True_Start[i, j, 1] = RecoTrackPrimaryParticleTruePositionTrackStart[j*4 + 1] Primary_True_Start[i, j, 2] = RecoTrackPrimaryParticleTruePositionTrackStart[j*4 + 2] + Primary_True_px = MomentumTrackStart[j*4 + 0] + Primary_True_py = MomentumTrackStart[j*4 + 1] + Primary_True_pz = MomentumTrackStart[j*4 + 2] + True_KE[i, j] = np.sqrt((Primary_True_px**2 + Primary_True_py**2 + Primary_True_pz**2 + 105.7**2) - 105.7) True_TrackLength[i, j] = Muon_TrueTrackLength Reco_TrackLength[i, j] = Reco_Track_Length[j] + True_Charge[i, j] = True_PDG[j] + Reco_Charge[i, j] = Reco_Track_Charge[j] if RecoTrackPrimaryParticleTruePositionTrackEnd[j*4 + 0] > -8000. and not EndPos.size == 0: if (EndPos[j*3 + 2] - StartPos[j*3 + 2]) > 890. and (RecoTrackPrimaryParticleTruePositionTrackEnd[j*4 + 2] - RecoTrackPrimaryParticleTruePositionTrackStart[j*4 + 2]) > 890.: Reco_End[i, j, 0] = EndPos[j*3 + 0] @@ -125,25 +148,30 @@ def draw_performance(out_dir, input_filename): True_TrackLength = True_TrackLength[boolean_True_TrackLength] boolean_Reco_TrackLength = (Reco_TrackLength != -9999.) Reco_TrackLength = Reco_TrackLength[boolean_Reco_TrackLength] + boolean_Reco_Charge = (Reco_Charge != -9999.) + Reco_Charge = Reco_Charge[boolean_Reco_Charge] + boolean_True_Charge = (True_Charge != -9999.) + True_Charge = True_Charge[boolean_True_Charge] + boolean_True_KE = (True_KE != -9999.) + True_KE = True_KE[boolean_True_KE] # flatten arrays - Reco_Start_x = Reco_Start[:, 0]#.flatten() - Reco_Start_y = Reco_Start[:, 1]#.flatten() - Reco_Start_z = Reco_Start[:, 2]#.flatten() - Reco_End_x = Reco_End[:, 0]#.flatten() - Reco_End_y = Reco_End[:, 1]#.flatten() - Reco_End_z = Reco_End[:, 2]#.flatten() - Primary_True_Start_x = Primary_True_Start[:, 0]#.flatten() - Primary_True_Start_y = Primary_True_Start[:, 1]#.flatten() - Primary_True_Start_z = Primary_True_Start[:, 2]#.flatten() - Primary_True_End_x = Primary_True_End[:, 0]#.flatten() - Primary_True_End_y = Primary_True_End[:, 1]#.flatten() - Primary_True_End_z = Primary_True_End[:, 2]#.flatten() - #True_TrackLength.flatten() - #Reco_TrackLength.flatten() + Reco_Start_x = Reco_Start[:, 0] + Reco_Start_y = Reco_Start[:, 1] + Reco_Start_z = Reco_Start[:, 2] + Reco_End_x = Reco_End[:, 0] + Reco_End_y = Reco_End[:, 1] + Reco_End_z = Reco_End[:, 2] + Primary_True_Start_x = Primary_True_Start[:, 0] + Primary_True_Start_y = Primary_True_Start[:, 1] + Primary_True_Start_z = Primary_True_Start[:, 2] + Primary_True_End_x = Primary_True_End[:, 0] + Primary_True_End_y = Primary_True_End[:, 1] + Primary_True_End_z = Primary_True_End[:, 2] # total number of events after filtering print("#events reconstruction: ", len(Reco_Start), "# events truth: ", len(Primary_True_Start)) + print("tracklength truth: ", len(True_TrackLength), "reco: ", len(Reco_TrackLength)) # subtract reconstruction from truth for all directions Diff_Start_x = Primary_True_Start_x - Reco_Start_x @@ -192,8 +220,8 @@ def draw_performance(out_dir, input_filename): mp.grid(True, linestyle = '--', alpha = 0.2) mp.xlabel('Difference truth – reconstruction [mm]') mp.ylabel('#') - mp.title('Start x', fontsize = 'x-large') - mp.legend(loc = 'best', fontsize = 'x-large', markerscale = 1.0, columnspacing = 0.5, handlelength = 0.8) + mp.title('Start x', fontsize = 'xx-large') + mp.legend(loc = 'best', fontsize = 'xx-large', markerscale = 1.0, columnspacing = 0.5, handlelength = 0.8) mp.savefig('%s/Difference_Start_X_hist.png' % out_dir, bbox_inches = 'tight') mp.close() @@ -204,8 +232,8 @@ def draw_performance(out_dir, input_filename): mp.grid(True, linestyle = '--', alpha = 0.2) mp.xlabel('Difference truth – reconstruction [mm]') mp.ylabel('#') - mp.title('Start y', fontsize = 'x-large') - mp.legend(loc = 'best', fontsize = 'x-large', markerscale = 1.0, columnspacing = 0.5, handlelength = 0.8) + mp.title('Start y', fontsize = 'xx-large') + mp.legend(loc = 'best', fontsize = 'xx-large', markerscale = 1.0, columnspacing = 0.5, handlelength = 0.8) mp.savefig('%s/Difference_Start_Y_hist.png' % out_dir, bbox_inches = 'tight') mp.close() @@ -216,8 +244,8 @@ def draw_performance(out_dir, input_filename): mp.grid(True, linestyle = '--', alpha = 0.2) mp.xlabel('Difference truth – reconstruction [mm]') mp.ylabel('#') - mp.title('Start z', fontsize = 'x-large') - mp.legend(loc = 'best', fontsize = 'x-large', markerscale = 1.0, columnspacing = 0.5, handlelength = 0.8) + mp.title('Start z', fontsize = 'xx-large') + mp.legend(loc = 'best', fontsize = 'xx-large', markerscale = 1.0, columnspacing = 0.5, handlelength = 0.8) mp.savefig('%s/Difference_Start_Z_hist.png' % out_dir, bbox_inches = 'tight') mp.close() @@ -228,8 +256,8 @@ def draw_performance(out_dir, input_filename): mp.grid(True, linestyle = '--', alpha = 0.2) mp.xlabel('Difference truth – reconstruction [mm]') mp.ylabel('#') - mp.title('End x', fontsize = 'x-large') - mp.legend(loc = 'best', fontsize = 'x-large', markerscale = 1.0, columnspacing = 0.5, handlelength = 0.8) + mp.title('End x', fontsize = 'xx-large') + mp.legend(loc = 'best', fontsize = 'xx-large', markerscale = 1.0, columnspacing = 0.5, handlelength = 0.8) mp.savefig('%s/Difference_End_X_hist.png' % out_dir, bbox_inches = 'tight') mp.close() @@ -240,8 +268,8 @@ def draw_performance(out_dir, input_filename): mp.grid(True, linestyle = '--', alpha = 0.2) mp.xlabel('Difference truth – reconstruction [mm]') mp.ylabel('#') - mp.title('End y', fontsize = 'x-large') - mp.legend(loc = 'best', fontsize = 'x-large', markerscale = 1.0, columnspacing = 0.5, handlelength = 0.8) + mp.title('End y', fontsize = 'xx-large') + mp.legend(loc = 'best', fontsize = 'xx-large', markerscale = 1.0, columnspacing = 0.5, handlelength = 0.8) mp.savefig('%s/Difference_End_Y_hist.png' % out_dir, bbox_inches = 'tight') mp.close() @@ -252,8 +280,8 @@ def draw_performance(out_dir, input_filename): mp.grid(True, linestyle = '--', alpha = 0.2) mp.xlabel('Difference truth – reconstruction [mm]') mp.ylabel('#') - mp.title('End z', fontsize = 'x-large') - mp.legend(loc = 'best', fontsize = 'x-large', markerscale = 1.0, columnspacing = 0.5, handlelength = 0.8) + mp.title('End z', fontsize = 'xx-large') + mp.legend(loc = 'best', fontsize = 'xx-large', markerscale = 1.0, columnspacing = 0.5, handlelength = 0.8) mp.savefig('%s/Difference_End_Z_hist.png' % out_dir, bbox_inches = 'tight') mp.close() @@ -279,6 +307,33 @@ def draw_performance(out_dir, input_filename): Primary_True_End_y_hist, Primary_True_End_y_bins = np.histogram(Primary_True_End_y, bins = Reco_Start_y_bins) Primary_True_End_z_hist, Primary_True_End_z_bins = np.histogram(Primary_True_End_z, bins = Reco_Start_z_bins) + # filter histograms for tracks with and without X hits for Xlayers == true + if (Xlayers): + start_helper = np.array([check_orientation(int(Reco_Start_z[i])) == 'kXBar' for i in range(len(Reco_Start_z))], dtype = bool) + Reco_Start_x_WX = Reco_Start_x[start_helper] + Reco_Start_y_WX = Reco_Start_y[start_helper] + Reco_Start_z_WX = Reco_Start_z[start_helper] + end_helper = np.array([check_orientation(int(Reco_End_z[i])) == 'kXBar' for i in range(len(Reco_End_z))], dtype = bool) + Reco_End_x_WX = Reco_End_x[end_helper] + Reco_End_y_WX = Reco_End_y[end_helper] + Reco_End_z_WX = Reco_End_z[end_helper] + + Reco_Start_x_WX_hist, Reco_Start_x_bins = np.histogram(Reco_Start_x_WX, bins = Reco_Start_x_bins) + Reco_Start_y_WX_hist, Reco_Start_y_bins = np.histogram(Reco_Start_y_WX, bins = Reco_Start_y_bins) + Reco_Start_z_WX_hist, Reco_Start_z_bins = np.histogram(Reco_Start_z_WX, bins = Reco_Start_z_bins) + + Reco_End_x_WX_hist, Reco_End_x_bins = np.histogram(Reco_End_x_WX, bins = Reco_End_x_bins) + Reco_End_y_WX_hist, Reco_End_y_bins = np.histogram(Reco_End_y_WX, bins = Reco_End_y_bins) + Reco_End_z_WX_hist, Reco_End_z_bins = np.histogram(Reco_End_z_WX, bins = Reco_End_z_bins) + + Reco_Start_x_WX_histX, Reco_Start_x_WX_histY = histogram_arr_handle(Reco_Start_x_WX_hist, Reco_Start_x_bins) + Reco_Start_y_WX_histX, Reco_Start_y_WX_histY = histogram_arr_handle(Reco_Start_y_WX_hist, Reco_Start_y_bins) + Reco_Start_z_WX_histX, Reco_Start_z_WX_histY = histogram_arr_handle(Reco_Start_z_WX_hist, Reco_Start_z_bins) + + Reco_End_x_WX_histX, Reco_End_x_WX_histY = histogram_arr_handle(Reco_End_x_WX_hist, Reco_End_x_bins) + Reco_End_y_WX_histX, Reco_End_y_WX_histY = histogram_arr_handle(Reco_End_y_WX_hist, Reco_End_y_bins) + Reco_End_z_WX_histX, Reco_End_z_WX_histY = histogram_arr_handle(Reco_End_z_WX_hist, Reco_End_z_bins) + # make the histograms usable Reco_Start_x_histX, Reco_Start_x_histY = histogram_arr_handle(Reco_Start_x_hist, Reco_Start_x_bins) Reco_Start_y_histX, Reco_Start_y_histY = histogram_arr_handle(Reco_Start_y_hist, Reco_Start_y_bins) @@ -296,79 +351,87 @@ def draw_performance(out_dir, input_filename): # now plot this mp.plot(Primary_True_Start_x_histX, Primary_True_Start_x_histY, color = blue_cbf, linestyle = '-.', linewidth = 2, label = 'truth') mp.plot(Reco_Start_x_histX, Reco_Start_x_histY, color = red_cbf, linestyle = ':', linewidth = 2, label = 'reco') + if (Xlayers): mp.plot(Reco_Start_x_WX_histX, Reco_Start_x_WX_histY, color = black_cbf, linestyle = '-', linewidth = 2, label = 'just X (reco)') mp.fill_between(Primary_True_Start_x_histX, 0, Primary_True_Start_x_histY, color = blue_cbf, alpha = 0.6, hatch = '//') mp.fill_between(Reco_Start_x_histX, 0, Reco_Start_x_histY, alpha = 0.6, hatch = '\\\\', color = red_cbf) mp.grid(True, linestyle = '--', alpha = 0.2) - mp.legend(loc = 'best', fontsize = 'x-large', markerscale = 1.0, columnspacing = 0.5, handlelength = 0.8) + mp.legend(loc = 'best', fontsize = 'xx-large', markerscale = 1.0, columnspacing = 0.5, handlelength = 0.8) mp.xlabel('Start in TMS [mm]') mp.ylabel('#') - mp.title('Start x', fontsize = 'x-large') + mp.title('Start x', fontsize = 'xx-large') mp.savefig('%s/Track_Start_X_hist.png' % out_dir, bbox_inches = 'tight') mp.close() mp.plot(Primary_True_Start_y_histX, Primary_True_Start_y_histY, color = blue_cbf, linestyle = '-.', linewidth = 2, label = 'truth') mp.plot(Reco_Start_y_histX, Reco_Start_y_histY, color = red_cbf, linestyle = ':', linewidth = 2, label = 'reco') + if (Xlayers): mp.plot(Reco_Start_y_WX_histX, Reco_Start_y_WX_histY, color = black_cbf, linestyle = '-', linewidth = 2, label = 'just X (reco)') mp.fill_between(Primary_True_Start_y_histX, 0, Primary_True_Start_y_histY, color = blue_cbf, alpha = 0.6, hatch = '//') mp.fill_between(Reco_Start_y_histX, 0, Reco_Start_y_histY, alpha = 0.6, hatch = '\\\\', color = red_cbf) mp.grid(True, linestyle = '--', alpha = 0.2) - mp.legend(loc = 'best', fontsize = 'x-large', markerscale = 1.0, columnspacing = 0.5, handlelength = 0.8) + mp.legend(loc = 'best', fontsize = 'xx-large', markerscale = 1.0, columnspacing = 0.5, handlelength = 0.8) mp.xlabel('Start in TMS [mm]') mp.ylabel('#') - mp.title('Start y', fontsize = 'x-large') + mp.title('Start y', fontsize = 'xx-large') mp.savefig('%s/Track_Start_Y_hist.png' % out_dir, bbox_inches = 'tight') mp.close() mp.plot(Primary_True_Start_z_histX, Primary_True_Start_z_histY, color = blue_cbf, linestyle = '-.', linewidth = 2, label = 'truth') mp.plot(Reco_Start_z_histX, Reco_Start_z_histY, color = red_cbf, linestyle = ':', linewidth = 2, label = 'reco') + if (Xlayers): mp.plot(Reco_Start_z_WX_histX, Reco_Start_z_WX_histY, color = black_cbf, linestyle = '-', linewidth = 2, label = 'just X (reco)') mp.fill_between(Primary_True_Start_z_histX, 0, Primary_True_Start_z_histY, color = blue_cbf, alpha = 0.6, hatch = '//') mp.fill_between(Reco_Start_z_histX, 0, Reco_Start_z_histY, alpha = 0.6, hatch = '\\\\', color = red_cbf) mp.grid(True, linestyle = '--', alpha = 0.2) - mp.legend(loc = 'best', fontsize = 'x-large', markerscale = 1.0, columnspacing = 0.5, handlelength = 0.8) + mp.legend(loc = 'best', fontsize = 'xx-large', markerscale = 1.0, columnspacing = 0.5, handlelength = 0.8) mp.xlabel('Start in TMS [mm]') mp.ylabel('#') - mp.title('Start z', fontsize = 'x-large') + mp.title('Start z', fontsize = 'xx-large') mp.savefig('%s/Track_Start_Z_hist.png' % out_dir, bbox_inches = 'tight') mp.close() mp.plot(Primary_True_End_x_histX, Primary_True_End_x_histY, color = blue_cbf, linestyle = '-.', linewidth = 2, label = 'truth') mp.plot(Reco_End_x_histX, Reco_End_x_histY, color = red_cbf, linestyle = ':', linewidth = 2, label = 'reco') + if (Xlayers): mp.plot(Reco_End_x_WX_histX, Reco_End_x_WX_histY, color = black_cbf, linestyle = '-', linewidth = 2, label = 'just X (reco)') mp.fill_between(Primary_True_End_x_histX, 0, Primary_True_End_x_histY, color = blue_cbf, alpha = 0.6, hatch = '//') mp.fill_between(Reco_End_x_histX, 0, Reco_End_x_histY, alpha = 0.6, hatch = '\\\\', color = red_cbf) + if (Xlayers): mp.fill_between(Reco_End_x_WX_histX, 0, Reco_End_x_WX_histY, color = black_cbf, alpha = 0.2) mp.grid(True, linestyle = '--', alpha = 0.2) - mp.legend(loc = 'best', fontsize = 'x-large', markerscale = 1.0, columnspacing = 0.5, handlelength = 0.8) + mp.legend(loc = 'best', fontsize = 'xx-large', markerscale = 1.0, columnspacing = 0.5, handlelength = 0.8) mp.xlabel('End in TMS [mm]') mp.ylabel('#') - mp.title('End x', fontsize = 'x-large') + mp.title('End x', fontsize = 'xx-large') mp.savefig('%s/Track_End_X_hist.png' % out_dir, bbox_inches = 'tight') mp.close() mp.plot(Primary_True_End_y_histX, Primary_True_End_y_histY, color = blue_cbf, linestyle = '-.', linewidth = 2, label = 'truth') mp.plot(Reco_End_y_histX, Reco_End_y_histY, color = red_cbf, linestyle = ':', linewidth = 2, label = 'reco') + if (Xlayers): mp.plot(Reco_End_y_WX_histX, Reco_End_y_WX_histY, color = black_cbf, linestyle = '-', linewidth = 2, label = 'just X (reco)') mp.fill_between(Primary_True_End_y_histX, 0, Primary_True_End_y_histY, color = blue_cbf, alpha = 0.6, hatch = '//') mp.fill_between(Reco_End_y_histX, 0, Reco_End_y_histY, alpha = 0.6, hatch = '\\\\', color = red_cbf) + if (Xlayers): mp.fill_between(Reco_End_y_WX_histX, 0, Reco_End_y_WX_histY, color = black_cbf, alpha = 0.2) mp.grid(True, linestyle = '--', alpha = 0.2) - mp.legend(loc = 'best', fontsize = 'x-large', markerscale = 1.0, columnspacing = 0.5, handlelength = 0.8) + mp.legend(loc = 'best', fontsize = 'xx-large', markerscale = 1.0, columnspacing = 0.5, handlelength = 0.8) mp.xlabel('End in TMS [mm]') mp.ylabel('#') - mp.title('End y', fontsize = 'x-large') + mp.title('End y', fontsize = 'xx-large') mp.savefig('%s/Track_End_Y_hist.png' % out_dir, bbox_inches = 'tight') mp.close() mp.plot(Primary_True_End_z_histX, Primary_True_End_z_histY, color = blue_cbf, linestyle = '-.', linewidth = 2, label = 'truth') mp.plot(Reco_End_z_histX, Reco_End_z_histY, color = red_cbf, linestyle = ':', linewidth = 2, label = 'reco') + if (Xlayers): mp.plot(Reco_End_z_WX_histX, Reco_End_z_WX_histY, color = black_cbf, linestyle = '-', linewidth = 2, label = 'just X (reco)') mp.fill_between(Primary_True_End_z_histX, 0, Primary_True_End_z_histY, color = blue_cbf, alpha = 0.6, hatch = '//') mp.fill_between(Reco_End_z_histX, 0, Reco_End_z_histY, alpha = 0.6, hatch = '\\\\', color = red_cbf) + if (Xlayers): mp.fill_between(Reco_End_z_WX_histX, 0, Reco_End_z_WX_histY, color = black_cbf, alpha = 0.2) mp.grid(True, linestyle = '--', alpha = 0.2) - mp.legend(loc = 'best', fontsize = 'x-large', markerscale = 1.0, columnspacing = 0.5, handlelength = 0.8) + mp.legend(loc = 'best', fontsize = 'xx-large', markerscale = 1.0, columnspacing = 0.5, handlelength = 0.8) mp.xlabel('End in TMS [mm]') mp.ylabel('#') - mp.title('End z', fontsize = 'x-large') + mp.title('End z', fontsize = 'xx-large') mp.savefig('%s/Track_End_Z_hist.png' % out_dir, bbox_inches = 'tight') mp.close() ### plot difference in dependence of hit position # create 2d histograms - #print(min(min(Primary_True_Start[:, 1]), min(Primary_True_End[:, 1])), max(max(Primary_True_Start[:, 1]), max(Primary_True_End[:, 1])), min(min(Reco_Start[:, 1]), min(Reco_End[:, 1])), max(max(Reco_Start[:, 1]), max(Reco_End[:, 1]))) dependence_Start_x_hist, dependence_Start_x_binsX, dependence_Start_x_binsY = np.histogram2d(Primary_True_Start_x, Diff_Start_x, bins = [Primary_True_Start_x_bins, Diff_Start_x_bins]) dependence_Start_y_hist, dependence_Start_y_binsX, dependence_Start_y_binsY = np.histogram2d(Primary_True_Start_y, Diff_Start_y, bins = [Primary_True_Start_y_bins, Diff_Start_y_bins]) @@ -382,7 +445,7 @@ def draw_performance(out_dir, input_filename): im = mp.pcolormesh(dependence_Start_x_binsX, dependence_Start_x_binsY, np.transpose(dependence_Start_x_hist), cmap = cmap); mp.xlabel('Start in TMS (X) [mm]') mp.ylabel('Difference truth - reco [mm]') - mp.title('Start x', fontsize = 'x-large') + mp.title('Start x', fontsize = 'xx-large') mp.colorbar(im); mp.savefig('%s/Start_X_real_and_difference.png' % out_dir, bbox_inches = 'tight'); mp.close(); @@ -390,7 +453,7 @@ def draw_performance(out_dir, input_filename): im = mp.pcolormesh(dependence_Start_y_binsX, dependence_Start_y_binsY, np.transpose(dependence_Start_y_hist), cmap = cmap); mp.xlabel('Start in TMS (Y) [mm]') mp.ylabel('Difference truth - reco [mm]') - mp.title('Start y', fontsize = 'x-large') + mp.title('Start y', fontsize = 'xx-large') mp.colorbar(im); mp.savefig('%s/Start_Y_real_and_difference.png' % out_dir, bbox_inches = 'tight'); mp.close(); @@ -398,7 +461,7 @@ def draw_performance(out_dir, input_filename): im = mp.pcolormesh(dependence_Start_z_binsX, dependence_Start_z_binsY, np.transpose(dependence_Start_z_hist), cmap = cmap); mp.xlabel('Start in TMS (Z) [mm]') mp.ylabel('Difference truth - reco [mm]') - mp.title('Start z', fontsize = 'x-large') + mp.title('Start z', fontsize = 'xx-large') mp.colorbar(im); mp.savefig('%s/Start_Z_real_and_difference.png' % out_dir, bbox_inches = 'tight'); mp.close(); @@ -406,7 +469,7 @@ def draw_performance(out_dir, input_filename): im = mp.pcolormesh(dependence_End_x_binsX, dependence_End_x_binsY, np.transpose(dependence_End_x_hist), cmap = cmap); mp.xlabel('End in TMS (X) [mm]') mp.ylabel('Difference truth - reco [mm]') - mp.title('End x', fontsize = 'x-large') + mp.title('End x', fontsize = 'xx-large') mp.colorbar(im); mp.savefig('%s/End_X_real_and_difference.png' % out_dir, bbox_inches = 'tight'); mp.close(); @@ -414,7 +477,7 @@ def draw_performance(out_dir, input_filename): im = mp.pcolormesh(dependence_End_y_binsX, dependence_End_y_binsY, np.transpose(dependence_End_y_hist), cmap = cmap); mp.xlabel('End in TMS (Y) [mm]') mp.ylabel('Difference truth - reco [mm]') - mp.title('End y', fontsize = 'x-large') + mp.title('End y', fontsize = 'xx-large') mp.colorbar(im); mp.savefig('%s/End_Y_real_and_difference.png' % out_dir, bbox_inches = 'tight'); mp.close(); @@ -422,7 +485,7 @@ def draw_performance(out_dir, input_filename): im = mp.pcolormesh(dependence_End_z_binsX, dependence_End_z_binsY, np.transpose(dependence_End_z_hist), cmap = cmap); mp.xlabel('End in TMS (Z) [mm]') mp.ylabel('Difference truth - reco [mm]') - mp.title('End z', fontsize = 'x-large') + mp.title('End z', fontsize = 'xx-large') mp.colorbar(im); mp.savefig('%s/End_Z_real_and_difference.png' % out_dir, bbox_inches = 'tight'); mp.close(); @@ -431,8 +494,11 @@ def draw_performance(out_dir, input_filename): # histogram # single histogram as well - reco_tracklength_hist, reco_tracklength_bins = np.histogram(Reco_TrackLength / 10, bins = 40, range = (50, 500))#(min(min(Reco_TrackLength / 10), min(True_TrackLength)), max(max(Reco_TrackLength / 10), max(True_TrackLength)))) - true_tracklength_hist, true_tracklength_bins = np.histogram(True_TrackLength, bins = reco_tracklength_bins) + reco_tracklength_hist, reco_tracklength_bins = np.histogram(True_TrackLength - Reco_TrackLength / 10, bins = 40)#, range = (min(min(True_TrackLength - Reco_TrackLength / 10), min(True_TrackLength)), max(max(True_TrackLength - Reco_TrackLength / 10), max(True_TrackLength))))#(50, 5000))# + true_tracklength_hist, true_tracklength_bins = np.histogram(True_TrackLength, bins = 40)#reco_tracklength_bins) + + #true_tracklength_0 = True_TrackLength[True_TrackLength == 0] + #print(len(true_tracklength_0)) reco_tracklength_histX, reco_tracklength_histY = histogram_arr_handle(reco_tracklength_hist, reco_tracklength_bins) true_tracklength_histX, true_tracklength_histY = histogram_arr_handle(true_tracklength_hist, true_tracklength_bins) @@ -441,22 +507,99 @@ def draw_performance(out_dir, input_filename): mp.plot(reco_tracklength_histX, reco_tracklength_histY, color = orange_cbf, linestyle = ':', linewidth = 2, label = 'reco') mp.fill_between(true_tracklength_histX, 0, true_tracklength_histY, alpha = 0.6, color = blue_cbf, hatch = '\\\\') mp.fill_between(reco_tracklength_histX, 0, reco_tracklength_histY, alpha = 0.6, color = orange_cbf, hatch = '//') - mp.xlabel('Reco track length [$\\frac{g}{cm^2}$]') + mp.xlabel('True track length - Reco track length [$\\frac{g}{cm^2}$]') mp.ylabel('#') - mp.legend(loc = 'best', fontsize = 'x-large', markerscale = 1.0, columnspacing = 0.5, handlelength = 0.8) + mp.legend(loc = 'best', fontsize = 'xx-large', markerscale = 1.0, columnspacing = 0.5, handlelength = 0.8) mp.grid(True, linestyle = '--', alpha = 0.2) mp.savefig('%s/Reco_tracklength.png' % out_dir, bbox_inches = 'tight') mp.close() - tracklength_hist, tracklength_binsX, tracklength_binsY = np.histogram2d(True_TrackLength, Reco_TrackLength / 10, bins = [true_tracklength_bins, reco_tracklength_bins]) + tracklength_hist, tracklength_binsX, tracklength_binsY = np.histogram2d(True_TrackLength, True_TrackLength - Reco_TrackLength / 10, bins = [true_tracklength_bins, reco_tracklength_bins]) + + cmap = cm.get_cmap('cividis'); im = mp.pcolormesh(tracklength_binsX, tracklength_binsY, np.transpose(tracklength_hist), cmap = cmap) mp.xlabel('True track length [$\\frac{g}{cm^2}$]') - mp.ylabel('Reco track length [$\\frac{g}{cm^2}$]') + mp.ylabel('True track length - Reco track length [$\\frac{g}{cm^2}$]') mp.colorbar(im) mp.savefig('%s/Tracklength.png' % out_dir, bbox_inches = 'tight') mp.close() + ### single slices for comparison of geometries + # low energy slice + # find all events with energy between lower and upper bound + combined_TrackLength = True_TrackLength - Reco_TrackLength / 10 + low_helper = (True_TrackLength > 150.) & (True_TrackLength < 250.) + high_helper = (True_TrackLength > 400.) + low_energy_slice = combined_TrackLength[low_helper] + high_energy_slice = combined_TrackLength[high_helper] + + # calculate mean and std of this slice + low_energy_mean = np.mean(low_energy_slice) + low_energy_std = np.std(low_energy_slice) + print('Low energy: ', low_energy_mean, ' ± ', low_energy_std) + + high_energy_mean = np.mean(high_energy_slice) + high_energy_std = np.std(high_energy_slice) + print('High energy: ', high_energy_mean, ' ± ', high_energy_std) + + # plot slice distribution + low_energy_slice_hist, low_energy_slice_bins = np.histogram(low_energy_slice, bins = reco_tracklength_bins) + high_energy_slice_hist, high_energy_slice_bins = np.histogram(high_energy_slice, bins = reco_tracklength_bins) + + low_energy_slice_histX, low_energy_slice_histY = histogram_arr_handle(low_energy_slice_hist, low_energy_slice_bins) + high_energy_slice_histX, high_energy_slice_histY = histogram_arr_handle(high_energy_slice_hist, high_energy_slice_bins) + + + mp.plot(low_energy_slice_histX, low_energy_slice_histY, color = blue_cbf, linestyle = '--', linewidth = 2) + mp.plot(high_energy_slice_histX, high_energy_slice_histY, color = orange_cbf, linestyle = ':', linewidth = 2) + mp.fill_between(low_energy_slice_histX, 0, low_energy_slice_histY, alpha = 0.6, color = blue_cbf, hatch = '//', label = 'low energy [150 - 250]') + mp.fill_between(high_energy_slice_histX, 0, high_energy_slice_histY, alpha = 0.6, color = orange_cbf, hatch = '\\\\', label = 'high energy [> 400]') + mp.xlabel('True track length - Reco track length [$\\frac{g}{cm^2}$]') + mp.ylabel('#') + mp.legend(loc = 'upper left', fontsize = 'xx-large', markerscale = 1.0, columnspacing = 0.5, handlelength = 0.8) + mp.grid(True, linestyle = '--', alpha = 0.2) + mp.savefig('%s/Low_energy_slice.png' % out_dir, bbox_inches = 'tight') + mp.close() + + + #reco_pure_length = np.sqrt((Reco_End[:, 0] - Reco_Start[:, 0])**2 + (Reco_End[:, 1] - Reco_Start[:, 1])**2 + (Reco_End[:, 2] - Reco_Start[:, 2])**2) + #true_pure_length = np.sqrt((Primary_True_End[:, 0] - Primary_True_Start[:, 0])**2 + (Primary_True_End[:, 1] - Primary_True_Start[:, 1])**2 + (Primary_True_End[:, 2] - Primary_True_Start[:, 2])**2) + + #reco_pure_length_hist, reco_pure_length_bins = np.histogram(reco_pure_length, bins = 40, range = (min(min(reco_pure_length), min(true_pure_length)), max(max(reco_pure_length), max(true_pure_length)))) + #true_pure_length_hist, true_pure_length_bins = np.histogram(true_pure_length, bins = reco_pure_length_bins) + + #reco_pure_length_histX, reco_pure_length_histY = histogram_arr_handle(reco_pure_length_hist, reco_pure_length_bins) + #true_pure_length_histX, true_pure_length_histY = histogram_arr_handle(true_pure_length_hist, true_pure_length_bins) + + #mp.plot(true_pure_length_histX, true_pure_length_histY, color = blue_cbf, linestyle = '--', linewidth = 2, label = 'truth') + #mp.plot(reco_pure_length_histX, reco_pure_length_histY, color = orange_cbf, linestyle = ':', linewidth = 2, label = 'reco') + #mp.fill_between(true_pure_length_histX, 0, true_pure_length_histY, alpha = 0.6, color = blue_cbf, hatch = '\\\\') + #mp.fill_between(reco_pure_length_histX, 0, reco_pure_length_histY, alpha = 0.6, color = orange_cbf, hatch = '//') + #mp.xlabel('Track length [mm]') + #mp.ylabel('#') + #mp.legend(loc = 'best', fontsize = 'xx-large', markerscale = 1.0, columnspacing = 0.5, handlelength = 0.8) + #mp.grid(True, linestyle = '--', alpha = 0.2) + #mp.savefig('%s/Reco_pure_length.png' % out_dir, bbox_inches = 'tight') + #mp.close(); + + """length_density_hist, length_density_binsX, length_density_binsY = np.histogram2d(reco_pure_length, Reco_TrackLength / 10, bins = [30, 30], range = [[min(reco_pure_length), 6500], [min(Reco_TrackLength / 10), max(Reco_TrackLength / 10)]]) + + im = mp.pcolormesh(length_density_binsX, length_density_binsY, np.transpose(length_density_hist), cmap = cmap) + mp.xlabel('Reco track length [mm]') + mp.ylabel('Reco track density length [$\\frac{g}{cm^2}$]') + mp.colorbar(im) + mp.savefig('%s/TrackLength_Density.png' % out_dir, bbox_inches = 'tight') + mp.close() + + length_density_true_hist, length_density_true_binsX, length_density_true_binsY = np.histogram2d(true_pure_length, True_TrackLength, bins = [30, 30], range = [[min(reco_pure_length), 6500], [min(Reco_TrackLength / 10), max(Reco_TrackLength / 10)]]) + + im = mp.pcolormesh(length_density_true_binsX, length_density_true_binsY, np.transpose(length_density_true_hist), cmap = cmap) + mp.xlabel('True track length [mm]') + mp.ylabel('True track density length [$\\frac{g}{cm^2}$]') + mp.colorbar(im) + mp.savefig('%s/TrackLength_Density_truth.png' % out_dir, bbox_inches = 'tight') + mp.close()""" ### Stopping vs exiting evaluation opposite_direction_counter = 0 @@ -478,7 +621,7 @@ def draw_performance(out_dir, input_filename): reco_stopping = False # check if line crosses detector volume in yz - if (true_slope * 18294 + true_intercept) < -2949 and Primary_True_End[i, 1] < (-2949 + 338.3): + if (true_slope * 18294 + true_intercept) < -2949 and Primary_True_End[i, 1] <= (-2939):# + 338.3): truth_exiting = True elif Primary_True_End[i, 2] < 17750: truth_stopping = True @@ -564,6 +707,147 @@ def draw_performance(out_dir, input_filename): mp.savefig('%s/Exiting_venn.png' % out_dir, bbox_inches = 'tight') mp.close() + + ### Charge ID + # Filter out all events without truth +/-13 as PDG + boolean_true_muon = (np.abs(True_Charge) == 13) + True_Charge = True_Charge[boolean_true_muon] + Reco_Charge = Reco_Charge[boolean_true_muon] + True_KE = True_KE[boolean_true_muon] + boolean_issues = (Reco_Charge == -9.99999999e+08) + not_identified = Reco_Charge[boolean_issues] + #True_KE = True_KE[Reco_Charge != -9.99999999e+08] + + # Counter for both positive, both negative, different reco/truth positive and negative + true_muons = len(True_Charge[True_Charge == 13]) + true_antimuons = len(True_Charge[True_Charge == -13]) + reco_muons = len(Reco_Charge[Reco_Charge == 13]) + reco_antimuons = len(Reco_Charge[Reco_Charge == -13]) + counter_true_positive = 0 #truth and reco agree on muon + counter_true_negative = 0 #truth and reco agree on antimuon + counter_false_positive = 0 #truth and reco disagree, reco -> muon + counter_false_negative = 0 #truth and reco disagree, reco -> antimuon + + # Prepare the KE arrays + Muons_KE = np.ones(len(True_KE), dtype = float) * -9999. + AMuons_KE = np.ones(len(True_KE), dtype = float) * -9999. + True_Muons_KE = True_KE[True_Charge == 13] + True_Antimuons_KE = True_KE[True_Charge == -13] + + # Compare sign of truth and reco + for i in range(len(True_Charge)): + if True_Charge[i] == Reco_Charge[i]: + if True_Charge[i] > 0: + counter_true_positive += 1 + Muons_KE[i] = True_KE[i] + else: + counter_true_negative += 1 + AMuons_KE[i] = True_KE[i] + else: + if Reco_Charge[i] > 0: + counter_false_positive += 1 + else: + counter_false_negative += 1 + + if (counter_true_positive + counter_true_negative + counter_false_positive + counter_false_negative) != len(True_Charge): + print("Counters don't add up") + print("Length: ", len(True_Charge)) + + print("Charge ID numbers") + print(" Not identified: ", len(not_identified)) + print(" True muons: ", counter_true_positive) + print(" True antimuons: ", counter_true_negative) + print(" False muons: ", counter_false_positive) + print(" False antimuons: ", counter_false_negative, " here the not identified go into!") + + #TODO look into 'edge' cases with region jumps + + # Performance evaluation + efficiency_muons = counter_true_positive / (counter_true_positive + (counter_false_negative - len(not_identified))) + efficiency_antimuons = counter_true_negative / (counter_true_negative + counter_false_positive) + purity_muons = counter_true_positive / (counter_true_positive + counter_false_positive) + purity_antimuons = counter_true_negative / (counter_true_negative + (counter_false_negative - len(not_identified))) + accuracy_anti_muons = (counter_true_positive + counter_true_negative) / (counter_true_positive + counter_true_negative + counter_false_positive + (counter_false_negative - len(not_identified))) + + print(" Muons (efficiency | purity): ", efficiency_muons, " | ", purity_muons) + print(" Antimuons (efficiency | purity): ", efficiency_antimuons, " | ", purity_antimuons) + print(" Accuracy (both): ", accuracy_anti_muons) + + # Plotting + # Total + v = venn2(subsets = (true_muons - counter_true_positive, reco_muons- counter_true_positive, counter_true_positive), set_labels = ('Truth', 'Reco'), set_colors = ('#03558e', '#e95125', '#f39e92'), alpha = 1) + c = venn2_circles(subsets = (true_muons - counter_true_positive, reco_muons - counter_true_positive, counter_true_positive), color = 'grey', linestyle = 'dashed', linewidth = 1.5) + v.get_label_by_id('10').set_text('%i' % true_muons) + v.get_label_by_id('10').set_color('lightgrey') + v.get_label_by_id('10').set_fontsize('x-large') + v.get_label_by_id('A').set_color('#03558e') + v.get_label_by_id('A').set_fontsize('xx-large') + v.get_patch_by_id('10').set_color('#03558e') + v.get_label_by_id('01').set_text('%i' % reco_muons) + v.get_label_by_id('01').set_color('black') + v.get_label_by_id('01').set_fontsize('x-large') + v.get_label_by_id('B').set_color('#e95125') + v.get_label_by_id('B').set_fontsize('xx-large') + v.get_patch_by_id('01').set_color('#e95125') + if counter_true_positive > 0: + v.get_label_by_id('11').set_text('%i' % counter_true_positive) + v.get_label_by_id('11').set_color('black') + v.get_label_by_id('11').set_fontsize('x-large') + v.get_patch_by_id('11').set_color('#f39e92') + mp.title('Charge ID $\\mu$') + mp.savefig('%s/Muons_venn.png' % out_dir, bbox_inches = 'tight') + mp.close() + + v = venn2(subsets = (true_antimuons - counter_true_negative, reco_antimuons - counter_true_negative, counter_true_negative), set_labels = ('Truth', 'Reco'), set_colors = ('#03558e', '#e95125', '#f39e92'), alpha = 1) + c = venn2_circles(subsets = (true_antimuons - counter_true_negative, reco_antimuons - counter_true_negative, counter_true_negative), color = 'grey', linestyle = 'dashed', linewidth = 1.5) + v.get_label_by_id('10').set_text('%i' % true_antimuons) + v.get_label_by_id('10').set_color('lightgrey') + v.get_label_by_id('10').set_fontsize('x-large') + v.get_label_by_id('A').set_color('#03558e') + v.get_label_by_id('A').set_fontsize('xx-large') + v.get_patch_by_id('10').set_color('#03558e') + v.get_label_by_id('01').set_text('%i' % reco_antimuons) + v.get_label_by_id('01').set_color('black') + v.get_label_by_id('01').set_fontsize('x-large') + v.get_label_by_id('B').set_color('#e95125') + v.get_label_by_id('B').set_fontsize('xx-large') + v.get_patch_by_id('01').set_color('#e95125') + if counter_true_negative > 0: + v.get_label_by_id('11').set_text('%i' % counter_true_negative) + v.get_label_by_id('11').set_color('black') + v.get_label_by_id('11').set_fontsize('x-large') + v.get_patch_by_id('11').set_color('#f39e92') + mp.title('Charge ID $\\bar{\\mu}$') + mp.savefig('%s/Antimuons_venn.png' % out_dir, bbox_inches = 'tight') + mp.close() + + # Energy dependent + Muons_KE = Muons_KE[Muons_KE != -9999.] + AMuons_KE = AMuons_KE[AMuons_KE != -9999.] + + muons_ke_hist, muons_ke_bins = np.histogram(Muons_KE, bins = 100, range = (0, 5000)) + amuons_ke_hist, amuons_ke_bins = np.histogram(AMuons_KE, bins = muons_ke_bins) + + true_muons_ke_hist, muons_ke_bins = np.histogram(True_Muons_KE, bins = muons_ke_bins) + true_antimuons_ke_hist, amuons_ke_bins = np.histogram(True_Antimuons_KE, bins = amuons_ke_bins) + + muons_ke_hist_x, muons_ke_hist_y = histogram_arr_handle(muons_ke_hist, muons_ke_bins) + amuons_ke_hist_x, amuons_ke_hist_y = histogram_arr_handle(amuons_ke_hist, amuons_ke_bins) + + true_muons_ke_hist_x, true_muons_ke_hist_y = histogram_arr_handle(true_muons_ke_hist, muons_ke_bins) + true_antimuons_ke_hist_x, true_antimuons_ke_hist_y = histogram_arr_handle(true_antimuons_ke_hist, amuons_ke_bins) + + mp.plot(muons_ke_hist_x, muons_ke_hist_y / true_muons_ke_hist_y, color = blue_cbf, linestyle = '--', linewidth = 2) + mp.plot(amuons_ke_hist_x, amuons_ke_hist_y / true_antimuons_ke_hist_y, color = orange_cbf, linestyle = '--', linewidth = 2) + mp.fill_between(muons_ke_hist_x, 0, muons_ke_hist_y / true_muons_ke_hist_y, color = blue_cbf, alpha = 0.6, hatch = '//', label = '$\\mu$ correct') + mp.fill_between(amuons_ke_hist_x, 0, amuons_ke_hist_y / true_antimuons_ke_hist_y, color = orange_cbf, alpha = 0.6, hatch = '\\\\', label = '$\\bar{\\mu}$ correct') + mp.xlabel('True Muon KE [MeV]') + mp.ylabel('Fraction') + mp.legend(loc = 'lower center', fontsize = 'xx-large', markerscale = 1.0, columnspacing = 0.5, handlelength = 0.8) + mp.grid(True, linestyle = '--', alpha = 0.2) + mp.savefig('%s/Charge_ID_KE.png' % out_dir, bbox_inches = 'tight') + mp.close() + return def eff_pur_acc(overlap, truth, reco): @@ -611,15 +895,72 @@ def histogram_arr_handle(bins, edges):#(u bins,v edges) return x_output, y_output +def check_orientation(hit_z): + return layer_dict["%s" % hit_z] + +first_z = 11368 +layer_dict = { "%s" % first_z : "kUBar" } + +def calculate_layers(Xlayers): + increment = 2 + if Xlayers: + increment = 3 + thin_layers = 39 + thick_layers = 61 + # Calculate the z position for each layer for the thin section + for i in range(thin_layers): + hit_z = first_z + i * 55 + # even layers + if not Xlayers: + if (((hit_z - first_z) / 55) % increment) == 0: + layer_dict.update({ "%s" % hit_z : "kUBar" }) + # odd layers + elif (((hit_z - first_z) / 55) % increment) == 1: + layer_dict.update({ "%s" % hit_z : "kVBar" }) + # x layers + if Xlayers: + # even layers + if (((hit_z - first_z) / 55) % 2) == 0: + layer_dict.update({ "%s" % hit_z : "kUBar" }) + # odd layers + elif (((hit_z - first_z) / 55) % 2) == 1: + layer_dict.update({ "%s" % hit_z : "kVBar" }) + if (((hit_z - first_z) / 55) % increment) == 0: + layer_dict.update({ "%s" % hit_z : "kXBar" }) + # Calculate the z position for each layers the thick section + start_thick = first_z + thin_layers * 55 + for i in range(thick_layers): + hit_z = start_thick + i * 80 + if not Xlayers: + # even layers + if (((hit_z - start_thick) / 80) % increment) == 0: + layer_dict.update({ "%s" % hit_z : "kVBar" }) + # odd layers + elif (((hit_z - start_thick) / 80) % increment) == 1: + layer_dict.update({ "%s" % hit_z : "kUBar" }) + # x layers + if Xlayers: + # even layers + if (((hit_z - start_thick) / 80) % 2) == 0: + layer_dict.update({ "%s" % hit_z : "kVBar" }) + # odd layers + elif (((hit_z - start_thick) / 80) % 2) == 1: + layer_dict.update({ "%s" % hit_z : "kUBar" }) + if (((hit_z - start_thick) / 80) % increment) == 0: + layer_dict.update({ "%s" % hit_z : "kXBar" }) if __name__ == "__main__": parser = argparse.ArgumentParser(description = "Draws spills.") parser.add_argument('--outdir', "-o", type = str, help = "The output dir. Will be made if it doesn't exist. Default = spills/", default = "spills") parser.add_argument('--input_filename', "-f", type = str, help = "The file with the events to draw.") + parser.add_argument('--Xlayers', "-X", help = "Does the geometry use X (90 degree orientated) scintillator layers? Yes -> --Xlayers, No -> --no-Xlayers", action = argparse.BooleanOptionalAction) args = parser.parse_args() + Xlayers = args.Xlayers + calculate_layers(Xlayers) + #print(layer_dict) out_dir = args.outdir input_filename = args.input_filename - draw_performance(out_dir, input_filename) + draw_performance(out_dir, input_filename, Xlayers)