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