From ce2875de18f352d4d4bd0855e559ed10fa8b0a8a Mon Sep 17 00:00:00 2001 From: Anushka Durg Date: Tue, 29 Oct 2024 13:59:18 +0100 Subject: [PATCH 1/5] added Ntuple file as input for runsim_SND Test1 dataprofile --- shipLHC/run_simSND.py | 29 ++++++++++++++++++++++++++++- 1 file changed, 28 insertions(+), 1 deletion(-) diff --git a/shipLHC/run_simSND.py b/shipLHC/run_simSND.py index 96249d5329..58348c414b 100644 --- a/shipLHC/run_simSND.py +++ b/shipLHC/run_simSND.py @@ -59,7 +59,7 @@ parser.add_argument("--EmuDet","--nuTargetActive",dest="nuTargetPassive",help="activate emulsiondetector", required=False,action="store_false") parser.add_argument("--NagoyaEmu","--useNagoyaEmulsions",dest="useNagoyaEmulsions",help="use bricks of 57 Nagoya emulsion films instead of 60 Slavich", required=False,action="store_true") parser.add_argument("-y", dest="year", help="specify the year to generate the respective TI18 detector setup", required=False, type=int, default=2024) - +parser.add_argument("--muonDataProfile", help="Path to the ROOT file with muon track data") options = parser.parse_args() # user hook @@ -166,6 +166,7 @@ def Exec(self,opt): primGen = ROOT.FairPrimaryGenerator() # -----Particle Gun----------------------- + if simEngine == "PG": myPgun = ROOT.FairBoxGenerator(options.pID,1) myPgun.SetPRange(options.Estart,options.Eend) @@ -191,7 +192,33 @@ def Exec(self,opt): f'({options.EVx},{options.EVy},{options.EVz})[cm × cm × cm] \n' f'with a uniform x-y spread of (Dx,Dy)=({options.Dx},{options.Dy})[cm × cm]' f' and {options.nZSlices} z slices in steps of {options.zSliceStep}[cm].') + ROOT.FairLogger.GetLogger().SetLogScreenLevel("WARNING") # otherwise stupid printout for each event +#________NTUPLE-------- +if options.muonDataProfile: + print(f"Using file: {options.muonDataProfile}") + file = ROOT.TFile(options.muonDataProfile) + ntuple = file.Get("muon_tracks") + primGen = ROOT.FairPrimaryGenerator() + z_plane = 250.0 * u.cm + pgun = ROOT.FairBoxGenerator(13, 1) # Particle ID 13 for muons, charge is irrelevant + pgun.SetPRange(30, 3500) # Set momentum range + n_entries = ntuple.GetEntries() # Number of entries to cycle through + if options.nEvents > n_entries: + print("Warning: The number of events requested exceeds the number of entries in the data profile. Cycling through entries.") + for i in range(options.nEvents): + ntuple.GetEntry(i % n_entries) # Cycle through entries + x_plane = ntuple.x + y_plane = ntuple.y + slope_xz = ntuple.slopes_xz + slope_yz = ntuple.slopes_yz + x0 = x_plane - slope_xz * z_plane + y0 = y_plane - slope_yz * z_plane + z0 = 0 + pgun.SetXYZ(x0 * u.cm, y0 * u.cm, z0 * u.cm) # Set the coordinates of each muon + primGen.AddGenerator(pgun) + run.SetGenerator(primGen) + ROOT.FairLogger.GetLogger().SetLogScreenLevel("WARNING") # -----muon DIS Background------------------------ if simEngine == "muonDIS": ut.checkFileExists(inputFile) From cbdb7c013ecc3e54c6fb4425144ddc5cb50698c0 Mon Sep 17 00:00:00 2001 From: Anushka Durg Date: Tue, 29 Oct 2024 13:59:18 +0100 Subject: [PATCH 2/5] added Ntuple file as input for runsim_SND Test1 dataprofile --- shipLHC/run_simSND.py | 33 ++++++++++++++++++++++++++++++--- 1 file changed, 30 insertions(+), 3 deletions(-) diff --git a/shipLHC/run_simSND.py b/shipLHC/run_simSND.py index 96249d5329..6a2f87ecf3 100644 --- a/shipLHC/run_simSND.py +++ b/shipLHC/run_simSND.py @@ -59,7 +59,7 @@ parser.add_argument("--EmuDet","--nuTargetActive",dest="nuTargetPassive",help="activate emulsiondetector", required=False,action="store_false") parser.add_argument("--NagoyaEmu","--useNagoyaEmulsions",dest="useNagoyaEmulsions",help="use bricks of 57 Nagoya emulsion films instead of 60 Slavich", required=False,action="store_true") parser.add_argument("-y", dest="year", help="specify the year to generate the respective TI18 detector setup", required=False, type=int, default=2024) - +parser.add_argument("--muonDataProfile", help="Path to the ROOT file with muon track data") options = parser.parse_args() # user hook @@ -85,7 +85,7 @@ def Exec(self,opt): if options.muonback: simEngine = "MuonBack" if options.nuage: simEngine = "Nuage" if options.mudis: simEngine = "muonDIS" - +if options.muonDataProfile: simEngine = "muonDataProfile" if options.inputFile: if options.inputFile == "none": options.inputFile = None inputFile = options.inputFile @@ -152,7 +152,7 @@ def Exec(self,opt): run.SetName(mcEngine) # Transport engine run.SetSink(ROOT.FairRootFileSink(outFile)) # Output file run.SetUserConfig("g4Config.C") # user configuration file default g4Config.C -rtdb = run.GetRuntimeDb() +rtdb = run.GetRuntimeDb() # add user task if userTask: userTask = MyTask() @@ -166,6 +166,7 @@ def Exec(self,opt): primGen = ROOT.FairPrimaryGenerator() # -----Particle Gun----------------------- + if simEngine == "PG": myPgun = ROOT.FairBoxGenerator(options.pID,1) myPgun.SetPRange(options.Estart,options.Eend) @@ -191,7 +192,33 @@ def Exec(self,opt): f'({options.EVx},{options.EVy},{options.EVz})[cm × cm × cm] \n' f'with a uniform x-y spread of (Dx,Dy)=({options.Dx},{options.Dy})[cm × cm]' f' and {options.nZSlices} z slices in steps of {options.zSliceStep}[cm].') + ROOT.FairLogger.GetLogger().SetLogScreenLevel("WARNING") # otherwise stupid printout for each event +#------NTUPLE------ +if options.muonDataProfile: + print(f"Using file: {options.muonDataProfile}") + file = ROOT.TFile(options.muonDataProfile) + ntuple = file.Get("muon_tracks") + primGen = ROOT.FairPrimaryGenerator() + z_plane = 250.0 * u.cm + n_entries = ntuple.GetEntries() # Number of entries to cycle through + if options.nEvents > n_entries: + print("Warning: The number of events requested exceeds the number of entries in the data profile. Cycling through entries.") + for i in range(options.nEvents): + ntuple.GetEntry(i % n_entries) # Cycle through entries if needed + x_plane = ntuple.x + y_plane = ntuple.y + slope_xz = ntuple.slopes_xz + slope_yz = ntuple.slopes_yz + x0 = x_plane - slope_xz * z_plane + y0 = y_plane - slope_yz * z_plane + z0 = 0 + pgun = ROOT.FairBoxGenerator(13, 1) #Particle ID 13 for muons, charge is irrelevent + pgun.SetPRange(30, 3500) # Set momentum range + pgun.SetXYZ(x0 * u.cm, y0 * u.cm, z0 * u.cm) # Set the coordinates of each muon + primGen.AddGenerator(pgun) + run.SetGenerator(primGen) + ROOT.FairLogger.GetLogger().SetLogScreenLevel("WARNING") # -----muon DIS Background------------------------ if simEngine == "muonDIS": ut.checkFileExists(inputFile) From 9c28ced0b20daec68eba28cda1a32b3e99459617 Mon Sep 17 00:00:00 2001 From: Anushka Date: Thu, 14 Nov 2024 17:31:50 +0100 Subject: [PATCH 3/5] file to help generate a root file for MC --- shipLHC/mc_create_root.py | 54 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 54 insertions(+) create mode 100644 shipLHC/mc_create_root.py diff --git a/shipLHC/mc_create_root.py b/shipLHC/mc_create_root.py new file mode 100644 index 0000000000..8c11f0b7f5 --- /dev/null +++ b/shipLHC/mc_create_root.py @@ -0,0 +1,54 @@ +#!/usr/bin/env python +import ROOT +import numpy as np +import os + +ROOT.gROOT.SetBatch(True) + +# Define the fixed z coordinate where we will extract x, y positions +z_plane = 250.0 # adjust as needed + +# Set up the TChain for the ROOT files +track_chain = ROOT.TChain("rawConv") # Ensure the TTree name matches + +# Directory containing ROOT files +root_files_dir = "/eos/experiment/sndlhc/users/sii/4705/" +file_pattern = "sndsw_raw-{}_{}_4705_muonReco.root" + +# Adding all ROOT files to the TChain +for part in range(10, 15): + for it in range(10): + file_path = os.path.join(root_files_dir, file_pattern.format(it, part)) + track_chain.Add(file_path) + +# Prepare output ROOT file and TTree +output_file = ROOT.TFile("extracted_muon_tracks.root", "RECREATE") +ntuple = ROOT.TNtuple("muon_tracks", "Extrapolated Muon Tracks", + "event:x:y:slopes_xz:slopes_yz") + +# Loop over each event and track +for i_event, event in enumerate(track_chain): + for track in event.Reco_MuonTracks: + if track.getTrackType() != 11 or track.getTrackFlag() != 1: + continue # Skip non-Scifi tracks or bad quality tracks + + # Extrapolate position and calculate slopes + pos = track.getStart() + mom = track.getTrackMom() + + # Calculate lamda for z-plane intersection + lam = (z_plane - pos.z()) / mom.z() + extrap_x = pos.x() + lam * mom.x() + extrap_y = pos.y() + lam * mom.y() + + # Slopes + slope_xz = mom.x() / mom.z() + slope_yz = mom.y() / mom.z() + + # Store in the ntuple + ntuple.Fill(i_event, extrap_x, extrap_y, slope_xz, slope_yz) + +# Write to the output file and close it +output_file.Write() +output_file.Close() +print("Extraction completed. Results are saved in 'extracted_muon_tracks.root'.") From ae0ed9b1e160a37057c4de4291b48bd0368697e2 Mon Sep 17 00:00:00 2001 From: Anushka Date: Fri, 15 Nov 2024 16:22:03 +0100 Subject: [PATCH 4/5] fixup! added Ntuple file as input for runsim_SND --- shipLHC/run_simSND.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/shipLHC/run_simSND.py b/shipLHC/run_simSND.py index 6a2f87ecf3..74a353b954 100644 --- a/shipLHC/run_simSND.py +++ b/shipLHC/run_simSND.py @@ -200,7 +200,6 @@ def Exec(self,opt): file = ROOT.TFile(options.muonDataProfile) ntuple = file.Get("muon_tracks") primGen = ROOT.FairPrimaryGenerator() - z_plane = 250.0 * u.cm n_entries = ntuple.GetEntries() # Number of entries to cycle through if options.nEvents > n_entries: print("Warning: The number of events requested exceeds the number of entries in the data profile. Cycling through entries.") @@ -208,16 +207,17 @@ def Exec(self,opt): ntuple.GetEntry(i % n_entries) # Cycle through entries if needed x_plane = ntuple.x y_plane = ntuple.y - slope_xz = ntuple.slopes_xz - slope_yz = ntuple.slopes_yz - x0 = x_plane - slope_xz * z_plane - y0 = y_plane - slope_yz * z_plane - z0 = 0 + z_plane = ntuple.z + theta = ntuple.theta + phi = ntuple.phi pgun = ROOT.FairBoxGenerator(13, 1) #Particle ID 13 for muons, charge is irrelevent - pgun.SetPRange(30, 3500) # Set momentum range - pgun.SetXYZ(x0 * u.cm, y0 * u.cm, z0 * u.cm) # Set the coordinates of each muon + pgun.SetThetaRange(theta, theta) + pgun.SetPhiRange(phi,phi) + pgun.SetMomentumRange(30,35000) + pgun.SetXYZ(x_plane * u.cm, y_plane * u.cm, z_plane * u.cm) # Set the coordinates of each muon primGen.AddGenerator(pgun) run.SetGenerator(primGen) + run.SetGenerator(primGen) ROOT.FairLogger.GetLogger().SetLogScreenLevel("WARNING") # -----muon DIS Background------------------------ if simEngine == "muonDIS": From c87a534916f95c8ebcb25f3678a29f60212660b9 Mon Sep 17 00:00:00 2001 From: Anushka Date: Fri, 15 Nov 2024 16:27:14 +0100 Subject: [PATCH 5/5] fixup! file to help generate a root file for MC --- shipLHC/mc_create_root.py | 56 ++++++++++++++++++++------------------- 1 file changed, 29 insertions(+), 27 deletions(-) diff --git a/shipLHC/mc_create_root.py b/shipLHC/mc_create_root.py index 8c11f0b7f5..097ae9104e 100644 --- a/shipLHC/mc_create_root.py +++ b/shipLHC/mc_create_root.py @@ -1,30 +1,32 @@ #!/usr/bin/env python import ROOT -import numpy as np -import os +import argparse ROOT.gROOT.SetBatch(True) -# Define the fixed z coordinate where we will extract x, y positions -z_plane = 250.0 # adjust as needed - -# Set up the TChain for the ROOT files +# Parse command-line arguments +parser = argparse.ArgumentParser(description="Extract muon track data at a specific z-plane.") +parser.add_argument("--zplane", type=float, required=True, + help="Fixed z-plane value for extrapolation.") +parser.add_argument("--inputfile", type=str, required=True, + help="Input ROOT file containing muon tracks.") +parser.add_argument("--outputfile", type=str, default="extracted_muon_tracks.root", + help="Output ROOT file to save the extracted data.") +args = parser.parse_args() + +# Get the z-plane value and file paths from user input +z_plane = args.zplane +input_file = args.inputfile +output_file_name = args.outputfile + +# Set up the TChain for the ROOT file track_chain = ROOT.TChain("rawConv") # Ensure the TTree name matches - -# Directory containing ROOT files -root_files_dir = "/eos/experiment/sndlhc/users/sii/4705/" -file_pattern = "sndsw_raw-{}_{}_4705_muonReco.root" - -# Adding all ROOT files to the TChain -for part in range(10, 15): - for it in range(10): - file_path = os.path.join(root_files_dir, file_pattern.format(it, part)) - track_chain.Add(file_path) +track_chain.Add(input_file) # Add the specified input file to the TChain # Prepare output ROOT file and TTree -output_file = ROOT.TFile("extracted_muon_tracks.root", "RECREATE") +output_file = ROOT.TFile(output_file_name, "RECREATE") ntuple = ROOT.TNtuple("muon_tracks", "Extrapolated Muon Tracks", - "event:x:y:slopes_xz:slopes_yz") + "event:x:y:z:theta:phi") # Loop over each event and track for i_event, event in enumerate(track_chain): @@ -32,23 +34,23 @@ if track.getTrackType() != 11 or track.getTrackFlag() != 1: continue # Skip non-Scifi tracks or bad quality tracks - # Extrapolate position and calculate slopes - pos = track.getStart() - mom = track.getTrackMom() + # Extract position and momentum + pos = track.getStart() # pos.x(), pos.y(), pos.z() + mom = track.getTrackMom() # mom.x(), mom.y(), mom.z() - # Calculate lamda for z-plane intersection + # Calculate lambda for z-plane intersection lam = (z_plane - pos.z()) / mom.z() extrap_x = pos.x() + lam * mom.x() extrap_y = pos.y() + lam * mom.y() - # Slopes - slope_xz = mom.x() / mom.z() - slope_yz = mom.y() / mom.z() + # Calculate theta and phi angles + theta = ROOT.TMath.ATan(ROOT.TMath.Sqrt((mom.x() ** 2 + mom.y() ** 2) / mom.z() ** 2)) + phi = ROOT.TMath.ATan2(mom.y(), mom.x()) # Store in the ntuple - ntuple.Fill(i_event, extrap_x, extrap_y, slope_xz, slope_yz) + ntuple.Fill(i_event, extrap_x, extrap_y, z_plane, theta, phi) # Write to the output file and close it output_file.Write() output_file.Close() -print("Extraction completed. Results are saved in 'extracted_muon_tracks.root'.") +print(f"Extraction completed. Results are saved in '{output_file_name}'. Used z-plane: {z_plane}")