1 #if !defined(__CINT__) || defined(__MAKECINT__)
10 #include "TGraphAsymmErrors.h"
19 #include "AliHFSystErr.h"
20 #include <Riostream.h>
24 ///////////////////////////////////////////////////////////////////////////////
26 // Macro to compute the Raa, taking as inputs the output of the corrected yields
27 // and the pp reference
29 // R_AB = [ ( dsigma/dpt )_AB / sigma_AB ] / <TAB> * ( dsigma/dpt )_pp
33 // 1. ppfile = reference pp in the same pt binning
34 // 2. ABfile = corrected AB yields
35 // 3. outfile = output file name
36 // 4. decay = decay as in HFSystErr class
37 // 5. sigmaABCINT1B = cross section for normalization (**USE THE SAME AS ON 2.**)
38 // 6. fdMethod = feed-down subtraction method (kNb, kfc)
39 // 7. cc = centrality class
40 // 8. Energy = colliding energy (k276,k55)
41 // 9. MinHypo = minimum energy loss hypothesis (Default 1./3.)
42 // 10. MaxHypo = maximum energy loss hypothesis (Default 3.0)
43 // 11. MaxRb = maximum Raa(b) hypothesis (Default 6.0, won't do anything)
44 // 12. RbHypo : flag to decide whether the Eloss hypothesis is Rb or Rb/Rc
45 // 13. CentralHypo = central energy loss hypothesis, DEFAULT TO 1.0
46 // 14. isRaavsEP = flag to compute the Raa IN/OUT of plane, divides the reference by 2.0
47 // 15. ScaledAndExtrapRef: flag to tag scaled+reference pp-scaled-data
49 // Complains to : Zaida Conesa del Valle
51 ///////////////////////////////////////////////////////////////////////////////
53 enum centrality{ kpp, k07half, kpPb0100, k010, k1020, k020, k2040, k2030, k3040, k4050, k3050, k5060, k4060, k6080, k4080, k5080, k80100 };
54 enum energy{ k276, k5dot023, k55 };
55 enum BFDSubtrMethod { kfc, kNb };
56 enum RaavsEP {kPhiIntegrated, kInPlane, kOutOfPlane};
58 Bool_t printout = false;
59 Double_t ptprintout = 1.5;
60 Double_t NormPPUnc = 0.035;
61 Double_t NormABUnc = 0.037;
62 Bool_t elossFDQuadSum = true;
64 //____________________________________________________________
65 Bool_t PbPbDataSyst(AliHFSystErr *syst, Double_t pt, Int_t cc, Double_t &dataSystUp, Double_t &dataSystDown);
67 //____________________________________________________________
68 Double_t ExtractFDSyst(Double_t total, Double_t fd) {
69 // total^2 = data^2 + fd^2
70 Double_t data2 = total*total - fd*fd ;
71 return TMath::Sqrt( data2 );
74 //____________________________________________________________
75 Int_t FindGraphBin(TGraphAsymmErrors *gr, Double_t pt)
78 Int_t npoints = gr->GetN();
79 for(Int_t i=0; i<=npoints; i++){
82 if ( TMath::Abs ( x - pt ) < 0.4 ) {
93 // R_AB = [ ( dsigma/dpt )_AB / sigma_AB ] / <TAB> * ( dsigma/dpt )_pp
96 //____________________________________________________________
97 void HFPtSpectrumRaa(const char *ppfile="HFPtSpectrum_D0Kpi_method2_rebinnedth_230311_newsigma.root",
98 const char *ABfile="HFPtSpectrum_D0Kpi_PbPbcuts_method2_rebinnedth_230311_newsigma.root",
99 const char *outfile="HFPtSpectrumRaa.root",
101 Double_t sigmaABCINT1B=54.e9,
102 Int_t fdMethod = kNb, Int_t cc=kpp, Int_t Energy=k276,
103 Double_t MinHypo=1./3., Double_t MaxHypo=3.0, Double_t MaxRb=6.0,
104 Bool_t isRbHypo=false, Double_t CentralHypo = 1.0,
105 Bool_t isUseTaaForRaa=true, const char *shadRbcFile="", Int_t nSigmaShad=3.0,
106 Int_t isRaavsEP=kPhiIntegrated, Bool_t isScaledAndExtrapRef=kFALSE)
109 gROOT->Macro("$ALICE_ROOT/PWGHF/vertexingHF/macros/LoadLibraries.C");
112 // Defining the TAB values for the given centrality class
114 Double_t Tab = 1., TabSyst = 0., A=207.2, B=207.2;
115 if ( Energy!=k276 && Energy!=k5dot023) {
116 printf("\n The Tab values for this cms energy have not yet been implemented, please do it ! \n");
122 // Values from Alberica's twiki:
123 // https://twiki.cern.ch/twiki/bin/viewauth/ALICE/CentStudies
124 if ( cc == k07half ) {
125 Tab = 24.81; TabSyst = 0.8037;
126 } else if ( cc == k010 ) {
127 Tab = 23.48; TabSyst = 0.97;
128 } else if ( cc == k1020 ) {
129 Tab = 14.4318; TabSyst = 0.5733;
130 } else if ( cc == k020 ) {
131 Tab = 18.93; TabSyst = 0.74;
132 } else if ( cc == k2040 ) {
133 Tab = 6.86; TabSyst = 0.28;
134 } else if ( cc == k2030 ) {
135 Tab = 8.73769; TabSyst = 0.370219;
136 } else if ( cc == k3040 ) {
137 Tab = 5.02755; TabSyst = 0.22099;
138 } else if ( cc == k4050 ) {
139 Tab = 2.68327; TabSyst = 0.137073;
140 } else if ( cc == k3050 ) {
141 Tab = 3.87011; TabSyst = 0.183847;
142 } else if ( cc == k4060 ) {
143 Tab = 2.00; TabSyst= 0.11;
144 } else if ( cc == k4080 ) {
145 Tab = 1.20451; TabSyst = 0.071843;
146 } else if ( cc == k5060 ) {
147 Tab = 1.32884; TabSyst = 0.0929536;
148 } else if ( cc == k6080 ) {
149 Tab = 0.419; TabSyst = 0.033;
150 } else if ( cc == k5080 ) {
151 Tab = 0.719; TabSyst = 0.054;
152 } else if ( cc == k80100 ){
153 Tab = 0.0690; TabSyst = 0.0062;
156 // pPb Glauber (A. Toia)
157 // https://twiki.cern.ch/twiki/bin/viewauth/ALICE/PACentStudies#Glauber_Calculations_with_sigma
158 else if( cc == kpPb0100 ){
159 Tab = 0.098334; TabSyst = 0.0070679;
164 // Reading the pp file
166 TFile * ppf = new TFile(ppfile,"read");
168 TGraphAsymmErrors * gSigmaPPSyst;
169 TGraphAsymmErrors * gSigmaPPSystData = (TGraphAsymmErrors*)ppf->Get("gScaledDataSystData");
170 TGraphAsymmErrors * gSigmaPPSystTheory = (TGraphAsymmErrors*)ppf->Get("gScaledDataSystExtrap");
171 TGraphAsymmErrors * gSigmaPPSystFeedDown = (TGraphAsymmErrors*)ppf->Get("gScaledDataSystFeedDown");
172 TH1I * hCombinedReferenceFlag;
173 TGraphAsymmErrors * gReferenceFdSyst;
174 if(isScaledAndExtrapRef){
175 hCombinedReferenceFlag = (TH1I*)ppf->Get("hCombinedReferenceFlag");
176 hSigmaPP = (TH1D*)ppf->Get("hReference");
177 gSigmaPPSyst = (TGraphAsymmErrors*)ppf->Get("gReferenceSyst");
178 gReferenceFdSyst = (TGraphAsymmErrors*)ppf->Get("gReferenceFdSyst");
180 hSigmaPP = (TH1D*)ppf->Get("fhScaledData");
181 gSigmaPPSyst = (TGraphAsymmErrors*)ppf->Get("gScaledData");
184 // Call the systematics uncertainty class for a given decay
185 AliHFSystErr *systematicsPP = new AliHFSystErr();
186 systematicsPP->Init(decay);
189 // Reading the AB collisions file
191 TFile * ABf = new TFile(ABfile,"read");
192 TH1D *hSigmaAB = (TH1D*)ABf->Get("histoSigmaCorr");
193 TH2D *hSigmaABRcb = (TH2D*)ABf->Get("histoSigmaCorrRcb");
194 TGraphAsymmErrors * gSigmaABSyst = (TGraphAsymmErrors*)ABf->Get("gSigmaCorr");
195 TGraphAsymmErrors * gSigmaABSystFeedDown = (TGraphAsymmErrors*)ABf->Get("gSigmaCorrConservative");
196 TNtuple * nSigmaAB = (TNtuple*)ABf->Get("fnSigma");
198 TH1D *hMassAB = (TH1D*)ABf->Get("hRECpt");
199 TH1D *hDirectEffptAB = (TH1D*)ABf->Get("hDirectEffpt");
200 TH1D *histofcAB = (TH1D*)ABf->Get("histofc");
202 TH1D* fhStatUncEffcSigmaAB = (TH1D*)ABf->Get("fhStatUncEffcSigma");
203 TH1D* fhStatUncEffbSigmaAB = (TH1D*)ABf->Get("fhStatUncEffbSigma");
204 TH1D* fhStatUncEffcFDAB = (TH1D*)ABf->Get("fhStatUncEffcFD");
205 TH1D* fhStatUncEffbFDAB = (TH1D*)ABf->Get("fhStatUncEffbFD");
207 TH1D* fhStatUncEffcSigmaAB_Raa = (TH1D*)fhStatUncEffcSigmaAB->Clone("fhStatUncEffcSigmaAB_Raa");
208 TH1D* fhStatUncEffbSigmaAB_Raa = (TH1D*)fhStatUncEffbSigmaAB->Clone("fhStatUncEffbSigmaAB_Raa");
209 TH1D* fhStatUncEffcFDAB_Raa = (TH1D*)fhStatUncEffcFDAB->Clone("fhStatUncEffcFDAB_Raa");
210 TH1D* fhStatUncEffbFDAB_Raa = (TH1D*)fhStatUncEffbFDAB->Clone("fhStatUncEffbFDAB_Raa");
211 fhStatUncEffcSigmaAB_Raa->Reset();
212 fhStatUncEffbSigmaAB_Raa->Reset();
213 fhStatUncEffcFDAB_Raa->Reset();
214 fhStatUncEffbFDAB_Raa->Reset();
215 fhStatUncEffcSigmaAB_Raa->SetName("fhStatUncEffcSigmaAB_Raa");
216 fhStatUncEffbSigmaAB_Raa->SetName("fhStatUncEffbSigmaAB_Raa");
217 fhStatUncEffcFDAB_Raa->SetName("fhStatUncEffcFDAB_Raa");
218 fhStatUncEffbFDAB_Raa->SetName("fhStatUncEffvFDAB_Raa");
222 // Call the systematics uncertainty class for a given decay
223 AliHFSystErr *systematicsAB = new AliHFSystErr();
224 systematicsAB->SetCollisionType(1);
225 if ( cc == k07half ) systematicsAB->SetCentrality("07half");
226 else if ( cc == k010 ) systematicsAB->SetCentrality("010");
227 else if ( cc == k1020 ) systematicsAB->SetCentrality("1020");
228 else if ( cc == k020 ) systematicsAB->SetCentrality("020");
229 else if ( cc == k2040 || cc == k2030 || cc == k3040 ) {
230 systematicsAB->SetCentrality("2040");
231 systematicsAB->SetIsPbPb2010EnergyScan(true);
233 else if ( cc == k4060 || cc == k4050 || cc == k5060 ) systematicsAB->SetCentrality("4060");
234 else if ( cc == k6080 || cc == k5080 ) systematicsAB->SetCentrality("6080");
235 else if ( cc == k4080 ) systematicsAB->SetCentrality("4080");
236 else if ( cc == k3050 ) {
237 if (isRaavsEP == kPhiIntegrated) systematicsAB->SetCentrality("4080");
238 else if (isRaavsEP == kInPlane) systematicsAB->SetCentrality("3050InPlane");
239 else if (isRaavsEP == kOutOfPlane) systematicsAB->SetCentrality("3050OutOfPlane");
242 else if ( cc == kpPb0100 ){
243 systematicsAB->SetCollisionType(2);
246 cout << " Systematics not yet implemented " << endl;
250 systematicsAB->Init(decay);
253 Int_t entries = nSigmaAB->GetEntries();
254 Float_t pt=0., signal=0., Rb=0., Rcb=0., fcAB=0., yieldAB=0., sigmaAB=0., statUncSigmaAB=0., sigmaABMin=0.,sigmaABMax=0.;
255 nSigmaAB->SetBranchAddress("pt",&pt);
256 nSigmaAB->SetBranchAddress("Signal",&signal);
257 if (fdMethod==kNb) nSigmaAB->SetBranchAddress("Rb",&Rb);
258 else if (fdMethod==kfc) nSigmaAB->SetBranchAddress("Rcb",&Rcb);
259 nSigmaAB->SetBranchAddress("fc",&fcAB);
260 nSigmaAB->SetBranchAddress("Yield",&yieldAB);
261 nSigmaAB->SetBranchAddress("Sigma",&sigmaAB);
262 nSigmaAB->SetBranchAddress("SigmaStatUnc",&statUncSigmaAB);
263 nSigmaAB->SetBranchAddress("SigmaMax",&sigmaABMax);
264 nSigmaAB->SetBranchAddress("SigmaMin",&sigmaABMin);
267 // define the binning
268 Int_t nbins = hSigmaAB->GetNbinsX();
269 Double_t binwidth = hSigmaAB->GetBinWidth(1);
270 Double_t *limits = new Double_t[nbins+1];
271 Double_t *binwidths = new Double_t[nbins];
273 for (Int_t i=1; i<=nbins; i++) {
274 binwidth = hSigmaAB->GetBinWidth(i);
275 xlow = hSigmaAB->GetBinLowEdge(i);
277 binwidths[i-1] = binwidth;
279 limits[nbins] = xlow + binwidth;
283 // Read the shadowing file if given as input
285 Double_t centralRbcShad[nbins+1], minRbcShad[nbins+1], maxRbcShad[nbins+1];
286 for(Int_t i=0; i<=nbins; i++) { centralRbcShad[i]=1.0; minRbcShad[i]=6.0; maxRbcShad[i]=0.0; }
287 Bool_t isShadHypothesis = false;
288 if( strcmp(shadRbcFile,"")!=0 ) {
289 isShadHypothesis = true;
290 cout<<endl<<">> Beware, using the shadowing prediction file with an "<<nSigmaShad<<"*sigma <<"<<endl<<endl;
291 TFile *fshad = new TFile(shadRbcFile,"read");
292 if(!fshad){ cout <<" >> Shadowing file not properly opened!!!"<<endl<<endl; return;}
293 // TH1D *hRbcShadCentral = (TH1D*)fshad->Get("hDfromBoverPromptD_Shadowing_central");
294 // TH1D *hRbcShadMin = (TH1D*)fshad->Get("hDfromBoverPromptD_Shadowing_upper");
295 // TH1D *hRbcShadMax = (TH1D*)fshad->Get("hDfromBoverPromptD_Shadowing_lower");
296 TH1D *hRbcShadCentral = (TH1D*)fshad->Get("hDfromBoverDfromc_L0");
297 TH1D *hRbcShadMin = (TH1D*)fshad->Get("hDfromBoverDfromc_L0");
298 TH1D *hRbcShadMax = (TH1D*)fshad->Get("hDfromBoverDfromc_L1");
299 if(!hRbcShadCentral || !hRbcShadMin || !hRbcShadMax) {
300 cout<< endl <<">> Shadowing input histograms are not ok !! "<<endl<<endl;
305 for(Int_t i=1; i<=nbins; i++) {
306 Double_t xpt = hSigmaAB->GetBinCenter(i);
308 centralRbcShad[i] = hRbcShadCentral->GetBinContent( hRbcShadCentral->FindBin(xpt) );
309 Double_t minValue0 = hRbcShadMin->GetBinContent( hRbcShadMin->FindBin(xpt) );
310 Double_t maxValue0 = hRbcShadMax->GetBinContent( hRbcShadMax->FindBin(xpt) );
311 Double_t arrayEl[3] = {minValue0,maxValue0, centralRbcShad[i]};
312 Double_t minValue = TMath::MinElement(3,arrayEl);
313 Double_t maxValue = TMath::MaxElement(3,arrayEl);
314 cout<<">> Shadowing pt="<<xpt<<" central="<<centralRbcShad[i]<<" min="<<minValue<<" max="<<maxValue<<endl;
315 if(minValue>centralRbcShad[i]){ minValue = centralRbcShad[i]; }
316 if(maxValue<centralRbcShad[i]){ maxValue = centralRbcShad[i]; }
317 minRbcShad[i] = centralRbcShad[i] - nSigmaShad*(centralRbcShad[i] - minValue);
318 maxRbcShad[i] = centralRbcShad[i] + nSigmaShad*(maxValue - centralRbcShad[i]);
319 cout<<">> Shadowing hypothesis pt="<<xpt<<" central="<<centralRbcShad[i]<<" min="<<minRbcShad[i]<<" max="<<maxRbcShad[i]<<endl;
324 // define the bins correspondence bw histos/files/graphs
327 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.);
328 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.);
329 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.);
330 Int_t nbinsHypo=800;//200;
331 Double_t *limitsHypo = new Double_t[nbinsHypo+1];
332 for(Int_t i=1; i<=nbinsHypo+1; i++) limitsHypo[i-1]= i*4./800.;
333 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);
334 TH2D *hRCharmVsRBeauty[nbins+1];
335 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);
336 TH2D *hRCharmVsElossHypo[nbins+1];
337 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);
339 TH1D *hRABEloss00= new TH1D("hRABEloss00","hRABEloss00",nbins,limits);
340 TH1D *hRABEloss05= new TH1D("hRABEloss05","hRABEloss05",nbins,limits);
341 TH1D *hRABEloss10= new TH1D("hRABEloss10","hRABEloss10",nbins,limits);
342 TH1D *hRABEloss15= new TH1D("hRABEloss15","hRABEloss15",nbins,limits);
343 TH1D *hRABEloss20= new TH1D("hRABEloss20","hRABEloss20",nbins,limits);
345 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.);
346 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.);
348 TH1D * hRABvsRbFDhigh_proj = new TH1D("hRABvsRbFDhigh_proj","hRABvsRbFDhigh_proj",nbins,limits);
349 TH1D * hRABvsRbFDlow_proj = new TH1D("hRABvsRbFDlow_proj","hRABvsRbFDlow_proj",nbins,limits);
351 TNtuple *ntupleRAB=0x0 ;
353 ntupleRAB = new TNtuple("ntupleRAB","ntupleRAB (Nb)","pt:TAB:sigmaPP:sigmaAB:invyieldAB:invyieldABFDHigh:invyieldABFDLow:RABCharm:RABCharmFDHigh:RABCharmFDLow:RABBeauty:fc",100000);
354 } else if (fdMethod==kfc) {
355 ntupleRAB = new TNtuple("ntupleRAB","ntupleRAB (fc)","pt:TAB:sigmaPP:sigmaAB:invyieldAB:invyieldABFDHigh:invyieldABFDLow:Rcb:RABCharm:RABCharmFDHigh:RABCharmFDLow:RABBeauty:RABBeautyFDHigh:RABBeautyFDLow:fc",100000);
357 if(!ntupleRAB) printf("ERROR: Wrong method option");
359 TH1D * hYieldABvsPt = new TH1D("hYieldABvsPt"," Yield_{AB}(c) vs p_{T} (no Eloss hypothesis); p_{T} [GeV/c] ; Yield_{charm} ",nbins,limits);
360 TH1D * hRABvsPt = new TH1D("hRABvsPt"," R_{AB}(c) vs p_{T} (no Eloss hypothesis); p_{T} [GeV/c] ; R_{charm} ",nbins,limits);
361 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);
362 TGraphAsymmErrors *gRAB_ElossHypothesis = new TGraphAsymmErrors(nbins+1);
363 gRAB_ElossHypothesis->SetNameTitle("gRAB_ElossHypothesis","RAB Eloss systematics");
364 TGraphAsymmErrors *gRAB_FeedDownSystematics = new TGraphAsymmErrors(nbins+1);
365 gRAB_FeedDownSystematics->SetNameTitle("gRAB_FeedDownSystematics","RAB Feed-Down systematics");
366 TGraphAsymmErrors *gRAB_fcFeedDownOnly = new TGraphAsymmErrors(nbins+1);
367 gRAB_fcFeedDownOnly->SetNameTitle("gRAB_fcFeedDownOnly","RAB fc Feed-Down Only");
368 TGraphAsymmErrors *gRAB_FeedDownSystematicsElossHypothesis = new TGraphAsymmErrors(nbins+1);
369 gRAB_FeedDownSystematicsElossHypothesis->SetNameTitle("gRAB_FeedDownSystematicsElossHypothesis","RAB Feed-Down systematics considering Eloss hypothesis");
370 TGraphAsymmErrors *gRAB_DataSystematics = new TGraphAsymmErrors(nbins+1);
371 gRAB_DataSystematics->SetNameTitle("gRAB_DataSystematics","RAB Measurement (no FD, no Eloss) systematics");
372 TGraphAsymmErrors *gRAB_DataSystematicsPP = new TGraphAsymmErrors(nbins+1);
373 gRAB_DataSystematicsPP->SetNameTitle("gRAB_DataSystematicsPP","RAB Measurement PP meas. systematics (data+scaling)");
374 TGraphAsymmErrors *gRAB_DataSystematicsAB = new TGraphAsymmErrors(nbins+1);
375 gRAB_DataSystematicsAB->SetNameTitle("gRAB_DataSystematicsAB","RAB Measurement AB (no FD, no Eloss, no PP data) systematics");
376 TGraphAsymmErrors *gRAB_GlobalSystematics = new TGraphAsymmErrors(nbins+1);
377 gRAB_GlobalSystematics->SetNameTitle("gRAB_GlobalSystematics","RAB Measurement global (data, FD, Eloss) systematics");
378 Double_t ElossMax[nbins+1], ElossMin[nbins+1];
379 for(Int_t i=0; i<=nbins; i++) { ElossMax[i]=0.; ElossMin[i]=6.; }
380 Double_t fcElossMax[nbins+1], fcElossMin[nbins+1];
381 for(Int_t i=0; i<=nbins; i++) { fcElossMax[i]=0.; fcElossMin[i]=6.; }
382 Double_t FDElossMax[nbins+1], FDElossMin[nbins+1];
383 for(Int_t i=0; i<=nbins; i++) { FDElossMax[i]=0.; FDElossMin[i]=6.; }
385 TGraphAsymmErrors *gRAB_Norm = new TGraphAsymmErrors(1);
386 gRAB_Norm->SetNameTitle("gRAB_Norm","RAB Normalization systematics (pp norm + Tab)");
387 Double_t normUnc = TMath::Sqrt ( NormPPUnc*NormPPUnc + (TabSyst/Tab)*(TabSyst/Tab) );
388 if(!isUseTaaForRaa) normUnc = TMath::Sqrt ( NormPPUnc*NormPPUnc + NormABUnc*NormABUnc );
389 gRAB_Norm->SetPoint(1,0.5,1.);
390 gRAB_Norm->SetPointError(1,0.25,0.25,normUnc,normUnc);
393 // R_AB = ( dN/dpt )_AB / <Ncoll_AB> * ( dN/dpt )_pp ; <Ncoll> = <Tab> * sigma_NN^inel
394 // R_AB = [ ( dsigma/dpt )_AB / sigma_AB ] / <TAB> * ( dsigma/dpt )_pp
396 Int_t istartPPfd=0, istartPPsyst=0, istartABfd=0, istartPPextr=0;
397 Double_t yPPh=0., yPPl=0., yABh=0., yABl=0.;
398 Double_t RaaCharm =0., RaaBeauty=0.;
399 Double_t RaaCharmFDhigh = 0., RaaCharmFDlow = 0.;
400 Double_t RaaBeautyFDhigh = 0., RaaBeautyFDlow = 0.;
401 Double_t systUp=0., systLow=0., systPPUp=0., systPPLow=0., systABUp=0., systABLow=0.;
404 // Search the central value of the energy loss hypothesis Rb = Rc (bin)
406 Double_t ElossCentral[nbins+1];
407 for(Int_t i=0; i<=nbins; i++) { ElossCentral[i]=0.; }
409 for(Int_t ientry=0; ientry<=entries; ientry++){
411 nSigmaAB->GetEntry(ientry);
412 // cout << " pt="<< pt<<" sigma-AB="<<sigmaAB<<endl;
413 if ( !(sigmaAB>0.) ) continue;
414 //if(decay==2 && pt<2.) continue;
416 // Compute RAB and the statistical uncertainty
417 Int_t hppbin = hSigmaPP->FindBin( pt );
418 Int_t hABbin = hSigmaAB->FindBin( pt );
419 Double_t sigmapp = hSigmaPP->GetBinContent( hppbin );
420 // cout << " pt="<< pt<<", sigma-pp="<< sigmapp<<endl;
421 if (isRaavsEP>0.) sigmapp = 0.5*sigmapp;
422 if ( !(sigmapp>0.) ) continue;
424 RaaCharm = ( sigmaAB / sigmaABCINT1B ) / ((Tab*1e3) * sigmapp *1e-12 ) ;
425 if(!isUseTaaForRaa) {
426 RaaCharm = ( sigmaAB ) / ( (A*B) * sigmapp ) ;
432 else if (fdMethod==kfc) {
433 RaaBeauty = ( RaaCharm / Rcb ) ;
436 Double_t ElossHypo = 0.;
437 if (fdMethod==kfc) { ElossHypo = 1. / Rcb; }
438 else { ElossHypo = 1. / (RaaCharm / RaaBeauty) ; }
439 if(isRbHypo) ElossHypo = RaaBeauty;
441 // If using shadowing hypothesis, change the central hypothesis too
442 if(isShadHypothesis) CentralHypo = centralRbcShad[hABbin];
444 // cout <<" pt "<< pt << " Raa charm " << RaaCharm << " Raa beauty " << RaaBeauty << " eloss hypo "<< ElossHypo<<endl;
446 // Find the bin for the central Eloss hypo
448 if( TMath::Abs( ElossHypo - CentralHypo ) < 0.075 ){
449 Double_t DeltaIni = TMath::Abs( ElossCentral[ hABbin ] - CentralHypo );
450 Double_t DeltaV = TMath::Abs( ElossHypo - CentralHypo );
451 // cout << " pt " << pt << " ECentral " << ElossCentral[ hABbin ] << " Ehypo "<< ElossHypo ;
452 if ( DeltaV < DeltaIni ) ElossCentral[ hABbin ] = ElossHypo;
453 // cout << " final ECentral " << ElossCentral[ hABbin ] << endl;
457 // Calculation of the Raa and its uncertainties
459 for(Int_t ientry=0; ientry<entries; ientry++){
461 nSigmaAB->GetEntry(ientry);
462 if ( !(sigmaAB>0.) ) continue;
463 // if ( pt<2 || pt>16) continue;
466 // Compute RAB and the statistical uncertainty
467 Int_t hppbin = hSigmaPP->FindBin( pt );
468 Double_t sigmapp = hSigmaPP->GetBinContent( hppbin );
469 if (isRaavsEP>0.) sigmapp = 0.5*sigmapp;
470 if ( !(sigmapp>0.) ) continue;
472 RaaCharm = ( sigmaAB / sigmaABCINT1B ) / ((Tab*1e3) * sigmapp *1e-12 );
473 if(!isUseTaaForRaa) {
474 RaaCharm = ( sigmaAB ) / ( (A*B) * sigmapp ) ;
477 // Flag to know if it is an scaled or extrapolated point of the pp reference
478 Bool_t isExtrapolatedBin = kFALSE;
479 if(isScaledAndExtrapRef) isExtrapolatedBin = hCombinedReferenceFlag->GetBinContent( hppbin );
481 istartPPsyst = FindGraphBin(gSigmaPPSyst,pt);
484 // FONLL Feed-Down systematics
487 if(!isExtrapolatedBin) istartPPfd = FindGraphBin(gSigmaPPSystFeedDown,pt);
489 istartABfd = FindGraphBin(gSigmaABSystFeedDown,pt);
491 // cout << " Starting bin for pp is "<< istartPPfd <<", for AA is "<<istartABfd << endl;
492 if(isExtrapolatedBin){
493 Int_t ibinfd = FindGraphBin(gReferenceFdSyst,pt);
494 yPPh = gReferenceFdSyst->GetErrorYhigh(ibinfd);
495 yPPl = gReferenceFdSyst->GetErrorYlow(ibinfd);
497 yPPh = gSigmaPPSystFeedDown->GetErrorYhigh(istartPPfd);
498 yPPl = gSigmaPPSystFeedDown->GetErrorYlow(istartPPfd);
505 yABh = gSigmaABSystFeedDown->GetErrorYhigh(istartABfd);
506 yABl = gSigmaABSystFeedDown->GetErrorYlow(istartABfd);
509 RaaCharmFDhigh = ( sigmaABMax / sigmaABCINT1B ) / ((Tab*1e3) * (sigmapp+yPPh) *1e-12 ) ;
510 RaaCharmFDlow = ( sigmaABMin / sigmaABCINT1B ) / ((Tab*1e3) * (sigmapp-yPPl) *1e-12 ) ;
511 if(printout && TMath::Abs(ptprintout-pt)<0.1 ) cout << endl<<" pt "<< pt << " Raa " << RaaCharm <<" high "<< RaaCharmFDhigh << " low "<< RaaCharmFDlow<<endl;
512 if(!isUseTaaForRaa) {
513 RaaCharmFDhigh = ( sigmaABMax ) / ( (A*B)* (sigmapp+yPPh) ) ;
514 RaaCharmFDlow = ( sigmaABMin ) / ( (A*B)* (sigmapp-yPPl) ) ;
520 RaaBeautyFDlow = Rb ;
521 RaaBeautyFDhigh = Rb ;
522 ntupleRAB->Fill( pt, Tab*1e3, sigmapp*1e-12, sigmaAB*1e-12, sigmaAB/sigmaABCINT1B,
523 sigmaABMax / sigmaABCINT1B, sigmaABMin / sigmaABCINT1B,
524 RaaCharm, RaaCharmFDhigh, RaaCharmFDlow, RaaBeauty, fcAB );
526 else if (fdMethod==kfc) {
527 RaaBeauty = ( RaaCharm / Rcb ) ;
528 RaaBeautyFDlow = ( RaaCharmFDlow / Rcb ) ;
529 RaaBeautyFDhigh = ( RaaCharmFDhigh / Rcb ) ;
530 hRABvsRcb->Fill( pt, RaaCharm, RaaBeauty );
531 ntupleRAB->Fill( pt, Tab*1e3, sigmapp*1e-12, sigmaAB*1e-12, sigmaAB/sigmaABCINT1B,
532 sigmaABMax / sigmaABCINT1B, sigmaABMin / sigmaABCINT1B,
533 Rcb, RaaCharm, RaaCharmFDhigh, RaaCharmFDlow, RaaBeauty, RaaBeautyFDhigh, RaaBeautyFDlow, fcAB );
535 hRABvsRb->Fill( pt, RaaCharm, RaaBeauty );
536 hRABvsRbFDlow->Fill( pt, RaaCharmFDlow, RaaBeautyFDlow );
537 hRABvsRbFDhigh->Fill( pt, RaaCharmFDhigh, RaaBeautyFDhigh );
538 if(printout && TMath::Abs(ptprintout-pt)<0.1) cout << " pt "<< pt << " Rb " << RaaBeauty <<" high "<< RaaBeautyFDhigh << " low "<< RaaBeautyFDlow <<endl;
540 hRABCharmVsRBeautyVsPt->Fill( pt, RaaBeauty, RaaCharm );
541 Int_t ptbin = hRABvsPt->FindBin( pt );
542 hRCharmVsRBeauty[ptbin]->Fill( RaaBeauty, RaaCharm );
543 hRCharmVsRBeauty[ptbin]->Fill( RaaBeautyFDlow, RaaCharmFDlow );
544 hRCharmVsRBeauty[ptbin]->Fill( RaaBeautyFDhigh, RaaCharmFDhigh );
548 if( TMath::Abs(Rcb-0.015)<0.009 ) hRABEloss00->Fill(pt,RaaCharm);
549 if( TMath::Abs(Rcb-0.5)<0.009 ) hRABEloss05->Fill(pt,RaaCharm);
550 if( TMath::Abs(Rcb-1.0)<0.009 ) {
551 hRABEloss10->Fill(pt,RaaCharm);
552 hRABvsRbFDhigh_proj->Fill(pt,RaaCharmFDhigh);
553 hRABvsRbFDlow_proj->Fill(pt,RaaCharmFDlow);
555 if( TMath::Abs(Rcb-1.5)<0.009 ) hRABEloss15->Fill(pt,RaaCharm);
556 if( TMath::Abs(Rcb-2.0)<0.009 ) hRABEloss20->Fill(pt,RaaCharm);
558 else if (fdMethod==kNb) {
559 if( TMath::Abs(RaaBeauty-0.015)<0.009 ) hRABEloss00->Fill(pt,RaaCharm);
560 if( TMath::Abs(RaaBeauty-0.5)<0.009 ) hRABEloss05->Fill(pt,RaaCharm);
561 if( TMath::Abs(RaaBeauty-1.0)<0.009 ) {
562 hRABEloss10->Fill(pt,RaaCharm);
563 hRABvsRbFDhigh_proj->Fill(pt,RaaCharmFDhigh);
564 hRABvsRbFDlow_proj->Fill(pt,RaaCharmFDlow);
566 if( TMath::Abs(RaaBeauty-1.5)<0.009 ) hRABEloss15->Fill(pt,RaaCharm);
567 if( TMath::Abs(RaaBeauty-2.0)<0.009 ) hRABEloss20->Fill(pt,RaaCharm);
571 Int_t hABbin = hMassAB->FindBin( pt );
572 if(isShadHypothesis) CentralHypo = centralRbcShad[hABbin];
574 if(printout && TMath::Abs(ptprintout-pt)<0.1)
575 if ( fdMethod==kNb && TMath::Abs(Rb -CentralHypo)< 0.05) {
576 cout << " pt "<< pt <<", at bin "<<hABbin<<endl;
577 cout<<" entries "<<entries<<", i="<<ientry<<", pt="<<pt<<", Rb="<<Rb<<", Tab="<<Tab<<", sigmaAB="<<sigmaAB<<", sigmapp="<<sigmapp<<", Raacharm="<<RaaCharm<<", RaaBeauty="<<RaaBeauty<<endl;
578 cout << " AB basis: mass "<< hMassAB->GetBinContent(hABbin)<<", eff "<< hDirectEffptAB->GetBinContent(hABbin)<<endl;
579 cout<<" FD low, err low AB "<< (sigmaAB-sigmaABMin)<<" err low PP "<< yPPl<<" Raacharm="<<RaaCharmFDlow<<", RaaBeauty="<<RaaBeautyFDlow<<endl;
580 cout<<" FD high, err high AB "<< (sigmaABMax-sigmaAB)<<" err high PP "<< yPPh<<" Raacharm="<<RaaCharmFDhigh<<", RaaBeauty="<<RaaBeautyFDhigh<<endl;
582 if(printout && TMath::Abs(ptprintout-pt)<0.1)
583 if ( fdMethod==kfc) if(TMath::Abs(Rcb -CentralHypo)< 0.05 ){
584 cout << " pt "<< pt <<", at bin "<<hABbin<<endl;
585 cout<<" entries "<<entries<<", i="<<ientry<<", pt="<<pt<<", Rcb="<<Rcb<<", Tab="<<Tab<<", sigmaAB="<<sigmaAB<<", sigmapp="<<sigmapp<<", Raacharm="<<RaaCharm<<", RaaBeauty="<<RaaBeauty<<endl;
586 cout << " AB basis: mass "<< hMassAB->GetBinContent(hABbin)<<", eff "<< hDirectEffptAB->GetBinContent(hABbin)<<", fc "<<histofcAB->GetBinContent(hABbin)<< endl;
587 cout<<" FD low, err low AB "<< (sigmaAB-sigmaABMin)<<" err low PP "<< yPPl<<" Raacharm="<<RaaCharmFDlow<<", RaaBeauty="<<RaaBeautyFDlow<<endl;
588 cout<<" FD high, err high AB "<< (sigmaABMax-sigmaAB)<<" err high PP "<< yPPh<<" Raacharm="<<RaaCharmFDhigh<<", RaaBeauty="<<RaaBeautyFDhigh<<endl;
593 // Fill in the global properties ?
595 Double_t ElossHypo = 0.;
596 if (fdMethod==kfc) { ElossHypo = 1./ Rcb; }
597 else { ElossHypo = 1. / (RaaCharm / RaaBeauty); }
598 if(isRbHypo) ElossHypo = RaaBeauty;
599 hRCharmVsElossHypo[ptbin]->Fill( ElossHypo, RaaCharm );
601 // If using shadowing hypothesis, change the limit hypothesis too
602 if(isShadHypothesis) {
603 MinHypo = minRbcShad[ hABbin ];
604 MaxHypo = maxRbcShad[ hABbin ];
607 // cout <<" pt "<< pt << " Raa charm " << RaaCharm << " Raa beauty " << RaaBeauty << " eloss hypo "<< ElossHypo
608 if(ientry==0) cout<<" pt"<< pt<< " ElossCentral "<< ElossCentral[hABbin] << " min-hypo "<<MinHypo << " max-hypo "<<MaxHypo<<endl;
611 // Fill in histos charm (null Eloss hypothesis)
613 Double_t minFdSyst = 0., maxFdSyst = 0.;
614 if ( ElossHypo == ElossCentral[ hABbin ] ) {
617 // Data stat uncertainty
619 Double_t sigmappStat = hSigmaPP->GetBinError( hppbin );
620 if (isRaavsEP>0.) sigmappStat = sigmappStat*0.5;
621 Int_t hRABbin = hRABvsPt->FindBin( pt );
622 Double_t stat = RaaCharm * TMath::Sqrt( (statUncSigmaAB/sigmaAB)*(statUncSigmaAB/sigmaAB) +
623 (sigmappStat/sigmapp)*(sigmappStat/sigmapp) ) ;
624 if ( RaaCharm==0 ) stat =0.;
626 hRABvsPt->SetBinContent( hRABbin, RaaCharm );
627 hRABvsPt->SetBinError( hRABbin, stat );
628 hYieldABvsPt->SetBinContent( hRABbin, sigmaAB/sigmaABCINT1B );
629 hYieldABvsPt->SetBinError( hRABbin, statUncSigmaAB/sigmaABCINT1B );
631 cout << "pt="<< pt<< " Raa " << RaaCharm << " stat unc. "<<
632 " sigma-pp "<< sigmapp <<" sigma-AB "<< sigmaAB<<endl;
633 if(printout && TMath::Abs(ptprintout-pt)<0.1) {
634 cout << " Raa " << RaaCharm << " stat unc. "<< stat << " is "<< stat/RaaCharm * 100. <<
635 "%, stat-pp "<< sigmappStat/sigmapp*100. <<"% stat-AB "<< statUncSigmaAB/sigmaAB*100.<<"%"<<endl;
638 Double_t errstatEff = fhStatUncEffcSigmaAB->GetBinError( hRABbin );
639 fhStatUncEffcSigmaAB_Raa->SetBinError( hRABbin, errstatEff*RaaCharm );
640 errstatEff = fhStatUncEffbSigmaAB->GetBinError( hRABbin );
641 fhStatUncEffbSigmaAB_Raa->SetBinError( hRABbin, errstatEff*RaaCharm );
642 errstatEff = fhStatUncEffcFDAB->GetBinError( hRABbin );
643 fhStatUncEffcFDAB_Raa->SetBinError( hRABbin, errstatEff*RaaCharm );
644 errstatEff = fhStatUncEffbFDAB->GetBinError( hRABbin );
645 fhStatUncEffbFDAB_Raa->SetBinError( hRABbin, errstatEff*RaaCharm );
651 // Data systematics (sigma syst-but FD + extrap) syst
654 // Data syst: a) Syst in p-p
656 Double_t ptwidth = hSigmaAB->GetBinWidth(hABbin) / 2. ;
658 if(!isExtrapolatedBin) istartPPextr = FindGraphBin(gSigmaPPSystTheory,pt);
660 Double_t dataPPUp=0., dataPPLow=0.;
661 if(isExtrapolatedBin) {
662 dataPPUp = gSigmaPPSyst->GetErrorYhigh(istartPPsyst);
663 dataPPLow = gSigmaPPSyst->GetErrorYlow(istartPPsyst);
665 systPPLow = dataPPLow;
667 dataPPUp = ExtractFDSyst( gSigmaPPSystData->GetErrorYhigh(istartPPextr), gSigmaPPSystFeedDown->GetErrorYhigh(istartPPfd) );
668 dataPPLow = ExtractFDSyst( gSigmaPPSystData->GetErrorYlow(istartPPextr), gSigmaPPSystFeedDown->GetErrorYlow(istartPPfd) );
669 systPPUp = TMath::Sqrt( dataPPUp*dataPPUp + gSigmaPPSystTheory->GetErrorYhigh(istartPPextr)*gSigmaPPSystTheory->GetErrorYhigh(istartPPextr) );
670 systPPLow = TMath::Sqrt( dataPPLow*dataPPLow + gSigmaPPSystTheory->GetErrorYlow(istartPPextr)*gSigmaPPSystTheory->GetErrorYlow(istartPPextr) );
673 dataPPUp = dataPPUp*0.5;
674 dataPPLow = dataPPLow*0.5;
675 if(isExtrapolatedBin) {
677 systPPLow = dataPPLow;
679 systPPUp = TMath::Sqrt( dataPPUp*dataPPUp + 0.5*gSigmaPPSystTheory->GetErrorYhigh(istartPPextr)*0.5*gSigmaPPSystTheory->GetErrorYhigh(istartPPextr) );
680 systPPLow = TMath::Sqrt( dataPPLow*dataPPLow + 0.5*gSigmaPPSystTheory->GetErrorYlow(istartPPextr)*0.5*gSigmaPPSystTheory->GetErrorYlow(istartPPextr) );
685 if(printout && TMath::Abs(ptprintout-pt)<0.1) {
686 cout << " pt : "<< pt<<" Syst-pp-data "<< dataPPUp/sigmapp << "%, ";
687 if(!isExtrapolatedBin){
688 if (isRaavsEP>0.) cout <<" extr unc + "<< 0.5*gSigmaPPSystTheory->GetErrorYhigh(istartPPextr)/sigmapp <<" - "<< 0.5*gSigmaPPSystTheory->GetErrorYlow(istartPPextr)/sigmapp <<" %";
689 else cout <<" extr unc + "<< gSigmaPPSystTheory->GetErrorYhigh(istartPPextr)/sigmapp <<" - "<< gSigmaPPSystTheory->GetErrorYlow(istartPPextr)/sigmapp <<" %";
695 // Data syst: b) Syst in PbPb
697 Double_t dataSystUp=0., dataSystDown=0.;
698 Bool_t PbPbDataSystOk = PbPbDataSyst(systematicsAB,pt,cc,dataSystUp,dataSystDown);
699 if (!PbPbDataSystOk) { cout <<" There is some issue with the PbPb data systematics, please check and rerun"<<endl; return; }
700 systABUp = sigmaAB * TMath::Sqrt( dataSystUp*dataSystUp +
701 (hDirectEffptAB->GetBinError(hABbin)/hDirectEffptAB->GetBinContent(hABbin))*(hDirectEffptAB->GetBinError(hABbin)/hDirectEffptAB->GetBinContent(hABbin)) );
703 systABLow = sigmaAB * TMath::Sqrt( dataSystDown*dataSystDown +
704 (hDirectEffptAB->GetBinError(hABbin)/hDirectEffptAB->GetBinContent(hABbin))*(hDirectEffptAB->GetBinError(hABbin)/hDirectEffptAB->GetBinContent(hABbin)) );
706 // Data syst : c) combine pp & PbPb
708 systLow = sigmapp>0. ?
709 RaaCharm * TMath::Sqrt( (systABLow/sigmaAB)*(systABLow/sigmaAB) + (systPPUp/sigmapp)*(systPPUp/sigmapp) )
712 systUp = sigmapp>0. ?
713 RaaCharm * TMath::Sqrt( (systABUp/sigmaAB)*(systABUp/sigmaAB) + (systPPLow/sigmapp)*(systPPLow/sigmapp) )
715 if ( RaaCharm==0 ) { systPPUp =0.; systPPLow =0.; }
718 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;
721 hRABvsPt_DataSystematics->SetBinContent( hRABbin, RaaCharm );
722 hRABvsPt_DataSystematics->SetBinError( hRABbin, systUp );
723 gRAB_DataSystematics->SetPoint( hABbin, pt, RaaCharm ); // i, x, y
724 gRAB_DataSystematics->SetPointError( hABbin, ptwidth, ptwidth, systLow, systUp );
725 gRAB_DataSystematics->SetPointEXlow(hABbin, 0.4); gRAB_DataSystematics->SetPointEXhigh(hABbin,0.4);
726 gRAB_DataSystematicsPP->SetPoint( hABbin, pt, RaaCharm ); // i, x, y
727 gRAB_DataSystematicsPP->SetPointError( hABbin, ptwidth, ptwidth, RaaCharm *(systPPUp/sigmapp), RaaCharm *systPPLow/sigmapp );
728 gRAB_DataSystematicsAB->SetPoint( hABbin, pt, RaaCharm ); // i, x, y
729 gRAB_DataSystematicsAB->SetPointError( hABbin, ptwidth, ptwidth, RaaCharm *systABLow/sigmaAB, RaaCharm *systABUp/sigmaAB );
733 // Feed-down Systematics
735 Double_t FDL=0., FDH=0.;
736 if ( RaaCharmFDhigh > RaaCharmFDlow ){
737 FDH = RaaCharmFDhigh; FDL = RaaCharmFDlow;
739 FDL = RaaCharmFDhigh; FDH = RaaCharmFDlow;
742 if(printout && TMath::Abs(ptprintout-pt)<0.1) cout<<" Raa "<<RaaCharm<<", Raa-fd-low "<<RaaCharmFDlow <<", Raa-fd-high "<<RaaCharmFDhigh <<endl;
743 maxFdSyst = TMath::Abs(FDH - RaaCharm);
744 minFdSyst = TMath::Abs(RaaCharm - FDL);
746 gRAB_FeedDownSystematics->SetPoint( hABbin, pt, RaaCharm ); // i, x, y
747 gRAB_FeedDownSystematics->SetPointError( hABbin, 0.3, 0.3, minFdSyst, maxFdSyst ); // i, x, y
748 gRAB_fcFeedDownOnly->SetPoint( hABbin, pt,fcAB );
749 gRAB_fcFeedDownOnly->SetPointError(hABbin, 0.3, 0.3, fcAB-(sigmaABMin/sigmaAB*fcAB), (sigmaABMax/sigmaAB*fcAB)-fcAB );
753 cout<<" FD syst +"<< maxFdSyst/RaaCharm <<" - "<<minFdSyst/RaaCharm<<endl;
754 cout<<" fc = "<<fcAB<<", ("<< sigmaABMax/sigmaAB * fcAB <<","<< sigmaABMin/sigmaAB * fcAB <<")"<<endl;
758 // Filling part of the Eloss scenarii information
761 gRAB_ElossHypothesis->SetPoint( hABbin, pt, RaaCharm ); // i, x, y
762 gRAB_ElossHypothesis->SetPointEXlow( hABbin, ptwidth);
763 gRAB_ElossHypothesis->SetPointEXhigh( hABbin, ptwidth);
764 gRAB_FeedDownSystematicsElossHypothesis->SetPoint( hABbin, pt, RaaCharm ); // i, x, y
765 gRAB_FeedDownSystematicsElossHypothesis->SetPointEXlow( hABbin, ptwidth);
766 gRAB_FeedDownSystematicsElossHypothesis->SetPointEXhigh( hABbin, ptwidth);
767 gRAB_GlobalSystematics->SetPoint( hABbin, pt, RaaCharm ); // i, x, y
768 gRAB_GlobalSystematics->SetPointEXlow(hABbin,0.4); gRAB_GlobalSystematics->SetPointEXhigh(hABbin,0.4);
773 // Filling Eloss scenarii information
775 // trick in case not fine enough Rb hypothesis to cope with the min/max range
776 // if( RaaCharm>0 && ( (ElossHypo >= MinHypo && ElossHypo <=MaxHypo) || ElossHypo == ElossCentral[ hABbin ] ) && RaaBeauty<=MaxRb ) {
777 // by default better not use it, to monitor when this happens (could affect results)
778 if( RaaCharm>0 && ElossHypo >= MinHypo && ElossHypo <=MaxHypo && RaaBeauty<=MaxRb ) {
780 Double_t Ehigh = ElossMax[ hABbin ] ;
781 Double_t Elow = ElossMin[ hABbin ] ;
782 if ( RaaCharm > Ehigh ) ElossMax[ hABbin ] = RaaCharm ;
783 if ( RaaCharm < Elow ) ElossMin[ hABbin ] = RaaCharm ;
784 if(printout && TMath::Abs(ptprintout-pt)<0.1) {
785 cout<<" Hypothesis " << ElossHypo << " sigma-AB "<< sigmaAB <<", Raa "<< RaaCharm <<", Raa Eloss max "<< ElossMax[hABbin] <<" Raa Eloss min "<< ElossMin[hABbin] << " Rb="<< RaaBeauty <<endl;
786 cout<<" Rb="<< RaaBeauty <<" max "<< RaaBeautyFDhigh <<" min "<< RaaBeautyFDlow <<endl;
788 Double_t fcEhigh = fcElossMax[ hABbin ] ;
789 Double_t fcElow = fcElossMin[ hABbin ] ;
790 if ( fcAB > fcEhigh ) fcElossMax[ hABbin ] = fcAB ;
791 if ( fcAB < fcElow ) fcElossMin[ hABbin ] = fcAB ;
792 Double_t FDEhigh = FDElossMax[ hABbin ];
793 Double_t FDEmin = FDElossMin[ hABbin ];
794 Double_t RFDhigh = RaaCharmFDhigh>RaaCharmFDlow ? RaaCharmFDhigh : RaaCharmFDlow;
795 Double_t RFDlow = RaaCharmFDlow<RaaCharmFDhigh ? RaaCharmFDlow : RaaCharmFDhigh;
796 if ( RFDhigh > FDEhigh ) FDElossMax[ hABbin ] = RFDhigh ;
797 if ( RFDlow < FDEmin ) FDElossMin[ hABbin ] = RFDlow ;
798 if(printout && TMath::Abs(ptprintout-pt)<0.1)
799 cout<<" Hypothesis " << ElossHypo << " sigma-AB "<< sigmaAB <<", Raa FD-max Eloss max "<< FDElossMax[hABbin] <<" Raa FD-min Eloss min "<< FDElossMin[hABbin] <<endl;
808 // Finish filling the y-uncertainties of the Eloss scenarii
809 for (Int_t ibin=0; ibin<=nbins; ibin++){
810 Double_t ipt=0., value =0.;
811 gRAB_ElossHypothesis->GetPoint(ibin,ipt,value);
814 // Uncertainty on Raa due to the Eloss hypothesis
815 Double_t elossYhigh = TMath::Abs( ElossMax[ibin] - value );
816 Double_t elossYlow = TMath::Abs( value - ElossMin[ibin] );
817 gRAB_ElossHypothesis->SetPointEYhigh(ibin, elossYhigh );
818 gRAB_ElossHypothesis->SetPointEYlow(ibin, elossYlow );
819 gRAB_ElossHypothesis->SetPointEXhigh(ibin, 0.2);
820 gRAB_ElossHypothesis->SetPointEXlow(ibin, 0.2);
821 cout << " pt "<< ipt << " Raa "<< value <<" max "<< ElossMax[ibin] << " min " <<ElossMin[ibin] <<endl;
822 cout<<" Eloss syst +"<< elossYhigh <<" - "<< elossYlow <<endl;
823 // cout << " fc max "<< fcElossMax[ibin] << " fc min " <<fcElossMin[ibin] <<endl;
825 // Uncertainty on Raa due to the FD unc & Eloss hypothesis
826 Double_t fdElossEYhigh = TMath::Abs( FDElossMax[ibin] - value );
827 Double_t fdElossEYlow = TMath::Abs( value - FDElossMin[ibin] );
829 Double_t fdEYhigh = gRAB_FeedDownSystematics->GetErrorYhigh(ibin);
830 fdElossEYhigh = TMath::Sqrt( elossYhigh*elossYhigh + fdEYhigh*fdEYhigh );
831 Double_t fdEYlow = gRAB_FeedDownSystematics->GetErrorYlow(ibin);
832 fdElossEYlow = TMath::Sqrt( elossYlow*elossYlow + fdEYlow*fdEYlow );
834 gRAB_FeedDownSystematicsElossHypothesis->SetPointEYhigh(ibin, fdElossEYhigh );
835 gRAB_FeedDownSystematicsElossHypothesis->SetPointEYlow(ibin, fdElossEYlow );
836 gRAB_FeedDownSystematicsElossHypothesis->SetPointEXhigh(ibin, 0.25);
837 gRAB_FeedDownSystematicsElossHypothesis->SetPointEXlow(ibin, 0.25);
838 cout<<" FD & Eloss syst +"<< fdElossEYhigh <<" - "<< fdElossEYlow
839 <<" = + "<< fdElossEYhigh/value <<" - "<< fdElossEYlow/value <<" %" <<endl;
841 // All uncertainty on Raa (FD unc & Eloss + data)
842 Double_t systdatal = gRAB_DataSystematics->GetErrorYlow(ibin);
843 Double_t systdatah = gRAB_DataSystematics->GetErrorYhigh(ibin);
844 Double_t systgbhUnc = TMath::Sqrt( systdatah*systdatah + fdElossEYhigh*fdElossEYhigh );
845 Double_t systgblUnc = TMath::Sqrt( systdatal*systdatal + fdElossEYlow*fdElossEYlow );
846 gRAB_GlobalSystematics->SetPointEYhigh(ibin,systgbhUnc);
847 gRAB_GlobalSystematics->SetPointEYlow(ibin,systgblUnc);
848 cout<<" Data syst +"<< systdatah <<" - "<< systdatal <<" = + "<< systdatah/value <<" - " << systdatal/value << " % "<<endl;
849 cout<<" Global syst +"<< systgbhUnc <<" - "<< systgblUnc << " = + "<< systgbhUnc/value <<" - "<< systgblUnc/value << " %" <<endl;
854 gROOT->SetStyle("Plain");
855 gStyle->SetPalette(1);
856 gStyle->SetOptStat(0);
859 TCanvas *cRABvsRb = new TCanvas("RABvsRb","RAB vs Rb");
860 hRABvsRb->Draw("colz");
863 TCanvas *cRABvsRbvsPt = new TCanvas("cRABvsRbvsPt","RAB vs Rb vs pt");
864 hRABCharmVsRBeautyVsPt->Draw("lego3z");
865 cRABvsRbvsPt->Update();
868 TCanvas *cRABvsRbFDl = new TCanvas("RABvsRbFDl","RAB vs Rb (FD low)");
869 hRABvsRbFDlow->Draw("cont4z");
870 cRABvsRbFDl->Update();
871 TCanvas *cRABvsRbFDh = new TCanvas("RABvsRbFDh","RAB vs Rb (FD high)");
872 hRABvsRbFDhigh->Draw("cont4z");
873 cRABvsRbFDh->Update();
875 TCanvas * cSigmaABptEloss = new TCanvas("cSigmaABptEloss","SigmaAB vs pt, Eloss hypothesis");
876 TH1D *hSigmaABEloss00= new TH1D("hSigmaABEloss00","hSigmaABEloss00",nbins,limits);
877 TH1D *hSigmaABEloss05= new TH1D("hSigmaABEloss05","hSigmaABEloss05",nbins,limits);
878 TH1D *hSigmaABEloss10= new TH1D("hSigmaABEloss10","hSigmaABEloss10",nbins,limits);
879 TH1D *hSigmaABEloss15= new TH1D("hSigmaABEloss15","hSigmaABEloss15",nbins,limits);
880 TH1D *hSigmaABEloss20= new TH1D("hSigmaABEloss20","hSigmaABEloss20",nbins,limits);
881 for (Int_t i=0; i<=nSigmaAB->GetEntriesFast(); i++) {
882 nSigmaAB->GetEntry(i);
884 if( TMath::Abs(Rcb-0.015)<0.009 ) hSigmaABEloss00->Fill(pt,sigmaAB);
885 if( TMath::Abs(Rcb-0.5)<0.009 ) hSigmaABEloss05->Fill(pt,sigmaAB);
886 if( TMath::Abs(Rcb-1.0)<0.009 ) hSigmaABEloss10->Fill(pt,sigmaAB);
887 if( TMath::Abs(Rcb-1.5)<0.009 ) hSigmaABEloss15->Fill(pt,sigmaAB);
888 if( TMath::Abs(Rcb-2.0)<0.009 ) hSigmaABEloss20->Fill(pt,sigmaAB);
890 else if (fdMethod==kNb) {
891 if( TMath::Abs(Rb-0.015)<0.009 ) hSigmaABEloss00->Fill(pt,sigmaAB);
892 if( TMath::Abs(Rb-0.5)<0.009 ) hSigmaABEloss05->Fill(pt,sigmaAB);
893 if( TMath::Abs(Rb-1.0)<0.009 ) hSigmaABEloss10->Fill(pt,sigmaAB);
894 if( TMath::Abs(Rb-1.5)<0.009 ) hSigmaABEloss15->Fill(pt,sigmaAB);
895 if( TMath::Abs(Rb-2.0)<0.009 ) hSigmaABEloss20->Fill(pt,sigmaAB);
898 hSigmaABEloss00->SetLineColor(2);
899 hSigmaABEloss05->SetLineColor(3);
900 hSigmaABEloss10->SetLineColor(4);
901 hSigmaABEloss15->SetLineColor(kMagenta+1);
902 hSigmaABEloss20->SetLineColor(kGreen+2);
903 hSigmaABEloss00->SetMarkerStyle(22);
904 hSigmaABEloss05->SetMarkerStyle(26);
905 hSigmaABEloss10->SetMarkerStyle(20);
906 hSigmaABEloss15->SetMarkerStyle(25);
907 hSigmaABEloss20->SetMarkerStyle(21);
909 hSigmaABEloss05->Draw("ph");
910 hSigmaABEloss10->Draw("phsame");
911 hSigmaABEloss15->Draw("phsame");
912 hSigmaABEloss20->Draw("phsame");
915 hSigmaABEloss20->Draw("p");
916 hSigmaABEloss00->Draw("phsame");
917 hSigmaABEloss05->Draw("phsame");
918 hSigmaABEloss10->Draw("phsame");
919 hSigmaABEloss15->Draw("phsame");
920 hSigmaABEloss20->Draw("phsame");
922 TLegend *legrcb = new TLegend(0.8,0.8,0.95,0.9);
923 legrcb->SetFillColor(0);
924 legrcb->AddEntry(hSigmaABEloss00,"Rc/b=0.0","lp");
925 legrcb->AddEntry(hSigmaABEloss05,"Rc/b=0.5","lp");
926 legrcb->AddEntry(hSigmaABEloss10,"Rc/b=1.0","lp");
927 legrcb->AddEntry(hSigmaABEloss15,"Rc/b=1.5","lp");
928 legrcb->AddEntry(hSigmaABEloss20,"Rc/b=2.0","lp");
930 cSigmaABptEloss->Update();
933 TCanvas * cRABptEloss = new TCanvas("cRABptEloss","RAB vs pt, Eloss hypothesis");
934 hRABEloss00->SetLineColor(2);
935 hRABEloss05->SetLineColor(3);
936 hRABEloss10->SetLineColor(4);
937 hRABEloss15->SetLineColor(kMagenta+1);
938 hRABEloss20->SetLineColor(kGreen+2);
939 hRABEloss00->SetMarkerStyle(22);
940 hRABEloss05->SetMarkerStyle(26);
941 hRABEloss10->SetMarkerStyle(20);
942 hRABEloss15->SetMarkerStyle(25);
943 hRABEloss20->SetMarkerStyle(21);
945 hRABEloss05->Draw("ph");
946 hRABEloss10->Draw("phsame");
947 hRABEloss15->Draw("phsame");
948 hRABEloss20->Draw("phsame");
951 hRABEloss20->Draw("p");
952 hRABEloss00->Draw("phsame");
953 hRABEloss05->Draw("phsame");
954 hRABEloss10->Draw("phsame");
955 hRABEloss15->Draw("phsame");
956 hRABEloss20->Draw("phsame");
958 legrcb = new TLegend(0.8,0.8,0.95,0.9);
959 legrcb->SetFillColor(0);
961 legrcb->AddEntry(hRABEloss00,"Rc/b=0.0","lp");
962 legrcb->AddEntry(hRABEloss05,"Rc/b=0.5","lp");
963 legrcb->AddEntry(hRABEloss10,"Rc/b=1.0","lp");
964 legrcb->AddEntry(hRABEloss15,"Rc/b=0.5","lp");
965 legrcb->AddEntry(hRABEloss20,"Rc/b=2.0","lp");
967 else if (fdMethod==kNb) {
968 legrcb->AddEntry(hRABEloss00,"Rb=0.0","lp");
969 legrcb->AddEntry(hRABEloss05,"Rb=0.5","lp");
970 legrcb->AddEntry(hRABEloss10,"Rb=1.0","lp");
971 legrcb->AddEntry(hRABEloss15,"Rb=0.5","lp");
972 legrcb->AddEntry(hRABEloss20,"Rb=2.0","lp");
975 cRABptEloss->Update();
978 TCanvas * cRABpt = new TCanvas("cRABpt","RAB vs pt, no hypothesis");
979 hRABEloss10->Draw("");
982 TCanvas * cRABptFDUnc = new TCanvas("cRABptFDUnc","RAB vs pt, FD Uncertainties");
983 hRABvsRbFDlow_proj->Draw("");
984 hRABEloss10->Draw("phsame");
985 hRABvsRbFDhigh_proj->SetLineColor(kMagenta+1);
986 hRABvsRbFDhigh_proj->Draw("same");
987 hRABvsRbFDlow_proj->SetLineColor(kGreen+2);
988 hRABvsRbFDlow_proj->Draw("same");
989 legrcb = new TLegend(0.8,0.8,0.95,0.9);
990 legrcb->SetFillColor(0);
991 legrcb->AddEntry(hRABEloss10,"FD Central","lp");
992 legrcb->AddEntry(hRABvsRbFDhigh_proj,"FD Upper unc.","l");
993 legrcb->AddEntry(hRABvsRbFDlow_proj,"FD Lower unc.","l");
995 cRABptFDUnc->Update();
997 TCanvas *RaaPlot = new TCanvas("RaaPlot","RAB vs pt, plot all");
998 RaaPlot->SetTopMargin(0.085);
999 RaaPlot->SetBottomMargin(0.1);
1000 RaaPlot->SetTickx();
1001 RaaPlot->SetTicky();
1002 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);
1003 hRaaCanvas->GetXaxis()->SetTitleSize(0.05);
1004 hRaaCanvas->GetXaxis()->SetTitleOffset(0.9);
1005 hRaaCanvas->GetYaxis()->SetTitleSize(0.05);
1006 hRaaCanvas->GetYaxis()->SetTitleOffset(0.9);
1008 gRAB_Norm->SetFillStyle(1001);
1009 gRAB_Norm->SetFillColor(kGray+2);
1010 gRAB_Norm->Draw("2");
1011 TLine *line = new TLine(0.0172415,1.0,40.,1.0);
1012 line->SetLineStyle(2);
1014 hRABvsPt->SetMarkerColor(kBlue);
1015 hRABvsPt->SetMarkerColor(kBlue);
1016 hRABvsPt->SetMarkerStyle(21);
1017 hRABvsPt->SetMarkerSize(1.1);
1018 hRABvsPt->SetLineWidth(2);
1019 hRABvsPt->Draw("psame");
1020 gRAB_DataSystematics->SetLineColor(kBlue);
1021 gRAB_DataSystematics->SetLineWidth(3);
1022 gRAB_DataSystematics->SetLineWidth(2);
1023 gRAB_DataSystematics->SetFillColor(kRed);
1024 gRAB_DataSystematics->SetFillStyle(0);
1025 gRAB_DataSystematics->Draw("2");
1026 gRAB_FeedDownSystematics->SetFillColor(kViolet+1);
1027 gRAB_FeedDownSystematics->SetFillStyle(1001);
1028 gRAB_FeedDownSystematics->Draw("2");
1029 gRAB_ElossHypothesis->SetLineColor(kMagenta-7);
1030 gRAB_ElossHypothesis->SetFillColor(kMagenta-7);
1031 gRAB_ElossHypothesis->SetFillStyle(1001);
1032 gRAB_ElossHypothesis->Draw("2");
1033 hRABvsPt->Draw("psame");
1034 gRAB_DataSystematics->Draw("2");
1035 legrcb = new TLegend(0.5517241,0.6504237,0.8520115,0.8728814,NULL,"brNDC");
1036 legrcb->SetBorderSize(0);
1037 legrcb->SetTextSize(0.03389831);
1038 legrcb->SetLineColor(1);
1039 legrcb->SetLineStyle(1);
1040 legrcb->SetLineWidth(1);
1041 legrcb->SetFillColor(0);
1042 legrcb->SetFillStyle(1001);
1043 if(cc==k020) legrcb->AddEntry(hRABvsPt,"R_{AA} 0-20% CC","pe");
1044 else if(cc==k4080) legrcb->AddEntry(hRABvsPt,"R_{AA} 40-80% CC","pe");
1045 else legrcb->AddEntry(hRABvsPt,"R_{AA} and stat. unc.","pe");
1046 legrcb->AddEntry(gRAB_DataSystematics,"Syst. from data","f");
1047 legrcb->AddEntry(gRAB_ElossHypothesis,"Syst. from R_{AA}(B)","f");
1048 legrcb->AddEntry(gRAB_FeedDownSystematics,"Syst. from B feed-down","f");
1051 TString system = "Pb-Pb #sqrt{s_{NN}}=2.76 TeV";
1052 if(cc==kpPb0100) system = "p-Pb #sqrt{s_{NN}}=5.023 TeV";
1053 if(decay==1) tc =new TLatex(0.18,0.82,Form("D^{0}, %s ",system.Data()));
1054 else if(decay==2) tc =new TLatex(0.18,0.82,Form("D^{+}, %s ",system.Data()));
1055 else if(decay==3) tc =new TLatex(0.18,0.82,Form("D^{*+}, %s ",system.Data()));
1056 else if(decay==4) tc =new TLatex(0.18,0.82,Form("D_{s}^{+}, %s ",system.Data()));
1057 else tc =new TLatex(0.18,0.82,Form("any (?) D meson, %s ",system.Data()));
1059 tc->SetTextSize(0.038);
1060 tc->SetTextFont(42);
1065 TCanvas *RaaPlotFDEloss = new TCanvas("RaaPlotFDEloss","RAB vs pt, plot FD & ElossUnc");
1066 RaaPlotFDEloss->SetTopMargin(0.085);
1067 RaaPlotFDEloss->SetBottomMargin(0.1);
1070 hRABvsPt->Draw("psame");
1071 gRAB_FeedDownSystematics->SetFillColor(kViolet+1);
1072 gRAB_FeedDownSystematics->SetFillStyle(1001);
1073 gRAB_FeedDownSystematics->Draw("2");
1074 gRAB_ElossHypothesis->SetLineColor(kMagenta-7);
1075 gRAB_ElossHypothesis->SetFillColor(kMagenta-7);
1076 gRAB_ElossHypothesis->SetFillStyle(1001);
1077 gRAB_ElossHypothesis->Draw("2");
1078 gRAB_FeedDownSystematicsElossHypothesis->SetLineColor(kBlack);
1079 gRAB_FeedDownSystematicsElossHypothesis->SetFillStyle(0);
1080 gRAB_FeedDownSystematicsElossHypothesis->SetFillColor(kViolet+1);
1081 gRAB_FeedDownSystematicsElossHypothesis->Draw("2");
1082 hRABvsPt->Draw("psame");
1083 legrcb = new TLegend(0.6,0.6,0.9,0.9);
1084 legrcb->SetBorderSize(0);
1085 legrcb->SetTextSize(0.03389831);
1086 legrcb->SetLineColor(1);
1087 legrcb->SetLineStyle(1);
1088 legrcb->SetLineWidth(1);
1089 legrcb->SetFillColor(0);
1090 legrcb->SetFillStyle(1001);
1091 legrcb->AddEntry(hRABvsPt,"R_{PbPb} and stat. unc.","pe");
1092 legrcb->AddEntry(gRAB_ElossHypothesis,"Energy loss syst.","f");
1093 legrcb->AddEntry(gRAB_FeedDownSystematics,"Feed down syst.","f");
1094 legrcb->AddEntry(gRAB_FeedDownSystematicsElossHypothesis,"Feed down & Eloss syst.","f");
1096 RaaPlotFDEloss->Update();
1099 TCanvas *RaaPlotGlob = new TCanvas("RaaPlotGlob","RAB vs pt, plot Global unc");
1100 RaaPlotGlob->SetTopMargin(0.085);
1101 RaaPlotGlob->SetBottomMargin(0.1);
1102 RaaPlotGlob->SetTickx();
1103 RaaPlotGlob->SetTicky();
1106 hRABvsPt->Draw("psame");
1107 gRAB_DataSystematics->Draw("2");
1108 gRAB_FeedDownSystematicsElossHypothesis->Draw("2");
1109 gRAB_GlobalSystematics->SetLineColor(kRed);
1110 gRAB_GlobalSystematics->SetLineWidth(2);
1111 gRAB_GlobalSystematics->SetFillColor(kRed);
1112 gRAB_GlobalSystematics->SetFillStyle(3002);
1113 gRAB_GlobalSystematics->Draw("2");
1114 hRABvsPt->Draw("psame");
1115 legrcb = new TLegend(0.6,0.6,0.9,0.9);
1116 legrcb->SetBorderSize(0);
1117 legrcb->SetTextSize(0.03389831);
1118 legrcb->SetLineColor(1);
1119 legrcb->SetLineStyle(1);
1120 legrcb->SetLineWidth(1);
1121 legrcb->SetFillColor(0);
1122 legrcb->SetFillStyle(1001);
1123 legrcb->AddEntry(hRABvsPt,"R_{PbPb} and stat. unc.","pe");
1124 legrcb->AddEntry(gRAB_DataSystematics,"Data syst.","f");
1125 legrcb->AddEntry(gRAB_FeedDownSystematicsElossHypothesis,"Feed down & Eloss syst.","f");
1126 legrcb->AddEntry(gRAB_GlobalSystematics,"Global syst.","f");
1128 RaaPlotGlob->Update();
1132 TCanvas *RaaPlotSimple = new TCanvas("RaaPlotSimple","RAB vs pt, plot Simple unc");
1133 RaaPlotSimple->SetTopMargin(0.085);
1134 RaaPlotSimple->SetBottomMargin(0.1);
1135 RaaPlotSimple->SetTickx();
1136 RaaPlotSimple->SetTicky();
1139 hRABvsPt->Draw("psame");
1140 gRAB_GlobalSystematics->SetLineColor(kBlue);
1141 gRAB_GlobalSystematics->SetLineWidth(2);
1142 gRAB_GlobalSystematics->SetFillStyle(0);
1143 gRAB_GlobalSystematics->Draw("2");
1144 gRAB_Norm->Draw("2");
1145 hRABvsPt->Draw("psame");
1146 legrcb = new TLegend(0.5991379,0.6949153,0.8534483,0.8559322,NULL,"brNDC");
1147 legrcb->SetBorderSize(0);
1148 legrcb->SetTextSize(0.03389831);
1149 legrcb->SetLineColor(1);
1150 legrcb->SetLineStyle(1);
1151 legrcb->SetLineWidth(1);
1152 legrcb->SetFillColor(0);
1153 legrcb->SetFillStyle(1001);
1154 if(cc==k020) legrcb->AddEntry(hRABvsPt,"R_{AA} 0-20% CC","pe");
1155 else if(cc==k4080) legrcb->AddEntry(hRABvsPt,"R_{AA} 40-80% CC","pe");
1156 else legrcb->AddEntry(hRABvsPt,"R_{AA} and stat. unc.","pe");
1157 legrcb->AddEntry(gRAB_GlobalSystematics,"Systematics","f");
1160 RaaPlotSimple->Update();
1163 TCanvas *c = new TCanvas("c","");
1164 systematicsAB->DrawErrors();
1167 TCanvas *cStatUnc = new TCanvas("cStatUnc","stat unc");
1168 cStatUnc->Divide(2,2);
1170 fhStatUncEffcSigmaAB_Raa->Draw("e");
1172 fhStatUncEffbSigmaAB_Raa->Draw("e");
1174 fhStatUncEffcFDAB_Raa->Draw("e");
1176 fhStatUncEffbFDAB_Raa->Draw("e");
1180 // Write the information to a file
1182 TFile * out = new TFile(outfile,"recreate");
1187 hRABCharmVsRBeautyVsPt->Write();
1188 for(Int_t j=0; j<=nbins; j++) hRCharmVsRBeauty[j]->Write();
1189 // for(Int_t j=0; j<=nbins; j++) hRCharmVsElossHypo[j]->Write();
1191 hRABvsPt_DataSystematics->Write();
1192 gRAB_ElossHypothesis->Write();
1193 gRAB_FeedDownSystematics->Write();
1194 gRAB_fcFeedDownOnly->Write();
1195 gRAB_DataSystematics->Write();
1196 gRAB_DataSystematicsPP->Write();
1197 gSigmaPPSystTheory->Write();
1198 gRAB_DataSystematicsAB->Write();
1200 gRAB_FeedDownSystematicsElossHypothesis->Write();
1201 gRAB_GlobalSystematics->Write();
1202 if(isScaledAndExtrapRef) hCombinedReferenceFlag->Write();
1208 //____________________________________________________________
1209 Bool_t PbPbDataSyst(AliHFSystErr *syst, Double_t pt, Int_t cc, Double_t &dataSystUp, Double_t &dataSystDown)
1212 Double_t err=0., errUp=1., errDown=1.;
1216 err = syst->GetTotalSystErr(pt)*syst->GetTotalSystErr(pt);
1220 // Apply an asymetric PID uncertainty for 2010 PbPb data only
1221 if( syst->GetRunNumber()==10 && syst->GetCollisionType()==1 ) {
1222 if( cc==k07half || cc==k020 || cc==k010 || cc==k1020 || cc==k2040 ) {
1223 if(pt<6) pidunc = 0.15;
1225 errUp = err + pidunc*pidunc - syst->GetPIDEffErr(pt)*syst->GetPIDEffErr(pt);
1228 else if ( cc==k3050 || cc==k4080 || cc==k4060 || cc==k6080 ){
1229 if(pt<3.1) pidunc = 0.10;
1231 errUp = err + pidunc*pidunc - syst->GetPIDEffErr(pt)*syst->GetPIDEffErr(pt);
1235 else { isOk = true; }
1237 dataSystUp = TMath::Sqrt(errUp);
1238 dataSystDown = TMath::Sqrt(errDown);