Merge branch 'mcgen' into 'master'

Add DIF to MC production scripts

See merge request faser/calypso!325

Control/CalypsoExample/Generation/python/faser_parser.py

@@ -47,6 +47,12 @@ def faser_pgparser():
     parser.add_argument("--zpos", default=None, type=float,
                         help="Specify z position of particles (in mm) (helpful to avoid FASERnu)")
 
+    parser.add_argument("--pidd1", default=None, type=int, 
+                        help="Specify PDG ID of daugther 1 for DIF generator")
+    parser.add_argument("--pidd2", default=None, type=int, 
+                        help="Specify PDG ID of daugther 2 for DIF generator")
+
+
     parser.add_argument("--sampler", default="log",
                         help="Specify energy sampling (log, lin, const, hist, hist2D)")
     parser.add_argument("--hist_name", default="log",

Control/CalypsoExample/Generation/scripts/faser_particlegun.py

@@ -101,8 +101,11 @@ if __name__ == '__main__':
     from AthenaCommon.PhysicalConstants import pi
 
     if isinstance(args.pid, list):
-        # Note args.pid is a list, must make this a set for ParticleGun
-        pidarg = set(args.pid)
+        # Note args.pid is a list, must make this a set for ParticleGun if have > 1 
+        if len(args.pid) > 1:
+            pidarg = set(args.pid)
+        else:
+            pidarg = args.pid[0]
     else:
         # Just pass a single value
         pidarg = args.pid
@@ -112,13 +115,32 @@ if __name__ == '__main__':
     # Create the simgun dictionary
     # Negative radius gives uniform sampling
     # Positive radius gives Gaussian sampling
-    sg_dict = {
-        "Generator" : "SingleParticle", 
-        "pid" : pidarg, "mass" : args.mass, 
-        "theta" :  PG.GaussianSampler(0, args.angle, oneside = True), 
-        "phi" : [0, 2*pi], "radius" : args.radius, 
-        "randomSeed" : args.outfile 
-    }
+
+    if args.pidd1:
+        if args.pidd2 is None:
+            args.pidd2 = -args.pidd1
+            
+        sg_dict = {
+            "Generator" : "DecayInFlight", 
+            "mother_pid" : pidarg,
+            "daughter1_pid" : args.pidd1,
+            "daughter2_pid" : args.pidd2,        
+            "mass" : args.mass, 
+            "theta" :  PG.GaussianSampler(0, args.angle, oneside = True) if args.angle is not None and args.angle != "None" else None,
+            "phi" : [0, 2*pi], "radius" : args.radius, 
+            "randomSeed" : args.outfile 
+            }
+
+    else:
+
+        sg_dict = {
+            "Generator" : "SingleParticle", 
+            "pid" : pidarg, "mass" : args.mass, 
+            "theta" :  PG.GaussianSampler(0, args.angle, oneside = True) if args.angle is not None and args.angle != "None" else None,
+            "phi" : [0, 2*pi], "radius" : args.radius, 
+            "randomSeed" : args.outfile 
+            }
+
 
     # -1000 is safely upstream of detector (to be checked)
     # Note zpos is in mm!
@@ -132,6 +154,12 @@ if __name__ == '__main__':
         sg_dict["energy"] = PG.LogSampler(args.minE*GeV, args.maxE*GeV)
     elif args.sampler == "const":
         sg_dict["energy"] = PG.ConstSampler(args.maxE*GeV)
+    elif args.sampler == "hist":
+        fname, hname = args.hist_name.split(":")
+        sg_dict["energy"] = PG.TH1Sampler(fname, hname)
+    elif args.sampler == "hist2D":
+        fname, hname = args.hist_name.split(":")
+        sg_dict["energy"] = PG.TH2Sampler(fname, hname)        
     else:
         print(f"Sampler {args.sampler} not known!")
         sys.exit(1)

Generators/DIFGenerator/python/DIFSampler.py

@@ -93,9 +93,8 @@ class DIFSampler(PG.ParticleSampler):
     # def __init__(self, daughter1_pid = 11, daughter2_pid = -11, mother_pid = None, mother_mom = PG.EThetaMPhiSampler(sqrt((1*TeV)**2 + (10*MeV)**2),0,10*MeV,0),mother_pos = CylinderSampler([0, 100**2],[0, 2*pi],[-1500, 0],0)):
     def __init__(self, daughter1_pid=13, daughter2_pid=-13, mother_pid=None,
                  mother_mom=PG.EThetaMPhiSampler(sqrt((1*TeV)**2 + (500*MeV)**2), [0, 0.0002], 500*MeV, [0, 2*pi]),
-                 my_z_position=-1500):
-                 # mother_pos=CylinderSampler([0, 100**2], [0, 2*pi], [-1500, 0], 0)):
-        self._mother_sampler = PG.ParticleSampler(pid = mother_pid, mom = mother_mom, pos=CylinderSampler([0, 100**2], [0, 2*pi], my_z_position, 0))
+                 mother_pos=CylinderSampler([0, 100**2], [0, 2*pi], [-1500, 0], 0)):
+        self._mother_sampler = PG.ParticleSampler(pid = mother_pid, mom = mother_mom, pos=mother_pos)
 
         self.daughter1 = self.particle(daughter1_pid)
         self.daughter2 = self.particle(daughter2_pid)

Generators/FaserParticleGun/python/FaserParticleGunConfig.py

