From 26b84127b59a3c7d1ff3a3bc4f3700e8bb29a3d9 Mon Sep 17 00:00:00 2001
From: Carl Gwilliam <gwilliam@hep.ph.liv.ac.uk>
Date: Wed, 24 May 2023 10:18:28 +0100
Subject: [PATCH] Update generator scripts for B-L and add grids for ALPS
 models

---
 .../ForeseeGenerator/python/Validate.py       |   4 +-
 .../share/generate_forsee_events.py           | 100 +++++++++++++++++-
 .../share/point_ALP_photon.json               |  14 +++
 Generators/ForeseeGenerator/share/points.json |  18 ++++
 .../ForeseeGenerator/share/points_ALP_W.json  |  15 +++
 .../share/points_ALP_photon.json              |  14 +++
 .../share/points_BminusL.json                 |  25 +++++
 .../ForeseeGenerator/share/points_extra.json  |  19 ++++
 .../ForeseeGenerator/share/points_extra2.json |  12 +++
 .../ForeseeGenerator/share/points_john.json   |  12 +++
 .../ForeseeGenerator/share/validate_grid.py   |   9 +-
 11 files changed, 232 insertions(+), 10 deletions(-)
 create mode 100644 Generators/ForeseeGenerator/share/point_ALP_photon.json
 create mode 100644 Generators/ForeseeGenerator/share/points.json
 create mode 100644 Generators/ForeseeGenerator/share/points_ALP_W.json
 create mode 100644 Generators/ForeseeGenerator/share/points_ALP_photon.json
 create mode 100644 Generators/ForeseeGenerator/share/points_BminusL.json
 create mode 100644 Generators/ForeseeGenerator/share/points_extra.json
 create mode 100644 Generators/ForeseeGenerator/share/points_extra2.json
 create mode 100644 Generators/ForeseeGenerator/share/points_john.json

diff --git a/Generators/ForeseeGenerator/python/Validate.py b/Generators/ForeseeGenerator/python/Validate.py
index 58f34f08b..b0bd5c262 100644
--- a/Generators/ForeseeGenerator/python/Validate.py
+++ b/Generators/ForeseeGenerator/python/Validate.py
@@ -109,13 +109,13 @@ class EvgenValidation(EvgenAnalysisAlg):
         # Vertex
         self.hists.add("Vtx_X_LLP", 50, -100, 100)
         self.hists.add("Vtx_Y_LLP", 50, -100, 100)
-        self.hists.add("Vtx_Z_LLP", 500, -5000, 0)
+        self.hists.add("Vtx_Z_LLP", 500, -5000, 5000)
         self.hists.add("Vtx_R_LLP", 20, 0, 200)        
         self.hists.add("Vtx_XY_LLP", 50, -100, 100, 50, -100, 100)
 
         self.hists.add("Vtx_X_All", 50, -100, 100)
         self.hists.add("Vtx_Y_All", 50, -100, 100)
-        self.hists.add("Vtx_Z_All", 500, -5000, 0)
+        self.hists.add("Vtx_Z_All", 500, -5000, 5000)
         self.hists.add("Vtx_R_All", 20, 0, 200)        
         self.hists.add("Vtx_XY_All", 50, -100, 100, 50, -100, 100)      
         
diff --git a/Generators/ForeseeGenerator/share/generate_forsee_events.py b/Generators/ForeseeGenerator/share/generate_forsee_events.py
index 6ef5d1d03..9190f70fa 100644
--- a/Generators/ForeseeGenerator/share/generate_forsee_events.py
+++ b/Generators/ForeseeGenerator/share/generate_forsee_events.py
@@ -1,4 +1,4 @@
-import os
+import os, sys
 
 import numpy as np
 import matplotlib.pyplot as plt
@@ -53,8 +53,17 @@ class ForeseeGenerator(object):
 #                                   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.065)**2)< 0.1",  
-                                  channels=[self.mode], distance=480, length=1.5 ,
+        # Generate 6.5 cm high to account for translation from ATLAS to FASER coord. system
+        # Relax r to 11 cm to acount for any mismatches with calypso and for duaghter entering if mother doesn't
+
+        if "ALP" in self.modelname:
+            # Since not relying on tracks then extend length to include spectrometer (magents + tracking stations)
+            length = 1.5 + 2 + 0.5 # m
+        else:
+            length = 1.5 # m
+
+        self.foresee.set_detector(selection="np.sqrt(x.x**2 + (x.y - 0.065)**2)< 0.11",  
+                                  channels=[self.mode], distance=480, length=length ,
                                   luminosity=1/1000.) # 1 pb-1        
 
 
