PWGGA/PHOSTasks/PHOS_pp_pi0/macros/MakeFinalSpectrum.C

   1 void MakeFinalSpectrum()
   2 {
   3   //-----------------------------------------------------------------------------
   4   // This macro takes the raw pi0 spectrum from LHC11a_FitResult_20130314.root,
   5   // correct it for feed down, then for efficiency,
   6   // adds all systematic errors and produces the production invariant spectrum
   7   //--
   8   // Last modification: 06.07.2012 by Yuri Kharlov
   9   //-----------------------------------------------------------------------------
  10
  11   TFile *f  = new TFile("LHC11a_FitResult_20130913.root");
  12   TH1D * nr1CB    = (TH1D*)f->Get("Mix_CB_yr1")  ;
  13   TH1D * nr1intCB = (TH1D*)f->Get("Mix_CB_yr1int")  ;
  14   TH1D * nr2CB    = (TH1D*)f->Get("Mix_CB_yr2")  ;
  15   TH1D * nr2intCB = (TH1D*)f->Get("Mix_CB_yr2int")  ;
  16
  17   TH1D * nr1GS    = (TH1D*)f->Get("Mix_yr1")  ;
  18   TH1D * nr1intGS = (TH1D*)f->Get("Mix_yr1int")  ;
  19   TH1D * nr2GS    = (TH1D*)f->Get("Mix_yr2")  ;
  20   TH1D * nr2intGS = (TH1D*)f->Get("Mix_yr2int")  ;
  21
  22   //Divide by bin width
  23   for(Int_t i=1;i<= nr1CB->GetNbinsX();i++){
  24     nr1CB   ->SetBinContent(i,nr1CB->GetBinContent(i)/nr1CB->GetXaxis()->GetBinWidth(i)) ;
  25     nr1CB   ->SetBinError  (i,nr1CB->GetBinError(i)/nr1CB->GetXaxis()->GetBinWidth(i)) ;
  26     nr1intCB->SetBinContent(i,nr1intCB->GetBinContent(i)/nr1intCB->GetXaxis()->GetBinWidth(i)) ;
  27     nr1intCB->SetBinError  (i,nr1intCB->GetBinError(i)/nr1intCB->GetXaxis()->GetBinWidth(i)) ;
  28     nr2CB   ->SetBinContent(i,nr2CB->GetBinContent(i)/nr2CB->GetXaxis()->GetBinWidth(i)) ;
  29     nr2CB   ->SetBinError  (i,nr2CB->GetBinError(i)/nr2CB->GetXaxis()->GetBinWidth(i)) ;
  30     nr2intCB->SetBinContent(i,nr2intCB->GetBinContent(i)/nr2intCB->GetXaxis()->GetBinWidth(i)) ;
  31     nr2intCB->SetBinError  (i,nr2intCB->GetBinError(i)/nr2intCB->GetXaxis()->GetBinWidth(i)) ;
  32
  33     nr1GS   ->SetBinContent(i,nr1GS->GetBinContent(i)/nr1GS->GetXaxis()->GetBinWidth(i)) ;
  34     nr1GS   ->SetBinError  (i,nr1GS->GetBinError(i)/nr1GS->GetXaxis()->GetBinWidth(i)) ;
  35     nr1intGS->SetBinContent(i,nr1intGS->GetBinContent(i)/nr1intGS->GetXaxis()->GetBinWidth(i)) ;
  36     nr1intGS->SetBinError  (i,nr1intGS->GetBinError(i)/nr1intGS->GetXaxis()->GetBinWidth(i)) ;
  37     nr2GS   ->SetBinContent(i,nr2GS->GetBinContent(i)/nr2GS->GetXaxis()->GetBinWidth(i)) ;
  38     nr2GS   ->SetBinError  (i,nr2GS->GetBinError(i)/nr2GS->GetXaxis()->GetBinWidth(i)) ;
  39     nr2intGS->SetBinContent(i,nr2intGS->GetBinContent(i)/nr2intGS->GetXaxis()->GetBinWidth(i)) ;
  40     nr2intGS->SetBinError  (i,nr2intGS->GetBinError(i)/nr2intGS->GetXaxis()->GetBinWidth(i)) ;
  41   }
  42
  43   // feed down correction
  44
  45   TF1 *fKaonContaminationToPi0 = new TF1("kaonCont","1./(1.-1.33*1.2*exp(-2.95-0.16*x))",0.,30.);
  46   nr1GS   ->Divide(fKaonContaminationToPi0);
  47   nr1intGS->Divide(fKaonContaminationToPi0);
  48   nr2GS   ->Divide(fKaonContaminationToPi0);
  49   nr2intGS->Divide(fKaonContaminationToPi0);
  50   nr1CB   ->Divide(fKaonContaminationToPi0);
  51   nr2CB   ->Divide(fKaonContaminationToPi0);
  52   nr1intCB->Divide(fKaonContaminationToPi0);
  53   nr2intCB->Divide(fKaonContaminationToPi0);
  54
  55   // SPD pileup correction
  56
  57   TF1 *fSPDpileup = new TF1("SPDpileup","0.988",0.,30.);
  58   nr1GS   ->Multiply(fSPDpileup);
  59   nr1intGS->Multiply(fSPDpileup);
  60   nr2GS   ->Multiply(fSPDpileup);
  61   nr2intGS->Multiply(fSPDpileup);
  62   nr1CB   ->Multiply(fSPDpileup);
  63   nr2CB   ->Multiply(fSPDpileup);
  64   nr1intCB->Multiply(fSPDpileup);
  65   nr2intCB->Multiply(fSPDpileup);
  66
  67   //correct for efficiency
  68   TFile *fEff = new TFile("Pi0_efficiency_LHC11a__20131029_Mall.root") ;
  69
  70   TF1 * effGS=fEff->Get("eff_Pi0_Gaus_2760GeV") ;
  71   TF1 * effCB=fEff->Get("eff_Pi0_CB_2760GeV") ;
  72   effGS   ->SetRange(0.,25.) ;
  73   effCB   ->SetRange(0.,25.) ;
  74
  75   nr1GS   ->Divide(effGS) ;
  76   nr1intGS->Divide(effGS) ;
  77   nr2GS   ->Divide(effGS) ;
  78   nr2intGS->Divide(effGS) ;
  79   nr1CB   ->Divide(effCB) ;
  80   nr2CB   ->Divide(effCB) ;
  81   nr1intCB->Divide(effCB) ;
  82   nr2intCB->Divide(effCB) ;
  83
  84   //make 1/pt
  85   for(Int_t i=1;i<=nr1CB->GetNbinsX();i++){
  86     Double_t pt = TMath::TwoPi()*nr1CB->GetXaxis()->GetBinCenter(i);
  87     nr1CB   ->SetBinContent(i,nr1CB   ->GetBinContent(i)/pt) ;
  88     nr1CB   ->SetBinError  (i,nr1CB   ->GetBinError(i)  /pt) ;
  89     nr1intCB->SetBinContent(i,nr1intCB->GetBinContent(i)/pt) ;
  90     nr1intCB->SetBinError  (i,nr1intCB->GetBinError(i)  /pt) ;
  91     nr2CB   ->SetBinContent(i,nr2CB   ->GetBinContent(i)/pt) ;
  92     nr2CB   ->SetBinError  (i,nr2CB   ->GetBinError(i)  /pt) ;
  93     nr2intCB->SetBinContent(i,nr2intCB->GetBinContent(i)/pt) ;
  94     nr2intCB->SetBinError  (i,nr2intCB->GetBinError(i)  /pt) ;
  95
  96     nr1GS   ->SetBinContent(i,nr1GS   ->GetBinContent(i)/pt) ;
  97     nr1GS   ->SetBinError  (i,nr1GS   ->GetBinError(i)  /pt) ;
  98     nr1intGS->SetBinContent(i,nr1intGS->GetBinContent(i)/pt) ;
  99     nr1intGS->SetBinError  (i,nr1intGS->GetBinError(i)  /pt) ;
 100     nr2GS   ->SetBinContent(i,nr2GS   ->GetBinContent(i)/pt) ;
 101     nr2GS   ->SetBinError  (i,nr2GS   ->GetBinError(i)  /pt) ;
 102     nr2intGS->SetBinContent(i,nr2intGS->GetBinContent(i)/pt) ;
 103     nr2intGS->SetBinError  (i,nr2intGS->GetBinError(i)  /pt) ;
 104   }
 105
 106   //For the final spectrum we take average of fits
 107   //with numerical integration of entries in signal
 108
 109   TH1D * hStat = (TH1D*)nr2intCB->Clone("hPi02760GeVStat") ;
 110   TH1D * hSys  = (TH1D*)hStat   ->Clone("hPi02760GeVSys") ;
 111   TH1D * hSys2 = (TH1D*)hStat   ->Clone("hPi02760GeVSysTypeB") ;
 112   TH1D * hSys3 = (TH1D*)hStat   ->Clone("hPi02760GeVSysTypeC") ;
 113   hStat->SetAxisRange(0.,14.9,"X");
 114   hSys ->SetAxisRange(0.,14.9,"X");
 115
 116   //For systematic error estimate take largest deviation
 117   //of integrated yeilds (note, they are efficiency corrected)
 118   for(Int_t i=1;i<=nr1CB->GetNbinsX();i++){
 119     Double_t mean= hStat->GetBinContent(i) ;
 120     Double_t dev = TMath::Max(
 121                    TMath::Max(TMath::Abs(nr1intCB->GetBinContent(i)-mean),
 122                               TMath::Abs(nr2intCB->GetBinContent(i)-mean)),
 123                    TMath::Max(TMath::Abs(nr1intGS->GetBinContent(i)-mean),
 124                               TMath::Abs(nr2intGS->GetBinContent(i)-mean))
 125                               );
 126     hSys ->SetBinError(i,dev) ;
 127     hSys2->SetBinError(i,dev) ;
 128   }
 129
 130   //Add other sys errors
 131   TF1 * globalE = new TF1("globalE","1.-((x+1.354)/(x*1.002+1.354))^6.18 ",1.,30.) ;
 132   TF1 * conv    = new TF1("conversion","0.035",0.,30.) ;
 133   TF1 * accept  = new TF1("accept"    ,"0.01" ,0.,30.) ;
 134   TF1 * pileup  = new TF1("pileup"    ,"0.004",0.,30.) ;
 135   TF1 * calib   = new TF1("calib"     ,"0.005",0.,30.) ;
 136   TF1 * modDiff = new TF1("modDiff"   ,"0.04",0.,30.) ;
 137   // TF1 * modDiff = new TF1("modDiff"   ,"16.9*exp(-4.5*x)+0.033",0.,30.) ;
 138   TF1 * tofCut  = new TF1("tofCut"    ,"0.0105" ,0.,30.) ;
 139
 140   //Borya's estimate of non-linearity (found for pp @ 7 TeV)
 141   TF1 * nonlin= new TF1("nl","0.015+7.38*exp(-x/0.24)",0.,30.) ;
 142
 143   //Draw sys errors
 144   TH1D * hRelSysRaw = (TH1D*)hSys->Clone("RelSysRaw") ;
 145   hRelSysRaw->SetTitle("Summary of systematic errors");
 146   for(Int_t i=1;i<=hSys->GetNbinsX();i++){
 147     Double_t mean= hSys->GetBinContent(i) ;
 148     Double_t a=hSys->GetBinError(i) ;
 149     if(mean>0)
 150       hRelSysRaw->SetBinContent(i,a/mean) ;
 151     else
 152       hRelSysRaw->SetBinContent(i,0.) ;
 153       hRelSysRaw->SetBinError(i,0.) ;
 154   }
 155
 156   //Add errors in sys errors
 157
 158   TH1D * hRelSysTot = (TH1D*)hSys->Clone("RelSys") ;
 159   hRelSysTot->SetTitle("Summary of systematic uncertainties");
 160
 161   for(Int_t i=1;i<=hSys->GetNbinsX();i++){
 162     Double_t pt   = hSys->GetXaxis()->GetBinCenter(i) ;
 163     Double_t mean = hSys->GetBinContent(i) ;
 164     Double_t a    = hSys->GetBinError(i) ;
 165     // Double_t b    = mean * hSysErrModules->GetBinContent(i) ;
 166     Double_t tot= mean*mean*nonlin ->Eval(pt)*nonlin ->Eval(pt)
 167                  +mean*mean*conv   ->Eval(pt)*conv   ->Eval(pt)
 168                  +mean*mean*accept ->Eval(pt)*accept ->Eval(pt)
 169                  +mean*mean*pileup ->Eval(pt)*pileup ->Eval(pt)
 170                  +mean*mean*calib  ->Eval(pt)*calib  ->Eval(pt)
 171                  +mean*mean*modDiff->Eval(pt)*modDiff->Eval(pt)
 172                  +mean*mean*tofCut ->Eval(pt)*tofCut ->Eval(pt)
 173                  +mean*mean*globalE->Eval(pt)*globalE->Eval(pt);
 174     Double_t raa= mean*mean*nonlin->Eval(pt)*nonlin->Eval(pt)
 175                  +mean*mean*pileup->Eval(pt)*pileup->Eval(pt)
 176                  +mean*mean*calib ->Eval(pt)*calib ->Eval(pt);
 177     hSys3->SetBinError(i,TMath::Sqrt(tot)) ;
 178     // hSys->SetBinError(i,TMath::Sqrt(tot + a*a + b*b)) ;
 179     hSys2->SetBinError(i,TMath::Sqrt(raa + a*a)) ;
 180     hSys ->SetBinError(i,TMath::Sqrt(tot + a*a)) ;
 181
 182     a = hSys->GetBinError(i) ;
 183     printf("i=%d, %g+-%g\n",i,mean,a);
 184     if(mean>0)
 185       hRelSysTot->SetBinContent(i,a/mean) ;
 186     else {
 187       hRelSysTot->SetBinContent(i,0.) ;
 188       hRelSysTot->SetBinError(i,0.) ;
 189     }
 190   }
 191
 192   // TFile *fSysErrModules = TFile::Open("PHOS_sysErr_modules.root");
 193   // TH1D* hSysErrModules = (TH1D*)fSysErrModules->Get("hSyserr");
 194
 195   gStyle->SetOptStat(0);
 196   TCanvas * c = new TCanvas("c","SysErrors") ;
 197   c->SetLogy();
 198   c->cd() ;
 199   gPad->SetRightMargin(0.02);
 200   gPad->SetTopMargin(0.07);
 201   hRelSysTot->SetAxisRange(0.8    ,11.9,"X");
 202   hRelSysTot->SetAxisRange(0.0031,0.45,"Y");
 203   hRelSysTot->GetYaxis()->SetMoreLogLabels();
 204   hRelSysTot->GetYaxis()->SetNoExponent();
 205   hRelSysTot->SetNdivisions(520,"X");
 206   hRelSysTot->SetLineColor(1) ;
 207   hRelSysTot->SetLineWidth(2) ;
 208   hRelSysRaw->SetLineColor(2) ;
 209   hRelSysRaw->SetLineWidth(2) ;
 210   globalE->SetLineColor(kGreen+2) ;
 211   nonlin ->SetLineColor(4) ;
 212   conv   ->SetLineColor(6) ;
 213   accept ->SetLineColor(kOrange) ;
 214   pileup ->SetLineColor(42);
 215   calib  ->SetLineColor(44);
 216   calib  ->SetLineStyle(2);
 217   modDiff->SetLineColor(52);
 218   modDiff->SetLineStyle(2);
 219   tofCut ->SetLineColor(53);
 220   tofCut ->SetLineStyle(3);
 221   // hSysErrModules->SetLineColor(kOrange+2);
 222   // hSysErrModules->SetLineStyle(2);
 223   hRelSysTot->SetXTitle("p_{T} (GeV/c)") ;
 224   hRelSysTot->SetYTitle("Rel.syst.error") ;
 225   hRelSysTot->GetYaxis()->SetTitleOffset(1.2) ;
 226   hRelSysTot->Draw("hist") ;
 227   hRelSysRaw->Draw("h same") ;
 228   globalE->Draw("same") ;
 229   nonlin ->Draw("same") ;
 230   conv   ->Draw("same") ;
 231   accept ->Draw("same") ;
 232   pileup ->Draw("same");
 233   calib  ->Draw("same");
 234   modDiff->Draw("same");
 235   tofCut ->Draw("same");
 236   // hSysErrModules->Draw("same ][");
 237   TLegend * l = new TLegend(0.57,0.62,0.85,0.925) ;
 238   l->SetFillColor(kWhite);
 239   l->SetBorderSize(0);
 240   l->AddEntry(hRelSysTot,"Total uncertainty","l") ;
 241   l->AddEntry(hRelSysRaw,"Raw extraction","l") ;
 242   l->AddEntry(conv   ,"Conversion","l") ;
 243   l->AddEntry(nonlin ,"Non-linearity","l") ;
 244   l->AddEntry(accept ,"Acceptance","l");
 245   l->AddEntry(pileup ,"Pileup","l");
 246   l->AddEntry(calib  ,"Rel.calib.","l");
 247   l->AddEntry(modDiff,"Per-module yield","l");
 248   l->AddEntry(tofCut ,"Timing cut","l");
 249   // l->AddEntry(hSysErrModules,"Intermodule spectra","l");
 250   l->AddEntry(globalE,"Global E scale","l") ;
 251   l->Draw() ;
 252
 253   c->Print("LHC11a_SysErrors.eps");
 254
 255   hStat->SetTitle("Normalized production #pi^{0} yield, pp @ 2.76 TeV, stat.err.");
 256   hStat->SetXTitle("p_{T}, GeV/c");
 257   hStat->SetYTitle("1/N_{ev} 1/(2#pi p_{T}) d^{2}N/dydp_{T} (GeV^{-2}c^{2})");
 258   hSys ->SetTitle("Normalized production #pi^{0} yield, pp @ 2.76 TeV, total syst.err.");
 259   hSys ->SetXTitle("p_{T}, GeV/c");
 260   hSys ->SetYTitle("1/N_{ev} 1/(2#pi p_{T}) d^{2}N/dydp_{T} (GeV^{-2}c^{2})");
 261   hSys2->SetTitle("Normalized production #pi^{0} yield, pp @ 2.76 TeV, syst.err. for R_{AA}");
 262   hSys3->SetTitle("Normalized production #pi^{0} yield, pp @ 2.76 TeV, apparatus syst.err.");
 263   hSys ->SetFillColor(kBlue-10) ;
 264   hSys ->SetNdivisions(520,"X");
 265
 266   hSysFinal = (TH1F*)hSys->Clone("hSysFinal");
 267   hSysFinal ->SetTitle("Production #pi^{0} yield, pp @ 2.76 TeV, 07.11.2013");
 268   hSysFinal ->GetYaxis()->SetTitleOffset(1.2);
 269   hStat->GetYaxis()->SetTitleOffset(1.2);
 270   hSysFinal ->SetAxisRange(0.8,9.9,"X");
 271
 272   TCanvas *c2 = new TCanvas("c2","Production spectrum");
 273   gPad->SetRightMargin(0.02);
 274   gPad->SetTopMargin(0.07);
 275   gPad->SetGridx();
 276   gPad->SetGridy();
 277   c2->SetLogy();
 278   hSysFinal->Draw("E2") ;
 279   hStat->SetMarkerStyle(20) ;
 280   hStat->SetMarkerColor(4) ;
 281   hStat->SetLineColor(4) ;
 282   hStat->Draw("same") ;
 283
 284   // //Apply bin width correction
 285   // TH1D * hBWcorr = BinWidthCorrection(hStat) ;
 286   // hSys->Divide(hBWcorr) ;
 287   // hSys2->Divide(hBWcorr) ;
 288   // hSys3->Divide(hBWcorr) ;
 289   // hStat->Divide(hBWcorr) ;
 290
 291   c2->Print("LHC11a_pi0Spectrum.eps");
 292
 293   TF1 *tsallis = FitTsallis(hStat,hSys);
 294
 295   TFile fout("PHOS_pi0_2760eV_noBinWidthCorr_20131112.root","recreate") ;
 296   hSys   ->Write() ;
 297   hSys2  ->Write() ;
 298   hSys3  ->Write() ;
 299   hStat  ->Write() ;
 300   tsallis->Write() ;
 301   fout.Close() ;
 302
 303
 304 }
 305
 306 //-----------------------------------------------------------------------------
 307 TF1 *FitTsallis(TH1D* hStat, TH1D* hSys)
 308 {
 309   TF1 * fit = new TF1("Tsalis","[0]/2./3.1415*([2]-1.)*([2]-2.)/([2]*[1]*([2]*[1]+0.135*([2]-2.)))*(1.+(sqrt(x*x+0.135*0.135)-0.135)/([2]*[1]))^-[2]",0.5,25.) ;
 310   fit->SetParameters(10.,0.2,8.) ;
 311   fit->SetLineColor(kBlack) ;
 312   fit->SetLineWidth(2) ;
 313   fit->SetParName(0,"dN/dy") ;
 314   fit->SetParName(1,"T") ;
 315   fit->SetParName(2,"n") ;
 316
 317   TH1D * hSysRatio  = (TH1D*)hSys ->Clone("RatioSys") ;
 318   TH1D * hStatRatio = (TH1D*)hStat->Clone("RatioStat") ;
 319
 320   TH1D * hsum = (TH1D*)hStat->Clone("sum") ;
 321   for(Int_t i=1; i<=hsum->GetNbinsX();i++){
 322     Double_t a=hSys->GetBinError(i) ;
 323     Double_t b=hStat->GetBinError(i) ;
 324     hsum->SetBinError(i,TMath::Sqrt(a*a+b*b)) ;
 325   }
 326   hsum->SetStats(1) ;
 327   hsum->Fit(fit,"Q") ;
 328   Double_t meanPt=fit->Moment(1,1.,10.) ;
 329   printf("<pt>=%f \n",meanPt) ;
 330
 331   hSysRatio ->Divide(fit) ;
 332   hStatRatio->Divide(fit) ;
 333
 334   return fit;
 335 }
 336 //-----------------------------------------------------------------------------
 337 TH1D * BinWidthCorrection(TH1D * h){
 338   //We apply bin width a-la PHENIX
 339   //Use Tsalis fit to calculate shift in y direction
 340
 341   TF1 * fit = new TF1("hag","[0]*(([2]*[1]+sqrt(x*x+0.135*0.135))/([2]*[1]+0.135))^-[2]",0.5,25.) ;
 342   fit->SetParameters(10.,0.2,8.) ;
 343   TCanvas * corr = new TCanvas("BWcorr","Bin width correction") ;
 344   Int_t col[6]={kRed,kOrange,kMagenta,kGreen,kCyan,kBlue} ;
 345   TH1D * hcorr[20] ;
 346   char key[55] ;
 347   Double_t rMax=10 ;
 348   Int_t iteration=0 ;
 349   TH1D * htmp = (TH1D*)h->Clone("tmp") ;
 350   while(iteration<6){
 351     printf(" Iteration %d: rMax=%f \n",iteration, rMax) ;
 352     htmp->Fit(fit,"N") ;
 353     sprintf(key,"Ineration%d",iteration) ;
 354     hcorr[iteration]=(TH1D*)h->Clone(key);
 355     rMax= 0;
 356     for(Int_t i=1;i<=h->GetNbinsX();i++){
 357       Double_t a=h->GetXaxis()->GetBinLowEdge(i) ;
 358       Double_t b=h->GetXaxis()->GetBinUpEdge(i) ;
 359       Double_t r=fit->Integral(a,b)/(b-a)/fit->Eval(0.5*(a+b)) ;
 360       hcorr[iteration]->SetBinContent(i,r) ;
 361       hcorr[iteration]->SetBinError(i,0.) ;
 362       if(rMax<r)rMax=r ;
 363     }
 364     delete htmp ;
 365     htmp = (TH1D*)h->Clone("tmp") ;
 366     htmp->Divide(hcorr[iteration]) ;
 367     corr->cd() ;
 368     hcorr[iteration]->SetLineColor(col[iteration]);
 369     if(iteration==0)
 370       hcorr[iteration]->Draw() ;
 371     else
 372       hcorr[iteration]->Draw("same") ;
 373     corr->Update() ;
 374     iteration++ ;
 375   }
 376
 377   hcorr[5]->SetTitle("Bin-width correction for #pi^{0} spectrum");
 378   hcorr[5]->SetYTitle("Bin-width corrected / uncorrected");
 379   TFile fout("PHOS_pi0_7TeV_BinWidthCorrection.root","recreate") ;
 380   hcorr[5]->Write();
 381   fout.Close();
 382
 383   return hcorr[5] ;
 384 }