Merge branch 'fasermc-dev' into 'master'

Consolodate s0010 script updates

See merge request faser/calypso!345

Control/CalypsoExample/Generation/python/faser_parser.py

@@ -44,6 +44,10 @@ def faser_pgparser():
                         help="Specify radius (in mm)")
     parser.add_argument("--angle", default=0.005, type=float_or_none,
                         help="Specify angular width (in Rad)")
+    parser.add_argument("--xpos", default=None, type=float,
+                        help="Specify x position of particles (in mm)")
+    parser.add_argument("--ypos", default=None, type=float,
+                        help="Specify y position of particles (in mm)")
     parser.add_argument("--zpos", default=None, type=float,
                         help="Specify z position of particles (in mm) (helpful to avoid FASERnu)")
 

Control/CalypsoExample/Generation/scripts/faser_particlegun.py

@@ -148,9 +148,15 @@ if __name__ == '__main__':
 
     # -1000 is safely upstream of detector (to be checked)
     # Note zpos is in mm!
-    if args.zpos:
+    if args.zpos != None:
         sg_dict["z"] = args.zpos
 
+    if args.xpos != None:
+        sg_dict["x"] = args.xpos
+
+    if args.ypos != None:
+        sg_dict["y"] = args.ypos
+
     # Determine energy sampling
     if args.sampler == "lin":
         sg_dict["energy"] = PG.UniformSampler(args.minE*GeV, args.maxE*GeV)
@@ -159,11 +165,28 @@ if __name__ == '__main__':
     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)
+        nargs = len(args.hist_name.split(":"))
+        if nargs == 2:
+            fname, hname = args.hist_name.split(":")
+            sg_dict["energy"] = PG.TH1Sampler(fname, hname)
+        elif nargs == 3:
+            fname, hname, scale = args.hist_name.split(":")
+            sg_dict["energy"] = PG.TH1Sampler(fname, hname, scale)
+        else:
+            print(f"Can't parse histogram {args.hist_name}!")
+            sys.exit(1)
+
     elif args.sampler == "hist2D":
-        fname, hname = args.hist_name.split(":")
-        sg_dict["energy"] = PG.TH2Sampler(fname, hname)        
+        nargs = len(args.hist_name.split(":"))
+        if nargs == 2:
+            fname, hname = args.hist_name.split(":")
+            sg_dict["energy"] = PG.TH2Sampler(fname, hname)
+        elif nargs == 4:
+            fname, hname, scalex, scaley = args.hist_name.split(":")
+            sg_dict["energy"] = PG.TH2Sampler(fname, hname, scalex, scaley)
+        else:
+            print(f"Can't parse histogram {args.hist_name}!")
+            sys.exit(1)
     else:
         print(f"Sampler {args.sampler} not known!")
         sys.exit(1)
@@ -255,8 +278,12 @@ if __name__ == '__main__':
 
     b = time.time()
     log.info("Run G4FaserAlg in " + str(b-a) + " seconds")
-#
-# Success should be 0
-#
-    sys.exit(not sc.isSuccess())
+
+# Signal errors
+if sc.isSuccess():
+    log.info("Execution succeeded")
+    sys.exit(0)
+else:
+    log.info("Execution failed, return 1")
+    sys.exit(1)
 

Control/CalypsoExample/Generation/scripts/submit_faser_particlegun.sh

@@ -185,9 +185,12 @@ cd "${config_file_stem}-${seg_str}"
 # Run job
 if [[ -z "$tag" ]]; then
     faser_particlegun.py  "--conf=$config_path" "--segment=$seg_str"
+    gen_code=$?
 else
     faser_particlegun.py  "--conf=$config_path" "--segment=$seg_str" "--tag=$tag"
+    gen_code=$?
 fi
+echo "Return code: $gen_code"
 #
 # Print out ending time
 date
