[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 <TColor.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.065);

  // 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
  //

 // 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.04); // 4% PDG2010

  // Tracking efficiency
  fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
  for(Int_t i=4;i<=20;i++) fTrackingEff->SetBinContent(i,0.03); // 3% (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.10);
  for(Int_t i=4;i<=20;i++) fRawYield->SetBinContent(i,0.055);  //5 to 10%

  // 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.05); // 5%
  fPIDEff->SetBinContent(3,0.13); // 13%
 

  // MC dN/dpt  (copied from D0 : will update later)
  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.12);
  for(Int_t i=4;i<=20;i++) fPartAntipart->SetBinContent(i,0.05);   //5 to 12%
  */
  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(TGraphAsymmErrors *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->SetAxisRange(2.,11.9,"X");
  hFrame->SetAxisRange(0.,0.5,"Y");
  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(3002);
    fNorm->Draw("same");
    leg->AddEntry(fNorm,"Normalization","f");
  }
  if(grErrFeeddown) {
    grErrFeeddown->SetFillColor(kTeal-8);
    grErrFeeddown->SetFillStyle(3001);
    grErrFeeddown->Draw("2");
    leg->AddEntry(grErrFeeddown,"Feed-down from B","f");
  }
  if(fTrackingEff) {
    fTrackingEff->SetFillColor(4);
    fTrackingEff->SetFillStyle(3006);
    fTrackingEff->Draw("same");
    leg->AddEntry(fTrackingEff,"Tracking efficiency","f");
  }
  if(fBR) {
    fBR->SetFillColor(6);
    fBR->SetFillStyle(3005);
    //fBR->SetFillStyle(3020);
    fBR->Draw("same");
    leg->AddEntry(fBR,"Branching ratio","f");
  }
  if(fRawYield) {
    Int_t ci;   // for color index setting
    ci = TColor::GetColor("#00cc00");
    fRawYield->SetLineColor(ci);
    //    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) {
    Int_t ci = TColor::GetColor("#9933ff");
    fMCPtShape->SetLineColor(ci);
    //    fMCPtShape->SetLineColor(8);
    fMCPtShape->SetLineWidth(3);
    fMCPtShape->Draw("same");
    leg->AddEntry(fMCPtShape,"MC p_{t} shape","l");
  }
  if(fPartAntipart) {
    Int_t ci = TColor::GetColor("#ff6600");
    fPartAntipart->SetLineColor(ci);
    //    fPartAntipart->SetLineColor(9);
    fPartAntipart->SetLineWidth(3);
    fPartAntipart->Draw("same");
    leg->AddEntry(fPartAntipart,"D = #bar{D}","l");
  }

  
  TH1F *hTotErr=new TH1F("hTotErr","",20,0,20);
  for(Int_t i=1;i<=20;i++) {
    Double_t pt=hTotErr->GetBinCenter(i);
    if(grErrFeeddown) {
      Double_t x=0., y=0., erryh = 0.;
      for(Int_t j=0; j<grErrFeeddown->GetN(); j++) {
	grErrFeeddown->GetPoint(j,x,y);
	if ( TMath::Abs(x - pt) < 0.5) erryh = grErrFeeddown->GetErrorYhigh(j);   
      }
      hTotErr->SetBinContent(i,GetTotalSystErr(pt,erryh));
    } else { 
      hTotErr->SetBinContent(i,GetTotalSystErr(pt));
    }
  }
  hTotErr->SetLineColor(1);
  hTotErr->SetLineWidth(3);
  hTotErr->Draw("same");
  leg->AddEntry(hTotErr,"Total (excl. norm.)","l");
  

  leg->Draw();

  cSystErr->SaveAs("RelativeSystematics.eps");

  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>
