PWGHF/vertexingHF/macros/HFPtSpectrumRaa.C

   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$ */
  24 ///////////////////////////////////////////////////////////////////////////////
  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 //
  51 ///////////////////////////////////////////////////////////////////////////////
  52
  53 enum centrality{ kpp, k07half, kpPb0100, k010, k1020, k020, k2040, k2030, k3040, k4050, k3050, k5060, k4060, k6080, k4080, k80100 };
  54 enum energy{ k276, k55 };
  55 enum BFDSubtrMethod { kfc, kNb };
  56 enum RaavsEP {kPhiIntegrated, kInPlane, kOutOfPlane};
  57
  58 Bool_t printout = false;
  59 Double_t NormPPUnc = 0.035;
  60 Bool_t elossFDQuadSum = true;
  61
  62 //____________________________________________________________
  63 Bool_t PbPbDataSyst(AliHFSystErr *syst, Double_t pt, Int_t cc, Double_t &dataSystUp, Double_t &dataSystDown);
  64
  65 //____________________________________________________________
  66 Double_t ExtractFDSyst(Double_t total, Double_t fd) {
  67   // total^2 = data^2 + fd^2
  68   Double_t data2 = total*total - fd*fd ;
  69   return TMath::Sqrt( data2 );
  70 }
  71
  72 //____________________________________________________________
  73 Int_t FindGraphBin(TGraphAsymmErrors *gr, Double_t pt)
  74 {
  75   Int_t istart =0;
  76   Int_t npoints = gr->GetN();
  77   for(Int_t i=0; i<=npoints; i++){
  78     Double_t x=0.,y=0.;
  79     gr->GetPoint(i,x,y);
  80     if ( TMath::Abs ( x - pt ) < 0.4 ) {
  81       istart = i;
  82       break;
  83     }
  84   }
  85   return istart;
  86 }
  87
  88
  89 //
  90 //
  91 // R_AB =  [ ( dsigma/dpt )_AB / sigma_AB ]  / <TAB> *  ( dsigma/dpt )_pp
  92 //
  93 //
  94 //____________________________________________________________
  95 void HFPtSpectrumRaa(const char *ppfile="HFPtSpectrum_D0Kpi_method2_rebinnedth_230311_newsigma.root",
  96                      const char *ABfile="HFPtSpectrum_D0Kpi_PbPbcuts_method2_rebinnedth_230311_newsigma.root",
  97                      const char *outfile="HFPtSpectrumRaa.root",
  98                      Int_t decay=1,
  99                      Double_t sigmaABCINT1B=54.e9,
 100                      Int_t fdMethod = kNb, Int_t cc=kpp, Int_t Energy=k276,
 101                      Double_t MinHypo=1./3., Double_t MaxHypo=3.0, Double_t MaxRb=6.0,
 102                      Bool_t isRbHypo=false, Double_t CentralHypo = 1.0,
 103                      Int_t isRaavsEP=kPhiIntegrated, Bool_t isScaledAndExtrapRef=kFALSE)
 104 {
 105
 106   gROOT->Macro("$ALICE_ROOT/PWGHF/vertexingHF/macros/LoadLibraries.C");
 107
 108   //
 109   // Defining the TAB values for the given centrality class
 110   //
 111   Double_t Tab = 1., TabSyst = 0.;
 112   if ( Energy!=k276 ) {
 113     printf("\n The Tab values for this cms energy have not yet been implemented, please do it ! \n");
 114     return;
 115   }
 116   if ( cc == kpp ){
 117     Tab = 1.;
 118   }
 119   // Values from Alberica's twiki:
 120   //   https://twiki.cern.ch/twiki/bin/viewauth/ALICE/CentStudies
 121   if ( cc == k07half ) {
 122     Tab = 24.81; TabSyst = 0.8037;
 123   } else if ( cc == k010 ) {
 124     Tab = 23.48; TabSyst = 0.97;
 125   } else if ( cc == k1020 ) {
 126     Tab = 14.4318; TabSyst = 0.5733;
 127   } else if ( cc == k020 ) {
 128     Tab = 18.93; TabSyst = 0.74;
 129   } else if ( cc == k2040 ) {
 130     Tab = 6.86; TabSyst = 0.28;
 131   } else if ( cc == k2030 ) {
 132     Tab = 8.73769; TabSyst = 0.370219;
 133   } else if ( cc == k3040 ) {
 134     Tab = 5.02755; TabSyst = 0.22099;
 135   } else if ( cc == k4050 ) {
 136     Tab = 2.68327; TabSyst = 0.137073;
 137   } else if ( cc == k3050 ) {
 138     Tab = 3.87011; TabSyst = 0.183847;
 139   } else if ( cc == k4060 ) {
 140     Tab = 2.00;  TabSyst= 0.11;
 141   } else if ( cc == k4080 ) {
 142     Tab = 1.20451; TabSyst = 0.071843;
 143   } else if ( cc == k5060 ) {
 144     Tab = 1.32884; TabSyst = 0.0929536;
 145   } else if ( cc == k6080 ) {
 146     Tab = 0.419; TabSyst = 0.033;
 147   } else if ( cc == k80100 ){
 148     Tab = 0.0690; TabSyst = 0.0062;
 149   }
 150
 151   // pPb Glauber (A. Toia)
 152   // https://twiki.cern.ch/twiki/bin/viewauth/ALICE/PACentStudies#Glauber_Calculations_with_sigma
 153   else if( cc == kpPb0100 ){
 154     tab = 0.098334; tabUnc = 0.0070679;
 155   }
 156
 157   //
 158   // Reading the pp file
 159   //
 160   TFile * ppf = new TFile(ppfile,"read");
 161   TH1D * hSigmaPP;// = (TH1D*)ppf->Get("fhScaledData");
 162   TGraphAsymmErrors * gSigmaPPSyst;// = (TGraphAsymmErrors*)ppf->Get("gScaledData");
 163   TGraphAsymmErrors * gSigmaPPSystData = (TGraphAsymmErrors*)ppf->Get("gScaledDataSystData");
 164   TGraphAsymmErrors * gSigmaPPSystTheory = (TGraphAsymmErrors*)ppf->Get("gScaledDataSystExtrap");
 165   TGraphAsymmErrors * gSigmaPPSystFeedDown = (TGraphAsymmErrors*)ppf->Get("gScaledDataSystFeedDown");
 166   TH1I * hCombinedReferenceFlag;
 167   TGraphAsymmErrors * gReferenceFdSyst;
 168   if(isScaledAndExtrapRef){
 169     hCombinedReferenceFlag = (TH1I*)ppf->Get("hCombinedReferenceFlag");
 170     hSigmaPP = (TH1D*)ppf->Get("hReference");
 171     gSigmaPPSyst = (TGraphAsymmErrors*)ppf->Get("gReferenceSyst");
 172     gReferenceFdSyst = (TGraphAsymmErrors*)ppf->Get("gReferenceFdSyst");
 173   } else {
 174     hSigmaPP = (TH1D*)ppf->Get("fhScaledData");
 175     gSigmaPPSyst = (TGraphAsymmErrors*)ppf->Get("gScaledData");
 176   }
 177
 178   // Call the systematics uncertainty class for a given decay
 179   AliHFSystErr *systematicsPP = new AliHFSystErr();
 180   systematicsPP->Init(decay);
 181
 182   //
 183   // Reading the AB collisions file
 184   //
 185   TFile * ABf = new TFile(ABfile,"read");
 186   TH1D *hSigmaAB = (TH1D*)ABf->Get("histoSigmaCorr");
 187   TH2D *hSigmaABRcb = (TH2D*)ABf->Get("histoSigmaCorrRcb");
 188   TGraphAsymmErrors * gSigmaABSyst = (TGraphAsymmErrors*)ABf->Get("gSigmaCorr");
 189   TGraphAsymmErrors * gSigmaABSystFeedDown = (TGraphAsymmErrors*)ABf->Get("gSigmaCorrConservative");
 190   TNtuple * nSigmaAB = (TNtuple*)ABf->Get("fnSigma");
 191   //
 192   TH1D *hMassAB = (TH1D*)ABf->Get("hRECpt");
 193   TH1D *hDirectEffptAB = (TH1D*)ABf->Get("hDirectEffpt");
 194   TH1D *histofcAB = (TH1D*)ABf->Get("histofc");
 195   //
 196   TH1D* fhStatUncEffcSigmaAB = (TH1D*)ABf->Get("fhStatUncEffcSigma");
 197   TH1D* fhStatUncEffbSigmaAB = (TH1D*)ABf->Get("fhStatUncEffbSigma");
 198   TH1D* fhStatUncEffcFDAB = (TH1D*)ABf->Get("fhStatUncEffcFD");
 199   TH1D* fhStatUncEffbFDAB = (TH1D*)ABf->Get("fhStatUncEffbFD");
 200   //
 201   TH1D* fhStatUncEffcSigmaAB_Raa = (TH1D*)fhStatUncEffcSigmaAB->Clone("fhStatUncEffcSigmaAB_Raa");
 202   TH1D* fhStatUncEffbSigmaAB_Raa = (TH1D*)fhStatUncEffbSigmaAB->Clone("fhStatUncEffbSigmaAB_Raa");
 203   TH1D* fhStatUncEffcFDAB_Raa = (TH1D*)fhStatUncEffcFDAB->Clone("fhStatUncEffcFDAB_Raa");
 204   TH1D* fhStatUncEffbFDAB_Raa = (TH1D*)fhStatUncEffbFDAB->Clone("fhStatUncEffbFDAB_Raa");
 205   fhStatUncEffcSigmaAB_Raa->Reset();
 206   fhStatUncEffbSigmaAB_Raa->Reset();
 207   fhStatUncEffcFDAB_Raa->Reset();
 208   fhStatUncEffbFDAB_Raa->Reset();
 209   fhStatUncEffcSigmaAB_Raa->SetName("fhStatUncEffcSigmaAB_Raa");
 210   fhStatUncEffbSigmaAB_Raa->SetName("fhStatUncEffbSigmaAB_Raa");
 211   fhStatUncEffcFDAB_Raa->SetName("fhStatUncEffcFDAB_Raa");
 212   fhStatUncEffbFDAB_Raa->SetName("fhStatUncEffvFDAB_Raa");
 213
 214
 215   //
 216   // Call the systematics uncertainty class for a given decay
 217   AliHFSystErr *systematicsAB = new AliHFSystErr();
 218   systematicsAB->SetCollisionType(1);
 219   if ( cc == k07half ) systematicsAB->SetCentrality("07half");
 220   else if ( cc == k010 ) systematicsAB->SetCentrality("010");
 221   else if ( cc == k1020 ) systematicsAB->SetCentrality("1020");
 222   else if ( cc == k020 ) systematicsAB->SetCentrality("020");
 223   else if ( cc == k2040 || cc == k2030 || cc == k3040 ) {
 224     systematicsAB->SetCentrality("2040");
 225     systematicsAB->SetIsPbPb2010EnergyScan(true);
 226   }
 227   else if ( cc == k4060 || cc == k4050 || cc == k5060 ) systematicsAB->SetCentrality("4060");
 228   else if ( cc == k6080 ) systematicsAB->SetCentrality("6080");
 229   else if ( cc == k4080 ) systematicsAB->SetCentrality("4080");
 230   else if ( cc == k3050 ) {
 231     if (isRaavsEP == kPhiIntegrated) systematicsAB->SetCentrality("4080");
 232     else if (isRaavsEP == kInPlane) systematicsAB->SetCentrality("3050InPlane");
 233     else if (isRaavsEP == kOutOfPlane) systematicsAB->SetCentrality("3050OutOfPlane");
 234   }
 235   else {
 236     cout << " Systematics not yet implemented " << endl;
 237     return;
 238   }
 239   //
 240   if ( cc == kpPb0100 ){
 241     systematics->SetCollisionType(0);
 242     cout <<endl<<" Beware pPb systematics not yet implemented, using pp at 7 TeV !! "<<endl<<endl;
 243   }
 244   //
 245   systematicsAB->Init(decay);
 246
 247   //
 248   Int_t entries = nSigmaAB->GetEntries();
 249   Float_t pt=0., signal=0., Rb=0., Rcb=0., fcAB=0., yieldAB=0., sigmaAB=0., statUncSigmaAB=0., sigmaABMin=0.,sigmaABMax=0.;
 250   nSigmaAB->SetBranchAddress("pt",&pt);
 251   nSigmaAB->SetBranchAddress("Signal",&signal);
 252   if (fdMethod==kNb) nSigmaAB->SetBranchAddress("Rb",&Rb);
 253   else if (fdMethod==kfc) nSigmaAB->SetBranchAddress("Rcb",&Rcb);
 254   nSigmaAB->SetBranchAddress("fc",&fcAB);
 255   nSigmaAB->SetBranchAddress("Yield",&yieldAB);
 256   nSigmaAB->SetBranchAddress("Sigma",&sigmaAB);
 257   nSigmaAB->SetBranchAddress("SigmaStatUnc",&statUncSigmaAB);
 258   nSigmaAB->SetBranchAddress("SigmaMax",&sigmaABMax);
 259   nSigmaAB->SetBranchAddress("SigmaMin",&sigmaABMin);
 260
 261
 262   // define the binning
 263   Int_t nbins = hSigmaAB->GetNbinsX();
 264   Double_t binwidth = hSigmaAB->GetBinWidth(1);
 265   Double_t *limits = new Double_t[nbins+1];
 266   Double_t *binwidths = new Double_t[nbins];
 267   Double_t xlow=0.;
 268   for (Int_t i=1; i<=nbins; i++) {
 269     binwidth = hSigmaAB->GetBinWidth(i);
 270     xlow = hSigmaAB->GetBinLowEdge(i);
 271     limits[i-1] = xlow;
 272     binwidths[i-1] = binwidth;
 273   }
 274   limits[nbins] = xlow + binwidth;
 275   //
 276   // define the bins correspondence bw histos/files/graphs
 277   //
 278   //
 279   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.);
 280   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.);
 281   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.);
 282   Int_t nbinsHypo=800;//200;
 283   Double_t *limitsHypo = new Double_t[nbinsHypo+1];
 284   for(Int_t i=1; i<=nbinsHypo+1; i++) limitsHypo[i-1]= i*4./800.;
 285   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);
 286   TH2D *hRCharmVsRBeauty[nbins+1];
 287   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);
 288   TH2D *hRCharmVsElossHypo[nbins+1];
 289   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);
 290   //
 291   TH1D *hRABEloss00= new TH1D("hRABEloss00","hRABEloss00",nbins,limits);
 292   TH1D *hRABEloss05= new TH1D("hRABEloss05","hRABEloss05",nbins,limits);
 293   TH1D *hRABEloss10= new TH1D("hRABEloss10","hRABEloss10",nbins,limits);
 294   TH1D *hRABEloss15= new TH1D("hRABEloss15","hRABEloss15",nbins,limits);
 295   TH1D *hRABEloss20= new TH1D("hRABEloss20","hRABEloss20",nbins,limits);
 296   //
 297   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.);
 298   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.);
 299   //
 300   TH1D * hRABvsRbFDhigh_proj = new TH1D("hRABvsRbFDhigh_proj","hRABvsRbFDhigh_proj",nbins,limits);
 301   TH1D * hRABvsRbFDlow_proj = new TH1D("hRABvsRbFDlow_proj","hRABvsRbFDlow_proj",nbins,limits);
 302   //
 303   TNtuple *ntupleRAB=0x0 ;
 304   if (fdMethod==kNb) {
 305     ntupleRAB = new TNtuple("ntupleRAB","ntupleRAB (Nb)","pt:TAB:sigmaPP:sigmaAB:invyieldAB:invyieldABFDHigh:invyieldABFDLow:RABCharm:RABCharmFDHigh:RABCharmFDLow:RABBeauty:fc");
 306   } else if (fdMethod==kfc) {
 307     ntupleRAB = new TNtuple("ntupleRAB","ntupleRAB (fc)","pt:TAB:sigmaPP:sigmaAB:invyieldAB:invyieldABFDHigh:invyieldABFDLow:Rcb:RABCharm:RABCharmFDHigh:RABCharmFDLow:RABBeauty:RABBeautyFDHigh:RABBeautyFDLow:fc");
 308   }
 309   if(!ntupleRAB) printf("ERROR: Wrong method option");
 310
 311   TH1D * hYieldABvsPt = new TH1D("hYieldABvsPt"," Yield_{AB}(c) vs p_{T} (no Eloss hypothesis); p_{T} [GeV/c] ; Yield_{charm} ",nbins,limits);
 312   TH1D * hRABvsPt = new TH1D("hRABvsPt"," R_{AB}(c) vs p_{T} (no Eloss hypothesis); p_{T} [GeV/c] ; R_{charm} ",nbins,limits);
 313   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);
 314   TGraphAsymmErrors *gRAB_ElossHypothesis = new TGraphAsymmErrors(nbins+1);
 315   gRAB_ElossHypothesis->SetNameTitle("gRAB_ElossHypothesis","RAB Eloss systematics");
 316   TGraphAsymmErrors *gRAB_FeedDownSystematics = new TGraphAsymmErrors(nbins+1);
 317   gRAB_FeedDownSystematics->SetNameTitle("gRAB_FeedDownSystematics","RAB Feed-Down systematics");
 318   TGraphAsymmErrors *gRAB_fcFeedDownOnly = new TGraphAsymmErrors(nbins+1);
 319   gRAB_fcFeedDownOnly->SetNameTitle("gRAB_fcFeedDownOnly","RAB fc Feed-Down Only");
 320   TGraphAsymmErrors *gRAB_FeedDownSystematicsElossHypothesis = new TGraphAsymmErrors(nbins+1);
 321   gRAB_FeedDownSystematicsElossHypothesis->SetNameTitle("gRAB_FeedDownSystematicsElossHypothesis","RAB Feed-Down systematics considering Eloss hypothesis");
 322   TGraphAsymmErrors *gRAB_DataSystematics = new TGraphAsymmErrors(nbins+1);
 323   gRAB_DataSystematics->SetNameTitle("gRAB_DataSystematics","RAB Measurement (no FD, no Eloss) systematics");
 324   TGraphAsymmErrors *gRAB_DataSystematicsPP = new TGraphAsymmErrors(nbins+1);
 325   gRAB_DataSystematicsPP->SetNameTitle("gRAB_DataSystematicsPP","RAB Measurement PP meas. systematics (data+scaling)");
 326   TGraphAsymmErrors *gRAB_DataSystematicsAB = new TGraphAsymmErrors(nbins+1);
 327   gRAB_DataSystematicsAB->SetNameTitle("gRAB_DataSystematicsAB","RAB Measurement AB (no FD, no Eloss, no PP data) systematics");
 328   TGraphAsymmErrors *gRAB_GlobalSystematics = new TGraphAsymmErrors(nbins+1);
 329   gRAB_GlobalSystematics->SetNameTitle("gRAB_GlobalSystematics","RAB Measurement global (data, FD, Eloss) systematics");
 330   Double_t ElossMax[nbins+1], ElossMin[nbins+1];
 331   for(Int_t i=0; i<=nbins; i++) { ElossMax[i]=0.; ElossMin[i]=1.; }
 332   Double_t fcElossMax[nbins+1], fcElossMin[nbins+1];
 333   for(Int_t i=0; i<=nbins; i++) { fcElossMax[i]=0.; fcElossMin[i]=1.; }
 334   Double_t FDElossMax[nbins+1], FDElossMin[nbins+1];
 335   for(Int_t i=0; i<=nbins; i++) { FDElossMax[i]=0.; FDElossMin[i]=1.; }
 336
 337   TGraphAsymmErrors *gRAB_Norm = new TGraphAsymmErrors(1);
 338   gRAB_Norm->SetNameTitle("gRAB_Norm","RAB Normalization systematics (pp norm + Tab)");
 339   Double_t normUnc = TMath::Sqrt ( NormPPUnc*NormPPUnc + (TabSyst/Tab)*(TabSyst/Tab) );
 340   gRAB_Norm->SetPoint(1,0.5,1.);
 341   gRAB_Norm->SetPointError(1,0.25,0.25,normUnc,normUnc);
 342
 343   //
 344   // R_AB =  ( dN/dpt )_AB  / <Ncoll_AB> *  ( dN/dpt )_pp ; <Ncoll> = <Tab> * sigma_NN^inel
 345   // R_AB =  [ ( dsigma/dpt )_AB / sigma_AB ]  / <TAB> *  ( dsigma/dpt )_pp
 346   //
 347   Int_t istartPPfd=0, istartPPsyst=0, istartABfd=0, istartPPextr=0;
 348   Double_t yPPh=0., yPPl=0., yABh=0., yABl=0.;
 349   Double_t RaaCharm =0., RaaBeauty=0.;
 350   Double_t RaaCharmFDhigh = 0., RaaCharmFDlow = 0.;
 351   Double_t RaaBeautyFDhigh = 0., RaaBeautyFDlow = 0.;
 352   Double_t systUp=0., systLow=0., systPPUp=0., systPPLow=0., systABUp=0., systABLow=0.;
 353   //
 354   //
 355   // Search the central value of the energy loss hypothesis Rb = Rc (bin)
 356   //
 357   Double_t ElossCentral[nbins+1];
 358   for(Int_t i=0; i<=nbins; i++) { ElossCentral[i]=0.; }
 359   //
 360   for(Int_t ientry=0; ientry<=entries; ientry++){
 361
 362     nSigmaAB->GetEntry(ientry);
 363     //    cout << " pt="<< pt<<" sigma-AB="<<sigmaAB<<endl;
 364     if ( !(sigmaAB>0.) ) continue;
 365
 366     // Compute RAB and the statistical uncertainty
 367     Int_t hppbin = hSigmaPP->FindBin( pt );
 368     Double_t sigmapp = hSigmaPP->GetBinContent( hppbin );
 369     //    cout << " pt="<< pt<<", sigma-pp="<< sigmapp<<endl;
 370     if (isRaavsEP>0.) sigmapp = 0.5*sigmapp;
 371     if ( !(sigmapp>0.) ) continue;
 372     RaaCharm =  ( sigmaAB / sigmaABCINT1B ) / ((Tab*1e3) * sigmapp *1e-12 ) ;
 373     if (fdMethod==kNb) {
 374       RaaBeauty = Rb ;
 375     }
 376     else if (fdMethod==kfc) {
 377       RaaBeauty = ( RaaCharm / Rcb ) ;
 378     }
 379
 380     Double_t ElossHypo = 0.;
 381     if (fdMethod==kfc) { ElossHypo = 1. / Rcb; }
 382     else  { ElossHypo = 1. / (RaaCharm / RaaBeauty) ; }
 383     if(isRbHypo) ElossHypo = RaaBeauty;
 384     //    cout <<" pt "<< pt << " Raa charm " << RaaCharm << " Raa beauty " << RaaBeauty << " eloss hypo "<< ElossHypo<<endl;
 385     //
 386     // Find the bin for the central Eloss hypo
 387     //
 388     if( TMath::Abs( ElossHypo - CentralHypo ) < 0.075 ){
 389       Int_t hABbin = hSigmaAB->FindBin( pt );
 390       Double_t DeltaIni = TMath::Abs( ElossCentral[ hABbin ] - CentralHypo );
 391       Double_t DeltaV = TMath::Abs( ElossHypo - CentralHypo );
 392       //      cout << " pt " << pt << " ECentral " << ElossCentral[ hABbin ] << " Ehypo "<< ElossHypo ;
 393       if ( DeltaV < DeltaIni ) ElossCentral[ hABbin ] = ElossHypo;
 394       //      cout << " final ECentral " << ElossCentral[ hABbin ] << endl;
 395     }
 396   }
 397   //
 398   // Calculation of the Raa and its uncertainties
 399   //
 400   for(Int_t ientry=0; ientry<entries; ientry++){
 401
 402     nSigmaAB->GetEntry(ientry);
 403     if ( !(sigmaAB>0.) ) continue;
 404     //    if ( pt<2 || pt>16) continue;
 405
 406
 407     // Compute RAB and the statistical uncertainty
 408     Int_t hppbin = hSigmaPP->FindBin( pt );
 409     Double_t sigmapp = hSigmaPP->GetBinContent( hppbin );
 410     if (isRaavsEP>0.) sigmapp = 0.5*sigmapp;
 411     if ( !(sigmapp>0.) ) continue;
 412     RaaCharm =  ( sigmaAB / sigmaABCINT1B ) / ((Tab*1e3) * sigmapp *1e-12 );
 413
 414     // Flag to know if it is an scaled or extrapolated point of the pp reference
 415     Bool_t isExtrapolatedBin = kFALSE;
 416     if(isScaledAndExtrapRef) isExtrapolatedBin = hCombinedReferenceFlag->GetBinContent( hppbin );
 417     istartPPsyst = -1;
 418     istartPPsyst = FindGraphBin(gSigmaPPSyst,pt);
 419
 420     //
 421     // FONLL Feed-Down systematics
 422     //
 423     istartPPfd = -1;
 424     if(!isExtrapolatedBin) istartPPfd = FindGraphBin(gSigmaPPSystFeedDown,pt);
 425     istartABfd = -1;
 426     istartABfd = FindGraphBin(gSigmaABSystFeedDown,pt);
 427
 428     //      cout << " Starting bin for pp is "<< istartPPfd <<", for AA is "<<istartABfd << endl;
 429     if(isExtrapolatedBin){
 430       Int_t ibinfd = FindGraphBin(gReferenceFdSyst,pt);
 431       yPPh = gReferenceFdSyst->GetErrorYhigh(ibinfd);
 432       yPPl = gReferenceFdSyst->GetErrorYlow(ibinfd);
 433     } else {
 434       yPPh = gSigmaPPSystFeedDown->GetErrorYhigh(istartPPfd);
 435       yPPl = gSigmaPPSystFeedDown->GetErrorYlow(istartPPfd);
 436     }
 437     if (isRaavsEP>0.) {
 438       yPPh = yPPh*0.5;
 439       yPPl = yPPl*0.5;
 440     }
 441
 442     yABh = gSigmaABSystFeedDown->GetErrorYhigh(istartABfd);
 443     yABl = gSigmaABSystFeedDown->GetErrorYlow(istartABfd);
 444
 445
 446     RaaCharmFDhigh = ( sigmaABMax / sigmaABCINT1B ) / ((Tab*1e3) * (sigmapp+yPPh) *1e-12 ) ;
 447     RaaCharmFDlow =  ( sigmaABMin / sigmaABCINT1B ) / ((Tab*1e3) * (sigmapp-yPPl) *1e-12 ) ;
 448     if(printout) cout << endl<<" pt "<< pt << " Raa " << RaaCharm <<" high "<< RaaCharmFDhigh << " low "<< RaaCharmFDlow<<endl;
 449
 450
 451     if (fdMethod==kNb) {
 452       RaaBeauty = Rb ;
 453       RaaBeautyFDlow = Rb ;
 454       RaaBeautyFDhigh = Rb ;
 455       ntupleRAB->Fill( pt, Tab*1e3, sigmapp*1e-12, sigmaAB*1e-12, sigmaAB/sigmaABCINT1B,
 456                        sigmaABMax / sigmaABCINT1B, sigmaABMin / sigmaABCINT1B,
 457                        RaaCharm, RaaCharmFDhigh, RaaCharmFDlow, RaaBeauty, fcAB );
 458     }
 459     else if (fdMethod==kfc) {
 460       RaaBeauty = ( RaaCharm / Rcb ) ;
 461       RaaBeautyFDlow = ( RaaCharmFDlow / Rcb ) ;
 462       RaaBeautyFDhigh = ( RaaCharmFDhigh / Rcb ) ;
 463       hRABvsRcb->Fill( pt, RaaCharm, RaaBeauty );
 464       ntupleRAB->Fill( pt, Tab*1e3, sigmapp*1e-12, sigmaAB*1e-12, sigmaAB/sigmaABCINT1B,
 465                        sigmaABMax / sigmaABCINT1B, sigmaABMin / sigmaABCINT1B,
 466                        Rcb, RaaCharm, RaaCharmFDhigh, RaaCharmFDlow, RaaBeauty, RaaBeautyFDhigh, RaaBeautyFDlow, fcAB );
 467     }
 468     hRABvsRb->Fill( pt, RaaCharm, RaaBeauty );
 469     hRABvsRbFDlow->Fill( pt, RaaCharmFDlow, RaaBeautyFDlow );
 470     hRABvsRbFDhigh->Fill( pt, RaaCharmFDhigh, RaaBeautyFDhigh );
 471     if(printout) cout << " pt "<< pt << " Rb " << RaaBeauty <<" high "<< RaaBeautyFDhigh << " low "<< RaaBeautyFDlow <<endl;
 472
 473     hRABCharmVsRBeautyVsPt->Fill( pt, RaaBeauty, RaaCharm );
 474     Int_t ptbin = hRABvsPt->FindBin( pt );
 475     hRCharmVsRBeauty[ptbin]->Fill( RaaBeauty, RaaCharm );
 476     hRCharmVsRBeauty[ptbin]->Fill( RaaBeautyFDlow, RaaCharmFDlow );
 477     hRCharmVsRBeauty[ptbin]->Fill( RaaBeautyFDhigh, RaaCharmFDhigh );
 478
 479
 480     if (fdMethod==kfc) {
 481       if( TMath::Abs(Rcb-0.015)<0.009 ) hRABEloss00->Fill(pt,RaaCharm);
 482       if( TMath::Abs(Rcb-0.5)<0.009 ) hRABEloss05->Fill(pt,RaaCharm);
 483       if( TMath::Abs(Rcb-1.0)<0.009 ) {
 484         hRABEloss10->Fill(pt,RaaCharm);
 485         hRABvsRbFDhigh_proj->Fill(pt,RaaCharmFDhigh);
 486         hRABvsRbFDlow_proj->Fill(pt,RaaCharmFDlow);
 487       }
 488       if( TMath::Abs(Rcb-1.5)<0.009 ) hRABEloss15->Fill(pt,RaaCharm);
 489       if( TMath::Abs(Rcb-2.0)<0.009 ) hRABEloss20->Fill(pt,RaaCharm);
 490     }
 491     else if (fdMethod==kNb) {
 492       if( TMath::Abs(RaaBeauty-0.015)<0.009 ) hRABEloss00->Fill(pt,RaaCharm);
 493       if( TMath::Abs(RaaBeauty-0.5)<0.009 ) hRABEloss05->Fill(pt,RaaCharm);
 494       if( TMath::Abs(RaaBeauty-1.0)<0.009 ) {
 495         hRABEloss10->Fill(pt,RaaCharm);
 496         hRABvsRbFDhigh_proj->Fill(pt,RaaCharmFDhigh);
 497         hRABvsRbFDlow_proj->Fill(pt,RaaCharmFDlow);
 498       }
 499       if( TMath::Abs(RaaBeauty-1.5)<0.009 ) hRABEloss15->Fill(pt,RaaCharm);
 500       if( TMath::Abs(RaaBeauty-2.0)<0.009 ) hRABEloss20->Fill(pt,RaaCharm);
 501     }
 502
 503
 504     Int_t hABbin = hMassAB->FindBin( pt );
 505     if(printout)
 506     if ( fdMethod==kNb && TMath::Abs(Rb -1.0)< 0.05) {
 507       cout << " pt "<< pt <<", at bin "<<hABbin<<endl;
 508       cout<<" entries "<<entries<<", i="<<ientry<<", pt="<<pt<<", Rb="<<Rb<<", Tab="<<Tab<<", sigmaAB="<<sigmaAB<<", sigmapp="<<sigmapp<<", Raacharm="<<RaaCharm<<", RaaBeauty="<<RaaBeauty<<endl;
 509       cout << "  AB  basis: mass "<< hMassAB->GetBinContent(hABbin)<<", eff "<< hDirectEffptAB->GetBinContent(hABbin)<<endl;
 510       cout<<"   FD low,  err low AB "<< (sigmaAB-sigmaABMin)<<"  err low PP "<< yPPl<<" Raacharm="<<RaaCharmFDlow<<", RaaBeauty="<<RaaBeautyFDlow<<endl;
 511       cout<<"   FD high, err high AB "<< (sigmaABMax-sigmaAB)<<"  err high PP "<< yPPh<<" Raacharm="<<RaaCharmFDhigh<<", RaaBeauty="<<RaaBeautyFDhigh<<endl;
 512     }
 513     if(printout)
 514     if ( fdMethod==kfc) if(TMath::Abs(Rcb -1.0)< 0.05 ){
 515       cout << " pt "<< pt <<", at bin "<<hABbin<<endl;
 516       cout<<" entries "<<entries<<", i="<<ientry<<", pt="<<pt<<", Rcb="<<Rcb<<", Tab="<<Tab<<", sigmaAB="<<sigmaAB<<", sigmapp="<<sigmapp<<", Raacharm="<<RaaCharm<<", RaaBeauty="<<RaaBeauty<<endl;
 517       cout << "  AB  basis: mass "<< hMassAB->GetBinContent(hABbin)<<", eff "<< hDirectEffptAB->GetBinContent(hABbin)<<", fc "<<histofcAB->GetBinContent(hABbin)<< endl;
 518       cout<<"   FD low,  err low AB "<< (sigmaAB-sigmaABMin)<<"  err low PP "<< yPPl<<" Raacharm="<<RaaCharmFDlow<<", RaaBeauty="<<RaaBeautyFDlow<<endl;
 519       cout<<"   FD high, err high AB "<< (sigmaABMax-sigmaAB)<<"  err high PP "<< yPPh<<" Raacharm="<<RaaCharmFDhigh<<", RaaBeauty="<<RaaBeautyFDhigh<<endl;
 520     }
 521
 522
 523     //
 524     // Fill in the global properties ?
 525     //
 526     Double_t ElossHypo = 0.;
 527     if (fdMethod==kfc) { ElossHypo = 1./ Rcb; }
 528     else  { ElossHypo = 1. / (RaaCharm / RaaBeauty); }
 529     if(isRbHypo) ElossHypo = RaaBeauty;
 530     hRCharmVsElossHypo[ptbin]->Fill( ElossHypo, RaaCharm );
 531
 532     //    cout <<" pt "<< pt << " Raa charm " << RaaCharm << " Raa beauty " << RaaBeauty << " eloss hypo "<< ElossHypo<<endl;
 533     //
 534     // Fill in histos charm (null Eloss hypothesis)
 535     //
 536     Double_t minFdSyst = 0., maxFdSyst = 0.;
 537     if ( ElossHypo == ElossCentral[ hABbin ] ) {
 538
 539       //
 540       // Data stat uncertainty
 541       //
 542       Double_t sigmappStat = hSigmaPP->GetBinError( hppbin );
 543       if (isRaavsEP>0.) sigmappStat = sigmappStat*0.5;
 544       Int_t hRABbin = hRABvsPt->FindBin( pt );
 545       Double_t stat = RaaCharm * TMath::Sqrt( (statUncSigmaAB/sigmaAB)*(statUncSigmaAB/sigmaAB) +
 546                                               (sigmappStat/sigmapp)*(sigmappStat/sigmapp) ) ;
 547       if ( RaaCharm==0 ) stat =0.;
 548       if ( RaaCharm>0 ) {
 549         hRABvsPt->SetBinContent( hRABbin, RaaCharm );
 550         hRABvsPt->SetBinError( hRABbin, stat );
 551         hYieldABvsPt->SetBinContent( hRABbin, sigmaAB/sigmaABCINT1B );
 552         hYieldABvsPt->SetBinError( hRABbin, statUncSigmaAB/sigmaABCINT1B );
 553
 554         if(printout)
 555           cout << " Raa " << RaaCharm << " stat unc. "<< stat << " is "<< stat/RaaCharm * 100. <<
 556             "%, stat-pp "<< sigmappStat/sigmapp*100. <<"% stat-AB "<< statUncSigmaAB/sigmaAB*100.<<"%"<<endl;
 557
 558         Double_t errstatEff = fhStatUncEffcSigmaAB->GetBinError( hRABbin );
 559         fhStatUncEffcSigmaAB_Raa->SetBinError( hRABbin, errstatEff*RaaCharm );
 560         errstatEff = fhStatUncEffbSigmaAB->GetBinError( hRABbin );
 561         fhStatUncEffbSigmaAB_Raa->SetBinError( hRABbin, errstatEff*RaaCharm );
 562         errstatEff = fhStatUncEffcFDAB->GetBinError( hRABbin );
 563         fhStatUncEffcFDAB_Raa->SetBinError( hRABbin, errstatEff*RaaCharm );
 564         errstatEff = fhStatUncEffbFDAB->GetBinError( hRABbin );
 565         fhStatUncEffbFDAB_Raa->SetBinError( hRABbin, errstatEff*RaaCharm );
 566       }
 567
 568
 569       //
 570       //
 571       // Data systematics (sigma syst-but FD + extrap) syst
 572       //
 573       //
 574       // Data syst: a) Syst in p-p
 575       //
 576       Double_t ptwidth = hSigmaAB->GetBinWidth(hABbin) / 2. ;
 577       istartPPextr = -1;
 578       if(!isExtrapolatedBin) istartPPextr = FindGraphBin(gSigmaPPSystTheory,pt);
 579
 580       Double_t dataPPUp=0., dataPPLow=0.;
 581       if(isExtrapolatedBin) {
 582         dataPPUp = gSigmaPPSyst->GetErrorYhigh(istartPPsyst);
 583         dataPPLow = gSigmaPPSyst->GetErrorYlow(istartPPsyst);
 584         systPPUp = dataPPUp;
 585         systPPLow = dataPPLow;
 586       } else {
 587         dataPPUp = ExtractFDSyst( gSigmaPPSystData->GetErrorYhigh(istartPPextr), gSigmaPPSystFeedDown->GetErrorYhigh(istartPPfd) );
 588         dataPPLow = ExtractFDSyst( gSigmaPPSystData->GetErrorYlow(istartPPextr), gSigmaPPSystFeedDown->GetErrorYlow(istartPPfd) );
 589         systPPUp = TMath::Sqrt( dataPPUp*dataPPUp + gSigmaPPSystTheory->GetErrorYhigh(istartPPextr)*gSigmaPPSystTheory->GetErrorYhigh(istartPPextr) );
 590         systPPLow = TMath::Sqrt( dataPPLow*dataPPLow + gSigmaPPSystTheory->GetErrorYlow(istartPPextr)*gSigmaPPSystTheory->GetErrorYlow(istartPPextr) );
 591       }
 592       if (isRaavsEP>0.) {
 593         dataPPUp = dataPPUp*0.5;
 594         dataPPLow = dataPPLow*0.5;
 595         if(isExtrapolatedBin) {
 596           systPPUp = dataPPUp;
 597           systPPLow = dataPPLow;
 598         } else {
 599           systPPUp = TMath::Sqrt( dataPPUp*dataPPUp + 0.5*gSigmaPPSystTheory->GetErrorYhigh(istartPPextr)*0.5*gSigmaPPSystTheory->GetErrorYhigh(istartPPextr) );
 600           systPPLow = TMath::Sqrt( dataPPLow*dataPPLow + 0.5*gSigmaPPSystTheory->GetErrorYlow(istartPPextr)*0.5*gSigmaPPSystTheory->GetErrorYlow(istartPPextr) );
 601         }
 602       }
 603
 604
 605       if(printout) {
 606         cout << " pt : "<< pt<<" Syst-pp-data "<< dataPPUp/sigmapp << "%, ";
 607         if(!isExtrapolatedBin){
 608           if (isRaavsEP>0.) cout <<" extr unc + "<< 0.5*gSigmaPPSystTheory->GetErrorYhigh(istartPPextr)/sigmapp <<" - "<< 0.5*gSigmaPPSystTheory->GetErrorYlow(istartPPextr)/sigmapp <<" %";
 609           else cout <<" extr unc + "<< gSigmaPPSystTheory->GetErrorYhigh(istartPPextr)/sigmapp <<" - "<< gSigmaPPSystTheory->GetErrorYlow(istartPPextr)/sigmapp <<" %";
 610         }
 611         cout << endl;
 612       }
 613
 614       //
 615       // Data syst: b) Syst in PbPb
 616       //
 617       Double_t dataSystUp=0., dataSystDown=0.;
 618       Bool_t PbPbDataSystOk = PbPbDataSyst(systematicsAB,pt,cc,dataSystUp,dataSystDown);
 619       if (!PbPbDataSystOk) { cout <<" There is some issue with the PbPb data systematics, please check and rerun"<<endl; return; }
 620       systABUp = sigmaAB * TMath::Sqrt( dataSystUp*dataSystUp +
 621                                         (hDirectEffptAB->GetBinError(hABbin)/hDirectEffptAB->GetBinContent(hABbin))*(hDirectEffptAB->GetBinError(hABbin)/hDirectEffptAB->GetBinContent(hABbin)) );
 622
 623       systABLow = sigmaAB * TMath::Sqrt( dataSystDown*dataSystDown +
 624                                         (hDirectEffptAB->GetBinError(hABbin)/hDirectEffptAB->GetBinContent(hABbin))*(hDirectEffptAB->GetBinError(hABbin)/hDirectEffptAB->GetBinContent(hABbin)) );
 625       //
 626       // Data syst : c) combine pp & PbPb
 627       //
 628       systLow = sigmapp>0. ?
 629         RaaCharm * TMath::Sqrt( (systABLow/sigmaAB)*(systABLow/sigmaAB) + (systPPUp/sigmapp)*(systPPUp/sigmapp) )
 630         : 0.;
 631
 632       systUp = sigmapp>0. ?
 633         RaaCharm * TMath::Sqrt( (systABUp/sigmaAB)*(systABUp/sigmaAB) + (systPPLow/sigmapp)*(systPPLow/sigmapp) )
 634         : 0.;
 635       if ( RaaCharm==0 ) { systPPUp =0.; systPPLow =0.; }
 636
 637       //      if(printout)
 638         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;
 639
 640       if ( RaaCharm>0 ) {
 641         hRABvsPt_DataSystematics->SetBinContent( hRABbin, RaaCharm );
 642         hRABvsPt_DataSystematics->SetBinError( hRABbin, systUp );
 643         gRAB_DataSystematics->SetPoint( hABbin, pt, RaaCharm ); // i, x, y
 644         gRAB_DataSystematics->SetPointError( hABbin, ptwidth, ptwidth, systLow, systUp );
 645         gRAB_DataSystematics->SetPointEXlow(hABbin, 0.4); gRAB_DataSystematics->SetPointEXhigh(hABbin,0.4);
 646         gRAB_DataSystematicsPP->SetPoint( hABbin, pt, RaaCharm ); // i, x, y
 647         gRAB_DataSystematicsPP->SetPointError( hABbin, ptwidth, ptwidth, RaaCharm *(systPPUp/sigmapp), RaaCharm *systPPLow/sigmapp );
 648         gRAB_DataSystematicsAB->SetPoint( hABbin, pt, RaaCharm ); // i, x, y
 649         gRAB_DataSystematicsAB->SetPointError( hABbin, ptwidth, ptwidth, RaaCharm *systABLow/sigmaAB, RaaCharm *systABUp/sigmaAB );
 650       }
 651
 652       //
 653       // Feed-down Systematics
 654       //
 655       Double_t FDL=0., FDH=0.;
 656       if ( RaaCharmFDhigh > RaaCharmFDlow ){
 657         FDH = RaaCharmFDhigh; FDL = RaaCharmFDlow;
 658       } else {
 659         FDL = RaaCharmFDhigh; FDH = RaaCharmFDlow;
 660       }
 661
 662       if(printout) cout<<" Raa "<<RaaCharm<<", Raa-fd-low "<<RaaCharmFDlow <<", Raa-fd-high "<<RaaCharmFDhigh <<endl;
 663       maxFdSyst = TMath::Abs(FDH - RaaCharm);
 664       minFdSyst = TMath::Abs(RaaCharm - FDL);
 665       if ( RaaCharm>0 ) {
 666         gRAB_FeedDownSystematics->SetPoint( hABbin, pt, RaaCharm ); // i, x, y
 667         gRAB_FeedDownSystematics->SetPointError( hABbin, 0.3, 0.3, minFdSyst, maxFdSyst ); // i, x, y
 668         gRAB_fcFeedDownOnly->SetPoint( hABbin, pt,fcAB );
 669         gRAB_fcFeedDownOnly->SetPointError(hABbin, 0.3, 0.3, fcAB-(sigmaABMin/sigmaAB*fcAB), (sigmaABMax/sigmaAB*fcAB)-fcAB );
 670       }
 671
 672       //      if(printout) {
 673         cout<<" FD syst  +"<< maxFdSyst/RaaCharm <<" - "<<minFdSyst/RaaCharm<<endl;
 674         cout<<"  fc = "<<fcAB<<", ("<< sigmaABMax/sigmaAB * fcAB <<","<< sigmaABMin/sigmaAB * fcAB <<")"<<endl;
 675         //      }
 676
 677       //
 678       // Filling part of the Eloss scenarii information
 679       //
 680       if(RaaCharm>0 ) {
 681         gRAB_ElossHypothesis->SetPoint( hABbin, pt, RaaCharm ); // i, x, y
 682         gRAB_ElossHypothesis->SetPointEXlow( hABbin, ptwidth);
 683         gRAB_ElossHypothesis->SetPointEXhigh( hABbin, ptwidth);
 684         gRAB_FeedDownSystematicsElossHypothesis->SetPoint( hABbin, pt, RaaCharm ); // i, x, y
 685         gRAB_FeedDownSystematicsElossHypothesis->SetPointEXlow( hABbin, ptwidth);
 686         gRAB_FeedDownSystematicsElossHypothesis->SetPointEXhigh( hABbin, ptwidth);
 687         gRAB_GlobalSystematics->SetPoint( hABbin, pt, RaaCharm ); // i, x, y
 688         gRAB_GlobalSystematics->SetPointEXlow(hABbin,0.4); gRAB_GlobalSystematics->SetPointEXhigh(hABbin,0.4);
 689       }
 690     }
 691
 692     //
 693     // Filling Eloss scenarii information
 694     //
 695     if( RaaCharm>0 && ElossHypo >= MinHypo && ElossHypo <=MaxHypo && RaaBeauty<=MaxRb ) {
 696       Double_t Ehigh =  ElossMax[ hABbin ] ;
 697       Double_t Elow =  ElossMin[ hABbin ] ;
 698       if ( RaaCharm > Ehigh ) ElossMax[ hABbin ] = RaaCharm ;
 699       if ( RaaCharm < Elow ) ElossMin[ hABbin ] = RaaCharm ;
 700       if(printout) {
 701         cout<<" Hypothesis " << ElossHypo << " sigma-AB "<< sigmaAB <<", Raa "<< RaaCharm <<", Raa Eloss max "<< ElossMax[hABbin] <<" Raa Eloss min "<< ElossMin[hABbin] << " Rb="<< RaaBeauty <<endl;
 702         cout<<"  Rb="<< RaaBeauty <<" max "<< RaaBeautyFDhigh <<" min "<< RaaBeautyFDlow <<endl;
 703       }
 704       Double_t fcEhigh =  fcElossMax[ hABbin ] ;
 705       Double_t fcElow =  fcElossMin[ hABbin ] ;
 706       if ( fcAB > fcEhigh ) fcElossMax[ hABbin ] = fcAB ;
 707       if ( fcAB < fcElow ) fcElossMin[ hABbin ] = fcAB ;
 708       Double_t FDEhigh = FDElossMax[ hABbin ];
 709       Double_t FDEmin = FDElossMin[ hABbin ];
 710       Double_t RFDhigh = RaaCharmFDhigh>RaaCharmFDlow ? RaaCharmFDhigh : RaaCharmFDlow;
 711       Double_t RFDlow = RaaCharmFDlow<RaaCharmFDhigh ? RaaCharmFDlow : RaaCharmFDhigh;
 712       if ( RFDhigh > FDEhigh ) FDElossMax[ hABbin ] = RFDhigh ;
 713       if ( RFDlow < FDEmin ) FDElossMin[ hABbin ] = RFDlow ;
 714       if(printout)
 715         cout<<" Hypothesis " << ElossHypo << " sigma-AB "<< sigmaAB <<", Raa FD-max Eloss max "<< FDElossMax[hABbin] <<" Raa FD-min Eloss min "<< FDElossMin[hABbin] <<endl;
 716     }
 717
 718
 719   }
 720   delete [] limits;
 721   delete [] binwidths;
 722
 723
 724   // Finish filling the y-uncertainties of the Eloss scenarii
 725   for (Int_t ibin=0; ibin<=nbins; ibin++){
 726     Double_t ipt=0., value =0.;
 727     gRAB_ElossHypothesis->GetPoint(ibin,ipt,value);
 728     if(ipt<=0) continue;
 729     //
 730     // Uncertainty on Raa due to the Eloss hypothesis
 731     Double_t elossYhigh = TMath::Abs( ElossMax[ibin] - value );
 732     Double_t elossYlow = TMath::Abs( value - ElossMin[ibin] );
 733     gRAB_ElossHypothesis->SetPointEYhigh(ibin, elossYhigh );
 734     gRAB_ElossHypothesis->SetPointEYlow(ibin, elossYlow );
 735     gRAB_ElossHypothesis->SetPointEXhigh(ibin, 0.2);
 736     gRAB_ElossHypothesis->SetPointEXlow(ibin, 0.2);
 737     cout << " pt "<< ipt << " Raa "<< value <<" max "<< ElossMax[ibin] << " min " <<ElossMin[ibin] <<endl;
 738     cout<<" Eloss syst  +"<< elossYhigh <<" - "<< elossYlow <<endl;
 739     //    cout << " fc max "<< fcElossMax[ibin] << " fc min " <<fcElossMin[ibin] <<endl;
 740     //
 741     // Uncertainty on Raa due to the FD unc & Eloss hypothesis
 742     Double_t fdElossEYhigh = TMath::Abs( FDElossMax[ibin] - value );
 743     Double_t fdElossEYlow = TMath::Abs( value - FDElossMin[ibin] );
 744     if(elossFDQuadSum){
 745       Double_t fdEYhigh = gRAB_FeedDownSystematics->GetErrorYhigh(ibin);
 746       fdElossEYhigh = TMath::Sqrt( elossYhigh*elossYhigh + fdEYhigh*fdEYhigh );
 747       Double_t fdEYlow = gRAB_FeedDownSystematics->GetErrorYlow(ibin);
 748       fdElossEYlow = TMath::Sqrt( elossYlow*elossYlow + fdEYlow*fdEYlow );
 749     }
 750     gRAB_FeedDownSystematicsElossHypothesis->SetPointEYhigh(ibin, fdElossEYhigh );
 751     gRAB_FeedDownSystematicsElossHypothesis->SetPointEYlow(ibin, fdElossEYlow );
 752     gRAB_FeedDownSystematicsElossHypothesis->SetPointEXhigh(ibin, 0.25);
 753     gRAB_FeedDownSystematicsElossHypothesis->SetPointEXlow(ibin, 0.25);
 754     cout<<" FD & Eloss syst  +"<< fdElossEYhigh <<" - "<< fdElossEYlow
 755         <<" = + "<< fdElossEYhigh/value <<" - "<< fdElossEYlow/value <<" %" <<endl;
 756     //
 757     // All uncertainty on Raa (FD unc & Eloss + data)
 758     Double_t systdatal = gRAB_DataSystematics->GetErrorYlow(ibin);
 759     Double_t systdatah = gRAB_DataSystematics->GetErrorYhigh(ibin);
 760     Double_t systgbhUnc = TMath::Sqrt( systdatah*systdatah + fdElossEYhigh*fdElossEYhigh );
 761     Double_t systgblUnc = TMath::Sqrt( systdatal*systdatal + fdElossEYlow*fdElossEYlow );
 762     gRAB_GlobalSystematics->SetPointEYhigh(ibin,systgbhUnc);
 763     gRAB_GlobalSystematics->SetPointEYlow(ibin,systgblUnc);
 764     cout<<" Data syst  +"<< systdatah <<" - "<<  systdatal <<" = + "<< systdatah/value <<" - " <<  systdatal/value << " % "<<endl;
 765     cout<<" Global syst  +"<< systgbhUnc <<" - "<<  systgblUnc << " = + "<< systgbhUnc/value <<" - "<<  systgblUnc/value << " %" <<endl;
 766     //
 767   }
 768
 769
 770   gROOT->SetStyle("Plain");
 771   gStyle->SetPalette(1);
 772   gStyle->SetOptStat(0);
 773
 774
 775   TCanvas *cRABvsRb = new TCanvas("RABvsRb","RAB vs Rb");
 776   hRABvsRb->Draw("colz");
 777   cRABvsRb->Update();
 778
 779   TCanvas *cRABvsRbvsPt = new TCanvas("cRABvsRbvsPt","RAB vs Rb vs pt");
 780   hRABCharmVsRBeautyVsPt->Draw("lego3z");
 781   cRABvsRbvsPt->Update();
 782
 783
 784   TCanvas *cRABvsRbFDl = new TCanvas("RABvsRbFDl","RAB vs Rb (FD low)");
 785   hRABvsRbFDlow->Draw("cont4z");
 786   cRABvsRbFDl->Update();
 787   TCanvas *cRABvsRbFDh = new TCanvas("RABvsRbFDh","RAB vs Rb (FD high)");
 788   hRABvsRbFDhigh->Draw("cont4z");
 789   cRABvsRbFDh->Update();
 790
 791   TCanvas * cSigmaABptEloss = new TCanvas("cSigmaABptEloss","SigmaAB vs pt, Eloss hypothesis");
 792   TH1D *hSigmaABEloss00= new TH1D("hSigmaABEloss00","hSigmaABEloss00",nbins,limits);
 793   TH1D *hSigmaABEloss05= new TH1D("hSigmaABEloss05","hSigmaABEloss05",nbins,limits);
 794   TH1D *hSigmaABEloss10= new TH1D("hSigmaABEloss10","hSigmaABEloss10",nbins,limits);
 795   TH1D *hSigmaABEloss15= new TH1D("hSigmaABEloss15","hSigmaABEloss15",nbins,limits);
 796   TH1D *hSigmaABEloss20= new TH1D("hSigmaABEloss20","hSigmaABEloss20",nbins,limits);
 797   for (Int_t i=0; i<=nSigmaAB->GetEntriesFast(); i++) {
 798     nSigmaAB->GetEntry(i);
 799     if (fdMethod==kfc) {
 800       if( TMath::Abs(Rcb-0.015)<0.009 ) hSigmaABEloss00->Fill(pt,sigmaAB);
 801       if( TMath::Abs(Rcb-0.5)<0.009 ) hSigmaABEloss05->Fill(pt,sigmaAB);
 802       if( TMath::Abs(Rcb-1.0)<0.009 ) hSigmaABEloss10->Fill(pt,sigmaAB);
 803       if( TMath::Abs(Rcb-1.5)<0.009 ) hSigmaABEloss15->Fill(pt,sigmaAB);
 804       if( TMath::Abs(Rcb-2.0)<0.009 ) hSigmaABEloss20->Fill(pt,sigmaAB);
 805     }
 806     else if (fdMethod==kNb) {
 807       if( TMath::Abs(Rb-0.015)<0.009 ) hSigmaABEloss00->Fill(pt,sigmaAB);
 808       if( TMath::Abs(Rb-0.5)<0.009 ) hSigmaABEloss05->Fill(pt,sigmaAB);
 809       if( TMath::Abs(Rb-1.0)<0.009 ) hSigmaABEloss10->Fill(pt,sigmaAB);
 810       if( TMath::Abs(Rb-1.5)<0.009 ) hSigmaABEloss15->Fill(pt,sigmaAB);
 811       if( TMath::Abs(Rb-2.0)<0.009 ) hSigmaABEloss20->Fill(pt,sigmaAB);
 812     }
 813   }
 814   hSigmaABEloss00->SetLineColor(2);
 815   hSigmaABEloss05->SetLineColor(3);
 816   hSigmaABEloss10->SetLineColor(4);
 817   hSigmaABEloss15->SetLineColor(kMagenta+1);
 818   hSigmaABEloss20->SetLineColor(kGreen+2);
 819   hSigmaABEloss00->SetMarkerStyle(22);
 820   hSigmaABEloss05->SetMarkerStyle(26);
 821   hSigmaABEloss10->SetMarkerStyle(20);
 822   hSigmaABEloss15->SetMarkerStyle(25);
 823   hSigmaABEloss20->SetMarkerStyle(21);
 824   if (fdMethod==kNb) {
 825     hSigmaABEloss05->Draw("ph");
 826     hSigmaABEloss10->Draw("phsame");
 827     hSigmaABEloss15->Draw("phsame");
 828     hSigmaABEloss20->Draw("phsame");
 829   }
 830   else {
 831     hSigmaABEloss20->Draw("p");
 832     hSigmaABEloss00->Draw("phsame");
 833     hSigmaABEloss05->Draw("phsame");
 834     hSigmaABEloss10->Draw("phsame");
 835     hSigmaABEloss15->Draw("phsame");
 836     hSigmaABEloss20->Draw("phsame");
 837   }
 838   TLegend *legrcb = new TLegend(0.8,0.8,0.95,0.9);
 839   legrcb->SetFillColor(0);
 840   legrcb->AddEntry(hSigmaABEloss00,"Rc/b=0.0","lp");
 841   legrcb->AddEntry(hSigmaABEloss05,"Rc/b=0.5","lp");
 842   legrcb->AddEntry(hSigmaABEloss10,"Rc/b=1.0","lp");
 843   legrcb->AddEntry(hSigmaABEloss15,"Rc/b=1.5","lp");
 844   legrcb->AddEntry(hSigmaABEloss20,"Rc/b=2.0","lp");
 845   legrcb->Draw();
 846   cSigmaABptEloss->Update();
 847
 848
 849   TCanvas * cRABptEloss = new TCanvas("cRABptEloss","RAB vs pt, Eloss hypothesis");
 850   hRABEloss00->SetLineColor(2);
 851   hRABEloss05->SetLineColor(3);
 852   hRABEloss10->SetLineColor(4);
 853   hRABEloss15->SetLineColor(kMagenta+1);
 854   hRABEloss20->SetLineColor(kGreen+2);
 855   hRABEloss00->SetMarkerStyle(22);
 856   hRABEloss05->SetMarkerStyle(26);
 857   hRABEloss10->SetMarkerStyle(20);
 858   hRABEloss15->SetMarkerStyle(25);
 859   hRABEloss20->SetMarkerStyle(21);
 860   if (fdMethod==kNb) {
 861     hRABEloss05->Draw("ph");
 862     hRABEloss10->Draw("phsame");
 863     hRABEloss15->Draw("phsame");
 864     hRABEloss20->Draw("phsame");
 865   }
 866   else {
 867     hRABEloss20->Draw("p");
 868     hRABEloss00->Draw("phsame");
 869     hRABEloss05->Draw("phsame");
 870     hRABEloss10->Draw("phsame");
 871     hRABEloss15->Draw("phsame");
 872     hRABEloss20->Draw("phsame");
 873   }
 874   legrcb = new TLegend(0.8,0.8,0.95,0.9);
 875   legrcb->SetFillColor(0);
 876   if (fdMethod==kfc) {
 877     legrcb->AddEntry(hRABEloss00,"Rc/b=0.0","lp");
 878     legrcb->AddEntry(hRABEloss05,"Rc/b=0.5","lp");
 879     legrcb->AddEntry(hRABEloss10,"Rc/b=1.0","lp");
 880     legrcb->AddEntry(hRABEloss15,"Rc/b=0.5","lp");
 881     legrcb->AddEntry(hRABEloss20,"Rc/b=2.0","lp");
 882   }
 883   else if (fdMethod==kNb) {
 884     legrcb->AddEntry(hRABEloss00,"Rb=0.0","lp");
 885     legrcb->AddEntry(hRABEloss05,"Rb=0.5","lp");
 886     legrcb->AddEntry(hRABEloss10,"Rb=1.0","lp");
 887     legrcb->AddEntry(hRABEloss15,"Rb=0.5","lp");
 888     legrcb->AddEntry(hRABEloss20,"Rb=2.0","lp");
 889   }
 890   legrcb->Draw();
 891   cRABptEloss->Update();
 892
 893
 894   TCanvas * cRABpt = new TCanvas("cRABpt","RAB vs pt, no hypothesis");
 895   hRABEloss10->Draw("");
 896   cRABpt->Update();
 897
 898   TCanvas * cRABptFDUnc = new TCanvas("cRABptFDUnc","RAB vs pt, FD Uncertainties");
 899   hRABvsRbFDlow_proj->Draw("");
 900   hRABEloss10->Draw("phsame");
 901   hRABvsRbFDhigh_proj->SetLineColor(kMagenta+1);
 902   hRABvsRbFDhigh_proj->Draw("same");
 903   hRABvsRbFDlow_proj->SetLineColor(kGreen+2);
 904   hRABvsRbFDlow_proj->Draw("same");
 905   legrcb = new TLegend(0.8,0.8,0.95,0.9);
 906   legrcb->SetFillColor(0);
 907   legrcb->AddEntry(hRABEloss10,"FD Central","lp");
 908   legrcb->AddEntry(hRABvsRbFDhigh_proj,"FD Upper unc.","l");
 909   legrcb->AddEntry(hRABvsRbFDlow_proj,"FD Lower unc.","l");
 910   legrcb->Draw();
 911   cRABptFDUnc->Update();
 912
 913   TCanvas *RaaPlot = new TCanvas("RaaPlot","RAB vs pt, plot all");
 914   RaaPlot->SetTopMargin(0.085);
 915   RaaPlot->SetBottomMargin(0.1);
 916   RaaPlot->SetTickx();
 917   RaaPlot->SetTicky();
 918   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.,2.0);
 919   hRaaCanvas->GetXaxis()->SetTitleSize(0.05);
 920   hRaaCanvas->GetXaxis()->SetTitleOffset(0.9);
 921   hRaaCanvas->GetYaxis()->SetTitleSize(0.05);
 922   hRaaCanvas->GetYaxis()->SetTitleOffset(0.9);
 923   hRaaCanvas->Draw();
 924   gRAB_Norm->SetFillStyle(1001);
 925   gRAB_Norm->SetFillColor(kGray+2);
 926   gRAB_Norm->Draw("2");
 927   TLine *line = new TLine(0.0172415,1.0,40.,1.0);
 928   line->SetLineStyle(2);
 929   line->Draw();
 930   hRABvsPt->SetMarkerColor(kBlue);
 931   hRABvsPt->SetMarkerColor(kBlue);
 932   hRABvsPt->SetMarkerStyle(21);
 933   hRABvsPt->SetMarkerSize(1.1);
 934   hRABvsPt->SetLineWidth(2);
 935   hRABvsPt->Draw("psame");
 936   gRAB_DataSystematics->SetLineColor(kBlue);
 937   gRAB_DataSystematics->SetLineWidth(3);
 938   gRAB_DataSystematics->SetLineWidth(2);
 939   gRAB_DataSystematics->SetFillColor(kRed);
 940   gRAB_DataSystematics->SetFillStyle(0);
 941   gRAB_DataSystematics->Draw("2");
 942   gRAB_FeedDownSystematics->SetFillColor(kViolet+1);
 943   gRAB_FeedDownSystematics->SetFillStyle(1001);
 944   gRAB_FeedDownSystematics->Draw("2");
 945   gRAB_ElossHypothesis->SetLineColor(kMagenta-7);
 946   gRAB_ElossHypothesis->SetFillColor(kMagenta-7);
 947   gRAB_ElossHypothesis->SetFillStyle(1001);
 948   gRAB_ElossHypothesis->Draw("2");
 949   hRABvsPt->Draw("psame");
 950   gRAB_DataSystematics->Draw("2");
 951   legrcb = new TLegend(0.5517241,0.6504237,0.8520115,0.8728814,NULL,"brNDC");
 952   legrcb->SetBorderSize(0);
 953   legrcb->SetTextSize(0.03389831);
 954   legrcb->SetLineColor(1);
 955   legrcb->SetLineStyle(1);
 956   legrcb->SetLineWidth(1);
 957   legrcb->SetFillColor(0);
 958   legrcb->SetFillStyle(1001);
 959   if(cc==k020) legrcb->AddEntry(hRABvsPt,"R_{AA} 0-20% CC","pe");
 960   else if(cc==k4080) legrcb->AddEntry(hRABvsPt,"R_{AA} 40-80% CC","pe");
 961   else legrcb->AddEntry(hRABvsPt,"R_{AA} and stat. unc.","pe");
 962   legrcb->AddEntry(gRAB_DataSystematics,"Syst. from data","f");
 963   legrcb->AddEntry(gRAB_ElossHypothesis,"Syst. from R_{AA}(B)","f");
 964   legrcb->AddEntry(gRAB_FeedDownSystematics,"Syst. from B feed-down","f");
 965   legrcb->Draw();
 966   TLatex* tc;
 967   if(decay==1) tc =new TLatex(0.18,0.82,"D^{0},  Pb-Pb    #sqrt{s_{NN}}=2.76 TeV");
 968   else if(decay==2) tc =new TLatex(0.18,0.82,"D^{+},  Pb-Pb    #sqrt{s_{NN}}=2.76 TeV");
 969   else if(decay==3) tc =new TLatex(0.18,0.82,"D^{*+},  Pb-Pb    #sqrt{s_{NN}}=2.76 TeV");
 970   else if(decay==4) tc =new TLatex(0.18,0.82,"D_{s}^{+},  Pb-Pb    #sqrt{s_{NN}}=2.76 TeV");
 971   else  tc =new TLatex(0.18,0.82,"any (?) D meson,  Pb-Pb    #sqrt{s_{NN}}=2.76 TeV");
 972   tc->SetNDC();
 973   tc->SetTextSize(0.038);
 974   tc->SetTextFont(42);
 975   tc->Draw();
 976   RaaPlot->Update();
 977
 978
 979   TCanvas *RaaPlotFDEloss = new TCanvas("RaaPlotFDEloss","RAB vs pt, plot FD & ElossUnc");
 980   RaaPlotFDEloss->SetTopMargin(0.085);
 981   RaaPlotFDEloss->SetBottomMargin(0.1);
 982   hRaaCanvas->Draw();
 983   line->Draw();
 984   hRABvsPt->Draw("psame");
 985   gRAB_FeedDownSystematics->SetFillColor(kViolet+1);
 986   gRAB_FeedDownSystematics->SetFillStyle(1001);
 987   gRAB_FeedDownSystematics->Draw("2");
 988   gRAB_ElossHypothesis->SetLineColor(kMagenta-7);
 989   gRAB_ElossHypothesis->SetFillColor(kMagenta-7);
 990   gRAB_ElossHypothesis->SetFillStyle(1001);
 991   gRAB_ElossHypothesis->Draw("2");
 992   gRAB_FeedDownSystematicsElossHypothesis->SetLineColor(kBlack);
 993   gRAB_FeedDownSystematicsElossHypothesis->SetFillStyle(0);
 994   gRAB_FeedDownSystematicsElossHypothesis->SetFillColor(kViolet+1);
 995   gRAB_FeedDownSystematicsElossHypothesis->Draw("2");
 996   hRABvsPt->Draw("psame");
 997   legrcb = new TLegend(0.6,0.6,0.9,0.9);
 998   legrcb->SetBorderSize(0);
 999   legrcb->SetTextSize(0.03389831);
