refining angle defs in 2-3 scattering; adding mtrx elements for 2nu2gamma interaction; minor track length fixes in FluxBremIsotropic

athompson-git · athompson-git · commit 8d5fa785271c · 2023-10-08T14:27:57.000-05:00
diff --git a/cross_section_mc.py b/cross_section_mc.py
@@ -252,8 +252,10 @@ class Scatter2to3MC:
     """
     2 -> 3 scattering for general masses in fixed target frame
     Based on Byckling, Kajantie ch. V.4
+    
+    ATTN: assumes particle A is at rest.
     """
-    def __init__(self, mtrx2: MatrixElement2to3, p1=LorentzVector(), p2=LorentzVector(), n_samples=1000):
+    def __init__(self, mtrx2: MatrixElement2to3, p1=LorentzVector(), p2=LorentzVector(), n_samples=100):
         self.mtrx2 = mtrx2
 
         self.ma = mtrx2.ma
@@ -274,6 +276,7 @@ def __init__(self, mtrx2: MatrixElement2to3, p1=LorentzVector(), p2=LorentzVecto
 
     def cayley_det(self, x, y, z, u, v, w):
         # Gramm-Schmidt det but with invariants
+        # Eq. 5.24 of Byckling, Kajantie
         return -0.5 * np.linalg.det([[0, 1, 1, 1, 1],
                                     [1, 0, v, x, z],
                                     [1, v, 0, u, y],
@@ -283,6 +286,9 @@ def cayley_det(self, x, y, z, u, v, w):
     def kallen(self, x, y, z):
         return np.power((x - y - z), 2) - 4*y*z
     
+    def flux_factor(self, s):
+        return power(2*sqrt(self.kallen(s, self.ma**2, self.mb**2))*(2*pi)**5, -1)
+    
     def s2MaxMin(self, s):
         return (sqrt(s) - self.m1)**2, (self.m2 + self.m3)**2
     
@@ -293,7 +299,7 @@ def t1MaxMin(self, s, s2):
             + np.nan_to_num(sqrt(self.kallen(s, self.ma**2, self.mb**2)*self.kallen(s, s2, self.m1**2))))/(2*s)
     
     def phase_space_heaviside(self, s, t1, s2):
-        return np.heaviside(-np.nan_to_num(self.cayley_det(s, t1, s2, self.ma**2, self.mb**2, self.m1**2)), 0.0) \
+        return np.heaviside(np.nan_to_num(-self.cayley_det(s, t1, s2, self.ma**2, self.mb**2, self.m1**2)), 0.0) \
             * np.heaviside(self.kallen(s, self.m1**2, s2), 0.0) \
             * np.heaviside(self.kallen(s2, self.m2**2, self.m3**2), 0.0)
     
@@ -310,11 +316,30 @@ def s1_from_angle(self, phib, s, t1, s2, t2):
                              [s-self.ma**2 + self.mb**2, s+s2-self.m1**2, 0.0]]) \
                 + 2*np.nan_to_num(sqrt(self.cayley_det(s,t1,s2,self.ma**2,self.mb**2,self.m1**2)*self.cayley_det(s2,t2,self.m3**2,t1,self.mb**2,self.m2**2)))*cos(phib))
 
+    def paR23(self, s, t1, s2):
+        return sqrt(((s + t1 - self.mb**2 - self.m1**2)/(2*sqrt(s2)))**2 - self.ma**2)
+    
+    def pbR23(self, s2, t1):
+        return sqrt(self.kallen(s2, self.mb**2, t1)/s2) / 2
+    
+    def p1R23(self, s, s2):
+        e1_cm = (s + self.m1**2 - s2)/(2*sqrt(s))
+        return np.sqrt(e1_cm**2 - self.m1**2)
+
+    def sinTheta_b1(self, s, t1, s2):
+        # R23 frame
+        return -4 * s2 * self.cayley_det(s, t1, s2, self.ma**2, self.mb**2, self.m1**2) \
+            / (self.kallen(s2, self.mb**2, t1)*self.kallen(s, s2, self.m1**2))
+    
+    def cosTheta_ab(self, s, s2, t1):
+        sinTheta_b1 = self.sinTheta_b1(s, t1, s2)
+        cosTheta_b1 = sqrt(1 - sinTheta_b1**2)
+        return (self.p1R23(s, s2)*self.pbR23(s2, t1)*cosTheta_b1 - self.pbR23(s2, t1)**2) / (self.pbR23(s2, t1)*self.paR23(s, t1, s2))
+
     def dsigma_ds2dt1dOmega3(self, s, s2, t1, cosThetab3, phib):
         t2 = self.t2_from_angle(cosThetab3, t1, s2)
         s1 = self.s1_from_angle(phib, s, t1, s2, t2)
