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