[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.                  *
 **************************************************************************/

/* $Id$ */

/////////////////////////////////////////////////////////////
//
// 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=1;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=1;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=1;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.10); // 10%

  // 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);
  for(Int_t i=3;i<=6;i++) fPartAntipart->SetBinContent(i,0.08);
  
  return;
}
//--------------------------------------------------------------------------
void AliHFSystErr::InitDplustoKpipi() {
  // 
  // D+->Kpipi syst errors. Responsible: R. Bala
  //

 // Normalization
  fNorm = new TH1F("fNorm","fNorm",20,0,20);
  for(Int_t i=1;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=1;i<=20;i++) fBR->SetBinContent(i,0.04); // 4% PDG2010

  // Tracking efficiency
  fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
  for(Int_t i=1;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.20);
  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: A. Grelli, Y. Wang
  //

 // Normalization
  fNorm = new TH1F("fNorm","fNorm",20,0,20);
  for(Int_t i=1;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=1;i<=20;i++) fBR->SetBinContent(i,0.015); // 1.5% PDG2010

  // Tracking efficiency
  fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
  fTrackingEff->SetBinContent(1,0.12);
  fTrackingEff->SetBinContent(2,0.08);
  fTrackingEff->SetBinContent(3,0.05);
  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,0.12);
  fRawYield->SetBinContent(3,0.09);
  fRawYield->SetBinContent(4,0.08);
  fRawYield->SetBinContent(5,0.06);
  for(Int_t i=5;i<=20;i++) fRawYield->SetBinContent(i,0.04);  //4%

  // 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.04); // 3%
  fPIDEff->SetBinContent(1,1); // 100%
  fPIDEff->SetBinContent(2,1); // 100%
  fPIDEff->SetBinContent(3,0.05); // 5%
  fPIDEff->SetBinContent(4,0.05); // 5%
  fPIDEff->SetBinContent(5,0.05); // 5%
 

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

  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",300,80,640,500);
  cSystErr->Range(0.20,-0.5,18.4,0.34);
  cSystErr->SetRightMargin(0.318);
  cSystErr->SetFillColor(0);

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

  TLegend *leg = new TLegend(0.69,0.44,0.98,0.86,NULL,"brNDC");
  leg->SetTextSize(0.03601695);
  leg->SetFillStyle(0);
  leg->SetBorderSize(0);
  
  TH1F *hTotErr=new TH1F("hTotErr","",20,0,20);
  Int_t nbins = fNorm->GetNbinsX();
  TGraphAsymmErrors *gTotErr = new TGraphAsymmErrors(nbins);
  for(Int_t i=1;i<=20;i++) {
    Double_t pt = hTotErr->GetBinCenter(i);
    Double_t ptwidth = hTotErr->GetBinWidth(i);

    if(grErrFeeddown) {
      Double_t x=0., y=0., errxl=0., errxh=0., erryl=0., erryh=0.;
      Double_t toterryl=0., toterryh=0.;
      for(Int_t j=0; j<grErrFeeddown->GetN(); j++) {
	grErrFeeddown->GetPoint(j,x,y);
	errxh = grErrFeeddown->GetErrorXhigh(j);
	errxl = grErrFeeddown->GetErrorXlow(j);
	if ( ( (x-errxl) <= pt) && ( (x+errxl) >= pt) ) {
	  erryh = grErrFeeddown->GetErrorYhigh(j);
	  erryl = grErrFeeddown->GetErrorYlow(j);
	}
      }
      if (erryl>=1e-3) toterryl = GetTotalSystErr(pt,erryl);
      else toterryl = GetTotalSystErr(pt);
      if (erryh>=1e-3) toterryh = GetTotalSystErr(pt,erryh);
      else toterryh = GetTotalSystErr(pt);

      hTotErr->SetBinContent(i,toterryh);
      gTotErr->SetPoint(i,pt,0.);
      gTotErr->SetPointError(i,ptwidth/2.,ptwidth/2.,toterryl,toterryh); // i, exl, exh, eyl, eyh
    }
    else {
      hTotErr->SetBinContent(i,GetTotalSystErr(pt));
      gTotErr->SetPoint(i,pt,0.);
      gTotErr->SetPointError(i,ptwidth/2.,ptwidth/2.,GetTotalSystErr(pt),GetTotalSystErr(pt)); // i, exl, exh, eyl, eyh
    }

  }
  gTotErr->SetLineColor(kBlack);
  gTotErr->SetFillColor(kRed);
  gTotErr->SetFillStyle(3002);
  gTotErr->Draw("2");
  leg->AddEntry(gTotErr,"Total (excl. norm.)","f");
