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