SND@LHC Software
Loading...
Searching...
No Matches
Functions | Variables
proton_bremsstrahlung Namespace Reference

Functions

 rhoFormFactor (m)
 
 energy (p, m)
 
 penaltyFactor (m)
 
 zeta (p, theta)
 
 pTransverse (p, theta)
 
 ptSquare (p, theta)
 
 H (p, theta, mDarkPhoton)
 
 wba (p, theta, mDarkPhoton, epsilon)
 
 sigma (s)
 
 es (p, mDarkPhoton)
 
 sigmaRatio (p, mDarkPhoton)
 
 dNdZdPtSquare (p, mDarkPhoton, theta, epsilon)
 
 dPt2dTheta (p, theta)
 
 dZdP (p, theta)
 
 dNdPdTheta (p, theta, mDarkPhoton, epsilon)
 
 pMin (mDarkPhoton)
 
 pMax (mDarkPhoton)
 
 prodRate (mDarkPhoton, epsilon, tmin=-0.5 *math.pi, tmax=0.5 *math.pi)
 
 normalisedProductionPDF (p, theta, mDarkPhoton, epsilon, norm)
 
 hProdPDF (mDarkPhoton, epsilon, norm, binsp, binstheta, tmin=-0.5 *math.pi, tmax=0.5 *math.pi, suffix="")
 

Variables

float mProton = 0.938272081
 
int protonEnergy = 400.
 
 protonMomentum = math.sqrt(protonEnergy*protonEnergy - mProton*mProton)
 

Function Documentation

◆ dNdPdTheta()

proton_bremsstrahlung.dNdPdTheta (   p,
  theta,
  mDarkPhoton,
  epsilon 
) 
 Differential A' rate per p.o.t. as a function of P and theta

Definition at line 125 of file proton_bremsstrahlung.py.

125def dNdPdTheta(p, theta, mDarkPhoton, epsilon):
126 """ Differential A' rate per p.o.t. as a function of P and theta """
127 diffRate = dNdZdPtSquare(p,mDarkPhoton,theta,epsilon) * dPt2dTheta(p,theta) * dZdP(p,theta)
128 return math.fabs(diffRate) # integrating in (-pi, pi)...
129
130

◆ dNdZdPtSquare()

proton_bremsstrahlung.dNdZdPtSquare (   p,
  mDarkPhoton,
  theta,
  epsilon 
) 
 Differential A' rate per p.o.t. as a function of Z and Pt^2

Definition at line 109 of file proton_bremsstrahlung.py.

109def dNdZdPtSquare(p, mDarkPhoton, theta, epsilon):
110 """ Differential A' rate per p.o.t. as a function of Z and Pt^2 """
111 return sigmaRatio(p,mDarkPhoton)*wba(p,theta,mDarkPhoton,epsilon)
112
113

◆ dPt2dTheta()

proton_bremsstrahlung.dPt2dTheta (   p,
  theta 
) 
 Jacobian Pt^2->theta

Definition at line 114 of file proton_bremsstrahlung.py.

114def dPt2dTheta(p, theta):
115 """ Jacobian Pt^2->theta """
116 z2 = pow(zeta(p,theta),2.)
117 return 2.*theta*z2*protonMomentum*protonMomentum
118
119

◆ dZdP()

proton_bremsstrahlung.dZdP (   p,
  theta 
) 
 Jacobian z->p

Definition at line 120 of file proton_bremsstrahlung.py.

120def dZdP(p, theta):
121 """ Jacobian z->p """
122 return 1./( protonMomentum* math.sqrt(theta*theta+1.) )
123
124

◆ energy()

proton_bremsstrahlung.energy (   p,
  m 
) 
 Compute energy from momentum and mass

Definition at line 34 of file proton_bremsstrahlung.py.

34def energy(p,m):
35 """ Compute energy from momentum and mass """
36 return math.sqrt(p*p + m*m)
37

◆ es()

proton_bremsstrahlung.es (   p,
  mDarkPhoton 
) 
 s(p,mA)

Definition at line 99 of file proton_bremsstrahlung.py.

99def es(p, mDarkPhoton):
100 """ s(p,mA) """
101 return 2.*mProton*(energy(protonMomentum,mProton)-energy(p,mDarkPhoton))
102
103

◆ H()

proton_bremsstrahlung.H (   p,
  theta,
  mDarkPhoton 
) 
 A kinematic term

