From 70066380c4b18c05fc35bd54042446b11dc1788d Mon Sep 17 00:00:00 2001 From: Jeffrey Kleykamp Date: Fri, 20 Sep 2024 00:10:26 -0500 Subject: [PATCH 1/9] Add 'make sanitize' which is super useful. Super. Useful. --- Makefile | 4 ++++ app/Makefile | 4 ++++ src/Makefile | 4 ++++ 3 files changed, 12 insertions(+) diff --git a/Makefile b/Makefile index 3434c42d..f5fa954b 100644 --- a/Makefile +++ b/Makefile @@ -18,3 +18,7 @@ clean: all: check-submodule make -j8 -C src make -j8 -C app + +sanitize: check-submodule + make -j8 -C src sanitize + make -j8 -C app sanitize diff --git a/app/Makefile b/app/Makefile index c5e47966..b786e785 100755 --- a/app/Makefile +++ b/app/Makefile @@ -59,6 +59,10 @@ MKDIR_P := mkdir -p all: directories $(EXE_DIR)/ConvertToTMSTree.exe $(EXE_DIR)/BetheBloch_Example.exe $(EXE_DIR)/DBSCAN_test.exe $(EXE_DIR)/DrawEvents.exe $(EXE_DIR)/TrackLengthTester.exe $(EXE_DIR)/ShootRay.exe $(EXE_DIR)/BField_tester.exe $(EXE_DIR)/CherryPickEvents.exe +# Sanitize target +sanitize: CXXFLAGS += -fsanitize=address +sanitize: all + directories: $(EXE_DIR) $(EXE_DIR): diff --git a/src/Makefile b/src/Makefile index b7e44e33..3675205a 100644 --- a/src/Makefile +++ b/src/Makefile @@ -73,6 +73,10 @@ endif all: directories libTMS_Prod +# Sanitize target +sanitize: CXXFLAGS := $(CXXFLAGS) -fsanitize=address +sanitize: all + directories: $(LIB_DIR) $(LIB_DIR): From e8223b48cae9553e8cac7c0a4662eb9a586a3d4a Mon Sep 17 00:00:00 2001 From: Jeffrey Kleykamp Date: Mon, 23 Sep 2024 13:29:26 -0500 Subject: [PATCH 2/9] Make sure accumulator is within bounds --- src/TMS_Reco.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/TMS_Reco.cpp b/src/TMS_Reco.cpp index d7fc8e74..9106f3aa 100644 --- a/src/TMS_Reco.cpp +++ b/src/TMS_Reco.cpp @@ -3794,8 +3794,9 @@ void TMS_TrackFinder::Accumulate(double xhit, double zhit) { std::cout << "cbin: " << c_bin << std::endl; } */ - // Fill the accumulator - Accumulator[i][c_bin]++; + // Fill the accumulator, but only within bounds + // We don't care about lines outside of bounds + if (i > 0 && c_bin > 0 && i < nSlope && c_bin < nIntercept) Accumulator[i][c_bin]++; } } From 2e4f48d52bbf42100bce30efb79330c07502ec28 Mon Sep 17 00:00:00 2001 From: Jeffrey Kleykamp Date: Mon, 23 Sep 2024 13:42:25 -0500 Subject: [PATCH 3/9] Make sure pnt+1 is within bounds --- src/TMS_Event.cpp | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/src/TMS_Event.cpp b/src/TMS_Event.cpp index 1060130f..409d40b7 100644 --- a/src/TMS_Event.cpp +++ b/src/TMS_Event.cpp @@ -16,7 +16,6 @@ TMS_Event::TMS_Event() { LightWeight = true; } - static bool TMS_TrueParticle_NotWorthSaving(TMS_TrueParticle tp) { if (tp.GetTrueVisibleEnergy() == 0 && !tp.IsPrimary()) return true; // Don't worry about really low visible energy @@ -1016,16 +1015,16 @@ double TMS_Event::GetMuonTrueTrackLength() { std::vector pos = (*it).GetPositionPoints(); int num = 0; - for (auto pnt = pos.begin(); pnt != pos.end(); ++pnt, ++num) { + for (auto pnt = pos.begin(); (pnt+1) != pos.end(); ++pnt, ++num) { auto nextpnt = *(pnt+1); TVector3 point1((*pnt).X(), (*pnt).Y(), (*pnt).Z()); //-200 TVector3 point2(nextpnt.X(), nextpnt.Y(), nextpnt.Z()); //-200 if (TMS_Geom::GetInstance().IsInsideTMS(point1) && TMS_Geom::GetInstance().IsInsideTMS(point2)) { - if ((point2-point1).Mag() > 100) { - continue; - } - double tracklength = TMS_Geom::GetInstance().GetTrackLength(point1, point2); - total += tracklength; + if ((point2-point1).Mag() > 100) { + continue; + } + double tracklength = TMS_Geom::GetInstance().GetTrackLength(point1, point2); + total += tracklength; } } } From aacf9fd3a6ece4a961ca40543b497159b59598c6 Mon Sep 17 00:00:00 2001 From: Jeffrey Kleykamp Date: Mon, 23 Sep 2024 14:06:58 -0500 Subject: [PATCH 4/9] Always make sure iterators are within bounds --- src/TMS_Reco.cpp | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/TMS_Reco.cpp b/src/TMS_Reco.cpp index 9106f3aa..a3a59de8 100644 --- a/src/TMS_Reco.cpp +++ b/src/TMS_Reco.cpp @@ -1260,13 +1260,13 @@ std::vector TMS_TrackFinder::TrackMatching3D() { // Handling cases of two hits in same plane to be matched // By adding a loop into these statements one could also take care of more than 2 hits in the same plane/with the same z coordinate - if (UTracks[itU].GetZ() == UTracks[itU - 1].GetZ()) { + if (itU > 0 && UTracks[itU].GetZ() == UTracks[itU - 1].GetZ()) { CalculateRecoY(UTracks[itU - 1], UTracks[itU - 1], VTracks[itV]); CalculateRecoX(UTracks[itU - 1], VTracks[itV], UTracks[itU - 1]); (aTrack.Hits).push_back(UTracks[itU]); // This adds the original hit if (itU > 0) --itU; // and this allows for the other hit then to be added with the next push_back statement } - if (VTracks[itV].GetZ() == VTracks[itV - 1].GetZ()) { + if (itV > 0 && VTracks[itV].GetZ() == VTracks[itV - 1].GetZ()) { CalculateRecoY(VTracks[itV - 1], UTracks[itU], VTracks[itV - 1]); CalculateRecoX(UTracks[itU], VTracks[itV - 1], VTracks[itV - 1]); (aTrack.Hits).push_back(VTracks[itV]); @@ -1284,20 +1284,20 @@ std::vector TMS_TrackFinder::TrackMatching3D() { CalculateRecoX(UTracks[itU], VTracks[itV], VTracks[itV]); // Handling cases of two hits in same plane - if (UTracks[itU].GetZ() == UTracks[itU - 1].GetZ()) { + if (itU > 0 && UTracks[itU].GetZ() == UTracks[itU - 1].GetZ()) { UTracks[itU - 1].SetRecoY(CompareY(UTracks[itU - 1], VTracks[itV], XTracks[itX])); CalculateRecoX(UTracks[itU - 1], VTracks[itV], UTracks[itU - 1]); (aTrack.Hits).push_back(UTracks[itU]); // This adds the original hit if (itU > 0) --itU; // and this allows for the other hit then to be aded with the next push_back statement } - if (VTracks[itV].GetZ() == VTracks[itV - 1].GetZ()) { + if (itV > 0 && VTracks[itV].GetZ() == VTracks[itV - 1].GetZ()) { VTracks[itV - 1].SetRecoY(CompareY(UTracks[itU], VTracks[itV - 1], XTracks[itX])); CalculateRecoX(UTracks[itU], VTracks[itV - 1], VTracks[itV - 1]); (aTrack.Hits).push_back(VTracks[itV]); if (itV > 0) --itV; } - if (XTracks[itX].GetZ() == XTracks[itX - 1].GetZ()) { + if (itX > 0 && XTracks[itX].GetZ() == XTracks[itX - 1].GetZ()) { XTracks[itX - 1].SetRecoY(XTracks[itX - 1].GetNotZ()); CalculateRecoX(UTracks[itU], VTracks[itV], XTracks[itX - 1]); @@ -1322,13 +1322,13 @@ std::vector TMS_TrackFinder::TrackMatching3D() { std::cout << "same in UV, not X" << std::endl; std::cout << "Hit U: " << UTracks[itU].GetRecoX() << " | " << UTracks[itU].GetRecoY() << " | " << UTracks[itU].GetZ() << " / V: " << VTracks[itV].GetRecoX() << " | " << VTracks[itV].GetRecoY() << " | " << VTracks[itV].GetZ() << " / X: " << XTracks[itX].GetNotZ() << " | " << XTracks[itX].GetZ() << std::endl; #endif - if (UTracks[itU].GetZ() == UTracks[itU - 1].GetZ()) { + if (itU > 0 && UTracks[itU].GetZ() == UTracks[itU - 1].GetZ()) { CalculateRecoY(UTracks[itU - 1], UTracks[itU - 1], VTracks[itV]); CalculateRecoX(UTracks[itU - 1], VTracks[itV], UTracks[itU - 1]); (aTrack.Hits).push_back(UTracks[itU]); // This adds the original hit if (itU > 0) --itU; // and this allows for the other hit then to be added with the next push_back statement } - if (VTracks[itV].GetZ() == VTracks[itV - 1].GetZ()) { + if (itV > 0 && VTracks[itV].GetZ() == VTracks[itV - 1].GetZ()) { CalculateRecoY(VTracks[itV - 1], UTracks[itU], VTracks[itV - 1]); CalculateRecoX(UTracks[itU], VTracks[itV - 1], VTracks[itV - 1]); (aTrack.Hits).push_back(VTracks[itV]); From f360cb6aebfd99cc67fe00e7c6562f4450db98de Mon Sep 17 00:00:00 2001 From: Jeffrey Kleykamp Date: Mon, 23 Sep 2024 15:06:04 -0500 Subject: [PATCH 5/9] Don't look outside range --- src/TMS_Kalman.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/TMS_Kalman.cpp b/src/TMS_Kalman.cpp index 7cc5d71a..1cd5ac53 100644 --- a/src/TMS_Kalman.cpp +++ b/src/TMS_Kalman.cpp @@ -56,7 +56,9 @@ TMS_Kalman::TMS_Kalman(std::vector &Candidates) : } } - const int N_LAYER_BACK = 10; + int N_LAYER_BACK = 10; + // Can't look back further than the first element + if (Candidates.size() < (unsigned)N_LAYER_BACK) N_LAYER_BACK = Candidates.size(); AverageXSlope = (Candidates[Candidates.size() - N_LAYER_BACK].GetRecoX() - Candidates.back().GetRecoX())/(Candidates[Candidates.size() - N_LAYER_BACK].GetZ() - Candidates.back().GetZ()); AverageYSlope = (Candidates.front().GetRecoY() - Candidates.back().GetRecoY())/(Candidates.front().GetZ() - Candidates.back().GetZ()); From bad17a32e3d55fc1fd28e6809e783d48924c2ca2 Mon Sep 17 00:00:00 2001 From: Jeffrey Kleykamp Date: Mon, 23 Sep 2024 15:08:24 -0500 Subject: [PATCH 6/9] Don't update material if density is unknown. See issue 160 --- src/Material.h | 2 ++ src/TMS_Kalman.cpp | 32 +++++++++++++++++++++++++------- 2 files changed, 27 insertions(+), 7 deletions(-) diff --git a/src/Material.h b/src/Material.h index d9b61281..f8b4feb9 100644 --- a/src/Material.h +++ b/src/Material.h @@ -71,6 +71,8 @@ class Material { fMaterialType = kAir; } else { fMaterialType = kUnknown; + std::cout<<"Made material with unknown fMaterialType, given density "< TMS_Const::TMS_End[0]) ) // point outside x region + if ( (PreviousState.x < TMS_Const::TMS_Start[0]) || (PreviousState.x > TMS_Const::TMS_End[0]) ) { // point outside x region std::cerr << "Predicted x value outside TMS: " << PreviousState.y << "\tTMS: [" << TMS_Const::TMS_Start[0] << ", "<< TMS_Const::TMS_End[0] << "]" << std::endl; - if ( (PreviousState.y < TMS_Const::TMS_Start[1]) || (PreviousState.y > TMS_Const::TMS_End[1]) ) // point outside y region + Node.PrintTrueReco(); + PreviousState.Print(); + CurrentState.Print(); + } + if ( (PreviousState.y < TMS_Const::TMS_Start[1]) || (PreviousState.y > TMS_Const::TMS_End[1]) ) { // point outside y region std::cerr << "Predicted y value outside TMS: " << PreviousState.y << "\tTMS: [" << TMS_Const::TMS_Start[1] << ", "<< TMS_Const::TMS_End[1] << "]" << std::endl; + Node.PrintTrueReco(); + PreviousState.Print(); + CurrentState.Print(); + } TVectorD UpdateVec = Transfer*(PreviousVec); @@ -233,9 +241,15 @@ void TMS_Kalman::Predict(TMS_KalmanNode &Node) { TotalLength += thickness; // Update the Bethe Bloch calculator to use this material - Material matter(density); - - Bethe.fMaterial = matter; + try { + Material matter(density); + Bethe.fMaterial = matter; + // Set the material for the multiple scattering + MSC.fMaterial = matter; + } + catch (std::invalid_argument const& ex) { + std::cout<<"Could not make a material using density "< TMS_Const::TMS_End[0]) ) // point outside x region { std::cerr << "lol x value outside TMS: " << CurrentState.y << "\tTMS: [" << TMS_Const::TMS_Start[0] << ", "<< TMS_Const::TMS_End[0] << "]" << std::endl; + Node.PrintTrueReco(); + PreviousState.Print(); + CurrentState.Print(); } if ( (CurrentState.y < TMS_Const::TMS_Start[1]) || (CurrentState.y > TMS_Const::TMS_End[1]) ) // point outside y region { std::cerr << "lol y value outside TMS: " << CurrentState.y << "\tTMS: [" << TMS_Const::TMS_Start[1] << ", "<< TMS_Const::TMS_End[1] << "]" << std::endl; + Node.PrintTrueReco(); + PreviousState.Print(); + CurrentState.Print(); } Node.SetRecoXY(CurrentState); From 1e62abdcbcb15fe414325259f26426008c7e27a2 Mon Sep 17 00:00:00 2001 From: Jeffrey Kleykamp Date: Tue, 24 Sep 2024 15:33:01 -0500 Subject: [PATCH 7/9] Make kalman stuff less verbose --- src/TMS_Kalman.cpp | 2 +- src/TMS_Reco.cpp | 11 ++++++----- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/src/TMS_Kalman.cpp b/src/TMS_Kalman.cpp index 8e042a01..005d3f03 100644 --- a/src/TMS_Kalman.cpp +++ b/src/TMS_Kalman.cpp @@ -314,7 +314,7 @@ void TMS_Kalman::Predict(TMS_KalmanNode &Node) { CovarianceMatrix(2,2) = 1.50; CovarianceMatrix(3,3) = 2.50; CovarianceMatrix(4,4) = 1.0; - std::cout << "Initialising covariance!" << std::endl; + if (Talk) std::cout << "Initialising covariance!" << std::endl; } Node.FillNoiseMatrix(); // Full the matrix for multiple scattering diff --git a/src/TMS_Reco.cpp b/src/TMS_Reco.cpp index a3a59de8..caa47640 100644 --- a/src/TMS_Reco.cpp +++ b/src/TMS_Reco.cpp @@ -733,19 +733,20 @@ void TMS_TrackFinder::FindTracks(TMS_Event &event) { KalmanFilter = TMS_Kalman(trk.Hits); kalman_reco_mom = KalmanFilter.GetMomentum(); - std::cout << "Kalman filter momentum: " << kalman_reco_mom << " MeV" << std::endl; + bool verbose_kalman = false; + if (verbose_kalman) std::cout << "Kalman filter momentum: " << kalman_reco_mom << " MeV" << std::endl; trk.SetMomentum(kalman_reco_mom); // Fill the momentum of the TMS_Track obj - std::cout << "Kalman filter start pos : " << KalmanFilter.Start[0] << ", " << KalmanFilter.Start[1] << ", " << KalmanFilter.Start[2] << std::endl; + if (verbose_kalman) std::cout << "Kalman filter start pos : " << KalmanFilter.Start[0] << ", " << KalmanFilter.Start[1] << ", " << KalmanFilter.Start[2] << std::endl; trk.SetStartPosition(KalmanFilter.Start[0], KalmanFilter.Start[1], KalmanFilter.Start[2]); // Fill the momentum of the TMS_Track obj - std::cout << "Kalman filter end pos : " << KalmanFilter.End[0] << ", " << KalmanFilter.End[1] << ", " << KalmanFilter.End[2] << std::endl; + if (verbose_kalman) std::cout << "Kalman filter end pos : " << KalmanFilter.End[0] << ", " << KalmanFilter.End[1] << ", " << KalmanFilter.End[2] << std::endl; trk.SetEndPosition(KalmanFilter.End[0], KalmanFilter.End[1], KalmanFilter.End[2]); // Fill the momentum of the TMS_Track obj - std::cout << "Kalman filter start dir : " << KalmanFilter.StartDirection[0] << ", " << KalmanFilter.StartDirection[1] << ", " << KalmanFilter.StartDirection[2] << std::endl; + if (verbose_kalman) std::cout << "Kalman filter start dir : " << KalmanFilter.StartDirection[0] << ", " << KalmanFilter.StartDirection[1] << ", " << KalmanFilter.StartDirection[2] << std::endl; trk.SetStartDirection(KalmanFilter.StartDirection[0], KalmanFilter.StartDirection[1], KalmanFilter.StartDirection[2]); // Fill the momentum of the TMS_Track obj - std::cout << "Kalman filter end dir : " << KalmanFilter.EndDirection[0] << ", " << KalmanFilter.EndDirection[1] << ", " << KalmanFilter.EndDirection[2] << std::endl; + if (verbose_kalman) std::cout << "Kalman filter end dir : " << KalmanFilter.EndDirection[0] << ", " << KalmanFilter.EndDirection[1] << ", " << KalmanFilter.EndDirection[2] << std::endl; trk.SetEndDirection(KalmanFilter.EndDirection[0], KalmanFilter.EndDirection[1], KalmanFilter.EndDirection[2]); // Fill the momentum of the TMS_Track obj trk.KalmanNodes = KalmanFilter.GetKalmanNodes(); // Fill the KalmanNodes of the TMS_Track From 4894e8f25a46aad13e0cfdfa12ade9a06cb6cf46 Mon Sep 17 00:00:00 2001 From: Jeffrey Kleykamp Date: Tue, 24 Sep 2024 16:08:04 -0500 Subject: [PATCH 8/9] Fix capitalized N typo --- src/TMS_TreeWriter.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/TMS_TreeWriter.cpp b/src/TMS_TreeWriter.cpp index c154145e..d312eabb 100644 --- a/src/TMS_TreeWriter.cpp +++ b/src/TMS_TreeWriter.cpp @@ -234,7 +234,7 @@ void TMS_TreeWriter::MakeBranches() { Reco_Tree->Branch("KalmanPos", RecoTrackKalmanPos, "TrackHitPos[nTracks][200][3]/F"); Reco_Tree->Branch("KalmanTruePos", RecoTrackKalmanTruePos, "TrackHitTruePos[nTracks][200][3]/F"); Reco_Tree->Branch("StartDirection", RecoTrackStartDirection, "StartDirection[nTracks][3]/F"); - Reco_Tree->Branch("EndDirectioN", RecoTrackEndDirection, "EndDirection[nTracks][3]/F"); + Reco_Tree->Branch("EndDirection", RecoTrackEndDirection, "EndDirection[nTracks][3]/F"); Reco_Tree->Branch("StartPos", RecoTrackStartPos, "StartPos[nTracks][3]/F"); Reco_Tree->Branch("EndPos", RecoTrackEndPos, "EndPos[nTracks][3]/F"); Reco_Tree->Branch("EnergyRange", RecoTrackEnergyRange, "EnergyRange[nTracks]/F"); From eff42e0b8d98c6ed665bde778d41db40dc5181d2 Mon Sep 17 00:00:00 2001 From: jdkio <68720233+jdkio@users.noreply.github.com> Date: Wed, 25 Sep 2024 10:35:31 -0400 Subject: [PATCH 9/9] Update TMS_Reco.cpp --- src/TMS_Reco.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/TMS_Reco.cpp b/src/TMS_Reco.cpp index caa47640..eeee2847 100644 --- a/src/TMS_Reco.cpp +++ b/src/TMS_Reco.cpp @@ -3797,7 +3797,7 @@ void TMS_TrackFinder::Accumulate(double xhit, double zhit) { */ // Fill the accumulator, but only within bounds // We don't care about lines outside of bounds - if (i > 0 && c_bin > 0 && i < nSlope && c_bin < nIntercept) Accumulator[i][c_bin]++; + if (i >= 0 && c_bin >= 0 && i < nSlope && c_bin < nIntercept) Accumulator[i][c_bin]++; } }