-        return pi * self.phase_space_heaviside(s, t1, s2) * self.mtrx2(s, s2, t1, s1, t2) \
-            * power(2*sqrt((s - (self.ma + self.mb)**2)*(s - (self.ma - self.mb)**2)), -1) \
+        return pi * self.phase_space_heaviside(s, t1, s2) * self.mtrx2(s, s2, t1, s1, t2) * self.flux_factor(s) \
             * np.nan_to_num(sqrt(self.kallen(s, self.ma**2, self.mb**2)/self.kallen(s2, self.m2**2, self.m3**2)))/(16*s2)
 
     def r23_velocity(self, s2, t1):
@@ -364,23 +389,76 @@ def simulate_particle1(self, s):
         t1_rnd = np.random.uniform(t1Min, t1Max, self.n_samples)
         cos_rnd = np.random.uniform(-1, 1, self.n_samples)
         phi_rnd = np.random.uniform(0.0, 2*pi, self.n_samples)
+
+        print("t1 min, max = ", t1Min, t1Max)
+        print("s2 min, max = ", s2Min, s2Max)
         
         mc_volume = (t1Max - t1Min)*(s2Max - s2Min)*4*pi / self.n_samples
         
         # Get particle 1's CM spectra
         for i in range(self.n_samples):
+            if self.phase_space_heaviside(s, t1_rnd[i], s2_rnd[i]) < 1.0:
+                continue
+
             weights = self.dsigma_ds2dt1dOmega3(s, s2_rnd[i], t1_rnd[i], cos_rnd[i], phi_rnd[i])
             e1_cm = (s + self.m1**2 - s2_rnd[i])/(2*sqrt(s))
             p1_cm = np.sqrt(e1_cm**2 - self.m1**2)
             theta_b1_cm = arcsin(sqrt(-4*s2_rnd[i]*self.cayley_det(s, t1_rnd[i], s2_rnd[i], self.ma**2, self.mb**2, self.m1**2) \
                 / (self.kallen(s2_rnd[i], self.mb**2, t1_rnd[i])*self.kallen(s, s2_rnd[i], self.m1**2))))
             lorentz_vector_p1_cm = LorentzVector(e1_cm, p1_cm*sin(theta_b1_cm), 0.0, p1_cm*cos(theta_b1_cm))
-            v_cm_to_lab = self.r23_velocity(s2_rnd[i], t1_rnd[i])
-            lorentz_vector_p1_lab = lorentz_boost(lorentz_vector_p1_cm, v_cm_to_lab)
+
+            # Boost to lab frame: use pa velocity assuming pa is at rest
+            pa_R23 = self.paR23(s, t1_rnd[i], s2_rnd[i])
+            ea_R23 = sqrt(pa_R23**2 + self.ma**2)
+            cosTheta_ab_R23 = self.cosTheta_ab(s, s2_rnd[i], t1_rnd[i])
+            va_R23_3vec = Vector3(-pa_R23*sqrt(1-cosTheta_ab_R23**2)/ea_R23, 0.0, -pa_R23*cosTheta_ab_R23/ea_R23)
+
+            #v_cm_to_lab = self.r23_velocity(s2_rnd[i], t1_rnd[i])
+            lorentz_vector_p1_lab = lorentz_boost(lorentz_vector_p1_cm, va_R23_3vec)
             self.p1_lab_4vectors.append(lorentz_vector_p1_lab)
             self.p1_lab_energies.append(lorentz_vector_p1_lab.p0)
+            #print("dsigma = {}, CM energy = {}, Lab energy = {}".format(weights, lorentz_vector_p1_cm.p0, lorentz_vector_p1_lab.p0))
             self.p1_lab_angles.append(lorentz_vector_p1_lab.theta())
             self.dsigma.append(weights * mc_volume)
