]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGHF/vertexingHF/macros/HFPtSpectrumRaa.C
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / macros / HFPtSpectrumRaa.C
index 6faa746058fcbe72fbb590339fe38cab09154fa2..7cb960cc208379353d7ad22d26f837a3636c0fd9 100644 (file)
 //      9. MinHypo  = minimum energy loss hypothesis (Default 1./3.)
 //      10. MaxHypo = maximum energy loss hypothesis (Default 3.0)
 //      11. MaxRb = maximum Raa(b) hypothesis (Default 6.0, won't do anything)
-//      11. kRaavsEP = flag to compute the Raa IN/OUT of plane, divides the reference by 2.0
+//      12. RbHypo : flag to decide whether the Eloss hypothesis is Rb or Rb/Rc
+//      13. CentralHypo = central energy loss hypothesis, DEFAULT TO 1.0
+//      14. isRaavsEP = flag to compute the Raa IN/OUT of plane, divides the reference by 2.0
+//      15. ScaledAndExtrapRef: flag to tag scaled+reference pp-scaled-data
 //
 //  Complains to : Zaida Conesa del Valle
 //
 ///////////////////////////////////////////////////////////////////////////////
 
-enum centrality{ kpp, k07half, k010, k1020, k020, k2040, k3050, k4060, k6080, k4080, k80100 };
-enum energy{ k276, k55 };
+enum centrality{ kpp, k07half, kpPb0100, k010, k1020, k020, k2040, k2030, k3040, k4050, k3050, k5060, k4060, k6080, k4080, k5080, k80100, kpPb020, kpPb2040, kpPb4060, kpPb60100 };
+enum centestimator{ kV0M, kV0A, kZNA };
+enum energy{ k276, k5dot023, k55 };
 enum BFDSubtrMethod { kfc, kNb };
+enum RaavsEP {kPhiIntegrated, kInPlane, kOutOfPlane};
 
 Bool_t printout = false;
+Double_t ptprintout = 1.5;
 Double_t NormPPUnc = 0.035;
+Double_t NormABUnc = 0.037;
 Bool_t elossFDQuadSum = true;
 
 //____________________________________________________________
@@ -65,6 +72,23 @@ Double_t ExtractFDSyst(Double_t total, Double_t fd) {
   return TMath::Sqrt( data2 );
 }
 
+//____________________________________________________________
+Int_t FindGraphBin(TGraphAsymmErrors *gr, Double_t pt)
+{
+  Int_t istart =0;
+  Int_t npoints = gr->GetN();
+  for(Int_t i=0; i<=npoints; i++){
+    Double_t x=0.,y=0.;
+    gr->GetPoint(i,x,y);
+    if ( TMath::Abs ( x - pt ) < 0.4 ) { 
+      istart = i; 
+      break;
+    }
+  }
+  return istart;
+}
+
+
 //
 //
 // R_AB =  [ ( dsigma/dpt )_AB / sigma_AB ]  / <TAB> *  ( dsigma/dpt )_pp
