AliPhysics  1168478 (1168478)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
HFPtSpectrumRaa.C
Go to the documentation of this file.
1 #if !defined(__CINT__) || defined(__MAKECINT__)
2 #include "TFile.h"
3 #include "TH1.h"
4 #include "TH1D.h"
5 #include "TH2.h"
6 #include "TH2D.h"
7 #include "TH3.h"
8 #include "TH3D.h"
9 #include "TNtuple.h"
10 #include "TGraphAsymmErrors.h"
11 #include "TMath.h"
12 #include "TCanvas.h"
13 #include "TLegend.h"
14 #include "TROOT.h"
15 #include "TStyle.h"
16 #include "TLine.h"
17 #include "TLatex.h"
18 
19 #include "AliHFSystErr.h"
20 #include <Riostream.h>
21 #endif
22 
23 /* $Id$ */
25 //
26 // Macro to compute the Raa, taking as inputs the output of the corrected yields
27 // and the pp reference
28 //
29 // R_AB = [ ( dsigma/dpt )_AB / sigma_AB ] / <TAB> * ( dsigma/dpt )_pp
30 //
31 //
32 // Parameters:
33 // 1. ppfile = reference pp in the same pt binning
34 // 2. ABfile = corrected AB yields
35 // 3. outfile = output file name
36 // 4. decay = decay as in HFSystErr class
37 // 5. sigmaABCINT1B = cross section for normalization (**USE THE SAME AS ON 2.**)
38 // 6. fdMethod = feed-down subtraction method (kNb, kfc)
39 // 7. cc = centrality class
40 // 8. Energy = colliding energy (k276,k55)
41 // 9. MinHypo = minimum energy loss hypothesis (Default 1./3.)
42 // 10. MaxHypo = maximum energy loss hypothesis (Default 3.0)
43 // 11. MaxRb = maximum Raa(b) hypothesis (Default 6.0, won't do anything)
44 // 12. RbHypo : flag to decide whether the Eloss hypothesis is Rb or Rb/Rc
45 // 13. CentralHypo = central energy loss hypothesis, DEFAULT TO 1.0
46 // 14. isRaavsEP = flag to compute the Raa IN/OUT of plane, divides the reference by 2.0
47 // 15. ScaledAndExtrapRef: flag to tag scaled+reference pp-scaled-data
48 //
49 // Complains to : Zaida Conesa del Valle
50 //
52 
55 enum energy{ k276, k5dot023, k55 };
60 
61 
62 Bool_t printout = false;
67 
68 //____________________________________________________________
69 Bool_t PbPbDataSyst(AliHFSystErr *syst, Double_t pt, Int_t cc, Double_t &dataSystUp, Double_t &dataSystDown);
70 
71 //____________________________________________________________
73  // total^2 = data^2 + fd^2
74  Double_t data2 = total*total - fd*fd ;
75  return TMath::Sqrt( data2 );
76 }
77 
78 //____________________________________________________________
80 {
81  Int_t istart =0;
82  Int_t npoints = gr->GetN();
83  for(Int_t i=0; i<=npoints; i++){
84  Double_t x=0.,y=0.;
85  gr->GetPoint(i,x,y);
86  if ( TMath::Abs ( x - pt ) < 0.4 ) {
87  istart = i;
88  break;
89  }
90  }
91  return istart;
92 }
93 
94 
95 //
96 //
97 // R_AB = [ ( dsigma/dpt )_AB / sigma_AB ] / <TAB> * ( dsigma/dpt )_pp
98 //
99 //
100 //____________________________________________________________
101 void HFPtSpectrumRaa(const char *ppfile="HFPtSpectrum_D0Kpi_method2_rebinnedth_230311_newsigma.root",
102  const char *ABfile="HFPtSpectrum_D0Kpi_PbPbcuts_method2_rebinnedth_230311_newsigma.root",
103  const char *outfile="HFPtSpectrumRaa.root",
104  Int_t decay=1,
105  Double_t sigmaABCINT1B=54.e9,
106  Int_t fdMethod = kNb, Int_t cc=kpp, Int_t Energy=k276,
107  Double_t MinHypo=1./3., Double_t MaxHypo=3.0, Double_t MaxRb=6.0,
108  Bool_t isRbHypo=false, Double_t CentralHypo = 1.0,
109  Int_t ccestimator = kV0M,
110  Bool_t isUseTaaForRaa=true, const char *shadRbcFile="", Int_t nSigmaShad=3.0,
111  Int_t isRaavsEP=kPhiIntegrated, Bool_t isScaledAndExtrapRef=kFALSE,
112  Int_t rapiditySlice=kdefault, Int_t analysisSpeciality=kTopological)
113 {
114 
115  gROOT->Macro("$ALICE_PHYSICS/PWGHF/vertexingHF/macros/LoadLibraries.C");
116 
117  //
118  // Defining the TAB values for the given centrality class
119  //
120  Double_t Tab = 1., TabSyst = 0., A=207.2, B=207.2;
121  if ( Energy!=k276 && Energy!=k5dot023) {
122  printf("\n The Tab values for this cms energy have not yet been implemented, please do it ! \n");
123  return;
124  }
125  if ( cc == kpp ){
126  Tab = 1.;
127  }
128  // Values from Alberica's twiki:
129  // https://twiki.cern.ch/twiki/bin/viewauth/ALICE/CentStudies
130  if( (ccestimator == kV0M) && (Energy==k276) ) {
131  if ( cc == k07half ) {
132  Tab = 24.81; TabSyst = 0.8037;
133  } else if ( cc == k010 ) {
134  Tab = 23.48; TabSyst = 0.97;
135  } else if ( cc == k1020 ) {
136  Tab = 14.4318; TabSyst = 0.5733;
137  } else if ( cc == k020 ) {
138  Tab = 18.93; TabSyst = 0.74;
139  } else if ( cc == k2040 ) {
140  Tab = 6.86; TabSyst = 0.28;
141  } else if ( cc == k2030 ) {
142  Tab = 8.73769; TabSyst = 0.370219;
143  } else if ( cc == k3040 ) {
144  Tab = 5.02755; TabSyst = 0.22099;
145  } else if ( cc == k4050 ) {
146  Tab = 2.68327; TabSyst = 0.137073;
147  } else if ( cc == k3050 ) {
148  Tab = 3.87011; TabSyst = 0.183847;
149  } else if ( cc == k4060 ) {
150  Tab = 2.00; TabSyst= 0.11;
151  } else if ( cc == k4080 ) {
152  Tab = 1.20451; TabSyst = 0.071843;
153  } else if ( cc == k5060 ) {
154  Tab = 1.32884; TabSyst = 0.0929536;
155  } else if ( cc == k6080 ) {
156  Tab = 0.419; TabSyst = 0.033;
157  } else if ( cc == k5080 ) {
158  Tab = 0.719; TabSyst = 0.054;
159  } else if ( cc == k80100 ){
160  Tab = 0.0690; TabSyst = 0.0062;
161  }
162  }
163  if( (ccestimator == kV0M) && (Energy==k5dot023) ) {
164  if ( cc == k3050 ) {
165  Tab = 3.76; TabSyst = 0.13;
166  }
167  }
168 
169 
170  // pPb Glauber (A. Toia)
171  // https://twiki.cern.ch/twiki/bin/viewauth/ALICE/PACentStudies#Glauber_Calculations_with_sigma
172  if( cc == kpPb0100 ){
173  Tab = 0.098334; TabSyst = 0.0070679;
174  A=207.2; B=1.;
175  }
176  else if( ccestimator == kV0A ){
177  if ( cc == kpPb020 ) {
178  Tab = 0.183; TabSyst = 0.006245;
179  } else if ( cc == kpPb2040 ) {
180  Tab = 0.134; TabSyst = 0.004899;
181  } else if ( cc == kpPb4060 ) {
182  Tab = 0.092; TabSyst = 0.004796;
183  } else if ( cc == kpPb60100 ) {
184  Tab = 0.041; TabSyst = 0.008832;
185  }
186  }
187  else if( ccestimator == kZNA ){
188  if ( cc== kpPb010 ){
189  Tab = 0.17; TabSyst = 0.01275;
190  else if ( cc == kpPb020 ) {
191  Tab = 0.164; TabSyst = 0.010724;
192  } else if ( cc == kpPb2040 ) {
193  Tab = 0.137; TabSyst = 0.005099;
194  } else if ( cc == kpPb4060 ) {
195  Tab = 0.1011; TabSyst = 0.006;
196  } else if ( cc == kpPb60100 ) {
197  Tab = 0.0459; TabSyst = 0.003162;
198  }
199  }
200  else if( ccestimator == kCL1 ){
201  if ( cc == kpPb020 ) {
202  Tab = 0.19; TabSyst = 0.007;
203  } else if ( cc == kpPb2040 ) {
204  Tab = 0.136; TabSyst = 0.005;
205  } else if ( cc == kpPb4060 ) {
206  Tab = 0.088; TabSyst = 0.005;
207  } else if ( cc == kpPb60100 ) {
208  Tab = 0.0369; TabSyst = 0.0085;
209  }
210  }
211 
212  //
213  // Reading the pp file
214  //
215  TFile * ppf = new TFile(ppfile,"read");
216  TH1D * hSigmaPP;
217  TGraphAsymmErrors * gSigmaPPSyst;
218  TGraphAsymmErrors * gSigmaPPSystData = (TGraphAsymmErrors*)ppf->Get("gScaledDataSystData");
219  TGraphAsymmErrors * gSigmaPPSystTheory = (TGraphAsymmErrors*)ppf->Get("gScaledDataSystExtrap");
220  TGraphAsymmErrors * gSigmaPPSystFeedDown = (TGraphAsymmErrors*)ppf->Get("gScaledDataSystFeedDown");
221  TH1I * hCombinedReferenceFlag=0x0;
222  TGraphAsymmErrors * gReferenceFdSyst=0x0;
223  if(isScaledAndExtrapRef){
224  hCombinedReferenceFlag = (TH1I*)ppf->Get("hCombinedReferenceFlag");
225  hSigmaPP = (TH1D*)ppf->Get("hReference");
226  gSigmaPPSyst = (TGraphAsymmErrors*)ppf->Get("gReferenceSyst");
227  gReferenceFdSyst = (TGraphAsymmErrors*)ppf->Get("gReferenceFdSyst");
228  } else {
229  hSigmaPP = (TH1D*)ppf->Get("fhScaledData");
230  gSigmaPPSyst = (TGraphAsymmErrors*)ppf->Get("gScaledData");
231  }
232  Double_t scalePPRefToMatchRapidityBin = 1.0;
233 
234 
235  // Call the systematics uncertainty class for a given decay
236  AliHFSystErr *systematicsPP = 0x0;
237  if(ppf->Get("AliHFSystErr")){
238  systematicsPP = (AliHFSystErr*)ppf->Get("AliHFSystErr");
239  printf(" --> AliHFSystErr object for pp reference read from file. Object title = %s\n",systematicsPP->GetTitle());
240  }else{
241  systematicsPP = new AliHFSystErr();
242  if(analysisSpeciality==kLowPt){
243  systematicsPP->SetIsLowPtAnalysis(true);
244  }
245  if(analysisSpeciality==kPP7TeVPass4){
246  systematicsPP->SetIsPass4Analysis(true);
247  }
248  systematicsPP->Init(decay);
249  printf(" --> AliHFSystErr object for pp reference created based on macro arguments. Object title = %s\n",systematicsPP->GetTitle());
250  }
251  //
252  // Reading the AB collisions file
253  //
254  TFile * ABf = new TFile(ABfile,"read");
255  TH1D *hSigmaAB = (TH1D*)ABf->Get("histoSigmaCorr");
256  // TH2D *hSigmaABRcb = (TH2D*)ABf->Get("histoSigmaCorrRcb");
257  // TGraphAsymmErrors * gSigmaABSyst = (TGraphAsymmErrors*)ABf->Get("gSigmaCorr");
258  TGraphAsymmErrors * gSigmaABSystFeedDown = (TGraphAsymmErrors*)ABf->Get("gSigmaCorrConservative");
259  TNtuple * nSigmaAB = (TNtuple*)ABf->Get("fnSigma");
260  //
261  TH1D *hMassAB = (TH1D*)ABf->Get("hRECpt");
262  TH1D *hDirectEffptAB = (TH1D*)ABf->Get("hDirectEffpt");
263  TH1D *histofcAB = (TH1D*)ABf->Get("histofc");
264  //
265  TH1D* fhStatUncEffcSigmaAB = (TH1D*)ABf->Get("fhStatUncEffcSigma");
266  TH1D* fhStatUncEffbSigmaAB = (TH1D*)ABf->Get("fhStatUncEffbSigma");
267  TH1D* fhStatUncEffcFDAB = (TH1D*)ABf->Get("fhStatUncEffcFD");
268  TH1D* fhStatUncEffbFDAB = (TH1D*)ABf->Get("fhStatUncEffbFD");
269  //
270  TH1D* fhStatUncEffcSigmaAB_Raa = (TH1D*)fhStatUncEffcSigmaAB->Clone("fhStatUncEffcSigmaAB_Raa");
271  TH1D* fhStatUncEffbSigmaAB_Raa = (TH1D*)fhStatUncEffbSigmaAB->Clone("fhStatUncEffbSigmaAB_Raa");
272  TH1D* fhStatUncEffcFDAB_Raa = (TH1D*)fhStatUncEffcFDAB->Clone("fhStatUncEffcFDAB_Raa");
273  TH1D* fhStatUncEffbFDAB_Raa = (TH1D*)fhStatUncEffbFDAB->Clone("fhStatUncEffbFDAB_Raa");
274  fhStatUncEffcSigmaAB_Raa->Reset();
275  fhStatUncEffbSigmaAB_Raa->Reset();
276  fhStatUncEffcFDAB_Raa->Reset();
277  fhStatUncEffbFDAB_Raa->Reset();
278  fhStatUncEffcSigmaAB_Raa->SetName("fhStatUncEffcSigmaAB_Raa");
279  fhStatUncEffbSigmaAB_Raa->SetName("fhStatUncEffbSigmaAB_Raa");
280  fhStatUncEffcFDAB_Raa->SetName("fhStatUncEffcFDAB_Raa");
281  fhStatUncEffbFDAB_Raa->SetName("fhStatUncEffvFDAB_Raa");
282 
283 
284  //
285  // Call the systematics uncertainty class for a given decay
286  AliHFSystErr *systematicsAB = 0x0;
287  if(ABf->Get("AliHFSystErr")){
288  systematicsAB=(AliHFSystErr*)ABf->Get("AliHFSystErr");
289  printf(" --> AliHFSystErr object for A-A read from HFPtSpectrum file. Object title = %s\n",systematicsAB->GetTitle());
290  }else{
291  systematicsAB = new AliHFSystErr();
292  systematicsAB->SetCollisionType(1);
293  systematicsAB->SetRunNumber(2016);// check this
294  if(Energy==k276){
295  if ( cc == k07half ) systematicsAB->SetCentrality("07half");
296  else if ( cc == k010 ) systematicsAB->SetCentrality("010");
297  else if ( cc == k1020 ) systematicsAB->SetCentrality("1020");
298  else if ( cc == k020 ) systematicsAB->SetCentrality("020");
299  else if ( cc == k2040 || cc == k2030 || cc == k3040 ) {
300  systematicsAB->SetCentrality("2040");
301  systematicsAB->SetIsPbPb2010EnergyScan(true);
302  }
303  else if ( cc == k4060 || cc == k4050 || cc == k5060 ) systematicsAB->SetCentrality("4060");
304  else if ( cc == k6080 || cc == k5080 ) systematicsAB->SetCentrality("6080");
305  else if ( cc == k4080 ) systematicsAB->SetCentrality("4080");
306  else if ( cc == k3050 ) {
307  if (isRaavsEP == kPhiIntegrated) systematicsAB->SetCentrality("4080");
308  else if (isRaavsEP == kInPlane) systematicsAB->SetCentrality("3050InPlane");
309  else if (isRaavsEP == kOutOfPlane) systematicsAB->SetCentrality("3050OutOfPlane");
310  }
311  } else if (Energy==k5dot023){
312  if ( cc == k3050 ){
313  systematicsAB->SetRunNumber(15);
314  systematicsAB->SetCentrality("3050");
315  }
316  }
317  //
318  else if ( cc == kpPb0100 || cc == kpPb010 || cc == kpPb020 || cc == kpPb2040 || cc == kpPb4060 || cc == kpPb60100 ) {
319  systematicsAB->SetCollisionType(2);
320  // Rapidity slices
321  if(rapiditySlice!=kdefault){
322  systematicsAB->SetIspPb2011RapidityScan(true);
323  TString rapidity="";
324  switch(rapiditySlice) {
325  case k08to04: rapidity="0804"; scalePPRefToMatchRapidityBin=(0.093+0.280)/1.0; break;
326  case k07to04: rapidity="0804"; scalePPRefToMatchRapidityBin=0.280/1.0; break;
327  case k04to01: rapidity="0401"; scalePPRefToMatchRapidityBin=0.284/1.0; break;
328  case k01to01: rapidity="0101"; scalePPRefToMatchRapidityBin=0.191/1.0; break;
329  case k01to04: rapidity="0104"; scalePPRefToMatchRapidityBin=0.288/1.0; break;
330  case k04to07: rapidity="0408"; scalePPRefToMatchRapidityBin=0.288/1.0; break;
331  case k04to08: rapidity="0408"; scalePPRefToMatchRapidityBin=(0.288+0.096)/1.0; break;
332  case k01to05: rapidity="0401"; scalePPRefToMatchRapidityBin=0.4; break;
333  }
334  systematicsAB->SetRapidity(rapidity);
335  }
336  // Centrality slices
337  if(ccestimator==kV0A) {
338  if(cc == kpPb020) systematicsAB->SetCentrality("020V0A");
339  else if(cc == kpPb2040) systematicsAB->SetCentrality("2040V0A");
340  else if(cc == kpPb4060) systematicsAB->SetCentrality("4060V0A");
341  else if(cc == kpPb60100) systematicsAB->SetCentrality("60100V0A");
342  } else if (ccestimator==kZNA) {
343  if(cc == kpPb010) systematicsAB->SetCentrality("010ZNA");
344  else if(cc == kpPb020) systematicsAB->SetCentrality("020ZNA");
345  else if(cc == kpPb2040) systematicsAB->SetCentrality("2040ZNA");
346  else if(cc == kpPb4060) systematicsAB->SetCentrality("4060ZNA");
347  else if(cc == kpPb60100) systematicsAB->SetCentrality("60100ZNA");
348  } else if (ccestimator==kCL1) {
349  if(cc == kpPb020) systematicsAB->SetCentrality("020CL1");
350  else if(cc == kpPb2040) systematicsAB->SetCentrality("2040CL1");
351  else if(cc == kpPb4060) systematicsAB->SetCentrality("4060CL1");
352  else if(cc == kpPb60100) systematicsAB->SetCentrality("60100CL1");
353  }else {
354  if(!(cc == kpPb0100)) {
355  cout <<" Error on the pPb options"<<endl;
356  return;
357  }
358  }
359  }
360  }
361  else {
362  cout << " Systematics not yet implemented " << endl;
363  return;
364  }
365  if(analysisSpeciality==kLowPt){
366  systematicsAB->SetIsLowPtAnalysis(true);
367  }
368  else if(analysisSpeciality==kBDT){
369  systematicsAB->SetIsBDTAnalysis(true);
370  }
371  //
372  systematicsAB->Init(decay);
373  printf(" --> AliHFSystErr object for A-A created based on macro arguments. Object title = %s\n",systematicsAB->GetTitle());
374  }
375  //
376  Int_t entries = nSigmaAB->GetEntries();
377  Float_t pt=0., signal=0., Rb=0., Rcb=0., fcAB=0., yieldAB=0., sigmaAB=0., statUncSigmaAB=0., sigmaABMin=0.,sigmaABMax=0.;
378  nSigmaAB->SetBranchAddress("pt",&pt);
379  nSigmaAB->SetBranchAddress("Signal",&signal);
380  if (fdMethod==kNb) nSigmaAB->SetBranchAddress("Rb",&Rb);
381  else if (fdMethod==kfc) nSigmaAB->SetBranchAddress("Rcb",&Rcb);
382  nSigmaAB->SetBranchAddress("fc",&fcAB);
383  nSigmaAB->SetBranchAddress("Yield",&yieldAB);
384  nSigmaAB->SetBranchAddress("Sigma",&sigmaAB);
385  nSigmaAB->SetBranchAddress("SigmaStatUnc",&statUncSigmaAB);
386  nSigmaAB->SetBranchAddress("SigmaMax",&sigmaABMax);
387  nSigmaAB->SetBranchAddress("SigmaMin",&sigmaABMin);
388 
389 
390  // define the binning
391  Int_t nbins = hSigmaAB->GetNbinsX();
392  Double_t binwidth = hSigmaAB->GetBinWidth(1);
393  Double_t *limits = new Double_t[nbins+1];
394  Double_t *binwidths = new Double_t[nbins];
395  Double_t xlow=0.;
396  for (Int_t i=1; i<=nbins; i++) {
397  binwidth = hSigmaAB->GetBinWidth(i);
398  xlow = hSigmaAB->GetBinLowEdge(i);
399  limits[i-1] = xlow;
400  binwidths[i-1] = binwidth;
401  }
402  limits[nbins] = xlow + binwidth;
403 
404 
405  //
406  // Read the shadowing file if given as input
407  //
408  Double_t centralRbcShad[nbins+1], minRbcShad[nbins+1], maxRbcShad[nbins+1];
409  for(Int_t i=0; i<=nbins; i++) { centralRbcShad[i]=1.0; minRbcShad[i]=6.0; maxRbcShad[i]=0.0; }
410  Bool_t isShadHypothesis = false;
411  if( strcmp(shadRbcFile,"")!=0 ) {
412  isShadHypothesis = true;
413  cout<<endl<<">> Beware, using the shadowing prediction file with an "<<nSigmaShad<<"*sigma <<"<<endl<<endl;
414  TFile *fshad = new TFile(shadRbcFile,"read");
415  if(!fshad){ cout <<" >> Shadowing file not properly opened!!!"<<endl<<endl; return;}
416  // TH1D *hRbcShadCentral = (TH1D*)fshad->Get("hDfromBoverPromptD_Shadowing_central");
417  // TH1D *hRbcShadMin = (TH1D*)fshad->Get("hDfromBoverPromptD_Shadowing_upper");
418  // TH1D *hRbcShadMax = (TH1D*)fshad->Get("hDfromBoverPromptD_Shadowing_lower");
419  TH1D *hRbcShadCentral = (TH1D*)fshad->Get("hDfromBoverDfromc_L0");
420  TH1D *hRbcShadMin = (TH1D*)fshad->Get("hDfromBoverDfromc_L0");
421  TH1D *hRbcShadMax = (TH1D*)fshad->Get("hDfromBoverDfromc_L1");
422  if(!hRbcShadCentral || !hRbcShadMin || !hRbcShadMax) {
423  cout<< endl <<">> Shadowing input histograms are not ok !! "<<endl<<endl;
424  return;
425  }
426  // nSigmaShad
427  // nSigmaShad
428  for(Int_t i=1; i<=nbins; i++) {
429  Double_t xpt = hSigmaAB->GetBinCenter(i);
430  if(xpt>24) xpt = 20;
431  centralRbcShad[i] = hRbcShadCentral->GetBinContent( hRbcShadCentral->FindBin(xpt) );
432  Double_t minValue0 = hRbcShadMin->GetBinContent( hRbcShadMin->FindBin(xpt) );
433  Double_t maxValue0 = hRbcShadMax->GetBinContent( hRbcShadMax->FindBin(xpt) );
434  Double_t arrayEl[3] = {minValue0,maxValue0, centralRbcShad[i]};
435  Double_t minValue = TMath::MinElement(3,arrayEl);
436  Double_t maxValue = TMath::MaxElement(3,arrayEl);
437  cout<<">> Shadowing pt="<<xpt<<" central="<<centralRbcShad[i]<<" min="<<minValue<<" max="<<maxValue<<endl;
438  if(minValue>centralRbcShad[i]){ minValue = centralRbcShad[i]; }
439  if(maxValue<centralRbcShad[i]){ maxValue = centralRbcShad[i]; }
440  minRbcShad[i] = centralRbcShad[i] - nSigmaShad*(centralRbcShad[i] - minValue);
441  maxRbcShad[i] = centralRbcShad[i] + nSigmaShad*(maxValue - centralRbcShad[i]);
442  cout<<">> Shadowing hypothesis pt="<<xpt<<" central="<<centralRbcShad[i]<<" min="<<minRbcShad[i]<<" max="<<maxRbcShad[i]<<endl;
443  }
444  }
445 
446  //
447  // define the bins correspondence bw histos/files/graphs
448  //
449  //
450  TH2D * hRABvsRcb = new TH2D("hRABvsRcb"," R_{AB}(c) vs Rcb Eloss hypothesis; p_{T} [GeV/c] ; R_{AB}(c) ; Rcb Eloss hypothesis ",nbins,limits,800,0.,4.);
451  TH2D * hRABvsRb = new TH2D("hRABvsRb"," R_{AB}(c) vs Rb Eloss hypothesis; p_{T} [GeV/c] ; R_{AB}(c) ; Rb Eloss hypothesis ",nbins,limits,800,0.,4.);
452  // TH2D * hRABBeautyvsRCharm = new TH2D("hRABBeautyvsRCharm"," R_{AB}(c) vs Rb Eloss hypothesis; p_{T} [GeV/c] ; R_{AB}(b) ; R_{AB}(c) ",nbins,limits,800,0.,4.);
453  Int_t nbinsHypo=800;//200;
454  Double_t *limitsHypo = new Double_t[nbinsHypo+1];
455  for(Int_t i=1; i<=nbinsHypo+1; i++) limitsHypo[i-1]= i*4./800.;
456  TH3D * hRABCharmVsRBeautyVsPt = new TH3D("hRABCharmVsRBeautyVsPt"," R_{AB}(c) vs Rb vs p_{T} Eloss hypothesis; p_{T} [GeV/c] ; R_{AB}(b) ; R_{AB}(c) ",nbins,limits,nbinsHypo,limitsHypo,nbinsHypo,limitsHypo);
457  TH2D *hRCharmVsRBeauty[nbins+1];
458  for(Int_t i=0; i<=nbins; i++) hRCharmVsRBeauty[i] = new TH2D(Form("hRCharmVsRBeauty_%i",i),Form("RAB(c) vs RAB(b) for pt bin %i ; R_{AB}(b) ; R_{AB}(c)",i),nbinsHypo,limitsHypo,nbinsHypo,limitsHypo);
459  TH2D *hRCharmVsElossHypo[nbins+1];
460  for(Int_t i=0; i<=nbins; i++) hRCharmVsElossHypo[i] = new TH2D(Form("hRCharmVsElossHypo_%i",i),Form("RAB(c) vs ElossHypo for pt bin %i ; Eloss Hypothesis (c/b) ; R_{AB}(c)",i),nbinsHypo,limitsHypo,nbinsHypo,limitsHypo);
461  //
462  TH1D *hRABEloss00= new TH1D("hRABEloss00","hRABEloss00",nbins,limits);
463  TH1D *hRABEloss05= new TH1D("hRABEloss05","hRABEloss05",nbins,limits);
464  TH1D *hRABEloss10= new TH1D("hRABEloss10","hRABEloss10",nbins,limits);
465  TH1D *hRABEloss15= new TH1D("hRABEloss15","hRABEloss15",nbins,limits);
466  TH1D *hRABEloss20= new TH1D("hRABEloss20","hRABEloss20",nbins,limits);
467  //
468  TH2D * hRABvsRbFDlow = new TH2D("hRABvsRbFDlow"," R_{AB}(c) vs Rb Eloss hypothesis (FD low); p_{T} [GeV/c] ; Rb Eloss hypothesis ; R_{AB}(c) ",nbins,limits,800,0.,4.);
469  TH2D * hRABvsRbFDhigh = new TH2D("hRABvsRbFDhigh"," R_{AB}(c) vs Rb Eloss hypothesis (FD high); p_{T} [GeV/c] ; Rb Eloss hypothesis ; R_{AB}(c) ",nbins,limits,800,0.,4.);
470  //
471  TH1D * hRABvsRbFDhigh_proj = new TH1D("hRABvsRbFDhigh_proj","hRABvsRbFDhigh_proj",nbins,limits);
472  TH1D * hRABvsRbFDlow_proj = new TH1D("hRABvsRbFDlow_proj","hRABvsRbFDlow_proj",nbins,limits);
473  //
474  TNtuple *ntupleRAB=0x0 ;
475  if (fdMethod==kNb) {
476  ntupleRAB = new TNtuple("ntupleRAB","ntupleRAB (Nb)","pt:TAB:sigmaPP:sigmaAB:invyieldAB:invyieldABFDHigh:invyieldABFDLow:RABCharm:RABCharmFDHigh:RABCharmFDLow:RABBeauty:fc",100000);
477  } else if (fdMethod==kfc) {
478  ntupleRAB = new TNtuple("ntupleRAB","ntupleRAB (fc)","pt:TAB:sigmaPP:sigmaAB:invyieldAB:invyieldABFDHigh:invyieldABFDLow:Rcb:RABCharm:RABCharmFDHigh:RABCharmFDLow:RABBeauty:RABBeautyFDHigh:RABBeautyFDLow:fc",100000);
479  }
480  if(!ntupleRAB) printf("ERROR: Wrong method option");
481 
482  TH1D * hYieldABvsPt = new TH1D("hYieldABvsPt"," Yield_{AB}(c) vs p_{T} (no Eloss hypothesis); p_{T} [GeV/c] ; Yield_{charm} ",nbins,limits);
483  TH1D * hRABvsPt = new TH1D("hRABvsPt"," R_{AB}(c) vs p_{T} (no Eloss hypothesis); p_{T} [GeV/c] ; R_{charm} ",nbins,limits);
484  TH1D * hRABvsPt_DataSystematics = new TH1D("hRABvsPt_DataSystematics"," Systematics of R_{AB} (c) vs p_{T} (no Eloss hypothesis); p_{T} [GeV/c] ; R_{charm} ",nbins,limits);
485  TGraphAsymmErrors *gRAB_ElossHypothesis = new TGraphAsymmErrors(nbins+1);
486  gRAB_ElossHypothesis->SetNameTitle("gRAB_ElossHypothesis","RAB Eloss systematics");
487  TGraphAsymmErrors *gRAB_FeedDownSystematics = new TGraphAsymmErrors(nbins+1);
488  gRAB_FeedDownSystematics->SetNameTitle("gRAB_FeedDownSystematics","RAB Feed-Down systematics");
489  TGraphAsymmErrors *gRAB_fcFeedDownOnly = new TGraphAsymmErrors(nbins+1);
490  gRAB_fcFeedDownOnly->SetNameTitle("gRAB_fcFeedDownOnly","RAB fc Feed-Down Only");
491  TGraphAsymmErrors *gRAB_FeedDownSystematicsElossHypothesis = new TGraphAsymmErrors(nbins+1);
492  gRAB_FeedDownSystematicsElossHypothesis->SetNameTitle("gRAB_FeedDownSystematicsElossHypothesis","RAB Feed-Down systematics considering Eloss hypothesis");
493  TGraphAsymmErrors *gRAB_DataSystematics = new TGraphAsymmErrors(nbins+1);
494  gRAB_DataSystematics->SetNameTitle("gRAB_DataSystematics","RAB Measurement (no FD, no Eloss) systematics");
495  TGraphAsymmErrors *gRAB_DataSystematicsPP = new TGraphAsymmErrors(nbins+1);
496  gRAB_DataSystematicsPP->SetNameTitle("gRAB_DataSystematicsPP","RAB Measurement PP meas. systematics (data+scaling)");
497  TGraphAsymmErrors *gRAB_DataSystematicsAB = new TGraphAsymmErrors(nbins+1);
498  gRAB_DataSystematicsAB->SetNameTitle("gRAB_DataSystematicsAB","RAB Measurement AB (no FD, no Eloss, no PP data) systematics");
499  TGraphAsymmErrors *gRAB_GlobalSystematics = new TGraphAsymmErrors(nbins+1);
500  gRAB_GlobalSystematics->SetNameTitle("gRAB_GlobalSystematics","RAB Measurement global (data, FD, Eloss) systematics");
501  Double_t ElossMax[nbins+1], ElossMin[nbins+1];
502  for(Int_t i=0; i<=nbins; i++) { ElossMax[i]=0.; ElossMin[i]=6.; }
503  Double_t fcElossMax[nbins+1], fcElossMin[nbins+1];
504  for(Int_t i=0; i<=nbins; i++) { fcElossMax[i]=0.; fcElossMin[i]=6.; }
505  Double_t FDElossMax[nbins+1], FDElossMin[nbins+1];
506  for(Int_t i=0; i<=nbins; i++) { FDElossMax[i]=0.; FDElossMin[i]=6.; }
507 
508  TGraphAsymmErrors *gRAB_Norm = new TGraphAsymmErrors(1);
509  gRAB_Norm->SetNameTitle("gRAB_Norm","RAB Normalization systematics (pp norm + Tab)");
510  Double_t normUnc = TMath::Sqrt ( NormPPUnc*NormPPUnc + (TabSyst/Tab)*(TabSyst/Tab) );
511  if(!isUseTaaForRaa) normUnc = TMath::Sqrt ( NormPPUnc*NormPPUnc + NormABUnc*NormABUnc );
512  gRAB_Norm->SetPoint(1,0.5,1.);
513  gRAB_Norm->SetPointError(1,0.25,0.25,normUnc,normUnc);
514 
515  //
516  // R_AB = ( dN/dpt )_AB / <Ncoll_AB> * ( dN/dpt )_pp ; <Ncoll> = <Tab> * sigma_NN^inel
517  // R_AB = [ ( dsigma/dpt )_AB / sigma_AB ] / <TAB> * ( dsigma/dpt )_pp
518  //
519  Int_t istartPPfd=0, istartPPsyst=0, istartABfd=0, istartPPextr=0;
520  Double_t yPPh=0., yPPl=0., yABh=0., yABl=0.;
521  Double_t RaaCharm =0., RaaBeauty=0.;
522  Double_t RaaCharmFDhigh = 0., RaaCharmFDlow = 0.;
523  Double_t RaaBeautyFDhigh = 0., RaaBeautyFDlow = 0.;
524  Double_t systUp=0., systLow=0., systPPUp=0., systPPLow=0., systABUp=0., systABLow=0.;
525  //
526  //
527  // Search the central value of the energy loss hypothesis Rb = Rc (bin)
528  //
529  Double_t ElossCentral[nbins+1];
530  for(Int_t i=0; i<=nbins; i++) { ElossCentral[i]=0.; }
531  //
532  for(Int_t ientry=0; ientry<=entries; ientry++){
533 
534  nSigmaAB->GetEntry(ientry);
535  // cout << " pt="<< pt<<" sigma-AB="<<sigmaAB<<endl;
536  if ( !(sigmaAB>0.) ) continue;
537  //if(decay==2 && pt<2.) continue;
538 
539  // Compute RAB and the statistical uncertainty
540  Int_t hppbin = hSigmaPP->FindBin( pt );
541  Int_t hABbin = hSigmaAB->FindBin( pt );
542  Double_t sigmapp = hSigmaPP->GetBinContent( hppbin );
543  sigmapp *= scalePPRefToMatchRapidityBin; // scale to the proper rapidity bin width
544  // cout << " pt="<< pt<<", sigma-pp="<< sigmapp<<endl;
545  if (isRaavsEP>0.) sigmapp = 0.5*sigmapp;
546  if ( !(sigmapp>0.) ) continue;
547 
548  RaaCharm = ( sigmaAB / sigmaABCINT1B ) / ((Tab*1e3) * sigmapp *1e-12 ) ;
549  if(!isUseTaaForRaa) {
550  RaaCharm = ( sigmaAB ) / ( (A*B) * sigmapp ) ;
551  }
552 
553  if (fdMethod==kNb) {
554  RaaBeauty = Rb ;
555  }
556  else if (fdMethod==kfc) {
557  RaaBeauty = ( RaaCharm / Rcb ) ;
558  }
559 
560  Double_t ElossHypo = 0.;
561  if (fdMethod==kfc) { ElossHypo = 1. / Rcb; }
562  else { ElossHypo = 1. / (RaaCharm / RaaBeauty) ; }
563  if(isRbHypo) ElossHypo = RaaBeauty;
564 
565  // If using shadowing hypothesis, change the central hypothesis too
566  if(isShadHypothesis) CentralHypo = centralRbcShad[hABbin];
567 
568  // cout <<" pt "<< pt << " Raa charm " << RaaCharm << " Raa beauty " << RaaBeauty << " eloss hypo "<< ElossHypo<<endl;
569  //
570  // Find the bin for the central Eloss hypo
571  //
572  if( TMath::Abs( ElossHypo - CentralHypo ) < 0.075 ){
573  Double_t DeltaIni = TMath::Abs( ElossCentral[ hABbin ] - CentralHypo );
574  Double_t DeltaV = TMath::Abs( ElossHypo - CentralHypo );
575  // cout << " pt " << pt << " ECentral " << ElossCentral[ hABbin ] << " Ehypo "<< ElossHypo ;
576  if ( DeltaV < DeltaIni ) ElossCentral[ hABbin ] = ElossHypo;
577  // cout << " final ECentral " << ElossCentral[ hABbin ] << endl;
578  }
579  }
580  //
581  // Calculation of the Raa and its uncertainties
582  //
583  for(Int_t ientry=0; ientry<entries; ientry++){
584 
585  nSigmaAB->GetEntry(ientry);
586  if ( !(sigmaAB>0.) ) continue;
587  // if ( pt<2 || pt>16) continue;
588 
589 
590  // Compute RAB and the statistical uncertainty
591  Int_t hppbin = hSigmaPP->FindBin( pt );
592  Double_t sigmapp = hSigmaPP->GetBinContent( hppbin );
593  if (isRaavsEP>0.) sigmapp = 0.5*sigmapp;
594  sigmapp *= scalePPRefToMatchRapidityBin; // scale to the proper rapidity bin width
595  if ( !(sigmapp>0.) ) continue;
596 
597  RaaCharm = ( sigmaAB / sigmaABCINT1B ) / ((Tab*1e3) * sigmapp *1e-12 );
598  if(!isUseTaaForRaa) {
599  RaaCharm = ( sigmaAB ) / ( (A*B) * sigmapp ) ;
600  }
601 
602  // Flag to know if it is an scaled or extrapolated point of the pp reference
603  Bool_t isExtrapolatedBin = kFALSE;
604  if(isScaledAndExtrapRef && hCombinedReferenceFlag) isExtrapolatedBin = hCombinedReferenceFlag->GetBinContent( hppbin );
605  istartPPsyst = -1;
606  istartPPsyst = FindGraphBin(gSigmaPPSyst,pt);
607 
608  //
609  // FONLL Feed-Down systematics
610  //
611  istartPPfd = -1;
612  if(!isExtrapolatedBin) istartPPfd = FindGraphBin(gSigmaPPSystFeedDown,pt);
613  istartABfd = -1;
614  istartABfd = FindGraphBin(gSigmaABSystFeedDown,pt);
615 
616  // cout << " Starting bin for pp is "<< istartPPfd <<", for AA is "<<istartABfd << endl;
617  if(isExtrapolatedBin){
618  if(gReferenceFdSyst){
619  Int_t ibinfd = FindGraphBin(gReferenceFdSyst,pt);
620  yPPh = gReferenceFdSyst->GetErrorYhigh(ibinfd);
621  yPPl = gReferenceFdSyst->GetErrorYlow(ibinfd);
622  }
623  } else {
624  yPPh = gSigmaPPSystFeedDown->GetErrorYhigh(istartPPfd);
625  yPPl = gSigmaPPSystFeedDown->GetErrorYlow(istartPPfd);
626  }
627  if (isRaavsEP>0.) {
628  yPPh = yPPh*0.5;
629  yPPl = yPPl*0.5;
630  }
631  yPPh *= scalePPRefToMatchRapidityBin; // scale to the proper rapidity bin width
632  yPPl *= scalePPRefToMatchRapidityBin; // scale to the proper rapidity bin width
633 
634  yABh = gSigmaABSystFeedDown->GetErrorYhigh(istartABfd);
635  yABl = gSigmaABSystFeedDown->GetErrorYlow(istartABfd);
636 
637 
638  RaaCharmFDhigh = ( sigmaABMax / sigmaABCINT1B ) / ((Tab*1e3) * (sigmapp+yPPh) *1e-12 ) ;
639  RaaCharmFDlow = ( sigmaABMin / sigmaABCINT1B ) / ((Tab*1e3) * (sigmapp-yPPl) *1e-12 ) ;
640  if(printout && TMath::Abs(ptprintout-pt)<0.1 ) cout << endl<<" pt "<< pt << " Raa " << RaaCharm <<" high "<< RaaCharmFDhigh << " low "<< RaaCharmFDlow<<endl;
641  if(!isUseTaaForRaa) {
642  RaaCharmFDhigh = ( sigmaABMax ) / ( (A*B)* (sigmapp+yPPh) ) ;
643  RaaCharmFDlow = ( sigmaABMin ) / ( (A*B)* (sigmapp-yPPl) ) ;
644  }
645 
646 
647  if (fdMethod==kNb) {
648  RaaBeauty = Rb ;
649  RaaBeautyFDlow = Rb ;
650  RaaBeautyFDhigh = Rb ;
651  ntupleRAB->Fill( pt, Tab*1e3, sigmapp*1e-12, sigmaAB*1e-12, sigmaAB/sigmaABCINT1B,
652  sigmaABMax / sigmaABCINT1B, sigmaABMin / sigmaABCINT1B,
653  RaaCharm, RaaCharmFDhigh, RaaCharmFDlow, RaaBeauty, fcAB );
654  }
655  else if (fdMethod==kfc) {
656  RaaBeauty = ( RaaCharm / Rcb ) ;
657  RaaBeautyFDlow = ( RaaCharmFDlow / Rcb ) ;
658  RaaBeautyFDhigh = ( RaaCharmFDhigh / Rcb ) ;
659  hRABvsRcb->Fill( pt, RaaCharm, RaaBeauty );
660  ntupleRAB->Fill( pt, Tab*1e3, sigmapp*1e-12, sigmaAB*1e-12, sigmaAB/sigmaABCINT1B,
661  sigmaABMax / sigmaABCINT1B, sigmaABMin / sigmaABCINT1B,
662  Rcb, RaaCharm, RaaCharmFDhigh, RaaCharmFDlow, RaaBeauty, RaaBeautyFDhigh, RaaBeautyFDlow, fcAB );
663  }
664  hRABvsRb->Fill( pt, RaaCharm, RaaBeauty );
665  hRABvsRbFDlow->Fill( pt, RaaCharmFDlow, RaaBeautyFDlow );
666  hRABvsRbFDhigh->Fill( pt, RaaCharmFDhigh, RaaBeautyFDhigh );
667  if(printout && TMath::Abs(ptprintout-pt)<0.1) cout << " pt "<< pt << " Rb " << RaaBeauty <<" high "<< RaaBeautyFDhigh << " low "<< RaaBeautyFDlow <<endl;
668 
669  hRABCharmVsRBeautyVsPt->Fill( pt, RaaBeauty, RaaCharm );
670  Int_t ptbin = hRABvsPt->FindBin( pt );
671  hRCharmVsRBeauty[ptbin]->Fill( RaaBeauty, RaaCharm );
672  hRCharmVsRBeauty[ptbin]->Fill( RaaBeautyFDlow, RaaCharmFDlow );
673  hRCharmVsRBeauty[ptbin]->Fill( RaaBeautyFDhigh, RaaCharmFDhigh );
674 
675 
676  if (fdMethod==kfc) {
677  if( TMath::Abs(Rcb-0.015)<0.009 ) hRABEloss00->Fill(pt,RaaCharm);
678  if( TMath::Abs(Rcb-0.5)<0.009 ) hRABEloss05->Fill(pt,RaaCharm);
679  if( TMath::Abs(Rcb-1.0)<0.009 ) {
680  hRABEloss10->Fill(pt,RaaCharm);
681  hRABvsRbFDhigh_proj->Fill(pt,RaaCharmFDhigh);
682  hRABvsRbFDlow_proj->Fill(pt,RaaCharmFDlow);
683  }
684  if( TMath::Abs(Rcb-1.5)<0.009 ) hRABEloss15->Fill(pt,RaaCharm);
685  if( TMath::Abs(Rcb-2.0)<0.009 ) hRABEloss20->Fill(pt,RaaCharm);
686  }
687  else if (fdMethod==kNb) {
688  if( TMath::Abs(RaaBeauty-0.015)<0.009 ) hRABEloss00->Fill(pt,RaaCharm);
689  if( TMath::Abs(RaaBeauty-0.5)<0.009 ) hRABEloss05->Fill(pt,RaaCharm);
690  if( TMath::Abs(RaaBeauty-1.0)<0.009 ) {
691  hRABEloss10->Fill(pt,RaaCharm);
692  hRABvsRbFDhigh_proj->Fill(pt,RaaCharmFDhigh);
693  hRABvsRbFDlow_proj->Fill(pt,RaaCharmFDlow);
694  }
695  if( TMath::Abs(RaaBeauty-1.5)<0.009 ) hRABEloss15->Fill(pt,RaaCharm);
696  if( TMath::Abs(RaaBeauty-2.0)<0.009 ) hRABEloss20->Fill(pt,RaaCharm);
697  }
698 
699 
700  Int_t hABbin = hMassAB->FindBin( pt );
701  if(isShadHypothesis) CentralHypo = centralRbcShad[hABbin];
702 
703  if(printout && TMath::Abs(ptprintout-pt)<0.1)
704  if ( fdMethod==kNb && TMath::Abs(Rb -CentralHypo)< 0.05) {
705  cout << " pt "<< pt <<", at bin "<<hABbin<<endl;
706  cout<<" entries "<<entries<<", i="<<ientry<<", pt="<<pt<<", Rb="<<Rb<<", Tab="<<Tab<<", sigmaAB="<<sigmaAB<<", sigmapp="<<sigmapp<<", Raacharm="<<RaaCharm<<", RaaBeauty="<<RaaBeauty<<endl;
707  cout << " AB basis: mass "<< hMassAB->GetBinContent(hABbin)<<", eff "<< hDirectEffptAB->GetBinContent(hABbin)<<endl;
708  cout<<" FD low, err low AB "<< (sigmaAB-sigmaABMin)<<" err low PP "<< yPPl<<" Raacharm="<<RaaCharmFDlow<<", RaaBeauty="<<RaaBeautyFDlow<<endl;
709  cout<<" FD high, err high AB "<< (sigmaABMax-sigmaAB)<<" err high PP "<< yPPh<<" Raacharm="<<RaaCharmFDhigh<<", RaaBeauty="<<RaaBeautyFDhigh<<endl;
710  }
711  if(printout && TMath::Abs(ptprintout-pt)<0.1)
712  if ( fdMethod==kfc) if(TMath::Abs(Rcb -CentralHypo)< 0.05 ){
713  cout << " pt "<< pt <<", at bin "<<hABbin<<endl;
714  cout<<" entries "<<entries<<", i="<<ientry<<", pt="<<pt<<", Rcb="<<Rcb<<", Tab="<<Tab<<", sigmaAB="<<sigmaAB<<", sigmapp="<<sigmapp<<", Raacharm="<<RaaCharm<<", RaaBeauty="<<RaaBeauty<<endl;
715  cout << " AB basis: mass "<< hMassAB->GetBinContent(hABbin)<<", eff "<< hDirectEffptAB->GetBinContent(hABbin)<<", fc "<<histofcAB->GetBinContent(hABbin)<< endl;
716  cout<<" FD low, err low AB "<< (sigmaAB-sigmaABMin)<<" err low PP "<< yPPl<<" Raacharm="<<RaaCharmFDlow<<", RaaBeauty="<<RaaBeautyFDlow<<endl;
717  cout<<" FD high, err high AB "<< (sigmaABMax-sigmaAB)<<" err high PP "<< yPPh<<" Raacharm="<<RaaCharmFDhigh<<", RaaBeauty="<<RaaBeautyFDhigh<<endl;
718  }
719 
720 
721  //
722  // Fill in the global properties ?
723  //
724  Double_t ElossHypo = 0.;
725  if (fdMethod==kfc) { ElossHypo = 1./ Rcb; }
726  else { ElossHypo = 1. / (RaaCharm / RaaBeauty); }
727  if(isRbHypo) ElossHypo = RaaBeauty;
728  hRCharmVsElossHypo[ptbin]->Fill( ElossHypo, RaaCharm );
729 
730  // If using shadowing hypothesis, change the limit hypothesis too
731  if(isShadHypothesis) {
732  MinHypo = minRbcShad[ hABbin ];
733  MaxHypo = maxRbcShad[ hABbin ];
734  }
735 
736  // cout <<" pt "<< pt << " Raa charm " << RaaCharm << " Raa beauty " << RaaBeauty << " eloss hypo "<< ElossHypo
737  if(ientry==0) cout<<" pt"<< pt<< " ElossCentral "<< ElossCentral[hABbin] << " min-hypo "<<MinHypo << " max-hypo "<<MaxHypo<<endl;
738 
739  //
740  // Fill in histos charm (null Eloss hypothesis)
741  //
742  Double_t minFdSyst = 0., maxFdSyst = 0.;
743  if ( ElossHypo == ElossCentral[ hABbin ] ) {
744 
745  //
746  // Data stat uncertainty
747  //
748  Double_t sigmappStat = hSigmaPP->GetBinError( hppbin );
749  if (isRaavsEP>0.) sigmappStat = sigmappStat*0.5;
750  sigmappStat *= scalePPRefToMatchRapidityBin; // scale to the proper rapidity bin width
751  Int_t hRABbin = hRABvsPt->FindBin( pt );
752  Double_t stat = RaaCharm * TMath::Sqrt( (statUncSigmaAB/sigmaAB)*(statUncSigmaAB/sigmaAB) +
753  (sigmappStat/sigmapp)*(sigmappStat/sigmapp) ) ;
754  if ( RaaCharm==0 ) stat =0.;
755  if ( RaaCharm>0 ) {
756  hRABvsPt->SetBinContent( hRABbin, RaaCharm );
757  hRABvsPt->SetBinError( hRABbin, stat );
758  hYieldABvsPt->SetBinContent( hRABbin, sigmaAB/sigmaABCINT1B );
759  hYieldABvsPt->SetBinError( hRABbin, statUncSigmaAB/sigmaABCINT1B );
760 
761  cout << "pt="<< pt<< " Raa " << RaaCharm << " stat unc. "<< stat <<
762  " sigma-pp "<< sigmapp <<" sigma-AB "<< sigmaAB<<endl;
763  if(printout && TMath::Abs(ptprintout-pt)<0.1) {
764  cout << " Raa " << RaaCharm << " stat unc. "<< stat << " is "<< stat/RaaCharm * 100. <<
765  "%, stat-pp "<< sigmappStat/sigmapp*100. <<"% stat-AB "<< statUncSigmaAB/sigmaAB*100.<<"%"<<endl;
766  }
767 
768  Double_t errstatEff = fhStatUncEffcSigmaAB->GetBinError( hRABbin );
769  fhStatUncEffcSigmaAB_Raa->SetBinError( hRABbin, errstatEff*RaaCharm );
770  errstatEff = fhStatUncEffbSigmaAB->GetBinError( hRABbin );
771  fhStatUncEffbSigmaAB_Raa->SetBinError( hRABbin, errstatEff*RaaCharm );
772  errstatEff = fhStatUncEffcFDAB->GetBinError( hRABbin );
773  fhStatUncEffcFDAB_Raa->SetBinError( hRABbin, errstatEff*RaaCharm );
774  errstatEff = fhStatUncEffbFDAB->GetBinError( hRABbin );
775  fhStatUncEffbFDAB_Raa->SetBinError( hRABbin, errstatEff*RaaCharm );
776  }
777 
778 
779  //
780  //
781  // Data systematics (sigma syst-but FD + extrap) syst
782  //
783  //
784  // Data syst: a) Syst in p-p
785  //
786  Double_t ptwidth = hSigmaAB->GetBinWidth(hABbin) / 2. ;
787  istartPPextr = -1;
788  if(!isExtrapolatedBin) istartPPextr = FindGraphBin(gSigmaPPSystTheory,pt);
789 
790  Double_t dataPPUp=0., dataPPLow=0.;
791  if(isExtrapolatedBin) {
792  dataPPUp = gSigmaPPSyst->GetErrorYhigh(istartPPsyst);
793  dataPPLow = gSigmaPPSyst->GetErrorYlow(istartPPsyst);
794  systPPUp = dataPPUp;
795  systPPLow = dataPPLow;
796  } else {
797  dataPPUp = ExtractFDSyst( gSigmaPPSystData->GetErrorYhigh(istartPPextr), gSigmaPPSystFeedDown->GetErrorYhigh(istartPPfd) );
798  dataPPLow = ExtractFDSyst( gSigmaPPSystData->GetErrorYlow(istartPPextr), gSigmaPPSystFeedDown->GetErrorYlow(istartPPfd) );
799  systPPUp = TMath::Sqrt( dataPPUp*dataPPUp + gSigmaPPSystTheory->GetErrorYhigh(istartPPextr)*gSigmaPPSystTheory->GetErrorYhigh(istartPPextr) );
800  systPPLow = TMath::Sqrt( dataPPLow*dataPPLow + gSigmaPPSystTheory->GetErrorYlow(istartPPextr)*gSigmaPPSystTheory->GetErrorYlow(istartPPextr) );
801  }
802  if (isRaavsEP>0.) {
803  dataPPUp = dataPPUp*0.5;
804  dataPPLow = dataPPLow*0.5;
805  if(isExtrapolatedBin) {
806  systPPUp = dataPPUp;
807  systPPLow = dataPPLow;
808  } else {
809  systPPUp = TMath::Sqrt( dataPPUp*dataPPUp + 0.5*gSigmaPPSystTheory->GetErrorYhigh(istartPPextr)*0.5*gSigmaPPSystTheory->GetErrorYhigh(istartPPextr) );
810  systPPLow = TMath::Sqrt( dataPPLow*dataPPLow + 0.5*gSigmaPPSystTheory->GetErrorYlow(istartPPextr)*0.5*gSigmaPPSystTheory->GetErrorYlow(istartPPextr) );
811  }
812  }
813  systPPUp *= scalePPRefToMatchRapidityBin; // scale to the proper rapidity bin width
814  systPPLow *= scalePPRefToMatchRapidityBin; // scale to the proper rapidity bin width
815 
816 
817  if(printout && TMath::Abs(ptprintout-pt)<0.1) {
818  cout << " pt : "<< pt<<" Syst-pp-data "<< dataPPUp/sigmapp << "%, ";
819  if(!isExtrapolatedBin){
820  if (isRaavsEP>0.) cout <<" extr unc + "<< 0.5*gSigmaPPSystTheory->GetErrorYhigh(istartPPextr)/sigmapp <<" - "<< 0.5*gSigmaPPSystTheory->GetErrorYlow(istartPPextr)/sigmapp <<" %";
821  else cout <<" extr unc + "<< (gSigmaPPSystTheory->GetErrorYhigh(istartPPextr)*scalePPRefToMatchRapidityBin)/sigmapp <<" - "<< (gSigmaPPSystTheory->GetErrorYlow(istartPPextr)*scalePPRefToMatchRapidityBin)/sigmapp <<" %";
822  }
823  cout << endl;
824  }
825 
826  //
827  // Data syst: b) Syst in PbPb
828  //
829  Double_t dataSystUp=0., dataSystDown=0.;
830  Bool_t PbPbDataSystOk = PbPbDataSyst(systematicsAB,pt,cc,dataSystUp,dataSystDown);
831  if (!PbPbDataSystOk) { cout <<" There is some issue with the PbPb data systematics, please check and rerun"<<endl; return; }
832  systABUp = sigmaAB * TMath::Sqrt( dataSystUp*dataSystUp +
833  (hDirectEffptAB->GetBinError(hABbin)/hDirectEffptAB->GetBinContent(hABbin))*(hDirectEffptAB->GetBinError(hABbin)/hDirectEffptAB->GetBinContent(hABbin)) );
834 
835  systABLow = sigmaAB * TMath::Sqrt( dataSystDown*dataSystDown +
836  (hDirectEffptAB->GetBinError(hABbin)/hDirectEffptAB->GetBinContent(hABbin))*(hDirectEffptAB->GetBinError(hABbin)/hDirectEffptAB->GetBinContent(hABbin)) );
837  //
838  // Data syst : c) combine pp & PbPb
839  //
840  systLow = sigmapp>0. ?
841  RaaCharm * TMath::Sqrt( (systABLow/sigmaAB)*(systABLow/sigmaAB) + (systPPUp/sigmapp)*(systPPUp/sigmapp) )
842  : 0.;
843 
844  systUp = sigmapp>0. ?
845  RaaCharm * TMath::Sqrt( (systABUp/sigmaAB)*(systABUp/sigmaAB) + (systPPLow/sigmapp)*(systPPLow/sigmapp) )
846  : 0.;
847  if ( RaaCharm==0 ) { systPPUp =0.; systPPLow =0.; }
848 
849  // if(printout)
850  cout << " Syst-pp-up "<< systPPUp/sigmapp <<"%, syst-pp-low "<< systPPLow/sigmapp <<"%, syst-AB-up "<<systABUp/sigmaAB<<"%, syst-AB-low "<<systABLow/sigmaAB<<"%, tot-syst-up "<<systUp/RaaCharm<<"%, tot-syst-low "<<systLow/RaaCharm<<"%"<<endl;
851 
852  if ( RaaCharm>0 ) {
853  hRABvsPt_DataSystematics->SetBinContent( hRABbin, RaaCharm );
854  hRABvsPt_DataSystematics->SetBinError( hRABbin, systUp );
855  gRAB_DataSystematics->SetPoint( hABbin, pt, RaaCharm ); // i, x, y
856  gRAB_DataSystematics->SetPointError( hABbin, ptwidth, ptwidth, systLow, systUp );
857  gRAB_DataSystematics->SetPointEXlow(hABbin, 0.4); gRAB_DataSystematics->SetPointEXhigh(hABbin,0.4);
858  gRAB_DataSystematicsPP->SetPoint( hABbin, pt, RaaCharm ); // i, x, y
859  gRAB_DataSystematicsPP->SetPointError( hABbin, ptwidth, ptwidth, RaaCharm *(systPPUp/sigmapp), RaaCharm *systPPLow/sigmapp );
860  gRAB_DataSystematicsAB->SetPoint( hABbin, pt, RaaCharm ); // i, x, y
861  gRAB_DataSystematicsAB->SetPointError( hABbin, ptwidth, ptwidth, RaaCharm *systABLow/sigmaAB, RaaCharm *systABUp/sigmaAB );
862  }
863 
864  //
865  // Feed-down Systematics
866  //
867  Double_t FDL=0., FDH=0.;
868  if ( RaaCharmFDhigh > RaaCharmFDlow ){
869  FDH = RaaCharmFDhigh; FDL = RaaCharmFDlow;
870  } else {
871  FDL = RaaCharmFDhigh; FDH = RaaCharmFDlow;
872  }
873 
874  if(printout && TMath::Abs(ptprintout-pt)<0.1) cout<<" Raa "<<RaaCharm<<", Raa-fd-low "<<RaaCharmFDlow <<", Raa-fd-high "<<RaaCharmFDhigh <<endl;
875  maxFdSyst = TMath::Abs(FDH - RaaCharm);
876  minFdSyst = TMath::Abs(RaaCharm - FDL);
877  if ( RaaCharm>0 ) {
878  gRAB_FeedDownSystematics->SetPoint( hABbin, pt, RaaCharm ); // i, x, y
879  gRAB_FeedDownSystematics->SetPointError( hABbin, 0.3, 0.3, minFdSyst, maxFdSyst ); // i, x, y
880  gRAB_fcFeedDownOnly->SetPoint( hABbin, pt,fcAB );
881  gRAB_fcFeedDownOnly->SetPointError(hABbin, 0.3, 0.3, fcAB-(sigmaABMin/sigmaAB*fcAB), (sigmaABMax/sigmaAB*fcAB)-fcAB );
882  }
883 
884  // if(printout) {
885  cout<<" FD syst +"<< maxFdSyst/RaaCharm <<" - "<<minFdSyst/RaaCharm<<endl;
886  cout<<" fc = "<<fcAB<<", ("<< sigmaABMax/sigmaAB * fcAB <<","<< sigmaABMin/sigmaAB * fcAB <<")"<<endl;
887  // }
888 
889  //
890  // Filling part of the Eloss scenarii information
891  //
892  if(RaaCharm>0 ) {
893  gRAB_ElossHypothesis->SetPoint( hABbin, pt, RaaCharm ); // i, x, y
894  gRAB_ElossHypothesis->SetPointEXlow( hABbin, ptwidth);
895  gRAB_ElossHypothesis->SetPointEXhigh( hABbin, ptwidth);
896  gRAB_FeedDownSystematicsElossHypothesis->SetPoint( hABbin, pt, RaaCharm ); // i, x, y
897  gRAB_FeedDownSystematicsElossHypothesis->SetPointEXlow( hABbin, ptwidth);
898  gRAB_FeedDownSystematicsElossHypothesis->SetPointEXhigh( hABbin, ptwidth);
899  gRAB_GlobalSystematics->SetPoint( hABbin, pt, RaaCharm ); // i, x, y
900  gRAB_GlobalSystematics->SetPointEXlow(hABbin,0.4); gRAB_GlobalSystematics->SetPointEXhigh(hABbin,0.4);
901  }
902  }
903 
904  //
905  // Filling Eloss scenarii information
906  //
907  // trick in case not fine enough Rb hypothesis to cope with the min/max range
908  // if( RaaCharm>0 && ( (ElossHypo >= MinHypo && ElossHypo <=MaxHypo) || ElossHypo == ElossCentral[ hABbin ] ) && RaaBeauty<=MaxRb ) {
909  // by default better not use it, to monitor when this happens (could affect results)
910  if( RaaCharm>0 && ElossHypo >= MinHypo && ElossHypo <=MaxHypo && RaaBeauty<=MaxRb ) {
911 
912  Double_t Ehigh = ElossMax[ hABbin ] ;
913  Double_t Elow = ElossMin[ hABbin ] ;
914  if ( RaaCharm > Ehigh ) ElossMax[ hABbin ] = RaaCharm ;
915  if ( RaaCharm < Elow ) ElossMin[ hABbin ] = RaaCharm ;
916  if(printout && TMath::Abs(ptprintout-pt)<0.1) {
917  cout<<" Hypothesis " << ElossHypo << " sigma-AB "<< sigmaAB <<", Raa "<< RaaCharm <<", Raa Eloss max "<< ElossMax[hABbin] <<" Raa Eloss min "<< ElossMin[hABbin] << " Rb="<< RaaBeauty <<endl;
918  cout<<" Rb="<< RaaBeauty <<" max "<< RaaBeautyFDhigh <<" min "<< RaaBeautyFDlow <<endl;
919  }
920  Double_t fcEhigh = fcElossMax[ hABbin ] ;
921  Double_t fcElow = fcElossMin[ hABbin ] ;
922  if ( fcAB > fcEhigh ) fcElossMax[ hABbin ] = fcAB ;
923  if ( fcAB < fcElow ) fcElossMin[ hABbin ] = fcAB ;
924  Double_t FDEhigh = FDElossMax[ hABbin ];
925  Double_t FDEmin = FDElossMin[ hABbin ];
926  Double_t RFDhigh = RaaCharmFDhigh>RaaCharmFDlow ? RaaCharmFDhigh : RaaCharmFDlow;
927  Double_t RFDlow = RaaCharmFDlow<RaaCharmFDhigh ? RaaCharmFDlow : RaaCharmFDhigh;
928  if ( RFDhigh > FDEhigh ) FDElossMax[ hABbin ] = RFDhigh ;
929  if ( RFDlow < FDEmin ) FDElossMin[ hABbin ] = RFDlow ;
930  if(printout && TMath::Abs(ptprintout-pt)<0.1)
931  cout<<" Hypothesis " << ElossHypo << " sigma-AB "<< sigmaAB <<", Raa FD-max Eloss max "<< FDElossMax[hABbin] <<" Raa FD-min Eloss min "<< FDElossMin[hABbin] <<endl;
932  }
933 
934 
935  }
936 
937 
938  // Finish filling the y-uncertainties of the Eloss scenarii
939  for (Int_t ibin=0; ibin<=nbins; ibin++){
940  Double_t ipt=0., value =0.;
941  gRAB_ElossHypothesis->GetPoint(ibin,ipt,value);
942  if(ipt<=0) continue;
943  //
944  // Uncertainty on Raa due to the Eloss hypothesis
945  Double_t elossYhigh = TMath::Abs( ElossMax[ibin] - value );
946  Double_t elossYlow = TMath::Abs( value - ElossMin[ibin] );
947  gRAB_ElossHypothesis->SetPointEYhigh(ibin, elossYhigh );
948  gRAB_ElossHypothesis->SetPointEYlow(ibin, elossYlow );
949  gRAB_ElossHypothesis->SetPointEXhigh(ibin, 0.2);
950  gRAB_ElossHypothesis->SetPointEXlow(ibin, 0.2);
951  cout << " pt "<< ipt << " Raa "<< value <<" max "<< ElossMax[ibin] << " min " <<ElossMin[ibin] <<endl;
952  cout<<" Eloss syst +"<< elossYhigh <<" - "<< elossYlow <<endl;
953  // cout << " fc max "<< fcElossMax[ibin] << " fc min " <<fcElossMin[ibin] <<endl;
954  //
955  // Uncertainty on Raa due to the FD unc & Eloss hypothesis
956  Double_t fdElossEYhigh = TMath::Abs( FDElossMax[ibin] - value );
957  Double_t fdElossEYlow = TMath::Abs( value - FDElossMin[ibin] );
958  if(elossFDQuadSum){
959  Double_t fdEYhigh = gRAB_FeedDownSystematics->GetErrorYhigh(ibin);
960  fdElossEYhigh = TMath::Sqrt( elossYhigh*elossYhigh + fdEYhigh*fdEYhigh );
961  Double_t fdEYlow = gRAB_FeedDownSystematics->GetErrorYlow(ibin);
962  fdElossEYlow = TMath::Sqrt( elossYlow*elossYlow + fdEYlow*fdEYlow );
963  }
964  gRAB_FeedDownSystematicsElossHypothesis->SetPointEYhigh(ibin, fdElossEYhigh );
965  gRAB_FeedDownSystematicsElossHypothesis->SetPointEYlow(ibin, fdElossEYlow );
966  gRAB_FeedDownSystematicsElossHypothesis->SetPointEXhigh(ibin, 0.25);
967  gRAB_FeedDownSystematicsElossHypothesis->SetPointEXlow(ibin, 0.25);
968  cout<<" FD & Eloss syst +"<< fdElossEYhigh <<" - "<< fdElossEYlow
969  <<" = + "<< fdElossEYhigh/value <<" - "<< fdElossEYlow/value <<" %" <<endl;
970  //
971  // All uncertainty on Raa (FD unc & Eloss + data)
972  Double_t systdatal = gRAB_DataSystematics->GetErrorYlow(ibin);
973  Double_t systdatah = gRAB_DataSystematics->GetErrorYhigh(ibin);
974  Double_t systgbhUnc = TMath::Sqrt( systdatah*systdatah + fdElossEYhigh*fdElossEYhigh );
975  Double_t systgblUnc = TMath::Sqrt( systdatal*systdatal + fdElossEYlow*fdElossEYlow );
976  gRAB_GlobalSystematics->SetPointEYhigh(ibin,systgbhUnc);
977  gRAB_GlobalSystematics->SetPointEYlow(ibin,systgblUnc);
978  cout<<" Data syst +"<< systdatah <<" - "<< systdatal <<" = + "<< systdatah/value <<" - " << systdatal/value << " % "<<endl;
979  cout<<" Global syst +"<< systgbhUnc <<" - "<< systgblUnc << " = + "<< systgbhUnc/value <<" - "<< systgblUnc/value << " %" <<endl;
980  //
981  }
982 
983  cout<<endl<<" Calculation finished, now drawing"<<endl<<endl;
984 
985 
986  gROOT->SetStyle("Plain");
987  gStyle->SetPalette(1);
988  gStyle->SetOptStat(0);
989 
990 
991  TCanvas *cRABvsRb = new TCanvas("RABvsRb","RAB vs Rb");
992  hRABvsRb->Draw("colz");
993  cRABvsRb->Update();
994 
995 // TCanvas *cRABvsRbvsPt = new TCanvas("cRABvsRbvsPt","RAB vs Rb vs pt");
996 // hRABCharmVsRBeautyVsPt->Draw("lego3z");
997 // cRABvsRbvsPt->Update();
998 
999 
1000  cout<< " Drawing feed-down contribution"<<endl;
1001  TCanvas *cRABvsRbFDl = new TCanvas("RABvsRbFDl","RAB vs Rb (FD low)");
1002  hRABvsRbFDlow->Draw("cont4z");
1003  cRABvsRbFDl->Update();
1004  TCanvas *cRABvsRbFDh = new TCanvas("RABvsRbFDh","RAB vs Rb (FD high)");
1005  hRABvsRbFDhigh->Draw("cont4z");
1006  cRABvsRbFDh->Update();
1007 
1008  TCanvas * cSigmaABptEloss = new TCanvas("cSigmaABptEloss","SigmaAB vs pt, Eloss hypothesis");
1009  TH1D *hSigmaABEloss00= new TH1D("hSigmaABEloss00","hSigmaABEloss00",nbins,limits);
1010  TH1D *hSigmaABEloss05= new TH1D("hSigmaABEloss05","hSigmaABEloss05",nbins,limits);
1011  TH1D *hSigmaABEloss10= new TH1D("hSigmaABEloss10","hSigmaABEloss10",nbins,limits);
1012  TH1D *hSigmaABEloss15= new TH1D("hSigmaABEloss15","hSigmaABEloss15",nbins,limits);
1013  TH1D *hSigmaABEloss20= new TH1D("hSigmaABEloss20","hSigmaABEloss20",nbins,limits);
1014 
1015  delete [] limits;
1016  delete [] binwidths;
1017 
1018  for (Int_t i=0; i<=nSigmaAB->GetEntriesFast(); i++) {
1019  nSigmaAB->GetEntry(i);
1020  if (fdMethod==kfc) {
1021  if( TMath::Abs(Rcb-0.015)<0.009 ) hSigmaABEloss00->Fill(pt,sigmaAB);
1022  if( TMath::Abs(Rcb-0.5)<0.009 ) hSigmaABEloss05->Fill(pt,sigmaAB);
1023  if( TMath::Abs(Rcb-1.0)<0.009 ) hSigmaABEloss10->Fill(pt,sigmaAB);
1024  if( TMath::Abs(Rcb-1.5)<0.009 ) hSigmaABEloss15->Fill(pt,sigmaAB);
1025  if( TMath::Abs(Rcb-2.0)<0.009 ) hSigmaABEloss20->Fill(pt,sigmaAB);
1026  }
1027  else if (fdMethod==kNb) {
1028  if( TMath::Abs(Rb-0.015)<0.009 ) hSigmaABEloss00->Fill(pt,sigmaAB);
1029  if( TMath::Abs(Rb-0.5)<0.009 ) hSigmaABEloss05->Fill(pt,sigmaAB);
1030  if( TMath::Abs(Rb-1.0)<0.009 ) hSigmaABEloss10->Fill(pt,sigmaAB);
1031  if( TMath::Abs(Rb-1.5)<0.009 ) hSigmaABEloss15->Fill(pt,sigmaAB);
1032  if( TMath::Abs(Rb-2.0)<0.009 ) hSigmaABEloss20->Fill(pt,sigmaAB);
1033  }
1034  }
1035  hSigmaABEloss00->SetLineColor(2);
1036  hSigmaABEloss05->SetLineColor(3);
1037  hSigmaABEloss10->SetLineColor(4);
1038  hSigmaABEloss15->SetLineColor(kMagenta+1);
1039  hSigmaABEloss20->SetLineColor(kGreen+2);
1040  hSigmaABEloss00->SetMarkerStyle(22);
1041  hSigmaABEloss05->SetMarkerStyle(26);
1042  hSigmaABEloss10->SetMarkerStyle(20);
1043  hSigmaABEloss15->SetMarkerStyle(25);
1044  hSigmaABEloss20->SetMarkerStyle(21);
1045  if (fdMethod==kNb) {
1046  hSigmaABEloss05->Draw("ph");
1047  hSigmaABEloss10->Draw("phsame");
1048  hSigmaABEloss15->Draw("phsame");
1049  hSigmaABEloss20->Draw("phsame");
1050  }
1051  else {
1052  hSigmaABEloss20->Draw("p");
1053  hSigmaABEloss00->Draw("phsame");
1054  hSigmaABEloss05->Draw("phsame");
1055  hSigmaABEloss10->Draw("phsame");
1056  hSigmaABEloss15->Draw("phsame");
1057  hSigmaABEloss20->Draw("phsame");
1058  }
1059  TLegend *legrcb = new TLegend(0.8,0.8,0.95,0.9);
1060  legrcb->SetFillColor(0);
1061  legrcb->AddEntry(hSigmaABEloss00,"Rc/b=0.0","lp");
1062  legrcb->AddEntry(hSigmaABEloss05,"Rc/b=0.5","lp");
1063  legrcb->AddEntry(hSigmaABEloss10,"Rc/b=1.0","lp");
1064  legrcb->AddEntry(hSigmaABEloss15,"Rc/b=1.5","lp");
1065  legrcb->AddEntry(hSigmaABEloss20,"Rc/b=2.0","lp");
1066  legrcb->Draw();
1067  cSigmaABptEloss->Update();
1068 
1069 
1070  TCanvas * cRABptEloss = new TCanvas("cRABptEloss","RAB vs pt, Eloss hypothesis");
1071  hRABEloss00->SetLineColor(2);
1072  hRABEloss05->SetLineColor(3);
1073  hRABEloss10->SetLineColor(4);
1074  hRABEloss15->SetLineColor(kMagenta+1);
1075  hRABEloss20->SetLineColor(kGreen+2);
1076  hRABEloss00->SetMarkerStyle(22);
1077  hRABEloss05->SetMarkerStyle(26);
1078  hRABEloss10->SetMarkerStyle(20);
1079  hRABEloss15->SetMarkerStyle(25);
1080  hRABEloss20->SetMarkerStyle(21);
1081  if (fdMethod==kNb) {
1082  hRABEloss05->Draw("ph");
1083  hRABEloss10->Draw("phsame");
1084  hRABEloss15->Draw("phsame");
1085  hRABEloss20->Draw("phsame");
1086  }
1087  else {
1088  hRABEloss20->Draw("p");
1089  hRABEloss00->Draw("phsame");
1090  hRABEloss05->Draw("phsame");
1091  hRABEloss10->Draw("phsame");
1092  hRABEloss15->Draw("phsame");
1093  hRABEloss20->Draw("phsame");
1094  }
1095  legrcb = new TLegend(0.8,0.8,0.95,0.9);
1096  legrcb->SetFillColor(0);
1097  if (fdMethod==kfc) {
1098  legrcb->AddEntry(hRABEloss00,"Rc/b=0.0","lp");
1099  legrcb->AddEntry(hRABEloss05,"Rc/b=0.5","lp");
1100  legrcb->AddEntry(hRABEloss10,"Rc/b=1.0","lp");
1101  legrcb->AddEntry(hRABEloss15,"Rc/b=0.5","lp");
1102  legrcb->AddEntry(hRABEloss20,"Rc/b=2.0","lp");
1103  }
1104  else if (fdMethod==kNb) {
1105  legrcb->AddEntry(hRABEloss00,"Rb=0.0","lp");
1106  legrcb->AddEntry(hRABEloss05,"Rb=0.5","lp");
1107  legrcb->AddEntry(hRABEloss10,"Rb=1.0","lp");
1108  legrcb->AddEntry(hRABEloss15,"Rb=0.5","lp");
1109  legrcb->AddEntry(hRABEloss20,"Rb=2.0","lp");
1110  }
1111  legrcb->Draw();
1112  cRABptEloss->Update();
1113 
1114 
1115  cout<< " Drawing summary results"<<endl;
1116  TCanvas * cRABpt = new TCanvas("cRABpt","RAB vs pt, no hypothesis");
1117  hRABEloss10->Draw("");
1118  cRABpt->Update();
1119 
1120  TCanvas * cRABptFDUnc = new TCanvas("cRABptFDUnc","RAB vs pt, FD Uncertainties");
1121  hRABvsRbFDlow_proj->Draw("");
1122  hRABEloss10->Draw("phsame");
1123  hRABvsRbFDhigh_proj->SetLineColor(kMagenta+1);
1124  hRABvsRbFDhigh_proj->Draw("same");
1125  hRABvsRbFDlow_proj->SetLineColor(kGreen+2);
1126  hRABvsRbFDlow_proj->Draw("same");
1127  legrcb = new TLegend(0.8,0.8,0.95,0.9);
1128  legrcb->SetFillColor(0);
1129  legrcb->AddEntry(hRABEloss10,"FD Central","lp");
1130  legrcb->AddEntry(hRABvsRbFDhigh_proj,"FD Upper unc.","l");
1131  legrcb->AddEntry(hRABvsRbFDlow_proj,"FD Lower unc.","l");
1132  legrcb->Draw();
1133  cRABptFDUnc->Update();
1134 
1135  TCanvas *RaaPlot = new TCanvas("RaaPlot","RAB vs pt, plot all");
1136  RaaPlot->SetTopMargin(0.085);
1137  RaaPlot->SetBottomMargin(0.1);
1138  RaaPlot->SetTickx();
1139  RaaPlot->SetTicky();
1140  TH2D *hRaaCanvas = new TH2D("hRaaCanvas"," R_{AB}(c) vs p_{T} (no Eloss hypothesis); p_{t} [GeV/c] ; R_{AA} prompt D",40,0.,40.,100,0.,3.0);
1141  hRaaCanvas->GetXaxis()->SetTitleSize(0.05);
1142  hRaaCanvas->GetXaxis()->SetTitleOffset(0.9);
1143  hRaaCanvas->GetYaxis()->SetTitleSize(0.05);
1144  hRaaCanvas->GetYaxis()->SetTitleOffset(0.9);
1145  hRaaCanvas->Draw();
1146  gRAB_Norm->SetFillStyle(1001);
1147  gRAB_Norm->SetFillColor(kGray+2);
1148  gRAB_Norm->Draw("2");
1149  TLine *line = new TLine(0.0172415,1.0,40.,1.0);
1150  line->SetLineStyle(2);
1151  line->Draw();
1152  hRABvsPt->SetMarkerColor(kBlue);
1153  hRABvsPt->SetMarkerColor(kBlue);
1154  hRABvsPt->SetMarkerStyle(21);
1155  hRABvsPt->SetMarkerSize(1.1);
1156  hRABvsPt->SetLineWidth(2);
1157  hRABvsPt->Draw("psame");
1158  gRAB_DataSystematics->SetLineColor(kBlue);
1159  gRAB_DataSystematics->SetLineWidth(3);
1160  gRAB_DataSystematics->SetLineWidth(2);
1161  gRAB_DataSystematics->SetFillColor(kRed);
1162  gRAB_DataSystematics->SetFillStyle(0);
1163  gRAB_DataSystematics->Draw("2");
1164  gRAB_FeedDownSystematics->SetFillColor(kViolet+1);
1165  gRAB_FeedDownSystematics->SetFillStyle(1001);
1166  gRAB_FeedDownSystematics->Draw("2");
1167  gRAB_ElossHypothesis->SetLineColor(kMagenta-7);
1168  gRAB_ElossHypothesis->SetFillColor(kMagenta-7);
1169  gRAB_ElossHypothesis->SetFillStyle(1001);
1170  gRAB_ElossHypothesis->Draw("2");
1171  hRABvsPt->Draw("psame");
1172  gRAB_DataSystematics->Draw("2");
1173  legrcb = new TLegend(0.5517241,0.6504237,0.8520115,0.8728814,NULL,"brNDC");
1174  legrcb->SetBorderSize(0);
1175  legrcb->SetTextSize(0.03389831);
1176  legrcb->SetLineColor(1);
1177  legrcb->SetLineStyle(1);
1178  legrcb->SetLineWidth(1);
1179  legrcb->SetFillColor(0);
1180  legrcb->SetFillStyle(1001);
1181  if(cc==k020) legrcb->AddEntry(hRABvsPt,"R_{AA} 0-20% CC","pe");
1182  else if(cc==k4080) legrcb->AddEntry(hRABvsPt,"R_{AA} 40-80% CC","pe");
1183  else legrcb->AddEntry(hRABvsPt,"R_{AA} and stat. unc.","pe");
1184  legrcb->AddEntry(gRAB_DataSystematics,"Syst. from data","f");
1185  legrcb->AddEntry(gRAB_ElossHypothesis,"Syst. from R_{AA}(B)","f");
1186  legrcb->AddEntry(gRAB_FeedDownSystematics,"Syst. from B feed-down","f");
1187  legrcb->Draw();
1188  TLatex* tc;
1189  TString system = "Pb-Pb #sqrt{s_{NN}}=2.76 TeV";
1190  if( cc==kpPb0100 || cc==kpPb020 || cc==kpPb2040 || cc==kpPb4060 || cc==kpPb60100 ) system = "p-Pb #sqrt{s_{NN}}=5.023 TeV";
1191  if(decay==1) tc =new TLatex(0.18,0.82,Form("D^{0}, %s ",system.Data()));
1192  else if(decay==2) tc =new TLatex(0.18,0.82,Form("D^{+}, %s ",system.Data()));
1193  else if(decay==3) tc =new TLatex(0.18,0.82,Form("D^{*+}, %s ",system.Data()));
1194  else if(decay==4) tc =new TLatex(0.18,0.82,Form("D_{s}^{+}, %s ",system.Data()));
1195  else tc =new TLatex(0.18,0.82,Form("any (?) D meson, %s ",system.Data()));
1196  tc->SetNDC();
1197  tc->SetTextSize(0.038);
1198  tc->SetTextFont(42);
1199  tc->Draw();
1200  RaaPlot->Update();
1201 
1202 
1203  TCanvas *RaaPlotFDEloss = new TCanvas("RaaPlotFDEloss","RAB vs pt, plot FD & ElossUnc");
1204  RaaPlotFDEloss->SetTopMargin(0.085);
1205  RaaPlotFDEloss->SetBottomMargin(0.1);
1206  hRaaCanvas->Draw();
1207  line->Draw();
1208  hRABvsPt->Draw("psame");
1209  gRAB_FeedDownSystematics->SetFillColor(kViolet+1);
1210  gRAB_FeedDownSystematics->SetFillStyle(1001);
1211  gRAB_FeedDownSystematics->Draw("2");
1212  gRAB_ElossHypothesis->SetLineColor(kMagenta-7);
1213  gRAB_ElossHypothesis->SetFillColor(kMagenta-7);
1214  gRAB_ElossHypothesis->SetFillStyle(1001);
1215  gRAB_ElossHypothesis->Draw("2");
1216  gRAB_FeedDownSystematicsElossHypothesis->SetLineColor(kBlack);
1217  gRAB_FeedDownSystematicsElossHypothesis->SetFillStyle(0);
1218  gRAB_FeedDownSystematicsElossHypothesis->SetFillColor(kViolet+1);
1219  gRAB_FeedDownSystematicsElossHypothesis->Draw("2");
1220  hRABvsPt->Draw("psame");
1221  legrcb = new TLegend(0.6,0.6,0.9,0.9);
1222  legrcb->SetBorderSize(0);
1223  legrcb->SetTextSize(0.03389831);
1224  legrcb->SetLineColor(1);
1225  legrcb->SetLineStyle(1);
1226  legrcb->SetLineWidth(1);
1227  legrcb->SetFillColor(0);
1228  legrcb->SetFillStyle(1001);
1229  legrcb->AddEntry(hRABvsPt,"R_{PbPb} and stat. unc.","pe");
1230  legrcb->AddEntry(gRAB_ElossHypothesis,"Energy loss syst.","f");
1231  legrcb->AddEntry(gRAB_FeedDownSystematics,"Feed down syst.","f");
1232  legrcb->AddEntry(gRAB_FeedDownSystematicsElossHypothesis,"Feed down & Eloss syst.","f");
1233  legrcb->Draw();
1234  RaaPlotFDEloss->Update();
1235 
1236 
1237  TCanvas *RaaPlotGlob = new TCanvas("RaaPlotGlob","RAB vs pt, plot Global unc");
1238  RaaPlotGlob->SetTopMargin(0.085);
1239  RaaPlotGlob->SetBottomMargin(0.1);
1240  RaaPlotGlob->SetTickx();
1241  RaaPlotGlob->SetTicky();
1242  hRaaCanvas->Draw();
1243  line->Draw();
1244  hRABvsPt->Draw("psame");
1245  gRAB_DataSystematics->Draw("2");
1246  gRAB_FeedDownSystematicsElossHypothesis->Draw("2");
1247  gRAB_GlobalSystematics->SetLineColor(kRed);
1248  gRAB_GlobalSystematics->SetLineWidth(2);
1249  gRAB_GlobalSystematics->SetFillColor(kRed);
1250  gRAB_GlobalSystematics->SetFillStyle(3002);
1251  gRAB_GlobalSystematics->Draw("2");
1252  hRABvsPt->Draw("psame");
1253  legrcb = new TLegend(0.6,0.6,0.9,0.9);
1254  legrcb->SetBorderSize(0);
1255  legrcb->SetTextSize(0.03389831);
1256  legrcb->SetLineColor(1);
1257  legrcb->SetLineStyle(1);
1258  legrcb->SetLineWidth(1);
1259  legrcb->SetFillColor(0);
1260  legrcb->SetFillStyle(1001);
1261  legrcb->AddEntry(hRABvsPt,"R_{PbPb} and stat. unc.","pe");
1262  legrcb->AddEntry(gRAB_DataSystematics,"Data syst.","f");
1263  legrcb->AddEntry(gRAB_FeedDownSystematicsElossHypothesis,"Feed down & Eloss syst.","f");
1264  legrcb->AddEntry(gRAB_GlobalSystematics,"Global syst.","f");
1265  legrcb->Draw();
1266  RaaPlotGlob->Update();
1267 
1268 
1269 
1270  TCanvas *RaaPlotSimple = new TCanvas("RaaPlotSimple","RAB vs pt, plot Simple unc");
1271  RaaPlotSimple->SetTopMargin(0.085);
1272  RaaPlotSimple->SetBottomMargin(0.1);
1273  RaaPlotSimple->SetTickx();
1274  RaaPlotSimple->SetTicky();
1275  hRaaCanvas->Draw();
1276  line->Draw();
1277  hRABvsPt->Draw("psame");
1278  gRAB_GlobalSystematics->SetLineColor(kBlue);
1279  gRAB_GlobalSystematics->SetLineWidth(2);
1280  gRAB_GlobalSystematics->SetFillStyle(0);
1281  gRAB_GlobalSystematics->Draw("2");
1282  gRAB_Norm->Draw("2");
1283  hRABvsPt->Draw("psame");
1284  legrcb = new TLegend(0.5991379,0.6949153,0.8534483,0.8559322,NULL,"brNDC");
1285  legrcb->SetBorderSize(0);
1286  legrcb->SetTextSize(0.03389831);
1287  legrcb->SetLineColor(1);
1288  legrcb->SetLineStyle(1);
1289  legrcb->SetLineWidth(1);
1290  legrcb->SetFillColor(0);
1291  legrcb->SetFillStyle(1001);
1292  if(cc==k020) legrcb->AddEntry(hRABvsPt,"R_{AA} 0-20% CC","pe");
1293  else if(cc==k4080) legrcb->AddEntry(hRABvsPt,"R_{AA} 40-80% CC","pe");
1294  else legrcb->AddEntry(hRABvsPt,"R_{AA} and stat. unc.","pe");
1295  legrcb->AddEntry(gRAB_GlobalSystematics,"Systematics","f");
1296  legrcb->Draw();
1297  tc->Draw();
1298  RaaPlotSimple->Update();
1299 
1300 
1301  TCanvas *c = new TCanvas("c","");
1302  systematicsAB->DrawErrors();
1303  c->Update();
1304 
1305  TCanvas *cStatUnc = new TCanvas("cStatUnc","stat unc");
1306  cStatUnc->Divide(2,2);
1307  cStatUnc->cd(1);
1308  fhStatUncEffcSigmaAB_Raa->Draw("e");
1309  cStatUnc->cd(2);
1310  fhStatUncEffbSigmaAB_Raa->Draw("e");
1311  cStatUnc->cd(3);
1312  fhStatUncEffcFDAB_Raa->Draw("e");
1313  cStatUnc->cd(4);
1314  fhStatUncEffbFDAB_Raa->Draw("e");
1315  cStatUnc->Update();
1316 
1317  //
1318  // Write the information to a file
1319  //
1320  cout<<endl<< " Save results in the output file"<<endl<<endl;
1321  TFile * out = new TFile(outfile,"recreate");
1322 
1323  ntupleRAB->Write();
1324  hRABvsRcb->Write();
1325  hRABvsRb->Write();
1326  hRABCharmVsRBeautyVsPt->Write();
1327  for(Int_t j=0; j<=nbins; j++) hRCharmVsRBeauty[j]->Write();
1328  // for(Int_t j=0; j<=nbins; j++) hRCharmVsElossHypo[j]->Write();
1329  hRABvsPt->Write();
1330  hRABvsPt_DataSystematics->Write();
1331  gRAB_ElossHypothesis->Write();
1332  gRAB_FeedDownSystematics->Write();
1333  gRAB_fcFeedDownOnly->Write();
1334  gRAB_DataSystematics->Write();
1335  gRAB_DataSystematicsPP->Write();
1336  gSigmaPPSystTheory->Write();
1337  gRAB_DataSystematicsAB->Write();
1338  gRAB_Norm->Write();
1339  gRAB_FeedDownSystematicsElossHypothesis->Write();
1340  gRAB_GlobalSystematics->Write();
1341  if(isScaledAndExtrapRef && hCombinedReferenceFlag) hCombinedReferenceFlag->Write();
1342  systematicsPP->SetName("AliHFSystErrPP");
1343  systematicsPP->Write();
1344  systematicsAB->SetName("AliHFSystErrAA");
1345  systematicsAB->Write();
1346  out->Write();
1347 
1348 }
1349 
1350 //____________________________________________________________
1351 Bool_t PbPbDataSyst(AliHFSystErr *syst, Double_t pt, Int_t cc, Double_t &dataSystUp, Double_t &dataSystDown)
1352 {
1353 
1354  Double_t err=0., errUp=1., errDown=1.;
1355  Bool_t isOk=false;
1356  Double_t pidunc=0.;
1357 
1358  err = syst->GetTotalSystErr(pt)*syst->GetTotalSystErr(pt);
1359  errUp = err ;
1360  errDown = err ;
1361 
1362  // Apply an asymetric PID uncertainty for 2010 PbPb data only
1363  if( syst->GetRunNumber()==10 && syst->GetCollisionType()==1 ) {
1364  if( cc==k07half || cc==k020 || cc==k010 || cc==k1020 || cc==k2040 ) {
1365  if(pt<6) pidunc = 0.15;
1366  else pidunc = 0.05;
1367  errUp = err + pidunc*pidunc - syst->GetPIDEffErr(pt)*syst->GetPIDEffErr(pt);
1368  isOk = true;
1369  }
1370  else if ( cc==k3050 || cc==k4080 || cc==k4060 || cc==k6080 ){
1371  if(pt<3.1) pidunc = 0.10;
1372  else pidunc = 0.05;
1373  errUp = err + pidunc*pidunc - syst->GetPIDEffErr(pt)*syst->GetPIDEffErr(pt);
1374  isOk = true;
1375  }
1376  }
1377  else { isOk = true; }
1378 
1379  dataSystUp = TMath::Sqrt(errUp);
1380  dataSystDown = TMath::Sqrt(errDown);
1381 
1382  return isOk;
1383 }
particularity
Definition: HFPtSpectrum.C:48
const Color_t cc[]
Definition: DrawKs.C:1
double Double_t
Definition: External.C:58
BFDSubtrMethod
Definition: HFPtSpectrum.C:45
energy
Definition: HFPtSpectrum.C:44
centrality
void SetCentrality(TString centrality)
Definition: AliHFSystErr.h:60
TCanvas * c
Definition: TestFitELoss.C:172
void SetIsPbPb2010EnergyScan(Bool_t flag)
Definition: AliHFSystErr.h:86
Bool_t elossFDQuadSum
Int_t GetRunNumber() const
Definition: AliHFSystErr.h:48
void SetIsBDTAnalysis(Bool_t flag)
Definition: AliHFSystErr.h:80
int Int_t
Definition: External.C:63
Definition: External.C:204
float Float_t
Definition: External.C:68
RaavsEP
Definition: HFPtSpectrum.C:46
centestimator
Definition: External.C:252
Definition: External.C:228
Definition: External.C:212
rapidity
Definition: HFPtSpectrum.C:47
void SetIsPass4Analysis(Bool_t flag)
Definition: AliHFSystErr.h:72
void HFPtSpectrumRaa(const char *ppfile="HFPtSpectrum_D0Kpi_method2_rebinnedth_230311_newsigma.root", const char *ABfile="HFPtSpectrum_D0Kpi_PbPbcuts_method2_rebinnedth_230311_newsigma.root", const char *outfile="HFPtSpectrumRaa.root", Int_t decay=1, Double_t sigmaABCINT1B=54.e9, Int_t fdMethod=kNb, Int_t cc=kpp, Int_t Energy=k276, Double_t MinHypo=1./3., Double_t MaxHypo=3.0, Double_t MaxRb=6.0, Bool_t isRbHypo=false, Double_t CentralHypo=1.0, Int_t ccestimator=kV0M, Bool_t isUseTaaForRaa=true, const char *shadRbcFile="", Int_t nSigmaShad=3.0, Int_t isRaavsEP=kPhiIntegrated, Bool_t isScaledAndExtrapRef=kFALSE, Int_t rapiditySlice=kdefault, Int_t analysisSpeciality=kTopological)
Double_t GetTotalSystErr(Double_t pt, Double_t feeddownErr=0) const
void SetIsLowPtAnalysis(Bool_t flag)
Definition: AliHFSystErr.h:68
Double_t ExtractFDSyst(Double_t total, Double_t fd)
Double_t ptprintout
decay
Definition: HFPtSpectrum.C:41
void Init(Int_t decay)
Function to initialize the variables/histograms.
void SetIspPb2011RapidityScan(Bool_t flag)
Definition: AliHFSystErr.h:96
Bool_t printout
void SetCollisionType(Int_t type)
Definition: AliHFSystErr.h:51
Double_t NormABUnc
Bool_t PbPbDataSyst(AliHFSystErr *syst, Double_t pt, Int_t cc, Double_t &dataSystUp, Double_t &dataSystDown)
const Int_t nbins
bool Bool_t
Definition: External.C:53
Int_t FindGraphBin(TGraphAsymmErrors *gr, Double_t pt)
void SetRapidity(TString rapidity)
Settings of rapidity ranges for pPb 0-100% CC.
Definition: AliHFSystErr.h:92
Int_t GetCollisionType() const
Definition: AliHFSystErr.h:57
void SetRunNumber(Int_t number)
Definition: AliHFSystErr.h:44
Double_t NormPPUnc
Double_t GetPIDEffErr(Double_t pt) const