1000   legrcb->SetLineColor(1);
1001   legrcb->SetLineStyle(1);
1002   legrcb->SetLineWidth(1);
1003   legrcb->SetFillColor(0);
1004   legrcb->SetFillStyle(1001);
1005   legrcb->AddEntry(hRABvsPt,"R_{PbPb} and stat. unc.","pe");
1006   legrcb->AddEntry(gRAB_ElossHypothesis,"Energy loss syst.","f");
1007   legrcb->AddEntry(gRAB_FeedDownSystematics,"Feed down syst.","f");
1008   legrcb->AddEntry(gRAB_FeedDownSystematicsElossHypothesis,"Feed down & Eloss syst.","f");
1009   legrcb->Draw();
1010   RaaPlotFDEloss->Update();
1011
1012
1013   TCanvas *RaaPlotGlob = new TCanvas("RaaPlotGlob","RAB vs pt, plot Global unc");
1014   RaaPlotGlob->SetTopMargin(0.085);
1015   RaaPlotGlob->SetBottomMargin(0.1);
1016   RaaPlotGlob->SetTickx();
1017   RaaPlotGlob->SetTicky();
1018   hRaaCanvas->Draw();
1019   line->Draw();
1020   hRABvsPt->Draw("psame");
1021   gRAB_DataSystematics->Draw("2");
1022   gRAB_FeedDownSystematicsElossHypothesis->Draw("2");
1023   gRAB_GlobalSystematics->SetLineColor(kRed);
1024   gRAB_GlobalSystematics->SetLineWidth(2);
1025   gRAB_GlobalSystematics->SetFillColor(kRed);
1026   gRAB_GlobalSystematics->SetFillStyle(3002);
1027   gRAB_GlobalSystematics->Draw("2");
1028   hRABvsPt->Draw("psame");
1029   legrcb = new TLegend(0.6,0.6,0.9,0.9);
1030   legrcb->SetBorderSize(0);
1031   legrcb->SetTextSize(0.03389831);
1032   legrcb->SetLineColor(1);
1033   legrcb->SetLineStyle(1);
1034   legrcb->SetLineWidth(1);
1035   legrcb->SetFillColor(0);
1036   legrcb->SetFillStyle(1001);
1037   legrcb->AddEntry(hRABvsPt,"R_{PbPb} and stat. unc.","pe");
1038   legrcb->AddEntry(gRAB_DataSystematics,"Data syst.","f");
1039   legrcb->AddEntry(gRAB_FeedDownSystematicsElossHypothesis,"Feed down & Eloss syst.","f");
1040   legrcb->AddEntry(gRAB_GlobalSystematics,"Global syst.","f");
1041   legrcb->Draw();
1042   RaaPlotGlob->Update();
1043
1044
1045
1046   TCanvas *RaaPlotSimple = new TCanvas("RaaPlotSimple","RAB vs pt, plot Simple unc");
1047   RaaPlotSimple->SetTopMargin(0.085);
1048   RaaPlotSimple->SetBottomMargin(0.1);
1049   RaaPlotSimple->SetTickx();
1050   RaaPlotSimple->SetTicky();
1051   hRaaCanvas->Draw();
1052   line->Draw();
1053   hRABvsPt->Draw("psame");
1054   gRAB_GlobalSystematics->SetLineColor(kBlue);
1055   gRAB_GlobalSystematics->SetLineWidth(2);
1056   gRAB_GlobalSystematics->SetFillStyle(0);
1057   gRAB_GlobalSystematics->Draw("2");
1058   gRAB_Norm->Draw("2");
1059   hRABvsPt->Draw("psame");
1060   legrcb = new TLegend(0.5991379,0.6949153,0.8534483,0.8559322,NULL,"brNDC");
1061   legrcb->SetBorderSize(0);
1062   legrcb->SetTextSize(0.03389831);
1063   legrcb->SetLineColor(1);
1064   legrcb->SetLineStyle(1);
1065   legrcb->SetLineWidth(1);
1066   legrcb->SetFillColor(0);
1067   legrcb->SetFillStyle(1001);
1068   if(cc==k020) legrcb->AddEntry(hRABvsPt,"R_{AA} 0-20% CC","pe");
1069   else if(cc==k4080) legrcb->AddEntry(hRABvsPt,"R_{AA} 40-80% CC","pe");
1070   else legrcb->AddEntry(hRABvsPt,"R_{AA} and stat. unc.","pe");
1071   legrcb->AddEntry(gRAB_GlobalSystematics,"Systematics","f");
1072   legrcb->Draw();
1073   tc->Draw();
1074   RaaPlotSimple->Update();
1075
1076
1077   TCanvas *c = new TCanvas("c","");
1078   systematicsAB->DrawErrors();
1079   c->Update();
1080
1081   TCanvas *cStatUnc = new TCanvas("cStatUnc","stat unc");
1082   cStatUnc->Divide(2,2);
1083   cStatUnc->cd(1);
1084   fhStatUncEffcSigmaAB_Raa->Draw("e");
1085   cStatUnc->cd(2);
1086   fhStatUncEffbSigmaAB_Raa->Draw("e");
1087   cStatUnc->cd(3);
1088   fhStatUncEffcFDAB_Raa->Draw("e");
1089   cStatUnc->cd(4);
1090   fhStatUncEffbFDAB_Raa->Draw("e");
1091   cStatUnc->Update();
1092
1093   //
1094   // Write the information to a file
1095   //
1096   TFile * out = new TFile(outfile,"recreate");
1097
1098   ntupleRAB->Write();
1099   hRABvsRcb->Write();
1100   hRABvsRb->Write();
1101   hRABCharmVsRBeautyVsPt->Write();
1102   for(Int_t j=0; j<=nbins; j++) hRCharmVsRBeauty[j]->Write();
1103   //  for(Int_t j=0; j<=nbins; j++) hRCharmVsElossHypo[j]->Write();
1104   hRABvsPt->Write();
1105   hRABvsPt_DataSystematics->Write();
1106   gRAB_ElossHypothesis->Write();
1107   gRAB_FeedDownSystematics->Write();
1108   gRAB_fcFeedDownOnly->Write();
1109   gRAB_DataSystematics->Write();
1110   gRAB_DataSystematicsPP->Write();
1111   gSigmaPPSystTheory->Write();
1112   gRAB_DataSystematicsAB->Write();
1113   gRAB_Norm->Write();
1114   gRAB_FeedDownSystematicsElossHypothesis->Write();
1115   gRAB_GlobalSystematics->Write();
1116   if(isScaledAndExtrapRef) hCombinedReferenceFlag->Write();
1117
1118   out->Write();
1119
1120 }
1121
1122 //____________________________________________________________
1123 Bool_t PbPbDataSyst(AliHFSystErr *syst, Double_t pt, Int_t cc, Double_t &dataSystUp, Double_t &dataSystDown)
1124 {
1125
1126   Double_t err=0., errUp=1., errDown=1.;
1127   Bool_t isOk=false;
1128   Double_t pidunc=0.;
1129
1130   err = syst->GetTotalSystErr(pt)*syst->GetTotalSystErr(pt);
1131   errDown = err ;
1132
1133   if( cc==k07half || cc==k020 || cc==k010 || cc==k1020 || cc==k2040 ) {
1134     if(pt<6) pidunc = 0.15;
1135     else pidunc = 0.05;
1136     errUp = err + pidunc*pidunc - syst->GetPIDEffErr(pt)*syst->GetPIDEffErr(pt);
1137     isOk = true;
1138   }
1139   else if ( cc==k3050 || cc==k4080 || cc==k4060 || cc==k6080 ){
1140     if(pt<3.1) pidunc = 0.10;
1141     else pidunc = 0.05;
1142     errUp = err + pidunc*pidunc - syst->GetPIDEffErr(pt)*syst->GetPIDEffErr(pt);
1143     isOk = true;
1144   }
1145
1146   dataSystUp = TMath::Sqrt(errUp);
1147   dataSystDown = TMath::Sqrt(errDown);
1148
1149   /// << ------------------------- PATCH FOR 2011 LHC11h data ------------------------------
1150   dataSystUp = dataSystDown;
1151
1152   return isOk;
1153 }