Incrementing class version.
[u/mrichter/AliRoot.git] / PWGHF / correlationHF / AliHFDmesonCorrAverage.cxx
CommitLineData
ec9cbe08 1/**************************************************************************
2 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16
17/* $Id: $ */
18
19/////////////////////////////////////////////////////////////
20// class to average D meson -hadron correlations
21//
22// Author: A. Rossi, andrea.rossi@cern.ch
23/////////////////////////////////////////////////////////////
24
25#include "AliHFDmesonCorrAverage.h"
26#include "AliHFDhadronCorrSystUnc.h"
27#include <TMath.h>
28#include <TH1D.h>
29#include <TCanvas.h>
30#include <TGraphAsymmErrors.h>
31
32
33//using std::cout;
34//using std::endl;
35
36ClassImp(AliHFDmesonCorrAverage)
37
38AliHFDmesonCorrAverage::AliHFDmesonCorrAverage() :
39TNamed(),
40 fSystDzero(0),
41 fSystDstar(0),
42 fSystDplus(0),
43 fSystDaverage(0),
44 fincludeDzero(kFALSE),
45 fincludeDstar(kFALSE),
46 fincludeDplus(kFALSE),
47 fmethod(),
48 fptminD(),
49 fptmaxD(),
50 fptminAsso(),
51 fptmaxAsso(),
52 fhDzero(0x0),
53 fhDstar(0x0),
54 fhDplus(0x0),
55 fhDaverage(0x0),
56 fgrTotSystAverage(0x0),
57 fgrFDSystAverage(0x0),
58 fgrNonFDSystAverage(0x0),
59 fweightsDzeroStat(0x0),
60 fweightsDstarStat(0x0),
61 fweightsDplusStat(0x0),
62 fweightsDzeroSystYield(0x0),
63 fweightsDstarSystYield(0x0),
64 fweightsDplusSystYield(0x0),
65 fweightsDzeroSystBkg(0x0),
66 fweightsDstarSystBkg(0x0),
67 fweightsDplusSystBkg(0x0),
68 fnbinsphi(),
69 fsys(-2),
acd5f443 70 fyear(-2),
71 fSystAlreadySet(kFALSE)
ec9cbe08 72{// standard constructor
73
74}
75
76
77AliHFDmesonCorrAverage::AliHFDmesonCorrAverage(const char* name) :
78 TNamed(name,name),
79 fSystDzero(0),
80 fSystDstar(0),
81 fSystDplus(0),
82 fSystDaverage(0),
83 fincludeDzero(kFALSE),
84 fincludeDstar(kFALSE),
85 fincludeDplus(kFALSE),
86 fmethod(),
87 fptminD(),
88 fptmaxD(),
89 fptminAsso(),
90 fptmaxAsso(),
91 fhDzero(0x0),
92 fhDstar(0x0),
93 fhDplus(0x0),
94 fhDaverage(0x0),
95 fgrTotSystAverage(0x0),
96 fgrFDSystAverage(0x0),
97 fgrNonFDSystAverage(0x0),
98 fweightsDzeroStat(0x0),
99 fweightsDstarStat(0x0),
100 fweightsDplusStat(0x0),
101 fweightsDzeroSystYield(0x0),
102 fweightsDstarSystYield(0x0),
103 fweightsDplusSystYield(0x0),
104 fweightsDzeroSystBkg(0x0),
105 fweightsDstarSystBkg(0x0),
106 fweightsDplusSystBkg(0x0),
107 fnbinsphi(),
108 fsys(-2),
acd5f443 109 fyear(-2),
110 fSystAlreadySet(kFALSE)
ec9cbe08 111{// default constructor
112
113}
114
115AliHFDmesonCorrAverage::~AliHFDmesonCorrAverage(){
116
117 delete fSystDzero;
118 delete fSystDstar;
119 delete fSystDplus;
120 delete fSystDaverage;
121 delete fhDzero;
122 delete fhDstar;
123 delete fhDplus;
124 delete fhDaverage;
125 delete fgrTotSystAverage;
126 delete fgrFDSystAverage;
127 delete fgrNonFDSystAverage;
128 delete fweightsDzeroStat;
129 delete fweightsDstarStat;
130 delete fweightsDplusStat;
131 delete fweightsDzeroSystYield;
132 delete fweightsDstarSystYield;
133 delete fweightsDplusSystYield;
134 delete fweightsDzeroSystBkg;
135 delete fweightsDstarSystBkg;
136 delete fweightsDplusSystBkg;
137
138}
139
140
141
142Bool_t AliHFDmesonCorrAverage::InitSystematicUncertainty(Int_t system,Int_t year){
acd5f443 143 if(!fSystAlreadySet){
144 if(system!=-1){
145 if(system==0){ //pp
146 if(fincludeDzero){
147 if(year==2010){
148 fSystDzero=new AliHFDhadronCorrSystUnc("fSystDzero");
149 fSystDzero->InitStandardUncertaintiesPP2010(0,(fptminD+fptmaxD)*0.5,fptminAsso);
150 fSystDzero->BuildTotalUncHisto();
151 fSystDzero->BuildTotalNonFDUncHisto();
152 fSystDzero->BuildGraphsUnc(fhDzero);
153 fSystDzero->BuildGraphsRelUnc();
154 }
155 else {
156 Printf("No values for syst unc foreseen for this dataset");
157 }
ec9cbe08 158 }
acd5f443 159
160 if(fincludeDstar){
161 if(year==2010){
162 fSystDstar=new AliHFDhadronCorrSystUnc("fSystDstar");
163 fSystDstar->InitStandardUncertaintiesPP2010(1,(fptminD+fptmaxD)*0.5,fptminAsso);
164 fSystDstar->BuildTotalUncHisto();
165 fSystDstar->BuildTotalNonFDUncHisto();
166 fSystDstar->BuildGraphsUnc(fhDstar);
167 fSystDstar->BuildGraphsRelUnc();
168 }
169 else {
170 Printf("No values for syst unc foreseen for this dataset");
171 }
ec9cbe08 172 }
acd5f443 173
174 if(fincludeDplus){
175 if(year==2010){
176 fSystDplus=new AliHFDhadronCorrSystUnc("fSystDplus");
177 fSystDplus->InitStandardUncertaintiesPP2010(2,(fptminD+fptmaxD)*0.5,fptminAsso);
178 fSystDplus->BuildTotalUncHisto();
179 fSystDplus->BuildTotalNonFDUncHisto();
180 fSystDplus->BuildGraphsUnc(fhDplus);
181 fSystDplus->BuildGraphsRelUnc();
182 }
183 else {
184 Printf("No values for syst unc foreseen for this dataset");
185 }
ec9cbe08 186 }
187 }
acd5f443 188 else if(system==1){ //p-Pb
189 if(fincludeDzero){
190 if(year==2013){
191 fSystDzero=new AliHFDhadronCorrSystUnc("fSystDzero");
192 fSystDzero->InitStandardUncertaintiesPPb2013(0,(fptminD+fptmaxD)*0.5,fptminAsso);
193 fSystDzero->BuildTotalUncHisto();
194 fSystDzero->BuildTotalNonFDUncHisto();
195 fSystDzero->BuildGraphsUnc(fhDzero);
196 fSystDzero->BuildGraphsRelUnc();
197 }
198 else {
199 Printf("No values for syst unc foreseen for this dataset");
200 }
ec9cbe08 201 }
acd5f443 202 if(fincludeDstar){
203 if(year==2013){
204 fSystDstar=new AliHFDhadronCorrSystUnc("fSystDstar");
205 fSystDstar->InitStandardUncertaintiesPPb2013(1,(fptminD+fptmaxD)*0.5,fptminAsso);
206 fSystDstar->BuildTotalUncHisto();
207 fSystDstar->BuildTotalNonFDUncHisto();
208 fSystDstar->BuildGraphsUnc(fhDstar);
209 fSystDstar->BuildGraphsRelUnc();
210 }
211 else {
212 Printf("No values for syst unc foreseen for this dataset");
213 }
ec9cbe08 214 }
acd5f443 215
216 if(fincludeDplus){
217 if(year==2013){
218 fSystDplus=new AliHFDhadronCorrSystUnc("fSystDplus");
219 fSystDplus->InitStandardUncertaintiesPPb2013(2,(fptminD+fptmaxD)*0.5,fptminAsso);
220 fSystDplus->BuildTotalUncHisto();
221 fSystDplus->BuildTotalNonFDUncHisto();
222 fSystDplus->BuildGraphsUnc(fhDplus);
223 fSystDplus->BuildGraphsRelUnc();
224 }
225 else {
226 Printf("No values for syst unc foreseen for this dataset");
227 }
ec9cbe08 228 }
229 }
acd5f443 230 else {
231 Printf("Cannot initiate syst uncertainties: wrong system selected");
232 return kFALSE;
ec9cbe08 233 }
234 }
ec9cbe08 235 }
236
237 fSystDaverage=new AliHFDhadronCorrSystUnc("fSystDaverage");
238 if(fincludeDzero)fSystDaverage->SetHistoTemplate(fSystDzero->GetHistoTemplate(),"fhDeltaPhiTemplateDaverage");
239 else if(fincludeDstar)fSystDaverage->SetHistoTemplate(fSystDstar->GetHistoTemplate(),"fhDeltaPhiTemplateDaverage");
240 else if(fincludeDplus)fSystDaverage->SetHistoTemplate(fSystDplus->GetHistoTemplate(),"fhDeltaPhiTemplateDaverage");
241
242 fSystDaverage->InitEmptyHistosFromTemplate();
243
244 return kTRUE;
245}
246
247
248void AliHFDmesonCorrAverage::InitAverageHisto(TH1D *h){
249 fhDaverage=(TH1D*)h->Clone("fhDaverage");
250 fhDaverage->SetXTitle("#Delta#varphi = #varphi_{assoc} - #varphi_{D}");
251 fhDaverage->SetYTitle("#frac{1}{N_{D}}#frac{dN^{h}}{d#Delta#varphi} (rad^{-1})");
252 fhDaverage->SetTitle("");
253 fhDaverage->SetMinimum(0);
254
255 fhDaverage->Reset(0);
256 fhDaverage->SetLineWidth(2);
257 fhDaverage->SetLineColor(kBlack);
258 fhDaverage->SetMarkerStyle(20);
259 fhDaverage->SetMarkerSize(1.2);
260 fhDaverage->SetMarkerColor(kBlack);
261
262}
263
264void AliHFDmesonCorrAverage::CalculateAverage(){
265 if(fmethod!=10&&fmethod!=20){
266 Printf("This option for the calculation of the average has not been implemented yet");
267 return;
268 }
269 // check that histos were set
270 if(fincludeDzero&&(!fhDzero)){
271 Printf("Dzero histo not set");
272 return;
273 }
274 if(fincludeDstar&&(!fhDstar)){
275 Printf("Dstar histo not set");
276 return;
277 }
278 if(fincludeDplus&&(!fhDplus)){
279 Printf("Dplus histo not set");
280 return;
281 }
282
283
284 // check that they have the same binning
285 if(fincludeDzero&&fincludeDstar){
286 if(fhDzero->GetNbinsX()!=fhDstar->GetNbinsX()){
287 Printf("Dzero and Dstar histos have different binning");
288 return ;
289 }
290 fnbinsphi=fhDzero->GetNbinsX();
291 for(Int_t j=1;j<=fnbinsphi;j++){
292 if(TMath::Abs(fhDzero->GetBinLowEdge(j)-fhDstar->GetBinLowEdge(j))>0.001){
293 Printf("Dzero and Dstar histos have different binning");
294 return;
295 }
296 }
297 InitAverageHisto(fhDzero);
298 }
299 if(fincludeDzero&&fincludeDplus){
300 if(fhDzero->GetNbinsX()!=fhDplus->GetNbinsX()){
301 Printf("Dzero and Dplus histos have different binning");
302 return ;
303 }
304 fnbinsphi=fhDzero->GetNbinsX();
305 for(Int_t j=1;j<=fnbinsphi;j++){
306 if(TMath::Abs(fhDzero->GetBinLowEdge(j)-fhDplus->GetBinLowEdge(j))>0.001){
307 Printf("Dzero and Dplus histos have different binning");
308 return;
309 }
310 }
311 if(!fhDaverage)InitAverageHisto(fhDzero);
312 }
313 if(fincludeDstar&&fincludeDplus){
314 if(fhDstar->GetNbinsX()!=fhDplus->GetNbinsX()){
315 Printf("Dstar and Dplus histos have different binning");
316 return ;
317 }
318 fnbinsphi=fhDstar->GetNbinsX();
319 for(Int_t j=1;j<=fnbinsphi;j++){
320 if(TMath::Abs(fhDstar->GetBinLowEdge(j)-fhDplus->GetBinLowEdge(j))>0.001){
321 Printf("Dplus and Dstar histos have different binning");
322 return;
323 }
324 }
325 if(!fhDaverage)InitAverageHisto(fhDstar);
326 }
327
328 if(!(fincludeDstar||fincludeDplus||fincludeDzero)){
329 Printf("No mesons choosen for average");
330 return ;
331 }
332 Int_t nmeson=0;
333 if(fincludeDzero)nmeson++;
334 if(fincludeDstar)nmeson++;
335 if(fincludeDplus)nmeson++;
336
337 printf("Settin syst \n");
acd5f443 338
ec9cbe08 339 if(!InitSystematicUncertainty(fsys,fyear)){
340 Printf("Cannot init syst uncertainties \n");
341 return;
342 }
ec9cbe08 343
344 SetWeights();
345 printf("Weights set \n");
346
347 for(Int_t j=1;j<=fnbinsphi;j++){
348 Double_t value=0.,errStatValue=0.,errSystYieldValue=0.,errSystBkgValue=0.,weight=0.;
349 Double_t systMCclosureMin=0.,systMCcorrectionsMin=0.,systFDmin=0.,systMCDefficiencyMin=0.,systSecContaminationMin=0.;
350 Double_t systMCclosureMax=0.,systMCcorrectionsMax=0.,systFDmax=0.,systMCDefficiencyMax=0.,systSecContaminationMax=0.;
351 Double_t sumweights=0;
352 Double_t errsyst;
353
354 if(fincludeDzero){
355 // stat error + yield unc + bkg subtr
356 weight=1./(1./fweightsDzeroStat[j-1]+1./fweightsDzeroSystYield[j-1]+1./fweightsDzeroSystBkg[j-1]);// need to do this way since we stored separately the stat and syst weight (=1/variance)
357 value+=fhDzero->GetBinContent(j)*weight;
acd5f443 358 errStatValue+=1./fweightsDzeroStat[j-1]*weight*weight;
359 errSystYieldValue+=1./fweightsDzeroSystYield[j-1]*weight*weight;
360 errSystBkgValue+=1./fweightsDzeroSystBkg[j-1]*weight*weight;
ec9cbe08 361 sumweights+=weight;
362
363 printf("Dzero the value is: %f, weight: %f \n",value, weight);
364 // MCcorrections : associated tracks, correlated
365 errsyst=TMath::Abs(fhDzero->GetBinContent(j)*fSystDzero->GetHistoMCcorrectionsMin()->GetBinContent(j));
366 systMCcorrectionsMin+=errsyst*weight;
367 errsyst=fhDzero->GetBinContent(j)*fSystDzero->GetHistoMCcorrectionsMax()->GetBinContent(j);
368 systMCcorrectionsMax+=errsyst*weight;
369
370 printf(" Dzero trackeff the min syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDzero->GetHistoMCcorrectionsMin()->GetBinContent(j), systMCcorrectionsMin);
371 printf(" Dzero trackeff the max syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDzero->GetHistoMCcorrectionsMax()->GetBinContent(j), systMCcorrectionsMax);
372
373 // MCDefficiency : uncorrelated
374 errsyst=fhDzero->GetBinContent(j)*fSystDzero->GetHistoMCDefficiencyMin()->GetBinContent(j);
375 systMCDefficiencyMin+=errsyst*errsyst*weight*weight;
376 errsyst=fhDzero->GetBinContent(j)*fSystDzero->GetHistoMCDefficiencyMax()->GetBinContent(j);
377 systMCDefficiencyMax+=errsyst*errsyst*weight*weight;
378
379 // SecContamination : correlated
380 errsyst=TMath::Abs(fhDzero->GetBinContent(j)*fSystDzero->GetHistoSecContaminationMin()->GetBinContent(j));
381 systSecContaminationMin+=errsyst*weight;
382 errsyst=fhDzero->GetBinContent(j)*fSystDzero->GetHistoSecContaminationMax()->GetBinContent(j);
383 systSecContaminationMax+=errsyst*weight;
384
385 // MC closure: fully correlated
386 errsyst=TMath::Abs(fhDzero->GetBinContent(j)*fSystDzero->GetHistoMCclosureTestMin()->GetBinContent(j));
387 systMCclosureMin+=errsyst*weight;
388 errsyst=fhDzero->GetBinContent(j)*fSystDzero->GetHistoMCclosureTestMax()->GetBinContent(j);
389 systMCclosureMax+=errsyst*weight;
390
391 printf(" Dzero Mcclosure the min syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDzero->GetHistoMCclosureTestMin()->GetBinContent(j), systMCclosureMin);
392 printf(" Dzero Mcclosure the max syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDzero->GetHistoMCclosureTestMax()->GetBinContent(j), systMCclosureMax);
393
394 // FD (assume full correlations)
395 errsyst=TMath::Abs(fhDzero->GetBinContent(j)*fSystDzero->GetHistoFDmin()->GetBinContent(j));
396 systFDmin+=errsyst*weight;
397 errsyst=fhDzero->GetBinContent(j)*fSystDzero->GetHistoFDmax()->GetBinContent(j);
398 systFDmax+=errsyst*weight;
399 printf(" Dzero feeddown the min syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDzero->GetHistoFDmin()->GetBinContent(j), systFDmin);
400 printf(" Dzero feeddown the max syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDzero->GetHistoFDmax()->GetBinContent(j), systFDmax);
401
402 printf("Dzero the value is: %f, weight: %f \n",value, weight);
403 }
404 if(fincludeDstar){
405
406 // stat error + yield unc + bkg subtr
407 weight=1./(1./fweightsDstarStat[j-1]+1./fweightsDstarSystYield[j-1]+1./fweightsDstarSystBkg[j-1]);// need to do this way since we stored separately the stat and syst weight (=1/variance)
408 value+=fhDstar->GetBinContent(j)*weight;
acd5f443 409 errStatValue+=1./fweightsDstarStat[j-1]*weight*weight;
410 errSystYieldValue+=1./fweightsDstarSystYield[j-1]*weight*weight;
411 errSystBkgValue+=1./fweightsDstarSystBkg[j-1]*weight*weight;
ec9cbe08 412 sumweights+=weight;
413
414 printf("Dstar the value is: %f, weight: %f \n",value, weight);
415
416 // MCcorrections: associated tracks, correlated
417 errsyst=TMath::Abs(fhDstar->GetBinContent(j)*fSystDstar->GetHistoMCcorrectionsMin()->GetBinContent(j));
418 systMCcorrectionsMin+=errsyst*weight;
419 errsyst=TMath::Abs(fhDstar->GetBinContent(j)*fSystDstar->GetHistoMCcorrectionsMax()->GetBinContent(j));
420 systMCcorrectionsMax+=errsyst*weight;
421
422 printf(" Dstar trackeff the min syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDstar->GetHistoMCcorrectionsMin()->GetBinContent(j), systMCcorrectionsMin);
423 printf(" Dstar trackeff the max syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDstar->GetHistoMCcorrectionsMax()->GetBinContent(j), systMCcorrectionsMax);
424
425 // MCDefficiency : uncorrelated
426 errsyst=fhDstar->GetBinContent(j)*fSystDstar->GetHistoMCDefficiencyMin()->GetBinContent(j);
427 systMCDefficiencyMin+=errsyst*errsyst*weight*weight;
428 errsyst=fhDstar->GetBinContent(j)*fSystDstar->GetHistoMCDefficiencyMax()->GetBinContent(j);
429 systMCDefficiencyMax+=errsyst*errsyst*weight*weight;
430
431 // SecContamination : correlated
432 errsyst=TMath::Abs(fhDstar->GetBinContent(j)*fSystDstar->GetHistoSecContaminationMin()->GetBinContent(j));
433 systSecContaminationMin+=errsyst*weight;
434 errsyst=fhDstar->GetBinContent(j)*fSystDstar->GetHistoSecContaminationMax()->GetBinContent(j);
435 systSecContaminationMax+=errsyst*weight;
436
437 // MC closure
438 errsyst=TMath::Abs(fhDstar->GetBinContent(j)*fSystDstar->GetHistoMCclosureTestMin()->GetBinContent(j));
439 systMCclosureMin+=errsyst*weight;
440 errsyst=fhDstar->GetBinContent(j)*fSystDstar->GetHistoMCclosureTestMax()->GetBinContent(j);
441 systMCclosureMax+=errsyst*weight;
442
443 printf(" Dstar Mcclosure the min syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDstar->GetHistoMCclosureTestMin()->GetBinContent(j), systMCclosureMin);
444 printf(" Dstar Mcclosure the max syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDstar->GetHistoMCclosureTestMax()->GetBinContent(j), systMCclosureMax);
445
446
447 // FD (assume full correlations)
448 errsyst=TMath::Abs(fhDstar->GetBinContent(j)*fSystDstar->GetHistoFDmin()->GetBinContent(j));
449 systFDmin+=errsyst*weight;
450 errsyst=fhDstar->GetBinContent(j)*fSystDstar->GetHistoFDmax()->GetBinContent(j);
451 systFDmax+=errsyst*weight;
452
453 printf(" Dstar feeddown the min syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDstar->GetHistoFDmin()->GetBinContent(j), systFDmin);
454 printf(" Dstar feeddown the max syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDstar->GetHistoFDmax()->GetBinContent(j), systFDmax);
455
456 printf("Dstar the value is: %f, weight: %f \n",value, weight);
457 }
458 if(fincludeDplus){
459 // stat error + yield unc + bkg subtr
460 weight=1./(1./fweightsDplusStat[j-1]+1./fweightsDplusSystYield[j-1]+1./fweightsDplusSystBkg[j-1]);// need to do this way since we stored separately the stat and syst weight (=1/variance)
461 value+=fhDplus->GetBinContent(j)*weight;
acd5f443 462 errStatValue+=1./fweightsDplusStat[j-1]*weight*weight;
463 errSystYieldValue+=1./fweightsDplusSystYield[j-1]*weight*weight;
464 errSystBkgValue+=1./fweightsDplusSystBkg[j-1]*weight*weight;
ec9cbe08 465 sumweights+=weight;
466
467 printf("Dplus the value is: %f, weight: %f \n",value, weight);
468
469 // MCcorrections : associated tracsk, correlated
470 errsyst=TMath::Abs(fhDplus->GetBinContent(j)*fSystDplus->GetHistoMCcorrectionsMin()->GetBinContent(j));
471 systMCcorrectionsMin+=errsyst*weight;
472 errsyst=fhDplus->GetBinContent(j)*fSystDplus->GetHistoMCcorrectionsMax()->GetBinContent(j);
473 systMCcorrectionsMax+=errsyst*weight;
474
475 printf(" Dplus trackeff the min syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDplus->GetHistoMCcorrectionsMin()->GetBinContent(j), systMCcorrectionsMin);
476 printf(" Dplus trackeff the max syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDplus->GetHistoMCcorrectionsMax()->GetBinContent(j), systMCcorrectionsMax);
477
478 // MCDefficiency : uncorrelated
479 errsyst=fhDplus->GetBinContent(j)*fSystDplus->GetHistoMCDefficiencyMin()->GetBinContent(j);
480 systMCDefficiencyMin+=errsyst*errsyst*weight*weight;
481 errsyst=fhDplus->GetBinContent(j)*fSystDplus->GetHistoMCDefficiencyMax()->GetBinContent(j);
482 systMCDefficiencyMax+=errsyst*errsyst*weight*weight;
483
484 printf(" Dplus cutvariat the min syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDplus->GetHistoMCDefficiencyMin()->GetBinContent(j), systMCDefficiencyMin);
485 printf(" Dplus cutvariat the max syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDplus->GetHistoMCDefficiencyMax()->GetBinContent(j), systMCDefficiencyMax);
486
487 // SecContamination : correlated
488 errsyst=TMath::Abs(fhDplus->GetBinContent(j)*fSystDplus->GetHistoSecContaminationMin()->GetBinContent(j));
489 systSecContaminationMin+=errsyst*weight;
490 errsyst=fhDplus->GetBinContent(j)*fSystDplus->GetHistoSecContaminationMax()->GetBinContent(j);
491 systSecContaminationMax+=errsyst*weight;
492
493 // MC closure
494 errsyst=TMath::Abs(fhDplus->GetBinContent(j)*fSystDplus->GetHistoMCclosureTestMin()->GetBinContent(j));
495 systMCclosureMin+=errsyst*weight;
496 errsyst=fhDplus->GetBinContent(j)*fSystDplus->GetHistoMCclosureTestMax()->GetBinContent(j);
497 systMCclosureMax+=errsyst*weight;
498
499 printf(" Dplus Mcclosure the min syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDplus->GetHistoMCclosureTestMin()->GetBinContent(j), systMCclosureMin);
500 printf(" Dplus Mcclosure the max syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDplus->GetHistoMCclosureTestMax()->GetBinContent(j), systMCclosureMax);
501
502 // FD (assume full correlations)
503 errsyst=TMath::Abs(fhDplus->GetBinContent(j)*fSystDplus->GetHistoFDmin()->GetBinContent(j));
504 systFDmin+=errsyst*weight;
505 errsyst=fhDplus->GetBinContent(j)*fSystDplus->GetHistoFDmax()->GetBinContent(j);
506 systFDmax+=errsyst*weight;
507
508 printf(" Dplus feeddown the min syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDplus->GetHistoFDmin()->GetBinContent(j), systMCclosureMin);
509 printf(" Dplus feeddown the max syst value is: %f (rel %f), sum: %f \n",errsyst,fSystDplus->GetHistoFDmax()->GetBinContent(j), systMCclosureMax);
510
511 printf("Dplus the value is: %f, weight: %f \n",value, weight);
512
513 }
514
515 value/=sumweights;
acd5f443 516 errStatValue/=sumweights*sumweights;// was: ((Double_t)nmeson*sumweights);
517 errSystYieldValue/=sumweights*sumweights;// was: ((Double_t)nmeson*sumweights);
518 errSystBkgValue/=sumweights*sumweights; // was: ((Double_t)nmeson*sumweights);
ec9cbe08 519 fhDaverage->SetBinContent(j,value);
520 fhDaverage->SetBinError(j,TMath::Sqrt(errStatValue));
521
522 // Settting unc on yield and back sub
523 fSystDaverage->GetHistoYieldUnc()->SetBinContent(j,TMath::Sqrt(errSystYieldValue)/TMath::Abs(value));
524 fSystDaverage->GetHistoBackSubUncMin()->SetBinContent(j,-TMath::Sqrt(errSystBkgValue)/TMath::Abs(value));
525 fSystDaverage->GetHistoBackSubUncMax()->SetBinContent(j,TMath::Sqrt(errSystBkgValue)/TMath::Abs(value));
526
527
528
529 // MCcorrections, associated tracks
530 if(systMCcorrectionsMin>=0.){
531 systMCcorrectionsMin=systMCcorrectionsMin/sumweights;
532 }
533 else {
534 Printf("Negative sum in error calculation \n");
535 return;
536 }
537
538 if(systMCcorrectionsMax>=0.){
539 systMCcorrectionsMax=systMCcorrectionsMax/sumweights;
540 }
541 else {
542 Printf("Negative sum in error calculation \n");
543 return;
544 }
545
546
547 // MCDefficiency
548 if(systMCDefficiencyMin>=0.){
549 systMCDefficiencyMin=TMath::Sqrt(systMCDefficiencyMin)/sumweights;
550 }
551 else {
552 Printf("Negative sum in error calculation \n");
553 return;
554 }
555
556 if(systMCDefficiencyMax>=0.){
557 systMCDefficiencyMax=TMath::Sqrt(systMCDefficiencyMax)/sumweights;
558 }
559 else {
560 Printf("Negative sum in error calculation \n");
561 return;
562 }
563
564 // SecContamination
565 if(systSecContaminationMin>=0.){
566 systSecContaminationMin=systSecContaminationMin/sumweights;
567 }
568 else {
569 Printf("Negative sum in error calculation \n");
570 return;
571 }
572
573 if(systSecContaminationMax>=0.){
574 systSecContaminationMax=systSecContaminationMax/sumweights;
575 }
576 else {
577 Printf("Negative sum in error calculation \n");
578 return;
579 }
580
581 // MCclosureTest
582 if(systMCclosureMin>=0.){
583 systMCclosureMin=systMCclosureMin/sumweights;
584 }
585 else {
586 Printf("Negative sum in error calculation \n");
587 return;
588 }
589
590 if(systMCclosureMax>=0.){
591 systMCclosureMax=systMCclosureMax/sumweights;
592 }
593 else {
594 Printf("Negative sum in error calculation \n");
595 return;
596 }
597
598
599
600 // Feed down
601 if(systFDmin>=0.){
602 systFDmin=systFDmin/sumweights;
603 }
604 else {
605 Printf("Negative sum in error calculation \n");
606 return;
607 }
608
609 if(systFDmax>=0.){
610 systFDmax=systFDmax/sumweights;
611 }
612 else {
613 Printf("Negative sum in error calculation \n");
614 return;
615 }
616
617 value=TMath::Abs(value);// protection to avoid flipping of min syst histo sign
618
619 fSystDaverage->GetHistoMCcorrectionsMin()->SetBinContent(j,-systMCcorrectionsMin/value);
620 fSystDaverage->GetHistoMCcorrectionsMax()->SetBinContent(j,systMCcorrectionsMax/value);
621
622 fSystDaverage->GetHistoMCDefficiencyMin()->SetBinContent(j,-systMCDefficiencyMin/value);
623 fSystDaverage->GetHistoMCDefficiencyMax()->SetBinContent(j,systMCDefficiencyMax/value);
624
625 fSystDaverage->GetHistoSecContaminationMin()->SetBinContent(j,-systSecContaminationMin/value);
626 fSystDaverage->GetHistoSecContaminationMax()->SetBinContent(j,systSecContaminationMax/value);
627
628 fSystDaverage->GetHistoMCclosureTestMin()->SetBinContent(j,-systMCclosureMin/value);
629 fSystDaverage->GetHistoMCclosureTestMax()->SetBinContent(j,systMCclosureMax/value);
630
631 fSystDaverage->GetHistoFDmin()->SetBinContent(j,-systFDmin/value);
632 fSystDaverage->GetHistoFDmax()->SetBinContent(j,systFDmax/value);
633
634 }
635 fSystDaverage->BuildTotalUncHisto();
636 fSystDaverage->BuildTotalNonFDUncHisto();
637 fSystDaverage->BuildGraphsUnc(fhDaverage);
638 fSystDaverage->BuildGraphsRelUnc();
639
640 return ;
641
642}
643
644
645TH1D* AliHFDmesonCorrAverage::ReflectHisto(TH1D *h){
646
647 TH1D *h2=new TH1D(Form("%sReflected",h->GetName()),Form("%sReflected",h->GetName()),h->GetNbinsX()/2.,0.,TMath::Pi());
648 for(Int_t j=1;j<=h->GetNbinsX();j++){
649 Double_t x=h->GetBinCenter(j);
650 Double_t y0=h->GetBinContent(j);
651 Double_t ey0=h->GetBinError(j);
652 Int_t j2;
653 if(x>0&&x<TMath::Pi())j2=h2->FindBin(x);
654 else if(x<0)j2=h2->FindBin(-1.*x);
655 else if(x>TMath::Pi())j2=h2->FindBin(2.*TMath::Pi()-x);
656 else printf("Point %d excluded \n",j);
657 Double_t y=h2->GetBinContent(j2);
658 Double_t ey=h2->GetBinError(j2);
659 h2->SetBinContent(j2,y+y0);
660 h2->SetBinError(j2,TMath::Sqrt(ey0*ey0+ey*ey));
661 }
662
663 return h2;
664}
665
666void AliHFDmesonCorrAverage::SetWeights(){
667
668
669 if(fincludeDzero){
670 TH1D *hYieldUnc=fSystDzero->GetHistoYieldUnc();
671 TH1D *hBkgUncMin=fSystDzero->GetHistoBackSubUncMin();
672 TH1D *hBkgUncMax=fSystDzero->GetHistoBackSubUncMax();
673 if(fweightsDzeroStat)delete fweightsDzeroStat;
674 if(fweightsDzeroSystYield)delete fweightsDzeroSystYield;
675 if(fweightsDzeroSystBkg)delete fweightsDzeroSystBkg;
676 fweightsDzeroStat=new Double_t[fnbinsphi];
677 fweightsDzeroSystYield=new Double_t[fnbinsphi];
678 fweightsDzeroSystBkg=new Double_t[fnbinsphi];
679 for(Int_t j=0;j<fnbinsphi;j++){
680 if(fmethod==10){
681 fweightsDzeroStat[j]=1./(fhDzero->GetBinError(j+1)*fhDzero->GetBinError(j+1));
682 fweightsDzeroSystYield[j]=1./(hYieldUnc->GetBinContent(j+1)*fhDzero->GetBinContent(j+1)*hYieldUnc->GetBinContent(j+1)*fhDzero->GetBinContent(j+1));
683 //for asymmetric uncertainties...
684 if(TMath::Abs(hBkgUncMin->GetBinContent(j+1)) > TMath::Abs(hBkgUncMax->GetBinContent(j+1))) fweightsDzeroSystBkg[j]=1./(hBkgUncMin->GetBinContent(j+1)*fhDzero->GetBinContent(j+1)*hBkgUncMin->GetBinContent(j+1)*fhDzero->GetBinContent(j+1));
685 else fweightsDzeroSystBkg[j]=1./(hBkgUncMax->GetBinContent(j+1)*fhDzero->GetBinContent(j+1)*hBkgUncMax->GetBinContent(j+1)*fhDzero->GetBinContent(j+1));
686 }
687 else{
688 Printf("This option for the calculation of the average has not been implemented yet");
689 return;
690 }
691 }
692 }
693
694 if(fincludeDstar){
695 TH1D *hYieldUnc=fSystDstar->GetHistoYieldUnc();
696 TH1D *hBkgUncMin=fSystDstar->GetHistoBackSubUncMin();
697 TH1D *hBkgUncMax=fSystDstar->GetHistoBackSubUncMax();
698 if(fweightsDstarStat)delete fweightsDstarStat;
699 if(fweightsDstarSystYield)delete fweightsDstarSystYield;
700 if(fweightsDstarSystBkg)delete fweightsDstarSystBkg;
701 fweightsDstarStat=new Double_t[fnbinsphi];
702 fweightsDstarSystYield=new Double_t[fnbinsphi];
703 fweightsDstarSystBkg=new Double_t[fnbinsphi];
704 for(Int_t j=0;j<fnbinsphi;j++){
705 if(fmethod==10){
706 fweightsDstarStat[j]=1./(fhDstar->GetBinError(j+1)*fhDstar->GetBinError(j+1));
707 fweightsDstarSystYield[j]=1./(hYieldUnc->GetBinContent(j+1)*fhDstar->GetBinContent(j+1)*hYieldUnc->GetBinContent(j+1)*fhDstar->GetBinContent(j+1));
708 //for asymmetric uncertainties...
709 if(TMath::Abs(hBkgUncMin->GetBinContent(j+1)) > TMath::Abs(hBkgUncMax->GetBinContent(j+1))) fweightsDstarSystBkg[j]=1./(hBkgUncMin->GetBinContent(j+1)*fhDstar->GetBinContent(j+1)*hBkgUncMin->GetBinContent(j+1)*fhDstar->GetBinContent(j+1));
710 else fweightsDstarSystBkg[j]=1./(hBkgUncMax->GetBinContent(j+1)*fhDstar->GetBinContent(j+1)*hBkgUncMax->GetBinContent(j+1)*fhDstar->GetBinContent(j+1));
711 }
712 else{
713 Printf("This option for the calculation of the average has not been implemented yet");
714 return;
715 }
716 }
717 }
718
719 if(fincludeDplus){
720 TH1D *hYieldUnc=fSystDplus->GetHistoYieldUnc();
721 TH1D *hBkgUncMin=fSystDplus->GetHistoBackSubUncMin();
722 TH1D *hBkgUncMax=fSystDplus->GetHistoBackSubUncMax();
723 if(fweightsDplusStat)delete fweightsDplusStat;
724 if(fweightsDplusSystYield)delete fweightsDplusSystYield;
725 if(fweightsDplusSystBkg)delete fweightsDplusSystBkg;
726 fweightsDplusStat=new Double_t[fnbinsphi];
727 fweightsDplusSystYield=new Double_t[fnbinsphi];
728 fweightsDplusSystBkg=new Double_t[fnbinsphi];
729 for(Int_t j=0;j<fnbinsphi;j++){
730 if(fmethod==10){
731 fweightsDplusStat[j]=1./(fhDplus->GetBinError(j+1)*fhDplus->GetBinError(j+1));
732 fweightsDplusSystYield[j]=1./(hYieldUnc->GetBinContent(j+1)*fhDplus->GetBinContent(j+1)*hYieldUnc->GetBinContent(j+1)*fhDplus->GetBinContent(j+1));
733 //for asymmetric uncertainties...
734 if(TMath::Abs(hBkgUncMin->GetBinContent(j+1)) > TMath::Abs(hBkgUncMax->GetBinContent(j+1))) fweightsDplusSystBkg[j]=1./(hBkgUncMin->GetBinContent(j+1)*fhDplus->GetBinContent(j+1)*hBkgUncMin->GetBinContent(j+1)*fhDplus->GetBinContent(j+1));
735 else fweightsDplusSystBkg[j]=1./(hBkgUncMax->GetBinContent(j+1)*fhDplus->GetBinContent(j+1)*hBkgUncMax->GetBinContent(j+1)*fhDplus->GetBinContent(j+1));
736 }
737 else{
738 Printf("This option for the calculation of the average has not been implemented yet");
739 return;
740 }
741 }
742 }
743
744
745
746
747}
748
749