@@ -9,7 +9,7 @@ import ParticleGun as PG
 from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
 from AthenaConfiguration.ComponentFactory import CompFactory
 # from AthenaCommon.Constants import VERBOSE, INFO
-from AthenaCommon.SystemOfUnits import TeV
+from AthenaCommon.SystemOfUnits import TeV, GeV, MeV
 from AthenaCommon.PhysicalConstants import pi
 
 ### add radial pos sampler ### with gaussian beam implemented
@@ -144,7 +144,31 @@ def FaserParticleGunDecayInFlightCfg(ConfigFlags, **kwargs) :
     pg = cfg.getPrimary()
 
     from DIFGenerator import DIFSampler
-    pg.sampler = DIFSampler(**kwargs)
+
+    if "radius" in kwargs:
+        kwargs["mother_pos"]  =  RadialPosSampler(x = kwargs.setdefault("x", 0.0), 
+                                                  y = kwargs.setdefault("y", 0.0), 
+                                                  z = kwargs.setdefault("z", -3750.0),
+                                                  r = kwargs.setdefault("radius", 1.0),
+                                                  t = kwargs.setdefault("t", 0.0) )
+    else:
+        from DIFGenerator.DIFSampler import CylinderSampler
+        kwargs["mother_pos"]  = CylinderSampler([0, 100**2], [0, 2*pi], [-1500, 0], 0)
+    
+    if not "mother_mom" in kwargs:
+        kwargs["mother_mom"] = PG.EThetaMPhiSampler(kwargs.setdefault("energy", ((1*TeV)**2 + (500*MeV)**2)**0.5),
+                                                    kwargs.setdefault("theta", [0, 0.0002]),
+                                                    kwargs.setdefault("mass", 500*MeV),
+                                                    kwargs.setdefault("phi", [0, 2*pi]) )                                  
+        
+
+        
+    pg.sampler = DIFSampler(kwargs.setdefault("daughter1_pid", 13),
+                            kwargs.setdefault("daughter2_pid", -13),
+                            kwargs.setdefault("mother_pid", None),
+                            kwargs["mother_mom"],
+                            kwargs["mother_pos"] )
+
 
     return cfg
 
@@ -190,9 +214,12 @@ def FaserParticleGunForeseeCfg(ConfigFlags, **kwargs) :
     return cfg
 
 def FaserParticleGunCfg(ConfigFlags) :
-    # generator = ConfigFlags.Sim.Gun.setdefault("Generator", "SingleParticle")
-    generator = ConfigFlags.Sim.Gun.setdefault("Generator", "DecayInFlight")
+    generator = ConfigFlags.Sim.Gun.setdefault("Generator", "SingleParticle")
+    # generator = ConfigFlags.Sim.Gun.setdefault("Generator", "DecayInFlight")
+
     kwargs = ConfigFlags.Sim.Gun
+    del kwargs["Generator"]
+    
     if generator == "SingleEcalParticle" :
         return FaserParticleGunSingleEcalParticleCfg(ConfigFlags, **kwargs)
     elif generator == "Cosmics" :

Generators/ForeseeGenerator/share/generate_forsee_events.py

@@ -47,13 +47,13 @@ class ForeseeGenerator(object):
         else:
             self.foresee = Foresee(path = self.path)
 
-        # Generate 6 cm high to account for translation from ATLAS to FASER coord. system
+        # Generate 6.5 cm high to account for translation from ATLAS to FASER coord. system
         # TODO: relax this a bit as daughters may enter even if mother doesn't
 #         self.foresee.set_detector(selection="np.sqrt(x.x**2 + x.y**2)< 0.1",  
 #                                   channels=[self.mode], distance=480, length=1.5 ,
 #                                   luminosity=1/1000.) # 1 pb-1        
 
-        self.foresee.set_detector(selection="np.sqrt(x.x**2 + (x.y - 0.06)**2)< 0.1",  
+        self.foresee.set_detector(selection="np.sqrt(x.x**2 + (x.y - 0.065)**2)< 0.1",  
                                   channels=[self.mode], distance=480, length=1.5 ,
                                   luminosity=1/1000.) # 1 pb-1        
 

Generators/ForeseeGenerator/share/validate_grid.py

 import json
 import subprocess as proc
-import os
+from glob import glob
+import os, sys
 
 grid = "../calypso/Generators/ForeseeGenerator/share/points.json"
 
 with open(grid) as f:
     data = json.load(f)
 
-name = "DarkPhotonGrid13p6"
+name = "DarkPhotonGrid13p6_65mm"
 path = f"/bundle/data/FASER/Carl_Output/Foresee/{name}/"
 ecom = 13.6
 pid1 = 11
@@ -15,7 +16,12 @@ pid2 = -11
 
 for coup, masses in data["samples"].items():
     for m in masses:
-        infile = f"{path}events_{ecom}TeV_m{m}GeV_c{float(coup):.1e}to_{pid1}_{pid2}_{name}.hepmc"
+        files = glob(f"{path}events_{ecom}TeV_m{m}GeV_c{float(coup):.1e}to_{pid1}_{pid2}_{name}*.hepmc")
+
+        if len(files) != 1:
+            continue
+        
+        infile = files[0]
         outfile = infile.replace(".hepmc", ".EVNT.pool.root")
         valfile = infile.replace(".hepmc", ".VAL.pool.root")
         valdir = infile.replace(".hepmc", "")