1
#include "DM_Particle.hpp"
2

3
#include <iostream>
4
#include <cmath>
5
#include <functional>
6

7
#include "Numerics_Functions.hpp"
8
#include "Physics_Functions.hpp"
9

10
//6. DM-nucleus scattering cross sections:
11
	//Nuclear Helm form factor
12 2
		double FormFactor_N(double q,double A,bool ldm)
13
		{
14 2
			if(ldm || q==0.0) return 1.0;
15
			else
16
			{
17 0
				double a = 0.52*fm;
18 0
				double c = (1.23*pow(A,1.0/3.0)-0.6)*fm;
19 0
				double s = 0.9*fm;
20 0
				double rn = sqrt(c*c+7.0/3.0*pow(M_PI*a,2.0)-5*s*s);
21 0
				double qr = q*rn;
22 0
				return 3.0*(sin(qr)/pow(qr,3.0)-cos(qr)/pow(qr,2.0))*exp(-q*q*s*s/2.0);
23
			}
24
		}
25
		
26
	//3. DM particle struct
27 2
		DM_Particle::DM_Particle()
28
		{
29 2
			mass=0.0;
30 2
			sigma_n=0.0;
31 2
			sigma_e = 0.0;
32 2
			ldm = true;
33 2
			formfactor=" ";
34 2
			ZorA = " ";
35 2
			screening=false;
36 2
			mMediator=0.0;
37
		}
38

39 2
		DM_Particle::DM_Particle(double mDM,double sn,double se,bool light,std::string ff,std::string za,bool scr,double mMed)
40
		{
41 2
			mass=mDM;
42 2
			sigma_n=sn;
43 2
			sigma_e = se;
44 2
			ldm = light;
45 2
			formfactor=ff;
46 2
			ZorA = za;
47 2
			screening=scr;
48 2
			mMediator=mMed;
49
		}
50
		//Atomic form factor (charge screening)
51 2
		double DM_Particle::FormFactor_A(double q,int Z) const
52
		{
53 2
			if(!screening) return 1.0;
54
			else
55
			{
56 0
				double a2q2 = pow(q*Thomas_Fermi_Radius(Z),2.0);
57 0
				return a2q2/(1.0+a2q2);
58
			}
59
		}
60 2
		void DM_Particle::Set_Mass(double m)
61
		{
62 2
			mass = m;
63
		}
64 2
		void DM_Particle::Set_Sigma_n(double s)
65
		{
66 2
			sigma_n = s;
67 2
			sigma_e = pow(Reduced_Mass(mass,mElectron)/Reduced_Mass(mass,mProton),2.0)*sigma_n;
68
		}
69 2
		void DM_Particle::Set_Sigma_e(double s)
70
		{
71 2
			sigma_e = s;
72 2
			sigma_n = pow(Reduced_Mass(mass,mProton)/Reduced_Mass(mass,mElectron),2.0)*sigma_e;
73
		}
74
	//DM form factor
75 2
		double DM_Particle::FormFactor(double q) const 
76
		{
77
			//Contact interactions
78 2
				if(formfactor=="Contact") return 1.0;
79
			//General dark photon
80 0
				else if(formfactor=="General")	return (qRef*qRef+mMediator*mMediator)/(q*q+mMediator*mMediator);
81
			//Long range interaction
82 0
				else if (formfactor=="Long-Range")	return qRef*qRef/q/q;
83
			//Electric dipole interaction
84 0
				else if (formfactor=="Electric-Dipole")		return qRef/q;
85
			//Error
86
				else
87
				{
88 0
					std::cerr <<"Error in FormFactor(): Form factor "<<formfactor <<"not recognized."<<endl;
89 0
					std::exit(EXIT_FAILURE);
90
				}
91
		}
92
	//Zero momentum transfer spin-independent cross-section
93 2
		double DM_Particle::sigmaSI(int Z,double A) const
94
		{
95 2
			double X=-1.0;
96 2
			double mAux=-1.0;
97 2
			if(ZorA=="Z") 
98
			{
99 2
				X=Z;
100 2
				mAux = mProton;
101
			}
102 2
			else if (ZorA=="A") 
103
			{
104 2
				X=A;
105 2
				mAux = mNucleon;
106
			}
107
			else 
108
			{
109 0
				std::cerr <<"Error in sigmaSI: ZorA = " <<ZorA <<" not recognized."<<endl;
110 0
				std::exit(EXIT_FAILURE);
111
			}
112 2
			return sigma_n*pow(Reduced_Mass(mass,NucleusMass(A)),2.0)/pow(Reduced_Mass(mass,mAux),2.0)*pow(X,2.0);
113
		}
114
	//Differential cross-sections
115 2
		double DM_Particle::dSdER(double ER,int Z,double A,double vDM) const
116
		{
117 2
			double mA = A*mNucleon;
118 2
			double ERmax = 2.0*pow(Reduced_Mass(mass,mA)*vDM,2.0)/mA;
119 2
			double q = sqrt(2.0*mA*ER);	
120 2
			if(q==0&&(formfactor=="Electric-Dipole"||formfactor=="Long-Range")&&screening) q=1e-15*keV; //Screening cause 1/q to cancel
121 2
			return sigmaSI(Z,A) / ERmax * pow(FormFactor_N(q,A,ldm),2.0) * pow(FormFactor(q),2.0) * pow(FormFactor_A(q,Z),2.0); 
122
		}
