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