@@ -68,6 +77,10 @@ class ForeseeGenerator(object):
             self.data = self.darkphoton()
         elif self.modelname == "ALP-W":
             self.data = self.alp_W()
+        elif self.modelname == "ALP-photon":
+            self.data = self.alp_photon()
+        elif self.modelname == "U(1)B-L":
+            self.data = self.bminusl()                        
         else:
             sys.exit(f"Unknown model {self.modelname}")
 
@@ -227,6 +240,84 @@ class ForeseeGenerator(object):
                 )   
 
         return self.decays()
+        
+
+    def alp_photon(self):
+
+        if self.version == 1:
+            sys.exit("ALP-photon model not available in this version of FORESEE")
+
+        self.nbinsample = 100 # resample bins to help with smoothness
+        
+        masses = [float(x) for x in ['{:0.4e}'.format(m) for m in np.logspace(-2,0,20+1)]]
+        self.model.add_production_direct(
+            label = "Prim",
+            energy = self.energy,
+            coupling_ref = 1,
+            masses=masses,
+            )
+
+        self.model.set_br_1d(
+            modes = ["gamma_gamma"],
+            finalstates = [[22, 22]],
+            filenames=["model/br/gamma_gamma.txt"]
+            )   
+        
+        
+        return self.decays()
+
+    def bminusl(self):
+
+        if self.version == 1:
+            sys.exit("U(1)B-L model not available in this version of FORESEE")
+
+        self.nbinsample = 100 # TODO
+
+        self.model.add_production_2bodydecay(
+            pid0 =  "111",
+            pid1 = "22",
+            br = "2.*0.99 * (coupling/0.303)**2 * pow(1.-pow(mass/self.masses('111'),2),3)",
+            generator = "EPOSLHC",
+            energy = self.energy,
+            nsample = 100)
+    
+        self.model.add_production_2bodydecay(
+            pid0 = "221",
+            pid1 = "22",
+            br = "2.*0.25*0.39 * (coupling/0.303)**2 * pow(1.-pow(mass/self.masses('221'),2),3)",
+            generator = "EPOSLHC",
+            energy = self.energy,
+            nsample = 100)
+
+
+        masses_brem = [ 
+            0.01  ,  0.0126,  0.0158,  0.02  ,  0.0251,  0.0316,  0.0398,
+            0.0501,  0.0631,  0.0794,  0.1   ,  0.1122,  0.1259,  0.1413,
+            0.1585,  0.1778,  0.1995,  0.2239,  0.2512,  0.2818,  0.3162,
+            0.3548,  0.3981,  0.4467,  0.5012,  0.5623,  0.6026,  0.631 ,
+            0.6457,  0.6607,  0.6761,  0.6918,  0.7079,  0.7244,  0.7413,
+            0.7586,  0.7762,  0.7943,  0.8128,  0.8318,  0.8511,  0.871 ,
+            0.8913,  0.912 ,  0.9333,  0.955 ,  0.9772,  1.    ,  1.122 ,
+            1.2589,  1.4125,  1.5849,  1.7783,  1.9953,  2.2387,  2.5119,
+            2.8184,  3.1623,  3.9811,  5.0119,  6.3096,  7.9433, 10.    
+            ]
+        
+        self.model.add_production_direct(
+            label = "Brem",
+            energy = self.energy,
+            condition = "1./0.3**2 * (p.pt<1)",
+            coupling_ref=1,
+            masses = masses_brem,
+            )
+            
+        self.model.set_br_1d(
+        modes = ["e_e", "mu_mu", "nu_nu", "Hadrons"],
+            finalstates = [[11, -11], [13, -13], [12, -12], None],
+            filenames=["model/br/e_e.txt", "model/br/mu_mu.txt", "model/br/nu_nu.txt", "model/br/Hadrons.txt"],
+            #filenames=[f"model/br/all.txt"] 
+            )
+        
+        return self.decays()
 
 
     def decays(self):
@@ -468,7 +559,8 @@ def production(args):
                 "Mass (GeV)" : m,
                 "Coupling" : f"{couplings[0]:e}",
                 "Cross-section (pb)" : sum(weights[0]), # / data["nevents"],