//   hTotErr->SetLineColor(1);
//   hTotErr->SetLineWidth(3);
//   hTotErr->Draw("same");
//   leg->AddEntry(hTotErr,"Total (excl. norm.)","l");
  

  fNorm->SetFillColor(1);
  fNorm->SetFillStyle(3002);
  //fNorm->Draw("same");
  //TH1F *hNormRefl = ReflectHisto(fNorm);
  //hNormRefl->Draw("same");
  leg->AddEntry(fNorm,"Normalization (10%)","");

  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");
    TH1F *hTrackingEffRefl = ReflectHisto(fTrackingEff);
    hTrackingEffRefl->Draw("same");
    leg->AddEntry(fTrackingEff,"Tracking efficiency","f");
  }
  if(fBR) {
    fBR->SetFillColor(6);
    fBR->SetFillStyle(3005);
    //fBR->SetFillStyle(3020);
    fBR->Draw("same");
    TH1F *hBRRefl = ReflectHisto(fBR);
    hBRRefl->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");
    TH1F *hRawYieldRefl = ReflectHisto(fRawYield);
    hRawYieldRefl->Draw("same");
    leg->AddEntry(fRawYield,"Yield extraction","l");
  }
  if(fCutsEff) {
    fCutsEff->SetLineColor(4);
    fCutsEff->SetLineWidth(3);
    fCutsEff->Draw("same");
    TH1F *hCutsEffRefl = ReflectHisto(fCutsEff);
    hCutsEffRefl->Draw("same");
    leg->AddEntry(fCutsEff,"Cuts efficiency","l");
  }
  if(fPIDEff) {
    fPIDEff->SetLineColor(7);
    fPIDEff->SetLineWidth(3);
    fPIDEff->Draw("same");
    TH1F *hPIDEffRefl = ReflectHisto(fPIDEff);
    hPIDEffRefl->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");
    TH1F *hMCPtShapeRefl = ReflectHisto(fMCPtShape);
    hMCPtShapeRefl->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");
    TH1F *hPartAntipartRefl = ReflectHisto(fPartAntipart);
    hPartAntipartRefl->Draw("same");
    leg->AddEntry(fPartAntipart,"D = #bar{D}","l");
  }


  leg->Draw();

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

  return;
}
//-------------------------------------------------------------------------
TH1F* AliHFSystErr::ReflectHisto(TH1F *hin) const {
  //
  // Clones and reflects histogram 
  // 
  TH1F *hout=(TH1F*)hin->Clone("hout");
  hout->Scale(-1.);

  return hout;
}
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
27de2dfb	16	/* $Id$ */
27de2dfb	17
0233abe6	18	/////////////////////////////////////////////////////////////
	19	//
	20	// Class to handle systematic errors for charm hadrons
	21	//
	22	// Usage:
	23	// AliHFSystEff syst(DECAY); // DECAY = 1 for D0, 2, for D+, 3 for D*
	24	// syst.DrawErrors(); // to see a plot of the error contributions
	25	// syst.GetTotalSystErr(pt); // to get the total err at pt
	26	//
	27	// Author: A.Dainese, andrea.dainese@pd.infn.it
	28	/////////////////////////////////////////////////////////////
	29
	30	#include <TStyle.h>
	31	#include <TGraphAsymmErrors.h>
	32	#include <TMath.h>
	33	#include <TCanvas.h>
	34	#include <TH2F.h>
	35	#include <TLegend.h>
