Updated normalization syst unc
[u/mrichter/AliRoot.git] / PWGHF / 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   case 4: // D+s->KKpi
134     if (fCollisionType==0) InitDstoKKpi2010pp();
135     else AliError("Not yet implemented");
136     break;
137     
138   default:
139     printf("Invalid decay type: %d\n",decay);
140     break;
141   }
142
143 }
144   
145 //--------------------------------------------------------------------------
146 void AliHFSystErr::InitD0toKpi2010pp() {
147   // 
148   // D0->Kpi syst errors. Responsible: A. Rossi
149   //   2010 pp sample
150   //
151
152   // Normalization
153   fNorm = new TH1F("fNorm","fNorm",24,0,24);
154   for(Int_t i=1;i<=24;i++) fNorm->SetBinContent(i,0.04); // 4% error on sigmaV0and
155
156   // Branching ratio 
157   fBR = new TH1F("fBR","fBR",24,0,24);
158   for(Int_t i=1;i<=24;i++) fBR->SetBinContent(i,0.012); // 1.2% PDG2010
159
160   // Tracking efficiency
161   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",24,0,24);
162   for(Int_t i=1;i<=24;i++) fTrackingEff->SetBinContent(i,0.08); // 8% (4% per track)
163
164   // Raw yield extraction
165   fRawYield = new TH1F("fRawYield","fRawYield",24,0,24);
166   fRawYield->SetBinContent(1,1);
167   fRawYield->SetBinContent(2,0.22);
168   fRawYield->SetBinContent(3,0.1);
169   for(Int_t i=4;i<=7;i++) fRawYield->SetBinContent(i,0.04);
170   for(Int_t i=8;i<=12;i++) fRawYield->SetBinContent(i,0.07);
171   for(Int_t i=13;i<=16;i++) fRawYield->SetBinContent(i,0.10);
172   for(Int_t i=17;i<=24;i++) fRawYield->SetBinContent(i,1);
173
174   // Cuts efficiency (from cuts variation)
175   fCutsEff = new TH1F("fCutsEff","fCutsEff",24,0,24);
176   for(Int_t i=1;i<=24;i++) fCutsEff->SetBinContent(i,0.10); // 10%
177
178   // PID efficiency (from PID/noPID)
179   fPIDEff = new TH1F("fPIDEff","fPIDEff",24,0,24);
180   for(Int_t i=1;i<=24;i++) fPIDEff->SetBinContent(i,0.03); // 3%
181   fPIDEff->SetBinContent(2,0.05); // 5%
182
183   // MC dN/dpt
184   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",24,0,24);
185   for(Int_t i=1;i<=24;i++) fMCPtShape->SetBinContent(i,0);
186   fMCPtShape->SetBinContent(1,0.03);
187   fMCPtShape->SetBinContent(2,0.03);
188
189   // particle-antiparticle
190   //  fPartAntipart = new TH1F("fPartAntipart","fPartAntipart",24,0,24);
191   //  fPartAntipart->SetBinContent(1,1); 
192   //  for(Int_t i=2;i<=24;i++) fPartAntipart->SetBinContent(i,0.05);
193   
194   return;
195 }
196 //--------------------------------------------------------------------------
197 void AliHFSystErr::InitD0toKpi2010PbPb020() {
198   // 
199   // D0->Kpi syst errors. Responsible: A. Rossi
200   //   2010 PbPb sample, 0-20 CC
201   //
202
203   // Normalization
204   fNorm = new TH1F("fNorm","fNorm",20,0,20);
205   for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.05); // TAA and pp norm
206
207   // Branching ratio 
208   fBR = new TH1F("fBR","fBR",20,0,20);
209   for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.012); // 1.2% PDG2010
210
211   // Tracking efficiency
212   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
213   for(Int_t i=1;i<=20;i++) fTrackingEff->SetBinContent(i,0.10);// Jacek, 5% per track
214
215   // Raw yield extraction
216   fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
217   fRawYield->SetBinContent(1,0);
218   fRawYield->SetBinContent(2,0);
219   fRawYield->SetBinContent(3,0.08);
220   for(Int_t i=4;i<=12;i++) fRawYield->SetBinContent(i,0.06);
221   for(Int_t i=13;i<=16;i++) fRawYield->SetBinContent(i,0.10);
222
223   // Cuts efficiency (from cuts variation)
224   fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
225   fCutsEff->SetBinContent(1,0.);
226   fCutsEff->SetBinContent(2,0.);
227   fCutsEff->SetBinContent(3,0.13);
228   fCutsEff->SetBinContent(4,0.11);  
229   for(Int_t i=5;i<=16;i++) fCutsEff->SetBinContent(i,0.10);
230   for(Int_t i=17;i<=20;i++) fCutsEff->SetBinContent(i,0.);
231
232   // PID efficiency (from PID/noPID)
233   fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
234   for(Int_t i=3;i<=16;i++) fPIDEff->SetBinContent(i,0.05);
235
236   // MC dN/dpt
237   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
238   for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,0.01);
239   fMCPtShape->SetBinContent(3,0.04);
240   fMCPtShape->SetBinContent(4,0.02);
241   for(Int_t i=13;i<=16;i++) fMCPtShape->SetBinContent(i,0.03); 
242
243 //   // particle-antiparticle
244 //   fPartAntipart = new TH1F("fPartAntipart","fPartAntipart",20,0,20);
245 //   for(Int_t i=3;i<=12;i++) fPartAntipart->SetBinContent(i,0.05);
246 //   fPartAntipart->SetBinContent(3,0.10);
247 //   fPartAntipart->SetBinContent(4,0.10);
248 //   fPartAntipart->SetBinContent(7,0.10);
249 //   fPartAntipart->SetBinContent(8,0.10);
250   
251   return;
252 }
253 //--------------------------------------------------------------------------
254 void AliHFSystErr::InitD0toKpi2010PbPb4080() {
255   // 
256   // D0->Kpi syst errors. Responsible: A. Rossi
257   //   2010 PbPb sample, 40-80 CC
258   //
259
260   // Normalization
261   fNorm = new TH1F("fNorm","fNorm",20,0,20);
262   for(Int_t i=1;i<=24;i++) fNorm->SetBinContent(i,0.07); // TAA and pp norm
263
264   // Branching ratio 
265   fBR = new TH1F("fBR","fBR",20,0,20);
266   for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.012); // 1.2% PDG2010
267
268   // Tracking efficiency
269   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
270   for(Int_t i=1;i<=20;i++) fTrackingEff->SetBinContent(i,0.10); // Jacek, 5% per track
271
272
273   // Raw yield extraction
274   fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
275   fRawYield->SetBinContent(1,0);
276   fRawYield->SetBinContent(2,0);
277   for(Int_t i=3;i<=16;i++) fRawYield->SetBinContent(i,0.05);
278   //for(Int_t i=13;i<=20;i++) fRawYield->SetBinContent(i,0);
279
280   // Cuts efficiency (from cuts variation)
281   fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
282   fCutsEff->SetBinContent(1,0.);
283   fCutsEff->SetBinContent(2,0.);
284   fCutsEff->SetBinContent(3,0.13);
285   fCutsEff->SetBinContent(4,0.11);  
286   for(Int_t i=5;i<=16;i++) fCutsEff->SetBinContent(i,0.10);
287   for(Int_t i=17;i<=20;i++) fCutsEff->SetBinContent(i,0.);
288
289   // PID efficiency (from PID/noPID)
290   fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
291 //   for(Int_t i=3;i<=6;i++) fPIDEff->SetBinContent(i,0.10);
292 //   for(Int_t i=7;i<=16;i++) fPIDEff->SetBinContent(i,0.05);
293   for(Int_t i=3;i<=16;i++) fPIDEff->SetBinContent(i,0.05);
294
295   // MC dN/dpt
296   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
297   for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,0.01);
298   //  fMCPtShape->SetBinContent(3,0.04); Not set for peripherals (Raa Vs pt is flat)
299   //  fMCPtShape->SetBinContent(4,0.02);
300   for(Int_t i=13;i<=16;i++) fMCPtShape->SetBinContent(i,0.03); 
301   
302   //   // particle-antiparticle
303   //   fPartAntipart = new TH1F("fPartAntipart","fPartAntipart",20,0,20);
304   //   for(Int_t i=3;i<=12;i++) fPartAntipart->SetBinContent(i,0.05);
305   
306   return;
307 }
308
309 //--------------------------------------------------------------------------
310 void AliHFSystErr::InitD0toKpi2010ppLowEn() {
311   // 
312   // D0->Kpi syst errors. Low energy run
313   //   2011 2.76 TeV pp sample
314   //
315   AliInfo(" Settings for D0 --> K pi, p-p collisions at 2.76 TeV"); 
316
317   // Normalization
318   fNorm = new TH1F("fNorm","fNorm",20,0,20);
319   for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.031); // 4% error on sigmaV0and
320
321   // Branching ratio 
322   fBR = new TH1F("fBR","fBR",20,0,20);
323   for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.012); // 1.2% PDG2010
324
325   // Tracking efficiency
326   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
327   for(Int_t i=1;i<=20;i++) fTrackingEff->SetBinContent(i,0.12); //10% (5% per track)
328
329   // Raw yield extraction
330   fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
331   fRawYield->SetBinContent(1,1);
332   for(Int_t i=1;i<=20;i++) fRawYield->SetBinContent(i,0.15);
333
334   // Cuts efficiency (from cuts variation)
335   fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
336   for(Int_t i=1;i<=20;i++) fCutsEff->SetBinContent(i,0.10); // 10% 
337   fCutsEff->SetBinContent(2,0.20);
338   for(Int_t i=7;i<=20;i++) fCutsEff->SetBinContent(i,0.15); // 10% 
339
340
341   // PID efficiency (from PID/noPID)
342   fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
343   for(Int_t i=1;i<=20;i++) fPIDEff->SetBinContent(i,0.15); // 10%
344   //  fPIDEff->SetBinContent(2,0.20);
345   for(Int_t i=7;i<=20;i++) fPIDEff->SetBinContent(i,0.05); // 10%
346
347   // MC dN/dpt
348   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
349   for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,0);
350 //   fMCPtShape->SetBinContent(1,0.03);
351 //   fMCPtShape->SetBinContent(2,0.03);
352
353 //   // particle-antiparticle
354 //   fPartAntipart = new TH1F("fPartAntipart","fPartAntipart",20,0,20);
355 //   fPartAntipart->SetBinContent(1,1);
356 //   fPartAntipart->SetBinContent(2,1);
357 //   for(Int_t i=3;i<=6;i++) fPartAntipart->SetBinContent(i,0.08);
358 //   for(Int_t i=1;i<=20;i++) fPartAntipart->SetBinContent(i,0.);
359   
360   return;
361 }
362
363 //--------------------------------------------------------------------------
364 void AliHFSystErr::InitDplustoKpipi2010pp() {
365   // 
366   // D+->Kpipi syst errors. Responsible: R. Bala
367   //  2010 pp sample
368   //
369
370
371 // Normalization
372   fNorm = new TH1F("fNorm","fNorm",24,0,24);
373   for(Int_t i=1;i<=24;i++) fNorm->SetBinContent(i,0.04); // 4% error on sigmaV0and
374
375   // Branching ratio 
376   fBR = new TH1F("fBR","fBR",24,0,24);
377   for(Int_t i=1;i<=24;i++) fBR->SetBinContent(i,0.021); // 2.1% PDG2010
378
379   // Tracking efficiency
380   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",24,0,24);
381   for(Int_t i=1;i<=24;i++) fTrackingEff->SetBinContent(i,0.12); // 12% (4% per track)
382
383
384   // Raw yield extraction
385   fRawYield = new TH1F("fRawYield","fRawYield",24,0,24);
386   fRawYield->SetBinContent(1,1);
387   fRawYield->SetBinContent(2,0.25);
388   fRawYield->SetBinContent(3,0.25);
389   fRawYield->SetBinContent(4,0.20);
390   fRawYield->SetBinContent(5,0.09);
391   fRawYield->SetBinContent(6,0.09);
392   for(Int_t i=7;i<=12;i++) fRawYield->SetBinContent(i,0.05);  
393   for(Int_t i=12;i<=24;i++) fRawYield->SetBinContent(i,0.10);  
394   
395   // Cuts efficiency (from cuts variation)
396   fCutsEff = new TH1F("fCutsEff","fCutsEff",24,0,24);
397   for(Int_t i=1;i<=24;i++) fCutsEff->SetBinContent(i,0.10); // 10%
398
399   // PID efficiency (from PID/noPID)
400   fPIDEff = new TH1F("fPIDEff","fPIDEff",24,0,24);
401   for(Int_t i=1;i<=24;i++) fPIDEff->SetBinContent(i,0.05); // 5%
402   fPIDEff->SetBinContent(1,0.15); // 15%
403   fPIDEff->SetBinContent(2,0.15); // 15%
404   fPIDEff->SetBinContent(3,0.15); // 15%
405   fPIDEff->SetBinContent(4,0.15); // 15%
406   for(Int_t i=12;i<=16;i++) fPIDEff->SetBinContent(i,0.10); // 5%
407
408   // MC dN/dpt  
409   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",24,0,24);
410   for(Int_t i=1;i<=24;i++) fMCPtShape->SetBinContent(i,0);
411   fMCPtShape->SetBinContent(1,0.03);
412   fMCPtShape->SetBinContent(2,0.03);
413
414
415   // particle-antiparticle
416   /*
417   fPartAntipart = new TH1F("fPartAntipart","fPartAntipart",20,0,20);
418   fPartAntipart->SetBinContent(1,1);
419   fPartAntipart->SetBinContent(2,1);
420   fPartAntipart->SetBinContent(3,0.12);
421   for(Int_t i=4;i<=20;i++) fPartAntipart->SetBinContent(i,0.05);   //5 to 12%
422   */
423   return;
424 }
425  
426 //--------------------------------------------------------------------------
427 void AliHFSystErr::InitDstoKKpi2010pp() {
428   // 
429   // D+s->KKpi syst errors. Responsible: G.M. Innocenti
430   //  2010 pp sample
431   //
432
433
434 // Normalization
435   fNorm = new TH1F("fNorm","fNorm",12,0,12);
436   for(Int_t i=1;i<=12;i++) fNorm->SetBinContent(i,0.07); // 7% error on sigmaV0and
437
438   // Branching ratio 
439   fBR = new TH1F("fBR","fBR",12,0,12);
440   for(Int_t i=1;i<=12;i++) fBR->SetBinContent(i,0.05); // 5% PDG2010
441
442   // Tracking efficiency
443   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",12,0,12);
444   for(Int_t i=1;i<=12;i++) fTrackingEff->SetBinContent(i,0.03); // 3% (1% per track)
445
446
447   // Raw yield extraction
448   fRawYield = new TH1F("fRawYield","fRawYield",12,0,12);
449   fRawYield->SetBinContent(1,1);
450   fRawYield->SetBinContent(2,1);
451   fRawYield->SetBinContent(3,0.20);
452   fRawYield->SetBinContent(4,0.20);
453   fRawYield->SetBinContent(5,0.15);
454   fRawYield->SetBinContent(6,0.15);
455   fRawYield->SetBinContent(7,0.15);
456   fRawYield->SetBinContent(8,0.15);
457   fRawYield->SetBinContent(9,0.20);
458   fRawYield->SetBinContent(10,0.20);
459   fRawYield->SetBinContent(11,0.20);
460   fRawYield->SetBinContent(12,0.20);
461   
462   // Cuts efficiency (from cuts variation)
463   fCutsEff = new TH1F("fCutsEff","fCutsEff",12,0,12);
464   for(Int_t i=1;i<=12;i++) fCutsEff->SetBinContent(i,0.15); // 15%
465
466   // PID efficiency (from PID/noPID)
467   fPIDEff = new TH1F("fPIDEff","fPIDEff",12,0,12);
468   for(Int_t i=1;i<=12;i++) fPIDEff->SetBinContent(i,0.07); // 7%
469
470   // MC dN/dpt  (copied from D0 : will update later)
471   //fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",12,0,12);
472   //for(Int_t i=1;i<=12;i++) fMCPtShape->SetBinContent(i,(Float_t)i*0.006);
473
474
475   // particle-antiparticle
476   /*
477   fPartAntipart = new TH1F("fPartAntipart","fPartAntipart",12,0,12);
478   fPartAntipart->SetBinContent(1,1);
479   fPartAntipart->SetBinContent(2,1);
480   fPartAntipart->SetBinContent(3,0.12);
481   for(Int_t i=4;i<=12;i++) fPartAntipart->SetBinContent(i,0.05);   //5 to 12%
482   */
483   return;
484 }
485    
486  
487 //--------------------------------------------------------------------------
488 void AliHFSystErr::InitDplustoKpipi2010PbPb020() {
489   // 
490   // D+->Kpipi syst errors. Responsible: ??
491   //  2010 PbPb sample, 0-20 CC
492   //
493
494  // Normalization
495   fNorm = new TH1F("fNorm","fNorm",20,0,20);
496   for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.05); // TAA and pp norm
497
498   // Branching ratio 
499   fBR = new TH1F("fBR","fBR",20,0,20);
500   for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.021); // 2.1% PDG2010
501
502   // Tracking efficiency
503   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
504   for(Int_t i=1;i<=20;i++) fTrackingEff->SetBinContent(i,0.15); // Jacek, 5% per track
505
506   // Raw yield extraction
507   fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
508   for(Int_t i=1;i<=20;i++) fRawYield->SetBinContent(i,.10);  //5 to 10%
509   // fRawYield->SetBinContent(5,0.23);
510   //fRawYield->SetBinContent(6,0.23);
511   fRawYield->SetBinContent(7,0.20);
512   fRawYield->SetBinContent(8,0.20);
513   fRawYield->SetBinContent(9,0.15);
514   fRawYield->SetBinContent(10,0.15);
515   fRawYield->SetBinContent(11,0.15);
516   fRawYield->SetBinContent(12,0.15);
517
518   // Cuts efficiency (from cuts variation)
519   fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
520   for(Int_t i=1;i<=20;i++) fCutsEff->SetBinContent(i,0.15); // 10%
521
522   // PID efficiency (from PID/noPID)
523   fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
524   for(Int_t i=1;i<=20;i++) fPIDEff->SetBinContent(i,0.05); // 5%
525
526   // MC dN/dpt  (2/2/2012)
527   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
528   for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,0.);
529   for(Int_t iBin=7; iBin<=8; iBin++) fMCPtShape->SetBinContent(iBin,0.01);
530   for(Int_t iBin=9; iBin<=12; iBin++) fMCPtShape->SetBinContent(iBin,0.05);
531   for(Int_t iBin=13; iBin<=16; iBin++) fMCPtShape->SetBinContent(iBin,0.05);
532
533
534   // particle-antiparticle
535   /*
536   fPartAntipart = new TH1F("fPartAntipart","fPartAntipart",20,0,20);
537   fPartAntipart->SetBinContent(1,1);
538   fPartAntipart->SetBinContent(2,1);
539   fPartAntipart->SetBinContent(3,0.12);
540   for(Int_t i=4;i<=20;i++) fPartAntipart->SetBinContent(i,0.05);   //5 to 12%
541   */
542
543   return;
544 }
545
546 //--------------------------------------------------------------------------
547 void AliHFSystErr::InitDplustoKpipi2010PbPb4080() {
548   // 
549   // D+->Kpipi syst errors. Responsible: ??
550   //  2010 PbPb sample, 40-80 CC
551   //
552   
553
554  // Normalization
555   fNorm = new TH1F("fNorm","fNorm",20,0,20);
556   for(Int_t i=1;i<=24;i++) fNorm->SetBinContent(i,0.07); // TAA and pp norm
557
558   // Branching ratio 
559   fBR = new TH1F("fBR","fBR",20,0,20);
560   for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.021); // 2.1% 
561
562   // Tracking efficiency
563   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
564   for(Int_t i=1;i<=20;i++) fTrackingEff->SetBinContent(i,0.15); // Jacek, 5% per track
565
566
567   // Raw yield extraction
568   fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
569   fRawYield->SetBinContent(1,1);
570   fRawYield->SetBinContent(2,1);
571   fRawYield->SetBinContent(3,1);
572   fRawYield->SetBinContent(4,0.15);
573   fRawYield->SetBinContent(5,0.05);
574   fRawYield->SetBinContent(6,0.05);
575   fRawYield->SetBinContent(7,0.15);
576   fRawYield->SetBinContent(8,0.15);
577   for(Int_t i=9;i<=12;i++) fRawYield->SetBinContent(i,0.15);
578   for(Int_t i=13;i<=20;i++) fRawYield->SetBinContent(i,1);  //5 to 10%
579
580   // Cuts efficiency (from cuts variation)
581   fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
582   for(Int_t i=1;i<=20;i++) fCutsEff->SetBinContent(i,0.10); // 10%
583
584   // PID efficiency (from PID/noPID)
585   fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
586   for(Int_t i=1;i<=20;i++) fPIDEff->SetBinContent(i,0.05); // 5%
587   fPIDEff->SetBinContent(3,0.13); // 13%
588  
589
590   // MC dN/dpt  (2/2/2012)
591   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
592   for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,0);
593   for(Int_t iBin=4; iBin<=8; iBin++) fMCPtShape->SetBinContent(iBin,0.01);
594   for(Int_t iBin=9; iBin<=12; iBin++) fMCPtShape->SetBinContent(iBin,0.03);
595   for(Int_t iBin=13; iBin<=16; iBin++) fMCPtShape->SetBinContent(iBin,0.03);
596
597
598   // particle-antiparticle
599   /*
600   fPartAntipart = new TH1F("fPartAntipart","fPartAntipart",20,0,20);
601   fPartAntipart->SetBinContent(1,1);
602   fPartAntipart->SetBinContent(2,1);
603   fPartAntipart->SetBinContent(3,0.12);
604   for(Int_t i=4;i<=20;i++) fPartAntipart->SetBinContent(i,0.05);   //5 to 12%
605   */
606   return;
607 }
608
609 //--------------------------------------------------------------------------
610 void AliHFSystErr::InitDplustoKpipi2010ppLowEn() {
611
612   // 
613   // D+->Kpipi syst errors. Responsible: R. Bala
614   //  2011 2.76 TeV pp sample
615   //
616   AliInfo(" Settings for D+ --> K pi pi p-p collisions at 2.76 TeV"); 
617
618   // Normalization
619   fNorm = new TH1F("fNorm","fNorm",20,0,20);
620   for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.031); // 10% error on sigmaV0and
621
622   // Branching ratio 
623   fBR = new TH1F("fBR","fBR",20,0,20);
624   for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.021); // 2.1% PDG2010
625
626   // Tracking efficiency
627   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
628   for(Int_t i=1;i<=20;i++) fTrackingEff->SetBinContent(i,0.18); // 3% (1% per track)
629
630   // Raw yield extraction
631   fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
632   fRawYield->SetBinContent(1,1);
633   fRawYield->SetBinContent(2,1);
634   for(Int_t i=3;i<=6;i++) fRawYield->SetBinContent(i,0.10);  //5 to 10%
635   fRawYield->SetBinContent(7,0.15);
636   fRawYield->SetBinContent(8,0.15); 
637   for(Int_t i=9;i<=20;i++) fRawYield->SetBinContent(i,0.055);  //5 to 10%
638
639   // Cuts efficiency (from cuts variation)
640   fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
641   for(Int_t i=1;i<=20;i++) fCutsEff->SetBinContent(i,0.15); // 10%
642
643   // PID efficiency (from PID/noPID)
644   fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
645   for(Int_t i=1;i<=20;i++) fPIDEff->SetBinContent(i,0.05); // 5%
646   fPIDEff->SetBinContent(3,0.10); // 13%
647   fPIDEff->SetBinContent(4,0.10); // 13%
648  
649   // MC dN/dpt  (copied from D0 : will update later)
650   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
651   for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,0.);
652   fMCPtShape->SetBinContent(1,0.03);
653   fMCPtShape->SetBinContent(2,0.03);
654
655   return;
656 }
657
658 //--------------------------------------------------------------------------
659 void AliHFSystErr::InitDstartoD0pi2010pp() {
660   // 
661   // D*+->D0pi syst errors. Responsible: A. Grelli, Y. Wang
662   //  2010 pp sample
663   //
664
665  // Normalization
666   fNorm = new TH1F("fNorm","fNorm",24,0,24);
667   for(Int_t i=1;i<=24;i++) fNorm->SetBinContent(i,0.04); // 4% error on sigmaV0and
668
669   // Branching ratio 
670   fBR = new TH1F("fBR","fBR",24,0,24);
671   for(Int_t i=1;i<=24;i++) fBR->SetBinContent(i,0.015); // 1.5% PDG2010
672
673   // Tracking efficiency
674   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",24,0,24);
675   fTrackingEff->SetBinContent(1,1.0);
676   fTrackingEff->SetBinContent(2,0.13); // 10% (ITSsa) \oplus 8% (4% per ITSTPC track)
677   fTrackingEff->SetBinContent(3,0.12);
678   fTrackingEff->SetBinContent(3,0.12);
679   for(Int_t i=4;i<=24;i++) fTrackingEff->SetBinContent(i,0.12); // 12% (4% per track)
680
681
682   // Raw yield extraction
683   fRawYield = new TH1F("fRawYield","fRawYield",24,0,24);
684   fRawYield->SetBinContent(1,1.0);
685   fRawYield->SetBinContent(2,0.10);
686   fRawYield->SetBinContent(3,0.04);
687   fRawYield->SetBinContent(4,0.03);
688   fRawYield->SetBinContent(5,0.03);
689   fRawYield->SetBinContent(6,0.05);
690   fRawYield->SetBinContent(7,0.05);
691   fRawYield->SetBinContent(8,0.05);
692   for(Int_t i=9;i<=12;i++) fRawYield->SetBinContent(i,0.04);  //4%
693   for(Int_t i=13;i<=16;i++) fRawYield->SetBinContent(i,0.09);  //4%
694   for(Int_t i=17;i<=24;i++) fRawYield->SetBinContent(i,0.2);  //4%
695
696   // Cuts efficiency (from cuts variation)
697   fCutsEff = new TH1F("fCutsEff","fCutsEff",24,0,24);
698   fCutsEff->SetBinContent(2,0.22);
699   for(Int_t i=3;i<=24;i++) fCutsEff->SetBinContent(i,0.10); // 10%
700
701   // PID efficiency (from PID/noPID)
702   fPIDEff = new TH1F("fPIDEff","fPIDEff",24,0,24);
703   for(Int_t i=1;i<=24;i++) fPIDEff->SetBinContent(i,0.04); // 3%
704  
705
706   // MC dN/dpt  (copied from D0 : will update later)
707   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",24,0,24);
708   for(Int_t i=1;i<=24;i++) fMCPtShape->SetBinContent(i,0);
709   fMCPtShape->SetBinContent(1,0.03);
710   fMCPtShape->SetBinContent(2,0.03);
711
712   return;
713
714
715 }
716 //--------------------------------------------------------------------------
717 void AliHFSystErr::InitDstartoD0pi2010ppLowEn() {
718
719   // 
720   // D+->Kpipi syst errors. Responsible: A. Grelli
721   //  2011 2.76 TeV pp sample
722   //
723   AliInfo(" Settings for D*+ --> D0 pi p-p collisions at 2.76 TeV"); 
724
725 // Normalization
726   fNorm = new TH1F("fNorm","fNorm",20,0,20);
727   for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.031); // 10% error on sigmaV0and
728
729   // Branching ratio 
730   fBR = new TH1F("fBR","fBR",20,0,20);
731   for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.015); // 1.5% PDG2010
732
733   // Tracking efficiency
734   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
735   for(Int_t i=1;i<=20;i++) fTrackingEff->SetBinContent(i,0.18); //10% (to be checked!!)
736
737   // Raw yield extraction
738   fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
739   fRawYield->SetBinContent(1,1);
740   fRawYield->SetBinContent(2,1);
741   fRawYield->SetBinContent(3,0.14);
742   fRawYield->SetBinContent(4,0.14);
743   fRawYield->SetBinContent(5,0.12);
744   fRawYield->SetBinContent(6,0.12);
745   fRawYield->SetBinContent(7,0.06);
746   fRawYield->SetBinContent(8,0.06);
747   fRawYield->SetBinContent(9,0.08);
748   fRawYield->SetBinContent(10,0.08);
749   fRawYield->SetBinContent(11,0.08);
750   fRawYield->SetBinContent(12,0.08);
751   for(Int_t i=9;i<=20;i++) fRawYield->SetBinContent(i,0.065);
752
753   // Cuts efficiency (from cuts variation)
754   fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
755   for(Int_t i=1;i<=20;i++) fCutsEff->SetBinContent(i,0.10);  
756   fCutsEff->SetBinContent(3,0.15);
757   fCutsEff->SetBinContent(4,0.15);
758   fCutsEff->SetBinContent(5,0.15);
759   fCutsEff->SetBinContent(6,0.15);
760   fCutsEff->SetBinContent(7,0.10);
761   fCutsEff->SetBinContent(8,0.10);
762   fCutsEff->SetBinContent(9,0.10);
763   fCutsEff->SetBinContent(10,0.10);
764   fCutsEff->SetBinContent(11,0.10);
765   fCutsEff->SetBinContent(12,0.10);
766
767   // PID efficiency (from PID/noPID)
768   fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
769   for(Int_t i=1;i<=20;i++) fPIDEff->SetBinContent(i,0.05); // 10%
770
771   // MC dN/dpt
772   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
773   for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,0.);
774   fMCPtShape->SetBinContent(1,0.03);
775   fMCPtShape->SetBinContent(2,0.03);
776
777   return;
778 }
779
780 //------------------------------------------------------------------------
781 void AliHFSystErr::InitDstartoD0pi2010PbPb020() {
782   // 
783   // D*+->D0pi syst errors. Responsible: A. Grelli
784   //  2010 PbPb sample, 0-20 CC
785   //
786
787   AliInfo(" Settings for D*+ --> D0pi Pb-Pb collisions at 2.76 TeV - 0-20 centrality - DUMMY"); 
788
789  // Normalization
790   fNorm = new TH1F("fNorm","fNorm",24,0,24);
791   for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.05); // TAA and pp norm
792
793   // Branching ratio 
794   fBR = new TH1F("fBR","fBR",24,0,24);
795   for(Int_t i=1;i<=24;i++) fBR->SetBinContent(i,0.015); // 1.5% PDG2010
796
797   // Tracking efficiency
798   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",24,0,24);;
799   for(Int_t i=1;i<=24;i++) fTrackingEff->SetBinContent(i,0.15); // Jacek, 5% per track
800
801
802   // Raw yield extraction
803   fRawYield = new TH1F("fRawYield","fRawYield",24,0,24);
804   for(Int_t i=1;i<=24;i++) fRawYield->SetBinContent(i,0.1);  //4%
805   fRawYield->SetBinContent(3,0.2);
806   fRawYield->SetBinContent(4,0.2);
807   fRawYield->SetBinContent(5,0.2);
808   fRawYield->SetBinContent(6,0.2);
809  
810   // Cuts efficiency (from cuts variation)
811   fCutsEff = new TH1F("fCutsEff","fCutsEff",24,0,24);
812   for(Int_t i=1;i<=24;i++) fCutsEff->SetBinContent(i,0.10); // 10%
813   fCutsEff->SetBinContent(4,0.15);
814   fCutsEff->SetBinContent(5,0.15);
815   fCutsEff->SetBinContent(6,0.15);
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.05); // 3%
820  
821
822   // MC dN/dpt  (from study on D* pt shape)
823   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",24,0,24);
824   for(Int_t i=1;i<=24;i++) fMCPtShape->SetBinContent(i,0.045);
825   fMCPtShape->SetBinContent(4,0.025);
826   fMCPtShape->SetBinContent(5,0.025);
827   fMCPtShape->SetBinContent(6,0.025);
828   fMCPtShape->SetBinContent(7,0.04);
829   fMCPtShape->SetBinContent(8,0.04);
830   fMCPtShape->SetBinContent(9,0.03);
831   fMCPtShape->SetBinContent(10,0.03);
832   fMCPtShape->SetBinContent(11,0.03);
833   fMCPtShape->SetBinContent(12,0.03);
834   
835   
836
837   return;
838
839 }
840
841 //-------------------------------------------------------------------------
842 void AliHFSystErr::InitDstartoD0pi2010PbPb2040() {
843   // 
844   // D*+->D0pi syst errors. Responsible: A. Grelli
845   //  2010 PbPb sample, 20-40 CC
846   //
847
848   AliInfo(" Settings for D*+ --> D0pi Pb-Pb collisions at 2.76 TeV - 20-40 centrality - DUMMY"); 
849
850  // Normalization
851   fNorm = new TH1F("fNorm","fNorm",24,0,24);
852   for(Int_t i=1;i<=24;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
853
854   // Branching ratio 
855   fBR = new TH1F("fBR","fBR",24,0,24);
856   for(Int_t i=1;i<=24;i++) fBR->SetBinContent(i,0.015); // 1.5% PDG2010
857
858   // Tracking efficiency
859   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",24,0,24);;
860   for(Int_t i=1;i<=24;i++) fTrackingEff->SetBinContent(i,0.15); // Jacek, 5% per track
861
862
863   // Raw yield extraction
864   fRawYield = new TH1F("fRawYield","fRawYield",24,0,24);
865   for(Int_t i=1;i<=24;i++) fRawYield->SetBinContent(i,0.15);  //4%
866  
867   // Cuts efficiency (from cuts variation)
868   fCutsEff = new TH1F("fCutsEff","fCutsEff",24,0,24);
869   for(Int_t i=1;i<=24;i++) fCutsEff->SetBinContent(i,0.10); // 10%
870
871   // PID efficiency (from PID/noPID)
872   fPIDEff = new TH1F("fPIDEff","fPIDEff",24,0,24);
873   for(Int_t i=1;i<=24;i++) fPIDEff->SetBinContent(i,0.04); // 3%
874  
875
876   // MC dN/dpt  (copied from D0 : will update later)
877   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",24,0,24);
878   for(Int_t i=1;i<=24;i++) fMCPtShape->SetBinContent(i,0.);
879   fMCPtShape->SetBinContent(1,0.03);
880   fMCPtShape->SetBinContent(2,0.03);
881
882   return;
883
884 }
885
886 //--------------------------------------------------------------------------
887 void AliHFSystErr::InitDstartoD0pi2010PbPb4080() {
888   // 
889   // D*+->D0pi syst errors. Responsible: A. Grelli
890   //  2010 PbPb sample, 40-80 CC
891   //
892
893   AliInfo(" Settings for D*+ --> D0pi Pb-Pb collisions at 2.76 TeV - 40-80 centrality - DUMMY"); 
894
895  // Normalization
896   fNorm = new TH1F("fNorm","fNorm",24,0,24);
897   for(Int_t i=1;i<=24;i++) fNorm->SetBinContent(i,0.07); // TAA and pp norm
898
899   // Branching ratio 
900   fBR = new TH1F("fBR","fBR",24,0,24);
901   for(Int_t i=1;i<=24;i++) fBR->SetBinContent(i,0.015); // 1.5% PDG2010
902
903   // Tracking efficiency
904   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",24,0,24);;
905   for(Int_t i=1;i<=24;i++) fTrackingEff->SetBinContent(i,0.15); // Jacek, 5% per track
906
907
908   // Raw yield extraction
909   fRawYield = new TH1F("fRawYield","fRawYield",24,0,24);
910   for(Int_t i=1;i<=24;i++) fRawYield->SetBinContent(i,0.2);  //4%
911   fRawYield->SetBinContent(1,1);
912   fRawYield->SetBinContent(2,0.15);
913   fRawYield->SetBinContent(3,0.15);
914   fRawYield->SetBinContent(4,0.15);
915   fRawYield->SetBinContent(5,0.15);
916   fRawYield->SetBinContent(6,0.10);
917   fRawYield->SetBinContent(7,0.10);
918   fRawYield->SetBinContent(8,0.10);
919   fRawYield->SetBinContent(9,0.11);
920   fRawYield->SetBinContent(10,0.11);
921   fRawYield->SetBinContent(11,0.11);
922   fRawYield->SetBinContent(12,0.11);
923   fRawYield->SetBinContent(13,0.08);
924   fRawYield->SetBinContent(14,0.08);
925   fRawYield->SetBinContent(15,0.08);
926   fRawYield->SetBinContent(16,0.08);
927
928
929   // Cuts efficiency (from cuts variation)
930   fCutsEff = new TH1F("fCutsEff","fCutsEff",24,0,24);
931   for(Int_t i=1;i<=24;i++) fCutsEff->SetBinContent(i,0.10); // 10%
932
933   // PID efficiency (from PID/noPID)
934   fPIDEff = new TH1F("fPIDEff","fPIDEff",24,0,24);
935   for(Int_t i=1;i<=24;i++) fPIDEff->SetBinContent(i,0.05); // 3%
936  
937
938   // MC dN/dpt  (copied from D0 : will update later)
939   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",24,0,24);
940   for(Int_t i=1;i<=24;i++) fMCPtShape->SetBinContent(i,0.01);
941   fMCPtShape->SetBinContent(2,0.05);
942   fMCPtShape->SetBinContent(3,0.05);
943   fMCPtShape->SetBinContent(4,0.05);
944   fMCPtShape->SetBinContent(5,0.04);
945   fMCPtShape->SetBinContent(6,0.02);
946   fMCPtShape->SetBinContent(7,0.04);
947   fMCPtShape->SetBinContent(8,0.04);
948  
949   return;
950
951 }
952
953 //--------------------------------------------------------------------------
954 void AliHFSystErr::InitD0toKpi2010PbPb010CentScan(){
955   // define errors for RAA vs. centrality
956   InitD0toKpi2010PbPb020();
957   for(Int_t i=7;i<=12;i++) fRawYield->SetBinContent(i,0.05);
958   for(Int_t i=7;i<=12;i++) fMCPtShape->SetBinContent(i,0.08); 
959 }
960 //--------------------------------------------------------------------------
961 void AliHFSystErr::InitD0toKpi2010PbPb1020CentScan(){
962   // define errors for RAA vs. centrality
963   InitD0toKpi2010PbPb020();
964   for(Int_t i=7;i<=12;i++) fRawYield->SetBinContent(i,0.05);
965   for(Int_t i=7;i<=12;i++) fMCPtShape->SetBinContent(i,0.08); 
966 }
967 //--------------------------------------------------------------------------
968 void AliHFSystErr::InitD0toKpi2010PbPb2040CentScan(){
969   // define errors for RAA vs. centrality
970   InitD0toKpi2010PbPb4080();
971   for(Int_t i=7;i<=12;i++) fRawYield->SetBinContent(i,0.05);
972   for(Int_t i=7;i<=12;i++) fMCPtShape->SetBinContent(i,0.08); 
973 }
974 //--------------------------------------------------------------------------
975 void AliHFSystErr::InitD0toKpi2010PbPb4060CentScan(){
976    // define errors for RAA vs. centrality
977   InitD0toKpi2010PbPb4080();
978   for(Int_t i=7;i<=12;i++) fRawYield->SetBinContent(i,0.06);
979   for(Int_t i=7;i<=12;i++) fMCPtShape->SetBinContent(i,0.08); 
980 }
981 //--------------------------------------------------------------------------
982 void AliHFSystErr::InitD0toKpi2010PbPb6080CentScan(){
983    // define errors for RAA vs. centrality
984   InitD0toKpi2010PbPb4080();
985   for(Int_t i=7;i<=12;i++) fRawYield->SetBinContent(i,0.08);
986   for(Int_t i=7;i<=12;i++) fMCPtShape->SetBinContent(i,0.08); 
987 }
988 //--------------------------------------------------------------------------
989 void AliHFSystErr::InitDplustoKpipi2010PbPb010CentScan(){
990   // define errors for RAA vs. centrality
991   InitDplustoKpipi2010PbPb020();
992   for(Int_t i=7;i<=12;i++) fRawYield->SetBinContent(i,0.18);
993   for(Int_t i=7;i<=12;i++) fMCPtShape->SetBinContent(i,0.09);
994
995 }
996 //--------------------------------------------------------------------------
997 void AliHFSystErr::InitDplustoKpipi2010PbPb1020CentScan(){
998   // define errors for RAA vs. centrality
999   InitDplustoKpipi2010PbPb020();
1000   for(Int_t i=7;i<=12;i++) fRawYield->SetBinContent(i,0.23);
1001   for(Int_t i=7;i<=12;i++) fMCPtShape->SetBinContent(i,0.08);
1002 }
1003 //--------------------------------------------------------------------------
1004 void AliHFSystErr::InitDplustoKpipi2010PbPb2040CentScan(){
1005   // define errors for RAA vs. centrality
1006   InitDplustoKpipi2010PbPb020();
1007   for(Int_t i=7;i<=12;i++) fRawYield->SetBinContent(i,0.08);
1008   for(Int_t i=7;i<=12;i++) fMCPtShape->SetBinContent(i,0.095);
1009 }
1010 //--------------------------------------------------------------------------
1011 void AliHFSystErr::InitDplustoKpipi2010PbPb4060CentScan(){
1012   // define errors for RAA vs. centrality
1013   InitDplustoKpipi2010PbPb4080();
1014   for(Int_t i=7;i<=12;i++) fRawYield->SetBinContent(i,0.08);
1015   for(Int_t i=7;i<=12;i++) fMCPtShape->SetBinContent(i,0.08);
1016 }
1017 //--------------------------------------------------------------------------
1018 void AliHFSystErr::InitDplustoKpipi2010PbPb6080CentScan(){
1019   // define errors for RAA vs. centrality
1020   InitDplustoKpipi2010PbPb4080();
1021   for(Int_t i=7;i<=12;i++) fRawYield->SetBinContent(i,0.15);
1022   for(Int_t i=7;i<=12;i++) fMCPtShape->SetBinContent(i,0.07);
1023 }
1024
1025 //--------------------------------------------------------------------------
1026 void AliHFSystErr::InitDstartoD0pi2010PbPb010CentScan(){
1027   // define errors for RAA vs. centrality
1028   InitDstartoD0pi2010PbPb020();
1029   for(Int_t i=7;i<=12;i++) fRawYield->SetBinContent(i,0.16); 
1030   for(Int_t i=7;i<=12;i++) fMCPtShape->SetBinContent(i,0.15);
1031 }
1032 //--------------------------------------------------------------------------
1033 void AliHFSystErr::InitDstartoD0pi2010PbPb1020CentScan(){
1034   // define errors for RAA vs. centrality
1035   InitDstartoD0pi2010PbPb020();
1036   for(Int_t i=7;i<=12;i++) fRawYield->SetBinContent(i,0.05); 
1037   for(Int_t i=7;i<=12;i++) fMCPtShape->SetBinContent(i,0.15);
1038 }
1039 //--------------------------------------------------------------------------
1040 void AliHFSystErr::InitDstartoD0pi2010PbPb2040CentScan(){
1041   // define errors for RAA vs. centrality
1042   InitDstartoD0pi2010PbPb2040();
1043   for(Int_t i=7;i<=12;i++) fRawYield->SetBinContent(i,0.10); 
1044   for(Int_t i=7;i<=12;i++) fMCPtShape->SetBinContent(i,0.08);
1045 }
1046 //--------------------------------------------------------------------------
1047 void AliHFSystErr::InitDstartoD0pi2010PbPb4060CentScan(){
1048   // define errors for RAA vs. centrality
1049   InitDstartoD0pi2010PbPb4080();
1050   for(Int_t i=7;i<=12;i++) fRawYield->SetBinContent(i,0.10); 
1051   for(Int_t i=7;i<=12;i++) fMCPtShape->SetBinContent(i,0.045);
1052 }
1053 //--------------------------------------------------------------------------
1054 void AliHFSystErr::InitDstartoD0pi2010PbPb6080CentScan(){
1055   // define errors for RAA vs. centrality
1056   InitDstartoD0pi2010PbPb4080();
1057   for(Int_t i=7;i<=12;i++) fRawYield->SetBinContent(i,0.10); 
1058   for(Int_t i=7;i<=12;i++) fMCPtShape->SetBinContent(i,0.045);
1059 }
1060
1061
1062 //--------------------------------------------------------------------------
1063 Double_t AliHFSystErr::GetCutsEffErr(Double_t pt) const {
1064   // 
1065   // Get error
1066   //
1067
1068   Int_t bin=fCutsEff->FindBin(pt);
1069
1070   return fCutsEff->GetBinContent(bin);
1071 }
1072 //--------------------------------------------------------------------------
1073 Double_t AliHFSystErr::GetMCPtShapeErr(Double_t pt) const {
1074   // 
1075   // Get error
1076   //
1077
1078   Int_t bin=fMCPtShape->FindBin(pt);
1079
1080   return fMCPtShape->GetBinContent(bin);
1081 }
1082 //--------------------------------------------------------------------------
1083 Double_t AliHFSystErr::GetSeleEffErr(Double_t pt) const {
1084   // 
1085   // Get error
1086   //
1087
1088   Double_t err=GetCutsEffErr(pt)*GetCutsEffErr(pt)+GetMCPtShapeErr(pt)*GetMCPtShapeErr(pt);
1089
1090   return TMath::Sqrt(err);
1091 }
1092 //--------------------------------------------------------------------------
1093 Double_t AliHFSystErr::GetPIDEffErr(Double_t pt) const {
1094   // 
1095   // Get error
1096   //
1097
1098   Int_t bin=fPIDEff->FindBin(pt);
1099
1100   return fPIDEff->GetBinContent(bin);
1101 }
1102 //--------------------------------------------------------------------------
1103 Double_t AliHFSystErr::GetTrackingEffErr(Double_t pt) const {
1104   // 
1105   // Get error
1106   //
1107
1108   Int_t bin=fTrackingEff->FindBin(pt);
1109
1110   return fTrackingEff->GetBinContent(bin);
1111 }
1112 //--------------------------------------------------------------------------
1113 Double_t AliHFSystErr::GetRawYieldErr(Double_t pt) const {
1114   // 
1115   // Get error
1116   //
1117
1118   Int_t bin=fRawYield->FindBin(pt);
1119
1120   return fRawYield->GetBinContent(bin);
1121 }
1122 //--------------------------------------------------------------------------
1123 Double_t AliHFSystErr::GetPartAntipartErr(Double_t pt) const {
1124   // 
1125   // Get error
1126   //
1127
1128   Int_t bin=fPartAntipart->FindBin(pt);
1129
1130   return fPartAntipart->GetBinContent(bin);
1131 }
1132 //--------------------------------------------------------------------------
1133 Double_t AliHFSystErr::GetTotalSystErr(Double_t pt,Double_t feeddownErr) const {
1134   // 
1135   // Get total syst error (except norm. error)
1136   //
1137
1138   Double_t err=0.;
1139
1140   if(fRawYield) err += GetRawYieldErr(pt)*GetRawYieldErr(pt);
1141   if(fTrackingEff) err += GetTrackingEffErr(pt)*GetTrackingEffErr(pt);
1142   //  if(fBR) err += GetBRErr()*GetBRErr();
1143   if(fCutsEff) err += GetCutsEffErr(pt)*GetCutsEffErr(pt);
1144   if(fPIDEff) err += GetPIDEffErr(pt)*GetPIDEffErr(pt);
1145   if(fMCPtShape) err += GetMCPtShapeErr(pt)*GetMCPtShapeErr(pt);
1146   if(fPartAntipart) err += GetPartAntipartErr(pt)*GetPartAntipartErr(pt);
1147
1148   err += feeddownErr*feeddownErr;
1149
1150   return TMath::Sqrt(err);
1151 }
1152 //---------------------------------------------------------------------------
1153 void AliHFSystErr::DrawErrors(TGraphAsymmErrors *grErrFeeddown) const {
1154   //
1155   // Draw errors
1156   //
1157   gStyle->SetOptStat(0);
1158
1159   TCanvas *cSystErr = new TCanvas("cSystErr","Systematic Errors",300,80,640,500);
1160   cSystErr->Range(0.20,-0.5,18.4,0.34);
1161   cSystErr->SetRightMargin(0.318);
1162   cSystErr->SetFillColor(0);
1163
1164   TH2F *hFrame = new TH2F("hFrame","Systematic errors; p_{t} [GeV/c]; Relative Error",20,0,20,100,-1,+1);
1165   hFrame->SetAxisRange(1.,15.9,"X");
1166   hFrame->SetAxisRange(-0.5,0.5,"Y");
1167   hFrame->Draw();
1168
1169   TLegend *leg = new TLegend(0.69,0.44,0.98,0.86,NULL,"brNDC");
1170   leg->SetTextSize(0.03601695);
1171   leg->SetFillStyle(0);
1172   leg->SetBorderSize(0);
1173   
1174   TH1F *hTotErr=new TH1F("hTotErr","",24,0,24);
1175   Int_t nbins = fNorm->GetNbinsX();
1176   TGraphAsymmErrors *gTotErr = new TGraphAsymmErrors(nbins);
1177   for(Int_t i=1;i<=24;i++) {
1178     Double_t pt = hTotErr->GetBinCenter(i);
1179     Double_t ptwidth = hTotErr->GetBinWidth(i);
1180
1181     if(grErrFeeddown) {
1182       Double_t x=0., y=0., errxl=0., errxh=0., erryl=0., erryh=0.;
1183       Double_t toterryl=0., toterryh=0.;
1184       for(Int_t j=0; j<grErrFeeddown->GetN(); j++) {
1185         grErrFeeddown->GetPoint(j,x,y);
1186         errxh = grErrFeeddown->GetErrorXhigh(j);
1187         errxl = grErrFeeddown->GetErrorXlow(j);
1188         if ( ( (x-errxl) <= pt) && ( (x+errxl) >= pt) ) {
1189           erryh = grErrFeeddown->GetErrorYhigh(j);
1190           erryl = grErrFeeddown->GetErrorYlow(j);
1191         }
1192       }
1193       if (erryl>=1e-3) toterryl = GetTotalSystErr(pt,erryl);
1194       else toterryl = GetTotalSystErr(pt);
1195       if (erryh>=1e-3) toterryh = GetTotalSystErr(pt,erryh);
1196       else toterryh = GetTotalSystErr(pt);
1197
1198       hTotErr->SetBinContent(i,toterryh);
1199       gTotErr->SetPoint(i,pt,0.);
1200       gTotErr->SetPointError(i,ptwidth/2.,ptwidth/2.,toterryl,toterryh); // i, exl, exh, eyl, eyh
1201     }
1202     else {
1203       hTotErr->SetBinContent(i,GetTotalSystErr(pt));
1204       gTotErr->SetPoint(i,pt,0.);
1205       gTotErr->SetPointError(i,ptwidth/2.,ptwidth/2.,GetTotalSystErr(pt),GetTotalSystErr(pt)); // i, exl, exh, eyl, eyh
1206     }
1207
1208   }
1209   gTotErr->SetLineColor(kBlack);
1210   gTotErr->SetFillColor(kRed);
1211   gTotErr->SetFillStyle(3002);
1212   gTotErr->Draw("2");
1213   leg->AddEntry(gTotErr,"Total (excl. norm.)","f");
1214 //   hTotErr->SetLineColor(1);
1215 //   hTotErr->SetLineWidth(3);
1216 //   hTotErr->Draw("same");
1217 //   leg->AddEntry(hTotErr,"Total (excl. norm.)","l");
1218   
1219
1220   fNorm->SetFillColor(1);
1221   fNorm->SetFillStyle(3002);
1222   //fNorm->Draw("same");
1223   //TH1F *hNormRefl = ReflectHisto(fNorm);
1224   //hNormRefl->Draw("same");
1225   leg->AddEntry(fNorm,"Normalization (10%)","");
1226
1227   if(grErrFeeddown) {
1228     grErrFeeddown->SetFillColor(kTeal-8);
1229     grErrFeeddown->SetFillStyle(3001);
1230     grErrFeeddown->Draw("2");
1231     leg->AddEntry(grErrFeeddown,"Feed-down from B","f");
1232   }
1233   if(fTrackingEff) {
1234     fTrackingEff->SetFillColor(4);
1235     fTrackingEff->SetFillStyle(3006);
1236     fTrackingEff->Draw("same");
1237     TH1F *hTrackingEffRefl = ReflectHisto(fTrackingEff);
1238     hTrackingEffRefl->Draw("same");
1239     leg->AddEntry(fTrackingEff,"Tracking efficiency","f");
1240   }
1241   if(fBR) {
1242     fBR->SetFillColor(6);
1243     fBR->SetFillStyle(3005);
1244     //fBR->SetFillStyle(3020);
1245     fBR->Draw("same");
1246     TH1F *hBRRefl = ReflectHisto(fBR);
1247     hBRRefl->Draw("same");
1248     leg->AddEntry(fBR,"Branching ratio","f");
1249   }
1250   if(fRawYield) {
1251     Int_t ci;   // for color index setting
1252     ci = TColor::GetColor("#00cc00");
1253     fRawYield->SetLineColor(ci);
1254     //    fRawYield->SetLineColor(3);
1255     fRawYield->SetLineWidth(3);
1256     fRawYield->Draw("same");
1257     TH1F *hRawYieldRefl = ReflectHisto(fRawYield);
1258     hRawYieldRefl->Draw("same");
1259     leg->AddEntry(fRawYield,"Yield extraction","l");
1260   }
1261   if(fCutsEff) {
1262     fCutsEff->SetLineColor(4);
1263     fCutsEff->SetLineWidth(3);
1264     fCutsEff->Draw("same");
1265     TH1F *hCutsEffRefl = ReflectHisto(fCutsEff);
1266     hCutsEffRefl->Draw("same");
1267     leg->AddEntry(fCutsEff,"Cuts efficiency","l");
1268   }
1269   if(fPIDEff) {
1270     fPIDEff->SetLineColor(7);
1271     fPIDEff->SetLineWidth(3);
1272     fPIDEff->Draw("same");
1273     TH1F *hPIDEffRefl = ReflectHisto(fPIDEff);
1274     hPIDEffRefl->Draw("same");
1275     leg->AddEntry(fPIDEff,"PID efficiency","l");
1276   }
1277   if(fMCPtShape) {
1278     Int_t ci = TColor::GetColor("#9933ff");
1279     fMCPtShape->SetLineColor(ci);
1280     //    fMCPtShape->SetLineColor(8);
1281     fMCPtShape->SetLineWidth(3);
1282     fMCPtShape->Draw("same");
1283     TH1F *hMCPtShapeRefl = ReflectHisto(fMCPtShape);
1284     hMCPtShapeRefl->Draw("same");
1285     leg->AddEntry(fMCPtShape,"MC p_{t} shape","l");
1286   }
1287   if(fPartAntipart) {
1288     Int_t ci = TColor::GetColor("#ff6600");
1289     fPartAntipart->SetLineColor(ci);
1290     //    fPartAntipart->SetLineColor(9);
1291     fPartAntipart->SetLineWidth(3);
1292     fPartAntipart->Draw("same");
1293     TH1F *hPartAntipartRefl = ReflectHisto(fPartAntipart);
1294     hPartAntipartRefl->Draw("same");
1295     leg->AddEntry(fPartAntipart,"D = #bar{D}","l");
1296   }
1297
1298
1299   leg->Draw();
1300
1301   cSystErr->SaveAs("RelativeSystematics.eps");
1302
1303   return;
1304 }
1305 //-------------------------------------------------------------------------
1306 TH1F* AliHFSystErr::ReflectHisto(TH1F *hin) const {
1307   //
1308   // Clones and reflects histogram 
1309   // 
1310   TH1F *hout=(TH1F*)hin->Clone("hout");
1311   hout->Scale(-1.);
1312
1313   return hout;
1314 }