AliPhysics  648edd6 (648edd6)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ComputeDmesonYield.C
Go to the documentation of this file.
1 #if !defined(__CINT__) || defined(__MAKECINT__)
2 #include "TFile.h"
3 #include "TNtuple.h"
4 #include "TString.h"
5 #include "TGraphAsymmErrors.h"
6 #include "TH1.h"
7 #include "TH2.h"
8 #include "TLine.h"
9 #include "TMarker.h"
10 #include "TStyle.h"
11 #include "TLatex.h"
12 #include "TCanvas.h"
13 #include "TMath.h"
14 #include "TLegend.h"
15 #include "AliHFSystErr.h"
16 #endif
17 
18 // input files
19 TString filnamPPref="~/alice/Charm/PbYield/2011/Final/ppref/D0Kpi_276TeV_FONLLExtrapolationAndExtrapolation_from7TeVAlicedata_combinedFD_141211_010312_160712_rebin_1_2_3_4_5_6_8_12_16_24.root";
20 
21 TString filnamSpectrumNb="~/alice/Charm/PbYield/2011/Final/010/Dzero/HFPtSpectrum_D0toKpi_010_Nb_160615.root";
22 TString filnamSpectrumFc="~/alice/Charm/PbYield/2011/Final/010/Dzero/HFPtSpectrum_D0toKpi_010_fc_160615.root";
23 TString filnamRaaNb="~/alice/Charm/PbYield/2011/Final/010/Dzero/HFPtSpectrumRaa_D0toKpi_010_Nb_160615.root";
24 TString filnamRaaFc="~/alice/Charm/PbYield/2011/Final/010/Dzero/HFPtSpectrumRaa_D0toKpi_010_fc_160615.root";
25 
26 //configuration
27 enum dspec{kDzero,kDplus,kDstar,kDs};
28 Int_t mesonSpecie=kDzero;
29 TString centrality="0-10";
30 const Int_t nPtBins=9;
31 Double_t binlim[nPtBins+1]={1.,2.,3.,4.,5.,6.,8.,12.,16.,24.};
32 Int_t method=2; // 1=fc; 2=Nb --> default=2
33 Int_t optErrFD=1; // 0=from histos, not combined;
34  // 1= from ntuple with Rbc hypo, not combining Nb and fc;
35  // 2= from ntuple with Rbc hypo, combining Nb and fc
36  // --> default=2; change to 1 only to produce files for D-meson ratios
37 Float_t centHypoFdOverPr=2.;
38 Float_t lowHypoFdOverPr=1.;
39 Float_t highHypoFdOverPr=3.;
40 Double_t normToCsec=1.; // put here the trigger cross section in ub; if ==1 per-event yield are computed
41 TString collSyst="Pb-Pb";
42 
43 
44 // Graphical styles
45 Bool_t draw[nPtBins]={1,1,0,1,0,1,0,0,1};
46 Int_t colors[nPtBins]={kGray+2,kMagenta+1,kMagenta,kBlue,kCyan,kGreen+2,kYellow+2,kOrange+1,kRed+1};
47 Int_t lstyle[nPtBins]={9,10,3,5,7,1,3,6,2};
48 
49 
50 Bool_t PbPbDataSyst(AliHFSystErr *syst, TH1D* heff, Double_t pt, Double_t &dataSystUp, Double_t &dataSystDown);
51 
52 
53 void ComputeDmesonYield(){
54 
55  TString mesName="Dzero";
56  Int_t mesCode=1;
57  Double_t brat=0.0393;
58  TString mesSymb="D^{0}";
59  Float_t ptForExtrap=36.;
60  if(mesonSpecie==kDplus){
61  mesName="Dplus";
62  mesCode=2;
63  brat=0.0946;
64  mesSymb="D^{+}";
65  ptForExtrap=24.;
66  }else if(mesonSpecie==kDstar){
67  mesName="Dstar";
68  mesCode=3;
69  brat=0.0393*0.677;
70  mesSymb="D^{*+}";
71  ptForExtrap=24.;
72  }else if(mesonSpecie==kDs){
73  mesName="Ds";
74  mesCode=4;
75  brat=0.0227;
76  mesSymb="D_{s}^{+}";
77  if(collSyst=="Pb-Pb"){
78  centHypoFdOverPr=1.;
79  lowHypoFdOverPr=1./3.;
80  highHypoFdOverPr=3.;
81  }
82  ptForExtrap=12.;
83  }
84 
85  TString centralityNoDash=centrality.Data();
86  centralityNoDash.ReplaceAll("-","");
87  TString filnamChi;
88  TString filnamCnt;
89  TString filnamCntS;
90  Int_t optCombineFDEloss=1; // 0=enevlope, 1= sum in quadrature
91  Bool_t correctForBR=kTRUE;
92  Bool_t showRcbSystNorm=kTRUE;
93 
94  if(method==1){
95  filnamChi=filnamSpectrumFc.Data();
96  filnamCnt=filnamRaaFc.Data();
97  filnamCntS=filnamRaaNb.Data();
98  }
99  else if(method==2){
100  filnamChi=filnamSpectrumNb.Data();
101  filnamCnt=filnamRaaNb.Data();
102  filnamCntS=filnamRaaFc.Data();
103  }
104 
105  Int_t colorSystFD=kGray+1;
106  Int_t colorSystRb=kOrange;
107  Int_t colorSystDatappC=kBlue+1;//kMagenta+1;
108  Int_t colorSystDataaaC=kRed+1;//kMagenta+1;
109  Int_t colorppC=kBlue+1;
110  Int_t coloraaC=kRed+1;
111  Int_t linestppC=1;
112  Int_t linestaaC=1;
113  Int_t linewippC=2;
114  Int_t linewiaaC=2;
115  Int_t markerppC=21;
116  Int_t markeraaC=23;
117  Int_t msizppC=1.2;
118  Int_t msizaaC=1.2;
119  Float_t sizesystdata=0.4;
120  Float_t sizesystfd=0.3;
121  Float_t sizesystrb=0.2;
122  Float_t sizesysttot=0.15;
123 
124 
125  // pp reference
126 
127  TFile *filPP=new TFile(filnamPPref.Data());
128  TH1D *hppRef = (TH1D*)filPP->Get("hReference");
129  if(!hppRef) hppRef = (TH1D*)filPP->Get("fhScaledData");
130  TH1D *hppRefSystData = (TH1D*)filPP->Get("fhScaledSystData");
131 
132  Float_t relsysterrLowDatapp[nPtBins],relsysterrHiDatapp[nPtBins];
133  TGraphAsymmErrors* gppRefSyst=(TGraphAsymmErrors*)filPP->Get("gReferenceSyst");
134  if(gppRefSyst){
135  for(Int_t i=0; i<gppRefSyst->GetN(); i++){
136  Double_t x,y;
137  gppRefSyst->GetPoint(i,x,y);
138  Int_t hBin=TMath::BinarySearch(nPtBins,binlim,x);
139  if(x>ptForExtrap){
140  relsysterrLowDatapp[hBin]=gppRefSyst->GetErrorYlow(i)/y;
141  relsysterrHiDatapp[hBin]=gppRefSyst->GetErrorYhigh(i)/y;
142  }
143  }
144  }
145  Float_t relstaterrpp[nPtBins];
146  for(Int_t ib=0; ib<nPtBins; ib++){
147  Int_t hBin=hppRef->FindBin(0.5*(binlim[ib]+binlim[ib+1]));
148  relstaterrpp[ib]=hppRef->GetBinError(hBin)/hppRef->GetBinContent(hBin);
149  printf("-- Relative stat errors PP: Bin %d(%f)--- histobin=%d Err %f\n",ib,0.5*(binlim[ib]+binlim[ib+1]),hBin,relstaterrpp[ib]);
150  Int_t hBin2=hppRefSystData->FindBin(0.5*(binlim[ib]+binlim[ib+1]));
151  if(hppRefSystData->GetBinCenter(hBin2)<ptForExtrap){
152  Float_t relsysterrDatapp=hppRefSystData->GetBinError(hBin2)/hppRefSystData->GetBinContent(hBin2);
153  relsysterrLowDatapp[ib]=relsysterrDatapp;
154  relsysterrHiDatapp[ib]=relsysterrDatapp;
155  }
156  printf(" ---- Check SYST err DATA PP Bin %d CS=%f Err+%f -%f\n",ib, hppRef->GetBinContent(hBin),relsysterrHiDatapp[ib],relsysterrLowDatapp[ib]);
157  }
158 
159  Float_t relsysterrLowEnScalpp[nPtBins];
160  Float_t relsysterrHiEnScalpp[nPtBins];
161  for(Int_t ib=0; ib<nPtBins; ib++){
162  relsysterrLowEnScalpp[ib]=0.;
163  relsysterrHiEnScalpp[ib]=0.;
164  }
165  TGraphAsymmErrors* gsystppEnSc=(TGraphAsymmErrors*)filPP->Get("gScaledDataSystExtrap");
166  for(Int_t i=0; i<gsystppEnSc->GetN(); i++){
167  Double_t x,y;
168  gsystppEnSc->GetPoint(i,x,y);
169  Int_t hBin=TMath::BinarySearch(nPtBins,binlim,x);
170  if(hBin>=0 && hBin<nPtBins && x<binlim[nPtBins]){
171  if(x<ptForExtrap){
172  relsysterrLowEnScalpp[hBin]=gsystppEnSc->GetErrorYlow(i)/y;
173  relsysterrHiEnScalpp[hBin]=gsystppEnSc->GetErrorYhigh(i)/y;
174  }
175  }
176  }
177 
178  TGraphAsymmErrors* gsystppFD=(TGraphAsymmErrors*)filPP->Get("gScaledDataSystFeedDown");
179  Float_t relsysterrLowFDpp[nPtBins];
180  Float_t relsysterrHiFDpp[nPtBins];
181  for(Int_t i=0; i<nPtBins; i++){
182  relsysterrHiFDpp[i]=0.;
183  relsysterrLowFDpp[i]=0.;
184  }
185  for(Int_t i=0; i<gsystppFD->GetN(); i++){
186  Double_t x,y;
187  gsystppFD->GetPoint(i,x,y);
188  Int_t hBin=TMath::BinarySearch(nPtBins,binlim,x);
189  if(hBin>=0 && hBin<nPtBins && x<binlim[nPtBins]){
190  relsysterrLowFDpp[hBin]=gsystppFD->GetErrorYlow(i)/y;
191  relsysterrHiFDpp[hBin]=gsystppFD->GetErrorYhigh(i)/y;
192  printf(" ---- Check syst err FD pp Bin %d Err+%f -%f\n",hBin,relsysterrHiFDpp[hBin],relsysterrLowFDpp[hBin]);
193  // printf("%d %f %f\n",hBin,relsysterrLowFDpp[hBin],relsysterrHiFDpp[hBin]);
194  }
195  }
196 
197  // A-A events
198 
199  printf("--- %s events ---\n",collSyst.Data());
200 
201  TFile *filChi= new TFile(filnamChi.Data());
202  TH1D *hSpC = (TH1D*)filChi->Get("hRECpt");
203  Float_t relstaterrPbPb[nPtBins];
204  for(Int_t ib=0; ib<nPtBins; ib++){
205  Int_t hBin=hSpC->FindBin(0.5*(binlim[ib]+binlim[ib+1]));
206  relstaterrPbPb[ib]=hSpC->GetBinError(hBin)/hSpC->GetBinContent(hBin);
207  printf("-- Relative stat errors AA from yield: Bin %d(%d) pt=%f Err %f\n",ib,hBin,hSpC->GetBinCenter(hBin),relstaterrPbPb[ib]);
208  }
209 
210  AliHFSystErr *systematicsABcent =(AliHFSystErr*)filChi->Get("AliHFSystErr");
211  if(!systematicsABcent){
212  printf("AliHFSysstErr generated on the fly\n");
213  systematicsABcent=new AliHFSystErr("AliHFSystErr","on the fly");
214  if(collSyst=="p-Pb"){
215  systematicsABcent->SetCollisionType(2);
216  systematicsABcent->SetRunNumber(16);
217  }else{
218  systematicsABcent->SetCollisionType(1);
219  }
220  systematicsABcent->SetCentrality(centralityNoDash.Data());
221  systematicsABcent->Init(mesCode);
222  }else{
223  printf("AliHFSystErr read from HFPtSpectrum output\n");
224  }
225  TH1D* heffC=(TH1D*)filChi->Get("hDirectEffpt");
226 
227 
228  TGraphAsymmErrors* gsystaaFDc=(TGraphAsymmErrors*)filChi->Get("gSigmaCorrConservative"); // FD due to scales
229  Float_t relsysterrLowFDaa[nPtBins];
230  Float_t relsysterrHiFDaa[nPtBins];
231  for(Int_t i=0; i<nPtBins; i++){
232  relsysterrLowFDaa[i]=0.;
233  relsysterrHiFDaa[i]=0.;
234  }
235  for(Int_t i=0; i<gsystaaFDc->GetN(); i++){
236  Double_t x,y;
237  gsystaaFDc->GetPoint(i,x,y);
238  Int_t hBin=TMath::BinarySearch(nPtBins,binlim,x);
239  if(hBin>=0 && hBin<nPtBins && x<binlim[nPtBins]){
240  relsysterrLowFDaa[hBin]=gsystaaFDc->GetErrorYlow(i)/y;
241  relsysterrHiFDaa[hBin]=gsystaaFDc->GetErrorYhigh(i)/y;
242  printf(" ---- Check syst err FD AA Bin %d Err+%f -%f\n",hBin,relsysterrHiFDaa[hBin],relsysterrLowFDaa[hBin]);
243  }
244  }
245 
246  TFile *filCnt=new TFile(filnamCnt.Data());
247  Float_t relstaterrPbPb2[nPtBins];
248  TH1D* hraaCcheck2=(TH1D*)filCnt->Get("hRABvsPt");
249  for(Int_t ib=0; ib<nPtBins; ib++){
250  Int_t hBin=hraaCcheck2->FindBin(0.5*(binlim[ib]+binlim[ib+1]));
251  Double_t aux=hraaCcheck2->GetBinError(hBin)/hraaCcheck2->GetBinContent(hBin);
252  relstaterrPbPb2[ib]=TMath::Sqrt(aux*aux-relstaterrpp[ib]*relstaterrpp[ib]);
253 
254  printf("-- Relative stat errors AA from RAA-PP: Bin %d(%f)---%d Err %f\n",ib,0.5*(binlim[ib]+binlim[ib+1]),hBin,relstaterrPbPb2[ib]);
255  }
256  AliHFSystErr *systematicsPP =(AliHFSystErr*)filCnt->Get("AliHFSystErrPP");
257 
258  TNtuple* ntC=(TNtuple*)filCnt->Get("ntupleRAB");
259  Float_t pt,TAB,sigmaPP,invyieldAB,RABCharm,RABBeauty;
260  Float_t invyieldABFDHigh,invyieldABFDLow,fprompt;
261  ntC->SetBranchAddress("pt",&pt);
262  ntC->SetBranchAddress("TAB",&TAB);
263  ntC->SetBranchAddress("sigmaPP",&sigmaPP);
264  ntC->SetBranchAddress("invyieldAB",&invyieldAB);
265  ntC->SetBranchAddress("RABCharm",&RABCharm);
266  ntC->SetBranchAddress("RABBeauty",&RABBeauty);
267  ntC->SetBranchAddress("invyieldABFDHigh",&invyieldABFDHigh);
268  ntC->SetBranchAddress("invyieldABFDLow",&invyieldABFDLow);
269  ntC->SetBranchAddress("fc",&fprompt);
270 
271  Float_t raac[nPtBins],invypp[nPtBins],invyPbPb[nPtBins],invyPbPbLo[nPtBins],invyPbPbHi[nPtBins],minval[nPtBins];
272  Float_t raacLowHyp[nPtBins],invyPbPbLowHyp[nPtBins],minvalLowHyp[nPtBins];
273  Float_t raacHigHyp[nPtBins],invyPbPbHigHyp[nPtBins],minvalHigHyp[nPtBins];
274  Float_t fPromptCent[nPtBins],fPromptHigHyp[nPtBins],fPromptLowHyp[nPtBins];
275  Float_t invyPbPbLoSingleSyst[nPtBins],invyPbPbHiSingleSyst[nPtBins];
276 
277  for(Int_t ib=0; ib<nPtBins; ib++){
278  minval[ib]=9999.;
279  invypp[ib]=-999.;
280  invyPbPb[ib]=-999.;
281  invyPbPbLo[ib]=-999.;
282  invyPbPbHi[ib]=-999.;
283  invyPbPbLoSingleSyst[ib]=-999.;
284  invyPbPbHiSingleSyst[ib]=-999.;
285  raac[ib]=-999.;
286  raacLowHyp[ib]=-999.;
287  invyPbPbLowHyp[ib]=-999.;
288  minvalLowHyp[ib]=9999.;
289  raacHigHyp[ib]=-999.;
290  invyPbPbHigHyp[ib]=-999.;
291  minvalHigHyp[ib]=9999.;
292  fPromptCent[ib]=-9999.;
293  fPromptHigHyp[ib]=-9999.;
294  fPromptLowHyp[ib]=-9999.;
295  }
296 
297  TGraph** gcb=new TGraph*[nPtBins];
298  TGraph** gcrbc=new TGraph*[nPtBins];
299  for(Int_t ib=0; ib<nPtBins; ib++){
300  gcb[ib]=new TGraph(0);
301  gcrbc[ib]=new TGraph(0);
302  }
303 
304  for(Int_t ie=0; ie<ntC->GetEntries(); ie++){
305  ntC->GetEvent(ie);
306  if(correctForBR){
307  invyieldAB/=brat;
308  invyieldABFDLow/=brat;
309  invyieldABFDHigh/=brat;
310  sigmaPP/=brat;
311  }
312  if(pt>binlim[nPtBins]) continue;
313  Int_t theBin=TMath::BinarySearch(nPtBins,binlim,(Double_t)pt);
314  if(theBin<0 || theBin>=nPtBins) continue;
315  Float_t rFdPr=RABBeauty/RABCharm;
316  Float_t dist=TMath::Abs(rFdPr-centHypoFdOverPr);
317  if(dist<minval[theBin]){
318  raac[theBin]=RABCharm;
319  minval[theBin]=dist;
320  invyPbPb[theBin]=invyieldAB*normToCsec;
321  invyPbPbLo[theBin]=invyieldABFDLow*normToCsec;
322  invyPbPbHi[theBin]=invyieldABFDHigh*normToCsec;
323  invypp[theBin]=TAB*sigmaPP;
324  if(collSyst=="p-Pb" && TMath::Abs(normToCsec-1.)>0.001) invypp[theBin]=208.*sigmaPP*1e6;
325  fPromptCent[theBin]=fprompt;
326  }
327  Float_t distLowHyp=TMath::Abs(rFdPr-lowHypoFdOverPr);
328  if(distLowHyp<minvalLowHyp[theBin]){
329  // LowHyp -> lower Raa(feeddown) -> less Fd to be subtracted -> higher fprompt -> higher prompt yield -> higher Raa(prompt)
330  raacLowHyp[theBin]=RABCharm;
331  invyPbPbLowHyp[theBin]=invyieldAB*normToCsec;
332  invyPbPbHiSingleSyst[theBin]=invyieldABFDHigh*normToCsec;
333  minvalLowHyp[theBin]=distLowHyp;
334  fPromptLowHyp[theBin]=fprompt;
335  }
336  Float_t distHigHyp=TMath::Abs(rFdPr-highHypoFdOverPr);
337  if(distHigHyp<minvalHigHyp[theBin]){
338  // HigHyp -> higher Raa(feeddown) -> more Fd to be subtracted -> lower fprompt -> lower prompt yield -> lower Raa(prompt)
339  raacHigHyp[theBin]=RABCharm;
340  invyPbPbHigHyp[theBin]=invyieldAB*normToCsec;
341  invyPbPbLoSingleSyst[theBin]=invyieldABFDLow*normToCsec;
342  minvalHigHyp[theBin]=distHigHyp;
343  fPromptHigHyp[theBin]=fprompt;
344  }
345  if(theBin<nPtBins && theBin>=0 && pt<binlim[nPtBins]){
346  gcb[theBin]->SetPoint(gcb[theBin]->GetN(),RABBeauty,RABCharm);
347  gcrbc[theBin]->SetPoint(gcrbc[theBin]->GetN(),rFdPr,RABCharm);
348  }
349  }
350 
351  TH1D* hcheck[nPtBins];
352  for(Int_t i=0; i<nPtBins; i++){
353  hcheck[i]=(TH1D*)filCnt->Get(Form("hRCharmVsElossHypo_%d",i+1));
354  }
355 
356  TMarker** mC=new TMarker*[nPtBins];
357  for(Int_t ib=0; ib<nPtBins; ib++){
358  mC[ib]=new TMarker(1.,raac[ib],20);
359  mC[ib]->SetMarkerSize(1.2);
360  if(showRcbSystNorm){
361  Double_t auxx,auxy;
362  for(Int_t ip=0; ip<gcrbc[ib]->GetN(); ip++){
363  gcrbc[ib]->GetPoint(ip,auxx,auxy);
364  auxy/=raac[ib];
365  gcrbc[ib]->SetPoint(ip,auxx,(auxy-1)*100.);
366  }
367  }
368  }
369 
370  TFile *filCntS=new TFile(filnamCntS.Data());
371  TNtuple* ntCS=(TNtuple*)filCntS->Get("ntupleRAB");
372  ntCS->SetBranchAddress("pt",&pt);
373  ntCS->SetBranchAddress("TAB",&TAB);
374  ntCS->SetBranchAddress("sigmaPP",&sigmaPP);
375  ntCS->SetBranchAddress("invyieldAB",&invyieldAB);
376  ntCS->SetBranchAddress("RABCharm",&RABCharm);
377  ntCS->SetBranchAddress("RABBeauty",&RABBeauty);
378  ntCS->SetBranchAddress("invyieldABFDHigh",&invyieldABFDHigh);
379  ntCS->SetBranchAddress("invyieldABFDLow",&invyieldABFDLow);
380  ntCS->SetBranchAddress("fc",&fprompt);
381 
382  Float_t invyPbPbS[nPtBins],invyPbPbLoS[nPtBins],invyPbPbHiS[nPtBins],minvalS[nPtBins];
383  Float_t invyPbPbLoSingleSystS[nPtBins],invyPbPbHiSingleSystS[nPtBins],minvalLowHypS[nPtBins],minvalHigHypS[nPtBins];
384  Float_t fPromptCentS[nPtBins],fPromptLowHypS[nPtBins],fPromptHigHypS[nPtBins];
385  for(Int_t ib=0; ib<nPtBins; ib++){
386  minvalS[ib]=9999.;
387  minvalLowHypS[ib]=9999.;
388  minvalHigHypS[ib]=9999.;
389  invyPbPbS[ib]=-999.;
390  invyPbPbLoS[ib]=-999.;
391  invyPbPbHiS[ib]=-999.;
392  invyPbPbLoSingleSystS[ib]=-999.;
393  invyPbPbHiSingleSystS[ib]=-999.;
394  fPromptLowHypS[ib]=-999.;
395  fPromptHigHypS[ib]=-999.;
396  }
397  for(Int_t ie=0; ie<ntCS->GetEntries(); ie++){
398  ntCS->GetEvent(ie);
399  if(correctForBR){
400  invyieldAB/=brat;
401  invyieldABFDLow/=brat;
402  invyieldABFDHigh/=brat;
403  sigmaPP/=brat;
404  }
405  if(pt>binlim[nPtBins]) continue;
406  Int_t theBin=TMath::BinarySearch(nPtBins,binlim,(Double_t)pt);
407  if(theBin<0 || theBin>=nPtBins) continue;
408  Float_t rFdPr=RABBeauty/RABCharm;
409  Float_t dist=TMath::Abs(rFdPr-centHypoFdOverPr);
410  if(dist<minvalS[theBin]){
411  minvalS[theBin]=dist;
412  invyPbPbS[theBin]=invyieldAB*normToCsec;
413  invyPbPbLoS[theBin]=invyieldABFDLow*normToCsec;
414  invyPbPbHiS[theBin]=invyieldABFDHigh*normToCsec;
415  fPromptCentS[theBin]=fprompt;
416  }
417  Float_t distLowHyp=TMath::Abs(rFdPr-lowHypoFdOverPr);
418  if(distLowHyp<minvalLowHypS[theBin]){
419  // LowHyp -> lower Raa(feeddown) -> less Fd to be subtracted -> higher fprompt -> higher prompt yield -> higher Raa(prompt)
420  invyPbPbHiSingleSystS[theBin]=invyieldABFDHigh*normToCsec;
421  minvalLowHypS[theBin]=distLowHyp;
422  fPromptLowHypS[theBin]=fprompt;
423  }
424  Float_t distHigHyp=TMath::Abs(rFdPr-highHypoFdOverPr);
425  if(distHigHyp<minvalHigHypS[theBin]){
426  // HigHyp -> higher Raa(feeddown) -> more Fd to be subtracted -> lower fprompt -> lower prompt yield -> lower Raa(prompt)
427  invyPbPbLoSingleSystS[theBin]=invyieldABFDLow*normToCsec;
428  minvalHigHypS[theBin]=distHigHyp;
429  fPromptHigHypS[theBin]=fprompt;
430  }
431  }
432 
433  TH1F* hfPromptCent=new TH1F("hfPromptCent"," ; p_{T} (Gev/c) ; f_{prompt}",nPtBins,binlim);
434  TH1F* hfPromptMinNb=new TH1F("hfPromptMinNb"," ; p_{T} (Gev/c) ; f_{prompt}",nPtBins,binlim);
435  TH1F* hfPromptMinfc=new TH1F("hfPromptMinfc"," ; p_{T} (Gev/c) ; f_{prompt}",nPtBins,binlim);
436  TH1F* hfPromptMaxNb=new TH1F("hfPromptMaxNb"," ; p_{T} (Gev/c) ; f_{prompt}",nPtBins,binlim);
437  TH1F* hfPromptMaxfc=new TH1F("hfPromptMaxfc"," ; p_{T} (Gev/c) ; f_{prompt}",nPtBins,binlim);
438  hfPromptCent->SetStats(0);
439  hfPromptMinNb->SetStats(0);
440  hfPromptMaxNb->SetStats(0);
441  hfPromptMinfc->SetStats(0);
442  hfPromptMaxfc->SetStats(0);
443 
444  for(Int_t ib=0; ib<nPtBins; ib++){
445  printf("Bin %d\n",ib);
446  printf(" fprompt central=%f --- Nb fpromptmin=%f fpromptmax=%f --- fc fpromptmin=%f fpromptmax=%f\n",fPromptCent[ib],fPromptHigHyp[ib],fPromptLowHyp[ib],fPromptHigHypS[ib],fPromptLowHypS[ib]);
447 
448  //add error from FONLL scale in quadrature
449  Double_t relerrFDscaleHigh = (invyPbPbHi[ib]-invyPbPb[ib])/invyPbPb[ib];
450  Double_t relerrFDscaleLow = (invyPbPb[ib]-invyPbPbLo[ib])/invyPbPb[ib];
451 
452  Double_t relerrElossHigh = (fPromptLowHyp[ib]-fPromptCent[ib])/fPromptCent[ib];
453  Double_t relerrElossLow = (fPromptCent[ib]-fPromptHigHyp[ib])/fPromptCent[ib];
454 
455  Double_t toterrFDhigh = TMath::Sqrt(relerrFDscaleHigh*relerrFDscaleHigh+relerrElossHigh*relerrElossHigh)*fPromptCent[ib];
456  Double_t toterrFDlow = TMath::Sqrt(relerrFDscaleLow*relerrFDscaleLow+relerrElossLow*relerrElossLow)*fPromptCent[ib];
457 
458  Double_t relerrFDscaleHighS = (invyPbPbHiS[ib]-invyPbPbS[ib])/invyPbPbS[ib];
459  Double_t relerrFDscaleLowS = (invyPbPbS[ib]-invyPbPbLoS[ib])/invyPbPbS[ib];
460 
461  Double_t relerrElossHighS = (fPromptLowHypS[ib]-fPromptCentS[ib])/fPromptCentS[ib];
462  Double_t relerrElossLowS = (fPromptCentS[ib]-fPromptHigHypS[ib])/fPromptCentS[ib];
463 
464  Double_t toterrFDhighS = TMath::Sqrt(relerrFDscaleHighS*relerrFDscaleHighS+relerrElossHighS*relerrElossHighS)*fPromptCentS[ib];
465  Double_t toterrFDlowS = TMath::Sqrt(relerrFDscaleLowS*relerrFDscaleLowS+relerrElossLowS*relerrElossLowS)*fPromptCentS[ib];
466 
467  hfPromptCent->SetBinContent(ib+1,fPromptCent[ib]);
468  hfPromptMinNb->SetBinContent(ib+1,fPromptCent[ib]-toterrFDlow);
469  hfPromptMaxNb->SetBinContent(ib+1,fPromptCent[ib]+toterrFDhigh);
470  hfPromptMinfc->SetBinContent(ib+1,fPromptCentS[ib]-toterrFDlowS);
471  hfPromptMaxfc->SetBinContent(ib+1,fPromptCentS[ib]+toterrFDhighS);
472 
473  printf("Bin %d\n",ib);
474  printf(" fprompt central=%f --- Nb fpromptmin=%f fpromptmax=%f --- fc fpromptmin=%f fpromptmax=%f\n",fPromptCent[ib],hfPromptMinNb->GetBinContent(ib+1),hfPromptMaxNb->GetBinContent(ib+1),hfPromptMinfc->GetBinContent(ib+1),hfPromptMaxfc->GetBinContent(ib+1));
475  }
476 
477  TH1F* hppC=new TH1F("hppC",Form("pp reference for %s%% CC",centrality.Data()),nPtBins,binlim);
478  TGraphAsymmErrors *gppCsystdata=new TGraphAsymmErrors(0);
479  gppCsystdata->SetName("gppCsystdata");
480  gppCsystdata->SetTitle(Form("Data Syst. Err. pp, scaled to %s%% CC",centrality.Data()));
481  TGraphAsymmErrors *gppCsystFD=new TGraphAsymmErrors(0);
482  gppCsystFD->SetName("gppCsystFD");
483  gppCsystFD->SetTitle(Form("B feed-down Syst. Err. pp, scaled to %s%% CC",centrality.Data()));
484  TH1F* hAAC=new TH1F("hAAC",Form("PbPb %s%% CC",centrality.Data()),nPtBins,binlim);
485  TGraphAsymmErrors *gaaCsystdata=new TGraphAsymmErrors(0);
486  gaaCsystdata->SetName("gaaCsystdata");
487  gaaCsystdata->SetTitle(Form("Data Syst. Err. PbPb %s%% CC",centrality.Data()));
488  TGraphAsymmErrors *gaaCsystFD=new TGraphAsymmErrors(0);
489  gaaCsystFD->SetName("gaaCsystFD");
490  gaaCsystFD->SetTitle(Form("B feed-down Syst. Err. PbPb %s%% CC",centrality.Data()));
491  TGraphAsymmErrors *gaaCsystRb=new TGraphAsymmErrors(0);
492  gaaCsystRb->SetName("gaaCsystRb");
493  gaaCsystRb->SetTitle(Form("Raa(B) Syst. Err. PbPb %s%% CC",centrality.Data()));
494  TGraphAsymmErrors *gaaCsystB=new TGraphAsymmErrors(0);
495  gaaCsystB->SetName("gaaCsystB");
496  gaaCsystB->SetTitle(Form("B Syst. Err. PbPb %s%% CC",centrality.Data()));
497  TGraphAsymmErrors *gaaCsystTot=new TGraphAsymmErrors(0);
498  gaaCsystTot->SetName("gaaCsystTot");
499  gaaCsystTot->SetTitle(Form("Tot Syst. Err. PbPb %s%% CC",centrality.Data()));
500  TH1F* hRAAC=new TH1F("hRAAC","",nPtBins,binlim);
501  TGraphAsymmErrors *graaC=new TGraphAsymmErrors(0);
502 
503  for(Int_t ib=0; ib<nPtBins; ib++){
504  printf("Bin %d\n",ib);
505  printf("raa(centHyp)=%f raa(lowHyp)=%f raa(highHyp)=%f RelDiff=%f\n",
506  raac[ib],raacLowHyp[ib],raacHigHyp[ib],(raacHigHyp[ib]-raacLowHyp[ib])/raac[ib]);
507  if(raac[ib]>0.){
508  Float_t relstaterrRaa=TMath::Sqrt(relstaterrpp[ib]*relstaterrpp[ib]+relstaterrPbPb2[ib]*relstaterrPbPb2[ib]);
509  hRAAC->SetBinContent(ib+1,raac[ib]);
510  hRAAC->SetBinError(ib+1,relstaterrRaa*raac[ib]);
511  Int_t nP=graaC->GetN();
512  graaC->SetPoint(nP,hRAAC->GetBinCenter(ib+1),raac[ib]);
513  graaC->SetPointEXlow(nP,hRAAC->GetBinCenter(ib+1)-hRAAC->GetBinLowEdge(ib+1));
514  graaC->SetPointEXhigh(nP,hRAAC->GetBinLowEdge(ib+2)-hRAAC->GetBinCenter(ib+1));
515  graaC->SetPointEYlow(nP,raac[ib]-raacHigHyp[ib]);
516  graaC->SetPointEYhigh(nP,raacLowHyp[ib]-raac[ib]);
517  }
518  if(invypp[ib]>0.){
519  hppC->SetBinContent(ib+1,invypp[ib]);
520  hppC->SetBinError(ib+1,relstaterrpp[ib]*invypp[ib]);
521  Int_t nP=gppCsystdata->GetN();
522  gppCsystdata->SetPoint(nP,hppC->GetBinCenter(ib+1),invypp[ib]);
523  gppCsystdata->SetPointEXlow(nP,sizesystdata);
524  gppCsystdata->SetPointEXhigh(nP,sizesystdata);
525  Double_t edatl=relsysterrLowDatapp[ib]*invypp[ib];
526  Double_t edath=relsysterrHiDatapp[ib]*invypp[ib];
527  Double_t eenscl=relsysterrLowEnScalpp[ib]*invypp[ib];
528  Double_t eensch=relsysterrHiEnScalpp[ib]*invypp[ib];
529  Double_t elow=TMath::Sqrt(edatl*edatl+eenscl*eenscl);
530  Double_t ehig=TMath::Sqrt(edath*edath+eensch*eensch);
531  gppCsystdata->SetPointEYlow(nP,elow);
532  gppCsystdata->SetPointEYhigh(nP,ehig);
533  nP=gppCsystFD->GetN();
534  gppCsystFD->SetPoint(nP,hppC->GetBinCenter(ib+1),invypp[ib]);
535  gppCsystFD->SetPointEXlow(nP,sizesystfd);
536  gppCsystFD->SetPointEXhigh(nP,sizesystfd);
537  gppCsystFD->SetPointEYlow(nP,relsysterrLowFDpp[ib]*invypp[ib]);
538  gppCsystFD->SetPointEYhigh(nP,relsysterrHiFDpp[ib]*invypp[ib]);
539  }
540  if(invyPbPb[ib]>0.){
541  printf("yielddAA(centHyp)=%f yielddAA(lowHyp)=%f yielddAA(highHyp)=%f RelDiff=%f\n",
542  invyPbPb[ib],invyPbPbLowHyp[ib],invyPbPbHigHyp[ib],(invyPbPbHigHyp[ib]-invyPbPbLowHyp[ib])/invyPbPb[ib]);
543  hAAC->SetBinContent(ib+1,invyPbPb[ib]);
544  hAAC->SetBinError(ib+1,relstaterrPbPb2[ib]*invyPbPb[ib]);
545  Int_t nP=gaaCsystdata->GetN();
546  gaaCsystdata->SetPoint(nP,hAAC->GetBinCenter(ib+1),invyPbPb[ib]);
547  gaaCsystdata->SetPointEXlow(nP,sizesystdata);
548  gaaCsystdata->SetPointEXhigh(nP,sizesystdata);
549  Double_t systundatalow,systundatahigh;
550  Bool_t isOk=PbPbDataSyst(systematicsABcent,heffC,hAAC->GetBinCenter(ib+1),systundatahigh,systundatalow);
551  printf("Check PID syst: %f %f %f\n",systematicsABcent->GetTotalSystErr(hAAC->GetBinCenter(ib+1)),systundatahigh,systundatalow);
552  systundatalow*=invyPbPb[ib];
553  systundatahigh*=invyPbPb[ib];
554  gaaCsystdata->SetPointEYlow(nP,systundatalow);
555  gaaCsystdata->SetPointEYhigh(nP,systundatahigh);
556 
557  nP=gaaCsystRb->GetN();
558  gaaCsystRb->SetPoint(nP,hAAC->GetBinCenter(ib+1),invyPbPb[ib]);
559  gaaCsystRb->SetPointEXlow(nP,sizesystrb);
560  gaaCsystRb->SetPointEXhigh(nP,sizesystrb);
561  gaaCsystRb->SetPointEYlow(nP,invyPbPb[ib]-invyPbPbHigHyp[ib]);
562  gaaCsystRb->SetPointEYhigh(nP,invyPbPbLowHyp[ib]-invyPbPb[ib]);
563  nP=gaaCsystFD->GetN();
564  gaaCsystFD->SetPoint(nP,hAAC->GetBinCenter(ib+1),invyPbPb[ib]);
565  gaaCsystFD->SetPointEXlow(nP,sizesystfd);
566  gaaCsystFD->SetPointEXhigh(nP,sizesystfd);
567  gaaCsystB->SetPoint(nP,hAAC->GetBinCenter(ib+1),invyPbPb[ib]);
568  gaaCsystB->SetPointEXlow(nP,sizesystfd);
569  gaaCsystB->SetPointEXhigh(nP,sizesystfd);
570  gaaCsystTot->SetPoint(nP,hAAC->GetBinCenter(ib+1),invyPbPb[ib]);
571  gaaCsystTot->SetPointEXlow(nP,sizesysttot);
572  gaaCsystTot->SetPointEXhigh(nP,sizesysttot);
573  if(optErrFD==0){
574  gaaCsystFD->SetPointEYlow(nP,relsysterrLowFDaa[ib]*invyPbPb[ib]);
575  gaaCsystFD->SetPointEYhigh(nP,relsysterrHiFDaa[ib]*invyPbPb[ib]);
576  }else{
577  Double_t minyield,maxyield;
578  Double_t minyieldsing,maxyieldsing;
579  if(optErrFD==1){
580  minyield=invyPbPbLo[ib];
581  maxyield=invyPbPbHi[ib];
582  minyieldsing=invyPbPbLoSingleSyst[ib];
583  maxyieldsing=invyPbPbHiSingleSyst[ib];
584  }else{
585  minyield=TMath::Min(invyPbPbLo[ib],invyPbPbLoS[ib]);
586  maxyield=TMath::Max(invyPbPbHi[ib],invyPbPbHiS[ib]);
587  minyieldsing=TMath::Min(invyPbPbLoSingleSyst[ib],invyPbPbLoSingleSystS[ib]);
588  maxyieldsing=TMath::Max(invyPbPbHiSingleSyst[ib],invyPbPbHiSingleSystS[ib]);
589  }
590  gaaCsystFD->SetPointEYlow(nP,invyPbPb[ib]-minyield);
591  gaaCsystFD->SetPointEYhigh(nP,maxyield-invyPbPb[ib]);
592  printf("Relative syst from FD (FONLL scales+masses)=-%f +%f\n",gaaCsystFD->GetErrorYlow(nP)/invyPbPb[ib],gaaCsystFD->GetErrorYhigh(nP)/invyPbPb[ib]);
593  Double_t systunBlow=0.;
594  Double_t systunBhigh=0.;
595  if(optCombineFDEloss==0){ // envelope
596  systunBlow=invyPbPb[ib]-minyieldsing;
597  systunBhigh=maxyieldsing-invyPbPb[ib];
598  }else{ // sum in quadrature
599  Double_t e1l=gaaCsystFD->GetErrorYlow(nP);
600  Double_t e1h=gaaCsystFD->GetErrorYhigh(nP);
601  Double_t e2l=gaaCsystRb->GetErrorYlow(nP);
602  Double_t e2h=gaaCsystRb->GetErrorYhigh(nP);
603  systunBlow=TMath::Sqrt(e1l*e1l+e2l*e2l);
604  systunBhigh=TMath::Sqrt(e1h*e1h+e2h*e2h);
605  }
606  gaaCsystB->SetPointEYlow(nP,systunBlow);
607  gaaCsystB->SetPointEYhigh(nP,systunBhigh);
608 
609  Double_t totSystlow=TMath::Sqrt(systunBlow*systunBlow+systundatalow*systundatalow);
610  Double_t totSysthigh=TMath::Sqrt(systunBhigh*systunBhigh+systundatahigh*systundatahigh);
611  gaaCsystTot->SetPointEYlow(nP,totSystlow);
612  gaaCsystTot->SetPointEYhigh(nP,totSysthigh);
613  }
614  }
615  }
616 
617  TH1F* hcheckRAAC=(TH1F*)hAAC->Clone("hcheckRAAC");
618  hcheckRAAC->Divide(hppC);
619 
620 
621  // Plots
622 
623  gStyle->SetOptTitle(0);
624 
625  TCanvas* c1=new TCanvas("c1","RcVsRb");
626  TLegend* leg= new TLegend(0.19,0.18,0.55,0.41);
627  leg->SetFillColor(0);
628  leg->SetFillStyle(0);
629  leg->SetBorderSize(0);
630  leg->SetTextFont(43);
631  leg->SetTextSize(28);
632  leg->SetMargin(0.32);
633  TLegendEntry* ent;
634  Bool_t first=kTRUE;
635 
636  for(Int_t ib=0; ib<nPtBins; ib++){
637  if(draw[ib] && gcb[ib]->GetN()>0){
638  gcb[ib]->SetLineColor(colors[ib]);
639  gcb[ib]->SetLineWidth(3);
640  gcb[ib]->SetLineStyle(lstyle[ib]);
641  if(first){
642  gcb[ib]->Draw("AL");
643  gcb[ib]->GetXaxis()->SetLimits(0.,1.);
644  gcb[ib]->GetXaxis()->SetTitle("#it{R}_{AA} feed-down");
645  gcb[ib]->GetYaxis()->SetTitle("#it{R}_{AA} prompt");
646  gcb[ib]->SetMinimum(0.);
647  gcb[ib]->SetMaximum(1.);
648  first=kFALSE;
649  }else{
650  gcb[ib]->Draw("lsame");
651  }
652  ent=leg->AddEntry(gcb[ib],Form("%d<#it{p}_{T}<%d GeV/#it{c}",(Int_t)binlim[ib],(Int_t)binlim[ib+1]),"L");
653  }
654  }
655  leg->Draw();
656 
657  first=kTRUE;
658 
659  TCanvas* c2=new TCanvas("c2x","RcVsRcb",700,700);
660  c2->SetTickx();
661  c2->SetTicky();
662  // c2->SetGridx();
663  // c2->SetGridy();
664  c2->SetLeftMargin(0.14);
665  c2->SetBottomMargin(0.13);
666  c2->SetTopMargin(0.035);
667  c2->SetRightMargin(0.045);
668  for(Int_t ib=0; ib<nPtBins; ib++){
669  if(draw[ib] && gcrbc[ib]->GetN()>0){
670  gcrbc[ib]->SetLineColor(colors[ib]);
671  gcrbc[ib]->SetLineWidth(3);
672  gcrbc[ib]->SetLineStyle(lstyle[ib]);
673  if(first){
674  gcrbc[ib]->Draw("AL");
675  gcrbc[ib]->GetXaxis()->SetLimits(lowHypoFdOverPr,highHypoFdOverPr);
676  gcrbc[ib]->GetXaxis()->SetTitle("Hypothesis on (#it{R}_{AA} feed-down)/(#it{R}_{AA} prompt)");
677  gcrbc[ib]->GetYaxis()->SetTitleOffset(1.2);
678  gcrbc[ib]->GetXaxis()->SetTitleOffset(1.2);
679  gcrbc[ib]->GetYaxis()->SetTitleFont(43);
680  gcrbc[ib]->GetXaxis()->SetTitleFont(43);
681  gcrbc[ib]->GetYaxis()->SetTitleSize(30);
682  gcrbc[ib]->GetXaxis()->SetTitleSize(30);
683  gcrbc[ib]->GetYaxis()->SetLabelFont(43);
684  gcrbc[ib]->GetXaxis()->SetLabelFont(43);
685  gcrbc[ib]->GetYaxis()->SetLabelSize(28);
686  gcrbc[ib]->GetXaxis()->SetLabelSize(28);
687  if(!showRcbSystNorm){
688  gcrbc[ib]->GetYaxis()->SetTitle("#it{R}_{AA} prompt");
689  gcrbc[ib]->SetMinimum(0.);
690  gcrbc[ib]->SetMaximum(0.8);
691  }else{
692  gcrbc[ib]->GetYaxis()->SetTitle("Relative variation of #it{R}_{AA} prompt (%)");
693  gcrbc[ib]->SetMinimum(-20);
694  gcrbc[ib]->SetMaximum(20);
695  }
696  first=kFALSE;
697  }else{
698  gcrbc[ib]->Draw("lsame");
699  }
700  if(!showRcbSystNorm){
701  mC[ib]->SetMarkerColor(colors[ib]);
702  mC[ib]->Draw("same");
703  }
704  }
705  }
706 // TBox* b=new TBox(0.5,0.,2.,1.);
707 // b->SetFillColor(kYellow);
708 // b->SetFillStyle(3003);
709  // b->Draw("same");
710  leg->Draw();
711  TLatex* tali0=new TLatex(0.18,0.89,"ALICE");
712  tali0->SetNDC();
713  tali0->SetTextFont(43);
714  tali0->SetTextSize(28);
715  tali0->Draw();
716  TLatex* tpbpb=new TLatex(0.41,0.89,Form("%s, %s%%",collSyst.Data(),centrality.Data()));
717  tpbpb->SetNDC();
718  tpbpb->SetTextFont(43);
719  tpbpb->SetTextSize(28);
720  tpbpb->Draw();
721  // TLatex* tdmes=new TLatex(0.63,0.85,Form("%s meson",mesSymb.Data()));
722  TLatex* tdmes=new TLatex(0.3,0.43,Form("%s meson",mesSymb.Data()));
723  tdmes->SetNDC();
724  tdmes->SetTextFont(43);
725  tdmes->SetTextSize(28);
726  tdmes->Draw();
727 
728  TLine* lin0=new TLine(lowHypoFdOverPr,0.,highHypoFdOverPr,0.);
729  lin0->SetLineStyle(2);
730  lin0->SetLineColor(kGray+1);
731  lin0->Draw();
732  if(method==2 && optErrFD==2){
733  c2->SaveAs(Form("%s-RcVsRcb_method%d_optErrFD%d_br%d.eps",mesName.Data(),method,optErrFD,correctForBR));
734  }
735 
736 
737  TCanvas* cfp=new TCanvas("cfp","fprompt",700,700);
738  gPad->SetTickx();
739  gPad->SetTicky();
740  hfPromptCent->SetMinimum(0);
741  hfPromptCent->SetMaximum(1.05);
742  hfPromptCent->SetLineWidth(3);
743  hfPromptCent->Draw();
744  hfPromptMinNb->SetLineWidth(2);
745  hfPromptMaxNb->SetLineWidth(2);
746  hfPromptMinNb->SetLineStyle(2);
747  hfPromptMaxNb->SetLineStyle(2);
748  hfPromptMinNb->Draw("same");
749  hfPromptMaxNb->Draw("same");
750  hfPromptMinfc->SetLineWidth(2);
751  hfPromptMaxfc->SetLineWidth(2);
752  hfPromptMinfc->SetLineStyle(2);
753  hfPromptMaxfc->SetLineStyle(2);
754  hfPromptMinfc->SetLineColor(2);
755  hfPromptMaxfc->SetLineColor(2);
756  if(optErrFD==2){
757  hfPromptMinfc->Draw("same");
758  hfPromptMaxfc->Draw("same");
759  TLatex* t1=new TLatex(0.17,0.15,"fc");
760  t1->SetNDC();
761  t1->SetTextColor(2);
762  t1->Draw();
763  }
764  TLatex* t2=new TLatex(0.17,0.2,"Nb");
765  t2->SetNDC();
766  t2->Draw();
767  cfp->SaveAs(Form("fprompt-%s.eps",mesName.Data()));
768 
769  hppC->SetMarkerStyle(markerppC);
770  hppC->SetMarkerSize(msizppC);
771  hppC->SetMarkerColor(colorppC);
772  hppC->SetLineColor(colorppC);
773  hppC->SetLineWidth(linewippC);
774  hppC->SetLineStyle(linestppC);
775  gppCsystdata->SetLineColor(colorSystDatappC);
776  gppCsystdata->SetLineWidth(2);
777  gppCsystdata->SetFillStyle(0);
778  gppCsystdata->SetFillColor(colorSystDatappC);
779  gppCsystFD->SetLineColor(colorSystFD);
780  gppCsystFD->SetFillColor(colorSystFD);
781  hAAC->SetMarkerStyle(markeraaC);
782  hAAC->SetMarkerSize(msizaaC);
783  hAAC->SetMarkerColor(coloraaC);
784  hAAC->SetLineColor(coloraaC);
785  hAAC->SetLineWidth(linewiaaC);
786  hAAC->SetLineStyle(linestaaC);
787  gaaCsystdata->SetLineColor(colorSystDataaaC);
788  gaaCsystdata->SetLineWidth(2);
789  gaaCsystdata->SetFillStyle(0);
790 
791  gaaCsystB->SetFillColor(colorSystRb);
792  gaaCsystRb->SetFillColor(colorSystRb);
793  gaaCsystFD->SetFillColor(colorSystFD);
794 
795  gaaCsystTot->SetLineColor(1);
796  gaaCsystTot->SetLineWidth(3);
797  gaaCsystTot->SetFillStyle(0);
798 
799 
800  Float_t ymax=4.7*hppC->GetMaximum();
801  Float_t ymin=999.;
802  Float_t y16C=hAAC->GetBinContent(hAAC->FindBin(14.));
803  if(y16C>0. && y16C<ymin) ymin=y16C;
804  Float_t y12C=hAAC->GetBinContent(hAAC->FindBin(10.));
805  if(y12C>0. && y12C<ymin) ymin=y12C;
806  ymin*=0.255;
807  ymax=19.6;
808  ymin=7E-7;
809  if(collSyst=="p-Pb" && TMath::Abs(normToCsec-1.)>0.001){
810  ymax=110000.;
811  ymin=1.1;
812  }
813 
814  TH2F *hempty=new TH2F("hempty","",100,0.,binlim[nPtBins]*1.02,100,ymin,ymax);
815  hempty->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
816  hempty->GetYaxis()->SetTitle("d#it{N}/d#it{p}_{T}_{ }|_{ |#it{y}|<0.5} (1/GeV/#it{c})");
817  hempty->GetYaxis()->SetTitleSize(0.05);
818  hempty->GetXaxis()->SetTitleSize(0.05);
819  hempty->GetYaxis()->SetTitleOffset(1.3);
820  hempty->GetYaxis()->SetTitleFont(42);
821  hempty->GetXaxis()->SetTitleFont(42);
822  hempty->GetYaxis()->SetLabelFont(42);
823  hempty->GetXaxis()->SetLabelFont(42);
824  hempty->SetStats(0);
825  TLatex* tdec =new TLatex(0.18,0.89,Form("%s",mesSymb.Data()));
826  tdec->SetNDC();
827  tdec->SetTextSize(0.038);
828  tdec->SetTextFont(42);
829 
830  gStyle->SetMarkerSize(1.8);
831 
832  TCanvas* c3=new TCanvas("c3","Yield1Pad",750,800);
833  c3->SetLeftMargin(0.15);
834  c3->SetRightMargin(0.07);
835  c3->SetBottomMargin(0.12);
836  c3->SetTopMargin(0.07);
837  c3->SetLogy();
838  c3->SetTicky();
839  hempty->Draw();
840  gppCsystFD->Draw("e2same");
841  gppCsystdata->DrawClone("e2same");
842  gaaCsystFD->Draw("E2same");
843  gaaCsystRb->Draw("E2same");
844  gaaCsystdata->DrawClone("e2same");
845  hppC->DrawClone("Psame");
846  hAAC->DrawClone("Psame");
847  // TLegend* legC= new TLegend(0.59,0.77,0.89,0.91);
848  TLegend* legC= new TLegend(0.53,0.55,0.89,0.69);
849  //TLegend* legC= new TLegend(0.59,0.37,0.89,0.51);
850  legC->SetHeader(Form("Centrality %s",centrality.Data()));
851  legC->SetFillColor(0);
852  legC->SetTextFont(42);
853  ent=legC->AddEntry(hppC,"p-p rescaled reference","PL");
854  ent=legC->AddEntry("","(#pm 6% norm. unc. not shown)","");
855  //ent->SetTextColor(hppC->GetLineColor());
856  ent=legC->AddEntry(hAAC,collSyst.Data(),"PL");
857  //ent->SetTextColor(hAAC->GetLineColor());
858  legC->Draw();
859  //TLegend* legSy= new TLegend(0.18,0.16,0.45,0.32);
860  //TLegend* legSy= new TLegend(0.62,0.57,0.89,0.75);
861  TLegend* legSy= new TLegend(0.566,0.715,0.89,0.86);
862  legSy->SetFillColor(0);
863  legSy->SetLineColor(kGray+2);
864  legSy->SetTextFont(42);
865  ent=legSy->AddEntry(gppCsystdata,"Syst. unc. from Data","F");
866  ent=legSy->AddEntry(gppCsystFD,"Syst. unc. from FONLL feed-down corr.","F");
867  ent=legSy->AddEntry(gaaCsystRb,"Syst. unc. from #it{R}_{AA} feed-down","F");
868  legSy->Draw();
869  tdec->Draw();
870  c3->SaveAs(Form("%s-Yields_1pad_method%d_optErrFD%d_br%d.eps",mesName.Data(),method,optErrFD,correctForBR));
871 
872 
873 
874  TH2F *hemptyr=new TH2F("hemptyr","",100,0.,binlim[nPtBins]*1.02,100,0.,1.7);
875  hemptyr->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
876  hemptyr->GetYaxis()->SetTitle(Form("#it{R}_{AA} (prompt %s)",mesSymb.Data()));
877  hemptyr->GetYaxis()->SetTitleOffset(1.2);
878  hemptyr->SetStats(0);
879 
880  TCanvas* c5=new TCanvas("c5","RAAcheck");
881  c5->SetGridy();
882  hRAAC->SetLineWidth(2);
883  hRAAC->SetMarkerStyle(20);
884  hRAAC->SetMarkerColor(1);
885  hRAAC->SetLineColor(1);
886  hRAAC->SetMinimum(0.);
887  hRAAC->SetMaximum(1.2);
888  hcheckRAAC->SetMarkerStyle(24);
889  hcheckRAAC->SetMarkerSize(1.2);
890  hcheckRAAC->SetMarkerColor(2);
891  hraaCcheck2->SetMarkerStyle(25);
892  hraaCcheck2->SetMarkerColor(6);
893  hraaCcheck2->SetMarkerSize(1.5);
894 
895  graaC->SetFillColor(kGray);
896  hemptyr->Draw();
897  graaC->Draw("E2same");
898  hRAAC->Draw("PEsame");
899  hraaCcheck2->Draw("PSAME");
900  hcheckRAAC->Draw("PSAME");
901  TLegend* legr=new TLegend(0.7,0.7,0.89,0.89);
902  ent=legr->AddEntry(hRAAC,Form("%s%%",centrality.Data()),"PL");
903  legr->SetFillColor(0);
904  legr->SetTextFont(42);
905  legr->Draw();
906  c5->Update();
907  // c5->SaveAs(Form("%s-RAA_check_method%d.eps",mesName.Data(),method));
908 
909 
910  // TCanvas* c6=new TCanvas("c6","Rccheck");
911  // c6->Divide(3,2);
912  // for(Int_t ib=0; ib<nPtBins; ib++){
913  // c6->cd(ib+1);
914  // if(gcrbc[ib]->GetN()>0){
915  // hcheck[ib]->Draw("col");
916  // gcrbc[ib]->SetLineColor(colors[ib]);
917  // gcrbc[ib]->Draw("lsame");
918  // }
919  // }
920 
921 
922  TString type="Yield";
923  if(TMath::Abs(normToCsec-1)>0.001) type="CrossSec";
924  TFile* outfil=new TFile(Form("%s%s_method%d_fd%d_br%d.root",mesName.Data(),type.Data(),method,optErrFD,correctForBR),"recreate");
925  hAAC->Write();
926  gppCsystFD->Write();
927  gppCsystdata->Write();
928  gaaCsystFD->Write();
929  gaaCsystRb->Write();
930  gaaCsystB->Write();
931  gaaCsystdata->Write();
932  gaaCsystTot->Write();
933  hppC->Write();
934  hfPromptCent->Write();
935  hfPromptMinNb->Write();
936  hfPromptMaxNb->Write();
937  hfPromptMinfc->Write();
938  hfPromptMaxfc->Write();
939  systematicsABcent->SetName("AliHFSystErrAA");
940  systematicsABcent->Write();
941  if(systematicsPP) systematicsPP->Write();
942  outfil->Close();
943 }
944 
945 
946 //____________________________________________________________
947 Bool_t PbPbDataSyst(AliHFSystErr *syst, TH1D* heff, Double_t pt, Double_t &dataSystUp, Double_t &dataSystDown)
948 {
949 
950  Double_t err = syst->GetTotalSystErr(pt)*syst->GetTotalSystErr(pt);
951  Int_t theBin=heff->FindBin(pt);
952  Double_t errRel=heff->GetBinError(theBin)/heff->GetBinContent(theBin);
953 
954  err += (errRel*errRel);
955 
956  Double_t errDown = err ;
957  Double_t errUp = err ;
958 
959 
960  dataSystUp = TMath::Sqrt(errUp);
961  dataSystDown = TMath::Sqrt(errDown);
962 
963  printf("Pt %f Bin %d Eff %f RelErrEff %f TotSyst +%f -%f\n",pt,theBin,heff->GetBinContent(theBin),errRel,dataSystUp,dataSystDown);
964 
965  return kTRUE;
966 }
967 
968 
kDstar
Definition: ComputeDmesonYield.C:27
ymax
const Double_t ymax
Definition: AddTaskCFDStar.C:7
Double_t
double Double_t
Definition: External.C:58
TH2F
Definition: External.C:236
filnamSpectrumFc
TString filnamSpectrumFc
Definition: ComputeDmesonYield.C:22
mesonSpecie
Int_t mesonSpecie
Definition: ComputeDmesonYield.C:28
filnamSpectrumNb
TString filnamSpectrumNb
Definition: ComputeDmesonYield.C:21
filnamRaaFc
TString filnamRaaFc
Definition: ComputeDmesonYield.C:24
centrality
centrality
Definition: AverageDmesonRaa.C:72
AliHFSystErr::SetCentrality
void SetCentrality(TString centrality)
Definition: AliHFSystErr.h:60
kDs
Definition: ComputeDmesonYield.C:27
lowHypoFdOverPr
Float_t lowHypoFdOverPr
Definition: ComputeDmesonYield.C:38
kDzero
Definition: ComputeDmesonYield.C:27
nPtBins
const Int_t nPtBins
Definition: ComputeDmesonYield.C:30
Int_t
int Int_t
Definition: External.C:63
kDplus
Definition: ComputeDmesonYield.C:27
Float_t
float Float_t
Definition: External.C:68
method
Int_t method
Definition: ComputeDmesonYield.C:32
TH1D
Definition: External.C:212
ComputeDmesonYield
void ComputeDmesonYield()
Definition: ComputeDmesonYield.C:53
AliHFSystErr::GetTotalSystErr
Double_t GetTotalSystErr(Double_t pt, Double_t feeddownErr=0) const
Definition: AliHFSystErr.cxx:7226
colors
Int_t colors[nPtBins]
Definition: ComputeDmesonYield.C:46
binlim
Double_t binlim[nPtBins+1]
Definition: ComputeDmesonYield.C:31
AliHFSystErr::Init
void Init(Int_t decay)
Function to initialize the variables/histograms.
Definition: AliHFSystErr.cxx:97
draw
Bool_t draw[nPtBins]
Definition: ComputeDmesonYield.C:45
normToCsec
Double_t normToCsec
Definition: ComputeDmesonYield.C:40
filnamPPref
TString filnamPPref
Definition: ComputeDmesonYield.C:19
AliHFSystErr.h
ymin
const Double_t ymin
Definition: AddTaskCFDStar.C:6
AliHFSystErr::SetCollisionType
void SetCollisionType(Int_t type)
Definition: AliHFSystErr.h:51
centHypoFdOverPr
Float_t centHypoFdOverPr
Definition: ComputeDmesonYield.C:37
highHypoFdOverPr
Float_t highHypoFdOverPr
Definition: ComputeDmesonYield.C:39
Bool_t
bool Bool_t
Definition: External.C:53
optErrFD
Int_t optErrFD
Definition: ComputeDmesonYield.C:33
dspec
dspec
Definition: ComputeDmesonYield.C:27
AliHFSystErr::SetRunNumber
void SetRunNumber(Int_t number)
Definition: AliHFSystErr.h:44
collSyst
TString collSyst
Definition: ComputeDmesonYield.C:41
lstyle
Int_t lstyle[nPtBins]
Definition: ComputeDmesonYield.C:47
PbPbDataSyst
Bool_t PbPbDataSyst(AliHFSystErr *syst, TH1D *heff, Double_t pt, Double_t &dataSystUp, Double_t &dataSystDown)
Definition: ComputeDmesonYield.C:947
filnamRaaNb
TString filnamRaaNb
Definition: ComputeDmesonYield.C:23