Definition at line 61 of file proton_bremsstrahlung.py.

61def H(p, theta, mDarkPhoton):
62 """ A kinematic term """
63 return ptSquare(p,theta) + (1.-zeta(p,theta))*mDarkPhoton*mDarkPhoton + pow(zeta(p,theta),2.)*mProton*mProton
64
65

◆ hProdPDF()

proton_bremsstrahlung.hProdPDF (   mDarkPhoton,
  epsilon,
  norm,
  binsp,
  binstheta,
  tmin = -0.5 * math.pi,
  tmax = 0.5 * math.pi,
  suffix = "" 
) 
 Histogram of the PDF for A' production in SHIP

Definition at line 165 of file proton_bremsstrahlung.py.

165def hProdPDF(mDarkPhoton, epsilon, norm, binsp, binstheta, tmin = -0.5 * math.pi, tmax = 0.5 * math.pi, suffix=""):
166 """ Histogram of the PDF for A' production in SHIP """
167 angles = np.linspace(tmin,tmax,binstheta).tolist()
168 anglestep = 2.*(tmax - tmin)/binstheta
169 momentumStep = (pMax(mDarkPhoton)-pMin(mDarkPhoton))/(binsp-1)
170 momenta = np.linspace(pMin(mDarkPhoton),pMax(mDarkPhoton),binsp,endpoint=False).tolist()
171 hPDF = r.TH2F("hPDF_eps%s_m%s"%(epsilon,mDarkPhoton) ,"hPDF_eps%s_m%s"%(epsilon,mDarkPhoton),
172 binsp,pMin(mDarkPhoton)-0.5*momentumStep,pMax(mDarkPhoton)-0.5*momentumStep,
173 binstheta,tmin-0.5*anglestep,tmax-0.5*anglestep)
174 hPDF.SetTitle("PDF for A' production (m_{A'}=%s GeV, #epsilon =%s)"%(mDarkPhoton,epsilon))
175 hPDF.GetXaxis().SetTitle("P_{A'} [GeV]")
176 hPDF.GetYaxis().SetTitle("#theta_{A'} [rad]")
177 hPDFtheta = r.TH1F("hPDFtheta_eps%s_m%s"%(epsilon,mDarkPhoton),
178 "hPDFtheta_eps%s_m%s"%(epsilon,mDarkPhoton),
179 binstheta,tmin-0.5*anglestep,tmax-0.5*anglestep)
180 hPDFp = r.TH1F("hPDFp_eps%s_m%s"%(epsilon,mDarkPhoton),
181 "hPDFp_eps%s_m%s"%(epsilon,mDarkPhoton),
182 binsp,pMin(mDarkPhoton)-0.5*momentumStep,pMax(mDarkPhoton)-0.5*momentumStep)
183 hPDFp.GetXaxis().SetTitle("P_{A'} [GeV]")
184 hPDFtheta.GetXaxis().SetTitle("#theta_{A'} [rad]")
185 for theta in angles:
186 for p in momenta:
187 w = normalisedProductionPDF(p,theta,mDarkPhoton,epsilon,norm)
188 hPDF.Fill(p,theta,w)
189 hPDFtheta.Fill(theta,w)
190 hPDFp.Fill(p,w)
191 hPdfFilename = sys.modules['__main__'].outputDir+"/ParaPhoton_eps%s_m%s%s.root"%(epsilon,mDarkPhoton,suffix)
192 outfile = r.TFile(hPdfFilename,"recreate")
193 #weight = hPDF.Integral("width")
194 #print "Weight = %3.3f"%weight
195 #hPDF.Scale(1./weight)
196 hPDF.Write()
197 hPDFp.Write()
198 hPDFtheta.Write()
199 outfile.Close()
200 del angles
201 del momenta
202 return hPDF
203
204

◆ normalisedProductionPDF()

proton_bremsstrahlung.normalisedProductionPDF (   p,
  theta,
  mDarkPhoton,
  epsilon,
  norm 
) 
 Probability density function for A' production in SHIP

Definition at line 160 of file proton_bremsstrahlung.py.

160def normalisedProductionPDF(p, theta, mDarkPhoton, epsilon, norm):
161 """ Probability density function for A' production in SHIP """
162 return (1. / norm) * dNdPdTheta(p, theta, mDarkPhoton, epsilon)
163
164

