PWG3/vertexingHF/AliHFSystErr.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2010, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 /////////////////////////////////////////////////////////////
  19 //
  20 // Class to handle systematic errors for charm hadrons
  21 //
  22 // Usage:
  23 // AliHFSystEff syst;           // DECAY = 1 for D0, 2, for D+, 3 for D*
  24 // syst.SetRunNumber(YEAR);     // YEAR = two last numbers of the year (is 10 for 2010)
  25 // syst.SetCollisionType(TYPE);  // TYPE =  0 is pp, 1 is PbPb
  26 // syst.SetCentrality(CENT);     // CENT is centrality, 0100 for MB, 020 (4080) for 0-20 (40-80) CC...
  27 // syst.Init(DECAY);             // DECAY = 1 for D0, 2, for D+, 3 for D*
  28 // syst.DrawErrors(); // to see a plot of the error contributions
  29 // syst.GetTotalSystErr(pt); // to get the total err at pt
  30 //
  31 // Author: A.Dainese, andrea.dainese@pd.infn.it
  32 /////////////////////////////////////////////////////////////
  33
  34 #include <TStyle.h>
  35 #include <TGraphAsymmErrors.h>
  36 #include <TMath.h>
  37 #include <TCanvas.h>
  38 #include <TH2F.h>
  39 #include <TLegend.h>
  40 #include <TColor.h>
  41
  42 #include "AliLog.h"
  43 #include "AliHFSystErr.h"
  44
  45
  46 ClassImp(AliHFSystErr)
  47
  48 //--------------------------------------------------------------------------
  49 AliHFSystErr::AliHFSystErr(const Char_t* name, const Char_t* title) :
  50 TNamed(name,title),
  51 fNorm(0),
  52 fRawYield(0),
  53 fTrackingEff(0),
  54 fBR(0),
  55 fCutsEff(0),
  56 fPIDEff(0),
  57 fMCPtShape(0),
  58 fPartAntipart(0),
  59 fRunNumber(10),
  60 fCollisionType(0),
  61 fCentralityClass(0100)
  62 {
  63   //
  64   // Default Constructor
  65   //
  66 }
  67
  68 //--------------------------------------------------------------------------
  69 AliHFSystErr::~AliHFSystErr() {
  70   //
  71   // Default Destructor
  72   //
  73
  74   if(fNorm)         { delete fNorm; fNorm=0; }
  75   if(fRawYield)     { delete fRawYield; fRawYield=0; }
  76   if(fTrackingEff)  { delete fTrackingEff; fTrackingEff=0; }
  77   if(fBR)           { delete fBR; fBR=0; }
  78   if(fCutsEff)      { delete fCutsEff; fCutsEff=0; }
  79   if(fPIDEff)       { delete fPIDEff; fPIDEff=0; }
  80   if(fMCPtShape)    { delete fMCPtShape; fMCPtShape=0; }
  81   if(fPartAntipart) { delete fPartAntipart; fPartAntipart=0; }
  82
  83 }
  84
  85 //--------------------------------------------------------------------------
  86 void AliHFSystErr::Init(Int_t decay){
  87   //
  88   // Variables/histos initialization
  89   //
  90
  91   if (fRunNumber!=10) {
  92     AliError("Only settings for 2010 are implemented so far");
  93   }
  94   if (fCentralityClass!=020 && fCentralityClass!=4080 && fCentralityClass!=0100){
  95     AliError("Only settings for MB2010 are implemented so far");
  96   }
  97
  98   switch(decay) {
  99   case 1: // D0->Kpi
 100     if (fCollisionType==0) InitD0toKpi2010pp();
 101     else if (fCollisionType==1) {
 102       if (fCentralityClass==020) InitD0toKpi2010PbPb020();
 103       else if (fCentralityClass==4080) InitD0toKpi2010PbPb4080();
 104       else AliError("Not yet implemented");
 105     }
 106     break;
 107   case 2: // D+->Kpipi
 108     if (fCollisionType==0) InitDplustoKpipi2010pp();
 109     else if (fCollisionType==1) {
 110       if (fCentralityClass==020) InitDplustoKpipi2010PbPb020();
 111       else if (fCentralityClass==4080) InitDplustoKpipi2010PbPb4080();
 112       else AliError("Not yet implemented");
 113     }
 114     break;
 115   case 3: // D*->D0pi
 116     if (fCollisionType==0) InitDstartoD0pi2010pp();
 117     else if (fCollisionType==1) {
 118       if (fCentralityClass==020) InitDstartoD0pi2010PbPb020();
 119       else if (fCentralityClass==4080) InitDstartoD0pi2010PbPb4080();
 120       else AliError("Not yet implemented");
 121     }
 122     break;
 123   default:
 124     printf("Invalid decay type: %d\n",decay);
 125     break;
 126   }
 127
 128 }
 129
 130 //--------------------------------------------------------------------------
 131 void AliHFSystErr::InitD0toKpi2010pp() {
 132   //
 133   // D0->Kpi syst errors. Responsible: A. Rossi
 134   //   2010 pp sample
 135   //
 136
 137   // Normalization
 138   fNorm = new TH1F("fNorm","fNorm",20,0,20);
 139   for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
 140
 141   // Branching ratio
 142   fBR = new TH1F("fBR","fBR",20,0,20);
 143   for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.012); // 1.2% PDG2010
 144
 145   // Tracking efficiency
 146   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
 147   for(Int_t i=1;i<=20;i++) fTrackingEff->SetBinContent(i,0.02); // 2% (1% per track)
 148
 149   // Raw yield extraction
 150   fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
 151   fRawYield->SetBinContent(1,1);
 152   fRawYield->SetBinContent(2,1);
 153   fRawYield->SetBinContent(3,0.15);
 154   for(Int_t i=4;i<=20;i++) fRawYield->SetBinContent(i,0.065);
 155
 156   // Cuts efficiency (from cuts variation)
 157   fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
 158   for(Int_t i=1;i<=20;i++) fCutsEff->SetBinContent(i,0.10); // 10%
 159
 160   // PID efficiency (from PID/noPID)
 161   fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
 162   for(Int_t i=1;i<=20;i++) fPIDEff->SetBinContent(i,0.03); // 3%
 163   fPIDEff->SetBinContent(4,0.10); // 10%
 164
 165   // MC dN/dpt
 166   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
 167   for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,(Float_t)i*0.006);
 168
 169   // particle-antiparticle
 170   fPartAntipart = new TH1F("fPartAntipart","fPartAntipart",20,0,20);
 171   fPartAntipart->SetBinContent(1,1);
 172   fPartAntipart->SetBinContent(2,1);
 173   for(Int_t i=3;i<=6;i++) fPartAntipart->SetBinContent(i,0.08);
 174
 175   return;
 176 }
 177 //--------------------------------------------------------------------------
 178 void AliHFSystErr::InitD0toKpi2010PbPb020() {
 179   //
 180   // D0->Kpi syst errors. Responsible: ??
 181   //   2010 PbPb sample, 0-20 CC
 182   //
 183
 184   // Normalization
 185   fNorm = new TH1F("fNorm","fNorm",20,0,20);
 186   for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
 187
 188   // Branching ratio
 189   fBR = new TH1F("fBR","fBR",20,0,20);
 190   for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.012); // 1.2% PDG2010
 191
 192   // Tracking efficiency
 193   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
 194   for(Int_t i=1;i<=20;i++) fTrackingEff->SetBinContent(i,0.02); // 2% (1% per track)
 195
 196   // Raw yield extraction
 197   fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
 198   fRawYield->SetBinContent(1,0);
 199   fRawYield->SetBinContent(2,0);
 200   fRawYield->SetBinContent(3,0);
 201   fRawYield->SetBinContent(4,0.27);
 202   fRawYield->SetBinContent(5,0.27);
 203   fRawYield->SetBinContent(6,0);
 204   fRawYield->SetBinContent(7,0.07);
 205   fRawYield->SetBinContent(8,0.07);
 206   for(Int_t i=9;i<=12;i++) fRawYield->SetBinContent(i,0.09);
 207   for(Int_t i=12;i<=20;i++) fRawYield->SetBinContent(i,0);
 208
 209   // Cuts efficiency (from cuts variation)
 210   fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
 211   for(Int_t i=1;i<=20;i++) fCutsEff->SetBinContent(i,0.10); // 10%
 212
 213   // PID efficiency (from PID/noPID)
 214   fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
 215   for(Int_t i=1;i<=20;i++) fPIDEff->SetBinContent(i,0.03); // 3%
 216   fPIDEff->SetBinContent(4,0.10); // 10%
 217
 218   // MC dN/dpt
 219   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
 220   for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,(Float_t)i*0.006);
 221
 222   // particle-antiparticle
 223   fPartAntipart = new TH1F("fPartAntipart","fPartAntipart",20,0,20);
 224   fPartAntipart->SetBinContent(1,1);
 225   fPartAntipart->SetBinContent(2,1);
 226   for(Int_t i=3;i<=6;i++) fPartAntipart->SetBinContent(i,0.08);
 227
 228   return;
 229 }
 230 //--------------------------------------------------------------------------
 231 void AliHFSystErr::InitD0toKpi2010PbPb4080() {
 232   //
 233   // D0->Kpi syst errors. Responsible: ??
 234   //   2010 PbPb sample, 40-80 CC
 235   //
 236
 237   // Normalization
 238   fNorm = new TH1F("fNorm","fNorm",20,0,20);
 239   for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
 240
 241   // Branching ratio
 242   fBR = new TH1F("fBR","fBR",20,0,20);
 243   for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.012); // 1.2% PDG2010
 244
 245   // Tracking efficiency
 246   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
 247   for(Int_t i=1;i<=20;i++) fTrackingEff->SetBinContent(i,0.02); // 2% (1% per track)
 248
 249   // Raw yield extraction
 250   fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
 251   fRawYield->SetBinContent(1,0);
 252   fRawYield->SetBinContent(2,0);
 253   fRawYield->SetBinContent(3,0);
 254   fRawYield->SetBinContent(4,0);
 255   fRawYield->SetBinContent(5,0.3);
 256   fRawYield->SetBinContent(6,0);
 257   fRawYield->SetBinContent(7,0.33);
 258   fRawYield->SetBinContent(8,0.33);
 259   for(Int_t i=9;i<=20;i++) fRawYield->SetBinContent(i,0);
 260
 261   // Cuts efficiency (from cuts variation)
 262   fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
 263   for(Int_t i=1;i<=20;i++) fCutsEff->SetBinContent(i,0.10); // 10%
 264
 265   // PID efficiency (from PID/noPID)
 266   fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
 267   for(Int_t i=1;i<=20;i++) fPIDEff->SetBinContent(i,0.03); // 3%
 268   fPIDEff->SetBinContent(4,0.10); // 10%
 269
 270   // MC dN/dpt
 271   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
 272   for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,(Float_t)i*0.006);
 273
 274   // particle-antiparticle
 275   fPartAntipart = new TH1F("fPartAntipart","fPartAntipart",20,0,20);
 276   fPartAntipart->SetBinContent(1,1);
 277   fPartAntipart->SetBinContent(2,1);
 278   for(Int_t i=3;i<=6;i++) fPartAntipart->SetBinContent(i,0.08);
 279
 280   return;
 281 }
 282 //--------------------------------------------------------------------------
 283 void AliHFSystErr::InitDplustoKpipi2010pp() {
 284   //
 285   // D+->Kpipi syst errors. Responsible: R. Bala
 286   //  2010 pp sample
 287   //
 288
 289  // Normalization
 290   fNorm = new TH1F("fNorm","fNorm",20,0,20);
 291   for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
 292
 293   // Branching ratio
 294   fBR = new TH1F("fBR","fBR",20,0,20);
 295   for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.04); // 4% PDG2010
 296
 297   // Tracking efficiency
 298   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
 299   for(Int_t i=1;i<=20;i++) fTrackingEff->SetBinContent(i,0.03); // 3% (1% per track)
 300
 301
 302   // Raw yield extraction
 303   fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
 304   fRawYield->SetBinContent(1,1);
 305   fRawYield->SetBinContent(2,1);
 306   fRawYield->SetBinContent(3,0.20);
 307   for(Int_t i=4;i<=20;i++) fRawYield->SetBinContent(i,0.055);  //5 to 10%
 308
 309   // Cuts efficiency (from cuts variation)
 310   fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
 311   for(Int_t i=1;i<=20;i++) fCutsEff->SetBinContent(i,0.10); // 10%
 312
 313   // PID efficiency (from PID/noPID)
 314   fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
 315   for(Int_t i=1;i<=20;i++) fPIDEff->SetBinContent(i,0.05); // 5%
 316   fPIDEff->SetBinContent(3,0.13); // 13%
 317
 318
 319   // MC dN/dpt  (copied from D0 : will update later)
 320   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
 321   for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,(Float_t)i*0.006);
 322
 323
 324   // particle-antiparticle
 325   /*
 326   fPartAntipart = new TH1F("fPartAntipart","fPartAntipart",20,0,20);
 327   fPartAntipart->SetBinContent(1,1);
 328   fPartAntipart->SetBinContent(2,1);
 329   fPartAntipart->SetBinContent(3,0.12);
 330   for(Int_t i=4;i<=20;i++) fPartAntipart->SetBinContent(i,0.05);   //5 to 12%
 331   */
 332   return;
 333 }
 334 //--------------------------------------------------------------------------
 335 void AliHFSystErr::InitDplustoKpipi2010PbPb020() {
 336   //
 337   // D+->Kpipi syst errors. Responsible: ??
 338   //  2010 PbPb sample, 0-20 CC
 339   //
 340
 341   AliInfo("Not yet implemented");
 342   return;
 343 }
 344
 345 //--------------------------------------------------------------------------
 346 void AliHFSystErr::InitDplustoKpipi2010PbPb4080() {
 347   //
 348   // D+->Kpipi syst errors. Responsible: ??
 349   //  2010 PbPb sample, 40-80 CC
 350   //
 351
 352   AliInfo("Not yet implemented");
 353   return;
 354 }
 355
 356 //--------------------------------------------------------------------------
 357 void AliHFSystErr::InitDstartoD0pi2010pp() {
 358   //
 359   // D*+->D0pi syst errors. Responsible: A. Grelli, Y. Wang
 360   //  2010 pp sample
 361   //
 362
 363  // Normalization
 364   fNorm = new TH1F("fNorm","fNorm",20,0,20);
 365   for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
 366
 367   // Branching ratio
 368   fBR = new TH1F("fBR","fBR",20,0,20);
 369   for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.015); // 1.5% PDG2010
 370
 371   // Tracking efficiency
 372   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
 373   fTrackingEff->SetBinContent(1,0.12);
 374   fTrackingEff->SetBinContent(2,0.08);
 375   fTrackingEff->SetBinContent(3,0.05);
 376   for(Int_t i=4;i<=20;i++) fTrackingEff->SetBinContent(i,0.03); // 3% (1% per track)
 377
 378
 379   // Raw yield extraction
 380   fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
 381   fRawYield->SetBinContent(1,1);
 382   fRawYield->SetBinContent(2,0.12);
 383   fRawYield->SetBinContent(3,0.09);
 384   fRawYield->SetBinContent(4,0.08);
 385   fRawYield->SetBinContent(5,0.06);
 386   for(Int_t i=5;i<=20;i++) fRawYield->SetBinContent(i,0.04);  //4%
 387
 388   // Cuts efficiency (from cuts variation)
 389   fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
 390   for(Int_t i=1;i<=20;i++) fCutsEff->SetBinContent(i,0.10); // 10%
 391
 392   // PID efficiency (from PID/noPID)
 393   fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
 394   for(Int_t i=1;i<=20;i++) fPIDEff->SetBinContent(i,0.04); // 3%
 395   fPIDEff->SetBinContent(1,1); // 100%
 396   fPIDEff->SetBinContent(2,1); // 100%
 397   fPIDEff->SetBinContent(3,0.05); // 5%
 398   fPIDEff->SetBinContent(4,0.05); // 5%
 399   fPIDEff->SetBinContent(5,0.05); // 5%
 400
 401
 402   // MC dN/dpt  (copied from D0 : will update later)
 403   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
 404   for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,(Float_t)i*0.006);
 405
 406   return;
 407 }
 408 //--------------------------------------------------------------------------
 409 void AliHFSystErr::InitDstartoD0pi2010PbPb020() {
 410   //
 411   // D*+->D0pi syst errors. Responsible: ??
 412   //  2010 PbPb sample, 0-20 CC
 413   //
 414
 415   AliInfo("Not yet implemented");
 416   return;
 417 }
 418 //--------------------------------------------------------------------------
 419 void AliHFSystErr::InitDstartoD0pi2010PbPb4080() {
 420   //
 421   // D*+->D0pi syst errors. Responsible: ??
 422   //  2010 PbPb sample, 40-80 CC
 423   //
 424
 425   AliInfo("Not yet implemented");
 426   return;
 427 }
 428 //--------------------------------------------------------------------------
 429 Double_t AliHFSystErr::GetCutsEffErr(Double_t pt) const {
 430   //
 431   // Get error
 432   //
 433
 434   Int_t bin=fCutsEff->FindBin(pt);
 435
 436   return fCutsEff->GetBinContent(bin);
 437 }
 438 //--------------------------------------------------------------------------
 439 Double_t AliHFSystErr::GetMCPtShapeErr(Double_t pt) const {
 440   //
 441   // Get error
 442   //
 443
 444   Int_t bin=fMCPtShape->FindBin(pt);
 445
 446   return fMCPtShape->GetBinContent(bin);
 447 }
 448 //--------------------------------------------------------------------------
 449 Double_t AliHFSystErr::GetSeleEffErr(Double_t pt) const {
 450   //
 451   // Get error
 452   //
 453
 454   Double_t err=GetCutsEffErr(pt)*GetCutsEffErr(pt)+GetMCPtShapeErr(pt)*GetMCPtShapeErr(pt);
 455
 456   return TMath::Sqrt(err);
 457 }
 458 //--------------------------------------------------------------------------
 459 Double_t AliHFSystErr::GetPIDEffErr(Double_t pt) const {
 460   //
 461   // Get error
 462   //
 463
 464   Int_t bin=fPIDEff->FindBin(pt);
 465
 466   return fPIDEff->GetBinContent(bin);
 467 }
 468 //--------------------------------------------------------------------------
 469 Double_t AliHFSystErr::GetTrackingEffErr(Double_t pt) const {
 470   //
 471   // Get error
 472   //
 473
 474   Int_t bin=fTrackingEff->FindBin(pt);
 475
 476   return fTrackingEff->GetBinContent(bin);
 477 }
 478 //--------------------------------------------------------------------------
 479 Double_t AliHFSystErr::GetRawYieldErr(Double_t pt) const {
 480   //
 481   // Get error
 482   //
 483
 484   Int_t bin=fRawYield->FindBin(pt);
 485
 486   return fRawYield->GetBinContent(bin);
 487 }
 488 //--------------------------------------------------------------------------
 489 Double_t AliHFSystErr::GetPartAntipartErr(Double_t pt) const {
 490   //
 491   // Get error
 492   //
 493
 494   Int_t bin=fPartAntipart->FindBin(pt);
 495
 496   return fPartAntipart->GetBinContent(bin);
 497 }
 498 //--------------------------------------------------------------------------
 499 Double_t AliHFSystErr::GetTotalSystErr(Double_t pt,Double_t feeddownErr) const {
 500   //
 501   // Get total syst error (except norm. error)
 502   //
 503
 504   Double_t err=0.;
 505
 506   if(fRawYield) err += GetRawYieldErr(pt)*GetRawYieldErr(pt);
 507   if(fTrackingEff) err += GetTrackingEffErr(pt)*GetTrackingEffErr(pt);
 508   if(fBR) err += GetBRErr()*GetBRErr();
 509   if(fCutsEff) err += GetCutsEffErr(pt)*GetCutsEffErr(pt);
 510   if(fPIDEff) err += GetPIDEffErr(pt)*GetPIDEffErr(pt);
 511   if(fMCPtShape) err += GetMCPtShapeErr(pt)*GetMCPtShapeErr(pt);
 512   if(fPartAntipart) err += GetPartAntipartErr(pt)*GetPartAntipartErr(pt);
 513
 514   err += feeddownErr*feeddownErr;
 515
 516   return TMath::Sqrt(err);
 517 }
 518 //---------------------------------------------------------------------------
 519 void AliHFSystErr::DrawErrors(TGraphAsymmErrors *grErrFeeddown) const {
 520   //
 521   // Draw errors
 522   //
 523   gStyle->SetOptStat(0);
 524
 525   TCanvas *cSystErr = new TCanvas("cSystErr","Systematic Errors",300,80,640,500);
 526   cSystErr->Range(0.20,-0.5,18.4,0.34);
 527   cSystErr->SetRightMargin(0.318);
 528   cSystErr->SetFillColor(0);
 529
 530   TH2F *hFrame = new TH2F("hFrame","Systematic errors; p_{t} [GeV/c]; Relative Error",20,0,20,100,-1,+1);
 531   hFrame->SetAxisRange(2.,11.9,"X");
 532   hFrame->SetAxisRange(-0.5,0.5,"Y");
 533   hFrame->Draw();
 534
 535   TLegend *leg = new TLegend(0.69,0.44,0.98,0.86,NULL,"brNDC");
 536   leg->SetTextSize(0.03601695);
 537   leg->SetFillStyle(0);
 538   leg->SetBorderSize(0);
 539
 540   TH1F *hTotErr=new TH1F("hTotErr","",20,0,20);
 541   Int_t nbins = fNorm->GetNbinsX();
 542   TGraphAsymmErrors *gTotErr = new TGraphAsymmErrors(nbins);
 543   for(Int_t i=1;i<=20;i++) {
 544     Double_t pt = hTotErr->GetBinCenter(i);
 545     Double_t ptwidth = hTotErr->GetBinWidth(i);
 546
 547     if(grErrFeeddown) {
 548       Double_t x=0., y=0., errxl=0., errxh=0., erryl=0., erryh=0.;
 549       Double_t toterryl=0., toterryh=0.;
 550       for(Int_t j=0; j<grErrFeeddown->GetN(); j++) {
 551         grErrFeeddown->GetPoint(j,x,y);
 552         errxh = grErrFeeddown->GetErrorXhigh(j);
 553         errxl = grErrFeeddown->GetErrorXlow(j);
 554         if ( ( (x-errxl) <= pt) && ( (x+errxl) >= pt) ) {
 555           erryh = grErrFeeddown->GetErrorYhigh(j);
 556           erryl = grErrFeeddown->GetErrorYlow(j);
 557         }
 558       }
 559       if (erryl>=1e-3) toterryl = GetTotalSystErr(pt,erryl);
 560       else toterryl = GetTotalSystErr(pt);
 561       if (erryh>=1e-3) toterryh = GetTotalSystErr(pt,erryh);
 562       else toterryh = GetTotalSystErr(pt);
 563
 564       hTotErr->SetBinContent(i,toterryh);
 565       gTotErr->SetPoint(i,pt,0.);
 566       gTotErr->SetPointError(i,ptwidth/2.,ptwidth/2.,toterryl,toterryh); // i, exl, exh, eyl, eyh
 567     }
 568     else {
 569       hTotErr->SetBinContent(i,GetTotalSystErr(pt));
 570       gTotErr->SetPoint(i,pt,0.);
 571       gTotErr->SetPointError(i,ptwidth/2.,ptwidth/2.,GetTotalSystErr(pt),GetTotalSystErr(pt)); // i, exl, exh, eyl, eyh
 572     }
 573
 574   }
 575   gTotErr->SetLineColor(kBlack);
 576   gTotErr->SetFillColor(kRed);
 577   gTotErr->SetFillStyle(3002);
 578   gTotErr->Draw("2");
 579   leg->AddEntry(gTotErr,"Total (excl. norm.)","f");
 580 //   hTotErr->SetLineColor(1);
 581 //   hTotErr->SetLineWidth(3);
 582 //   hTotErr->Draw("same");
 583 //   leg->AddEntry(hTotErr,"Total (excl. norm.)","l");
 584
 585
 586   fNorm->SetFillColor(1);
 587   fNorm->SetFillStyle(3002);
 588   //fNorm->Draw("same");
 589   //TH1F *hNormRefl = ReflectHisto(fNorm);
 590   //hNormRefl->Draw("same");
 591   leg->AddEntry(fNorm,"Normalization (10%)","");
 592
 593   if(grErrFeeddown) {
 594     grErrFeeddown->SetFillColor(kTeal-8);
 595     grErrFeeddown->SetFillStyle(3001);
 596     grErrFeeddown->Draw("2");
 597     leg->AddEntry(grErrFeeddown,"Feed-down from B","f");
 598   }
 599   if(fTrackingEff) {
 600     fTrackingEff->SetFillColor(4);
 601     fTrackingEff->SetFillStyle(3006);
 602     fTrackingEff->Draw("same");
 603     TH1F *hTrackingEffRefl = ReflectHisto(fTrackingEff);
 604     hTrackingEffRefl->Draw("same");
 605     leg->AddEntry(fTrackingEff,"Tracking efficiency","f");
 606   }
 607   if(fBR) {
 608     fBR->SetFillColor(6);
 609     fBR->SetFillStyle(3005);
 610     //fBR->SetFillStyle(3020);
 611     fBR->Draw("same");
 612     TH1F *hBRRefl = ReflectHisto(fBR);
 613     hBRRefl->Draw("same");
 614     leg->AddEntry(fBR,"Branching ratio","f");
 615   }
 616   if(fRawYield) {
 617     Int_t ci;   // for color index setting
 618     ci = TColor::GetColor("#00cc00");
 619     fRawYield->SetLineColor(ci);
 620     //    fRawYield->SetLineColor(3);
 621     fRawYield->SetLineWidth(3);
 622     fRawYield->Draw("same");
 623     TH1F *hRawYieldRefl = ReflectHisto(fRawYield);
 624     hRawYieldRefl->Draw("same");
 625     leg->AddEntry(fRawYield,"Yield extraction","l");
 626   }
 627   if(fCutsEff) {
 628     fCutsEff->SetLineColor(4);
 629     fCutsEff->SetLineWidth(3);
 630     fCutsEff->Draw("same");
 631     TH1F *hCutsEffRefl = ReflectHisto(fCutsEff);
 632     hCutsEffRefl->Draw("same");
 633     leg->AddEntry(fCutsEff,"Cuts efficiency","l");
 634   }
 635   if(fPIDEff) {
 636     fPIDEff->SetLineColor(7);
 637     fPIDEff->SetLineWidth(3);
 638     fPIDEff->Draw("same");
 639     TH1F *hPIDEffRefl = ReflectHisto(fPIDEff);
 640     hPIDEffRefl->Draw("same");
 641     leg->AddEntry(fPIDEff,"PID efficiency","l");
 642   }
 643   if(fMCPtShape) {
 644     Int_t ci = TColor::GetColor("#9933ff");
 645     fMCPtShape->SetLineColor(ci);
 646     //    fMCPtShape->SetLineColor(8);
 647     fMCPtShape->SetLineWidth(3);
 648     fMCPtShape->Draw("same");
 649     TH1F *hMCPtShapeRefl = ReflectHisto(fMCPtShape);
 650     hMCPtShapeRefl->Draw("same");
 651     leg->AddEntry(fMCPtShape,"MC p_{t} shape","l");
 652   }
 653   if(fPartAntipart) {
 654     Int_t ci = TColor::GetColor("#ff6600");
 655     fPartAntipart->SetLineColor(ci);
 656     //    fPartAntipart->SetLineColor(9);
 657     fPartAntipart->SetLineWidth(3);
 658     fPartAntipart->Draw("same");
 659     TH1F *hPartAntipartRefl = ReflectHisto(fPartAntipart);
 660     hPartAntipartRefl->Draw("same");
 661     leg->AddEntry(fPartAntipart,"D = #bar{D}","l");
 662   }
 663
 664
 665   leg->Draw();
 666
 667   cSystErr->SaveAs("RelativeSystematics.eps");
 668
 669   return;
 670 }
 671 //-------------------------------------------------------------------------
 672 TH1F* AliHFSystErr::ReflectHisto(TH1F *hin) const {
 673   //
 674   // Clones and reflects histogram
 675   //
 676   TH1F *hout=(TH1F*)hin->Clone("hout");
 677   hout->Scale(-1.);
 678
 679   return hout;
 680 }