-                "ECom (TeV)" : float(args.Ecom)
+                "ECom (TeV)" : float(args.Ecom),
+                "NEvents" : float(data["nevents"])                
                 }
 
             print(xsect)
diff --git a/Generators/ForeseeGenerator/share/point_ALP_photon.json b/Generators/ForeseeGenerator/share/point_ALP_photon.json
new file mode 100644
index 000000000..986bbc376
--- /dev/null
+++ b/Generators/ForeseeGenerator/share/point_ALP_photon.json
@@ -0,0 +1,14 @@
+{
+"name" : "ALPPhotontGrid13p6_65mm",
+"model" : "ALP-photon",
+"pid1" : 22,
+"pid2" : 22,
+"nevents" : 50000,
+"nrun" : 110100,
+"samples" : {
+	"7e-05" : [0.06],
+	"1.5e-04" : [0.06],
+	"3e-04" : [0.06]
+   }
+}
+
diff --git a/Generators/ForeseeGenerator/share/points.json b/Generators/ForeseeGenerator/share/points.json
new file mode 100644
index 000000000..ac3df72f7
--- /dev/null
+++ b/Generators/ForeseeGenerator/share/points.json
@@ -0,0 +1,18 @@
+{
+"name" : "DarkPhotonGrid13p6_65mm",
+"model" : "DarkPhoton",
+"pid1" : 11,
+"pid2" : -11,
+"nevents" : 50000,
+"nrun" : 110021,
+"samples" : {
+       "1e-06" : [0.01, 0.0251, 0.0501, 0.0631, 0.0794],
+       "5e-06" : [0.01, 0.0251, 0.0501, 0.0631, 0.0794, 0.1, 0.1122],
+       "1e-05" : [0.01, 0.0251, 0.0501, 0.0631, 0.0794, 0.1, 0.1122],
+       "3e-05" : [0.01, 0.0251, 0.0501, 0.0631, 0.0794, 0.1, 0.1122],
+       "5e-05" : [0.01, 0.0251, 0.0501, 0.0631],
+       "1e-04" : [0.01, 0.0251, 0.0501],
+       "3e-04" : [0.01, 0.0251]
+   }
+}
+
diff --git a/Generators/ForeseeGenerator/share/points_ALP_W.json b/Generators/ForeseeGenerator/share/points_ALP_W.json
new file mode 100644
index 000000000..8afa2bfa2
--- /dev/null
+++ b/Generators/ForeseeGenerator/share/points_ALP_W.json
@@ -0,0 +1,15 @@
+{
+"name" : "ALPWGrid13p6_65mm_110mm_3p5m_extra",
+"model" : "ALP-W",
+"pid1" : 22,
+"pid2" : 22,
+"nevents" : 50000,
+"nrun" : 110206,
+"samples" : {
+	"4e-04" : [0.1],
+	"1e-04" : [0.2],
+	"5e-05" : [0.1],
+	"2e-05" : [0.35]
+   }
+}
+
diff --git a/Generators/ForeseeGenerator/share/points_ALP_photon.json b/Generators/ForeseeGenerator/share/points_ALP_photon.json
new file mode 100644
index 000000000..dd31b37b9
--- /dev/null
+++ b/Generators/ForeseeGenerator/share/points_ALP_photon.json
@@ -0,0 +1,14 @@
+{
+"name" : "ALPPhotonGrid13p6_65mm_110mm_3p5m",
+"model" : "ALP-photon",
+"pid1" : 22,
+"pid2" : 22,
+"nevents" : 50000,
+"nrun" : 110103,
+"samples" : {
+	"7e-05" : [0.06],
+	"1.5e-04" : [0.06],
+	"3e-04" : [0.06]
+   }
+}
+
diff --git a/Generators/ForeseeGenerator/share/points_BminusL.json b/Generators/ForeseeGenerator/share/points_BminusL.json
new file mode 100644
index 000000000..c0574ccbc
--- /dev/null
+++ b/Generators/ForeseeGenerator/share/points_BminusL.json
@@ -0,0 +1,25 @@
+{
+"name" : "U1B-LGrid13p6_65mm_110mm",
+"model" : "U(1)B-L",
+"pid1" : 11,
+"pid2" : -11,
+"nevents" : 50000,
+"nrun" : 110300,
+"samples" : {
+	"4e-5"   : [0.01],
+	"3e-5"   : [0.01],
+	"2.4e-5" :       [0.014],
+	"2e-5"   :       [0.014],
+	"1.4e-5" :              [0.02],
+	"1.1e-5" :              [0.02, 0.024],
+	"1e-5"   :                    [0.024],
+	"8e-6"   :                           [0.03],
+	"7e-6"   :                           [0.03],
+	"6e-6"   :                           [0.03],
+	"5e-6"   :                           [0.03, 0.034],
+	"4e-6"   : [0.01, 0.014, 0.02,        0.03, 0.034, 0.04],
+	"3e-6  " : [0.01, 0.014, 0.02, 0.024, 0.03, 0.034, 0.04],
+	"2.4e-6" : [                   0.024, 0.03]
+   }
+}
+
diff --git a/Generators/ForeseeGenerator/share/points_extra.json b/Generators/ForeseeGenerator/share/points_extra.json
new file mode 100644
index 000000000..2557f41bf
--- /dev/null
+++ b/Generators/ForeseeGenerator/share/points_extra.json
@@ -0,0 +1,19 @@
+{
+"name" : "DarkPhotonGrid13p6_65mm_extra",
+"model" : "DarkPhoton",
+"pid1" : 11,
+"pid2" : -11,
+"nevents" : 50000,
+"nrun" : 110056,
+"samples" : {
+       "3e-04" : [0.0126, 0.0158],
+       "1e-04" : [0.0126, 0.0158],
+       "5e-05" : [        0.0158, 0.0316, 0.0358],
+       "3e-05" : [                0.0316, 0.0358],
+       "2e-05" : [                                              0.0501, 0.0631, 0.0794],
+       "1e-05" : [0.0126, 0.0158, 0.0316, 0.0358],
+       "5e-06" : [0.0126, 0.0158, 0.0316, 0.0358]
+       "3e-06" : [0.01, 0.0126, 0.0158, 0.0251, 0.0316, 0.0358, 0.0501, 0.0631, 0.0794],
+   }
+}
+
diff --git a/Generators/ForeseeGenerator/share/points_extra2.json b/Generators/ForeseeGenerator/share/points_extra2.json
new file mode 100644
index 000000000..2649af218
--- /dev/null
+++ b/Generators/ForeseeGenerator/share/points_extra2.json
@@ -0,0 +1,12 @@
+{
+"name" : "DarkPhotonGrid13p6_65mm_extra2",
+"model" : "DarkPhoton",
+"pid1" : 11,
+"pid2" : -11,
+"nevents" : 50000,
+"nrun" : 110090,
+"samples" : {
+	"1.5e-04" : [0.01, 0.0126, 0.0158, 0.0251]
+   }
+}
+
diff --git a/Generators/ForeseeGenerator/share/points_john.json b/Generators/ForeseeGenerator/share/points_john.json
new file mode 100644
index 000000000..1aee546ca
--- /dev/null
+++ b/Generators/ForeseeGenerator/share/points_john.json
@@ -0,0 +1,12 @@
+{
+"name" : "DarkPhotonGrid13p6b",
+"model" : "DarkPhoton",
+"pid1" : 11,
+"pid2" : -11,
+"nevents" : 50000,
+"nrun" : 110056,
+"samples" : {
+	"3e-06" : [0.0398]
+   }
+}
+
diff --git a/Generators/ForeseeGenerator/share/validate_grid.py b/Generators/ForeseeGenerator/share/validate_grid.py
index 88edb1194..7e13d5779 100644
--- a/Generators/ForeseeGenerator/share/validate_grid.py
+++ b/Generators/ForeseeGenerator/share/validate_grid.py
@@ -3,19 +3,20 @@ import subprocess as proc
 from glob import glob
 import os, sys
 
-grid = "../calypso/Generators/ForeseeGenerator/share/points.json"
+grid = "../calypso/Generators/ForeseeGenerator/share/points_ALP_W.json"
 
 with open(grid) as f:
     data = json.load(f)
 
-name = "DarkPhotonGrid13p6_65mm"
+name = data["name"]
 path = f"/bundle/data/FASER/Carl_Output/Foresee/{name}/"
 ecom = 13.6
-pid1 = 11
-pid2 = -11
+pid1 = data["pid1"]
+pid2 = data["pid2"]
 
 for coup, masses in data["samples"].items():
     for m in masses:
+        print(coup, m, 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:
-- 
GitLab