5ad8ad40	36	#include <TColor.h>
0233abe6	37
	38	#include "AliHFSystErr.h"
	39
	40
	41	ClassImp(AliHFSystErr)
	42
	43	//--------------------------------------------------------------------------
	44	AliHFSystErr::AliHFSystErr(const Char_t* name, const Char_t* title) :
	45	TNamed(name,title),
	46	fNorm(0),
	47	fRawYield(0),
	48	fTrackingEff(0),
	49	fBR(0),
	50	fCutsEff(0),
	51	fPIDEff(0),
	52	fMCPtShape(0),
	53	fPartAntipart(0)
	54	{
	55	//
	56	// Default Constructor
	57	//
	58	}
	59	//--------------------------------------------------------------------------
	60	AliHFSystErr::AliHFSystErr(Int_t decay,const Char_t* name, const Char_t* title) :
	61	TNamed(name,title),
	62	fNorm(0),
	63	fRawYield(0),
	64	fTrackingEff(0),
	65	fBR(0),
	66	fCutsEff(0),
	67	fPIDEff(0),
	68	fMCPtShape(0),
	69	fPartAntipart(0)
	70	{
	71	//
	72	// Standard Constructor
	73	//
	74
	75	switch(decay) {
	76	case 1: // D0->Kpi
	77	InitD0toKpi();
	78	break;
	79	case 2: // D+->Kpipi
	80	InitDplustoKpipi();
	81	break;
	82	case 3: // D*->D0pi
	83	InitDstartoD0pi();
	84	break;
	85	default:
	86	printf("Invalid decay type: %d\n",decay);
	87	break;
	88	}
	89
	90
	91	}
	92	//--------------------------------------------------------------------------
	93	AliHFSystErr::~AliHFSystErr() {
	94	//
	95	// Default Destructor
	96	//
	97
	98	if(fNorm) { delete fNorm; fNorm=0; }
	99	if(fRawYield) { delete fRawYield; fRawYield=0; }
	100	if(fTrackingEff) { delete fTrackingEff; fTrackingEff=0; }
101	if(fBR) { delete fBR; fBR=0; }
102	if(fCutsEff) { delete fCutsEff; fCutsEff=0; }
103	if(fPIDEff) { delete fPIDEff; fPIDEff=0; }
104	if(fMCPtShape) { delete fMCPtShape; fMCPtShape=0; }
105	if(fPartAntipart) { delete fPartAntipart; fPartAntipart=0; }
106
107	}
108	//--------------------------------------------------------------------------
109	void AliHFSystErr::InitD0toKpi() {
110	//
111	// D0->Kpi syst errors. Responsible: A. Rossi
112	//
113
114	// Normalization
115	fNorm = new TH1F("fNorm","fNorm",20,0,20);
844405bc	116	for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
0233abe6	117
	118	// Branching ratio
	119	fBR = new TH1F("fBR","fBR",20,0,20);
844405bc	120	for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.012); // 1.2% PDG2010
0233abe6	121
	122	// Tracking efficiency
	123	fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
844405bc	124	for(Int_t i=1;i<=20;i++) fTrackingEff->SetBinContent(i,0.02); // 2% (1% per track)
0233abe6	125
	126	// Raw yield extraction
	127	fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
	128	fRawYield->SetBinContent(1,1);
	129	fRawYield->SetBinContent(2,1);
	130	fRawYield->SetBinContent(3,0.15);
5ad8ad40	131	for(Int_t i=4;i<=20;i++) fRawYield->SetBinContent(i,0.065);
0233abe6	132
	133	// Cuts efficiency (from cuts variation)
	134	fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
	135	for(Int_t i=1;i<=20;i++) fCutsEff->SetBinContent(i,0.10); // 10%
	136
	137	// PID efficiency (from PID/noPID)
	138	fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
	139	for(Int_t i=1;i<=20;i++) fPIDEff->SetBinContent(i,0.03); // 3%
844405bc	140	fPIDEff->SetBinContent(4,0.10); // 10%
0233abe6	141
	142	// MC dN/dpt
	143	fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
	144	for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,(Float_t)i*0.006);
	145
	146	// particle-antiparticle
	147	fPartAntipart = new TH1F("fPartAntipart","fPartAntipart",20,0,20);
	148	fPartAntipart->SetBinContent(1,1);
	149	fPartAntipart->SetBinContent(2,1);
844405bc	150	for(Int_t i=3;i<=6;i++) fPartAntipart->SetBinContent(i,0.08);
0233abe6	151
	152	return;
	153	}
	154	//--------------------------------------------------------------------------
	155	void AliHFSystErr::InitDplustoKpipi() {
	156	//
	157	// D+->Kpipi syst errors. Responsible: R. Bala
	158	//
	159
626bcf57	160	// Normalization
626bcf57	161	fNorm = new TH1F("fNorm","fNorm",20,0,20);
844405bc	162	for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
626bcf57	163
	164	// Branching ratio
	165	fBR = new TH1F("fBR","fBR",20,0,20);
844405bc	166	for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.04); // 4% PDG2010
626bcf57	167
	168	// Tracking efficiency
	169	fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
