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