Update (Andrea)
[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.05); // 5%
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
386 // Normalization
387   fNorm = new TH1F("fNorm","fNorm",24,0,24);
388   for(Int_t i=1;i<=24;i++) fNorm->SetBinContent(i,0.07); // 10% error on sigmaV0and
389
390   // Branching ratio 
391   fBR = new TH1F("fBR","fBR",24,0,24);
392   for(Int_t i=1;i<=24;i++) fBR->SetBinContent(i,0.04); // 4% PDG2010
393
394   // Tracking efficiency
395   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",24,0,24);
396   for(Int_t i=1;i<=24;i++) fTrackingEff->SetBinContent(i,0.03); // 3% (1% per track)
397
398
399   // Raw yield extraction
400   fRawYield = new TH1F("fRawYield","fRawYield",24,0,24);
401   fRawYield->SetBinContent(1,0.25);
402   fRawYield->SetBinContent(2,0.25);
403   fRawYield->SetBinContent(3,0.25);
404   fRawYield->SetBinContent(4,0.25);
405   fRawYield->SetBinContent(5,0.09);
406   fRawYield->SetBinContent(6,0.09);
407   for(Int_t i=7;i<=12;i++) fRawYield->SetBinContent(i,0.05);  //5 to 10%
408   for(Int_t i=12;i<=24;i++) fRawYield->SetBinContent(i,0.10);  // 15%
409   
410   // Cuts efficiency (from cuts variation)
411   fCutsEff = new TH1F("fCutsEff","fCutsEff",24,0,24);
412   for(Int_t i=1;i<=24;i++) fCutsEff->SetBinContent(i,0.10); // 10%
413
414   // PID efficiency (from PID/noPID)
415   fPIDEff = new TH1F("fPIDEff","fPIDEff",24,0,24);
416   for(Int_t i=1;i<=24;i++) fPIDEff->SetBinContent(i,0.05); // 5%
417   fPIDEff->SetBinContent(1,0.15); // 15%
418   fPIDEff->SetBinContent(2,0.15); // 15%
419   fPIDEff->SetBinContent(3,0.15); // 15%
420   fPIDEff->SetBinContent(4,0.15); // 15%
421   for(Int_t i=12;i<=16;i++) fPIDEff->SetBinContent(i,0.10); // 5%
422
423   // MC dN/dpt  
424   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",24,0,24);
425   for(Int_t i=1;i<=24;i++) fMCPtShape->SetBinContent(i,0);
426   fMCPtShape->SetBinContent(1,0.03);
427   fMCPtShape->SetBinContent(2,0.03);
428
429
430   // particle-antiparticle
431   /*
432   fPartAntipart = new TH1F("fPartAntipart","fPartAntipart",20,0,20);
433   fPartAntipart->SetBinContent(1,1);
434   fPartAntipart->SetBinContent(2,1);
435   fPartAntipart->SetBinContent(3,0.12);
436   for(Int_t i=4;i<=20;i++) fPartAntipart->SetBinContent(i,0.05);   //5 to 12%
437   */
438   return;
439 }
440  
441 //--------------------------------------------------------------------------
442 void AliHFSystErr::InitDplustoKpipi2010PbPb020() {
443   // 
444   // D+->Kpipi syst errors. Responsible: ??
445   //  2010 PbPb sample, 0-20 CC
446   //
447
448  // Normalization
449   fNorm = new TH1F("fNorm","fNorm",20,0,20);
450   for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
451
452   // Branching ratio 
453   fBR = new TH1F("fBR","fBR",20,0,20);
454   for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.04); // 4% PDG2010
455
456   // Tracking efficiency
457   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
458   for(Int_t i=1;i<=20;i++) fTrackingEff->SetBinContent(i,i,0.075+0.005*(Float_t)i); 
459
460
461   // Raw yield extraction
462   fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
463   for(Int_t i=1;i<=20;i++) fRawYield->SetBinContent(i,1);  //5 to 10%
464 //   for(Int_t i=6;i<=7;i++) fRawYield->SetBinContent(i,0.32);
465 //   for(Int_t i=7;i<=9;i++) fRawYield->SetBinContent(i,0.32);
466 //   for(Int_t i=9;i<=12;i++) fRawYield->SetBinContent(i,0.05);//31./156.);
467   fRawYield->SetBinContent(6,0.12);
468   fRawYield->SetBinContent(7,0.17);
469   fRawYield->SetBinContent(8,0.17);
470   fRawYield->SetBinContent(9,0.08);
471   fRawYield->SetBinContent(10,0.08);
472   fRawYield->SetBinContent(11,0.08);
473   fRawYield->SetBinContent(12,0.08);
474
475   // Cuts efficiency (from cuts variation)
476   fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
477   for(Int_t i=1;i<=20;i++) fCutsEff->SetBinContent(i,0.15); // 10%
478
479   // PID efficiency (from PID/noPID)
480   fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
481   for(Int_t i=1;i<=20;i++) fPIDEff->SetBinContent(i,0.05); // 5%
482   fPIDEff->SetBinContent(3,0.13); // 13%
483  
484
485   // MC dN/dpt  (copied from D0 : will update later)
486   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
487   for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,0.);
488   fMCPtShape->SetBinContent(1,0.03);
489   fMCPtShape->SetBinContent(2,0.03);
490
491
492   // particle-antiparticle
493   /*
494   fPartAntipart = new TH1F("fPartAntipart","fPartAntipart",20,0,20);
495   fPartAntipart->SetBinContent(1,1);
496   fPartAntipart->SetBinContent(2,1);
497   fPartAntipart->SetBinContent(3,0.12);
498   for(Int_t i=4;i<=20;i++) fPartAntipart->SetBinContent(i,0.05);   //5 to 12%
499   */
500   return;
501 }
502
503 //--------------------------------------------------------------------------
504 void AliHFSystErr::InitDplustoKpipi2010PbPb4080() {
505   // 
506   // D+->Kpipi syst errors. Responsible: ??
507   //  2010 PbPb sample, 40-80 CC
508   //
509   
510
511  // Normalization
512   fNorm = new TH1F("fNorm","fNorm",20,0,20);
513   for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
514
515   // Branching ratio 
516   fBR = new TH1F("fBR","fBR",20,0,20);
517   for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.04); // 4% PDG2010
518
519   // Tracking efficiency
520   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
521   for(Int_t i=1;i<=20;i++) fTrackingEff->SetBinContent(i,0.03); // 3% (1% per track)
522
523
524   // Raw yield extraction
525   fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
526   fRawYield->SetBinContent(1,1);
527   fRawYield->SetBinContent(2,1);
528   fRawYield->SetBinContent(3,1);
529   fRawYield->SetBinContent(4,1.);
530   fRawYield->SetBinContent(5,4./23.);
531   fRawYield->SetBinContent(6,4./23.);
532   fRawYield->SetBinContent(7,5./20.);
533   fRawYield->SetBinContent(8,5./20.);
534   for(Int_t i=9;i<=12;i++) fRawYield->SetBinContent(i,1./16.);
535   for(Int_t i=13;i<=20;i++) fRawYield->SetBinContent(i,1);  //5 to 10%
536
537   // Cuts efficiency (from cuts variation)
538   fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
539   for(Int_t i=1;i<=20;i++) fCutsEff->SetBinContent(i,0.10); // 10%
540
541   // PID efficiency (from PID/noPID)
542   fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
543   for(Int_t i=1;i<=20;i++) fPIDEff->SetBinContent(i,0.05); // 5%
544   fPIDEff->SetBinContent(3,0.13); // 13%
545  
546
547   // MC dN/dpt  (copied from D0 : will update later)
548   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
549   for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,0);
550   fMCPtShape->SetBinContent(1,0.03);
551   fMCPtShape->SetBinContent(2,0.03);
552
553
554   // particle-antiparticle
555   /*
556   fPartAntipart = new TH1F("fPartAntipart","fPartAntipart",20,0,20);
557   fPartAntipart->SetBinContent(1,1);
558   fPartAntipart->SetBinContent(2,1);
559   fPartAntipart->SetBinContent(3,0.12);
560   for(Int_t i=4;i<=20;i++) fPartAntipart->SetBinContent(i,0.05);   //5 to 12%
561   */
562   return;
563 }
564
565 //--------------------------------------------------------------------------
566 void AliHFSystErr::InitDplustoKpipi2010ppLowEn() {
567
568   // 
569   // D+->Kpipi syst errors. Responsible: R. Bala
570   //  2011 2.76 TeV pp sample
571   //
572   AliInfo(" Settings for D+ --> K pi pi p-p collisions at 2.76 TeV"); 
573
574   // Normalization
575   fNorm = new TH1F("fNorm","fNorm",20,0,20);
576   for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
577
578   // Branching ratio 
579   fBR = new TH1F("fBR","fBR",20,0,20);
580   for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.04); // 4% PDG2010
581
582   // Tracking efficiency
583   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
584   for(Int_t i=1;i<=20;i++) fTrackingEff->SetBinContent(i,0.15); // 3% (1% per track)
585
586   // Raw yield extraction
587   fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
588   fRawYield->SetBinContent(1,1);
589   fRawYield->SetBinContent(2,1);
590   fRawYield->SetBinContent(3,0.20);
591   fRawYield->SetBinContent(4,0.20);
592   for(Int_t i=5;i<=20;i++) fRawYield->SetBinContent(i,0.055);  //5 to 10%
593
594   // Cuts efficiency (from cuts variation)
595   fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
596   for(Int_t i=1;i<=20;i++) fCutsEff->SetBinContent(i,0.10); // 10%
597
598   // PID efficiency (from PID/noPID)
599   fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
600   for(Int_t i=1;i<=20;i++) fPIDEff->SetBinContent(i,0.05); // 5%
601   fPIDEff->SetBinContent(3,0.15); // 13%
602   fPIDEff->SetBinContent(4,0.15); // 13%
603  
604   // MC dN/dpt  (copied from D0 : will update later)
605   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
606   for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,0.);
607   fMCPtShape->SetBinContent(1,0.03);
608   fMCPtShape->SetBinContent(2,0.03);
609
610   return;
611 }
612
613 //--------------------------------------------------------------------------
614 void AliHFSystErr::InitDstartoD0pi2010pp() {
615   // 
616   // D*+->D0pi syst errors. Responsible: A. Grelli, Y. Wang
617   //  2010 pp sample
618   //
619
620  // Normalization
621   fNorm = new TH1F("fNorm","fNorm",24,0,24);
622   for(Int_t i=1;i<=24;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
623
624   // Branching ratio 
625   fBR = new TH1F("fBR","fBR",24,0,24);
626   for(Int_t i=1;i<=24;i++) fBR->SetBinContent(i,0.015); // 1.5% PDG2010
627
628   // Tracking efficiency
629   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",24,0,24);
630   fTrackingEff->SetBinContent(1,1.0);
631   fTrackingEff->SetBinContent(2,0.12);
632   fTrackingEff->SetBinContent(3,0.08);
633   fTrackingEff->SetBinContent(3,0.05);
634   for(Int_t i=4;i<=24;i++) fTrackingEff->SetBinContent(i,0.03); // 3% (1% per track)
635
636
637   // Raw yield extraction
638   fRawYield = new TH1F("fRawYield","fRawYield",24,0,24);
639   fRawYield->SetBinContent(1,1.0);
640   fRawYield->SetBinContent(2,0.09);
641   fRawYield->SetBinContent(3,0.04);
642   fRawYield->SetBinContent(4,0.02);
643   fRawYield->SetBinContent(5,0.03);
644   fRawYield->SetBinContent(6,0.06);
645   fRawYield->SetBinContent(7,0.04);
646   fRawYield->SetBinContent(8,0.11);
647   for(Int_t i=5;i<=24;i++) fRawYield->SetBinContent(i,0.03);  //4%
648   fRawYield->SetBinContent(13,0.09);
649   fRawYield->SetBinContent(14,0.09);
650   fRawYield->SetBinContent(15,0.09);
651   fRawYield->SetBinContent(16,0.09);
652   fRawYield->SetBinContent(17,0.24);
653   fRawYield->SetBinContent(18,0.24);
654   fRawYield->SetBinContent(19,0.24);
655   fRawYield->SetBinContent(20,0.24);
656   fRawYield->SetBinContent(21,0.24);
657   fRawYield->SetBinContent(22,0.24);
658   fRawYield->SetBinContent(23,0.24);
659   fRawYield->SetBinContent(24,0.24);
660
661   // Cuts efficiency (from cuts variation)
662   fCutsEff = new TH1F("fCutsEff","fCutsEff",24,0,24);
663   fCutsEff->SetBinContent(2,0.22);
664   for(Int_t i=3;i<=24;i++) fCutsEff->SetBinContent(i,0.10); // 10%
665
666   // PID efficiency (from PID/noPID)
667   fPIDEff = new TH1F("fPIDEff","fPIDEff",24,0,24);
668   for(Int_t i=1;i<=24;i++) fPIDEff->SetBinContent(i,0.04); // 3%
669  
670
671   // MC dN/dpt  (copied from D0 : will update later)
672   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",24,0,24);
673   for(Int_t i=1;i<=24;i++) fMCPtShape->SetBinContent(i,0);
674   fMCPtShape->SetBinContent(1,0.03);
675   fMCPtShape->SetBinContent(2,0.03);
676
677   return;
678
679
680 }
681 //--------------------------------------------------------------------------
682 void AliHFSystErr::InitDstartoD0pi2010ppLowEn() {
683
684   // 
685   // D+->Kpipi syst errors. Responsible: A. Grelli
686   //  2011 2.76 TeV pp sample
687   //
688   AliInfo(" Settings for D*+ --> D0 pi p-p collisions at 2.76 TeV"); 
689
690 // Normalization
691   fNorm = new TH1F("fNorm","fNorm",20,0,20);
692   for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
693
694   // Branching ratio 
695   fBR = new TH1F("fBR","fBR",20,0,20);
696   for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.015); // 1.5% PDG2010
697
698   // Tracking efficiency
699   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
700   for(Int_t i=1;i<=20;i++) fTrackingEff->SetBinContent(i,0.15); //10% (to be checked!!)
701
702   // Raw yield extraction
703   fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
704   fRawYield->SetBinContent(1,1);
705   fRawYield->SetBinContent(2,1);
706   fRawYield->SetBinContent(3,0.2);
707   fRawYield->SetBinContent(4,0.2);
708   fRawYield->SetBinContent(5,0.08);
709   fRawYield->SetBinContent(6,0.08);
710   fRawYield->SetBinContent(7,0.1);
711   fRawYield->SetBinContent(8,0.1);
712   fRawYield->SetBinContent(9,0.2);
713   fRawYield->SetBinContent(10,0.2);
714   for(Int_t i=9;i<=20;i++) fRawYield->SetBinContent(i,0.065);
715
716   // Cuts efficiency (from cuts variation)
717   fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
718   for(Int_t i=1;i<=20;i++) fCutsEff->SetBinContent(i,0.18);  
719   fCutsEff->SetBinContent(3,0.2);
720   fCutsEff->SetBinContent(4,0.2);
721   fCutsEff->SetBinContent(5,0.1);
722   fCutsEff->SetBinContent(6,0.1);
723   fCutsEff->SetBinContent(7,0.12);
724   fCutsEff->SetBinContent(8,0.12);
725   fCutsEff->SetBinContent(9,0.2);
726   fCutsEff->SetBinContent(10,0.2);
727   fCutsEff->SetBinContent(11,0.2);
728   fCutsEff->SetBinContent(12,0.2);
729
730   // PID efficiency (from PID/noPID)
731   fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
732   for(Int_t i=1;i<=20;i++) fPIDEff->SetBinContent(i,0.05); // 10%
733
734   // MC dN/dpt
735   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
736   for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,0.);
737   fMCPtShape->SetBinContent(1,0.03);
738   fMCPtShape->SetBinContent(2,0.03);
739
740
741   return;
742 }
743 //------------------------------------------------------------------------
744 void AliHFSystErr::InitDstartoD0pi2010PbPb020() {
745   // 
746   // D*+->D0pi syst errors. Responsible: A. Grelli
747   //  2010 PbPb sample, 0-20 CC
748   //
749
750   AliInfo(" Settings for D*+ --> D0pi Pb-Pb collisions at 2.76 TeV - 0-20 centrality - DUMMY"); 
751
752  // Normalization
753   fNorm = new TH1F("fNorm","fNorm",24,0,24);
754   for(Int_t i=1;i<=24;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
755
756   // Branching ratio 
757   fBR = new TH1F("fBR","fBR",24,0,24);
758   for(Int_t i=1;i<=24;i++) fBR->SetBinContent(i,0.015); // 1.5% PDG2010
759
760   // Tracking efficiency
761   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",24,0,24);;
762   for(Int_t i=1;i<=24;i++) fTrackingEff->SetBinContent(i,0.12); // 3% (1% per track)
763
764
765   // Raw yield extraction
766   fRawYield = new TH1F("fRawYield","fRawYield",24,0,24);
767   for(Int_t i=1;i<=24;i++) fRawYield->SetBinContent(i,0.15);  //4%
768  
769   // Cuts efficiency (from cuts variation)
770   fCutsEff = new TH1F("fCutsEff","fCutsEff",24,0,24);
771   for(Int_t i=1;i<=24;i++) fCutsEff->SetBinContent(i,0.10); // 10%
772
773   // PID efficiency (from PID/noPID)
774   fPIDEff = new TH1F("fPIDEff","fPIDEff",24,0,24);
775   for(Int_t i=1;i<=24;i++) fPIDEff->SetBinContent(i,0.04); // 3%
776  
777
778   // MC dN/dpt  (copied from D0 : will update later)
779   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",24,0,24);
780   for(Int_t i=1;i<=24;i++) fMCPtShape->SetBinContent(i,0.);
781   fMCPtShape->SetBinContent(1,0.03);
782   fMCPtShape->SetBinContent(2,0.03);
783
784   return;
785
786 }
787 //-------------------------------------------------------------------------
788 void AliHFSystErr::InitDstartoD0pi2010PbPb2040() {
789   // 
790   // D*+->D0pi syst errors. Responsible: A. Grelli
791   //  2010 PbPb sample, 20-40 CC
792   //
793
794   AliInfo(" Settings for D*+ --> D0pi Pb-Pb collisions at 2.76 TeV - 20-40 centrality - DUMMY"); 
795
796  // Normalization
797   fNorm = new TH1F("fNorm","fNorm",24,0,24);
798   for(Int_t i=1;i<=24;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
799
800   // Branching ratio 
801   fBR = new TH1F("fBR","fBR",24,0,24);
802   for(Int_t i=1;i<=24;i++) fBR->SetBinContent(i,0.015); // 1.5% PDG2010
803
804   // Tracking efficiency
805   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",24,0,24);;
806   for(Int_t i=1;i<=24;i++) fTrackingEff->SetBinContent(i,0.12); // 3% (1% per track)
807
808
809   // Raw yield extraction
810   fRawYield = new TH1F("fRawYield","fRawYield",24,0,24);
811   for(Int_t i=1;i<=24;i++) fRawYield->SetBinContent(i,0.15);  //4%
812  
813   // Cuts efficiency (from cuts variation)
814   fCutsEff = new TH1F("fCutsEff","fCutsEff",24,0,24);
815   for(Int_t i=1;i<=24;i++) fCutsEff->SetBinContent(i,0.10); // 10%
816
817   // PID efficiency (from PID/noPID)
818   fPIDEff = new TH1F("fPIDEff","fPIDEff",24,0,24);
819   for(Int_t i=1;i<=24;i++) fPIDEff->SetBinContent(i,0.04); // 3%
820  
821
822   // MC dN/dpt  (copied from D0 : will update later)
823   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",24,0,24);
824   for(Int_t i=1;i<=24;i++) fMCPtShape->SetBinContent(i,0.);
825   fMCPtShape->SetBinContent(1,0.03);
826   fMCPtShape->SetBinContent(2,0.03);
827
828   return;
829
830 }
831 //--------------------------------------------------------------------------
832 void AliHFSystErr::InitDstartoD0pi2010PbPb4080() {
833   // 
834   // D*+->D0pi syst errors. Responsible: A. Grelli
835   //  2010 PbPb sample, 40-80 CC
836   //
837
838   AliInfo(" Settings for D*+ --> D0pi Pb-Pb collisions at 2.76 TeV - 40-80 centrality - DUMMY"); 
839
840  // Normalization
841   fNorm = new TH1F("fNorm","fNorm",24,0,24);
842   for(Int_t i=1;i<=24;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
843
844   // Branching ratio 
845   fBR = new TH1F("fBR","fBR",24,0,24);
846   for(Int_t i=1;i<=24;i++) fBR->SetBinContent(i,0.015); // 1.5% PDG2010
847
848   // Tracking efficiency
849   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",24,0,24);;
850   for(Int_t i=1;i<=24;i++) fTrackingEff->SetBinContent(i,0.12); // 3% (1% per track)
851
852
853   // Raw yield extraction
854   fRawYield = new TH1F("fRawYield","fRawYield",24,0,24);
855   for(Int_t i=1;i<=24;i++) fRawYield->SetBinContent(i,0.15);  //4%
856  
857   // Cuts efficiency (from cuts variation)
858   fCutsEff = new TH1F("fCutsEff","fCutsEff",24,0,24);
859   for(Int_t i=1;i<=24;i++) fCutsEff->SetBinContent(i,0.10); // 10%
860
861   // PID efficiency (from PID/noPID)
862   fPIDEff = new TH1F("fPIDEff","fPIDEff",24,0,24);
863   for(Int_t i=1;i<=24;i++) fPIDEff->SetBinContent(i,0.04); // 3%
864  
865
866   // MC dN/dpt  (copied from D0 : will update later)
867   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",24,0,24);
868   for(Int_t i=1;i<=24;i++) fMCPtShape->SetBinContent(i,0.);
869   fMCPtShape->SetBinContent(1,0.03);
870   fMCPtShape->SetBinContent(2,0.03);
871
872   return;
873
874 }
875 //--------------------------------------------------------------------------
876 Double_t AliHFSystErr::GetCutsEffErr(Double_t pt) const {
877   // 
878   // Get error
879   //
880
881   Int_t bin=fCutsEff->FindBin(pt);
882
883   return fCutsEff->GetBinContent(bin);
884 }
885 //--------------------------------------------------------------------------
886 Double_t AliHFSystErr::GetMCPtShapeErr(Double_t pt) const {
887   // 
888   // Get error
889   //
890
891   Int_t bin=fMCPtShape->FindBin(pt);
892
893   return fMCPtShape->GetBinContent(bin);
894 }
895 //--------------------------------------------------------------------------
896 Double_t AliHFSystErr::GetSeleEffErr(Double_t pt) const {
897   // 
898   // Get error
899   //
900
901   Double_t err=GetCutsEffErr(pt)*GetCutsEffErr(pt)+GetMCPtShapeErr(pt)*GetMCPtShapeErr(pt);
902
903   return TMath::Sqrt(err);
904 }
905 //--------------------------------------------------------------------------
906 Double_t AliHFSystErr::GetPIDEffErr(Double_t pt) const {
907   // 
908   // Get error
909   //
910
911   Int_t bin=fPIDEff->FindBin(pt);
912
913   return fPIDEff->GetBinContent(bin);
914 }
915 //--------------------------------------------------------------------------
916 Double_t AliHFSystErr::GetTrackingEffErr(Double_t pt) const {
917   // 
918   // Get error
919   //
920
921   Int_t bin=fTrackingEff->FindBin(pt);
922
923   return fTrackingEff->GetBinContent(bin);
924 }
925 //--------------------------------------------------------------------------
926 Double_t AliHFSystErr::GetRawYieldErr(Double_t pt) const {
927   // 
928   // Get error
929   //
930
931   Int_t bin=fRawYield->FindBin(pt);
932
933   return fRawYield->GetBinContent(bin);
934 }
935 //--------------------------------------------------------------------------
936 Double_t AliHFSystErr::GetPartAntipartErr(Double_t pt) const {
937   // 
938   // Get error
939   //
940
941   Int_t bin=fPartAntipart->FindBin(pt);
942
943   return fPartAntipart->GetBinContent(bin);
944 }
945 //--------------------------------------------------------------------------
946 Double_t AliHFSystErr::GetTotalSystErr(Double_t pt,Double_t feeddownErr) const {
947   // 
948   // Get total syst error (except norm. error)
949   //
950
951   Double_t err=0.;
952
953   if(fRawYield) err += GetRawYieldErr(pt)*GetRawYieldErr(pt);
954   if(fTrackingEff) err += GetTrackingEffErr(pt)*GetTrackingEffErr(pt);
955   if(fBR) err += GetBRErr()*GetBRErr();
956   if(fCutsEff) err += GetCutsEffErr(pt)*GetCutsEffErr(pt);
957   if(fPIDEff) err += GetPIDEffErr(pt)*GetPIDEffErr(pt);
958   if(fMCPtShape) err += GetMCPtShapeErr(pt)*GetMCPtShapeErr(pt);
959   if(fPartAntipart) err += GetPartAntipartErr(pt)*GetPartAntipartErr(pt);
960
961   err += feeddownErr*feeddownErr;
962
963   return TMath::Sqrt(err);
964 }
965 //---------------------------------------------------------------------------
966 void AliHFSystErr::DrawErrors(TGraphAsymmErrors *grErrFeeddown) const {
967   //
968   // Draw errors
969   //
970   gStyle->SetOptStat(0);
971
972   TCanvas *cSystErr = new TCanvas("cSystErr","Systematic Errors",300,80,640,500);
973   cSystErr->Range(0.20,-0.5,18.4,0.34);
974   cSystErr->SetRightMargin(0.318);
975   cSystErr->SetFillColor(0);
976
977   TH2F *hFrame = new TH2F("hFrame","Systematic errors; p_{t} [GeV/c]; Relative Error",20,0,20,100,-1,+1);
978   hFrame->SetAxisRange(1.,15.9,"X");
979   hFrame->SetAxisRange(-0.5,0.5,"Y");
980   hFrame->Draw();
981
982   TLegend *leg = new TLegend(0.69,0.44,0.98,0.86,NULL,"brNDC");
983   leg->SetTextSize(0.03601695);
984   leg->SetFillStyle(0);
985   leg->SetBorderSize(0);
986   
987   TH1F *hTotErr=new TH1F("hTotErr","",24,0,24);
988   Int_t nbins = fNorm->GetNbinsX();
989   TGraphAsymmErrors *gTotErr = new TGraphAsymmErrors(nbins);
990   for(Int_t i=1;i<=24;i++) {
991     Double_t pt = hTotErr->GetBinCenter(i);
992     Double_t ptwidth = hTotErr->GetBinWidth(i);
993
994     if(grErrFeeddown) {
995       Double_t x=0., y=0., errxl=0., errxh=0., erryl=0., erryh=0.;
996       Double_t toterryl=0., toterryh=0.;
997       for(Int_t j=0; j<grErrFeeddown->GetN(); j++) {
998         grErrFeeddown->GetPoint(j,x,y);
999         errxh = grErrFeeddown->GetErrorXhigh(j);
1000         errxl = grErrFeeddown->GetErrorXlow(j);
1001         if ( ( (x-errxl) <= pt) && ( (x+errxl) >= pt) ) {
1002           erryh = grErrFeeddown->GetErrorYhigh(j);
1003           erryl = grErrFeeddown->GetErrorYlow(j);
1004         }
1005       }
1006       if (erryl>=1e-3) toterryl = GetTotalSystErr(pt,erryl);
1007       else toterryl = GetTotalSystErr(pt);
1008       if (erryh>=1e-3) toterryh = GetTotalSystErr(pt,erryh);
1009       else toterryh = GetTotalSystErr(pt);
1010
1011       hTotErr->SetBinContent(i,toterryh);
1012       gTotErr->SetPoint(i,pt,0.);
1013       gTotErr->SetPointError(i,ptwidth/2.,ptwidth/2.,toterryl,toterryh); // i, exl, exh, eyl, eyh
1014     }
1015     else {
1016       hTotErr->SetBinContent(i,GetTotalSystErr(pt));
1017       gTotErr->SetPoint(i,pt,0.);
1018       gTotErr->SetPointError(i,ptwidth/2.,ptwidth/2.,GetTotalSystErr(pt),GetTotalSystErr(pt)); // i, exl, exh, eyl, eyh
1019     }
1020
1021   }
1022   gTotErr->SetLineColor(kBlack);
1023   gTotErr->SetFillColor(kRed);
1024   gTotErr->SetFillStyle(3002);
1025   gTotErr->Draw("2");
1026   leg->AddEntry(gTotErr,"Total (excl. norm.)","f");
1027 //   hTotErr->SetLineColor(1);
1028 //   hTotErr->SetLineWidth(3);
1029 //   hTotErr->Draw("same");
1030 //   leg->AddEntry(hTotErr,"Total (excl. norm.)","l");
1031   
1032
1033   fNorm->SetFillColor(1);
1034   fNorm->SetFillStyle(3002);
1035   //fNorm->Draw("same");
1036   //TH1F *hNormRefl = ReflectHisto(fNorm);
1037   //hNormRefl->Draw("same");
1038   leg->AddEntry(fNorm,"Normalization (10%)","");
1039
1040   if(grErrFeeddown) {
1041     grErrFeeddown->SetFillColor(kTeal-8);
1042     grErrFeeddown->SetFillStyle(3001);
1043     grErrFeeddown->Draw("2");
1044     leg->AddEntry(grErrFeeddown,"Feed-down from B","f");
1045   }
1046   if(fTrackingEff) {
1047     fTrackingEff->SetFillColor(4);
1048     fTrackingEff->SetFillStyle(3006);
1049     fTrackingEff->Draw("same");
1050     TH1F *hTrackingEffRefl = ReflectHisto(fTrackingEff);
1051     hTrackingEffRefl->Draw("same");
1052     leg->AddEntry(fTrackingEff,"Tracking efficiency","f");
1053   }
1054   if(fBR) {
1055     fBR->SetFillColor(6);
1056     fBR->SetFillStyle(3005);
1057     //fBR->SetFillStyle(3020);
1058     fBR->Draw("same");
1059     TH1F *hBRRefl = ReflectHisto(fBR);
1060     hBRRefl->Draw("same");
1061     leg->AddEntry(fBR,"Branching ratio","f");
1062   }
1063   if(fRawYield) {
1064     Int_t ci;   // for color index setting
1065     ci = TColor::GetColor("#00cc00");
1066     fRawYield->SetLineColor(ci);
1067     //    fRawYield->SetLineColor(3);
1068     fRawYield->SetLineWidth(3);
1069     fRawYield->Draw("same");
1070     TH1F *hRawYieldRefl = ReflectHisto(fRawYield);
1071     hRawYieldRefl->Draw("same");
1072     leg->AddEntry(fRawYield,"Yield extraction","l");
1073   }
1074   if(fCutsEff) {
1075     fCutsEff->SetLineColor(4);
1076     fCutsEff->SetLineWidth(3);
1077     fCutsEff->Draw("same");
1078     TH1F *hCutsEffRefl = ReflectHisto(fCutsEff);
1079     hCutsEffRefl->Draw("same");
1080     leg->AddEntry(fCutsEff,"Cuts efficiency","l");
1081   }
1082   if(fPIDEff) {
1083     fPIDEff->SetLineColor(7);
1084     fPIDEff->SetLineWidth(3);
1085     fPIDEff->Draw("same");
1086     TH1F *hPIDEffRefl = ReflectHisto(fPIDEff);
1087     hPIDEffRefl->Draw("same");
1088     leg->AddEntry(fPIDEff,"PID efficiency","l");
1089   }
1090   if(fMCPtShape) {
1091     Int_t ci = TColor::GetColor("#9933ff");
1092     fMCPtShape->SetLineColor(ci);
1093     //    fMCPtShape->SetLineColor(8);
1094     fMCPtShape->SetLineWidth(3);
1095     fMCPtShape->Draw("same");
1096     TH1F *hMCPtShapeRefl = ReflectHisto(fMCPtShape);
1097     hMCPtShapeRefl->Draw("same");
1098     leg->AddEntry(fMCPtShape,"MC p_{t} shape","l");
1099   }
1100   if(fPartAntipart) {
1101     Int_t ci = TColor::GetColor("#ff6600");
1102     fPartAntipart->SetLineColor(ci);
1103     //    fPartAntipart->SetLineColor(9);
1104     fPartAntipart->SetLineWidth(3);
1105     fPartAntipart->Draw("same");
1106     TH1F *hPartAntipartRefl = ReflectHisto(fPartAntipart);
1107     hPartAntipartRefl->Draw("same");
1108     leg->AddEntry(fPartAntipart,"D = #bar{D}","l");
1109   }
1110
1111
1112   leg->Draw();
1113
1114   cSystErr->SaveAs("RelativeSystematics.eps");
1115
1116   return;
1117 }
1118 //-------------------------------------------------------------------------
1119 TH1F* AliHFSystErr::ReflectHisto(TH1F *hin) const {
1120   //
1121   // Clones and reflects histogram 
1122   // 
1123   TH1F *hout=(TH1F*)hin->Clone("hout");
1124   hout->Scale(-1.);
1125
1126   return hout;
1127 }