844405bc	170	for(Int_t i=1;i<=20;i++) fTrackingEff->SetBinContent(i,0.03); // 3% (1% per track)
626bcf57	171
	172
	173	// Raw yield extraction
	174	fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
	175	fRawYield->SetBinContent(1,1);
	176	fRawYield->SetBinContent(2,1);
5c73c2f5	177	fRawYield->SetBinContent(3,0.20);
5ad8ad40	178	for(Int_t i=4;i<=20;i++) fRawYield->SetBinContent(i,0.055); //5 to 10%
626bcf57	179
	180	// Cuts efficiency (from cuts variation)
	181	fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
	182	for(Int_t i=1;i<=20;i++) fCutsEff->SetBinContent(i,0.10); // 10%
	183
	184	// PID efficiency (from PID/noPID)
	185	fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
	186	for(Int_t i=1;i<=20;i++) fPIDEff->SetBinContent(i,0.05); // 5%
	187	fPIDEff->SetBinContent(3,0.13); // 13%
	188
	189
	190	// MC dN/dpt (copied from D0 : will update later)
	191	fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
	192	for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,(Float_t)i*0.006);
	193
	194
	195	// particle-antiparticle
	196	/*
	197	fPartAntipart = new TH1F("fPartAntipart","fPartAntipart",20,0,20);
	198	fPartAntipart->SetBinContent(1,1);
	199	fPartAntipart->SetBinContent(2,1);
	200	fPartAntipart->SetBinContent(3,0.12);
	201	for(Int_t i=4;i<=20;i++) fPartAntipart->SetBinContent(i,0.05); //5 to 12%
	202	*/
0233abe6	203	return;
	204	}
	205	//--------------------------------------------------------------------------
	206	void AliHFSystErr::InitDstartoD0pi() {
	207	//
2cf98d0d	208	// D*+->D0pi syst errors. Responsible: A. Grelli, Y. Wang
0233abe6	209	//
0233abe6	210
2cf98d0d	211	// Normalization
	212	fNorm = new TH1F("fNorm","fNorm",20,0,20);
	213	for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
	214
	215	// Branching ratio
	216	fBR = new TH1F("fBR","fBR",20,0,20);
8accac87	217	for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.015); // 1.5% PDG2010
2cf98d0d	218
	219	// Tracking efficiency
	220	fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
	221	fTrackingEff->SetBinContent(1,0.12);
	222	fTrackingEff->SetBinContent(2,0.08);
8accac87	223	fTrackingEff->SetBinContent(3,0.05);
2cf98d0d	224	for(Int_t i=4;i<=20;i++) fTrackingEff->SetBinContent(i,0.03); // 3% (1% per track)
	225
	226
	227	// Raw yield extraction
	228	fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
8accac87	229	fRawYield->SetBinContent(1,1);
	230	fRawYield->SetBinContent(2,0.12);
	231	fRawYield->SetBinContent(3,0.09);
	232	fRawYield->SetBinContent(4,0.08);
	233	fRawYield->SetBinContent(5,0.06);
	234	for(Int_t i=5;i<=20;i++) fRawYield->SetBinContent(i,0.04); //4%
2cf98d0d	235
	236	// Cuts efficiency (from cuts variation)
	237	fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
	238	for(Int_t i=1;i<=20;i++) fCutsEff->SetBinContent(i,0.10); // 10%
	239
	240	// PID efficiency (from PID/noPID)
	241	fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
8accac87	242	for(Int_t i=1;i<=20;i++) fPIDEff->SetBinContent(i,0.04); // 3%
	243	fPIDEff->SetBinContent(1,1); // 100%
	244	fPIDEff->SetBinContent(2,1); // 100%
	245	fPIDEff->SetBinContent(3,0.05); // 5%
	246	fPIDEff->SetBinContent(4,0.05); // 5%
	247	fPIDEff->SetBinContent(5,0.05); // 5%
