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