@@ -77,7 +101,11 @@ void HFPtSpectrumRaa(const char *ppfile="HFPtSpectrum_D0Kpi_method2_rebinnedth_2
                     Int_t decay=1,
                     Double_t sigmaABCINT1B=54.e9,
                     Int_t fdMethod = kNb, Int_t cc=kpp, Int_t Energy=k276,
-                    Double_t MinHypo=1./3., Double_t MaxHypo=3.0, Double_t MaxRb=6.0, Bool_t kRaavsEP=kFALSE)
+                    Double_t MinHypo=1./3., Double_t MaxHypo=3.0, Double_t MaxRb=6.0,
+                    Bool_t isRbHypo=false, Double_t CentralHypo = 1.0,
+                    Int_t ccestimator = kV0M,
+                    Bool_t isUseTaaForRaa=true, const char *shadRbcFile="", Int_t nSigmaShad=3.0,
+                    Int_t isRaavsEP=kPhiIntegrated, Bool_t isScaledAndExtrapRef=kFALSE)
 {
 
   gROOT->Macro("$ALICE_ROOT/PWGHF/vertexingHF/macros/LoadLibraries.C");
@@ -85,8 +113,8 @@ void HFPtSpectrumRaa(const char *ppfile="HFPtSpectrum_D0Kpi_method2_rebinnedth_2
   //
   // Defining the TAB values for the given centrality class
   //
-  Double_t Tab = 1., TabSyst = 0.;
-  if ( Energy!=k276 ) {
+  Double_t Tab = 1., TabSyst = 0., A=207.2, B=207.2;
+  if ( Energy!=k276 && Energy!=k5dot023) {
     printf("\n The Tab values for this cms energy have not yet been implemented, please do it ! \n");
     return;
   }
@@ -95,38 +123,89 @@ void HFPtSpectrumRaa(const char *ppfile="HFPtSpectrum_D0Kpi_method2_rebinnedth_2
   }
   // Values from Alberica's twiki:
   //   https://twiki.cern.ch/twiki/bin/viewauth/ALICE/CentStudies
-  if ( cc == k07half ) {
-    Tab = 24.81; TabSyst = 0.8037;
-  } else if ( cc == k010 ) {
-    Tab = 23.48; TabSyst = 0.97;
-  } else if ( cc == k1020 ) {
-    Tab = 14.4318; TabSyst = 0.573289;
-  } else if ( cc == k020 ) {
-    Tab = 18.93; TabSyst = 0.74;
-  } else if ( cc == k2040 ) {
-    Tab = 6.86; TabSyst = 0.28;
-  } else if ( cc == k3050 ) {
-    Tab = 3.87011; TabSyst = 0.183847;
-  } else if ( cc == k4060 ) {
-    Tab = 2.00; TabSyst = 0.11;
-  } else if ( cc == k6080 ) {
-    Tab = 0.419; TabSyst = 0.033;
-  } else if ( cc == k4080 ) {
-    Tab = 1.20451; TabSyst = 0.071843;
-  } else if ( cc == k80100 ){
-    Tab = 0.0690; TabSyst = 0.0062;
+  if( ccestimator == kV0M ) {
+    if ( cc == k07half ) {
+      Tab = 24.81; TabSyst = 0.8037;
+    } else if ( cc == k010 ) {
+      Tab = 23.48; TabSyst = 0.97;
+    } else if ( cc == k1020 ) {
+      Tab = 14.4318; TabSyst = 0.5733;
+    } else if ( cc == k020 ) {
+      Tab = 18.93; TabSyst = 0.74;
+    } else if ( cc == k2040 ) {
+      Tab = 6.86; TabSyst = 0.28;
+    } else if ( cc == k2030 ) {
+      Tab = 8.73769; TabSyst = 0.370219;
+    } else if ( cc == k3040 ) {
+      Tab = 5.02755; TabSyst = 0.22099;
+    } else if ( cc == k4050 ) {
+      Tab = 2.68327; TabSyst = 0.137073;
+    } else if ( cc == k3050 ) {
+      Tab = 3.87011; TabSyst = 0.183847;
+    } else if ( cc == k4060 ) {
+      Tab = 2.00;  TabSyst= 0.11;
+    } else if ( cc == k4080 ) {
+      Tab = 1.20451; TabSyst = 0.071843;
+    } else if ( cc == k5060 ) {
+      Tab = 1.32884; TabSyst = 0.0929536;
+    } else if ( cc == k6080 ) {
+      Tab = 0.419; TabSyst = 0.033;
+    } else if ( cc == k5080 ) {
+      Tab = 0.719; TabSyst = 0.054;
+    } else if ( cc == k80100 ){
+      Tab = 0.0690; TabSyst = 0.0062;
+    }
   }
 
+  // pPb Glauber (A. Toia)
+  // https://twiki.cern.ch/twiki/bin/viewauth/ALICE/PACentStudies#Glauber_Calculations_with_sigma
+  if( cc == kpPb0100 ){
+    Tab = 0.098334; TabSyst = 0.0070679;
+    A=207.2; B=1.;
+  }
+  else if( ccestimator == kV0A ){
+    if ( cc == kpPb020 ) {
+      Tab = 0.183; TabSyst = 0.006245;
+    } else if ( cc == kpPb2040 ) {
+      Tab = 0.134; TabSyst = 0.004899;
+    } else if ( cc == kpPb4060 ) {
+      Tab = 0.092; TabSyst = 0.004796;
+    } else if ( cc == kpPb60100 ) {
+      Tab = 0.041; TabSyst = 0.008832;
+    }
+  }
+  else if( ccestimator == kZNA ){
+    if ( cc == kpPb020 ) {
+      Tab = 0.164; TabSyst = 0.010724;
+    } else if ( cc == kpPb2040 ) {
+      Tab = 0.137; TabSyst = 0.005099;
+    } else if ( cc == kpPb4060 ) {
+      Tab = 0.1011; TabSyst = 0.006;
+    } else if ( cc == kpPb60100 ) {
+      Tab = 0.0459; TabSyst = 0.003162;
+    }
+  }
 
   //
   // Reading the pp file 
   //
   TFile * ppf = new TFile(ppfile,"read");
-  TH1D * hSigmaPP = (TH1D*)ppf->Get("fhScaledData");
-  TGraphAsymmErrors * gSigmaPPSyst = (TGraphAsymmErrors*)ppf->Get("gScaledData");
+  TH1D * hSigmaPP;
+  TGraphAsymmErrors * gSigmaPPSyst;
   TGraphAsymmErrors * gSigmaPPSystData = (TGraphAsymmErrors*)ppf->Get("gScaledDataSystData");
   TGraphAsymmErrors * gSigmaPPSystTheory = (TGraphAsymmErrors*)ppf->Get("gScaledDataSystExtrap");
   TGraphAsymmErrors * gSigmaPPSystFeedDown = (TGraphAsymmErrors*)ppf->Get("gScaledDataSystFeedDown");
+  TH1I * hCombinedReferenceFlag;
+  TGraphAsymmErrors * gReferenceFdSyst;
+  if(isScaledAndExtrapRef){
+    hCombinedReferenceFlag = (TH1I*)ppf->Get("hCombinedReferenceFlag");
+    hSigmaPP = (TH1D*)ppf->Get("hReference");
+    gSigmaPPSyst = (TGraphAsymmErrors*)ppf->Get("gReferenceSyst");
+    gReferenceFdSyst = (TGraphAsymmErrors*)ppf->Get("gReferenceFdSyst");
+  } else { 
+    hSigmaPP = (TH1D*)ppf->Get("fhScaledData");
+    gSigmaPPSyst = (TGraphAsymmErrors*)ppf->Get("gScaledData");
+  }
   
   // Call the systematics uncertainty class for a given decay
   AliHFSystErr *systematicsPP = new AliHFSystErr();
@@ -137,8 +216,8 @@ void HFPtSpectrumRaa(const char *ppfile="HFPtSpectrum_D0Kpi_method2_rebinnedth_2
   // 
   TFile * ABf = new TFile(ABfile,"read");
   TH1D *hSigmaAB = (TH1D*)ABf->Get("histoSigmaCorr");
-  TH2D *hSigmaABRcb = (TH2D*)ABf->Get("histoSigmaCorrRcb");
-  TGraphAsymmErrors * gSigmaABSyst = (TGraphAsymmErrors*)ABf->Get("gSigmaCorr");
+  //  TH2D *hSigmaABRcb = (TH2D*)ABf->Get("histoSigmaCorrRcb");
+  //  TGraphAsymmErrors * gSigmaABSyst = (TGraphAsymmErrors*)ABf->Get("gSigmaCorr");
   TGraphAsymmErrors * gSigmaABSystFeedDown = (TGraphAsymmErrors*)ABf->Get("gSigmaCorrConservative");
   TNtuple * nSigmaAB = (TNtuple*)ABf->Get("fnSigma");
   //
@@ -164,28 +243,54 @@ void HFPtSpectrumRaa(const char *ppfile="HFPtSpectrum_D0Kpi_method2_rebinnedth_2
   fhStatUncEffcFDAB_Raa->SetName("fhStatUncEffcFDAB_Raa");
   fhStatUncEffbFDAB_Raa->SetName("fhStatUncEffvFDAB_Raa");
 
+
   //
   // Call the systematics uncertainty class for a given decay
   AliHFSystErr *systematicsAB = new AliHFSystErr();
   systematicsAB->SetCollisionType(1);
-  if ( cc == k07half ) systematicsAB->SetCentrality("010");
+  if ( cc == k07half ) systematicsAB->SetCentrality("07half");
   else if ( cc == k010 ) systematicsAB->SetCentrality("010");
   else if ( cc == k1020 ) systematicsAB->SetCentrality("1020");
   else if ( cc == k020 ) systematicsAB->SetCentrality("020");
-  else if ( cc == k2040 ) {
+  else if ( cc == k2040 || cc == k2030 || cc == k3040 ) {
     systematicsAB->SetCentrality("2040");
     systematicsAB->SetIsPbPb2010EnergyScan(true);
   }
-  else if ( cc == k4060 ) systematicsAB->SetCentrality("4060");
-  else if ( cc == k6080 ) systematicsAB->SetCentrality("6080");
+  else if ( cc == k4060 || cc == k4050 || cc == k5060 ) systematicsAB->SetCentrality("4060");
+  else if ( cc == k6080 || cc == k5080 ) systematicsAB->SetCentrality("6080");
   else if ( cc == k4080 ) systematicsAB->SetCentrality("4080");
-  else if ( cc == k3050 ) systematicsAB->SetCentrality("4080");
+  else if ( cc == k3050 ) {
+    if (isRaavsEP == kPhiIntegrated) systematicsAB->SetCentrality("4080");
+    else if (isRaavsEP == kInPlane) systematicsAB->SetCentrality("3050InPlane");
+    else if (isRaavsEP == kOutOfPlane) systematicsAB->SetCentrality("3050OutOfPlane");
+  }
+  //
+  else if ( cc == kpPb0100 || cc == kpPb020 || cc == kpPb2040 || cc == kpPb4060 || cc == kpPb60100 ) {
+    systematicsAB->SetCollisionType(2);
+    if(ccestimator==kV0A) {
+      if(cc == kpPb020) systematicsAB->SetCentrality("020V0A");
+      else if(cc == kpPb2040) systematicsAB->SetCentrality("2040V0A");
+      else if(cc == kpPb4060) systematicsAB->SetCentrality("4060V0A");
+      else if(cc == kpPb60100) systematicsAB->SetCentrality("60100V0A");
+    } else if (ccestimator==kZNA) {
+      if(cc == kpPb020) systematicsAB->SetCentrality("020ZNA");
+      else if(cc == kpPb2040) systematicsAB->SetCentrality("2040ZNA");
+      else if(cc == kpPb4060) systematicsAB->SetCentrality("4060ZNA");
+      else if(cc == kpPb60100) systematicsAB->SetCentrality("60100ZNA");
+    } else {
+      if(!(cc == kpPb0100)) {
+       cout <<" Error on the pPb options"<<endl;
+       return;
+      }
+    }
+  }
   else { 
     cout << " Systematics not yet implemented " << endl;
     return;
-  } 
-  if(decay!=4) systematicsAB->Init(decay);
-  else systematicsAB->Init(2);
+  }
+  //
+  systematicsAB->Init(decay);
+
   //
   Int_t entries = nSigmaAB->GetEntries();
   Float_t pt=0., signal=0., Rb=0., Rcb=0., fcAB=0., yieldAB=0., sigmaAB=0., statUncSigmaAB=0., sigmaABMin=0.,sigmaABMax=0.;
@@ -214,20 +319,63 @@ void HFPtSpectrumRaa(const char *ppfile="HFPtSpectrum_D0Kpi_method2_rebinnedth_2
     binwidths[i-1] = binwidth;
   }
   limits[nbins] = xlow + binwidth;
+
+
+  //
+  // Read the shadowing file if given as input
+  //
+  Double_t centralRbcShad[nbins+1], minRbcShad[nbins+1], maxRbcShad[nbins+1];
+  for(Int_t i=0; i<=nbins; i++) { centralRbcShad[i]=1.0; minRbcShad[i]=6.0; maxRbcShad[i]=0.0; }
+  Bool_t isShadHypothesis = false;
+  if( strcmp(shadRbcFile,"")!=0 ) {
+    isShadHypothesis = true;
+    cout<<endl<<">>  Beware, using the shadowing prediction file with an "<<nSigmaShad<<"*sigma <<"<<endl<<endl;
+    TFile *fshad = new TFile(shadRbcFile,"read");
+    if(!fshad){ cout <<" >> Shadowing file not properly opened!!!"<<endl<<endl; return;}
+    // TH1D *hRbcShadCentral = (TH1D*)fshad->Get("hDfromBoverPromptD_Shadowing_central");
+    // TH1D *hRbcShadMin = (TH1D*)fshad->Get("hDfromBoverPromptD_Shadowing_upper");
+    // TH1D *hRbcShadMax = (TH1D*)fshad->Get("hDfromBoverPromptD_Shadowing_lower");
+    TH1D *hRbcShadCentral = (TH1D*)fshad->Get("hDfromBoverDfromc_L0");
+    TH1D *hRbcShadMin = (TH1D*)fshad->Get("hDfromBoverDfromc_L0");
+    TH1D *hRbcShadMax = (TH1D*)fshad->Get("hDfromBoverDfromc_L1");
+    if(!hRbcShadCentral || !hRbcShadMin || !hRbcShadMax) {
+      cout<< endl <<">> Shadowing input histograms are not ok !! "<<endl<<endl;
+      return;
+    }
+    //       nSigmaShad
+    //    nSigmaShad
+    for(Int_t i=1; i<=nbins; i++) {
+      Double_t xpt = hSigmaAB->GetBinCenter(i);
+      if(xpt>24) xpt = 20;
+      centralRbcShad[i] = hRbcShadCentral->GetBinContent( hRbcShadCentral->FindBin(xpt) );
+      Double_t minValue0 = hRbcShadMin->GetBinContent( hRbcShadMin->FindBin(xpt) );
+      Double_t maxValue0 = hRbcShadMax->GetBinContent( hRbcShadMax->FindBin(xpt) );
+      Double_t arrayEl[3] = {minValue0,maxValue0, centralRbcShad[i]};
+      Double_t minValue = TMath::MinElement(3,arrayEl);
+      Double_t maxValue = TMath::MaxElement(3,arrayEl);
+      cout<<">> Shadowing pt="<<xpt<<"  central="<<centralRbcShad[i]<<"  min="<<minValue<<"  max="<<maxValue<<endl;
+      if(minValue>centralRbcShad[i]){ minValue = centralRbcShad[i]; }
+      if(maxValue<centralRbcShad[i]){ maxValue = centralRbcShad[i]; }
+      minRbcShad[i] = centralRbcShad[i] - nSigmaShad*(centralRbcShad[i] - minValue);
+      maxRbcShad[i] = centralRbcShad[i] + nSigmaShad*(maxValue - centralRbcShad[i]);
+      cout<<">> Shadowing hypothesis pt="<<xpt<<"  central="<<centralRbcShad[i]<<"  min="<<minRbcShad[i]<<"  max="<<maxRbcShad[i]<<endl;
+    }
+  }
+
   //
   // define the bins correspondence bw histos/files/graphs
   //
   //
   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.);
   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.);
-  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.);
+  //  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.);
   Int_t nbinsHypo=800;//200;
   Double_t *limitsHypo = new Double_t[nbinsHypo+1];
   for(Int_t i=1; i<=nbinsHypo+1; i++) limitsHypo[i-1]= i*4./800.;
   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);
-  TH2D *hRCharmVsRBeauty[nbins];
+  TH2D *hRCharmVsRBeauty[nbins+1];
   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);
-  TH2D *hRCharmVsElossHypo[nbins];
+  TH2D *hRCharmVsElossHypo[nbins+1];
   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);
   //
   TH1D *hRABEloss00= new TH1D("hRABEloss00","hRABEloss00",nbins,limits);