2cf98d0d	248
	249
	250	// MC dN/dpt (copied from D0 : will update later)
	251	fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
	252	for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,(Float_t)i*0.006);
	253
0233abe6	254	return;
	255	}
	256	//--------------------------------------------------------------------------
	257	Double_t AliHFSystErr::GetCutsEffErr(Double_t pt) const {
	258	//
	259	// Get error
	260	//
	261
	262	Int_t bin=fCutsEff->FindBin(pt);
	263
	264	return fCutsEff->GetBinContent(bin);
	265	}
	266	//--------------------------------------------------------------------------
	267	Double_t AliHFSystErr::GetMCPtShapeErr(Double_t pt) const {
	268	//
	269	// Get error
	270	//
	271
	272	Int_t bin=fMCPtShape->FindBin(pt);
	273
	274	return fMCPtShape->GetBinContent(bin);
	275	}
	276	//--------------------------------------------------------------------------
	277	Double_t AliHFSystErr::GetSeleEffErr(Double_t pt) const {
	278	//
	279	// Get error
	280	//
	281
	282	Double_t err=GetCutsEffErr(pt)GetCutsEffErr(pt)+GetMCPtShapeErr(pt)GetMCPtShapeErr(pt);
	283
	284	return TMath::Sqrt(err);
	285	}
	286	//--------------------------------------------------------------------------
	287	Double_t AliHFSystErr::GetPIDEffErr(Double_t pt) const {
	288	//
	289	// Get error
	290	//
	291
	292	Int_t bin=fPIDEff->FindBin(pt);
	293
	294	return fPIDEff->GetBinContent(bin);
	295	}
	296	//--------------------------------------------------------------------------
	297	Double_t AliHFSystErr::GetTrackingEffErr(Double_t pt) const {
	298	//
	299	// Get error
	300	//
	301
	302	Int_t bin=fTrackingEff->FindBin(pt);
	303
	304	return fTrackingEff->GetBinContent(bin);
	305	}
	306	//--------------------------------------------------------------------------
	307	Double_t AliHFSystErr::GetRawYieldErr(Double_t pt) const {
	308	//
	309	// Get error
	310	//
	311
	312	Int_t bin=fRawYield->FindBin(pt);
	313
	314	return fRawYield->GetBinContent(bin);
	315	}
	316	//--------------------------------------------------------------------------
	317	Double_t AliHFSystErr::GetPartAntipartErr(Double_t pt) const {
318	//
319	// Get error
320	//
321
322	Int_t bin=fPartAntipart->FindBin(pt);
323
324	return fPartAntipart->GetBinContent(bin);
325	}
326	//--------------------------------------------------------------------------
327	Double_t AliHFSystErr::GetTotalSystErr(Double_t pt,Double_t feeddownErr) const {
328	//
329	// Get total syst error (except norm. error)
330	//
331
332	Double_t err=0.;
333
334	if(fRawYield) err += GetRawYieldErr(pt)*GetRawYieldErr(pt);
335	if(fTrackingEff) err += GetTrackingEffErr(pt)*GetTrackingEffErr(pt);
336	if(fBR) err += GetBRErr()*GetBRErr();
337	if(fCutsEff) err += GetCutsEffErr(pt)*GetCutsEffErr(pt);
338	if(fPIDEff) err += GetPIDEffErr(pt)*GetPIDEffErr(pt);
339	if(fMCPtShape) err += GetMCPtShapeErr(pt)*GetMCPtShapeErr(pt);
340	if(fPartAntipart) err += GetPartAntipartErr(pt)*GetPartAntipartErr(pt);
341
342	err += feeddownErr*feeddownErr;
343
344	return TMath::Sqrt(err);
345	}
346	//---------------------------------------------------------------------------
5ad8ad40	347	void AliHFSystErr::DrawErrors(TGraphAsymmErrors *grErrFeeddown) const {
0233abe6	348	//
	349	// Draw errors
	350	//
	351	gStyle->SetOptStat(0);
	352
a07118d0	353	TCanvas *cSystErr = new TCanvas("cSystErr","Systematic Errors",300,80,640,500);
	354	cSystErr->Range(0.20,-0.5,18.4,0.34);
	355	cSystErr->SetRightMargin(0.318);
0233abe6	356	cSystErr->SetFillColor(0);
0233abe6	357
2ee3afe2	358	TH2F *hFrame = new TH2F("hFrame","Systematic errors; p_{t} [GeV/c]; Relative Error",20,0,20,100,-1,+1);
5ad8ad40	359	hFrame->SetAxisRange(2.,11.9,"X");
a07118d0	360	hFrame->SetAxisRange(-0.5,0.5,"Y");
0233abe6	361	hFrame->Draw();
0233abe6	362
a07118d0	363	TLegend *leg = new TLegend(0.69,0.44,0.98,0.86,NULL,"brNDC");
a07118d0	364	leg->SetTextSize(0.03601695);
0233abe6	365	leg->SetFillStyle(0);
0233abe6	366	leg->SetBorderSize(0);
2ee3afe2	367
	368	TH1F *hTotErr=new TH1F("hTotErr","",20,0,20);
	369	Int_t nbins = fNorm->GetNbinsX();
	370	TGraphAsymmErrors *gTotErr = new TGraphAsymmErrors(nbins);
	371	for(Int_t i=1;i<=20;i++) {
	372	Double_t pt = hTotErr->GetBinCenter(i);
	373	Double_t ptwidth = hTotErr->GetBinWidth(i);
	374
	375	if(grErrFeeddown) {
	376	Double_t x=0., y=0., errxl=0., errxh=0., erryl=0., erryh=0.;
	377	Double_t toterryl=0., toterryh=0.;
	378	for(Int_t j=0; j<grErrFeeddown->GetN(); j++) {
	379	grErrFeeddown->GetPoint(j,x,y);
	380	errxh = grErrFeeddown->GetErrorXhigh(j);
	381	errxl = grErrFeeddown->GetErrorXlow(j);
a07118d0	382	if ( ( (x-errxl) <= pt) && ( (x+errxl) >= pt) ) {
2ee3afe2	383	erryh = grErrFeeddown->GetErrorYhigh(j);
	384	erryl = grErrFeeddown->GetErrorYlow(j);
	385	}
	386	}
	387	if (erryl>=1e-3) toterryl = GetTotalSystErr(pt,erryl);
	388	else toterryl = GetTotalSystErr(pt);
	389	if (erryh>=1e-3) toterryh = GetTotalSystErr(pt,erryh);
	390	else toterryh = GetTotalSystErr(pt);
	391
	392	hTotErr->SetBinContent(i,toterryh);
	393	gTotErr->SetPoint(i,pt,0.);
	394	gTotErr->SetPointError(i,ptwidth/2.,ptwidth/2.,toterryl,toterryh); // i, exl, exh, eyl, eyh
	395	}
	396	else {
	397	hTotErr->SetBinContent(i,GetTotalSystErr(pt));
	398	gTotErr->SetPoint(i,pt,0.);
	399	gTotErr->SetPointError(i,ptwidth/2.,ptwidth/2.,GetTotalSystErr(pt),GetTotalSystErr(pt)); // i, exl, exh, eyl, eyh
	400	}
	401
	402	}
f968ef30	403	gTotErr->SetLineColor(kBlack);
	404	gTotErr->SetFillColor(kRed);
	405	gTotErr->SetFillStyle(3002);
2ee3afe2	406	gTotErr->Draw("2");
	407	leg->AddEntry(gTotErr,"Total (excl. norm.)","f");
	408	// hTotErr->SetLineColor(1);
	409	// hTotErr->SetLineWidth(3);
	410	// hTotErr->Draw("same");
	411	// leg->AddEntry(hTotErr,"Total (excl. norm.)","l");
	412
0233abe6	413
e373b9f0	414	fNorm->SetFillColor(1);
	415	fNorm->SetFillStyle(3002);
	416	//fNorm->Draw("same");
	417	//TH1F *hNormRefl = ReflectHisto(fNorm);
	418	//hNormRefl->Draw("same");
	419	leg->AddEntry(fNorm,"Normalization (10%)","");
	420
5ad8ad40	421	if(grErrFeeddown) {
	422	grErrFeeddown->SetFillColor(kTeal-8);
	423	grErrFeeddown->SetFillStyle(3001);
	424	grErrFeeddown->Draw("2");
	425	leg->AddEntry(grErrFeeddown,"Feed-down from B","f");
	426	}
0233abe6	427	if(fTrackingEff) {
5ad8ad40	428	fTrackingEff->SetFillColor(4);
5ad8ad40	429	fTrackingEff->SetFillStyle(3006);
0233abe6	430	fTrackingEff->Draw("same");
f968ef30	431	TH1F *hTrackingEffRefl = ReflectHisto(fTrackingEff);
f968ef30	432	hTrackingEffRefl->Draw("same");
0233abe6	433	leg->AddEntry(fTrackingEff,"Tracking efficiency","f");
	434	}
	435	if(fBR) {
	436	fBR->SetFillColor(6);
	437	fBR->SetFillStyle(3005);
5ad8ad40	438	//fBR->SetFillStyle(3020);
0233abe6	439	fBR->Draw("same");
f968ef30	440	TH1F *hBRRefl = ReflectHisto(fBR);
f968ef30	441	hBRRefl->Draw("same");
0233abe6	442	leg->AddEntry(fBR,"Branching ratio","f");
	443	}
	444	if(fRawYield) {
5ad8ad40	445	Int_t ci; // for color index setting
	446	ci = TColor::GetColor("#00cc00");
	447	fRawYield->SetLineColor(ci);
	448	// fRawYield->SetLineColor(3);
0233abe6	449	fRawYield->SetLineWidth(3);
0233abe6	450	fRawYield->Draw("same");
f968ef30	451	TH1F *hRawYieldRefl = ReflectHisto(fRawYield);
f968ef30	452	hRawYieldRefl->Draw("same");
844405bc	453	leg->AddEntry(fRawYield,"Yield extraction","l");
0233abe6	454	}
	455	if(fCutsEff) {
	456	fCutsEff->SetLineColor(4);
	457	fCutsEff->SetLineWidth(3);
	458	fCutsEff->Draw("same");
f968ef30	459	TH1F *hCutsEffRefl = ReflectHisto(fCutsEff);
f968ef30	460	hCutsEffRefl->Draw("same");
0233abe6	461	leg->AddEntry(fCutsEff,"Cuts efficiency","l");
	462	}
	463	if(fPIDEff) {
	464	fPIDEff->SetLineColor(7);
	465	fPIDEff->SetLineWidth(3);
	466	fPIDEff->Draw("same");
f968ef30	467	TH1F *hPIDEffRefl = ReflectHisto(fPIDEff);
f968ef30	468	hPIDEffRefl->Draw("same");
0233abe6	469	leg->AddEntry(fPIDEff,"PID efficiency","l");
	470	}
	471	if(fMCPtShape) {
5ad8ad40	472	Int_t ci = TColor::GetColor("#9933ff");
	473	fMCPtShape->SetLineColor(ci);
	474	// fMCPtShape->SetLineColor(8);
0233abe6	475	fMCPtShape->SetLineWidth(3);
0233abe6	476	fMCPtShape->Draw("same");
f968ef30	477	TH1F *hMCPtShapeRefl = ReflectHisto(fMCPtShape);
f968ef30	478	hMCPtShapeRefl->Draw("same");
0233abe6	479	leg->AddEntry(fMCPtShape,"MC p_{t} shape","l");
	480	}
	481	if(fPartAntipart) {
5ad8ad40	482	Int_t ci = TColor::GetColor("#ff6600");
	483	fPartAntipart->SetLineColor(ci);
	484	// fPartAntipart->SetLineColor(9);
0233abe6	485	fPartAntipart->SetLineWidth(3);
0233abe6	486	fPartAntipart->Draw("same");
f968ef30	487	TH1F *hPartAntipartRefl = ReflectHisto(fPartAntipart);
f968ef30	488	hPartAntipartRefl->Draw("same");
0233abe6	489	leg->AddEntry(fPartAntipart,"D = #bar{D}","l");
0233abe6	490	}
0233abe6	491
0233abe6	492
	493	leg->Draw();
	494
5ad8ad40	495	cSystErr->SaveAs("RelativeSystematics.eps");
5ad8ad40	496
0233abe6	497	return;
0233abe6	498	}
f968ef30	499	//-------------------------------------------------------------------------
	500	TH1F* AliHFSystErr::ReflectHisto(TH1F *hin) const {
	501	//
	502	// Clones and reflects histogram
	503	//
	504	TH1F hout=(TH1F)hin->Clone("hout");
	505	hout->Scale(-1.);
	506
	507	return hout;
	508	}