PWGLF/RESONANCES/extra/Plot_XiStarResults.C

   1 //gROOT->Reset();
   2 #include <math.h>
   3 #include <time.h>
   4 #include <stdio.h>
   5 #include <stdlib.h>
   6 #include <Riostream.h>
   7 #include "RooFit.h"
   8 #include "RooVoigtian.h"
   9 #include "RooAbsReal.h"
  10 #include "RooRealVar.h"
  11 #include "RooComplex.h"
  12 #include "RooMath.h"
  13
  14 #include "TVector2.h"
  15 #include "TFile.h"
  16 #include "TString.h"
  17 #include "TF1.h"
  18 #include "TH1.h"
  19 #include "TH2.h"
  20 #include "TH3.h"
  21 #include "TProfile.h"
  22 #include "TProfile2D.h"
  23 #include "TMath.h"
  24 #include "TText.h"
  25 #include "TRandom3.h"
  26 #include "TArray.h"
  27 #include "TLegend.h"
  28 #include "TStyle.h"
  29 #include "TMinuit.h"
  30
  31 #include "TCanvas.h"
  32 #include "TLatex.h"
  33 #include "TImage.h"
  34 #include "TGaxis.h"
  35 #include "TGraphErrors.h"
  36 #include "TPaveStats.h"
  37
  38 using namespace std;
  39
  40
  41 double myLevyPtFunc(Double_t *x, Double_t *par);
  42 double myMtExpFunc(Double_t *x, Double_t *par);
  43 double myPtExpFunc(Double_t *x, Double_t *par);
  44 double myBoltzFunc(Double_t *x, Double_t *par);
  45 double myPowerLawPtFunc(Double_t *x, Double_t *par);
  46 double myPowerLawFunc(Double_t *x, Double_t *par);
  47 double myBoltsBWFunc(Double_t *x, Double_t *par);
  48 double IntegrandBG(const double *x, const double *p);
  49
  50 void Plot_XiStarResults(){
  51
  52   TGaxis::SetMaxDigits(3);
  53   gStyle->SetOptStat(0);
  54   gStyle->SetOptFit(01111);
  55
  56
  57   TFile *filesCutVar[13];
  58   TH1D *Spectra_CutVar[13];
  59   TF1 *LevyFits_CutVar[13];
  60   TH1D *Yields_CutVar[13];
  61   TH1D *Eff_CutVar[13];
  62   for(int i=0; i<13; i++){
  63     TString *name=new TString("CutVariation_");
  64     *name += i;
  65     name->Append(".root");
  66     filesCutVar[i]=new TFile(name->Data());
  67     Spectra_CutVar[i] = (TH1D*)filesCutVar[i]->Get("h_Spectrum");
  68     Spectra_CutVar[i]->SetDirectory(0);
  69     LevyFits_CutVar[i] = (TF1*)filesCutVar[i]->Get("myLevy");
  70     Yields_CutVar[i] = (TH1D*)filesCutVar[i]->Get("RawYields");
  71     Eff_CutVar[i] = (TH1D*)filesCutVar[i]->Get("Efficiency");
  72   }
  73   //
  74   TFile *filesFitVar[4];
  75   TH1D *Spectra_FitVar[4];
  76   TF1 *LevyFits_FitVar[4];
  77   TH1D *Yields_FitVar[4];
  78   for(int i=0; i<4; i++){
  79     TString *name=new TString("FitVariation_");
  80     *name += i+1;
  81     name->Append(".root");
  82     filesFitVar[i]=new TFile(name->Data());
  83     Spectra_FitVar[i] = (TH1D*)filesFitVar[i]->Get("h_Spectrum");
  84     Spectra_FitVar[i]->SetDirectory(0);
  85     LevyFits_FitVar[i] = (TF1*)filesFitVar[i]->Get("myLevy");
  86     Yields_FitVar[i] = (TH1D*)filesFitVar[i]->Get("RawYields");
  87   }
  88
  89   double FVavg=0;
  90   for(int bin=2; bin<=9; bin++){
  91     // FV 1 has largest spectrum point variation
  92     FVavg += fabs(Spectra_FitVar[0]->GetBinContent(bin)-Spectra_FitVar[1]->GetBinContent(bin))/Spectra_FitVar[0]->GetBinContent(bin);
  93   }
  94   cout<<"Average point-by-point spectrum variation due to fitting: "<<FVavg/8<<endl;
  95
  96   double TrackVavg=0;
  97   for(int bin=2; bin<=9; bin++){
  98     // CutVar 5 has largest spectrum point variation in the "tracking systematic" category
  99     TrackVavg += fabs(Spectra_CutVar[0]->GetBinContent(bin)-Spectra_CutVar[5]->GetBinContent(bin))/Spectra_CutVar[0]->GetBinContent(bin);
 100   }
 101   cout<<"Average point-by-point spectrum variation due to Tracking: "<<TrackVavg/8<<endl;
 102
 103   double TopologVavg=0;
 104   for(int bin=2; bin<=9; bin++){
 105     // CutVar 8 has largest spectrum point variation in the "topological systematic" category
 106     TopologVavg += fabs(Spectra_CutVar[0]->GetBinContent(bin)-Spectra_CutVar[8]->GetBinContent(bin))/Spectra_CutVar[0]->GetBinContent(bin);
 107   }
 108   cout<<"Average point-by-point spectrum variation due to Topological Cuts: "<<TopologVavg/8<<endl;
 109
 110   ////////////////////////////////////
 111   TH1D *FinalSpectrum = (TH1D*)Spectra_CutVar[0]->Clone();
 112   FinalSpectrum->SetMarkerSize(1.);
 113   double StatError[8]={0};
 114   for(int ptbin=1; ptbin<=8; ptbin++) StatError[ptbin-1] = FinalSpectrum->GetBinError(ptbin+1);// add 1 since first bin is a dummy
 115
 116
 117   TH1D *PureSysFromCutVar[12];
 118   double MeanSysCutVar[12]={0};
 119   double sigmaSysCutVar[12]={0};
 120   for(int i=0; i<12; i++) {
 121     TString *name=new TString("PureSysFromCutVar");
 122     *name += i+1;
 123     PureSysFromCutVar[i] = new TH1D(name->Data(),"",8,0.5,8.5);
 124     PureSysFromCutVar[i]->GetXaxis()->SetTitle("pt bin");
 125     PureSysFromCutVar[i]->GetYaxis()->SetTitle("|y_{s}-y_{v}|/y_{s}");
 126   }
 127
 128   //////////////////////////////////////
 129   // Sys error assignment for pt bins to be used for mean value fits
 130
 131   int CutVarLow=1, CutVarHigh=13;
 132   for(int i=CutVarLow; i<CutVarHigh; i++){// 1 to 13
 133     for(int ptbin=2; ptbin<=9; ptbin++){
 134       double SysVar = fabs(FinalSpectrum->GetBinContent(ptbin) - Spectra_CutVar[i]->GetBinContent(ptbin))/FinalSpectrum->GetBinContent(ptbin);
 135       double sigma = sqrt(fabs(pow(FinalSpectrum->GetBinError(ptbin),2) - pow(Spectra_CutVar[i]->GetBinError(ptbin),2)));// stat part
 136       PureSysFromCutVar[i-1]->SetBinContent(ptbin-1, SysVar);
 137       PureSysFromCutVar[i-1]->SetBinError(ptbin-1, sigma/FinalSpectrum->GetBinContent(ptbin));
 138     }
 139     PureSysFromCutVar[i-1]->Fit("pol0","IMEQ","",0,9);
 140     MeanSysCutVar[i-1] = pol0->GetParameter(0);
 141     sigmaSysCutVar[i-1] = pol0->GetParError(0);
 142     cout<<"Cut Var "<<i<<"  Mean Sys = "<<MeanSysCutVar[i-1]<<"   += "<<pol0->GetParError(0)<<endl;
 143     //cout<<"Chi2/NDF = "<<pol0->GetChisquare()/pol0->GetNDF()<<endl;
 144     //TPaveStats *sb2 = (TPaveStats*)(PureSysFromCutVar[i-1]->GetListOfFunctions()->FindObject("stats"));
 145     //sb2->SetX1NDC(.2);
 146   }
 147
 148
 149   double SysError[9]={0};
 150   // Cut Variations
 151   //int CutVarN=12;
 152   for(int i=CutVarLow; i<CutVarHigh; i++){// 1 to 13
 153     for(int ptbin=1; ptbin<=9; ptbin++){
 154       if(MeanSysCutVar[i-1] > 0 && MeanSysCutVar[i-1]>sigmaSysCutVar[i-1]) {
 155         SysError[ptbin-1] += pow(MeanSysCutVar[i-1]*FinalSpectrum->GetBinContent(ptbin),2) - pow(sigmaSysCutVar[i-1]*FinalSpectrum->GetBinContent(ptbin),2);
 156       }
 157     }
 158   }
 159   // Fit Variations
 160   for(int i=2; i<3; i++){// now only use i=2 (max deviation).  was 0 to 4
 161     for(int ptbin=1; ptbin<=9; ptbin++){
 162       SysError[ptbin-1] += pow(FinalSpectrum->GetBinContent(ptbin) - Spectra_FitVar[i]->GetBinContent(ptbin),2);
 163     }
 164   }
 165
 166   // Add Stat errors and set Final Error
 167   for(int ptbin=1; ptbin<=9; ptbin++){
 168     double TotError = pow(FinalSpectrum->GetBinError(ptbin),2);// Stat
 169     TotError += SysError[ptbin-1];// Add on Sys
 170     FinalSpectrum->SetBinError(ptbin, sqrt(TotError));
 171     cout<<"Spectrum ptbin:"<<ptbin-1<<"  "<<FinalSpectrum->GetBinContent(ptbin)<<"  +- "<<FinalSpectrum->GetBinError(ptbin)<<endl;
 172   }
 173
 174   for(int ptbin=1; ptbin<=9; ptbin++) SysError[ptbin-1] = sqrt(SysError[ptbin-1]);
 175
 176   // print efficiencies
 177   for(int ptbin=1; ptbin<=8; ptbin++){
 178     cout<<"Efficiency ptbin:"<<ptbin<<"  "<<Eff_CutVar[0]->GetBinContent(ptbin)<<"  +- "<<Eff_CutVar[0]->GetBinError(ptbin)<<endl;
 179   }
 180
 181   TCanvas *can = new TCanvas("can","can",13,34,700,500);
 182   can->Range(-1.25,-0.2625,11.25,2.3625);
 183   can->SetFillColor(10);
 184   can->SetBorderMode(0);
 185   can->SetBorderSize(2);
 186   //can->SetGridx();
 187   //can->SetGridy();
 188   can->SetFrameFillColor(0);
 189   can->SetFrameBorderMode(0);
 190   can->SetFrameBorderMode(0);
 191
 192   TLegend *legend = new TLegend(.1,.7,.35,.9,NULL,"brNDC");
 193   legend->SetBorderSize(1);
 194   legend->SetTextSize(.04);// small .03; large .036
 195   //legend->SetLineColor(0);
 196   legend->SetFillColor(0);
 197   TLegend *legend2 = new TLegend(.35,.62,.9,.72,NULL,"brNDC");
 198   legend2->SetBorderSize(1);
 199   legend2->SetTextSize(.04);// small .03; large .036
 200   //legend2->SetLineColor(0);
 201   legend2->SetFillColor(0);
 202
 203
 204   double minFitpoint = 0.8; double maxFitpoint = 5.6;
 205
 206
 207
 208   TF1 *myLevy=new TF1("myLevy",myLevyPtFunc,0,8.5,3);
 209   myLevy->SetParName(0,"dN/dy");
 210   myLevy->SetParName(1,"C");
 211   myLevy->SetParName(2,"n");
 212   myLevy->SetParameter(0,.008);
 213   myLevy->SetParameter(1,.3);
 214   myLevy->SetParameter(2,15);
 215   myLevy->SetParLimits(0,.001,.2);
 216   myLevy->SetParLimits(1,.1,1);
 217   myLevy->SetParLimits(2,1,500);
 218   myLevy->SetLineColor(1);
 219   TF1 *myPtExp = new TF1("myPtExp",myPtExpFunc,0,8.5,2);
 220   myPtExp->SetParName(0,"dN/dy");
 221   myPtExp->SetParName(1,"C");
 222   myPtExp->SetParameter(0,.003);
 223   myPtExp->SetParameter(1,.3);
 224   myPtExp->SetParLimits(0,0.0001,0.1);
 225   myPtExp->SetParLimits(1,.1,1);
 226   myPtExp->SetLineColor(2);
 227   TF1 *myMtExp = new TF1("myMtExp",myMtExpFunc,0,8.5,2);
 228   myMtExp->SetParName(0,"dN/dy");
 229   myMtExp->SetParName(1,"C");
 230   myMtExp->SetParameter(0,.003);
 231   myMtExp->SetParameter(1,.3);
 232   myMtExp->SetParLimits(0,0.0001,0.1);
 233   myMtExp->SetParLimits(1,.1,1);
 234   myMtExp->SetLineColor(4);
 235   TF1 *myBoltz = new TF1("myBoltz",myBoltzFunc,0,8.5,2);
 236   myBoltz->SetParName(0,"dN/dy");
 237   myBoltz->SetParName(1,"C");
 238   myBoltz->SetParameter(0,.003);
 239   myBoltz->SetParameter(1,.3);
 240   myBoltz->SetParLimits(0,0.0001,0.1);
 241   myBoltz->SetParLimits(1,.1,1);
 242   myBoltz->SetLineColor(6);
 243   TF1 *myPowerLawPt=new TF1("myPowerLawPt",myPowerLawPtFunc,0,8.5,3);
 244   myPowerLawPt->SetParName(0,"dN/dy");
 245   myPowerLawPt->SetParName(1,"pt0");
 246   myPowerLawPt->SetParName(2,"n");
 247   myPowerLawPt->SetParameter(0,.008);
 248   myPowerLawPt->SetParameter(1,50);
 249   myPowerLawPt->SetParameter(2,50);
 250   myPowerLawPt->SetParLimits(0,.001,.01);
 251   myPowerLawPt->SetParLimits(1,1,500);
 252   myPowerLawPt->SetParLimits(2,1,500);
 253   myPowerLawPt->SetLineColor(7);
 254   TF1 *myPowerLaw=new TF1("myPowerLaw",myPowerLawFunc,0,8.5,3);
 255   myPowerLaw->SetParName(0,"dN/dy");
 256   myPowerLaw->SetParName(1,"pt0");
 257   myPowerLaw->SetParName(2,"n");
 258   myPowerLaw->SetParameter(0,.008);
 259   myPowerLaw->SetParameter(1,50);
 260   myPowerLaw->SetParameter(2,50);
 261   myPowerLaw->SetParLimits(0,.001,.01);
 262   myPowerLaw->SetParLimits(1,1,500);
 263   myPowerLaw->SetParLimits(2,1,500);
 264   myPowerLaw->SetLineColor(8);
 265   TF1 *myBGBlastWave=new TF1("myBGBlastWave",myBoltsBWFunc,0,8.5,5);
 266   myBGBlastWave->SetParName(0,"mass");
 267   myBGBlastWave->SetParName(1,"beta");
 268   myBGBlastWave->SetParName(2,"C");
 269   myBGBlastWave->SetParName(3,"n");
 270   myBGBlastWave->SetParName(4,"Norm");
 271   myBGBlastWave->FixParameter(0,1.5318);
 272   myBGBlastWave->SetParameter(1,.3);
 273   myBGBlastWave->SetParameter(2,.3);
 274   myBGBlastWave->SetParameter(3,5);
 275   myBGBlastWave->SetParameter(4,.3);
 276   myBGBlastWave->SetParLimits(1,.001,.9);
 277   myBGBlastWave->SetParLimits(2,.1,200);
 278   myBGBlastWave->SetParLimits(3,.1,200);
 279   myBGBlastWave->SetLineColor(9);
 280
 281
 282   FinalSpectrum->Draw();
 283   legend->AddEntry(FinalSpectrum,"(#Xi(1530) + cc)/2","p");
 284   //legend->AddEntry(FinalSpectrum,"(#Xi(1530) + cc)/2: Current Results","p");
 285   FinalSpectrum->Fit(myLevy,"IME","",minFitpoint,maxFitpoint);
 286   FinalSpectrum->Fit(myPtExp,"IMEN","",minFitpoint,maxFitpoint);
 287   FinalSpectrum->Fit(myMtExp,"IMEN","",minFitpoint,maxFitpoint);
 288   FinalSpectrum->Fit(myBoltz,"IMEN","",minFitpoint,maxFitpoint);
 289   FinalSpectrum->Fit(myPowerLaw,"IMEN","",minFitpoint,maxFitpoint);
 290   FinalSpectrum->Fit(myPowerLawPt,"IMEN","",minFitpoint,maxFitpoint);
 291   //FinalSpectrum->Fit(myBGBlastWave,"IME","",minFitpoint,maxFitpoint);
 292   myLevy->Draw("same");
 293   /*myPtExp->Draw("same");
 294   myMtExp->Draw("same");
 295   myBoltz->Draw("same");
 296   myPowerLaw->Draw("same");
 297   myPowerLawPt->Draw("same");*/
 298   //myBGBlastWave->Draw("same");
 299   //
 300   legend->AddEntry(myLevy,"Levy","l");
 301   /*legend->AddEntry(myPtExp,"Pt-exp","l");
 302   legend->AddEntry(myMtExp,"Mt-exp","l");
 303   legend->AddEntry(myBoltz,"Pt-Boltzmann","l");
 304   legend->AddEntry(myPowerLaw,"PowerLaw","l");
 305   legend->AddEntry(myPowerLawPt,"Pt-PowerLaw","l");
 306   */
 307   double AvgLowPtExtrap = myLevy->Integral(0,0.8) + myPtExp->Integral(0,0.8) + myPowerLawPt->Integral(0,0.8) + myMtExp->Integral(0,0.8);
 308   AvgLowPtExtrap /= 4.;
 309   cout<<"Levy Extrap % = "<<myLevy->Integral(0,0.8)/myLevy->Integral(0,10.)<<"  Avg Extrap % = "<<AvgLowPtExtrap/myLevy->Integral(0,10.)<<endl;
 310   //double LowPtExtrapUncertaintyPercentage = 0.17;// Levy underestimation of Pythia
 311   //LowPtExtrapUncertaintyPercentage -= (AvgLowPtExtrap - myLevy->Integral(0,0.8))/myLevy->Integral(0,10.);
 312   double High_lowptExtrap = (myPtExp->Integral(0,0.8) + myPowerLawPt->Integral(0,0.8))/2.;
 313   double Low_lowptExtrap = myMtExp->Integral(0,0.8);
 314   double High_LowPtExtrapUncertaintyPercentage = (High_lowptExtrap - myLevy->Integral(0,0.8))/myLevy->Integral(0,10.);
 315   double Low_LowPtExtrapUncertaintyPercentage = (myLevy->Integral(0,0.8) - Low_lowptExtrap)/myLevy->Integral(0,10.);
 316   cout<<"high low-pt uncertainty = "<<High_LowPtExtrapUncertaintyPercentage<<endl;
 317   cout<<"low low-pt uncertainty = "<<Low_LowPtExtrapUncertaintyPercentage<<endl;
 318
 319
 320   //PureSysFromCutVar->Draw();
 321
 322
 323   double dNdY = myLevy->Integral(0,10);
 324   double dNdY_exp = myMtExp->Integral(0,10);
 325   double dNdY_Stat = LevyFits_CutVar[0]->GetParError(0);
 326   double dNdY_Covered = myLevy->Integral(0.8,5.6);
 327   double dNdY_Covered_exp = myMtExp->Integral(0.8,5.6);
 328   double dNdY_DataPoints = 0, dNdY_DataPoints_e=0;
 329   for(int ii=0; ii<8; ii++){// to explicitly integrate dNdY without fit (Christina's suggestion)
 330     dNdY_DataPoints += Spectra_CutVar[0]->GetBinContent(ii+2)*Spectra_CutVar[0]->GetBinWidth(ii+2);
 331     dNdY_DataPoints_e += pow(Spectra_CutVar[0]->GetBinError(ii+2)*Spectra_CutVar[0]->GetBinWidth(ii+2),2);
 332   }
 333   dNdY_DataPoints_e = sqrt(dNdY_DataPoints_e);
 334   double dNdYSysPlus = 0;
 335   double dNdYSysMinus = 0;
 336   double CSys = 0;
 337   double nSys = 0;
 338   //////////////////////////////////////
 339   // Sys errors of Levy parameters
 340   double sigma_dNdY_cv[12]={0};
 341   double sigma_C_cv[12]={0};
 342   double sigma_n_cv[12]={0};
 343   double SysVar_dNdY_cv[12]={0};
 344   double SysVar_C_cv[12]={0};
 345   double SysVar_n_cv[12]={0};
 346
 347   for(int cv=1; cv<13; cv++){
 348     if(MeanSysCutVar[cv-1] > 0 && MeanSysCutVar[cv-1] > sigmaSysCutVar[cv-1]) {// significant deviation
 349       sigma_dNdY_cv[cv-1] = sqrt(fabs(pow(LevyFits_CutVar[0]->GetParError(0),2) - pow(LevyFits_CutVar[cv]->GetParError(0),2)));
 350       sigma_C_cv[cv-1] = sqrt(fabs(pow(LevyFits_CutVar[0]->GetParError(1),2) - pow(LevyFits_CutVar[cv]->GetParError(1),2)));
 351       sigma_n_cv[cv-1] = sqrt(fabs(pow(LevyFits_CutVar[0]->GetParError(2),2) - pow(LevyFits_CutVar[cv]->GetParError(2),2)));
 352       SysVar_dNdY_cv[cv-1] = fabs(LevyFits_CutVar[0]->GetParameter(0)-LevyFits_CutVar[cv]->GetParameter(0));
 353       SysVar_C_cv[cv-1] = fabs(LevyFits_CutVar[0]->GetParameter(1)-LevyFits_CutVar[cv]->GetParameter(1));
 354       SysVar_n_cv[cv-1] = fabs(LevyFits_CutVar[0]->GetParameter(2)-LevyFits_CutVar[cv]->GetParameter(2));
 355
 356       if(SysVar_dNdY_cv[cv-1] > sigma_dNdY_cv[cv-1]){
 357         dNdYSysPlus += pow(SysVar_dNdY_cv[cv-1],2)-pow(sigma_dNdY_cv[i-1],2);
 358         dNdYSysMinus += pow(SysVar_dNdY_cv[cv-1],2)-pow(sigma_dNdY_cv[i-1],2);
 359       }
 360       if(SysVar_C_cv[cv-1] > sigma_C_cv[cv-1]){
 361         CSys += pow(SysVar_C_cv[cv-1],2)-pow(sigma_C_cv[cv-1],2);
 362       }
 363       if(SysVar_n_cv[cv-1] > sigma_n_cv[cv-1]){
 364         nSys += pow(SysVar_n_cv[cv-1],2)-pow(sigma_n_cv[cv-1],2);
 365       }
 366
 367     }
 368
 369   }
 370
 371
 372
 373   double sigma_dNdY_fv[4]={0};
 374   double sigma_C_fv[4]={0};
 375   double sigma_n_fv[4]={0};
 376   double SysVar_dNdY_fv[4]={0};
 377   double SysVar_C_fv[4]={0};
 378   double SysVar_n_fv[4]={0};
 379   for(int fv=0; fv<4; fv++){
 380     sigma_dNdY_fv[fv] = sqrt(fabs(pow(LevyFits_CutVar[0]->GetParError(0),2) - pow(LevyFits_FitVar[fv]->GetParError(0),2)));
 381     SysVar_dNdY_fv[fv] = sqrt(pow(LevyFits_CutVar[0]->GetParameter(0)-LevyFits_FitVar[fv]->GetParameter(0),2));
 382     sigma_C_fv[fv] = sqrt(fabs(pow(LevyFits_CutVar[0]->GetParError(1),2) - pow(LevyFits_FitVar[fv]->GetParError(1),2)));
 383     SysVar_C_fv[fv] = sqrt(pow(LevyFits_CutVar[0]->GetParameter(1)-LevyFits_FitVar[fv]->GetParameter(1),2));
 384     sigma_n_fv[fv] = sqrt(fabs(pow(LevyFits_CutVar[0]->GetParError(2),2) - pow(LevyFits_FitVar[fv]->GetParError(2),2)));
 385     SysVar_n_fv[fv] = sqrt(pow(LevyFits_CutVar[0]->GetParameter(2)-LevyFits_FitVar[fv]->GetParameter(2),2));
 386
 387     if(fv==2){//max case for dN/dy
 388       //dNdYSysPlus += pow(LevyFits_CutVar[0]->GetParameter(0)-LevyFits_FitVar[fv]->GetParameter(0),2) - pow(LevyFits_CutVar[0]->GetParError(0)-LevyFits_FitVar[fv]->GetParError(0),2);
 389       //dNdYSysMinus += pow(LevyFits_CutVar[0]->GetParameter(0)-LevyFits_FitVar[fv]->GetParameter(0),2) - pow(LevyFits_CutVar[0]->GetParError(0)-LevyFits_FitVar[fv]->GetParError(0),2);
 390       dNdYSysPlus += pow(SysVar_dNdY_fv[fv],2) - pow(sigma_dNdY_fv[fv],2);
 391       dNdYSysMinus += pow(SysVar_dNdY_fv[fv],2) - pow(sigma_dNdY_fv[fv],2);
 392     }
 393     //if(fv==2) CSys += pow(LevyFits_CutVar[0]->GetParameter(1)-LevyFits_FitVar[fv]->GetParameter(1),2) - pow(LevyFits_CutVar[0]->GetParError(1)-LevyFits_FitVar[fv]->GetParError(1),2);
 394     //if(fv==1) nSys += pow(LevyFits_CutVar[0]->GetParameter(2)-LevyFits_FitVar[fv]->GetParameter(2),2) - pow(LevyFits_CutVar[0]->GetParError(2)-LevyFits_FitVar[fv]->GetParError(2),2);
 395     if(fv==2) CSys += pow(SysVar_C_fv[fv],2) - pow(sigma_C_fv[fv],2);
 396     if(fv==1) nSys += pow(SysVar_n_fv[fv],2) - pow(sigma_n_fv[fv],2);
 397   }
 398
 399
 400   double dNdYSys_forRatio = dNdYSysPlus;
 401   dNdYSysPlus += pow(0.062*dNdY,2);// INEL Norm uncertainty
 402   dNdYSysPlus += pow(0.04*dNdY,2);// Material Budget uncertainty
 403   dNdYSysPlus += pow(0.01*dNdY,2);// Geant3/Fluka correction uncertainty
 404
 405   dNdYSysPlus += pow(High_LowPtExtrapUncertaintyPercentage * dNdY,2);// extrapolation uncertainty (was LowPtExtrapUncertaintyPercentage * fabs(dNdY-dNdY_Covered))
 406   dNdYSys_forRatio += pow(High_LowPtExtrapUncertaintyPercentage * dNdY,2);// extrapolation uncertainty
 407   dNdYSysPlus = sqrt(dNdYSysPlus);
 408   dNdYSys_forRatio = sqrt(dNdYSys_forRatio);
 409   //
 410   dNdYSysMinus += pow(0.03*dNdY,2);// INEL Norm uncertainty
 411   dNdYSysMinus += pow(0.04*dNdY,2);// Material Budget uncertainty
 412   dNdYSysMinus += pow(0.01*dNdY,2);// Geant3/Fluka correction uncertainty
 413   dNdYSysMinus += pow(Low_LowPtExtrapUncertaintyPercentage * dNdY,2);// extrapolation uncertainty (40% from Pythia study).
 414   dNdYSysMinus = sqrt(dNdYSysMinus);
 415
 416   CSys = sqrt(CSys);
 417   nSys = sqrt(nSys);
 418
 419   cout<<"dN/dy = "<<dNdY<<"  , Covered range = "<<dNdY_Covered<<"  , DataPoint Sum = "<<dNdY_DataPoints<<"  , stat = "<<dNdY_DataPoints_e<<endl;
 420   cout<<"dN/dy stat = "<<dNdY_Stat<<"   +sys = "<<dNdYSysPlus<<"  -sys = "<<dNdYSysMinus<<endl;
 421
 422   cout<<"C = "<<myLevy->GetParameter(1)<<"  , stat = "<<LevyFits_CutVar[0]->GetParError(1)<<"  , sys = "<<CSys<<endl;
 423   cout<<"n = "<<myLevy->GetParameter(2)<<"  , stat = "<<LevyFits_CutVar[0]->GetParError(2)<<"  , sys = "<<nSys<<endl;
 424
 425   // mean pt
 426   double ptEdges[9]={0.8, 1.2, 1.6, 2.0, 2.4, 3.2, 4.0, 4.8, 5.6};
 427   double meanpt_def=0, meanpt_Sys1=0, meanpt_Sys2=0, meanpt_Sys3=0;// Sys1 = Spectrum only, Sys2 = bin centers, Sys3 = low-pt extrapolation
 428   double yieldSum_def=0, yieldSum_Sys1=0, yieldSum_Sys2=0, yieldSum_Sys3=0;
 429   double meanpt_statError=0;
 430   meanpt_def += myLevy->Moment(1,0.,0.8) * myLevy->Integral(0,0.8);// low-pt extrapolation
 431   yieldSum_def += myLevy->Integral(0,0.8);
 432   meanpt_Sys1 += meanpt_def; meanpt_Sys2 += meanpt_def;
 433   yieldSum_Sys1 += yieldSum_def; yieldSum_Sys2 += yieldSum_def;
 434   meanpt_Sys3 += fabs(myLevy->Integral(0,0.8) - High_LowPtExtrapUncertaintyPercentage*dNdY) * myLevy->Moment(1,0.,0.8);// low-pt extrapolation uncertainty
 435   yieldSum_Sys3 += fabs(myLevy->Integral(0,0.8) - High_LowPtExtrapUncertaintyPercentage*dNdY);// low-pt extrapolation uncertainty
 436   meanpt_statError += myLevy->Moment(1,0.,0.8) * myLevy->Integral(0,0.8) * myLevy->GetParError(0)/myLevy->GetParameter(0);
 437
 438
 439   double StatErrorPtBins=0;
 440   for(int bin=1; bin<=8; bin++){
 441     double bin_width = fabs(ptEdges[bin-1]-ptEdges[bin]);
 442     meanpt_def += myLevy->Moment(1,ptEdges[bin-1],ptEdges[bin]) * FinalSpectrum->GetBinContent(bin+1)*bin_width;
 443     yieldSum_def += FinalSpectrum->GetBinContent(bin+1)*bin_width;
 444     meanpt_Sys1 += myLevy->Moment(1,ptEdges[bin-1],ptEdges[bin]) * myLevy->Integral(ptEdges[bin-1],ptEdges[bin]);
 445     yieldSum_Sys1 += myLevy->Integral(ptEdges[bin-1],ptEdges[bin]);
 446     meanpt_Sys2 += (ptEdges[bin]+ptEdges[bin-1])/2. * FinalSpectrum->GetBinContent(bin+1)*bin_width;
 447     yieldSum_Sys2 += FinalSpectrum->GetBinContent(bin+1)*bin_width;
 448     meanpt_Sys3 += myLevy->Moment(1,ptEdges[bin-1],ptEdges[bin]) * FinalSpectrum->GetBinContent(bin+1)*bin_width;
 449     yieldSum_Sys3 += FinalSpectrum->GetBinContent(bin+1)*bin_width;
 450     //
 451     StatErrorPtBins += pow(myLevy->Moment(1,ptEdges[bin-1],ptEdges[bin]) * Spectra_CutVar[0]->GetBinError(bin+1)*bin_width,2);
 452   }
 453
 454   meanpt_def += myLevy->Moment(1,5.6,10.) * myLevy->Integral(5.6,10.);// high-pt extrapolation
 455   yieldSum_def += myLevy->Integral(5.6,10.);
 456   meanpt_Sys1 += myLevy->Moment(1,5.6,10.) * myLevy->Integral(5.6,10.);// high-pt extrapolation
 457   yieldSum_Sys1 += myLevy->Integral(5.6,10.);
 458   meanpt_Sys2 += myLevy->Moment(1,5.6,10.) * myLevy->Integral(5.6,10.);// high-pt extrapolation
 459   yieldSum_Sys2 += myLevy->Integral(5.6,10.);
 460   meanpt_Sys3 += fabs(myLevy->Integral(5.6,10.) - High_LowPtExtrapUncertaintyPercentage*myLevy->Integral(5.6,10.))*myLevy->Moment(1,5.6,10.);// high-pt extrapolation
 461   yieldSum_Sys3 += fabs(myLevy->Integral(5.6,10.) - High_LowPtExtrapUncertaintyPercentage*myLevy->Integral(5.6,10.));
 462   //
 463   meanpt_statError += sqrt(StatErrorPtBins);
 464   meanpt_statError += myLevy->Moment(1,5.6,10.) * myLevy->Integral(5.6,10.) * myLevy->GetParError(0)/myLevy->GetParameter(0);// high-pt extrapolation
 465   meanpt_statError /= yieldSum_def;
 466
 467   meanpt_def /= yieldSum_def;
 468   meanpt_Sys1 /= yieldSum_Sys1;
 469   meanpt_Sys2 /= yieldSum_Sys2;
 470   meanpt_Sys3 /= yieldSum_Sys3;
 471
 472   double meanpt_Sys = sqrt(pow(meanpt_def-meanpt_Sys1,2)+pow(meanpt_def-meanpt_Sys2,2)+pow(meanpt_def-meanpt_Sys3,2));
 473
 474   cout<<"<pt> = "<<meanpt_def<<"  , stat = "<<meanpt_statError<<"  , sys = "<<meanpt_Sys<<endl;
 475
 476   /*double meanpt = myLevy->Moment(1,0.,10.);// 0.,30 for full range
 477   double meanpt_CoveredLevy = myLevy->Moment(1,0.8,5.6);
 478   double mPt_Stat = 0;
 479   double mPt_Sys = 0;
 480   double mPt_StatSys = 0;
 481   double base_C = myLevy->GetParameter(1);
 482   double base_n = myLevy->GetParameter(2);
 483
 484
 485   // Sys
 486   myLevy->SetParameter(1, base_C + myLevy->GetParError(1));
 487   mPt_StatSys += pow(myLevy->Moment(1,0.,10)-meanpt,2);
 488   //
 489   //
 490   myLevy->SetParameter(2, base_n + myLevy->GetParError(2));
 491   mPt_StatSys += pow(myLevy->Moment(1,0.,10)-meanpt,2);
 492   mPt_StatSys = sqrt(mPt_StatSys);
 493   //
 494   //
 495   // Stat
 496   myLevy->SetParameter(1, base_C + LevyFits_CutVar[0]->GetParError(1));
 497   mPt_Stat += pow(myLevy->Moment(1,0.,10)-meanpt,2);
 498   //
 499   //
 500   myLevy->SetParameter(2, base_n + LevyFits_CutVar[0]->GetParError(2));
 501   mPt_Stat += pow(myLevy->Moment(1,0.,10)-meanpt,2);
 502   mPt_Stat = sqrt(mPt_Stat);
 503   //
 504   myLevy->SetParameter(1, base_C); myLevy->SetParameter(1, base_n);
 505
 506   cout<<"Levy mean pt = "<<meanpt<<"  , stat = "<<mPt_Stat<<"  , sys = "<<sqrt( pow(mPt_StatSys,2)-pow(mPt_Stat,2) )<<endl;
 507   cout<<"Levy mean pt in covered range = "<<meanpt_CoveredLevy<<endl;
 508   */
 509
 510   double XidNdY = 0.0079;// (Xi^+ + Xi^-)/2. Published
 511   double XiStarToXi = dNdY/XidNdY;
 512   double XiStarToXi_Stat = sqrt(pow(dNdY_Stat/XidNdY,2) + pow(dNdY*0.0001/pow(XidNdY,2),2));// 0.0001 is stat error of Xi^{\pm}
 513   double XiStartoXi_Sys = dNdYSys_forRatio/XidNdY;
 514
 515
 516   cout<<"Xi*/Xi = "<<XiStarToXi<<"  , stat = "<<XiStarToXi_Stat<<"  , sys = "<<XiStartoXi_Sys<<endl;
 517   cout<<endl;
 518   cout<<endl;
 519
 520   //Old QM2011 values for (Xi* + cc)/2
 521   double Old_XiStarPoints[8]={0.001259, 0.0008586, 0.0006316, 0.0003604, 1.856e-4, 5.4601e-5, 2.0187e-5, 7.6559e-6};
 522   double Old_XiStarPoints_e[8]={2.44e-4, 1.224e-4, 7.136e-5, 3.85e-5, 1.59e-5, 7.49e-6, 3.51e-6, 1.822e-6};
 523   double pt_points[8]={1.0, 1.4, 1.8, 2.2, 2.8, 3.6, 4.4, 5.2};
 524   double pt_points_e[8]={.2, .2, .2, .2, .4, .4, .4, .4};
 525   TGraphErrors *gr_OldXiStar = new TGraphErrors(8, pt_points, Old_XiStarPoints, pt_points_e, Old_XiStarPoints_e);
 526   gr_OldXiStar->SetMarkerStyle(21);
 527   gr_OldXiStar->SetMarkerColor(2);
 528   gr_OldXiStar->SetLineColor(2);
 529   double NewOldRatio[8]={0};
 530   for(int ii=0; ii<8; ii++){
 531     NewOldRatio[ii] = FinalSpectrum->GetBinContent(ii+2)/Old_XiStarPoints[ii];
 532     //cout<<FinalSpectrum->GetBinError(ii+2)/FinalSpectrum->GetBinContent(ii+2)<<endl;
 533   }
 534   TGraph *gr_NewOldRatio = new TGraph(8, pt_points, NewOldRatio);
 535   gr_NewOldRatio->SetMarkerStyle(20);
 536   gr_NewOldRatio->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 537   gr_NewOldRatio->GetYaxis()->SetTitle("New/Old");
 538   gr_NewOldRatio->SetTitle("New Spectrum Divided by QM11 Spectrum");
 539   //gr_NewOldRatio->Draw("AP");
 540
 541   //gr_OldXiStar->Draw("P");
 542   //legend->AddEntry(gr_OldXiStar,"(#Xi(1530) + cc)/2: QM11 Results","p");
 543   // To Fit Old points
 544   /*TH1D *h_OldXiStar = (TH1D*)h_Xispectrum->Clone();
 545   h_OldXiStar->Add(h_OldXiStar,-1);
 546   for(int i=0; i<8; i++){
 547     h_OldXiStar->SetBinContent(i+2, Old_XiStarPoints[i]);
 548     h_OldXiStar->SetBinError(i+2, Old_XiStarPoints_e[i]);
 549   }
 550   h_OldXiStar->SetBinContent(1,+100);// first bin has no data
 551   h_OldXiStar->Fit(myLevy,"IME","",minFitpoint,maxFitpoint);
 552   */
 553   //legend->Draw("same");
 554
 555   double StatPercentage[8]={};
 556   double SysPercentage[8]={};
 557   double TotErrPercentage[8]={};
 558   for(int ii=0; ii<8; ii++){
 559     StatPercentage[ii] = StatError[ii]/FinalSpectrum->GetBinContent(ii+2);
 560     SysPercentage[ii] = SysError[ii+1]/FinalSpectrum->GetBinContent(ii+2);
 561     TotErrPercentage[ii] = sqrt(pow(StatError[ii],2)+pow(SysError[ii+1],2))/FinalSpectrum->GetBinContent(ii+2);
 562   }
 563   TGraph *gr_StatPercentage = new TGraph(8, pt_points, StatPercentage);
 564   gr_StatPercentage->SetMarkerStyle(25);
 565   gr_StatPercentage->SetMarkerColor(4);
 566   gr_StatPercentage->SetLineColor(4);
 567   gr_StatPercentage->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 568   gr_StatPercentage->GetYaxis()->SetTitle("Uncertainty (fraction)");
 569   gr_StatPercentage->SetTitle("");
 570   gr_StatPercentage->SetMinimum(0);
 571   gr_StatPercentage->SetMaximum(0.2);// 0.35 to compare with Enrico
 572   //gr_StatPercentage->Draw("APC");
 573   TGraph *gr_SysPercentage = new TGraph(8, pt_points, SysPercentage);
 574   gr_SysPercentage->SetMarkerStyle(25);
 575   gr_SysPercentage->SetMarkerColor(2);
 576   gr_SysPercentage->SetLineColor(2);
 577   //gr_SysPercentage->Draw("PC");
 578   TGraph *gr_TotErrPercentage = new TGraph(8, pt_points, TotErrPercentage);
 579   gr_TotErrPercentage->SetMarkerStyle(25);
 580   gr_TotErrPercentage->SetMarkerColor(1);
 581   gr_TotErrPercentage->SetLineColor(1);
 582   //gr_TotErrPercentage->Draw("PC");
 583   //legend->AddEntry(gr_StatPercentage,"Stat. uncertainty","p");
 584   //legend->AddEntry(gr_SysPercentage,"Sys. uncertainty","p");
 585   //legend->AddEntry(gr_TotErrPercentage,"#sqrt{Stat^{2}+Sys^{2}}","p");
 586   legend->Draw("same");
 587
 588   ////////////////////////////////
 589   // Cut Var
 590   for(int var=1; var<13; var++){// 0 to 4 for FitVar,  1 to 13 for CutVar
 591     TString *Varname = new TString();
 592     if(var==1) Varname->Append("Nclusters TPC");
 593     if(var==2) Varname->Append("DCA PV proton");
 594     if(var==3) Varname->Append("DCA PV 1st pion");
 595     if(var==4) Varname->Append("DCA PV 2nd pion");
 596     if(var==5) Varname->Append("DCA PV 3rd pion");
 597     if(var==6) Varname->Append("DCA PV 4th pion");
 598     if(var==7) Varname->Append("DCA proton-pion");
 599     if(var==8) Varname->Append("DCA Lambda-pion");
 600     if(var==9) Varname->Append("Decay length xy Lambda");
 601     if(var==10) Varname->Append("Decay length xy Xi");
 602     if(var==11) Varname->Append("Cos PointingAngle Lambda");
 603     if(var==12) Varname->Append("Cos PointingAngle Xi");
 604
 605     cout<<Varname->Data()<<" & ";
 606     cout.precision(2);
 607     for(int par=0; par<3; par++) {
 608       cout<<100.*(LevyFits_CutVar[var]->GetParameter(par) - LevyFits_CutVar[0]->GetParameter(par))/LevyFits_CutVar[0]->GetParameter(par)<<"(";
 609       if(par==0) cout<<100.*sigma_dNdY_cv[var-1]/LevyFits_CutVar[0]->GetParameter(par)<<")"<<" ";
 610       if(par==1) cout<<100.*sigma_C_cv[var-1]/LevyFits_CutVar[0]->GetParameter(par)<<")"<<" ";
 611       if(par==2) cout<<100.*sigma_n_cv[var-1]/LevyFits_CutVar[0]->GetParameter(par)<<")"<<" ";
 612
 613       if((par+1)!=3) cout<<"& ";
 614       else cout<<" \\\\\ \\hline";
 615     }
 616     cout<<endl;
 617   }
 618   cout<<endl;
 619
 620   //////////////////////////////
 621   // Fit Var
 622   for(int var=0; var<4; var++){// 0 to 4 for FitVar,  1 to 13 for CutVar
 623     TString *Varname = new TString();
 624     if(var==0) Varname->Append("Pure Voigtian");
 625     if(var==1) Varname->Append("Pure Voigtian w/o bkg subtraction");
 626     if(var==2) Varname->Append("bkg norm to the right");
 627     if(var==3) Varname->Append("1/2 bin-counting range");
 628
 629     cout<<Varname->Data()<<" & ";
 630     cout.precision(2);
 631     for(int par=0; par<3; par++) {
 632       cout<<100.*(LevyFits_FitVar[var]->GetParameter(par) - LevyFits_CutVar[0]->GetParameter(par))/(LevyFits_CutVar[0]->GetParameter(par))<<"(";
 633       if(par==0) cout<<100.*sigma_dNdY_fv[var]/LevyFits_CutVar[0]->GetParameter(par)<<")"<<" ";
 634       if(par==1) cout<<100.*sigma_C_fv[var]/LevyFits_CutVar[0]->GetParameter(par)<<")"<<" ";
 635       if(par==2) cout<<100.*sigma_n_fv[var]/LevyFits_CutVar[0]->GetParameter(par)<<")"<<" ";
 636
 637       if((par+1)!=3) cout<<"& ";
 638       else cout<<" \\\\\ \\hline";
 639     }
 640     cout<<endl;
 641   }
 642   cout<<endl;
 643
 644   ///////////////////////////////////////////
 645
 646   /*
 647   cout<<"Yields"<<endl;
 648   // print yields
 649   for(int var=1; var<=13; var++){
 650     TString *Varname = new TString();
 651     if(var==1) Varname->Append("Standard ");
 652     else {Varname->Append("Cut Var "); *Varname += var-1;}
 653
 654     cout<<Varname->Data()<<" & ";
 655     cout.precision(4);// 4 or 6
 656     for(int ptbin=1; ptbin<=8; ptbin++){
 657       //cout<<Yields_CutVar[var-1]->GetBinContent(ptbin)/Eff_CutVar[var-1]->GetBinContent(ptbin)<<" ";
 658       cout<<Eff_CutVar[var-1]->GetBinContent(ptbin)<<" ";
 659       if((ptbin)!=8) cout<<"& ";
 660       else cout<<" \\\\\ \\hline";
 661     }
 662     cout<<endl;
 663   }
 664   */
 665   /* for(int var=1; var<5; var++){
 666     TString *Varname = new TString("Fit Var ");
 667     *Varname += var;
 668     cout<<Varname->Data()<<" & ";
 669     cout.precision(6);
 670     for(int ptbin=1; ptbin<=8; ptbin++){
 671       cout<<Yields_FitVar[var-1]->GetBinContent(ptbin)/Eff_CutVar[0]->GetBinContent(ptbin)<<" ";
 672       if((ptbin)!=8) cout<<"& ";
 673       else cout<<" \\\\\ \\hline";
 674     }
 675     cout<<endl;
 676     }*/
 677   cout<<endl;
 678
 679   //TString *name=new TString("../AnaNotes/XiStar/XiStar_yield_ErrorsCutVar");
 680   //*name += CutVarN;
 681   /*TString *name=new TString("../AnaNotes/XiStar/XiStar_CutVarSys");
 682   *name += CutVarLow;
 683   name->Append(".png");
 684   c1->SaveAs(name->Data());*/
 685   //TString *name=new TString("../AnaNotes/XiStar/XiStar_CutVarSummary");
 686   //*name += CutVarLow;
 687   //name->Append(".png");
 688   //can->SaveAs(name->Data());
 689
 690 }
 691
 692
 693 //________________________________________________________________________
 694 double myLevyPtFunc(Double_t *x, Double_t *par)
 695 {
 696   Double_t lMass=0;
 697   lMass = 1.5318; //Xi* mass
 698
 699
 700   Double_t ldNdy  = par[0]; // dN/dy
 701   Double_t l2pi   = 2*TMath::Pi(); // 2pi (cancels in return statement)
 702   Double_t lTemp = par[1]; // Temperature
 703   Double_t lPower = par[2]; // power=n
 704
 705   Double_t lBigCoef = ((lPower-1)*(lPower-2)) / (l2pi*lPower*lTemp*(lPower*lTemp+lMass*(lPower-2)));
 706   Double_t lInPower = 1 + (TMath::Sqrt(x[0]*x[0]+lMass*lMass)-lMass) / (lPower*lTemp);
 707
 708   return l2pi * ldNdy * x[0] * lBigCoef * TMath::Power(lInPower,(-1)*lPower);
 709 }
 710 //________________________________________________________________________
 711 double myMtExpFunc(Double_t *x, Double_t *par)
 712 {
 713   Double_t lMass=0;
 714   lMass = 1.5318; //Xi* mass
 715
 716   return 2*TMath::Pi()*x[0]*par[0]*exp(-TMath::Sqrt(x[0]*x[0]+lMass*lMass)/par[1]);
 717 }
 718 //________________________________________________________________________
 719 double myPtExpFunc(Double_t *x, Double_t *par)
 720 {
 721   return 2*TMath::Pi()*x[0]*par[0]*exp(-x[0]/par[1]);
 722 }
 723 //________________________________________________________________________
 724 double myBoltzFunc(Double_t *x, Double_t *par)
 725 {
 726   Double_t lMass=0;
 727   lMass = 1.5318; //Xi* mass
 728   Double_t Mt = sqrt(pow(x[0],2)+pow(lMass,2));
 729   return par[0]*x[0]*Mt*exp(-Mt/par[1]);
 730
 731 }
 732 //________________________________________________________________________
 733 double myPowerLawPtFunc(Double_t *x, Double_t *par)
 734 {
 735   return ( x[0]*par[0]*pow(1+x[0]/par[1], -par[2]) );
 736 }
 737 //________________________________________________________________________
 738 double myPowerLawFunc(Double_t *x, Double_t *par)
 739 {
 740   return ( par[0]*pow(1+x[0]/par[1], -par[2]) );
 741 }
 742 //________________________________________________________________________
 743 double myBoltsBWFunc(Double_t *x, Double_t *par)
 744 {
 745   double pT = x[0];
 746   double mass    = par[0];
 747   double beta    = par[1];
 748   double temp    = par[2];
 749   double n       = par[3];
 750   double norm    = par[4];
 751
 752   static TF1 * fIntBG = 0;
 753   if(!fIntBG)
 754     fIntBG = new TF1 ("fIntBG", IntegrandBG, 0, 1, 5);
 755
 756   fIntBG->SetParameters(mass, pT, beta, temp,n);
 757   double result = fIntBG->Integral(0,1);
 758   return result*norm;
 759 }
 760 //________________________________________________________________________
 761 double IntegrandBG(const double * x, const double* par){
 762   double x0 = x[0];
 763
 764   double mass     = par[0];
 765   double pT       = par[1];
 766   double beta_max = par[2];
 767   double temp     = par[3];
 768   Double_t n      = par[4];
 769
 770   // Keep beta within reasonable limits
 771   Double_t beta = beta_max * TMath::Power(x0, n);
 772   if (beta > 0.9999999999999999) beta = 0.9999999999999999;
 773
 774   double mT      = TMath::Sqrt(mass*mass+pT*pT);
 775
 776   double rho0   = TMath::ATanH(beta);
 777   double arg00 = pT*TMath::SinH(rho0)/temp;
 778   if (arg00 > 700.) arg00 = 700.; // avoid FPE
 779   double arg01 = mT*TMath::CosH(rho0)/temp;
 780   double f0 = x0*mT*TMath::BesselI0(arg00)*TMath::BesselK1(arg01);
 781
 782   //  printf("r=%f, pt=%f, beta_max=%f, temp=%f, n=%f, mt=%f, beta=%f, rho=%f, argI0=%f, argK1=%f\n", x0, pT, beta_max, temp, n, mT, beta, rho0, arg00, arg01);
 783
 784   return f0;
 785 }