From 6d9d3ed5e0014830ee946a7e41851792f91988eb Mon Sep 17 00:00:00 2001 From: anupama-reghunath Date: Thu, 28 Nov 2024 21:58:25 +0100 Subject: [PATCH 1/7] analysis_toolkit;Added pre- selection function as a start --- python/analysis_toolkit.py | 241 +++++++++++++++++++++++++++++++++++++ 1 file changed, 241 insertions(+) create mode 100644 python/analysis_toolkit.py diff --git a/python/analysis_toolkit.py b/python/analysis_toolkit.py new file mode 100644 index 000000000..e63de8d15 --- /dev/null +++ b/python/analysis_toolkit.py @@ -0,0 +1,241 @@ +#!/usr/bin/env python +"""Toolkit for Analysis.""" + +import numpy as np +import ROOT +import shipunit as u +import yaml +from rootpyPickler import Unpickler +from ShipGeoConfig import AttrDict + + +class selection_check: + """Class to perform various selection checks on the candidate.""" + + def __init__(self, g): + """Initialize the selection_check class with geometry and configuration.""" + geo_file = ROOT.TFile.Open(g, "read") + self.geometry_manager = geo_file.FAIRGeom + unpickler = Unpickler(geo_file) + self.ship_geo = unpickler.load("ShipGeo") + + fairship = ROOT.gSystem.Getenv("FAIRSHIP") + + if self.ship_geo.DecayVolumeMedium == "helium": + with open(fairship + "/geometry/veto_config_helium.yaml", "r") as file: + config = yaml.safe_load(file) + self.veto_geo = AttrDict(config) + self.veto_geo.z0 + if self.ship_geo.DecayVolumeMedium == "vacuums": + with open(fairship + "/geometry/veto_config_vacuums.yaml", "r") as file: + config = yaml.safe_load(file) + self.veto_geo = AttrDict(config) + self.veto_geo.z0 + + def access_event(self, tree): + """Access event data.""" + self.tree = tree + + def define_candidate_time(self, candidate): + """Calculate time associated with the candidate decay vertex using strawtubes MCPoint info.""" + t0 = self.tree.ShipEventHeader.GetEventTime() + candidate_pos = ROOT.TVector3() + candidate.GetVertex(candidate_pos) + + d1, d2 = candidate.GetDaughter(0), candidate.GetDaughter(1) + d1_mc, d2_mc = self.tree.fitTrack2MC[d1], self.tree.fitTrack2MC[d2] + + time_vtx_from_strawhits = [] + + for hit in self.tree.strawtubesPoint: + if not ( + int(str(hit.GetDetectorID())[:1]) == 1 + or int(str(hit.GetDetectorID())[:1]) == 2 + ): + continue + + if not (hit.GetTrackID() == d1_mc or hit.GetTrackID() == d2_mc): + continue + + t_straw = hit.GetTime() + + dist = np.sqrt( + (candidate_pos.X() - hit.GetX()) ** 2 + + (candidate_pos.Y() - hit.GetY()) ** 2 + + (candidate_pos.Z() - hit.GetZ()) ** 2 + ) # distance to the vertex #in cm + + d_mom = self.tree.MCTrack[hit.GetTrackID()].GetP() / u.GeV + mass = self.tree.MCTrack[hit.GetTrackID()].GetMass() + v = u.c_light * d_mom / np.sqrt(d_mom**2 + (mass) ** 2) + + t_vertex = t_straw - (dist / v) + + time_vtx_from_strawhits.append(t_vertex) + + t_vtx = np.average(time_vtx_from_strawhits) + t0 + + return t_vtx # units in ns + + def impact_parameter(self, candidate): + """Calculate the impact parameter of the candidate relative to (0,0,target.z0).""" + candidate_pos = ROOT.TVector3() + candidate.GetVertex(candidate_pos) + + candidate_mom = ROOT.TLorentzVector() + candidate.Momentum(candidate_mom) + target_point = ROOT.TVector3(0, 0, self.ship_geo.target.z0) + + projection_factor = 0 + if hasattr(candidate_mom, "P"): + P = candidate_mom.P() + else: + P = candidate_mom.Mag() + for i in range(3): + projection_factor += ( + candidate_mom(i) / P * (target_point(i) - candidate_pos(i)) + ) + + dist = 0 + for i in range(3): + dist += ( + target_point(i) + - candidate_pos(i) + - projection_factor * candidate_mom(i) / P + ) ** 2 + dist = ROOT.TMath.Sqrt(dist) + + return dist # in cm + + def dist_to_innerwall(self, candidate): + """Calculate the minimum distance(in XY plane) of the candidate decay vertex to the inner wall of the decay vessel. If outside the decay volume, or if distance > 100cm,Return 0.""" + candidate_pos = ROOT.TVector3() + candidate.GetVertex(candidate_pos) + position = (candidate_pos.X(), candidate_pos.Y(), candidate_pos.Z()) + + nsteps = 8 + dalpha = 2 * ROOT.TMath.Pi() / nsteps + min_distance = float("inf") + + node = self.geometry_manager.FindNode(*position) + if not node: + return 0 # is outside the decay volume + + # Loop over directions in the XY plane + for n in range(nsteps): + alpha = n * dalpha + direction = ( + ROOT.TMath.Sin(alpha), + ROOT.TMath.Cos(alpha), + 0.0, + ) # Direction vector in XY plane + self.geometry_manager.InitTrack(*position, *direction) + if not self.geometry_manager.FindNextBoundary(): + continue + # Get the distance to the boundary and update the minimum distance + distance = self.geometry_manager.GetStep() + min_distance = min(min_distance, distance) + + return min_distance if min_distance < 100 * u.m else 0 + + def dist_to_entrance_vessel(self, candidate): + """Calculate the distance of the candidate decay vertex to the entrance of the decay vessel.""" + candidate_pos = ROOT.TVector3() + candidate.GetVertex(candidate_pos) + return candidate_pos.Z() - self.veto_geo.z0 + + def nDOF(self, candidate): + """Return the number of degrees of freedom (nDOF) for the particle's daughter tracks.""" + flag = True + nmeas = [] + t1, t2 = candidate.GetDaughter(0), candidate.GetDaughter(1) + ndf_cut = 25 # cut on nDOF + for tr in [t1, t2]: + fit_status = self.tree.FitTracks[tr].getFitStatus() + nmeas.append( + int(round(fit_status.getNdf())) + ) # nmeas.append(fit_status.getNdf()) + if nmeas[-1] <= ndf_cut: + flag = False # dof<=25=False + return nmeas, flag + + def daughtermomentum(self, candidate): + """Return the momentum(Mag) of the particle's daughter tracks.""" + flag = True + daughter_mom = [] + t1, t2 = candidate.GetDaughter(0), candidate.GetDaughter(1) + for trD in [t1, t2]: + x = self.tree.FitTracks[trD] + xx = x.getFittedState() + daughter_mom.append((xx.getMom().Mag())) + if daughter_mom[-1] < 1 * u.GeV: + flag = False + return daughter_mom, flag + + def invariant_mass(self, candidate): + """Invariant mass of the candidate.""" + return candidate.GetMass() + + def DOCA(self, candidate): + """Distance of Closest Approach.""" + return candidate.GetDoca() + + def is_in_fiducial(self, candidate): + """Check if the candidate decay vertex is within the Fiducial Volume.""" + candidate_pos = ROOT.TVector3() + candidate.GetVertex(candidate_pos) + + if candidate_pos.Z() > self.ship_geo.TrackStation1.z: + return False + if candidate_pos.Z() < self.veto_geo.z0: + return False + + # if self.dist2InnerWall(candidate)<=5*u.cm: return False + + vertex_node = ROOT.gGeoManager.FindNode( + candidate_pos.X(), candidate_pos.Y(), candidate_pos.Z() + ) + vertex_elem = vertex_node.GetVolume().GetName() + if not vertex_elem.startswith("DecayVacuum_"): + return False + return True + + def chi2nDOF(self, candidate): + """Return the reduced chi^2 of the particle's daughter tracks.""" + t1, t2 = candidate.GetDaughter(0), candidate.GetDaughter(1) + # cut on chi2nDOF + flag = True + chi2ndf_cut = 5 + chi2ndf = [] + for tr in [t1, t2]: + fit_status = self.tree.FitTracks[tr].getFitStatus() + chi2ndf.append(fit_status.getChi2() / fit_status.getNdf()) + if chi2ndf[-1] >= chi2ndf_cut: + flag = False + return chi2ndf, flag + + def preselection_cut(self, candidate, IP_cut=250): + """Umbrella method to apply the pre-selection cuts on the candidate.""" + if len(self.tree.Particles) != 1: + return False + if self.dist_to_innerwall(candidate) <= 5 * u.cm: + return False + if not (self.is_in_fiducial(candidate)): + return False + if self.dist_to_entrance_vessel(candidate) < 100 * u.cm: + return False + if self.impact_parameter(candidate) > IP_cut: + return False + nmeas, flag = self.nDOF(candidate) + if not (flag): + return False + chi2ndf, flag = self.chi2nDOF(candidate) + if not (flag): + return False + d_mom, flag = self.daughtermomentum(candidate) + if not (flag): + return False + if self.DOCA(candidate) >= 1 * u.cm: + return False + + return True From 7d891d669497e9a37dc8a66c1183e11d19676d04 Mon Sep 17 00:00:00 2001 From: anupama-reghunath Date: Thu, 28 Nov 2024 23:33:26 +0100 Subject: [PATCH 2/7] add example script for analysis_toolkit --- macro/analysis_example.py | 41 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 41 insertions(+) create mode 100644 macro/analysis_example.py diff --git a/macro/analysis_example.py b/macro/analysis_example.py new file mode 100644 index 000000000..95b85503f --- /dev/null +++ b/macro/analysis_example.py @@ -0,0 +1,41 @@ +#!/usr/bin/env python +"""Example script for usage of the analysis_toolkit for signal selection .""" + +import analysis_toolkit +import ROOT + +geo_file = "geofile_full.conical.Pythia8-TGeant4.root" +input_file = "ship.conical.Pythia8-TGeant4_rec.root" + +f = ROOT.TFile.Open(input_file) +tree = f.cbmsim + +selection = analysis_toolkit.selection_check(geo_file) + +for event_nr, event in enumerate(tree): + if event_nr > 200: + break + + selection.access_event(event) + + if len(event.Particles) == 0: + continue + + for candidate_id_in_event, signal in enumerate(tree.Particles): + print(f"Event:{event_nr} Candidate_index: {candidate_id_in_event}") + + print(f"\t vertex time: {selection.define_candidate_time(signal)} ns") + print(f"\t Impact Parameter: {selection.impact_parameter(signal)} cm") + print(f"\t is within Fiducial volume: {selection.is_in_fiducial(signal)}") + print(f"\t\tDist2InnerWall: {selection.dist_to_innerwall(signal)} cm") + print(f"\t\tDist2EntranceLid: {selection.dist_to_entrance_vessel(signal)} cm") + print(f"\t Daughter momentum>1 GeV: {selection.daughtermomentum(signal)} GeV") + print(f"\t Invariant mass: {selection.invariant_mass(signal)} GeV") + print(f"\t Degrees of Freedom ([d1,d2],Cut(>25)): {selection.nDOF(signal)}") + print(f"\t Reduced Chi^2 ([d1,d2],Cut(<5): {selection.chi2nDOF(signal)}") + print(f"\t DOCA: {selection.DOCA(signal)} cm") + + print( + f"\t Preselection Cut passed:{selection.preselection_cut(signal,IP_cut=250)}" + ) + print("------------------------------------------------------------") From e8254a8798a991563b167d6197f60688800145e6 Mon Sep 17 00:00:00 2001 From: anupama-reghunath Date: Fri, 29 Nov 2024 14:37:27 +0100 Subject: [PATCH 3/7] tabulate pre-selection parameters --- macro/analysis_example.py | 111 ++++++++++++++++++++------ python/analysis_toolkit.py | 154 ++++++++++++++++++++++++++++--------- 2 files changed, 204 insertions(+), 61 deletions(-) diff --git a/macro/analysis_example.py b/macro/analysis_example.py index 95b85503f..1499e859b 100644 --- a/macro/analysis_example.py +++ b/macro/analysis_example.py @@ -1,41 +1,106 @@ #!/usr/bin/env python """Example script for usage of the analysis_toolkit for signal selection .""" +from argparse import ArgumentParser + import analysis_toolkit import ROOT +import rootUtils as ut + +parser = ArgumentParser(description=__doc__) +parser.add_argument("--path", help="Path to simulation file", default="./") +options = parser.parse_args() -geo_file = "geofile_full.conical.Pythia8-TGeant4.root" -input_file = "ship.conical.Pythia8-TGeant4_rec.root" -f = ROOT.TFile.Open(input_file) +geo_file = options.path + "/geofile_full.conical.Pythia8-TGeant4.root" +input_file = options.path + "/ship.conical.Pythia8-TGeant4_rec.root" + +f = ROOT.TFile.Open(input_file, "read") tree = f.cbmsim selection = analysis_toolkit.selection_check(geo_file) -for event_nr, event in enumerate(tree): - if event_nr > 200: - break +hist_dict = {} +ut.bookHist(hist_dict, "event_weight", " Event weight; ; ", 100, 100, 100) +ut.bookHist( + hist_dict, + "candidate_time", + " candidate time @ decay vertex; ns ; ", + 200, + 100, + 100, +) +ut.bookHist(hist_dict, "impact_parameter", "Impact Parameter;cm;", 200, 100, 100) +ut.bookHist(hist_dict, "dist_to_innerwall", "Distance to Innerwall;cm;", 200, 100, 100) +ut.bookHist( + hist_dict, + "dist_to_vesselentrance", + "Distance to Decay Vessel Entrance;cm;", + 200, + 100, + 100, +) +ut.bookHist(hist_dict, "inv_mass", "Invariant mass ;GeV;", 200, 100, 100) +ut.bookHist(hist_dict, "DOCA", "Distance of Closest Approach ;cm;", 200, 100, 100) +ut.bookHist( + hist_dict, + "len_Particles", + " len(tree.Particles) ;Number of candidates per event;", + 200, + 100, + 100, +) +ut.bookHist( + hist_dict, + "d_mom", + "momentum of Daughters ;d1 (GeV);d2 (GeV)", + 200, + 100, + 100, + 200, + 100, + 100, +) +ut.bookHist(hist_dict, "nDOF", "nDOF ;d1;d2", 200, 100, 100, 200, 100, 100) +ut.bookHist(hist_dict, "chi2nDOF", "chi2nDOF ;d1;d2", 200, 100, 100, 200, 100, 100) + + +for event_nr, event in enumerate(tree): selection.access_event(event) if len(event.Particles) == 0: continue - for candidate_id_in_event, signal in enumerate(tree.Particles): - print(f"Event:{event_nr} Candidate_index: {candidate_id_in_event}") - - print(f"\t vertex time: {selection.define_candidate_time(signal)} ns") - print(f"\t Impact Parameter: {selection.impact_parameter(signal)} cm") - print(f"\t is within Fiducial volume: {selection.is_in_fiducial(signal)}") - print(f"\t\tDist2InnerWall: {selection.dist_to_innerwall(signal)} cm") - print(f"\t\tDist2EntranceLid: {selection.dist_to_entrance_vessel(signal)} cm") - print(f"\t Daughter momentum>1 GeV: {selection.daughtermomentum(signal)} GeV") - print(f"\t Invariant mass: {selection.invariant_mass(signal)} GeV") - print(f"\t Degrees of Freedom ([d1,d2],Cut(>25)): {selection.nDOF(signal)}") - print(f"\t Reduced Chi^2 ([d1,d2],Cut(<5): {selection.chi2nDOF(signal)}") - print(f"\t DOCA: {selection.DOCA(signal)} cm") - - print( - f"\t Preselection Cut passed:{selection.preselection_cut(signal,IP_cut=250)}" + event_weight = event.MCTrack[2].GetWeight() + + hist_dict["event_weight"].Fill(event_weight) + + for candidate_id_in_event, signal in enumerate(event.Particles): + hist_dict["candidate_time"].Fill(selection.define_candidate_time(signal)) + hist_dict["impact_parameter"].Fill(selection.impact_parameter(signal)) + hist_dict["dist_to_innerwall"].Fill(selection.dist_to_innerwall(signal)) + hist_dict["dist_to_vesselentrance"].Fill( + selection.dist_to_vesselentrance(signal) + ) + hist_dict["DOCA"].Fill(selection.DOCA(signal)) + hist_dict["inv_mass"].Fill(selection.invariant_mass(signal)) + hist_dict["len_Particles"].Fill(len(tree.Particles)) + hist_dict["d_mom"].Fill(*selection.daughtermomentum(signal)) + hist_dict["nDOF"].Fill(*selection.nDOF(signal)) + hist_dict["chi2nDOF"].Fill(*selection.chi2nDOF(signal)) + + IP_cut = 250 + preselection_flag = selection.preselection_cut( + signal, IP_cut=IP_cut, show_table=True ) - print("------------------------------------------------------------") + + if preselection_flag: + print( + f"Event:{event_nr} Candidate_index: {candidate_id_in_event} <--passes the pre-selection\n\n" + ) + else: + print(f"Event:{event_nr} Candidate_index: {candidate_id_in_event} \n\n") + + +ut.writeHists(hist_dict, "preselectionparameters.root") diff --git a/python/analysis_toolkit.py b/python/analysis_toolkit.py index e63de8d15..894f7e5e9 100644 --- a/python/analysis_toolkit.py +++ b/python/analysis_toolkit.py @@ -7,6 +7,7 @@ import yaml from rootpyPickler import Unpickler from ShipGeoConfig import AttrDict +from tabulate import tabulate class selection_check: @@ -63,7 +64,7 @@ def define_candidate_time(self, candidate): (candidate_pos.X() - hit.GetX()) ** 2 + (candidate_pos.Y() - hit.GetY()) ** 2 + (candidate_pos.Z() - hit.GetZ()) ** 2 - ) # distance to the vertex #in cm + ) # distance to the vertex in cm d_mom = self.tree.MCTrack[hit.GetTrackID()].GetP() / u.GeV mass = self.tree.MCTrack[hit.GetTrackID()].GetMass() @@ -138,7 +139,7 @@ def dist_to_innerwall(self, candidate): return min_distance if min_distance < 100 * u.m else 0 - def dist_to_entrance_vessel(self, candidate): + def dist_to_vesselentrance(self, candidate): """Calculate the distance of the candidate decay vertex to the entrance of the decay vessel.""" candidate_pos = ROOT.TVector3() candidate.GetVertex(candidate_pos) @@ -146,31 +147,27 @@ def dist_to_entrance_vessel(self, candidate): def nDOF(self, candidate): """Return the number of degrees of freedom (nDOF) for the particle's daughter tracks.""" - flag = True nmeas = [] t1, t2 = candidate.GetDaughter(0), candidate.GetDaughter(1) - ndf_cut = 25 # cut on nDOF + for tr in [t1, t2]: fit_status = self.tree.FitTracks[tr].getFitStatus() nmeas.append( int(round(fit_status.getNdf())) ) # nmeas.append(fit_status.getNdf()) - if nmeas[-1] <= ndf_cut: - flag = False # dof<=25=False - return nmeas, flag + + return np.array(nmeas) def daughtermomentum(self, candidate): """Return the momentum(Mag) of the particle's daughter tracks.""" - flag = True daughter_mom = [] t1, t2 = candidate.GetDaughter(0), candidate.GetDaughter(1) for trD in [t1, t2]: x = self.tree.FitTracks[trD] xx = x.getFittedState() daughter_mom.append((xx.getMom().Mag())) - if daughter_mom[-1] < 1 * u.GeV: - flag = False - return daughter_mom, flag + + return np.array(daughter_mom) def invariant_mass(self, candidate): """Invariant mass of the candidate.""" @@ -203,39 +200,120 @@ def is_in_fiducial(self, candidate): def chi2nDOF(self, candidate): """Return the reduced chi^2 of the particle's daughter tracks.""" t1, t2 = candidate.GetDaughter(0), candidate.GetDaughter(1) - # cut on chi2nDOF - flag = True - chi2ndf_cut = 5 + chi2ndf = [] for tr in [t1, t2]: fit_status = self.tree.FitTracks[tr].getFitStatus() chi2ndf.append(fit_status.getChi2() / fit_status.getNdf()) - if chi2ndf[-1] >= chi2ndf_cut: - flag = False - return chi2ndf, flag - def preselection_cut(self, candidate, IP_cut=250): + return np.array(chi2ndf) + + def preselection_cut(self, candidate, IP_cut=250, show_table=True): """Umbrella method to apply the pre-selection cuts on the candidate.""" + flag = True + if len(self.tree.Particles) != 1: - return False - if self.dist_to_innerwall(candidate) <= 5 * u.cm: - return False + flag = False if not (self.is_in_fiducial(candidate)): - return False - if self.dist_to_entrance_vessel(candidate) < 100 * u.cm: - return False - if self.impact_parameter(candidate) > IP_cut: - return False - nmeas, flag = self.nDOF(candidate) - if not (flag): - return False - chi2ndf, flag = self.chi2nDOF(candidate) - if not (flag): - return False - d_mom, flag = self.daughtermomentum(candidate) - if not (flag): - return False + flag = False + if self.dist_to_innerwall(candidate) <= 5 * u.cm: + flag = False + if self.dist_to_vesselentrance(candidate) <= 100 * u.cm: + flag = False + if self.impact_parameter(candidate) >= IP_cut * u.cm: + flag = False if self.DOCA(candidate) >= 1 * u.cm: - return False - - return True + flag = False + if np.any(self.nDOF(candidate) <= 25 * u.cm): + flag = False + if np.any(self.chi2nDOF(candidate) >= 5 * u.cm): + flag = False + if np.any(self.daughtermomentum(candidate) <= 1 * u.GeV): + flag = False + + if show_table: + table = [ + [ + "Number of candidates in event", + len(self.tree.Particles), + "==1", + len(self.tree.Particles) == 1, + ], + [ + "Time @ decay vertex (ns)", + self.define_candidate_time(candidate), + "", + "", + ], + [ + "Impact Parameter (cm)", + self.impact_parameter(candidate), + f"IP < {IP_cut*u.cm} cm", + self.impact_parameter(candidate) < IP_cut * u.cm, + ], + [ + "DOCA (cm)", + self.DOCA(candidate), + "DOCA < 1 cm", + self.DOCA(candidate) < 1 * u.cm, + ], + [ + "Is within Fiducial Volume?", + self.is_in_fiducial(candidate), + "True", + self.is_in_fiducial(candidate), + ], + [ + "Dist2InnerWall (cm)", + self.dist_to_innerwall(candidate), + "> 5 cm", + self.dist_to_innerwall(candidate) > 5 * u.cm, + ], + [ + "Dist2VesselEntrance (cm)", + self.dist_to_vesselentrance(candidate), + "> 100 cm", + self.dist_to_vesselentrance(candidate) > 100 * u.cm, + ], + ["Invariant Mass (GeV)", self.invariant_mass(candidate), "", ""], + [ + "Daughter Momentum [d1, d2] (GeV)", + self.daughtermomentum(candidate), + "> 1 GeV", + np.all(self.daughtermomentum(candidate) > 1 * u.GeV), + ], + [ + "Degrees of Freedom [d1, d2]", + self.nDOF(candidate), + "> 25", + np.all(self.nDOF(candidate) > 25), + ], + [ + "Reduced Chi^2 [d1, d2]", + self.chi2nDOF(candidate), + "< 5", + np.all(self.chi2nDOF(candidate) < 5), + ], + ["\033[1mPre-selection passed:\033[0m", "", "", flag], + ] + + for row in table: + row[3] = ( + f"\033[1;32m{row[3]}\033[0m" + if row[3] + else f"\033[1;31m{row[3]}\033[0m" + ) # Green for True, Red for False + + print( + tabulate( + table, + headers=[ + "Parameter", + "Value", + "Pre-selection cut", + "Pre-selection Check", + ], + tablefmt="grid", + ) + ) + return flag From 46e3b26453ffbd3135b82caae333de05f7fc36b1 Mon Sep 17 00:00:00 2001 From: anupama-reghunath Date: Thu, 5 Dec 2024 11:07:47 +0100 Subject: [PATCH 4/7] cleanup --- macro/analysis_example.py | 174 +++++++++++++++++++------------------ python/analysis_toolkit.py | 11 ++- 2 files changed, 97 insertions(+), 88 deletions(-) diff --git a/macro/analysis_example.py b/macro/analysis_example.py index 1499e859b..af4a7a889 100644 --- a/macro/analysis_example.py +++ b/macro/analysis_example.py @@ -12,95 +12,101 @@ options = parser.parse_args() -geo_file = options.path + "/geofile_full.conical.Pythia8-TGeant4.root" -input_file = options.path + "/ship.conical.Pythia8-TGeant4_rec.root" - -f = ROOT.TFile.Open(input_file, "read") +f = ROOT.TFile.Open(options.path + "/ship.conical.Pythia8-TGeant4_rec.root", "read") tree = f.cbmsim +geo_file = ROOT.TFile.Open( + options.path + "/geofile_full.conical.Pythia8-TGeant4.root", "read" +) + selection = analysis_toolkit.selection_check(geo_file) -hist_dict = {} -ut.bookHist(hist_dict, "event_weight", " Event weight; ; ", 100, 100, 100) -ut.bookHist( - hist_dict, - "candidate_time", - " candidate time @ decay vertex; ns ; ", - 200, - 100, - 100, -) -ut.bookHist(hist_dict, "impact_parameter", "Impact Parameter;cm;", 200, 100, 100) -ut.bookHist(hist_dict, "dist_to_innerwall", "Distance to Innerwall;cm;", 200, 100, 100) -ut.bookHist( - hist_dict, - "dist_to_vesselentrance", - "Distance to Decay Vessel Entrance;cm;", - 200, - 100, - 100, -) -ut.bookHist(hist_dict, "inv_mass", "Invariant mass ;GeV;", 200, 100, 100) -ut.bookHist(hist_dict, "DOCA", "Distance of Closest Approach ;cm;", 200, 100, 100) -ut.bookHist( - hist_dict, - "len_Particles", - " len(tree.Particles) ;Number of candidates per event;", - 200, - 100, - 100, -) -ut.bookHist( - hist_dict, - "d_mom", - "momentum of Daughters ;d1 (GeV);d2 (GeV)", - 200, - 100, - 100, - 200, - 100, - 100, -) -ut.bookHist(hist_dict, "nDOF", "nDOF ;d1;d2", 200, 100, 100, 200, 100, 100) -ut.bookHist(hist_dict, "chi2nDOF", "chi2nDOF ;d1;d2", 200, 100, 100, 200, 100, 100) - - -for event_nr, event in enumerate(tree): - selection.access_event(event) - - if len(event.Particles) == 0: - continue - - event_weight = event.MCTrack[2].GetWeight() - - hist_dict["event_weight"].Fill(event_weight) - - for candidate_id_in_event, signal in enumerate(event.Particles): - hist_dict["candidate_time"].Fill(selection.define_candidate_time(signal)) - hist_dict["impact_parameter"].Fill(selection.impact_parameter(signal)) - hist_dict["dist_to_innerwall"].Fill(selection.dist_to_innerwall(signal)) - hist_dict["dist_to_vesselentrance"].Fill( - selection.dist_to_vesselentrance(signal) - ) - hist_dict["DOCA"].Fill(selection.DOCA(signal)) - hist_dict["inv_mass"].Fill(selection.invariant_mass(signal)) - hist_dict["len_Particles"].Fill(len(tree.Particles)) - hist_dict["d_mom"].Fill(*selection.daughtermomentum(signal)) - hist_dict["nDOF"].Fill(*selection.nDOF(signal)) - hist_dict["chi2nDOF"].Fill(*selection.chi2nDOF(signal)) - - IP_cut = 250 - preselection_flag = selection.preselection_cut( - signal, IP_cut=IP_cut, show_table=True - ) - - if preselection_flag: - print( - f"Event:{event_nr} Candidate_index: {candidate_id_in_event} <--passes the pre-selection\n\n" +def Main_function(): + """Sample function to analyse the pre-selection parameters.""" + hist_dict = {} + + ut.bookHist(hist_dict, "event_weight", " Event weight; ; ", 100, 100, 100) + ut.bookHist( + hist_dict, + "candidate_time", + " candidate time @ decay vertex; ns ; ", + 200, + 100, + 100, + ) + ut.bookHist(hist_dict, "impact_parameter", "Impact Parameter;cm;", 200, 100, 100) + ut.bookHist( + hist_dict, "dist_to_innerwall", "Distance to Innerwall;cm;", 200, 100, 100 + ) + ut.bookHist( + hist_dict, + "dist_to_vesselentrance", + "Distance to Decay Vessel Entrance;cm;", + 200, + 100, + 100, + ) + ut.bookHist(hist_dict, "inv_mass", "Invariant mass ;GeV;", 200, 100, 100) + ut.bookHist(hist_dict, "DOCA", "Distance of Closest Approach ;cm;", 200, 100, 100) + ut.bookHist( + hist_dict, + "len_Particles", + " len(tree.Particles) ;Number of candidates per event;", + 200, + 100, + 100, + ) + ut.bookHist( + hist_dict, + "d_mom", + "momentum of Daughters ;d1 (GeV);d2 (GeV)", + 200, + 100, + 100, + 200, + 100, + 100, + ) + ut.bookHist(hist_dict, "nDOF", "nDOF ;d1;d2", 200, 100, 100, 200, 100, 100) + ut.bookHist(hist_dict, "chi2nDOF", "chi2nDOF ;d1;d2", 200, 100, 100, 200, 100, 100) + + for event_nr, event in enumerate(tree): + selection.access_event(event) + if len(event.Particles) == 0: + continue + + event_weight = event.MCTrack[2].GetWeight() + + hist_dict["event_weight"].Fill(event_weight) + + for candidate_id_in_event, signal in enumerate(event.Particles): + hist_dict["candidate_time"].Fill(selection.define_candidate_time(signal)) + hist_dict["impact_parameter"].Fill(selection.impact_parameter(signal)) + hist_dict["dist_to_innerwall"].Fill(selection.dist_to_innerwall(signal)) + hist_dict["dist_to_vesselentrance"].Fill( + selection.dist_to_vesselentrance(signal) + ) + hist_dict["DOCA"].Fill(selection.DOCA(signal)) + hist_dict["inv_mass"].Fill(selection.invariant_mass(signal)) + hist_dict["len_Particles"].Fill(len(tree.Particles)) + hist_dict["d_mom"].Fill(*selection.daughtermomentum(signal)) + hist_dict["nDOF"].Fill(*selection.nDOF(signal)) + hist_dict["chi2nDOF"].Fill(*selection.chi2nDOF(signal)) + + IP_cut = 250 + preselection_flag = selection.preselection_cut( + signal, IP_cut=IP_cut, show_table=False ) - else: - print(f"Event:{event_nr} Candidate_index: {candidate_id_in_event} \n\n") + + if preselection_flag: + print( + f"Event:{event_nr} Candidate_index: {candidate_id_in_event} <--passes the pre-selection\n\n" + ) + else: + print(f"Event:{event_nr} Candidate_index: {candidate_id_in_event} \n\n") + + ut.writeHists(hist_dict, "preselectionparameters.root") -ut.writeHists(hist_dict, "preselectionparameters.root") +Main_function() diff --git a/python/analysis_toolkit.py b/python/analysis_toolkit.py index 894f7e5e9..df99594c6 100644 --- a/python/analysis_toolkit.py +++ b/python/analysis_toolkit.py @@ -13,9 +13,8 @@ class selection_check: """Class to perform various selection checks on the candidate.""" - def __init__(self, g): + def __init__(self, geo_file): """Initialize the selection_check class with geometry and configuration.""" - geo_file = ROOT.TFile.Open(g, "read") self.geometry_manager = geo_file.FAIRGeom unpickler = Unpickler(geo_file) self.ship_geo = unpickler.load("ShipGeo") @@ -208,8 +207,12 @@ def chi2nDOF(self, candidate): return np.array(chi2ndf) - def preselection_cut(self, candidate, IP_cut=250, show_table=True): - """Umbrella method to apply the pre-selection cuts on the candidate.""" + def preselection_cut(self, candidate, IP_cut=250, show_table=False): + """ + Umbrella method to apply the pre-selection cuts on the candidate. + + show_table=True tabulates the pre-selection parameters. + """ flag = True if len(self.tree.Particles) != 1: From 715979daa90ae760f100348ac82388c59e217b2e Mon Sep 17 00:00:00 2001 From: anupama-reghunath Date: Wed, 11 Dec 2024 22:43:46 +0100 Subject: [PATCH 5/7] change main function name --- macro/analysis_example.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/macro/analysis_example.py b/macro/analysis_example.py index af4a7a889..ee427635f 100644 --- a/macro/analysis_example.py +++ b/macro/analysis_example.py @@ -22,7 +22,7 @@ selection = analysis_toolkit.selection_check(geo_file) -def Main_function(): +def main(): """Sample function to analyse the pre-selection parameters.""" hist_dict = {} @@ -109,4 +109,5 @@ def Main_function(): ut.writeHists(hist_dict, "preselectionparameters.root") -Main_function() +if __name__ == "__main__": + main() From 6af146c1111014ef7e45caf098a5bafbc5da8c83 Mon Sep 17 00:00:00 2001 From: anupama-reghunath Date: Tue, 17 Dec 2024 23:50:54 +0100 Subject: [PATCH 6/7] Shift scripts to experimental / example directory --- {macro => examples}/analysis_example.py | 2 +- python/experimental/__init__.py | 5 +++++ python/{ => experimental}/analysis_toolkit.py | 5 ++--- 3 files changed, 8 insertions(+), 4 deletions(-) rename {macro => examples}/analysis_example.py (98%) create mode 100644 python/experimental/__init__.py rename python/{ => experimental}/analysis_toolkit.py (98%) diff --git a/macro/analysis_example.py b/examples/analysis_example.py similarity index 98% rename from macro/analysis_example.py rename to examples/analysis_example.py index ee427635f..3929c9751 100644 --- a/macro/analysis_example.py +++ b/examples/analysis_example.py @@ -3,9 +3,9 @@ from argparse import ArgumentParser -import analysis_toolkit import ROOT import rootUtils as ut +from experimental import analysis_toolkit parser = ArgumentParser(description=__doc__) parser.add_argument("--path", help="Path to simulation file", default="./") diff --git a/python/experimental/__init__.py b/python/experimental/__init__.py new file mode 100644 index 000000000..4821959d6 --- /dev/null +++ b/python/experimental/__init__.py @@ -0,0 +1,5 @@ +""" +This is the 'experimental' package, which contains scripts that are still being tested and may undergo significant changes. + +Suggestions & feedback are highly appreciated as this is a work-in-progress. +""" diff --git a/python/analysis_toolkit.py b/python/experimental/analysis_toolkit.py similarity index 98% rename from python/analysis_toolkit.py rename to python/experimental/analysis_toolkit.py index df99594c6..436d8ca37 100644 --- a/python/analysis_toolkit.py +++ b/python/experimental/analysis_toolkit.py @@ -30,7 +30,6 @@ def __init__(self, geo_file): with open(fairship + "/geometry/veto_config_vacuums.yaml", "r") as file: config = yaml.safe_load(file) self.veto_geo = AttrDict(config) - self.veto_geo.z0 def access_event(self, tree): """Access event data.""" @@ -227,9 +226,9 @@ def preselection_cut(self, candidate, IP_cut=250, show_table=False): flag = False if self.DOCA(candidate) >= 1 * u.cm: flag = False - if np.any(self.nDOF(candidate) <= 25 * u.cm): + if np.any(self.nDOF(candidate) <= 25): flag = False - if np.any(self.chi2nDOF(candidate) >= 5 * u.cm): + if np.any(self.chi2nDOF(candidate) >= 5): flag = False if np.any(self.daughtermomentum(candidate) <= 1 * u.GeV): flag = False From f5aa6c0c386dadc1f7abe2fe1a70822225577518 Mon Sep 17 00:00:00 2001 From: anupama-reghunath Date: Wed, 18 Dec 2024 09:09:28 +0100 Subject: [PATCH 7/7] Add ChangeLOG entry for analysis_toolkit prototype --- CHANGELOG.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index b5daa4285..3deab4014 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -14,6 +14,9 @@ it in future. ### Added +* New analysis toolkit prototype added as part of the 'experimental' package. +* Simple analysis example script now available in 'examples/' + ### Fixed * Use ConstructedAt instead of remove pythonization for TClonesArray