Updates in PbPb cuts (Andrea Rossi)
[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 fIsLowEnergy(false)
63 {
64   //
65   // Default Constructor
66   //
67 }
68
69 //--------------------------------------------------------------------------
70 AliHFSystErr::~AliHFSystErr() {
71   //  
72   // Default Destructor
73   //
74   /*
75
76   if(fNorm)         { delete fNorm; fNorm=0; }
77   if(fRawYield)     { delete fRawYield; fRawYield=0; }
78   if(fTrackingEff)  { delete fTrackingEff; fTrackingEff=0; }
79   if(fBR)           { delete fBR; fBR=0; }
80   if(fCutsEff)      { delete fCutsEff; fCutsEff=0; }
81   if(fPIDEff)       { delete fPIDEff; fPIDEff=0; }
82   if(fMCPtShape)    { delete fMCPtShape; fMCPtShape=0; }
83   if(fPartAntipart) { delete fPartAntipart; fPartAntipart=0; }
84   */
85 }
86
87 //--------------------------------------------------------------------------
88 void AliHFSystErr::Init(Int_t decay){
89   //
90   // Variables/histos initialization
91   //
92
93   if (fRunNumber!=10 && fIsLowEnergy==false) { 
94     AliError("Only settings for 2010 and the low energy runs are implemented so far");
95   }
96   if (fCentralityClass!="020" && fCentralityClass!="4080" && fCentralityClass!="0100"){
97     AliError("Only settings for MB2010 are implemented so far");
98   }
99
100   switch(decay) {
101   case 1: // D0->Kpi
102     if (fCollisionType==0) {
103       if (fIsLowEnergy) InitD0toKpi2010ppLowEn();
104       else InitD0toKpi2010pp();
105     } else if (fCollisionType==1) {
106       if (fCentralityClass=="020") InitD0toKpi2010PbPb020();
107       else if (fCentralityClass=="4080") InitD0toKpi2010PbPb4080();
108       else AliError("Not yet implemented");
109     }
110     //    else if (fCollisionType==2) InitD0toKpi2010ppLowEn();
111     break;
112   case 2: // D+->Kpipi
113     if (fCollisionType==0) {
114       if (fIsLowEnergy) InitDplustoKpipi2010ppLowEn();
115       else InitDplustoKpipi2010pp();
116     } else if (fCollisionType==1) {
117       if (fCentralityClass=="020") InitDplustoKpipi2010PbPb020();
118       else if (fCentralityClass=="4080") InitDplustoKpipi2010PbPb4080();
119       else AliError("Not yet implemented");
120     }
121     break;
122   case 3: // D*->D0pi
123     if (fCollisionType==0) {
124       if(fIsLowEnergy)  InitDstartoD0pi2010ppLowEn();
125       else InitDstartoD0pi2010pp();
126     }else if (fCollisionType==1) {
127       if (fCentralityClass=="020")  InitDstartoD0pi2010PbPb020();
128       if (fCentralityClass=="2040") InitDstartoD0pi2010PbPb2040();
129       if (fCentralityClass=="4080") InitDstartoD0pi2010PbPb4080();
130       if (fCentralityClass!="4080" && fCentralityClass!="2040" && fCentralityClass!="020")  AliError("Not yet implemented");
131     }
132     break;
133   default:
134     printf("Invalid decay type: %d\n",decay);
135     break;
136   }
137
138 }
139   
140 //--------------------------------------------------------------------------
141 void AliHFSystErr::InitD0toKpi2010pp() {
142   // 
143   // D0->Kpi syst errors. Responsible: A. Rossi
144   //   2010 pp sample
145   //
146
147   // Normalization
148   fNorm = new TH1F("fNorm","fNorm",24,0,24);
149   for(Int_t i=1;i<=24;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
150
151   // Branching ratio 
152   fBR = new TH1F("fBR","fBR",24,0,24);
153   for(Int_t i=1;i<=24;i++) fBR->SetBinContent(i,0.012); // 1.2% PDG2010
154
155   // Tracking efficiency
156   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",24,0,24);
157   for(Int_t i=1;i<=24;i++) fTrackingEff->SetBinContent(i,0.02); // 2% (1% per track)
158
159   // Raw yield extraction
160   fRawYield = new TH1F("fRawYield","fRawYield",24,0,24);
161   fRawYield->SetBinContent(1,1);
162   fRawYield->SetBinContent(2,0.2); 
163   for(Int_t i=3;i<=16;i++) fRawYield->SetBinContent(i,0.05);
164   for(Int_t i=13;i<=16;i++) fRawYield->SetBinContent(i,0.10);
165   for(Int_t i=17;i<=24;i++) fRawYield->SetBinContent(i,1);
166
167   // Cuts efficiency (from cuts variation)
168   fCutsEff = new TH1F("fCutsEff","fCutsEff",24,0,24);
169   for(Int_t i=1;i<=24;i++) fCutsEff->SetBinContent(i,0.10); // 10%
170
171   // PID efficiency (from PID/noPID)
172   fPIDEff = new TH1F("fPIDEff","fPIDEff",24,0,24);
173   for(Int_t i=1;i<=24;i++) fPIDEff->SetBinContent(i,0.03); // 3%
174   fPIDEff->SetBinContent(2,0.10); // 10%
175
176   // MC dN/dpt
177   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",24,0,24);
178   for(Int_t i=1;i<=24;i++) fMCPtShape->SetBinContent(i,0);
179   fMCPtShape->SetBinContent(1,0.03);
180   fMCPtShape->SetBinContent(2,0.03);
181
182   // particle-antiparticle
183   //  fPartAntipart = new TH1F("fPartAntipart","fPartAntipart",24,0,24);
184   //  fPartAntipart->SetBinContent(1,1); 
185   //  for(Int_t i=2;i<=24;i++) fPartAntipart->SetBinContent(i,0.05);
186   
187   return;
188 }
189 //--------------------------------------------------------------------------
190 void AliHFSystErr::InitD0toKpi2010PbPb020() {
191   // 
192   // D0->Kpi syst errors. Responsible: A. Rossi
193   //   2010 PbPb sample, 0-20 CC
194   //
195
196   // Normalization
197   fNorm = new TH1F("fNorm","fNorm",20,0,20);
198   for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
199
200   // Branching ratio 
201   fBR = new TH1F("fBR","fBR",20,0,20);
202   for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.012); // 1.2% PDG2010
203
204   // Tracking efficiency
205   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
206   for(Int_t i=2;i<=12;i++) fTrackingEff->SetBinContent(i,0.05+0.005*(Float_t)i);
207  
208
209   // Raw yield extraction
210   fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
211   fRawYield->SetBinContent(1,0);
212   fRawYield->SetBinContent(2,0);
213   fRawYield->SetBinContent(3,0.05);
214   fRawYield->SetBinContent(4,0.05);
215   fRawYield->SetBinContent(5,0.10);
216   fRawYield->SetBinContent(6,0.10);
217   for(Int_t i=7;i<=12;i++) fRawYield->SetBinContent(i,0.05);
218  // for(Int_t i=13;i<=20;i++) fRawYield->SetBinContent(i,0);
219
220   // Cuts efficiency (from cuts variation)
221   fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
222   fCutsEff->SetBinContent(1,0.);
223   fCutsEff->SetBinContent(2,0.);
224   fCutsEff->SetBinContent(3,0.18);
225   fCutsEff->SetBinContent(4,0.18);
226   fCutsEff->SetBinContent(5,0.18);
227   fCutsEff->SetBinContent(6,0.15);
228   fCutsEff->SetBinContent(7,0.15);
229   fCutsEff->SetBinContent(8,0.15);
230 fCutsEff->SetBinContent(9,0.15);
231  fCutsEff->SetBinContent(10,0.15);
232   fCutsEff->SetBinContent(11,0.15);
233   fCutsEff->SetBinContent(12,0.15);
234  // for(Int_t i=13;i<=20;i++) fCutsEff->SetBinContent(i,0.);
235
236   // PID efficiency (from PID/noPID)
237   fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
238   for(Int_t i=3;i<=12;i++) fPIDEff->SetBinContent(i,0.05);
239   fPIDEff->SetBinContent(5,0.10);
240   fPIDEff->SetBinContent(6,0.10);
241
242   // MC dN/dpt
243   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
244   for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,0);
245   fMCPtShape->SetBinContent(1,0.03);
246   fMCPtShape->SetBinContent(2,0.03);
247
248
249   // particle-antiparticle
250   fPartAntipart = new TH1F("fPartAntipart","fPartAntipart",20,0,20);
251   for(Int_t i=3;i<=12;i++) fPartAntipart->SetBinContent(i,0.05);
252   fPartAntipart->SetBinContent(3,0.10);
253   fPartAntipart->SetBinContent(4,0.10);
254   fPartAntipart->SetBinContent(7,0.10);
255   fPartAntipart->SetBinContent(8,0.10);
256   
257   return;
258 }
259 //--------------------------------------------------------------------------
260 void AliHFSystErr::InitD0toKpi2010PbPb4080() {
261   // 
262   // D0->Kpi syst errors. Responsible: A. Rossi
263   //   2010 PbPb sample, 40-80 CC
264   //
265
266   // Normalization
267   fNorm = new TH1F("fNorm","fNorm",20,0,20);
268   for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
269
270   // Branching ratio 
271   fBR = new TH1F("fBR","fBR",20,0,20);
272   for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.012); // 1.2% PDG2010
273
274   // Tracking efficiency
275   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
276   for(Int_t i=2;i<=12;i++) fTrackingEff->SetBinContent(i,0.5*(0.05+0.005*(Float_t)i));
277
278
279   // Raw yield extraction
280   fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
281   fRawYield->SetBinContent(1,0);
282   fRawYield->SetBinContent(2,0);
283   for(Int_t i=3;i<=12;i++) fRawYield->SetBinContent(i,0.05);
284   //for(Int_t i=13;i<=20;i++) fRawYield->SetBinContent(i,0);
285
286   // Cuts efficiency (from cuts variation)
287   fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
288   fCutsEff->SetBinContent(1,0.);
289   fCutsEff->SetBinContent(2,0.);
290   fCutsEff->SetBinContent(3,0.15);
291   fCutsEff->SetBinContent(4,0.15);
292   fCutsEff->SetBinContent(5,0.15);
293   fCutsEff->SetBinContent(6,0.15);
294   fCutsEff->SetBinContent(7,0.15);
295   fCutsEff->SetBinContent(8,0.15);
296   fCutsEff->SetBinContent(9,0.15);
297   fCutsEff->SetBinContent(10,0.15);
298   fCutsEff->SetBinContent(11,0.15);
299   fCutsEff->SetBinContent(12,0.15);
300  // for(Int_t i=13;i<=20;i++) fCutsEff->SetBinContent(i,0.);
301
302   // PID efficiency (from PID/noPID)
303   fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
304   for(Int_t i=3;i<=12;i++) fPIDEff->SetBinContent(i,0.05);
305
306   // MC dN/dpt
307   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
308   for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,0);
309   fMCPtShape->SetBinContent(1,0.03);
310   fMCPtShape->SetBinContent(2,0.03);
311
312   // particle-antiparticle
313   fPartAntipart = new TH1F("fPartAntipart","fPartAntipart",20,0,20);
314   for(Int_t i=3;i<=12;i++) fPartAntipart->SetBinContent(i,0.05);
315   
316   return;
317 }
318
319 //--------------------------------------------------------------------------
320 void AliHFSystErr::InitD0toKpi2010ppLowEn() {
321   // 
322   // D0->Kpi syst errors. Low energy run
323   //   2011 2.76 TeV pp sample
324   //
325   AliInfo(" Settings for D0 --> K pi, p-p collisions at 2.76 TeV"); 
326
327   // Normalization
328   fNorm = new TH1F("fNorm","fNorm",20,0,20);
329   for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
330
331   // Branching ratio 
332   fBR = new TH1F("fBR","fBR",20,0,20);
333   for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.012); // 1.2% PDG2010
334
335   // Tracking efficiency
336   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
337   for(Int_t i=1;i<=20;i++) fTrackingEff->SetBinContent(i,0.10); //10% (5% per track)
338
339   // Raw yield extraction
340   fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
341   fRawYield->SetBinContent(1,1);
342   fRawYield->SetBinContent(2,1);
343   fRawYield->SetBinContent(3,0.2);
344   fRawYield->SetBinContent(4,0.2);
345   fRawYield->SetBinContent(5,0.1);
346   fRawYield->SetBinContent(6,0.1);
347   fRawYield->SetBinContent(7,0.2);
348   fRawYield->SetBinContent(8,0.2);
349   for(Int_t i=9;i<=20;i++) fRawYield->SetBinContent(i,0.065);
350
351   // Cuts efficiency (from cuts variation)
352   fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
353   for(Int_t i=1;i<=20;i++) fCutsEff->SetBinContent(i,0.18); // 10% 
354
355   // PID efficiency (from PID/noPID)
356   fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
357   for(Int_t i=1;i<=20;i++) fPIDEff->SetBinContent(i,0.20); // 10%
358   fPIDEff->SetBinContent(3,0.25); // 10%
359   fPIDEff->SetBinContent(4,0.25); // 10%
360   fPIDEff->SetBinContent(7,0.25); // 10%
361   fPIDEff->SetBinContent(8,0.25); // 10%
362
363   // MC dN/dpt
364   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
365   for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,0);
366   fMCPtShape->SetBinContent(1,0.03);
367   fMCPtShape->SetBinContent(2,0.03);
368
369   // particle-antiparticle
370   fPartAntipart = new TH1F("fPartAntipart","fPartAntipart",20,0,20);
371   fPartAntipart->SetBinContent(1,1);
372   fPartAntipart->SetBinContent(2,1);
373   for(Int_t i=3;i<=6;i++) fPartAntipart->SetBinContent(i,0.08);
374   
375   return;
376 }
377
378 //--------------------------------------------------------------------------
379 void AliHFSystErr::InitDplustoKpipi2010pp() {
380   // 
381   // D+->Kpipi syst errors. Responsible: R. Bala
382   //  2010 pp sample
383   //
384
385  // Normalization
386   fNorm = new TH1F("fNorm","fNorm",20,0,20);
387   for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
388
389   // Branching ratio 
390   fBR = new TH1F("fBR","fBR",20,0,20);
391   for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.04); // 4% PDG2010
392
393   // Tracking efficiency
394   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
395   for(Int_t i=1;i<=20;i++) fTrackingEff->SetBinContent(i,0.03); // 3% (1% per track)
396
397
398   // Raw yield extraction
399   fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
400   fRawYield->SetBinContent(1,1);
401   fRawYield->SetBinContent(2,1);
402   fRawYield->SetBinContent(3,0.20);
403   for(Int_t i=4;i<=20;i++) fRawYield->SetBinContent(i,0.055);  //5 to 10%
404
405   // Cuts efficiency (from cuts variation)
406   fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
407   for(Int_t i=1;i<=20;i++) fCutsEff->SetBinContent(i,0.10); // 10%
408
409   // PID efficiency (from PID/noPID)
410   fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
411   for(Int_t i=1;i<=20;i++) fPIDEff->SetBinContent(i,0.05); // 5%
412   fPIDEff->SetBinContent(3,0.13); // 13%
413  
414
415   // MC dN/dpt  (copied from D0 : will update later)
416   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
417   for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,0);
418   fMCPtShape->SetBinContent(1,0.03);
419   fMCPtShape->SetBinContent(2,0.03);
420
421
422   // particle-antiparticle
423   /*
424   fPartAntipart = new TH1F("fPartAntipart","fPartAntipart",20,0,20);
425   fPartAntipart->SetBinContent(1,1);
426   fPartAntipart->SetBinContent(2,1);
427   fPartAntipart->SetBinContent(3,0.12);
428   for(Int_t i=4;i<=20;i++) fPartAntipart->SetBinContent(i,0.05);   //5 to 12%
429   */
430   return;
431 }
432 //--------------------------------------------------------------------------
433 void AliHFSystErr::InitDplustoKpipi2010PbPb020() {
434   // 
435   // D+->Kpipi syst errors. Responsible: ??
436   //  2010 PbPb sample, 0-20 CC
437   //
438
439  // Normalization
440   fNorm = new TH1F("fNorm","fNorm",20,0,20);
441   for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
442
443   // Branching ratio 
444   fBR = new TH1F("fBR","fBR",20,0,20);
445   for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.04); // 4% PDG2010
446
447   // Tracking efficiency
448   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
449   for(Int_t i=1;i<=20;i++) fTrackingEff->SetBinContent(i,i,0.075+0.005*(Float_t)i); 
450
451
452   // Raw yield extraction
453   fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
454   for(Int_t i=1;i<=20;i++) fRawYield->SetBinContent(i,1);  //5 to 10%
455 //   for(Int_t i=6;i<=7;i++) fRawYield->SetBinContent(i,0.32);
456 //   for(Int_t i=7;i<=9;i++) fRawYield->SetBinContent(i,0.32);
457 //   for(Int_t i=9;i<=12;i++) fRawYield->SetBinContent(i,0.05);//31./156.);
458   fRawYield->SetBinContent(6,0.12);
459   fRawYield->SetBinContent(7,0.17);
460   fRawYield->SetBinContent(8,0.17);
461   fRawYield->SetBinContent(9,0.08);
462   fRawYield->SetBinContent(10,0.08);
463   fRawYield->SetBinContent(11,0.08);
464   fRawYield->SetBinContent(12,0.08);
465
466   // Cuts efficiency (from cuts variation)
467   fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
468   for(Int_t i=1;i<=20;i++) fCutsEff->SetBinContent(i,0.15); // 10%
469
470   // PID efficiency (from PID/noPID)
471   fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
472   for(Int_t i=1;i<=20;i++) fPIDEff->SetBinContent(i,0.05); // 5%
473   fPIDEff->SetBinContent(3,0.13); // 13%
474  
475
476   // MC dN/dpt  (copied from D0 : will update later)
477   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
478   for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,0.);
479   fMCPtShape->SetBinContent(1,0.03);
480   fMCPtShape->SetBinContent(2,0.03);
481
482
483   // particle-antiparticle
484   /*
485   fPartAntipart = new TH1F("fPartAntipart","fPartAntipart",20,0,20);
486   fPartAntipart->SetBinContent(1,1);
487   fPartAntipart->SetBinContent(2,1);
488   fPartAntipart->SetBinContent(3,0.12);
489   for(Int_t i=4;i<=20;i++) fPartAntipart->SetBinContent(i,0.05);   //5 to 12%
490   */
491   return;
492 }
493
494 //--------------------------------------------------------------------------
495 void AliHFSystErr::InitDplustoKpipi2010PbPb4080() {
496   // 
497   // D+->Kpipi syst errors. Responsible: ??
498   //  2010 PbPb sample, 40-80 CC
499   //
500   
501
502  // Normalization
503   fNorm = new TH1F("fNorm","fNorm",20,0,20);
504   for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
505
506   // Branching ratio 
507   fBR = new TH1F("fBR","fBR",20,0,20);
508   for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.04); // 4% PDG2010
509
510   // Tracking efficiency
511   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
512   for(Int_t i=1;i<=20;i++) fTrackingEff->SetBinContent(i,0.03); // 3% (1% per track)
513
514
515   // Raw yield extraction
516   fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
517   fRawYield->SetBinContent(1,1);
518   fRawYield->SetBinContent(2,1);
519   fRawYield->SetBinContent(3,1);
520   fRawYield->SetBinContent(4,1.);
521   fRawYield->SetBinContent(5,4./23.);
522   fRawYield->SetBinContent(6,4./23.);
523   fRawYield->SetBinContent(7,5./20.);
524   fRawYield->SetBinContent(8,5./20.);
525   for(Int_t i=9;i<=12;i++) fRawYield->SetBinContent(i,1./16.);
526   for(Int_t i=13;i<=20;i++) fRawYield->SetBinContent(i,1);  //5 to 10%
527
528   // Cuts efficiency (from cuts variation)
529   fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
530   for(Int_t i=1;i<=20;i++) fCutsEff->SetBinContent(i,0.10); // 10%
531
532   // PID efficiency (from PID/noPID)
533   fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
534   for(Int_t i=1;i<=20;i++) fPIDEff->SetBinContent(i,0.05); // 5%
535   fPIDEff->SetBinContent(3,0.13); // 13%
536  
537
538   // MC dN/dpt  (copied from D0 : will update later)
539   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
540   for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,0);
541   fMCPtShape->SetBinContent(1,0.03);
542   fMCPtShape->SetBinContent(2,0.03);
543
544
545   // particle-antiparticle
546   /*
547   fPartAntipart = new TH1F("fPartAntipart","fPartAntipart",20,0,20);
548   fPartAntipart->SetBinContent(1,1);
549   fPartAntipart->SetBinContent(2,1);
550   fPartAntipart->SetBinContent(3,0.12);
551   for(Int_t i=4;i<=20;i++) fPartAntipart->SetBinContent(i,0.05);   //5 to 12%
552   */
553   return;
554 }
555
556 //--------------------------------------------------------------------------
557 void AliHFSystErr::InitDplustoKpipi2010ppLowEn() {
558
559   // 
560   // D+->Kpipi syst errors. Responsible: R. Bala
561   //  2011 2.76 TeV pp sample
562   //
563   AliInfo(" Settings for D+ --> K pi pi p-p collisions at 2.76 TeV"); 
564
565   // Normalization
566   fNorm = new TH1F("fNorm","fNorm",20,0,20);
567   for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
568
569   // Branching ratio 
570   fBR = new TH1F("fBR","fBR",20,0,20);
571   for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.04); // 4% PDG2010
572
573   // Tracking efficiency
574   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
575   for(Int_t i=1;i<=20;i++) fTrackingEff->SetBinContent(i,0.15); // 3% (1% per track)
576
577   // Raw yield extraction
578   fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
579   fRawYield->SetBinContent(1,1);
580   fRawYield->SetBinContent(2,1);
581   fRawYield->SetBinContent(3,0.20);
582   fRawYield->SetBinContent(4,0.20);
583   for(Int_t i=5;i<=20;i++) fRawYield->SetBinContent(i,0.055);  //5 to 10%
584
585   // Cuts efficiency (from cuts variation)
586   fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
587   for(Int_t i=1;i<=20;i++) fCutsEff->SetBinContent(i,0.10); // 10%
588
589   // PID efficiency (from PID/noPID)
590   fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
591   for(Int_t i=1;i<=20;i++) fPIDEff->SetBinContent(i,0.05); // 5%
592   fPIDEff->SetBinContent(3,0.15); // 13%
593   fPIDEff->SetBinContent(4,0.15); // 13%
594  
595   // MC dN/dpt  (copied from D0 : will update later)
596   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
597   for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,0.);
598   fMCPtShape->SetBinContent(1,0.03);
599   fMCPtShape->SetBinContent(2,0.03);
600
601   return;
602 }
603
604 //--------------------------------------------------------------------------
605 void AliHFSystErr::InitDstartoD0pi2010pp() {
606   // 
607   // D*+->D0pi syst errors. Responsible: A. Grelli, Y. Wang
608   //  2010 pp sample
609   //
610
611  // Normalization
612   fNorm = new TH1F("fNorm","fNorm",24,0,24);
613   for(Int_t i=1;i<=24;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
614
615   // Branching ratio 
616   fBR = new TH1F("fBR","fBR",24,0,24);
617   for(Int_t i=1;i<=24;i++) fBR->SetBinContent(i,0.015); // 1.5% PDG2010
618
619   // Tracking efficiency
620   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",24,0,24);
621   fTrackingEff->SetBinContent(1,1.0);
622   fTrackingEff->SetBinContent(2,0.12);
623   fTrackingEff->SetBinContent(3,0.08);
624   fTrackingEff->SetBinContent(3,0.05);
625   for(Int_t i=4;i<=24;i++) fTrackingEff->SetBinContent(i,0.03); // 3% (1% per track)
626
627
628   // Raw yield extraction
629   fRawYield = new TH1F("fRawYield","fRawYield",24,0,24);
630   fRawYield->SetBinContent(1,1.0);
631   fRawYield->SetBinContent(2,0.09);
632   fRawYield->SetBinContent(3,0.04);
633   fRawYield->SetBinContent(4,0.02);
634   fRawYield->SetBinContent(5,0.03);
635   fRawYield->SetBinContent(6,0.06);
636   fRawYield->SetBinContent(7,0.04);
637   fRawYield->SetBinContent(8,0.11);
638   for(Int_t i=5;i<=24;i++) fRawYield->SetBinContent(i,0.03);  //4%
639   fRawYield->SetBinContent(13,0.09);
640   fRawYield->SetBinContent(14,0.09);
641   fRawYield->SetBinContent(15,0.09);
642   fRawYield->SetBinContent(16,0.09);
643   fRawYield->SetBinContent(17,0.24);
644   fRawYield->SetBinContent(18,0.24);
645   fRawYield->SetBinContent(19,0.24);
646   fRawYield->SetBinContent(20,0.24);
647   fRawYield->SetBinContent(21,0.24);
648   fRawYield->SetBinContent(22,0.24);
649   fRawYield->SetBinContent(23,0.24);
650   fRawYield->SetBinContent(24,0.24);
651
652   // Cuts efficiency (from cuts variation)
653   fCutsEff = new TH1F("fCutsEff","fCutsEff",24,0,24);
654   fCutsEff->SetBinContent(2,0.22);
655   for(Int_t i=3;i<=24;i++) fCutsEff->SetBinContent(i,0.10); // 10%
656
657   // PID efficiency (from PID/noPID)
658   fPIDEff = new TH1F("fPIDEff","fPIDEff",24,0,24);
659   for(Int_t i=1;i<=24;i++) fPIDEff->SetBinContent(i,0.04); // 3%
660  
661
662   // MC dN/dpt  (copied from D0 : will update later)
663   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",24,0,24);
664   for(Int_t i=1;i<=24;i++) fMCPtShape->SetBinContent(i,0);
665   fMCPtShape->SetBinContent(1,0.03);
666   fMCPtShape->SetBinContent(2,0.03);
667
668   return;
669
670
671 }
672 //--------------------------------------------------------------------------
673 void AliHFSystErr::InitDstartoD0pi2010ppLowEn() {
674
675   // 
676   // D+->Kpipi syst errors. Responsible: A. Grelli
677   //  2011 2.76 TeV pp sample
678   //
679   AliInfo(" Settings for D*+ --> D0 pi p-p collisions at 2.76 TeV"); 
680
681 // Normalization
682   fNorm = new TH1F("fNorm","fNorm",20,0,20);
683   for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
684
685   // Branching ratio 
686   fBR = new TH1F("fBR","fBR",20,0,20);
687   for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.015); // 1.5% PDG2010
688
689   // Tracking efficiency
690   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
691   for(Int_t i=1;i<=20;i++) fTrackingEff->SetBinContent(i,0.15); //10% (to be checked!!)
692
693   // Raw yield extraction
694   fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
695   fRawYield->SetBinContent(1,1);
696   fRawYield->SetBinContent(2,1);
697   fRawYield->SetBinContent(3,0.2);
698   fRawYield->SetBinContent(4,0.2);
699   fRawYield->SetBinContent(5,0.08);
700   fRawYield->SetBinContent(6,0.08);
701   fRawYield->SetBinContent(7,0.1);
702   fRawYield->SetBinContent(8,0.1);
703   fRawYield->SetBinContent(9,0.2);
704   fRawYield->SetBinContent(10,0.2);
705   for(Int_t i=9;i<=20;i++) fRawYield->SetBinContent(i,0.065);
706
707   // Cuts efficiency (from cuts variation)
708   fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
709   for(Int_t i=1;i<=20;i++) fCutsEff->SetBinContent(i,0.18);  
710   fCutsEff->SetBinContent(3,0.2);
711   fCutsEff->SetBinContent(4,0.2);
712   fCutsEff->SetBinContent(5,0.1);
713   fCutsEff->SetBinContent(6,0.1);
714   fCutsEff->SetBinContent(7,0.12);
715   fCutsEff->SetBinContent(8,0.12);
716   fCutsEff->SetBinContent(9,0.2);
717   fCutsEff->SetBinContent(10,0.2);
718
719   // PID efficiency (from PID/noPID)
720   fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
721   for(Int_t i=1;i<=20;i++) fPIDEff->SetBinContent(i,0.05); // 10%
722
723   // MC dN/dpt
724   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
725   for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,0.);
726   fMCPtShape->SetBinContent(1,0.03);
727   fMCPtShape->SetBinContent(2,0.03);
728
729
730   return;
731 }
732 //------------------------------------------------------------------------
733 void AliHFSystErr::InitDstartoD0pi2010PbPb020() {
734   // 
735   // D*+->D0pi syst errors. Responsible: A. Grelli
736   //  2010 PbPb sample, 0-20 CC
737   //
738
739   AliInfo(" Settings for D*+ --> D0pi Pb-Pb collisions at 2.76 TeV - 0-20 centrality - DUMMY"); 
740
741  // Normalization
742   fNorm = new TH1F("fNorm","fNorm",24,0,24);
743   for(Int_t i=1;i<=24;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
744
745   // Branching ratio 
746   fBR = new TH1F("fBR","fBR",24,0,24);
747   for(Int_t i=1;i<=24;i++) fBR->SetBinContent(i,0.015); // 1.5% PDG2010
748
749   // Tracking efficiency
750   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",24,0,24);;
751   for(Int_t i=1;i<=24;i++) fTrackingEff->SetBinContent(i,0.12); // 3% (1% per track)
752
753
754   // Raw yield extraction
755   fRawYield = new TH1F("fRawYield","fRawYield",24,0,24);
756   for(Int_t i=1;i<=24;i++) fRawYield->SetBinContent(i,0.15);  //4%
757  
758   // Cuts efficiency (from cuts variation)
759   fCutsEff = new TH1F("fCutsEff","fCutsEff",24,0,24);
760   for(Int_t i=1;i<=24;i++) fCutsEff->SetBinContent(i,0.10); // 10%
761
762   // PID efficiency (from PID/noPID)
763   fPIDEff = new TH1F("fPIDEff","fPIDEff",24,0,24);
764   for(Int_t i=1;i<=24;i++) fPIDEff->SetBinContent(i,0.04); // 3%
765  
766
767   // MC dN/dpt  (copied from D0 : will update later)
768   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",24,0,24);
769   for(Int_t i=1;i<=24;i++) fMCPtShape->SetBinContent(i,0.);
770   fMCPtShape->SetBinContent(1,0.03);
771   fMCPtShape->SetBinContent(2,0.03);
772
773   return;
774
775 }
776 //-------------------------------------------------------------------------
777 void AliHFSystErr::InitDstartoD0pi2010PbPb2040() {
778   // 
779   // D*+->D0pi syst errors. Responsible: A. Grelli
780   //  2010 PbPb sample, 20-40 CC
781   //
782
783   AliInfo(" Settings for D*+ --> D0pi Pb-Pb collisions at 2.76 TeV - 20-40 centrality - DUMMY"); 
784
785  // Normalization
786   fNorm = new TH1F("fNorm","fNorm",24,0,24);
787   for(Int_t i=1;i<=24;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
788
789   // Branching ratio 
790   fBR = new TH1F("fBR","fBR",24,0,24);
791   for(Int_t i=1;i<=24;i++) fBR->SetBinContent(i,0.015); // 1.5% PDG2010
792
793   // Tracking efficiency
794   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",24,0,24);;
795   for(Int_t i=1;i<=24;i++) fTrackingEff->SetBinContent(i,0.12); // 3% (1% per track)
796
797
798   // Raw yield extraction
799   fRawYield = new TH1F("fRawYield","fRawYield",24,0,24);
800   for(Int_t i=1;i<=24;i++) fRawYield->SetBinContent(i,0.15);  //4%
801  
802   // Cuts efficiency (from cuts variation)
803   fCutsEff = new TH1F("fCutsEff","fCutsEff",24,0,24);
804   for(Int_t i=1;i<=24;i++) fCutsEff->SetBinContent(i,0.10); // 10%
805
806   // PID efficiency (from PID/noPID)
807   fPIDEff = new TH1F("fPIDEff","fPIDEff",24,0,24);
808   for(Int_t i=1;i<=24;i++) fPIDEff->SetBinContent(i,0.04); // 3%
809  
810
811   // MC dN/dpt  (copied from D0 : will update later)
812   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",24,0,24);
813   for(Int_t i=1;i<=24;i++) fMCPtShape->SetBinContent(i,0.);
814   fMCPtShape->SetBinContent(1,0.03);
815   fMCPtShape->SetBinContent(2,0.03);
816
817   return;
818
819 }
820 //--------------------------------------------------------------------------
821 void AliHFSystErr::InitDstartoD0pi2010PbPb4080() {
822   // 
823   // D*+->D0pi syst errors. Responsible: A. Grelli
824   //  2010 PbPb sample, 40-80 CC
825   //
826
827   AliInfo(" Settings for D*+ --> D0pi Pb-Pb collisions at 2.76 TeV - 40-80 centrality - DUMMY"); 
828
829  // Normalization
830   fNorm = new TH1F("fNorm","fNorm",24,0,24);
831   for(Int_t i=1;i<=24;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
832
833   // Branching ratio 
834   fBR = new TH1F("fBR","fBR",24,0,24);
835   for(Int_t i=1;i<=24;i++) fBR->SetBinContent(i,0.015); // 1.5% PDG2010
836
837   // Tracking efficiency
838   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",24,0,24);;
839   for(Int_t i=1;i<=24;i++) fTrackingEff->SetBinContent(i,0.12); // 3% (1% per track)
840
841
842   // Raw yield extraction
843   fRawYield = new TH1F("fRawYield","fRawYield",24,0,24);
844   for(Int_t i=1;i<=24;i++) fRawYield->SetBinContent(i,0.15);  //4%
845  
846   // Cuts efficiency (from cuts variation)
847   fCutsEff = new TH1F("fCutsEff","fCutsEff",24,0,24);
848   for(Int_t i=1;i<=24;i++) fCutsEff->SetBinContent(i,0.10); // 10%
849
850   // PID efficiency (from PID/noPID)
851   fPIDEff = new TH1F("fPIDEff","fPIDEff",24,0,24);
852   for(Int_t i=1;i<=24;i++) fPIDEff->SetBinContent(i,0.04); // 3%
853  
854
855   // MC dN/dpt  (copied from D0 : will update later)
856   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",24,0,24);
857   for(Int_t i=1;i<=24;i++) fMCPtShape->SetBinContent(i,0.);
858   fMCPtShape->SetBinContent(1,0.03);
859   fMCPtShape->SetBinContent(2,0.03);
860
861   return;
862
863 }
864 //--------------------------------------------------------------------------
865 Double_t AliHFSystErr::GetCutsEffErr(Double_t pt) const {
866   // 
867   // Get error
868   //
869
870   Int_t bin=fCutsEff->FindBin(pt);
871
872   return fCutsEff->GetBinContent(bin);
873 }
874 //--------------------------------------------------------------------------
875 Double_t AliHFSystErr::GetMCPtShapeErr(Double_t pt) const {
876   // 
877   // Get error
878   //
879
880   Int_t bin=fMCPtShape->FindBin(pt);
881
882   return fMCPtShape->GetBinContent(bin);
883 }
884 //--------------------------------------------------------------------------
885 Double_t AliHFSystErr::GetSeleEffErr(Double_t pt) const {
886   // 
887   // Get error
888   //
889
890   Double_t err=GetCutsEffErr(pt)*GetCutsEffErr(pt)+GetMCPtShapeErr(pt)*GetMCPtShapeErr(pt);
891
892   return TMath::Sqrt(err);
893 }
894 //--------------------------------------------------------------------------
895 Double_t AliHFSystErr::GetPIDEffErr(Double_t pt) const {
896   // 
897   // Get error
898   //
899
900   Int_t bin=fPIDEff->FindBin(pt);
901
902   return fPIDEff->GetBinContent(bin);
903 }
904 //--------------------------------------------------------------------------
905 Double_t AliHFSystErr::GetTrackingEffErr(Double_t pt) const {
906   // 
907   // Get error
908   //
909
910   Int_t bin=fTrackingEff->FindBin(pt);
911
912   return fTrackingEff->GetBinContent(bin);
913 }
914 //--------------------------------------------------------------------------
915 Double_t AliHFSystErr::GetRawYieldErr(Double_t pt) const {
916   // 
917   // Get error
918   //
919
920   Int_t bin=fRawYield->FindBin(pt);
921
922   return fRawYield->GetBinContent(bin);
923 }
924 //--------------------------------------------------------------------------
925 Double_t AliHFSystErr::GetPartAntipartErr(Double_t pt) const {
926   // 
927   // Get error
928   //
929
930   Int_t bin=fPartAntipart->FindBin(pt);
931
932   return fPartAntipart->GetBinContent(bin);
933 }
934 //--------------------------------------------------------------------------
935 Double_t AliHFSystErr::GetTotalSystErr(Double_t pt,Double_t feeddownErr) const {
936   // 
937   // Get total syst error (except norm. error)
938   //
939
940   Double_t err=0.;
941
942   if(fRawYield) err += GetRawYieldErr(pt)*GetRawYieldErr(pt);
943   if(fTrackingEff) err += GetTrackingEffErr(pt)*GetTrackingEffErr(pt);
944   if(fBR) err += GetBRErr()*GetBRErr();
945   if(fCutsEff) err += GetCutsEffErr(pt)*GetCutsEffErr(pt);
946   if(fPIDEff) err += GetPIDEffErr(pt)*GetPIDEffErr(pt);
947   if(fMCPtShape) err += GetMCPtShapeErr(pt)*GetMCPtShapeErr(pt);
948   if(fPartAntipart) err += GetPartAntipartErr(pt)*GetPartAntipartErr(pt);
949
950   err += feeddownErr*feeddownErr;
951
952   return TMath::Sqrt(err);
953 }
954 //---------------------------------------------------------------------------
955 void AliHFSystErr::DrawErrors(TGraphAsymmErrors *grErrFeeddown) const {
956   //
957   // Draw errors
958   //
959   gStyle->SetOptStat(0);
960
961   TCanvas *cSystErr = new TCanvas("cSystErr","Systematic Errors",300,80,640,500);
962   cSystErr->Range(0.20,-0.5,18.4,0.34);
963   cSystErr->SetRightMargin(0.318);
964   cSystErr->SetFillColor(0);
965
966   TH2F *hFrame = new TH2F("hFrame","Systematic errors; p_{t} [GeV/c]; Relative Error",20,0,20,100,-1,+1);
967   hFrame->SetAxisRange(1.,15.9,"X");
968   hFrame->SetAxisRange(-0.5,0.5,"Y");
969   hFrame->Draw();
970
971   TLegend *leg = new TLegend(0.69,0.44,0.98,0.86,NULL,"brNDC");
972   leg->SetTextSize(0.03601695);
973   leg->SetFillStyle(0);
974   leg->SetBorderSize(0);
975   
976   TH1F *hTotErr=new TH1F("hTotErr","",24,0,24);
977   Int_t nbins = fNorm->GetNbinsX();
978   TGraphAsymmErrors *gTotErr = new TGraphAsymmErrors(nbins);
979   for(Int_t i=1;i<=24;i++) {
980     Double_t pt = hTotErr->GetBinCenter(i);
981     Double_t ptwidth = hTotErr->GetBinWidth(i);
982
983     if(grErrFeeddown) {
984       Double_t x=0., y=0., errxl=0., errxh=0., erryl=0., erryh=0.;
985       Double_t toterryl=0., toterryh=0.;
986       for(Int_t j=0; j<grErrFeeddown->GetN(); j++) {
987         grErrFeeddown->GetPoint(j,x,y);
988         errxh = grErrFeeddown->GetErrorXhigh(j);
989         errxl = grErrFeeddown->GetErrorXlow(j);
990         if ( ( (x-errxl) <= pt) && ( (x+errxl) >= pt) ) {
991           erryh = grErrFeeddown->GetErrorYhigh(j);
992           erryl = grErrFeeddown->GetErrorYlow(j);
993         }
994       }
995       if (erryl>=1e-3) toterryl = GetTotalSystErr(pt,erryl);
996       else toterryl = GetTotalSystErr(pt);
997       if (erryh>=1e-3) toterryh = GetTotalSystErr(pt,erryh);
998       else toterryh = GetTotalSystErr(pt);
999
1000       hTotErr->SetBinContent(i,toterryh);
1001       gTotErr->SetPoint(i,pt,0.);
1002       gTotErr->SetPointError(i,ptwidth/2.,ptwidth/2.,toterryl,toterryh); // i, exl, exh, eyl, eyh
1003     }
1004     else {
1005       hTotErr->SetBinContent(i,GetTotalSystErr(pt));
1006       gTotErr->SetPoint(i,pt,0.);
1007       gTotErr->SetPointError(i,ptwidth/2.,ptwidth/2.,GetTotalSystErr(pt),GetTotalSystErr(pt)); // i, exl, exh, eyl, eyh
1008     }
1009
1010   }
1011   gTotErr->SetLineColor(kBlack);
1012   gTotErr->SetFillColor(kRed);
1013   gTotErr->SetFillStyle(3002);
1014   gTotErr->Draw("2");
1015   leg->AddEntry(gTotErr,"Total (excl. norm.)","f");
1016 //   hTotErr->SetLineColor(1);
1017 //   hTotErr->SetLineWidth(3);
1018 //   hTotErr->Draw("same");
1019 //   leg->AddEntry(hTotErr,"Total (excl. norm.)","l");
1020   
1021
1022   fNorm->SetFillColor(1);
1023   fNorm->SetFillStyle(3002);
1024   //fNorm->Draw("same");
1025   //TH1F *hNormRefl = ReflectHisto(fNorm);
1026   //hNormRefl->Draw("same");
1027   leg->AddEntry(fNorm,"Normalization (10%)","");
1028
1029   if(grErrFeeddown) {
1030     grErrFeeddown->SetFillColor(kTeal-8);
1031     grErrFeeddown->SetFillStyle(3001);
1032     grErrFeeddown->Draw("2");
1033     leg->AddEntry(grErrFeeddown,"Feed-down from B","f");
1034   }
1035   if(fTrackingEff) {
1036     fTrackingEff->SetFillColor(4);
1037     fTrackingEff->SetFillStyle(3006);
1038     fTrackingEff->Draw("same");
1039     TH1F *hTrackingEffRefl = ReflectHisto(fTrackingEff);
1040     hTrackingEffRefl->Draw("same");
1041     leg->AddEntry(fTrackingEff,"Tracking efficiency","f");
1042   }
1043   if(fBR) {
1044     fBR->SetFillColor(6);
1045     fBR->SetFillStyle(3005);
1046     //fBR->SetFillStyle(3020);
1047     fBR->Draw("same");
1048     TH1F *hBRRefl = ReflectHisto(fBR);
1049     hBRRefl->Draw("same");
1050     leg->AddEntry(fBR,"Branching ratio","f");
1051   }
1052   if(fRawYield) {
1053     Int_t ci;   // for color index setting
1054     ci = TColor::GetColor("#00cc00");
1055     fRawYield->SetLineColor(ci);
1056     //    fRawYield->SetLineColor(3);
1057     fRawYield->SetLineWidth(3);
1058     fRawYield->Draw("same");
1059     TH1F *hRawYieldRefl = ReflectHisto(fRawYield);
1060     hRawYieldRefl->Draw("same");
1061     leg->AddEntry(fRawYield,"Yield extraction","l");
1062   }
1063   if(fCutsEff) {
1064     fCutsEff->SetLineColor(4);
1065     fCutsEff->SetLineWidth(3);
1066     fCutsEff->Draw("same");
1067     TH1F *hCutsEffRefl = ReflectHisto(fCutsEff);
1068     hCutsEffRefl->Draw("same");
1069     leg->AddEntry(fCutsEff,"Cuts efficiency","l");
1070   }
1071   if(fPIDEff) {
1072     fPIDEff->SetLineColor(7);
1073     fPIDEff->SetLineWidth(3);
1074     fPIDEff->Draw("same");
1075     TH1F *hPIDEffRefl = ReflectHisto(fPIDEff);
1076     hPIDEffRefl->Draw("same");
1077     leg->AddEntry(fPIDEff,"PID efficiency","l");
1078   }
1079   if(fMCPtShape) {
1080     Int_t ci = TColor::GetColor("#9933ff");
1081     fMCPtShape->SetLineColor(ci);
1082     //    fMCPtShape->SetLineColor(8);
1083     fMCPtShape->SetLineWidth(3);
1084     fMCPtShape->Draw("same");
1085     TH1F *hMCPtShapeRefl = ReflectHisto(fMCPtShape);
1086     hMCPtShapeRefl->Draw("same");
1087     leg->AddEntry(fMCPtShape,"MC p_{t} shape","l");
1088   }
1089   if(fPartAntipart) {
1090     Int_t ci = TColor::GetColor("#ff6600");
1091     fPartAntipart->SetLineColor(ci);
1092     //    fPartAntipart->SetLineColor(9);
1093     fPartAntipart->SetLineWidth(3);
1094     fPartAntipart->Draw("same");
1095     TH1F *hPartAntipartRefl = ReflectHisto(fPartAntipart);
1096     hPartAntipartRefl->Draw("same");
1097     leg->AddEntry(fPartAntipart,"D = #bar{D}","l");
1098   }
1099
1100
1101   leg->Draw();
1102
1103   cSystErr->SaveAs("RelativeSystematics.eps");
1104
1105   return;
1106 }
1107 //-------------------------------------------------------------------------
1108 TH1F* AliHFSystErr::ReflectHisto(TH1F *hin) const {
1109   //
1110   // Clones and reflects histogram 
1111   // 
1112   TH1F *hout=(TH1F*)hin->Clone("hout");
1113   hout->Scale(-1.);
1114
1115   return hout;
1116 }