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 fIsLowEnergy(false)
  63 {
  64   //
  65   // Default Constructor
  66   //
  67 }
  68
  69 //--------------------------------------------------------------------------
  70 AliHFSystErr::~AliHFSystErr() {
  71   //
  72   // Default Destructor
  73   //
  74   /*
  75
  76   if(fNorm)         { delete fNorm; fNorm=0; }
  77   if(fRawYield)     { delete fRawYield; fRawYield=0; }
  78   if(fTrackingEff)  { delete fTrackingEff; fTrackingEff=0; }
  79   if(fBR)           { delete fBR; fBR=0; }
  80   if(fCutsEff)      { delete fCutsEff; fCutsEff=0; }
  81   if(fPIDEff)       { delete fPIDEff; fPIDEff=0; }
  82   if(fMCPtShape)    { delete fMCPtShape; fMCPtShape=0; }
  83   if(fPartAntipart) { delete fPartAntipart; fPartAntipart=0; }
  84   */
  85 }
  86
  87 //--------------------------------------------------------------------------
  88 void AliHFSystErr::Init(Int_t decay){
  89   //
  90   // Variables/histos initialization
  91   //
  92
  93   if (fRunNumber!=10 && fIsLowEnergy==false) {
  94     AliError("Only settings for 2010 and the low energy runs are implemented so far");
  95   }
  96   if (fCentralityClass!="020" && fCentralityClass!="4080" && fCentralityClass!="0100"){
  97     AliError("Only settings for MB2010 are implemented so far");
  98   }
  99
 100   switch(decay) {
 101   case 1: // D0->Kpi
 102     if (fCollisionType==0) {
 103       if (fIsLowEnergy) InitD0toKpi2010ppLowEn();
 104       else InitD0toKpi2010pp();
 105     } else if (fCollisionType==1) {
 106       if (fCentralityClass=="020") InitD0toKpi2010PbPb020();
 107       else if (fCentralityClass=="4080") InitD0toKpi2010PbPb4080();
 108       else AliError("Not yet implemented");
 109     }
 110     //    else if (fCollisionType==2) InitD0toKpi2010ppLowEn();
 111     break;
 112   case 2: // D+->Kpipi
 113     if (fCollisionType==0) {
 114       if (fIsLowEnergy) InitDplustoKpipi2010ppLowEn();
 115       else InitDplustoKpipi2010pp();
 116     } else if (fCollisionType==1) {
 117       if (fCentralityClass=="020") InitDplustoKpipi2010PbPb020();
 118       else if (fCentralityClass=="4080") InitDplustoKpipi2010PbPb4080();
 119       else AliError("Not yet implemented");
 120     }
 121     break;
 122   case 3: // D*->D0pi
 123     if (fCollisionType==0) {
 124       if(fIsLowEnergy)  InitDstartoD0pi2010ppLowEn();
 125       else InitDstartoD0pi2010pp();
 126     }else if (fCollisionType==1) {
 127       if (fCentralityClass=="020")  InitDstartoD0pi2010PbPb020();
 128       if (fCentralityClass=="2040") InitDstartoD0pi2010PbPb2040();
 129       if (fCentralityClass=="4080") InitDstartoD0pi2010PbPb4080();
 130       if (fCentralityClass!="4080" && fCentralityClass!="2040" && fCentralityClass!="020")  AliError("Not yet implemented");
 131     }
 132     break;
 133   default:
 134     printf("Invalid decay type: %d\n",decay);
 135     break;
 136   }
 137
 138 }
 139
 140 //--------------------------------------------------------------------------
 141 void AliHFSystErr::InitD0toKpi2010pp() {
 142   //
 143   // D0->Kpi syst errors. Responsible: A. Rossi
 144   //   2010 pp sample
 145   //
 146
 147   // Normalization
 148   fNorm = new TH1F("fNorm","fNorm",24,0,24);
 149   for(Int_t i=1;i<=24;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
 150
 151   // Branching ratio
 152   fBR = new TH1F("fBR","fBR",24,0,24);
 153   for(Int_t i=1;i<=24;i++) fBR->SetBinContent(i,0.012); // 1.2% PDG2010
 154
 155   // Tracking efficiency
 156   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",24,0,24);
 157   for(Int_t i=1;i<=24;i++) fTrackingEff->SetBinContent(i,0.02); // 2% (1% per track)
 158
 159   // Raw yield extraction
 160   fRawYield = new TH1F("fRawYield","fRawYield",24,0,24);
 161   fRawYield->SetBinContent(1,1);
 162   fRawYield->SetBinContent(2,0.2);
 163   for(Int_t i=3;i<=16;i++) fRawYield->SetBinContent(i,0.05);
 164   for(Int_t i=13;i<=16;i++) fRawYield->SetBinContent(i,0.10);
 165   for(Int_t i=17;i<=24;i++) fRawYield->SetBinContent(i,1);
 166
 167   // Cuts efficiency (from cuts variation)
 168   fCutsEff = new TH1F("fCutsEff","fCutsEff",24,0,24);
 169   for(Int_t i=1;i<=24;i++) fCutsEff->SetBinContent(i,0.10); // 10%
 170
 171   // PID efficiency (from PID/noPID)
 172   fPIDEff = new TH1F("fPIDEff","fPIDEff",24,0,24);
 173   for(Int_t i=1;i<=24;i++) fPIDEff->SetBinContent(i,0.03); // 3%
 174   fPIDEff->SetBinContent(2,0.10); // 10%
 175
 176   // MC dN/dpt
 177   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",24,0,24);
 178   for(Int_t i=1;i<=24;i++) fMCPtShape->SetBinContent(i,0);
 179   fMCPtShape->SetBinContent(1,0.03);
 180   fMCPtShape->SetBinContent(2,0.03);
 181
 182   // particle-antiparticle
 183   //  fPartAntipart = new TH1F("fPartAntipart","fPartAntipart",24,0,24);
 184   //  fPartAntipart->SetBinContent(1,1);
 185   //  for(Int_t i=2;i<=24;i++) fPartAntipart->SetBinContent(i,0.05);
 186
 187   return;
 188 }
 189 //--------------------------------------------------------------------------
 190 void AliHFSystErr::InitD0toKpi2010PbPb020() {
 191   //
 192   // D0->Kpi syst errors. Responsible: A. Rossi
 193   //   2010 PbPb sample, 0-20 CC
 194   //
 195
 196   // Normalization
 197   fNorm = new TH1F("fNorm","fNorm",20,0,20);
 198   for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
 199
 200   // Branching ratio
 201   fBR = new TH1F("fBR","fBR",20,0,20);
 202   for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.012); // 1.2% PDG2010
 203
 204   // Tracking efficiency
 205   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
 206   for(Int_t i=2;i<=12;i++) fTrackingEff->SetBinContent(i,0.05+0.005*(Float_t)i);
 207
 208
 209   // Raw yield extraction
 210   fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
 211   fRawYield->SetBinContent(1,0);
 212   fRawYield->SetBinContent(2,0);
 213   fRawYield->SetBinContent(3,0.05);
 214   fRawYield->SetBinContent(4,0.05);
 215   fRawYield->SetBinContent(5,0.10);
 216   fRawYield->SetBinContent(6,0.10);
 217   for(Int_t i=7;i<=12;i++) fRawYield->SetBinContent(i,0.05);
 218  // for(Int_t i=13;i<=20;i++) fRawYield->SetBinContent(i,0);
 219
 220   // Cuts efficiency (from cuts variation)
 221   fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
 222   fCutsEff->SetBinContent(1,0.);
 223   fCutsEff->SetBinContent(2,0.);
 224   fCutsEff->SetBinContent(3,0.18);
 225   fCutsEff->SetBinContent(4,0.18);
 226   fCutsEff->SetBinContent(5,0.18);
 227   fCutsEff->SetBinContent(6,0.15);
 228   fCutsEff->SetBinContent(7,0.15);
 229   fCutsEff->SetBinContent(8,0.15);
 230 fCutsEff->SetBinContent(9,0.15);
 231  fCutsEff->SetBinContent(10,0.15);
 232   fCutsEff->SetBinContent(11,0.15);
 233   fCutsEff->SetBinContent(12,0.15);
 234  // for(Int_t i=13;i<=20;i++) fCutsEff->SetBinContent(i,0.);
 235
 236   // PID efficiency (from PID/noPID)
 237   fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
 238   for(Int_t i=3;i<=12;i++) fPIDEff->SetBinContent(i,0.05);
 239   fPIDEff->SetBinContent(5,0.10);
 240   fPIDEff->SetBinContent(6,0.10);
 241
 242   // MC dN/dpt
 243   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
 244   for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,0);
 245   fMCPtShape->SetBinContent(1,0.03);
 246   fMCPtShape->SetBinContent(2,0.03);
 247
 248
 249   // particle-antiparticle
 250   fPartAntipart = new TH1F("fPartAntipart","fPartAntipart",20,0,20);
 251   for(Int_t i=3;i<=12;i++) fPartAntipart->SetBinContent(i,0.05);
 252   fPartAntipart->SetBinContent(3,0.10);
 253   fPartAntipart->SetBinContent(4,0.10);
 254   fPartAntipart->SetBinContent(7,0.10);
 255   fPartAntipart->SetBinContent(8,0.10);
 256
 257   return;
 258 }
 259 //--------------------------------------------------------------------------
 260 void AliHFSystErr::InitD0toKpi2010PbPb4080() {
 261   //
 262   // D0->Kpi syst errors. Responsible: A. Rossi
 263   //   2010 PbPb sample, 40-80 CC
 264   //
 265
 266   // Normalization
 267   fNorm = new TH1F("fNorm","fNorm",20,0,20);
 268   for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
 269
 270   // Branching ratio
 271   fBR = new TH1F("fBR","fBR",20,0,20);
 272   for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.012); // 1.2% PDG2010
 273
 274   // Tracking efficiency
 275   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
 276   for(Int_t i=2;i<=12;i++) fTrackingEff->SetBinContent(i,0.5*(0.05+0.005*(Float_t)i));
 277
 278
 279   // Raw yield extraction
 280   fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
 281   fRawYield->SetBinContent(1,0);
 282   fRawYield->SetBinContent(2,0);
 283   for(Int_t i=3;i<=12;i++) fRawYield->SetBinContent(i,0.05);
 284   //for(Int_t i=13;i<=20;i++) fRawYield->SetBinContent(i,0);
 285
 286   // Cuts efficiency (from cuts variation)
 287   fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
 288   fCutsEff->SetBinContent(1,0.);
 289   fCutsEff->SetBinContent(2,0.);
 290   fCutsEff->SetBinContent(3,0.15);
 291   fCutsEff->SetBinContent(4,0.15);
 292   fCutsEff->SetBinContent(5,0.15);
 293   fCutsEff->SetBinContent(6,0.15);
 294   fCutsEff->SetBinContent(7,0.15);
 295   fCutsEff->SetBinContent(8,0.15);
 296   fCutsEff->SetBinContent(9,0.15);
 297   fCutsEff->SetBinContent(10,0.15);
 298   fCutsEff->SetBinContent(11,0.15);
 299   fCutsEff->SetBinContent(12,0.15);
 300  // for(Int_t i=13;i<=20;i++) fCutsEff->SetBinContent(i,0.);
 301
 302   // PID efficiency (from PID/noPID)
 303   fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
 304   for(Int_t i=3;i<=12;i++) fPIDEff->SetBinContent(i,0.05);
 305
 306   // MC dN/dpt
 307   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
 308   for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,0);
 309   fMCPtShape->SetBinContent(1,0.03);
 310   fMCPtShape->SetBinContent(2,0.03);
 311
 312   // particle-antiparticle
 313   fPartAntipart = new TH1F("fPartAntipart","fPartAntipart",20,0,20);
 314   for(Int_t i=3;i<=12;i++) fPartAntipart->SetBinContent(i,0.05);
 315
 316   return;
 317 }
 318
 319 //--------------------------------------------------------------------------
 320 void AliHFSystErr::InitD0toKpi2010ppLowEn() {
 321   //
 322   // D0->Kpi syst errors. Low energy run
 323   //   2011 2.76 TeV pp sample
 324   //
 325   AliInfo(" Settings for D0 --> K pi, p-p collisions at 2.76 TeV");
 326
 327   // Normalization
 328   fNorm = new TH1F("fNorm","fNorm",20,0,20);
 329   for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
 330
 331   // Branching ratio
 332   fBR = new TH1F("fBR","fBR",20,0,20);
 333   for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.012); // 1.2% PDG2010
 334
 335   // Tracking efficiency
 336   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
 337   for(Int_t i=1;i<=20;i++) fTrackingEff->SetBinContent(i,0.10); //10% (5% per track)
 338
 339   // Raw yield extraction
 340   fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
 341   fRawYield->SetBinContent(1,1);
 342   fRawYield->SetBinContent(2,1);
 343   fRawYield->SetBinContent(3,0.2);
 344   fRawYield->SetBinContent(4,0.2);
 345   fRawYield->SetBinContent(5,0.1);
 346   fRawYield->SetBinContent(6,0.1);
 347   fRawYield->SetBinContent(7,0.2);
 348   fRawYield->SetBinContent(8,0.2);
 349   for(Int_t i=9;i<=20;i++) fRawYield->SetBinContent(i,0.065);
 350
 351   // Cuts efficiency (from cuts variation)
 352   fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
 353   for(Int_t i=1;i<=20;i++) fCutsEff->SetBinContent(i,0.18); // 10%
 354
 355   // PID efficiency (from PID/noPID)
 356   fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
 357   for(Int_t i=1;i<=20;i++) fPIDEff->SetBinContent(i,0.20); // 10%
 358   fPIDEff->SetBinContent(3,0.25); // 10%
 359   fPIDEff->SetBinContent(4,0.25); // 10%
 360   fPIDEff->SetBinContent(7,0.25); // 10%
 361   fPIDEff->SetBinContent(8,0.25); // 10%
 362
 363   // MC dN/dpt
 364   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
 365   for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,0);
 366   fMCPtShape->SetBinContent(1,0.03);
 367   fMCPtShape->SetBinContent(2,0.03);
 368
 369   // particle-antiparticle
 370   fPartAntipart = new TH1F("fPartAntipart","fPartAntipart",20,0,20);
 371   fPartAntipart->SetBinContent(1,1);
 372   fPartAntipart->SetBinContent(2,1);
 373   for(Int_t i=3;i<=6;i++) fPartAntipart->SetBinContent(i,0.08);
 374
 375   return;
 376 }
 377
 378 //--------------------------------------------------------------------------
 379 void AliHFSystErr::InitDplustoKpipi2010pp() {
 380   //
 381   // D+->Kpipi syst errors. Responsible: R. Bala
 382   //  2010 pp sample
 383   //
 384
 385  // Normalization
 386   fNorm = new TH1F("fNorm","fNorm",20,0,20);
 387   for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
 388
 389   // Branching ratio
 390   fBR = new TH1F("fBR","fBR",20,0,20);
 391   for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.04); // 4% PDG2010
 392
 393   // Tracking efficiency
 394   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
 395   for(Int_t i=1;i<=20;i++) fTrackingEff->SetBinContent(i,0.03); // 3% (1% per track)
 396
 397
 398   // Raw yield extraction
 399   fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
 400   fRawYield->SetBinContent(1,1);
 401   fRawYield->SetBinContent(2,1);
 402   fRawYield->SetBinContent(3,0.20);
 403   for(Int_t i=4;i<=20;i++) fRawYield->SetBinContent(i,0.055);  //5 to 10%
 404
 405   // Cuts efficiency (from cuts variation)
 406   fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
 407   for(Int_t i=1;i<=20;i++) fCutsEff->SetBinContent(i,0.10); // 10%
 408
 409   // PID efficiency (from PID/noPID)
 410   fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
 411   for(Int_t i=1;i<=20;i++) fPIDEff->SetBinContent(i,0.05); // 5%
 412   fPIDEff->SetBinContent(3,0.13); // 13%
 413
 414
 415   // MC dN/dpt  (copied from D0 : will update later)
 416   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
 417   for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,0);
 418   fMCPtShape->SetBinContent(1,0.03);
 419   fMCPtShape->SetBinContent(2,0.03);
 420
 421
 422   // particle-antiparticle
 423   /*
 424   fPartAntipart = new TH1F("fPartAntipart","fPartAntipart",20,0,20);
 425   fPartAntipart->SetBinContent(1,1);
 426   fPartAntipart->SetBinContent(2,1);
 427   fPartAntipart->SetBinContent(3,0.12);
 428   for(Int_t i=4;i<=20;i++) fPartAntipart->SetBinContent(i,0.05);   //5 to 12%
 429   */
 430   return;
 431 }
 432 //--------------------------------------------------------------------------
 433 void AliHFSystErr::InitDplustoKpipi2010PbPb020() {
 434   //
 435   // D+->Kpipi syst errors. Responsible: ??
 436   //  2010 PbPb sample, 0-20 CC
 437   //
 438
 439  // Normalization
 440   fNorm = new TH1F("fNorm","fNorm",20,0,20);
 441   for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
 442
 443   // Branching ratio
 444   fBR = new TH1F("fBR","fBR",20,0,20);
 445   for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.04); // 4% PDG2010
 446
 447   // Tracking efficiency
 448   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
 449   for(Int_t i=1;i<=20;i++) fTrackingEff->SetBinContent(i,i,0.075+0.005*(Float_t)i);
 450
 451
 452   // Raw yield extraction
 453   fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
 454   for(Int_t i=1;i<=20;i++) fRawYield->SetBinContent(i,1);  //5 to 10%
 455 //   for(Int_t i=6;i<=7;i++) fRawYield->SetBinContent(i,0.32);
 456 //   for(Int_t i=7;i<=9;i++) fRawYield->SetBinContent(i,0.32);
 457 //   for(Int_t i=9;i<=12;i++) fRawYield->SetBinContent(i,0.05);//31./156.);
 458   fRawYield->SetBinContent(6,0.12);
 459   fRawYield->SetBinContent(7,0.17);
 460   fRawYield->SetBinContent(8,0.17);
 461   fRawYield->SetBinContent(9,0.08);
 462   fRawYield->SetBinContent(10,0.08);
 463   fRawYield->SetBinContent(11,0.08);
 464   fRawYield->SetBinContent(12,0.08);
 465
 466   // Cuts efficiency (from cuts variation)
 467   fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
 468   for(Int_t i=1;i<=20;i++) fCutsEff->SetBinContent(i,0.15); // 10%
 469
 470   // PID efficiency (from PID/noPID)
 471   fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
 472   for(Int_t i=1;i<=20;i++) fPIDEff->SetBinContent(i,0.05); // 5%
 473   fPIDEff->SetBinContent(3,0.13); // 13%
 474
 475
 476   // MC dN/dpt  (copied from D0 : will update later)
 477   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
 478   for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,0.);
 479   fMCPtShape->SetBinContent(1,0.03);
 480   fMCPtShape->SetBinContent(2,0.03);
 481
 482
 483   // particle-antiparticle
 484   /*
 485   fPartAntipart = new TH1F("fPartAntipart","fPartAntipart",20,0,20);
 486   fPartAntipart->SetBinContent(1,1);
 487   fPartAntipart->SetBinContent(2,1);
 488   fPartAntipart->SetBinContent(3,0.12);
 489   for(Int_t i=4;i<=20;i++) fPartAntipart->SetBinContent(i,0.05);   //5 to 12%
 490   */
 491   return;
 492 }
 493
 494 //--------------------------------------------------------------------------
 495 void AliHFSystErr::InitDplustoKpipi2010PbPb4080() {
 496   //
 497   // D+->Kpipi syst errors. Responsible: ??
 498   //  2010 PbPb sample, 40-80 CC
 499   //
 500
 501
 502  // Normalization
 503   fNorm = new TH1F("fNorm","fNorm",20,0,20);
 504   for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
 505
 506   // Branching ratio
 507   fBR = new TH1F("fBR","fBR",20,0,20);
 508   for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.04); // 4% PDG2010
 509
 510   // Tracking efficiency
 511   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
 512   for(Int_t i=1;i<=20;i++) fTrackingEff->SetBinContent(i,0.03); // 3% (1% per track)
 513
 514
 515   // Raw yield extraction
 516   fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
 517   fRawYield->SetBinContent(1,1);
 518   fRawYield->SetBinContent(2,1);
 519   fRawYield->SetBinContent(3,1);
 520   fRawYield->SetBinContent(4,1.);
 521   fRawYield->SetBinContent(5,4./23.);
 522   fRawYield->SetBinContent(6,4./23.);
 523   fRawYield->SetBinContent(7,5./20.);
 524   fRawYield->SetBinContent(8,5./20.);
 525   for(Int_t i=9;i<=12;i++) fRawYield->SetBinContent(i,1./16.);
 526   for(Int_t i=13;i<=20;i++) fRawYield->SetBinContent(i,1);  //5 to 10%
 527
 528   // Cuts efficiency (from cuts variation)
 529   fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
 530   for(Int_t i=1;i<=20;i++) fCutsEff->SetBinContent(i,0.10); // 10%
 531
 532   // PID efficiency (from PID/noPID)
 533   fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
 534   for(Int_t i=1;i<=20;i++) fPIDEff->SetBinContent(i,0.05); // 5%
 535   fPIDEff->SetBinContent(3,0.13); // 13%
 536
 537
 538   // MC dN/dpt  (copied from D0 : will update later)
 539   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
 540   for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,0);
 541   fMCPtShape->SetBinContent(1,0.03);
 542   fMCPtShape->SetBinContent(2,0.03);
 543
 544
 545   // particle-antiparticle
 546   /*
 547   fPartAntipart = new TH1F("fPartAntipart","fPartAntipart",20,0,20);
 548   fPartAntipart->SetBinContent(1,1);
 549   fPartAntipart->SetBinContent(2,1);
 550   fPartAntipart->SetBinContent(3,0.12);
 551   for(Int_t i=4;i<=20;i++) fPartAntipart->SetBinContent(i,0.05);   //5 to 12%
 552   */
 553   return;
 554 }
 555
 556 //--------------------------------------------------------------------------
 557 void AliHFSystErr::InitDplustoKpipi2010ppLowEn() {
 558
 559   //
 560   // D+->Kpipi syst errors. Responsible: R. Bala
 561   //  2011 2.76 TeV pp sample
 562   //
 563   AliInfo(" Settings for D+ --> K pi pi p-p collisions at 2.76 TeV");
 564
 565   // Normalization
 566   fNorm = new TH1F("fNorm","fNorm",20,0,20);
 567   for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
 568
 569   // Branching ratio
 570   fBR = new TH1F("fBR","fBR",20,0,20);
 571   for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.04); // 4% PDG2010
 572
 573   // Tracking efficiency
 574   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
 575   for(Int_t i=1;i<=20;i++) fTrackingEff->SetBinContent(i,0.15); // 3% (1% per track)
 576
 577   // Raw yield extraction
 578   fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
 579   fRawYield->SetBinContent(1,1);
 580   fRawYield->SetBinContent(2,1);
 581   fRawYield->SetBinContent(3,0.20);
 582   fRawYield->SetBinContent(4,0.20);
 583   for(Int_t i=5;i<=20;i++) fRawYield->SetBinContent(i,0.055);  //5 to 10%
 584
 585   // Cuts efficiency (from cuts variation)
 586   fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
 587   for(Int_t i=1;i<=20;i++) fCutsEff->SetBinContent(i,0.10); // 10%
 588
 589   // PID efficiency (from PID/noPID)
 590   fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
 591   for(Int_t i=1;i<=20;i++) fPIDEff->SetBinContent(i,0.05); // 5%
 592   fPIDEff->SetBinContent(3,0.15); // 13%
 593   fPIDEff->SetBinContent(4,0.15); // 13%
 594
 595   // MC dN/dpt  (copied from D0 : will update later)
 596   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
 597   for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,0.);
 598   fMCPtShape->SetBinContent(1,0.03);
 599   fMCPtShape->SetBinContent(2,0.03);
 600
 601   return;
 602 }
 603
 604 //--------------------------------------------------------------------------
 605 void AliHFSystErr::InitDstartoD0pi2010pp() {
 606   //
 607   // D*+->D0pi syst errors. Responsible: A. Grelli, Y. Wang
 608   //  2010 pp sample
 609   //
 610
 611  // Normalization
 612   fNorm = new TH1F("fNorm","fNorm",24,0,24);
 613   for(Int_t i=1;i<=24;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
 614
 615   // Branching ratio
 616   fBR = new TH1F("fBR","fBR",24,0,24);
 617   for(Int_t i=1;i<=24;i++) fBR->SetBinContent(i,0.015); // 1.5% PDG2010
 618
 619   // Tracking efficiency
 620   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",24,0,24);
 621   fTrackingEff->SetBinContent(1,1.0);
 622   fTrackingEff->SetBinContent(2,0.12);
 623   fTrackingEff->SetBinContent(3,0.08);
 624   fTrackingEff->SetBinContent(3,0.05);
 625   for(Int_t i=4;i<=24;i++) fTrackingEff->SetBinContent(i,0.03); // 3% (1% per track)
 626
 627
 628   // Raw yield extraction
 629   fRawYield = new TH1F("fRawYield","fRawYield",24,0,24);
 630   fRawYield->SetBinContent(1,1.0);
 631   fRawYield->SetBinContent(2,0.09);
 632   fRawYield->SetBinContent(3,0.04);
 633   fRawYield->SetBinContent(4,0.02);
 634   fRawYield->SetBinContent(5,0.03);
 635   fRawYield->SetBinContent(6,0.06);
 636   fRawYield->SetBinContent(7,0.04);
 637   fRawYield->SetBinContent(8,0.11);
 638   for(Int_t i=5;i<=24;i++) fRawYield->SetBinContent(i,0.03);  //4%
 639   fRawYield->SetBinContent(13,0.09);
 640   fRawYield->SetBinContent(14,0.09);
 641   fRawYield->SetBinContent(15,0.09);
 642   fRawYield->SetBinContent(16,0.09);
 643   fRawYield->SetBinContent(17,0.24);
 644   fRawYield->SetBinContent(18,0.24);
 645   fRawYield->SetBinContent(19,0.24);
 646   fRawYield->SetBinContent(20,0.24);
 647   fRawYield->SetBinContent(21,0.24);
 648   fRawYield->SetBinContent(22,0.24);
 649   fRawYield->SetBinContent(23,0.24);
 650   fRawYield->SetBinContent(24,0.24);
 651
 652   // Cuts efficiency (from cuts variation)
 653   fCutsEff = new TH1F("fCutsEff","fCutsEff",24,0,24);
 654   fCutsEff->SetBinContent(2,0.22);
 655   for(Int_t i=3;i<=24;i++) fCutsEff->SetBinContent(i,0.10); // 10%
 656
 657   // PID efficiency (from PID/noPID)
 658   fPIDEff = new TH1F("fPIDEff","fPIDEff",24,0,24);
 659   for(Int_t i=1;i<=24;i++) fPIDEff->SetBinContent(i,0.04); // 3%
 660
 661
 662   // MC dN/dpt  (copied from D0 : will update later)
 663   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",24,0,24);
 664   for(Int_t i=1;i<=24;i++) fMCPtShape->SetBinContent(i,0);
 665   fMCPtShape->SetBinContent(1,0.03);
 666   fMCPtShape->SetBinContent(2,0.03);
 667
 668   return;
 669
 670
 671 }
 672 //--------------------------------------------------------------------------
 673 void AliHFSystErr::InitDstartoD0pi2010ppLowEn() {
 674
 675   //
 676   // D+->Kpipi syst errors. Responsible: A. Grelli
 677   //  2011 2.76 TeV pp sample
 678   //
 679   AliInfo(" Settings for D*+ --> D0 pi p-p collisions at 2.76 TeV");
 680
 681 // Normalization
 682   fNorm = new TH1F("fNorm","fNorm",20,0,20);
 683   for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
 684
 685   // Branching ratio
 686   fBR = new TH1F("fBR","fBR",20,0,20);
 687   for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.015); // 1.5% PDG2010
 688
 689   // Tracking efficiency
 690   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
 691   for(Int_t i=1;i<=20;i++) fTrackingEff->SetBinContent(i,0.15); //10% (to be checked!!)
 692
 693   // Raw yield extraction
 694   fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
 695   fRawYield->SetBinContent(1,1);
 696   fRawYield->SetBinContent(2,1);
 697   fRawYield->SetBinContent(3,0.2);
 698   fRawYield->SetBinContent(4,0.2);
 699   fRawYield->SetBinContent(5,0.08);
 700   fRawYield->SetBinContent(6,0.08);
 701   fRawYield->SetBinContent(7,0.1);
 702   fRawYield->SetBinContent(8,0.1);
 703   fRawYield->SetBinContent(9,0.2);
 704   fRawYield->SetBinContent(10,0.2);
 705   for(Int_t i=9;i<=20;i++) fRawYield->SetBinContent(i,0.065);
 706
 707   // Cuts efficiency (from cuts variation)
 708   fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
 709   for(Int_t i=1;i<=20;i++) fCutsEff->SetBinContent(i,0.18);
 710   fCutsEff->SetBinContent(3,0.2);
 711   fCutsEff->SetBinContent(4,0.2);
 712   fCutsEff->SetBinContent(5,0.1);
 713   fCutsEff->SetBinContent(6,0.1);
 714   fCutsEff->SetBinContent(7,0.12);
 715   fCutsEff->SetBinContent(8,0.12);
 716   fCutsEff->SetBinContent(9,0.2);
 717   fCutsEff->SetBinContent(10,0.2);
 718
 719   // PID efficiency (from PID/noPID)
 720   fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
 721   for(Int_t i=1;i<=20;i++) fPIDEff->SetBinContent(i,0.05); // 10%
 722
 723   // MC dN/dpt
 724   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
 725   for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,0.);
 726   fMCPtShape->SetBinContent(1,0.03);
 727   fMCPtShape->SetBinContent(2,0.03);
 728
 729
 730   return;
 731 }
 732 //------------------------------------------------------------------------
 733 void AliHFSystErr::InitDstartoD0pi2010PbPb020() {
 734   //
 735   // D*+->D0pi syst errors. Responsible: A. Grelli
 736   //  2010 PbPb sample, 0-20 CC
 737   //
 738
 739   AliInfo(" Settings for D*+ --> D0pi Pb-Pb collisions at 2.76 TeV - 0-20 centrality - DUMMY");
 740
 741  // Normalization
 742   fNorm = new TH1F("fNorm","fNorm",24,0,24);
 743   for(Int_t i=1;i<=24;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
 744
 745   // Branching ratio
 746   fBR = new TH1F("fBR","fBR",24,0,24);
 747   for(Int_t i=1;i<=24;i++) fBR->SetBinContent(i,0.015); // 1.5% PDG2010
 748
 749   // Tracking efficiency
 750   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",24,0,24);;
 751   for(Int_t i=1;i<=24;i++) fTrackingEff->SetBinContent(i,0.12); // 3% (1% per track)
 752
 753
 754   // Raw yield extraction
 755   fRawYield = new TH1F("fRawYield","fRawYield",24,0,24);
 756   for(Int_t i=1;i<=24;i++) fRawYield->SetBinContent(i,0.15);  //4%
 757
 758   // Cuts efficiency (from cuts variation)
 759   fCutsEff = new TH1F("fCutsEff","fCutsEff",24,0,24);
 760   for(Int_t i=1;i<=24;i++) fCutsEff->SetBinContent(i,0.10); // 10%
 761
 762   // PID efficiency (from PID/noPID)
 763   fPIDEff = new TH1F("fPIDEff","fPIDEff",24,0,24);
 764   for(Int_t i=1;i<=24;i++) fPIDEff->SetBinContent(i,0.04); // 3%
 765
 766
 767   // MC dN/dpt  (copied from D0 : will update later)
 768   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",24,0,24);
 769   for(Int_t i=1;i<=24;i++) fMCPtShape->SetBinContent(i,0.);
 770   fMCPtShape->SetBinContent(1,0.03);
 771   fMCPtShape->SetBinContent(2,0.03);
 772
 773   return;
 774
 775 }
 776 //-------------------------------------------------------------------------
 777 void AliHFSystErr::InitDstartoD0pi2010PbPb2040() {
 778   //
 779   // D*+->D0pi syst errors. Responsible: A. Grelli
 780   //  2010 PbPb sample, 20-40 CC
 781   //
 782
 783   AliInfo(" Settings for D*+ --> D0pi Pb-Pb collisions at 2.76 TeV - 20-40 centrality - DUMMY");
 784
 785  // Normalization
 786   fNorm = new TH1F("fNorm","fNorm",24,0,24);
 787   for(Int_t i=1;i<=24;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
 788
 789   // Branching ratio
 790   fBR = new TH1F("fBR","fBR",24,0,24);
 791   for(Int_t i=1;i<=24;i++) fBR->SetBinContent(i,0.015); // 1.5% PDG2010
 792
 793   // Tracking efficiency
 794   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",24,0,24);;
 795   for(Int_t i=1;i<=24;i++) fTrackingEff->SetBinContent(i,0.12); // 3% (1% per track)
 796
 797
 798   // Raw yield extraction
 799   fRawYield = new TH1F("fRawYield","fRawYield",24,0,24);
 800   for(Int_t i=1;i<=24;i++) fRawYield->SetBinContent(i,0.15);  //4%
 801
 802   // Cuts efficiency (from cuts variation)
 803   fCutsEff = new TH1F("fCutsEff","fCutsEff",24,0,24);
 804   for(Int_t i=1;i<=24;i++) fCutsEff->SetBinContent(i,0.10); // 10%
 805
 806   // PID efficiency (from PID/noPID)
 807   fPIDEff = new TH1F("fPIDEff","fPIDEff",24,0,24);
 808   for(Int_t i=1;i<=24;i++) fPIDEff->SetBinContent(i,0.04); // 3%
 809
 810
 811   // MC dN/dpt  (copied from D0 : will update later)
 812   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",24,0,24);
 813   for(Int_t i=1;i<=24;i++) fMCPtShape->SetBinContent(i,0.);
 814   fMCPtShape->SetBinContent(1,0.03);
 815   fMCPtShape->SetBinContent(2,0.03);
 816
 817   return;
 818
 819 }
 820 //--------------------------------------------------------------------------
 821 void AliHFSystErr::InitDstartoD0pi2010PbPb4080() {
 822   //
 823   // D*+->D0pi syst errors. Responsible: A. Grelli
 824   //  2010 PbPb sample, 40-80 CC
 825   //
 826
 827   AliInfo(" Settings for D*+ --> D0pi Pb-Pb collisions at 2.76 TeV - 40-80 centrality - DUMMY");
 828
 829  // Normalization
 830   fNorm = new TH1F("fNorm","fNorm",24,0,24);
 831   for(Int_t i=1;i<=24;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
 832
 833   // Branching ratio
 834   fBR = new TH1F("fBR","fBR",24,0,24);
 835   for(Int_t i=1;i<=24;i++) fBR->SetBinContent(i,0.015); // 1.5% PDG2010
 836
 837   // Tracking efficiency
 838   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",24,0,24);;
 839   for(Int_t i=1;i<=24;i++) fTrackingEff->SetBinContent(i,0.12); // 3% (1% per track)
 840
 841
 842   // Raw yield extraction
 843   fRawYield = new TH1F("fRawYield","fRawYield",24,0,24);
 844   for(Int_t i=1;i<=24;i++) fRawYield->SetBinContent(i,0.15);  //4%
 845
 846   // Cuts efficiency (from cuts variation)
 847   fCutsEff = new TH1F("fCutsEff","fCutsEff",24,0,24);
 848   for(Int_t i=1;i<=24;i++) fCutsEff->SetBinContent(i,0.10); // 10%
 849
 850   // PID efficiency (from PID/noPID)
 851   fPIDEff = new TH1F("fPIDEff","fPIDEff",24,0,24);
 852   for(Int_t i=1;i<=24;i++) fPIDEff->SetBinContent(i,0.04); // 3%
 853
 854
 855   // MC dN/dpt  (copied from D0 : will update later)
 856   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",24,0,24);
 857   for(Int_t i=1;i<=24;i++) fMCPtShape->SetBinContent(i,0.);
 858   fMCPtShape->SetBinContent(1,0.03);
 859   fMCPtShape->SetBinContent(2,0.03);
 860
 861   return;
 862
 863 }
 864 //--------------------------------------------------------------------------
 865 Double_t AliHFSystErr::GetCutsEffErr(Double_t pt) const {
 866   //
 867   // Get error
 868   //
 869
 870   Int_t bin=fCutsEff->FindBin(pt);
 871
 872   return fCutsEff->GetBinContent(bin);
 873 }
 874 //--------------------------------------------------------------------------
 875 Double_t AliHFSystErr::GetMCPtShapeErr(Double_t pt) const {
 876   //
 877   // Get error
 878   //
 879
 880   Int_t bin=fMCPtShape->FindBin(pt);
 881
 882   return fMCPtShape->GetBinContent(bin);
 883 }
 884 //--------------------------------------------------------------------------
 885 Double_t AliHFSystErr::GetSeleEffErr(Double_t pt) const {
 886   //
 887   // Get error
 888   //
 889
 890   Double_t err=GetCutsEffErr(pt)*GetCutsEffErr(pt)+GetMCPtShapeErr(pt)*GetMCPtShapeErr(pt);
 891
 892   return TMath::Sqrt(err);
 893 }
 894 //--------------------------------------------------------------------------
 895 Double_t AliHFSystErr::GetPIDEffErr(Double_t pt) const {
 896   //
 897   // Get error
 898   //
 899
 900   Int_t bin=fPIDEff->FindBin(pt);
 901
 902   return fPIDEff->GetBinContent(bin);
 903 }
 904 //--------------------------------------------------------------------------
 905 Double_t AliHFSystErr::GetTrackingEffErr(Double_t pt) const {
 906   //
 907   // Get error
 908   //
 909
 910   Int_t bin=fTrackingEff->FindBin(pt);
 911
 912   return fTrackingEff->GetBinContent(bin);
 913 }
 914 //--------------------------------------------------------------------------
 915 Double_t AliHFSystErr::GetRawYieldErr(Double_t pt) const {
 916   //
 917   // Get error
 918   //
 919
 920   Int_t bin=fRawYield->FindBin(pt);
 921
 922   return fRawYield->GetBinContent(bin);
 923 }
 924 //--------------------------------------------------------------------------
 925 Double_t AliHFSystErr::GetPartAntipartErr(Double_t pt) const {
 926   //
 927   // Get error
 928   //
 929
 930   Int_t bin=fPartAntipart->FindBin(pt);
 931
 932   return fPartAntipart->GetBinContent(bin);
 933 }
 934 //--------------------------------------------------------------------------
 935 Double_t AliHFSystErr::GetTotalSystErr(Double_t pt,Double_t feeddownErr) const {
 936   //
 937   // Get total syst error (except norm. error)
 938   //
 939
 940   Double_t err=0.;
 941
 942   if(fRawYield) err += GetRawYieldErr(pt)*GetRawYieldErr(pt);
 943   if(fTrackingEff) err += GetTrackingEffErr(pt)*GetTrackingEffErr(pt);
 944   if(fBR) err += GetBRErr()*GetBRErr();
 945   if(fCutsEff) err += GetCutsEffErr(pt)*GetCutsEffErr(pt);
 946   if(fPIDEff) err += GetPIDEffErr(pt)*GetPIDEffErr(pt);
 947   if(fMCPtShape) err += GetMCPtShapeErr(pt)*GetMCPtShapeErr(pt);
 948   if(fPartAntipart) err += GetPartAntipartErr(pt)*GetPartAntipartErr(pt);
 949
 950   err += feeddownErr*feeddownErr;
 951
 952   return TMath::Sqrt(err);
 953 }
 954 //---------------------------------------------------------------------------
 955 void AliHFSystErr::DrawErrors(TGraphAsymmErrors *grErrFeeddown) const {
 956   //
 957   // Draw errors
 958   //
 959   gStyle->SetOptStat(0);
 960
 961   TCanvas *cSystErr = new TCanvas("cSystErr","Systematic Errors",300,80,640,500);
 962   cSystErr->Range(0.20,-0.5,18.4,0.34);
 963   cSystErr->SetRightMargin(0.318);
 964   cSystErr->SetFillColor(0);
 965
 966   TH2F *hFrame = new TH2F("hFrame","Systematic errors; p_{t} [GeV/c]; Relative Error",20,0,20,100,-1,+1);
 967   hFrame->SetAxisRange(1.,15.9,"X");
 968   hFrame->SetAxisRange(-0.5,0.5,"Y");
 969   hFrame->Draw();
 970
 971   TLegend *leg = new TLegend(0.69,0.44,0.98,0.86,NULL,"brNDC");
 972   leg->SetTextSize(0.03601695);
 973   leg->SetFillStyle(0);
 974   leg->SetBorderSize(0);
 975
 976   TH1F *hTotErr=new TH1F("hTotErr","",24,0,24);
 977   Int_t nbins = fNorm->GetNbinsX();
 978   TGraphAsymmErrors *gTotErr = new TGraphAsymmErrors(nbins);
 979   for(Int_t i=1;i<=24;i++) {
 980     Double_t pt = hTotErr->GetBinCenter(i);
 981     Double_t ptwidth = hTotErr->GetBinWidth(i);
 982
 983     if(grErrFeeddown) {
 984       Double_t x=0., y=0., errxl=0., errxh=0., erryl=0., erryh=0.;
 985       Double_t toterryl=0., toterryh=0.;
 986       for(Int_t j=0; j<grErrFeeddown->GetN(); j++) {
 987         grErrFeeddown->GetPoint(j,x,y);
 988         errxh = grErrFeeddown->GetErrorXhigh(j);
 989         errxl = grErrFeeddown->GetErrorXlow(j);
 990         if ( ( (x-errxl) <= pt) && ( (x+errxl) >= pt) ) {
 991           erryh = grErrFeeddown->GetErrorYhigh(j);
 992           erryl = grErrFeeddown->GetErrorYlow(j);
 993         }
 994       }
 995       if (erryl>=1e-3) toterryl = GetTotalSystErr(pt,erryl);
 996       else toterryl = GetTotalSystErr(pt);
 997       if (erryh>=1e-3) toterryh = GetTotalSystErr(pt,erryh);
 998       else toterryh = GetTotalSystErr(pt);
 999
1000       hTotErr->SetBinContent(i,toterryh);
1001       gTotErr->SetPoint(i,pt,0.);
1002       gTotErr->SetPointError(i,ptwidth/2.,ptwidth/2.,toterryl,toterryh); // i, exl, exh, eyl, eyh
1003     }
1004     else {
1005       hTotErr->SetBinContent(i,GetTotalSystErr(pt));
1006       gTotErr->SetPoint(i,pt,0.);
1007       gTotErr->SetPointError(i,ptwidth/2.,ptwidth/2.,GetTotalSystErr(pt),GetTotalSystErr(pt)); // i, exl, exh, eyl, eyh
1008     }
1009
1010   }
1011   gTotErr->SetLineColor(kBlack);
1012   gTotErr->SetFillColor(kRed);
1013   gTotErr->SetFillStyle(3002);
1014   gTotErr->Draw("2");
1015   leg->AddEntry(gTotErr,"Total (excl. norm.)","f");
1016 //   hTotErr->SetLineColor(1);
1017 //   hTotErr->SetLineWidth(3);
1018 //   hTotErr->Draw("same");
1019 //   leg->AddEntry(hTotErr,"Total (excl. norm.)","l");
1020
1021
1022   fNorm->SetFillColor(1);
1023   fNorm->SetFillStyle(3002);
1024   //fNorm->Draw("same");
1025   //TH1F *hNormRefl = ReflectHisto(fNorm);
1026   //hNormRefl->Draw("same");
1027   leg->AddEntry(fNorm,"Normalization (10%)","");
1028
1029   if(grErrFeeddown) {
1030     grErrFeeddown->SetFillColor(kTeal-8);
1031     grErrFeeddown->SetFillStyle(3001);
1032     grErrFeeddown->Draw("2");
1033     leg->AddEntry(grErrFeeddown,"Feed-down from B","f");
1034   }
1035   if(fTrackingEff) {
1036     fTrackingEff->SetFillColor(4);
1037     fTrackingEff->SetFillStyle(3006);
1038     fTrackingEff->Draw("same");
1039     TH1F *hTrackingEffRefl = ReflectHisto(fTrackingEff);
1040     hTrackingEffRefl->Draw("same");
1041     leg->AddEntry(fTrackingEff,"Tracking efficiency","f");
1042   }
1043   if(fBR) {
1044     fBR->SetFillColor(6);
1045     fBR->SetFillStyle(3005);
1046     //fBR->SetFillStyle(3020);
1047     fBR->Draw("same");
1048     TH1F *hBRRefl = ReflectHisto(fBR);
1049     hBRRefl->Draw("same");
1050     leg->AddEntry(fBR,"Branching ratio","f");
1051   }
1052   if(fRawYield) {
1053     Int_t ci;   // for color index setting
1054     ci = TColor::GetColor("#00cc00");
1055     fRawYield->SetLineColor(ci);
1056     //    fRawYield->SetLineColor(3);
1057     fRawYield->SetLineWidth(3);
1058     fRawYield->Draw("same");
1059     TH1F *hRawYieldRefl = ReflectHisto(fRawYield);
1060     hRawYieldRefl->Draw("same");
1061     leg->AddEntry(fRawYield,"Yield extraction","l");
1062   }
1063   if(fCutsEff) {
1064     fCutsEff->SetLineColor(4);
1065     fCutsEff->SetLineWidth(3);
1066     fCutsEff->Draw("same");
1067     TH1F *hCutsEffRefl = ReflectHisto(fCutsEff);
1068     hCutsEffRefl->Draw("same");
1069     leg->AddEntry(fCutsEff,"Cuts efficiency","l");
1070   }
1071   if(fPIDEff) {
1072     fPIDEff->SetLineColor(7);
1073     fPIDEff->SetLineWidth(3);
1074     fPIDEff->Draw("same");
1075     TH1F *hPIDEffRefl = ReflectHisto(fPIDEff);
1076     hPIDEffRefl->Draw("same");
1077     leg->AddEntry(fPIDEff,"PID efficiency","l");
1078   }
1079   if(fMCPtShape) {
1080     Int_t ci = TColor::GetColor("#9933ff");
1081     fMCPtShape->SetLineColor(ci);
1082     //    fMCPtShape->SetLineColor(8);
1083     fMCPtShape->SetLineWidth(3);
1084     fMCPtShape->Draw("same");
1085     TH1F *hMCPtShapeRefl = ReflectHisto(fMCPtShape);
1086     hMCPtShapeRefl->Draw("same");
1087     leg->AddEntry(fMCPtShape,"MC p_{t} shape","l");
1088   }
1089   if(fPartAntipart) {
1090     Int_t ci = TColor::GetColor("#ff6600");
1091     fPartAntipart->SetLineColor(ci);
1092     //    fPartAntipart->SetLineColor(9);
1093     fPartAntipart->SetLineWidth(3);
1094     fPartAntipart->Draw("same");
1095     TH1F *hPartAntipartRefl = ReflectHisto(fPartAntipart);
1096     hPartAntipartRefl->Draw("same");
1097     leg->AddEntry(fPartAntipart,"D = #bar{D}","l");
1098   }
1099
1100
1101   leg->Draw();
1102
1103   cSystErr->SaveAs("RelativeSystematics.eps");
1104
1105   return;
1106 }
1107 //-------------------------------------------------------------------------
1108 TH1F* AliHFSystErr::ReflectHisto(TH1F *hin) const {
1109   //
1110   // Clones and reflects histogram
1111   //
1112   TH1F *hout=(TH1F*)hin->Clone("hout");
1113   hout->Scale(-1.);
1114
1115   return hout;
1116 }