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