1 /**************************************************************************
2 * Copyright(c) 1998-2010, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 /////////////////////////////////////////////////////////////
18 // Class to handle systematic errors for charm hadrons
21 // AliHFSystEff syst(DECAY); // DECAY = 1 for D0, 2, for D+, 3 for D*
22 // syst.DrawErrors(); // to see a plot of the error contributions
23 // syst.GetTotalSystErr(pt); // to get the total err at pt
25 // Author: A.Dainese, andrea.dainese@pd.infn.it
26 /////////////////////////////////////////////////////////////
29 #include <TGraphAsymmErrors.h>
36 #include "AliHFSystErr.h"
39 ClassImp(AliHFSystErr)
41 //--------------------------------------------------------------------------
42 AliHFSystErr::AliHFSystErr(const Char_t* name, const Char_t* title) :
54 // Default Constructor
57 //--------------------------------------------------------------------------
58 AliHFSystErr::AliHFSystErr(Int_t decay,const Char_t* name, const Char_t* title) :
70 // Standard Constructor
84 printf("Invalid decay type: %d\n",decay);
90 //--------------------------------------------------------------------------
91 AliHFSystErr::~AliHFSystErr() {
96 if(fNorm) { delete fNorm; fNorm=0; }
97 if(fRawYield) { delete fRawYield; fRawYield=0; }
98 if(fTrackingEff) { delete fTrackingEff; fTrackingEff=0; }
99 if(fBR) { delete fBR; fBR=0; }
100 if(fCutsEff) { delete fCutsEff; fCutsEff=0; }
101 if(fPIDEff) { delete fPIDEff; fPIDEff=0; }
102 if(fMCPtShape) { delete fMCPtShape; fMCPtShape=0; }
103 if(fPartAntipart) { delete fPartAntipart; fPartAntipart=0; }
106 //--------------------------------------------------------------------------
107 void AliHFSystErr::InitD0toKpi() {
109 // D0->Kpi syst errors. Responsible: A. Rossi
113 fNorm = new TH1F("fNorm","fNorm",20,0,20);
114 for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
117 fBR = new TH1F("fBR","fBR",20,0,20);
118 for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.012); // 1.2% PDG2010
120 // Tracking efficiency
121 fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
122 for(Int_t i=1;i<=20;i++) fTrackingEff->SetBinContent(i,0.02); // 2% (1% per track)
124 // Raw yield extraction
125 fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
126 fRawYield->SetBinContent(1,1);
127 fRawYield->SetBinContent(2,1);
128 fRawYield->SetBinContent(3,0.15);
129 for(Int_t i=4;i<=20;i++) fRawYield->SetBinContent(i,0.065);
131 // Cuts efficiency (from cuts variation)
132 fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
133 for(Int_t i=1;i<=20;i++) fCutsEff->SetBinContent(i,0.10); // 10%
135 // PID efficiency (from PID/noPID)
136 fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
137 for(Int_t i=1;i<=20;i++) fPIDEff->SetBinContent(i,0.03); // 3%
138 fPIDEff->SetBinContent(4,0.10); // 10%
141 fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
142 for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,(Float_t)i*0.006);
144 // particle-antiparticle
145 fPartAntipart = new TH1F("fPartAntipart","fPartAntipart",20,0,20);
146 fPartAntipart->SetBinContent(1,1);
147 fPartAntipart->SetBinContent(2,1);
148 for(Int_t i=3;i<=6;i++) fPartAntipart->SetBinContent(i,0.08);
152 //--------------------------------------------------------------------------
153 void AliHFSystErr::InitDplustoKpipi() {
155 // D+->Kpipi syst errors. Responsible: R. Bala
159 fNorm = new TH1F("fNorm","fNorm",20,0,20);
160 for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
163 fBR = new TH1F("fBR","fBR",20,0,20);
164 for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.04); // 4% PDG2010
166 // Tracking efficiency
167 fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
168 for(Int_t i=1;i<=20;i++) fTrackingEff->SetBinContent(i,0.03); // 3% (1% per track)
171 // Raw yield extraction
172 fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
173 fRawYield->SetBinContent(1,1);
174 fRawYield->SetBinContent(2,1);
175 fRawYield->SetBinContent(3,0.20);
176 for(Int_t i=4;i<=20;i++) fRawYield->SetBinContent(i,0.055); //5 to 10%
178 // Cuts efficiency (from cuts variation)
179 fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
180 for(Int_t i=1;i<=20;i++) fCutsEff->SetBinContent(i,0.10); // 10%
182 // PID efficiency (from PID/noPID)
183 fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
184 for(Int_t i=1;i<=20;i++) fPIDEff->SetBinContent(i,0.05); // 5%
185 fPIDEff->SetBinContent(3,0.13); // 13%
188 // MC dN/dpt (copied from D0 : will update later)
189 fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
190 for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,(Float_t)i*0.006);
193 // particle-antiparticle
195 fPartAntipart = new TH1F("fPartAntipart","fPartAntipart",20,0,20);
196 fPartAntipart->SetBinContent(1,1);
197 fPartAntipart->SetBinContent(2,1);
198 fPartAntipart->SetBinContent(3,0.12);
199 for(Int_t i=4;i<=20;i++) fPartAntipart->SetBinContent(i,0.05); //5 to 12%
203 //--------------------------------------------------------------------------
204 void AliHFSystErr::InitDstartoD0pi() {
206 // D*+->D0pi syst errors. Responsible: A. Grelli, Y. Wang
210 fNorm = new TH1F("fNorm","fNorm",20,0,20);
211 for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
214 fBR = new TH1F("fBR","fBR",20,0,20);
215 for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.07); // 7% PDG2010
217 // Tracking efficiency
218 fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
219 fTrackingEff->SetBinContent(1,0.12);
220 fTrackingEff->SetBinContent(2,0.08);
221 fTrackingEff->SetBinContent(3,0.032);
222 for(Int_t i=4;i<=20;i++) fTrackingEff->SetBinContent(i,0.03); // 3% (1% per track)
225 // Raw yield extraction
226 fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
227 fRawYield->SetBinContent(1,0.19);
228 fRawYield->SetBinContent(2,0.10);
229 fRawYield->SetBinContent(3,0.10);
230 fRawYield->SetBinContent(4,0.13);
231 fRawYield->SetBinContent(5,0.09);
232 for(Int_t i=5;i<=20;i++) fRawYield->SetBinContent(i,0.08); //8%
234 // Cuts efficiency (from cuts variation)
235 fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
236 for(Int_t i=1;i<=20;i++) fCutsEff->SetBinContent(i,0.10); // 10%
238 // PID efficiency (from PID/noPID)
239 fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
240 for(Int_t i=1;i<=20;i++) fPIDEff->SetBinContent(i,0.03); // 3%
241 fPIDEff->SetBinContent(4,0.10); // 10%
242 fPIDEff->SetBinContent(5,0.06); // 6%
245 // MC dN/dpt (copied from D0 : will update later)
246 fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
247 for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,(Float_t)i*0.006);
251 //--------------------------------------------------------------------------
252 Double_t AliHFSystErr::GetCutsEffErr(Double_t pt) const {
257 Int_t bin=fCutsEff->FindBin(pt);
259 return fCutsEff->GetBinContent(bin);
261 //--------------------------------------------------------------------------
262 Double_t AliHFSystErr::GetMCPtShapeErr(Double_t pt) const {
267 Int_t bin=fMCPtShape->FindBin(pt);
269 return fMCPtShape->GetBinContent(bin);
271 //--------------------------------------------------------------------------
272 Double_t AliHFSystErr::GetSeleEffErr(Double_t pt) const {
277 Double_t err=GetCutsEffErr(pt)*GetCutsEffErr(pt)+GetMCPtShapeErr(pt)*GetMCPtShapeErr(pt);
279 return TMath::Sqrt(err);
281 //--------------------------------------------------------------------------
282 Double_t AliHFSystErr::GetPIDEffErr(Double_t pt) const {
287 Int_t bin=fPIDEff->FindBin(pt);
289 return fPIDEff->GetBinContent(bin);
291 //--------------------------------------------------------------------------
292 Double_t AliHFSystErr::GetTrackingEffErr(Double_t pt) const {
297 Int_t bin=fTrackingEff->FindBin(pt);
299 return fTrackingEff->GetBinContent(bin);
301 //--------------------------------------------------------------------------
302 Double_t AliHFSystErr::GetRawYieldErr(Double_t pt) const {
307 Int_t bin=fRawYield->FindBin(pt);
309 return fRawYield->GetBinContent(bin);
311 //--------------------------------------------------------------------------
312 Double_t AliHFSystErr::GetPartAntipartErr(Double_t pt) const {
317 Int_t bin=fPartAntipart->FindBin(pt);
319 return fPartAntipart->GetBinContent(bin);
321 //--------------------------------------------------------------------------
322 Double_t AliHFSystErr::GetTotalSystErr(Double_t pt,Double_t feeddownErr) const {
324 // Get total syst error (except norm. error)
329 if(fRawYield) err += GetRawYieldErr(pt)*GetRawYieldErr(pt);
330 if(fTrackingEff) err += GetTrackingEffErr(pt)*GetTrackingEffErr(pt);
331 if(fBR) err += GetBRErr()*GetBRErr();
332 if(fCutsEff) err += GetCutsEffErr(pt)*GetCutsEffErr(pt);
333 if(fPIDEff) err += GetPIDEffErr(pt)*GetPIDEffErr(pt);
334 if(fMCPtShape) err += GetMCPtShapeErr(pt)*GetMCPtShapeErr(pt);
335 if(fPartAntipart) err += GetPartAntipartErr(pt)*GetPartAntipartErr(pt);
337 err += feeddownErr*feeddownErr;
339 return TMath::Sqrt(err);
341 //---------------------------------------------------------------------------
342 void AliHFSystErr::DrawErrors(TGraphAsymmErrors *grErrFeeddown) const {
346 gStyle->SetOptStat(0);
348 TCanvas *cSystErr = new TCanvas("cSystErr","Systematic Errors",300,80,640,500);
349 cSystErr->Range(0.20,-0.5,18.4,0.34);
350 cSystErr->SetRightMargin(0.318);
351 cSystErr->SetFillColor(0);
353 TH2F *hFrame = new TH2F("hFrame","Systematic errors; p_{t} [GeV/c]; Relative Error",20,0,20,100,-1,+1);
354 hFrame->SetAxisRange(2.,11.9,"X");
355 hFrame->SetAxisRange(-0.5,0.5,"Y");
358 TLegend *leg = new TLegend(0.69,0.44,0.98,0.86,NULL,"brNDC");
359 leg->SetTextSize(0.03601695);
360 leg->SetFillStyle(0);
361 leg->SetBorderSize(0);
363 TH1F *hTotErr=new TH1F("hTotErr","",20,0,20);
364 Int_t nbins = fNorm->GetNbinsX();
365 TGraphAsymmErrors *gTotErr = new TGraphAsymmErrors(nbins);
366 for(Int_t i=1;i<=20;i++) {
367 Double_t pt = hTotErr->GetBinCenter(i);
368 Double_t ptwidth = hTotErr->GetBinWidth(i);
371 Double_t x=0., y=0., errxl=0., errxh=0., erryl=0., erryh=0.;
372 Double_t toterryl=0., toterryh=0.;
373 for(Int_t j=0; j<grErrFeeddown->GetN(); j++) {
374 grErrFeeddown->GetPoint(j,x,y);
375 errxh = grErrFeeddown->GetErrorXhigh(j);
376 errxl = grErrFeeddown->GetErrorXlow(j);
377 if ( ( (x-errxl) <= pt) && ( (x+errxl) >= pt) ) {
378 erryh = grErrFeeddown->GetErrorYhigh(j);
379 erryl = grErrFeeddown->GetErrorYlow(j);
382 if (erryl>=1e-3) toterryl = GetTotalSystErr(pt,erryl);
383 else toterryl = GetTotalSystErr(pt);
384 if (erryh>=1e-3) toterryh = GetTotalSystErr(pt,erryh);
385 else toterryh = GetTotalSystErr(pt);
387 hTotErr->SetBinContent(i,toterryh);
388 gTotErr->SetPoint(i,pt,0.);
389 gTotErr->SetPointError(i,ptwidth/2.,ptwidth/2.,toterryl,toterryh); // i, exl, exh, eyl, eyh
392 hTotErr->SetBinContent(i,GetTotalSystErr(pt));
393 gTotErr->SetPoint(i,pt,0.);
394 gTotErr->SetPointError(i,ptwidth/2.,ptwidth/2.,GetTotalSystErr(pt),GetTotalSystErr(pt)); // i, exl, exh, eyl, eyh
398 gTotErr->SetLineColor(kBlack);
399 gTotErr->SetFillColor(kRed);
400 gTotErr->SetFillStyle(3002);
402 leg->AddEntry(gTotErr,"Total (excl. norm.)","f");
403 // hTotErr->SetLineColor(1);
404 // hTotErr->SetLineWidth(3);
405 // hTotErr->Draw("same");
406 // leg->AddEntry(hTotErr,"Total (excl. norm.)","l");
410 fNorm->SetFillColor(1);
411 fNorm->SetFillStyle(3002);
412 //fNorm->Draw("same");
413 //TH1F *hNormRefl = ReflectHisto(fNorm);
414 //hNormRefl->Draw("same");
415 leg->AddEntry(fNorm,"Normalization (10%)","");
418 grErrFeeddown->SetFillColor(kTeal-8);
419 grErrFeeddown->SetFillStyle(3001);
420 grErrFeeddown->Draw("2");
421 leg->AddEntry(grErrFeeddown,"Feed-down from B","f");
424 fTrackingEff->SetFillColor(4);
425 fTrackingEff->SetFillStyle(3006);
426 fTrackingEff->Draw("same");
427 TH1F *hTrackingEffRefl = ReflectHisto(fTrackingEff);
428 hTrackingEffRefl->Draw("same");
429 leg->AddEntry(fTrackingEff,"Tracking efficiency","f");
432 fBR->SetFillColor(6);
433 fBR->SetFillStyle(3005);
434 //fBR->SetFillStyle(3020);
436 TH1F *hBRRefl = ReflectHisto(fBR);
437 hBRRefl->Draw("same");
438 leg->AddEntry(fBR,"Branching ratio","f");
441 Int_t ci; // for color index setting
442 ci = TColor::GetColor("#00cc00");
443 fRawYield->SetLineColor(ci);
444 // fRawYield->SetLineColor(3);
445 fRawYield->SetLineWidth(3);
446 fRawYield->Draw("same");
447 TH1F *hRawYieldRefl = ReflectHisto(fRawYield);
448 hRawYieldRefl->Draw("same");
449 leg->AddEntry(fRawYield,"Yield extraction","l");
452 fCutsEff->SetLineColor(4);
453 fCutsEff->SetLineWidth(3);
454 fCutsEff->Draw("same");
455 TH1F *hCutsEffRefl = ReflectHisto(fCutsEff);
456 hCutsEffRefl->Draw("same");
457 leg->AddEntry(fCutsEff,"Cuts efficiency","l");
460 fPIDEff->SetLineColor(7);
461 fPIDEff->SetLineWidth(3);
462 fPIDEff->Draw("same");
463 TH1F *hPIDEffRefl = ReflectHisto(fPIDEff);
464 hPIDEffRefl->Draw("same");
465 leg->AddEntry(fPIDEff,"PID efficiency","l");
468 Int_t ci = TColor::GetColor("#9933ff");
469 fMCPtShape->SetLineColor(ci);
470 // fMCPtShape->SetLineColor(8);
471 fMCPtShape->SetLineWidth(3);
472 fMCPtShape->Draw("same");
473 TH1F *hMCPtShapeRefl = ReflectHisto(fMCPtShape);
474 hMCPtShapeRefl->Draw("same");
475 leg->AddEntry(fMCPtShape,"MC p_{t} shape","l");
478 Int_t ci = TColor::GetColor("#ff6600");
479 fPartAntipart->SetLineColor(ci);
480 // fPartAntipart->SetLineColor(9);
481 fPartAntipart->SetLineWidth(3);
482 fPartAntipart->Draw("same");
483 TH1F *hPartAntipartRefl = ReflectHisto(fPartAntipart);
484 hPartAntipartRefl->Draw("same");
485 leg->AddEntry(fPartAntipart,"D = #bar{D}","l");
491 cSystErr->SaveAs("RelativeSystematics.eps");
495 //-------------------------------------------------------------------------
496 TH1F* AliHFSystErr::ReflectHisto(TH1F *hin) const {
498 // Clones and reflects histogram
500 TH1F *hout=(TH1F*)hin->Clone("hout");