+    
+    def simulate_particle1_phase_space_sampled(self, s):
+        s2Max, s2Min = self.s2MaxMin(s)
+
+        s2_grid = np.linspace(s2Min, s2Max, self.n_samples)
+        s2_points = (s2_grid[1:] + s2_grid[:-1])/2
+        delta_s2 = s2_grid[1] - s2_grid[0]
+
+        cos_rnd = np.random.uniform(-1, 1, self.n_samples)
+        phi_rnd = np.random.uniform(0.0, 2*pi, self.n_samples)
+        
+        # Get particle 1's CM spectra
+        for s2 in s2_points:
+            t1Max, t1Min = self.t1MaxMin(s, s2)
+            t1_rnd = np.random.uniform(t1Min, t1Max, self.n_samples)
+            mc_volume = (t1Max - t1Min)*delta_s2*4*pi / self.n_samples
+            for i, t1 in enumerate(t1_rnd):
+                if self.phase_space_heaviside(s, t1, s2) < 1.0:
+                    continue
+
+                weights = self.dsigma_ds2dt1dOmega3(s, s2, t1, cos_rnd[i], phi_rnd[i])
+                e1_cm = (s + self.m1**2 - s2)/(2*sqrt(s))
+                p1_cm = np.sqrt(e1_cm**2 - self.m1**2)
+                theta_b1_cm = arcsin(sqrt(-4*s2*self.cayley_det(s, t1, s2, self.ma**2, self.mb**2, self.m1**2) \
+                    / (self.kallen(s2, self.mb**2, t1)*self.kallen(s, s2, self.m1**2))))
+                lorentz_vector_p1_cm = LorentzVector(e1_cm, p1_cm*sin(theta_b1_cm), 0.0, p1_cm*cos(theta_b1_cm))
+
+                # Boost to lab frame: use pa velocity assuming pa is at rest
+                pa_R23 = self.paR23(s, t1, s2)
+                ea_R23 = sqrt(pa_R23**2 + self.ma**2)
+                cosTheta_ab_R23 = self.cosTheta_ab(s, s2, t1)
+                va_R23_3vec = Vector3(-pa_R23*sqrt(1-cosTheta_ab_R23**2)/ea_R23, 0.0, -pa_R23*cosTheta_ab_R23/ea_R23)
+
+                lorentz_vector_p1_lab = lorentz_boost(lorentz_vector_p1_cm, va_R23_3vec)
+                self.p1_lab_4vectors.append(lorentz_vector_p1_lab)
+                self.p1_lab_energies.append(lorentz_vector_p1_lab.p0)
+                self.p1_lab_angles.append(lorentz_vector_p1_lab.theta())
+                self.dsigma.append(weights * mc_volume)
+
 
 
 
diff --git a/fluxes.py b/fluxes.py
@@ -293,22 +293,23 @@ def simulate_single(self, electron):
         self.axion_energy.extend(ea_rnd)
         self.axion_flux.extend(el_wgt * diff_br)
 
-    def simulate(self):
+    def simulate(self, use_track_length=False):
         self.axion_energy = []
         self.axion_flux = []
         self.scatter_axion_weight = []
         self.decay_axion_weight = []
 
-        """
-        ep_min = max((self.ma**2 - M_E**2)/(2*M_E), M_E)
-        for i, el in enumerate(self.electron_flux):
-            t_depth = np.random.uniform(0.0, 5.0)
-            new_energy = np.random.uniform(ep_min, el[0])
-            flux_weight = self.electron_flux_attenuated(t_depth, el[0], new_energy) * 5.0 * (el[0] - ep_min)
-            self.simulate_single([new_energy, flux_weight])
-        """
-        for i, el in enumerate(self.electron_flux):
-            self.simulate_single(el)
+        if use_track_length:
+            ep_min = max((self.ma**2 - M_E**2)/(2*M_E), M_E)
+            for i, el in enumerate(self.electron_flux):
+                t_depth = np.random.uniform(0.01, 5.0)
+                new_energy = np.random.uniform(ep_min, el[0])
+                flux_weight = self.electron_flux_attenuated(t_depth, el[0], new_energy) * 5.0 * (el[0] - ep_min)
+                self.simulate_single([new_energy, el[1] * flux_weight])
+            
+        else:
+            for i, el in enumerate(self.electron_flux):
+                self.simulate_single(el)
 
     def propagate(self, new_coupling=None):
         if new_coupling is not None:
