1
#include "DD_Nucleus.hpp"
2

3
#include <cmath>
4

5
#include "Physics_Functions.hpp"
6
#include "DD_General.hpp"
7

8
//1. Velocity necessary for a DM particle of mass mChi to recoil a nucleus of mass number A with a given recoil energy.
9 2
	double vMinimal_N(double RecoilEnergy,double mChi,double A)
10
	{
11 2
		return sqrt(mNucleon*A*RecoilEnergy/2.0/pow(Reduced_Mass(mChi,A*mNucleon),2.0));
12
	}
13
//2. Maximum recoil energy a DM particle of mass mChi and speed v can cause on a nucleus of mass number A
14 2
	double ERMax(double v,double mChi,double A)
15
	{
16 2
		return 2.0*v*v*pow(Reduced_Mass(mChi,A*mNucleon),2.0)/A/mNucleon;
17
	}
18

19
//3. Direct detection recoil spectrum
20 2
	double dRdER(double ER,const DM_Particle& DM,double X,int Z,double A)
21
	{
22 2
		double mA = A*mNucleon;
23 2
		double nDM = rhoDM/DM.mass;
24 2
		double v = 1.0; //Cancels in the following expression
25 2
		return X/mA*nDM*(v*v*DM.dSdER(ER,Z,A,v))*EtaFunction(vMinimal_N(ER,DM.mass,A),vEarth);
26
	}
27 2
	Interpolation dRdER(const DM_Particle& DM,double Emin,double Emax,double X,int Z,double A,const std::vector<double>& attenuation,const std::vector<DataPoint> &speeddata)
28
	{
29
		//0. maximal energy
30 2
			double ERmax = ERMax((vesc+vEarth),DM.mass,A);
31 2
			ERmax = std::min(ERmax,Emax); //only calculate points where dRdER !=0 for the interpolation
32
		// Check if the DM is able to cause recoils
33 2
			double vCutoff=vMinimal_N(Emin,DM.mass,A);
34 2
			if(vCutoff>(vesc+vEarth))
35
			{
36 2
				int points = 10;
37 2
				double dER = (Emax-Emin)/(points-1.0);
38 2
				std::vector<std::vector<double>> interpol_list;
39 2
				for(int i = 0;i<points;i++)
40
				{
41 2
					double ER = Emin+ i*dER;
42 2
					interpol_list.push_back(std::vector<double> {ER,0.0});
43
				}
44 2
				return Interpolation(interpol_list);
45
			}
46
		//1. Prefactor
47 2
			double mA = A * mNucleon;
48 2
			double prefactor = X/mA*rhoDM/DM.mass*attenuation[0];
49
		//2. Find the Kernel density estimate for the speed distribution function
50 2
			Interpolation kde = Perform_KDE(speeddata,vCutoff,(vesc+vEarth));
51
		//3. Create list to interpolate
52 2
			int points = 200;
53 2
			double dER = (ERmax-Emin)/(points-1.0);
54 2
			std::vector<std::vector<double>> interpol_list;
55 2
			for(int i = 0;i<points;i++)
56
			{
57 2
				double ER = Emin+ i*dER;
58
				//3.1 Create integrand.
59 2
					std::function<double(double)> integrand = [ER,DM,Z,A,&kde](double v)
60
					{
61 2
						return v*kde(v)*DM.dSdER(ER,Z,A,v);
62 2
					};
63
				//3.2 Integrate.
64 2
					double vMin = vMinimal_N(ER,DM.mass,A);
65
					double integral;
66 2
					if(vMin>=(vesc+vEarth))
67
					{
68 2
						integral = 0.0;
69
					}
70
					else
71
					{
72
				  		//Integrate
73
				  		//NOTE: Gives precision warnings regularly due to a too small choice of epsilon.
74 2
				  			double epsilon = Find_Epsilon(integrand,vMin,(vesc+vEarth),1.0e-3);
75 2
							integral = Integrate(integrand,vMin,(vesc+vEarth),epsilon);
76
	
77
					}
78
				//3.3 Append to list 
79 2
					interpol_list.push_back(std::vector<double> {ER,integral});
80
			}
81
			//3.4 If ERMax < Emax, fill the rest of the interval with zeros, so that we have a interpolation on the full interval.
82 2
				if(ERmax<Emax)
83
				{
84 2
					points = 5;
85 2
					dER = (Emax-ERmax)/(points-1);
86 2
					for(int i = 1;i<points;i++)
87
					{
88 2
						double ER = ERmax + i*dER;
89 2
						interpol_list.push_back(std::vector<double> {ER,0.0});
90
					}
91
				}
92
		//4. Interpolate and include prefactor.
93 2
			Interpolation drder(interpol_list);
94 2
			drder.Multiply(prefactor);
95
		
96 2
		return drder;
97
	}

Read our documentation on viewing source code .

Loading