@@ -198,25 +201,51 @@ export EOS_MGM_URL=root://eospublic.cern.ch
 #
 if ! [ -z "$outdest" ]
 then
+    echo "Output directory:"
     ls -l
-    echo "copy *-HITS.root to $outdest"
-    mkdir -p $outdest
-    eos cp *-HITS.root ${outdest}/ || true
+    thefile=`ls *-HITS.root`
+    if [ $? -eq 0 ]; then
+	echo "copy $thefile to $outdest"
+	eos mkdir -p $outdest
+	eos cp $thefile ${outdest}/${thefile} || true
+
+	# Check that it worked
+	eos ls ${outdest}/${thefile} > /dev/null
+	if [ $? -eq 0 ]; then
+	    echo "file $thefile copied to $outdest"
+	    copy_code=0
+	else
+	    echo "didnt find $thefile in $outdest !"
+	    copy_code=1
+	fi
+    else
+	echo "ls *-xAOD.root returned nothing!"
+	copy_code=1
+    fi
+
 fi
 #
 # Also copy log file
 if ! [ -z "$logdest" ]
 then
     cd ..
+    echo "Working directory:"
     ls -l
     echo "copy $logfile to $logdest"
-    mkdir -p $logdest
+    eos mkdir -p $logdest
     eos cp $logfile $logdest/$logfile
 elif ! [ -z "$outdest" ]
 then 
     cd ..
     ls -l
     echo "copy $logfile to $outdest"
-    mkdir -p $outdest
+    eos mkdir -p $outdest
     eos cp $logfile $outdest/$logfile
 fi
+
+# Make sure to return an error is calypso failed
+if [ $gen_code -ne 0 ] || [ $copy_code -ne 0 ]; then
+  exit 1
+else
+  exit 0
+fi

Generators/ParticleGun/python/histsampling.py

@@ -14,6 +14,7 @@ def load_hist(*args):
     """
     Load a histogram from a filename/TFile and histo name. If a single arg is
     provided, it has to be a histo object and will be cloned before return.
+
     """
     h = None
     if len(args) == 1 and issubclass(type(args[0]), ROOT.TH1):
@@ -25,6 +26,7 @@ def load_hist(*args):
             #f.Close()
         elif type(args[0]) is ROOT.TFile and type(args[1]) is str:
             h = args[0].Get(args[1]).Clone()
+
     if h is None:
         raise Exception("Error in histogram loading from " + args)
     else:
@@ -77,7 +79,7 @@ def get_random_bin(globalbins, cheights):
     raise Exception("Sample fell outside range of cumulative distribution?!?!")
 
 
-def get_random_x(h, globalbins, cheights, globalbin_to_axisbin):
+def get_random_x(h, globalbins, cheights, globalbin_to_axisbin, scale=1.):
     """
     Choose a random bin via get_random_bin, then pick a uniform random x
     point in that bin (without any attempt at estimating the in-bin distribution).
@@ -86,10 +88,10 @@ def get_random_x(h, globalbins, cheights, globalbin_to_axisbin):
     axisids = globalbin_to_axisbin.get(irand)
     assert axisids is not None
     xrand = random.uniform(h.GetXaxis().GetBinLowEdge(axisids[0]), h.GetXaxis().GetBinUpEdge(axisids[0]))
-    return xrand
+    return scale * xrand
 
 
-def get_random_xy(h2, globalbins, cheights, globalbin_to_axisbin):
+def get_random_xy(h2, globalbins, cheights, globalbin_to_axisbin, xscale=1., yscale=1.):
     """
     Choose a random bin via get_random_bin, then pick a uniform random x,y
     point in that bin (without any attempt at estimating the in-bin distribution).
@@ -99,34 +101,71 @@ def get_random_xy(h2, globalbins, cheights, globalbin_to_axisbin):
     assert axisids is not None
     xrand = random.uniform(h2.GetXaxis().GetBinLowEdge(axisids[0]), h2.GetXaxis().GetBinUpEdge(axisids[0]))
     yrand = random.uniform(h2.GetYaxis().GetBinLowEdge(axisids[1]), h2.GetYaxis().GetBinUpEdge(axisids[1]))
