Attachment 'WIMP_functions.C'
Download 1 #ifndef WIMP_FUNCTIONS_h
2 #define WIMP_FUNCTIONS_h
3
4 #include <iostream>
5 #include <cstdio>
6 #include <TMath.h>
7 #include <TSpline.h>
8 #include <TH2.h>
9
10 #ifndef FUNCTIONS_h
11 //#include "functions.h"
12 #endif
13 #include "Functions.C"
14 using namespace std;
15
16 // ##################################################################
17 // DECLARATIONS
18 // ##################################################################
19
20 int LoadWIMPParameters(FILE* f);
21 double GetWIMP_PolFactor(double kappa_eRpL, double kappa_eRpR,double kappa_eLpL,
22 double kappa_eLpR,double Pe, double Pp);
23 double DiffWIMPCrossSection(double ecms, double E, double theta);
24 double CoreDiffWIMPCrossSection(double ecms, double E, double theta);
25 double GetSignalWeight(double ecms, double E, double theta, double XSec_bg);
26 double GetCoreSignalWeight(double ecms, double E, double theta, double XSec_bg);
27 double TotalWIMPXSec(double Emin, double Emax, double PTCut, double coscut);
28 double TotalCoreWIMPXSec(double Emin, double Emax, double PTCut, double coscut);
29
30 // ##################################################################
31 // WIMP PARAMETERS
32 // ##################################################################
33
34 double _mass;
35 double _alpha_em;
36 double _BR;
37 double _sigma_an;
38 double _annihilator;
39 double _spin;
40 double _kappa_eRpL;
41 double _kappa_eRpR;
42 double _kappa_eLpL;
43 double _kappa_eLpR;
44 double _lambda;
45
46
47 double _BR_eff;
48 double _WIMP_PolFactor;
49
50
51
52
53 TSpline3* spline_eLpR;
54 TSpline3* spline_eRpL;
55
56
57 // ##################################################################
58 // DEFINITIONS
59 // ##################################################################
60
61
62 int LoadWIMPParameters(FILE* f){
63
64 char thisParName[16];
65 float thisValue;
66
67 while (fscanf (f, "%s %f ", thisParName, &thisValue) != EOF) {
68 // JUST PRINTOUT
69 // printf ("%s %.3f\n",thisParName,thisValue);
70
71 if (!memcmp (thisParName, "wimpMass", strlen(thisParName)))
72 _mass = thisValue;
73 if (!memcmp (thisParName, "alpha_em", strlen(thisParName)))
74 _alpha_em = thisValue;
75 if (!memcmp (thisParName, "sigma_an", strlen(thisParName)))
76 _sigma_an = thisValue;
77 if (!memcmp (thisParName, "wimpSpin", strlen(thisParName)))
78 _spin = thisValue;
79 if (!memcmp (thisParName, "Annihilator", strlen(thisParName)))
80 _annihilator = thisValue;
81 if (!memcmp (thisParName, "BR_ep", strlen(thisParName)))
82 _BR = thisValue;
83 if (!memcmp (thisParName, "kappa_eRpL", strlen(thisParName)))
84 _kappa_eRpL = thisValue;
85 if (!memcmp (thisParName, "kappa_eRpR", strlen(thisParName)))
86 _kappa_eRpR = thisValue;
87 if (!memcmp (thisParName, "kappa_eLpL", strlen(thisParName)))
88 _kappa_eLpL = thisValue;
89 if (!memcmp (thisParName, "kappa_eLpR", strlen(thisParName)))
90 _kappa_eLpR = thisValue;
91 }
92 return 0;
93 }
94
95
96
97 // #######################################################################
98
99 double GetWIMP_PolFactor(double kappa_eRpL, double kappa_eRpR, double kappa_eLpL, double kappa_eLpR,double Pe, double Pp){
100
101 double Polfactor = (0.25*(1.0+Pe) * ((1.0-Pp)*kappa_eRpL + (1.0+Pp)*kappa_eRpR)+
102 0.25*(1.0-Pe) * ((1.0-Pp)*kappa_eLpL + (1.0+Pp)*kappa_eLpR));
103 // printf("WIMP polarisation factor: %.3f\n",Polfactor);
104 return Polfactor;
105 }
106
107
108 // #################################################################
109
110 double DiffWIMPCrossSection(double ecms, double E, double theta){
111 const double PI = TMath::Pi();
112 double diffXSec(0.0);
113
114 double E_max = ((0.5*ecms)*(0.5*ecms)-(_mass*_mass))/(0.5*ecms);
115 if (E > E_max ){
116 diffXSec = 0.0;
117 return diffXSec;
118 }
119
120
121 double x = 2.0*E/ecms;
122 diffXSec = (_alpha_em * _BR_eff*_sigma_an/(16.0*PI))*
123 ((1.0+(1.0-x)*(1.0-x))/x)*1.0/(TMath::Sin(theta)*TMath::Sin(theta))*
124 (pow(2.0,2.0*_annihilator)*(2.0*_spin+1.0)*(2.0*_spin+1.0))*
125 pow((1.0-(4.0*_mass*_mass/((1.0-x)*ecms*ecms))),(0.5+_annihilator));
126 return diffXSec;
127 }
128
129
130 // #################################################################
131
132 double WIMPdsigmadx(double ecms, double E, double cos_cut){
133 const double PI = TMath::Pi();
134 double diffXSec(0.0);
135
136 double E_max = ((0.5*ecms)*(0.5*ecms)-(_mass*_mass))/(0.5*ecms);
137 if (E > E_max){
138 diffXSec = 0.0;
139 return diffXSec;
140 }
141
142 double theta = TMath::ACos(cos_cut);
143 double theta1 = theta;
144 double theta2 = PI - theta;
145 double integration_theta = TMath::Log(TMath::Tan(theta2/2.0)) - TMath::Log(TMath::Tan(theta1/2.0 ));
146
147
148 double x = 2.0*E/ecms;
149 diffXSec = (_alpha_em * _BR_eff*_sigma_an/(16.0*PI))*
150 ((1.0+(1.0-x)*(1.0-x))/x)*integration_theta*
151 (pow(2.0,2.0*_annihilator)*(2.0*_spin+1.0)*(2.0*_spin+1.0))*
152 pow((1.0-(4.0*_mass*_mass/((1.0-x)*ecms*ecms))),(0.5+_annihilator));
153 return diffXSec;
154 }
155
156 // #################################################################
157
158 TH1D* DsigmadE_WIMP(double ecms, double coscut){
159
160 TH1D* h = new TH1D("DsigmadE_WIMP","DsigmadE_WIMP;E[GeV];d#sigma /dE [fb]",300,0.0,300.0);
161 double E;
162
163 for (int i(1); i < 301; i++){
164 E = h->GetXaxis()->GetBinCenter(i);
165 double xsec = WIMPdsigmadx(ecms, E,coscut);
166 xsec *= 1000*2/500.0; // conversion pb->fb and x->E
167 h->SetBinContent(i,xsec);
168 }
169 return h;
170 }
171
172 // #################################################################
173
174 void DsigmadE_WIMP(TH1D* h, double ecms, double coscut){
175
176 double E;
177
178 for (int i(1); i < 301; i++){
179 E = h->GetXaxis()->GetBinCenter(i);
180 double xsec = WIMPdsigmadx(ecms, E,coscut);
181 xsec *= 1000*2/500.0; // conversion pb->fb and x->E
182 h->SetBinContent(i,xsec);
183 }
184 return;
185 }
186
187
188 // #################################################################
189
190 double CoreDiffWIMPCrossSection(double ecms, double E, double theta){ // has to be multiplied with kappa_eff and (2*S+1)^2
191 const double PI = TMath::Pi();
192 double diffXSec(0.0);
193
194 double E_max = ((0.5*ecms)*(0.5*ecms)-(_mass*_mass))/(0.5*ecms);
195 if (E > E_max){ // Only kinematic cut?
196 diffXSec = 0.0;
197 return diffXSec;
198 }
199
200
201 double x = 2.0*E/ecms;
202 diffXSec = (_alpha_em * _sigma_an/(16.0*PI))*
203 ((1.0+(1.0-x)*(1.0-x))/x)*1.0/(TMath::Sin(theta)*TMath::Sin(theta))*
204 (pow(2.0,2.0*_annihilator))*
205 pow((1.0-(4.0*_mass*_mass/((1.0-x)*ecms*ecms))),(0.5+_annihilator));
206
207 return diffXSec;
208 }
209
210 // #################################################################
211
212 double GetSignalWeight(double ecms, double E, double theta, double XSec_bg){
213
214 double weight(0.0);
215
216 double wimp_xsec = DiffWIMPCrossSection(ecms,E,theta);
217 weight = wimp_xsec/XSec_bg;
218
219
220 return weight;
221 }
222
223 // #################################################################
224
225 double GetCoreSignalWeight(double ecms, double E, double theta, double XSec_bg){
226
227 double weight(0.0);
228
229 double wimp_xsec = CoreDiffWIMPCrossSection(ecms,E,theta);
230 weight = wimp_xsec/XSec_bg;
231
232
233 return weight;
234 }
235
236 // #################################################################
237
238 double TotalWIMPXSec(double Emin, double Emax, double PTCut, double coscut){
239 double XSec(.0);
240
241 double dxs(0.0);
242
243 double x;
244 double cos;
245
246 double E;
247 double thet;
248 double pt;
249
250 TH2D* xsec = new TH2D("xsec","XSec;x;cos(#Theta)",1000,0.0,1.0,1000,-1.0,1.0);
251 for (int i(1); i < 1001; i++){
252 x = xsec->GetXaxis()->GetBinCenter(i);
253 // E = run_CMSEnergy*x/2.0;
254 E = CMSEnergy*x/2.0;
255 if (E < Emin || E > Emax) continue;
256
257 for (int j(1); j < 1001; j++){
258
259 cos = xsec->GetYaxis()->GetBinCenter(j);
260 thet = TMath::ACos(cos);
261 pt = E * TMath::Sin(thet);
262 if (TMath::Abs(cos) > coscut || pt < PTCut) continue;
263
264
265 dxs = DiffWIMPCrossSection(CMSEnergy,E,thet);
266 // dxs = DiffWIMPCrossSection(run_CMSEnergy,E,thet);
267 xsec->SetBinContent(i,j,dxs);
268 }
269 }
270 XSec = xsec->Integral("WIDTH");
271 delete xsec;
272 return XSec;
273 }
274
275
276 // #################################################################
277
278 double TotalCoreWIMPXSec(double Emin, double Emax, double PTCut, double coscut){
279 double XSec(.0);
280
281 double dxs(0.0);
282
283 double x;
284 double cos;
285
286 double E;
287 double thet;
288 double pt;
289
290 TH2D* xsec = new TH2D("xsec","XSec;x;cos(#Theta)",1000,0.0,1.0,1000,-1.0,1.0);
291 for (int i(1); i < 1001; i++){
292 x = xsec->GetXaxis()->GetBinCenter(i);
293 // E = run_CMSEnergy*x/2.0;
294 E = CMSEnergy*x/2.0;
295 if (E < Emin || E > Emax) continue;
296
297 for (int j(1); j < 1001; j++){
298
299 cos = xsec->GetYaxis()->GetBinCenter(j);
300 thet = TMath::ACos(cos);
301 pt = E * TMath::Sin(thet);
302 if (TMath::Abs(cos) > coscut || pt < PTCut) continue;
303
304 dxs = CoreDiffWIMPCrossSection(CMSEnergy,E,thet);
305 // dxs = CoreDiffWIMPCrossSection(run_CMSEnergy,E,thet);
306 xsec->SetBinContent(i,j,dxs);
307 }
308 }
309 XSec = xsec->Integral("WIDTH");
310 delete xsec;
311 return XSec;
312 }
313
314
315 double FitShape(double* e,double* par){
316
317 // const double PI = TMath::Pi();
318 double fitval;
319
320 // par[0]: norm, par[1]: annihilator, par[2]: mass, par[3]: ecms
321
322 double x = e[0]*2.0/par[3];
323 if (( 4.0*par[2]*par[2]/((1.0-x)*par[3]*par[3]) ) > 1.0) {
324 fitval = 0.0;
325 return fitval;
326 }
327
328 fitval = par[0]*
329 //(_alpha_em * _sigma_an/(16.0*PI))*
330 ((1.0+(1.0-x)*(1.0-x))/x)*
331 (pow(2.0,2.0*par[1]))*
332 pow((1.0-(4.0*par[2]*par[2]/((1.0-x)*par[3]*par[3]))),(0.5+par[1]));
333 return fitval;
334 }
335
336 double fitcont(double* e,double* par){
337
338 double fitval;
339 fitval = par[0] * (PolarisationWeights[0] * spline_eLpR->Eval(e[0]) + PolarisationWeights[1] * spline_eRpL->Eval(e[0]));
340 return fitval;
341 }
342
343
344 double fitfunc(double* e,double* par){
345 double fitval;
346 fitval = FitShape(e,par) + fitcont(e,&par[4]);
347 return fitval;
348 }
349
350
351 // Calculation cross secrion using diff operrators
352 // #################################################################
353
354 double WIMPdsigmadx_new(double ecms, double E, double cos_cut){
355 const double PI = TMath::Pi();
356 double diffXSec(0.0);
357
358 double E_max = ((0.5*ecms)*(0.5*ecms)-(_mass*_mass))/(0.5*ecms);
359
360
361 if ( E > E_max ){
362 diffXSec = 0.0;
363 return diffXSec;
364
365 }
366
367 double theta = TMath::ACos(cos_cut);
368 double theta1 = theta;
369 double theta2 = PI - theta;
370 double integration_theta2 = TMath::Cos(theta1)-TMath::Cos(theta2);
371 double integration_theta1 = TMath::Log(TMath::Tan(theta2/2.0)) - TMath::Log(TMath::Tan(theta1/2.0 ));
372 // if(_mass==100)
373 // cout<<"Integral Theta1 = "<<integration_theta1<<"Integral Theta2 = "<<integration_theta2<<endl;
374
375 double z = 2.0*E/ecms;
376 double mu=_mass/ecms;
377
378 /* //Vector Integration theta
379 diffXSec = (_alpha_em/(12.0*PI*PI))*(ecms*ecms/pow(_lambda,4))*
380 (1/z)*sqrt((1-z-4*mu*mu)/(1-z))*
381 (1-z+2*mu*mu)*
382 (4*(1-z)*integration_theta1+z*z*
383 (2*integration_theta1-integration_theta2));
384 */
385
386 //Axial-vector
387 /* diffXSec = (_alpha_em/(12.0*PI*PI))*(ecms*ecms/pow(_lambda,4))*
388 (1/z)*pow(((1-z-4*mu*mu)/(1-z)),1.5)*(1-z)*
389 (4*(1-z)*integration_theta1+z*z*
390 (2*integration_theta1-integration_theta2));
391
392 */
393
394 //Scalar, s-chanel
395 diffXSec = (_alpha_em/(8.0*PI*PI))*(ecms*ecms/pow(_lambda,4))*
396 (1/z)*integration_theta1*
397 pow(((1-z-4*mu*mu)/(1-z)),1.5)*
398 (1-z)*(2*(1-z)+z*z);
399
400
401
402
403 return diffXSec;
404 }
405
406 // #################################################################
407
408 void DsigmadE_WIMP_new(TH1D* h, double ecms, double coscut){
409
410 double E,counter;
411
412 for (int i(1); i < 301; i++){
413 E = h->GetXaxis()->GetBinCenter(i);
414 double xsec = WIMPdsigmadx_new(ecms, E,coscut);
415 xsec *= 2/500.0; // conversion z->E,
416 xsec *= 0.38938*pow(10.0,12.0); // conversion 1/Gev^2 to fb;
417 if(_mass==1&&i>0&&i<211){
418 cout<<"Cross section = "<<xsec<<" Value mass = "<<_mass<<" Spin = "<<_spin<<" Lambda = "<<_lambda<<" Photon Energy = "<<E<<endl;
419 counter+=xsec;
420 // if(i==210)
421 cout<<"Counter = "<<counter<<" E = "<<E<<endl;
422 }
423
424 h->SetBinContent(i,xsec);
425 }
426 return;
427 }
428
429 // #######################################################################
430
431 double Sigma_PolFactor(double sigma_eRpL, double sigma_eRpR, double sigma_eLpL, double sigma_eLpR,double Pe, double Pp){
432
433 double Polfactor = (0.25*(1.0+Pe) * ((1.0+Pp)*sigma_eRpR + (1.0-Pp)*sigma_eRpL)+
434 0.25*(1.0-Pe) * ((1.0+Pp)*sigma_eLpR + (1.0-Pp)*sigma_eLpL));
435 // printf("WIMP polarisation factor: %.3f\n",Polfactor);
436 return Polfactor;
437 }
438
439
440 //double FitShape1(double* e,double* par){
441 //
442 // const double PI = TMath::Pi();
443 // double fitval;
444 //
445 // // par[0]: norm, par[1]: annihilator, par[2]: mass, par[3]: ecms
446 //
447 // char str[256];
448 //
449 // sprintf(str,"((1.0+(1.0-x)*(1.0-x))/x)* %.10f * TMath::Power((1.0-(4.0*%.10f/((1.0-x)*%.10f))),(0.5+%.10f))",(pow(2.0,2.0*par[1])),par[2]*par[2],par[3]*par[3],par[1]);
450 //
451 // TF1* f = new TF1("f1",str);
452 //
453 // double integral = f->Integral(par[5]*2.0/par[3],par[6]*2.0/par[3]);
454 // cout << "INT: "<< integral << endl;
455 // par[0] = par[4]/integral;
456 //
457 // double x = e[0]*2.0/par[3];
458 // if (( 4.0*par[2]*par[2]/((1.0-x)*par[3]*par[3]) ) > 1.0) {
459 // fitval = 0.0;
460 // return fitval;
461 // }
462 //
463 // fitval = par[0]*
464 // //(_alpha_em * _sigma_an/(16.0*PI))*
465 // ((1.0+(1.0-x)*(1.0-x))/x)*
466 // (pow(2.0,2.0*par[1]))*
467 // pow((1.0-(4.0*par[2]*par[2]/((1.0-x)*par[3]*par[3]))),(0.5+par[1]));
468 //
469 //
470 // return fitval;
471 //}
472
473 #endif
474
Attached Files
To refer to attachments on a page, use attachment:filename, as shown below in the list of files. Do NOT use the URL of the [get] link, since this is subject to change and can break easily.You are not allowed to attach a file to this page.