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.
  • [get | view] (2012-01-30 13:28:08, 4.3 KB) [[attachment:CalibrateData.C]]
  • [get | view] (2012-01-30 13:27:40, 12.0 KB) [[attachment:CalibrationFunctions.C]]
  • [get | view] (2014-06-20 09:25:03, 5.2 KB) [[attachment:CreateDataBackground.C]]
  • [get | view] (2014-06-20 09:21:28, 22.3 KB) [[attachment:CreateDataForSignalAndBackground.C]]
  • [get | view] (2014-06-20 09:24:20, 0.2 KB) [[attachment:CreateDataForSignalAndBackground.py]]
  • [get | view] (2012-01-30 12:32:43, 12.0 KB) [[attachment:FindMergeAngle.C]]
  • [get | view] (2012-01-30 13:22:33, 9.8 KB) [[attachment:PreprocessData.C]]
  • [get | view] (2012-01-30 13:51:59, 16.8 KB) [[attachment:Selection.C]]
  • [get | view] (2014-06-20 09:26:52, 13.2 KB) [[attachment:TSysLimit_searches.C]]
  • [get | view] (2014-06-20 09:28:25, 14.9 KB) [[attachment:TSysLimit_searches_onebin.C]]
  • [get | view] (2014-06-20 09:25:07, 12.9 KB) [[attachment:WIMP_functions.C]]
  • [get | view] (2014-06-20 12:52:10, 1.3 KB) [[attachment:runLambdaSearch.sh]]
  • [get | view] (2014-06-20 12:52:16, 0.2 KB) [[attachment:submitter.sh]]
 All files | Selected Files: delete move to page copy to page

You are not allowed to attach a file to this page.