Changes for PbPb (Zaida)
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliHFSystErr.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2010, 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 /* $Id$ */
17
18 /////////////////////////////////////////////////////////////
19 //
20 // Class to handle systematic errors for charm hadrons
21 //
22 // Usage:
23 // AliHFSystEff syst;           // DECAY = 1 for D0, 2, for D+, 3 for D* 
24 // syst.SetRunNumber(YEAR);     // YEAR = two last numbers of the year (is 10 for 2010)
25 // syst.SetCollisionType(TYPE);  // TYPE =  0 is pp, 1 is PbPb
26 // syst.SetCentrality(CENT);     // CENT is centrality, 0100 for MB, 020 (4080) for 0-20 (40-80) CC...
27 // syst.Init(DECAY);             // DECAY = 1 for D0, 2, for D+, 3 for D* 
28 // syst.DrawErrors(); // to see a plot of the error contributions
29 // syst.GetTotalSystErr(pt); // to get the total err at pt 
30 //
31 // Author: A.Dainese, andrea.dainese@pd.infn.it
32 /////////////////////////////////////////////////////////////
33
34 #include <TStyle.h>
35 #include <TGraphAsymmErrors.h>
36 #include <TMath.h>
37 #include <TCanvas.h>
38 #include <TH2F.h>
39 #include <TLegend.h>
40 #include <TColor.h>
41
42 #include "AliLog.h"
43 #include "AliHFSystErr.h"
44
45
46 ClassImp(AliHFSystErr)
47
48 //--------------------------------------------------------------------------
49 AliHFSystErr::AliHFSystErr(const Char_t* name, const Char_t* title) : 
50 TNamed(name,title),
51 fNorm(0),
52 fRawYield(0),
53 fTrackingEff(0),
54 fBR(0),
55 fCutsEff(0),
56 fPIDEff(0),
57 fMCPtShape(0),
58 fPartAntipart(0),
59 fRunNumber(10),
60 fCollisionType(0),
61 fCentralityClass(0100)
62 {
63   //
64   // Default Constructor
65   //
66 }
67
68 //--------------------------------------------------------------------------
69 AliHFSystErr::~AliHFSystErr() {
70   //  
71   // Default Destructor
72   //
73
74   if(fNorm)         { delete fNorm; fNorm=0; }
75   if(fRawYield)     { delete fRawYield; fRawYield=0; }
76   if(fTrackingEff)  { delete fTrackingEff; fTrackingEff=0; }
77   if(fBR)           { delete fBR; fBR=0; }
78   if(fCutsEff)      { delete fCutsEff; fCutsEff=0; }
79   if(fPIDEff)       { delete fPIDEff; fPIDEff=0; }
80   if(fMCPtShape)    { delete fMCPtShape; fMCPtShape=0; }
81   if(fPartAntipart) { delete fPartAntipart; fPartAntipart=0; }
82
83 }
84
85 //--------------------------------------------------------------------------
86 void AliHFSystErr::Init(Int_t decay){
87   //
88   // Variables/histos initialization
89   //
90
91   if (fRunNumber!=10) { 
92     AliError("Only settings for 2010 are implemented so far");
93   }
94   if (fCentralityClass!=020 && fCentralityClass!=4080 && fCentralityClass!=0100){
95     AliError("Only settings for MB2010 are implemented so far");
96   }
97
98   switch(decay) {
99   case 1: // D0->Kpi
100     if (fCollisionType==0) InitD0toKpi2010pp();
101     else if (fCollisionType==1) {
102       if (fCentralityClass==020) InitD0toKpi2010PbPb020();
103       else if (fCentralityClass==4080) InitD0toKpi2010PbPb4080();
104       else AliError("Not yet implemented");
105     }
106     break;
107   case 2: // D+->Kpipi
108     if (fCollisionType==0) InitDplustoKpipi2010pp();
109     else if (fCollisionType==1) {
110       if (fCentralityClass==020) InitDplustoKpipi2010PbPb020();
111       else if (fCentralityClass==4080) InitDplustoKpipi2010PbPb4080();
112       else AliError("Not yet implemented");
113     }
114     break;
115   case 3: // D*->D0pi
116     if (fCollisionType==0) InitDstartoD0pi2010pp();
117     else if (fCollisionType==1) {
118       if (fCentralityClass==020) InitDstartoD0pi2010PbPb020();
119       else if (fCentralityClass==4080) InitDstartoD0pi2010PbPb4080();
120       else AliError("Not yet implemented");
121     }
122     break;
123   default:
124     printf("Invalid decay type: %d\n",decay);
125     break;
126   }
127
128 }
129   
130 //--------------------------------------------------------------------------
131 void AliHFSystErr::InitD0toKpi2010pp() {
132   // 
133   // D0->Kpi syst errors. Responsible: A. Rossi
134   //   2010 pp sample
135   //
136
137   // Normalization
138   fNorm = new TH1F("fNorm","fNorm",20,0,20);
139   for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
140
141   // Branching ratio 
142   fBR = new TH1F("fBR","fBR",20,0,20);
143   for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.012); // 1.2% PDG2010
144
145   // Tracking efficiency
146   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
147   for(Int_t i=1;i<=20;i++) fTrackingEff->SetBinContent(i,0.02); // 2% (1% per track)
148
149   // Raw yield extraction
150   fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
151   fRawYield->SetBinContent(1,1);
152   fRawYield->SetBinContent(2,1);
153   fRawYield->SetBinContent(3,0.15);
154   for(Int_t i=4;i<=20;i++) fRawYield->SetBinContent(i,0.065);
155
156   // Cuts efficiency (from cuts variation)
157   fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
158   for(Int_t i=1;i<=20;i++) fCutsEff->SetBinContent(i,0.10); // 10%
159
160   // PID efficiency (from PID/noPID)
161   fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
162   for(Int_t i=1;i<=20;i++) fPIDEff->SetBinContent(i,0.03); // 3%
163   fPIDEff->SetBinContent(4,0.10); // 10%
164
165   // MC dN/dpt
166   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
167   for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,(Float_t)i*0.006);
168
169   // particle-antiparticle
170   fPartAntipart = new TH1F("fPartAntipart","fPartAntipart",20,0,20);
171   fPartAntipart->SetBinContent(1,1);
172   fPartAntipart->SetBinContent(2,1);
173   for(Int_t i=3;i<=6;i++) fPartAntipart->SetBinContent(i,0.08);
174   
175   return;
176 }
177 //--------------------------------------------------------------------------
178 void AliHFSystErr::InitD0toKpi2010PbPb020() {
179   // 
180   // D0->Kpi syst errors. Responsible: ??
181   //   2010 PbPb sample, 0-20 CC
182   //
183
184   // Normalization
185   fNorm = new TH1F("fNorm","fNorm",20,0,20);
186   for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
187
188   // Branching ratio 
189   fBR = new TH1F("fBR","fBR",20,0,20);
190   for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.012); // 1.2% PDG2010
191
192   // Tracking efficiency
193   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
194   for(Int_t i=1;i<=20;i++) fTrackingEff->SetBinContent(i,0.02); // 2% (1% per track)
195
196   // Raw yield extraction
197   fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
198   fRawYield->SetBinContent(1,0);
199   fRawYield->SetBinContent(2,0);
200   fRawYield->SetBinContent(3,0);
201   fRawYield->SetBinContent(4,0.27);
202   fRawYield->SetBinContent(5,0.27);
203   fRawYield->SetBinContent(6,0);
204   fRawYield->SetBinContent(7,0.07);
205   fRawYield->SetBinContent(8,0.07);
206   for(Int_t i=9;i<=12;i++) fRawYield->SetBinContent(i,0.09);
207   for(Int_t i=12;i<=20;i++) fRawYield->SetBinContent(i,0);
208
209   // Cuts efficiency (from cuts variation)
210   fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
211   for(Int_t i=1;i<=20;i++) fCutsEff->SetBinContent(i,0.10); // 10%
212
213   // PID efficiency (from PID/noPID)
214   fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
215   for(Int_t i=1;i<=20;i++) fPIDEff->SetBinContent(i,0.03); // 3%
216   fPIDEff->SetBinContent(4,0.10); // 10%
217
218   // MC dN/dpt
219   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
220   for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,(Float_t)i*0.006);
221
222   // particle-antiparticle
223   fPartAntipart = new TH1F("fPartAntipart","fPartAntipart",20,0,20);
224   fPartAntipart->SetBinContent(1,1);
225   fPartAntipart->SetBinContent(2,1);
226   for(Int_t i=3;i<=6;i++) fPartAntipart->SetBinContent(i,0.08);
227   
228   return;
229 }
230 //--------------------------------------------------------------------------
231 void AliHFSystErr::InitD0toKpi2010PbPb4080() {
232   // 
233   // D0->Kpi syst errors. Responsible: ??
234   //   2010 PbPb sample, 40-80 CC
235   //
236
237   // Normalization
238   fNorm = new TH1F("fNorm","fNorm",20,0,20);
239   for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
240
241   // Branching ratio 
242   fBR = new TH1F("fBR","fBR",20,0,20);
243   for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.012); // 1.2% PDG2010
244
245   // Tracking efficiency
246   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
247   for(Int_t i=1;i<=20;i++) fTrackingEff->SetBinContent(i,0.02); // 2% (1% per track)
248
249   // Raw yield extraction
250   fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
251   fRawYield->SetBinContent(1,0);
252   fRawYield->SetBinContent(2,0);
253   fRawYield->SetBinContent(3,0);
254   fRawYield->SetBinContent(4,0);
255   fRawYield->SetBinContent(5,0.3);
256   fRawYield->SetBinContent(6,0);
257   fRawYield->SetBinContent(7,0.33);
258   fRawYield->SetBinContent(8,0.33);
259   for(Int_t i=9;i<=20;i++) fRawYield->SetBinContent(i,0);
260
261   // Cuts efficiency (from cuts variation)
262   fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
263   for(Int_t i=1;i<=20;i++) fCutsEff->SetBinContent(i,0.10); // 10%
264
265   // PID efficiency (from PID/noPID)
266   fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
267   for(Int_t i=1;i<=20;i++) fPIDEff->SetBinContent(i,0.03); // 3%
268   fPIDEff->SetBinContent(4,0.10); // 10%
269
270   // MC dN/dpt
271   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
272   for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,(Float_t)i*0.006);
273
274   // particle-antiparticle
275   fPartAntipart = new TH1F("fPartAntipart","fPartAntipart",20,0,20);
276   fPartAntipart->SetBinContent(1,1);
277   fPartAntipart->SetBinContent(2,1);
278   for(Int_t i=3;i<=6;i++) fPartAntipart->SetBinContent(i,0.08);
279   
280   return;
281 }
282 //--------------------------------------------------------------------------
283 void AliHFSystErr::InitDplustoKpipi2010pp() {
284   // 
285   // D+->Kpipi syst errors. Responsible: R. Bala
286   //  2010 pp sample
287   //
288
289  // Normalization
290   fNorm = new TH1F("fNorm","fNorm",20,0,20);
291   for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
292
293   // Branching ratio 
294   fBR = new TH1F("fBR","fBR",20,0,20);
295   for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.04); // 4% PDG2010
296
297   // Tracking efficiency
298   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
299   for(Int_t i=1;i<=20;i++) fTrackingEff->SetBinContent(i,0.03); // 3% (1% per track)
300
301
302   // Raw yield extraction
303   fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
304   fRawYield->SetBinContent(1,1);
305   fRawYield->SetBinContent(2,1);
306   fRawYield->SetBinContent(3,0.20);
307   for(Int_t i=4;i<=20;i++) fRawYield->SetBinContent(i,0.055);  //5 to 10%
308
309   // Cuts efficiency (from cuts variation)
310   fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
311   for(Int_t i=1;i<=20;i++) fCutsEff->SetBinContent(i,0.10); // 10%
312
313   // PID efficiency (from PID/noPID)
314   fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
315   for(Int_t i=1;i<=20;i++) fPIDEff->SetBinContent(i,0.05); // 5%
316   fPIDEff->SetBinContent(3,0.13); // 13%
317  
318
319   // MC dN/dpt  (copied from D0 : will update later)
320   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
321   for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,(Float_t)i*0.006);
322
323
324   // particle-antiparticle
325   /*
326   fPartAntipart = new TH1F("fPartAntipart","fPartAntipart",20,0,20);
327   fPartAntipart->SetBinContent(1,1);
328   fPartAntipart->SetBinContent(2,1);
329   fPartAntipart->SetBinContent(3,0.12);
330   for(Int_t i=4;i<=20;i++) fPartAntipart->SetBinContent(i,0.05);   //5 to 12%
331   */
332   return;
333 }
334 //--------------------------------------------------------------------------
335 void AliHFSystErr::InitDplustoKpipi2010PbPb020() {
336   // 
337   // D+->Kpipi syst errors. Responsible: ??
338   //  2010 PbPb sample, 0-20 CC
339   //
340   
341   AliInfo("Not yet implemented");
342   return;
343 }
344
345 //--------------------------------------------------------------------------
346 void AliHFSystErr::InitDplustoKpipi2010PbPb4080() {
347   // 
348   // D+->Kpipi syst errors. Responsible: ??
349   //  2010 PbPb sample, 40-80 CC
350   //
351   
352   AliInfo("Not yet implemented");
353   return;
354 }
355
356 //--------------------------------------------------------------------------
357 void AliHFSystErr::InitDstartoD0pi2010pp() {
358   // 
359   // D*+->D0pi syst errors. Responsible: A. Grelli, Y. Wang
360   //  2010 pp sample
361   //
362
363  // Normalization
364   fNorm = new TH1F("fNorm","fNorm",20,0,20);
365   for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
366
367   // Branching ratio 
368   fBR = new TH1F("fBR","fBR",20,0,20);
369   for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.015); // 1.5% PDG2010
370
371   // Tracking efficiency
372   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
373   fTrackingEff->SetBinContent(1,0.12);
374   fTrackingEff->SetBinContent(2,0.08);
375   fTrackingEff->SetBinContent(3,0.05);
376   for(Int_t i=4;i<=20;i++) fTrackingEff->SetBinContent(i,0.03); // 3% (1% per track)
377
378
379   // Raw yield extraction
380   fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
381   fRawYield->SetBinContent(1,1);
382   fRawYield->SetBinContent(2,0.12);
383   fRawYield->SetBinContent(3,0.09);
384   fRawYield->SetBinContent(4,0.08);
385   fRawYield->SetBinContent(5,0.06);
386   for(Int_t i=5;i<=20;i++) fRawYield->SetBinContent(i,0.04);  //4%
387
388   // Cuts efficiency (from cuts variation)
389   fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
390   for(Int_t i=1;i<=20;i++) fCutsEff->SetBinContent(i,0.10); // 10%
391
392   // PID efficiency (from PID/noPID)
393   fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
394   for(Int_t i=1;i<=20;i++) fPIDEff->SetBinContent(i,0.04); // 3%
395   fPIDEff->SetBinContent(1,1); // 100%
396   fPIDEff->SetBinContent(2,1); // 100%
397   fPIDEff->SetBinContent(3,0.05); // 5%
398   fPIDEff->SetBinContent(4,0.05); // 5%
399   fPIDEff->SetBinContent(5,0.05); // 5%
400  
401
402   // MC dN/dpt  (copied from D0 : will update later)
403   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
404   for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,(Float_t)i*0.006);
405
406   return;
407 }
408 //--------------------------------------------------------------------------
409 void AliHFSystErr::InitDstartoD0pi2010PbPb020() {
410   // 
411   // D*+->D0pi syst errors. Responsible: ??
412   //  2010 PbPb sample, 0-20 CC
413   //
414
415   AliInfo("Not yet implemented");
416   return;
417 }
418 //--------------------------------------------------------------------------
419 void AliHFSystErr::InitDstartoD0pi2010PbPb4080() {
420   // 
421   // D*+->D0pi syst errors. Responsible: ??
422   //  2010 PbPb sample, 40-80 CC
423   //
424
425   AliInfo("Not yet implemented");
426   return;
427 }
428 //--------------------------------------------------------------------------
429 Double_t AliHFSystErr::GetCutsEffErr(Double_t pt) const {
430   // 
431   // Get error
432   //
433
434   Int_t bin=fCutsEff->FindBin(pt);
435
436   return fCutsEff->GetBinContent(bin);
437 }
438 //--------------------------------------------------------------------------
439 Double_t AliHFSystErr::GetMCPtShapeErr(Double_t pt) const {
440   // 
441   // Get error
442   //
443
444   Int_t bin=fMCPtShape->FindBin(pt);
445
446   return fMCPtShape->GetBinContent(bin);
447 }
448 //--------------------------------------------------------------------------
449 Double_t AliHFSystErr::GetSeleEffErr(Double_t pt) const {
450   // 
451   // Get error
452   //
453
454   Double_t err=GetCutsEffErr(pt)*GetCutsEffErr(pt)+GetMCPtShapeErr(pt)*GetMCPtShapeErr(pt);
455
456   return TMath::Sqrt(err);
457 }
458 //--------------------------------------------------------------------------
459 Double_t AliHFSystErr::GetPIDEffErr(Double_t pt) const {
460   // 
461   // Get error
462   //
463
464   Int_t bin=fPIDEff->FindBin(pt);
465
466   return fPIDEff->GetBinContent(bin);
467 }
468 //--------------------------------------------------------------------------
469 Double_t AliHFSystErr::GetTrackingEffErr(Double_t pt) const {
470   // 
471   // Get error
472   //
473
474   Int_t bin=fTrackingEff->FindBin(pt);
475
476   return fTrackingEff->GetBinContent(bin);
477 }
478 //--------------------------------------------------------------------------
479 Double_t AliHFSystErr::GetRawYieldErr(Double_t pt) const {
480   // 
481   // Get error
482   //
483
484   Int_t bin=fRawYield->FindBin(pt);
485
486   return fRawYield->GetBinContent(bin);
487 }
488 //--------------------------------------------------------------------------
489 Double_t AliHFSystErr::GetPartAntipartErr(Double_t pt) const {
490   // 
491   // Get error
492   //
493
494   Int_t bin=fPartAntipart->FindBin(pt);
495
496   return fPartAntipart->GetBinContent(bin);
497 }
498 //--------------------------------------------------------------------------
499 Double_t AliHFSystErr::GetTotalSystErr(Double_t pt,Double_t feeddownErr) const {
500   // 
501   // Get total syst error (except norm. error)
502   //
503
504   Double_t err=0.;
505
506   if(fRawYield) err += GetRawYieldErr(pt)*GetRawYieldErr(pt);
507   if(fTrackingEff) err += GetTrackingEffErr(pt)*GetTrackingEffErr(pt);
508   if(fBR) err += GetBRErr()*GetBRErr();
509   if(fCutsEff) err += GetCutsEffErr(pt)*GetCutsEffErr(pt);
510   if(fPIDEff) err += GetPIDEffErr(pt)*GetPIDEffErr(pt);
511   if(fMCPtShape) err += GetMCPtShapeErr(pt)*GetMCPtShapeErr(pt);
512   if(fPartAntipart) err += GetPartAntipartErr(pt)*GetPartAntipartErr(pt);
513
514   err += feeddownErr*feeddownErr;
515
516   return TMath::Sqrt(err);
517 }
518 //---------------------------------------------------------------------------
519 void AliHFSystErr::DrawErrors(TGraphAsymmErrors *grErrFeeddown) const {
520   //
521   // Draw errors
522   //
523   gStyle->SetOptStat(0);
524
525   TCanvas *cSystErr = new TCanvas("cSystErr","Systematic Errors",300,80,640,500);
526   cSystErr->Range(0.20,-0.5,18.4,0.34);
527   cSystErr->SetRightMargin(0.318);
528   cSystErr->SetFillColor(0);
529
530   TH2F *hFrame = new TH2F("hFrame","Systematic errors; p_{t} [GeV/c]; Relative Error",20,0,20,100,-1,+1);
531   hFrame->SetAxisRange(2.,11.9,"X");
532   hFrame->SetAxisRange(-0.5,0.5,"Y");
533   hFrame->Draw();
534
535   TLegend *leg = new TLegend(0.69,0.44,0.98,0.86,NULL,"brNDC");
536   leg->SetTextSize(0.03601695);
537   leg->SetFillStyle(0);
538   leg->SetBorderSize(0);
539   
540   TH1F *hTotErr=new TH1F("hTotErr","",20,0,20);
541   Int_t nbins = fNorm->GetNbinsX();
542   TGraphAsymmErrors *gTotErr = new TGraphAsymmErrors(nbins);
543   for(Int_t i=1;i<=20;i++) {
544     Double_t pt = hTotErr->GetBinCenter(i);
545     Double_t ptwidth = hTotErr->GetBinWidth(i);
546
547     if(grErrFeeddown) {
548       Double_t x=0., y=0., errxl=0., errxh=0., erryl=0., erryh=0.;
549       Double_t toterryl=0., toterryh=0.;
550       for(Int_t j=0; j<grErrFeeddown->GetN(); j++) {
551         grErrFeeddown->GetPoint(j,x,y);
552         errxh = grErrFeeddown->GetErrorXhigh(j);
553         errxl = grErrFeeddown->GetErrorXlow(j);
554         if ( ( (x-errxl) <= pt) && ( (x+errxl) >= pt) ) {
555           erryh = grErrFeeddown->GetErrorYhigh(j);
556           erryl = grErrFeeddown->GetErrorYlow(j);
557         }
558       }
559       if (erryl>=1e-3) toterryl = GetTotalSystErr(pt,erryl);
560       else toterryl = GetTotalSystErr(pt);
561       if (erryh>=1e-3) toterryh = GetTotalSystErr(pt,erryh);
562       else toterryh = GetTotalSystErr(pt);
563
564       hTotErr->SetBinContent(i,toterryh);
565       gTotErr->SetPoint(i,pt,0.);
566       gTotErr->SetPointError(i,ptwidth/2.,ptwidth/2.,toterryl,toterryh); // i, exl, exh, eyl, eyh
567     }
568     else {
569       hTotErr->SetBinContent(i,GetTotalSystErr(pt));
570       gTotErr->SetPoint(i,pt,0.);
571       gTotErr->SetPointError(i,ptwidth/2.,ptwidth/2.,GetTotalSystErr(pt),GetTotalSystErr(pt)); // i, exl, exh, eyl, eyh
572     }
573
574   }
575   gTotErr->SetLineColor(kBlack);
576   gTotErr->SetFillColor(kRed);
577   gTotErr->SetFillStyle(3002);
578   gTotErr->Draw("2");
579   leg->AddEntry(gTotErr,"Total (excl. norm.)","f");
580 //   hTotErr->SetLineColor(1);
581 //   hTotErr->SetLineWidth(3);
582 //   hTotErr->Draw("same");
583 //   leg->AddEntry(hTotErr,"Total (excl. norm.)","l");
584   
585
586   fNorm->SetFillColor(1);
587   fNorm->SetFillStyle(3002);
588   //fNorm->Draw("same");
589   //TH1F *hNormRefl = ReflectHisto(fNorm);
590   //hNormRefl->Draw("same");
591   leg->AddEntry(fNorm,"Normalization (10%)","");
592
593   if(grErrFeeddown) {
594     grErrFeeddown->SetFillColor(kTeal-8);
595     grErrFeeddown->SetFillStyle(3001);
596     grErrFeeddown->Draw("2");
597     leg->AddEntry(grErrFeeddown,"Feed-down from B","f");
598   }
599   if(fTrackingEff) {
600     fTrackingEff->SetFillColor(4);
601     fTrackingEff->SetFillStyle(3006);
602     fTrackingEff->Draw("same");
603     TH1F *hTrackingEffRefl = ReflectHisto(fTrackingEff);
604     hTrackingEffRefl->Draw("same");
605     leg->AddEntry(fTrackingEff,"Tracking efficiency","f");
606   }
607   if(fBR) {
608     fBR->SetFillColor(6);
609     fBR->SetFillStyle(3005);
610     //fBR->SetFillStyle(3020);
611     fBR->Draw("same");
612     TH1F *hBRRefl = ReflectHisto(fBR);
613     hBRRefl->Draw("same");
614     leg->AddEntry(fBR,"Branching ratio","f");
615   }
616   if(fRawYield) {
617     Int_t ci;   // for color index setting
618     ci = TColor::GetColor("#00cc00");
619     fRawYield->SetLineColor(ci);
620     //    fRawYield->SetLineColor(3);
621     fRawYield->SetLineWidth(3);
622     fRawYield->Draw("same");
623     TH1F *hRawYieldRefl = ReflectHisto(fRawYield);
624     hRawYieldRefl->Draw("same");
625     leg->AddEntry(fRawYield,"Yield extraction","l");
626   }
627   if(fCutsEff) {
628     fCutsEff->SetLineColor(4);
629     fCutsEff->SetLineWidth(3);
630     fCutsEff->Draw("same");
631     TH1F *hCutsEffRefl = ReflectHisto(fCutsEff);
632     hCutsEffRefl->Draw("same");
633     leg->AddEntry(fCutsEff,"Cuts efficiency","l");
634   }
635   if(fPIDEff) {
636     fPIDEff->SetLineColor(7);
637     fPIDEff->SetLineWidth(3);
638     fPIDEff->Draw("same");
639     TH1F *hPIDEffRefl = ReflectHisto(fPIDEff);
640     hPIDEffRefl->Draw("same");
641     leg->AddEntry(fPIDEff,"PID efficiency","l");
642   }
643   if(fMCPtShape) {
644     Int_t ci = TColor::GetColor("#9933ff");
645     fMCPtShape->SetLineColor(ci);
646     //    fMCPtShape->SetLineColor(8);
647     fMCPtShape->SetLineWidth(3);
648     fMCPtShape->Draw("same");
649     TH1F *hMCPtShapeRefl = ReflectHisto(fMCPtShape);
650     hMCPtShapeRefl->Draw("same");
651     leg->AddEntry(fMCPtShape,"MC p_{t} shape","l");
652   }
653   if(fPartAntipart) {
654     Int_t ci = TColor::GetColor("#ff6600");
655     fPartAntipart->SetLineColor(ci);
656     //    fPartAntipart->SetLineColor(9);
657     fPartAntipart->SetLineWidth(3);
658     fPartAntipart->Draw("same");
659     TH1F *hPartAntipartRefl = ReflectHisto(fPartAntipart);
660     hPartAntipartRefl->Draw("same");
661     leg->AddEntry(fPartAntipart,"D = #bar{D}","l");
662   }
663
664
665   leg->Draw();
666
667   cSystErr->SaveAs("RelativeSystematics.eps");
668
669   return;
670 }
671 //-------------------------------------------------------------------------
672 TH1F* AliHFSystErr::ReflectHisto(TH1F *hin) const {
673   //
674   // Clones and reflects histogram 
675   // 
676   TH1F *hout=(TH1F*)hin->Clone("hout");
677   hout->Scale(-1.);
678
679   return hout;
680 }