diff --git a/matrix_element.py b/matrix_element.py
@@ -746,4 +746,44 @@ def __call__(self, s, s2, t1, s1, t2, gae=1.0):
         m = self.m2
         M = self.ma
         prefactor = self.ffp(t1) * self.ff2(sqrt(abs(t1))) * power(gae * 4*pi*ALPHA, 2.0) / power(t2 * (m**2-t2) * (m**2+mA**2-s2-t2+t1), 2)
-        return prefactor * 4*(m**4*(2*M**4*(mA**2+t1)+2*M**2*(3*mA**4+mA**2*(3*t1-2*(s+s2))-t1*(2*s+s2)+s2**2)+2*(mA**3-mA*s)**2+t1*(mA**4+2*mA**2*(s2-s)+2*s**2-2*s*s2+s2**2)+2*t1**2*(s-s2)+t1**3)+m**2*(2*M**4*((2*mA**2+t1)*(mA**2-s2+t1)-2*t2*(mA**2+t1))+2*M**2*((mA**2-s2+t1)*(3*mA**4-2*mA**2*(s+s1+s2-t1)-t1*(2*s+s2)+s2**2)+2*t2*(-2*mA**4+mA**2*(2*s+s2-2*t1)+2*s*t1+s2*(t1-s2)))+mA**6*(t1-4*s1)+mA**4*(4*s*(s1+t2)+4*s1*s2-4*s1*t1-s2*t1-4*s2*t2+3*t1**2+2*t1*t2)+mA**2*(2*s**2*(t1-2*t2)-2*s*(2*s1+t1)*(s2-t1)+4*s*s2*t2+t1*(s2**2-4*s2*(t1+t2)+3*t1**2))+t1*(2*s**2+2*s*(t1-s2)+(s2-t1)**2)*(-s2+t1-2*t2))+2*M**4*(mA**2-s2+t1-t2)*(mA**4-mA**2*(s2-t1+t2)-t1*t2)-2*M**2*(2*mA**2*s1*(mA**2-s2+t1)**2+t2*(mA**2-s2+t1)*(mA**4-2*mA**2*(s+s1)-t1*(2*s+s2)+s2**2)-t2**2*(mA**4+mA**2*(t1-2*s)-t1*(2*s+s2)+s2**2))+mA**2*(2*s1**2+2*s1*t1+t1**2)*(mA**2-s2+t1)**2-t2*(mA**2-s2+t1)*(mA**4*t1+2*mA**2*(2*s*s1-2*s1*s2+2*s1*t1+t1**2)+t1*(2*s**2+2*s*(t1-s2)+(s2-t1)**2))+t2**2*(mA**4*t1+2*mA**2*(t1*(s-s2)+(s-s2)**2+t1**2)+t1*(2*s**2+2*s*(t1-s2)+(s2-t1)**2)))
+        return prefactor * 4*(m**4*(2*M**4*(mA**2+t1)+2*M**2*(3*mA**4+mA**2*(3*t1-2*(s+s2))-t1*(2*s+s2)+s2**2)+2*(mA**3-mA*s)**2+t1*(mA**4+2*mA**2*(s2-s)+2*s**2-2*s*s2+s2**2)+2*t1**2*(s-s2)+t1**3)+m**2*(2*M**4*((2*mA**2+t1)*(mA**2-s2+t1)-2*t2*(mA**2+t1))+2*M**2*((mA**2-s2+t1)*(3*mA**4-2*mA**2*(s+s1+s2-t1)-t1*(2*s+s2)+s2**2)+2*t2*(-2*mA**4+mA**2*(2*s+s2-2*t1)+2*s*t1+s2*(t1-s2)))+mA**6*(t1-4*s1)+mA**4*(4*s*(s1+t2)+4*s1*s2-4*s1*t1-s2*t1-4*s2*t2+3*t1**2+2*t1*t2)+mA**2*(2*s**2*(t1-2*t2)-2*s*(2*s1+t1)*(s2-t1)+4*s*s2*t2+t1*(s2**2-4*s2*(t1+t2)+3*t1**2))+t1*(2*s**2+2*s*(t1-s2)+(s2-t1)**2)*(-s2+t1-2*t2))+2*M**4*(mA**2-s2+t1-t2)*(mA**4-mA**2*(s2-t1+t2)-t1*t2)-2*M**2*(2*mA**2*s1*(mA**2-s2+t1)**2+t2*(mA**2-s2+t1)*(mA**4-2*mA**2*(s+s1)-t1*(2*s+s2)+s2**2)-t2**2*(mA**4+mA**2*(t1-2*s)-t1*(2*s+s2)+s2**2))+mA**2*(2*s1**2+2*s1*t1+t1**2)*(mA**2-s2+t1)**2-t2*(mA**2-s2+t1)*(mA**4*t1+2*mA**2*(2*s*s1-2*s1*s2+2*s1*t1+t1**2)+t1*(2*s**2+2*s*(t1-s2)+(s2-t1)**2))+t2**2*(mA**4*t1+2*mA**2*(t1*(s-s2)+(s-s2)**2+t1**2)+t1*(2*s**2+2*s*(t1-s2)+(s2-t1)**2)))
+
+
+
+
+class M2DMFourFermi(MatrixElement2):
+    def __init__(self, mDM=1.0e3):
+        super().__init__(M_P, mDM, mDM, M_P)
+        self.mPsi = mDM
+
+    def __call__(self, s, t, lamByMx=0.1):
+        return power(lamByMx, 4) * ((s+t)**2 - 4*(M_P**2 - self.mPsi**2)**2)
+
+
+
+class M2NeutrinoBrem(MatrixElement2to3):
+    def __init__(self, mf=1.0, mi=0.0, MTarget=M_P, z=6, LamScale=1.0e6):
+        super().__init__(MTarget, mi, 0.0, mf, MTarget)
+        self.mf = mf
+        self.mi = mi
+        self.M = MTarget
+        self.ff = AtomicElasticFF(z)
+        self.Lam = LamScale
+    """
+    def __call__(self, s, s2, t1, s1, t2):
+        prefactor = self.ff(sqrt(abs(-2*self.M**2+s-s1+t2)))*(2*pi*ALPHA)/(self.Lam**6 * power(-2*self.M**2+s-s1+t2, 2))
+        return prefactor*abs((self.M**2+2*self.mf*self.mi-s2+t1-t2)*\
+                             (1/8*(s1-self.mf**2)*(-(self.mf**2+s-s1-s2)*(self.M**2+self.mi**2-s+s1+t1-t2)+(self.M**2-t1)\
+                                                   *(2*self.M**2+self.mf**2+self.mi**2-s2-t2)+4*self.M**2*(self.M**2-s+s2-t1)) \
+                                                    +1/8*(-self.M**2+s-s2+t1)*((self.mf**2+s-s1-s2)*(self.M**2+self.mi**2-s+s1+t1-t2) \
+                                                                               -(self.M**2-t1)*(2*self.M**2+self.mf**2+self.mi**2-s2-t2)) \
+                                                                                +1/4*self.mi*(self.M**2-t1)*(self.mf**2+s-s1-s2)+1/4*(self.M**2-t1)\
+                                                                                    *(self.M**2-s2+t1-t2)*(self.mf**2+s-s1-s2)+1/4*self.mf**2 \
+                                                                                        *(self.M**2-t1)*(self.mf**2+s-s1-s2)+1/4*self.M**2*(self.mf**2-s1)**2+1/4*self.M**2*(self.M**2-s+s2-t1)**2))
+    """
+    def __call__(self, s, s2, t1, s1, t2):
+        prefactor = self.ff(sqrt(abs(-2*self.M**2+s-s1+t2)))*(2*pi*ALPHA)/(self.Lam**6 * power(-2*self.M**2+s-s1+t2, 2))
+        return prefactor*abs((self.M**2 - s2 + t1 - t2)*(2*self.M**2 * (2*self.mf**2*(2*s - 2*(s1 + s2) + t1) + 2*self.mf**4 + 2*s**2 - \
+                                                        4*s*(s1 + s2) + 3*s*t1 + 2*s1**2 + 4*s1*s2 - 3*s1*t1 + 2*s2**2 - \
+                                                        2*s2*t1 + 2*t1**2 + t1*t2) - self.M**4 * (4*self.mf**2 + 5*s - 5*s1 - 4*s2 + 8*t1 + t2) + \
+                                                            4*self.M**6 - (s - s1 + t2)*((self.mf**2 + s - s1 - s2)**2 + t1**2)))
