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