◆ penaltyFactor()

proton_bremsstrahlung.penaltyFactor (   m) 
 Penalty factor for high masses - dipole form factor in the proton-A' vertex 
 m in GeV

Definition at line 38 of file proton_bremsstrahlung.py.

38def penaltyFactor(m):
39 """ Penalty factor for high masses - dipole form factor in the proton-A' vertex """
40 """ m in GeV """
41 if m*m>0.71:
42 return math.pow(m*m/0.71,-4)
43 else:
44 return 1
45

◆ pMax()

proton_bremsstrahlung.pMax (   mDarkPhoton)

Definition at line 135 of file proton_bremsstrahlung.py.

135def pMax(mDarkPhoton):
136 #return min(0.86*protonMomentum, math.sqrt( (energy(protonMomentum,mProton)**2. - mDarkPhoton**2.) - mDarkPhoton**2.))
137 return math.sqrt( (energy(protonMomentum,mProton)**2. - mDarkPhoton**2.) - mDarkPhoton**2.)
138
139

◆ pMin()

proton_bremsstrahlung.pMin (   mDarkPhoton)

Definition at line 131 of file proton_bremsstrahlung.py.

131def pMin(mDarkPhoton):
132 return max(0.14*protonMomentum, mDarkPhoton)
133
134

◆ prodRate()

proton_bremsstrahlung.prodRate (   mDarkPhoton,
  epsilon,
  tmin = -0.5 * math.pi,
  tmax = 0.5 * math.pi 
) 
 dNdPdTheta integrated over p and theta

Definition at line 140 of file proton_bremsstrahlung.py.

140def prodRate(mDarkPhoton, epsilon, tmin = -0.5 * math.pi, tmax = 0.5 * math.pi):
141 """ dNdPdTheta integrated over p and theta """
142 integral = dblquad( dNdPdTheta, # integrand
143 tmin, tmax, # theta boundaries (2nd argument of integrand)
144 lambda x: pMin(mDarkPhoton), lambda x: pMax(mDarkPhoton), # p boundaries (1st argument of integrand)
145 args=(mDarkPhoton, epsilon) ) # extra parameters to pass to integrand
146 return integral[0]
147
148# total production rate of A'
149#norm = prodRate(1.1,3.e-7) #mDarkPhoton,epsilon)
150# number of A' produced
151# numDarkPhotons = int(math.floor(norm*protonFlux))
152#
153# print
154# print "Epsilon \t %s"%epsilon
155# print "mDarkPhoton \t %s"%mDarkPhoton
156#print "A' production rate per p.o.t: \t %.8g"%norm
157# print "Number of A' produced in SHiP: \t %.8g"%numDarkPhotons
158
159

◆ pTransverse()

proton_bremsstrahlung.pTransverse (   p,
  theta 
) 
 Paraphoton transverse momentum in the lab frame

Definition at line 51 of file proton_bremsstrahlung.py.

51def pTransverse(p, theta):
52 """ Paraphoton transverse momentum in the lab frame """
53 return protonMomentum*theta*zeta(p,theta)
54
55

◆ ptSquare()

proton_bremsstrahlung.ptSquare (   p,
  theta 
) 
 Square paraphoton transverse momentum in the lab frame

Definition at line 56 of file proton_bremsstrahlung.py.

56def ptSquare(p, theta):
57 """ Square paraphoton transverse momentum in the lab frame """
58 return pow(pTransverse(p,theta), 2.)
59
60

◆ rhoFormFactor()

proton_bremsstrahlung.rhoFormFactor (   m) 
 From https://arxiv.org/abs/0910.5589

Definition at line 16 of file proton_bremsstrahlung.py.

16def rhoFormFactor(m):
17 """ From https://arxiv.org/abs/0910.5589 """
18 #constants from the code from Inar: https://github.com/pgdeniverville/BdNMC/blob/master/src/Proton_Brem_Distribution.cpp
19 f1ra = 0.6165340033101271
20 f1rb = 0.22320420111672623
21 f1rc = -0.33973820442685326
22 f1wa = 1.0117544786579074
23 f1wb = -0.8816565944110686
24 f1wc = 0.3699021157531611
25 f1prho = f1ra*0.77**2/(0.77**2-m**2-0.77*0.15j)
26 f1prhop = f1rb*1.25**2/(1.25**2-m**2-1.25*0.3j)
27 f1prhopp = f1rc*1.45**2/(1.45**2-m**2-1.45*0.5j)
28 f1pomega = f1wa*0.77**2/(0.77**2-m**2-0.77*0.0085j)
29 f1pomegap = f1wb*1.25**2/(1.25**2-m**2-1.25*0.3j)
30 f1pomegapp = f1wc*1.45**2/(1.45**2-m**2-1.45*0.5j)
31 return abs(f1prho+f1prhop+f1prhopp+f1pomega+f1pomegap+f1pomegapp)
32
33# useful functions

