AliPhysics  251aa1e (251aa1e)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AverageDmesonRaa.C
Go to the documentation of this file.
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 //------------------------------------------------------------------------------------------
24 //
25 // MACRO TO COMPUTE THE AVERAGE D0, D+ and D*+ meson RAA
26 //
27 // Different options considered as weights (relative stat, relative stat+uncorr-syst,
28 // relative stat+global-syst, absolute stat) as option of the macro
29 //
30 // Output: merged RAA root file, plus a tex file with the results
31 //
32 // Usage:
33 // 1. First set the average pt bining (variables nbins and ptbinlimits on top of the macro)
34 // 2. Load the libraries and compile the macro
35 // 3. Macro parameters:
36 // a: D0 RAA file, b: D0 pp reference file
37 // d: D+ RAA file, d: D+ pp reference file
38 // e: D*+ RAA file, f: D*+ pp reference file
39 // g: output filename
40 // h: average option (weights)
41 // i: centrality range, j: centrality estimator
42 // k: flag in case the pp reference has the split separated uncertainties
43 // l,m,n: flag in case D0/D+/D*+ pp reference has some extrapolated bin
44 //
45 //
46 // Author: Z. Conesa del Valle
47 //
48 //------------------------------------------------------------------------------------------
49 
50 
52 // Compilation instructions
53 //
54 // gSystem->SetIncludePath("-I. -I$ROOTSYS/include -I$ALICE_ROOT -I$ALICE_ROOT/include -I$ALICE_PHYSICS -I$ALICE_PHYSICS/include -I$ALICE_ROOT/ITS -I$ALICE_ROOT/TPC -I$ALICE_ROOT/CONTAINERS -I$ALICE_ROOT/STEER/STEER -I$ALICE_ROOT/STEER/STEERBase -I$ALICE_ROOT/STEER/ESD -I$ALICE_ROOT/STEER/AOD -I$ALICE_ROOT/TRD -I$ALICE_ROOT/macros -I$ALICE_ROOT/ANALYSIS -I$ALICE_PHYSICS/OADB -I$ALICE_PHYSICS/PWGHF -I$ALICE_PHYSICS/PWGHF/base -I$ALICE_PHYSICS/PWGHF/vertexingHF -I$ALICE_PHYSICS/PWG/FLOW/Base -I$ALICE_PHYSICS/PWG/FLOW/Tasks -I$ALICE_PHYSICS/PWG -g");
55 // .x $ALICE_PHYSICS/PWGHF/vertexingHF/macros/LoadLibraries.C
56 // .L $ALICE_PHYSICS/PWGHF/vertexingHF/macros/AverageDmesonRaa.C+
57 //
59 
60 //const Int_t nbins=10;
61 //Double_t ptbinlimits[nbins+1]={1.,2.,3.,4.,5.,6.,7.,8.,12.,16.,24.};
62 //Double_t ptbinlimits[nbins+1]={1.,2.,3.,4.,5.,6.,8.,12.,16.,24.,36.};
63 const Int_t nbins=7;
64 Double_t ptbinlimits[nbins+1]={1.,2.,4.,6.,8.,12.,16.,24.};
66 
67 //const Int_t ptmaxPPRefData[3] = { 16., 24., 24. };
68 
70 
73 
74 //____________________________________________________________
76 {
77  Int_t istart =0;
78  Int_t npoints = gr->GetN();
79  for(Int_t i=0; i<=npoints; i++){
80  Double_t x=0.,y=0.;
81  gr->GetPoint(i,x,y);
82  if ( TMath::Abs ( x - pt ) < 0.4 ) {
83  istart = i;
84  break;
85  }
86  }
87  return istart;
88 }
89 //____________________________________________________________
91 {
92  Int_t istart =-1;
93  Int_t npoints = gr->GetN();
94  Double_t x=0.,y=0.;
95  for(Int_t i=0; i<=npoints; i++){
96  gr->GetPoint(i,x,y);
97  if ( TMath::Abs ( x - pt ) < 0.4 ) {
98  istart = i;
99  break;
100  }
101  }
102 
103  uncLow=0.; uncHigh=0.;
104 
105  if(istart>-1){
106  uncLow = gr->GetErrorYlow(istart)/y;
107  uncHigh = gr->GetErrorYhigh(istart)/y;
108  // cout<<" pt="<<pt<<" -"<<uncLow<<" +"<<uncHigh<<endl;
109  }
110 }
111 
112 //____________________________________________________________
113 Double_t GetWeight(Int_t averageoption, Double_t pt,
114  TH1D* hRaa, Double_t raaSystLow, Double_t raaSystHigh,
115  Double_t ppSystRawYieldCutVar, Double_t ppSystRawYieldCutVarPid,
116  Double_t ABSystRawYieldCutVar, Double_t ABSystRawYieldCutVarPid)
117 {
118  Double_t weight=1.0;
119  Int_t hbin = hRaa->FindBin(pt);
120 
121  Double_t stat = hRaa->GetBinError(hbin);
122  Double_t relativeStat = stat / hRaa->GetBinContent(hbin);
123  Double_t weightStat=0.;
124  Double_t relativeSyst = 0.;
125  if(averageoption==kRelativeStatUnc) {
126  weightStat = relativeStat;
127  }
128  else if(averageoption==kAbsoluteStatUnc) {
129  weightStat = stat;
130  }
131  else if(averageoption==kRelativeStatUncorrWoPidSyst) {
132  relativeSyst = TMath::Sqrt( ppSystRawYieldCutVar*ppSystRawYieldCutVar + ABSystRawYieldCutVar*ABSystRawYieldCutVar );
133  }
134  else if(averageoption==kRelativeStatUncorrWPidSyst) {
135  relativeSyst = TMath::Sqrt( ppSystRawYieldCutVarPid*ppSystRawYieldCutVarPid + ABSystRawYieldCutVar*ABSystRawYieldCutVarPid );
136  }
137  else if(averageoption==kRelativeStatGlobalSyst) {
138  relativeSyst = raaSystHigh>raaSystLow ? raaSystHigh : raaSystLow;
139  }
140 
141  // weight = TMath::Sqrt( relativeStat*relativeStat + relativeSyst*relativeSyst );
142  weight = TMath::Sqrt( weightStat*weightStat + relativeSyst*relativeSyst );
143  // cout<< endl<<" rel stat "<< relativeStat<<" rel syst "<< relativeSyst<<" weight="<<(1.0/(weight*weight))<<endl;
144 
145  return (1.0/(weight*weight));
146 }
147 
148 
149 //____________________________________________________________
150 void AverageDmesonRaa( const char* fD0Raa="", const char* fD0ppRef="",
151  const char* fDplusRaa="", const char* fDplusppRef="",
152  const char* fDstarRaa="", const char* fDstarppRef="",
153  const char* outfile="", Int_t averageOption=kRelativeStatUnc, Int_t cc=kpPb0100, Int_t ccestimator=kV0M,
154  Bool_t isReadAllPPUnc=false, Bool_t isPPRefExtrapD0=false, Bool_t isPPRefExtrapDplus=false, Bool_t isPPRefExtrapDstar=false)
155 {
156 
157  FILE *resultFile;
158  TString foutname = "Average";
159  if(fD0Raa) foutname += "Dzero";
160  if(fDplusRaa) foutname += "Dplus";
161  if(fDstarRaa) foutname += "Dstar";
162  if(averageOption==kRelativeStatUnc) foutname+= "_RelStatUncWeight";
163  else if(averageOption==kAbsoluteStatUnc) foutname+= "_AbsStatUncWeight";
164  else if(averageOption==kRelativeStatUncorrWoPidSyst) foutname+= "_RelStatUncorrWeight";
165  TDatime d;
166  TString ndate = Form("%02d%02d%04d",d.GetDay(),d.GetMonth(),d.GetYear());
167  resultFile = fopen( Form("%s_result_%s.txt",foutname.Data(),ndate.Data()),"w");
168  fprintf(resultFile,"Ptmin (GeV/c) Ptmax (GeV/c) Raa(Daverage) +-(stat) +(syst) - (syst) \n\n");
169 
170 
171  // Define here all needed histograms/graphs to be retrieved
172  TH1D *hDmesonRaa[3];
173  TH1I *hCombinedReferenceFlag[3];
174  TGraphAsymmErrors *gDataSystematicsPP[3], *gDataSystematicsAB[3];
175  TGraphAsymmErrors *gScalingUncPP[3];
176  TGraphAsymmErrors *gRABFeedDownSystematicsElossHypothesis[3];
177  TGraphAsymmErrors *gRAB_GlobalSystematics[3];
178  TH1D *hDmesonPPRef[3], *hDmesonPPYieldExtrUnc[3], *hDmesonPPCutVarUnc[3], *hDmesonPPIDPUnc[3], *hDmesonPPMCPtUnc[3];
179 
180  // Define here all output histograms/graphs
181  TH1D *hDmesonAverageRAB = new TH1D("hDmesonAverageRAB","D meson average Raa ; p_{T} (GeV/c)",nbins,ptbinlimits);
182  TGraphAsymmErrors *gRABNorm;
183  TGraphAsymmErrors *gRAB_DmesonAverage_GlobalSystematics = new TGraphAsymmErrors(0);
184  TGraphAsymmErrors *gRAB_DmesonAverage_FeedDownSystematicsElossHypothesis = new TGraphAsymmErrors(0);
185  TGraphAsymmErrors *gRAB_DmesonAverage_ScalingSystematicsPP = new TGraphAsymmErrors(0);
186  TGraphAsymmErrors *gRAB_DmesonAverage_DataSystematicsPP = new TGraphAsymmErrors(0);
187  TGraphAsymmErrors *gRAB_DmesonAverage_DataSystematicsAB = new TGraphAsymmErrors(0);
188  TGraphAsymmErrors *gRAB_DmesonAverage_TrackingSystematicsPP = new TGraphAsymmErrors(0);
189  TGraphAsymmErrors *gRAB_DmesonAverage_TrackingSystematicsAB = new TGraphAsymmErrors(0);
190  gRAB_DmesonAverage_GlobalSystematics->SetNameTitle("gRAB_DmesonAverage_GlobalSystematics","DmesonAverage GlobalSystematics");
191  gRAB_DmesonAverage_FeedDownSystematicsElossHypothesis->SetNameTitle("gRAB_DmesonAverage_FeedDownSystematicsElossHypothesis","DmesonAverage FeedDownSystematicsElossHypothesis");
192  gRAB_DmesonAverage_ScalingSystematicsPP->SetNameTitle("gRAB_DmesonAverage_ScalingSystematicsPP","DmesonAverage Scaling uncertainty PP");
193  gRAB_DmesonAverage_DataSystematicsPP->SetNameTitle("gRAB_DmesonAverage_DataSystematicsPP","DmesonRaaAverage DataSystematicsPP - tracking uncertainty PP");
194  gRAB_DmesonAverage_DataSystematicsAB->SetNameTitle("gRAB_DmesonAverage_DataSystematicsAB","DmesonRaaAverage DataSystematicsAB - tracking uncertainty AB");
195  gRAB_DmesonAverage_TrackingSystematicsPP->SetNameTitle("gRAB_DmesonAverage_TrackingSystematicsPP","DmesonRaaAverage tracking uncertainty PP");
196  gRAB_DmesonAverage_TrackingSystematicsAB->SetNameTitle("gRAB_DmesonAverage_TrackingSystematicsAB","DmesonRaaAverage tracking uncertainty AB");
197 
198 
199  const char *filenamesRaa[3] = { fD0Raa, fDplusRaa, fDstarRaa };
200  const char *filenamesReference[3] = { fD0ppRef, fDplusppRef, fDstarppRef };
201  const char *filenamesSuffixes[3] = { "Dzero", "Dplus", "Dstar" };
202  Bool_t isDmeson[3] = { true, true, true };
203 
204  if(strcmp(filenamesRaa[0],"")==0) { cout<<" Dzero not set, error"<<endl; return; }
205 
206  //
207  // Get Raa file histos and graphs
208  //
209  for(Int_t j=0; j<3; j++) {
210  if(strcmp(filenamesRaa[j],"")==0) { isDmeson[j]=false; continue; }
211  cout<<" Reading file "<<filenamesRaa[j]<<"..."<<endl;
212  TFile *fDRaa = TFile::Open(filenamesRaa[j],"read");
213  if(!fDRaa){ cout<<" Error on file !!!!!"<<filenamesRaa[j]<<endl; return; }
214  hDmesonRaa[j] = (TH1D*)fDRaa->Get("hRABvsPt");
215  hDmesonRaa[j]->SetName(Form("%s_%s",hDmesonRaa[j]->GetName(),filenamesSuffixes[j]));
216 
217  // cout<< hDmesonRaa[j]<< " bins="<< hDmesonRaa[j]->GetNbinsX()<<", value bin 1 ="<<hDmesonRaa[j]->GetBinContent(1)<<endl;
218 
219  if(j==0) gRABNorm = (TGraphAsymmErrors*)fDRaa->Get("gRAB_Norm");
220  //
221  gDataSystematicsPP[j] = (TGraphAsymmErrors*)fDRaa->Get("gRAB_DataSystematicsPP");
222  gDataSystematicsPP[j]->SetName(Form("%s_%s",gDataSystematicsPP[j]->GetName(),filenamesSuffixes[j]));
223  gDataSystematicsAB[j] = (TGraphAsymmErrors*)fDRaa->Get("gRAB_DataSystematicsAB");
224  gDataSystematicsAB[j]->SetName(Form("%s_%s",gDataSystematicsAB[j]->GetName(),filenamesSuffixes[j]));
225  gRABFeedDownSystematicsElossHypothesis[j] = (TGraphAsymmErrors*)fDRaa->Get("gRAB_FeedDownSystematicsElossHypothesis");
226  gRABFeedDownSystematicsElossHypothesis[j]->SetName(Form("%s_%s",gRABFeedDownSystematicsElossHypothesis[j]->GetName(),filenamesSuffixes[j]));
227  gRAB_GlobalSystematics[j] = (TGraphAsymmErrors*)fDRaa->Get("gRAB_GlobalSystematics");
228  gRAB_GlobalSystematics[j]->SetName(Form("%s_%s",gRAB_GlobalSystematics[j]->GetName(),filenamesSuffixes[j]));
229  }
230  //
231  // Get pp-reference file histos and graphs
232  //
233  const char *pprefhgnames[6] = { "fhScaledData","gScaledDataSystExtrap",
234  "fhScaledSystRebinYieldExtraction","fhScaledSystRebinCutVariation","fhScaledSystRebinPIDUnc","fhScaledSystRebinMCPt"};
235  Bool_t isPPRefExtrap[3] = { isPPRefExtrapD0, isPPRefExtrapDplus, isPPRefExtrapDstar };
236  for(Int_t j=0; j<3; j++) {
237  if(strcmp(filenamesReference[j],"")==0) { isDmeson[j]=false; continue; }
238  cout<<" Reading file "<<filenamesReference[j]<<"..."<<endl;
239  TFile *fRef = TFile::Open(filenamesReference[j],"read");
240  gScalingUncPP[j] = (TGraphAsymmErrors*)fRef->Get(pprefhgnames[1]);
241  gScalingUncPP[j]->SetName(Form("%s_%s",gScalingUncPP[j]->GetName(),filenamesSuffixes[j]));
242  if(isPPRefExtrap[j]) {
243  hCombinedReferenceFlag[j] = (TH1I*)fRef->Get("hCombinedReferenceFlag");
244  hCombinedReferenceFlag[j]->SetName(Form("%s_%s",hCombinedReferenceFlag[j]->GetName(),filenamesSuffixes[j]));
245  }
246  if(isReadAllPPUnc){
247  const char*hname="fhScaledData";
248  if(isPPRefExtrap[j]) hname="hReference";
249  hDmesonPPRef[j] = (TH1D*)fRef->Get(hname);
250  hDmesonPPRef[j]->SetName(Form("%s_%s",hDmesonPPRef[j]->GetName(),filenamesSuffixes[j]));
251  hDmesonPPYieldExtrUnc[j] = (TH1D*)fRef->Get(pprefhgnames[2]);
252  hDmesonPPYieldExtrUnc[j]->SetName(Form("%s_%s",hDmesonPPYieldExtrUnc[j]->GetName(),filenamesSuffixes[j]));
253  hDmesonPPCutVarUnc[j] = (TH1D*)fRef->Get(pprefhgnames[3]);
254  hDmesonPPCutVarUnc[j]->SetName(Form("%s_%s",hDmesonPPCutVarUnc[j]->GetName(),filenamesSuffixes[j]));
255  hDmesonPPIDPUnc[j] = (TH1D*)fRef->Get(pprefhgnames[4]);
256  hDmesonPPIDPUnc[j]->SetName(Form("%s_%s",hDmesonPPIDPUnc[j]->GetName(),filenamesSuffixes[j]));
257  hDmesonPPMCPtUnc[j] = (TH1D*)fRef->Get(pprefhgnames[5]);
258  hDmesonPPMCPtUnc[j]->SetName(Form("%s_%s",hDmesonPPMCPtUnc[j]->GetName(),filenamesSuffixes[j]));
259  }
260  }
261 
262 
263  cout<<" Getting instances of AliHFSystErr... "<<endl;
264  AliHFSystErr *ppSyst[3];
265  Double_t ppTracking[3][nbins], ppSystRawYield[3][nbins], ppSystCutVar[3][nbins], ppSystPid[3][nbins];
266  ppSyst[0] = new AliHFSystErr("ppSyst_Dzero","ppSyst_Dzero");
267  ppSyst[1] = new AliHFSystErr("ppSyst_Dplus","ppSyst_Dplus");
268  ppSyst[2] = new AliHFSystErr("ppSyst_Dstar","ppSyst_Dstar");
269  for(Int_t j=0; j<3; j++) {
270  if(!isDmeson[j]) { delete ppSyst[j]; continue; }
271  ppSyst[j]->Init(j+1);
272  for(Int_t ipt=0; ipt<nbins; ipt++) {
273  Double_t ptval = ptbinlimits[ipt] + (ptbinlimits[ipt+1]-ptbinlimits[ipt])/2.;
274  ppTracking[j][ipt]=0.; ppSystRawYield[j][ipt]=0; ppSystCutVar[j][ipt]=0; ppSystPid[j][ipt]=0.;
275  ppTracking[j][ipt] =ppSyst[j]->GetTrackingEffErr(ptval);
276  ppSystRawYield[j][ipt]=ppSyst[j]->GetRawYieldErr(ptval);
277  ppSystCutVar[j][ipt] =ppSyst[j]->GetCutsEffErr(ptval);
278  ppSystPid[j][ipt] =ppSyst[j]->GetPIDEffErr(ptval);
279  }
280  delete ppSyst[j];
281  }
282 
283  AliHFSystErr *ABSyst[3];
284  Double_t ABTracking[3][nbins], ABSystRawYield[3][nbins], ABSystCutVar[3][nbins], ABSystPid[3][nbins];
285  ABSyst[0] = new AliHFSystErr("ABSyst_Dzero","ABSyst_Dzero");
286  ABSyst[1] = new AliHFSystErr("ABSyst_Dplus","ABSyst_Dplus");
287  ABSyst[2] = new AliHFSystErr("ABSyst_Dstar","ABSyst_Dstar");
288  //
289  for(Int_t j=0; j<3; j++) {
290  // PbPb systematics
291  ABSyst[j]->SetCollisionType(1); // PbPb by default
292  if ( cc == k010 ) ABSyst[j]->SetCentrality("010");
293  else if ( cc == k1020 ) ABSyst[j]->SetCentrality("1020");
294  else if ( cc == k2040 || cc == k2030 || cc == k3040 ) {
295  ABSyst[j]->SetCentrality("2040");
296  ABSyst[j]->SetIsPbPb2010EnergyScan(true);
297  }
298  else if ( cc == k3050 ) ABSyst[j]->SetCentrality("3050");
299  else if ( cc == k4060 || cc == k4050 || cc == k5060 ) ABSyst[j]->SetCentrality("4060");
300  else if ( cc == k6080 || cc == k5080 ) ABSyst[j]->SetCentrality("6080");
301  else if ( cc == k4080 ) ABSyst[j]->SetCentrality("4080");
302  // Going to pPb systematics
303  else if ( cc == kpPb0100 || cc == kpPb020 || cc == kpPb2040 || cc == kpPb4060 || cc == kpPb60100 ) {
304  ABSyst[j]->SetCollisionType(2);
305  if(ccestimator==kV0A) {
306  if(cc == kpPb020) ABSyst[j]->SetCentrality("020V0A");
307  else if(cc == kpPb2040) ABSyst[j]->SetCentrality("2040V0A");
308  else if(cc == kpPb4060) ABSyst[j]->SetCentrality("4060V0A");
309  else if(cc == kpPb60100) ABSyst[j]->SetCentrality("60100V0A");
310  } else if (ccestimator==kZNA) {
311  if(cc == kpPb020) ABSyst[j]->SetCentrality("020ZNA");
312  else if(cc == kpPb2040) ABSyst[j]->SetCentrality("2040ZNA");
313  else if(cc == kpPb4060) ABSyst[j]->SetCentrality("4060ZNA");
314  else if(cc == kpPb60100) ABSyst[j]->SetCentrality("60100ZNA");
315  } else if (ccestimator==kCL1) {
316  if(cc == kpPb020) ABSyst[j]->SetCentrality("020CL1");
317  else if(cc == kpPb2040) ABSyst[j]->SetCentrality("2040CL1");
318  else if(cc == kpPb4060) ABSyst[j]->SetCentrality("4060CL1");
319  else if(cc == kpPb60100) ABSyst[j]->SetCentrality("60100CL1");
320  } else {
321  if(!(cc == kpPb0100)) {
322  cout <<" Error on the pPb options"<<endl;
323  return;
324  }
325  }
326  }
327  //
328  }
329  for(Int_t j=0; j<3; j++) {
330  if(!isDmeson[j]) { delete ABSyst[j]; continue; }
331  ABSyst[j]->Init(j+1);
332  for(Int_t ipt=0; ipt<nbins; ipt++) {
333  Double_t ptval = ptbinlimits[ipt] + (ptbinlimits[ipt+1]-ptbinlimits[ipt])/2.;
334  ABTracking[j][ipt]=0.; ABSystRawYield[j][ipt]=0; ABSystCutVar[j][ipt]=0; ABSystPid[j][ipt]=0.;
335  ABTracking[j][ipt] =ABSyst[j]->GetTrackingEffErr(ptval);
336  ABSystRawYield[j][ipt]=ABSyst[j]->GetRawYieldErr(ptval);
337  ABSystCutVar[j][ipt] =ABSyst[j]->GetCutsEffErr(ptval);
338  ABSystPid[j][ipt] =ABSyst[j]->GetPIDEffErr(ptval);
339  }
340  delete ABSyst[j];
341  }
342 
343  //
344  // Loop per pt bin
345  //
346  for(Int_t ipt=0; ipt<nbins; ipt++) {
347 
348  cout<<" Calculation for pt bin ("<<ptbinlimits[ipt]<<","<<ptbinlimits[ipt+1]<<")"<<endl;
349 
350  Double_t ptval = ptbinlimits[ipt] + (ptbinlimits[ipt+1]-ptbinlimits[ipt])/2.;
351  Double_t RaaDmeson[3]={0.,0.,0.};
352  Double_t RaaDmesonStat[3]={0.,0.,0.};
353  Double_t RaaDmesonSystLow[3]={0.,0.,0.};
354  Double_t RaaDmesonSystHigh[3]={0.,0.,0.};
355  Double_t weight[3]={0.,0.,0.};
356  Double_t ppSystLow[3]={0.,0.,0.};
357  Double_t ppSystHigh[3]={0.,0.,0.};
358  Double_t ppSystUncorrLow[3]={0.,0.,0.};
359  Double_t ppSystUncorrHigh[3]={0.,0.,0.};
360  // Double_t ppTracking[3]={0.,0.,0.};
361  Double_t ScalingLow[3]={0.,0.,0.};
362  Double_t ScalingHigh[3]={0.,0.,0.};
363  Double_t ppSystRawYieldCutVar[3]={0.,0.,0.};
364  Double_t ppSystRawYieldCutVarPid[3]={0.,0.,0.};
365  Double_t ABSystLow[3]={0.,0.,0.};
366  Double_t ABSystHigh[3]={0.,0.,0.};
367  Double_t ABSystUncorrLow[3]={0.,0.,0.};
368  Double_t ABSystUncorrHigh[3]={0.,0.,0.};
369  Double_t ABSystRawYieldCutVar[3]={0.,0.,0.};
370  Double_t ABSystRawYieldCutVarPid[3]={0.,0.,0.};
371  // Double_t ABTracking[3]={0.,0.,0.};
372  Double_t RabFdElossLow[3]={0.,0.,0.};
373  Double_t RabFdElossHigh[3]={0.,0.,0.};
374  Double_t RabGlobalLow[3]={0.,0.,0.};
375  Double_t RabGlobalHigh[3]={0.,0.,0.};
376 
377  Double_t average=0., averageStat=0.;
378  Double_t weightTot=0.;
379  Double_t ppTrackingAv=0., ABTrackingAv=0.;
380  Double_t ppDataSystAvLow=0., ppDataSystAvHigh=0.;
381  Double_t ABDataSystAvLow=0., ABDataSystAvHigh=0.;
382  Double_t scalingLowAv=0., scalingHighAv=0.;
383  Double_t raaSystUncorrLow=0., raaSystUncorrHigh=0.;
384  Double_t raabeautyLow=0., raabeautyHigh=0.;
385 
386  Int_t histoBin=-1;
387 
388  // Get tracking uncertainties and raw yield and cut-variation and pid-systematics
389  if(isDebug) cout<<" Retrieving tracking + rawyield systematics"<<endl;
390  for(Int_t j=0; j<3; j++) {
391  if(!isDmeson[j]) continue;
392  ppSystRawYieldCutVar[j] = TMath::Sqrt( ppSystRawYield[j][ipt]*ppSystRawYield[j][ipt]
393  + ppSystCutVar[j][ipt]*ppSystCutVar[j][ipt] );
394  ppSystRawYieldCutVarPid[j] = TMath::Sqrt( ppSystRawYield[j][ipt]*ppSystRawYield[j][ipt]
395  + ppSystCutVar[j][ipt]*ppSystCutVar[j][ipt]
396  + ppSystPid[j][ipt]*ppSystPid[j][ipt] );
397  ABSystRawYieldCutVar[j] = TMath::Sqrt( ABSystRawYield[j][ipt]*ABSystRawYield[j][ipt]
398  + ABSystCutVar[j][ipt]*ABSystCutVar[j][ipt] );
399  ABSystRawYieldCutVarPid[j] = TMath::Sqrt( ABSystRawYield[j][ipt]*ABSystRawYield[j][ipt]
400  + ABSystCutVar[j][ipt]*ABSystCutVar[j][ipt]
401  + ABSystPid[j][ipt]*ABSystPid[j][ipt] );
402  if(isDebug) cout<<" j="<<j<<" pt="<< ptval<<" ppref unc RY+CV="<<ppSystRawYieldCutVar[j]<<" RY+CV+PID="<<ppSystRawYieldCutVarPid[j]<<endl;
403  if(isDebug) cout<<" j="<<j<<" pt="<< ptval<<" AB unc RY+CV="<<ABSystRawYieldCutVar[j]<<" RY+CV+PID="<<ABSystRawYieldCutVarPid[j]<<endl;
404  }
405 
406  if(isReadAllPPUnc){
407  if(isDebug) cout<<" Retrieving all pp reference systematics from the rebinned file"<<endl;
408  for(Int_t j=0; j<3; j++) {
409  if(!isDmeson[j]) continue;
410  Int_t ibin = hDmesonPPRef[j]->FindBin(ptval);
411  Double_t ppval = hDmesonPPRef[j]->GetBinContent(ibin);
412  Double_t rawyield = hDmesonPPYieldExtrUnc[j]->GetBinContent( hDmesonPPYieldExtrUnc[j]->FindBin(ptval) )/ppval;
413  Double_t cutvar = hDmesonPPCutVarUnc[j]->GetBinContent( hDmesonPPCutVarUnc[j]->FindBin(ptval) )/ppval;
414  Double_t pid = hDmesonPPIDPUnc[j]->GetBinContent( hDmesonPPIDPUnc[j]->FindBin(ptval) )/ppval;
415  ppSystRawYieldCutVar[j] = TMath::Sqrt( rawyield*rawyield + cutvar*cutvar );
416  ppSystRawYieldCutVarPid[j] = TMath::Sqrt( rawyield*rawyield + cutvar*cutvar + pid*pid );
417  if(isDebug) cout<<"redo j="<<j<<" pt="<< ptval<<" ppref unc RY+CV="<<ppSystRawYieldCutVar[j]<<" RY+CV+PID="<<ppSystRawYieldCutVarPid[j]<<endl;
418  }
419  }
420  // Check for the pp reference systematics for extrapolated pt bins
421  for(Int_t j=0; j<3; j++) {
422  if(isPPRefExtrap[j]){
423  // if(ptval>ptmaxPPRefData[j]) {
424  Int_t ippbin = hCombinedReferenceFlag[j]->FindBin(ptval);
425  Bool_t flag = hCombinedReferenceFlag[j]->GetBinContent(ippbin);
426  if(!flag) continue;
427  // cout<<" pp ref flag="<<flag<<" >>> ";
428  // Get pp reference relative systematics
429  Double_t ppSystTotLow=0., ppSystTotHigh=0.;
430  FindGraphRelativeUnc(gDataSystematicsPP[j],ptval,ppSystTotLow,ppSystTotHigh);
431  ppSystRawYieldCutVar[j] = ppSystTotLow > ppSystTotHigh ? ppSystTotLow : ppSystTotHigh ;
432  ppSystRawYieldCutVarPid[j] = ppSystRawYieldCutVar[j];
433  }
434  }
435  //
436  // Loop per meson to get the Raa values and uncertainties for the given pt bin
437  //
438  if(isDebug) cout<<" Retrieving all Raa values and uncertainties"<<endl;
439  for(Int_t j=0; j<3; j++) {
440  // Get value, stat unc and weight
441  if(!isDmeson[j]) continue;
442  if(!hDmesonRaa[j]) continue;
443  histoBin = hDmesonRaa[j]->FindBin(ptval);
444  RaaDmeson[j] = hDmesonRaa[j]->GetBinContent( histoBin );
445  if (RaaDmeson[j]<=0) continue;
446  RaaDmesonStat[j] = hDmesonRaa[j]->GetBinError( histoBin );
447  // Get global systematics
448  FindGraphRelativeUnc(gRAB_GlobalSystematics[j],ptval,RaaDmesonSystLow[j],RaaDmesonSystHigh[j]);
449  // Evaluate the weight
450  weight[j] = GetWeight(averageOption,ptval,hDmesonRaa[j],
451  RaaDmesonSystLow[j],RaaDmesonSystHigh[j],
452  ppSystRawYieldCutVar[j],ppSystRawYieldCutVarPid[j],
453  ABSystRawYieldCutVar[j],ABSystRawYieldCutVarPid[j]);
454  cout<<" raa "<<filenamesSuffixes[j]<<" meson = "<<RaaDmeson[j]<<" +-"<<RaaDmesonStat[j]<<"(stat) -> (weight="<<weight[j]<<") ,";
455  // Get pp reference relative systematics
456  FindGraphRelativeUnc(gDataSystematicsPP[j],ptval,ppSystLow[j],ppSystHigh[j]);
457  // Get pp-extrapolation relative uncertainty
458  FindGraphRelativeUnc(gScalingUncPP[j],ptval,ScalingHigh[j],ScalingLow[j]); // exchanging low-high bc has oposite influence on Raa
459  if(isPPRefExtrap[j]){
460  Int_t ippbin = hCombinedReferenceFlag[j]->FindBin(ptval);
461  Bool_t flag = hCombinedReferenceFlag[j]->GetBinContent(ippbin);
462  if(isDebug) cout<< " bin="<<j<<" pp ref flag on? "<<flag<<endl;
463  if(flag){ ScalingHigh[j]=0.; ScalingLow[j]=0.; ppTracking[j][ipt]=0.; }
464  }
465  // Get pp reference systematics minus tracking systematics minus extrapolation uncertainties
466  ppSystUncorrLow[j] = TMath::Sqrt( ppSystLow[j]*ppSystLow[j]
467  - ScalingLow[j]*ScalingLow[j]
468  - ppTracking[j][ipt]*ppTracking[j][ipt] );
469  ppSystUncorrHigh[j] = TMath::Sqrt( ppSystHigh[j]*ppSystHigh[j]
470  - ScalingHigh[j]*ScalingHigh[j]
471  - ppTracking[j][ipt]*ppTracking[j][ipt] );
472  // Get AB relative systematics
473  FindGraphRelativeUnc(gDataSystematicsAB[j],ptval,ABSystLow[j],ABSystHigh[j]);
474  // Get AB relative systematics minus tracking systematics
475  ABSystUncorrLow[j] = TMath::Sqrt( ABSystLow[j]*ABSystLow[j]
476  - ABTracking[j][ipt]*ABTracking[j][ipt] );
477  ABSystUncorrHigh[j] = TMath::Sqrt( ABSystHigh[j]*ABSystHigh[j]
478  - ABTracking[j][ipt]*ABTracking[j][ipt] );
479  // Get Feed-Down and Eloss relative uncertainties on the Raa
480  FindGraphRelativeUnc(gRABFeedDownSystematicsElossHypothesis[j],ptval,RabFdElossLow[j],RabFdElossHigh[j]);
481  //
482  // Check with global Raa uncertainties
483  FindGraphRelativeUnc(gRAB_GlobalSystematics[j],ptval,RabGlobalLow[j],RabGlobalHigh[j]);
484  Double_t testLow = TMath::Sqrt( RabFdElossLow[j]*RabFdElossLow[j] + ABSystLow[j]*ABSystLow[j]
485  + ppSystLow[j]*ppSystLow[j] + ScalingLow[j]*ScalingLow[j] );
486  Double_t testHigh = TMath::Sqrt( RabFdElossHigh[j]*RabFdElossHigh[j] + ABSystHigh[j]*ABSystHigh[j]
487  + ppSystHigh[j]*ppSystHigh[j] + ScalingHigh[j]*ScalingHigh[j] );
488  if (TMath::Abs( testLow - RabGlobalLow[j] ) > 0.015) {
489  cout << endl<<" >>>> Error on the global Raa uncertainties low : test-sum = "<< testLow<<", global = "<< RabGlobalLow[j]<<" ppref="<<ppSystLow[j]<<endl;
490  }
491  if (TMath::Abs( testHigh - RabGlobalHigh[j] ) > 0.015) {
492  cout << endl<<" >>>> Error on the global Raa uncertainties high : test-sum = "<< testHigh<<", global = "<< RabGlobalHigh[j]<<" ppref="<<ppSystHigh[j]<<endl<<endl;
493  }
494  //
495  histoBin = -1;
496  }
497  cout<<endl;
498 
499  //
500  // Evaluate Dmeson average
501  //
502  if(isDebug) cout<<" Evaluating the average"<<endl;
503  for(Int_t j=0; j<3; j++){
504  if(!isDmeson[j]) continue;
505  if( !(RaaDmeson[j]>0.) ) continue;
506  weightTot += weight[j];
507  // weighted average
508  average += RaaDmeson[j]*weight[j];
509  // stat absolute uncertainty (uncorrelated) : sum in quadrature
510  averageStat += (RaaDmesonStat[j]*weight[j])*(RaaDmesonStat[j]*weight[j]);
511  // pp tracking relative uncertainty (correlated) : linear sum
512  ppTrackingAv += ppTracking[j][ipt]*RaaDmeson[j]*weight[j];
513  // AB tracking relative uncertainty (correlated) : linear sum
514  ABTrackingAv += ABTracking[j][ipt]*RaaDmeson[j]*weight[j];
515  // pp scaling relative uncertainty (correlated) : linear sum
516  scalingLowAv += ScalingLow[j]*RaaDmeson[j]*weight[j];
517  scalingHighAv += ScalingHigh[j]*RaaDmeson[j]*weight[j];
518  // Get pp and AB relative uncorrelated systematics : sum in quadrature
519  raaSystUncorrLow += (ppSystUncorrLow[j]*RaaDmeson[j]*weight[j])*(ppSystUncorrLow[j]*RaaDmeson[j]*weight[j])
520  + (ABSystUncorrLow[j]*RaaDmeson[j]*weight[j])*(ABSystUncorrLow[j]*RaaDmeson[j]*weight[j]);
521  raaSystUncorrHigh += (ppSystUncorrHigh[j]*RaaDmeson[j]*weight[j])*(ppSystUncorrHigh[j]*RaaDmeson[j]*weight[j])
522  + (ABSystUncorrHigh[j]*RaaDmeson[j]*weight[j])*(ABSystUncorrHigh[j]*RaaDmeson[j]*weight[j]);
523  ppDataSystAvLow += (ppSystUncorrLow[j]*RaaDmeson[j]*weight[j])*(ppSystUncorrLow[j]*RaaDmeson[j]*weight[j]);
524  ABDataSystAvLow += (ABSystUncorrLow[j]*RaaDmeson[j]*weight[j])*(ABSystUncorrLow[j]*RaaDmeson[j]*weight[j]);
525  ppDataSystAvHigh += (ppSystUncorrHigh[j]*RaaDmeson[j]*weight[j])*(ppSystUncorrHigh[j]*RaaDmeson[j]*weight[j]);
526  ABDataSystAvHigh += (ABSystUncorrHigh[j]*RaaDmeson[j]*weight[j])*(ABSystUncorrHigh[j]*RaaDmeson[j]*weight[j]);
527  // Beauty uncertainties: evaluate Raa average for the upper / lower bands
528  raabeautyLow += (1-RabFdElossLow[j])*RaaDmeson[j]*weight[j];
529  raabeautyHigh += (1+RabFdElossHigh[j])*RaaDmeson[j]*weight[j];
530  }
531  average /= weightTot;
532  averageStat = TMath::Sqrt(averageStat)/weightTot;
533  ppTrackingAv /= weightTot;
534  ABTrackingAv /= weightTot;
535  scalingLowAv /= weightTot;
536  scalingHighAv /= weightTot;
537  raaSystUncorrLow = TMath::Sqrt(raaSystUncorrLow)/weightTot;
538  raaSystUncorrHigh = TMath::Sqrt(raaSystUncorrHigh)/weightTot;
539  ppDataSystAvLow = TMath::Sqrt(ppDataSystAvLow)/weightTot;
540  ppDataSystAvHigh = TMath::Sqrt(ppDataSystAvHigh)/weightTot;
541  ABDataSystAvLow = TMath::Sqrt(ABDataSystAvLow)/weightTot;
542  ABDataSystAvHigh = TMath::Sqrt(ABDataSystAvHigh)/weightTot;
543 
544  // finalization beauty uncertainties
545  raabeautyLow /= weightTot;
546  raabeautyHigh /= weightTot;
547  Double_t RaaBeauty[3] = { average, raabeautyLow, raabeautyHigh };
548  Double_t beautyUncLow = average-TMath::MinElement(3,RaaBeauty);
549  Double_t beautyUncHigh = TMath::MaxElement(3,RaaBeauty)-average;
550 
551  // finalization global uncertainties
552  Double_t totalUncLow = TMath::Sqrt( raaSystUncorrLow*raaSystUncorrLow
553  + ppTrackingAv*ppTrackingAv + ABTrackingAv*ABTrackingAv
554  + scalingLowAv*scalingLowAv
555  + beautyUncLow*beautyUncLow );
556  Double_t totalUncHigh = TMath::Sqrt( raaSystUncorrHigh*raaSystUncorrHigh
557  + ppTrackingAv*ppTrackingAv + ABTrackingAv*ABTrackingAv
558  + scalingHighAv*scalingHighAv
559  + beautyUncHigh*beautyUncHigh );
560  if(isDebug) cout<<" Raa="<<average<<" +-"<<averageStat<<"(stat) +"<<ppDataSystAvHigh<<" -"<<ppDataSystAvLow<<" (pp-data) +-"<<ppTrackingAv<<" (pp-track) +"<<ABDataSystAvHigh<<" -"<<ABDataSystAvLow<<" (ab-data) +-"<<ABTrackingAv<<" (ab-track) +"<<scalingHighAv<<" -"<<scalingLowAv<<" (scal) +"<<beautyUncHigh<<" -"<<beautyUncLow<<" (fd)"<<endl;
561  //
562  // Fill output histos/graphs
563  //
564  histoBin = hDmesonAverageRAB->FindBin(ptval);
565  hDmesonAverageRAB->SetBinContent(histoBin,average);
566  hDmesonAverageRAB->SetBinError(histoBin,averageStat);
567  Double_t ept = hDmesonAverageRAB->GetBinWidth(histoBin)/2.;
568  ept=0.3;
569  gRAB_DmesonAverage_GlobalSystematics->SetPoint(ipt,ptval,average);
570  gRAB_DmesonAverage_GlobalSystematics->SetPointError(ipt,ept,ept,totalUncLow,totalUncHigh);
571  ept=0.2;
572  gRAB_DmesonAverage_FeedDownSystematicsElossHypothesis->SetPoint(ipt,ptval,average);
573  gRAB_DmesonAverage_FeedDownSystematicsElossHypothesis->SetPointError(ipt,ept,ept,beautyUncLow,beautyUncHigh);
574  //
575  ept=0.1;
576  gRAB_DmesonAverage_TrackingSystematicsPP->SetPoint(ipt,ptval,average);
577  gRAB_DmesonAverage_TrackingSystematicsPP->SetPointError(ipt,ept,ept,ppTrackingAv,ppTrackingAv);
578  gRAB_DmesonAverage_TrackingSystematicsAB->SetPoint(ipt,ptval,average);
579  gRAB_DmesonAverage_TrackingSystematicsAB->SetPointError(ipt,ept,ept,ABTrackingAv,ABTrackingAv);
580  gRAB_DmesonAverage_ScalingSystematicsPP->SetPoint(ipt,ptval,average);
581  gRAB_DmesonAverage_ScalingSystematicsPP->SetPointError(ipt,ept,ept,scalingLowAv,scalingHighAv);
582  gRAB_DmesonAverage_DataSystematicsPP->SetPoint(ipt,ptval,average);
583  gRAB_DmesonAverage_DataSystematicsPP->SetPointError(ipt,ept,ept,ppDataSystAvLow,ppDataSystAvHigh);
584  gRAB_DmesonAverage_DataSystematicsAB->SetPoint(ipt,ptval,average);
585  gRAB_DmesonAverage_DataSystematicsAB->SetPointError(ipt,ept,ept,ABDataSystAvLow,ABDataSystAvHigh);
586  histoBin = -1;
587  //
588  // Printout
589  cout<< " pt min (GeV/c), pt max (GeV/c), Raa(Daverage), +- (stat), + (syst) , - (syst) "<<endl;
590  cout<< ptbinlimits[ipt] <<" "<< ptbinlimits[ipt+1]<< " "<< average<< " "<< averageStat<< " "<< totalUncHigh<<" "<<totalUncLow<<endl;
591  fprintf(resultFile,"%02.0f %02.0f %2.2f %2.2f %2.2f %2.2f\n",ptbinlimits[ipt],ptbinlimits[ipt+1],average,averageStat,totalUncHigh,totalUncLow);
592  } // end loop on pt bins
593 
594  fclose(resultFile);
595 
596  //
597  // Now can start drawing
598  //
599  TCanvas *cAvCheck = new TCanvas("cAvCheck","Average Dmeson check");
600  hDmesonAverageRAB->SetLineColor(kBlack);
601  hDmesonAverageRAB->SetMarkerStyle(20);
602  hDmesonAverageRAB->SetMarkerColor(kBlack);
603  hDmesonAverageRAB->Draw("e");
604  for(Int_t j=0; j<3; j++) {
605  if(!isDmeson[j]) continue;
606  hDmesonRaa[j]->SetLineColor(kBlack);
607  hDmesonRaa[j]->SetMarkerColor(2+j);
608  hDmesonRaa[j]->SetMarkerStyle(21+j);
609  gRAB_GlobalSystematics[j]->SetFillStyle(0);
610  gRAB_GlobalSystematics[j]->SetLineColor(2+j);
611  gRAB_GlobalSystematics[j]->Draw("2");
612  hDmesonRaa[j]->Draw("esame");
613  }
614  gRAB_DmesonAverage_GlobalSystematics->SetFillStyle(0);
615  gRAB_DmesonAverage_GlobalSystematics->SetLineColor(kBlack);
616  gRAB_DmesonAverage_GlobalSystematics->Draw("2");
617  hDmesonAverageRAB->Draw("esame");
618  cAvCheck->Update();
619 
620  TCanvas *cAv = new TCanvas("cAv","Average Dmeson");
621  hDmesonAverageRAB->Draw("e");
622  gRAB_DmesonAverage_FeedDownSystematicsElossHypothesis->SetFillStyle(1001);
623  gRAB_DmesonAverage_FeedDownSystematicsElossHypothesis->SetFillColor(kMagenta-7);
624  gRAB_DmesonAverage_FeedDownSystematicsElossHypothesis->Draw("2");
625  gRAB_DmesonAverage_GlobalSystematics->Draw("2");
626  hDmesonAverageRAB->Draw("esame");
627  cAv->Update();
628 
629  //
630  // Now can start saving the output
631  //
632  TFile *fout = new TFile(outfile,"recreate");
633  hDmesonAverageRAB->Write();
634  gRABNorm->Write();
635  gRAB_DmesonAverage_GlobalSystematics->Write();
636  gRAB_DmesonAverage_FeedDownSystematicsElossHypothesis->Write();
637  gRAB_DmesonAverage_TrackingSystematicsPP->Write();
638  gRAB_DmesonAverage_TrackingSystematicsAB->Write();
639  gRAB_DmesonAverage_ScalingSystematicsPP->Write();
640  gRAB_DmesonAverage_DataSystematicsPP->Write();
641  gRAB_DmesonAverage_DataSystematicsAB->Write();
642  fout->Write();
643 
644 }
const Color_t cc[]
Definition: DrawKs.C:1
double Double_t
Definition: External.C:58
Double_t GetWeight(Int_t averageoption, Double_t pt, TH1D *hRaa, Double_t raaSystLow, Double_t raaSystHigh, Double_t ppSystRawYieldCutVar, Double_t ppSystRawYieldCutVarPid, Double_t ABSystRawYieldCutVar, Double_t ABSystRawYieldCutVarPid)
Double_t ptbinlimits[nbins+1]
void FindGraphRelativeUnc(TGraphAsymmErrors *gr, Double_t pt, Double_t &uncLow, Double_t &uncHigh)
centrality
void SetCentrality(TString centrality)
Definition: AliHFSystErr.h:60
void SetIsPbPb2010EnergyScan(Bool_t flag)
Definition: AliHFSystErr.h:82
Double_t GetCutsEffErr(Double_t pt) const
Int_t FindGraphBin(TGraphAsymmErrors *gr, Double_t pt)
Bool_t isDebug
int Int_t
Definition: External.C:63
Definition: External.C:204
centestimator
Double_t GetTrackingEffErr(Double_t pt) const
Definition: External.C:212
void AverageDmesonRaa(const char *fD0Raa="", const char *fD0ppRef="", const char *fDplusRaa="", const char *fDplusppRef="", const char *fDstarRaa="", const char *fDstarppRef="", const char *outfile="", Int_t averageOption=kRelativeStatUnc, Int_t cc=kpPb0100, Int_t ccestimator=kV0M, Bool_t isReadAllPPUnc=false, Bool_t isPPRefExtrapD0=false, Bool_t isPPRefExtrapDplus=false, Bool_t isPPRefExtrapDstar=false)
void Init(Int_t decay)
Function to initialize the variables/histograms.
Double_t GetRawYieldErr(Double_t pt) const
void SetCollisionType(Int_t type)
Definition: AliHFSystErr.h:51
AverageOption
const Int_t nbins
bool Bool_t
Definition: External.C:53
Double_t GetPIDEffErr(Double_t pt) const