Showing 15 of 76 files from the diff.
Newly tracked file
src/DD_Nucleus.cpp created.
Newly tracked file
src/DD_General.cpp created.
Newly tracked file
src/DM_Particle.cpp created.
Newly tracked file
src/DD_Electron.cpp created.
Other files ignored by Codecov
Plots.nb has changed.
Makefile has changed.
bin/README.md was deleted.
.gitignore has changed.
README.md has changed.
build/README.md has changed.
.travis.yml has changed.

@@ -0,0 +1,97 @@
Loading
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 +
	double vMinimal_N(double RecoilEnergy,double mChi,double A)
10 +
	{
11 +
		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 +
	double ERMax(double v,double mChi,double A)
15 +
	{
16 +
		return 2.0*v*v*pow(Reduced_Mass(mChi,A*mNucleon),2.0)/A/mNucleon;
17 +
	}
18 +
19 +
//3. Direct detection recoil spectrum
20 +
	double dRdER(double ER,const DM_Particle& DM,double X,int Z,double A)
21 +
	{
22 +
		double mA = A*mNucleon;
23 +
		double nDM = rhoDM/DM.mass;
24 +
		double v = 1.0; //Cancels in the following expression
25 +
		return X/mA*nDM*(v*v*DM.dSdER(ER,Z,A,v))*EtaFunction(vMinimal_N(ER,DM.mass,A),vEarth);
26 +
	}
27 +
	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 +
			double ERmax = ERMax((vesc+vEarth),DM.mass,A);
31 +
			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 +
			double vCutoff=vMinimal_N(Emin,DM.mass,A);
34 +
			if(vCutoff>(vesc+vEarth))
35 +
			{
36 +
				int points = 10;
37 +
				double dER = (Emax-Emin)/(points-1.0);
38 +
				std::vector<std::vector<double>> interpol_list;
39 +
				for(int i = 0;i<points;i++)
40 +
				{
41 +
					double ER = Emin+ i*dER;
42 +
					interpol_list.push_back(std::vector<double> {ER,0.0});
43 +
				}
44 +
				return Interpolation(interpol_list);
45 +
			}
46 +
		//1. Prefactor
47 +
			double mA = A * mNucleon;
48 +
			double prefactor = X/mA*rhoDM/DM.mass*attenuation[0];
49 +
		//2. Find the Kernel density estimate for the speed distribution function
50 +
			Interpolation kde = Perform_KDE(speeddata,vCutoff,(vesc+vEarth));
51 +
		//3. Create list to interpolate
52 +
			int points = 200;
53 +
			double dER = (ERmax-Emin)/(points-1.0);
54 +
			std::vector<std::vector<double>> interpol_list;
55 +
			for(int i = 0;i<points;i++)
56 +
			{
57 +
				double ER = Emin+ i*dER;
58 +
				//3.1 Create integrand.
59 +
					std::function<double(double)> integrand = [ER,DM,Z,A,&kde](double v)
60 +
					{
61 +
						return v*kde(v)*DM.dSdER(ER,Z,A,v);
62 +
					};
63 +
				//3.2 Integrate.
64 +
					double vMin = vMinimal_N(ER,DM.mass,A);
65 +
					double integral;
66 +
					if(vMin>=(vesc+vEarth))
67 +
					{
68 +
						integral = 0.0;
69 +
					}
70 +
					else
71 +
					{
72 +
				  		//Integrate
73 +
				  		//NOTE: Gives precision warnings regularly due to a too small choice of epsilon.
74 +
				  			double epsilon = Find_Epsilon(integrand,vMin,(vesc+vEarth),1.0e-3);
75 +
							integral = Integrate(integrand,vMin,(vesc+vEarth),epsilon);
76 +
	
77 +
					}
78 +
				//3.3 Append to list 
79 +
					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 +
				if(ERmax<Emax)
83 +
				{
84 +
					points = 5;
85 +
					dER = (Emax-ERmax)/(points-1);
86 +
					for(int i = 1;i<points;i++)
87 +
					{
88 +
						double ER = ERmax + i*dER;
89 +
						interpol_list.push_back(std::vector<double> {ER,0.0});
90 +
					}
91 +
				}
92 +
		//4. Interpolate and include prefactor.
93 +
			Interpolation drder(interpol_list);
94 +
			drder.Multiply(prefactor);
95 +
		
96 +
		return drder;
97 +
	}

@@ -0,0 +1,57 @@
Loading
1 +
//Contains functions related with direct detection rates for DM-Electron experiments.
2 +
#ifndef __DD_Electron_hpp_
3 +
#define __DD_Electron_hpp_
4 +
5 +
#include <string>
6 +
#include <vector>
7 +
8 +
#include "DM_Particle.hpp"
9 +
#include "Numerics_Functions.hpp"
10 +
11 +
//1. Import crystal/ionization form factor
12 +
	extern void Import_FormFactor(const std::string& target);
13 +
14 +
//2. Semiconductors
15 +
16 +
	//Minimal velocity
17 +
	extern double vMinimal_e(double q,double Ee,double mDM);
18 +
19 +
	//Recoil spectra
20 +
	extern double dRdEe(double Ee,const DM_Particle& DM,const std::string& target,double efficiency);
21 +
	extern Interpolation dRdEe(const DM_Particle& DM,const std::string& target,double efficiency,const std::vector<double>& attenuation,const std::vector<DataPoint> &speeddata,double vCutoff);
22 +
	//Total event rates
23 +
24 +
	extern double TotalRate(const DM_Particle& DM,int Qthreshold,double efficiency,const std::string& target);
25 +
	extern double TotalRate(const DM_Particle& DM,int Qthreshold,double efficiency,const std::string& target,const std::vector<double>& attenuation,const std::vector<DataPoint> &speeddata,double vCutoff);
26 +
27 +
//3. Atomic shell struct for liquid noble targets.
28 +
	struct Shell 
29 +
	{
30 +
		std::string name;
31 +
		double Ebinding;
32 +
		int neSecondary;
33 +
		std::vector< std::vector<double> > FormFactor;
34 +
		double logkMin,logkMax,logqMin,logqMax, dlogk,dlogq;
35 +
		int nk,nq;
36 +
			//Constructors
37 +
				Shell(std::string nm,std::string filepath,double Eb,int Ns,double lnkMin,double lnkMax,int Nk,double lnqMin,double lnqMax,int Nq);
38 +
	};
39 +
	extern std::vector<Shell> Shells;
40 +
41 +
//4. Liquid noble gas experiments
42 +
43 +
	//Recoil energy spectrum
44 +
	extern double dRdlog10Ee(double Ee,const DM_Particle& DM,const std::string& target,const Shell& shell);
45 +
46 +
	//Spectrum number of electrons
47 +
	extern double PDFne(const std::string& target,unsigned int ne,double Ee,const Shell& shell);
48 +
	extern double dRdne(unsigned int ne,const DM_Particle& DM,const std::string& target,const Shell& shell);
49 +
	extern std::vector<double> dRdne(const DM_Particle& DM,const std::string& target,const std::vector<Shell>& shells,const std::vector<double> &attenuation,const std::vector<DataPoint> &speeddata,double vCutoff);
50 +
	
51 +
	//Spectrum number of PE
52 +
	extern double dRdnPE(unsigned int nPE,const DM_Particle& DM,const std::string& target,double muPE,double sigPE,const std::vector<Shell>& shells);
53 +
	extern std::vector<double> dRdnPE(double muPE,double sigPE,const std::vector<double>& dRdn);
54 +
	extern double dNdnPE(unsigned int nPE,const DM_Particle& DM,const std::string& target,const std::string& detector,const std::vector<Shell>& shells);
55 +
	extern std::vector<double> dNdnPE(const std::string& detector,const std::vector<double>& dRdnpe);
56 +
57 +
#endif

@@ -0,0 +1,585 @@
Loading
1 +
#include "Input_Parameters.hpp"
2 +
3 +
#include <iostream>
4 +
#include <fstream>
5 +
#include <libconfig.h++>
6 +
#include <sys/types.h> // required for stat.h
7 +
#include <sys/stat.h>	//required to create a folder
8 +
9 +
#include "DD_Electron.hpp"
10 +
11 +
using namespace libconfig;
12 +
13 +
14 +
//1. Global variables and input parameters
15 +
	//Version
16 +
		std::string version = "1.1";
17 +
	//Input variables
18 +
		//Simulation ID
19 +
			std::string ID="default";
20 +
		//Statistical parameter
21 +
			unsigned int SampleSize=0;
22 +
			//Certainty Level
23 +
			double CL=-1.0;
24 +
			//Importance sampling:
25 +
			double IS_Angle=0.0;	//Scattering Angle
26 +
			double IS_MFP=0.0;	//MFP
27 +
			//Importance splitting
28 +
			bool GIS=false;
29 +
			double GIS_Splits=0.0;
30 +
			double GIS_Kappa=0.0;
31 +
32 +
			//Parameter Scan
33 +
			double mMin=0.0;
34 +
			double mMax=0.0;
35 +
			double dm=0.0;
36 +
			int Masses=0;
37 +
			double dSigma=0.0;
38 +
		//Layer structure
39 +
		bool Atmosphere = false;
40 +
		double Altitude = 0.0;
41 +
		int Atmo_Layers=0;
42 +
		std::vector<Layer> Layers;
43 +
		//Dark Matter
44 +
		DM_Particle DM;
45 +
		//Detector
46 +
		std::string Detector="default";
47 +
		//DM-Electron scattering with semiconductors
48 +
		std::string DMe_target="default";
49 +
		double DMe_exposure=0.0;
50 +
		int DMe_threshold=0.0;
51 +
		double DMe_efficiency=1.0;
52 +
		unsigned int DMe_events=0;
53 +
		//Detector Location
54 +
		double Detector_Depth=-1.0;
55 +
		int Detector_Index=-1;
56 +
57 +
//2. Read in the input parameters
58 +
	void Copy_Config_File(const char *inputfile)
59 +
	{
60 +
		std::ifstream    inFile;
61 +
		std::ofstream outFile;
62 +
		inFile.open(inputfile);
63 +
    	outFile.open("../results/"+ID+"/"+ID+".cfg");
64 +
		outFile << inFile.rdbuf();
65 +
		inFile.close();
66 +
		outFile.close();
67 +
	}
68 +
69 +
	void Read_Config_File(const char *inputfile,int rank)
70 +
	{
71 +
		std::string line="----------";
72 +
		Config cfg;
73 +
		// Read the file. If there is an error, report it and exit.
74 +
			try
75 +
			{		
76 +
				cfg.readFile(inputfile);
77 +
				if(rank==0) cout <<"Config file:\t\t"<<inputfile<<endl<<line<<endl;
78 +
			}
79 +
			catch(const FileIOException &fioex)
80 +
			{
81 +
				std::cerr << "I/O error while reading configuration file." << std::endl;
82 +
				exit(EXIT_FAILURE);
83 +
			}
84 +
			catch(const ParseException &pex)
85 +
			{
86 +
				std::cerr << "Configurate file parse error at " << pex.getFile() << ":" << pex.getLine() << " - " << pex.getError() << std::endl;
87 +
				exit(EXIT_FAILURE);
88 +
			}
89 +
		//Simulation ID
90 +
			try
91 +
			{
92 +
				ID = cfg.lookup("simID").c_str();
93 +
				if(rank==0) cout <<"\tSimulation ID:\t\t" <<ID <<endl;
94 +
			}
95 +
			catch(const SettingNotFoundException &nfex)
96 +
			{
97 +
				cerr << "No 'simID' setting in configuration file." << endl;
98 +
				exit(EXIT_FAILURE);
99 +
			}
100 +
		//Create folder for results
101 +
			if(rank==0) 
102 +
			{
103 +
				cout <<endl<<"\tCreating folder /results/"+ID <<"." <<endl;
104 +
				std::string sPath = "../results/"+ID;
105 +
				mode_t nMode = 0733; // UNIX style permissions
106 +
				int nError = 0;
107 +
				#if defined(_WIN32)
108 +
				  nError = _mkdir(sPath.c_str()); // can be used on Windows
109 +
				#else 
110 +
				  nError = mkdir(sPath.c_str(),nMode); // can be used on non-Windows
111 +
				#endif
112 +
				if (nError != 0) 
113 +
				{
114 +
					cout <<"\tThe folder already exists, data will be overwritten."<<endl;
115 +
				}
116 +
			}
117 +
		//Copy cfg file into this folder
118 +
			Copy_Config_File(inputfile);
119 +
120 +
		//Initial SampleSize
121 +
			if(rank == 0) cout <<endl<<"\tStatistics:"<<endl;
122 +
123 +
			try
124 +
			{
125 +
				SampleSize = cfg.lookup("samplesize");
126 +
				if(rank==0) cout <<"\t\tMinimal sample size:\t" <<SampleSize <<endl;
127 +
			}
128 +
			catch(const SettingNotFoundException &nfex)
129 +
			{
130 +
				cerr << "No 'samplesize' setting in configuration file." << endl;
131 +
				exit(EXIT_FAILURE);
132 +
			}
133 +
		//Certainty Level
134 +
			try
135 +
			{
136 +
				CL = cfg.lookup("cl");
137 +
				if(rank==0) cout <<"\t\tCertainty level:\t" <<CL <<endl;
138 +
			}
139 +
			catch(const SettingNotFoundException &nfex)
140 +
			{
141 +
				cerr << "No 'cl' setting in configuration file." << endl;
142 +
				exit(EXIT_FAILURE);
143 +
			}
144 +
		//Halo parameter
145 +
			if(rank == 0) cout <<endl<<"\tDM halo parameter:"<<endl;
146 +
			//DM energy density
147 +
			try
148 +
			{
149 +
				rhoDM = cfg.lookup("rhoDM");
150 +
				rhoDM *=GeV/cm/cm/cm;
151 +
				if(rank==0) cout <<"\t\trhoDM\t\t" <<InUnits(rhoDM,GeV/cm/cm/cm)<<" GeV/cm^3" <<endl;
152 +
153 +
			}
154 +
			catch(const SettingNotFoundException &nfex)
155 +
			{
156 +
				cerr << "No 'rhoDM' setting in configuration file." << endl;
157 +
				exit(EXIT_FAILURE);
158 +
			}
159 +
			//Velocity dispersion
160 +
			try
161 +
			{
162 +
				v0 = cfg.lookup("v0");
163 +
				v0 *=km/sec;
164 +
				if(rank==0) cout <<"\t\tv0:\t\t" <<InUnits(v0,km/sec)<<" km/sec" <<endl;
165 +
166 +
			}
167 +
			catch(const SettingNotFoundException &nfex)
168 +
			{
169 +
				cerr << "No 'v0' setting in configuration file." << endl;
170 +
				exit(EXIT_FAILURE);
171 +
			}
172 +
			//Earth velocity
173 +
			try
174 +
			{
175 +
				vEarth = cfg.lookup("vEarth");
176 +
				vEarth *=km/sec;
177 +
				if(rank==0) cout <<"\t\tvE:\t\t" <<InUnits(vEarth,km/sec)<<" km/sec" <<endl;
178 +
179 +
			}
180 +
			catch(const SettingNotFoundException &nfex)
181 +
			{
182 +
				cerr << "No 'vEarth' setting in configuration file." << endl;
183 +
				exit(EXIT_FAILURE);
184 +
			}
185 +
			//Galactic escape velocity
186 +
			try
187 +
			{
188 +
				vesc = cfg.lookup("vEscape");
189 +
				vesc *=km/sec;
190 +
				if(rank==0) cout <<"\t\tvEscape:\t" <<InUnits(vesc,km/sec)<<" km/sec" <<endl;
191 +
192 +
			}
193 +
			catch(const SettingNotFoundException &nfex)
194 +
			{
195 +
				cerr << "No 'vEscape' setting in configuration file." << endl;
196 +
				exit(EXIT_FAILURE);
197 +
			}
198 +
199 +
		//Variation reduction
200 +
			if(rank == 0) cout <<endl<<"\tVariation reduction techniques:"<<endl;
201 +
		//Importance sampling parameter for the scattering angle:
202 +
			try
203 +
			{
204 +
				IS_Angle = cfg.lookup("is_angle");
205 +
			}
206 +
			catch(const SettingNotFoundException &nfex)
207 +
			{
208 +
				cerr << "No 'is_angle' setting in configuration file." << endl;
209 +
				exit(EXIT_FAILURE);
210 +
			}
211 +
		//Importance sampling parameter for the mean free path:
212 +
			try
213 +
			{
214 +
				IS_MFP = cfg.lookup("is_mfp");
215 +
			}
216 +
			catch(const SettingNotFoundException &nfex)
217 +
			{
218 +
				cerr << "No 'is_mfp' setting in configuration file." << endl;
219 +
				exit(EXIT_FAILURE);
220 +
			}
221 +
			if(rank==0 && (IS_Angle>0.0||IS_MFP>0.0))
222 +
			{
223 +
				cout <<"\t\tImportance Sampling\t[x]"<<endl;
224 +
				cout <<"\t\t\tIS (angle):\t" <<IS_Angle <<endl;
225 +
				cout <<"\t\t\tIS (MFP):\t" <<IS_MFP <<endl;
226 +
			}
227 +
			else if(rank ==0) cout <<"\t\tImportance Sampling\t[ ]"<<endl;
228 +
		//Importance splitting layer:
229 +
			try
230 +
			{
231 +
				GIS = cfg.lookup("splitting");
232 +
			}
233 +
			catch(const SettingNotFoundException &nfex)
234 +
			{
235 +
				cerr << "No 'splitting' setting in configuration file." << endl;
236 +
				exit(EXIT_FAILURE);
237 +
			}
238 +
		//Importance splits at boundary:
239 +
			try
240 +
			{
241 +
				GIS_Splits = cfg.lookup("splits");
242 +
				if(GIS_Splits<=1&&GIS)
243 +
				{
244 +
					cerr << "Error: Option 'splits' in configuration file has to be larger than 1 for active GIS." << endl;
245 +
					exit(EXIT_FAILURE);
246 +
				}
247 +
			}
248 +
			catch(const SettingNotFoundException &nfex)
249 +
			{
250 +
				cerr << "No 'splits' setting in configuration file." << endl;
251 +
				exit(EXIT_FAILURE);
252 +
			}
253 +
		//Coefficients
254 +
			try
255 +
			{
256 +
				GIS_Kappa = cfg.lookup("kappa");
257 +
			}
258 +
			catch(const SettingNotFoundException &nfex)
259 +
			{
260 +
				cerr << "No 'kappa' setting in configuration file." << endl;
261 +
				exit(EXIT_FAILURE);
262 +
			}
263 +
			if(rank==0 &&GIS)
264 +
			{
265 +
				cout <<"\t\tImportance Splitting\t[x]"<<endl;
266 +
				cout <<"\t\t\tSplits:\t\t" <<GIS_Splits <<endl;
267 +
				cout <<"\t\t\tKappa:\t\t" <<GIS_Kappa <<endl<<endl;
268 +
			}
269 +
			else if(rank ==0) cout <<"\t\tImportance Splitting\t[ ]"<<endl<<endl;
270 +
271 +
		//Interaction parameter
272 +
			if(rank == 0) cout <<"\tDM interactions:"<<endl;
273 +
			//De-/Activate light DM option
274 +
			bool LDM;
275 +
			try
276 +
			{
277 +
				LDM = cfg.lookup("LDM");
278 +
				if(rank==0) cout <<"\t\tLight DM:\t\t" <<((LDM)? "[x]" : "[ ]")<<endl;
279 +
			}
280 +
			catch(const SettingNotFoundException &nfex)
281 +
			{
282 +
				cerr << "No 'LDM' setting in configuration file." << endl;
283 +
				exit(EXIT_FAILURE);
284 +
			}
285 +
			//Form factor F_DM
286 +
			std::string FF_DM;
287 +
			try
288 +
			{
289 +
				FF_DM = cfg.lookup("DM_FormFactor").c_str();
290 +
				if(rank==0) cout <<"\t\tForm factor F_DM:\t" <<FF_DM <<endl;
291 +
			}
292 +
			catch(const SettingNotFoundException &nfex)
293 +
			{
294 +
				cerr << "No 'DM_FormFactor' setting in configuration file." << endl;
295 +
				exit(EXIT_FAILURE);
296 +
			}
297 +
			//Coupling to Z or A
298 +
			std::string ZorA;
299 +
			try
300 +
			{
301 +
				ZorA = cfg.lookup("ZorA").c_str();
302 +
				if(rank==0) cout <<"\t\tCoupling to:\t\t" <<ZorA <<endl;
303 +
			}
304 +
			catch(const SettingNotFoundException &nfex)
305 +
			{
306 +
				cerr << "No 'ZorA' setting in configuration file." << endl;
307 +
				exit(EXIT_FAILURE);
308 +
			}
309 +
			//Screening
310 +
			bool Screening;
311 +
			if(FF_DM=="Electric-Dipole"||FF_DM=="Long-Range") Screening = true;
312 +
			else
313 +
			{
314 +
				try
315 +
				{
316 +
					Screening = cfg.lookup("Screening");
317 +
				}
318 +
				catch(const SettingNotFoundException &nfex)
319 +
				{
320 +
					cerr << "No 'Screening' setting in configuration file." << endl;
321 +
					exit(EXIT_FAILURE);
322 +
				}
323 +
			}
324 +
			if(rank==0) cout <<"\t\tCharge screening:\t" <<((Screening)? "[x]" : "[ ]") <<endl;
325 +
			//Mediator mass
326 +
			double mMediator;
327 +
			if(FF_DM=="General")
328 +
			{
329 +
				try
330 +
				{
331 +
					mMediator = cfg.lookup("mMediator");
332 +
					mMediator*=MeV;
333 +
					if(rank==0) cout <<"\tmMediator:" <<mMediator <<endl;
334 +
				}
335 +
				catch(const SettingNotFoundException &nfex)
336 +
				{
337 +
					cerr << "No 'mMediator' setting in configuration file." << endl;
338 +
					exit(EXIT_FAILURE);
339 +
				}
340 +
			}
341 +
			else
342 +
			{
343 +
				mMediator=0.0;
344 +
			}
345 +
			DM = DM_Particle(0.0,0.0,0.0,LDM,FF_DM,ZorA,Screening,mMediator);
346 +
			if(rank==0) cout <<endl;
347 +
		//Experiment
348 +
			if(rank==0) cout <<"\tExperiment:" <<endl;
349 +
			try
350 +
			{
351 +
				Detector = cfg.lookup("experiment").c_str();
352 +
				if(rank==0) cout <<"\t\tDetector:\t" <<Detector <<endl;
353 +
			}
354 +
			catch(const SettingNotFoundException &nfex)
355 +
			{
356 +
				cerr << "No 'experiment' setting in configuration file." << endl;
357 +
				exit(EXIT_FAILURE);
358 +
			}
359 +
		//DM-Electron
360 +
			if(Detector == "Semiconductor")
361 +
			{
362 +
				
363 +
				try
364 +
				{
365 +
					DMe_target = cfg.lookup("target").c_str();
366 +
					if(DMe_target!="Si"&&DMe_target!="Ge")
367 +
					{
368 +
						std::cerr <<"Error in cfg file. Semiconductor target " <<DMe_target<<" not recognized."<<endl;
369 +
						std::exit(EXIT_FAILURE);
370 +
					}
371 +
					//Import form factor
372 +
					else
373 +
					{
374 +
						Import_FormFactor(DMe_target);
375 +
					}
376 +
					DMe_threshold = cfg.lookup("threshold");
377 +
					DMe_exposure = cfg.lookup("exposure");
378 +
					DMe_efficiency = cfg.lookup("efficiency");
379 +
					DMe_events = cfg.lookup("events");
380 +
					if(rank==0) cout <<"\t\tTarget:\t\t" <<DMe_target <<endl
381 +
									<<"\t\tThreshold:\t" <<DMe_threshold <<endl
382 +
								<<"\t\tExposure[g yr]:\t" <<DMe_exposure <<endl
383 +
								<<"\t\tEfficiency:\t" <<DMe_efficiency <<endl
384 +
								<<"\t\tEvents:\t\t" <<DMe_events <<endl<<endl;
385 +
					DMe_exposure*=gram*year;
386 +
				}
387 +
				catch(const SettingNotFoundException &nfex)
388 +
				{
389 +
					cerr << "No 'Semiconductor' setting in configuration file." << endl;
390 +
					exit(EXIT_FAILURE);
391 +
				}	
392 +
			}
393 +
			else if(Detector=="SENSEI"||Detector=="SENSEI-surface")
394 +
			{
395 +
				
396 +
					DMe_target="Si";
397 +
					Import_FormFactor(DMe_target);//Import form factor
398 +
					DMe_threshold = 1;
399 +
					DMe_exposure = (Detector=="SENSEI-surface")?	0.07*gram*456*minute : 0.246*gram*day;
400 +
					DMe_efficiency = 1.0;
401 +
			}
402 +
			else if(Detector=="SuperCDMS"||Detector=="DAMIC-M")
403 +
			{
404 +
				
405 +
					DMe_target="Si";
406 +
					Import_FormFactor(DMe_target);//Import form factor
407 +
					DMe_threshold = 1;
408 +
					DMe_exposure = (Detector=="SuperCDMS")? 0.487*gram*day : 1.0*kg*year;
409 +
					DMe_efficiency = (Detector=="SuperCDMS")? 0.9545 : 1.0;
410 +
			}
411 +
			else if(Detector == "XENON10e"||Detector=="XENON100e")
412 +
			{
413 +
				//Import form factor
414 +
				DMe_target="Xe";
415 +
				// Import_FFion();	
416 +
				Import_FormFactor(DMe_target);
417 +
			}
418 +
			else if(Detector == "DarkSide-50")
419 +
			{
420 +
				//Import form factor
421 +
				DMe_target="Ar";
422 +
				// Import_FFion();	
423 +
				Import_FormFactor(DMe_target);
424 +
			}
425 +
		//Parameter scan
426 +
			if(rank==0) cout <<"\tParameter scan:" <<endl;
427 +
			//mMin
428 +
				try
429 +
				{
430 +
					mMin = cfg.lookup("mMin");
431 +
					if(rank==0) cout <<"\t\tmMin [GeV]:\t" <<mMin<<endl;
432 +
				}
433 +
				catch(const SettingNotFoundException &nfex)
434 +
				{
435 +
					cerr << "No 'mMin' setting in configuration file." << endl;
436 +
					exit(EXIT_FAILURE);
437 +
				}
438 +
			//mMax
439 +
				try
440 +
				{
441 +
					mMax = cfg.lookup("mMax");
442 +
					if(rank==0) cout <<"\t\tmMax [GeV]:\t" <<mMax <<endl;
443 +
				}
444 +
				catch(const SettingNotFoundException &nfex)
445 +
				{
446 +
					cerr << "No 'mMax' setting in configuration file." << endl;
447 +
					exit(EXIT_FAILURE);
448 +
				}
449 +
			//masses
450 +
				try
451 +
				{
452 +
					Masses = cfg.lookup("masses");
453 +
					if(rank==0) cout <<"\t\tMass steps:\t" <<Masses <<endl;
454 +
				}
455 +
				catch(const SettingNotFoundException &nfex)
456 +
				{
457 +
					cerr << "No 'masses' setting in configuration file." << endl;
458 +
					exit(EXIT_FAILURE);
459 +
				}
460 +
				//masses
461 +
				try
462 +
				{
463 +
					dSigma = cfg.lookup("dSigma");
464 +
					if(rank==0) cout <<"\t\tCS stepsize:\t" <<dSigma <<endl;
465 +
				}
466 +
				catch(const SettingNotFoundException &nfex)
467 +
				{
468 +
					cerr << "No 'dSigma' setting in configuration file." << endl;
469 +
					exit(EXIT_FAILURE);
470 +
				}
471 +
			//Number of steps
472 +
				if(mMin==mMax) Masses=1;
473 +
				if(Masses==1) 	dm=0;
474 +
				else 	dm = (log10(mMax)-log10(mMin))/(Masses-1.0);
475 +
		//Layer Structure
476 +
				if(rank==0) cout <<endl<<"\tShielding Layers:"<<endl <<endl;
477 +
			//Atmosphere
478 +
				try
479 +
				{
480 +
					Atmosphere = cfg.lookup("atmosphere");
481 +
				}
482 +
				catch(const SettingNotFoundException &nfex)
483 +
				{
484 +
					cerr << "No 'atmosphere' setting in configuration file." << endl;
485 +
					exit(EXIT_FAILURE);
486 +
				}
487 +
				if(Atmosphere)
488 +
				{
489 +
					if(rank==0) cout <<"\tAtmosphere\t[x]"<<endl;
490 +
					try
491 +
					{
492 +
						Atmo_Layers = cfg.lookup("atmo_layers");
493 +
						if(rank==0) cout <<"\t\tAtm. layers:\t" <<Atmo_Layers <<endl;
494 +
					}
495 +
					catch(const SettingNotFoundException &nfex)
496 +
					{
497 +
						cerr << "No 'atmo_layers' setting in configuration file." << endl;
498 +
						exit(EXIT_FAILURE);
499 +
					}
500 +
					try
501 +
					{
502 +
						Altitude = cfg.lookup("altitude");
503 +
						if(rank==0) cout <<"\t\tAltitude[m]:\t" <<Altitude <<endl;
504 +
						if(Altitude>=86000.0||Altitude<0.0)
505 +
						{
506 +
							cerr << "Error: Altitude must be chosen between 0.0 and 86000.0." << endl;
507 +
							exit(EXIT_FAILURE);
508 +
						}
509 +
						Altitude*=meter;
510 +
					}
511 +
					catch(const SettingNotFoundException &nfex)
512 +
					{
513 +
						cerr << "No 'altitude' setting in configuration file." << endl;
514 +
						exit(EXIT_FAILURE);
515 +
					}
516 +
					Layers = Atmospheric_Layers(Atmo_Layers,Altitude);
517 +
				
518 +
				}
519 +
				else
520 +
				{
521 +
					if(rank==0) cout <<"\tAtmosphere\t[ ]"<<endl;
522 +
				}
523 +
			//User defined layers
524 +
			double layer_depth=(Layers.size()==0)? 0.0 : Layers.back().depth+Layers.back().thickness;
525 +
			const Setting& root = cfg.getRoot();
526 +
			try
527 +
			{
528 +
				const Setting &layers = root["layers"];
529 +
				int count=layers.getLength();
530 +
				if(rank==0) cout <<endl<<"\tLayers:\t\t\t" <<count<<endl;
531 +
			 	if(rank==0) cout<<"\t"<<line<<endl
532 +
			 				<<"\tLayer 0:\t\tSpace"<<endl;
533 +
				for(int i=0;i<count;i++)
534 +
				{
535 +
					const Setting &layer =layers[i];
536 +
					std::string layer_name=layer.lookup("name").c_str();
537 +
					double layer_density=layer.lookup("density");
538 +
					double layer_thickness=layer.lookup("thickness");
539 +
					if(layer_thickness==0.0||layer_density==0.0) continue;
540 +
					//Units
541 +
						layer_density*=gram/cm/cm/cm;
542 +
						layer_thickness*=meter;
543 +
					//Composition
544 +
						std::vector<std::vector<double>> layer_composition;
545 +
						int element_count=layer.lookup("composition").getLength();
546 +
						for(int j=0;j<element_count;j++)
547 +
						{
548 +
							std::vector<double> element;
549 +
							for(int k=0;k<3;k++) element.push_back(layer.lookup("composition")[j][k]);
550 +
							layer_composition.push_back(element);	
551 +
						}
552 +
					//Construct Layer
553 +
						unsigned int layer_index=Layers.size()+1;
554 +
						Layer next_layer(layer_name,layer_index,layer_density,layer_thickness,layer_depth,layer_composition);
555 +
						Layers.push_back(next_layer);
556 +
						
557 +
					//Increase depth by the layer's thickness
558 +
						layer_depth+=layer_thickness;
559 +
				}
560 +
				for(unsigned int i =0;i<Layers.size();i++)
561 +
				{
562 +
					if(rank==0) 
563 +
					{
564 +
						cout<<"\t"<<line<<endl;
565 +
						Layers[i].Print_Summary();
566 +
					}
567 +
				}
568 +
				if(rank==0) cout<<"\t"<<line<<endl;
569 +
			}
570 +
			catch(const SettingNotFoundException &nfex)
571 +
			{
572 +
			    // Ignore.
573 +
			}
574 +
		//Derived Quantities
575 +
			//Halo
576 +
			Nesc=M_PI*v0*v0*(sqrt(M_PI)*v0*erf(vesc/v0)-2*vesc*exp(-vesc*vesc/v0/v0));
577 +
			//Detector
578 +
			Detector_Depth = layer_depth;
579 +
			Detector_Index = Layers.size()+1;
580 +
			if(rank==0) cout 	<<"\tLayer " <<Detector_Index<<":\t\tDetector" <<endl
581 +
								<<"\tType:\t\t\t" <<Detector <<endl
582 +
			 					<<"\tDepth:\t\t\t" <<InUnits(Detector_Depth,meter)<<"m" <<endl
583 +
			 					<<"\t"<<line<<endl
584 +
			 					<<line<<endl;
585 +
	}

@@ -0,0 +1,1244 @@
Loading
1 +
#include "DD_General.hpp"
2 +
3 +
#include <iostream>
4 +
#include <fstream>
5 +
#include <cmath>
6 +
#include <cstdlib>
7 +
#include <functional>
8 +
#include <algorithm>
9 +
#include <chrono>
10 +
#include <numeric> //for std::accumulate
11 +
12 +
#include "Physics_Functions.hpp"
13 +
#include "Input_Parameters.hpp"
14 +
#include "DD_Nucleus.hpp"
15 +
#include "DD_Electron.hpp"
16 +
17 +
using namespace std::placeholders;
18 +
19 +
//1. Attenuation
20 +
	std::vector<double> Compute_Attenuation(unsigned long long int particles, const std::vector<DataPoint> &data,double vCutoff)
21 +
	{
22 +
		//1. We only simulate particles with v>vMin, hence we have to find out how many particles we neglected.
23 +
			//Define integrand and epsilon
24 +
				auto integrand = std::bind(SpeedDistribution,_1,vEarth);
25 +
				double epsilon = 0.00001;
26 +
			//Integrate
27 +
				double Simulated_Fraction = Integrate(integrand,vCutoff,vEarth+vesc,epsilon);
28 +
			//True number of particles
29 +
				double Ntrue = 1.0*particles/Simulated_Fraction;
30 +
		//2. Attenuation is the sum of the weights
31 +
			std::vector<double> att;
32 +
			double Weight_Sum = 0.0;
33 +
			double Weight2_Sum = 0.0;
34 +
			for(unsigned int i=0;i<data.size();i++)
35 +
			{
36 +
				Weight_Sum+=data[i].weight;
37 +
				Weight2_Sum+=data[i].weight*data[i].weight;
38 +
			}
39 +
			//Compute attenuation
40 +
				att.push_back(Weight_Sum/Ntrue);
41 +
				att.push_back(sqrt(Weight2_Sum)/Ntrue);
42 +
		return att;
43 +
	}
44 +
45 +
//2. Eta-Function: Integral of f(v)/v with lower integration limit v_Min
46 +
	//Analytic
47 +
	double EtaFunction(double vMin,double vE)
48 +
	{
49 +
		double xMin = vMin/v0;
50 +
		double xEsc = vesc/v0;
51 +
		double xE=vE/v0;
52 +
		if(xMin>xE+xEsc) return 0.0;
53 +
		else if(fabs(xMin-xE-xEsc)<1e-8) return 0.0;
54 +
		else if(xMin>abs(xE-xEsc)) return pow(M_PI,1.5)*v0*v0/2.0/Nesc/xE*(erf(xEsc)-erf(xMin-xE)-2.0/sqrt(M_PI)*(xE+xEsc-xMin)*exp(-xEsc*xEsc));
55 +
		else if (xEsc>xE) return pow(M_PI,1.5)*v0*v0/2.0/Nesc/xE*(erf(xMin+xE)-erf(xMin-xE)-4.0/sqrt(M_PI)*xE*exp(-xEsc*xEsc));
56 +
		else return 1.0/v0/xE;
57 +
	}
58 +
	
59 +
	//MC based
60 +
	Interpolation EtaFunction(const std::vector<double>& attenuation,const std::vector<DataPoint> &speeddata,double vCutoff,double vE)
61 +
	{
62 +
		Interpolation kde = Perform_KDE(speeddata,vCutoff,(vesc+vEarth));
63 +
		kde.Multiply(attenuation[0]);
64 +
		int points =200;
65 +
		double dv=(vesc+vE-vCutoff)/(points-1.0);
66 +
		std::vector<std::vector<double>> interpol_list;
67 +
		for(int i=0;i<199;i++)
68 +
		{
69 +
			double vMin = vCutoff+i*dv;
70 +
			std::function<double(double)> integrand = [&kde](double v)
71 +
			{
72 +
				return kde(v)/v;
73 +
			};
74 +
			double epsilon=1e-2*(vesc+vE-vMin)*(integrand(vMin)+integrand((vMin+vEarth+vesc)/2.0))/2.0;
75 +
			double value=Integrate(integrand,vMin,vesc+vE,epsilon);
76 +
			interpol_list.push_back(std::vector<double> {vMin,value});
77 +
		}
78 +
		interpol_list.push_back(std::vector<double> {(vesc+vE),0.0});
79 +
		interpol_list.push_back(std::vector<double> {1.0,0.0});
80 +
		interpol_list.push_back(std::vector<double> {10.0,0.0});
81 +
82 +
		Interpolation eta(interpol_list);
83 +
		return eta;
84 +
	}
85 +
86 +
//3. Compute general recoil spectrum
87 +
	Interpolation Compute_Spectrum(const DM_Particle &DM,const std::string& experiment)
88 +
	{
89 +
		if(experiment == "DAMIC")
90 +
		{
91 +
			//1. Detector input
92 +
				// double efficiency = 0.75;
93 +
				double efficiency = 1.0;
94 +
				int Z = 14;
95 +
				double A=28.0;
96 +
				double X=1.0;
97 +
				double Ethr = 0.55*keV;
98 +
				double Emax = 7.0*keV;
99 +
			//2. Tabulate and interpolate
100 +
				//2.1 Find maximum ER, i.e. the domain of the spectrum
101 +
					double ERmax = ERMax((vesc+vEarth),DM.mass,A);
102 +
					ERmax = std::min(ERmax,Emax);
103 +
				//2.2 Tabulate the spectrum
104 +
					int points = 250;
105 +
					double dER = (ERmax-Ethr)/(points-1.0);
106 +
					std::vector<std::vector<double>> interpol_list;
107 +
					for(int i = 0;i<points;i++)
108 +
					{
109 +
						double ER = Ethr + i*dER;
110 +
						double dR = efficiency*dRdER(ER,DM,X,Z,A);
111 +
						interpol_list.push_back(std::vector<double> {ER,dR});
112 +
					}
113 +
				//2.3 if the maximum value is not Emax we append one more point, to avoid problems for the integration later on.
114 +
					if(ERmax<Emax)interpol_list.push_back(std::vector<double> {Emax,0.0});
115 +
				//2.4 Interpolate and return
116 +
					Interpolation spectrum(interpol_list);
117 +
					return spectrum;
118 +
				
119 +
		}
120 +
		else if(experiment=="XENON1T")
121 +
		{
122 +
			//1. Detector input
123 +
				double efficiency = 0.82;
124 +
				int Z = 54;
125 +
				double A=131.0;
126 +
				double X=1.0;
127 +
				double Ethr = 5*keV;
128 +
				double Emax = 40.0*keV;
129 +
			//2. Tabulate and interpolate
130 +
				//2.1 Find maximum ER, i.e. the domain of the spectrum
131 +
					double ERmax = ERMax((vesc+vEarth),DM.mass,A);
132 +
					ERmax = std::min(ERmax,Emax);
133 +
				//2.2 Tabulate the spectrum
134 +
					int points = 250;
135 +
					double dER = (ERmax-Ethr)/(points-1.0);
136 +
					std::vector<std::vector<double>> interpol_list;
137 +
					for(int i = 0;i<points;i++)
138 +
					{
139 +
						double ER = Ethr + i*dER;
140 +
						double dR = efficiency*dRdER(ER,DM,X,Z,A);
141 +
						interpol_list.push_back(std::vector<double> {ER,dR});
142 +
					}
143 +
				//2.3 if the maximum value is not Emax we append one more point, to avoid problems for the integration later on.
144 +
					if(ERmax<Emax)interpol_list.push_back(std::vector<double> {Emax,0.0});
145 +
				//2.4 Interpolate and return
146 +
					Interpolation spectrum(interpol_list);
147 +
					return spectrum;
148 +
		}
149 +
		else if(experiment=="CRESST-II")
150 +
		{
151 +
			//Detector Input
152 +
				int Z[3]	= {8,20,74};
153 +
				double A[3] = {16.,40.,184.};
154 +
				double X[3] = {0.22,0.14,0.64};				
155 +
				double Ethr = 307*eV;
156 +
				double resolution = Ethr/5.0;
157 +
				double Emax = 40*keV;
158 +
			//1. Import efficiencies
159 +
				Interpolation eff_O("../detectors/CRESST-II/Lise_eff_AR_O.dat",keV);
160 +
				Interpolation eff_Ca("../detectors/CRESST-II/Lise_eff_AR_Ca.dat",keV);
161 +
				Interpolation eff_W("../detectors/CRESST-II/Lise_eff_AR_W.dat",keV);
162 +
163 +
			//2. Tabulate dR/dE
164 +
				//2.1 Find the maximum E for which we have contributions.
165 +
					std::vector<double> ERmax_aux = {ERMax((vesc+vEarth),DM.mass,A[0]),ERMax((vesc+vEarth),DM.mass,A[1]),ERMax((vesc+vEarth),DM.mass,A[2])};
166 +
					double ERmax = *std::max_element(ERmax_aux.begin(),ERmax_aux.end());
167 +
					if((ERmax+6.0*resolution)<Emax) Emax =  (ERmax+6.0*resolution);
168 +
				//2.2 Tabulate dRdE
169 +
					std::vector<std::vector<double>> table;
170 +
					int steps = 200;
171 +
					double dE = (Emax - Ethr)/(steps-1.0);
172 +
173 +
					for(int i =0; i<steps; i++)
174 +
					{
175 +
						double E = Ethr+ i*dE;
176 +
						//2.2.1 Find minimum and maximum ER contributing to dR/dE(E):
177 +
							std::vector<double> aux = {E-6.0*resolution,Ethr-3.0*resolution};
178 +
							double eMin = *std::max_element(aux.begin(),aux.end());
179 +
							double eMax = E+6.0*resolution;
180 +
						//2.2.2 Define integrand
181 +
							std::function<double(double)> integrand = [E,X,Z,A,DM,resolution,&eff_O,&eff_Ca,&eff_W] (double ER)
182 +
							{
183 +
								double sum = eff_O(E)*dRdER(ER,DM,X[0],Z[0],A[0])+eff_Ca(E)*dRdER(ER,DM,X[1],Z[1],A[1])+eff_W(E)*dRdER(ER,DM,X[2],Z[2],A[2]);
184 +
								return PDF_Gauss(E,ER,resolution)*sum;
185 +
							};
186 +
187 +
						//2.2.3 Integrate and append to list
188 +
							// double epsilon = Find_Epsilon(integrand,ERmin,ERmax,1e-5);
189 +
							double epsilon = 1e-3*(eMax-eMin)*integrand(eMin);
190 +
							double  dRdE = Integrate(integrand, eMin,eMax,epsilon);
191 +
							table.push_back(std::vector<double> {E,dRdE});
192 +
193 +
194 +
					}
195 +
196 +
				//2.3 If we lowered Emax we have to put a final value, to not mess with the numerical integration later on.
197 +
				if(Emax<40*keV)
198 +
				{
199 +
					table.push_back(std::vector<double> {Emax+resolution,0.0});
200 +
					table.push_back(std::vector<double> {40*keV+6.0*resolution,0.0});
201 +
				} 
202 +
			//3. Interpolate and return spectrum
203 +
				Interpolation spectrum(table);
204 +
				return spectrum;
205 +
		}
206 +
		else if(experiment=="CRESST-surface")
207 +
		{
208 +
			//Detector Input
209 +
				double Z[2] = {8,13};
210 +
				double A[2] = {16.,27.};
211 +
				double X[2] = {0.47,0.53};				
212 +
				double Ethr = 19.7*eV;
213 +
				double resolution = 3.74*eV;
214 +
				double Emax = 600*eV;
215 +
			//1. Import efficiencies
216 +
				//not necessary here
217 +
			//2. Tabulate dR/dE
218 +
				//2.1 Find the maximum E for which we have contributions.
219 +
					std::vector<double> ERmax_aux = {ERMax((vesc+vEarth),DM.mass,A[0]),ERMax((vesc+vEarth),DM.mass,A[1])};
220 +
					double ERmax = *std::max_element(ERmax_aux.begin(),ERmax_aux.end());
221 +
					if((ERmax+6.0*resolution)<Emax) Emax =  (ERmax+6.0*resolution);
222 +
				//2.2 Tabulate dRdE
223 +
					std::vector<std::vector<double>> table;
224 +
					int steps = 200;
225 +
					double dE = (Emax - Ethr)/(steps-1.0);
226 +
					for(int i =0; i<steps; i++)
227 +
					{
228 +
						double E = Ethr+ i*dE;
229 +
						//2.2.1 Find minimum and maximum ER contributing to dR/dE(E):
230 +
							std::vector<double> aux = {E-6.0*resolution,Ethr-3.0*resolution};
231 +
							double eMin = *std::max_element(aux.begin(),aux.end());
232 +
							double eMax = E+6.0*resolution;
233 +
						//2.2.2 Define integrand
234 +
							std::function<double(double)> integrand = [E,X,Z,A,DM,resolution] (double ER)
235 +
							{
236 +
								double sum = dRdER(ER,DM,X[0],Z[0],A[0])+dRdER(ER,DM,X[1],Z[1],A[1]);
237 +
								return PDF_Gauss(E,ER,resolution)*sum;
238 +
							};
239 +
						//2.2.3 Integrate and append to list
240 +
							double epsilon = 1e-3*(eMax-eMin)*integrand(eMin);
241 +
							double  dRdE = Integrate(integrand, eMin,eMax,epsilon);
242 +
							table.push_back(std::vector<double> {E,dRdE});
243 +
244 +
					}
245 +
				//2.3 If we lowered Emax we have to put a final value, to not mess with the numerical integration later on.
246 +
				if(Emax<600*eV)
247 +
				{
248 +
					table.push_back(std::vector<double> {Emax+resolution,0.0});
249 +
					table.push_back(std::vector<double> {600*eV+6.0*resolution,0.0});
250 +
				} 
251 +
			//3. Interpolate and return spectrum
252 +
				Interpolation spectrum(table);
253 +
				return spectrum;
254 +
		}
255 +
		else if(experiment == "Semiconductor"||experiment=="SENSEI"||experiment == "SENSEI-surface"||experiment=="SuperCDMS"||experiment=="DAMIC-M")
256 +
		{
257 +
			//Tabulate and interpolate
258 +
				//1. Domain of the spectrum
259 +
					double Eemax = 50*eV;
260 +
					double Eemin = 0.1*eV;
261 +
				//2. Tabulate the spectrum
262 +
					int points = 250;
263 +
					double dEe = (Eemax-Eemin)/(points-1.0);
264 +
					std::vector<std::vector<double>> interpol_list;
265 +
					for(int i = 0;i<points;i++)
266 +
					{
267 +
						double Ee = Eemin + i*dEe;
268 +
						double dR = dRdEe(Ee,DM,DMe_target,DMe_efficiency);
269 +
						interpol_list.push_back(std::vector<double> {Ee,dR});
270 +
					}
271 +
				//3. Interpolate and return
272 +
					Interpolation spectrum(interpol_list);
273 +
					return spectrum;
274 +
		}
275 +
		else if(experiment=="XENON10e"||experiment=="XENON100e")
276 +
		{
277 +
			//Tabulate and interpolate
278 +
				//1. Tabulate the PE spectrum
279 +
					int points = 250;
280 +
					std::vector<std::vector<double>> interpol_list;
281 +
					for(int nPE = 1;nPE<=points;nPE++)
282 +
					{
283 +
						double dN = dNdnPE(nPE,DM,DMe_target,experiment,Shells);
284 +
						interpol_list.push_back(std::vector<double> {1.0*nPE,dN});
285 +
					}
286 +
				//2. Interpolate and return
287 +
					Interpolation spectrum(interpol_list);
288 +
					return spectrum;
289 +
		}
290 +
		else if(experiment=="DarkSide-50")
291 +
		{
292 +
			//1. Tabulate the ne spectrum
293 +
				std::vector<std::vector<double>> interpol_list;
294 +
				for(int ne=1;ne<=10;ne++)
295 +
				{
296 +
					double dR=0.0;
297 +
					for(unsigned int i=0;i<Shells.size();i++)
298 +
					{						
299 +
						dR+= dRdne(ne,DM,DMe_target,Shells[i]);
300 +
					} 
301 +
					interpol_list.push_back(std::vector<double> {1.0*ne,dR});
302 +
				}
303 +
			//2. Interpolate and return
304 +
				Interpolation spectrum(interpol_list);
305 +
				return spectrum;
306 +
		}
307 +
		else
308 +
		{
309 +
			//Error
310 +
				cerr <<"Error in Compute_Spectrum(): Experiment " <<experiment <<" not recognized."<<endl;
311 +
				std::exit(EXIT_FAILURE);
312 +
		}
313 +
	}
314 +
	Interpolation Compute_Spectrum(const DM_Particle &DM,const std::string& experiment,const std::vector<DataPoint> &data,const std::vector<double>& attenuation)
315 +
	{
316 +
		if(experiment == "DAMIC")
317 +
		{
318 +
			//Detector input
319 +
				double efficiency = 1.0;
320 +
				int Z = 14;
321 +
				double A=28.0;
322 +
				double X=1.0;
323 +
				double Ethr = 0.55*keV;
324 +
				double Emax = 7.0*keV;
325 +
			//spectrum:
326 +
				Interpolation spectrum = dRdER(DM,Ethr,Emax,X,Z,A,attenuation,data);
327 +
				spectrum.Multiply(efficiency);
328 +
				return spectrum;
329 +
				
330 +
		}
331 +
		else if(experiment=="XENON1T")
332 +
		{
333 +
			//Detector input
334 +
				double efficiency = 0.82;
335 +
				int Z = 54;
336 +
				double A=131.0;
337 +
				double X=1.0;
338 +
				double Ethr = 5*keV; 
339 +
				double Emax = 40.0*keV;
340 +
			//spectrum:
341 +
				Interpolation spectrum = dRdER(DM,Ethr,Emax,X,Z,A,attenuation,data);
342 +
				spectrum.Multiply(efficiency);
343 +
				return spectrum;
344 +
		}
345 +
		else if(experiment=="CRESST-II")
346 +
		{
347 +
348 +
			//0. Detector Parameter
349 +
				int Z[3]	= {8,20,74};
350 +
				double A[3] = {16.,40.,184.};
351 +
				double X[3] = {0.22,0.14,0.64};				
352 +
				double Ethr = 307*eV;
353 +
				double resolution = Ethr/5.0;
354 +
				double Emax = 40*keV;
355 +
			//1. Import efficiencies
356 +
				Interpolation eff_O("../detectors/CRESST-II/Lise_eff_AR_O.dat",keV);
357 +
				Interpolation eff_Ca("../detectors/CRESST-II/Lise_eff_AR_Ca.dat",keV);
358 +
				Interpolation eff_W("../detectors/CRESST-II/Lise_eff_AR_W.dat",keV);
359 +
			//2. Tabulate and interpolate dR/dE
360 +
					
361 +
				//2.1 Interpolate the 3 recoil spectra dR/dER
362 +
					double ERmin = Ethr-3.0*resolution;
363 +
					double ERmax = Emax + 6.0*resolution;
364 +
					Interpolation dRdER_O = dRdER(DM, ERmin,ERmax, X[0],Z[0],A[0], attenuation,data);
365 +
					Interpolation dRdER_Ca = dRdER( DM,ERmin,ERmax, X[1],Z[1],A[1], attenuation,data);
366 +
					Interpolation dRdER_W = dRdER(DM, ERmin,ERmax, X[2],Z[2],A[2],attenuation,data);
367 +
368 +
				//2.2 Find the minimum and maximum ER for which we have contributions.
369 +
					std::vector<double> ERmax_aux = {ERMax((vesc+vEarth),DM.mass,A[0]),ERMax((vesc+vEarth),DM.mass,A[1]),ERMax((vesc+vEarth),DM.mass,A[2])};
370 +
					ERmax = *std::max_element(ERmax_aux.begin(),ERmax_aux.end());
371 +
					if((ERmax+6.0*resolution)<Emax) Emax =  (ERmax+6.0*resolution);
372 +
				//2.3 Tabulate
373 +
					std::vector<std::vector<double>> table;
374 +
					int steps = 200;
375 +
					double dE = (Emax - Ethr)/(steps-1.0);
376 +
					for(int i =0; i<steps; i++)
377 +
					{
378 +
						double E = Ethr+ i*dE;
379 +
						//2.3.1 Find minimum and maximum ER contributing to dR/dE(E):
380 +
							std::vector<double> aux = {E-6.0*resolution,Ethr-3.0*resolution};
381 +
							double eMin = *std::max_element(aux.begin(),aux.end());
382 +
							double eMax = E+6.0*resolution;
383 +
						//2.3.2 Define integrand
384 +
							std::function<double(double)> integrand = [E,resolution,&dRdER_O,&dRdER_Ca,&dRdER_W,&eff_O,&eff_Ca,&eff_W] (double ER)
385 +
							{
386 +
								double sum = eff_O(E)*dRdER_O(ER)+eff_Ca(E)*dRdER_Ca(ER)+eff_W(E)*dRdER_W(ER);
387 +
								return PDF_Gauss(E,ER,resolution)*sum;
388 +
							};
389 +
390 +
						//2.3.3 Integrate and append to list							
391 +
							double epsilon = 1e-3*(eMax-eMin)*(dRdER_O(eMin)+dRdER_Ca(eMin)+dRdER_W(eMin));
392 +
							double dRdE=0.0;
393 +
							if(epsilon>0) //epsilon can become negative for reasons of double precision
394 +
							{
395 +
								dRdE = Integrate(integrand, eMin,eMax,epsilon);
396 +
							}
397 +
							table.push_back(std::vector<double> {E,dRdE});
398 +
399 +
					}
400 +
				//2.4 If we lowered Emax we have to put a final value, to not mess with the numerical integration later on.
401 +
					
402 +
					if(Emax<40*keV)
403 +
					{
404 +
						table.push_back(std::vector<double> {Emax+resolution,0.0});
405 +
						table.push_back(std::vector<double> {40*keV+6.0*resolution,0.0});
406 +
					} 
407 +
			//3. Interpolate and return spectrum
408 +
				Interpolation spectrum(table);
409 +
				return spectrum;
410 +
		}
411 +
		else if(experiment=="CRESST-surface")
412 +
		{
413 +
			//0. Detector Parameter
414 +
				int Z[2]	= {8,13};
415 +
				double A[2] = {16.,27.};
416 +
				double X[2] = {0.47,0.53};				
417 +
				double Ethr = 19.7*eV;
418 +
				double resolution = 3.74*eV;
419 +
				double Emax = 600*eV;
420 +
			//1. Import efficiencies
421 +
				//not necessary here
422 +
			//2. Tabulate and interpolate dR/dE
423 +
					
424 +
				//2.1 Interpolate the 3 recoil spectra dR/dER
425 +
					double ERmin = Ethr-3.0*resolution;
426 +
					double ERmax = Emax + 6.0*resolution;
427 +
					Interpolation dRdER_O = 	dRdER(DM,ERmin,ERmax,X[0],Z[0],A[0],attenuation,data);
428 +
					Interpolation dRdER_Al = 	dRdER(DM,ERmin,ERmax,X[1],Z[1],A[1],attenuation,data);
429 +
430 +
				//2.2 Find the minimum and maximum ER for which we have contributions.
431 +
					std::vector<double> ERmax_aux = {ERMax((vesc+vEarth),DM.mass,A[0]),ERMax((vesc+vEarth),DM.mass,A[1])};
432 +
					ERmax = *std::max_element(ERmax_aux.begin(),ERmax_aux.end());
433 +
					if((ERmax+6.0*resolution)<Emax) Emax =  (ERmax+6.0*resolution);
434 +
				//2.3 Tabulate
435 +
					std::vector<std::vector<double>> table;
436 +
					int steps = 200;
437 +
					double dE = (Emax - Ethr)/(steps-1.0);
438 +
					for(int i =0; i<steps; i++)
439 +
					{
440 +
						double E = Ethr+ i*dE;
441 +
						//2.3.1 Find minimum and maximum ER contributing to dR/dE(E):
442 +
							std::vector<double> aux = {E-6.0*resolution,Ethr-3.0*resolution};
443 +
							double eMin = *std::max_element(aux.begin(),aux.end());
444 +
							double eMax = E+6.0*resolution;
445 +
						//2.3.2 Define integrand
446 +
							std::function<double(double)> integrand = [E,resolution,&dRdER_O,&dRdER_Al] (double ER)
447 +
							{
448 +
								double sum = dRdER_O(ER)+dRdER_Al(ER);
449 +
								return PDF_Gauss(E,ER,resolution)*sum;
450 +
							};
451 +
452 +
						//2.3.3 Integrate and append to list
453 +
							double epsilon = 1e-3*(eMax-eMin)*(dRdER_O(eMin)+dRdER_Al(eMin));
454 +
							double dRdE=0.0;
455 +
							if(epsilon>0) //epsilon can become negative for reasons of double precision
456 +
							{
457 +
								dRdE = Integrate(integrand, eMin,eMax,epsilon);
458 +
							}
459 +
							table.push_back(std::vector<double> {E,dRdE});
460 +
					}
461 +
				//2.4 If we lowered Emax we have to put a final value, to not mess with the numerical integration later on.
462 +
					
463 +
					if(Emax<600*eV)
464 +
					{
465 +
						table.push_back(std::vector<double> {Emax+resolution,0.0});
466 +
						table.push_back(std::vector<double> {600*eV+6.0*resolution,0.0});
467 +
					} 
468 +
			//3. Interpolate and return spectrum
469 +
				Interpolation spectrum(table);
470 +
				return spectrum;
471 +
		}
472 +
		else if(experiment=="Semiconductor"||experiment=="SENSEI"||experiment=="SENSEI-surface"||experiment=="XENON10e"||experiment=="XENON100e"||experiment=="DarkSide-50"||experiment=="SuperCDMS"||experiment=="DAMIC-M")
473 +
		{
474 +
			//not necessary
475 +
			return Interpolation();
476 +
		}
477 +
		else
478 +
		{
479 +
			cerr <<"Error in Compute_Spectrum(): Experiment " <<experiment <<" not recognized."<<endl;
480 +
			std::exit(EXIT_FAILURE);
481 +
		}		
482 +
	}
483 +
484 +
485 +
//4. Likelihood
486 +
	double Likelihood(const DM_Particle& DM,const std::string& experiment,const std::function<double(double)>& spectrum, double &N)
487 +
	{
488 +
		double llh=-1.0;
489 +
		//DAMIC
490 +
		if(experiment=="DAMIC")
491 +
		{
492 +
			double exposure = 0.107*kg*day;
493 +
			double Ethreshold = 0.55*keV;
494 +
			double Emax = 7*keV;
495 +
			int nEvents = 106;
496 +
			N = N_Signals(spectrum,Ethreshold,Emax,exposure);
497 +
			llh = CDF_Poisson(N,nEvents);
498 +
		}
499 +
		else if(experiment=="XENON1T")
500 +
		{
501 +
			double exposure = 34.2*day*1042*kg;
502 +
			double Ethreshold = 5*keV;
503 +
			double Emax = 40*keV;
504 +
			N = N_Signals(spectrum,Ethreshold,Emax,exposure);
505 +
			llh = CDF_Poisson(N,0);
506 +
		}
507 +
		else if(experiment=="CRESST-II")
508 +
		{
509 +
			//0. Detector parameters
510 +
				double Ethr = 307*eV;
511 +
				double Emax = 40*keV;
512 +
				double exposure = 52.15*kg*day;
513 +
514 +
			//1. Import and sort events
515 +
				std::vector<double> events = Read_List("../detectors/CRESST-II/Lise_AR.dat",keV);
516 +
				events.push_back(Ethr);
517 +
				events.push_back(Emax);
518 +
				std::sort(events.begin(),events.end());
519 +
			//2. Integrate spectrum inbetween the events
520 +
				std::vector<double> x;
521 +
				for(unsigned int i = 0;i < (events.size()-1); i++)
522 +
				{
523 +
					//2.1 The gap
524 +
					double E1 = events[i];
525 +
					double E2 = events[i+1];
526 +
					//2.2 Integrate the gap
527 +
						double epsilon = Find_Epsilon(spectrum,E1,E2,1e-2);
528 +
						double xGap = exposure*Integrate(spectrum, E1,E2,epsilon);
529 +
					//2.3 Append the x value to the list
530 +
					x.push_back(xGap);
531 +
				}
532 +
			//3. Find xMax
533 +
				double xMax = *std::max_element(x.begin(),x.end());
534 +
			//4. Return likelihood.
535 +
				N = std::accumulate(x.begin(),x.end(),0.0);
536 +
				llh = 1.0 - MaxGap_C0(xMax,N);
537 +
		}
538 +
		else if(experiment=="CRESST-surface")
539 +
		{
540 +
			//0. Detector parameters
541 +
				double Ethr = 19.7*eV;
542 +
				double Emax = 600*eV;
543 +
				double exposure = 0.046*gram*day;
544 +
545 +
			//1. Import and sort events
546 +
				std::vector<double> events = Read_List("../detectors/CRESST-surface/data.txt",keV);
547 +
				events.push_back(Ethr);
548 +
				events.push_back(Emax);
549 +
				std::sort(events.begin(),events.end());
550 +
			//2. Integrate spectrum inbetween the events
551 +
				std::vector<double> x;
552 +
				for(unsigned int i = 0;i < (events.size()-1); i++)
553 +
				{
554 +
					//2.1 The gap
555 +
					double E1 = events[i];
556 +
					double E2 = events[i+1];
557 +
					//2.2 Integrate the gap
558 +
						double epsilon = Find_Epsilon(spectrum,E1,E2,1e-2);
559 +
						double xGap = exposure*Integrate(spectrum, E1,E2,epsilon);
560 +
					//2.3 Append the x value to the list
561 +
					x.push_back(xGap);
562 +
				}
563 +
			//3. Find xMax
564 +
				double xMax = *std::max_element(x.begin(),x.end());
565 +
			//4. Return likelihood.
566 +
				N = std::accumulate(x.begin(),x.end(),0.0);
567 +
				llh = 1.0 - MaxGap_C0(xMax,N);
568 +
		}
569 +
		else if(experiment=="Semiconductor")
570 +
		{
571 +
			N = DMe_exposure*TotalRate(DM,DMe_threshold,DMe_efficiency,DMe_target);
572 +
			llh = CDF_Poisson(N,DMe_events);
573 +
574 +
		}
575 +
		else if(experiment=="SENSEI"||experiment=="SENSEI-surface"||experiment=="SuperCDMS"||experiment=="DAMIC-M")
576 +
		{
577 +
			//Efficiency and number of events
578 +
				std::vector<double> efficiency;
579 +
				std::vector<int> BinData;
580 +
				if(experiment=="SENSEI")
581 +
				{
582 +
					efficiency={1.0,0.62,0.48};
583 +
					BinData = {8516,87,0};
584 +
585 +
				}
586 +
				else if(experiment=="SENSEI-surface")
587 +
				{
588 +
					efficiency={0.668,0.41,0.32,0.27,0.24};
589 +
					BinData={140302,4676,131,1,0};
590 +
				}
591 +
				else if(experiment=="SuperCDMS")
592 +
				{
593 +
					//all points within 2sigma, taken from fig. 3 of [arXiv:1804.10697]
594 +
					efficiency={0.88,0.91,0.91,0.91,0.91,0.91};
595 +
					BinData={53000, 400, 74, 18, 7, 14};
596 +
				}
597 +
				else if(experiment=="DAMIC-M")
598 +
				{
599 +
					efficiency={1.0,1.0,1.0,1.0,1.0,1.0};
600 +
					BinData={100000,0,0,0,0,0};
601 +
				}
602 +
			//Compute rates
603 +
				std::vector<double> rates;
604 +
				for(unsigned int bin=1;bin<=BinData.size()+1;bin++) rates.push_back(TotalRate(DM,bin,DMe_efficiency,DMe_target));
605 +
			//Compute number of events and likelihood per bin
606 +
				std::vector<double> events;
607 +
				std::vector<double> llhs;
608 +
				for(unsigned int bin=(DMe_threshold-1);bin<BinData.size();bin++) 
609 +
				{
610 +
					events.push_back(efficiency[bin]*DMe_exposure*(rates[bin]-rates[bin+1]));
611 +
					llhs.push_back(CDF_Poisson(events[bin],BinData[bin]));
612 +
					N+=events[bin];
613 +
				}
614 +
			//Find minimum value
615 +
				llh = *std::min_element(llhs.begin(),llhs.end());
616 +
		}
617 +
		else if(experiment=="XENON10e")
618 +
		{
619 +
			//Data is binned in 7 bins [14,41),[41,68),[68,95),[95,122),[122,149),[149,176),[176,203)
620 +
				std::vector<int> data={126, 60, 12, 3, 2, 0, 2};
621 +
				std::vector<unsigned int> bins={14,41,68,95,122,149,176,203};
622 +
				std::vector<double> llhs;
623 +
			//Compute events in each bin.
624 +
				for(unsigned int bin=0;bin<data.size();bin++)
625 +
				{
626 +
					double Nbin=0;
627 +
					for(unsigned int PE= bins[bin];PE<bins[bin+1];PE++)
628 +
					{
629 +
						Nbin+=spectrum(PE);
630 +
					}
631 +
					N+=Nbin;
632 +
					double l=CDF_Poisson(Nbin,data[bin]);
633 +
					llhs.push_back(l);
634 +
635 +
				}
636 +
			llh = *std::min_element(llhs.begin(),llhs.end());
637 +
		}
638 +
		else if(experiment=="XENON100e")
639 +
		{
640 +
			//Data is binned in 3 bins [80,90),[90,110),[110,130)
641 +
				std::vector<int> data={794, 1218, 924, 776, 669, 630, 528, 488, 433, 387};
642 +
				std::vector<unsigned int> bins={80, 90, 110,130,150,170,190,210,230,250,270};
643 +
				std::vector<double> llhs;		
644 +
			//Compute events in each bin.
645 +
				//We only use the first 3 bins.
646 +
				for(unsigned int bin=0;bin<3;bin++)
647 +
				{
648 +
					double Nbin=0;
649 +
					for(unsigned int PE= bins[bin];PE<bins[bin+1];PE++)
650 +
					{
651 +
						Nbin+=dNdnPE(PE,DM,DMe_target,experiment,Shells);
652 +
					}
653 +
					N+=Nbin;
654 +
					double l=CDF_Poisson(Nbin,data[bin]);
655 +
					llhs.push_back(l);
656 +
				}
657 +
			llh = *std::min_element(llhs.begin(),llhs.end());
658 +
		}
659 +
		else if(experiment=="DarkSide-50")
660 +
		{
661 +
			std::vector<int>data={118643, 219893, 6131, 673, 252, 227, 198, 199, 189, 247, 230, 261,249, 329, 336, 349, 351, 352, 384, 411, 405, 461, 460, 436, 500, 546, 538, 536, 556, 583, 573, 630, 603, 635, 639, 682, 736, 755, 804, 811, 809, 882, 934, 935, 871, 965, 946, 1072, 997, 1060};
662 +
			std::vector<double> llhs;
663 +
			//Threshold and exposure
664 +
				double exposure = 6786.0*kg*day;
665 +
				int nThreshold=3;
666 +
			//Compute likelihood for each bin.
667 +
				for(unsigned int ne=nThreshold;ne<data.size();ne++)
668 +
				{
669 +
					double Ntot=0.0;
670 +
					//Add all shells
671 +
					for(unsigned int shell=0;shell<Shells.size();shell++)
672 +
					{					
673 +
						Ntot+=exposure/2.0 * (dRdne(ne,DM,DMe_target,Shells[shell])+dRdne(ne+1,DM,DMe_target,Shells[shell]));
674 +
					
675 +
					}
676 +
					N+=Ntot;
677 +
					double l=CDF_Poisson(Ntot,data[ne-1]);
678 +
					llhs.push_back(l);
679 +
				}
680 +
			llh = *std::min_element(llhs.begin(),llhs.end());
681 +
682 +
		}
683 +
		else
684 +
		{
685 +
			std::cerr <<"Error in Likelihood(): Experiment " <<experiment <<" not recognized."<<endl;
686 +
			std::exit(EXIT_FAILURE);
687 +
		}
688 +
		return llh;
689 +
	}
690 +
691 +
	std::vector<double> Likelihood(const DM_Particle& DM,const std::string& experiment,const std::function<double(double)>& spectrum,const std::vector<double>& attenuation,const std::vector<DataPoint> &data,double vCutoff,std::vector<double> &NMC)
692 +
	{
693 +
		std::vector<double> llh;
694 +
		//DAMIC
695 +
		if(experiment=="DAMIC")
696 +
		{
697 +
			double exposure = 0.107*kg*day;
698 +
			double Ethreshold = 0.55*keV;
699 +
			double Emax = 7*keV;
700 +
			double nEvents=106;
701 +
			double n = N_Signals(spectrum,Ethreshold,Emax,exposure);
702 +
			NMC.push_back(n);
703 +
			NMC.push_back(attenuation[1]/attenuation[0]*n);
704 +
			llh.push_back( CDF_Poisson(n,nEvents) );
705 +
			if(llh[0]==0.0) llh.push_back(0.0);
706 +
			else llh.push_back( NMC[1]*exp(-NMC[0]) * pow(NMC[0],nEvents)/Factorial(nEvents) );
707 +
		}
708 +
		else if(experiment=="XENON1T")
709 +
		{
710 +
			double exposure = 34.2*day*1042*kg;
711 +
			double Ethreshold = 5*keV;
712 +
			double Emax = 40*keV;
713 +
			double n = N_Signals(spectrum,Ethreshold,Emax,exposure);
714 +
			NMC.push_back(n);
715 +
			NMC.push_back(attenuation[1]/attenuation[0]*n);
716 +
			llh.push_back( CDF_Poisson(n,0) );
717 +
			llh.push_back( NMC[1]*exp(-NMC[0]) );
718 +
		}
719 +
		else if(experiment=="CRESST-II")
720 +
		{
721 +
			//0. Detector parameters
722 +
				double Ethr = 307*eV;
723 +
				double Emax = 40*keV;
724 +
				double exposure = 52.15*kg*day;
725 +
726 +
			//1. Import and sort events
727 +
				std::vector<double> events = Read_List("../detectors/CRESST-II/Lise_AR.dat",keV);
728 +
				events.push_back(Ethr);
729 +
				events.push_back(Emax);
730 +
				std::sort(events.begin(),events.end());
731 +
			//2. Integrate spectrum inbetween the events
732 +
				std::vector<double> x;
733 +
				for(unsigned int i = 0;i < (events.size()-1); i++)
734 +
				{
735 +
					//2.1 The gap
736 +
					double E1 = events[i];
737 +
					double E2 = events[i+1];
738 +
					//2.2 Integrate the gap
739 +
						double epsilon = Find_Epsilon(spectrum,E1,E2,1e-2);
740 +
						double xGap = exposure*Integrate(spectrum, E1,E2,epsilon);
741 +
					//2.3 Append the x value to the list
742 +
					x.push_back(xGap);
743 +
				}
744 +
			//3. Find xMax
745 +
				double xMax = *std::max_element(x.begin(),x.end());
746 +
			//4. Return likelihood.
747 +
				double N = std::accumulate(x.begin(),x.end(),0.0);
748 +
				double error = attenuation[1]/attenuation[0]*N;
749 +
				NMC.push_back( N );
750 +
				NMC.push_back( error );
751 +
				llh.push_back(1.0 - MaxGap_C0(xMax,N));
752 +
				llh.push_back(MaxGap_C0_error(xMax,N,error) );
753 +
		}
754 +
		else if(experiment=="CRESST-surface")
755 +
		{
756 +
			//0. Detector parameters
757 +
				double Ethr = 19.7*eV;
758 +
				double Emax = 600*eV;
759 +
				double exposure = 0.046*gram*day;
760 +
761 +
			//1. Import and sort events
762 +
				std::vector<double> events = Read_List("../detectors/CRESST-surface/data.txt",keV);
763 +
				events.push_back(Ethr);
764 +
				events.push_back(Emax);
765 +
				std::sort(events.begin(),events.end());
766 +
			//2. Integrate spectrum inbetween the events
767 +
				std::vector<double> x;
768 +
				for(unsigned int i = 0;i < (events.size()-1); i++)
769 +
				{
770 +
					//2.1 The gap
771 +
					double E1 = events[i];
772 +
					double E2 = events[i+1];
773 +
					//2.2 Integrate the gap
774 +
						double epsilon = Find_Epsilon(spectrum,E1,E2,1e-2);
775 +
						double xGap = exposure*Integrate(spectrum, E1,E2,epsilon);
776 +
					//2.3 Append the x value to the list
777 +
					x.push_back(xGap);
778 +
				}
779 +
			//3. Find xMax
780 +
				double xMax = *std::max_element(x.begin(),x.end());
781 +
			//4. Return likelihood.
782 +
				double N = std::accumulate(x.begin(),x.end(),0.0);
783 +
				double error = attenuation[1]/attenuation[0]*N;
784 +
				NMC.push_back( N );
785 +
				NMC.push_back( error );
786 +
				llh.push_back(1.0 - MaxGap_C0(xMax,N));
787 +
				llh.push_back(MaxGap_C0_error(xMax,N,error) );
788 +
		}
789 +
		else if(experiment=="Semiconductor")
790 +
		{
791 +
			double N = DMe_exposure*TotalRate(DM,DMe_threshold,DMe_efficiency,DMe_target,attenuation,data,vCutoff);
792 +
			double error = attenuation[1]/attenuation[0]*N;
793 +
			NMC.push_back( N );
794 +
			NMC.push_back( error );
795 +
			llh.push_back(CDF_Poisson(N,DMe_events));
796 +
			if(llh[0]==0.0) llh.push_back(0.0);
797 +
			else llh.push_back( NMC[1]*PDF_Poisson(NMC[0],DMe_events) );
798 +
799 +
		}
800 +
		else if(experiment=="SENSEI" || experiment=="SENSEI-surface"||experiment=="SuperCDMS"||experiment=="DAMIC-M")
801 +
		{
802 +
			//Efficiency and number of events per bin
803 +
				std::vector<double> efficiency;
804 +
				std::vector<int> BinData;
805 +
				if(experiment=="SENSEI")
806 +
				{
807 +
					efficiency={1.0,0.62,0.48};
808 +
					BinData = {8516,87,0};
809 +
				}
810 +
				else if(experiment=="SENSEI-surface")
811 +
				{
812 +
					efficiency={0.668,0.41,0.32,0.27,0.24};
813 +
					BinData={140302,4676,131,1,0};
814 +
				}
815 +
				else if(experiment=="SuperCDMS")
816 +
				{
817 +
					//all points within 2sigma, taken from fig 3 of 1804.10697
818 +
					efficiency={0.88,0.91,0.91,0.91,0.91,0.91};
819 +
					BinData={53000, 400, 74, 18, 7, 14};
820 +
				}
821 +
				else if(experiment=="DAMIC-M")
822 +
				{
823 +
					efficiency={1.0,1.0,1.0,1.0,1.0,1.0};
824 +
					BinData={100000,0,0,0,0,0};
825 +
				}
826 +
			//Compute spectrum
827 +
				Interpolation spectrum = dRdEe(DM,DMe_target,DMe_efficiency,attenuation,data,vCutoff);
828 +
			//Compute rates
829 +
				std::vector<double> rates;
830 +
				double Egap = 1.11*eV;
831 +
				double eps = 3.6*eV;
832 +
				for(unsigned int bin=DMe_threshold;bin<BinData.size()+1;bin++) 
833 +
				{
834 +
					double E1= eps*(bin-1.0)+Egap;
835 +
					double E2= eps*(bin)+Egap;
836 +
					double epsilon = 1.0e-5*(E2-E1)*spectrum(E1);
837 +
					double rate = Integrate(spectrum,E1,E2,epsilon);
838 +
					rates.push_back(rate);
839 +
				}
840 +
			//Compute number of events and likelihood per bin, and total number of events
841 +
				std::vector<double> events;
842 +
				std::vector<double> llhs;
843 +
				double N = 0;
844 +
				for(unsigned int bin=(DMe_threshold-1);bin<BinData.size();bin++) 
845 +
				{
846 +
					events.push_back(efficiency[bin]*DMe_exposure*rates[bin]);
847 +
					llhs.push_back(CDF_Poisson(events[bin],BinData[bin]));
848 +
					N+=events[bin];
849 +
				}
850 +
			//Total number of events
851 +
				double error = attenuation[1]/attenuation[0]*N;
852 +
				NMC.push_back( N );
853 +
				NMC.push_back( error );
854 +
			//Find minimum value
855 +
				llh.push_back( *std::min_element(llhs.begin(),llhs.end()) );
856 +
				llh.push_back(0.0);
857 +
858 +
		}
859 +
		else if(experiment=="XENON10e")
860 +
		{
861 +
			std::vector<double> dRdn=dRdne(DM,DMe_target,Shells,attenuation,data,vCutoff);
862 +
			double muPE = 27.0;
863 +
			double sigPE=6.7;
864 +
			std::vector<double> dRdnpe=dRdnPE(muPE,sigPE,dRdn);
865 +
			std::vector<double> dNdnpe=dNdnPE(experiment,dRdnpe);
866 +
			double N=0.0;
867 +
			//Data is binned in 7 bins [14,41),[41,68),[68,95),[95,122),[122,149),[149,176),[176,203)
868 +
				std::vector<int> events={126, 60, 12, 3, 2, 0, 2};
869 +
				std::vector<unsigned int> bins={14,41,68,95,122,149,176,203};
870 +
				std::vector<double> llhs;
871 +
872 +
			//Compute events in each bin.
873 +
				for(unsigned int bin=0;bin<events.size();bin++)
874 +
				{
875 +
					double Nbin=0;
876 +
					for(unsigned int PE= bins[bin];PE<bins[bin+1];PE++)
877 +
					{
878 +
						Nbin+=dNdnpe[PE-1];						
879 +
					}
880 +
					N+=Nbin;
881 +
					double l=CDF_Poisson(Nbin,events[bin]);
882 +
					llhs.push_back(l);
883 +
884 +
				}
885 +
			double error = attenuation[1]/attenuation[0]*N;
886 +
			double llh_min = *std::min_element(llhs.begin(),llhs.end());
887 +
888 +
			NMC.push_back( N );
889 +
			NMC.push_back( error );
890 +
			llh.push_back(llh_min);
891 +
			if(llh[0]==0.0) llh.push_back(0.0);
892 +
			else llh.push_back( NMC[1]*exp(-NMC[0]) );
893 +
894 +
		}
895 +
		else if(experiment=="XENON100e")
896 +
		{
897 +
			std::vector<double> dRdn=dRdne(DM,DMe_target,Shells,attenuation,data,vCutoff);
898 +
			double muPE = 19.7;
899 +
			double sigPE=6.2;
900 +
			std::vector<double> dRdnpe=dRdnPE(muPE,sigPE,dRdn);
901 +
			std::vector<double> dNdnpe=dNdnPE(experiment,dRdnpe);
902 +
			double N=0.0;
903 +
			//Data is binned in 3 bins [80,90),[90,110),[110,130)
904 +
				std::vector<int> events={794, 1218, 924, 776, 669, 630, 528, 488, 433, 387};
905 +
				std::vector<unsigned int> bins={80, 90, 110,130,150,170,190,210,230,250,270};
906 +
				std::vector<double> llhs;
907 +
908 +
			//Compute events in each bin for the first 3 bins.
909 +
				for(unsigned int bin=0;bin<3;bin++)
910 +
				{
911 +
					double Nbin=0;
912 +
					for(unsigned int PE= bins[bin];PE<bins[bin+1];PE++)
913 +
					{
914 +
						Nbin+=dNdnpe[PE-1];
915 +
					}
916 +
					N+=Nbin;
917 +
					double l=CDF_Poisson(Nbin,events[bin]);
918 +
					llhs.push_back(l);
919 +
920 +
				}
921 +
			double error = attenuation[1]/attenuation[0]*N;
922 +
			double llh_min = *std::min_element(llhs.begin(),llhs.end());
923 +
924 +
			NMC.push_back( N );
925 +
			NMC.push_back( error );
926 +
			llh.push_back(llh_min);
927 +
			if(llh[0]==0.0) llh.push_back(0.0);
928 +
			else llh.push_back( NMC[1]*exp(-NMC[0]) );
929 +
		}
930 +
		else if(experiment=="DarkSide-50")
931 +
		{
932 +
			//Experiment's parameter
933 +
				std::vector<int>events={118643, 219893, 6131, 673, 252, 227, 198, 199, 189, 247, 230, 261,249, 329, 336, 349, 351, 352, 384, 411, 405, 461, 460, 436, 500, 546, 538, 536, 556, 583, 573, 630, 603, 635, 639, 682, 736, 755, 804, 811, 809, 882, 934, 935, 871, 965, 946, 1072, 997, 1060};
934 +
				double exposure = 6786.0*kg*day;
935 +
				int nThreshold=3;
936 +
			//Spectrum
937 +
				std::vector<double> dRdn=dRdne(DM,DMe_target,Shells,attenuation,data,vCutoff);
938 +
			//Compute likelihood for each bin.
939 +
				std::vector<double> llhs;
940 +
				double Ntot=0.0;
941 +
				for(unsigned int ne=nThreshold;ne<events.size();ne++)
942 +
				{
943 +
					// double N=exposure * dRdn[ne-1];
944 +
					double N=exposure/2.0 * (dRdn[ne-1]+dRdn[ne]);
945 +
					Ntot+=N;
946 +
					double l=CDF_Poisson(N,events[ne-1]);
947 +
					llhs.push_back(l);
948 +
				}
949 +
			//Find the minimum likelihood.	
950 +
				double llh_min = *std::min_element(llhs.begin(),llhs.end());
951 +
				double error = attenuation[1]/attenuation[0]*Ntot;
952 +
953 +
			NMC.push_back( Ntot );
954 +
			NMC.push_back( error );
955 +
			llh.push_back(llh_min);
956 +
			if(llh[0]==0.0) llh.push_back(0.0);
957 +
			else llh.push_back( NMC[1]*exp(-NMC[0]) );
958 +
		}
959 +
		else
960 +
		{
961 +
			cerr <<"Error in Likelihood(): Experiment " <<experiment <<" not recognized."<<endl;
962 +
			std::exit(EXIT_FAILURE);
963 +
		}
964 +
		return llh;
965 +
	}
966 +
967 +
//5. Total number of events
968 +
	double N_Signals(const std::function<double(double)>& spectrum,double E1,double E2,double exposure)
969 +
	{
970 +
		double epsilon = Find_Epsilon(spectrum,E1,E2,1e-4);
971 +
		return exposure*Integrate(spectrum, E1,E2,epsilon);
972 +
	}
973 +
974 +
//6. Standard constraint
975 +
	 double Upper_Bound(const DM_Particle& DM,const std::string& experiment,double CL)
976 +
	{
977 +
	 	double sigma=-1.0;
978 +
979 +
		//DAMIC
980 +
		if(experiment=="DAMIC")
981 +
		{
982 +
			double N_Limit = Inv_CDF_Poisson(106,1.0-CL);
983 +
			double exposure = 0.107*kg*day;
984 +
			double Ethreshold = 0.55*keV;
985 +
			double Emax = 7*keV;
986 +
			DM_Particle dm=DM;
987 +
			dm.Set_Sigma_n(pb);
988 +
			std::function <double(double)> spectrum = Compute_Spectrum(dm,experiment);
989 +
			double N = N_Signals(spectrum,Ethreshold,Emax,exposure);
990 +
			sigma = N_Limit/N*pb;
991 +
		}
992 +
		else if(experiment=="XENON1T")
993 +
		{
994 +
			double N_Limit = -log(1.0-CL);
995 +
			double exposure = 34.2*day*1042*kg;
996 +
			double Ethreshold = 5*keV;
997 +
			double Emax = 40*keV;
998 +
			DM_Particle dm=DM;
999 +
			dm.Set_Sigma_n(pb);
1000 +
			std::function <double(double)> spectrum = Compute_Spectrum(dm,experiment);
1001 +
			double N = N_Signals(spectrum,Ethreshold,Emax,exposure);
1002 +
			sigma = N_Limit/N*pb;
1003 +
		}
1004 +
		else if(experiment == "CRESST-II")
1005 +
		{
1006 +
			//Find critical cross-section using bisection
1007 +
				double n;
1008 +
				double s1=1e-45*cm*cm;
1009 +
				double s2=1e-30*cm*cm;
1010 +
				double s3 = pow(10.0, log10(s1*s2)/2.0);
1011 +
				double e = 0.01;
1012 +
				double llh=0.0;
1013 +
				while(fabs(llh-(1.0-CL))>e)
1014 +
				{
1015 +
					DM_Particle dm=DM;
1016 +
					dm.Set_Sigma_n(s3);
1017 +
					Interpolation spectrum = Compute_Spectrum(dm,experiment);
1018 +
					llh = Likelihood(dm,experiment,spectrum,n);
1019 +
					if(llh>(1.0-CL))
1020 +
					{
1021 +
						s1 = s3;
1022 +
					}
1023 +
					else
1024 +
					{
1025 +
						s2 = s3;
1026 +
					}
1027 +
					s3 = pow(10.0, log10(s1*s2)/2.0);
1028 +
				}
1029 +
				sigma=s3;
1030 +
		}
1031 +
		else if(experiment == "CRESST-surface")
1032 +
		{
1033 +
			//Find critical cross-section using bisection
1034 +
				double n;
1035 +
				double s1=1e-36*cm*cm;
1036 +
				double s2=1e-28*cm*cm;
1037 +
				double s3 = pow(10.0, log10(s1*s2)/2.0);
1038 +
				double e = 0.001;
1039 +
				double llh=0.0;
1040 +
				while(fabs(llh-(1.0-CL))>e)
1041 +
				{
1042 +
					DM_Particle dm=DM;
1043 +
					dm.Set_Sigma_n(s3);
1044 +
					Interpolation spectrum = Compute_Spectrum(dm,experiment);
1045 +
					llh = Likelihood(dm,experiment,spectrum,n);
1046 +
					if(llh>(1.0-CL))
1047 +
					{
1048 +
						s1 = s3;
1049 +
					}
1050 +
					else
1051 +
					{
1052 +
						s2 = s3;
1053 +
					}
1054 +
					s3 = pow(10.0, log10(s1*s2)/2.0);
1055 +
				}
1056 +
				sigma=s3;
1057 +
		}
1058 +
		else if(experiment == "Semiconductor")
1059 +
		{
1060 +
			double N_Limit =Inv_CDF_Poisson(DMe_events,1.0-CL);
1061 +
			DM_Particle dm=DM;
1062 +
			dm.Set_Sigma_e(pb);
1063 +
			double N = DMe_exposure*TotalRate(dm,DMe_threshold,DMe_efficiency,DMe_target);
1064 +
			sigma = N_Limit/N*pb;
1065 +
		}
1066 +
		else if(experiment =="SENSEI-surface"||experiment=="SENSEI"||experiment=="SuperCDMS"||experiment=="DAMIC-M")
1067 +
		{
1068 +
			//Efficiency and number of events per bin
1069 +
				std::vector<double> efficiency;
1070 +
				std::vector<int> BinData;
1071 +
				if(experiment=="SENSEI")
1072 +
				{
1073 +
					efficiency={1.0,0.62,0.48};
1074 +
					BinData = {8516,87,0};
1075 +
1076 +
				}
1077 +
				else if(experiment=="SENSEI-surface")
1078 +
				{
1079 +
					efficiency={0.668,0.41,0.32,0.27,0.24};
1080 +
					BinData={140302,4676,131,1,0};
1081 +
				}
1082 +
				else if(experiment=="SuperCDMS")
1083 +
				{
1084 +
					// all points within 2sigma, taken from fig 3 of 1804.10697
1085 +
					efficiency={0.88,0.91,0.91,0.91,0.91,0.91};
1086 +
					BinData={53000, 400, 74, 18, 7, 14};
1087 +
				}
1088 +
				else if(experiment=="DAMIC-M")
1089 +
				{
1090 +
					efficiency={1.0,1.0,1.0,1.0,1.0,1.0};
1091 +
					BinData={100000,0,0,0,0,0};
1092 +
				}
1093 +
			//Compute rates for 1 pb
1094 +
				std::vector<double> rates;
1095 +
				DM_Particle dm=DM;
1096 +
				dm.Set_Sigma_e(pb);
1097 +
				for(unsigned int bin=1;bin<=(BinData.size()+1);bin++) 
1098 +
				{					
1099 +
					rates.push_back(TotalRate(dm,bin,DMe_efficiency,DMe_target));
1100 +
				}
1101 +
			//Compute number of events and limit per bin
1102 +
				std::vector<double>sigmas;
1103 +
				for(unsigned int bin=(DMe_threshold-1);bin<BinData.size();bin++)
1104 +
				{
1105 +
					double N_Limit=Inv_CDF_Poisson(BinData[bin],1.0-CL);
1106 +
					double Nbin =efficiency[bin]*DMe_exposure*(rates[bin]-rates[bin+1]);
1107 +
					double sig = N_Limit/Nbin*pb;
1108 +
					sigmas.push_back(sig);
1109 +
				}
1110 +
			//Find minimum value
1111 +
				sigma = *std::min_element(sigmas.begin(),sigmas.end());
1112 +
		}
1113 +
		else if(experiment == "XENON10e")
1114 +
		{
1115 +
			//Data is binned in 7 bins [14,41),[41,68),[68,95),[95,122),[122,149),[149,176),[176,203)
1116 +
				std::vector<int> data={126, 60, 12, 3, 2, 0, 2};
1117 +
				std::vector<unsigned int> bins={14,41,68,95,122,149,176,203};
1118 +
1119 +
				std::vector<double>sigmas;
1120 +
			//Compute events in each bin with 1pb.
1121 +
				DM_Particle dm=DM;
1122 +
				dm.Set_Sigma_e(pb);
1123 +
				for(unsigned int bin=0;bin<data.size();bin++)
1124 +
				{
1125 +
					double N_Limit=Inv_CDF_Poisson(data[bin],1.0-CL);
1126 +
					double Nbin=0;
1127 +
					for(unsigned int PE= bins[bin];PE<bins[bin+1];PE++)
1128 +
					{
1129 +
						Nbin+=dNdnPE(PE,dm,DMe_target,experiment,Shells);
1130 +
					}
1131 +
					
1132 +
					sigmas.push_back(N_Limit/Nbin*pb);
1133 +
				}
1134 +
			sigma = *std::min_element(sigmas.begin(),sigmas.end());
1135 +
		}
1136 +
		else if(experiment == "XENON100e")
1137 +
		{
1138 +
			//Data is binned in 3 bins [80,90),[90,110),[110,130)
1139 +
				std::vector<int> data={794, 1218, 924, 776, 669, 630, 528, 488, 433, 387};
1140 +
				std::vector<unsigned int> bins={80, 90, 110,130,150,170,190,210,230,250,270};
1141 +
				
1142 +
				//Inverting the Poisson CDF is not straight forward (look at inverse regularized gamma functions at some point)
1143 +
					std::vector<double> Nl;
1144 +
					for(unsigned int j=0;j<data.size();j++) Nl.push_back(Inv_CDF_Poisson(data[j],1.0-CL));
1145 +
				
1146 +
				std::vector<double>sigmas;
1147 +
			//Compute events in each bin with 1pb.
1148 +
				//We only use the first 3 bins.
1149 +
				DM_Particle dm=DM;
1150 +
				dm.Set_Sigma_e(pb);
1151 +
				for(unsigned int bin=0;bin<3;bin++)
1152 +
				{
1153 +
					double N_Limit=Nl[bin];
1154 +
					double Nbin=0;
1155 +
					for(unsigned int PE= bins[bin];PE<bins[bin+1];PE++)
1156 +
					{
1157 +
						Nbin+=dNdnPE(PE,dm,DMe_target,experiment,Shells);
1158 +
					}
1159 +
					sigmas.push_back(N_Limit/Nbin*pb);
1160 +
				}
1161 +
			sigma = *std::min_element(sigmas.begin(),sigmas.end());
1162 +
		}
1163 +
		else if (experiment=="DarkSide-50")
1164 +
		{
1165 +
			//Threshold and exposure
1166 +
				double exposure = 6786.0*kg*day;
1167 +
				int nThreshold=3;
1168 +
			//Data
1169 +
				std::vector<int>data={118643, 219893, 6131, 673, 252, 227, 198, 199, 189, 247, 230, 261,249, 329, 336, 349, 351, 352, 384, 411, 405, 461, 460, 436, 500, 546, 538, 536, 556, 583, 573, 630, 603, 635, 639, 682, 736, 755, 804, 811, 809, 882, 934, 935, 871, 965, 946, 1072, 997, 1060};
1170 +
			//Inverting the Poisson CDF
1171 +
				std::vector<double> Nl;
1172 +
				for(unsigned int j=0;j<data.size();j++) Nl.push_back(Inv_CDF_Poisson(data[j],1.0-CL));
1173 +
1174 +
			std::vector<double>sigmas;
1175 +
			//Compute events and limit in each bin for 1pb.
1176 +
				DM_Particle dm=DM;
1177 +
				dm.Set_Sigma_e(pb);
1178 +
				for(unsigned int ne=nThreshold;ne<data.size();ne++)
1179 +
				{
1180 +
					double N_Limit = Nl[ne-1];
1181 +
					double Nbin=0.0;
1182 +
					//Add all shells
1183 +
					for(unsigned int shell=0;shell<Shells.size();shell++)
1184 +
					{					
1185 +
						Nbin += 1.0/2.0*exposure*(dRdne(ne,dm,DMe_target,Shells[shell])+dRdne(ne+1,dm,DMe_target,Shells[shell]));
1186 +
					}
1187 +
					if(Nbin>0)
1188 +
					{
1189 +
						sigmas.push_back(N_Limit/Nbin*pb);
1190 +
					}
1191 +
					else
1192 +
					{
1193 +
						break;
1194 +
					}
1195 +
					
1196 +
					
1197 +
				}
1198 +
			//Take the minimal value
1199 +
				sigma = *std::min_element(sigmas.begin(),sigmas.end());
1200 +
		}
1201 +
		else
1202 +
		{
1203 +
			std::cerr <<"Error in Likelihood(): Experiment " <<experiment <<" not recognized."<<endl;
1204 +
			std::exit(EXIT_FAILURE);
1205 +
		}
1206 +
		return sigma;
1207 +
	}
1208 +
1209 +
	std::vector<std::vector<double>> Limit_Curve(DM_Particle DM,const std::string& experiment,const std::string& simID,double mMin,double mMax,int Nm,double CL,int rank)
1210 +
	{
1211 +
		if(rank == 0)cout <<"Compute analytic limits."<<endl;
1212 +
		std::vector<std::vector<double>> output;
1213 +
		//Time
1214 +
			std::chrono::high_resolution_clock::time_point tS,tE;
1215 +
		  	tS = std::chrono::high_resolution_clock::now();
1216 +
1217 +
		//Output file
1218 +
			ofstream f;
1219 +
			if(rank == 0) f.open("../results/"+simID+"/Limits_Analytic.txt");
1220 +
		//Mass scan
1221 +
			double dlogm = (Nm==1)? 0 : (log10(mMax)-log10(mMin))/(Nm-1);
1222 +
			for(int i=0;i<Nm;i++)
1223 +
			{
1224 +
				double mDM = pow(10.0,log10(mMin)+i*dlogm);
1225 +
				DM.Set_Mass(mDM);
1226 +
				double sigma = Upper_Bound(DM,experiment,CL);
1227 +
				output.push_back(std::vector<double> {mDM,sigma});
1228 +
				if(rank==0)
1229 +
				{
1230 +
					f<<mDM <<"\t" <<InUnits(sigma,cm*cm)<<endl;
1231 +
					cout <<"\r                                                     \r" 
1232 +
					<<Round(mDM) <<" GeV\t" <<Round(InUnits(sigma,cm*cm))<<"cm^2\t("<<floor(100.0*i/Nm)<<"\%)"<<flush;
1233 +
				}
1234 +
			}
1235 +
		//Finish
1236 +
			tE= std::chrono::high_resolution_clock::now();
1237 +
			double duration =1e-6*std::chrono::duration_cast<std::chrono::microseconds>( tE - tS ).count();
1238 +
			if(rank==0)
1239 +
			{
1240 +
				cout <<"\rDone.  ("<<Round(duration)<<"sec)                                 "<<endl;
1241 +
				f.close();
1242 +
			}	
1243 +
			return output;
1244 +
	}

@@ -1,35 +1,34 @@
Loading
1 +
//Contains the actual trajectory simulation function.
1 2
#ifndef __Trajector_Simulation_hpp_
2 3
#define __Trajector_Simulation_hpp_
3 4
4 -
#include <Eigen/Geometry>
5 5
#include <functional>
6 6
#include <random>
7 7
8 -
#include "Physical_Parameters.hpp"
9 -
#include "General_Utilities.hpp"
8 +
#include "Physics_Functions.hpp"
9 +
#include "Simulation_Essentials.hpp"
10 10
11 -
#include "Layer_Class.hpp"
12 -
#include "RN_Generators.hpp"
13 -
#include "Trajectory_Class.hpp"
11 +
//Cut off speed for a given experiment
12 +
	extern double Cutoff_Speed(const std::string& experiment,double mDM,int rank=0);
14 13
15 14
//Initial Condition Generator
16 -
	extern Event InitialCondition(double tIni,Eigen::Vector3d& xIni,std::mt19937& PRNG,double vCutoff = 0.0);
17 -
18 -
//Propagate a particle until it scatters, reaches the detector or gets back into space
19 -
	extern void Free_Propagation(Event& event,double& weight,std::vector<std::vector<double>>& prob,std::mt19937& PRNG,double logXi = -1.0);
20 -
21 -
//Find the scattering nucleus species:
22 -
	extern std::vector<double> ScatterNucleus(std::vector<std::vector<double>>& prob,std::mt19937& PRNG);
23 -
24 -
//Perform the scattering.
25 -
	extern double PDF_ScatterAngle(double cosalpha,double A,double mDM,double vDM,double deltaIS);
26 -
	extern void Scatter(std::vector<std::vector<double>>& prob,Event& event, double& weight,double mDM,std::mt19937& PRNG);
15 +
	extern Event InitialCondition(double tIni,Vector3D xIni,std::mt19937& PRNG,double vCutoff);
27 16
28 17
//Simulate one particle track:
29 -
	extern Trajectory ParticleTrack(double mDM,double sigman0,Event& IniCondi, double vcut, std::mt19937& PRNG,bool ldm);
30 -
31 -
//Generate velocity data
32 -
	extern std::vector<DataPoint> Simulate_Trajectories(int SampleSize,double mDM, double sigma,double vCutoff,double& hMax,long long unsigned int& ParticleCounter,std::mt19937& PRNG,bool ldm,int rank = 0);
33 -
34 -
18 +
	struct Result
19 +
	{
20 +
		bool success;
21 +
		Event final_event;
22 +
		long unsigned int nScattering;
23 +
		double max_horizontal_distance;
24 +
		double weight;
25 +
26 +
		//Constructors
27 +
		Result();
28 +
		Result(double w);
29 +
30 +
		void Summary();
31 +
	};
32 +
	extern std::vector<Result> Simulate_Trajectory(DM_Particle& DM, Event IC, double vMin,std::mt19937& PRNG,Result result=Result(),std::string split_ID ="Contact");
33 +
	
35 34
#endif
36 35
imilarity index 64%
37 36
ename from results/CRESST-II/CRESST-II_Limits.txt
38 37
ename to results/CRESST-II/Limits_MC.txt

@@ -0,0 +1,917 @@
Loading
1 +
//Disclaimer:
2 +
//Some of the function implemenations were made with the help of the 
3 +
//"Numerical Recipes 3rd Edition: The Art of Scientific Computing"
4 +
//by William H. Press, Saul A. Teukolsky, William T. Vetterling, Brian P. Flannery
5 +
6 +
#include "Numerics_Functions.hpp"
7 +
8 +
#include <iostream>
9 +
#include <fstream>
10 +
#include <cmath>
11 +
#include <cstdlib>
12 +
#include <algorithm>
13 +
14 +
using namespace std;
15 +
16 +
//1. Special functions
17 +
	//1.1 Simple functions
18 +
		int Sign(double arg)
19 +
		{
20 +
			if(arg>0.0) 		return 1;
21 +
			else if(arg==0.0)	return 0;
22 +
			else 				return -1;
23 +
		}
24 +
		double Sign(double x, double y)
25 +
		{
26 +
			if(Sign(x) == Sign(y)) return x;
27 +
			else return -1.0*x;
28 +
		}
29 +
		double StepFunction(double x)
30 +
		{
31 +
			if(x>=0) 	return 1.0;
32 +
			else		return 0.0;
33 +
		}
34 +
		//This function should only be used to 
35 +
		double Round(double N,unsigned int digits)
36 +
		{
37 +
			if(N==0) return 0;
38 +
			if(digits>5)
39 +
			{
40 +
				std::cerr <<"Error in Round(): Significant digits > 5."<<std::endl;
41 +
				std::exit(EXIT_FAILURE);
42 +
			}
43 +
			//Make the argument a positive number.
44 +
				double sign = Sign(N);
45 +
				N*=sign;
46 +
			//Cut of the decimal power
47 +
				double DecimalPower = floor( log10(N) );
48 +
				double prefactor = N*pow(10,-DecimalPower);
49 +
			//Cut off all unneeded decimal points, which might cause rounding error
50 +
				prefactor = (int) floor(0.5+1.0* pow(10,digits)*prefactor);
51 +
				prefactor/=10;
52 +
			//Round the prefactor with k significant digits
53 +
				prefactor = floor(0.5 +1.0*prefactor);
54 +
				prefactor/= pow(10,digits-1);
55 +
			//Final output
56 +
				return sign*prefactor*pow(10,DecimalPower);
57 +
		}
58 +
59 +
	//1.2 Gamma Functions from "Numerical Recipes"
60 +
		//Factorial
61 +
			std::vector<double> FactorialList= {1.0};
62 +
			double Factorial(unsigned int n)
63 +
			{
64 +
				if (n>170)
65 +
				{
66 +
67 +
					std::cerr <<"Error in Factorial: Overflow for " <<n <<"!."<<std::endl;
68 +
					std::exit(EXIT_FAILURE);
69 +
				}
70 +
				else if(n<FactorialList.size()) return FactorialList[n];
71 +
				else
72 +
				{
73 +
					while(FactorialList.size()<=n)
74 +
					{
75 +
						FactorialList.push_back( FactorialList.back()*FactorialList.size());
76 +
					}
77 +
					return FactorialList.back();
78 +
				}
79 +
			}
80 +
		//Binomial coefficients
81 +
			double BinomialCoefficient(int n,int k)
82 +
			{
83 +
				if(k<0 ||n<0 )
84 +
				{
85 +
					std::cerr <<"Warning in BinomialCoefficient(): negative arguments. Return 0." <<std::endl;
86 +
					std::exit(EXIT_FAILURE);
87 +
				}
88 +
				else if(n<k)	return 0;
89 +
				else if (n>170)
90 +
				{
91 +
					return floor(0.5+exp(GammaLn(n+1.0)-GammaLn(k+1.0)-GammaLn(n-k+1.0)));
92 +
				}
93 +
				else	return floor(0.5+Factorial(n)/Factorial(k)/Factorial(n-k));
94 +
			}
95 +
		//Logarithmic gamma function
96 +
			double GammaLn(double x)
97 +
			{
98 +
				double cof[14] = {57.1562356658629235,-59.5979603554754912, 14.1360979747417471,-0.491913816097620199,.339946499848118887e-4, .465236289270485756e-4,-.983744753048795646e-4,.158088703224912494e-3, -.210264441724104883e-3,.217439618115212643e-3,-.164318106536763890e-3, .844182239838527433e-4,-.261908384015814087e-4,.368991826595316234e-5};
99 +
				if(x<=0)
100 +
				{
101 +
					std::cerr<<"Error in GammaLn(x): x<=0."<<std::endl;
102 +
					std::exit(EXIT_FAILURE);
103 +
				}
104 +
				double sum = 0.999999999999997092;
105 +
				double y = x;
106 +
				double tmp = x+671.0/128.0;
107 +
				tmp = (x+0.5)*log(tmp)-tmp;
108 +
				for(int j=0;j<14;j++)
109 +
				{
110 +
					sum+=cof[j]/++y;
111 +
				}
112 +
				return tmp+log(2.5066282746310005*sum/x);
113 +
			}
114 +
		//Q(x,a) via integration;
115 +
			double GammaQint(double x,double a)
116 +
			{
117 +
				//Compute P(x,a)
118 +
					double gammaP;
119 +
					double gln=GammaLn(a);
120 +
					//How far to integrate N sqrt(a) around the peak at a-1:
121 +
						double N = 10;
122 +
						double tPeak = a-1.0;
123 +
						double tMin = std::max(0.0, tPeak - N*sqrt(a));
124 +
						double tMax = tPeak + N*sqrt(a);
125 +
		 			//If x lies far away from the peak:
126 +
						if(x>tMax) gammaP = 1.0;
127 +
						else if (x<tMin ) gammaP = 0.0;
128 +
					//Numerical integration
129 +
						else
130 +
						{
131 +
							
132 +
							//integrand
133 +
								std::function<double(double)> integrand = [a,gln] (double t)
134 +
								{
135 +
									return exp(-gln-t+log(t)*(a-1.0));
136 +
								};
137 +
								if(x<tMin) tMin = 0.0;
138 +
							//Precision
139 +
								double eps = Find_Epsilon(integrand,tMin,x,1e-5);
140 +
							//Integrate
141 +
								gammaP = Integrate(integrand,tMin,x,eps);
142 +
						}
143 +
				
144 +
				return 1.0-gammaP;
145 +
			}
146 +
		//Series expansion of P(x,a)
147 +
			double GammaPser(double x,double a)
148 +
			{
149 +
				double eps = std::numeric_limits<double>::epsilon();
150 +
				double sum,del,ap;
151 +
				double gln=GammaLn(a);
152 +
				ap = a;
153 +
				del = sum = 1.0/a;
154 +
				while(fabs(del)>fabs(sum)*eps)
155 +
				{
156 +
					ap++;
157 +
					del*=x/ap;
158 +
					sum+=del;
159 +
				}
160 +
				return sum*exp(-x+a*log(x)-gln);
161 +
			}
162 +
		//Continued fraction representation of Q(x,a)
163 +
			double GammaQcf(double x,double a)
164 +
			{
165 +
				//Precision
166 +
					double eps = std::numeric_limits<double>::epsilon();
167 +
					double FPMIN = std::numeric_limits<double>::min()/eps;
168 +
				double del=0.0;
169 +
				double gln=GammaLn(a);
170 +
				double b = x+1.0-a;
171 +
				double c=1.0/FPMIN;
172 +
				double d = 1.0/b;
173 +
				double h=d;
174 +
				int i=1;
175 +
				while(fabs(del-1.0)>eps)
176 +
				{
177 +
					double an = -1.0*i*(i-a);
178 +
					b+=2.0;
179 +
					d=an*d+b;
180 +
					if(fabs(d)<FPMIN) d=FPMIN;
181 +
					c=b+an/c;
182 +
					if(fabs(c)<FPMIN) c =FPMIN;
183 +
					d=1.0/d;
184 +
					del=d*c;
185 +
					h*=del;
186 +
				}
187 +
				return exp(-x+a*log(x)-gln)*h;
188 +
			}
189 +
		//Final function using different methods for different parts of the domain
190 +
			double GammaQ(double x,double a)
191 +
			{
192 +
				double aMax=100.0;
193 +
				if(x<0.0||a<=0.0)
194 +
				{
195 +
					std::cerr <<"Error in GammaQ("<<x<<","<<a<<"): Invalid arguments."<<std::endl;
196 +
					std::exit(EXIT_FAILURE);
197 +
				}
198 +
				else if(x==0)return 1.0;
199 +
				else if (a>aMax) return GammaQint(x,a);
200 +
				else if (x<a+1.0) return 1.0-GammaPser(x,a);
201 +
				else return GammaQcf(x,a);
202 +
			}
203 +
			double GammaP(double x,double a)
204 +
			{
205 +
				return 1.0-GammaQ(x,a);
206 +
			}
207 +
			
208 +
		//Inverse incomplete gamma function. 
209 +
			//Solves P(x,a)=p for x.
210 +
				double Inv_GammaP(double p,double a)
211 +
				{
212 +
					
213 +
					//Check the arguments
214 +
						if(a<=0.0)
215 +
						{
216 +
							std::cerr <<"Error in Inv_GammaP(): a must be positive."<<std::endl;
217 +
							std::exit(EXIT_FAILURE);
218 +
						}
219 +
						if(p>=1.0) return std::max(100.0,a+100.*sqrt(a));
220 +
						if(p<=0.0) return 0.0;
221 +
					//Parameter
222 +
						double x;
223 +
						double gln=GammaLn(a);
224 +
						double a1=a-1.0;
225 +
						double lna1 = log(a1);
226 +
						double afac = exp(a1*(lna1-1.0)-gln);
227 +
					//Initial guess 1
228 +
						if(a>1.0)
229 +
						{
230 +
							double pp = (p<0.5)? p : 1.0-p;
231 +
							double t = sqrt(-2.0*log(pp));
232 +
							x = (2.30753+t*0.27061)/(1.+t*(0.99229+t*0.04481)) - t;
233 +
							if(p<0.5) x=-x;
234 +
							x = std::max(1.0e-3,a*pow(1.0-1.0/(9.*a)-x/(3.*sqrt(a)),3.0));
235 +
236 +
						}
237 +
					//Initial guess 2
238 +
						else
239 +
						{
240 +
							double t = 1.0 - a*(0.253+a*0.12);
241 +
					        if (p < t) x = pow(p/t,1./a);
242 +
					        else x = 1.-log(1.-(p-t)/(1.-t));
243 +
						}
244 +
					//Halley's method
245 +
						double EPS = 1.0e-8;
246 +
						for(int i=0;i<12;i++)
247 +
						{
248 +
							if(x<=0.0) return 0.0;
249 +
							double error = GammaP(x,a) - p;
250 +
							double t;
251 +
							if(a>1.0) t = afac*exp(-(x-a1)+a1*(log(x)-lna1));
252 +
							else t = exp(-x+a1*log(x)-gln);
253 +
							double u = error/t;
254 +
							x-= (t = u/(1.-0.5*std::min(1.,u*((a-1.)/x - 1))));
255 +
							if (x <= 0.) x = 0.5*(x + t);
256 +
							if (fabs(t) < EPS*x ) break;
257 +
						}
258 +
					return x;
259 +
				}
260 +
			//Solves Q(x,a)=p for x.
261 +
				double Inv_GammaQ(double q,double a)
262 +
				{
263 +
					return Inv_GammaP(1.0-q,a);
264 +
				}
265 +
266 +
267 +
268 +
	//1.3 Statistics
269 +
		double PDF_Gauss(double x,double mu,double sigma)
270 +
		{
271 +
			return 1.0/sqrt(2.0*M_PI)/sigma*exp(-pow((x-mu)/sigma,2.0)/2.0);
272 +
		}
273 +
		double PDF_Binomial(int mu,double p,int trials)
274 +
		{
275 +
			return BinomialCoefficient(mu,trials)*pow(p,trials)*pow(1.0-p,(mu-trials));
276 +
		}
277 +
		double PDF_Poisson(double expected_events, unsigned int events)
278 +
		{
279 +
			if(expected_events==0&&events==0) return 1.0;
280 +
			else if(expected_events==0&&events>0) return 0.0;
281 +
			else
282 +
			{
283 +
				double sum=events*log(expected_events) -expected_events;
284 +
				for(unsigned int i=2;i<=events;i++)sum-= log(i);
285 +
				return exp(sum);
286 +
			}
287 +
		}
288 +
		double CDF_Poisson(double expectation_value,unsigned int observed_events)
289 +
		{
290 +
			double gq = GammaQ(expectation_value,observed_events+1);
291 +
			if (gq>=0) return gq;
292 +
			else return 0.0;
293 +
		}
294 +
		double Inv_CDF_Poisson(unsigned int observed_events, double cdf)
295 +
		{
296 +
			if(observed_events==0) return (-1.0)*log(cdf);
297 +
			else return Inv_GammaQ(cdf,observed_events+1); 	
298 +
		}
299 +
		//Maximum Gap a la Yellin
300 +
			double MaxGap_C0(double x,double mu)
301 +
			{
302 +
				if(x==mu) return 1.0-exp(-mu);
303 +
				else
304 +
				{
305 +
					int m = mu/x;
306 +
					double sum=0.0;
307 +
					for(int k=0;k<=m;k++) 
308 +
					{
309 +
						double term = pow(k*x-mu,k)/Factorial(k)*exp(-k*x)*(1.0+k/(mu-k*x));
310 +
						sum+= term ;
311 +
						if (fabs(term)<1e-20) break;
312 +
					}
313 +
					return sum;
314 +
				}
315 +
			}
316 +
			double MaxGap_C0_error(double x,double mu,double mu_error)
317 +
			{
318 +
				if(x==mu) return mu_error*exp(-mu);
319 +
				else
320 +
				{
321 +
					int m = mu/x;
322 +
					double sum=0.0;
323 +
					for(int k=0;k<=m;k++) 
324 +
					{
325 +
						double term = k*pow(k*x-mu,(-2.0+k))/Factorial(k)*exp(-k*x)*(1.0+k*(x-1.0)-mu);
326 +
						sum+= term ;
327 +
						if (fabs(term)<1e-20) break;
328 +
					}
329 +
					return fabs(sum*mu_error);
330 +
				}
331 +
			}
332 +
333 +
//2. Numerical Integration via Adaptive Simpson Method
334 +
	//Function to return a reasonable precision.
335 +
		double Find_Epsilon(std::function<double(double)> func, double a,double b,double precision)
336 +
		{
337 +
			double c = (a + b)/2;
338 +
			double h = b - a;                                                                  
339 +
				double fa = func(a);
340 +
				double fb = func(b);
341 +
				double fc = func(c);                                                           
342 +
				double S = (h/6)*(fa + 4*fc + fb);
343 +
				double epsilon = precision*S;
344 +
				return epsilon;
345 +
		}
346 +
	//Recursive functions for one-dimensional integration
347 +
		double Integrate(std::function<double(double)> func, double a,double b, double epsilon,int maxRecursionDepth)
348 +
		{
349 +
			int sign = +1;
350 +
			if(a==b) return 0.0;
351 +
			else if(a>b)
352 +
			{
353 +
				double aux = a;
354 +
				a = b;
355 +
				b=aux;
356 +
				sign = -1;
357 +
			}
358 +
			double c = (a + b)/2;
359 +
			double h = b - a;                                                                  
360 +
				double fa = func(a);
361 +
				double fb = func(b);
362 +
				double fc = func(c); 
363 +
				double S = (h/6)*(fa + 4*fc + fb);
364 +
				bool warning = false;
365 +
				double result =   auxAdaptiveSimpson(func, a, b, fabs(epsilon), S, fa, fb, fc, maxRecursionDepth,warning); 
366 +
				if(warning)   
367 +
				{
368 +
					std::cout <<"Warning in Integrate(): Numerical integration on the interval ("<<a<<","<<b<<") did not converge to the desired precision." <<std::endl;
369 +
					std::cout <<"\tDesired precision: " <<fabs(epsilon) <<" Result: " <<result<<std::endl;
370 +
				}
371 +
			if(std::isnan(result)) std::cout <<"Warning in Integrate(): Result is nan."<<std::endl;
372 +
			else if(std::isinf(result)) std::cout <<"Warning in Integrate(): Result is inf."<<std::endl;
373 +
			return sign*result;
374 +
		}
375 +
		double auxAdaptiveSimpson(std::function<double(double)> func, double a,double b, double epsilon,double S,double fa,double fb,double fc,int bottom,bool &warning)
376 +
		{
377 +
			double c = (a+b)/2;
378 +
			double h = b-a;
379 +
			double d = (a+c)/2;
380 +
			double e = (b+c)/2;
381 +
			double fd=func(d);
382 +
			double fe=func(e);
383 +
			double Sleft = (h/12)*(fa + 4*fd + fc);                                                           
384 +
			double Sright = (h/12)*(fc + 4*fe + fb);                                                          
385 +
			double S2 = Sleft + Sright; 
386 +
			if (bottom <= 0 || abs(S2 - S) <= 15*epsilon)//15 due to error analysis 
387 +
			{
388 +
				if(bottom<=0&&abs(S2 - S) > 15*epsilon) warning=true;
389 +
				return S2 + (S2 - S)/15;  
390 +
			}                                         
391 +
			else
392 +
			{
393 +
				return auxAdaptiveSimpson(func, a, c, epsilon/2, Sleft,  fa, fc, fd, bottom-1,warning)+auxAdaptiveSimpson(func, c, b, epsilon/2, Sright, fc, fb, fe, bottom-1,warning); 
394 +
			}    
395 +
		}
396 +
397 +
//3. Finding roots
398 +
	//Root finding with Ridder's method
399 +
	double Find_Root(std::function<double(double)>& func,double xLeft, double xRight,double xAccuracy)
400 +
	{
401 +
		const int Max_Iterations = 50;
402 +
		//1. Check if xLeft<xRight, otherwise swap.
403 +
			if(xLeft>xRight)
404 +
			{
405 +
				double temp = xLeft;
406 +
				xLeft = xRight;
407 +
				xRight = temp;
408 +
			}
409 +
		//2. Compute functions at boundary
410 +
			double fLeft = func(xLeft);
411 +
			double fRight = func(xRight);
412 +
		
413 +
		//3. Check if xLeft and xRight bracket a root or already yield a root. Also check for NaN's.
414 +
			if(std::isnan(fLeft)||std::isnan(fRight))
415 +
			{
416 +
				std::cerr <<"Error in Find_Root(): Function returns nan at the brackets."<<std::endl;
417 +
				std::exit(EXIT_FAILURE);
418 +
			}
419 +
			else if(fLeft*fRight>=0.0)
420 +
			{
421 +
				if(fLeft==0) return xLeft;
422 +
				else if(fRight==0) return xRight;
423 +
				else
424 +
				{
425 +
					std::cerr<<"Error in Find_Root(): f(xLeft)*f(xRight) = ("<<func(xLeft)<<")*("<<func(xRight) <<")>0.0"<<std::endl;
426 +
					std::exit(EXIT_FAILURE);
427 +
				}
428 +
			}
429 +
		//4. Ridder's method
430 +
			else
431 +
			{
432 +
				double x1 = xLeft;
433 +
				double x2 = xRight;
434 +
				double f1 = fLeft;
435 +
				double f2 = fRight;
436 +
				double result = -9.9e99;
437 +
				for(int i=0;i<Max_Iterations;i++)
438 +
				{
439 +
					//Mid point
440 +
						double x3 = (x1+x2)/2.0;
441 +
442 +
						double f3 = func(x3);
443 +
					//New point
444 +
						double x4 = x3 +(x3-x1) * Sign(f1-f2)*f3/sqrt(f3*f3-f1*f2);
445 +
					//Check if we found the root
446 +
						if(fabs(x4-result)<xAccuracy) return x4;
447 +
					//Prepare next iteration
448 +
						result = x4;
449 +
						double f4 = func(x4);
450 +
						if(f4==0.0) return result;
451 +
						//a) x3 and x4 bracket the root
452 +
						if(Sign(f3,f4) != f3)
453 +
						{
454 +
							x1 = x3;
455 +
							f1 = f3;
456 +
							x2 = x4;
457 +
							f2 = f4;
458 +
						}
459 +
						//b) x1 and x4 bracket the root
460 +
						else if(Sign(f1,f4) != f1)
461 +
						{
462 +
							x2 = x4;
463 +
							f2 = f4;
464 +
465 +
						}
466 +
						//c) x2 and x4 bracket the root
467 +
						else if(Sign(f2,f4) != f2)
468 +
						{
469 +
							x1 = x4;
470 +
							f1 = f4;
471 +
						}
472 +
						else
473 +
						{
474 +
							std::cerr <<"Error in Find_Root(). Ridder's method does not reach the root."<<std::endl;
475 +
							std::exit(EXIT_FAILURE);
476 +
						}
477 +
				}
478 +
				std::cout <<"Warning in Find_Root(): Iterations exceed the maximum. Final value f("<<result<<")="<<func(result)<<std::endl;
479 +
				return result;
480 +
			}
481 +
	}
482 +
483 +
484 +
//4. One-dimensional interpolation of tabulated data using Steffen splines 
485 +
	//Computation of the coefficients
486 +
		void Interpolation::Compute_Steffen_Coefficients(std::vector<std::vector<double>>& data, std::vector<double> &a,std::vector<double> &b,std::vector<double> &c,std::vector<double> &d)
487 +
		{
488 +
			unsigned int N= data.size();
489 +
			//Compute the Steffen coefficients for the interpolation
490 +
				//1. h and s.
491 +
					double h[N-1],s[N-1];
492 +
					for(unsigned int i=0;i<N-1;i++)
493 +
					{
494 +
						double x_i = data[i][0];
495 +
						double x_ip1 = data[i+1][0];
496 +
						
497 +
						double y_i = data[i][1];
498 +
						double y_ip1 = data[i+1][1];
499 +
						h[i] = x_ip1-x_i;
500 +
						s[i]= (y_ip1-y_i)/h[i];
501 +
					}
502 +
				//2. p and dy
503 +
					double dy[N],p[N];
504 +
					for(unsigned int i=0;i<N;i++)
505 +
					{
506 +
						//First point
507 +
							if(i==0)
508 +
							{
509 +
								p[i] = s[i]*(1.0+h[i]/(h[i]+h[i+1]))-s[i+1]*h[i]/(h[i]+h[i+1]);
510 +
								dy[i] = (Sign(p[i])+Sign(s[i]))*std::min(1.0*abs(s[i]),0.5*abs(p[i]));
511 +
							}
512 +
						//Last point
513 +
							else if(i==N-1)
514 +
							{
515 +
								p[i] = s[i-1]*(1.0+h[i-1]/(h[i-1]+h[i-2]))-s[i-2]*h[i-1]/(h[i-1]+h[i-2]);
516 +
								dy[i] = (Sign(p[i])+Sign(s[i-1]))*std::min(1.0*abs(s[i-1]),0.5*abs(p[i]));
517 +
							}
518 +
						//Points in the middle
519 +
							else
520 +
							{
521 +
								p[i] = (s[i-1]*h[i]+s[i]*h[i-1])/(h[i-1]+h[i]);
522 +
								dy[i] = (Sign(s[i-1])+Sign(s[i]))*std::min(1.0*abs(p[i])/2.0 ,std::min(1.0*abs(s[i]) ,1.0*abs(s[i-1]) ) );
523 +
							}
524 +
					}
525 +
				//3. a,b,c, and d
526 +
					for(unsigned int i=0;i<N-1;i++)
527 +
					{
528 +
						a.push_back((dy[i]+dy[i+1]-2.0*s[i])/pow(h[i],2.0));
529 +
						b.push_back((3.0*s[i]-2.0*dy[i]-dy[i+1])/h[i]);
530 +
						c.push_back(dy[i]);
531 +
						d.push_back(data[i][1]);
532 +
						if(std::isnan(a.back())||std::isnan(b.back())||std::isnan(c.back())||std::isnan(d.back())) 
533 +
						{
534 +
							std::cout <<"Warning: Steffen coefficients in interpolation are NAN."<<std::endl;
535 +
						}
536 +
					}
537 +
		}
538 +
		unsigned int Interpolation::Bisection(double x,int jLeft,int jRight)
539 +
		{
540 +
			while((jRight-jLeft)>1)
541 +
			{
542 +
				int jm = (jRight+jLeft) >>1 ;
543 +
				if (x >= TabulatedData[jm][0])	jLeft=jm;
544 +
				else jRight=jm;
545 +
			}
546 +
			return jLeft;
547 +
		}
548 +
		unsigned int Interpolation::Hunt(double x)
549 +
		{
550 +
551 +
			// 1. Hunting phase returns jd,ju which bracket j
552 +
				int dj = 1;
553 +
				int jd;
554 +
				unsigned int ju;
555 +
				//Hunt up
556 +
				if(x>TabulatedData[jLast][0])
557 +
				{
558 +
					jd = jLast;
559 +
					ju = jd+dj;
560 +
					while(x>TabulatedData[ju][0])
561 +
					{
562 +
						jd = ju;
563 +
						ju += dj;
564 +
						//Check if we ran off the range:
565 +
						if(ju > N_Data-1)
566 +
						{
567 +
							ju = N_Data-1;
568 +
							break;
569 +
						}
570 +
						else
571 +
						{
572 +
							dj += dj;
573 +
						}
574 +
					}
575 +
				}
576 +
				//Hunt down
577 +
				else if (x<TabulatedData[jLast][0])
578 +
				{
579 +
					ju = jLast;
580 +
					jd = ju-dj;
581 +
					while(x<TabulatedData[jd][0])
582 +
					{
583 +
						ju = jd;
584 +
						jd -= dj;
585 +
						//Check if we ran off the range:
586 +
						if(jd < 0)
587 +
						{
588 +
							jd = 0;
589 +
							break;
590 +
						}
591 +
						else
592 +
						{
593 +
							dj += dj;
594 +
						}
595 +
					}
596 +
				}
597 +
				else
598 +
				{
599 +
					return jLast;
600 +
				}
601 +
			// 2. Bisection phase
602 +
				if((ju-jd)>1) jd = Bisection(x,jd,ju);
603 +
			return jd;
604 +
		}
605 +
		// Find j such that list[j]<x<list[j+1]
606 +
		unsigned int Interpolation::Locate(double x)
607 +
		{
608 +
			if( ((xDomain[0]-x)>0.0) || ((x-xDomain[1])>0.0) )
609 +
			{
610 +
				printf("\nError in Interpolation::Locate(): x = %e lies outside the domain [%e,%e].\n\n",x,xDomain[0],xDomain[1]);
611 +
				std::exit(EXIT_FAILURE);
612 +
			}
613 +
			else
614 +
			{
615 +
				//Use Bisection() or the Hunt method, depending of the last calls were correlated.
616 +
				unsigned int j = corr ? Hunt(x): Bisection(x,0,N_Data-1);
617 +
				//Check if the points are still correlated.
618 +
				corr=(fabs((j-jLast))<10);
619 +
620 +
				jLast = j;
621 +
				return j;
622 +
			}
623 +
		}
624 +
	//Constructors
625 +
		Interpolation::Interpolation()
626 +
		{
627 +
			
628 +
			//Generate some data to interpolate zero
629 +
				std::vector<std::vector<double>> data;
630 +
				data.push_back(std::vector<double> {-1.0,0.0});
631 +
				data.push_back(std::vector<double> {0.0,0.0});
632 +
				data.push_back(std::vector<double> {+1.0,0.0});
633 +
			//Define members
634 +
				preFactor = 1.0;
635 +
				TabulatedData = data;
636 +
				N_Data = data.size();
637 +
				xDomain = {-1.0,1.0};
638 +
				jLast = 0;
639 +
				corr=false;
640 +
			//Compute coefficients
641 +
				Compute_Steffen_Coefficients(data,a,b,c,d);
642 +
643 +
		}
644 +
		Interpolation::Interpolation(const std::string& filename,double dim1,double dim2)
645 +
		{
646 +
			//1. Import data.
647 +
				std::ifstream inputfile;
648 +
				inputfile.open(filename);
649 +
				if (inputfile.good())
650 +
				{
651 +
			        while (!inputfile.eof())
652 +
			        {
653 +
			        	double x,y;
654 +
			            inputfile >>x;
655 +
			            inputfile >>y;
656 +
			            TabulatedData.push_back(std::vector<double> {x*dim1,y*dim2});
657 +
			        }
658 +
			        // Close the file.
659 +
			        inputfile.close();
660 +
		   		}
661 +
		   		else
662 +
		   		{
663 +
		        	std::cerr << "Error in Interpolation(" <<filename<<"): File does not exist."<<std::endl;
664 +
		    		std::exit(EXIT_FAILURE);
665 +
		    	}
666 +
			//2. Sort data.
667 +
				std::sort(TabulatedData.begin(), TabulatedData.end());
668 +
			//3. Interpolate data.
669 +
				preFactor = 1.0;
670 +
				N_Data = TabulatedData.size();
671 +
				jLast = 0;
672 +
				corr=false;
673 +
				//3.1 Find the function's domain:
674 +
					std::vector <double> x;
675 +
					for(unsigned i = 0; i < TabulatedData.size(); i++) x.push_back(TabulatedData[i][0]);
676 +
					xDomain.push_back( *min_element(x.begin(),x.end()));
677 +
					xDomain.push_back( *max_element(x.begin(),x.end()));
678 +
				//3.2 Compute the Steffen coefficients for the interpolation
679 +
					Compute_Steffen_Coefficients(TabulatedData,a,b,c,d);
680 +
681 +
		}
682 +
		Interpolation::Interpolation(const std::vector<std::vector<double>>& data,double dim1,double dim2)
683 +
		{
684 +
			preFactor = 1.0;
685 +
			for(unsigned int i =0;i<data.size();i++)
686 +
			{
687 +
				std::vector<double> aux ={data[i][0]*dim1,data[i][1]*dim2};
688 +
				TabulatedData.push_back(aux);
689 +
			}
690 +
			// TabulatedData = data;
691 +
			N_Data = TabulatedData.size();
692 +
			jLast = 0;
693 +
			corr=false;
694 +
			//Find the function's domain:
695 +
				std::vector <double> x;
696 +
				for(unsigned i = 0; i < TabulatedData.size(); i++) x.push_back(TabulatedData[i][0]);
697 +
				xDomain.push_back( *min_element(x.begin(),x.end()));
698 +
				xDomain.push_back( *max_element(x.begin(),x.end()));
699 +
			//Compute the Steffen coefficients for the interpolation
700 +
				Compute_Steffen_Coefficients(TabulatedData,a,b,c,d);
701 +
		}
702 +
	//Return values
703 +
		std::vector<std::vector<double>> Interpolation::Return_Data()
704 +
		{
705 +
			return TabulatedData;
706 +
		}
707 +
		std::vector<double> Interpolation::Return_Domain()
708 +
		{
709 +
			return xDomain;
710 +
		}
711 +
		double Interpolation::Return_Prefactor()
712 +
		{
713 +
			return preFactor;
714 +
		}
715 +
		std::vector<std::vector<double>> Interpolation::Return_Coefficients()
716 +
		{
717 +
			std::vector<std::vector<double>> output;
718 +
			for (unsigned i = 0 ; i < a.size() ; i++)
719 +
			{
720 +
				std::vector<double> aux {a[i],b[i],c[i],d[i]};
721 +
				output.push_back(aux);
722 +
			}
723 +
			return output;
724 +
		}
725 +
	//Set values
726 +
		void Interpolation::Set_Prefactor(double factor)
727 +
		{
728 +
			preFactor=factor;
729 +
		}
730 +
	//Interpolation
731 +
		double Interpolation::Interpolate(double x)
732 +
		{
733 +
			int j = Locate(x);
734 +
			double x_j = TabulatedData[j][0];
735 +
			return preFactor*(a[j]*pow(x-x_j,3.0)+b[j]*pow(x-x_j,2.0)+c[j]*(x-x_j)+d[j]);
736 +
		}
737 +
	//Multiply by a factor
738 +
		void Interpolation::Multiply(double factor)
739 +
		{
740 +
			preFactor*=factor;
741 +
		}
742 +
	//Save function
743 +
	    void Interpolation::Save_Function(std::string filename,unsigned int points)
744 +
	    {
745 +
	    	std::ofstream f;
746 +
	    	f.open(filename);
747 +
	    	double dx = (xDomain[1]-xDomain[0])/(points-1);
748 +
	    	for(unsigned int i = 0;i < points; i++)
749 +
	 		{
750 +
	 			double x = xDomain[0] + i * dx;
751 +
	 			f <<x <<"\t" <<Interpolate(x)<<std::endl;
752 +
	 		}
753 +
	    	f.close();
754 +
	    }
755 +
756 +
//5. Struct for weighted data points
757 +
	//Constructors
758 +
	    DataPoint::DataPoint(double v, double w)
759 +
	    {
760 +
	    	value=v;
761 +
	    	weight=w;
762 +
	    }
763 +
	    DataPoint::DataPoint(double v)
764 +
	    {
765 +
	    	value=v;
766 +
	    	weight=1.0;
767 +
	    }
768 +
	    DataPoint::DataPoint()
769 +
	    {
770 +
	    	value=0.0;
771 +
	    	weight=1.0;
772 +
	    }
773 +
	//Operators
774 +
	    bool operator <(const DataPoint& lhs,const DataPoint& rhs)
775 +
	    {
776 +
	    	return (lhs.value<rhs.value);
777 +
	    }
778 +
		bool operator >(const DataPoint& lhs,const DataPoint& rhs)
779 +
		{
780 +
			return (lhs.value>rhs.value);
781 +
		}
782 +
		bool operator ==(const DataPoint& lhs,const DataPoint& rhs)
783 +
		{
784 +
			return (lhs.value==rhs.value);
785 +
		}
786 +
		std::ostream& operator <<(std::ostream &output,const DataPoint& dp)
787 +
		{
788 +
			output <<dp.value<<"\t" <<dp.weight;
789 +
			return output;
790 +
		}
791 +
	//Weighted mean
792 +
		std::vector<double> Weighted_Average(std::vector<DataPoint>& data)
793 +
		{
794 +
			double sum =0.0;
795 +
			double wsum=0.0;
796 +
			long unsigned int N=data.size();
797 +
			//1. Average
798 +
				for(unsigned int i=0;i<N;i++)
799 +
				{
800 +
					sum+=data[i].weight*data[i].value;
801 +
					wsum+=data[i].weight;
802 +
				}
803 +
				double Average = sum/wsum;
804 +
				double wAverage= wsum/N;
805 +
			//2. Standard error (Cochran)
806 +
				double sum1=0.0,sum2=0.0,sum3=0.0;
807 +
				for(unsigned int i=0;i<data.size();i++)
808 +
				{
809 +
					sum1+=pow((data[i].weight*data[i].value-Average*wAverage),2.0);
810 +
					sum2+=(data[i].weight-wAverage)*(data[i].weight*data[i].value-Average*wAverage);
811 +
					sum3+=pow(data[i].weight-wAverage,2.0);
812 +
				}
813 +
				double SE=N/(N-1.0)/wsum/wsum*(sum1-2.0*Average*sum2+pow(Average,2.0)*sum3);
814 +
			//3. Return result
815 +
				return std::vector<double> {Average,sqrt(SE)};
816 +
		}
817 +
818 +
//6. Kernel density estimation
819 +
	double Gaussian_Kernel(double x)
820 +
	{
821 +
		return 1.0/sqrt(2.0*M_PI)*exp(-x*x/2.0);
822 +
	}
823 +
	Interpolation Perform_KDE(std::vector<DataPoint> data,double xMin,double xMax,double bw)
824 +
	{
825 +
		//Data points
826 +
			unsigned int N_Data = data.size();
827 +
			double Weight_Sum=0.0;
828 +
		//Count the weights
829 +
			for(unsigned int i = 0 ; i<N_Data ; i++) Weight_Sum += data[i].weight;
830 +
831 +
		//1. Bandwidth selection:
832 +
		//If the bandwidth is not set manually, we find it here.
833 +
			if(bw==0)
834 +
			{
835 +
				//1.1 Average.
836 +
					double AverageSum=0.0;
837 +
					for(unsigned int i = 0 ; i<N_Data ; i++)
838 +
					{
839 +
						AverageSum+=data[i].weight*data[i].value;
840 +
					}
841 +
					double Average = AverageSum/Weight_Sum;
842 +
				//1.2 Standard deviation
843 +
					double Variance=0.0;
844 +
					for(unsigned int i = 0 ; i<N_Data ; i++)
845 +
					{
846 +
						Variance += data[i].weight * pow(data[i].value - Average,2.0) / Weight_Sum;
847 +
					}
848 +
				//1.3 Bandwidth with rule-of-thumb estimator
849 +
					bw = sqrt(Variance)*pow(4.0/3.0/N_Data,0.2);
850 +
			}
851 +
			//Sort data:
852 +
				std::sort(data.begin(), data.end());
853 +
			// Pseudo data
854 +
				unsigned int N_PseudoData = N_Data/3.0;
855 +
		//2. Perform and tabulate the KDE
856 +
			int points = 100;
857 +
			double dx = (xMax-xMin)/(points-1);
858 +
			std::vector<std::vector<double>> Interpol_List;
859 +
			for(int j =0;j<points;j++)
860 +
			{
861 +
				double x = xMin + j*dx;
862 +
				double kde=0.0;
863 +
				for(unsigned int i = 0;i<N_Data;i++)
864 +
				{
865 +
					//Vanilla Gauss KDE
866 +
						kde+= data[i].weight*Gaussian_Kernel((x-data[i].value)/bw);
867 +
					//Cowling and Hall pseudo data method
868 +
						if(i<N_PseudoData)
869 +
						{
870 +
							double xPseudo = 4.0*xMin-6.0*data[i].value+4.0*data[2*i].value-data[3*i].value;
871 +
							// double wPseudo = (6.0*data[i].weight+4.0*data[2*i].weight+data[3*i].weight)/9.0;
872 +
							double wPseudo = (data[i].weight+data[2*i].weight+data[3*i].weight)/3.0;
873 +
							kde+= wPseudo*Gaussian_Kernel((x-xPseudo)/bw);
874 +
						}
875 +
					//Reflection method
876 +
						// double xRefl = 2.0*xMin-data[i];
877 +
						// kde+= Gaussian_Kernel((x-data[i])/bw)+Gaussian_Kernel((x-xRefl)/bw);
878 +
						// kde+=weights[i]*(Gaussian_Kernel((x-data[i])/bw)+Gaussian_Kernel((x-xRefl)/bw));
879 +
				}				
880 +
				kde/= bw*Weight_Sum;
881 +
				Interpol_List.push_back( std::vector<double> {x,kde});
882 +
			}
883 +
		//3. Interpolate the list.
884 +
			Interpolation result(Interpol_List);
885 +
		//4. Check normalization/ re-normalize.
886 +
			double norm = Integrate(result,xMin,xMax,0.0001);
887 +
			result.Multiply(1.0/norm);
888 +
		//5. Output
889 +
			return result;
890 +
	}
891 +
892 +
//7. Read in table into vector
893 +
	std::vector<double> Read_List(std::string filepath,double dimension)
894 +
	{
895 +
		std::vector<double> list;
896 +
		std::ifstream inputfile;
897 +
		inputfile.open(filepath);
898 +
		if (inputfile.good())
899 +
		{
900 +
901 +
	        while (!inputfile.eof())
902 +
	        {
903 +
	        	double x;
904 +
	            inputfile >>x;
905 +
	            list.push_back(x*dimension);
906 +
	        }
907 +
	        // Close the file.
908 +
	        inputfile.close();
909 +
   		}
910 +
   		else
911 +
   		{
912 +
        	std::cerr << "Error in Read_List(" <<filepath<<"): File does not exist."<<std::endl;
913 +
        	std::exit(EXIT_FAILURE);
914 +
    	}
915 +
916 +
		return list;
917 +
	}

@@ -0,0 +1,89 @@
Loading
1 +
#include "Physics_Functions.hpp"
2 +
3 +
#include <cmath>
4 +
#include <iostream>
5 +
#include <functional>
6 +
#include <cstdlib>
7 +
8 +
#include "Numerics_Functions.hpp"
9 +
10 +
//1. Units:
11 +
	//Energy
12 +
	const double GeV=1.0;
13 +
	const double eV=1.0E-9*GeV;
14 +
	const double keV=1.0E-6*GeV;
15 +
	const double MeV=1.0E-3*GeV;
16 +
	const double TeV=1.0E3*GeV;
17 +
	const double erg=gram*pow(cm/sec,2);
18 +
	const double Joule=kg*pow(meter/sec,2);
19 +
	//Mass
20 +
	const double gram=5.617977528089887E23*GeV;
21 +
	const double kg=1e3*gram;
22 +
	//Length
23 +
	const double cm=5.068e13/GeV;
24 +
	const double meter=100*cm;
25 +
	const double km=1000*meter;
26 +
	const double fm=1e-15*meter;
27 +
	const double pb=1e-36*pow(cm,2);
28 +
	const double parsec=3.0857e16*meter;
29 +
	const double kpc=1e3*parsec;
30 +
	const double Mpc=1e6*parsec;
31 +
	const double a0 = 5.29177e-11*meter; //Bohr radius
32 +
	//Time
33 +
	const double sec=299792458*meter;
34 +
	const double minute=60*sec;
35 +
	const double hour=60*minute;
36 +
	const double day=24*hour;
37 +
	const double year=365.24*day;
38 +
	//Temperature
39 +
	const double Kelvin=8.62E-14*GeV;
40 +
	//Angle
41 +
	const double deg=M_PI/180.0;
42 +
	
43 +
44 +
//2. Specific Parameters:
45 +
	//Masses
46 +
	const double mPlanck= 1.2209E19*GeV;
47 +
	const double mProton=0.938*GeV ;
48 +
	const double mElectron= 0.511*MeV;
49 +
	const double mNucleon = 0.932*GeV;
50 +
	//Coupling constants
51 +
	const double aEM = 1.0/137.035999139;
52 +
	const double GNewton=pow(mPlanck,-2);
53 +
	const double GFermi = 1.16637e-5/GeV/GeV;
54 +
	//Geographic Parameters
55 +
	const double mEarth=5.972E24*kg;
56 +
	const double rEarth=6371*km;
57 +
	const double rhoEarth=5.51*gram*pow(cm,-3);
58 +
	const double rhoCrust=2.7*gram*pow(cm,-3);
59 +
	//Solar Parameters
60 +
	const double mSun=1.989E30*kg;
61 +
	const double rSun=6.957E8*meter;
62 +
	//Dark Matter Halo Parameters (are defined in the cfg file)
63 +
	double v0=0.0;
64 +
	double vesc=0.0;
65 +
	double vEarth =0.0;
66 +
	double rhoDM=0.0;
67 +
	double Nesc=0.0;
68 +
	//Dark matter form factor reference momentum transfer
69 +
	const double qRef = aEM*mElectron;
70 +
//3. Unit Conversion
71 +
	double InUnits(double quantity, double dimension)
72 +
	{
73 +
		return quantity/dimension;
74 +
	}
75 +
76 +
//5. Simple Physics functions
77 +
	double Reduced_Mass(double m1,double m2)
78 +
	{
79 +
		return m1*m2/(m1+m2);
80 +
	}
81 +
	double NucleusMass(double A)
82 +
	{
83 +
		return A*mNucleon;
84 +
	}
85 +
	double Thomas_Fermi_Radius(int Z)
86 +
	{
87 +
		return pow(9*M_PI*M_PI/2.0/Z,1.0/3.0)/4.0*a0;
88 +
	}
89 +
			

@@ -0,0 +1,216 @@
Loading
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 +
		double FormFactor_N(double q,double A,bool ldm)
13 +
		{
14 +
			if(ldm || q==0.0) return 1.0;
15 +
			else
16 +
			{
17 +
				double a = 0.52*fm;
18 +
				double c = (1.23*pow(A,1.0/3.0)-0.6)*fm;
19 +
				double s = 0.9*fm;
20 +
				double rn = sqrt(c*c+7.0/3.0*pow(M_PI*a,2.0)-5*s*s);
21 +
				double qr = q*rn;
22 +
				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 +
		DM_Particle::DM_Particle()
28 +
		{
29 +
			mass=0.0;
30 +
			sigma_n=0.0;
31 +
			sigma_e = 0.0;
32 +
			ldm = true;
33 +
			formfactor=" ";
34 +
			ZorA = " ";
35 +
			screening=false;
36 +
			mMediator=0.0;
37 +
		}
38 +
39 +
		DM_Particle::DM_Particle(double mDM,double sn,double se,bool light,std::string ff,std::string za,bool scr,double mMed)
40 +
		{
41 +
			mass=mDM;
42 +
			sigma_n=sn;
43 +
			sigma_e = se;
44 +
			ldm = light;
45 +
			formfactor=ff;
46 +
			ZorA = za;
47 +
			screening=scr;
48 +
			mMediator=mMed;
49 +
		}
50 +
		//Atomic form factor (charge screening)
51 +
		double DM_Particle::FormFactor_A(double q,int Z) const
52 +
		{
53 +
			if(!screening) return 1.0;
54 +
			else
55 +
			{
56 +
				double a2q2 = pow(q*Thomas_Fermi_Radius(Z),2.0);
57 +
				return a2q2/(1.0+a2q2);
58 +
			}
59 +
		}
60 +
		void DM_Particle::Set_Mass(double m)
61 +
		{
62 +
			mass = m;
63 +
		}
64 +
		void DM_Particle::Set_Sigma_n(double s)
65 +
		{
66 +
			sigma_n = s;
67 +
			sigma_e = pow(Reduced_Mass(mass,mElectron)/Reduced_Mass(mass,mProton),2.0)*sigma_n;
68 +
		}
69 +
		void DM_Particle::Set_Sigma_e(double s)
70 +
		{
71 +
			sigma_e = s;
72 +
			sigma_n = pow(Reduced_Mass(mass,mProton)/Reduced_Mass(mass,mElectron),2.0)*sigma_e;
73 +
		}
74 +
	//DM form factor
75 +
		double DM_Particle::FormFactor(double q) const 
76 +
		{
77 +
			//Contact interactions
78 +
				if(formfactor=="Contact") return 1.0;
79 +
			//General dark photon
80 +
				else if(formfactor=="General")	return (qRef*qRef+mMediator*mMediator)/(q*q+mMediator*mMediator);
81 +
			//Long range interaction
82 +
				else if (formfactor=="Long-Range")	return qRef*qRef/q/q;
83 +
			//Electric dipole interaction
84 +
				else if (formfactor=="Electric-Dipole")		return qRef/q;
85 +
			//Error
86 +
				else
87 +
				{
88 +
					std::cerr <<"Error in FormFactor(): Form factor "<<formfactor <<"not recognized."<<endl;
89 +
					std::exit(EXIT_FAILURE);
90 +
				}
91 +
		}
92 +
	//Zero momentum transfer spin-independent cross-section
93 +
		double DM_Particle::sigmaSI(int Z,double A) const
94 +
		{
95 +
			double X=-1.0;
96 +
			double mAux=-1.0;
97 +
			if(ZorA=="Z") 
98 +
			{
99 +
				X=Z;
100 +
				mAux = mProton;
101 +
			}
102 +
			else if (ZorA=="A") 
103 +
			{
104 +
				X=A;
105 +
				mAux = mNucleon;
106 +
			}
107 +
			else 
108 +
			{
109 +
				std::cerr <<"Error in sigmaSI: ZorA = " <<ZorA <<" not recognized."<<endl;
110 +
				std::exit(EXIT_FAILURE);
111 +
			}
112 +
			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 +
		double DM_Particle::dSdER(double ER,int Z,double A,double vDM) const
116 +
		{
117 +
			double mA = A*mNucleon;
118 +
			double ERmax = 2.0*pow(Reduced_Mass(mass,mA)*vDM,2.0)/mA;
119 +
			double q = sqrt(2.0*mA*ER);	
120 +
			if(q==0&&(formfactor=="Electric-Dipole"||formfactor=="Long-Range")&&screening) q=1e-15*keV; //Screening cause 1/q to cancel
121 +
			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 +
		double DM_Particle::dSdq2(double q,int Z,double A,double vDM) const
124 +
		{
125 +
			double mA = A*mNucleon;
126 +
			double qMax2 = 4.0*pow(Reduced_Mass(mass,mA)*vDM,2.0);
127 +
			if(q==0&&(formfactor=="Electric-Dipole"||formfactor=="Long-Range")&&screening) q=1e-15*keV; //Screening cause 1/q to cancel
128 +
			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 +
		double DM_Particle::Sigma_Tot(int Z,double A,double vDM) const
131 +
		{
132 +
			if(!ldm)
133 +
			{
134 +
				double mA = A*mNucleon;
135 +
				double ERmax = 2.0*pow(Reduced_Mass(mass,mA)*vDM,2.0)/mA;
136 +
				double ERmin = 0.0;
137 +
				//integrate the diff cross section
138 +
					auto dodER = std::bind(&DM_Particle::dSdER,this,std::placeholders::_1,Z,A,vDM);
139 +
					//Numerical integration
140 +
			  			double epsilon = Find_Epsilon(dodER,ERmin,ERmax,1e-5);
141 +
			  			double integral =Integrate(dodER,ERmin,ERmax,epsilon);
142 +
					return integral;
143 +
			}
144 +
			else
145 +
			{
146 +
				double result = sigmaSI(Z,A);
147 +
				
148 +
				double q2max = pow(2.0*Reduced_Mass(mass,A*mNucleon)*vDM,2.0);
149 +
				double a2=pow(Thomas_Fermi_Radius(Z),2.0);
150 +
				//General interaction
151 +
					if(formfactor == "General")
152 +
					{
153 +
						double m2 = pow(mMediator,2.0);
154 +
						if(screening)
155 +
						{
156 +
							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 +
							result*= pow(qRef*qRef+m2,2.0)/m2/(m2+q2max);
161 +
						}
162 +
					}
163 +
				//Contact interaction
164 +
					else if(formfactor == "Contact" && screening)
165 +
					{
166 +
						result *=(1.0+1.0/(1+a2*q2max)-2.0/a2/q2max*log1p(a2*q2max));
167 +
					}
168 +
				//Electric dipole interaction
169 +
					else if(formfactor == "Electric-Dipole")
170 +
					{
171 +
						result *= pow(qRef,2.0)/q2max*(log1p(a2*q2max)-a2*q2max/(1+a2*q2max));
172 +
					}
173 +
				//Long range interaction
174 +
					else if(formfactor== "Long-Range")
175 +
					{
176 +
						result *= pow(a2,2.0)*pow(qRef,4.0)/(1+a2*q2max);
177 +
					}
178 +
179 +
				return result;
180 +
			}
181 +
		}
182 +
183 +
	//Find the stopping cross section
184 +
		double Stopping_CrossSection(const std::string& detector,double sigma,double mDM)
185 +
		{
186 +
			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 +
			else return sigma;
188 +
		}
189 +
190 +
//7. DM speed distribution
191 +
	double SpeedDistribution(double v,double vEarth)
192 +
	{
193 +
		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 +
 	double Average_Speed(double vEarth,double vMin)
196 +
 	{
197 +
 		//1. integrand
198 +
	 		std::function<double(double)> integrand = [vEarth] (double v)
199 +
	 		{
200 +
	 			return v*SpeedDistribution(v,vEarth);
201 +
	 		};
202 +
	 	//2. integrate
203 +
	 		double vMean = Integrate(integrand,vMin,(vesc+vEarth),1e-6);
204 +
	 	//3. Renormalize
205 +
	 		if(vMin>0.0)
206 +
	 		{
207 +
	 			
208 +
	 			std::function<double(double)> integrand2 = [vEarth] (double v)
209 +
		 		{
210 +
		 			return SpeedDistribution(v,vEarth);
211 +
		 		};
212 +
		 		double norm=Integrate(integrand2,vMin,vesc+vEarth,1e-6);
213 +
		 		vMean/=norm;
214 +
	 		}
215 +
	 		return vMean;
216 +
 	}

@@ -0,0 +1,477 @@
Loading
1 +
#include "DD_Electron.hpp"
2 +
3 +
#include <iostream>
4 +
#include <fstream>
5 +
#include <cmath>
6 +
#include <cstdlib>
7 +
8 +
#include "Physics_Functions.hpp"
9 +
#include "Input_Parameters.hpp"
10 +
#include "DD_General.hpp"
11 +
12 +
//1. Import crystal/ionization form factor
13 +
	//Crystal form factor
14 +
	double FFcrystalSi[900][500];
15 +
	double FFcrystalGe[900][500];
16 +
17 +
	//Ionization form factors for different atomic shells
18 +
	std::vector<Shell> Shells;
19 +
20 +
	//Import tables
21 +
	void Import_FormFactor(const std::string& target)
22 +
	{
23 +
		//a) Semiconductor
24 +
			if(target=="Si"||target=="Ge")
25 +
			{
26 +
				//Prefactors:
27 +
				double dE = 0.1*eV;
28 +
				double wk = 2.0/137.0;
29 +
				double prefactor = (target=="Si")? 2.0*eV: 1.8*eV;
30 +
				
31 +
				//Import data:
32 +
				ifstream f;
33 +
				f.open("../detectors/Semiconductors/C."+target+"137.dat");
34 +
				for(int Ei=1;Ei<=500;Ei++)
35 +
				{
36 +
					for(int qi=1;qi<=900;qi++)
37 +
					{
38 +
						if(target=="Si")
39 +
						{
40 +
							f >> FFcrystalSi[qi-1][Ei-1];
41 +
							FFcrystalSi[qi-1][Ei-1]*=prefactor*qi/dE*wk/4.0;
42 +
						}
43 +
						else
44 +
						{
45 +
							f >> FFcrystalGe[qi-1][Ei-1]; 
46 +
							FFcrystalGe[qi-1][Ei-1]*=prefactor*qi/dE*wk/4.0;
47 +
48 +
						}
49 +
					}
50 +
				}
51 +
				f.close();
52 +
			}
53 +
		//b) Xenon
54 +
			else if(target=="Xe")
55 +
			{
56 +
				Shells.push_back( Shell ("5p","../detectors/FF_Xenon/FF_5p.txt",12.4*eV,0,-2.4,2.4,97,-1,4,101) );
57 +
				Shells.push_back( Shell ("5s","../detectors/FF_Xenon/FF_5s.txt",25.7*eV,0,-2.4,1.8,43,-1,4,51) );
58 +
				Shells.push_back( Shell ("4d","../detectors/FF_Xenon/FF_4d.txt",75.6*eV,4,-2.4,1.75,84,-1,4,101) );
59 +
				Shells.push_back( Shell ("4p","../detectors/FF_Xenon/FF_4p.txt",163.5*eV,6,-2.4,1.65,36,-1,4,101) );
60 +
				Shells.push_back( Shell ("4s","../detectors/FF_Xenon/FF_4s.txt",213.8*eV,3,-2.4,1.8,43,-1,4,51) );
61 +
			}
62 +
		//c) Argon
63 +
			else if(target=="Ar")
64 +
			{
65 +
				Shells.push_back( Shell ("3p32","../detectors/FF_Argon/FF_3p32.txt",	15.7*eV,0,-2.4,1.8,14,-1,4,51) );
66 +
				Shells.push_back( Shell ("3p12","../detectors/FF_Argon/FF_3p12.txt",	15.9*eV,0,-2.4,1.8,14,-1,4,51) );
67 +
				Shells.push_back( Shell ("3s"	,"../detectors/FF_Argon/FF_3s.txt",		29.3*eV,0,-2.4,1.8,43,-1,4,51) );
68 +
			}
69 +
		//d) Error
70 +
			else
71 +
			{
72 +
				cerr<<"Error in Import_Formfactor(): Target "<<target<<" not recognized."<<endl;
73 +
				std::exit(EXIT_FAILURE);
74 +
			}
75 +
	}
76 +
77 +
78 +
//2. Semiconductors
79 +
	//Minimal velocity
80 +
	double vMinimal_e(double q,double Ee,double mDM)
81 +
	{
82 +
		return (Ee/q+q/2.0/mDM);
83 +
	}
84 +
85 +
	//Recoil spectra
86 +
	double dRdEe(double Ee,const DM_Particle& DM,const std::string& target,double efficiency)
87 +
	{
88 +
		double dE = 0.1*eV;
89 +
		double dq = 0.02*aEM*mElectron;
90 +
		
91 +
		//Target specific parameters:
92 +
		double Mcell;
93 +
		if(target == "Si") 			Mcell = 2.0*28.08*mNucleon;
94 +
		else if (target == "Ge")	Mcell = 2.0*72.6*mNucleon;
95 +
		else
96 +
		{
97 +
			std::cerr <<"Error in dRdEe(): target " <<target <<" not recognized."<<endl;
98 +
			std::exit(EXIT_FAILURE);
99 +
		}
100 +
		double (&FFcrystal)[900][500] = (target=="Si") ? FFcrystalSi : FFcrystalGe;
101 +
102 +
		double prefactor = efficiency*rhoDM/DM.mass/Mcell*DM.sigma_e*aEM*pow(mElectron/Reduced_Mass(DM.mass,mElectron),2.0);
103 +
		
104 +
		//Integrate by summing over the tabulated form factors
105 +
		int Ei = floor(Ee/dE)-1;
106 +
		double sum=0.0;
107 +
		for(int qi=0;qi<900;qi++) sum+=1.0/(qi+1)/(qi+1)/dq*EtaFunction(vMinimal_e((qi+1)*dq,Ee,DM.mass),vEarth)*pow(DM.FormFactor((qi+1)*dq),2.0)*FFcrystal[qi][Ei];
108 +
		
109 +
		return prefactor*sum;
110 +
	}
111 +
	Interpolation dRdEe(const DM_Particle& DM,const std::string& target,double efficiency,const std::vector<double>& attenuation,const std::vector<DataPoint> &speeddata,double vCutoff)
112 +
	{
113 +
		double dE = 0.1*eV;
114 +
		double dq = 0.02*aEM*mElectron;
115 +
		//Target specific parameters
116 +
		double Mcell,Egap,eps;
117 +
		if(target == "Si")
118 +
		{
119 +
			Mcell = 2.0*28.08*mNucleon;
120 +
			Egap = 1.11*eV;
121 +
			eps = 3.6*eV;
122 +
		}			
123 +
		else if (target == "Ge")	
124 +
		{
125 +
			Mcell = 2.0*72.6*mNucleon;
126 +
			Egap = 0.67*eV;
127 +
			eps = 2.9*eV;
128 +
		}
129 +
		else
130 +
		{
131 +
			std::cerr <<"Error in dRdEe(): Target " <<target <<" not recognized."<<endl;
132 +
			std::exit(EXIT_FAILURE);
133 +
		}
134 +
		double (&FFcrystal)[900][500] = (target=="Si") ? FFcrystalSi : FFcrystalGe;
135 +
136 +
		//MC eta function
137 +
		Interpolation etaMC = EtaFunction(attenuation,speeddata,vCutoff,vEarth);
138 +
139 +
		double prefactor = efficiency*rhoDM/DM.mass/Mcell*DM.sigma_e*aEM*pow(mElectron/Reduced_Mass(DM.mass,mElectron),2.0);
140 +
141 +
		//Interpolation
142 +
		int points =500;
143 +
		std::vector<std::vector<double>> interpol_list;
144 +
		double Ethr= eps*(DMe_threshold-1)+Egap;
145 +
		for(int Ei=Ethr/dE;Ei<points;Ei++)
146 +
		{
147 +
			double sum=0.0;
148 +
			for(int qi=0;qi<900;qi++)
149 +
			{
150 +
				sum+=1.0/(qi+1)/(qi+1)/dq*etaMC(vMinimal_e((qi+1)*dq,(Ei+1)*dE,DM.mass))*pow(DM.FormFactor((qi+1)*dq),2.0)*FFcrystal[qi][Ei];
151 +
			}
152 +
			interpol_list.push_back(std::vector<double> {(Ei)*dE,prefactor*sum});
153 +
		}
154 +
		Interpolation spectrum(interpol_list);
155 +
156 +
		return spectrum;
157 +
	}
158 +
159 +
	//Total event rates
160 +
	double TotalRate(const DM_Particle& DM,int Qthreshold,double efficiency,const std::string& target)
161 +
	{
162 +
		double dE = 0.1*eV;
163 +
		double dq = 0.02*aEM*mElectron;
164 +
165 +
		//Target specific parameters
166 +
		double Mcell,Egap,eps;
167 +
		if(target == "Si")
168 +
		{
169 +
			Mcell = 2.0*28.08*mNucleon;
170 +
			Egap = 1.11*eV;
171 +
			eps = 3.6*eV;
172 +
		}			
173 +
		else if (target == "Ge")	
174 +
		{
175 +
			Mcell = 2.0*72.6*mNucleon;
176 +
			Egap = 0.67*eV;
177 +
			eps = 2.9*eV;
178 +
		}
179 +
		else
180 +
		{
181 +
			std::cerr <<"Error in TotalRate(): Target " <<target <<" not recognized."<<endl;
182 +
			std::exit(EXIT_FAILURE);
183 +
		}
184 +
		double Ethr= eps*(Qthreshold-1)+Egap;
185 +
		double (&FFcrystal)[900][500] = (target=="Si") ? FFcrystalSi : FFcrystalGe;
186 +
187 +
		//Integrating over the form factors by summing the tabulated values.
188 +
		double prefactor =efficiency*rhoDM/DM.mass/Mcell*DM.sigma_e*aEM*pow(mElectron/Reduced_Mass(DM.mass,mElectron),2.0);
189 +
		double sum=0.0;
190 +
		for(int Ei=(Ethr/dE);Ei<500;Ei++)
191 +
		{
192 +
			for(int qi=0;qi<900;qi++) sum+=dE/(qi+1)/(qi+1)/dq*EtaFunction(vMinimal_e((qi+1)*dq,(Ei+1)*dE,DM.mass),vEarth)*pow(DM.FormFactor((qi+1)*dq),2.0)*FFcrystal[qi][Ei];
193 +
		}
194 +
195 +
		return prefactor*sum;
196 +
	}
197 +
	double TotalRate(const DM_Particle& DM,int Qthreshold,double efficiency,const std::string& target,const std::vector<double>& attenuation,const std::vector<DataPoint> &speeddata,double vCutoff)
198 +
	{
199 +
		double dE = 0.1*eV;
200 +
		double dq = 0.02*aEM*mElectron;
201 +
		//Target specific parameters
202 +
		double Mcell,Egap,eps;
203 +
		if(target == "Si")
204 +
		{
205 +
			Mcell = 2.0*28.08*mNucleon;
206 +
			Egap = 1.11*eV;
207 +
			eps = 3.6*eV;
208 +
		}			
209 +
		else if (target == "Ge")	
210 +
		{
211 +
			Mcell = 2.0*72.6*mNucleon;
212 +
			Egap = 0.67*eV;
213 +
			eps = 2.9*eV;
214 +
		}
215 +
		else
216 +
		{
217 +
			std::cerr <<"Error in TotalRate(): Target " <<target <<" not recognized."<<endl;
218 +
			std::exit(EXIT_FAILURE);
219 +
		}
220 +
221 +
		//Energy threshold
222 +
		double Ethr= eps*(Qthreshold-1)+Egap;
223 +
		
224 +
		//Crystal form factor
225 +
		double (&FFcrystal)[900][500] = (target=="Si") ? FFcrystalSi : FFcrystalGe;
226 +
		
227 +
		//MC eta function
228 +
		Interpolation etaMC = EtaFunction(attenuation,speeddata,vCutoff,vEarth);
229 +
		
230 +
		//Integrating over the form factors by summing the tabulated values.
231 +
		double prefactor = efficiency*rhoDM/DM.mass/Mcell*DM.sigma_e*aEM*pow(mElectron/Reduced_Mass(DM.mass,mElectron),2.0);
232 +
		double sum=0.0;
233 +
		for(int Ei=(Ethr/dE);Ei<500;Ei++)
234 +
		{
235 +
			for(int qi=0;qi<900;qi++) sum+=dE/(qi+1)/(qi+1)/dq*etaMC(vMinimal_e((qi+1)*dq,(Ei+1)*dE,DM.mass))*pow(DM.FormFactor((qi+1)*dq),2.0)*FFcrystal[qi][Ei];
236 +
		}
237 +
		
238 +
		return prefactor*sum;
239 +
	}
240 +
241 +
242 +
//3. Atomic shell struct	
243 +
	//Constructor
244 +
	Shell::Shell(std::string nm,std::string filepath,double Eb,int Ns,double lnkMin,double lnkMax,int Nk,double lnqMin,double lnqMax,int Nq)
245 +
	{
246 +
		name = nm;
247 +
		Ebinding = Eb;
248 +
		neSecondary	=Ns;
249 +
		//Ionization form factor
250 +
			//Data range in k and q (logarithmically in units of aEM mElectron)
251 +
			logkMin=lnkMin;
252 +
			logkMax=lnkMax;
253 +
			logqMin=lnqMin;
254 +
			logqMax=lnqMax;
255 +
			nk=Nk;
256 +
			nq=Nq;
257 +
			dlogk=(logkMax-logkMin)/(nk-1);
258 +
			dlogq=(logqMax-logqMin)/(nq-1);
259 +
260 +
			//Import form factor table
261 +
			ifstream inputfile;
262 +
			inputfile.open(filepath);
263 +
			if (inputfile.good())
264 +
			{
265 +
				for(int ki=0;ki<Nk;ki++)