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