diff --git a/shipLHC/mc_create_root.py b/shipLHC/mc_create_root.py new file mode 100644 index 0000000000..097ae9104e --- /dev/null +++ b/shipLHC/mc_create_root.py @@ -0,0 +1,56 @@ +#!/usr/bin/env python +import ROOT +import argparse + +ROOT.gROOT.SetBatch(True) + +# 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 +track_chain.Add(input_file) # Add the specified input file to the TChain + +# Prepare output ROOT file and TTree +output_file = ROOT.TFile(output_file_name, "RECREATE") +ntuple = ROOT.TNtuple("muon_tracks", "Extrapolated Muon Tracks", + "event:x:y:z:theta:phi") + +# 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 + + # Extract position and momentum + pos = track.getStart() # pos.x(), pos.y(), pos.z() + mom = track.getTrackMom() # mom.x(), mom.y(), mom.z() + + # 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() + + # 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, z_plane, theta, phi) + +# Write to the output file and close it +output_file.Write() +output_file.Close() +print(f"Extraction completed. Results are saved in '{output_file_name}'. Used z-plane: {z_plane}") diff --git a/shipLHC/run_simSND.py b/shipLHC/run_simSND.py index 96249d5329..49055d0e29 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,32 @@ 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() + 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 + z_plane = ntuple.z + theta = ntuple.theta + phi = ntuple.phi + pgun = ROOT.FairBoxGenerator(13, 1) #Particle ID 13 for muons, charge is irrelevent + 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) + ROOT.FairLogger.GetLogger().SetLogScreenLevel("WARNING") # -----muon DIS Background------------------------ if simEngine == "muonDIS": ut.checkFileExists(inputFile)