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