@@ -244,9 +392,9 @@ void HFPtSpectrumRaa(const char *ppfile="HFPtSpectrum_D0Kpi_method2_rebinnedth_2
   //
   TNtuple *ntupleRAB=0x0 ;
   if (fdMethod==kNb) {
-    ntupleRAB = new TNtuple("ntupleRAB","ntupleRAB (Nb)","pt:TAB:sigmaPP:sigmaAB:invyieldAB:invyieldABFDHigh:invyieldABFDLow:RABCharm:RABCharmFDHigh:RABCharmFDLow:RABBeauty");
+    ntupleRAB = new TNtuple("ntupleRAB","ntupleRAB (Nb)","pt:TAB:sigmaPP:sigmaAB:invyieldAB:invyieldABFDHigh:invyieldABFDLow:RABCharm:RABCharmFDHigh:RABCharmFDLow:RABBeauty:fc",100000);
   } else if (fdMethod==kfc) {
-    ntupleRAB = new TNtuple("ntupleRAB","ntupleRAB (fc)","pt:TAB:sigmaPP:sigmaAB:invyieldAB:invyieldABFDHigh:invyieldABFDLow:Rcb:RABCharm:RABCharmFDHigh:RABCharmFDLow:RABBeauty:RABBeautyFDHigh:RABBeautyFDLow");
+    ntupleRAB = new TNtuple("ntupleRAB","ntupleRAB (fc)","pt:TAB:sigmaPP:sigmaAB:invyieldAB:invyieldABFDHigh:invyieldABFDLow:Rcb:RABCharm:RABCharmFDHigh:RABCharmFDLow:RABBeauty:RABBeautyFDHigh:RABBeautyFDLow:fc",100000);
   }
   if(!ntupleRAB) printf("ERROR: Wrong method option");
 
@@ -269,16 +417,17 @@ void HFPtSpectrumRaa(const char *ppfile="HFPtSpectrum_D0Kpi_method2_rebinnedth_2
   gRAB_DataSystematicsAB->SetNameTitle("gRAB_DataSystematicsAB","RAB Measurement AB (no FD, no Eloss, no PP data) systematics");
   TGraphAsymmErrors *gRAB_GlobalSystematics = new TGraphAsymmErrors(nbins+1);
   gRAB_GlobalSystematics->SetNameTitle("gRAB_GlobalSystematics","RAB Measurement global (data, FD, Eloss) systematics");
-  Double_t ElossMax[nbins], ElossMin[nbins];
-  for(Int_t i=0; i<=nbins; i++) { ElossMax[i]=0.; ElossMin[i]=1.; }
-  Double_t fcElossMax[nbins], fcElossMin[nbins];
-  for(Int_t i=0; i<=nbins; i++) { fcElossMax[i]=0.; fcElossMin[i]=1.; }
-  Double_t FDElossMax[nbins], FDElossMin[nbins];
-  for(Int_t i=0; i<=nbins; i++) { FDElossMax[i]=0.; FDElossMin[i]=1.; }
+  Double_t ElossMax[nbins+1], ElossMin[nbins+1];
+  for(Int_t i=0; i<=nbins; i++) { ElossMax[i]=0.; ElossMin[i]=6.; }
+  Double_t fcElossMax[nbins+1], fcElossMin[nbins+1];
+  for(Int_t i=0; i<=nbins; i++) { fcElossMax[i]=0.; fcElossMin[i]=6.; }
+  Double_t FDElossMax[nbins+1], FDElossMin[nbins+1];
+  for(Int_t i=0; i<=nbins; i++) { FDElossMax[i]=0.; FDElossMin[i]=6.; }
 
   TGraphAsymmErrors *gRAB_Norm = new TGraphAsymmErrors(1);
   gRAB_Norm->SetNameTitle("gRAB_Norm","RAB Normalization systematics (pp norm + Tab)");
   Double_t normUnc = TMath::Sqrt ( NormPPUnc*NormPPUnc + (TabSyst/Tab)*(TabSyst/Tab) );
+  if(!isUseTaaForRaa) normUnc = TMath::Sqrt ( NormPPUnc*NormPPUnc + NormABUnc*NormABUnc );
   gRAB_Norm->SetPoint(1,0.5,1.);
   gRAB_Norm->SetPointError(1,0.25,0.25,normUnc,normUnc);
 
@@ -286,8 +435,7 @@ void HFPtSpectrumRaa(const char *ppfile="HFPtSpectrum_D0Kpi_method2_rebinnedth_2
   // R_AB =  ( dN/dpt )_AB  / <Ncoll_AB> *  ( dN/dpt )_pp ; <Ncoll> = <Tab> * sigma_NN^inel
   // R_AB =  [ ( dsigma/dpt )_AB / sigma_AB ]  / <TAB> *  ( dsigma/dpt )_pp
   //
-  Int_t istartPPfd=0, istartABfd=0, istartPPextr=0;
-  Double_t xPP=0., yPP=0., xAB=0., yAB=0.;
+  Int_t istartPPfd=0, istartPPsyst=0, istartABfd=0, istartPPextr=0;
   Double_t yPPh=0., yPPl=0., yABh=0., yABl=0.;
   Double_t RaaCharm =0., RaaBeauty=0.;
   Double_t RaaCharmFDhigh = 0., RaaCharmFDlow = 0.;
@@ -297,20 +445,29 @@ void HFPtSpectrumRaa(const char *ppfile="HFPtSpectrum_D0Kpi_method2_rebinnedth_2
   //
   // Search the central value of the energy loss hypothesis Rb = Rc (bin)
   //
-  Double_t ElossCentral[nbins];
+  Double_t ElossCentral[nbins+1];
   for(Int_t i=0; i<=nbins; i++) { ElossCentral[i]=0.; }
   //
   for(Int_t ientry=0; ientry<=entries; ientry++){
 
     nSigmaAB->GetEntry(ientry);
+    //    cout << " pt="<< pt<<" sigma-AB="<<sigmaAB<<endl;
     if ( !(sigmaAB>0.) ) continue;
+    //if(decay==2 && pt<2.) continue;
 
     // Compute RAB and the statistical uncertainty
     Int_t hppbin = hSigmaPP->FindBin( pt );
+    Int_t hABbin = hSigmaAB->FindBin( pt );
     Double_t sigmapp = hSigmaPP->GetBinContent( hppbin );
-    if (kRaavsEP) sigmapp = 0.5*sigmapp;
+    //    cout << " pt="<< pt<<", sigma-pp="<< sigmapp<<endl;
+    if (isRaavsEP>0.) sigmapp = 0.5*sigmapp;
     if ( !(sigmapp>0.) ) continue;
+
     RaaCharm =  ( sigmaAB / sigmaABCINT1B ) / ((Tab*1e3) * sigmapp *1e-12 ) ;
+    if(!isUseTaaForRaa) { 
+      RaaCharm =  ( sigmaAB ) / ( (A*B) * sigmapp ) ;
+    }
+
     if (fdMethod==kNb) {
       RaaBeauty = Rb ; 
     }
@@ -321,20 +478,22 @@ void HFPtSpectrumRaa(const char *ppfile="HFPtSpectrum_D0Kpi_method2_rebinnedth_2
     Double_t ElossHypo = 0.;
     if (fdMethod==kfc) { ElossHypo = 1. / Rcb; }
     else  { ElossHypo = 1. / (RaaCharm / RaaBeauty) ; }
+    if(isRbHypo) ElossHypo = RaaBeauty;
+
+    // If using shadowing hypothesis, change the central hypothesis too
+    if(isShadHypothesis) CentralHypo = centralRbcShad[hABbin];
 
     //    cout <<" pt "<< pt << " Raa charm " << RaaCharm << " Raa beauty " << RaaBeauty << " eloss hypo "<< ElossHypo<<endl; 
     //
     // Find the bin for the central Eloss hypo
     //
-    if( TMath::Abs( ElossHypo - 1.0 ) < 0.075 ){
-      Int_t hABbin = hSigmaAB->FindBin( pt );
-      Double_t DeltaIni = TMath::Abs( ElossCentral[ hABbin ] - 1.0 );
-      Double_t DeltaV = TMath::Abs( ElossHypo - 1.0 );
-      cout << " pt " << pt << " ECentral " << ElossCentral[ hABbin ] << " Ehypo "<< ElossHypo ;
+    if( TMath::Abs( ElossHypo - CentralHypo ) < 0.075 ){
+      Double_t DeltaIni = TMath::Abs( ElossCentral[ hABbin ] - CentralHypo );
+      Double_t DeltaV = TMath::Abs( ElossHypo - CentralHypo );
+      //      cout << " pt " << pt << " ECentral " << ElossCentral[ hABbin ] << " Ehypo "<< ElossHypo ;
       if ( DeltaV < DeltaIni ) ElossCentral[ hABbin ] = ElossHypo;
-      cout << " final ECentral " << ElossCentral[ hABbin ] << endl;
+      //      cout << " final ECentral " << ElossCentral[ hABbin ] << endl;
     }
-
   }
   //
   // Calculation of the Raa and its uncertainties
@@ -345,44 +504,57 @@ void HFPtSpectrumRaa(const char *ppfile="HFPtSpectrum_D0Kpi_method2_rebinnedth_2
     if ( !(sigmaAB>0.) ) continue;
     //    if ( pt<2 || pt>16) continue;
 
+
     // Compute RAB and the statistical uncertainty
     Int_t hppbin = hSigmaPP->FindBin( pt );
     Double_t sigmapp = hSigmaPP->GetBinContent( hppbin );
-    if (kRaavsEP) sigmapp = 0.5*sigmapp;
+    if (isRaavsEP>0.) sigmapp = 0.5*sigmapp;
     if ( !(sigmapp>0.) ) continue;
+
     RaaCharm =  ( sigmaAB / sigmaABCINT1B ) / ((Tab*1e3) * sigmapp *1e-12 );
+    if(!isUseTaaForRaa) {
+      RaaCharm =  ( sigmaAB ) / ( (A*B) * sigmapp ) ;
+    }
+
+    // Flag to know if it is an scaled or extrapolated point of the pp reference
+    Bool_t isExtrapolatedBin = kFALSE;
+    if(isScaledAndExtrapRef) isExtrapolatedBin = hCombinedReferenceFlag->GetBinContent( hppbin );
+    istartPPsyst = -1;
+    istartPPsyst = FindGraphBin(gSigmaPPSyst,pt);
 
     //
     // FONLL Feed-Down systematics
     //
-    Int_t n = gSigmaPPSystFeedDown->GetN();
-    for(Int_t j=1; j<=n; j++){
-      gSigmaPPSystFeedDown->GetPoint(j,xPP,yPP);
-      if ( TMath::Abs ( xPP -pt ) < 0.4 ) { 
-       istartPPfd = j; 
-       break;
-      }
+    istartPPfd = -1;
+    if(!isExtrapolatedBin) istartPPfd = FindGraphBin(gSigmaPPSystFeedDown,pt);
+    istartABfd = -1;
+    istartABfd = FindGraphBin(gSigmaABSystFeedDown,pt);
+
+    //      cout << " Starting bin for pp is "<< istartPPfd <<", for AA is "<<istartABfd << endl;
+    if(isExtrapolatedBin){
+      Int_t ibinfd = FindGraphBin(gReferenceFdSyst,pt);
+      yPPh = gReferenceFdSyst->GetErrorYhigh(ibinfd);
+      yPPl = gReferenceFdSyst->GetErrorYlow(ibinfd);
+    } else { 
+      yPPh = gSigmaPPSystFeedDown->GetErrorYhigh(istartPPfd);
+      yPPl = gSigmaPPSystFeedDown->GetErrorYlow(istartPPfd);
     }
-    n = gSigmaABSystFeedDown->GetN();
-    for(Int_t j=1; j<=n; j++){
-      gSigmaABSystFeedDown->GetPoint(j,xAB,yAB);
-      if ( TMath::Abs ( xAB -pt ) < 0.4 ) { 
-       istartABfd = j; 
-       break;
-      }
+    if (isRaavsEP>0.) {
+      yPPh = yPPh*0.5;
+      yPPl = yPPl*0.5;
     }
-    //      cout << " Starting bin for pp is "<< istartPPfd <<", for AA is "<<istartABfd << endl;
-    yPPh = gSigmaPPSystFeedDown->GetErrorYhigh(istartPPfd);
-    if (kRaavsEP) yPPh = yPPh*0.5;
-    yPPl = gSigmaPPSystFeedDown->GetErrorYlow(istartPPfd);
-    if (kRaavsEP) yPPl = yPPl*0.5;
+
     yABh = gSigmaABSystFeedDown->GetErrorYhigh(istartABfd);
     yABl = gSigmaABSystFeedDown->GetErrorYlow(istartABfd);
-    
+
 
     RaaCharmFDhigh = ( sigmaABMax / sigmaABCINT1B ) / ((Tab*1e3) * (sigmapp+yPPh) *1e-12 ) ;
     RaaCharmFDlow =  ( sigmaABMin / sigmaABCINT1B ) / ((Tab*1e3) * (sigmapp-yPPl) *1e-12 ) ;
-    if(printout) cout << endl<<" pt "<< pt << " Raa " << RaaCharm <<" high "<< RaaCharmFDhigh << " low "<< RaaCharmFDlow<<endl;
+    if(printout && TMath::Abs(ptprintout-pt)<0.1 ) cout << endl<<" pt "<< pt << " Raa " << RaaCharm <<" high "<< RaaCharmFDhigh << " low "<< RaaCharmFDlow<<endl;
+    if(!isUseTaaForRaa) {
+      RaaCharmFDhigh = ( sigmaABMax ) / ( (A*B)* (sigmapp+yPPh) ) ;
+      RaaCharmFDlow =  ( sigmaABMin ) / ( (A*B)* (sigmapp-yPPl) ) ;
+    }
 
 
     if (fdMethod==kNb) {
@@ -391,7 +563,7 @@ void HFPtSpectrumRaa(const char *ppfile="HFPtSpectrum_D0Kpi_method2_rebinnedth_2
       RaaBeautyFDhigh = Rb ;
       ntupleRAB->Fill( pt, Tab*1e3, sigmapp*1e-12, sigmaAB*1e-12, sigmaAB/sigmaABCINT1B,
                       sigmaABMax / sigmaABCINT1B, sigmaABMin / sigmaABCINT1B,
-                      RaaCharm, RaaCharmFDhigh, RaaCharmFDlow, RaaBeauty);
+                      RaaCharm, RaaCharmFDhigh, RaaCharmFDlow, RaaBeauty, fcAB );
     }
     else if (fdMethod==kfc) {
       RaaBeauty = ( RaaCharm / Rcb ) ;
@@ -400,12 +572,12 @@ void HFPtSpectrumRaa(const char *ppfile="HFPtSpectrum_D0Kpi_method2_rebinnedth_2
       hRABvsRcb->Fill( pt, RaaCharm, RaaBeauty );
       ntupleRAB->Fill( pt, Tab*1e3, sigmapp*1e-12, sigmaAB*1e-12, sigmaAB/sigmaABCINT1B,
                       sigmaABMax / sigmaABCINT1B, sigmaABMin / sigmaABCINT1B,
-                      Rcb, RaaCharm, RaaCharmFDhigh, RaaCharmFDlow, RaaBeauty, RaaBeautyFDhigh, RaaBeautyFDlow);
+                      Rcb, RaaCharm, RaaCharmFDhigh, RaaCharmFDlow, RaaBeauty, RaaBeautyFDhigh, RaaBeautyFDlow, fcAB );
     }
     hRABvsRb->Fill( pt, RaaCharm, RaaBeauty );
     hRABvsRbFDlow->Fill( pt, RaaCharmFDlow, RaaBeautyFDlow );
     hRABvsRbFDhigh->Fill( pt, RaaCharmFDhigh, RaaBeautyFDhigh );
-    if(printout) cout << " pt "<< pt << " Rb " << RaaBeauty <<" high "<< RaaBeautyFDhigh << " low "<< RaaBeautyFDlow <<endl;
+    if(printout && TMath::Abs(ptprintout-pt)<0.1) cout << " pt "<< pt << " Rb " << RaaBeauty <<" high "<< RaaBeautyFDhigh << " low "<< RaaBeautyFDlow <<endl;
 
     hRABCharmVsRBeautyVsPt->Fill( pt, RaaBeauty, RaaCharm );
     Int_t ptbin = hRABvsPt->FindBin( pt );
@@ -439,16 +611,18 @@ void HFPtSpectrumRaa(const char *ppfile="HFPtSpectrum_D0Kpi_method2_rebinnedth_2
 
 
     Int_t hABbin = hMassAB->FindBin( pt );
-    if(printout)
-    if ( fdMethod==kNb && TMath::Abs(Rb -1.0)< 0.05) {
+    if(isShadHypothesis) CentralHypo = centralRbcShad[hABbin];
+
+    if(printout && TMath::Abs(ptprintout-pt)<0.1)
+    if ( fdMethod==kNb && TMath::Abs(Rb -CentralHypo)< 0.05) {
       cout << " pt "<< pt <<", at bin "<<hABbin<<endl;
       cout<<" entries "<<entries<<", i="<<ientry<<", pt="<<pt<<", Rb="<<Rb<<", Tab="<<Tab<<", sigmaAB="<<sigmaAB<<", sigmapp="<<sigmapp<<", Raacharm="<<RaaCharm<<", RaaBeauty="<<RaaBeauty<<endl;
       cout << "  AB  basis: mass "<< hMassAB->GetBinContent(hABbin)<<", eff "<< hDirectEffptAB->GetBinContent(hABbin)<<endl;
       cout<<"   FD low,  err low AB "<< (sigmaAB-sigmaABMin)<<"  err low PP "<< yPPl<<" Raacharm="<<RaaCharmFDlow<<", RaaBeauty="<<RaaBeautyFDlow<<endl;
       cout<<"   FD high, err high AB "<< (sigmaABMax-sigmaAB)<<"  err high PP "<< yPPh<<" Raacharm="<<RaaCharmFDhigh<<", RaaBeauty="<<RaaBeautyFDhigh<<endl;
     }
-    if(printout)
-    if ( fdMethod==kfc) if(TMath::Abs(Rcb -1.0)< 0.05 ){
+    if(printout && TMath::Abs(ptprintout-pt)<0.1)
+    if ( fdMethod==kfc) if(TMath::Abs(Rcb -CentralHypo)< 0.05 ){
       cout << " pt "<< pt <<", at bin "<<hABbin<<endl;
       cout<<" entries "<<entries<<", i="<<ientry<<", pt="<<pt<<", Rcb="<<Rcb<<", Tab="<<Tab<<", sigmaAB="<<sigmaAB<<", sigmapp="<<sigmapp<<", Raacharm="<<RaaCharm<<", RaaBeauty="<<RaaBeauty<<endl;
       cout << "  AB  basis: mass "<< hMassAB->GetBinContent(hABbin)<<", eff "<< hDirectEffptAB->GetBinContent(hABbin)<<", fc "<<histofcAB->GetBinContent(hABbin)<< endl;       
@@ -463,9 +637,18 @@ void HFPtSpectrumRaa(const char *ppfile="HFPtSpectrum_D0Kpi_method2_rebinnedth_2
     Double_t ElossHypo = 0.;
     if (fdMethod==kfc) { ElossHypo = 1./ Rcb; }
     else  { ElossHypo = 1. / (RaaCharm / RaaBeauty); }
+    if(isRbHypo) ElossHypo = RaaBeauty;
     hRCharmVsElossHypo[ptbin]->Fill( ElossHypo, RaaCharm );
 
-    //    cout <<" pt "<< pt << " Raa charm " << RaaCharm << " Raa beauty " << RaaBeauty << " eloss hypo "<< ElossHypo<<endl; 
+   // If using shadowing hypothesis, change the limit hypothesis too
+    if(isShadHypothesis) {
+      MinHypo = minRbcShad[ hABbin ];
+      MaxHypo = maxRbcShad[ hABbin ];
+    }
+
+    //    cout <<" pt "<< pt << " Raa charm " << RaaCharm << " Raa beauty " << RaaBeauty << " eloss hypo "<< ElossHypo
+    if(ientry==0) cout<<" pt"<< pt<< " ElossCentral "<< ElossCentral[hABbin] << " min-hypo "<<MinHypo << " max-hypo "<<MaxHypo<<endl;
+
     //
     // Fill in histos charm (null Eloss hypothesis)
     //
@@ -476,7 +659,7 @@ void HFPtSpectrumRaa(const char *ppfile="HFPtSpectrum_D0Kpi_method2_rebinnedth_2
       // Data stat uncertainty
       //
       Double_t sigmappStat = hSigmaPP->GetBinError( hppbin );
-      if (kRaavsEP) sigmappStat = sigmappStat*0.5;
+      if (isRaavsEP>0.) sigmappStat = sigmappStat*0.5;
       Int_t hRABbin = hRABvsPt->FindBin( pt );
       Double_t stat = RaaCharm * TMath::Sqrt( (statUncSigmaAB/sigmaAB)*(statUncSigmaAB/sigmaAB) + 
                                              (sigmappStat/sigmapp)*(sigmappStat/sigmapp) ) ;
@@ -487,9 +670,12 @@ void HFPtSpectrumRaa(const char *ppfile="HFPtSpectrum_D0Kpi_method2_rebinnedth_2
        hYieldABvsPt->SetBinContent( hRABbin, sigmaAB/sigmaABCINT1B );
        hYieldABvsPt->SetBinError( hRABbin, statUncSigmaAB/sigmaABCINT1B );
        
-       if(printout) 
+       cout << "pt="<< pt<< " Raa " << RaaCharm << " stat unc. "<<
+         " sigma-pp "<< sigmapp <<" sigma-AB "<< sigmaAB<<endl;
+       if(printout && TMath::Abs(ptprintout-pt)<0.1) {
          cout << " Raa " << RaaCharm << " stat unc. "<< stat << " is "<< stat/RaaCharm * 100. <<
            "%, stat-pp "<< sigmappStat/sigmapp*100. <<"% stat-AB "<< statUncSigmaAB/sigmaAB*100.<<"%"<<endl;
+       }
 
        Double_t errstatEff = fhStatUncEffcSigmaAB->GetBinError( hRABbin );
        fhStatUncEffcSigmaAB_Raa->SetBinError( hRABbin, errstatEff*RaaCharm );
@@ -510,29 +696,41 @@ void HFPtSpectrumRaa(const char *ppfile="HFPtSpectrum_D0Kpi_method2_rebinnedth_2
       // Data syst: a) Syst in p-p 
       //
       Double_t ptwidth = hSigmaAB->GetBinWidth(hABbin) / 2. ;
-      n = gSigmaPPSystTheory->GetN();
-      for(Int_t j=1; j<=n; j++){
-       gSigmaPPSystTheory->GetPoint(j,xPP,yPP);
-       if ( TMath::Abs ( xPP -pt ) < 0.4 ) { 
-         istartPPextr = j; 
-         break;
+      istartPPextr = -1;
+      if(!isExtrapolatedBin) istartPPextr = FindGraphBin(gSigmaPPSystTheory,pt);
+
+      Double_t dataPPUp=0., dataPPLow=0.;
+      if(isExtrapolatedBin) {
+       dataPPUp = gSigmaPPSyst->GetErrorYhigh(istartPPsyst);
+       dataPPLow = gSigmaPPSyst->GetErrorYlow(istartPPsyst);
+       systPPUp = dataPPUp;
+       systPPLow = dataPPLow;
+      } else { 
+       dataPPUp = ExtractFDSyst( gSigmaPPSystData->GetErrorYhigh(istartPPextr), gSigmaPPSystFeedDown->GetErrorYhigh(istartPPfd) );
+       dataPPLow = ExtractFDSyst( gSigmaPPSystData->GetErrorYlow(istartPPextr), gSigmaPPSystFeedDown->GetErrorYlow(istartPPfd) );
+       systPPUp = TMath::Sqrt( dataPPUp*dataPPUp + gSigmaPPSystTheory->GetErrorYhigh(istartPPextr)*gSigmaPPSystTheory->GetErrorYhigh(istartPPextr) );
+       systPPLow = TMath::Sqrt( dataPPLow*dataPPLow + gSigmaPPSystTheory->GetErrorYlow(istartPPextr)*gSigmaPPSystTheory->GetErrorYlow(istartPPextr) );
+      }
+      if (isRaavsEP>0.) {
+       dataPPUp = dataPPUp*0.5;
+       dataPPLow = dataPPLow*0.5;
+       if(isExtrapolatedBin) {
+         systPPUp = dataPPUp;
+         systPPLow = dataPPLow;
+       } else {  
+         systPPUp = TMath::Sqrt( dataPPUp*dataPPUp + 0.5*gSigmaPPSystTheory->GetErrorYhigh(istartPPextr)*0.5*gSigmaPPSystTheory->GetErrorYhigh(istartPPextr) );
+         systPPLow = TMath::Sqrt( dataPPLow*dataPPLow + 0.5*gSigmaPPSystTheory->GetErrorYlow(istartPPextr)*0.5*gSigmaPPSystTheory->GetErrorYlow(istartPPextr) );
        }
       }
 
-      Double_t dataPPUp = ExtractFDSyst( gSigmaPPSystData->GetErrorYhigh(istartPPextr), gSigmaPPSystFeedDown->GetErrorYhigh(istartPPfd) );
-      if (kRaavsEP) dataPPUp = ExtractFDSyst( 0.5*gSigmaPPSystData->GetErrorYhigh(istartPPextr), 0.5*gSigmaPPSystFeedDown->GetErrorYhigh(istartPPfd) );
-      systPPUp = TMath::Sqrt( dataPPUp*dataPPUp + gSigmaPPSystTheory->GetErrorYhigh(istartPPextr)*gSigmaPPSystTheory->GetErrorYhigh(istartPPextr) );
-      if (kRaavsEP) systPPUp = TMath::Sqrt( dataPPUp*dataPPUp + 0.5*gSigmaPPSystTheory->GetErrorYhigh(istartPPextr)*0.5*gSigmaPPSystTheory->GetErrorYhigh(istartPPextr) );
-  
-      Double_t dataPPLow = ExtractFDSyst( gSigmaPPSystData->GetErrorYlow(istartPPextr), gSigmaPPSystFeedDown->GetErrorYlow(istartPPfd) );
-      if (kRaavsEP) dataPPLow = ExtractFDSyst( 0.5*gSigmaPPSystData->GetErrorYlow(istartPPextr), 0.5*gSigmaPPSystFeedDown->GetErrorYlow(istartPPfd) );
-      systPPLow = TMath::Sqrt( dataPPLow*dataPPLow + gSigmaPPSystTheory->GetErrorYlow(istartPPextr)*gSigmaPPSystTheory->GetErrorYlow(istartPPextr) );
-      if (kRaavsEP) systPPLow = TMath::Sqrt( dataPPLow*dataPPLow + 0.5*gSigmaPPSystTheory->GetErrorYlow(istartPPextr)*0.5*gSigmaPPSystTheory->GetErrorYlow(istartPPextr) );
-
 
-      if(printout) {
-       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;
-       else cout << " pt : "<< pt<<" Syst-pp-data "<< dataPPUp/sigmapp << "%, extr unc + "<< gSigmaPPSystTheory->GetErrorYhigh(istartPPextr)/sigmapp <<" - "<< gSigmaPPSystTheory->GetErrorYlow(istartPPextr)/sigmapp <<" %"<<endl;
+      if(printout && TMath::Abs(ptprintout-pt)<0.1) {
+       cout << " pt : "<< pt<<" Syst-pp-data "<< dataPPUp/sigmapp << "%, ";
+       if(!isExtrapolatedBin){
+         if (isRaavsEP>0.) cout <<" extr unc + "<< 0.5*gSigmaPPSystTheory->GetErrorYhigh(istartPPextr)/sigmapp <<" - "<< 0.5*gSigmaPPSystTheory->GetErrorYlow(istartPPextr)/sigmapp <<" %";
+         else cout <<" extr unc + "<< gSigmaPPSystTheory->GetErrorYhigh(istartPPextr)/sigmapp <<" - "<< gSigmaPPSystTheory->GetErrorYlow(istartPPextr)/sigmapp <<" %";
+       }
+       cout << endl;
       }
 
       //
@@ -550,11 +748,11 @@ void HFPtSpectrumRaa(const char *ppfile="HFPtSpectrum_D0Kpi_method2_rebinnedth_2
       // Data syst : c) combine pp & PbPb
       //
       systLow = sigmapp>0. ? 
-       RaaCharm * TMath::Sqrt( (systABLow/sigmaAB)*(systABLow/sigmaAB) + (systPPUp/sigmapp)*(systPPUp/sigmapp) )//+ (TabSyst/Tab)*(TabSyst/Tab) )
+       RaaCharm * TMath::Sqrt( (systABLow/sigmaAB)*(systABLow/sigmaAB) + (systPPUp/sigmapp)*(systPPUp/sigmapp) )
        : 0.;
 
       systUp = sigmapp>0. ? 
-       RaaCharm * TMath::Sqrt( (systABUp/sigmaAB)*(systABUp/sigmaAB) + (systPPLow/sigmapp)*(systPPLow/sigmapp) )//+ (TabSyst/Tab)*(TabSyst/Tab) )
+       RaaCharm * TMath::Sqrt( (systABUp/sigmaAB)*(systABUp/sigmaAB) + (systPPLow/sigmapp)*(systPPLow/sigmapp) )
        : 0.;
       if ( RaaCharm==0 ) { systPPUp =0.; systPPLow =0.; }
       
@@ -583,7 +781,7 @@ void HFPtSpectrumRaa(const char *ppfile="HFPtSpectrum_D0Kpi_method2_rebinnedth_2
        FDL = RaaCharmFDhigh; FDH = RaaCharmFDlow;
       } 
       
-      if(printout) cout<<" Raa "<<RaaCharm<<", Raa-fd-low "<<RaaCharmFDlow <<", Raa-fd-high "<<RaaCharmFDhigh <<endl;
+      if(printout && TMath::Abs(ptprintout-pt)<0.1) cout<<" Raa "<<RaaCharm<<", Raa-fd-low "<<RaaCharmFDlow <<", Raa-fd-high "<<RaaCharmFDhigh <<endl;
       maxFdSyst = TMath::Abs(FDH - RaaCharm);
       minFdSyst = TMath::Abs(RaaCharm - FDL);
       if ( RaaCharm>0 ) {
@@ -616,12 +814,16 @@ void HFPtSpectrumRaa(const char *ppfile="HFPtSpectrum_D0Kpi_method2_rebinnedth_2
     //
     // Filling Eloss scenarii information
     //
+    //    trick in case not fine enough Rb hypothesis to cope with the min/max range
+    //    if( RaaCharm>0 && ( (ElossHypo >= MinHypo && ElossHypo <=MaxHypo) || ElossHypo == ElossCentral[ hABbin ] ) && RaaBeauty<=MaxRb ) {
+    //      by default better not use it, to monitor when this happens (could affect results)
     if( RaaCharm>0 && ElossHypo >= MinHypo && ElossHypo <=MaxHypo && RaaBeauty<=MaxRb ) {
+
       Double_t Ehigh =  ElossMax[ hABbin ] ;
       Double_t Elow =  ElossMin[ hABbin ] ;
       if ( RaaCharm > Ehigh ) ElossMax[ hABbin ] = RaaCharm ;
       if ( RaaCharm < Elow ) ElossMin[ hABbin ] = RaaCharm ;
-      if(printout) {
+      if(printout && TMath::Abs(ptprintout-pt)<0.1) {
        cout<<" Hypothesis " << ElossHypo << " sigma-AB "<< sigmaAB <<", Raa "<< RaaCharm <<", Raa Eloss max "<< ElossMax[hABbin] <<" Raa Eloss min "<< ElossMin[hABbin] << " Rb="<< RaaBeauty <<endl;
        cout<<"  Rb="<< RaaBeauty <<" max "<< RaaBeautyFDhigh <<" min "<< RaaBeautyFDlow <<endl;
       }
@@ -635,7 +837,7 @@ void HFPtSpectrumRaa(const char *ppfile="HFPtSpectrum_D0Kpi_method2_rebinnedth_2
       Double_t RFDlow = RaaCharmFDlow<RaaCharmFDhigh ? RaaCharmFDlow : RaaCharmFDhigh;
       if ( RFDhigh > FDEhigh ) FDElossMax[ hABbin ] = RFDhigh ;
       if ( RFDlow < FDEmin ) FDElossMin[ hABbin ] = RFDlow ;
-      if(printout) 
+      if(printout && TMath::Abs(ptprintout-pt)<0.1
        cout<<" Hypothesis " << ElossHypo << " sigma-AB "<< sigmaAB <<", Raa FD-max Eloss max "<< FDElossMax[hABbin] <<" Raa FD-min Eloss min "<< FDElossMin[hABbin] <<endl;
     }
 
@@ -839,7 +1041,7 @@ void HFPtSpectrumRaa(const char *ppfile="HFPtSpectrum_D0Kpi_method2_rebinnedth_2
   RaaPlot->SetBottomMargin(0.1);
   RaaPlot->SetTickx();
   RaaPlot->SetTicky();
-  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);
+  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);
   hRaaCanvas->GetXaxis()->SetTitleSize(0.05);
   hRaaCanvas->GetXaxis()->SetTitleOffset(0.9);
   hRaaCanvas->GetYaxis()->SetTitleSize(0.05);
@@ -848,7 +1050,7 @@ void HFPtSpectrumRaa(const char *ppfile="HFPtSpectrum_D0Kpi_method2_rebinnedth_2
   gRAB_Norm->SetFillStyle(1001);
   gRAB_Norm->SetFillColor(kGray+2);
   gRAB_Norm->Draw("2");
-  TLine *line = new TLine(0.0172415,1.0,25.,1.0);
+  TLine *line = new TLine(0.0172415,1.0,40.,1.0);
   line->SetLineStyle(2);
   line->Draw();
   hRABvsPt->SetMarkerColor(kBlue);
@@ -888,11 +1090,13 @@ void HFPtSpectrumRaa(const char *ppfile="HFPtSpectrum_D0Kpi_method2_rebinnedth_2
   legrcb->AddEntry(gRAB_FeedDownSystematics,"Syst. from B feed-down","f");
   legrcb->Draw();
   TLatex* tc;
-  if(decay==1) tc =new TLatex(0.18,0.82,"D^{0},  Pb-Pb    #sqrt{s_{NN}}=2.76 TeV");
-  else if(decay==2) tc =new TLatex(0.18,0.82,"D^{+},  Pb-Pb    #sqrt{s_{NN}}=2.76 TeV");
-  else if(decay==3) tc =new TLatex(0.18,0.82,"D^{*+},  Pb-Pb    #sqrt{s_{NN}}=2.76 TeV");
-  else if(decay==4) tc =new TLatex(0.18,0.82,"D_{s}^{+},  Pb-Pb    #sqrt{s_{NN}}=2.76 TeV");
-  else  tc =new TLatex(0.18,0.82,"any (?) D meson,  Pb-Pb    #sqrt{s_{NN}}=2.76 TeV");
+  TString system = "Pb-Pb   #sqrt{s_{NN}}=2.76 TeV";
+  if( cc==kpPb0100 || cc==kpPb020 || cc==kpPb2040 || cc==kpPb4060 || cc==kpPb60100 ) system = "p-Pb   #sqrt{s_{NN}}=5.023 TeV";
+  if(decay==1) tc =new TLatex(0.18,0.82,Form("D^{0},  %s ",system.Data()));
+  else if(decay==2) tc =new TLatex(0.18,0.82,Form("D^{+},  %s ",system.Data()));
+  else if(decay==3) tc =new TLatex(0.18,0.82,Form("D^{*+},  %s ",system.Data()));
+  else if(decay==4) tc =new TLatex(0.18,0.82,Form("D_{s}^{+},  %s ",system.Data()));
+  else  tc =new TLatex(0.18,0.82,Form("any (?) D meson,  %s ",system.Data()));
   tc->SetNDC();
   tc->SetTextSize(0.038);
   tc->SetTextFont(42);
@@ -1037,6 +1241,7 @@ void HFPtSpectrumRaa(const char *ppfile="HFPtSpectrum_D0Kpi_method2_rebinnedth_2
   gRAB_Norm->Write();
   gRAB_FeedDownSystematicsElossHypothesis->Write();
   gRAB_GlobalSystematics->Write();
+  if(isScaledAndExtrapRef) hCombinedReferenceFlag->Write();
 
   out->Write();
 
@@ -1051,20 +1256,25 @@ Bool_t PbPbDataSyst(AliHFSystErr *syst, Double_t pt, Int_t cc, Double_t &dataSys
   Double_t pidunc=0.;
 
   err = syst->GetTotalSystErr(pt)*syst->GetTotalSystErr(pt);
+  errUp = err ;
   errDown = err ;
 
-  if( cc==k07half || cc==k020 || cc==k010 || cc==k1020 || cc==k2040 ) {
-    if(pt<6) pidunc = 0.15;
-    else pidunc = 0.05;
-    errUp = err + pidunc*pidunc - syst->GetPIDEffErr(pt)*syst->GetPIDEffErr(pt);
-    isOk = true;
-  }
-  else if ( cc==k3050 || cc==k4080 || cc==k4060 || cc==k6080 ){
-    if(pt<3.1) pidunc = 0.10;
-    else pidunc = 0.05;
-    errUp = err + pidunc*pidunc - syst->GetPIDEffErr(pt)*syst->GetPIDEffErr(pt);
-    isOk = true;
+  // Apply an asymetric PID uncertainty for 2010 PbPb data only
+  if( syst->GetRunNumber()==10 && syst->GetCollisionType()==1 ) {
+    if( cc==k07half || cc==k020 || cc==k010 || cc==k1020 || cc==k2040 ) {
+      if(pt<6) pidunc = 0.15;
+      else pidunc = 0.05;
+      errUp = err + pidunc*pidunc - syst->GetPIDEffErr(pt)*syst->GetPIDEffErr(pt);
+      isOk = true;
+    }
+    else if ( cc==k3050 || cc==k4080 || cc==k4060 || cc==k6080 ){
+      if(pt<3.1) pidunc = 0.10;
+      else pidunc = 0.05;
+      errUp = err + pidunc*pidunc - syst->GetPIDEffErr(pt)*syst->GetPIDEffErr(pt);
+      isOk = true;
+    }
   }
+  else { isOk = true; }
 
   dataSystUp = TMath::Sqrt(errUp);
   dataSystDown = TMath::Sqrt(errDown);