123 0
		double DM_Particle::dSdq2(double q,int Z,double A,double vDM) const
124
		{
125 0
			double mA = A*mNucleon;
126 0
			double qMax2 = 4.0*pow(Reduced_Mass(mass,mA)*vDM,2.0);
127 0
			if(q==0&&(formfactor=="Electric-Dipole"||formfactor=="Long-Range")&&screening) q=1e-15*keV; //Screening cause 1/q to cancel
128 0
			return sigmaSI(Z,A) / qMax2 * pow(FormFactor_N(q,A,ldm),2.0) * pow(FormFactor(q),2.0)* pow(FormFactor_A(q,Z),2.0); 
129
		}
130 2
		double DM_Particle::Sigma_Tot(int Z,double A,double vDM) const
131
		{
132 2
			if(!ldm)
133
			{
134 0
				double mA = A*mNucleon;
135 0
				double ERmax = 2.0*pow(Reduced_Mass(mass,mA)*vDM,2.0)/mA;
136 0
				double ERmin = 0.0;
137
				//integrate the diff cross section
138 0
					auto dodER = std::bind(&DM_Particle::dSdER,this,std::placeholders::_1,Z,A,vDM);
139
					//Numerical integration
140 0
			  			double epsilon = Find_Epsilon(dodER,ERmin,ERmax,1e-5);
141 0
			  			double integral =Integrate(dodER,ERmin,ERmax,epsilon);
142 0
					return integral;
143
			}
144
			else
145
			{
146 2
				double result = sigmaSI(Z,A);
147
				
148 2
				double q2max = pow(2.0*Reduced_Mass(mass,A*mNucleon)*vDM,2.0);
149 2
				double a2=pow(Thomas_Fermi_Radius(Z),2.0);
150
				//General interaction
151 2
					if(formfactor == "General")
152
					{
153 0
						double m2 = pow(mMediator,2.0);
154 0
						if(screening)
155
						{
156 0
							result*=(pow(a2,2.0)*pow(m2 + qRef*qRef,2.0)*(((a2*m2-1.0)*q2max*(q2max + m2*(2.0 + a2*q2max)))/((m2 + q2max)*(1.0 + a2*q2max)) + 2.0*m2*log((m2 + q2max)/(m2 + a2*m2*q2max))))/(pow(a2*m2-1.0,3.0)*q2max);
157
						}
158
						else
159
						{
160 0
							result*= pow(qRef*qRef+m2,2.0)/m2/(m2+q2max);
161
						}
162
					}
163
				//Contact interaction
164 2
					else if(formfactor == "Contact" && screening)
165
					{
166 0
						result *=(1.0+1.0/(1+a2*q2max)-2.0/a2/q2max*log1p(a2*q2max));
167
					}
168
				//Electric dipole interaction
169 2
					else if(formfactor == "Electric-Dipole")
170
					{
171 0
						result *= pow(qRef,2.0)/q2max*(log1p(a2*q2max)-a2*q2max/(1+a2*q2max));
172
					}
173
				//Long range interaction
174 2
					else if(formfactor== "Long-Range")
175
					{
176 0
						result *= pow(a2,2.0)*pow(qRef,4.0)/(1+a2*q2max);
177
					}
178

179 2
				return result;
180
			}
181
		}
182

183
	//Find the stopping cross section
184 2
		double Stopping_CrossSection(const std::string& detector,double sigma,double mDM)
185
		{
186 2
			if(detector == "Semiconductor"||detector=="SENSEI-surface"||detector=="SENSEI"||detector=="SuperCDMS"||detector=="DAMIC-M"||detector == "XENON10e"||detector == "XENON100e"||detector =="DarkSide-50")return pow(Reduced_Mass(mDM,mProton)/Reduced_Mass(mDM,mElectron),2.0)*sigma;
187 2
			else return sigma;
188
		}
189

190
//7. DM speed distribution
191 2
	double SpeedDistribution(double v,double vEarth)
192
	{
193 2
		return M_PI*v*v0*v0/Nesc/vEarth*(2*exp(-(v*v+vEarth*vEarth)/v0/v0)*sinh(2*v*vEarth/v0/v0)+(exp(-pow(v+vEarth,2.0)/v0/v0)-exp(-vesc*vesc/v0/v0))*StepFunction(abs(v+vEarth)-vesc)-(exp(-pow(v-vEarth,2.0)/v0/v0)-exp(-vesc*vesc/v0/v0))*StepFunction(abs(v-vEarth)-vesc) );
194
	}
195 2
 	double Average_Speed(double vEarth,double vMin)
196
 	{
197
 		//1. integrand
198 2
	 		std::function<double(double)> integrand = [vEarth] (double v)
199
	 		{
200 2
	 			return v*SpeedDistribution(v,vEarth);
201 2
	 		};
202
	 	//2. integrate
203 2
	 		double vMean = Integrate(integrand,vMin,(vesc+vEarth),1e-6);
204
	 	//3. Renormalize
205 2
	 		if(vMin>0.0)
206
	 		{
207
	 			
208 2
	 			std::function<double(double)> integrand2 = [vEarth] (double v)
209
		 		{
210
		 			return SpeedDistribution(v,vEarth);
211 2
		 		};
212 2
		 		double norm=Integrate(integrand2,vMin,vesc+vEarth,1e-6);
213 2
		 		vMean/=norm;
214
	 		}
215 2
	 		return vMean;
216
 	}

Read our documentation on viewing source code .

Loading