From 1cc554a0541f788443eea46b20ba50471bdef05d Mon Sep 17 00:00:00 2001 From: Jeffrey Kleykamp Date: Wed, 24 Apr 2024 15:43:36 -0500 Subject: [PATCH 01/30] Fix make_hists to work with single files --- scripts/Reco/make_hists.py | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/scripts/Reco/make_hists.py b/scripts/Reco/make_hists.py index 6fec32b5..5abb8351 100644 --- a/scripts/Reco/make_hists.py +++ b/scripts/Reco/make_hists.py @@ -574,18 +574,18 @@ def validate_then_run(args): if nfiles == 0: raise ValueError(f"Did not find any files in {indir}") print(f"Found {nfiles} files in {indir}") - if inlist != "": - # In this case the user specified a text file with the full paths - with open(inlist) as f: - file_data = f.read() - files_to_use = file_data.splitlines() - nfiles = len(files_to_use) - if nfiles == 0: - raise ValueError(f"Did not find any files in {inlist}") - print(f"Found {nfiles} files in {inlist}") - if infile != "": - # In this case, the user specified exactly one file. Usually they'd hadd many files together. - files_to_use = [infile] + if inlist != "": + # In this case the user specified a text file with the full paths + with open(inlist) as f: + file_data = f.read() + files_to_use = file_data.splitlines() + nfiles = len(files_to_use) + if nfiles == 0: + raise ValueError(f"Did not find any files in {inlist}") + print(f"Found {nfiles} files in {inlist}") + if infile != "": + # In this case, the user specified exactly one file. Usually they'd hadd many files together. + files_to_use = [infile] outdir = args.outdir if outdir == "": From df9d5f0525babba150de712d01d84629fdba1520 Mon Sep 17 00:00:00 2001 From: Jeffrey Kleykamp Date: Wed, 24 Apr 2024 15:52:05 -0500 Subject: [PATCH 02/30] Fix energy term of momentum vectors, at least particles in mass list --- src/TMS_TreeWriter.cpp | 17 +++++++++++----- src/TMS_TrueParticle.cpp | 15 ++++++++------ src/TMS_TrueParticle.h | 42 ++++++++++++++++++++++++---------------- 3 files changed, 46 insertions(+), 28 deletions(-) diff --git a/src/TMS_TreeWriter.cpp b/src/TMS_TreeWriter.cpp index 8582bc79..32db45db 100644 --- a/src/TMS_TreeWriter.cpp +++ b/src/TMS_TreeWriter.cpp @@ -312,12 +312,19 @@ void TMS_TreeWriter::MakeBranches() { } static void setMomentum(float *branch, TVector3 momentum, double energy = -9999) { - branch[0] = momentum.X(); - branch[1] = momentum.Y(); - branch[2] = momentum.Z(); + branch[0] = momentum.Px(); + branch[1] = momentum.Py(); + branch[2] = momentum.Pz(); branch[3] = energy; } +static void setMomentum(float *branch, TLorentzVector momentum) { + branch[0] = momentum.Px(); + branch[1] = momentum.Py(); + branch[2] = momentum.Pz(); + branch[3] = momentum.E(); +} + static void setPosition(float *branch, TLorentzVector position) { branch[0] = position.X(); branch[1] = position.Y(); @@ -420,7 +427,7 @@ void TMS_TreeWriter::Fill(TMS_Event &event) { setMomentum(BirthMomentum[index], (*it).GetBirthMomentum(), (*it).GetBirthEnergy()); setPosition(BirthPosition[index], (*it).GetBirthPosition()); - setMomentum(DeathMomentum[index], (*it).GetDeathMomentum()); + setMomentum(DeathMomentum[index], (*it).GetDeathMomentum(), (*it).GetDeathEnergy()); setPosition(DeathPosition[index], (*it).GetDeathPosition()); setMomentum(MomentumZIsLArEnd[index], (*it).GetMomentumZIsLArEnd()); @@ -1125,7 +1132,7 @@ void TMS_TreeWriter::Fill(TMS_Event &event) { std::cout<<"Error: Found true_primary_particle_index < 0. There should be at least one particle creating energy (assuming no dark noise) but instead the index is: "<= TrueParticles.size()) { - std::cout<<"Error: Found true_primary_particle_index >= TrueParticles.size(). This shouldn't happen since each index should point to a true particle"<= TrueParticles.size(). This shouldn't happen since each index should point to a true particle: "< Date: Mon, 29 Apr 2024 20:52:17 -0500 Subject: [PATCH 03/30] Add track length function that uses vectors --- src/TMS_Geom.h | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/src/TMS_Geom.h b/src/TMS_Geom.h index d32359a5..8c600f94 100644 --- a/src/TMS_Geom.h +++ b/src/TMS_Geom.h @@ -359,6 +359,29 @@ class TMS_Geom { TotalPathLength = Scale(TotalPathLength); return TotalPathLength; }; + + // Walks through all the pairs of nodes and returns the total track length + double GetTrackLength(std::vector nodes, bool ignore_y = false) { + double out = 0; + // for unsigned ints, 0-1 = MAX_UNSIGNED_INT, so need to verify that size > 0 + if (nodes.size() > 1) { + // Loop over all pairs of vectors + for (size_t i = 0; i < nodes.size() - 1; i++) { + auto p1 = nodes.at(i); + auto p2 = nodes.at(i+1); + if (ignore_y) { + TVector3 p1_without_y(p1); + TVector3 p2_without_y(p2); + p1_without_y.SetY(-200); + p2_without_y.SetY(-200); + p1 = p1_without_y; + p2 = p2_without_y; + } + out += GetTrackLength(p1, p2); + } + } + return out; + } std::vector > GetNodes(const TVector3 &point1_temp, const TVector3 &point2_temp) { // Make vectors have geometry scale From c1d2d4cf35540b2a8fbe6b7aa8dec2410a2a757d Mon Sep 17 00:00:00 2001 From: Jeffrey Kleykamp Date: Tue, 30 Apr 2024 21:15:17 -0500 Subject: [PATCH 04/30] Fix confusion about geant4 vs genie particles --- src/TMS_TreeWriter.cpp | 48 ++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 46 insertions(+), 2 deletions(-) diff --git a/src/TMS_TreeWriter.cpp b/src/TMS_TreeWriter.cpp index 32db45db..7a5a93a5 100644 --- a/src/TMS_TreeWriter.cpp +++ b/src/TMS_TreeWriter.cpp @@ -410,6 +410,7 @@ void TMS_TreeWriter::Fill(TMS_Event &event) { nTrueParticles = TrueParticles.size(); + if (nTrueParticles > __TMS_MAX_TRUE_PARTICLES__) nTrueParticles = __TMS_MAX_TRUE_PARTICLES__; for (auto it = TrueParticles.begin(); it != TrueParticles.end(); ++it) { int index = it - TrueParticles.begin(); @@ -1125,14 +1126,22 @@ void TMS_TreeWriter::Fill(TMS_Event &event) { true_secondary_particle_index = particle_info.indices[1]; } // Now for the primary index, find the true starting and ending momentum and position - // TODO add back if (true_primary_particle_index < 0) { // Do nothing, this means we didn't find a true particle associated with a reco track // This can't happen unless dark noise existed which is currently doesn't std::cout<<"Error: Found true_primary_particle_index < 0. There should be at least one particle creating energy (assuming no dark noise) but instead the index is: "<= TrueParticles.size()) { - std::cout<<"Error: Found true_primary_particle_index >= TrueParticles.size(). This shouldn't happen since each index should point to a true particle: "<= TrueParticles.size()) { + true_secondary_particle_index = -999999999; + } + RecoTrackTrueVisibleEnergy[itTrack] = total_true_visible_energy; RecoTrackPrimaryParticleIndex[itTrack] = true_primary_particle_index; RecoTrackPrimaryParticleTrueVisibleEnergy[itTrack] = true_primary_visible_energy; @@ -1338,6 +1352,36 @@ void TMS_TreeWriter::Clear() { RecoTrackEnergyDeposit[i] = -999; RecoTrackLength[i] = -999; } + + RecoTrackN = 0; + for (int i = 0; i < __TMS_MAX_LINES__; ++i) { + RecoTrackTrueVisibleEnergy[i] = -999; + RecoTrackPrimaryParticleIndex[i] = -999; + RecoTrackPrimaryParticleTrueVisibleEnergy[i] = -999; + RecoTrackSecondaryParticleIndex[i] = -999; + RecoTrackSecondaryParticleTrueVisibleEnergy[i] = -999; + for (int j = 0; j < 4; ++j) { + RecoTrackPrimaryParticleTrueMomentumTrackStart[i][j] = -999; + RecoTrackPrimaryParticleTruePositionTrackStart[i][j] = -999; + RecoTrackPrimaryParticleTrueMomentumTrackEnd[i][j] = -999; + RecoTrackPrimaryParticleTruePositionTrackEnd[i][j] = -999; + } + } + + nTrueParticles = 0; + for (int i = 0; i < __TMS_MAX_TRUE_PARTICLES__; ++i) { + VertexID[i] = -999; + Parent[i] = -999; + TrackId[i] = -999; + PDG[i] = -999; + TrueVisibleEnergy[i] = -999; + for (int j = 0; j < 4; ++j) { + BirthMomentum[i][j] = -999; + BirthPosition[i][j] = -999; + DeathMomentum[i][j] = -999; + DeathPosition[i][j] = -999; + } + } } From 61c367aab70f95b8dd265939819ba6543677d21d Mon Sep 17 00:00:00 2001 From: Jeffrey Kleykamp Date: Tue, 30 Apr 2024 22:14:41 -0500 Subject: [PATCH 05/30] Add the concept of interesting particles to save info about --- src/TMS_Event.cpp | 52 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 52 insertions(+) diff --git a/src/TMS_Event.cpp b/src/TMS_Event.cpp index 6a791422..e5006f58 100644 --- a/src/TMS_Event.cpp +++ b/src/TMS_Event.cpp @@ -1,5 +1,6 @@ #include "TMS_Event.h" #include "TMS_Readout_Manager.h" +#include "TDatabasePDG.h" #include // Initialise the event counter to 0 @@ -25,6 +26,7 @@ TMS_Event::TMS_Event(TG4Event &event, bool FillEvent) { bool TMSOnly = false; bool TMSLArOnly = false; bool OnlyPrimary = true; + bool OnlyPrimaryOrInteresting = false; bool LightWeight = TMS_Manager::GetInstance().Get_LightWeight_Truth(); /* if (LightWeight) { @@ -35,6 +37,8 @@ TMS_Event::TMS_Event(TG4Event &event, bool FillEvent) { } */ + TDatabasePDG *database = TDatabasePDG::Instance(); + // Save down the event number EventNumber = EventCounter; generator = std::default_random_engine(7890 + EventNumber); @@ -47,6 +51,12 @@ TMS_Event::TMS_Event(TG4Event &event, bool FillEvent) { //CheckIntegrity(); int vtxcounter = 0; + int nPrimary = 0; + int nInteresting = 0; + int nTotal = 0; + int nCharged = 0; + int nHighMomentum = 0; + int nChargedAndLowMomentum = 0; // Loop over the primary vertices for (TG4PrimaryVertexContainer::iterator it = event.Primaries.begin(); it != event.Primaries.end(); ++it) { //std::cout<<"for each event.Primaries "< 1000000000) { + // Numbers above 1000000000 are nuclei, and so aren't in the database + // They are charged though, but unlikely to have enough momentum to go far + isCharged = true; + } else { + auto particle = database->GetParticle(PDGcode); + if (!particle) { + std::cout<<"Warning: Couldn't find pdg code "<Charge()) > 0.2; + } + } + if (traj.Points.size() > 0) { + TVector3 initial_momentum = traj.Points[0].GetMomentum(); + if (initial_momentum.Mag() > 50) isHighMomentum = true; + //if (isCharged && isHighMomentum && !isPrimary) { + // std::cout<<"Found interesting non-primary particle "< mapping_track_to_vertex_id; int vertex_index = 0; From b82046449b79eb6471e780b9045eb626ab2b86dc Mon Sep 17 00:00:00 2001 From: Jeffrey Kleykamp Date: Tue, 30 Apr 2024 22:28:59 -0500 Subject: [PATCH 06/30] Add all pdg masses with function --- src/TMS_Constants.h | 23 +++++++++++++++-------- src/TMS_TrueParticle.h | 9 +-------- 2 files changed, 16 insertions(+), 16 deletions(-) diff --git a/src/TMS_Constants.h b/src/TMS_Constants.h index 84f39a73..212b2205 100644 --- a/src/TMS_Constants.h +++ b/src/TMS_Constants.h @@ -2,6 +2,8 @@ #define _TMSCONSTANTS_H_SEEN_ #include +#include "TDatabasePDG.h" + // Hidden from us inside edep-sim EDepSimRooTrackerKinematicsGenerator.hh // Number of maximum particles in an edep-sim event @@ -11,14 +13,19 @@ // Constants namespace TMS_KinConst { - const double mass_mu = 105.6583755; // Muon mass in MeV/c2 - const double mass_e = 0.510998950; // electron mass in MeV/c2 - const double mass_tau = 1776.86; // tau mass in MeV/c2 - const double mass_pic = 139.57039; // charged pion mass in MeV/c2 - const double mass_pi0 = 134.9768; // neutral pion mass in MeV/c2 - - const double mass_proton = 938.27208816; // proton mass in MeV/c2 - const double mass_neutron = 939.56542052; // neutron mass in MeV/c2 + static double GetMass(int pdg) { + TDatabasePDG *database = TDatabasePDG::Instance(); + double out = 0; + // The only particles not in the database are nuclei, and zero is a fine guess for their mass + auto particle = database->GetParticle(pdg); + if (particle) { + out = particle->Mass() * 1000; // Convert from GeV to MeV + } + return out; + } + + const double mass_mu = GetMass(13); // Used in Kalman filter + } // General constants for the TMS diff --git a/src/TMS_TrueParticle.h b/src/TMS_TrueParticle.h index 9c9a0234..d52469c2 100644 --- a/src/TMS_TrueParticle.h +++ b/src/TMS_TrueParticle.h @@ -145,14 +145,7 @@ class TMS_TrueParticle { TLorentzVector GetMomentumLeavingLAr() { return GetMomentumLeaving(TMS_Geom::StaticIsInsideLAr); }; double GetEnergyFromMomentum(TVector3 momentum) { - double mass = 0; - if (abs(PDG) == 13) mass = TMS_KinConst::mass_mu; - else if (abs(PDG) == 11) mass = TMS_KinConst::mass_e; - else if (abs(PDG) == 15) mass = TMS_KinConst::mass_tau; - else if (abs(PDG) == 211) mass = TMS_KinConst::mass_pic; - else if (abs(PDG) == 111) mass = TMS_KinConst::mass_pi0; - else if (abs(PDG) == 2212) mass = TMS_KinConst::mass_proton; - else if (abs(PDG) == 2112) mass = TMS_KinConst::mass_neutron; + double mass = TMS_KinConst::GetMass(PDG); return sqrt(momentum.Mag2()+mass*mass); } From 39be02f28ed67b586f15cf2d0c39ec6ac82d3f73 Mon Sep 17 00:00:00 2001 From: Jeffrey Kleykamp Date: Wed, 1 May 2024 15:56:01 -0500 Subject: [PATCH 07/30] Add various true track length measures --- src/TMS_Geom.h | 2 ++ src/TMS_TreeWriter.cpp | 66 ++++++++++++++++++++++++++++++++++++++-- src/TMS_TreeWriter.h | 11 +++++++ src/TMS_TrueParticle.cpp | 19 ++++++++++++ src/TMS_TrueParticle.h | 1 + 5 files changed, 97 insertions(+), 2 deletions(-) diff --git a/src/TMS_Geom.h b/src/TMS_Geom.h index 8c600f94..f13c2625 100644 --- a/src/TMS_Geom.h +++ b/src/TMS_Geom.h @@ -70,6 +70,8 @@ class TMS_Geom { bool IsInsideTMSFirstTwoModules(TVector3 position) const { return IsInsideBox(position, GetStartOfTMS(), GetEndOfTMSFirstTwoModules()); }; static bool StaticIsInsideTMSFirstTwoModules(TVector3 position) { return TMS_Geom::GetInstance().IsInsideTMSFirstTwoModules(position); }; + bool IsInsideReasonableSize(TVector3 position) const { return IsInsideBox(position, TVector3(-10000, -10000, 3000), TVector3(10000, 10000, 20000)); }; + static bool StaticIsInsideReasonableSize(TVector3 position) { return TMS_Geom::GetInstance().IsInsideReasonableSize(position); }; // Get the geometry diff --git a/src/TMS_TreeWriter.cpp b/src/TMS_TreeWriter.cpp index 7a5a93a5..137aef10 100644 --- a/src/TMS_TreeWriter.cpp +++ b/src/TMS_TreeWriter.cpp @@ -275,13 +275,33 @@ void TMS_TreeWriter::MakeBranches() { "RecoTrackPrimaryParticleTrueMomentumTrackEnd[RecoTrackN][4]/F"); Truth_Info->Branch("RecoTrackPrimaryParticleTruePositionTrackEnd", RecoTrackPrimaryParticleTruePositionTrackEnd, "RecoTrackPrimaryParticleTruePositionTrackEnd[RecoTrackN][4]/F"); - + + Truth_Info->Branch("RecoTrackPrimaryParticleTrueTrackLengthAsMeasured", + RecoTrackPrimaryParticleTrueTrackLengthAsMeasured, + "RecoTrackPrimaryParticleTrueTrackLengthAsMeasured[RecoTrackN]/F"); + Truth_Info->Branch("RecoTrackPrimaryParticleTrueTrackLengthAsMeasuredIgnoreY", + RecoTrackPrimaryParticleTrueTrackLengthAsMeasuredIgnoreY, + "RecoTrackPrimaryParticleTrueTrackLengthAsMeasuredIgnoreY[RecoTrackN]/F"); + Truth_Info->Branch("RecoTrackPrimaryParticleTrueTrackLength", RecoTrackPrimaryParticleTrueTrackLength, + "RecoTrackPrimaryParticleTrueTrackLength[RecoTrackN]/F"); + Truth_Info->Branch("RecoTrackPrimaryParticleTrueTrackLengthIgnoreY", RecoTrackPrimaryParticleTrueTrackLengthIgnoreY, + "RecoTrackPrimaryParticleTrueTrackLengthIgnoreY[RecoTrackN]/F"); + Truth_Info->Branch("RecoTrackPrimaryParticleTrueTrackLengthInTMS", RecoTrackPrimaryParticleTrueTrackLengthInTMS, + "RecoTrackPrimaryParticleTrueTrackLengthInTMS[RecoTrackN]/F"); + Truth_Info->Branch("RecoTrackPrimaryParticleTrueTrackLengthInTMSIgnoreY", + RecoTrackPrimaryParticleTrueTrackLengthInTMSIgnoreY, + "RecoTrackPrimaryParticleTrueTrackLengthInTMSIgnoreY[RecoTrackN]/F"); + Truth_Info->Branch("nTrueParticles", &nTrueParticles, "nTrueParticles/I"); Truth_Info->Branch("VertexID", VertexID, "VertexID[nTrueParticles]/I"); Truth_Info->Branch("Parent", Parent, "Parent[nTrueParticles]/I"); Truth_Info->Branch("TrackId", TrackId, "TrackId[nTrueParticles]/I"); Truth_Info->Branch("PDG", PDG, "PDG[nTrueParticles]/I"); Truth_Info->Branch("TrueVisibleEnergy", TrueVisibleEnergy, "TrueVisibleEnergy[nTrueParticles]/F"); + Truth_Info->Branch("TruePathLength", TruePathLength, "TruePathLength[nTrueParticles]/F"); + Truth_Info->Branch("TruePathLengthIgnoreY", TruePathLengthIgnoreY, "TruePathLengthIgnoreY[nTrueParticles]/F"); + Truth_Info->Branch("TruePathLengthInTMS", TruePathLengthInTMS, "TruePathLengthInTMS[nTrueParticles]/F"); + Truth_Info->Branch("TruePathLengthInTMSIgnoreY", TruePathLengthInTMSIgnoreY, "TruePathLengthInTMSIgnoreY[nTrueParticles]/F"); Truth_Info->Branch("BirthMomentum", BirthMomentum, "BirthMomentum[nTrueParticles][4]/F"); Truth_Info->Branch("BirthPosition", BirthPosition, "BirthPosition[nTrueParticles][4]/F"); @@ -430,7 +450,15 @@ void TMS_TreeWriter::Fill(TMS_Event &event) { setMomentum(DeathMomentum[index], (*it).GetDeathMomentum(), (*it).GetDeathEnergy()); setPosition(DeathPosition[index], (*it).GetDeathPosition()); - + + TruePathLength[index] = TMS_Geom::GetInstance().GetTrackLength((*it).GetPositionPoints(BirthPosition[index][2], DeathPosition[index][2])); + TruePathLengthIgnoreY[index] = + TMS_Geom::GetInstance().GetTrackLength((*it).GetPositionPoints(BirthPosition[index][2], DeathPosition[index][2]), true); + TruePathLengthInTMS[index] = + TMS_Geom::GetInstance().GetTrackLength((*it).GetPositionPoints(BirthPosition[index][2], DeathPosition[index][2], true)); + TruePathLengthInTMSIgnoreY[index] = + TMS_Geom::GetInstance().GetTrackLength((*it).GetPositionPoints(BirthPosition[index][2], DeathPosition[index][2], true), true); + setMomentum(MomentumZIsLArEnd[index], (*it).GetMomentumZIsLArEnd()); setPosition(PositionZIsLArEnd[index], (*it).GetPositionZIsLArEnd()); @@ -1110,6 +1138,11 @@ void TMS_TreeWriter::Fill(TMS_Event &event) { } // Now fill truth info + if (itTrack >= __TMS_MAX_LINES__) { + std::cout<<"Warning: RecoTrackN < __TMS_MAX_LINES__. If this happens often, increase __TMS_MAX_LINES__"< TMS_TrueParticle::GetPositionPoints(double z_start, double z_end, bool onlyInsideTMS) { + std::vector out; + for (size_t i = 0; i < GetPositionPoints().size(); i++) { + double z = GetPositionPoints()[i].Z(); + if (z >= z_start && z <= z_end) { + double x = GetPositionPoints()[i].X(); + double y = GetPositionPoints()[i].Y(); + TVector3 point(x, y, z); + bool add = false; + if (onlyInsideTMS) { + add = TMS_Geom::StaticIsInsideTMS(point); + } + else add = TMS_Geom::StaticIsInsideReasonableSize(point); + if (add) out.push_back(point); + } + } + return out; +} + TLorentzVector TMS_TrueParticle::GetPositionAtZ(double z, double max_z_dist) { TLorentzVector out(-99999999, -99999999, -99999999, -99999999); double z_dist_found = 999999; diff --git a/src/TMS_TrueParticle.h b/src/TMS_TrueParticle.h index d52469c2..69ee98ea 100644 --- a/src/TMS_TrueParticle.h +++ b/src/TMS_TrueParticle.h @@ -127,6 +127,7 @@ class TMS_TrueParticle { TLorentzVector GetPositionEnteringLAr() { return GetPositionEntering(TMS_Geom::StaticIsInsideLAr); }; TLorentzVector GetPositionLeavingLAr() { return GetPositionLeaving(TMS_Geom::StaticIsInsideLAr); }; + std::vector GetPositionPoints(double z_start, double z_end, bool onlyInsideTMS = false); TLorentzVector GetMomentumAtZ(double z, double max_z_dist = 220); // About 2 planes in either direction is the max z distance we'll tolerate, 110mm / thick plane TLorentzVector GetMomentumZIsLArEnd() { return GetMomentumAtZ(TMS_Geom::GetInstance().GetZEndOfLAr()); }; From c9ca334b49cbba0b9b8bc15423c70165985ccdb5 Mon Sep 17 00:00:00 2001 From: Jeffrey Kleykamp Date: Thu, 2 May 2024 22:11:04 -0500 Subject: [PATCH 08/30] Change TMS fiducial to bars only --- src/TMS_Constants.h | 5 +++++ src/TMS_Geom.h | 12 ++++++------ 2 files changed, 11 insertions(+), 6 deletions(-) diff --git a/src/TMS_Constants.h b/src/TMS_Constants.h index 212b2205..f0a9026f 100644 --- a/src/TMS_Constants.h +++ b/src/TMS_Constants.h @@ -62,9 +62,14 @@ namespace TMS_Const { const double LAr_End_Exact[] = {3478.48, 829.282, 9135.88}; // More exact locations of bars + // This seems to contain the steel as well const double TMS_Start_Exact[] = {-3520, -3864, TMS_Thin_Start}; const double TMS_End_Exact[] = {3520, 1159, TMS_Thick_End}; + // Plot TrueHitX,Y,Z and zoom in to see where the last hits are + const double TMS_Start_Bars_Only[] = {-3350, 240, TMS_Thin_Start}; + const double TMS_End_Bars_Only[] = {3350, -2950, TMS_Thick_End}; + // Gap for TMS region that is thin iron layer (mm) const double TMS_Thin_gap = 55; // Gap for TMS region that is thick iron layer (mm) diff --git a/src/TMS_Geom.h b/src/TMS_Geom.h index f13c2625..18166c6d 100644 --- a/src/TMS_Geom.h +++ b/src/TMS_Geom.h @@ -39,13 +39,13 @@ class TMS_Geom { inline double GetYEndOfLAr() const { return TMS_Const::LAr_End_Exact[1] - fiducial_volume_cut; }; inline double GetZEndOfLAr() const { return TMS_Const::LAr_End_Exact[2] - fiducial_volume_cut; }; inline TVector3 GetEndOfLAr() const { return TVector3(GetXEndOfLAr(), GetYEndOfLAr(), GetZEndOfLAr()); }; - inline double GetXStartOfTMS() const { return TMS_Const::TMS_Start_Exact[0] + fiducial_volume_cut;; }; - inline double GetYStartOfTMS() const { return TMS_Const::TMS_Start_Exact[1] + fiducial_volume_cut;; }; - inline double GetZStartOfTMS() const { return TMS_Const::TMS_Thin_Start; }; + inline double GetXStartOfTMS() const { return TMS_Const::TMS_Start_Bars_Only[0]; }; + inline double GetYStartOfTMS() const { return TMS_Const::TMS_Start_Bars_Only[1]; }; + inline double GetZStartOfTMS() const { return TMS_Const::TMS_Start_Bars_Only[2]; }; inline TVector3 GetStartOfTMS() const { return TVector3(GetXStartOfTMS(), GetYStartOfTMS(), GetZStartOfTMS()); }; - inline double GetXEndOfTMS() const { return TMS_Const::TMS_End_Exact[0] - fiducial_volume_cut; }; - inline double GetYEndOfTMS() const { return TMS_Const::TMS_End_Exact[1] - fiducial_volume_cut; }; - inline double GetZEndOfTMS() const { return TMS_Const::TMS_Thick_End - fiducial_volume_cut; }; + inline double GetXEndOfTMS() const { return TMS_Const::TMS_End_Bars_Only[0]; }; + inline double GetYEndOfTMS() const { return TMS_Const::TMS_End_Bars_Only[1]; }; + inline double GetZEndOfTMS() const { return TMS_Const::TMS_End_Bars_Only[2]; }; inline TVector3 GetEndOfTMS() const { return TVector3(GetXEndOfTMS(), GetYEndOfTMS(), GetZEndOfTMS()); }; inline double GetZEndOfTMSThin() const { return TMS_Const::TMS_Thick_Start; }; inline TVector3 GetEndOfTMSThin() const { return TVector3(GetXEndOfTMS(), GetYEndOfTMS(), GetZEndOfTMSThin()); }; From 9076cf14356f51fb49034cd9981adb819da0aba0 Mon Sep 17 00:00:00 2001 From: Jeffrey Kleykamp Date: Thu, 2 May 2024 22:20:39 -0500 Subject: [PATCH 09/30] Make LAr fiducial a config option --- config/TMS_Default_Config.toml | 1 + src/TMS_Geom.h | 2 +- src/TMS_Manager.cpp | 1 + src/TMS_Manager.h | 2 ++ 4 files changed, 5 insertions(+), 1 deletion(-) diff --git a/config/TMS_Default_Config.toml b/config/TMS_Default_Config.toml index 0fc6796e..f8b7f230 100644 --- a/config/TMS_Default_Config.toml +++ b/config/TMS_Default_Config.toml @@ -106,6 +106,7 @@ # LightWeight reduces run time signficantly because we only care about TMS objects (no LAr) [Truth] LightWeight = false + LArFiducialCut = 200.0 # mm # Draw PDF of "event display". Slows down reco considerably [Applications] diff --git a/src/TMS_Geom.h b/src/TMS_Geom.h index 18166c6d..3d2a9776 100644 --- a/src/TMS_Geom.h +++ b/src/TMS_Geom.h @@ -30,7 +30,7 @@ class TMS_Geom { } // Positions are fidicual volume - const double fiducial_volume_cut = 500; + const double fiducial_volume_cut = TMS_Manager::GetInstance().Get_LightWeight_Truth(); inline double GetXStartOfLAr() const { return TMS_Const::LAr_Start_Exact[0] + fiducial_volume_cut; }; inline double GetYStartOfLAr() const { return TMS_Const::LAr_Start_Exact[1] + fiducial_volume_cut; }; inline double GetZStartOfLAr() const { return TMS_Const::LAr_Start_Exact[2] + fiducial_volume_cut; }; diff --git a/src/TMS_Manager.cpp b/src/TMS_Manager.cpp index e9346035..44a11cce 100644 --- a/src/TMS_Manager.cpp +++ b/src/TMS_Manager.cpp @@ -70,6 +70,7 @@ TMS_Manager::TMS_Manager() { _RECO_CALIBRATION_EnergyCalibration = toml::find (data, "Recon", "Calibration", "EnergyCalibration"); _TRUTH_LIGHTWEIGHT = toml::find (data, "Truth", "LightWeight"); + _TRUTH_LAR_FIDUCIAL_CUT = toml::find (data, "Truth", "LArFiducialCut"); _APPLICATIONS_DrawPDF = toml::find (data, "Applications", "DrawPDF"); diff --git a/src/TMS_Manager.h b/src/TMS_Manager.h index 9b013f6f..d10b098d 100644 --- a/src/TMS_Manager.h +++ b/src/TMS_Manager.h @@ -62,6 +62,7 @@ class TMS_Manager { bool Get_Reco_Clustering() { return _RECO_CLUSTERING; }; bool Get_LightWeight_Truth() { return _TRUTH_LIGHTWEIGHT; }; + double Get_LArFiducialCut() { return _TRUTH_LAR_FIDUCIAL_CUT; }; bool Get_Reco_TIME_RunTimeSlicer() { return _RECO_TIME_RunTimeSlicer; }; bool Get_Reco_TIME_RunSimpleTimeSlicer() { return _RECO_TIME_RunSimpleTimeSlicer; }; @@ -130,6 +131,7 @@ class TMS_Manager { // Lightweight trajectory saving (ignore small trajectories and gammas) bool _TRUTH_LIGHTWEIGHT; + double _TRUTH_LAR_FIDUCIAL_CUT; bool _RECO_TIME_RunTimeSlicer; bool _RECO_TIME_RunSimpleTimeSlicer; From b77207c7051201e53f14d7c870662ec9577bd82a Mon Sep 17 00:00:00 2001 From: Jeffrey Kleykamp Date: Mon, 13 May 2024 14:35:29 -0500 Subject: [PATCH 10/30] Add Truth_Spill which has one entry per spill, instead of per slice --- app/ConvertToTMSTree.cpp | 1 + src/TMS_TreeWriter.cpp | 138 +++++++++++++++++++++++++++++++++++++++ src/TMS_TreeWriter.h | 3 + 3 files changed, 142 insertions(+) diff --git a/app/ConvertToTMSTree.cpp b/app/ConvertToTMSTree.cpp index a7fb52b3..f454b39f 100644 --- a/app/ConvertToTMSTree.cpp +++ b/app/ConvertToTMSTree.cpp @@ -124,6 +124,7 @@ bool ConvertToTMSTree(std::string filename, std::string output_filename) { // Save det sim information TMS_ReadoutTreeWriter::GetWriter().Fill(tms_event); + TMS_TreeWriter::GetWriter().FillSpill(tms_event); int nslices = TMS_TimeSlicer::GetSlicer().RunTimeSlicer(tms_event); diff --git a/src/TMS_TreeWriter.cpp b/src/TMS_TreeWriter.cpp index 137aef10..698b4046 100644 --- a/src/TMS_TreeWriter.cpp +++ b/src/TMS_TreeWriter.cpp @@ -50,6 +50,10 @@ TMS_TreeWriter::TMS_TreeWriter() { Truth_Info->SetDirectory(Output); Truth_Info->SetAutoSave(__TMS_AUTOSAVE__); + Truth_Spill = new TTree("Truth_Spill", "Truth_Spill"); + Truth_Spill->SetDirectory(Output); + Truth_Spill->SetAutoSave(__TMS_AUTOSAVE__); + // Make a metadata branch //TTree *MetaData = new TTree("Meta_Data", "Meta_Data"); @@ -329,6 +333,55 @@ void TMS_TreeWriter::MakeBranches() { Truth_Info->Branch("PositionZIsTMSStart", PositionZIsTMSStart, "PositionZIsTMSStart[nTrueParticles][4]/F"); Truth_Info->Branch("MomentumZIsTMSEnd", MomentumZIsTMSEnd, "MomentumZIsTMSEnd[nTrueParticles][4]/F"); Truth_Info->Branch("PositionZIsTMSEnd", PositionZIsTMSEnd, "PositionZIsTMSEnd[nTrueParticles][4]/F"); + + + // Truth information saved per spill only + Truth_Spill->Branch("EventNo", &EventNo, "EventNo/I"); + Truth_Spill->Branch("IsCC", &IsCC, "IsCC/O"); + Truth_Spill->Branch("Interaction", &Reaction); + Truth_Spill->Branch("NeutrinoPDG", &NeutrinoPDG, "NeutrinoPDG/I"); + Truth_Spill->Branch("NeutrinoP4", NeutrinoP4, "NeutrinoP4[4]/F"); + Truth_Spill->Branch("NeutrinoX4", NeutrinoX4, "NeutrinoX4[4]/F"); + + + + Truth_Spill->Branch("nTrueParticles", &nTrueParticles, "nTrueParticles/I"); + Truth_Spill->Branch("VertexID", VertexID, "VertexID[nTrueParticles]/I"); + Truth_Spill->Branch("Parent", Parent, "Parent[nTrueParticles]/I"); + Truth_Spill->Branch("TrackId", TrackId, "TrackId[nTrueParticles]/I"); + Truth_Spill->Branch("PDG", PDG, "PDG[nTrueParticles]/I"); + Truth_Spill->Branch("TrueVisibleEnergy", TrueVisibleEnergy, "TrueVisibleEnergy[nTrueParticles]/F"); + Truth_Spill->Branch("TruePathLength", TruePathLength, "TruePathLength[nTrueParticles]/F"); + Truth_Spill->Branch("TruePathLengthIgnoreY", TruePathLengthIgnoreY, "TruePathLengthIgnoreY[nTrueParticles]/F"); + Truth_Spill->Branch("TruePathLengthInTMS", TruePathLengthInTMS, "TruePathLengthInTMS[nTrueParticles]/F"); + Truth_Spill->Branch("TruePathLengthInTMSIgnoreY", TruePathLengthInTMSIgnoreY, "TruePathLengthInTMSIgnoreY[nTrueParticles]/F"); + + Truth_Spill->Branch("BirthMomentum", BirthMomentum, "BirthMomentum[nTrueParticles][4]/F"); + Truth_Spill->Branch("BirthPosition", BirthPosition, "BirthPosition[nTrueParticles][4]/F"); + Truth_Spill->Branch("DeathMomentum", DeathMomentum, "DeathMomentum[nTrueParticles][4]/F"); + Truth_Spill->Branch("DeathPosition", DeathPosition, "DeathPosition[nTrueParticles][4]/F"); + + // IsInside-based start/end + Truth_Spill->Branch("MomentumLArStart", MomentumLArStart, "MomentumLArStart[nTrueParticles][4]/F"); + Truth_Spill->Branch("PositionLArStart", PositionLArStart, "PositionLArStart[nTrueParticles][4]/F"); + Truth_Spill->Branch("MomentumLArEnd", MomentumLArEnd, "MomentumLArEnd[nTrueParticles][4]/F"); + Truth_Spill->Branch("PositionLArEnd", PositionLArEnd, "PositionLArEnd[nTrueParticles][4]/F"); + Truth_Spill->Branch("MomentumTMSStart", MomentumTMSStart, "MomentumTMSStart[nTrueParticles][4]/F"); + Truth_Spill->Branch("PositionTMSStart", PositionTMSStart, "PositionTMSStart[nTrueParticles][4]/F"); + Truth_Spill->Branch("MomentumTMSFirstTwoModulesEnd", MomentumTMSFirstTwoModulesEnd, "MomentumTMSFirstTwoModulesEnd[nTrueParticles][4]/F"); + Truth_Spill->Branch("PositionTMSFirstTwoModulesEnd", PositionTMSFirstTwoModulesEnd, "PositionTMSFirstTwoModulesEnd[nTrueParticles][4]/F"); + Truth_Spill->Branch("MomentumTMSThinEnd", MomentumTMSThinEnd, "MomentumTMSThinEnd[nTrueParticles][4]/F"); + Truth_Spill->Branch("PositionTMSThinEnd", PositionTMSThinEnd, "PositionTMSThinEnd[nTrueParticles][4]/F"); + Truth_Spill->Branch("MomentumTMSEnd", MomentumTMSEnd, "MomentumTMSEnd[nTrueParticles][4]/F"); + Truth_Spill->Branch("PositionTMSEnd", PositionTMSEnd, "PositionTMSEnd[nTrueParticles][4]/F"); + + // Z-based start/end + Truth_Spill->Branch("MomentumZIsLArEnd", MomentumZIsLArEnd, "MomentumZIsLArEnd[nTrueParticles][4]/F"); + Truth_Spill->Branch("PositionZIsLArEnd", PositionZIsLArEnd, "PositionZIsLArEnd[nTrueParticles][4]/F"); + Truth_Spill->Branch("MomentumZIsTMSStart", MomentumZIsTMSStart, "MomentumZIsTMSStart[nTrueParticles][4]/F"); + Truth_Spill->Branch("PositionZIsTMSStart", PositionZIsTMSStart, "PositionZIsTMSStart[nTrueParticles][4]/F"); + Truth_Spill->Branch("MomentumZIsTMSEnd", MomentumZIsTMSEnd, "MomentumZIsTMSEnd[nTrueParticles][4]/F"); + Truth_Spill->Branch("PositionZIsTMSEnd", PositionZIsTMSEnd, "PositionZIsTMSEnd[nTrueParticles][4]/F"); } static void setMomentum(float *branch, TVector3 momentum, double energy = -9999) { @@ -1227,6 +1280,91 @@ void TMS_TreeWriter::Fill(TMS_Event &event) { Truth_Info->Fill(); } +void TMS_TreeWriter::FillSpill(TMS_Event &event) { + // Clear old info + Clear(); + + // Fill the truth info + EventNo = event.GetEventNumber(); + SliceNo = event.GetSliceNumber(); + SpillNo = event.GetSpillNumber(); + Reaction = event.GetReaction(); + + NeutrinoPDG = event.GetNeutrinoPDG(); + NeutrinoP4[0] = event.GetNeutrinoP4().X(); + NeutrinoP4[1] = event.GetNeutrinoP4().Y(); + NeutrinoP4[2] = event.GetNeutrinoP4().Z(); + NeutrinoP4[3] = event.GetNeutrinoP4().T(); + NeutrinoX4[0] = event.GetNeutrinoX4().X(); + NeutrinoX4[1] = event.GetNeutrinoX4().Y(); + NeutrinoX4[2] = event.GetNeutrinoX4().Z(); + NeutrinoX4[3] = event.GetNeutrinoX4().T(); + NeutrinoP4[0] = event.GetNeutrinoP4().X(); + IsCC = (event.GetReaction().find("[CC]") != std::string::npos); + + std::vector TrueParticles = event.GetTrueParticles(); + nTrueParticles = TrueParticles.size(); + if (nTrueParticles > __TMS_MAX_TRUE_PARTICLES__) nTrueParticles = __TMS_MAX_TRUE_PARTICLES__; + for (auto it = TrueParticles.begin(); it != TrueParticles.end(); ++it) { + int index = it - TrueParticles.begin(); + + if (index >= __TMS_MAX_TRUE_PARTICLES__) { + std::cerr<<"WARNING: Found more particles than __TMS_MAX_TRUE_PARTICLES__. Stopping loop early. If this happens often, increase the max"<Fill(); +} + // Reset the variables void TMS_TreeWriter::Clear() { diff --git a/src/TMS_TreeWriter.h b/src/TMS_TreeWriter.h index 82db3a23..9d7a6922 100644 --- a/src/TMS_TreeWriter.h +++ b/src/TMS_TreeWriter.h @@ -30,12 +30,14 @@ class TMS_TreeWriter { } void Fill(TMS_Event &event); + void FillSpill(TMS_Event &event); void Write() { Output->cd(); Branch_Lines->Write(); Reco_Tree->Write(); Truth_Info->Write(); + Truth_Spill->Write(); std::cout << "TMS_TreeWriter wrote output to " << Output->GetName() << std::endl; Output->Close(); } @@ -62,6 +64,7 @@ class TMS_TreeWriter { TTree* Branch_Lines; // The TTree TTree* Reco_Tree; // The TTree TTree* Truth_Info; // Truth info + TTree* Truth_Spill; // Truth spill void Clear(); void MakeBranches(); // Make the output branches From 7586403e186cff6bf6a476e99de7ce8e634bcd2e Mon Sep 17 00:00:00 2001 From: Jeffrey Kleykamp Date: Mon, 13 May 2024 15:19:51 -0500 Subject: [PATCH 11/30] Add index mapping to Truth_Info and reco trees --- app/ConvertToTMSTree.cpp | 5 ++++- src/TMS_TreeWriter.cpp | 6 +++++- src/TMS_TreeWriter.h | 5 ++++- 3 files changed, 13 insertions(+), 3 deletions(-) diff --git a/app/ConvertToTMSTree.cpp b/app/ConvertToTMSTree.cpp index f454b39f..7d135f9f 100644 --- a/app/ConvertToTMSTree.cpp +++ b/app/ConvertToTMSTree.cpp @@ -74,6 +74,7 @@ bool ConvertToTMSTree(std::string filename, std::string output_filename) { Timer.Start(); int i = 0; + int truth_info_entry_number = 0; //N_entries = 500; // Do we overlay events @@ -124,10 +125,12 @@ bool ConvertToTMSTree(std::string filename, std::string output_filename) { // Save det sim information TMS_ReadoutTreeWriter::GetWriter().Fill(tms_event); - TMS_TreeWriter::GetWriter().FillSpill(tms_event); int nslices = TMS_TimeSlicer::GetSlicer().RunTimeSlicer(tms_event); + TMS_TreeWriter::GetWriter().FillSpill(tms_event, truth_info_entry_number, nslices); + truth_info_entry_number += nslices; + // Could save per spill info here //std::cout<<"Ran time slicer"<Branch("EventNo", &EventNo, "EventNo/I"); Truth_Spill->Branch("IsCC", &IsCC, "IsCC/O"); Truth_Spill->Branch("Interaction", &Reaction); + Truth_Spill->Branch("TruthInfoIndex", &TruthInfoIndex, "TruthInfoIndex/I"); + Truth_Spill->Branch("TruthInfoNSlices", &TruthInfoNSlices, "TruthInfoNSlices/I"); Truth_Spill->Branch("NeutrinoPDG", &NeutrinoPDG, "NeutrinoPDG/I"); Truth_Spill->Branch("NeutrinoP4", NeutrinoP4, "NeutrinoP4[4]/F"); Truth_Spill->Branch("NeutrinoX4", NeutrinoX4, "NeutrinoX4[4]/F"); @@ -1280,7 +1282,7 @@ void TMS_TreeWriter::Fill(TMS_Event &event) { Truth_Info->Fill(); } -void TMS_TreeWriter::FillSpill(TMS_Event &event) { +void TMS_TreeWriter::FillSpill(TMS_Event &event, int truth_info_entry_number, int truth_info_n_slices) { // Clear old info Clear(); @@ -1289,6 +1291,8 @@ void TMS_TreeWriter::FillSpill(TMS_Event &event) { SliceNo = event.GetSliceNumber(); SpillNo = event.GetSpillNumber(); Reaction = event.GetReaction(); + TruthInfoIndex = truth_info_entry_number; + TruthInfoNSlices = truth_info_n_slices; NeutrinoPDG = event.GetNeutrinoPDG(); NeutrinoP4[0] = event.GetNeutrinoP4().X(); diff --git a/src/TMS_TreeWriter.h b/src/TMS_TreeWriter.h index 9d7a6922..c7fbfbf7 100644 --- a/src/TMS_TreeWriter.h +++ b/src/TMS_TreeWriter.h @@ -30,7 +30,7 @@ class TMS_TreeWriter { } void Fill(TMS_Event &event); - void FillSpill(TMS_Event &event); + void FillSpill(TMS_Event &event, int truth_info_entry_number, int truth_info_n_slices); void Write() { Output->cd(); @@ -316,6 +316,9 @@ class TMS_TreeWriter { float RecoTrackPrimaryParticleTrueTrackLengthIgnoreY[__TMS_MAX_LINES__]; float RecoTrackPrimaryParticleTrueTrackLengthInTMS[__TMS_MAX_LINES__]; float RecoTrackPrimaryParticleTrueTrackLengthInTMSIgnoreY[__TMS_MAX_LINES__]; + + int TruthInfoIndex; + int TruthInfoNSlices; }; From d4d0c96a4190d607a24a650eb21c9356fdba6181 Mon Sep 17 00:00:00 2001 From: Jeffrey Kleykamp Date: Mon, 13 May 2024 23:03:41 -0500 Subject: [PATCH 12/30] Fix neutrino X4 units and add to FillTruthFromGRooTracker --- app/CherryPickEvents.cpp | 5 ++++- app/ConvertToTMSTree.cpp | 14 +++++++++++--- app/DrawEvents.cpp | 5 ++++- src/TMS_Event.cpp | 21 +++++++++++---------- src/TMS_Event.h | 5 +++-- 5 files changed, 33 insertions(+), 17 deletions(-) diff --git a/app/CherryPickEvents.cpp b/app/CherryPickEvents.cpp index 03f8a222..32dbbee9 100644 --- a/app/CherryPickEvents.cpp +++ b/app/CherryPickEvents.cpp @@ -57,11 +57,14 @@ int main(int argc, char** argv) { // Get the true neutrino vector from the gRooTracker object int StdHepPdg[__EDEP_SIM_MAX_PART__]; double StdHepP4[__EDEP_SIM_MAX_PART__][4]; + double StdHepX4[__EDEP_SIM_MAX_PART__][4]; gRoo->SetBranchStatus("*", false); gRoo->SetBranchStatus("StdHepPdg", true); gRoo->SetBranchStatus("StdHepP4", true); + gRoo->SetBranchStatus("StdHepX4", true); gRoo->SetBranchAddress("StdHepPdg", StdHepPdg); gRoo->SetBranchAddress("StdHepP4", StdHepP4); + gRoo->SetBranchAddress("StdHepX4", StdHepX4); // Get the event TG4Event *event = NULL; @@ -183,7 +186,7 @@ int main(int argc, char** argv) { } TMS_Event tms_event = TMS_Event(*event, true); - tms_event.FillTruthFromGRooTracker(StdHepPdg, StdHepP4); + tms_event.FillTruthFromGRooTracker(StdHepPdg, StdHepP4, StdHepX4); int pdg = tms_event.GetNeutrinoPDG(); double enu = tms_event.GetNeutrinoP4().E(); diff --git a/app/ConvertToTMSTree.cpp b/app/ConvertToTMSTree.cpp index 7d135f9f..3051bec5 100644 --- a/app/ConvertToTMSTree.cpp +++ b/app/ConvertToTMSTree.cpp @@ -101,7 +101,7 @@ bool ConvertToTMSTree(std::string filename, std::string output_filename) { // Fill up truth information from the GRooTracker object if (gRoo){ - tms_event.FillTruthFromGRooTracker(StdHepPdg, StdHepP4); + tms_event.FillTruthFromGRooTracker(StdHepPdg, StdHepP4, StdHepX4); } // Keep filling up the vector and move on to the next event @@ -128,6 +128,15 @@ bool ConvertToTMSTree(std::string filename, std::string output_filename) { int nslices = TMS_TimeSlicer::GetSlicer().RunTimeSlicer(tms_event); + // Check if this is not pileup + if (event->Primaries.size() == 1 && tms_event.GetNVertices() == 1) { + // Fill the info of the one and only true vertex in the spill + auto primary_vertex = event->Primaries[0]; + int interaction_number = primary_vertex.GetInteractionNumber(); + gRoo->GetEntry(interaction_number); + tms_event.FillTruthFromGRooTracker(StdHepPdg, StdHepP4, StdHepX4); + } + TMS_TreeWriter::GetWriter().FillSpill(tms_event, truth_info_entry_number, nslices); truth_info_entry_number += nslices; @@ -163,8 +172,7 @@ bool ConvertToTMSTree(std::string filename, std::string output_filename) { auto primary_vertex = event->Primaries[primary_vertex_id]; int interaction_number = primary_vertex.GetInteractionNumber(); gRoo->GetEntry(interaction_number); - tms_event_slice.FillTruthFromGRooTracker(StdHepPdg, StdHepP4); - tms_event_slice.FillAdditionalTruthFromGRooTracker(StdHepX4); + tms_event_slice.FillTruthFromGRooTracker(StdHepPdg, StdHepP4, StdHepX4); // And the lepton info int lepton_index = -1; int current_index = 0; diff --git a/app/DrawEvents.cpp b/app/DrawEvents.cpp index c690e22e..7de01320 100644 --- a/app/DrawEvents.cpp +++ b/app/DrawEvents.cpp @@ -56,12 +56,15 @@ int main(int argc, char** argv) { // Get the true neutrino vector from the gRooTracker object int StdHepPdg[__EDEP_SIM_MAX_PART__]; double StdHepP4[__EDEP_SIM_MAX_PART__][4]; + double StdHepX4[__EDEP_SIM_MAX_PART__][4]; if (gRoo) { gRoo->SetBranchStatus("*", false); gRoo->SetBranchStatus("StdHepPdg", true); gRoo->SetBranchStatus("StdHepP4", true); + gRoo->SetBranchStatus("StdHepX4", true); gRoo->SetBranchAddress("StdHepPdg", StdHepPdg); gRoo->SetBranchAddress("StdHepP4", StdHepP4); + gRoo->SetBranchStatus("StdHepX4", true); } // Get the event @@ -187,7 +190,7 @@ int main(int argc, char** argv) { TMS_Event tms_event = TMS_Event(*event, true); if (gRoo){ - tms_event.FillTruthFromGRooTracker(StdHepPdg, StdHepP4); + tms_event.FillTruthFromGRooTracker(StdHepPdg, StdHepP4, StdHepX4); } int pdg = tms_event.GetNeutrinoPDG(); diff --git a/src/TMS_Event.cpp b/src/TMS_Event.cpp index e5006f58..560ea343 100644 --- a/src/TMS_Event.cpp +++ b/src/TMS_Event.cpp @@ -11,6 +11,7 @@ TMS_Event::TMS_Event() { SliceNumber = 0; SpillNumber = -999; nTrueTrajectories = -999; + nVertices = -999; VertexIdOfMostEnergyInEvent = -999; LightWeight = true; } @@ -46,6 +47,7 @@ TMS_Event::TMS_Event(TG4Event &event, bool FillEvent) { SpillNumber = EventCounter; NSlices = 1; // By default there's at least one VertexIdOfMostEnergyInEvent = -999; + nVertices = 0; // Check the integrity of the event //CheckIntegrity(); @@ -218,6 +220,7 @@ TMS_Event::TMS_Event(TG4Event &event, bool FillEvent) { vtxcounter++; } // End if (FillEvent) } // End loop over the primary vertices, for (TG4PrimaryVertexContainer::iterator it + nVertices = vtxcounter; //std::cout<<"N total: "<> GetDeadChannelPositions() { return ChannelPositions; }; std::vector> GetDeadChannelTimes() { return DeadChannelTimes; }; @@ -119,6 +119,7 @@ class TMS_Event { // The number of true trajectories right out of edep-sim // No energy cuts, or number of deposits etc checked int nTrueTrajectories; + int nVertices; std::string Reaction; From c1b8bae83283819ae1e28a4bf5cca439fba0f82e Mon Sep 17 00:00:00 2001 From: Jeffrey Kleykamp Date: Tue, 14 May 2024 20:44:42 -0500 Subject: [PATCH 13/30] Add fiducial as config --- config/TMS_Default_Config.toml | 19 ++++++++++++++++++- src/TMS_Constants.h | 4 ++-- src/TMS_Manager.cpp | 13 +++++++++++++ src/TMS_Manager.h | 26 ++++++++++++++++++++++++++ 4 files changed, 59 insertions(+), 3 deletions(-) diff --git a/config/TMS_Default_Config.toml b/config/TMS_Default_Config.toml index f8b7f230..37627e54 100644 --- a/config/TMS_Default_Config.toml +++ b/config/TMS_Default_Config.toml @@ -100,7 +100,24 @@ [Recon.Calibration] EnergyCalibration = 0.08856 # MeV / PE - + +[Fiducial] + [Fiducial.TMS.Start] + X = -3350.0 # mm + Y = -2950.0 # mm + Z = 11362.0 # mm + [Fiducial.TMS.End] + X = 3350.0 # mm + Y = 240.0 # mm + Z = 18294.0 # mm + [Fiducial.LAr.Start] + X = -3478.48 # mm + Y = -2166.71 # mm + Z = 4179.24 # mm + [Fiducial.LAr.End] + X = 3478.48 # mm + Y = 829.282 # mm + Z = 9135.88 # mm # Options for processing and saving truth information # LightWeight reduces run time signficantly because we only care about TMS objects (no LAr) diff --git a/src/TMS_Constants.h b/src/TMS_Constants.h index f0a9026f..803daf54 100644 --- a/src/TMS_Constants.h +++ b/src/TMS_Constants.h @@ -67,8 +67,8 @@ namespace TMS_Const { const double TMS_End_Exact[] = {3520, 1159, TMS_Thick_End}; // Plot TrueHitX,Y,Z and zoom in to see where the last hits are - const double TMS_Start_Bars_Only[] = {-3350, 240, TMS_Thin_Start}; - const double TMS_End_Bars_Only[] = {3350, -2950, TMS_Thick_End}; + const double TMS_Start_Bars_Only[] = {-3350, -2950, TMS_Thin_Start}; + const double TMS_End_Bars_Only[] = {3350, 240, TMS_Thick_End}; // Gap for TMS region that is thin iron layer (mm) const double TMS_Thin_gap = 55; diff --git a/src/TMS_Manager.cpp b/src/TMS_Manager.cpp index 44a11cce..a37272ee 100644 --- a/src/TMS_Manager.cpp +++ b/src/TMS_Manager.cpp @@ -68,6 +68,19 @@ TMS_Manager::TMS_Manager() { _RECO_CLUSTERING = toml::find (data, "Recon", "Clustering"); _RECO_CALIBRATION_EnergyCalibration = toml::find (data, "Recon", "Calibration", "EnergyCalibration"); + + _FIDUCIAL_TMS_START_X = toml::find(data, "Fiducial", "TMS", "Start", "X"); + _FIDUCIAL_TMS_START_Y = toml::find(data, "Fiducial", "TMS", "Start", "Y"); + _FIDUCIAL_TMS_START_Z = toml::find(data, "Fiducial", "TMS", "Start", "Z"); + _FIDUCIAL_TMS_END_X = toml::find(data, "Fiducial", "TMS", "End", "X"); + _FIDUCIAL_TMS_END_Y = toml::find(data, "Fiducial", "TMS", "End", "Y"); + _FIDUCIAL_TMS_END_Z = toml::find(data, "Fiducial", "TMS", "End", "Z"); + _FIDUCIAL_LAR_START_X = toml::find(data, "Fiducial", "LAr", "Start", "X"); + _FIDUCIAL_LAR_START_Y = toml::find(data, "Fiducial", "LAr", "Start", "Y"); + _FIDUCIAL_LAR_START_Z = toml::find(data, "Fiducial", "LAr", "Start", "Z"); + _FIDUCIAL_LAR_END_X = toml::find(data, "Fiducial", "LAr", "End", "X"); + _FIDUCIAL_LAR_END_Y = toml::find(data, "Fiducial", "LAr", "End", "Y"); + _FIDUCIAL_LAR_END_Z = toml::find(data, "Fiducial", "LAr", "End", "Z"); _TRUTH_LIGHTWEIGHT = toml::find (data, "Truth", "LightWeight"); _TRUTH_LAR_FIDUCIAL_CUT = toml::find (data, "Truth", "LArFiducialCut"); diff --git a/src/TMS_Manager.h b/src/TMS_Manager.h index d10b098d..9f4554ed 100644 --- a/src/TMS_Manager.h +++ b/src/TMS_Manager.h @@ -74,6 +74,19 @@ class TMS_Manager { double Get_RECO_TIME_TimeSlicerMaxTime() { return _RECO_TIME_TimeSlicerMaxTime; }; double Get_RECO_CALIBRATION_EnergyCalibration() { return _RECO_CALIBRATION_EnergyCalibration; }; + + double Get_FIDUCIAL_TMS_START_X() { return _FIDUCIAL_TMS_START_X; }; + double Get_FIDUCIAL_TMS_START_Y() { return _FIDUCIAL_TMS_START_Y; }; + double Get_FIDUCIAL_TMS_START_Z() { return _FIDUCIAL_TMS_START_Z; }; + double Get_FIDUCIAL_TMS_END_X() { return _FIDUCIAL_TMS_END_X; }; + double Get_FIDUCIAL_TMS_END_Y() { return _FIDUCIAL_TMS_END_Y; }; + double Get_FIDUCIAL_TMS_END_Z() { return _FIDUCIAL_TMS_END_Z; }; + double Get_FIDUCIAL_LAR_START_X() { return _FIDUCIAL_LAR_START_X; }; + double Get_FIDUCIAL_LAR_START_Y() { return _FIDUCIAL_LAR_START_Y; }; + double Get_FIDUCIAL_LAR_START_Z() { return _FIDUCIAL_LAR_START_Z; }; + double Get_FIDUCIAL_LAR_END_X() { return _FIDUCIAL_LAR_END_X; }; + double Get_FIDUCIAL_LAR_END_Y() { return _FIDUCIAL_LAR_END_Y; }; + double Get_FIDUCIAL_LAR_END_Z() { return _FIDUCIAL_LAR_END_Z; }; bool Get_DrawPDF() { return _APPLICATIONS_DrawPDF; }; int Get_MaximumNEvents() { return _APPLICATIONS_MaximumNEvents; }; @@ -143,6 +156,19 @@ class TMS_Manager { double _RECO_TIME_TimeSlicerMaxTime; double _RECO_CALIBRATION_EnergyCalibration; + + double _FIDUCIAL_TMS_START_X; + double _FIDUCIAL_TMS_START_Y; + double _FIDUCIAL_TMS_START_Z; + double _FIDUCIAL_TMS_END_X; + double _FIDUCIAL_TMS_END_Y; + double _FIDUCIAL_TMS_END_Z; + double _FIDUCIAL_LAR_START_X; + double _FIDUCIAL_LAR_START_Y; + double _FIDUCIAL_LAR_START_Z; + double _FIDUCIAL_LAR_END_X; + double _FIDUCIAL_LAR_END_Y; + double _FIDUCIAL_LAR_END_Z; bool _APPLICATIONS_DrawPDF; int _APPLICATIONS_MaximumNEvents; From b1c79d8ea7fd07c3d15b2ef3e6163033aa881797 Mon Sep 17 00:00:00 2001 From: Jeffrey Kleykamp Date: Tue, 14 May 2024 20:56:06 -0500 Subject: [PATCH 14/30] Change functions to use config file fiducials --- config/TMS_Default_Config.toml | 8 ++++---- src/TMS_Geom.h | 17 +++++++++++++---- 2 files changed, 17 insertions(+), 8 deletions(-) diff --git a/config/TMS_Default_Config.toml b/config/TMS_Default_Config.toml index 37627e54..85564976 100644 --- a/config/TMS_Default_Config.toml +++ b/config/TMS_Default_Config.toml @@ -103,12 +103,12 @@ [Fiducial] [Fiducial.TMS.Start] - X = -3350.0 # mm - Y = -2950.0 # mm + X = -3300.0 # mm + Y = -2850.0 # mm Z = 11362.0 # mm [Fiducial.TMS.End] - X = 3350.0 # mm - Y = 240.0 # mm + X = 3300.0 # mm + Y = 160 # mm Z = 18294.0 # mm [Fiducial.LAr.Start] X = -3478.48 # mm diff --git a/src/TMS_Geom.h b/src/TMS_Geom.h index 3d2a9776..9ec2a29a 100644 --- a/src/TMS_Geom.h +++ b/src/TMS_Geom.h @@ -51,6 +51,15 @@ class TMS_Geom { inline TVector3 GetEndOfTMSThin() const { return TVector3(GetXEndOfTMS(), GetYEndOfTMS(), GetZEndOfTMSThin()); }; inline double GetZEndOfTMSFirstTwoModules() const { return GetZStartOfTMS() + 110; }; // module 2 - module 0 = 11cm inline TVector3 GetEndOfTMSFirstTwoModules() const { return TVector3(GetXEndOfTMS(), GetYEndOfTMS(), GetZEndOfTMSFirstTwoModules()); }; + + inline TVector3 GetStartOfTMSFiducial() const { return TVector3(TMS_Manager::GetInstance().Get_FIDUCIAL_TMS_START_X(), + TMS_Manager::GetInstance().Get_FIDUCIAL_TMS_START_Y(), TMS_Manager::GetInstance().Get_FIDUCIAL_TMS_START_Y()); }; + inline TVector3 GetEndOfTMSFiducial() const { return TVector3(TMS_Manager::GetInstance().Get_FIDUCIAL_TMS_END_X(), + TMS_Manager::GetInstance().Get_FIDUCIAL_TMS_END_Y(), TMS_Manager::GetInstance().Get_FIDUCIAL_TMS_END_Y()); }; + inline TVector3 GetStartOfLArFiducial() const { return TVector3(TMS_Manager::GetInstance().Get_FIDUCIAL_LAR_START_X(), + TMS_Manager::GetInstance().Get_FIDUCIAL_LAR_START_Y(), TMS_Manager::GetInstance().Get_FIDUCIAL_LAR_START_Y()); }; + inline TVector3 GetEndOfLArFiducial() const { return TVector3(TMS_Manager::GetInstance().Get_FIDUCIAL_LAR_END_X(), + TMS_Manager::GetInstance().Get_FIDUCIAL_LAR_END_Y(), TMS_Manager::GetInstance().Get_FIDUCIAL_LAR_END_Y()); }; bool IsInsideBox(TVector3 position, TVector3 start, TVector3 end) const { if (position.X() < start.X()) return false; @@ -61,13 +70,13 @@ class TMS_Geom { if (position.Z() > end.Z()) return false; return true; }; - bool IsInsideLAr(TVector3 position) const { return IsInsideBox(position, GetStartOfLAr(), GetEndOfLAr()); }; + bool IsInsideLAr(TVector3 position) const { return IsInsideBox(position, GetStartOfLArFiducial(), GetEndOfLArFiducial()); }; static bool StaticIsInsideLAr(TVector3 position) { return TMS_Geom::GetInstance().IsInsideLAr(position); }; - bool IsInsideTMS(TVector3 position) const { return IsInsideBox(position, GetStartOfTMS(), GetEndOfTMS()); }; + bool IsInsideTMS(TVector3 position) const { return IsInsideBox(position, GetStartOfTMSFiducial(), GetEndOfTMSFiducial()); }; static bool StaticIsInsideTMS(TVector3 position) { return TMS_Geom::GetInstance().IsInsideTMS(position); }; - bool IsInsideTMSThin(TVector3 position) const { return IsInsideBox(position, GetStartOfTMS(), GetEndOfTMSThin()); }; + bool IsInsideTMSThin(TVector3 position) const { return IsInsideBox(position, GetStartOfTMSFiducial(), GetEndOfTMSThin()); }; static bool StaticIsInsideTMSThin(TVector3 position) { return TMS_Geom::GetInstance().IsInsideTMSThin(position); }; - bool IsInsideTMSFirstTwoModules(TVector3 position) const { return IsInsideBox(position, GetStartOfTMS(), GetEndOfTMSFirstTwoModules()); }; + bool IsInsideTMSFirstTwoModules(TVector3 position) const { return IsInsideBox(position, GetStartOfTMSFiducial(), GetEndOfTMSFirstTwoModules()); }; static bool StaticIsInsideTMSFirstTwoModules(TVector3 position) { return TMS_Geom::GetInstance().IsInsideTMSFirstTwoModules(position); }; bool IsInsideReasonableSize(TVector3 position) const { return IsInsideBox(position, TVector3(-10000, -10000, 3000), TVector3(10000, 10000, 20000)); }; From b308caf1a24bf85948ab044c1dae67b1e1cd4243 Mon Sep 17 00:00:00 2001 From: Jeffrey Kleykamp Date: Tue, 14 May 2024 21:06:05 -0500 Subject: [PATCH 15/30] Add inside TMS mass function, etc --- config/TMS_Default_Config.toml | 2 +- src/TMS_Geom.h | 24 +++++++++++++++++------- src/TMS_Manager.cpp | 1 - src/TMS_Manager.h | 2 -- 4 files changed, 18 insertions(+), 11 deletions(-) diff --git a/config/TMS_Default_Config.toml b/config/TMS_Default_Config.toml index 85564976..46e2bf8a 100644 --- a/config/TMS_Default_Config.toml +++ b/config/TMS_Default_Config.toml @@ -108,7 +108,7 @@ Z = 11362.0 # mm [Fiducial.TMS.End] X = 3300.0 # mm - Y = 160 # mm + Y = 160.0 # mm Z = 18294.0 # mm [Fiducial.LAr.Start] X = -3478.48 # mm diff --git a/src/TMS_Geom.h b/src/TMS_Geom.h index 9ec2a29a..a4a084ef 100644 --- a/src/TMS_Geom.h +++ b/src/TMS_Geom.h @@ -30,14 +30,13 @@ class TMS_Geom { } // Positions are fidicual volume - const double fiducial_volume_cut = TMS_Manager::GetInstance().Get_LightWeight_Truth(); - inline double GetXStartOfLAr() const { return TMS_Const::LAr_Start_Exact[0] + fiducial_volume_cut; }; - inline double GetYStartOfLAr() const { return TMS_Const::LAr_Start_Exact[1] + fiducial_volume_cut; }; - inline double GetZStartOfLAr() const { return TMS_Const::LAr_Start_Exact[2] + fiducial_volume_cut; }; + inline double GetXStartOfLAr() const { return TMS_Const::LAr_Start_Exact[0]; }; + inline double GetYStartOfLAr() const { return TMS_Const::LAr_Start_Exact[1]; }; + inline double GetZStartOfLAr() const { return TMS_Const::LAr_Start_Exact[2]; }; inline TVector3 GetStartOfLAr() const { return TVector3(GetXStartOfLAr(), GetYStartOfLAr(), GetZStartOfLAr()); }; - inline double GetXEndOfLAr() const { return TMS_Const::LAr_End_Exact[0] - fiducial_volume_cut; }; - inline double GetYEndOfLAr() const { return TMS_Const::LAr_End_Exact[1] - fiducial_volume_cut; }; - inline double GetZEndOfLAr() const { return TMS_Const::LAr_End_Exact[2] - fiducial_volume_cut; }; + inline double GetXEndOfLAr() const { return TMS_Const::LAr_End_Exact[0]; }; + inline double GetYEndOfLAr() const { return TMS_Const::LAr_End_Exact[1]; }; + inline double GetZEndOfLAr() const { return TMS_Const::LAr_End_Exact[2]; }; inline TVector3 GetEndOfLAr() const { return TVector3(GetXEndOfLAr(), GetYEndOfLAr(), GetZEndOfLAr()); }; inline double GetXStartOfTMS() const { return TMS_Const::TMS_Start_Bars_Only[0]; }; inline double GetYStartOfTMS() const { return TMS_Const::TMS_Start_Bars_Only[1]; }; @@ -47,6 +46,14 @@ class TMS_Geom { inline double GetYEndOfTMS() const { return TMS_Const::TMS_End_Bars_Only[1]; }; inline double GetZEndOfTMS() const { return TMS_Const::TMS_End_Bars_Only[2]; }; inline TVector3 GetEndOfTMS() const { return TVector3(GetXEndOfTMS(), GetYEndOfTMS(), GetZEndOfTMS()); }; + inline double GetXStartOfTMSMass() const { return TMS_Const::TMS_Start_Exact[0]; }; + inline double GetYStartOfTMSMass() const { return TMS_Const::TMS_Start_Exact[1]; }; + inline double GetZStartOfTMSMass() const { return TMS_Const::TMS_Start_Exact[2]; }; + inline TVector3 GetStartOfTMSMass() const { return TVector3(GetXStartOfTMS(), GetYStartOfTMSMass(), GetZStartOfTMSMass()); }; + inline double GetXEndOfTMSMass() const { return TMS_Const::TMS_End_Exact[0]; }; + inline double GetYEndOfTMSMass() const { return TMS_Const::TMS_End_Exact[1]; }; + inline double GetZEndOfTMSMass() const { return TMS_Const::TMS_End_Exact[2]; }; + inline TVector3 GetEndOfTMSMass() const { return TVector3(GetXEndOfTMS(), GetYEndOfTMSMass(), GetZEndOfTMSMass()); }; inline double GetZEndOfTMSThin() const { return TMS_Const::TMS_Thick_Start; }; inline TVector3 GetEndOfTMSThin() const { return TVector3(GetXEndOfTMS(), GetYEndOfTMS(), GetZEndOfTMSThin()); }; inline double GetZEndOfTMSFirstTwoModules() const { return GetZStartOfTMS() + 110; }; // module 2 - module 0 = 11cm @@ -78,10 +85,13 @@ class TMS_Geom { static bool StaticIsInsideTMSThin(TVector3 position) { return TMS_Geom::GetInstance().IsInsideTMSThin(position); }; bool IsInsideTMSFirstTwoModules(TVector3 position) const { return IsInsideBox(position, GetStartOfTMSFiducial(), GetEndOfTMSFirstTwoModules()); }; static bool StaticIsInsideTMSFirstTwoModules(TVector3 position) { return TMS_Geom::GetInstance().IsInsideTMSFirstTwoModules(position); }; + bool IsInsideTMSMass(TVector3 position) const { return IsInsideBox(position, GetStartOfTMSMass(), GetEndOfTMSMass()); }; + static bool StaticIsInsideTMSMass(TVector3 position) { return TMS_Geom::GetInstance().IsInsideTMSMass(position); }; bool IsInsideReasonableSize(TVector3 position) const { return IsInsideBox(position, TVector3(-10000, -10000, 3000), TVector3(10000, 10000, 20000)); }; static bool StaticIsInsideReasonableSize(TVector3 position) { return TMS_Geom::GetInstance().IsInsideReasonableSize(position); }; + // Get the geometry TGeoManager* GetGeometry() { diff --git a/src/TMS_Manager.cpp b/src/TMS_Manager.cpp index a37272ee..fa4a9755 100644 --- a/src/TMS_Manager.cpp +++ b/src/TMS_Manager.cpp @@ -83,7 +83,6 @@ TMS_Manager::TMS_Manager() { _FIDUCIAL_LAR_END_Z = toml::find(data, "Fiducial", "LAr", "End", "Z"); _TRUTH_LIGHTWEIGHT = toml::find (data, "Truth", "LightWeight"); - _TRUTH_LAR_FIDUCIAL_CUT = toml::find (data, "Truth", "LArFiducialCut"); _APPLICATIONS_DrawPDF = toml::find (data, "Applications", "DrawPDF"); diff --git a/src/TMS_Manager.h b/src/TMS_Manager.h index 9f4554ed..df1ebb0b 100644 --- a/src/TMS_Manager.h +++ b/src/TMS_Manager.h @@ -62,7 +62,6 @@ class TMS_Manager { bool Get_Reco_Clustering() { return _RECO_CLUSTERING; }; bool Get_LightWeight_Truth() { return _TRUTH_LIGHTWEIGHT; }; - double Get_LArFiducialCut() { return _TRUTH_LAR_FIDUCIAL_CUT; }; bool Get_Reco_TIME_RunTimeSlicer() { return _RECO_TIME_RunTimeSlicer; }; bool Get_Reco_TIME_RunSimpleTimeSlicer() { return _RECO_TIME_RunSimpleTimeSlicer; }; @@ -144,7 +143,6 @@ class TMS_Manager { // Lightweight trajectory saving (ignore small trajectories and gammas) bool _TRUTH_LIGHTWEIGHT; - double _TRUTH_LAR_FIDUCIAL_CUT; bool _RECO_TIME_RunTimeSlicer; bool _RECO_TIME_RunSimpleTimeSlicer; From ab490576005a9b993ea339c2a825507737c75302 Mon Sep 17 00:00:00 2001 From: Jeffrey Kleykamp Date: Tue, 21 May 2024 15:05:21 -0500 Subject: [PATCH 16/30] Fix X4 -> EvtVtx --- app/ConvertToTMSTree.cpp | 14 ++++++++------ src/TMS_Event.cpp | 12 ++++++------ src/TMS_Event.h | 4 ++-- 3 files changed, 16 insertions(+), 14 deletions(-) diff --git a/app/ConvertToTMSTree.cpp b/app/ConvertToTMSTree.cpp index 3051bec5..e6b7dcde 100644 --- a/app/ConvertToTMSTree.cpp +++ b/app/ConvertToTMSTree.cpp @@ -43,17 +43,19 @@ bool ConvertToTMSTree(std::string filename, std::string output_filename) { TG4Event *event = NULL; events->SetBranchAddress("Event", &event); // Get the true neutrino vector from the gRooTracker object + // Nice list of vars to avoid future confusion (internet search for genie groottracker) + // https://internal.dunescience.org/doxygen/read__t2k__rootracker_8C.html int StdHepPdg[__EDEP_SIM_MAX_PART__]; double StdHepP4[__EDEP_SIM_MAX_PART__][4]; - double StdHepX4[__EDEP_SIM_MAX_PART__][4]; + double EvtVtx[__EDEP_SIM_MAX_PART__][4]; if (gRoo){ gRoo->SetBranchStatus("*", false); gRoo->SetBranchStatus("StdHepPdg", true); gRoo->SetBranchStatus("StdHepP4", true); - gRoo->SetBranchStatus("StdHepX4", true); + gRoo->SetBranchStatus("EvtVtx", true); gRoo->SetBranchAddress("StdHepPdg", StdHepPdg); gRoo->SetBranchAddress("StdHepP4", StdHepP4); - gRoo->SetBranchAddress("StdHepX4", StdHepX4); + gRoo->SetBranchAddress("EvtVtx", EvtVtx); } // The global manager TMS_Manager::GetInstance().SetFileName(filename); @@ -101,7 +103,7 @@ bool ConvertToTMSTree(std::string filename, std::string output_filename) { // Fill up truth information from the GRooTracker object if (gRoo){ - tms_event.FillTruthFromGRooTracker(StdHepPdg, StdHepP4, StdHepX4); + tms_event.FillTruthFromGRooTracker(StdHepPdg, StdHepP4, EvtVtx); } // Keep filling up the vector and move on to the next event @@ -134,7 +136,7 @@ bool ConvertToTMSTree(std::string filename, std::string output_filename) { auto primary_vertex = event->Primaries[0]; int interaction_number = primary_vertex.GetInteractionNumber(); gRoo->GetEntry(interaction_number); - tms_event.FillTruthFromGRooTracker(StdHepPdg, StdHepP4, StdHepX4); + tms_event.FillTruthFromGRooTracker(StdHepPdg, StdHepP4, EvtVtx); } TMS_TreeWriter::GetWriter().FillSpill(tms_event, truth_info_entry_number, nslices); @@ -172,7 +174,7 @@ bool ConvertToTMSTree(std::string filename, std::string output_filename) { auto primary_vertex = event->Primaries[primary_vertex_id]; int interaction_number = primary_vertex.GetInteractionNumber(); gRoo->GetEntry(interaction_number); - tms_event_slice.FillTruthFromGRooTracker(StdHepPdg, StdHepP4, StdHepX4); + tms_event_slice.FillTruthFromGRooTracker(StdHepPdg, StdHepP4, EvtVtx); // And the lepton info int lepton_index = -1; int current_index = 0; diff --git a/src/TMS_Event.cpp b/src/TMS_Event.cpp index 560ea343..1b400afd 100644 --- a/src/TMS_Event.cpp +++ b/src/TMS_Event.cpp @@ -843,17 +843,17 @@ void TMS_Event::AddEvent(TMS_Event &Other_Event) { // For now just fill the true neutrino // But shows how you can easily make a vector of rootracker particles for the TMS_Event to carry along -void TMS_Event::FillTruthFromGRooTracker(int pdg[__EDEP_SIM_MAX_PART__], double p4[__EDEP_SIM_MAX_PART__][4], double x4[__EDEP_SIM_MAX_PART__][4]) { +void TMS_Event::FillTruthFromGRooTracker(int pdg[__EDEP_SIM_MAX_PART__], double p4[__EDEP_SIM_MAX_PART__][4], + double vtx[__EDEP_SIM_MAX_PART__][4]) { TrueNeutrino.first.SetX(p4[0][0]); TrueNeutrino.first.SetY(p4[0][1]); TrueNeutrino.first.SetZ(p4[0][2]); TrueNeutrino.first.SetT(p4[0][3]); TrueNeutrino.second = pdg[0]; - // Save in mm - TrueNeutrinoPosition.SetX(x4[0][0] * 1000); - TrueNeutrinoPosition.SetY(x4[0][1] * 1000); - TrueNeutrinoPosition.SetZ(x4[0][2] * 1000); - TrueNeutrinoPosition.SetT(x4[0][3] * 1000); + TrueNeutrinoPosition.SetX(vtx[0][0]); + TrueNeutrinoPosition.SetY(vtx[0][1]); + TrueNeutrinoPosition.SetZ(vtx[0][2]); + TrueNeutrinoPosition.SetT(vtx[0][3]); } void TMS_Event::FillTrueLeptonInfo(int pdg, TLorentzVector position, TLorentzVector momentum) { diff --git a/src/TMS_Event.h b/src/TMS_Event.h index 1f904553..39b07e52 100644 --- a/src/TMS_Event.h +++ b/src/TMS_Event.h @@ -51,11 +51,11 @@ class TMS_Event { std::string GetReaction() { return Reaction; }; // Include some truth metadata, like process, energy, lepton momentum - void FillTruthFromGRooTracker(int pdg[100], double p4[100][4], double x4[100][4]); + void FillTruthFromGRooTracker(int pdg[100], double p4[100][4], double vtx[100][4]); int GetNeutrinoPDG() { return TrueNeutrino.second; }; TLorentzVector GetNeutrinoP4() { return TrueNeutrino.first; }; - TLorentzVector GetNeutrinoX4() { return TrueNeutrinoPosition; }; + TLorentzVector GetNeutrinoVtx() { return TrueNeutrinoPosition; }; int GetLeptonPDG() { return TrueLeptonPDG; }; TLorentzVector GetLeptonX4() { return TrueLeptonPosition; }; TLorentzVector GetLeptonP4() { return TrueLeptonMomentum; }; From d1a62ff023b6820896267921d187e868d29b43c1 Mon Sep 17 00:00:00 2001 From: Jeffrey Kleykamp Date: Tue, 21 May 2024 15:06:33 -0500 Subject: [PATCH 17/30] Add TMS_Spill, which has one entry per spill instead of per slice. Also add various vars related to whether interaction started in TMS, etc --- src/TMS_TreeWriter.cpp | 193 ++++++++++++++++++++------------------- src/TMS_TreeWriter.h | 16 ++++ src/TMS_TrueParticle.cpp | 12 +++ src/TMS_TrueParticle.h | 2 + 4 files changed, 127 insertions(+), 96 deletions(-) diff --git a/src/TMS_TreeWriter.cpp b/src/TMS_TreeWriter.cpp index 8883234c..0c477ffe 100644 --- a/src/TMS_TreeWriter.cpp +++ b/src/TMS_TreeWriter.cpp @@ -238,14 +238,12 @@ void TMS_TreeWriter::MakeBranches() { Reco_Tree->Branch("Length", RecoTrackLength, "Length[nTracks]/F"); + MakeTruthBranches(Truth_Info); + MakeTruthBranches(Truth_Spill); + // Truth information - Truth_Info->Branch("EventNo", &EventNo, "EventNo/I"); - Truth_Info->Branch("IsCC", &IsCC, "IsCC/O"); Truth_Info->Branch("nParticles", &nParticles, "nParticles/I"); - Truth_Info->Branch("Interaction", &Reaction); Truth_Info->Branch("NeutrinoPDG", &NeutrinoPDG, "NeutrinoPDG/I"); - Truth_Info->Branch("NeutrinoP4", NeutrinoP4, "NeutrinoP4[4]/F"); - Truth_Info->Branch("NeutrinoX4", NeutrinoX4, "NeutrinoX4[4]/F"); Truth_Info->Branch("LeptonPDG", &LeptonPDG, "LeptonPDG/I"); Truth_Info->Branch("LeptonP4", LeptonP4, "LeptonP4[4]/F"); Truth_Info->Branch("LeptonX4", LeptonX4, "LeptonX4[4]/F"); @@ -295,95 +293,69 @@ void TMS_TreeWriter::MakeBranches() { Truth_Info->Branch("RecoTrackPrimaryParticleTrueTrackLengthInTMSIgnoreY", RecoTrackPrimaryParticleTrueTrackLengthInTMSIgnoreY, "RecoTrackPrimaryParticleTrueTrackLengthInTMSIgnoreY[RecoTrackN]/F"); +} - Truth_Info->Branch("nTrueParticles", &nTrueParticles, "nTrueParticles/I"); - Truth_Info->Branch("VertexID", VertexID, "VertexID[nTrueParticles]/I"); - Truth_Info->Branch("Parent", Parent, "Parent[nTrueParticles]/I"); - Truth_Info->Branch("TrackId", TrackId, "TrackId[nTrueParticles]/I"); - Truth_Info->Branch("PDG", PDG, "PDG[nTrueParticles]/I"); - Truth_Info->Branch("TrueVisibleEnergy", TrueVisibleEnergy, "TrueVisibleEnergy[nTrueParticles]/F"); - Truth_Info->Branch("TruePathLength", TruePathLength, "TruePathLength[nTrueParticles]/F"); - Truth_Info->Branch("TruePathLengthIgnoreY", TruePathLengthIgnoreY, "TruePathLengthIgnoreY[nTrueParticles]/F"); - Truth_Info->Branch("TruePathLengthInTMS", TruePathLengthInTMS, "TruePathLengthInTMS[nTrueParticles]/F"); - Truth_Info->Branch("TruePathLengthInTMSIgnoreY", TruePathLengthInTMSIgnoreY, "TruePathLengthInTMSIgnoreY[nTrueParticles]/F"); - - Truth_Info->Branch("BirthMomentum", BirthMomentum, "BirthMomentum[nTrueParticles][4]/F"); - Truth_Info->Branch("BirthPosition", BirthPosition, "BirthPosition[nTrueParticles][4]/F"); - Truth_Info->Branch("DeathMomentum", DeathMomentum, "DeathMomentum[nTrueParticles][4]/F"); - Truth_Info->Branch("DeathPosition", DeathPosition, "DeathPosition[nTrueParticles][4]/F"); - - // IsInside-based start/end - Truth_Info->Branch("MomentumLArStart", MomentumLArStart, "MomentumLArStart[nTrueParticles][4]/F"); - Truth_Info->Branch("PositionLArStart", PositionLArStart, "PositionLArStart[nTrueParticles][4]/F"); - Truth_Info->Branch("MomentumLArEnd", MomentumLArEnd, "MomentumLArEnd[nTrueParticles][4]/F"); - Truth_Info->Branch("PositionLArEnd", PositionLArEnd, "PositionLArEnd[nTrueParticles][4]/F"); - Truth_Info->Branch("MomentumTMSStart", MomentumTMSStart, "MomentumTMSStart[nTrueParticles][4]/F"); - Truth_Info->Branch("PositionTMSStart", PositionTMSStart, "PositionTMSStart[nTrueParticles][4]/F"); - Truth_Info->Branch("MomentumTMSFirstTwoModulesEnd", MomentumTMSFirstTwoModulesEnd, "MomentumTMSFirstTwoModulesEnd[nTrueParticles][4]/F"); - Truth_Info->Branch("PositionTMSFirstTwoModulesEnd", PositionTMSFirstTwoModulesEnd, "PositionTMSFirstTwoModulesEnd[nTrueParticles][4]/F"); - Truth_Info->Branch("MomentumTMSThinEnd", MomentumTMSThinEnd, "MomentumTMSThinEnd[nTrueParticles][4]/F"); - Truth_Info->Branch("PositionTMSThinEnd", PositionTMSThinEnd, "PositionTMSThinEnd[nTrueParticles][4]/F"); - Truth_Info->Branch("MomentumTMSEnd", MomentumTMSEnd, "MomentumTMSEnd[nTrueParticles][4]/F"); - Truth_Info->Branch("PositionTMSEnd", PositionTMSEnd, "PositionTMSEnd[nTrueParticles][4]/F"); - - // Z-based start/end - Truth_Info->Branch("MomentumZIsLArEnd", MomentumZIsLArEnd, "MomentumZIsLArEnd[nTrueParticles][4]/F"); - Truth_Info->Branch("PositionZIsLArEnd", PositionZIsLArEnd, "PositionZIsLArEnd[nTrueParticles][4]/F"); - Truth_Info->Branch("MomentumZIsTMSStart", MomentumZIsTMSStart, "MomentumZIsTMSStart[nTrueParticles][4]/F"); - Truth_Info->Branch("PositionZIsTMSStart", PositionZIsTMSStart, "PositionZIsTMSStart[nTrueParticles][4]/F"); - Truth_Info->Branch("MomentumZIsTMSEnd", MomentumZIsTMSEnd, "MomentumZIsTMSEnd[nTrueParticles][4]/F"); - Truth_Info->Branch("PositionZIsTMSEnd", PositionZIsTMSEnd, "PositionZIsTMSEnd[nTrueParticles][4]/F"); - - +void TMS_TreeWriter::MakeTruthBranches(TTree* truth) { // Truth information saved per spill only - Truth_Spill->Branch("EventNo", &EventNo, "EventNo/I"); - Truth_Spill->Branch("IsCC", &IsCC, "IsCC/O"); - Truth_Spill->Branch("Interaction", &Reaction); - Truth_Spill->Branch("TruthInfoIndex", &TruthInfoIndex, "TruthInfoIndex/I"); - Truth_Spill->Branch("TruthInfoNSlices", &TruthInfoNSlices, "TruthInfoNSlices/I"); - Truth_Spill->Branch("NeutrinoPDG", &NeutrinoPDG, "NeutrinoPDG/I"); - Truth_Spill->Branch("NeutrinoP4", NeutrinoP4, "NeutrinoP4[4]/F"); - Truth_Spill->Branch("NeutrinoX4", NeutrinoX4, "NeutrinoX4[4]/F"); - + truth->Branch("EventNo", &EventNo, "EventNo/I"); + truth->Branch("IsCC", &IsCC, "IsCC/O"); + truth->Branch("Interaction", &Reaction); + truth->Branch("TruthInfoIndex", &TruthInfoIndex, "TruthInfoIndex/I"); + truth->Branch("TruthInfoNSlices", &TruthInfoNSlices, "TruthInfoNSlices/I"); + truth->Branch("nPrimaryVertices", &nPrimaryVertices, "nPrimaryVertices/I"); + truth->Branch("HasPileup", &HasPileup, "HasPileup/O"); + truth->Branch("NeutrinoPDG", &NeutrinoPDG, "NeutrinoPDG/I"); + truth->Branch("NeutrinoP4", NeutrinoP4, "NeutrinoP4[4]/F"); + truth->Branch("NeutrinoX4", NeutrinoX4, "NeutrinoX4[4]/F"); + + truth->Branch("nTrueParticles", &nTrueParticles, "nTrueParticles/I"); + truth->Branch("VertexID", VertexID, "VertexID[nTrueParticles]/I"); + truth->Branch("Parent", Parent, "Parent[nTrueParticles]/I"); + truth->Branch("TrackId", TrackId, "TrackId[nTrueParticles]/I"); + truth->Branch("PDG", PDG, "PDG[nTrueParticles]/I"); + truth->Branch("TrueVisibleEnergy", TrueVisibleEnergy, "TrueVisibleEnergy[nTrueParticles]/F"); + truth->Branch("TruePathLength", TruePathLength, "TruePathLength[nTrueParticles]/F"); + truth->Branch("TruePathLengthIgnoreY", TruePathLengthIgnoreY, "TruePathLengthIgnoreY[nTrueParticles]/F"); + truth->Branch("TruePathLengthInTMS", TruePathLengthInTMS, "TruePathLengthInTMS[nTrueParticles]/F"); + truth->Branch("TruePathLengthInTMSIgnoreY", TruePathLengthInTMSIgnoreY, "TruePathLengthInTMSIgnoreY[nTrueParticles]/F"); - - Truth_Spill->Branch("nTrueParticles", &nTrueParticles, "nTrueParticles/I"); - Truth_Spill->Branch("VertexID", VertexID, "VertexID[nTrueParticles]/I"); - Truth_Spill->Branch("Parent", Parent, "Parent[nTrueParticles]/I"); - Truth_Spill->Branch("TrackId", TrackId, "TrackId[nTrueParticles]/I"); - Truth_Spill->Branch("PDG", PDG, "PDG[nTrueParticles]/I"); - Truth_Spill->Branch("TrueVisibleEnergy", TrueVisibleEnergy, "TrueVisibleEnergy[nTrueParticles]/F"); - Truth_Spill->Branch("TruePathLength", TruePathLength, "TruePathLength[nTrueParticles]/F"); - Truth_Spill->Branch("TruePathLengthIgnoreY", TruePathLengthIgnoreY, "TruePathLengthIgnoreY[nTrueParticles]/F"); - Truth_Spill->Branch("TruePathLengthInTMS", TruePathLengthInTMS, "TruePathLengthInTMS[nTrueParticles]/F"); - Truth_Spill->Branch("TruePathLengthInTMSIgnoreY", TruePathLengthInTMSIgnoreY, "TruePathLengthInTMSIgnoreY[nTrueParticles]/F"); + truth->Branch("InteractionTMSFiducial", &InteractionTMSFiducial, "InteractionTMSFiducial/O"); + truth->Branch("InteractionTMSFirstTwoModules", &InteractionTMSFirstTwoModules, "InteractionTMSFirstTwoModules/O"); + truth->Branch("InteractionTMSThin", &InteractionTMSThin, "InteractionTMSThin/O"); + truth->Branch("InteractionLArFiducial", &InteractionLArFiducial, "InteractionLArFiducial/O"); + truth->Branch("TMSFiducialStart", TMSFiducialStart, "TMSFiducialStart[nTrueParticles]/O"); + truth->Branch("TMSFiducialTouch", TMSFiducialTouch, "TMSFiducialTouch[nTrueParticles]/O"); + truth->Branch("TMSFiducialEnd", TMSFiducialEnd, "TMSFiducialEnd[nTrueParticles]/O"); + truth->Branch("LArFiducialStart", LArFiducialStart, "LArFiducialStart[nTrueParticles]/O"); + truth->Branch("LArFiducialTouch", LArFiducialTouch, "LArFiducialTouch[nTrueParticles]/O"); + truth->Branch("LArFiducialEnd", LArFiducialEnd, "LArFiducialEnd[nTrueParticles]/O"); - Truth_Spill->Branch("BirthMomentum", BirthMomentum, "BirthMomentum[nTrueParticles][4]/F"); - Truth_Spill->Branch("BirthPosition", BirthPosition, "BirthPosition[nTrueParticles][4]/F"); - Truth_Spill->Branch("DeathMomentum", DeathMomentum, "DeathMomentum[nTrueParticles][4]/F"); - Truth_Spill->Branch("DeathPosition", DeathPosition, "DeathPosition[nTrueParticles][4]/F"); + truth->Branch("BirthMomentum", BirthMomentum, "BirthMomentum[nTrueParticles][4]/F"); + truth->Branch("BirthPosition", BirthPosition, "BirthPosition[nTrueParticles][4]/F"); + truth->Branch("DeathMomentum", DeathMomentum, "DeathMomentum[nTrueParticles][4]/F"); + truth->Branch("DeathPosition", DeathPosition, "DeathPosition[nTrueParticles][4]/F"); // IsInside-based start/end - Truth_Spill->Branch("MomentumLArStart", MomentumLArStart, "MomentumLArStart[nTrueParticles][4]/F"); - Truth_Spill->Branch("PositionLArStart", PositionLArStart, "PositionLArStart[nTrueParticles][4]/F"); - Truth_Spill->Branch("MomentumLArEnd", MomentumLArEnd, "MomentumLArEnd[nTrueParticles][4]/F"); - Truth_Spill->Branch("PositionLArEnd", PositionLArEnd, "PositionLArEnd[nTrueParticles][4]/F"); - Truth_Spill->Branch("MomentumTMSStart", MomentumTMSStart, "MomentumTMSStart[nTrueParticles][4]/F"); - Truth_Spill->Branch("PositionTMSStart", PositionTMSStart, "PositionTMSStart[nTrueParticles][4]/F"); - Truth_Spill->Branch("MomentumTMSFirstTwoModulesEnd", MomentumTMSFirstTwoModulesEnd, "MomentumTMSFirstTwoModulesEnd[nTrueParticles][4]/F"); - Truth_Spill->Branch("PositionTMSFirstTwoModulesEnd", PositionTMSFirstTwoModulesEnd, "PositionTMSFirstTwoModulesEnd[nTrueParticles][4]/F"); - Truth_Spill->Branch("MomentumTMSThinEnd", MomentumTMSThinEnd, "MomentumTMSThinEnd[nTrueParticles][4]/F"); - Truth_Spill->Branch("PositionTMSThinEnd", PositionTMSThinEnd, "PositionTMSThinEnd[nTrueParticles][4]/F"); - Truth_Spill->Branch("MomentumTMSEnd", MomentumTMSEnd, "MomentumTMSEnd[nTrueParticles][4]/F"); - Truth_Spill->Branch("PositionTMSEnd", PositionTMSEnd, "PositionTMSEnd[nTrueParticles][4]/F"); + truth->Branch("MomentumLArStart", MomentumLArStart, "MomentumLArStart[nTrueParticles][4]/F"); + truth->Branch("PositionLArStart", PositionLArStart, "PositionLArStart[nTrueParticles][4]/F"); + truth->Branch("MomentumLArEnd", MomentumLArEnd, "MomentumLArEnd[nTrueParticles][4]/F"); + truth->Branch("PositionLArEnd", PositionLArEnd, "PositionLArEnd[nTrueParticles][4]/F"); + truth->Branch("MomentumTMSStart", MomentumTMSStart, "MomentumTMSStart[nTrueParticles][4]/F"); + truth->Branch("PositionTMSStart", PositionTMSStart, "PositionTMSStart[nTrueParticles][4]/F"); + truth->Branch("MomentumTMSFirstTwoModulesEnd", MomentumTMSFirstTwoModulesEnd, "MomentumTMSFirstTwoModulesEnd[nTrueParticles][4]/F"); + truth->Branch("PositionTMSFirstTwoModulesEnd", PositionTMSFirstTwoModulesEnd, "PositionTMSFirstTwoModulesEnd[nTrueParticles][4]/F"); + truth->Branch("MomentumTMSThinEnd", MomentumTMSThinEnd, "MomentumTMSThinEnd[nTrueParticles][4]/F"); + truth->Branch("PositionTMSThinEnd", PositionTMSThinEnd, "PositionTMSThinEnd[nTrueParticles][4]/F"); + truth->Branch("MomentumTMSEnd", MomentumTMSEnd, "MomentumTMSEnd[nTrueParticles][4]/F"); + truth->Branch("PositionTMSEnd", PositionTMSEnd, "PositionTMSEnd[nTrueParticles][4]/F"); // Z-based start/end - Truth_Spill->Branch("MomentumZIsLArEnd", MomentumZIsLArEnd, "MomentumZIsLArEnd[nTrueParticles][4]/F"); - Truth_Spill->Branch("PositionZIsLArEnd", PositionZIsLArEnd, "PositionZIsLArEnd[nTrueParticles][4]/F"); - Truth_Spill->Branch("MomentumZIsTMSStart", MomentumZIsTMSStart, "MomentumZIsTMSStart[nTrueParticles][4]/F"); - Truth_Spill->Branch("PositionZIsTMSStart", PositionZIsTMSStart, "PositionZIsTMSStart[nTrueParticles][4]/F"); - Truth_Spill->Branch("MomentumZIsTMSEnd", MomentumZIsTMSEnd, "MomentumZIsTMSEnd[nTrueParticles][4]/F"); - Truth_Spill->Branch("PositionZIsTMSEnd", PositionZIsTMSEnd, "PositionZIsTMSEnd[nTrueParticles][4]/F"); + truth->Branch("MomentumZIsLArEnd", MomentumZIsLArEnd, "MomentumZIsLArEnd[nTrueParticles][4]/F"); + truth->Branch("PositionZIsLArEnd", PositionZIsLArEnd, "PositionZIsLArEnd[nTrueParticles][4]/F"); + truth->Branch("MomentumZIsTMSStart", MomentumZIsTMSStart, "MomentumZIsTMSStart[nTrueParticles][4]/F"); + truth->Branch("PositionZIsTMSStart", PositionZIsTMSStart, "PositionZIsTMSStart[nTrueParticles][4]/F"); + truth->Branch("MomentumZIsTMSEnd", MomentumZIsTMSEnd, "MomentumZIsTMSEnd[nTrueParticles][4]/F"); + truth->Branch("PositionZIsTMSEnd", PositionZIsTMSEnd, "PositionZIsTMSEnd[nTrueParticles][4]/F"); } static void setMomentum(float *branch, TVector3 momentum, double energy = -9999) { @@ -427,10 +399,10 @@ void TMS_TreeWriter::Fill(TMS_Event &event) { NeutrinoP4[1] = event.GetNeutrinoP4().Y(); NeutrinoP4[2] = event.GetNeutrinoP4().Z(); NeutrinoP4[3] = event.GetNeutrinoP4().T(); - NeutrinoX4[0] = event.GetNeutrinoX4().X(); - NeutrinoX4[1] = event.GetNeutrinoX4().Y(); - NeutrinoX4[2] = event.GetNeutrinoX4().Z(); - NeutrinoX4[3] = event.GetNeutrinoX4().T(); + NeutrinoX4[0] = event.GetNeutrinoVtx().X(); + NeutrinoX4[1] = event.GetNeutrinoVtx().Y(); + NeutrinoX4[2] = event.GetNeutrinoVtx().Z(); + NeutrinoX4[3] = event.GetNeutrinoVtx().T(); NeutrinoP4[0] = event.GetNeutrinoP4().X(); IsCC = (event.GetReaction().find("[CC]") != std::string::npos); // Lepton info @@ -482,7 +454,11 @@ void TMS_TreeWriter::Fill(TMS_Event &event) { Muon_Death[3] = (*it).GetDeathPosition().T(); } - + TVector3 interaction_location = event.GetNeutrinoVtx().Vect(); + InteractionTMSFiducial = TMS_Geom::GetInstance().IsInsideTMS(interaction_location); + InteractionTMSFirstTwoModules = TMS_Geom::GetInstance().IsInsideTMSFirstTwoModules(interaction_location); + InteractionTMSThin = TMS_Geom::GetInstance().IsInsideTMSThin(interaction_location); + InteractionLArFiducial = TMS_Geom::GetInstance().IsInsideLAr(interaction_location); nTrueParticles = TrueParticles.size(); if (nTrueParticles > __TMS_MAX_TRUE_PARTICLES__) nTrueParticles = __TMS_MAX_TRUE_PARTICLES__; @@ -499,6 +475,15 @@ void TMS_TreeWriter::Fill(TMS_Event &event) { TrackId[index] = (*it).GetTrackId(); PDG[index] = (*it).GetPDG(); TrueVisibleEnergy[index] = (*it).GetTrueVisibleEnergy(); + + TVector3 location_birth = (*it).GetBirthPosition().Vect(); + TVector3 location_death = (*it).GetDeathPosition().Vect(); + TMSFiducialStart[index] = TMS_Geom::GetInstance().IsInsideTMS(location_birth); + TMSFiducialTouch[index] = (*it).EntersVolume(TMS_Geom::StaticIsInsideTMS); + TMSFiducialEnd[index] = TMS_Geom::GetInstance().IsInsideTMS(location_death); + LArFiducialStart[index] = TMS_Geom::GetInstance().IsInsideLAr(location_birth); + LArFiducialTouch[index] = (*it).EntersVolume(TMS_Geom::StaticIsInsideLAr); + LArFiducialEnd[index] = TMS_Geom::GetInstance().IsInsideLAr(location_death); setMomentum(BirthMomentum[index], (*it).GetBirthMomentum(), (*it).GetBirthEnergy()); setPosition(BirthPosition[index], (*it).GetBirthPosition()); @@ -1293,18 +1278,25 @@ void TMS_TreeWriter::FillSpill(TMS_Event &event, int truth_info_entry_number, in Reaction = event.GetReaction(); TruthInfoIndex = truth_info_entry_number; TruthInfoNSlices = truth_info_n_slices; + HasPileup = event.GetNVertices() != 1; + nPrimaryVertices = event.GetNVertices(); NeutrinoPDG = event.GetNeutrinoPDG(); NeutrinoP4[0] = event.GetNeutrinoP4().X(); NeutrinoP4[1] = event.GetNeutrinoP4().Y(); NeutrinoP4[2] = event.GetNeutrinoP4().Z(); NeutrinoP4[3] = event.GetNeutrinoP4().T(); - NeutrinoX4[0] = event.GetNeutrinoX4().X(); - NeutrinoX4[1] = event.GetNeutrinoX4().Y(); - NeutrinoX4[2] = event.GetNeutrinoX4().Z(); - NeutrinoX4[3] = event.GetNeutrinoX4().T(); - NeutrinoP4[0] = event.GetNeutrinoP4().X(); + NeutrinoX4[0] = event.GetNeutrinoVtx().X(); + NeutrinoX4[1] = event.GetNeutrinoVtx().Y(); + NeutrinoX4[2] = event.GetNeutrinoVtx().Z(); + NeutrinoX4[3] = event.GetNeutrinoVtx().T(); IsCC = (event.GetReaction().find("[CC]") != std::string::npos); + + TVector3 interaction_location = event.GetNeutrinoVtx().Vect(); + InteractionTMSFiducial = TMS_Geom::GetInstance().IsInsideTMS(interaction_location); + InteractionTMSFirstTwoModules = TMS_Geom::GetInstance().IsInsideTMSFirstTwoModules(interaction_location); + InteractionTMSThin = TMS_Geom::GetInstance().IsInsideTMSThin(interaction_location); + InteractionLArFiducial = TMS_Geom::GetInstance().IsInsideLAr(interaction_location); std::vector TrueParticles = event.GetTrueParticles(); nTrueParticles = TrueParticles.size(); @@ -1322,6 +1314,15 @@ void TMS_TreeWriter::FillSpill(TMS_Event &event, int truth_info_entry_number, in TrackId[index] = (*it).GetTrackId(); PDG[index] = (*it).GetPDG(); TrueVisibleEnergy[index] = (*it).GetTrueVisibleEnergy(); + + TVector3 location_birth = (*it).GetBirthPosition().Vect(); + TVector3 location_death = (*it).GetDeathPosition().Vect(); + TMSFiducialStart[index] = TMS_Geom::GetInstance().IsInsideTMS(location_birth); + TMSFiducialTouch[index] = (*it).EntersVolume(TMS_Geom::StaticIsInsideTMS); + TMSFiducialEnd[index] = TMS_Geom::GetInstance().IsInsideTMS(location_death); + LArFiducialStart[index] = TMS_Geom::GetInstance().IsInsideLAr(location_birth); + LArFiducialTouch[index] = (*it).EntersVolume(TMS_Geom::StaticIsInsideLAr); + LArFiducialEnd[index] = TMS_Geom::GetInstance().IsInsideLAr(location_death); setMomentum(BirthMomentum[index], (*it).GetBirthMomentum(), (*it).GetBirthEnergy()); setPosition(BirthPosition[index], (*it).GetBirthPosition()); diff --git a/src/TMS_TreeWriter.h b/src/TMS_TreeWriter.h index c7fbfbf7..248f08a5 100644 --- a/src/TMS_TreeWriter.h +++ b/src/TMS_TreeWriter.h @@ -68,6 +68,7 @@ class TMS_TreeWriter { void Clear(); void MakeBranches(); // Make the output branches + void MakeTruthBranches(TTree* truth); // Make the output branches // The variables int EventNo; @@ -265,6 +266,18 @@ class TMS_TreeWriter { float TruePathLengthInTMS[__TMS_MAX_TRUE_PARTICLES__]; float TruePathLengthInTMSIgnoreY[__TMS_MAX_TRUE_PARTICLES__]; + // Flags for easy use + bool InteractionTMSFiducial; + bool InteractionTMSFirstTwoModules; + bool InteractionTMSThin; + bool InteractionLArFiducial; + bool TMSFiducialStart[__TMS_MAX_TRUE_PARTICLES__]; + bool TMSFiducialTouch[__TMS_MAX_TRUE_PARTICLES__]; + bool TMSFiducialEnd[__TMS_MAX_TRUE_PARTICLES__]; + bool LArFiducialStart[__TMS_MAX_TRUE_PARTICLES__]; + bool LArFiducialTouch[__TMS_MAX_TRUE_PARTICLES__]; + bool LArFiducialEnd[__TMS_MAX_TRUE_PARTICLES__]; + float BirthMomentum[__TMS_MAX_TRUE_PARTICLES__][4]; float BirthPosition[__TMS_MAX_TRUE_PARTICLES__][4]; @@ -319,6 +332,9 @@ class TMS_TreeWriter { int TruthInfoIndex; int TruthInfoNSlices; + + int nPrimaryVertices; + bool HasPileup; }; diff --git a/src/TMS_TrueParticle.cpp b/src/TMS_TrueParticle.cpp index 2a9e274f..3c88ea57 100644 --- a/src/TMS_TrueParticle.cpp +++ b/src/TMS_TrueParticle.cpp @@ -154,3 +154,15 @@ TLorentzVector TMS_TrueParticle::GetMomentumLeaving(IsInsideFunctionType isInsid double energy = GetEnergyFromMomentum(out); return TLorentzVector(out.Px(), out.Py(), out.Pz(), energy); } + +bool TMS_TrueParticle::EntersVolume(IsInsideFunctionType isInside) { + bool out = false; + for (size_t i = 0; i < GetPositionPoints().size(); i++) { + if (isInside(GetPositionPoints()[i].Vect())) { + // Found a single example, we can end our search + out = true; + break; + } + } + return out; +} diff --git a/src/TMS_TrueParticle.h b/src/TMS_TrueParticle.h index 69ee98ea..b16dc634 100644 --- a/src/TMS_TrueParticle.h +++ b/src/TMS_TrueParticle.h @@ -145,6 +145,8 @@ class TMS_TrueParticle { TLorentzVector GetMomentumEnteringLAr() { return GetMomentumEntering(TMS_Geom::StaticIsInsideLAr); }; TLorentzVector GetMomentumLeavingLAr() { return GetMomentumLeaving(TMS_Geom::StaticIsInsideLAr); }; + bool EntersVolume(IsInsideFunctionType isInside); + double GetEnergyFromMomentum(TVector3 momentum) { double mass = TMS_KinConst::GetMass(PDG); return sqrt(momentum.Mag2()+mass*mass); From 215872e950cef7d3d86550f720abe60c33453e4a Mon Sep 17 00:00:00 2001 From: Jeffrey Kleykamp Date: Tue, 21 May 2024 16:27:54 -0500 Subject: [PATCH 18/30] Add additional truth info about reco tracks --- src/TMS_TreeWriter.cpp | 108 ++++++++++++++++++++++++++++++++++++++--- src/TMS_TreeWriter.h | 33 ++++++++++++- 2 files changed, 131 insertions(+), 10 deletions(-) diff --git a/src/TMS_TreeWriter.cpp b/src/TMS_TreeWriter.cpp index 0c477ffe..2a5ff297 100644 --- a/src/TMS_TreeWriter.cpp +++ b/src/TMS_TreeWriter.cpp @@ -284,15 +284,54 @@ void TMS_TreeWriter::MakeBranches() { Truth_Info->Branch("RecoTrackPrimaryParticleTrueTrackLengthAsMeasuredIgnoreY", RecoTrackPrimaryParticleTrueTrackLengthAsMeasuredIgnoreY, "RecoTrackPrimaryParticleTrueTrackLengthAsMeasuredIgnoreY[RecoTrackN]/F"); + Truth_Info->Branch("RecoTrackPrimaryParticleTrueTrackLengthRecoStart", + RecoTrackPrimaryParticleTrueTrackLengthRecoStart, + "RecoTrackPrimaryParticleTrueTrackLengthRecoStart[RecoTrackN]/F"); + Truth_Info->Branch("RecoTrackPrimaryParticleTrueTrackLengthRecoStartIgnoreY", + RecoTrackPrimaryParticleTrueTrackLengthRecoStartIgnoreY, + "RecoTrackPrimaryParticleTrueTrackLengthRecoStartIgnoreY[RecoTrackN]/F"); + Truth_Info->Branch("RecoTrackPrimaryParticleTrueTrackLengthInTMSIgnoreY", + RecoTrackPrimaryParticleTrueTrackLengthInTMSIgnoreY, + "RecoTrackPrimaryParticleTrueTrackLengthInTMSIgnoreY[RecoTrackN]/F"); + + Truth_Info->Branch("RecoTrackPrimaryParticlePDG", &RecoTrackPrimaryParticlePDG, "RecoTrackPrimaryParticlePDG[RecoTrackN]/I"); + Truth_Info->Branch("RecoTrackPrimaryParticleTrueMomentum", RecoTrackPrimaryParticleTrueMomentum, + "RecoTrackPrimaryParticleTrueMomentum[RecoTrackN][4]/F"); + Truth_Info->Branch("RecoTrackPrimaryParticleTruePositionStart", RecoTrackPrimaryParticleTruePositionStart, + "RecoTrackPrimaryParticleTruePositionStart[RecoTrackN][4]/F"); + Truth_Info->Branch("RecoTrackPrimaryParticleTruePositionEnd", RecoTrackPrimaryParticleTruePositionEnd, + "RecoTrackPrimaryParticleTruePositionEnd[RecoTrackN][4]/F"); Truth_Info->Branch("RecoTrackPrimaryParticleTrueTrackLength", RecoTrackPrimaryParticleTrueTrackLength, "RecoTrackPrimaryParticleTrueTrackLength[RecoTrackN]/F"); Truth_Info->Branch("RecoTrackPrimaryParticleTrueTrackLengthIgnoreY", RecoTrackPrimaryParticleTrueTrackLengthIgnoreY, "RecoTrackPrimaryParticleTrueTrackLengthIgnoreY[RecoTrackN]/F"); Truth_Info->Branch("RecoTrackPrimaryParticleTrueTrackLengthInTMS", RecoTrackPrimaryParticleTrueTrackLengthInTMS, "RecoTrackPrimaryParticleTrueTrackLengthInTMS[RecoTrackN]/F"); - Truth_Info->Branch("RecoTrackPrimaryParticleTrueTrackLengthInTMSIgnoreY", - RecoTrackPrimaryParticleTrueTrackLengthInTMSIgnoreY, - "RecoTrackPrimaryParticleTrueTrackLengthInTMSIgnoreY[RecoTrackN]/F"); + Truth_Info->Branch("RecoTrackPrimaryParticleTrueMomentumEnteringTMS", RecoTrackPrimaryParticleTrueMomentumEnteringTMS, + "RecoTrackPrimaryParticleTrueMomentumEnteringTMS[RecoTrackN][4]/F"); + Truth_Info->Branch("RecoTrackPrimaryParticleTrueMomentumLeavingTMS", RecoTrackPrimaryParticleTrueMomentumLeavingTMS, + "RecoTrackPrimaryParticleTrueMomentumLeavingTMS[RecoTrackN][4]/F"); + + Truth_Info->Branch("RecoTrackPrimaryParticleTMSFiducialStart", RecoTrackPrimaryParticleTMSFiducialStart, + "RecoTrackPrimaryParticleTMSFiducialStart[RecoTrackN]/O"); + Truth_Info->Branch("RecoTrackPrimaryParticleTMSFiducialTouch", RecoTrackPrimaryParticleTMSFiducialTouch, + "RecoTrackPrimaryParticleTMSFiducialTouch[RecoTrackN]/O"); + Truth_Info->Branch("RecoTrackPrimaryParticleTMSFiducialEnd", RecoTrackPrimaryParticleTMSFiducialEnd, + "RecoTrackPrimaryParticleTMSFiducialEnd[RecoTrackN]/O"); + Truth_Info->Branch("RecoTrackPrimaryParticleLArFiducialStart", RecoTrackPrimaryParticleLArFiducialStart, + "RecoTrackPrimaryParticleLArFiducialStart[RecoTrackN]/O"); + Truth_Info->Branch("RecoTrackPrimaryParticleLArFiducialTouch", RecoTrackPrimaryParticleLArFiducialTouch, + "RecoTrackPrimaryParticleLArFiducialTouch[RecoTrackN]/O"); + Truth_Info->Branch("RecoTrackPrimaryParticleLArFiducialEnd", RecoTrackPrimaryParticleLArFiducialEnd, + "RecoTrackPrimaryParticleLArFiducialEnd[RecoTrackN]/O"); + + Truth_Info->Branch("RecoTrackSecondaryParticlePDG", &RecoTrackSecondaryParticlePDG, "RecoTrackSecondaryParticlePDG[RecoTrackN]/I"); + Truth_Info->Branch("RecoTrackSecondaryParticleTrueMomentum", RecoTrackSecondaryParticleTrueMomentum, + "RecoTrackSecondaryParticleTrueMomentum[RecoTrackN][4]/F"); + Truth_Info->Branch("RecoTrackSecondaryParticleTruePositionStart", RecoTrackSecondaryParticleTruePositionStart, + "RecoTrackSecondaryParticleTruePositionStart[RecoTrackN][4]/F"); + Truth_Info->Branch("RecoTrackSecondaryParticleTruePositionEnd", RecoTrackSecondaryParticleTruePositionEnd, + "RecoTrackSecondaryParticleTruePositionEnd[RecoTrackN][4]/F"); } void TMS_TreeWriter::MakeTruthBranches(TTree* truth) { @@ -1214,7 +1253,7 @@ void TMS_TreeWriter::Fill(TMS_Event &event) { // Another idea is to save pdg information regardless // We clear the vector elements in the clear function so there's no reason to override it here // But we do need to clear the index so we don't grab the wrong element - true_primary_particle_index = -999999999; + true_primary_particle_index = -888888888; } else { TMS_TrueParticle tp = TrueParticles[true_primary_particle_index]; @@ -1238,15 +1277,42 @@ void TMS_TreeWriter::Fill(TMS_Event &event) { // Or maybe crazy slowdown in sand? // So let's go slightly beyond end of TMS const double LARGE_Z = TMS_Const::TMS_Thick_End + 1000; - RecoTrackPrimaryParticleTrueTrackLength[itTrack] = + const double SMALL_Z = TMS_Const::LAr_Start_Exact[2] - 1000; + RecoTrackPrimaryParticleTrueTrackLengthRecoStart[itTrack] = TMS_Geom::GetInstance().GetTrackLength(tp.GetPositionPoints(start_z, LARGE_Z)); - RecoTrackPrimaryParticleTrueTrackLengthIgnoreY[itTrack] = + RecoTrackPrimaryParticleTrueTrackLengthRecoStartIgnoreY[itTrack] = TMS_Geom::GetInstance().GetTrackLength(tp.GetPositionPoints(start_z, LARGE_Z), true); // Again from true start to particle's end within tms RecoTrackPrimaryParticleTrueTrackLengthInTMS[itTrack] = TMS_Geom::GetInstance().GetTrackLength(tp.GetPositionPoints(start_z, LARGE_Z, true)); RecoTrackPrimaryParticleTrueTrackLengthInTMSIgnoreY[itTrack] = TMS_Geom::GetInstance().GetTrackLength(tp.GetPositionPoints(start_z, LARGE_Z, true), true); + + RecoTrackPrimaryParticleTrueTrackLength[itTrack] = + TMS_Geom::GetInstance().GetTrackLength(tp.GetPositionPoints(SMALL_Z, LARGE_Z)); + RecoTrackPrimaryParticleTrueTrackLengthIgnoreY[itTrack] = + TMS_Geom::GetInstance().GetTrackLength(tp.GetPositionPoints(SMALL_Z, LARGE_Z), true); + + RecoTrackPrimaryParticlePDG[itTrack] = tp.GetPDG(); + setMomentum(RecoTrackPrimaryParticleTrueMomentum[itTrack], tp.GetBirthMomentum()); + setPosition(RecoTrackPrimaryParticleTruePositionStart[itTrack], tp.GetBirthPosition()); + setMomentum(RecoTrackPrimaryParticleTruePositionEnd[itTrack], tp.GetDeathPosition()); + + setMomentum(RecoTrackPrimaryParticleTrueMomentumEnteringTMS[itTrack], tp.GetMomentumEnteringTMS()); + setPosition(RecoTrackPrimaryParticleTruePositionEnteringTMS[itTrack], tp.GetPositionEnteringTMS()); + setMomentum(RecoTrackPrimaryParticleTrueMomentumLeavingTMS[itTrack], tp.GetMomentumLeavingTMS()); + setPosition(RecoTrackPrimaryParticleTruePositionLeavingTMS[itTrack], tp.GetPositionLeavingTMS()); + setMomentum(RecoTrackPrimaryParticleTrueMomentumLeavingLAr[itTrack], tp.GetMomentumLeavingLAr()); + setPosition(RecoTrackPrimaryParticleTruePositionLeavingLAr[itTrack], tp.GetPositionLeavingLAr()); + + TVector3 location_birth = tp.GetBirthPosition().Vect(); + TVector3 location_death = tp.GetDeathPosition().Vect(); + RecoTrackPrimaryParticleTMSFiducialStart[itTrack] = TMS_Geom::GetInstance().IsInsideTMS(location_birth); + RecoTrackPrimaryParticleTMSFiducialTouch[itTrack] = tp.EntersVolume(TMS_Geom::StaticIsInsideTMS); + RecoTrackPrimaryParticleTMSFiducialEnd[itTrack] = TMS_Geom::GetInstance().IsInsideTMS(location_death); + RecoTrackPrimaryParticleLArFiducialStart[itTrack] = TMS_Geom::GetInstance().IsInsideLAr(location_birth); + RecoTrackPrimaryParticleLArFiducialTouch[itTrack] = tp.EntersVolume(TMS_Geom::StaticIsInsideLAr); + RecoTrackPrimaryParticleLArFiducialEnd[itTrack] = TMS_Geom::GetInstance().IsInsideLAr(location_death); } } @@ -1254,6 +1320,13 @@ void TMS_TreeWriter::Fill(TMS_Event &event) { if (true_secondary_particle_index < 0 || (size_t)true_secondary_particle_index >= TrueParticles.size()) { true_secondary_particle_index = -999999999; } + else { + TMS_TrueParticle tp = TrueParticles[true_secondary_particle_index]; + RecoTrackSecondaryParticlePDG[itTrack] = tp.GetPDG(); + setMomentum(RecoTrackSecondaryParticleTrueMomentum[itTrack], tp.GetBirthMomentum()); + setPosition(RecoTrackSecondaryParticleTruePositionStart[itTrack], tp.GetBirthPosition()); + setMomentum(RecoTrackSecondaryParticleTruePositionEnd[itTrack], tp.GetDeathPosition()); + } RecoTrackTrueVisibleEnergy[itTrack] = total_true_visible_energy; RecoTrackPrimaryParticleIndex[itTrack] = true_primary_particle_index; @@ -1560,16 +1633,35 @@ void TMS_TreeWriter::Clear() { RecoTrackPrimaryParticleTrueTrackLengthAsMeasured[i] = -999; RecoTrackPrimaryParticleTrueTrackLengthAsMeasuredIgnoreY[i] = -999; - RecoTrackPrimaryParticleTrueTrackLength[i] = -999; - RecoTrackPrimaryParticleTrueTrackLengthIgnoreY[i] = -999; + RecoTrackPrimaryParticleTrueTrackLengthRecoStart[i] = -999; + RecoTrackPrimaryParticleTrueTrackLengthRecoStartIgnoreY[i] = -999; RecoTrackPrimaryParticleTrueTrackLengthInTMS[i] = -999; RecoTrackPrimaryParticleTrueTrackLengthInTMSIgnoreY[i] = -999; + + RecoTrackPrimaryParticlePDG[i] = -999; + RecoTrackSecondaryParticlePDG[i] = -999; + RecoTrackPrimaryParticleTrueTrackLength[i] = -999; + RecoTrackPrimaryParticleTrueTrackLengthIgnoreY[i] = -999; for (int j = 0; j < 4; ++j) { RecoTrackPrimaryParticleTrueMomentumTrackStart[i][j] = -999; RecoTrackPrimaryParticleTruePositionTrackStart[i][j] = -999; RecoTrackPrimaryParticleTrueMomentumTrackEnd[i][j] = -999; RecoTrackPrimaryParticleTruePositionTrackEnd[i][j] = -999; + + RecoTrackPrimaryParticleTrueMomentum[i][j] = -999; + RecoTrackPrimaryParticleTruePositionStart[i][j] = -999; + RecoTrackPrimaryParticleTruePositionEnd[i][j] = -999; + RecoTrackPrimaryParticleTrueMomentumEnteringTMS[i][j] = -999; + RecoTrackPrimaryParticleTruePositionEnteringTMS[i][j] = -999; + RecoTrackPrimaryParticleTrueMomentumLeavingTMS[i][j] = -999; + RecoTrackPrimaryParticleTruePositionLeavingTMS[i][j] = -999; + RecoTrackPrimaryParticleTrueMomentumLeavingLAr[i][j] = -999; + RecoTrackPrimaryParticleTruePositionLeavingLAr[i][j] = -999; + + RecoTrackSecondaryParticleTrueMomentum[i][j] = -999; + RecoTrackSecondaryParticleTruePositionStart[i][j] = -999; + RecoTrackSecondaryParticleTruePositionEnd[i][j] = -999; } } diff --git a/src/TMS_TreeWriter.h b/src/TMS_TreeWriter.h index 248f08a5..853facb9 100644 --- a/src/TMS_TreeWriter.h +++ b/src/TMS_TreeWriter.h @@ -313,10 +313,39 @@ class TMS_TreeWriter { int RecoTrackN; float RecoTrackTrueVisibleEnergy[__TMS_MAX_LINES__]; + // deprecated, with pileup we can't guarentee a 1-1 relationship int RecoTrackPrimaryParticleIndex[__TMS_MAX_LINES__]; float RecoTrackPrimaryParticleTrueVisibleEnergy[__TMS_MAX_LINES__]; + // deprecated, with pileup we can't guarentee a 1-1 relationship int RecoTrackSecondaryParticleIndex[__TMS_MAX_LINES__]; float RecoTrackSecondaryParticleTrueVisibleEnergy[__TMS_MAX_LINES__]; + + // Save truth info about primary particle + int RecoTrackPrimaryParticlePDG[__TMS_MAX_LINES__]; + float RecoTrackPrimaryParticleTrueMomentum[__TMS_MAX_LINES__][4]; + float RecoTrackPrimaryParticleTruePositionStart[__TMS_MAX_LINES__][4]; + float RecoTrackPrimaryParticleTruePositionEnd[__TMS_MAX_LINES__][4]; + float RecoTrackPrimaryParticleTrueTrackLength[__TMS_MAX_LINES__]; + float RecoTrackPrimaryParticleTrueTrackLengthIgnoreY[__TMS_MAX_LINES__]; + float RecoTrackPrimaryParticleTrueMomentumEnteringTMS[__TMS_MAX_LINES__][4]; + float RecoTrackPrimaryParticleTruePositionEnteringTMS[__TMS_MAX_LINES__][4]; + float RecoTrackPrimaryParticleTrueMomentumLeavingTMS[__TMS_MAX_LINES__][4]; + float RecoTrackPrimaryParticleTruePositionLeavingTMS[__TMS_MAX_LINES__][4]; + float RecoTrackPrimaryParticleTrueMomentumLeavingLAr[__TMS_MAX_LINES__][4]; + float RecoTrackPrimaryParticleTruePositionLeavingLAr[__TMS_MAX_LINES__][4]; + + bool RecoTrackPrimaryParticleTMSFiducialStart[__TMS_MAX_TRUE_PARTICLES__]; + bool RecoTrackPrimaryParticleTMSFiducialTouch[__TMS_MAX_TRUE_PARTICLES__]; + bool RecoTrackPrimaryParticleTMSFiducialEnd[__TMS_MAX_TRUE_PARTICLES__]; + bool RecoTrackPrimaryParticleLArFiducialStart[__TMS_MAX_TRUE_PARTICLES__]; + bool RecoTrackPrimaryParticleLArFiducialTouch[__TMS_MAX_TRUE_PARTICLES__]; + bool RecoTrackPrimaryParticleLArFiducialEnd[__TMS_MAX_TRUE_PARTICLES__]; + + // Save truth info about secondary particle + int RecoTrackSecondaryParticlePDG[__TMS_MAX_LINES__]; + float RecoTrackSecondaryParticleTrueMomentum[__TMS_MAX_LINES__][4]; + float RecoTrackSecondaryParticleTruePositionStart[__TMS_MAX_LINES__][4]; + float RecoTrackSecondaryParticleTruePositionEnd[__TMS_MAX_LINES__][4]; float RecoTrackPrimaryParticleTrueMomentumTrackStart[__TMS_MAX_LINES__][4]; float RecoTrackPrimaryParticleTruePositionTrackStart[__TMS_MAX_LINES__][4]; @@ -325,8 +354,8 @@ class TMS_TreeWriter { float RecoTrackPrimaryParticleTrueTrackLengthAsMeasured[__TMS_MAX_LINES__]; float RecoTrackPrimaryParticleTrueTrackLengthAsMeasuredIgnoreY[__TMS_MAX_LINES__]; - float RecoTrackPrimaryParticleTrueTrackLength[__TMS_MAX_LINES__]; - float RecoTrackPrimaryParticleTrueTrackLengthIgnoreY[__TMS_MAX_LINES__]; + float RecoTrackPrimaryParticleTrueTrackLengthRecoStart[__TMS_MAX_LINES__]; + float RecoTrackPrimaryParticleTrueTrackLengthRecoStartIgnoreY[__TMS_MAX_LINES__]; float RecoTrackPrimaryParticleTrueTrackLengthInTMS[__TMS_MAX_LINES__]; float RecoTrackPrimaryParticleTrueTrackLengthInTMSIgnoreY[__TMS_MAX_LINES__]; From af81140ca7d2eb73db552311a482d6d47b3c00a0 Mon Sep 17 00:00:00 2001 From: Jeffrey Kleykamp Date: Tue, 28 May 2024 17:31:56 -0500 Subject: [PATCH 19/30] Add IsPrimary to TMS_TrueParticle --- src/TMS_TrueParticle.cpp | 22 ++++++++++++---------- src/TMS_TrueParticle.h | 3 ++- 2 files changed, 14 insertions(+), 11 deletions(-) diff --git a/src/TMS_TrueParticle.cpp b/src/TMS_TrueParticle.cpp index 3c88ea57..1eeefd35 100644 --- a/src/TMS_TrueParticle.cpp +++ b/src/TMS_TrueParticle.cpp @@ -1,6 +1,6 @@ #include "TMS_TrueParticle.h" -void TMS_TrueParticle::Print() { +void TMS_TrueParticle::Print(bool small) { std::cout << "Printing TMS_TrueParticle class: " << std::endl; std::cout << " Parent: " << Parent << std::endl; std::cout << " TrackId: " << TrackId << std::endl; @@ -18,15 +18,17 @@ void TMS_TrueParticle::Print() { std::cout << " Size of trajectory points vector: " << PositionPoints.size() << std::endl; std::cout << " Size of momentum at trajectory points vector: " << MomentumPoints.size() << std::endl; - std::cout << " Position of trajectory: " << std::endl; - int it = 0; - for (auto i = PositionPoints.begin(); i != PositionPoints.end(); ++i,++it) { - std::cout << " Point " << it+1 << "/" << PositionPoints.size() << std::endl; - std::cout << " Position: " << std::endl; - (*i).Print(); - std::cout << " Momentum: " << std::endl; - MomentumPoints[it].Print(); - std::cout << " GEANT4 Process, Subprocess: " << Process[it] << ", " << Subprocess[it] << std::endl; + if (!small) { + std::cout << " Position of trajectory: " << std::endl; + int it = 0; + for (auto i = PositionPoints.begin(); i != PositionPoints.end(); ++i,++it) { + std::cout << " Point " << it+1 << "/" << PositionPoints.size() << std::endl; + std::cout << " Position: " << std::endl; + (*i).Print(); + std::cout << " Momentum: " << std::endl; + MomentumPoints[it].Print(); + std::cout << " GEANT4 Process, Subprocess: " << Process[it] << ", " << Subprocess[it] << std::endl; + } } } diff --git a/src/TMS_TrueParticle.h b/src/TMS_TrueParticle.h index b16dc634..d3022430 100644 --- a/src/TMS_TrueParticle.h +++ b/src/TMS_TrueParticle.h @@ -58,7 +58,7 @@ class TMS_TrueParticle { } // Print - void Print(); + void Print(bool small = false); // Add a point in for this true particle void AddPoint(TLorentzVector &Position, TVector3 &Momentum) { @@ -83,6 +83,7 @@ class TMS_TrueParticle { int GetPDG() { return PDG; }; int GetParent() { return Parent; }; + bool IsPrimary() { return Parent < 0; }; int GetTrackId() { return TrackId; }; int GetVertexID() { return VertexID; }; double GetTrueVisibleEnergy() { return TrueVisibleEnergy; }; From cc88ed55572eb255213a642e43eab198b61d916b Mon Sep 17 00:00:00 2001 From: Jeffrey Kleykamp Date: Tue, 28 May 2024 17:34:16 -0500 Subject: [PATCH 20/30] Add primary particles vs non-primary particles. Also save only particle info if primary or deposited energy in TMS --- src/TMS_Event.cpp | 37 +++++++++++++++++++++++++++---- src/TMS_Event.h | 4 ++++ src/TMS_TreeWriter.cpp | 50 +++++++++++++++++++++++++++--------------- src/TMS_TreeWriter.h | 5 +++++ src/TMS_Utils.cpp | 8 +++---- src/TMS_Utils.h | 10 ++++++++- 6 files changed, 87 insertions(+), 27 deletions(-) diff --git a/src/TMS_Event.cpp b/src/TMS_Event.cpp index 1b400afd..82388afe 100644 --- a/src/TMS_Event.cpp +++ b/src/TMS_Event.cpp @@ -16,6 +16,12 @@ TMS_Event::TMS_Event() { LightWeight = true; } + +static bool TMS_TrueParticle_NotWorthSaving(TMS_TrueParticle tp) { + if (tp.GetTrueVisibleEnergy() == 0 && !tp.IsPrimary()) return true; + else return false; +}; + // Start the relatively tedious process of converting into TMS products! // Can also use FillEvent = false to get a simple meta data extractor TMS_Event::TMS_Event(TG4Event &event, bool FillEvent) { @@ -26,8 +32,9 @@ TMS_Event::TMS_Event(TG4Event &event, bool FillEvent) { bool OnlyMuon = false; bool TMSOnly = false; bool TMSLArOnly = false; - bool OnlyPrimary = true; + bool OnlyPrimary = false; bool OnlyPrimaryOrInteresting = false; + bool OnlyPrimaryOrVisibleEnergy = true; bool LightWeight = TMS_Manager::GetInstance().Get_LightWeight_Truth(); /* if (LightWeight) { @@ -124,12 +131,12 @@ TMS_Event::TMS_Event(TG4Event &event, bool FillEvent) { } if (traj.Points.size() > 0) { TVector3 initial_momentum = traj.Points[0].GetMomentum(); - if (initial_momentum.Mag() > 50) isHighMomentum = true; + if (initial_momentum.Mag() > 5) isHighMomentum = true; //if (isCharged && isHighMomentum && !isPrimary) { // std::cout<<"Found interesting non-primary particle "<Branch("RecoTrackPrimaryParticlePDG", &RecoTrackPrimaryParticlePDG, "RecoTrackPrimaryParticlePDG[RecoTrackN]/I"); + Truth_Info->Branch("RecoTrackPrimaryParticleIsPrimary", &RecoTrackPrimaryParticleIsPrimary, "RecoTrackPrimaryParticleIsPrimary[RecoTrackN]/O"); Truth_Info->Branch("RecoTrackPrimaryParticleTrueMomentum", RecoTrackPrimaryParticleTrueMomentum, "RecoTrackPrimaryParticleTrueMomentum[RecoTrackN][4]/F"); Truth_Info->Branch("RecoTrackPrimaryParticleTruePositionStart", RecoTrackPrimaryParticleTruePositionStart, @@ -326,6 +327,7 @@ void TMS_TreeWriter::MakeBranches() { "RecoTrackPrimaryParticleLArFiducialEnd[RecoTrackN]/O"); Truth_Info->Branch("RecoTrackSecondaryParticlePDG", &RecoTrackSecondaryParticlePDG, "RecoTrackSecondaryParticlePDG[RecoTrackN]/I"); + Truth_Info->Branch("RecoTrackSecondaryParticleIsPrimary", &RecoTrackSecondaryParticleIsPrimary, "RecoTrackSecondaryParticleIsPrimary[RecoTrackN]/O"); Truth_Info->Branch("RecoTrackSecondaryParticleTrueMomentum", RecoTrackSecondaryParticleTrueMomentum, "RecoTrackSecondaryParticleTrueMomentum[RecoTrackN][4]/F"); Truth_Info->Branch("RecoTrackSecondaryParticleTruePositionStart", RecoTrackSecondaryParticleTruePositionStart, @@ -348,10 +350,13 @@ void TMS_TreeWriter::MakeTruthBranches(TTree* truth) { truth->Branch("NeutrinoX4", NeutrinoX4, "NeutrinoX4[4]/F"); truth->Branch("nTrueParticles", &nTrueParticles, "nTrueParticles/I"); + truth->Branch("nTruePrimaryParticles", &nTruePrimaryParticles, "nTruePrimaryParticles/I"); + truth->Branch("nTrueForgottenParticles", &nTrueForgottenParticles, "nTrueForgottenParticles/I"); truth->Branch("VertexID", VertexID, "VertexID[nTrueParticles]/I"); truth->Branch("Parent", Parent, "Parent[nTrueParticles]/I"); truth->Branch("TrackId", TrackId, "TrackId[nTrueParticles]/I"); truth->Branch("PDG", PDG, "PDG[nTrueParticles]/I"); + truth->Branch("IsPrimary", IsPrimary, "IsPrimary[nTrueParticles]/O"); truth->Branch("TrueVisibleEnergy", TrueVisibleEnergy, "TrueVisibleEnergy[nTrueParticles]/F"); truth->Branch("TruePathLength", TruePathLength, "TruePathLength[nTrueParticles]/F"); truth->Branch("TruePathLengthIgnoreY", TruePathLengthIgnoreY, "TruePathLengthIgnoreY[nTrueParticles]/F"); @@ -500,12 +505,15 @@ void TMS_TreeWriter::Fill(TMS_Event &event) { InteractionLArFiducial = TMS_Geom::GetInstance().IsInsideLAr(interaction_location); nTrueParticles = TrueParticles.size(); + nTruePrimaryParticles = 0; + nTrueForgottenParticles = event.GetNTrueForgottenParticles(); if (nTrueParticles > __TMS_MAX_TRUE_PARTICLES__) nTrueParticles = __TMS_MAX_TRUE_PARTICLES__; for (auto it = TrueParticles.begin(); it != TrueParticles.end(); ++it) { int index = it - TrueParticles.begin(); if (index >= __TMS_MAX_TRUE_PARTICLES__) { std::cerr<<"WARNING: Found more particles than __TMS_MAX_TRUE_PARTICLES__. Stopping loop early. If this happens often, increase the max"< 0) { true_primary_visible_energy = particle_info.energies[0]; - true_primary_particle_index = particle_info.indices[0]; - } - if (particle_info.energies.size() > 1) { - true_secondary_visible_energy = particle_info.energies[1]; - true_secondary_particle_index = particle_info.indices[1]; + //std::cout<<"checking for primary particle trackid"<= TrueParticles.size()) { - // This can happen if TMS_Event.OnlyPrimary is false - // A primary particle is a genie particle, while secondary particles are particles created by - // interactions in geant4 simulation. - // todo save info of "special" secondary particles, such as secondary muons or michel electrons - // We don't care about the potentially 1000s of low-energy particles, but do care about - // ones long enough to make a track - // Another idea is to save pdg information regardless - // We clear the vector elements in the clear function so there's no reason to override it here - // But we do need to clear the index so we don't grab the wrong element - true_primary_particle_index = -888888888; + // This can happen if TMS_Event.OnlyPrimaryOrVisibleEnergy is false + std::cout<<"Warning: Found a true_primary_particle_index >= TrueParticles.size() case. "<= "< 1) { + true_secondary_visible_energy = particle_info.energies[1]; + true_secondary_particle_index = event.GetTrueParticleIndex(particle_info.trackids[1]); + } if (true_secondary_particle_index < 0 || (size_t)true_secondary_particle_index >= TrueParticles.size()) { true_secondary_particle_index = -999999999; } else { + if (itTrack < __TMS_MAX_LINES__) { TMS_TrueParticle tp = TrueParticles[true_secondary_particle_index]; RecoTrackSecondaryParticlePDG[itTrack] = tp.GetPDG(); + RecoTrackSecondaryParticleIsPrimary[itTrack] = tp.IsPrimary(); setMomentum(RecoTrackSecondaryParticleTrueMomentum[itTrack], tp.GetBirthMomentum()); setPosition(RecoTrackSecondaryParticleTruePositionStart[itTrack], tp.GetBirthPosition()); setMomentum(RecoTrackSecondaryParticleTruePositionEnd[itTrack], tp.GetDeathPosition()); + } } RecoTrackTrueVisibleEnergy[itTrack] = total_true_visible_energy; @@ -1373,6 +1381,8 @@ void TMS_TreeWriter::FillSpill(TMS_Event &event, int truth_info_entry_number, in std::vector TrueParticles = event.GetTrueParticles(); nTrueParticles = TrueParticles.size(); + nTruePrimaryParticles = 0; + nTrueForgottenParticles = event.GetNTrueForgottenParticles(); if (nTrueParticles > __TMS_MAX_TRUE_PARTICLES__) nTrueParticles = __TMS_MAX_TRUE_PARTICLES__; for (auto it = TrueParticles.begin(); it != TrueParticles.end(); ++it) { int index = it - TrueParticles.begin(); @@ -1386,6 +1396,8 @@ void TMS_TreeWriter::FillSpill(TMS_Event &event, int truth_info_entry_number, in Parent[index] = (*it).GetParent(); TrackId[index] = (*it).GetTrackId(); PDG[index] = (*it).GetPDG(); + IsPrimary[index] = (*it).IsPrimary(); + if ((*it).IsPrimary()) nTruePrimaryParticles += 1; TrueVisibleEnergy[index] = (*it).GetTrueVisibleEnergy(); TVector3 location_birth = (*it).GetBirthPosition().Vect(); @@ -1447,7 +1459,7 @@ void TMS_TreeWriter::FillSpill(TMS_Event &event, int truth_info_entry_number, in void TMS_TreeWriter::Clear() { // Reset truth information - EventNo = nParticles = NeutrinoPDG = LeptonPDG = Muon_TrueKE = Muon_TrueTrackLength = VertexIdOfMostEnergyInEvent = -999; + EventNo = nParticles = nTrueParticles = nTruePrimaryParticles = NeutrinoPDG = LeptonPDG = Muon_TrueKE = Muon_TrueTrackLength = VertexIdOfMostEnergyInEvent = -999; VertexIdOfMostEnergyInEvent = VisibleEnergyFromUVertexInSlice = TotalVisibleEnergyFromVertex = VisibleEnergyFromVVerticesInSlice = -999; Reaction = ""; IsCC = false; @@ -1640,6 +1652,8 @@ void TMS_TreeWriter::Clear() { RecoTrackPrimaryParticlePDG[i] = -999; RecoTrackSecondaryParticlePDG[i] = -999; + RecoTrackPrimaryParticleIsPrimary[i] = false; + RecoTrackSecondaryParticleIsPrimary[i] = false; RecoTrackPrimaryParticleTrueTrackLength[i] = -999; RecoTrackPrimaryParticleTrueTrackLengthIgnoreY[i] = -999; diff --git a/src/TMS_TreeWriter.h b/src/TMS_TreeWriter.h index 853facb9..51a86e34 100644 --- a/src/TMS_TreeWriter.h +++ b/src/TMS_TreeWriter.h @@ -255,11 +255,14 @@ class TMS_TreeWriter { bool IsCC; int nTrueParticles; + int nTruePrimaryParticles; + int nTrueForgottenParticles; int VertexID[__TMS_MAX_TRUE_PARTICLES__]; int Parent[__TMS_MAX_TRUE_PARTICLES__]; int TrackId[__TMS_MAX_TRUE_PARTICLES__]; int PDG[__TMS_MAX_TRUE_PARTICLES__]; + bool IsPrimary[__TMS_MAX_TRUE_PARTICLES__]; float TrueVisibleEnergy[__TMS_MAX_TRUE_PARTICLES__]; float TruePathLength[__TMS_MAX_TRUE_PARTICLES__]; float TruePathLengthIgnoreY[__TMS_MAX_TRUE_PARTICLES__]; @@ -322,6 +325,7 @@ class TMS_TreeWriter { // Save truth info about primary particle int RecoTrackPrimaryParticlePDG[__TMS_MAX_LINES__]; + int RecoTrackPrimaryParticleIsPrimary[__TMS_MAX_LINES__]; float RecoTrackPrimaryParticleTrueMomentum[__TMS_MAX_LINES__][4]; float RecoTrackPrimaryParticleTruePositionStart[__TMS_MAX_LINES__][4]; float RecoTrackPrimaryParticleTruePositionEnd[__TMS_MAX_LINES__][4]; @@ -343,6 +347,7 @@ class TMS_TreeWriter { // Save truth info about secondary particle int RecoTrackSecondaryParticlePDG[__TMS_MAX_LINES__]; + int RecoTrackSecondaryParticleIsPrimary[__TMS_MAX_LINES__]; float RecoTrackSecondaryParticleTrueMomentum[__TMS_MAX_LINES__][4]; float RecoTrackSecondaryParticleTruePositionStart[__TMS_MAX_LINES__][4]; float RecoTrackSecondaryParticleTruePositionEnd[__TMS_MAX_LINES__][4]; diff --git a/src/TMS_Utils.cpp b/src/TMS_Utils.cpp index c5748f83..9536059a 100644 --- a/src/TMS_Utils.cpp +++ b/src/TMS_Utils.cpp @@ -144,17 +144,17 @@ namespace TMS_Utils { std::cout<<"Starting with n hits: "< indices; + std::vector trackids; std::vector energies; + + void Print() { + std::cout<<"ParticleInfo.total_energy = "< Date: Tue, 28 May 2024 20:00:04 -0500 Subject: [PATCH 21/30] Add back setup.sh for SL7 users. Also adjust Makefiles --- app/Makefile | 15 ++++++++++++--- setup.sh | 51 +++++++++++++++++++++++++++++++++++++++++++++++++++ src/Makefile | 15 ++++++++++++--- 3 files changed, 75 insertions(+), 6 deletions(-) create mode 100644 setup.sh diff --git a/app/Makefile b/app/Makefile index 09eaecff..c5e47966 100755 --- a/app/Makefile +++ b/app/Makefile @@ -6,9 +6,6 @@ LD_SHARED=g++ ROOT_LIBS = $(shell root-config --evelibs) ROOT_INC = $(shell root-config --cflags) -GEANT4_LIBS = $(shell geant4-config --libs) -GEANT4_INC = $(shell geant4-config --prefix)/include/ - ifndef CLHEP_INC # If not defined, normal build. CLHEP_INC = $(shell clhep-config --prefix)/include CLHEP_LIBS = $(shell clhep-config --libs) @@ -18,8 +15,20 @@ else # github CI passes CLHEP_INC as an environment var, use this to set lib dir endif endif +# Check if geant4-config is available +GEANT4_CONFIG := $(shell command -v geant4-config 2>/dev/null) + +# Set GEANT4 and EDEP variables based on the availability of geant4-config +ifneq ($(GEANT4_CONFIG),) +GEANT4_LIBS = $(shell geant4-config --libs) +GEANT4_INC = $(shell geant4-config --prefix)/include/ EDEP_LIBS = -L $(EDEPSIM_LIB) -ledepsim_io $(CLHEP_LIBS) EDEP_INC = -I$(EDEPSIM_INC) -I$(GEANT4_INC) -I$(CLHEP_INC) +else +EDEP_LIBS = -L $(EDEPSIM_LIB) -ledepsim_io +# Need CLHEP for unit conversion, can probably just move this into a standalone header to remove dependency? +EDEP_INC = -I$(EDEPSIM_INC) -I$(CLHEP_INC) +endif # The TMS includes TMS_INC = -I../src diff --git a/setup.sh b/setup.sh new file mode 100644 index 00000000..2b04396a --- /dev/null +++ b/setup.sh @@ -0,0 +1,51 @@ +#!/bin/bash + +# Example of setup needed for SSRI dependencies +# Here we use our own version of edep-sim, but you should be able to use one set up by ups + +MY_SSRI_DIR=/dune/app/users/cwret/SSRI + +# Copied from Gavin's SSRI setup +source /cvmfs/dune.opensciencegrid.org/products/dune/setup_dune.sh + +#setup dk2nu v01_05_01b -q e15:prof +#setup genie v2_12_10c -q e15:prof +#setup genie_xsec v2_12_10 -q DefaultPlusValenciaMEC +#setup genie_phyopt v2_12_10 -q dkcharmtau +#setup geant4 v4_10_3_p01b -q e15:prof +#setup ifdhc + +#setup root v6_22_08b -f Linux64bit+3.10-2.17 -q e20:p383b:prof +setup cmake v3_17_3 -f Linux64bit+3.10-2.17 + +# Official duneanaobj UPS products doesn't include TMS info yet... +#setup duneanaobj v01_00_00 -f Linux64bit+3.10-2.17 -q debug:e17:gv2 + +# Setup custom UPS products, for CAFana format development +#source /dune/app/users/cwret/ups/setup.sh +#setup duneanaobj v02_01_01 -q e20:gv3:prof +# Flag for enabling linking against DUNEANAOBJ +#if [ ! -z "${DUNEANAOBJ_INC}" ]; then + #export DUNEANAOBJ_ENABLED=1 + #echo "Enabling DUNEANAOBJ..." +#fi + +# Need CLHEP for unit conversion in edep-sim info +setup clhep v2_4_4_1 -f Linux64bit+3.10-2.17 -q debug:e20 +# Specific requirements for duneanaobj +setup root v6_12_06a -f Linux64bit+3.10-2.17 -q e15:prof + +# Set up edep-sim +setup edepsim v3_2_0 -f Linux64bit+3.10-2.17 -q e20:prof + +# Add TMS execs and library directory to env +export TMS_DIR=${PWD} +export PATH=${PATH}:${TMS_DIR}/bin +export LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${TMS_DIR}/lib + +# Debugging options passed to makefile +#export DEBUG=1 +#export VERBOSE=1 + + +echo "Setup TMS environment" diff --git a/src/Makefile b/src/Makefile index f36ad8d7..f6d97aef 100644 --- a/src/Makefile +++ b/src/Makefile @@ -6,9 +6,6 @@ LD_SHARED=g++ ROOT_LIBS = $(shell root-config --evelibs) ROOT_INC = $(shell root-config --cflags) -GEANT4_LIBS = $(shell geant4-config --libs) -GEANT4_INC = $(shell geant4-config --prefix)/include/ - ifndef CLHEP_INC # If not defined, normal build. CLHEP_INC = $(shell clhep-config --prefix)/include CLHEP_LIBS = $(shell clhep-config --libs) @@ -16,8 +13,20 @@ else # github CI passes CLHEP_INC as an environment var, use this to set lib dir CLHEP_LIBS = -L$(CLHEP_INC)/../../lib64 endif +# Check if geant4-config is available +GEANT4_CONFIG := $(shell command -v geant4-config 2>/dev/null) + +# Set GEANT4 and EDEP variables based on the availability of geant4-config +ifneq ($(GEANT4_CONFIG),) +GEANT4_LIBS = $(shell geant4-config --libs) +GEANT4_INC = $(shell geant4-config --prefix)/include/ EDEP_LIBS = -L $(EDEPSIM_LIB) -ledepsim_io $(CLHEP_LIBS) EDEP_INC = -I$(EDEPSIM_INC) -I$(GEANT4_INC) -I$(CLHEP_INC) +else +EDEP_LIBS = -L $(EDEPSIM_LIB) -ledepsim_io +# Need CLHEP for unit conversion, can probably just move this into a standalone header to remove dependency? +EDEP_INC = -I$(EDEPSIM_INC) -I$(CLHEP_INC) +endif TOML_INC = -I../toml11/ From a481df314108b559570dbcc4994506a57614bffa Mon Sep 17 00:00:00 2001 From: Jeffrey Kleykamp Date: Tue, 28 May 2024 21:15:16 -0500 Subject: [PATCH 22/30] Fix typo in TMS_Geom --- src/TMS_Geom.h | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/TMS_Geom.h b/src/TMS_Geom.h index a4a084ef..249eecf7 100644 --- a/src/TMS_Geom.h +++ b/src/TMS_Geom.h @@ -60,13 +60,13 @@ class TMS_Geom { inline TVector3 GetEndOfTMSFirstTwoModules() const { return TVector3(GetXEndOfTMS(), GetYEndOfTMS(), GetZEndOfTMSFirstTwoModules()); }; inline TVector3 GetStartOfTMSFiducial() const { return TVector3(TMS_Manager::GetInstance().Get_FIDUCIAL_TMS_START_X(), - TMS_Manager::GetInstance().Get_FIDUCIAL_TMS_START_Y(), TMS_Manager::GetInstance().Get_FIDUCIAL_TMS_START_Y()); }; + TMS_Manager::GetInstance().Get_FIDUCIAL_TMS_START_Y(), TMS_Manager::GetInstance().Get_FIDUCIAL_TMS_START_Z()); }; inline TVector3 GetEndOfTMSFiducial() const { return TVector3(TMS_Manager::GetInstance().Get_FIDUCIAL_TMS_END_X(), - TMS_Manager::GetInstance().Get_FIDUCIAL_TMS_END_Y(), TMS_Manager::GetInstance().Get_FIDUCIAL_TMS_END_Y()); }; + TMS_Manager::GetInstance().Get_FIDUCIAL_TMS_END_Y(), TMS_Manager::GetInstance().Get_FIDUCIAL_TMS_END_Z()); }; inline TVector3 GetStartOfLArFiducial() const { return TVector3(TMS_Manager::GetInstance().Get_FIDUCIAL_LAR_START_X(), - TMS_Manager::GetInstance().Get_FIDUCIAL_LAR_START_Y(), TMS_Manager::GetInstance().Get_FIDUCIAL_LAR_START_Y()); }; + TMS_Manager::GetInstance().Get_FIDUCIAL_LAR_START_Y(), TMS_Manager::GetInstance().Get_FIDUCIAL_LAR_START_Z()); }; inline TVector3 GetEndOfLArFiducial() const { return TVector3(TMS_Manager::GetInstance().Get_FIDUCIAL_LAR_END_X(), - TMS_Manager::GetInstance().Get_FIDUCIAL_LAR_END_Y(), TMS_Manager::GetInstance().Get_FIDUCIAL_LAR_END_Y()); }; + TMS_Manager::GetInstance().Get_FIDUCIAL_LAR_END_Y(), TMS_Manager::GetInstance().Get_FIDUCIAL_LAR_END_Z()); }; bool IsInsideBox(TVector3 position, TVector3 start, TVector3 end) const { if (position.X() < start.X()) return false; From 089e41ab3caee3b387137c5bd635bd44b86f4868 Mon Sep 17 00:00:00 2001 From: Jeffrey Kleykamp Date: Tue, 28 May 2024 22:33:14 -0500 Subject: [PATCH 23/30] Add all vars to clearing function --- src/TMS_TreeWriter.cpp | 383 +++++++++++++++++++++++------------------ src/TMS_TreeWriter.h | 12 +- 2 files changed, 217 insertions(+), 178 deletions(-) diff --git a/src/TMS_TreeWriter.cpp b/src/TMS_TreeWriter.cpp index 900fbcf3..8f8259d3 100644 --- a/src/TMS_TreeWriter.cpp +++ b/src/TMS_TreeWriter.cpp @@ -243,7 +243,6 @@ void TMS_TreeWriter::MakeBranches() { // Truth information Truth_Info->Branch("nParticles", &nParticles, "nParticles/I"); - Truth_Info->Branch("NeutrinoPDG", &NeutrinoPDG, "NeutrinoPDG/I"); Truth_Info->Branch("LeptonPDG", &LeptonPDG, "LeptonPDG/I"); Truth_Info->Branch("LeptonP4", LeptonP4, "LeptonP4[4]/F"); Truth_Info->Branch("LeptonX4", LeptonX4, "LeptonX4[4]/F"); @@ -1458,239 +1457,279 @@ void TMS_TreeWriter::FillSpill(TMS_Event &event, int truth_info_entry_number, in // Reset the variables void TMS_TreeWriter::Clear() { + const float DEFAULT_CLEARING_FLOAT = -999999999; + // Reset truth information - EventNo = nParticles = nTrueParticles = nTruePrimaryParticles = NeutrinoPDG = LeptonPDG = Muon_TrueKE = Muon_TrueTrackLength = VertexIdOfMostEnergyInEvent = -999; - VertexIdOfMostEnergyInEvent = VisibleEnergyFromUVertexInSlice = TotalVisibleEnergyFromVertex = VisibleEnergyFromVVerticesInSlice = -999; + EventNo = nParticles = nTrueParticles = nTruePrimaryParticles = NeutrinoPDG = LeptonPDG = Muon_TrueKE = Muon_TrueTrackLength = VertexIdOfMostEnergyInEvent = DEFAULT_CLEARING_FLOAT; + VertexIdOfMostEnergyInEvent = VisibleEnergyFromUVertexInSlice = TotalVisibleEnergyFromVertex = VisibleEnergyFromVVerticesInSlice = DEFAULT_CLEARING_FLOAT; Reaction = ""; IsCC = false; for (int i = 0; i < 4; ++i) { - MuonP4[i]=-999; - Muon_Vertex[i]=-999; - Muon_Death[i]=-999; - NeutrinoP4[i]=-999; - NeutrinoX4[i]=-999; - LeptonP4[i]=-999; - LeptonX4[i]=-999; + MuonP4[i]=DEFAULT_CLEARING_FLOAT; + Muon_Vertex[i]=DEFAULT_CLEARING_FLOAT; + Muon_Death[i]=DEFAULT_CLEARING_FLOAT; + NeutrinoP4[i]=DEFAULT_CLEARING_FLOAT; + NeutrinoX4[i]=DEFAULT_CLEARING_FLOAT; + LeptonP4[i]=DEFAULT_CLEARING_FLOAT; + LeptonX4[i]=DEFAULT_CLEARING_FLOAT; } // Reset line information TMSStart = false; - nLinesU = -999; - nLinesV = -999; - nLinesX = -999; + nLinesU = DEFAULT_CLEARING_FLOAT; + nLinesV = DEFAULT_CLEARING_FLOAT; + nLinesX = DEFAULT_CLEARING_FLOAT; for (int i = 0; i < __TMS_MAX_LINES__; ++i) { - SlopeU[i] = -999; - SlopeV[i] = -999; - SlopeX[i] = -999; - InterceptU[i] = -999; - InterceptV[i] = -999; - InterceptX[i] = -999; - Slope_DownstreamU[i] = -999; - Slope_DownstreamV[i] = -999; - Slope_DownstreamX[i] = -999; - Intercept_DownstreamU[i] = -999; - Intercept_DownstreamV[i] = -999; - Intercept_DownstreamX[i] = -999; - Slope_UpstreamU[i] = -999; - Slope_UpstreamV[i] = -999; - Slope_UpstreamX[i] = -999; - Intercept_UpstreamU[i] = -999; - Intercept_UpstreamV[i] = -999; - Intercept_UpstreamX[i] = -999; - - DirectionZU[i] = -999; - DirectionXU[i] = -999; - DirectionZU_Upstream[i] = -999; - DirectionXU_Upstream[i] = -999; - DirectionZU_Downstream[i] = -999; - DirectionXU_Downstream[i] = -999; - - DirectionZV[i] = -999; - DirectionXV[i] = -999; - DirectionZV_Upstream[i] = -999; - DirectionXV_Upstream[i] = -999; - DirectionZV_Downstream[i] = -999; - DirectionXV_Downstream[i] = -999; - - DirectionZX[i] = -999; - DirectionYX[i] = -999; - DirectionZX_Upstream[i] = -999; - DirectionYX_Upstream[i] = -999; - DirectionZX_Downstream[i] = -999; - DirectionYX_Downstream[i] = -999; - - OccupancyU[i] = -999; - OccupancyV[i] = -999; - OccupancyX[i] = -999; - TrackLengthU[i] = -999; - TrackLengthV[i] = -999; - TrackLengthX[i] = -999; - TotalTrackEnergyU[i] = -999; - TotalTrackEnergyV[i] = -999; - TotalTrackEnergyX[i] = -999; - FirstPlaneU[i] = -999; - FirstPlaneV[i] = -999; - FirstPlaneX[i] = -999; - LastPlaneU[i] = -999; - LastPlaneV[i] = -999; - LastPlaneX[i] = -999; - nHitsInTrackU[i] = -999; - nHitsInTrackV[i] = -999; - nHitsInTrackX[i] = -999; + SlopeU[i] = DEFAULT_CLEARING_FLOAT; + SlopeV[i] = DEFAULT_CLEARING_FLOAT; + SlopeX[i] = DEFAULT_CLEARING_FLOAT; + InterceptU[i] = DEFAULT_CLEARING_FLOAT; + InterceptV[i] = DEFAULT_CLEARING_FLOAT; + InterceptX[i] = DEFAULT_CLEARING_FLOAT; + Slope_DownstreamU[i] = DEFAULT_CLEARING_FLOAT; + Slope_DownstreamV[i] = DEFAULT_CLEARING_FLOAT; + Slope_DownstreamX[i] = DEFAULT_CLEARING_FLOAT; + Intercept_DownstreamU[i] = DEFAULT_CLEARING_FLOAT; + Intercept_DownstreamV[i] = DEFAULT_CLEARING_FLOAT; + Intercept_DownstreamX[i] = DEFAULT_CLEARING_FLOAT; + Slope_UpstreamU[i] = DEFAULT_CLEARING_FLOAT; + Slope_UpstreamV[i] = DEFAULT_CLEARING_FLOAT; + Slope_UpstreamX[i] = DEFAULT_CLEARING_FLOAT; + Intercept_UpstreamU[i] = DEFAULT_CLEARING_FLOAT; + Intercept_UpstreamV[i] = DEFAULT_CLEARING_FLOAT; + Intercept_UpstreamX[i] = DEFAULT_CLEARING_FLOAT; + + DirectionZU[i] = DEFAULT_CLEARING_FLOAT; + DirectionXU[i] = DEFAULT_CLEARING_FLOAT; + DirectionZU_Upstream[i] = DEFAULT_CLEARING_FLOAT; + DirectionXU_Upstream[i] = DEFAULT_CLEARING_FLOAT; + DirectionZU_Downstream[i] = DEFAULT_CLEARING_FLOAT; + DirectionXU_Downstream[i] = DEFAULT_CLEARING_FLOAT; + + DirectionZV[i] = DEFAULT_CLEARING_FLOAT; + DirectionXV[i] = DEFAULT_CLEARING_FLOAT; + DirectionZV_Upstream[i] = DEFAULT_CLEARING_FLOAT; + DirectionXV_Upstream[i] = DEFAULT_CLEARING_FLOAT; + DirectionZV_Downstream[i] = DEFAULT_CLEARING_FLOAT; + DirectionXV_Downstream[i] = DEFAULT_CLEARING_FLOAT; + + DirectionZX[i] = DEFAULT_CLEARING_FLOAT; + DirectionYX[i] = DEFAULT_CLEARING_FLOAT; + DirectionZX_Upstream[i] = DEFAULT_CLEARING_FLOAT; + DirectionYX_Upstream[i] = DEFAULT_CLEARING_FLOAT; + DirectionZX_Downstream[i] = DEFAULT_CLEARING_FLOAT; + DirectionYX_Downstream[i] = DEFAULT_CLEARING_FLOAT; + + OccupancyU[i] = DEFAULT_CLEARING_FLOAT; + OccupancyV[i] = DEFAULT_CLEARING_FLOAT; + OccupancyX[i] = DEFAULT_CLEARING_FLOAT; + TrackLengthU[i] = DEFAULT_CLEARING_FLOAT; + TrackLengthV[i] = DEFAULT_CLEARING_FLOAT; + TrackLengthX[i] = DEFAULT_CLEARING_FLOAT; + TotalTrackEnergyU[i] = DEFAULT_CLEARING_FLOAT; + TotalTrackEnergyV[i] = DEFAULT_CLEARING_FLOAT; + TotalTrackEnergyX[i] = DEFAULT_CLEARING_FLOAT; + FirstPlaneU[i] = DEFAULT_CLEARING_FLOAT; + FirstPlaneV[i] = DEFAULT_CLEARING_FLOAT; + FirstPlaneX[i] = DEFAULT_CLEARING_FLOAT; + LastPlaneU[i] = DEFAULT_CLEARING_FLOAT; + LastPlaneV[i] = DEFAULT_CLEARING_FLOAT; + LastPlaneX[i] = DEFAULT_CLEARING_FLOAT; + nHitsInTrackU[i] = DEFAULT_CLEARING_FLOAT; + nHitsInTrackV[i] = DEFAULT_CLEARING_FLOAT; + nHitsInTrackX[i] = DEFAULT_CLEARING_FLOAT; TrackStoppingU[i] = false; TrackStoppingV[i] = false; TrackStoppingX[i] = false; } -/* Occupancy3D[i] = -999; - TrackLength3D[i] = -999; - TotalTrackEnergy3D[i] = -999; - FirstPlane3D[i] = -999; - LastPlane3D[i] = -999; - nHitsInTrack3D[i] = -999; +/* Occupancy3D[i] = DEFAULT_CLEARING_FLOAT; + TrackLength3D[i] = DEFAULT_CLEARING_FLOAT; + TotalTrackEnergy3D[i] = DEFAULT_CLEARING_FLOAT; + FirstPlane3D[i] = DEFAULT_CLEARING_FLOAT; + LastPlane3D[i] = DEFAULT_CLEARING_FLOAT; + nHitsInTrack3D[i] = DEFAULT_CLEARING_FLOAT; TrackStopping3D[i] = false; for (int j = 0; j < 2; ++j) { - FirstHitOne[i][j] = -999; - FirstHitOther[i][j] = -999; - LastHitOne[i][j] = -999; - LastHitOther[i][j] = -999; + FirstHitOne[i][j] = DEFAULT_CLEARING_FLOAT; + FirstHitOther[i][j] = DEFAULT_CLEARING_FLOAT; + LastHitOne[i][j] = DEFAULT_CLEARING_FLOAT; + LastHitOther[i][j] = DEFAULT_CLEARING_FLOAT; } for (int j = 0; j < 3; ++j) { - FirstHit3D[i][j] = -999; - LastHit3D[i][j] = -999; + FirstHit3D[i][j] = DEFAULT_CLEARING_FLOAT; + LastHit3D[i][j] = DEFAULT_CLEARING_FLOAT; } for (int j = 0; j < __TMS_MAX_LINE_HITS__; ++j) { - TrackHitEnergyOne[i][j]=-999; - TrackHitEnergyOther[i][j]=-999; - TrackHitPosOne[i][j][0]=-999; - TrackHitPosOne[i][j][1]=-999; - TrackHitPosOther[i][j][0]=-999; - TrackHitPosOther[i][j][1]=-999; - - TrackHitEnergy3D[i][j] = -999; - TrackHitPos3D[i][j][0] = -999; - TrackHitPos3D[i][j][1] = -999; - TrackHitPos3D[i][j][2] = -999; + TrackHitEnergyOne[i][j]=DEFAULT_CLEARING_FLOAT; + TrackHitEnergyOther[i][j]=DEFAULT_CLEARING_FLOAT; + TrackHitPosOne[i][j][0]=DEFAULT_CLEARING_FLOAT; + TrackHitPosOne[i][j][1]=DEFAULT_CLEARING_FLOAT; + TrackHitPosOther[i][j][0]=DEFAULT_CLEARING_FLOAT; + TrackHitPosOther[i][j][1]=DEFAULT_CLEARING_FLOAT; + + TrackHitEnergy3D[i][j] = DEFAULT_CLEARING_FLOAT; + TrackHitPos3D[i][j][0] = DEFAULT_CLEARING_FLOAT; + TrackHitPos3D[i][j][1] = DEFAULT_CLEARING_FLOAT; + TrackHitPos3D[i][j][2] = DEFAULT_CLEARING_FLOAT; } }*/ // Reset hit information - nHits = -999; + nHits = DEFAULT_CLEARING_FLOAT; for (int i = 0; i < __TMS_MAX_HITS__; ++i) { - for (int j = 0; j < 4; ++j) RecoHitPos[i][j] = -999; - RecoHitEnergy[i] = -999; + for (int j = 0; j < 4; ++j) RecoHitPos[i][j] = DEFAULT_CLEARING_FLOAT; + RecoHitEnergy[i] = DEFAULT_CLEARING_FLOAT; } // Reset Cluster info - nClustersU = -999; - nClustersV = -999; - nClustersX = -999; + nClustersU = DEFAULT_CLEARING_FLOAT; + nClustersV = DEFAULT_CLEARING_FLOAT; + nClustersX = DEFAULT_CLEARING_FLOAT; for (int i = 0; i < __TMS_MAX_CLUSTERS__; ++i) { - ClusterEnergyU[i] = -999; - ClusterEnergyV[i] = -999; - ClusterEnergyX[i] = -999; - nHitsInClusterU[i] = -999; - nHitsInClusterV[i] = -999; - nHitsInClusterX[i] = -999; + ClusterEnergyU[i] = DEFAULT_CLEARING_FLOAT; + ClusterEnergyV[i] = DEFAULT_CLEARING_FLOAT; + ClusterEnergyX[i] = DEFAULT_CLEARING_FLOAT; + nHitsInClusterU[i] = DEFAULT_CLEARING_FLOAT; + nHitsInClusterV[i] = DEFAULT_CLEARING_FLOAT; + nHitsInClusterX[i] = DEFAULT_CLEARING_FLOAT; for (int j = 0; j < 2; ++j) { - ClusterPosMeanU[i][j] = -999; - ClusterPosStdDevU[i][j] = -999; + ClusterPosMeanU[i][j] = DEFAULT_CLEARING_FLOAT; + ClusterPosStdDevU[i][j] = DEFAULT_CLEARING_FLOAT; - ClusterPosMeanV[i][j] = -999; - ClusterPosStdDevV[i][j] = -999; + ClusterPosMeanV[i][j] = DEFAULT_CLEARING_FLOAT; + ClusterPosStdDevV[i][j] = DEFAULT_CLEARING_FLOAT; - ClusterPosMeanX[i][j] = -999; - ClusterPosStdDevX[i][j] = -999; + ClusterPosMeanX[i][j] = DEFAULT_CLEARING_FLOAT; + ClusterPosStdDevX[i][j] = DEFAULT_CLEARING_FLOAT; } for (int j = 0; j < __TMS_MAX_LINE_HITS__; ++j) { - ClusterHitPosU[i][j][0] = -999; - ClusterHitPosU[i][j][1] = -999; - ClusterHitEnergyU[i][j] = -999; + ClusterHitPosU[i][j][0] = DEFAULT_CLEARING_FLOAT; + ClusterHitPosU[i][j][1] = DEFAULT_CLEARING_FLOAT; + ClusterHitEnergyU[i][j] = DEFAULT_CLEARING_FLOAT; - ClusterHitPosV[i][j][0] = -999; - ClusterHitPosV[i][j][1] = -999; - ClusterHitEnergyV[i][j] = -999; + ClusterHitPosV[i][j][0] = DEFAULT_CLEARING_FLOAT; + ClusterHitPosV[i][j][1] = DEFAULT_CLEARING_FLOAT; + ClusterHitEnergyV[i][j] = DEFAULT_CLEARING_FLOAT; - ClusterHitPosX[i][j][0] = -999; - ClusterHitPosX[i][j][1] = -999; - ClusterHitEnergyX[i][j] = -999; + ClusterHitPosX[i][j][0] = DEFAULT_CLEARING_FLOAT; + ClusterHitPosX[i][j][1] = DEFAULT_CLEARING_FLOAT; + ClusterHitEnergyX[i][j] = DEFAULT_CLEARING_FLOAT; } } // Reset track information - nTracks = -999; + nTracks = DEFAULT_CLEARING_FLOAT; for (int i = 0; i < __TMS_MAX_TRACKS__; ++i) { for (int j = 0; j < 3; ++j) { - RecoTrackStartPos[i][j] = -999; - RecoTrackDirection[i][j] = -999; - RecoTrackEndPos[i][j] = -999; + RecoTrackStartPos[i][j] = DEFAULT_CLEARING_FLOAT; + RecoTrackDirection[i][j] = DEFAULT_CLEARING_FLOAT; + RecoTrackEndPos[i][j] = DEFAULT_CLEARING_FLOAT; } for (int k = 0; k < __TMS_MAX_LINE_HITS__; ++k) { - RecoTrackHitPos[i][k][0] = -999; - RecoTrackHitPos[i][k][1] = -999; - RecoTrackHitPos[i][k][2] = -999; + RecoTrackHitPos[i][k][0] = DEFAULT_CLEARING_FLOAT; + RecoTrackHitPos[i][k][1] = DEFAULT_CLEARING_FLOAT; + RecoTrackHitPos[i][k][2] = DEFAULT_CLEARING_FLOAT; } - RecoTrackEnergyRange[i] = -999; - RecoTrackEnergyDeposit[i] = -999; - RecoTrackLength[i] = -999; + RecoTrackEnergyRange[i] = DEFAULT_CLEARING_FLOAT; + RecoTrackEnergyDeposit[i] = DEFAULT_CLEARING_FLOAT; + RecoTrackLength[i] = DEFAULT_CLEARING_FLOAT; } RecoTrackN = 0; for (int i = 0; i < __TMS_MAX_LINES__; ++i) { - RecoTrackTrueVisibleEnergy[i] = -999; - RecoTrackPrimaryParticleIndex[i] = -999; - RecoTrackPrimaryParticleTrueVisibleEnergy[i] = -999; - RecoTrackSecondaryParticleIndex[i] = -999; - RecoTrackSecondaryParticleTrueVisibleEnergy[i] = -999; - - RecoTrackPrimaryParticleTrueTrackLengthAsMeasured[i] = -999; - RecoTrackPrimaryParticleTrueTrackLengthAsMeasuredIgnoreY[i] = -999; - RecoTrackPrimaryParticleTrueTrackLengthRecoStart[i] = -999; - RecoTrackPrimaryParticleTrueTrackLengthRecoStartIgnoreY[i] = -999; - RecoTrackPrimaryParticleTrueTrackLengthInTMS[i] = -999; - RecoTrackPrimaryParticleTrueTrackLengthInTMSIgnoreY[i] = -999; + RecoTrackTrueVisibleEnergy[i] = DEFAULT_CLEARING_FLOAT; + RecoTrackPrimaryParticleIndex[i] = DEFAULT_CLEARING_FLOAT; + RecoTrackPrimaryParticleTrueVisibleEnergy[i] = DEFAULT_CLEARING_FLOAT; + RecoTrackSecondaryParticleIndex[i] = DEFAULT_CLEARING_FLOAT; + RecoTrackSecondaryParticleTrueVisibleEnergy[i] = DEFAULT_CLEARING_FLOAT; + + RecoTrackPrimaryParticleTrueTrackLengthAsMeasured[i] = DEFAULT_CLEARING_FLOAT; + RecoTrackPrimaryParticleTrueTrackLengthAsMeasuredIgnoreY[i] = DEFAULT_CLEARING_FLOAT; + RecoTrackPrimaryParticleTrueTrackLengthRecoStart[i] = DEFAULT_CLEARING_FLOAT; + RecoTrackPrimaryParticleTrueTrackLengthRecoStartIgnoreY[i] = DEFAULT_CLEARING_FLOAT; + RecoTrackPrimaryParticleTrueTrackLengthInTMS[i] = DEFAULT_CLEARING_FLOAT; + RecoTrackPrimaryParticleTrueTrackLengthInTMSIgnoreY[i] = DEFAULT_CLEARING_FLOAT; - RecoTrackPrimaryParticlePDG[i] = -999; - RecoTrackSecondaryParticlePDG[i] = -999; + RecoTrackPrimaryParticlePDG[i] = DEFAULT_CLEARING_FLOAT; + RecoTrackSecondaryParticlePDG[i] = DEFAULT_CLEARING_FLOAT; RecoTrackPrimaryParticleIsPrimary[i] = false; RecoTrackSecondaryParticleIsPrimary[i] = false; - RecoTrackPrimaryParticleTrueTrackLength[i] = -999; - RecoTrackPrimaryParticleTrueTrackLengthIgnoreY[i] = -999; + RecoTrackPrimaryParticleTrueTrackLength[i] = DEFAULT_CLEARING_FLOAT; + RecoTrackPrimaryParticleTrueTrackLengthIgnoreY[i] = DEFAULT_CLEARING_FLOAT; + + RecoTrackPrimaryParticleTMSFiducialStart[i] = false; + RecoTrackPrimaryParticleTMSFiducialTouch[i] = false; + RecoTrackPrimaryParticleTMSFiducialEnd[i] = false; + RecoTrackPrimaryParticleLArFiducialStart[i] = false; + RecoTrackPrimaryParticleLArFiducialTouch[i] = false; + RecoTrackPrimaryParticleLArFiducialEnd[i] = false; for (int j = 0; j < 4; ++j) { - RecoTrackPrimaryParticleTrueMomentumTrackStart[i][j] = -999; - RecoTrackPrimaryParticleTruePositionTrackStart[i][j] = -999; - RecoTrackPrimaryParticleTrueMomentumTrackEnd[i][j] = -999; - RecoTrackPrimaryParticleTruePositionTrackEnd[i][j] = -999; + RecoTrackPrimaryParticleTrueMomentumTrackStart[i][j] = DEFAULT_CLEARING_FLOAT; + RecoTrackPrimaryParticleTruePositionTrackStart[i][j] = DEFAULT_CLEARING_FLOAT; + RecoTrackPrimaryParticleTrueMomentumTrackEnd[i][j] = DEFAULT_CLEARING_FLOAT; + RecoTrackPrimaryParticleTruePositionTrackEnd[i][j] = DEFAULT_CLEARING_FLOAT; - RecoTrackPrimaryParticleTrueMomentum[i][j] = -999; - RecoTrackPrimaryParticleTruePositionStart[i][j] = -999; - RecoTrackPrimaryParticleTruePositionEnd[i][j] = -999; - RecoTrackPrimaryParticleTrueMomentumEnteringTMS[i][j] = -999; - RecoTrackPrimaryParticleTruePositionEnteringTMS[i][j] = -999; - RecoTrackPrimaryParticleTrueMomentumLeavingTMS[i][j] = -999; - RecoTrackPrimaryParticleTruePositionLeavingTMS[i][j] = -999; - RecoTrackPrimaryParticleTrueMomentumLeavingLAr[i][j] = -999; - RecoTrackPrimaryParticleTruePositionLeavingLAr[i][j] = -999; + RecoTrackPrimaryParticleTrueMomentum[i][j] = DEFAULT_CLEARING_FLOAT; + RecoTrackPrimaryParticleTruePositionStart[i][j] = DEFAULT_CLEARING_FLOAT; + RecoTrackPrimaryParticleTruePositionEnd[i][j] = DEFAULT_CLEARING_FLOAT; + RecoTrackPrimaryParticleTrueMomentumEnteringTMS[i][j] = DEFAULT_CLEARING_FLOAT; + RecoTrackPrimaryParticleTruePositionEnteringTMS[i][j] = DEFAULT_CLEARING_FLOAT; + RecoTrackPrimaryParticleTrueMomentumLeavingTMS[i][j] = DEFAULT_CLEARING_FLOAT; + RecoTrackPrimaryParticleTruePositionLeavingTMS[i][j] = DEFAULT_CLEARING_FLOAT; + RecoTrackPrimaryParticleTrueMomentumLeavingLAr[i][j] = DEFAULT_CLEARING_FLOAT; + RecoTrackPrimaryParticleTruePositionLeavingLAr[i][j] = DEFAULT_CLEARING_FLOAT; - RecoTrackSecondaryParticleTrueMomentum[i][j] = -999; - RecoTrackSecondaryParticleTruePositionStart[i][j] = -999; - RecoTrackSecondaryParticleTruePositionEnd[i][j] = -999; + RecoTrackSecondaryParticleTrueMomentum[i][j] = DEFAULT_CLEARING_FLOAT; + RecoTrackSecondaryParticleTruePositionStart[i][j] = DEFAULT_CLEARING_FLOAT; + RecoTrackSecondaryParticleTruePositionEnd[i][j] = DEFAULT_CLEARING_FLOAT; } } nTrueParticles = 0; for (int i = 0; i < __TMS_MAX_TRUE_PARTICLES__; ++i) { - VertexID[i] = -999; - Parent[i] = -999; - TrackId[i] = -999; - PDG[i] = -999; - TrueVisibleEnergy[i] = -999; + VertexID[i] = DEFAULT_CLEARING_FLOAT; + Parent[i] = DEFAULT_CLEARING_FLOAT; + TrackId[i] = DEFAULT_CLEARING_FLOAT; + PDG[i] = DEFAULT_CLEARING_FLOAT; + IsPrimary[i] = false; + TrueVisibleEnergy[i] = DEFAULT_CLEARING_FLOAT; + TruePathLength[i] = DEFAULT_CLEARING_FLOAT; + TruePathLengthIgnoreY[i] = DEFAULT_CLEARING_FLOAT; + TruePathLengthInTMS[i] = DEFAULT_CLEARING_FLOAT; + TruePathLengthInTMSIgnoreY[i] = DEFAULT_CLEARING_FLOAT; + + TMSFiducialStart[i] = false; + TMSFiducialTouch[i] = false; + TMSFiducialEnd[i] = false; + LArFiducialStart[i] = false; + LArFiducialTouch[i] = false; + LArFiducialEnd[i] = false; + for (int j = 0; j < 4; ++j) { - BirthMomentum[i][j] = -999; - BirthPosition[i][j] = -999; - DeathMomentum[i][j] = -999; - DeathPosition[i][j] = -999; + BirthMomentum[i][j] = DEFAULT_CLEARING_FLOAT; + BirthPosition[i][j] = DEFAULT_CLEARING_FLOAT; + DeathMomentum[i][j] = DEFAULT_CLEARING_FLOAT; + DeathPosition[i][j] = DEFAULT_CLEARING_FLOAT; + MomentumZIsLArEnd[i][j] = DEFAULT_CLEARING_FLOAT; + PositionZIsLArEnd[i][j] = DEFAULT_CLEARING_FLOAT; + MomentumLArStart[i][j] = DEFAULT_CLEARING_FLOAT; + PositionLArStart[i][j] = DEFAULT_CLEARING_FLOAT; + MomentumLArEnd[i][j] = DEFAULT_CLEARING_FLOAT; + PositionLArEnd[i][j] = DEFAULT_CLEARING_FLOAT; + MomentumZIsTMSStart[i][j] = DEFAULT_CLEARING_FLOAT; + PositionZIsTMSStart[i][j] = DEFAULT_CLEARING_FLOAT; + MomentumZIsTMSEnd[i][j] = DEFAULT_CLEARING_FLOAT; + PositionZIsTMSEnd[i][j] = DEFAULT_CLEARING_FLOAT; + MomentumTMSStart[i][j] = DEFAULT_CLEARING_FLOAT; + PositionTMSStart[i][j] = DEFAULT_CLEARING_FLOAT; + MomentumTMSFirstTwoModulesEnd[i][j] = DEFAULT_CLEARING_FLOAT; + PositionTMSFirstTwoModulesEnd[i][j] = DEFAULT_CLEARING_FLOAT; + MomentumTMSThinEnd[i][j] = DEFAULT_CLEARING_FLOAT; + PositionTMSThinEnd[i][j] = DEFAULT_CLEARING_FLOAT; + MomentumTMSEnd[i][j] = DEFAULT_CLEARING_FLOAT; + PositionTMSEnd[i][j] = DEFAULT_CLEARING_FLOAT; } } diff --git a/src/TMS_TreeWriter.h b/src/TMS_TreeWriter.h index 51a86e34..10c1e1d5 100644 --- a/src/TMS_TreeWriter.h +++ b/src/TMS_TreeWriter.h @@ -338,12 +338,12 @@ class TMS_TreeWriter { float RecoTrackPrimaryParticleTrueMomentumLeavingLAr[__TMS_MAX_LINES__][4]; float RecoTrackPrimaryParticleTruePositionLeavingLAr[__TMS_MAX_LINES__][4]; - bool RecoTrackPrimaryParticleTMSFiducialStart[__TMS_MAX_TRUE_PARTICLES__]; - bool RecoTrackPrimaryParticleTMSFiducialTouch[__TMS_MAX_TRUE_PARTICLES__]; - bool RecoTrackPrimaryParticleTMSFiducialEnd[__TMS_MAX_TRUE_PARTICLES__]; - bool RecoTrackPrimaryParticleLArFiducialStart[__TMS_MAX_TRUE_PARTICLES__]; - bool RecoTrackPrimaryParticleLArFiducialTouch[__TMS_MAX_TRUE_PARTICLES__]; - bool RecoTrackPrimaryParticleLArFiducialEnd[__TMS_MAX_TRUE_PARTICLES__]; + bool RecoTrackPrimaryParticleTMSFiducialStart[__TMS_MAX_LINES__]; + bool RecoTrackPrimaryParticleTMSFiducialTouch[__TMS_MAX_LINES__]; + bool RecoTrackPrimaryParticleTMSFiducialEnd[__TMS_MAX_LINES__]; + bool RecoTrackPrimaryParticleLArFiducialStart[__TMS_MAX_LINES__]; + bool RecoTrackPrimaryParticleLArFiducialTouch[__TMS_MAX_LINES__]; + bool RecoTrackPrimaryParticleLArFiducialEnd[__TMS_MAX_LINES__]; // Save truth info about secondary particle int RecoTrackSecondaryParticlePDG[__TMS_MAX_LINES__]; From c72cbc5e8065fb3de8d028ce91b14393e02b9d26 Mon Sep 17 00:00:00 2001 From: Jeffrey Kleykamp Date: Thu, 30 May 2024 11:06:32 -0500 Subject: [PATCH 24/30] Fix GetTrackLength bug by avoiding geometry bug --- src/TMS_Geom.h | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/src/TMS_Geom.h b/src/TMS_Geom.h index 249eecf7..05e1442c 100644 --- a/src/TMS_Geom.h +++ b/src/TMS_Geom.h @@ -353,6 +353,12 @@ class TMS_Geom { TVector3 point1 = Unscale(point1_temp); TVector3 point2 = Unscale(point2_temp); + // Root's geometry has a bug that if you call a point outside the world (or maybe inside sand?) + // then from then on we get no materials. Not sure what's going on, but we can avoid it + // by completely avoiding locations outside "reasonable size" + // Checked and GetMaterials returns a vector of size zero once the bug is triggered + if (!IsInsideReasonableSize(point1) || !IsInsideReasonableSize(point2)) return 0; + // First get the collection of materials between point1 and point2 std::vector > Materials = GetMaterials(point1, point2); @@ -400,6 +406,15 @@ class TMS_Geom { } out += GetTrackLength(p1, p2); } + double distance = (nodes.front() - nodes.back()).Mag(); + if (distance > 0.1 && out < 0.01) { + std::cout<<"Found track length < 0.01 with distance="< Date: Thu, 30 May 2024 21:53:44 -0500 Subject: [PATCH 25/30] Make default track length avg of U and V. Save reco track true hit info --- src/TMS_TreeWriter.cpp | 21 ++++++++++++++++++++- src/TMS_TreeWriter.h | 3 +++ 2 files changed, 23 insertions(+), 1 deletion(-) diff --git a/src/TMS_TreeWriter.cpp b/src/TMS_TreeWriter.cpp index 8f8259d3..db22f1d7 100644 --- a/src/TMS_TreeWriter.cpp +++ b/src/TMS_TreeWriter.cpp @@ -277,6 +277,10 @@ void TMS_TreeWriter::MakeBranches() { Truth_Info->Branch("RecoTrackPrimaryParticleTruePositionTrackEnd", RecoTrackPrimaryParticleTruePositionTrackEnd, "RecoTrackPrimaryParticleTruePositionTrackEnd[RecoTrackN][4]/F"); + Truth_Info->Branch("RecoTrackNHits", RecoTrackNHits, "RecoTrackNHits[RecoTrackN]/I"); + Truth_Info->Branch("RecoTrackTrueHitPosition", RecoTrackTrueHitPosition, + "RecoTrackTrueHitPosition[RecoTrackN][200][4]/F"); + Truth_Info->Branch("RecoTrackPrimaryParticleTrueTrackLengthAsMeasured", RecoTrackPrimaryParticleTrueTrackLengthAsMeasured, "RecoTrackPrimaryParticleTrueTrackLengthAsMeasured[RecoTrackN]/F"); @@ -1198,7 +1202,7 @@ void TMS_TreeWriter::Fill(TMS_Event &event) { nHitsIn3DTrack[itTrack] = (int) RecoTrack->Hits.size(); // Do we need to cast it? idk // std::cout << "TreeWriter number of hits: " << nHitsIn3DTrack[itTrack] << std::endl; RecoTrackEnergyRange[itTrack] = RecoTrack->EnergyRange; - RecoTrackLength[itTrack] = RecoTrack->Length; + RecoTrackLength[itTrack] = 0.5 * (TrackLengthU[itTrack] + TrackLengthV[itTrack]); // RecoTrack->Length;, 2d is better estimate than 3d because of y jumps RecoTrackEnergyDeposit[itTrack] = RecoTrack->EnergyDeposit; for (int j = 0; j < 3; j++) { RecoTrackStartPos[itTrack][j] = RecoTrack->Start[j]; @@ -1265,6 +1269,17 @@ void TMS_TreeWriter::Fill(TMS_Event &event) { setPosition(RecoTrackPrimaryParticleTruePositionTrackStart[itTrack], tp.GetPositionAtZ(start_z, max_z_distance)); setMomentum(RecoTrackPrimaryParticleTrueMomentumTrackEnd[itTrack], tp.GetMomentumAtZ(end_z, max_z_distance)); setPosition(RecoTrackPrimaryParticleTruePositionTrackEnd[itTrack], tp.GetPositionAtZ(end_z, max_z_distance)); + + // TODO needs fixing + //if ( RecoTrack->Hits.size() != RecoTrack->nHits) std::cout<<"N hits mismatch: "<< RecoTrack->Hits.size() << " vs "<< RecoTrack->nHits<Hits.size(); + for (size_t h = 0; h < RecoTrack->Hits.size() && h < __TMS_MAX_LINE_HITS__; h++) { + auto& hit = RecoTrack->Hits[h]; + RecoTrackTrueHitPosition[itTrack][h][0] = hit.GetTrueHit().GetX(); + RecoTrackTrueHitPosition[itTrack][h][1] = hit.GetTrueHit().GetY(); + RecoTrackTrueHitPosition[itTrack][h][2] = hit.GetTrueHit().GetZ(); + RecoTrackTrueHitPosition[itTrack][h][3] = hit.GetTrueHit().GetT(); + } // Now calulate the true track length from true start to true end RecoTrackPrimaryParticleTrueTrackLengthAsMeasured[itTrack] = @@ -1684,6 +1699,10 @@ void TMS_TreeWriter::Clear() { RecoTrackSecondaryParticleTrueMomentum[i][j] = DEFAULT_CLEARING_FLOAT; RecoTrackSecondaryParticleTruePositionStart[i][j] = DEFAULT_CLEARING_FLOAT; RecoTrackSecondaryParticleTruePositionEnd[i][j] = DEFAULT_CLEARING_FLOAT; + + for (int h = 0; h < __TMS_MAX_LINE_HITS__; h++) { + RecoTrackTrueHitPosition[i][h][j] = DEFAULT_CLEARING_FLOAT; + } } } diff --git a/src/TMS_TreeWriter.h b/src/TMS_TreeWriter.h index 10c1e1d5..769fb3f8 100644 --- a/src/TMS_TreeWriter.h +++ b/src/TMS_TreeWriter.h @@ -338,6 +338,9 @@ class TMS_TreeWriter { float RecoTrackPrimaryParticleTrueMomentumLeavingLAr[__TMS_MAX_LINES__][4]; float RecoTrackPrimaryParticleTruePositionLeavingLAr[__TMS_MAX_LINES__][4]; + int RecoTrackNHits[__TMS_MAX_LINES__]; + float RecoTrackTrueHitPosition[__TMS_MAX_LINES__][__TMS_MAX_LINE_HITS__][4]; + bool RecoTrackPrimaryParticleTMSFiducialStart[__TMS_MAX_LINES__]; bool RecoTrackPrimaryParticleTMSFiducialTouch[__TMS_MAX_LINES__]; bool RecoTrackPrimaryParticleTMSFiducialEnd[__TMS_MAX_LINES__]; From 55e5ec293a1872d583c96dad18b06d7a6d8a2eac Mon Sep 17 00:00:00 2001 From: Jeffrey Kleykamp Date: Thu, 30 May 2024 22:14:42 -0500 Subject: [PATCH 26/30] Add validation scripts --- scripts/Validation/Line_Candidates.h | 547 ++++++++++++++++++ scripts/Validation/Makefile | 24 + scripts/Validation/Plot_True_Track_Length.cpp | 313 ++++++++++ scripts/Validation/Reco_Tracks.cpp | 236 ++++++++ scripts/Validation/Reco_Tree.h | 203 +++++++ scripts/Validation/Truth_Info.h | 495 ++++++++++++++++ scripts/Validation/simply_draw_everything.py | 53 ++ 7 files changed, 1871 insertions(+) create mode 100644 scripts/Validation/Line_Candidates.h create mode 100644 scripts/Validation/Makefile create mode 100644 scripts/Validation/Plot_True_Track_Length.cpp create mode 100644 scripts/Validation/Reco_Tracks.cpp create mode 100644 scripts/Validation/Reco_Tree.h create mode 100644 scripts/Validation/Truth_Info.h create mode 100644 scripts/Validation/simply_draw_everything.py diff --git a/scripts/Validation/Line_Candidates.h b/scripts/Validation/Line_Candidates.h new file mode 100644 index 00000000..d3494b02 --- /dev/null +++ b/scripts/Validation/Line_Candidates.h @@ -0,0 +1,547 @@ +////////////////////////////////////////////////////////// +// This class has been automatically generated on +// Tue May 28 22:24:07 2024 by ROOT version 6.22/08 +// from TTree Line_Candidates/Line_Candidates +// found on file: neutrino.0_1698176146.edep_TMS_RecoCandidates_Hough_Cluster1.root +////////////////////////////////////////////////////////// + +#ifndef Line_Candidates_h +#define Line_Candidates_h + +#include +#include +#include + +// Header file for the classes stored in the TTree if any. + +class Line_Candidates { +public : + TTree *fChain; //!pointer to the analyzed TTree or TChain + Int_t fCurrent; //!current Tree number in a TChain + +// Fixed size dimensions of array or collections stored in the TTree if any. + + // Declaration of leaf types + Int_t EventNo; + Int_t SliceNo; + Int_t SpillNo; + Int_t nLinesU; + Int_t nLinesV; + Int_t nLinesX; + Int_t nLines3D; + Float_t SlopeU[4]; //[nLinesU] + Float_t InterceptU[4]; //[nLinesU] + Float_t Slope_DownstreamU[4]; //[nLinesU] + Float_t Intercept_DownstreamU[4]; //[nLinesU] + Float_t Slope_UpstreamU[4]; //[nLinesU] + Float_t Intercept_UpstreamU[4]; //[nLinesU] + Float_t SlopeV[5]; //[nLinesV] + Float_t InterceptV[5]; //[nLinesV] + Float_t Slope_DownstreamV[5]; //[nLinesV] + Float_t Intercept_DownstreamV[5]; //[nLinesV] + Float_t Slope_UpstreamV[5]; //[nLinesV] + Float_t Intercept_UpstreamV[5]; //[nLinesV] + Float_t SlopeX[1]; //[nLinesX] + Float_t InterceptX[1]; //[nLinesX] + Float_t Slope_DownstreamX[1]; //[nLinesX] + Float_t Intercept_DownstreamX[1]; //[nLinesX] + Float_t Slope_UpstreamX[1]; //[nLinesX] + Float_t Intercept_UpstreamX[1]; //[nLinesX] + Float_t DirectionZU[4]; //[nLinesU] + Float_t DirectionZV[5]; //[nLinesV] + Float_t DirectionZX[1]; //[nLinesX] + Float_t DirectionXU[4]; //[nLinesU] + Float_t DirectionXV[5]; //[nLinesV] + Float_t DirectionYX[1]; //[nLinesX] + Float_t DirectionZU_Downstream[4]; //[nLinesU] + Float_t DirectionXU_Downstream[4]; //[nLinesU] + Float_t DirectionZU_Upstream[4]; //[nLinesU] + Float_t DirectionXU_Upstream[4]; //[nLinesU] + Float_t DirectionZV_Downstream[5]; //[nLinesV] + Float_t DirectionXV_Downstream[5]; //[nLinesV] + Float_t DirectionZV_Upstream[5]; //[nLinesV] + Float_t DirectionXV_Upstream[5]; //[nLinesV] + Float_t DirectionZX_Downstream[1]; //[nLinesX] + Float_t DirectionYX_Downstream[1]; //[nLinesX] + Float_t DirectionZX_Upstream[1]; //[nLinesX] + Float_t DirectionYX_Upstream[1]; //[nLinesX] + Float_t FirstHoughHitU[4][2]; //[nLinesU] + Float_t FirstHoughHitV[5][2]; //[nLinesV] + Float_t FirstHoughHitX[1][2]; //[nLinesX] + Float_t LastHoughHitU[4][2]; //[nLinesU] + Float_t LastHoughHitV[5][2]; //[nLinesV] + Float_t LastHoughHitX[1][2]; //[nLinesX] + Float_t FirstHoughHitTimeU[4]; //[nLinesU] + Float_t FirstHoughHitTimeV[5]; //[nLinesV] + Float_t FirstHoughHitTimeX[1]; //[nLinesX] + Float_t LastHoughHitTimeU[4]; //[nLinesU] + Float_t LastHoughHitTimeV[5]; //[nLinesV] + Float_t LastHoughHitTImeX[1]; //[nLinesX] + Float_t HoughEarliestHitTimeU[4]; //[nLinesU] + Float_t HoughEarliestHitTimeV[5]; //[nLinesV] + Float_t HoughEarliestHitTimeX[1]; //[nLinesX] + Float_t HoughLatestHitTimeU[4]; //[nLinesU] + Float_t HoughLatestHitTimeV[5]; //[nLinesV] + Float_t HoughLatestHitTimeX[1]; //[nLinesX] + Int_t FirstHoughPlaneU[4]; //[nLinesU] + Int_t FirstHoughPlaneV[5]; //[nLinesV] + Int_t FirstHoughPlaneX[1]; //[nLinesX] + Int_t LastHoughPlaneU[4]; //[nLinesU] + Int_t LastHoughPlaneV[5]; //[nLinesV] + Int_t LastHoughPlaneX[1]; //[nLinesX] + Bool_t TMSStart; + Float_t TMSStartTime; + Float_t OccupancyU[4]; //[nLinesU] + Float_t OccupancyV[5]; //[nLinesV] + Float_t OccupancyX[1]; //[nLinesX] + Float_t TrackLengthU[4]; //[nLinesU] + Float_t TrackLengthV[5]; //[nLinesV] + Float_t TrackLengthX[1]; //[nLinesX] + Float_t TotalTrackEnergyU[4]; //[nLinesU] + Float_t TotalTrackEnergyV[5]; //[nLinesV] + Float_t TotalTrackEnergyX[1]; //[nLinesX] + Bool_t TrackStoppingU[4]; //[nLinesU] + Bool_t TrackStoppingV[5]; //[nLinesV] + Bool_t TrackStoppingX[1]; //[nLinesX] + Int_t nHitsInTrackU[4]; //[nLinesU] + Int_t nHitsInTrackV[5]; //[nLinesV] + Int_t nHitsInTrackX[1]; //[nLinesX] + Float_t TrackHitEnergyU[10][200]; + Float_t TrackHitEnergyV[10][200]; + Float_t TrackHitEnergyX[10][200]; + Float_t TrackHitPosU[10][200][2]; + Float_t TrackHitPosV[10][200][2]; + Float_t TrackHitPosX[10][200][2]; + Float_t TrackHitTimeU[10][200]; + Float_t TrackHitTimeV[10][200]; + Float_t TrackHitTimeX[10][200]; + Int_t nClustersU; + Int_t nClustersV; + Int_t nClusterX; + Float_t ClusterEnergyU[10]; //[nClustersU] + Float_t ClusterEnergyV[9]; //[nClustersV] + Float_t ClusterEnergyX[1]; //[nClustersX] + Float_t ClusterTimeU[10]; //[nClustersU] + Float_t ClusterTimeV[9]; //[nClustersV] + Float_t ClusterTimeX[1]; //[nClustersX] + Float_t ClusterPosMeanU[25][2]; + Float_t ClusterPosMeanV[25][2]; + Float_t ClusterPosMeanX[25][2]; + Float_t ClusterPosStdDevU[25][2]; + Float_t ClusterPosStdDevV[25][2]; + Float_t ClusterPosStdDevX[25][2]; + Int_t nHitsInClusterU[10]; //[nClustersU] + Int_t nHitsInClusterV[9]; //[nClustersV] + Int_t nHitsInClusterX[1]; //[nClustersX] + Float_t ClusterHitPosU[25][200][2]; + Float_t ClusterHitPosV[25][200][2]; + Float_t ClusterHitPosX[25][200][2]; + Float_t ClusterHitEnergyU[25][200]; + Float_t ClusterHitEnergyV[25][200]; + Float_t ClusterHitEnergyX[25][200]; + Float_t ClusterHitTimeU[25][200]; + Float_t ClusterHitTimeV[25][200]; + Float_t ClusterHitTimeX[25][200]; + Int_t ClusterHitSliceU[25][200]; + Int_t ClusterHitSliceV[25][200]; + Int_t ClusterHitSliceX[25][200]; + Int_t nHits; + Float_t RecoHitPos[201][4]; //[nHits] + Float_t RecoHitEnergy[201]; //[nHits] + Int_t RecoHitSlice[201]; //[nHits] + + // List of branches + TBranch *b_EventNo; //! + TBranch *b_SliceNo; //! + TBranch *b_SpillNo; //! + TBranch *b_nLinesU; //! + TBranch *b_nLinesV; //! + TBranch *b_nLinesX; //! + TBranch *b_nLines3D; //! + TBranch *b_SlopeU; //! + TBranch *b_InterceptU; //! + TBranch *b_Slope_DownstreamU; //! + TBranch *b_Intercept_DownstreamU; //! + TBranch *b_Slope_UpstreamU; //! + TBranch *b_Intercept_UpstreamU; //! + TBranch *b_SlopeV; //! + TBranch *b_InterceptV; //! + TBranch *b_Slope_DownstreamV; //! + TBranch *b_Intercept_DownstreamV; //! + TBranch *b_Slope_UpstreamV; //! + TBranch *b_Intercept_UpstreamV; //! + TBranch *b_SlopeX; //! + TBranch *b_InterceptX; //! + TBranch *b_Slope_DownstreamX; //! + TBranch *b_Intercept_DownstreamX; //! + TBranch *b_Slope_UpstreamX; //! + TBranch *b_Intercept_UpstreamX; //! + TBranch *b_DirectionZU; //! + TBranch *b_DirectionZV; //! + TBranch *b_DirectionZX; //! + TBranch *b_DirectionXU; //! + TBranch *b_DirectionXV; //! + TBranch *b_DirectionYX; //! + TBranch *b_DirectionZU_Downstream; //! + TBranch *b_DirectionXU_Downstream; //! + TBranch *b_DirectionZU_Upstream; //! + TBranch *b_DirectionXU_Upstream; //! + TBranch *b_DirectionZV_Downstream; //! + TBranch *b_DirectionXV_Downstream; //! + TBranch *b_DirectionZV_Upstream; //! + TBranch *b_DirectionXV_Upstream; //! + TBranch *b_DirectionZX_Downstream; //! + TBranch *b_DirectionYX_Downstream; //! + TBranch *b_DirectionZX_Upstream; //! + TBranch *b_DirectionYX_Upstream; //! + TBranch *b_FirstHoughHitU; //! + TBranch *b_FirstHoughHitV; //! + TBranch *b_FirstHoughHitX; //! + TBranch *b_LastHoughHitU; //! + TBranch *b_LastHoughHitV; //! + TBranch *b_LastHoughHitX; //! + TBranch *b_FirstHoughHitTimeU; //! + TBranch *b_FirstHoughHitTimeV; //! + TBranch *b_FirstHoughHitTimeX; //! + TBranch *b_LastHoughHitTimeU; //! + TBranch *b_LastHoughHitTimeV; //! + TBranch *b_LastHoughHitTImeX; //! + TBranch *b_HoughEarliestHitTimeU; //! + TBranch *b_HoughEarliestHitTimeV; //! + TBranch *b_HoughEarliestHitTimeX; //! + TBranch *b_HoughLatestHitTimeU; //! + TBranch *b_HoughLatestHitTimeV; //! + TBranch *b_HoughLatestHitTimeX; //! + TBranch *b_FirstHoughPlaneU; //! + TBranch *b_FirstHoughPlaneV; //! + TBranch *b_FirstHoughPlaneX; //! + TBranch *b_LastHoughPlaneU; //! + TBranch *b_LastHoughPlaneV; //! + TBranch *b_LastHoughPlaneX; //! + TBranch *b_TMSStart; //! + TBranch *b_TMSStartTime; //! + TBranch *b_OccupancyU; //! + TBranch *b_OccupancyV; //! + TBranch *b_OccupancyX; //! + TBranch *b_TrackLengthU; //! + TBranch *b_TrackLengthV; //! + TBranch *b_TrackLengthX; //! + TBranch *b_TotalTrackEnergyU; //! + TBranch *b_TotalTrackEnergyV; //! + TBranch *b_TotalTrackEnergyX; //! + TBranch *b_TrackStoppingU; //! + TBranch *b_TrackStoppingV; //! + TBranch *b_TrackStoppingX; //! + TBranch *b_nHitsInTrackU; //! + TBranch *b_nHitsInTrackV; //! + TBranch *b_nHitsInTrackX; //! + TBranch *b_TrackHitEnergyU; //! + TBranch *b_TrackHitEnergyV; //! + TBranch *b_TrackHitEnergyX; //! + TBranch *b_TrackHitPosU; //! + TBranch *b_TrackHitPosV; //! + TBranch *b_TrackHitPosX; //! + TBranch *b_TrackHitTimeU; //! + TBranch *b_TrackHitTimeV; //! + TBranch *b_TrackHitTimeX; //! + TBranch *b_nClustersU; //! + TBranch *b_nClustersV; //! + TBranch *b_nClustersX; //! + TBranch *b_ClusterEnergyU; //! + TBranch *b_ClusterEnergyV; //! + TBranch *b_ClusterEnergyX; //! + TBranch *b_ClusterTimeU; //! + TBranch *b_ClusterTimeV; //! + TBranch *b_ClusterTimeX; //! + TBranch *b_ClusterPosMeanU; //! + TBranch *b_ClusterPosMeanV; //! + TBranch *b_ClusterPosMeanX; //! + TBranch *b_ClusterPosStdDevU; //! + TBranch *b_ClusterPosStdDevV; //! + TBranch *b_ClusterPosStdDevX; //! + TBranch *b_nHitsInClusterU; //! + TBranch *b_nHitsInClusterV; //! + TBranch *b_nHitsInClusterX; //! + TBranch *b_ClusterHitPosU; //! + TBranch *b_ClusterHitPosV; //! + TBranch *b_ClusterHitPosX; //! + TBranch *b_ClusterHitEnergyU; //! + TBranch *b_ClusterHitEnergyV; //! + TBranch *b_ClusterHitEnergyX; //! + TBranch *b_ClusterHitTimeU; //! + TBranch *b_ClusterHitTimeV; //! + TBranch *b_ClusterHitTimeX; //! + TBranch *b_ClusterHitSliceU; //! + TBranch *b_ClusterHitSliceV; //! + TBranch *b_ClusterHitSliceX; //! + TBranch *b_nHits; //! + TBranch *b_RecoHitPos; //! + TBranch *b_RecoHitEnergy; //! + TBranch *b_RecoHitSlice; //! + + Line_Candidates(TTree *tree=0); + virtual ~Line_Candidates(); + virtual Int_t Cut(Long64_t entry); + virtual Int_t GetEntry(Long64_t entry); + virtual Long64_t LoadTree(Long64_t entry); + virtual void Init(TTree *tree); + virtual void Loop(); + virtual Bool_t Notify(); + virtual void Show(Long64_t entry = -1); +}; + +#endif + +#define Line_Candidates_cxx +#ifdef Line_Candidates_cxx +Line_Candidates::Line_Candidates(TTree *tree) : fChain(0) +{ +// if parameter tree is not specified (or zero), connect the file +// used to generate this class and read the Tree. + if (tree == 0) { + TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("neutrino.0_1698176146.edep_TMS_RecoCandidates_Hough_Cluster1.root"); + if (!f || !f->IsOpen()) { + f = new TFile("neutrino.0_1698176146.edep_TMS_RecoCandidates_Hough_Cluster1.root"); + } + f->GetObject("Line_Candidates",tree); + + } + Init(tree); +} + +Line_Candidates::~Line_Candidates() +{ + if (!fChain) return; +} + +Int_t Line_Candidates::GetEntry(Long64_t entry) +{ +// Read contents of entry. + if (!fChain) return 0; + return fChain->GetEntry(entry); +} +Long64_t Line_Candidates::LoadTree(Long64_t entry) +{ +// Set the environment to read one entry + if (!fChain) return -5; + Long64_t centry = fChain->LoadTree(entry); + if (centry < 0) return centry; + if (fChain->GetTreeNumber() != fCurrent) { + fCurrent = fChain->GetTreeNumber(); + Notify(); + } + return centry; +} + +void Line_Candidates::Init(TTree *tree) +{ + // The Init() function is called when the selector needs to initialize + // a new tree or chain. Typically here the branch addresses and branch + // pointers of the tree will be set. + // It is normally not necessary to make changes to the generated + // code, but the routine can be extended by the user if needed. + // Init() will be called many times when running on PROOF + // (once per file to be processed). + + // Set branch addresses and branch pointers + if (!tree) return; + fChain = tree; + fCurrent = -1; + fChain->SetMakeClass(1); + + fChain->SetBranchAddress("EventNo", &EventNo, &b_EventNo); + fChain->SetBranchAddress("SliceNo", &SliceNo, &b_SliceNo); + fChain->SetBranchAddress("SpillNo", &SpillNo, &b_SpillNo); + fChain->SetBranchAddress("nLinesU", &nLinesU, &b_nLinesU); + fChain->SetBranchAddress("nLinesV", &nLinesV, &b_nLinesV); + fChain->SetBranchAddress("nLinesX", &nLinesX, &b_nLinesX); + fChain->SetBranchAddress("nLines3D", &nLines3D, &b_nLines3D); + fChain->SetBranchAddress("SlopeU", SlopeU, &b_SlopeU); + fChain->SetBranchAddress("InterceptU", InterceptU, &b_InterceptU); + fChain->SetBranchAddress("Slope_DownstreamU", Slope_DownstreamU, &b_Slope_DownstreamU); + fChain->SetBranchAddress("Intercept_DownstreamU", Intercept_DownstreamU, &b_Intercept_DownstreamU); + fChain->SetBranchAddress("Slope_UpstreamU", Slope_UpstreamU, &b_Slope_UpstreamU); + fChain->SetBranchAddress("Intercept_UpstreamU", Intercept_UpstreamU, &b_Intercept_UpstreamU); + fChain->SetBranchAddress("SlopeV", SlopeV, &b_SlopeV); + fChain->SetBranchAddress("InterceptV", InterceptV, &b_InterceptV); + fChain->SetBranchAddress("Slope_DownstreamV", Slope_DownstreamV, &b_Slope_DownstreamV); + fChain->SetBranchAddress("Intercept_DownstreamV", Intercept_DownstreamV, &b_Intercept_DownstreamV); + fChain->SetBranchAddress("Slope_UpstreamV", Slope_UpstreamV, &b_Slope_UpstreamV); + fChain->SetBranchAddress("Intercept_UpstreamV", Intercept_UpstreamV, &b_Intercept_UpstreamV); + fChain->SetBranchAddress("SlopeX", &SlopeX, &b_SlopeX); + fChain->SetBranchAddress("InterceptX", &InterceptX, &b_InterceptX); + fChain->SetBranchAddress("Slope_DownstreamX", &Slope_DownstreamX, &b_Slope_DownstreamX); + fChain->SetBranchAddress("Intercept_DownstreamX", &Intercept_DownstreamX, &b_Intercept_DownstreamX); + fChain->SetBranchAddress("Slope_UpstreamX", &Slope_UpstreamX, &b_Slope_UpstreamX); + fChain->SetBranchAddress("Intercept_UpstreamX", &Intercept_UpstreamX, &b_Intercept_UpstreamX); + fChain->SetBranchAddress("DirectionZU", DirectionZU, &b_DirectionZU); + fChain->SetBranchAddress("DirectionZV", DirectionZV, &b_DirectionZV); + fChain->SetBranchAddress("DirectionZX", &DirectionZX, &b_DirectionZX); + fChain->SetBranchAddress("DirectionXU", DirectionXU, &b_DirectionXU); + fChain->SetBranchAddress("DirectionXV", DirectionXV, &b_DirectionXV); + fChain->SetBranchAddress("DirectionYX", &DirectionYX, &b_DirectionYX); + fChain->SetBranchAddress("DirectionZU_Downstream", DirectionZU_Downstream, &b_DirectionZU_Downstream); + fChain->SetBranchAddress("DirectionXU_Downstream", DirectionXU_Downstream, &b_DirectionXU_Downstream); + fChain->SetBranchAddress("DirectionZU_Upstream", DirectionZU_Upstream, &b_DirectionZU_Upstream); + fChain->SetBranchAddress("DirectionXU_Upstream", DirectionXU_Upstream, &b_DirectionXU_Upstream); + fChain->SetBranchAddress("DirectionZV_Downstream", DirectionZV_Downstream, &b_DirectionZV_Downstream); + fChain->SetBranchAddress("DirectionXV_Downstream", DirectionXV_Downstream, &b_DirectionXV_Downstream); + fChain->SetBranchAddress("DirectionZV_Upstream", DirectionZV_Upstream, &b_DirectionZV_Upstream); + fChain->SetBranchAddress("DirectionXV_Upstream", DirectionXV_Upstream, &b_DirectionXV_Upstream); + fChain->SetBranchAddress("DirectionZX_Downstream", &DirectionZX_Downstream, &b_DirectionZX_Downstream); + fChain->SetBranchAddress("DirectionYX_Downstream", &DirectionYX_Downstream, &b_DirectionYX_Downstream); + fChain->SetBranchAddress("DirectionZX_Upstream", &DirectionZX_Upstream, &b_DirectionZX_Upstream); + fChain->SetBranchAddress("DirectionYX_Upstream", &DirectionYX_Upstream, &b_DirectionYX_Upstream); + fChain->SetBranchAddress("FirstHoughHitU", FirstHoughHitU, &b_FirstHoughHitU); + fChain->SetBranchAddress("FirstHoughHitV", FirstHoughHitV, &b_FirstHoughHitV); + fChain->SetBranchAddress("FirstHoughHitX", &FirstHoughHitX, &b_FirstHoughHitX); + fChain->SetBranchAddress("LastHoughHitU", LastHoughHitU, &b_LastHoughHitU); + fChain->SetBranchAddress("LastHoughHitV", LastHoughHitV, &b_LastHoughHitV); + fChain->SetBranchAddress("LastHoughHitX", &LastHoughHitX, &b_LastHoughHitX); + fChain->SetBranchAddress("FirstHoughHitTimeU", FirstHoughHitTimeU, &b_FirstHoughHitTimeU); + fChain->SetBranchAddress("FirstHoughHitTimeV", FirstHoughHitTimeV, &b_FirstHoughHitTimeV); + fChain->SetBranchAddress("FirstHoughHitTimeX", &FirstHoughHitTimeX, &b_FirstHoughHitTimeX); + fChain->SetBranchAddress("LastHoughHitTimeU", LastHoughHitTimeU, &b_LastHoughHitTimeU); + fChain->SetBranchAddress("LastHoughHitTimeV", LastHoughHitTimeV, &b_LastHoughHitTimeV); + fChain->SetBranchAddress("LastHoughHitTImeX", &LastHoughHitTImeX, &b_LastHoughHitTImeX); + fChain->SetBranchAddress("HoughEarliestHitTimeU", HoughEarliestHitTimeU, &b_HoughEarliestHitTimeU); + fChain->SetBranchAddress("HoughEarliestHitTimeV", HoughEarliestHitTimeV, &b_HoughEarliestHitTimeV); + fChain->SetBranchAddress("HoughEarliestHitTimeX", &HoughEarliestHitTimeX, &b_HoughEarliestHitTimeX); + fChain->SetBranchAddress("HoughLatestHitTimeU", HoughLatestHitTimeU, &b_HoughLatestHitTimeU); + fChain->SetBranchAddress("HoughLatestHitTimeV", HoughLatestHitTimeV, &b_HoughLatestHitTimeV); + fChain->SetBranchAddress("HoughLatestHitTimeX", &HoughLatestHitTimeX, &b_HoughLatestHitTimeX); + fChain->SetBranchAddress("FirstHoughPlaneU", FirstHoughPlaneU, &b_FirstHoughPlaneU); + fChain->SetBranchAddress("FirstHoughPlaneV", FirstHoughPlaneV, &b_FirstHoughPlaneV); + fChain->SetBranchAddress("FirstHoughPlaneX", &FirstHoughPlaneX, &b_FirstHoughPlaneX); + fChain->SetBranchAddress("LastHoughPlaneU", LastHoughPlaneU, &b_LastHoughPlaneU); + fChain->SetBranchAddress("LastHoughPlaneV", LastHoughPlaneV, &b_LastHoughPlaneV); + fChain->SetBranchAddress("LastHoughPlaneX", &LastHoughPlaneX, &b_LastHoughPlaneX); + fChain->SetBranchAddress("TMSStart", &TMSStart, &b_TMSStart); + fChain->SetBranchAddress("TMSStartTime", &TMSStartTime, &b_TMSStartTime); + fChain->SetBranchAddress("OccupancyU", OccupancyU, &b_OccupancyU); + fChain->SetBranchAddress("OccupancyV", OccupancyV, &b_OccupancyV); + fChain->SetBranchAddress("OccupancyX", &OccupancyX, &b_OccupancyX); + fChain->SetBranchAddress("TrackLengthU", TrackLengthU, &b_TrackLengthU); + fChain->SetBranchAddress("TrackLengthV", TrackLengthV, &b_TrackLengthV); + fChain->SetBranchAddress("TrackLengthX", &TrackLengthX, &b_TrackLengthX); + fChain->SetBranchAddress("TotalTrackEnergyU", TotalTrackEnergyU, &b_TotalTrackEnergyU); + fChain->SetBranchAddress("TotalTrackEnergyV", TotalTrackEnergyV, &b_TotalTrackEnergyV); + fChain->SetBranchAddress("TotalTrackEnergyX", &TotalTrackEnergyX, &b_TotalTrackEnergyX); + fChain->SetBranchAddress("TrackStoppingU", TrackStoppingU, &b_TrackStoppingU); + fChain->SetBranchAddress("TrackStoppingV", TrackStoppingV, &b_TrackStoppingV); + fChain->SetBranchAddress("TrackStoppingX", &TrackStoppingX, &b_TrackStoppingX); + fChain->SetBranchAddress("nHitsInTrackU", nHitsInTrackU, &b_nHitsInTrackU); + fChain->SetBranchAddress("nHitsInTrackV", nHitsInTrackV, &b_nHitsInTrackV); + fChain->SetBranchAddress("nHitsInTrackX", &nHitsInTrackX, &b_nHitsInTrackX); + fChain->SetBranchAddress("TrackHitEnergyU", TrackHitEnergyU, &b_TrackHitEnergyU); + fChain->SetBranchAddress("TrackHitEnergyV", TrackHitEnergyV, &b_TrackHitEnergyV); + fChain->SetBranchAddress("TrackHitEnergyX", TrackHitEnergyX, &b_TrackHitEnergyX); + fChain->SetBranchAddress("TrackHitPosU", TrackHitPosU, &b_TrackHitPosU); + fChain->SetBranchAddress("TrackHitPosV", TrackHitPosV, &b_TrackHitPosV); + fChain->SetBranchAddress("TrackHitPosX", TrackHitPosX, &b_TrackHitPosX); + fChain->SetBranchAddress("TrackHitTimeU", TrackHitTimeU, &b_TrackHitTimeU); + fChain->SetBranchAddress("TrackHitTimeV", TrackHitTimeV, &b_TrackHitTimeV); + fChain->SetBranchAddress("TrackHitTimeX", TrackHitTimeX, &b_TrackHitTimeX); + fChain->SetBranchAddress("nClustersU", &nClustersU, &b_nClustersU); + fChain->SetBranchAddress("nClustersV", &nClustersV, &b_nClustersV); + fChain->SetBranchAddress("nClusterX", &nClusterX, &b_nClustersX); + fChain->SetBranchAddress("ClusterEnergyU", ClusterEnergyU, &b_ClusterEnergyU); + fChain->SetBranchAddress("ClusterEnergyV", ClusterEnergyV, &b_ClusterEnergyV); + fChain->SetBranchAddress("ClusterEnergyX", &ClusterEnergyX, &b_ClusterEnergyX); + fChain->SetBranchAddress("ClusterTimeU", ClusterTimeU, &b_ClusterTimeU); + fChain->SetBranchAddress("ClusterTimeV", ClusterTimeV, &b_ClusterTimeV); + fChain->SetBranchAddress("ClusterTimeX", &ClusterTimeX, &b_ClusterTimeX); + fChain->SetBranchAddress("ClusterPosMeanU", ClusterPosMeanU, &b_ClusterPosMeanU); + fChain->SetBranchAddress("ClusterPosMeanV", ClusterPosMeanV, &b_ClusterPosMeanV); + fChain->SetBranchAddress("ClusterPosMeanX", ClusterPosMeanX, &b_ClusterPosMeanX); + fChain->SetBranchAddress("ClusterPosStdDevU", ClusterPosStdDevU, &b_ClusterPosStdDevU); + fChain->SetBranchAddress("ClusterPosStdDevV", ClusterPosStdDevV, &b_ClusterPosStdDevV); + fChain->SetBranchAddress("ClusterPosStdDevX", ClusterPosStdDevX, &b_ClusterPosStdDevX); + fChain->SetBranchAddress("nHitsInClusterU", nHitsInClusterU, &b_nHitsInClusterU); + fChain->SetBranchAddress("nHitsInClusterV", nHitsInClusterV, &b_nHitsInClusterV); + fChain->SetBranchAddress("nHitsInClusterX", &nHitsInClusterX, &b_nHitsInClusterX); + fChain->SetBranchAddress("ClusterHitPosU", ClusterHitPosU, &b_ClusterHitPosU); + fChain->SetBranchAddress("ClusterHitPosV", ClusterHitPosV, &b_ClusterHitPosV); + fChain->SetBranchAddress("ClusterHitPosX", ClusterHitPosX, &b_ClusterHitPosX); + fChain->SetBranchAddress("ClusterHitEnergyU", ClusterHitEnergyU, &b_ClusterHitEnergyU); + fChain->SetBranchAddress("ClusterHitEnergyV", ClusterHitEnergyV, &b_ClusterHitEnergyV); + fChain->SetBranchAddress("ClusterHitEnergyX", ClusterHitEnergyX, &b_ClusterHitEnergyX); + fChain->SetBranchAddress("ClusterHitTimeU", ClusterHitTimeU, &b_ClusterHitTimeU); + fChain->SetBranchAddress("ClusterHitTimeV", ClusterHitTimeV, &b_ClusterHitTimeV); + fChain->SetBranchAddress("ClusterHitTimeX", ClusterHitTimeX, &b_ClusterHitTimeX); + fChain->SetBranchAddress("ClusterHitSliceU", ClusterHitSliceU, &b_ClusterHitSliceU); + fChain->SetBranchAddress("ClusterHitSliceV", ClusterHitSliceV, &b_ClusterHitSliceV); + fChain->SetBranchAddress("ClusterHitSliceX", ClusterHitSliceX, &b_ClusterHitSliceX); + fChain->SetBranchAddress("nHits", &nHits, &b_nHits); + fChain->SetBranchAddress("RecoHitPos", RecoHitPos, &b_RecoHitPos); + fChain->SetBranchAddress("RecoHitEnergy", RecoHitEnergy, &b_RecoHitEnergy); + fChain->SetBranchAddress("RecoHitSlice", RecoHitSlice, &b_RecoHitSlice); + Notify(); +} + +Bool_t Line_Candidates::Notify() +{ + // The Notify() function is called when a new file is opened. This + // can be either for a new TTree in a TChain or when when a new TTree + // is started when using PROOF. It is normally not necessary to make changes + // to the generated code, but the routine can be extended by the + // user if needed. The return value is currently not used. + + return kTRUE; +} + +void Line_Candidates::Show(Long64_t entry) +{ +// Print contents of entry. +// If entry is not specified, print current entry + if (!fChain) return; + fChain->Show(entry); +} +Int_t Line_Candidates::Cut(Long64_t entry) +{ +// This function may be called from Loop. +// returns 1 if entry is accepted. +// returns -1 otherwise. + entry += 1; + return 1; +} + + +void Line_Candidates::Loop() +{ +// In a ROOT session, you can do: +// root> .L Line_Candidates.C +// root> Line_Candidates t +// root> t.GetEntry(12); // Fill t data members with entry number 12 +// root> t.Show(); // Show values of entry 12 +// root> t.Show(16); // Read and show values of entry 16 +// root> t.Loop(); // Loop on all entries +// + +// This is the loop skeleton where: +// jentry is the global entry number in the chain +// ientry is the entry number in the current Tree +// Note that the argument to GetEntry must be: +// jentry for TChain::GetEntry +// ientry for TTree::GetEntry and TBranch::GetEntry +// +// To read only selected branches, Insert statements like: +// METHOD1: +// fChain->SetBranchStatus("*",0); // disable all branches +// fChain->SetBranchStatus("branchname",1); // activate branchname +// METHOD2: replace line +// fChain->GetEntry(jentry); //read all branches +//by b_branchname->GetEntry(ientry); //read only this branch + if (fChain == 0) return; + + Long64_t nentries = fChain->GetEntriesFast(); + + Long64_t nbytes = 0, nb = 0; + for (Long64_t jentry=0; jentryGetEntry(jentry); nbytes += nb; + // if (Cut(ientry) < 0) continue; + } +} +#endif // #ifdef Line_Candidates_cxx diff --git a/scripts/Validation/Makefile b/scripts/Validation/Makefile new file mode 100644 index 00000000..a9da6260 --- /dev/null +++ b/scripts/Validation/Makefile @@ -0,0 +1,24 @@ +CXX = g++ +ROOTCONFIG := $(shell which root-config) +ROOTCFLAGS := $(shell $(ROOTCONFIG) --cflags) +ROOTLIBS := $(shell $(ROOTCONFIG) --libs) +CXXFLAGS := -g -Wall -Wextra -DNDEBUG -O0 + +SRCS = $(wildcard *.cpp) +OBJS = $(SRCS:.cpp=.o) +EXEC = $(SRCS:.cpp=) + +.PHONY: all clean + +all: $(EXEC) + +#$(EXEC): $(OBJS) +#$(EXEC): %: %.o +$(EXEC): %: %.o + $(CXX) $*.o -o $@ $(ROOTLIBS) + + %.o: %.cpp + $(CXX) $(ROOTCFLAGS) $(CXXFLAGS) -c $< -o $@ + +clean: + rm -f $(OBJS) $(EXEC) diff --git a/scripts/Validation/Plot_True_Track_Length.cpp b/scripts/Validation/Plot_True_Track_Length.cpp new file mode 100644 index 00000000..d18b2377 --- /dev/null +++ b/scripts/Validation/Plot_True_Track_Length.cpp @@ -0,0 +1,313 @@ +#include +#include +#include +#include +#include + +// Root specific +#include +#include +#include +#include + +#include "Truth_Info.h" +#include "Reco_Tree.h" +#include "Line_Candidates.h" + +bool isTMSContained(TVector3 position, bool thin_only = false) { + bool out = true; + // Z positions of the first hits of the TMS + const double TMS_Thin_Start = 11362; + // Where do we transition to the thick region (first layer of scintillator before the change) + const double TMS_Thick_Start = 13500; + // Where does the thick region end + const double TMS_Thick_End = 18294; + const double TMS_Start_Bars_Only[] = {-3350, 240, TMS_Thin_Start}; + const double TMS_End_Bars_Only[] = {3350, -2950, TMS_Thick_End}; + if (position.X() < TMS_Start_Bars_Only[0]) out = false; + if (position.X() > TMS_End_Bars_Only[0]) out = false; + if (position.Y() > TMS_Start_Bars_Only[1]) out = false; + if (position.Y() < TMS_End_Bars_Only[1]) out = false; + if (position.Z() < TMS_Start_Bars_Only[2]) out = false; + if (thin_only) { if (position.Z() > TMS_Thick_Start) out = false; } + else { if (position.Z() > TMS_End_Bars_Only[2]) out = false; } + //std::cout<= 0 && particle_index < truth.nTrueParticles) pdg = truth.PDG[particle_index]; + //else std::cout<<"Found unusual RecoTrackPrimaryParticleIndex: "< 1) std::cout<<"ratio: "< " << std::endl; + return 1; + } + + // Extract input filename and number of events from command line arguments + std::string inputFilename = argv[1]; + int numEvents = -1; + if (argc > 2) numEvents = atoi(argv[2]); + + // Load the tree and make the Truth_Info object + TFile TF(inputFilename.c_str()); + TTree* truth = (TTree*) TF.Get("Truth_Info"); + TTree* reco = (TTree*) TF.Get("Reco_Tree"); + TTree* line_candidates = (TTree*) TF.Get("Line_Candidates"); + bool missing_ttree = false; + if (!truth) { + std::string message = inputFilename + " doesn't contain Truth_Info"; + std::cerr< +#include +#include +#include +#include + +// Root specific +#include +#include +#include + +#include "Truth_Info.h" +#include "Reco_Tree.h" + +bool isTMSContained(TVector3 position, bool thin_only = false) { + bool out = true; + // Z positions of the first hits of the TMS + const double TMS_Thin_Start = 11362; + // Where do we transition to the thick region (first layer of scintillator before the change) + const double TMS_Thick_Start = 13500; + // Where does the thick region end + const double TMS_Thick_End = 18294; + const double TMS_Start_Bars_Only[] = {-3350, 240, TMS_Thin_Start}; + const double TMS_End_Bars_Only[] = {3350, -2950, TMS_Thick_End}; + if (position.X() < TMS_Start_Bars_Only[0]) out = false; + if (position.X() > TMS_End_Bars_Only[0]) out = false; + if (position.Y() > TMS_Start_Bars_Only[1]) out = false; + if (position.Y() < TMS_End_Bars_Only[1]) out = false; + if (position.Z() < TMS_Start_Bars_Only[2]) out = false; + if (thin_only) { if (position.Z() > TMS_Thick_Start) out = false; } + else { if (position.Z() > TMS_End_Bars_Only[2]) out = false; } + //std::cout<= 0 && particle_index < truth.nTrueParticles) pdg = truth.PDG[particle_index]; + //else std::cout<<"Found unusual RecoTrackPrimaryParticleIndex: "< " << std::endl; + return 1; + } + + // Extract input filename and number of events from command line arguments + std::string inputFilename = argv[1]; + int numEvents = -1; + if (argc > 2) numEvents = atoi(argv[2]); + + // Load the tree and make the Truth_Info object + TFile TF(inputFilename.c_str()); + TTree* truth = (TTree*) TF.Get("Truth_Info"); + TTree* reco = (TTree*) TF.Get("Reco_Tree"); + bool missing_ttree = false; + if (!truth) { + std::string message = inputFilename + " doesn't contain Truth_Info"; + std::cerr< +#include +#include + +// Header file for the classes stored in the TTree if any. + +class Reco_Tree { +public : + TTree *fChain; //!pointer to the analyzed TTree or TChain + Int_t fCurrent; //!current Tree number in a TChain + +// Fixed size dimensions of array or collections stored in the TTree if any. + + // Declaration of leaf types + Int_t EventNo; + Int_t SliceNo; + Int_t SpillNo; + Int_t nTracks; + Int_t nHits[20]; //[nTracks] + Float_t TrackHitPos[20][200][3]; //[nTracks] + Float_t StartPos[20][3]; //[nTracks] + Float_t Direction[20][3]; //[nTracks] + Float_t EndPos[20][3]; //[nTracks] + Float_t EnergyRange[20]; //[nTracks] + Float_t EnergyDeposit[20]; //[nTracks] + Float_t Length[20]; //[nTracks] + + // List of branches + TBranch *b_EventNo; //! + TBranch *b_SliceNo; //! + TBranch *b_SpillNo; //! + TBranch *b_nTracks; //! + TBranch *b_nHits; //! + TBranch *b_TrackHitPos; //! + TBranch *b_StartPos; //! + TBranch *b_Direction; //! + TBranch *b_EndPos; //! + TBranch *b_EnergyRange; //! + TBranch *b_EnergyDeposit; //! + TBranch *b_Length; //! + + Reco_Tree(TTree *tree=0); + virtual ~Reco_Tree(); + virtual Int_t Cut(Long64_t entry); + virtual Int_t GetEntry(Long64_t entry); + virtual Long64_t LoadTree(Long64_t entry); + virtual void Init(TTree *tree); + virtual void Loop(); + virtual Bool_t Notify(); + virtual void Show(Long64_t entry = -1); + virtual Long64_t GetEntriesFast() { return fChain->GetEntriesFast(); }; +}; + +#endif + +#define Reco_Tree_cxx +#ifdef Reco_Tree_cxx +Reco_Tree::Reco_Tree(TTree *tree) : fChain(0) +{ +// if parameter tree is not specified (or zero), connect the file +// used to generate this class and read the Tree. + if (tree == 0) { + TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("/pnfs/dune/persistent/users/kleykamp/tmsreco_combined_files/2024-04-19_bfield_0p0T.tmsreco.root"); + if (!f || !f->IsOpen()) { + f = new TFile("/pnfs/dune/persistent/users/kleykamp/tmsreco_combined_files/2024-04-19_bfield_0p0T.tmsreco.root"); + } + f->GetObject("Reco_Tree",tree); + + } + Init(tree); +} + +Reco_Tree::~Reco_Tree() +{ + if (!fChain) return; + //delete fChain->GetCurrentFile(); +} + +Int_t Reco_Tree::GetEntry(Long64_t entry) +{ +// Read contents of entry. + if (!fChain) return 0; + return fChain->GetEntry(entry); +} +Long64_t Reco_Tree::LoadTree(Long64_t entry) +{ +// Set the environment to read one entry + if (!fChain) return -5; + Long64_t centry = fChain->LoadTree(entry); + if (centry < 0) return centry; + if (fChain->GetTreeNumber() != fCurrent) { + fCurrent = fChain->GetTreeNumber(); + Notify(); + } + return centry; +} + +void Reco_Tree::Init(TTree *tree) +{ + // The Init() function is called when the selector needs to initialize + // a new tree or chain. Typically here the branch addresses and branch + // pointers of the tree will be set. + // It is normally not necessary to make changes to the generated + // code, but the routine can be extended by the user if needed. + // Init() will be called many times when running on PROOF + // (once per file to be processed). + + // Set branch addresses and branch pointers + if (!tree) return; + fChain = tree; + fCurrent = -1; + fChain->SetMakeClass(1); + + fChain->SetBranchAddress("EventNo", &EventNo, &b_EventNo); + fChain->SetBranchAddress("SliceNo", &SliceNo, &b_SliceNo); + fChain->SetBranchAddress("SpillNo", &SpillNo, &b_SpillNo); + fChain->SetBranchAddress("nTracks", &nTracks, &b_nTracks); + fChain->SetBranchAddress("nHits", nHits, &b_nHits); + fChain->SetBranchAddress("TrackHitPos", TrackHitPos, &b_TrackHitPos); + fChain->SetBranchAddress("StartPos", StartPos, &b_StartPos); + fChain->SetBranchAddress("Direction", Direction, &b_Direction); + fChain->SetBranchAddress("EndPos", EndPos, &b_EndPos); + fChain->SetBranchAddress("EnergyRange", EnergyRange, &b_EnergyRange); + fChain->SetBranchAddress("EnergyDeposit", EnergyDeposit, &b_EnergyDeposit); + fChain->SetBranchAddress("Length", Length, &b_Length); + Notify(); +} + +Bool_t Reco_Tree::Notify() +{ + // The Notify() function is called when a new file is opened. This + // can be either for a new TTree in a TChain or when when a new TTree + // is started when using PROOF. It is normally not necessary to make changes + // to the generated code, but the routine can be extended by the + // user if needed. The return value is currently not used. + + return kTRUE; +} + +void Reco_Tree::Show(Long64_t entry) +{ +// Print contents of entry. +// If entry is not specified, print current entry + if (!fChain) return; + fChain->Show(entry); +} +Int_t Reco_Tree::Cut(Long64_t entry) +{ +// This function may be called from Loop. +// returns 1 if entry is accepted. +// returns -1 otherwise. + entry += 1; + return 1; +} + +void Reco_Tree::Loop() +{ +// In a ROOT session, you can do: +// root> .L Reco_Tree.C +// root> Reco_Tree t +// root> t.GetEntry(12); // Fill t data members with entry number 12 +// root> t.Show(); // Show values of entry 12 +// root> t.Show(16); // Read and show values of entry 16 +// root> t.Loop(); // Loop on all entries +// + +// This is the loop skeleton where: +// jentry is the global entry number in the chain +// ientry is the entry number in the current Tree +// Note that the argument to GetEntry must be: +// jentry for TChain::GetEntry +// ientry for TTree::GetEntry and TBranch::GetEntry +// +// To read only selected branches, Insert statements like: +// METHOD1: +// fChain->SetBranchStatus("*",0); // disable all branches +// fChain->SetBranchStatus("branchname",1); // activate branchname +// METHOD2: replace line +// fChain->GetEntry(jentry); //read all branches +//by b_branchname->GetEntry(ientry); //read only this branch + if (fChain == 0) return; + + Long64_t nentries = fChain->GetEntriesFast(); + + Long64_t nbytes = 0, nb = 0; + for (Long64_t jentry=0; jentryGetEntry(jentry); nbytes += nb; + // if (Cut(ientry) < 0) continue; + } +} +#endif // #ifdef Reco_Tree_cxx diff --git a/scripts/Validation/Truth_Info.h b/scripts/Validation/Truth_Info.h new file mode 100644 index 00000000..fc830c85 --- /dev/null +++ b/scripts/Validation/Truth_Info.h @@ -0,0 +1,495 @@ +////////////////////////////////////////////////////////// +// This class has been automatically generated on +// Thu May 2 20:30:27 2024 by ROOT version 6.22/08 +// from TTree Truth_Info/Truth_Info +// found on file: ../../neutrino.0_1698176146.edep_TMS_RecoCandidates_Hough_Cluster1.root +////////////////////////////////////////////////////////// + +#ifndef Truth_Info_h +#define Truth_Info_h + +#include +#include +#include + +// Header file for the classes stored in the TTree if any. +#include "string" + +class Truth_Info { +public : + TTree *fChain; //!pointer to the analyzed TTree or TChain + Int_t fCurrent; //!current Tree number in a TChain + +// Fixed size dimensions of array or collections stored in the TTree if any. + + // Declaration of leaf types + Int_t EventNo; + Bool_t IsCC; + std::string *Interaction; + Int_t TruthInfoIndex; + Int_t TruthInfoNSlices; + Int_t nPrimaryVertices; + Bool_t HasPileup; + Int_t NeutrinoPDG; + Float_t NeutrinoP4[4]; + Float_t NeutrinoX4[4]; + Int_t nTrueParticles; + Int_t nTruePrimaryParticles; + Int_t nTrueForgottenParticles; + Int_t VertexID[96]; //[nTrueParticles] + Int_t Parent[96]; //[nTrueParticles] + Int_t TrackId[96]; //[nTrueParticles] + Int_t PDG[96]; //[nTrueParticles] + Bool_t IsPrimary[96]; //[nTrueParticles] + Float_t TrueVisibleEnergy[96]; //[nTrueParticles] + Float_t TruePathLength[96]; //[nTrueParticles] + Float_t TruePathLengthIgnoreY[96]; //[nTrueParticles] + Float_t TruePathLengthInTMS[96]; //[nTrueParticles] + Float_t TruePathLengthInTMSIgnoreY[96]; //[nTrueParticles] + Bool_t InteractionTMSFiducial; + Bool_t InteractionTMSFirstTwoModules; + Bool_t InteractionTMSThin; + Bool_t InteractionLArFiducial; + Bool_t TMSFiducialStart[96]; //[nTrueParticles] + Bool_t TMSFiducialTouch[96]; //[nTrueParticles] + Bool_t TMSFiducialEnd[96]; //[nTrueParticles] + Bool_t LArFiducialStart[96]; //[nTrueParticles] + Bool_t LArFiducialTouch[96]; //[nTrueParticles] + Bool_t LArFiducialEnd[96]; //[nTrueParticles] + Float_t BirthMomentum[96][4]; //[nTrueParticles] + Float_t BirthPosition[96][4]; //[nTrueParticles] + Float_t DeathMomentum[96][4]; //[nTrueParticles] + Float_t DeathPosition[96][4]; //[nTrueParticles] + Float_t MomentumLArStart[96][4]; //[nTrueParticles] + Float_t PositionLArStart[96][4]; //[nTrueParticles] + Float_t MomentumLArEnd[96][4]; //[nTrueParticles] + Float_t PositionLArEnd[96][4]; //[nTrueParticles] + Float_t MomentumTMSStart[96][4]; //[nTrueParticles] + Float_t PositionTMSStart[96][4]; //[nTrueParticles] + Float_t MomentumTMSFirstTwoModulesEnd[96][4]; //[nTrueParticles] + Float_t PositionTMSFirstTwoModulesEnd[96][4]; //[nTrueParticles] + Float_t MomentumTMSThinEnd[96][4]; //[nTrueParticles] + Float_t PositionTMSThinEnd[96][4]; //[nTrueParticles] + Float_t MomentumTMSEnd[96][4]; //[nTrueParticles] + Float_t PositionTMSEnd[96][4]; //[nTrueParticles] + Float_t MomentumZIsLArEnd[96][4]; //[nTrueParticles] + Float_t PositionZIsLArEnd[96][4]; //[nTrueParticles] + Float_t MomentumZIsTMSStart[96][4]; //[nTrueParticles] + Float_t PositionZIsTMSStart[96][4]; //[nTrueParticles] + Float_t MomentumZIsTMSEnd[96][4]; //[nTrueParticles] + Float_t PositionZIsTMSEnd[96][4]; //[nTrueParticles] + Int_t nParticles; + Int_t LeptonPDG; + Float_t LeptonP4[4]; + Float_t LeptonX4[4]; + Float_t MuonP4[4]; + Float_t Muon_Vertex[4]; + Float_t Muon_Death[4]; + Float_t Muon_TrueKE; + Float_t Muon_TrueTrackLength; + Int_t VertexIdOfMostEnergyInEvent; + Float_t VisibleEnergyFromUVertexInSlice; + Float_t TotalVisibleEnergyFromVertex; + Float_t VisibleEnergyFromVVerticesInSlice; + Float_t VertexVisibleEnergyFractionInSlice; + Float_t PrimaryVertexVisibleEnergyFraction; + Int_t RecoTrackN; + Float_t RecoTrackTrueVisibleEnergy[8]; //[RecoTrackN] + Int_t RecoTrackPrimaryParticleIndex[8]; //[RecoTrackN] + Float_t RecoTrackPrimaryParticleTrueVisibleEnergy[8]; //[RecoTrackN] + Int_t RecoTrackSecondaryParticleIndex[8]; //[RecoTrackN] + Float_t RecoTrackSecondaryParticleTrueVisibleEnergy[8]; //[RecoTrackN] + Float_t RecoTrackPrimaryParticleTrueMomentumTrackStart[8][4]; //[RecoTrackN] + Float_t RecoTrackPrimaryParticleTruePositionTrackStart[8][4]; //[RecoTrackN] + Float_t RecoTrackPrimaryParticleTrueMomentumTrackEnd[8][4]; //[RecoTrackN] + Float_t RecoTrackPrimaryParticleTruePositionTrackEnd[8][4]; //[RecoTrackN] + Int_t RecoTrackNHits[8]; //[RecoTrackN] + Float_t RecoTrackTrueHitPosition[8][300][4]; //[RecoTrackN] + Float_t RecoTrackPrimaryParticleTrueTrackLengthAsMeasured[8]; //[RecoTrackN] + Float_t RecoTrackPrimaryParticleTrueTrackLengthAsMeasuredIgnoreY[8]; //[RecoTrackN] + Float_t RecoTrackPrimaryParticleTrueTrackLengthRecoStart[8]; //[RecoTrackN] + Float_t RecoTrackPrimaryParticleTrueTrackLengthRecoStartIgnoreY[8]; //[RecoTrackN] + Float_t RecoTrackPrimaryParticleTrueTrackLengthInTMSIgnoreY[8]; //[RecoTrackN] + Int_t RecoTrackPrimaryParticlePDG[8]; //[RecoTrackN] + Bool_t RecoTrackPrimaryParticleIsPrimary[8]; //[RecoTrackN] + Float_t RecoTrackPrimaryParticleTrueMomentum[8][4]; //[RecoTrackN] + Float_t RecoTrackPrimaryParticleTruePositionStart[8][4]; //[RecoTrackN] + Float_t RecoTrackPrimaryParticleTruePositionEnd[8][4]; //[RecoTrackN] + Float_t RecoTrackPrimaryParticleTrueTrackLength[8]; //[RecoTrackN] + Float_t RecoTrackPrimaryParticleTrueTrackLengthIgnoreY[8]; //[RecoTrackN] + Float_t RecoTrackPrimaryParticleTrueTrackLengthInTMS[8]; //[RecoTrackN] + Float_t RecoTrackPrimaryParticleTrueMomentumEnteringTMS[8][4]; //[RecoTrackN] + Float_t RecoTrackPrimaryParticleTrueMomentumLeavingTMS[8][4]; //[RecoTrackN] + Bool_t RecoTrackPrimaryParticleTMSFiducialStart[8]; //[RecoTrackN] + Bool_t RecoTrackPrimaryParticleTMSFiducialTouch[8]; //[RecoTrackN] + Bool_t RecoTrackPrimaryParticleTMSFiducialEnd[8]; //[RecoTrackN] + Bool_t RecoTrackPrimaryParticleLArFiducialStart[8]; //[RecoTrackN] + Bool_t RecoTrackPrimaryParticleLArFiducialTouch[8]; //[RecoTrackN] + Bool_t RecoTrackPrimaryParticleLArFiducialEnd[8]; //[RecoTrackN] + Int_t RecoTrackSecondaryParticlePDG[8]; //[RecoTrackN] + Bool_t RecoTrackSecondaryParticleIsPrimary[8]; //[RecoTrackN] + Float_t RecoTrackSecondaryParticleTrueMomentum[8][4]; //[RecoTrackN] + Float_t RecoTrackSecondaryParticleTruePositionStart[8][4]; //[RecoTrackN] + Float_t RecoTrackSecondaryParticleTruePositionEnd[8][4]; //[RecoTrackN] + + // List of branches + TBranch *b_EventNo; //! + TBranch *b_IsCC; //! + TBranch *b_Interaction; //! + TBranch *b_TruthInfoIndex; //! + TBranch *b_TruthInfoNSlices; //! + TBranch *b_nPrimaryVertices; //! + TBranch *b_HasPileup; //! + TBranch *b_NeutrinoPDG; //! + TBranch *b_NeutrinoP4; //! + TBranch *b_NeutrinoX4; //! + TBranch *b_nTrueParticles; //! + TBranch *b_nTruePrimaryParticles; //! + TBranch *b_nTrueForgottenParticles; //! + TBranch *b_VertexID; //! + TBranch *b_Parent; //! + TBranch *b_TrackId; //! + TBranch *b_PDG; //! + TBranch *b_IsPrimary; //! + TBranch *b_TrueVisibleEnergy; //! + TBranch *b_TruePathLength; //! + TBranch *b_TruePathLengthIgnoreY; //! + TBranch *b_TruePathLengthInTMS; //! + TBranch *b_TruePathLengthInTMSIgnoreY; //! + TBranch *b_InteractionTMSFiducial; //! + TBranch *b_InteractionTMSFirstTwoModules; //! + TBranch *b_InteractionTMSThin; //! + TBranch *b_InteractionLArFiducial; //! + TBranch *b_TMSFiducialStart; //! + TBranch *b_TMSFiducialTouch; //! + TBranch *b_TMSFiducialEnd; //! + TBranch *b_LArFiducialStart; //! + TBranch *b_LArFiducialTouch; //! + TBranch *b_LArFiducialEnd; //! + TBranch *b_BirthMomentum; //! + TBranch *b_BirthPosition; //! + TBranch *b_DeathMomentum; //! + TBranch *b_DeathPosition; //! + TBranch *b_MomentumLArStart; //! + TBranch *b_PositionLArStart; //! + TBranch *b_MomentumLArEnd; //! + TBranch *b_PositionLArEnd; //! + TBranch *b_MomentumTMSStart; //! + TBranch *b_PositionTMSStart; //! + TBranch *b_MomentumTMSFirstTwoModulesEnd; //! + TBranch *b_PositionTMSFirstTwoModulesEnd; //! + TBranch *b_MomentumTMSThinEnd; //! + TBranch *b_PositionTMSThinEnd; //! + TBranch *b_MomentumTMSEnd; //! + TBranch *b_PositionTMSEnd; //! + TBranch *b_MomentumZIsLArEnd; //! + TBranch *b_PositionZIsLArEnd; //! + TBranch *b_MomentumZIsTMSStart; //! + TBranch *b_PositionZIsTMSStart; //! + TBranch *b_MomentumZIsTMSEnd; //! + TBranch *b_PositionZIsTMSEnd; //! + TBranch *b_nParticles; //! + TBranch *b_LeptonPDG; //! + TBranch *b_LeptonP4; //! + TBranch *b_LeptonX4; //! + TBranch *b_MuonP4; //! + TBranch *b_Muon_Vertex; //! + TBranch *b_Muon_Death; //! + TBranch *b_Muon_TrueKE; //! + TBranch *b_Muon_TrueTrackLength; //! + TBranch *b_VertexIdOfMostEnergyInEvent; //! + TBranch *b_VisibleEnergyFromUVertexInSlice; //! + TBranch *b_TotalVisibleEnergyFromVertex; //! + TBranch *b_VisibleEnergyFromVVerticesInSlice; //! + TBranch *b_VertexVisibleEnergyFractionInSlice; //! + TBranch *b_PrimaryVertexVisibleEnergyFraction; //! + TBranch *b_RecoTrackN; //! + TBranch *b_RecoTrackTrueVisibleEnergy; //! + TBranch *b_RecoTrackPrimaryParticleIndex; //! + TBranch *b_RecoTrackPrimaryParticleTrueVisibleEnergy; //! + TBranch *b_RecoTrackSecondaryParticleIndex; //! + TBranch *b_RecoTrackSecondaryParticleTrueVisibleEnergy; //! + TBranch *b_RecoTrackPrimaryParticleTrueMomentumTrackStart; //! + TBranch *b_RecoTrackPrimaryParticleTruePositionTrackStart; //! + TBranch *b_RecoTrackPrimaryParticleTrueMomentumTrackEnd; //! + TBranch *b_RecoTrackPrimaryParticleTruePositionTrackEnd; //! + TBranch *b_RecoTrackNHits; //! + TBranch *b_RecoTrackTrueHitPosition; //! + TBranch *b_RecoTrackPrimaryParticleTrueTrackLengthAsMeasured; //! + TBranch *b_RecoTrackPrimaryParticleTrueTrackLengthAsMeasuredIgnoreY; //! + TBranch *b_RecoTrackPrimaryParticleTrueTrackLengthRecoStart; //! + TBranch *b_RecoTrackPrimaryParticleTrueTrackLengthRecoStartIgnoreY; //! + TBranch *b_RecoTrackPrimaryParticleTrueTrackLengthInTMSIgnoreY; //! + TBranch *b_RecoTrackPrimaryParticlePDG; //! + TBranch *b_RecoTrackPrimaryParticleIsPrimary; //! + TBranch *b_RecoTrackPrimaryParticleTrueMomentum; //! + TBranch *b_RecoTrackPrimaryParticleTruePositionStart; //! + TBranch *b_RecoTrackPrimaryParticleTruePositionEnd; //! + TBranch *b_RecoTrackPrimaryParticleTrueTrackLength; //! + TBranch *b_RecoTrackPrimaryParticleTrueTrackLengthIgnoreY; //! + TBranch *b_RecoTrackPrimaryParticleTrueTrackLengthInTMS; //! + TBranch *b_RecoTrackPrimaryParticleTrueMomentumEnteringTMS; //! + TBranch *b_RecoTrackPrimaryParticleTrueMomentumLeavingTMS; //! + TBranch *b_RecoTrackPrimaryParticleTMSFiducialStart; //! + TBranch *b_RecoTrackPrimaryParticleTMSFiducialTouch; //! + TBranch *b_RecoTrackPrimaryParticleTMSFiducialEnd; //! + TBranch *b_RecoTrackPrimaryParticleLArFiducialStart; //! + TBranch *b_RecoTrackPrimaryParticleLArFiducialTouch; //! + TBranch *b_RecoTrackPrimaryParticleLArFiducialEnd; //! + TBranch *b_RecoTrackSecondaryParticlePDG; //! + TBranch *b_RecoTrackSecondaryParticleIsPrimary; //! + TBranch *b_RecoTrackSecondaryParticleTrueMomentum; //! + TBranch *b_RecoTrackSecondaryParticleTruePositionStart; //! + TBranch *b_RecoTrackSecondaryParticleTruePositionEnd; //! + + Truth_Info(TTree *tree=0); + virtual ~Truth_Info(); + virtual Int_t Cut(Long64_t entry); + virtual Int_t GetEntry(Long64_t entry); + virtual Long64_t LoadTree(Long64_t entry); + virtual void Init(TTree *tree); + virtual void Loop(); + virtual Bool_t Notify(); + virtual void Show(Long64_t entry = -1); + virtual Long64_t GetEntriesFast() { return fChain->GetEntriesFast(); }; +}; + +#endif + +#define Truth_Info_cxx +#ifdef Truth_Info_cxx +Truth_Info::Truth_Info(TTree *tree) : fChain(0) +{ +// if parameter tree is not specified (or zero), connect the file +// used to generate this class and read the Tree. + if (tree == 0) { + TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("../../neutrino.0_1698176146.edep_TMS_RecoCandidates_Hough_Cluster1.root"); + if (!f || !f->IsOpen()) { + f = new TFile("../../neutrino.0_1698176146.edep_TMS_RecoCandidates_Hough_Cluster1.root"); + } + f->GetObject("Truth_Info",tree); + + } + Init(tree); +} + +Truth_Info::~Truth_Info() +{ + if (!fChain) return; + //delete fChain->GetCurrentFile(); +} + +Int_t Truth_Info::GetEntry(Long64_t entry) +{ +// Read contents of entry. + if (!fChain) return 0; + return fChain->GetEntry(entry); +} +Long64_t Truth_Info::LoadTree(Long64_t entry) +{ +// Set the environment to read one entry + if (!fChain) return -5; + Long64_t centry = fChain->LoadTree(entry); + if (centry < 0) return centry; + if (fChain->GetTreeNumber() != fCurrent) { + fCurrent = fChain->GetTreeNumber(); + Notify(); + } + return centry; +} + +void Truth_Info::Init(TTree *tree) +{ + // The Init() function is called when the selector needs to initialize + // a new tree or chain. Typically here the branch addresses and branch + // pointers of the tree will be set. + // It is normally not necessary to make changes to the generated + // code, but the routine can be extended by the user if needed. + // Init() will be called many times when running on PROOF + // (once per file to be processed). + + // Set object pointer + Interaction = 0; + // Set branch addresses and branch pointers + if (!tree) return; + fChain = tree; + fCurrent = -1; + fChain->SetMakeClass(1); + + fChain->SetBranchAddress("EventNo", &EventNo, &b_EventNo); + fChain->SetBranchAddress("IsCC", &IsCC, &b_IsCC); + fChain->SetBranchAddress("Interaction", &Interaction, &b_Interaction); + fChain->SetBranchAddress("TruthInfoIndex", &TruthInfoIndex, &b_TruthInfoIndex); + fChain->SetBranchAddress("TruthInfoNSlices", &TruthInfoNSlices, &b_TruthInfoNSlices); + fChain->SetBranchAddress("nPrimaryVertices", &nPrimaryVertices, &b_nPrimaryVertices); + fChain->SetBranchAddress("HasPileup", &HasPileup, &b_HasPileup); + fChain->SetBranchAddress("NeutrinoPDG", &NeutrinoPDG, &b_NeutrinoPDG); + fChain->SetBranchAddress("NeutrinoP4", NeutrinoP4, &b_NeutrinoP4); + fChain->SetBranchAddress("NeutrinoX4", NeutrinoX4, &b_NeutrinoX4); + fChain->SetBranchAddress("nTrueParticles", &nTrueParticles, &b_nTrueParticles); + fChain->SetBranchAddress("nTruePrimaryParticles", &nTruePrimaryParticles, &b_nTruePrimaryParticles); + fChain->SetBranchAddress("nTrueForgottenParticles", &nTrueForgottenParticles, &b_nTrueForgottenParticles); + fChain->SetBranchAddress("VertexID", VertexID, &b_VertexID); + fChain->SetBranchAddress("Parent", Parent, &b_Parent); + fChain->SetBranchAddress("TrackId", TrackId, &b_TrackId); + fChain->SetBranchAddress("PDG", PDG, &b_PDG); + fChain->SetBranchAddress("IsPrimary", IsPrimary, &b_IsPrimary); + fChain->SetBranchAddress("TrueVisibleEnergy", TrueVisibleEnergy, &b_TrueVisibleEnergy); + fChain->SetBranchAddress("TruePathLength", TruePathLength, &b_TruePathLength); + fChain->SetBranchAddress("TruePathLengthIgnoreY", TruePathLengthIgnoreY, &b_TruePathLengthIgnoreY); + fChain->SetBranchAddress("TruePathLengthInTMS", TruePathLengthInTMS, &b_TruePathLengthInTMS); + fChain->SetBranchAddress("TruePathLengthInTMSIgnoreY", TruePathLengthInTMSIgnoreY, &b_TruePathLengthInTMSIgnoreY); + fChain->SetBranchAddress("InteractionTMSFiducial", &InteractionTMSFiducial, &b_InteractionTMSFiducial); + fChain->SetBranchAddress("InteractionTMSFirstTwoModules", &InteractionTMSFirstTwoModules, &b_InteractionTMSFirstTwoModules); + fChain->SetBranchAddress("InteractionTMSThin", &InteractionTMSThin, &b_InteractionTMSThin); + fChain->SetBranchAddress("InteractionLArFiducial", &InteractionLArFiducial, &b_InteractionLArFiducial); + fChain->SetBranchAddress("TMSFiducialStart", TMSFiducialStart, &b_TMSFiducialStart); + fChain->SetBranchAddress("TMSFiducialTouch", TMSFiducialTouch, &b_TMSFiducialTouch); + fChain->SetBranchAddress("TMSFiducialEnd", TMSFiducialEnd, &b_TMSFiducialEnd); + fChain->SetBranchAddress("LArFiducialStart", LArFiducialStart, &b_LArFiducialStart); + fChain->SetBranchAddress("LArFiducialTouch", LArFiducialTouch, &b_LArFiducialTouch); + fChain->SetBranchAddress("LArFiducialEnd", LArFiducialEnd, &b_LArFiducialEnd); + fChain->SetBranchAddress("BirthMomentum", BirthMomentum, &b_BirthMomentum); + fChain->SetBranchAddress("BirthPosition", BirthPosition, &b_BirthPosition); + fChain->SetBranchAddress("DeathMomentum", DeathMomentum, &b_DeathMomentum); + fChain->SetBranchAddress("DeathPosition", DeathPosition, &b_DeathPosition); + fChain->SetBranchAddress("MomentumLArStart", MomentumLArStart, &b_MomentumLArStart); + fChain->SetBranchAddress("PositionLArStart", PositionLArStart, &b_PositionLArStart); + fChain->SetBranchAddress("MomentumLArEnd", MomentumLArEnd, &b_MomentumLArEnd); + fChain->SetBranchAddress("PositionLArEnd", PositionLArEnd, &b_PositionLArEnd); + fChain->SetBranchAddress("MomentumTMSStart", MomentumTMSStart, &b_MomentumTMSStart); + fChain->SetBranchAddress("PositionTMSStart", PositionTMSStart, &b_PositionTMSStart); + fChain->SetBranchAddress("MomentumTMSFirstTwoModulesEnd", MomentumTMSFirstTwoModulesEnd, &b_MomentumTMSFirstTwoModulesEnd); + fChain->SetBranchAddress("PositionTMSFirstTwoModulesEnd", PositionTMSFirstTwoModulesEnd, &b_PositionTMSFirstTwoModulesEnd); + fChain->SetBranchAddress("MomentumTMSThinEnd", MomentumTMSThinEnd, &b_MomentumTMSThinEnd); + fChain->SetBranchAddress("PositionTMSThinEnd", PositionTMSThinEnd, &b_PositionTMSThinEnd); + fChain->SetBranchAddress("MomentumTMSEnd", MomentumTMSEnd, &b_MomentumTMSEnd); + fChain->SetBranchAddress("PositionTMSEnd", PositionTMSEnd, &b_PositionTMSEnd); + fChain->SetBranchAddress("MomentumZIsLArEnd", MomentumZIsLArEnd, &b_MomentumZIsLArEnd); + fChain->SetBranchAddress("PositionZIsLArEnd", PositionZIsLArEnd, &b_PositionZIsLArEnd); + fChain->SetBranchAddress("MomentumZIsTMSStart", MomentumZIsTMSStart, &b_MomentumZIsTMSStart); + fChain->SetBranchAddress("PositionZIsTMSStart", PositionZIsTMSStart, &b_PositionZIsTMSStart); + fChain->SetBranchAddress("MomentumZIsTMSEnd", MomentumZIsTMSEnd, &b_MomentumZIsTMSEnd); + fChain->SetBranchAddress("PositionZIsTMSEnd", PositionZIsTMSEnd, &b_PositionZIsTMSEnd); + fChain->SetBranchAddress("nParticles", &nParticles, &b_nParticles); + fChain->SetBranchAddress("LeptonPDG", &LeptonPDG, &b_LeptonPDG); + fChain->SetBranchAddress("LeptonP4", LeptonP4, &b_LeptonP4); + fChain->SetBranchAddress("LeptonX4", LeptonX4, &b_LeptonX4); + fChain->SetBranchAddress("MuonP4", MuonP4, &b_MuonP4); + fChain->SetBranchAddress("Muon_Vertex", Muon_Vertex, &b_Muon_Vertex); + fChain->SetBranchAddress("Muon_Death", Muon_Death, &b_Muon_Death); + fChain->SetBranchAddress("Muon_TrueKE", &Muon_TrueKE, &b_Muon_TrueKE); + fChain->SetBranchAddress("Muon_TrueTrackLength", &Muon_TrueTrackLength, &b_Muon_TrueTrackLength); + fChain->SetBranchAddress("VertexIdOfMostEnergyInEvent", &VertexIdOfMostEnergyInEvent, &b_VertexIdOfMostEnergyInEvent); + fChain->SetBranchAddress("VisibleEnergyFromUVertexInSlice", &VisibleEnergyFromUVertexInSlice, &b_VisibleEnergyFromUVertexInSlice); + fChain->SetBranchAddress("TotalVisibleEnergyFromVertex", &TotalVisibleEnergyFromVertex, &b_TotalVisibleEnergyFromVertex); + fChain->SetBranchAddress("VisibleEnergyFromVVerticesInSlice", &VisibleEnergyFromVVerticesInSlice, &b_VisibleEnergyFromVVerticesInSlice); + fChain->SetBranchAddress("VertexVisibleEnergyFractionInSlice", &VertexVisibleEnergyFractionInSlice, &b_VertexVisibleEnergyFractionInSlice); + fChain->SetBranchAddress("PrimaryVertexVisibleEnergyFraction", &PrimaryVertexVisibleEnergyFraction, &b_PrimaryVertexVisibleEnergyFraction); + fChain->SetBranchAddress("RecoTrackN", &RecoTrackN, &b_RecoTrackN); + fChain->SetBranchAddress("RecoTrackTrueVisibleEnergy", RecoTrackTrueVisibleEnergy, &b_RecoTrackTrueVisibleEnergy); + fChain->SetBranchAddress("RecoTrackPrimaryParticleIndex", RecoTrackPrimaryParticleIndex, &b_RecoTrackPrimaryParticleIndex); + fChain->SetBranchAddress("RecoTrackPrimaryParticleTrueVisibleEnergy", RecoTrackPrimaryParticleTrueVisibleEnergy, &b_RecoTrackPrimaryParticleTrueVisibleEnergy); + fChain->SetBranchAddress("RecoTrackSecondaryParticleIndex", RecoTrackSecondaryParticleIndex, &b_RecoTrackSecondaryParticleIndex); + fChain->SetBranchAddress("RecoTrackSecondaryParticleTrueVisibleEnergy", RecoTrackSecondaryParticleTrueVisibleEnergy, &b_RecoTrackSecondaryParticleTrueVisibleEnergy); + fChain->SetBranchAddress("RecoTrackPrimaryParticleTrueMomentumTrackStart", RecoTrackPrimaryParticleTrueMomentumTrackStart, &b_RecoTrackPrimaryParticleTrueMomentumTrackStart); + fChain->SetBranchAddress("RecoTrackPrimaryParticleTruePositionTrackStart", RecoTrackPrimaryParticleTruePositionTrackStart, &b_RecoTrackPrimaryParticleTruePositionTrackStart); + fChain->SetBranchAddress("RecoTrackPrimaryParticleTrueMomentumTrackEnd", RecoTrackPrimaryParticleTrueMomentumTrackEnd, &b_RecoTrackPrimaryParticleTrueMomentumTrackEnd); + fChain->SetBranchAddress("RecoTrackPrimaryParticleTruePositionTrackEnd", RecoTrackPrimaryParticleTruePositionTrackEnd, &b_RecoTrackPrimaryParticleTruePositionTrackEnd); + fChain->SetBranchAddress("RecoTrackNHits", RecoTrackNHits, &b_RecoTrackNHits); + fChain->SetBranchAddress("RecoTrackTrueHitPosition", RecoTrackTrueHitPosition, &b_RecoTrackTrueHitPosition); + fChain->SetBranchAddress("RecoTrackPrimaryParticleTrueTrackLengthAsMeasured", RecoTrackPrimaryParticleTrueTrackLengthAsMeasured, &b_RecoTrackPrimaryParticleTrueTrackLengthAsMeasured); + fChain->SetBranchAddress("RecoTrackPrimaryParticleTrueTrackLengthAsMeasuredIgnoreY", RecoTrackPrimaryParticleTrueTrackLengthAsMeasuredIgnoreY, &b_RecoTrackPrimaryParticleTrueTrackLengthAsMeasuredIgnoreY); + fChain->SetBranchAddress("RecoTrackPrimaryParticleTrueTrackLengthRecoStart", RecoTrackPrimaryParticleTrueTrackLengthRecoStart, &b_RecoTrackPrimaryParticleTrueTrackLengthRecoStart); + fChain->SetBranchAddress("RecoTrackPrimaryParticleTrueTrackLengthRecoStartIgnoreY", RecoTrackPrimaryParticleTrueTrackLengthRecoStartIgnoreY, &b_RecoTrackPrimaryParticleTrueTrackLengthRecoStartIgnoreY); + fChain->SetBranchAddress("RecoTrackPrimaryParticleTrueTrackLengthInTMSIgnoreY", RecoTrackPrimaryParticleTrueTrackLengthInTMSIgnoreY, &b_RecoTrackPrimaryParticleTrueTrackLengthInTMSIgnoreY); + fChain->SetBranchAddress("RecoTrackPrimaryParticlePDG", RecoTrackPrimaryParticlePDG, &b_RecoTrackPrimaryParticlePDG); + fChain->SetBranchAddress("RecoTrackPrimaryParticleIsPrimary", RecoTrackPrimaryParticleIsPrimary, &b_RecoTrackPrimaryParticleIsPrimary); + fChain->SetBranchAddress("RecoTrackPrimaryParticleTrueMomentum", RecoTrackPrimaryParticleTrueMomentum, &b_RecoTrackPrimaryParticleTrueMomentum); + fChain->SetBranchAddress("RecoTrackPrimaryParticleTruePositionStart", RecoTrackPrimaryParticleTruePositionStart, &b_RecoTrackPrimaryParticleTruePositionStart); + fChain->SetBranchAddress("RecoTrackPrimaryParticleTruePositionEnd", RecoTrackPrimaryParticleTruePositionEnd, &b_RecoTrackPrimaryParticleTruePositionEnd); + fChain->SetBranchAddress("RecoTrackPrimaryParticleTrueTrackLength", RecoTrackPrimaryParticleTrueTrackLength, &b_RecoTrackPrimaryParticleTrueTrackLength); + fChain->SetBranchAddress("RecoTrackPrimaryParticleTrueTrackLengthIgnoreY", RecoTrackPrimaryParticleTrueTrackLengthIgnoreY, &b_RecoTrackPrimaryParticleTrueTrackLengthIgnoreY); + fChain->SetBranchAddress("RecoTrackPrimaryParticleTrueTrackLengthInTMS", RecoTrackPrimaryParticleTrueTrackLengthInTMS, &b_RecoTrackPrimaryParticleTrueTrackLengthInTMS); + fChain->SetBranchAddress("RecoTrackPrimaryParticleTrueMomentumEnteringTMS", RecoTrackPrimaryParticleTrueMomentumEnteringTMS, &b_RecoTrackPrimaryParticleTrueMomentumEnteringTMS); + fChain->SetBranchAddress("RecoTrackPrimaryParticleTrueMomentumLeavingTMS", RecoTrackPrimaryParticleTrueMomentumLeavingTMS, &b_RecoTrackPrimaryParticleTrueMomentumLeavingTMS); + fChain->SetBranchAddress("RecoTrackPrimaryParticleTMSFiducialStart", RecoTrackPrimaryParticleTMSFiducialStart, &b_RecoTrackPrimaryParticleTMSFiducialStart); + fChain->SetBranchAddress("RecoTrackPrimaryParticleTMSFiducialTouch", RecoTrackPrimaryParticleTMSFiducialTouch, &b_RecoTrackPrimaryParticleTMSFiducialTouch); + fChain->SetBranchAddress("RecoTrackPrimaryParticleTMSFiducialEnd", RecoTrackPrimaryParticleTMSFiducialEnd, &b_RecoTrackPrimaryParticleTMSFiducialEnd); + fChain->SetBranchAddress("RecoTrackPrimaryParticleLArFiducialStart", RecoTrackPrimaryParticleLArFiducialStart, &b_RecoTrackPrimaryParticleLArFiducialStart); + fChain->SetBranchAddress("RecoTrackPrimaryParticleLArFiducialTouch", RecoTrackPrimaryParticleLArFiducialTouch, &b_RecoTrackPrimaryParticleLArFiducialTouch); + fChain->SetBranchAddress("RecoTrackPrimaryParticleLArFiducialEnd", RecoTrackPrimaryParticleLArFiducialEnd, &b_RecoTrackPrimaryParticleLArFiducialEnd); + fChain->SetBranchAddress("RecoTrackSecondaryParticlePDG", RecoTrackSecondaryParticlePDG, &b_RecoTrackSecondaryParticlePDG); + fChain->SetBranchAddress("RecoTrackSecondaryParticleIsPrimary", RecoTrackSecondaryParticleIsPrimary, &b_RecoTrackSecondaryParticleIsPrimary); + fChain->SetBranchAddress("RecoTrackSecondaryParticleTrueMomentum", RecoTrackSecondaryParticleTrueMomentum, &b_RecoTrackSecondaryParticleTrueMomentum); + fChain->SetBranchAddress("RecoTrackSecondaryParticleTruePositionStart", RecoTrackSecondaryParticleTruePositionStart, &b_RecoTrackSecondaryParticleTruePositionStart); + fChain->SetBranchAddress("RecoTrackSecondaryParticleTruePositionEnd", RecoTrackSecondaryParticleTruePositionEnd, &b_RecoTrackSecondaryParticleTruePositionEnd); + Notify(); +} + +Bool_t Truth_Info::Notify() +{ + // The Notify() function is called when a new file is opened. This + // can be either for a new TTree in a TChain or when when a new TTree + // is started when using PROOF. It is normally not necessary to make changes + // to the generated code, but the routine can be extended by the + // user if needed. The return value is currently not used. + + return kTRUE; +} + +void Truth_Info::Show(Long64_t entry) +{ +// Print contents of entry. +// If entry is not specified, print current entry + if (!fChain) return; + fChain->Show(entry); +} +Int_t Truth_Info::Cut(Long64_t entry) +{ +// This function may be called from Loop. +// returns 1 if entry is accepted. +// returns -1 otherwise. + entry += 1; + return 1; +} + + +void Truth_Info::Loop() +{ +// In a ROOT session, you can do: +// root> .L Truth_Info.C +// root> Truth_Info t +// root> t.GetEntry(12); // Fill t data members with entry number 12 +// root> t.Show(); // Show values of entry 12 +// root> t.Show(16); // Read and show values of entry 16 +// root> t.Loop(); // Loop on all entries +// + +// This is the loop skeleton where: +// jentry is the global entry number in the chain +// ientry is the entry number in the current Tree +// Note that the argument to GetEntry must be: +// jentry for TChain::GetEntry +// ientry for TTree::GetEntry and TBranch::GetEntry +// +// To read only selected branches, Insert statements like: +// METHOD1: +// fChain->SetBranchStatus("*",0); // disable all branches +// fChain->SetBranchStatus("branchname",1); // activate branchname +// METHOD2: replace line +// fChain->GetEntry(jentry); //read all branches +//by b_branchname->GetEntry(ientry); //read only this branch + if (fChain == 0) return; + + Long64_t nentries = fChain->GetEntriesFast(); + + Long64_t nbytes = 0, nb = 0; + for (Long64_t jentry=0; jentryGetEntry(jentry); nbytes += nb; + // if (Cut(ientry) < 0) continue; + } +} +#endif // #ifdef Truth_Info_cxx diff --git a/scripts/Validation/simply_draw_everything.py b/scripts/Validation/simply_draw_everything.py new file mode 100644 index 00000000..e782bdac --- /dev/null +++ b/scripts/Validation/simply_draw_everything.py @@ -0,0 +1,53 @@ +import os + +import ROOT +ROOT.gROOT.SetBatch(True) + +import sys + +def draw_histograms(input_file): + # Open the input ROOT file + root_file = ROOT.TFile.Open(input_file) + + # Create a directory to save images + output_dir = os.path.splitext(input_file)[0] + "_images" + os.makedirs(output_dir, exist_ok=True) + + canvas = ROOT.TCanvas("canvas", "canvas", 800, 600) + + # Loop over all keys in the ROOT file + for key in root_file.GetListOfKeys(): + obj = key.ReadObj() + if isinstance(obj, ROOT.TH2): + # For 2D histograms, draw with "colz" option and save as png + obj.GetYaxis().SetTitleOffset(1.4) + obj.GetZaxis().SetTitleOffset(0.5) + obj.Draw("colz") + canvas.Print(os.path.join(output_dir, obj.GetName() + ".png")) + elif isinstance(obj, ROOT.TH1): + # For 1D histograms, draw and save as png + obj.Draw() + canvas.Print(os.path.join(output_dir, obj.GetName() + ".png")) + + # Close the input ROOT file + root_file.Close() + + + +if __name__ == "__main__": + # Check if input file is provided as argument + if len(sys.argv) != 2: + print("Usage: python script.py ") + sys.exit(1) + + input_file = sys.argv[1] + if not os.path.isfile(input_file): + print("Error: Input file does not exist!") + sys.exit(1) + + # Initialize ROOT + ROOT.gROOT.SetBatch(True) # Prevent ROOT from trying to open X11 windows + ROOT.gStyle.SetOptStat(0) # Hide statistics box in histograms + + # Call function to draw histograms + draw_histograms(input_file) From 4e7f80d1b202cd8874780a256752aa813746bc7a Mon Sep 17 00:00:00 2001 From: Jeffrey Kleykamp Date: Thu, 30 May 2024 22:15:51 -0500 Subject: [PATCH 27/30] Add LArSlicer.C --- scripts/TTreeSlicer/LArSlicer.C | 78 +++++++++++++++++++++++++++++++++ 1 file changed, 78 insertions(+) create mode 100644 scripts/TTreeSlicer/LArSlicer.C diff --git a/scripts/TTreeSlicer/LArSlicer.C b/scripts/TTreeSlicer/LArSlicer.C new file mode 100644 index 00000000..c0206df9 --- /dev/null +++ b/scripts/TTreeSlicer/LArSlicer.C @@ -0,0 +1,78 @@ +#include +#include +#include +#include +#include +#include + +void selectEvents(const char* inputFile, const char* outputFilename, int nmax) { + // Open the input file + TFile *file = TFile::Open(inputFile); + if (!file || file->IsZombie()) { + printf("Error: Unable to open input file %s\n", inputFile); + return; + } + + // Get the input tree + TTree *inputTree = (TTree*)file->Get("Truth_Info"); + if (!inputTree) { + printf("Error: Unable to retrieve input Truth_Info\n"); + file->Close(); + return; + } + + // Create a new file + TFile *outputFile = new TFile(outputFilename, "RECREATE"); + if (!outputFile || outputFile->IsZombie()) { + printf("Error: Unable to create output file %s\n", outputFilename); + file->Close(); + return; + } + + // Clone the input tree + TTree *outputTree = inputTree->CloneTree(0); // Create an empty tree with the same branches + + // Define the cut condition + // Example: Select events where positionX, positionY, and positionZ are within a certain range + const double LAr_Start_Exact[] = {-3478.48, -2166.71, 4179.24}; + const double LAr_End_Exact[] = {3478.48, 829.282, 9135.88}; + Double_t minX = LAr_Start_Exact[0], maxX = LAr_End_Exact[0]; + Double_t minY = LAr_Start_Exact[1], maxY = LAr_End_Exact[1]; + Double_t minZ = LAr_Start_Exact[2], maxZ = LAr_End_Exact[2]; + + // Set up TTreeReader for input tree + //TTreeReader reader(inputTree); + //inputTree->SetBranchStatus("*", true); + //TTreeReaderArray X4(reader, "LeptonX4"); + + float X4[4]; + inputTree->SetBranchAddress("LeptonX4", &X4); + + // Loop over events + long n = 0; + while (n < inputTree->GetEntries() && (nmax < 0 || n < nmax)) { + inputTree->GetEntry(n); + // Access the position components + Double_t posX = X4[0]; + Double_t posY = X4[1]; + Double_t posZ = X4[2]; + // Apply the cut condition + if (posX >= minX && posX <= maxX && + posY >= minY && posY <= maxY && + posZ >= minZ && posZ <= maxZ) { + outputTree->Fill(); // Fill the output tree with selected events + } + n += 1; + } + + // Write the output tree to the output file + outputTree->Write(); + + // Clean up + outputFile->Close(); + file->Close(); +} + +void LArSlicer(const char* inputFile, const char* outputFile, int nmax = -1) { + selectEvents(inputFile, outputFile, nmax); +} From ae6a3d1db775c73ef793c28a5011129f1268ab43 Mon Sep 17 00:00:00 2001 From: Jeffrey Kleykamp Date: Mon, 10 Jun 2024 15:16:26 -0500 Subject: [PATCH 28/30] Fix issue with GetMaterials getting stuck with tiny step sizes, leading to memory hold --- src/TMS_Geom.h | 25 +++++++++++++++++++++++-- 1 file changed, 23 insertions(+), 2 deletions(-) diff --git a/src/TMS_Geom.h b/src/TMS_Geom.h index 05e1442c..3b57d0ea 100644 --- a/src/TMS_Geom.h +++ b/src/TMS_Geom.h @@ -16,8 +16,8 @@ #include "TMS_Constants.h" #define __GEOM_LARGE_STEP__ 1E10 -#define __GEOM_SMALL_STEP__ 1E-5 -#define __GEOM_TINY_STEP__ 1E-10 +#define __GEOM_SMALL_STEP__ 1E-4 +#define __GEOM_TINY_STEP__ 1E-5 // Define the TMS geometry singleton class TMS_Geom { @@ -92,6 +92,21 @@ class TMS_Geom { static bool StaticIsInsideReasonableSize(TVector3 position) { return TMS_Geom::GetInstance().IsInsideReasonableSize(position); }; + std::string GetNameOfDetector(const TVector3 &point) { + std::string out = ""; + if (out == "" && IsInsideLAr(point)) out += "LAr"; + if (out == "" && IsInsideTMS(point)) out += "TMS"; + if (out == "" && IsInsideTMSMass(point)) out += "TMS Mass"; + if (out == "" && StaticIsInsideReasonableSize(point)) out += "Reasonable Size"; + if (out == "") out += "Other"; + out += " ("; + out += std::to_string(point.X()) + ", "; + out += std::to_string(point.Y()) + ", "; + out += std::to_string(point.Z()); + out += ")"; + return out; + }; + // Get the geometry TGeoManager* GetGeometry() { @@ -227,6 +242,7 @@ class TMS_Geom { geom->FindNode(); // How big was the step in this material step = geom->GetStep(); + //std::cout<<"Current total: "< 0 if (nodes.size() > 1) { // Loop over all pairs of vectors for (size_t i = 0; i < nodes.size() - 1; i++) { auto p1 = nodes.at(i); auto p2 = nodes.at(i+1); + //std::cout<<"p1 is in "< Date: Mon, 10 Jun 2024 16:37:49 -0500 Subject: [PATCH 29/30] Fix TMS thick end to be past reco hits --- config/TMS_Default_Config.toml | 2 +- src/TMS_Constants.h | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/config/TMS_Default_Config.toml b/config/TMS_Default_Config.toml index ec2ef7de..7316f164 100644 --- a/config/TMS_Default_Config.toml +++ b/config/TMS_Default_Config.toml @@ -111,7 +111,7 @@ [Fiducial.TMS.End] X = 3300.0 # mm Y = 160.0 # mm - Z = 18294.0 # mm + Z = 18314.0 # mm [Fiducial.LAr.Start] X = -3478.48 # mm Y = -2166.71 # mm diff --git a/src/TMS_Constants.h b/src/TMS_Constants.h index 803daf54..4c0d4397 100644 --- a/src/TMS_Constants.h +++ b/src/TMS_Constants.h @@ -50,7 +50,7 @@ namespace TMS_Const { // Where do we transition to the thick region (first layer of scintillator before the change) const double TMS_Thick_Start = 13500; // Where does the thick region end - const double TMS_Thick_End = 18294; + const double TMS_Thick_End = 18314; // Approximate starting and end positions of TMS detector in geometry for plotting hits, in {x,y,z} // in mm by default! From 58243b132e9d847e49613b3f95e8ef30eebeed9d Mon Sep 17 00:00:00 2001 From: Liam O'Sullivan Date: Tue, 2 Jul 2024 07:45:11 -0500 Subject: [PATCH 30/30] stashing --- src/TMS_TreeWriter.cpp | 1 + src/TMS_TrueParticle.cpp | 12 ++++++++---- 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/src/TMS_TreeWriter.cpp b/src/TMS_TreeWriter.cpp index db22f1d7..0928aebb 100644 --- a/src/TMS_TreeWriter.cpp +++ b/src/TMS_TreeWriter.cpp @@ -1267,6 +1267,7 @@ void TMS_TreeWriter::Fill(TMS_Event &event) { if (itTrack < __TMS_MAX_LINES__) { setMomentum(RecoTrackPrimaryParticleTrueMomentumTrackStart[itTrack], tp.GetMomentumAtZ(start_z, max_z_distance)); setPosition(RecoTrackPrimaryParticleTruePositionTrackStart[itTrack], tp.GetPositionAtZ(start_z, max_z_distance)); + //std::cout << "Setting tp shite: " << tp.GetPositionAtZ(start_z, max_z_distance).X() << " " << tp.GetPositionAtZ(start_z, max_z_distance).Y() << " " << tp.GetPositionAtZ(start_z, max_z_distance).Z() << ",\t" << start_z << " " << max_z_distance << std::endl; setMomentum(RecoTrackPrimaryParticleTrueMomentumTrackEnd[itTrack], tp.GetMomentumAtZ(end_z, max_z_distance)); setPosition(RecoTrackPrimaryParticleTruePositionTrackEnd[itTrack], tp.GetPositionAtZ(end_z, max_z_distance)); diff --git a/src/TMS_TrueParticle.cpp b/src/TMS_TrueParticle.cpp index 1eeefd35..5ed1a323 100644 --- a/src/TMS_TrueParticle.cpp +++ b/src/TMS_TrueParticle.cpp @@ -71,16 +71,20 @@ std::vector TMS_TrueParticle::GetPositionPoints(double z_start, double TLorentzVector TMS_TrueParticle::GetPositionAtZ(double z, double max_z_dist) { TLorentzVector out(-99999999, -99999999, -99999999, -99999999); double z_dist_found = 999999; + //double prev_z = 0; for (size_t i = 0; i < GetPositionPoints().size(); i++) { - double distance = abs(GetPositionPoints()[i].Z() - z); - if (distance <= max_z_dist) { + //std::cout << "pos: " << GetPositionPoints()[i].X() << ", "<< GetPositionPoints()[i].Y() << ", "<< GetPositionPoints()[i].Z() << " " << GetPositionPoints()[i].Z() - prev_z << std::endl; + //prev_z = GetPositionPoints()[i].Z(); + double distance = GetPositionPoints()[i].Z() - z; + if (abs(distance) <= max_z_dist) { // Found a candidate - if (distance < z_dist_found) { + if (abs(distance) < z_dist_found) { out = GetPositionPoints()[i]; z_dist_found = distance; } } - } + } + out.SetZ( out.Z() - z_dist_found ); return out; }