-    return xrand, yrand
+    return (xscale*xrand), (yscale*yrand)
 
 
 class TH1(object):
     "Minimal wrapper for ROOT TH1, for sampling consistency and easy loading"
 
     def __init__(self, *args):
-        self.th1 = load_hist(*args)
+        """ Args can be variable length, but these are now kwargs, so order matters
+
+        fname - filename
+        hname - histogram name
+        xscale - scaling factor for x axis variable
+        """
+        if len(args) > 2:
+            self.th1 = load_hist(*(args[0:2]))
+        else:
+            self.th1 = load_hist(*args)
+
         self.globalbins, self.globalbin_to_axisbin, self.cheights = None, None, None
 
+        if len(args) >= 3:
+            self.xscale = float(args[2])
+        else:
+            self.xscale = 1.
+
     def GetRandom(self):
         "A GetRandom that works for TH1s and uses Python random numbers"
         if self.globalbins is None or self.globalbin_to_axisbin is None or self.cheights is None:
             self.globalbins, self.globalbin_to_axisbin, self.cheights = get_sampling_vars(self.th1)
-        return get_random_x(self.th1, self.globalbins, self.cheights, self.globalbin_to_axisbin)
+        return get_random_x(self.th1, self.globalbins, self.cheights, self.globalbin_to_axisbin, self.xscale)
 
 
 class TH2(object):
     "Minimal wrapper for ROOT TH2, for easy loading and to allow 2D sampling"
 
     def __init__(self, *args):
-        self.th2 = load_hist(*args)
+        """ Args can be variable length, but these are now kwargs, so order matters
+
+        fname - filename
+        hname - histogram name
+        xscale - scaling factor for x axis variable
+        yscale - scaling factor for y axis variable
+        """
+        if len(args) > 2:
+            self.th2 = load_hist(*(args[0:2]))
+        else:
+            self.th2 = load_hist(*args)
+
         self.globalbins, self.globalbin_to_axisbin, self.cheights = None, None, None
 
+        if len(args) >= 3:
+            self.xscale = float(args[2])
+        else:
+            self.xscale = 1.
+
+        if len(args) >= 4:
+            self.yscale = float(args[3])
+        else:
+            self.yscale = 1.
+
+
     def GetRandom(self):
         "A GetRandom that works for TH2s"
         if self.globalbins is None or self.globalbin_to_axisbin is None or self.cheights is None:
             self.globalbins, self.globalbin_to_axisbin, self.cheights = get_sampling_vars(self.th2)
-        return get_random_xy(self.th2, self.globalbins, self.cheights, self.globalbin_to_axisbin)
+        return get_random_xy(self.th2, self.globalbins, self.cheights, self.globalbin_to_axisbin, self.xscale, self.yscale)
 
 

Generators/ParticleGun/python/samplers.py

@@ -633,14 +633,23 @@ class EThetaMPhiSampler(MomSampler):
         pt = p sin(theta)
         """
 
-        if self._theta is None:
-            e,theta = self.energy()
-        else:            
-            e = self.energy()
-            theta = self.theta()
-            
         m = self.mass()
-        p = math.sqrt( e**2 - m**2 )
+
+        count = 0
+        e = -1
+        while (e < m and count < 5):
+            count += 1
+            if self._theta is None:
+                e,theta = self.energy()
+            else:            
+                e = self.energy()
+                theta = self.theta()
+
+        try:            
+            p = math.sqrt( e**2 - m**2 )
+        except Exception:
+            raise Exception(f"Error generating E: {e} m: {m}!")
+
         pz = p * math.cos(theta)
         pt = p * math.sin(theta)
         phi = self.phi()