◆ sigma()

proton_bremsstrahlung.sigma (   s) 
 Parametrisation of sigma(s)

Definition at line 84 of file proton_bremsstrahlung.py.

84def sigma(s): # s in GeV^2 ---> sigma in mb
85 """ Parametrisation of sigma(s) """
86 a1 = 35.45
87 a2 = 0.308
88 a3 = 28.94
89 a4 = 33.34
90 a5 = 0.545
91 a6 = 0.458
92 a7 = 42.53
93 p1 = a2*pow(math.log(s/a3),2.)
94 p2 = a4*pow((1./s),a5)
95 p3 = a7*pow((1./s),a6)
96 return a1 + p1 - p2 + p3
97
98

◆ sigmaRatio()

proton_bremsstrahlung.sigmaRatio (   p,
  mDarkPhoton 
) 
 sigma(s') / sigma(s)

Definition at line 104 of file proton_bremsstrahlung.py.

104def sigmaRatio(p, mDarkPhoton):
105 """ sigma(s') / sigma(s) """
106 return sigma(es(p,mDarkPhoton)) / sigma(2.*mProton*energy(protonMomentum,mProton))
107
108

◆ wba()

proton_bremsstrahlung.wba (   p,
  theta,
  mDarkPhoton,
  epsilon 
) 
 Cross section weighting function in the Fermi-Weizsaeker-Williams approximation

Definition at line 66 of file proton_bremsstrahlung.py.

66def wba(p, theta, mDarkPhoton, epsilon):
67 """ Cross section weighting function in the Fermi-Weizsaeker-Williams approximation """
68 const = epsilon*epsilon*alphaQED / (2.*math.pi*H(p,theta,mDarkPhoton))
69
70 h2 = pow(H(p,theta,mDarkPhoton),2.)
71 oneMinusZSquare = pow(1.-zeta(p,theta),2.)
72 mp2 = mProton*mProton
73 mA2 = mDarkPhoton*mDarkPhoton
74
75 p1 = (1. + oneMinusZSquare) / zeta(p,theta)
76 p2 = ( 2. * zeta(p,theta) * (1.-zeta(p,theta)) * ( (2.*mp2 + mA2)/ H(p,theta,mDarkPhoton)
77 - pow(zeta(p,theta),2.)*2.*mp2*mp2/h2 ) )
78 #p3 = 2.*zeta(p,theta)*(1.-zeta(p,theta))*(zeta(p,theta)+oneMinusZSquare)*mp2*mA2/h2
79 p3 = 2.*zeta(p,theta)*(1.-zeta(p,theta))*(1+oneMinusZSquare)*mp2*mA2/h2
80 p4 = 2.*zeta(p,theta)*oneMinusZSquare*mA2*mA2/h2
81 return const*(p1-p2+p3+p4)
82
83

◆ zeta()

proton_bremsstrahlung.zeta (   p,
  theta 
) 
 Fraction of the proton momentum carried away by the paraphoton in the beam direction

Definition at line 46 of file proton_bremsstrahlung.py.

46def zeta(p, theta):
47 """ Fraction of the proton momentum carried away by the paraphoton in the beam direction """
48 return p / (protonMomentum * math.sqrt(theta*theta + 1.))
49
50

Variable Documentation

◆ mProton

float proton_bremsstrahlung.mProton = 0.938272081

Definition at line 11 of file proton_bremsstrahlung.py.

◆ protonEnergy

int proton_bremsstrahlung.protonEnergy = 400.

Definition at line 12 of file proton_bremsstrahlung.py.

◆ protonMomentum

proton_bremsstrahlung.protonMomentum = math.sqrt(protonEnergy*protonEnergy - mProton*mProton)

Definition at line 13 of file proton_bremsstrahlung.py.