PWG3/vertexingHF/macros/CombineFeedDownMCSubtractionMethodsUncertainties.C

   1 #include "TFile.h"
   2 #include "TH1.h"
   3 #include "TH1D.h"
   4 #include "TH2.h"
   5 #include "TH2F.h"
   6 #include "TGraphAsymmErrors.h"
   7 #include "TCanvas.h"
   8 #include "TLegend.h"
   9 #include "TMath.h"
  10 #include "TROOT.h"
  11 #include "TStyle.h"
  12
  13 #include "AliHFSystErr.h"
  14
  15 #include <Riostream.h>
  16
  17 //_________________________________________________________________________________________
  18 //
  19 //  Macro to combine the the MonteCarlo B feed-down subtraction uncertainties
  20 //
  21 //   Take as input the output files from the HFPtSpectrum class
  22 //    from both fc & Nb subtraction methods and combine the uncertainties.
  23 //   The final central value is set as the one from the Nb-method.
  24 //   The final uncertainties are defined as the envelope of both fc & Nb
  25 //      uncertainties with respect to the new central-value.
  26 //   The final global uncertainties are also defined and a preliminary drawing done.
  27 //
  28 //
  29 //   Usage parameters:
  30 //      1. HFPtSpectrum fc subtraction file
  31 //      2. HFPtSpectrum Nb subtraction file
  32 //      3. Output file name
  33 //      4. FONLL theoretical predictions file to draw on top
  34 //      5. Decay channel as defined in the AliHFSystErr class
  35 //
  36 //_________________________________________________________________________________________
  37
  38 enum centrality{ kpp7, kpp276, k010, k020, k2040, k4060, k6080, k4080, k80100 };
  39
  40 void CombineFeedDownMCSubtractionMethodsUncertainties(const char *fcfilename="HFPtSpectrum_D0Kpi_method1_221110_newnorm.root",
  41                                                       const char *nbfilename="HFPtSpectrum_D0Kpi_method2_221110_newnorm.root",
  42                                                       const char *outfilename="HFPtSpectrum_D0Kpi_combinedFD.root",
  43                                                       const char *thfilename="D0DplusDstarPredictions_y05.root",
  44                                                       Int_t decay=1, Int_t centrality=kpp7)
  45 {
  46
  47   //
  48   // Get fc file inputs
  49   TFile * fcfile = new TFile(fcfilename,"read");
  50   TH1D * histoSigmaCorrFc = (TH1D*)fcfile->Get("histoSigmaCorr");
  51   histoSigmaCorrFc->SetNameTitle("histoSigmaCorrFc","histoSigmaCorrFc");
  52   TGraphAsymmErrors * gSigmaCorrFc = (TGraphAsymmErrors*)fcfile->Get("gSigmaCorr");
  53   gSigmaCorrFc->SetNameTitle("gSigmaCorrFc","gSigmaCorrFc");
  54   TGraphAsymmErrors * gSigmaCorrConservativeFc = (TGraphAsymmErrors*)fcfile->Get("gSigmaCorrConservative");
  55   gSigmaCorrConservativeFc->SetNameTitle("gSigmaCorrConservativeFc","Cross section (fc prompt fraction)");
  56   TGraphAsymmErrors * gFcConservativeFc = (TGraphAsymmErrors*)fcfile->Get("gFcConservative");
  57   gFcConservativeFc->SetNameTitle("gFcConservativeFc","fc prompt fraction");
  58
  59   //
  60   // Get Nb file inputs
  61   TFile * nbfile = new TFile(nbfilename,"read");
  62   TH1D * histoSigmaCorrNb = (TH1D*)nbfile->Get("histoSigmaCorr");
  63   histoSigmaCorrNb->SetNameTitle("histoSigmaCorrNb","histoSigmaCorrNb");
  64   TGraphAsymmErrors * gSigmaCorrNb = (TGraphAsymmErrors*)nbfile->Get("gSigmaCorr");
  65   gSigmaCorrNb->SetNameTitle("gSigmaCorrNb","gSigmaCorrNb");
  66   TGraphAsymmErrors * gSigmaCorrConservativeNb = (TGraphAsymmErrors*)nbfile->Get("gSigmaCorrConservative");
  67   gSigmaCorrConservativeNb->SetNameTitle("gSigmaCorrConservativeNb","Cross section (Nb prompt fraction)");
  68   TGraphAsymmErrors * gFcConservativeNb = (TGraphAsymmErrors*)nbfile->Get("gFcConservative");
  69   gFcConservativeNb->SetNameTitle("gFcConservativeNb","Nb prompt fraction");
  70
  71   //
  72   // Get the predictions input
  73   TFile *thfile = new TFile(thfilename,"read");
  74   TGraphAsymmErrors * thD0KpifromBprediction = (TGraphAsymmErrors*)thfile->Get("D0Kpiprediction");
  75   TGraphAsymmErrors * thDpluskpipiprediction = (TGraphAsymmErrors*)thfile->Get("Dpluskpipiprediction");
  76   TGraphAsymmErrors * thDstarD0piprediction = (TGraphAsymmErrors*)thfile->Get("DstarD0piprediction");
  77   thD0KpifromBprediction->SetLineColor(4);
  78   thD0KpifromBprediction->SetFillColor(kAzure+9);
  79   thDpluskpipiprediction->SetLineColor(4);
  80   thDpluskpipiprediction->SetFillColor(kAzure+9);
  81   thDstarD0piprediction->SetLineColor(4);
  82   thDstarD0piprediction->SetFillColor(kAzure+9);
  83
  84   //
  85   // Get the spectra bins & limits
  86   Int_t nbins = histoSigmaCorrFc->GetNbinsX();
  87   Double_t *limits = new Double_t[nbins+1];
  88   Double_t xlow=0., binwidth=0.;
  89   for (Int_t i=1; i<=nbins; i++) {
  90     binwidth = histoSigmaCorrFc->GetBinWidth(i);
  91     xlow = histoSigmaCorrFc->GetBinLowEdge(i);
  92     limits[i-1] = xlow;
  93   }
  94   limits[nbins] = xlow + binwidth;
  95
  96
  97   //
  98   // Define a new histogram with the real-data reconstructed spectrum binning
  99   //   they will be filled with central value equal to the Nb result
 100   //   and uncertainties taken from the envelope of the result uncertainties
 101   //   The systematical unc. (but FD) will also be re-calculated
 102   //
 103   TH1D * histoSigmaCorr = new TH1D("histoSigmaCorr","corrected cross-section (combined fc and Nb MC feed-down subtraction)",nbins,limits);
 104   TGraphAsymmErrors * gSigmaCorr = new TGraphAsymmErrors(nbins+1);
 105   gSigmaCorr->SetNameTitle("gSigmaCorr","gSigmaCorr (combined fc and Nb MC FD)");
 106   TGraphAsymmErrors * gFcCorrConservative = new TGraphAsymmErrors(nbins+1);
 107   gFcCorrConservative->SetNameTitle("gFcCorrConservative","Combined prompt fraction");
 108   TGraphAsymmErrors * gSigmaCorrConservative = new TGraphAsymmErrors(nbins+1);
 109   gSigmaCorrConservative->SetNameTitle("gSigmaCorrConservative","Cross section (combined prompt fraction)");
 110   TGraphAsymmErrors * gSigmaCorrConservativePC = new TGraphAsymmErrors(nbins+1);
 111   gSigmaCorrConservativePC->SetNameTitle("gSigmaCorrConservativePC","Conservative gSigmaCorr (combined fc and Nb MC FD) in percentages [for drawing with AliHFSystErr]");
 112
 113   //
 114   // Call the systematics uncertainty class for a given decay
 115   //   will help to compute the systematical unc. (but FD)
 116   AliHFSystErr systematics;
 117   if( centrality==kpp276 ) {
 118     systematics.SetIsLowEnergy(true);
 119   } else if( centrality!=kpp7 )  {
 120     systematics.SetCollisionType(1);
 121     if ( centrality == k020 ) {
 122       systematics.SetCentrality("020");
 123     }
 124     else if ( centrality == k4080 ) {
 125       systematics.SetCentrality("4080");
 126     }
 127     else {
 128       cout << " Systematics not yet implemented " << endl;
 129       return;
 130     }
 131   }
 132   else { systematics.SetCollisionType(0); }
 133   systematics.Init(decay);
 134
 135   //
 136   // Loop on all the bins to do the calculations
 137   //
 138   Double_t pt=0., average = 0., averageStatUnc=0., avErrx=0., avErryl=0., avErryh=0., avErryfdl=0., avErryfdh=0.;
 139   Double_t avErrylPC=0., avErryhPC=0., avErryfdlPC=0., avErryfdhPC=0.;
 140   Double_t valFc = 0., valFcErrstat=0., valFcErrx=0., valFcErryl=0., valFcErryh=0., valFcErryfdl=0., valFcErryfdh=0.;
 141   Double_t valNb = 0., valNbErrstat=0., valNbErrx=0., valNbErryl=0., valNbErryh=0., valNbErryfdl=0., valNbErryfdh=0.;
 142   Double_t corrfd = 0., corrfdl=0., corrfdh=0.;
 143   //
 144   for(Int_t ibin=1; ibin<=nbins; ibin++){
 145
 146     // Get input values from fc method
 147     valFc = histoSigmaCorrFc->GetBinContent(ibin);
 148     pt = histoSigmaCorrFc->GetBinCenter(ibin);
 149     valFcErrstat = histoSigmaCorrFc->GetBinError(ibin);
 150     Double_t value =0., ptt=0.;
 151     gSigmaCorrConservativeFc->GetPoint(ibin,ptt,value);
 152     if (value<=0.) continue;
 153     if ( TMath::Abs(valFc-value)>0.1 || TMath::Abs(pt-ptt)>0.1 )
 154       cout << "Hey you ! There might be a problem with the fc input file, please, have a look !" << endl;
 155     valFcErrx = gSigmaCorrFc->GetErrorXlow(ibin);
 156     valFcErryl = gSigmaCorrFc->GetErrorYlow(ibin);
 157     valFcErryh = gSigmaCorrFc->GetErrorYhigh(ibin);
 158     valFcErryfdl = TMath::Abs( gSigmaCorrConservativeFc->GetErrorYlow(ibin) );
 159     valFcErryfdh = TMath::Abs( gSigmaCorrConservativeFc->GetErrorYhigh(ibin) );
 160     Double_t valfdFc = 0., x=0.;
 161     gFcConservativeFc->GetPoint(ibin,x,valfdFc);
 162     Double_t valfdFch = gFcConservativeFc->GetErrorYhigh(ibin);
 163     Double_t valfdFcl = gFcConservativeFc->GetErrorYlow(ibin);
 164
 165     // Get input values from Nb method
 166     valNb = histoSigmaCorrNb->GetBinContent(ibin);
 167     pt = histoSigmaCorrNb->GetBinCenter(ibin);
 168     valNbErrstat = histoSigmaCorrNb->GetBinError(ibin);
 169     gSigmaCorrConservativeNb->GetPoint(ibin,ptt,value);
 170     if ( TMath::Abs(valNb-value)>0.1 || TMath::Abs(pt-ptt)>0.1 )
 171       cout << "Hey you ! There might be a problem with the Nb input file, please, have a look !" << endl;
 172     valNbErrx = gSigmaCorrNb->GetErrorXlow(ibin);
 173     valNbErryl = gSigmaCorrNb->GetErrorYlow(ibin);
 174     valNbErryh = gSigmaCorrNb->GetErrorYhigh(ibin);
 175     valNbErryfdl = gSigmaCorrConservativeNb->GetErrorYlow(ibin);
 176     valNbErryfdh = gSigmaCorrConservativeNb->GetErrorYhigh(ibin);
 177     Double_t valfdNb = 0.;
 178     gFcConservativeNb->GetPoint(ibin,x,valfdNb);
 179     Double_t valfdNbh = gFcConservativeNb->GetErrorYhigh(ibin);
 180     Double_t valfdNbl = gFcConservativeNb->GetErrorYlow(ibin);
 181
 182
 183     // Compute the FD combined value
 184     //    average = valNb
 185     average = valNb ;
 186     corrfd = valfdNb;
 187     avErrx = valFcErrx;
 188     if ( TMath::Abs( valFcErrx - valNbErrx ) > 0.1 )
 189       cout << "Hey you ! There might be consistency problem with the fc & Nb input files, please, have a look !" << endl;
 190     averageStatUnc = valNbErrstat ;
 191 //     cout << " pt=" << pt << ", average="<<average<<endl;
 192 //     cout << "   stat unc (pc)=" << averageStatUnc/average << ", stat-fc (pc)="<<(valFcErrstat/valFc) << ", stat-Nb (pc)="<<(valNbErrstat/valNb)<<endl;
 193
 194     // now estimate the new feed-down combined uncertainties
 195     Double_t minimum[2] = { (valFc - valFcErryfdl), (valNb - valNbErryfdl) };
 196     Double_t maximum[2] = { (valFc + valFcErryfdh), (valNb + valNbErryfdh) };
 197     avErryfdl = average - TMath::MinElement(2,minimum);
 198     avErryfdh = TMath::MaxElement(2,maximum) - average;
 199     avErryfdlPC = avErryfdl / average ; // in percentage
 200     avErryfdhPC = avErryfdh / average ; // in percentage
 201 //     cout << " fc : val " << valFc << " + " << valFcErryfdh <<" - " << valFcErryfdl <<endl;
 202 //     cout << " Nb : val " << valNb << " + " << valNbErryfdh <<" - " << valNbErryfdl <<endl;
 203 //     cout << " fc  & Nb: val " << average << " + " << avErryfdh <<" - " << avErryfdl <<endl;
 204     Double_t minimumfc[2] = { (valfdNb - valfdNbl), (valfdFc - valfdFcl) };
 205     Double_t maximumfc[2] = { (valfdNb + valfdNbh), (valfdFc + valfdFch) };
 206     corrfdl = corrfd - TMath::MinElement(2,minimumfc);
 207     corrfdh = TMath::MaxElement(2,maximumfc) - corrfd;
 208
 209
 210     // compute the global systematics
 211     avErrylPC = systematics.GetTotalSystErr(pt,avErryfdlPC); // in percentage
 212     avErryhPC = systematics.GetTotalSystErr(pt,avErryfdhPC); // in percentage
 213     avErryl = avErrylPC * average ;
 214     avErryh = avErryhPC * average ;
 215 //     cout << "   syst av-l="<<avErryl<<", av-h="<<avErryh<<endl;
 216 //     cout << "     fd-l-pc="<<avErryfdlPC<<", fd-h-pc="<<avErryfdhPC<<", syst err(no fd)-pc="<<systematics.GetTotalSystErr(pt)<<", av-l-pc="<<avErrylPC<<", av-h-pc="<<avErryhPC<<endl;
 217
 218     // fill in the histos and TGraphs
 219     //   fill them only when for non empty bins
 220     if ( average > 0.1 ) {
 221       histoSigmaCorr->SetBinContent(ibin,average);
 222       histoSigmaCorr->SetBinError(ibin,averageStatUnc);
 223       gSigmaCorr->SetPoint(ibin,pt,average);
 224       gSigmaCorr->SetPointError(ibin,valFcErrx,valFcErrx,avErryl,avErryh);
 225       gSigmaCorrConservative->SetPoint(ibin,pt,average);
 226       gSigmaCorrConservative->SetPointError(ibin,valFcErrx,valFcErrx,avErryfdl,avErryfdh);
 227       gSigmaCorrConservativePC->SetPoint(ibin,pt,0.);
 228       gSigmaCorrConservativePC->SetPointError(ibin,valFcErrx,valFcErrx,avErryfdlPC,avErryfdhPC);
 229       gFcCorrConservative->SetPoint(ibin,pt,corrfd);
 230       gFcCorrConservative->SetPointError(ibin,valFcErrx,valFcErrx,corrfdl,corrfdh);
 231     }
 232
 233   }
 234
 235
 236   gROOT->SetStyle("Plain");
 237   gStyle->SetOptTitle(0);
 238
 239   //
 240   // Plot the results
 241   TH2F *histo2Draw = new TH2F("histo2Draw","histo2 (for drawing)",100,0,20.,100,1e3,5e7);
 242   histo2Draw->SetStats(0);
 243   histo2Draw->GetXaxis()->SetTitle("p_{T}  [GeV]");
 244   histo2Draw->GetXaxis()->SetTitleSize(0.05);
 245   histo2Draw->GetXaxis()->SetTitleOffset(0.95);
 246   histo2Draw->GetYaxis()->SetTitle("#frac{1}{BR} #times #frac{d#sigma}{dp_{T}} |_{|y|<0.5}");
 247   histo2Draw->GetYaxis()->SetTitleSize(0.05);
 248   //
 249   TCanvas *combinefdunc = new TCanvas("combinefdunc","show the FD results combination");
 250   //
 251   histo2Draw->Draw();
 252   //
 253   histoSigmaCorrFc->SetMarkerStyle(20);
 254   histoSigmaCorrFc->SetMarkerColor(kGreen+2);
 255   histoSigmaCorrFc->SetLineColor(kGreen+2);
 256   histoSigmaCorrFc->Draw("esame");
 257   gSigmaCorrConservativeFc->SetMarkerStyle(20);
 258   gSigmaCorrConservativeFc->SetMarkerColor(kGreen+2);
 259   gSigmaCorrConservativeFc->SetLineColor(kGreen+2);
 260   gSigmaCorrConservativeFc->SetFillStyle(3004);//2);
 261   gSigmaCorrConservativeFc->SetFillColor(kGreen);
 262   gSigmaCorrConservativeFc->Draw("2[]same");
 263   //
 264   histoSigmaCorrNb->SetMarkerStyle(25);
 265   histoSigmaCorrNb->SetMarkerColor(kViolet+5);
 266   histoSigmaCorrNb->SetLineColor(kViolet+5);
 267   histoSigmaCorrNb->Draw("esame");
 268   gSigmaCorrConservativeNb->SetMarkerStyle(25);
 269   gSigmaCorrConservativeNb->SetMarkerColor(kOrange+7);//kViolet+5);
 270   gSigmaCorrConservativeNb->SetLineColor(kOrange+7);//kOrange+7);//kViolet+5);
 271   gSigmaCorrConservativeNb->SetFillStyle(3018);//02);
 272   gSigmaCorrConservativeNb->SetFillColor(kMagenta);
 273   gSigmaCorrConservativeNb->Draw("2[]same");
 274   //
 275   gSigmaCorrConservative->SetLineColor(kRed);
 276   gSigmaCorrConservative->SetLineWidth(2);
 277   gSigmaCorrConservative->SetFillColor(kRed);
 278   gSigmaCorrConservative->SetFillStyle(0);
 279   gSigmaCorrConservative->Draw("2");
 280   histoSigmaCorr->SetMarkerColor(kRed);
 281   histoSigmaCorr->Draw("esame");
 282   //
 283   //
 284   TLegend* leg=combinefdunc->BuildLegend();
 285   leg->SetFillStyle(0);
 286   combinefdunc->SetLogy();
 287   combinefdunc->Update();
 288
 289   TCanvas *combinefcunc = new TCanvas("combinefcunc","show the fc FD results combination");
 290   //
 291   TH2F *histo3Draw = new TH2F("histo3Draw","histo3 (for drawing)",100,0,20.,10,0.,1.);
 292   histo3Draw->SetStats(0);
 293   histo3Draw->GetXaxis()->SetTitle("p_{T}  [GeV]");
 294   histo3Draw->GetXaxis()->SetTitleSize(0.05);
 295   histo3Draw->GetXaxis()->SetTitleOffset(0.95);
 296   histo3Draw->GetYaxis()->SetTitle("Prompt fraction of the raw yields");
 297   histo3Draw->GetYaxis()->SetTitleSize(0.05);
 298   histo3Draw->Draw();
 299   //
 300   gFcConservativeFc->SetMarkerStyle(20);
 301   gFcConservativeFc->SetMarkerColor(kGreen+2);
 302   gFcConservativeFc->SetLineColor(kGreen+2);
 303   gFcConservativeFc->SetFillStyle(3004);
 304   gFcConservativeFc->SetFillColor(kGreen);
 305   gFcConservativeFc->Draw("2P");
 306   //
 307   gFcConservativeNb ->SetMarkerStyle(25);
 308   gFcConservativeNb ->SetMarkerSize(1.3);
 309   gFcConservativeNb->SetMarkerColor(kOrange+7);//kViolet+5);
 310   gFcConservativeNb->SetLineColor(kOrange+7);//kViolet+5);
 311   gFcConservativeNb->SetFillStyle(3018);
 312   gFcConservativeNb->SetFillColor(kMagenta);
 313   gFcConservativeNb->Draw("2P");
 314   //
 315   gFcCorrConservative->SetMarkerStyle(21);
 316   gFcCorrConservative->SetLineColor(kRed);
 317   gFcCorrConservative->SetLineWidth(2);
 318   gFcCorrConservative->SetFillColor(kRed);
 319   gFcCorrConservative->SetFillStyle(0);
 320   gFcCorrConservative->Draw("2P");
 321   //
 322   leg=combinefcunc->BuildLegend();
 323   leg->SetFillStyle(0);
 324   //
 325   combinefcunc->Update();
 326
 327   //
 328   // Plot the results
 329   TCanvas *finalresults = new TCanvas("finalresults","show all combined results");
 330   //
 331   if ( decay==1 ) {
 332     thD0KpifromBprediction->SetLineColor(kGreen+2);
 333     thD0KpifromBprediction->SetLineWidth(3);
 334     thD0KpifromBprediction->SetFillColor(kGreen-6);
 335     thD0KpifromBprediction->Draw("3CA");
 336     thD0KpifromBprediction->Draw("CX");
 337   }
 338   else if ( decay==2 ) {
 339     thDpluskpipiprediction->SetLineColor(kGreen+2);
 340     thDpluskpipiprediction->SetLineWidth(3);
 341     thDpluskpipiprediction->SetFillColor(kGreen-6);
 342     thDpluskpipiprediction->Draw("3CA");
 343     thDpluskpipiprediction->Draw("CX");
 344   }
 345   else if ( decay==3 ) {
 346     thDstarD0piprediction->SetLineColor(kGreen+2);
 347     thDstarD0piprediction->SetLineWidth(3);
 348     thDstarD0piprediction->SetFillColor(kGreen-6);
 349     thDstarD0piprediction->Draw("3CA");
 350     thDstarD0piprediction->Draw("CX");
 351   }
 352   //
 353   gSigmaCorr->SetLineColor(kRed);
 354   gSigmaCorr->SetLineWidth(1);
 355   gSigmaCorr->SetFillColor(kRed);
 356   gSigmaCorr->SetFillStyle(0);
 357   gSigmaCorr->Draw("2");
 358   histoSigmaCorr->SetMarkerStyle(21);
 359   histoSigmaCorr->SetMarkerColor(kRed);
 360   histoSigmaCorr->Draw("esame");
 361   //
 362   leg = new TLegend(0.7,0.75,0.87,0.5);
 363   leg->SetBorderSize(0);
 364   leg->SetLineColor(0);
 365   leg->SetFillColor(0);
 366   leg->SetTextFont(42);
 367   if ( decay==1 ) leg->AddEntry(thD0KpifromBprediction,"FONLL ","fl");
 368   else if ( decay==2 ) leg->AddEntry(thDpluskpipiprediction,"FONLL ","fl");
 369   else if ( decay==3 ) leg->AddEntry(thDstarD0piprediction,"FONLL ","fl");
 370   leg->AddEntry(histoSigmaCorr,"data stat. unc.","pl");
 371   leg->AddEntry(gSigmaCorr,"data syst. unc.","f");
 372   leg->Draw();
 373   //
 374   finalresults->SetLogy();
 375   finalresults->Update();
 376
 377
 378   //
 379   // Draw all the systematics independently
 380   systematics.DrawErrors(gSigmaCorrConservativePC);
 381
 382
 383   // Write the output to a file
 384   TFile * out = new TFile(outfilename,"recreate");
 385   histoSigmaCorr->Write();
 386   gSigmaCorr->Write();
 387   gSigmaCorrConservative->Write();
 388   gSigmaCorrConservativePC->Write();
 389   gFcCorrConservative->Write();
 390   out->Write();
 391
 392 }