[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliHFSystErr.cxx

/**************************************************************************
 * Copyright(c) 1998-2010, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/////////////////////////////////////////////////////////////
//
// Class to handle systematic errors for charm hadrons
//
// Usage:
// AliHFSystEff syst(DECAY); // DECAY = 1 for D0, 2, for D+, 3 for D* 
// syst.DrawErrors(); // to see a plot of the error contributions
// syst.GetTotalSystErr(pt); // to get the total err at pt 
//
// Author: A.Dainese, andrea.dainese@pd.infn.it
/////////////////////////////////////////////////////////////

#include <TStyle.h>
#include <TGraphAsymmErrors.h>
#include <TMath.h>
#include <TCanvas.h>
#include <TH2F.h>
#include <TLegend.h>

#include "AliHFSystErr.h"


ClassImp(AliHFSystErr)

//--------------------------------------------------------------------------
AliHFSystErr::AliHFSystErr(const Char_t* name, const Char_t* title) : 
TNamed(name,title),
fNorm(0),
fRawYield(0),
fTrackingEff(0),
fBR(0),
fCutsEff(0),
fPIDEff(0),
fMCPtShape(0),
fPartAntipart(0)
{
  //
  // Default Constructor
  //
}
//--------------------------------------------------------------------------
AliHFSystErr::AliHFSystErr(Int_t decay,const Char_t* name, const Char_t* title) : 
TNamed(name,title),
fNorm(0),
fRawYield(0),
fTrackingEff(0),
fBR(0),
fCutsEff(0),
fPIDEff(0),
fMCPtShape(0),
fPartAntipart(0)
{
  //
  // Standard Constructor
  //

  switch(decay) {
  case 1: // D0->Kpi
    InitD0toKpi();
    break;
  case 2: // D+->Kpipi
    InitDplustoKpipi();
    break;
  case 3: // D*->D0pi
    InitDstartoD0pi();
    break;
  default:
    printf("Invalid decay type: %d\n",decay);
    break;
  }


}
//--------------------------------------------------------------------------
AliHFSystErr::~AliHFSystErr() {
  //  
  // Default Destructor
  //

  if(fNorm)         { delete fNorm; fNorm=0; }
  if(fRawYield)     { delete fRawYield; fRawYield=0; }
  if(fTrackingEff)  { delete fTrackingEff; fTrackingEff=0; }
  if(fBR)           { delete fBR; fBR=0; }
  if(fCutsEff)      { delete fCutsEff; fCutsEff=0; }
  if(fPIDEff)       { delete fPIDEff; fPIDEff=0; }
  if(fMCPtShape)    { delete fMCPtShape; fMCPtShape=0; }
  if(fPartAntipart) { delete fPartAntipart; fPartAntipart=0; }

}
//--------------------------------------------------------------------------
void AliHFSystErr::InitD0toKpi() {
  // 
  // D0->Kpi syst errors. Responsible: A. Rossi
  //
  
  // Normalization
  fNorm = new TH1F("fNorm","fNorm",20,0,20);
  for(Int_t i=4;i<=20;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and

  // Branching ratio 
  fBR = new TH1F("fBR","fBR",20,0,20);
  for(Int_t i=4;i<=20;i++) fBR->SetBinContent(i,0.012); // 1.2% PDG2010

  // Tracking efficiency
  fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
  for(Int_t i=4;i<=20;i++) fTrackingEff->SetBinContent(i,0.02); // 2% (1% per track)

  // Raw yield extraction
  fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
  fRawYield->SetBinContent(1,1);
  fRawYield->SetBinContent(2,1);
  fRawYield->SetBinContent(3,0.15);
  for(Int_t i=4;i<=20;i++) fRawYield->SetBinContent(i,0.06);

  // Cuts efficiency (from cuts variation)
  fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
  for(Int_t i=1;i<=20;i++) fCutsEff->SetBinContent(i,0.10); // 10%

  // PID efficiency (from PID/noPID)
  fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
  for(Int_t i=1;i<=20;i++) fPIDEff->SetBinContent(i,0.03); // 3%
  fPIDEff->SetBinContent(4,0.15); // 15%

  // MC dN/dpt
  fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
  for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,(Float_t)i*0.006);

  // particle-antiparticle
  fPartAntipart = new TH1F("fPartAntipart","fPartAntipart",20,0,20);
  fPartAntipart->SetBinContent(1,1);
  fPartAntipart->SetBinContent(2,1);
  fPartAntipart->SetBinContent(3,0.20);
  for(Int_t i=4;i<=20;i++) fPartAntipart->SetBinContent(i,0.05);
  
  return;
}
//--------------------------------------------------------------------------
void AliHFSystErr::InitDplustoKpipi() {
  // 
  // D+->Kpipi syst errors. Responsible: R. Bala
  //

  return;
}
//--------------------------------------------------------------------------
void AliHFSystErr::InitDstartoD0pi() {
  // 
  // D*+->D0pi syst errors. Responsible: tbd
  //

  return;
}
//--------------------------------------------------------------------------
Double_t AliHFSystErr::GetCutsEffErr(Double_t pt) const {
  // 
  // Get error
  //

  Int_t bin=fCutsEff->FindBin(pt);

  return fCutsEff->GetBinContent(bin);
}
//--------------------------------------------------------------------------
Double_t AliHFSystErr::GetMCPtShapeErr(Double_t pt) const {
  // 
  // Get error
  //

  Int_t bin=fMCPtShape->FindBin(pt);

  return fMCPtShape->GetBinContent(bin);
}
//--------------------------------------------------------------------------
Double_t AliHFSystErr::GetSeleEffErr(Double_t pt) const {
  // 
  // Get error
  //

  Double_t err=GetCutsEffErr(pt)*GetCutsEffErr(pt)+GetMCPtShapeErr(pt)*GetMCPtShapeErr(pt);

  return TMath::Sqrt(err);
}
//--------------------------------------------------------------------------
Double_t AliHFSystErr::GetPIDEffErr(Double_t pt) const {
  // 
  // Get error
  //

  Int_t bin=fPIDEff->FindBin(pt);

  return fPIDEff->GetBinContent(bin);
}
//--------------------------------------------------------------------------
Double_t AliHFSystErr::GetTrackingEffErr(Double_t pt) const {
  // 
  // Get error
  //

  Int_t bin=fTrackingEff->FindBin(pt);

  return fTrackingEff->GetBinContent(bin);
}
//--------------------------------------------------------------------------
Double_t AliHFSystErr::GetRawYieldErr(Double_t pt) const {
  // 
  // Get error
  //

  Int_t bin=fRawYield->FindBin(pt);

  return fRawYield->GetBinContent(bin);
}
//--------------------------------------------------------------------------
Double_t AliHFSystErr::GetPartAntipartErr(Double_t pt) const {
  // 
  // Get error
  //

  Int_t bin=fPartAntipart->FindBin(pt);

  return fPartAntipart->GetBinContent(bin);
}
//--------------------------------------------------------------------------
Double_t AliHFSystErr::GetTotalSystErr(Double_t pt,Double_t feeddownErr) const {
  // 
  // Get total syst error (except norm. error)
  //

  Double_t err=0.;

  if(fRawYield) err += GetRawYieldErr(pt)*GetRawYieldErr(pt);
  if(fTrackingEff) err += GetTrackingEffErr(pt)*GetTrackingEffErr(pt);
  if(fBR) err += GetBRErr()*GetBRErr();
  if(fCutsEff) err += GetCutsEffErr(pt)*GetCutsEffErr(pt);
  if(fPIDEff) err += GetPIDEffErr(pt)*GetPIDEffErr(pt);
  if(fMCPtShape) err += GetMCPtShapeErr(pt)*GetMCPtShapeErr(pt);
  if(fPartAntipart) err += GetPartAntipartErr(pt)*GetPartAntipartErr(pt);

  err += feeddownErr*feeddownErr;

  return TMath::Sqrt(err);
}
//---------------------------------------------------------------------------
void AliHFSystErr::DrawErrors(TGraphErrors *grErrFeeddown) const {
  //
  // Draw errors
  //
  gStyle->SetOptStat(0);

  TCanvas *cSystErr = new TCanvas("cSystErr","Systematic Errors",0,0,500,500);
  cSystErr->SetFillColor(0);

  TH2F *hFrame = new TH2F("hFrame","Systematic errors; p_{t} [GeV/c]; Relative Error",20,0,20,100,0,+1);
  hFrame->Draw();

  TLegend *leg=new TLegend(0.5,0.5,0.9,0.9);
  leg->SetFillStyle(0);
  leg->SetBorderSize(0);

  if(fNorm) {
    fNorm->SetFillColor(1);
    fNorm->SetFillStyle(3001);
    fNorm->Draw("same");
    leg->AddEntry(fNorm,"Normalization","f");
  }
  if(fTrackingEff) {
    fTrackingEff->SetFillColor(2);
    fTrackingEff->SetFillStyle(3004);
    fTrackingEff->Draw("same");
    leg->AddEntry(fTrackingEff,"Tracking efficiency","f");
  }
  if(fBR) {
    fBR->SetFillColor(6);
    fBR->SetFillStyle(3005);
    fBR->Draw("same");
    leg->AddEntry(fBR,"Branching ratio","f");
  }
  if(fRawYield) {
    fRawYield->SetLineColor(3);
    fRawYield->SetLineWidth(3);
    fRawYield->Draw("same");
    leg->AddEntry(fRawYield,"Inv. mass analysis","l");
  }
  if(fCutsEff) {
    fCutsEff->SetLineColor(4);
    fCutsEff->SetLineWidth(3);
    fCutsEff->Draw("same");
    leg->AddEntry(fCutsEff,"Cuts efficiency","l");
  }
  if(fPIDEff) {
    fPIDEff->SetLineColor(7);
    fPIDEff->SetLineWidth(3);
    fPIDEff->Draw("same");
    leg->AddEntry(fPIDEff,"PID efficiency","l");
  }
  if(fMCPtShape) {
    fMCPtShape->SetLineColor(8);
    fMCPtShape->SetLineWidth(3);
    fMCPtShape->Draw("same");
    leg->AddEntry(fMCPtShape,"MC p_{t} shape","l");
  }
  if(fPartAntipart) {
    fPartAntipart->SetLineColor(9);
    fPartAntipart->SetLineWidth(3);
    fPartAntipart->Draw("same");
    leg->AddEntry(fPartAntipart,"D = #bar{D}","l");
  }
  if(grErrFeeddown) {
    grErrFeeddown->SetFillColor(11);    
    grErrFeeddown->SetFillStyle(3001);
    grErrFeeddown->Draw("f");
    leg->AddEntry(grErrFeeddown,"Feed-down from B","f");
  }

  
  TH1F *hTotErr=new TH1F("hTotErr","",20,0,20);
  for(Int_t i=1;i<=20;i++) {
    Double_t pt=hTotErr->GetBinCenter(i);
    hTotErr->SetBinContent(i,GetTotalSystErr(pt));
  }
  hTotErr->SetLineColor(1);
  hTotErr->SetLineWidth(3);
  hTotErr->Draw("same");
  leg->AddEntry(hTotErr,"Total (excl. norm.)","l");
  

  leg->Draw();

  return;
}
Commit	Line	Data
0233abe6	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	/////////////////////////////////////////////////////////////
	17	//
	18	// Class to handle systematic errors for charm hadrons
	19	//
	20	// Usage:
	21	// AliHFSystEff syst(DECAY); // DECAY = 1 for D0, 2, for D+, 3 for D*
	22	// syst.DrawErrors(); // to see a plot of the error contributions
	23	// syst.GetTotalSystErr(pt); // to get the total err at pt
	24	//
	25	// Author: A.Dainese, andrea.dainese@pd.infn.it
	26	/////////////////////////////////////////////////////////////
	27
	28	#include <TStyle.h>
	29	#include <TGraphAsymmErrors.h>
	30	#include <TMath.h>
	31	#include <TCanvas.h>
	32	#include <TH2F.h>
	33	#include <TLegend.h>
	34
	35	#include "AliHFSystErr.h"
	36
	37
	38	ClassImp(AliHFSystErr)
	39
	40	//--------------------------------------------------------------------------
	41	AliHFSystErr::AliHFSystErr(const Char_t* name, const Char_t* title) :
	42	TNamed(name,title),
	43	fNorm(0),
	44	fRawYield(0),
	45	fTrackingEff(0),
	46	fBR(0),
	47	fCutsEff(0),
	48	fPIDEff(0),
	49	fMCPtShape(0),
	50	fPartAntipart(0)
	51	{
	52	//
	53	// Default Constructor
	54	//
	55	}
	56	//--------------------------------------------------------------------------
	57	AliHFSystErr::AliHFSystErr(Int_t decay,const Char_t* name, const Char_t* title) :
	58	TNamed(name,title),
	59	fNorm(0),
	60	fRawYield(0),
	61	fTrackingEff(0),
	62	fBR(0),
	63	fCutsEff(0),
	64	fPIDEff(0),
65	fMCPtShape(0),
66	fPartAntipart(0)
67	{
68	//
69	// Standard Constructor
70	//
71
72	switch(decay) {
73	case 1: // D0->Kpi
74	InitD0toKpi();
75	break;
76	case 2: // D+->Kpipi
77	InitDplustoKpipi();
78	break;
79	case 3: // D*->D0pi
80	InitDstartoD0pi();
81	break;
82	default:
83	printf("Invalid decay type: %d\n",decay);
84	break;
85	}
86
87
88	}
89	//--------------------------------------------------------------------------
90	AliHFSystErr::~AliHFSystErr() {
91	//
92	// Default Destructor
93	//
94
95	if(fNorm) { delete fNorm; fNorm=0; }
96	if(fRawYield) { delete fRawYield; fRawYield=0; }
97	if(fTrackingEff) { delete fTrackingEff; fTrackingEff=0; }
98	if(fBR) { delete fBR; fBR=0; }
99	if(fCutsEff) { delete fCutsEff; fCutsEff=0; }
100	if(fPIDEff) { delete fPIDEff; fPIDEff=0; }
101	if(fMCPtShape) { delete fMCPtShape; fMCPtShape=0; }
102	if(fPartAntipart) { delete fPartAntipart; fPartAntipart=0; }
103
104	}
105	//--------------------------------------------------------------------------
106	void AliHFSystErr::InitD0toKpi() {
107	//
108	// D0->Kpi syst errors. Responsible: A. Rossi
109	//
110
111	// Normalization
112	fNorm = new TH1F("fNorm","fNorm",20,0,20);
113	for(Int_t i=4;i<=20;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
114
115	// Branching ratio
116	fBR = new TH1F("fBR","fBR",20,0,20);
117	for(Int_t i=4;i<=20;i++) fBR->SetBinContent(i,0.012); // 1.2% PDG2010
118
119	// Tracking efficiency
120	fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
121	for(Int_t i=4;i<=20;i++) fTrackingEff->SetBinContent(i,0.02); // 2% (1% per track)
122
123	// Raw yield extraction
124	fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
125	fRawYield->SetBinContent(1,1);
126	fRawYield->SetBinContent(2,1);
127	fRawYield->SetBinContent(3,0.15);
128	for(Int_t i=4;i<=20;i++) fRawYield->SetBinContent(i,0.06);
129
130	// Cuts efficiency (from cuts variation)
131	fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
132	for(Int_t i=1;i<=20;i++) fCutsEff->SetBinContent(i,0.10); // 10%
133
134	// PID efficiency (from PID/noPID)
135	fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
136	for(Int_t i=1;i<=20;i++) fPIDEff->SetBinContent(i,0.03); // 3%
137	fPIDEff->SetBinContent(4,0.15); // 15%
138
139	// MC dN/dpt
140	fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
141	for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,(Float_t)i*0.006);
142
143	// particle-antiparticle
144	fPartAntipart = new TH1F("fPartAntipart","fPartAntipart",20,0,20);
145	fPartAntipart->SetBinContent(1,1);
146	fPartAntipart->SetBinContent(2,1);
147	fPartAntipart->SetBinContent(3,0.20);
148	for(Int_t i=4;i<=20;i++) fPartAntipart->SetBinContent(i,0.05);
149
150	return;
151	}
152	//--------------------------------------------------------------------------
153	void AliHFSystErr::InitDplustoKpipi() {
154	//
155	// D+->Kpipi syst errors. Responsible: R. Bala
156	//
157
158	return;
159	}
160	//--------------------------------------------------------------------------
161	void AliHFSystErr::InitDstartoD0pi() {
162	//
163	// D*+->D0pi syst errors. Responsible: tbd
164	//
165
166	return;
167	}
168	//--------------------------------------------------------------------------
169	Double_t AliHFSystErr::GetCutsEffErr(Double_t pt) const {
170	//
171	// Get error
172	//
173
174	Int_t bin=fCutsEff->FindBin(pt);
175
176	return fCutsEff->GetBinContent(bin);
177	}
178	//--------------------------------------------------------------------------
179	Double_t AliHFSystErr::GetMCPtShapeErr(Double_t pt) const {
180	//
181	// Get error
182	//
183
184	Int_t bin=fMCPtShape->FindBin(pt);
185
186	return fMCPtShape->GetBinContent(bin);
187	}
188	//--------------------------------------------------------------------------
189	Double_t AliHFSystErr::GetSeleEffErr(Double_t pt) const {
190	//
191	// Get error
192	//
193
194	Double_t err=GetCutsEffErr(pt)GetCutsEffErr(pt)+GetMCPtShapeErr(pt)GetMCPtShapeErr(pt);
195
196	return TMath::Sqrt(err);
197	}
198	//--------------------------------------------------------------------------
199	Double_t AliHFSystErr::GetPIDEffErr(Double_t pt) const {
200	//
201	// Get error
202	//
203
204	Int_t bin=fPIDEff->FindBin(pt);
205
206	return fPIDEff->GetBinContent(bin);
207	}
208	//--------------------------------------------------------------------------
209	Double_t AliHFSystErr::GetTrackingEffErr(Double_t pt) const {
210	//
211	// Get error
212	//
213
214	Int_t bin=fTrackingEff->FindBin(pt);
215
216	return fTrackingEff->GetBinContent(bin);
217	}
218	//--------------------------------------------------------------------------
219	Double_t AliHFSystErr::GetRawYieldErr(Double_t pt) const {
220	//
221	// Get error
222	//
223
224	Int_t bin=fRawYield->FindBin(pt);
225
226	return fRawYield->GetBinContent(bin);
227	}
228	//--------------------------------------------------------------------------
229	Double_t AliHFSystErr::GetPartAntipartErr(Double_t pt) const {
230	//
231	// Get error
232	//
233
234	Int_t bin=fPartAntipart->FindBin(pt);
235
236	return fPartAntipart->GetBinContent(bin);
237	}
238	//--------------------------------------------------------------------------
239	Double_t AliHFSystErr::GetTotalSystErr(Double_t pt,Double_t feeddownErr) const {
240	//
241	// Get total syst error (except norm. error)
242	//
243
244	Double_t err=0.;
245
246	if(fRawYield) err += GetRawYieldErr(pt)*GetRawYieldErr(pt);
247	if(fTrackingEff) err += GetTrackingEffErr(pt)*GetTrackingEffErr(pt);
248	if(fBR) err += GetBRErr()*GetBRErr();
249	if(fCutsEff) err += GetCutsEffErr(pt)*GetCutsEffErr(pt);
250	if(fPIDEff) err += GetPIDEffErr(pt)*GetPIDEffErr(pt);
251	if(fMCPtShape) err += GetMCPtShapeErr(pt)*GetMCPtShapeErr(pt);
252	if(fPartAntipart) err += GetPartAntipartErr(pt)*GetPartAntipartErr(pt);
253
254	err += feeddownErr*feeddownErr;
255
256	return TMath::Sqrt(err);
257	}
258	//---------------------------------------------------------------------------
259	void AliHFSystErr::DrawErrors(TGraphErrors *grErrFeeddown) const {
260	//
261	// Draw errors
262	//
263	gStyle->SetOptStat(0);
264
265	TCanvas *cSystErr = new TCanvas("cSystErr","Systematic Errors",0,0,500,500);
266	cSystErr->SetFillColor(0);
267
268	TH2F *hFrame = new TH2F("hFrame","Systematic errors; p_{t} [GeV/c]; Relative Error",20,0,20,100,0,+1);
269	hFrame->Draw();
270
271	TLegend *leg=new TLegend(0.5,0.5,0.9,0.9);
272	leg->SetFillStyle(0);
273	leg->SetBorderSize(0);
274
275	if(fNorm) {
276	fNorm->SetFillColor(1);
277	fNorm->SetFillStyle(3001);
278	fNorm->Draw("same");
279	leg->AddEntry(fNorm,"Normalization","f");
280	}
281	if(fTrackingEff) {
282	fTrackingEff->SetFillColor(2);
283	fTrackingEff->SetFillStyle(3004);
284	fTrackingEff->Draw("same");
285	leg->AddEntry(fTrackingEff,"Tracking efficiency","f");
286	}
287	if(fBR) {
288	fBR->SetFillColor(6);
289	fBR->SetFillStyle(3005);
290	fBR->Draw("same");
291	leg->AddEntry(fBR,"Branching ratio","f");
292	}
293	if(fRawYield) {
294	fRawYield->SetLineColor(3);
295	fRawYield->SetLineWidth(3);
296	fRawYield->Draw("same");
297	leg->AddEntry(fRawYield,"Inv. mass analysis","l");
298	}
299	if(fCutsEff) {
300	fCutsEff->SetLineColor(4);
301	fCutsEff->SetLineWidth(3);
302	fCutsEff->Draw("same");
303	leg->AddEntry(fCutsEff,"Cuts efficiency","l");
304	}
305	if(fPIDEff) {
306	fPIDEff->SetLineColor(7);
307	fPIDEff->SetLineWidth(3);
308	fPIDEff->Draw("same");
309	leg->AddEntry(fPIDEff,"PID efficiency","l");
310	}
311	if(fMCPtShape) {
312	fMCPtShape->SetLineColor(8);
313	fMCPtShape->SetLineWidth(3);
314	fMCPtShape->Draw("same");
315	leg->AddEntry(fMCPtShape,"MC p_{t} shape","l");
316	}
317	if(fPartAntipart) {
318	fPartAntipart->SetLineColor(9);
319	fPartAntipart->SetLineWidth(3);
320	fPartAntipart->Draw("same");
321	leg->AddEntry(fPartAntipart,"D = #bar{D}","l");
322	}
323	if(grErrFeeddown) {
324	grErrFeeddown->SetFillColor(11);
325	grErrFeeddown->SetFillStyle(3001);
326	grErrFeeddown->Draw("f");
327	leg->AddEntry(grErrFeeddown,"Feed-down from B","f");
328	}
329
330
331	TH1F *hTotErr=new TH1F("hTotErr","",20,0,20);
332	for(Int_t i=1;i<=20;i++) {
333	Double_t pt=hTotErr->GetBinCenter(i);
334	hTotErr->SetBinContent(i,GetTotalSystErr(pt));
335	}
336	hTotErr->SetLineColor(1);
337	hTotErr->SetLineWidth(3);
338	hTotErr->Draw("same");
339	leg->AddEntry(hTotErr,"Total (excl. norm.)","l");
340
341
342	leg->Draw();
343
344	return;
345	}