Changes in run_simSND

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}")

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)