5ad8ad40	34	#include <TColor.h>
0233abe6	35
	36	#include "AliHFSystErr.h"
	37
	38
	39	ClassImp(AliHFSystErr)
	40
	41	//--------------------------------------------------------------------------
	42	AliHFSystErr::AliHFSystErr(const Char_t* name, const Char_t* title) :
	43	TNamed(name,title),
	44	fNorm(0),
	45	fRawYield(0),
	46	fTrackingEff(0),
	47	fBR(0),
	48	fCutsEff(0),
	49	fPIDEff(0),
	50	fMCPtShape(0),
	51	fPartAntipart(0)
	52	{
	53	//
	54	// Default Constructor
	55	//
	56	}
	57	//--------------------------------------------------------------------------
	58	AliHFSystErr::AliHFSystErr(Int_t decay,const Char_t* name, const Char_t* title) :
	59	TNamed(name,title),
	60	fNorm(0),
	61	fRawYield(0),
	62	fTrackingEff(0),
	63	fBR(0),
	64	fCutsEff(0),
	65	fPIDEff(0),
	66	fMCPtShape(0),
	67	fPartAntipart(0)
	68	{
	69	//
	70	// Standard Constructor
	71	//
	72
	73	switch(decay) {
	74	case 1: // D0->Kpi
	75	InitD0toKpi();
	76	break;
	77	case 2: // D+->Kpipi
	78	InitDplustoKpipi();
	79	break;
	80	case 3: // D*->D0pi
	81	InitDstartoD0pi();
	82	break;
	83	default:
	84	printf("Invalid decay type: %d\n",decay);
	85	break;
	86	}
	87
	88
	89	}
	90	//--------------------------------------------------------------------------
	91	AliHFSystErr::~AliHFSystErr() {
	92	//
	93	// Default Destructor
	94	//
	95
	96	if(fNorm) { delete fNorm; fNorm=0; }
	97	if(fRawYield) { delete fRawYield; fRawYield=0; }
	98	if(fTrackingEff) { delete fTrackingEff; fTrackingEff=0; }
99	if(fBR) { delete fBR; fBR=0; }
100	if(fCutsEff) { delete fCutsEff; fCutsEff=0; }
101	if(fPIDEff) { delete fPIDEff; fPIDEff=0; }
102	if(fMCPtShape) { delete fMCPtShape; fMCPtShape=0; }
103	if(fPartAntipart) { delete fPartAntipart; fPartAntipart=0; }
104
105	}
106	//--------------------------------------------------------------------------
107	void AliHFSystErr::InitD0toKpi() {
108	//
109	// D0->Kpi syst errors. Responsible: A. Rossi
110	//
111
112	// Normalization
113	fNorm = new TH1F("fNorm","fNorm",20,0,20);
114	for(Int_t i=4;i<=20;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
115
116	// Branching ratio
117	fBR = new TH1F("fBR","fBR",20,0,20);
118	for(Int_t i=4;i<=20;i++) fBR->SetBinContent(i,0.012); // 1.2% PDG2010
119
120	// Tracking efficiency
121	fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
122	for(Int_t i=4;i<=20;i++) fTrackingEff->SetBinContent(i,0.02); // 2% (1% per track)
123
124	// Raw yield extraction
125	fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
126	fRawYield->SetBinContent(1,1);
127	fRawYield->SetBinContent(2,1);
128	fRawYield->SetBinContent(3,0.15);
5ad8ad40	129	for(Int_t i=4;i<=20;i++) fRawYield->SetBinContent(i,0.065);
0233abe6	130
	131	// Cuts efficiency (from cuts variation)
	132	fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
	133	for(Int_t i=1;i<=20;i++) fCutsEff->SetBinContent(i,0.10); // 10%
	134
	135	// PID efficiency (from PID/noPID)
	136	fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
	137	for(Int_t i=1;i<=20;i++) fPIDEff->SetBinContent(i,0.03); // 3%
	138	fPIDEff->SetBinContent(4,0.15); // 15%
	139
	140	// MC dN/dpt
	141	fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
	142	for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,(Float_t)i*0.006);
	143
	144	// particle-antiparticle
	145	fPartAntipart = new TH1F("fPartAntipart","fPartAntipart",20,0,20);
	146	fPartAntipart->SetBinContent(1,1);
	147	fPartAntipart->SetBinContent(2,1);
	148	fPartAntipart->SetBinContent(3,0.20);
	149	for(Int_t i=4;i<=20;i++) fPartAntipart->SetBinContent(i,0.05);
	150
	151	return;
	152	}
	153	//--------------------------------------------------------------------------
	154	void AliHFSystErr::InitDplustoKpipi() {
	155	//
	156	// D+->Kpipi syst errors. Responsible: R. Bala
	157	//
	158
626bcf57	159	// Normalization
	160	fNorm = new TH1F("fNorm","fNorm",20,0,20);
	161	for(Int_t i=4;i<=20;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
	162
	163	// Branching ratio
	164	fBR = new TH1F("fBR","fBR",20,0,20);
	165	for(Int_t i=4;i<=20;i++) fBR->SetBinContent(i,0.04); // 4% PDG2010
	166
	167	// Tracking efficiency
	168	fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
	169	for(Int_t i=4;i<=20;i++) fTrackingEff->SetBinContent(i,0.03); // 3% (1% per track)
	170
	171
	172	// Raw yield extraction
	173	fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
	174	fRawYield->SetBinContent(1,1);
	175	fRawYield->SetBinContent(2,1);
	176	fRawYield->SetBinContent(3,0.10);
5ad8ad40	177	for(Int_t i=4;i<=20;i++) fRawYield->SetBinContent(i,0.055); //5 to 10%
626bcf57	178
	179	// Cuts efficiency (from cuts variation)
	180	fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
	181	for(Int_t i=1;i<=20;i++) fCutsEff->SetBinContent(i,0.10); // 10%
	182
	183	// PID efficiency (from PID/noPID)
	184	fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
	185	for(Int_t i=1;i<=20;i++) fPIDEff->SetBinContent(i,0.05); // 5%
	186	fPIDEff->SetBinContent(3,0.13); // 13%
	187
	188
	189	// MC dN/dpt (copied from D0 : will update later)
	190	fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
	191	for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,(Float_t)i*0.006);
	192
	193
	194	// particle-antiparticle
	195	/*
	196	fPartAntipart = new TH1F("fPartAntipart","fPartAntipart",20,0,20);
	197	fPartAntipart->SetBinContent(1,1);
	198	fPartAntipart->SetBinContent(2,1);
	199	fPartAntipart->SetBinContent(3,0.12);
	200	for(Int_t i=4;i<=20;i++) fPartAntipart->SetBinContent(i,0.05); //5 to 12%
	201	*/
0233abe6	202	return;
	203	}
	204	//--------------------------------------------------------------------------
	205	void AliHFSystErr::InitDstartoD0pi() {
	206	//
	207	// D*+->D0pi syst errors. Responsible: tbd
	208	//
	209
	210	return;
	211	}
	212	//--------------------------------------------------------------------------
	213	Double_t AliHFSystErr::GetCutsEffErr(Double_t pt) const {
	214	//
	215	// Get error
	216	//
	217
	218	Int_t bin=fCutsEff->FindBin(pt);
	219
	220	return fCutsEff->GetBinContent(bin);
	221	}
	222	//--------------------------------------------------------------------------
	223	Double_t AliHFSystErr::GetMCPtShapeErr(Double_t pt) const {
	224	//
	225	// Get error
	226	//
	227
	228	Int_t bin=fMCPtShape->FindBin(pt);
	229
	230	return fMCPtShape->GetBinContent(bin);
	231	}
	232	//--------------------------------------------------------------------------
	233	Double_t AliHFSystErr::GetSeleEffErr(Double_t pt) const {
	234	//
	235	// Get error
	236	//
	237
	238	Double_t err=GetCutsEffErr(pt)GetCutsEffErr(pt)+GetMCPtShapeErr(pt)GetMCPtShapeErr(pt);
	239
	240	return TMath::Sqrt(err);
	241	}
	242	//--------------------------------------------------------------------------
	243	Double_t AliHFSystErr::GetPIDEffErr(Double_t pt) const {
	244	//
	245	// Get error
	246	//
	247
	248	Int_t bin=fPIDEff->FindBin(pt);
	249
	250	return fPIDEff->GetBinContent(bin);
	251	}
	252	//--------------------------------------------------------------------------
	253	Double_t AliHFSystErr::GetTrackingEffErr(Double_t pt) const {
	254	//
	255	// Get error
	256	//
	257
	258	Int_t bin=fTrackingEff->FindBin(pt);
	259
	260	return fTrackingEff->GetBinContent(bin);
	261	}
	262	//--------------------------------------------------------------------------
	263	Double_t AliHFSystErr::GetRawYieldErr(Double_t pt) const {
	264	//
	265	// Get error
266	//
267
268	Int_t bin=fRawYield->FindBin(pt);
269
270	return fRawYield->GetBinContent(bin);
271	}
272	//--------------------------------------------------------------------------
273	Double_t AliHFSystErr::GetPartAntipartErr(Double_t pt) const {
274	//
275	// Get error
276	//
277
278	Int_t bin=fPartAntipart->FindBin(pt);
279
280	return fPartAntipart->GetBinContent(bin);
281	}
282	//--------------------------------------------------------------------------
283	Double_t AliHFSystErr::GetTotalSystErr(Double_t pt,Double_t feeddownErr) const {
284	//
285	// Get total syst error (except norm. error)
286	//
287
288	Double_t err=0.;
289
290	if(fRawYield) err += GetRawYieldErr(pt)*GetRawYieldErr(pt);
291	if(fTrackingEff) err += GetTrackingEffErr(pt)*GetTrackingEffErr(pt);
292	if(fBR) err += GetBRErr()*GetBRErr();
293	if(fCutsEff) err += GetCutsEffErr(pt)*GetCutsEffErr(pt);
294	if(fPIDEff) err += GetPIDEffErr(pt)*GetPIDEffErr(pt);
295	if(fMCPtShape) err += GetMCPtShapeErr(pt)*GetMCPtShapeErr(pt);
296	if(fPartAntipart) err += GetPartAntipartErr(pt)*GetPartAntipartErr(pt);
297
298	err += feeddownErr*feeddownErr;
299
300	return TMath::Sqrt(err);
301	}
302	//---------------------------------------------------------------------------
5ad8ad40	303	void AliHFSystErr::DrawErrors(TGraphAsymmErrors *grErrFeeddown) const {
0233abe6	304	//
	305	// Draw errors
	306	//
	307	gStyle->SetOptStat(0);
	308
	309	TCanvas *cSystErr = new TCanvas("cSystErr","Systematic Errors",0,0,500,500);
	310	cSystErr->SetFillColor(0);
	311
	312	TH2F *hFrame = new TH2F("hFrame","Systematic errors; p_{t} [GeV/c]; Relative Error",20,0,20,100,0,+1);
5ad8ad40	313	hFrame->SetAxisRange(2.,11.9,"X");
5ad8ad40	314	hFrame->SetAxisRange(0.,0.5,"Y");
0233abe6	315	hFrame->Draw();
	316
	317	TLegend *leg=new TLegend(0.5,0.5,0.9,0.9);
	318	leg->SetFillStyle(0);
	319	leg->SetBorderSize(0);
	320
	321	if(fNorm) {
	322	fNorm->SetFillColor(1);
5ad8ad40	323	fNorm->SetFillStyle(3002);
0233abe6	324	fNorm->Draw("same");
	325	leg->AddEntry(fNorm,"Normalization","f");
	326	}
5ad8ad40	327	if(grErrFeeddown) {
	328	grErrFeeddown->SetFillColor(kTeal-8);
	329	grErrFeeddown->SetFillStyle(3001);
	330	grErrFeeddown->Draw("2");
	331	leg->AddEntry(grErrFeeddown,"Feed-down from B","f");
	332	}
0233abe6	333	if(fTrackingEff) {
5ad8ad40	334	fTrackingEff->SetFillColor(4);
5ad8ad40	335	fTrackingEff->SetFillStyle(3006);
0233abe6	336	fTrackingEff->Draw("same");
	337	leg->AddEntry(fTrackingEff,"Tracking efficiency","f");
	338	}
	339	if(fBR) {
	340	fBR->SetFillColor(6);
	341	fBR->SetFillStyle(3005);
5ad8ad40	342	//fBR->SetFillStyle(3020);
0233abe6	343	fBR->Draw("same");
	344	leg->AddEntry(fBR,"Branching ratio","f");
	345	}
	346	if(fRawYield) {
5ad8ad40	347	Int_t ci; // for color index setting
	348	ci = TColor::GetColor("#00cc00");
	349	fRawYield->SetLineColor(ci);
	350	// fRawYield->SetLineColor(3);
0233abe6	351	fRawYield->SetLineWidth(3);
	352	fRawYield->Draw("same");
	353	leg->AddEntry(fRawYield,"Inv. mass analysis","l");
	354	}
	355	if(fCutsEff) {
	356	fCutsEff->SetLineColor(4);
	357	fCutsEff->SetLineWidth(3);
	358	fCutsEff->Draw("same");
	359	leg->AddEntry(fCutsEff,"Cuts efficiency","l");
	360	}
	361	if(fPIDEff) {
	362	fPIDEff->SetLineColor(7);
	363	fPIDEff->SetLineWidth(3);
	364	fPIDEff->Draw("same");
	365	leg->AddEntry(fPIDEff,"PID efficiency","l");
	366	}
	367	if(fMCPtShape) {
5ad8ad40	368	Int_t ci = TColor::GetColor("#9933ff");
	369	fMCPtShape->SetLineColor(ci);
	370	// fMCPtShape->SetLineColor(8);
0233abe6	371	fMCPtShape->SetLineWidth(3);
	372	fMCPtShape->Draw("same");
	373	leg->AddEntry(fMCPtShape,"MC p_{t} shape","l");
	374	}
	375	if(fPartAntipart) {
5ad8ad40	376	Int_t ci = TColor::GetColor("#ff6600");
	377	fPartAntipart->SetLineColor(ci);
	378	// fPartAntipart->SetLineColor(9);
0233abe6	379	fPartAntipart->SetLineWidth(3);
	380	fPartAntipart->Draw("same");
	381	leg->AddEntry(fPartAntipart,"D = #bar{D}","l");
	382	}
0233abe6	383
	384
	385	TH1F *hTotErr=new TH1F("hTotErr","",20,0,20);
	386	for(Int_t i=1;i<=20;i++) {
	387	Double_t pt=hTotErr->GetBinCenter(i);
5ad8ad40	388	if(grErrFeeddown) {
	389	Double_t x=0., y=0., erryh = 0.;
	390	for(Int_t j=0; j<grErrFeeddown->GetN(); j++) {
	391	grErrFeeddown->GetPoint(j,x,y);
	392	if ( TMath::Abs(x - pt) < 0.5) erryh = grErrFeeddown->GetErrorYhigh(j);
	393	}
	394	hTotErr->SetBinContent(i,GetTotalSystErr(pt,erryh));
	395	} else {
	396	hTotErr->SetBinContent(i,GetTotalSystErr(pt));
	397	}
0233abe6	398	}
	399	hTotErr->SetLineColor(1);
	400	hTotErr->SetLineWidth(3);
	401	hTotErr->Draw("same");
	402	leg->AddEntry(hTotErr,"Total (excl. norm.)","l");
	403
	404
	405	leg->Draw();
	406
5ad8ad40	407	cSystErr->SaveAs("RelativeSystematics.eps");
5ad8ad40	408
0233abe6	409	return;
0233abe6	410	}