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: A. Rossi
 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=2;i<=12;i++) fTrackingEff->SetBinContent(i,0.05+0.005*(Float_t)i);
 195
 196
 197   // Raw yield extraction
 198   fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
 199   fRawYield->SetBinContent(1,0);
 200   fRawYield->SetBinContent(2,0);
 201   fRawYield->SetBinContent(3,0.05);
 202   fRawYield->SetBinContent(4,0.05);
 203   fRawYield->SetBinContent(5,0.10);
 204   fRawYield->SetBinContent(6,0.10);
 205   for(Int_t i=7;i<=12;i++) fRawYield->SetBinContent(i,0.05);
 206   for(Int_t i=13;i<=20;i++) fRawYield->SetBinContent(i,0);
 207
 208   // Cuts efficiency (from cuts variation)
 209   fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
 210   fCutsEff->SetBinContent(1,0.);
 211   fCutsEff->SetBinContent(2,0.);
 212   fCutsEff->SetBinContent(3,0.25);
 213   fCutsEff->SetBinContent(4,0.18);
 214   fCutsEff->SetBinContent(5,0.18);
 215   fCutsEff->SetBinContent(6,0.15);
 216   fCutsEff->SetBinContent(7,0.15);
 217   fCutsEff->SetBinContent(8,0.15);
 218   fCutsEff->SetBinContent(9,0.15);
 219   fCutsEff->SetBinContent(10,0.15);
 220   fCutsEff->SetBinContent(11,0.15);
 221   fCutsEff->SetBinContent(12,0.15);
 222   for(Int_t i=13;i<=20;i++) fCutsEff->SetBinContent(i,0.);
 223
 224   // PID efficiency (from PID/noPID)
 225   fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
 226   for(Int_t i=3;i<=12;i++) fPIDEff->SetBinContent(i,0.05);
 227   fPIDEff->SetBinContent(5,0.10);
 228   fPIDEff->SetBinContent(6,0.10);
 229
 230   // MC dN/dpt
 231   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
 232   for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,(Float_t)i*0.006);
 233
 234   // particle-antiparticle
 235   fPartAntipart = new TH1F("fPartAntipart","fPartAntipart",20,0,20);
 236   for(Int_t i=3;i<=12;i++) fPartAntipart->SetBinContent(i,0.05);
 237   fPartAntipart->SetBinContent(4,0.10);
 238   fPartAntipart->SetBinContent(7,0.10);
 239   fPartAntipart->SetBinContent(8,0.10);
 240
 241   return;
 242 }
 243 //--------------------------------------------------------------------------
 244 void AliHFSystErr::InitD0toKpi2010PbPb4080() {
 245   //
 246   // D0->Kpi syst errors. Responsible: A. Rossi
 247   //   2010 PbPb sample, 40-80 CC
 248   //
 249
 250   // Normalization
 251   fNorm = new TH1F("fNorm","fNorm",20,0,20);
 252   for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
 253
 254   // Branching ratio
 255   fBR = new TH1F("fBR","fBR",20,0,20);
 256   for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.012); // 1.2% PDG2010
 257
 258   // Tracking efficiency
 259   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
 260   for(Int_t i=2;i<=12;i++) fTrackingEff->SetBinContent(i,0.5*(0.05+0.005*(Float_t)i));
 261
 262
 263   // Raw yield extraction
 264   fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
 265   fRawYield->SetBinContent(1,0);
 266   fRawYield->SetBinContent(2,0);
 267   for(Int_t i=3;i<=12;i++) fRawYield->SetBinContent(i,0.05);
 268   for(Int_t i=13;i<=20;i++) fRawYield->SetBinContent(i,0);
 269
 270   // Cuts efficiency (from cuts variation)
 271   fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
 272   fCutsEff->SetBinContent(1,0.);
 273   fCutsEff->SetBinContent(2,0.);
 274   fCutsEff->SetBinContent(3,0.23);
 275   fCutsEff->SetBinContent(4,0.15);
 276   fCutsEff->SetBinContent(5,0.15);
 277   fCutsEff->SetBinContent(6,0.15);
 278   fCutsEff->SetBinContent(7,0.15);
 279   fCutsEff->SetBinContent(8,0.15);
 280   fCutsEff->SetBinContent(9,0.15);
 281   fCutsEff->SetBinContent(10,0.15);
 282   fCutsEff->SetBinContent(11,0.15);
 283   fCutsEff->SetBinContent(12,0.15);
 284   for(Int_t i=13;i<=20;i++) fCutsEff->SetBinContent(i,0.);
 285
 286   // PID efficiency (from PID/noPID)
 287   fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
 288   for(Int_t i=3;i<=12;i++) fPIDEff->SetBinContent(i,0.05);
 289
 290   // MC dN/dpt
 291   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
 292   for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,(Float_t)i*0.006);
 293
 294   // particle-antiparticle
 295   fPartAntipart = new TH1F("fPartAntipart","fPartAntipart",20,0,20);
 296   for(Int_t i=3;i<=12;i++) fPartAntipart->SetBinContent(i,0.05);
 297
 298   return;
 299 }
 300 //--------------------------------------------------------------------------
 301 void AliHFSystErr::InitDplustoKpipi2010pp() {
 302   //
 303   // D+->Kpipi syst errors. Responsible: R. Bala
 304   //  2010 pp sample
 305   //
 306
 307  // Normalization
 308   fNorm = new TH1F("fNorm","fNorm",20,0,20);
 309   for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
 310
 311   // Branching ratio
 312   fBR = new TH1F("fBR","fBR",20,0,20);
 313   for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.04); // 4% PDG2010
 314
 315   // Tracking efficiency
 316   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
 317   for(Int_t i=1;i<=20;i++) fTrackingEff->SetBinContent(i,0.03); // 3% (1% per track)
 318
 319
 320   // Raw yield extraction
 321   fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
 322   fRawYield->SetBinContent(1,1);
 323   fRawYield->SetBinContent(2,1);
 324   fRawYield->SetBinContent(3,0.20);
 325   for(Int_t i=4;i<=20;i++) fRawYield->SetBinContent(i,0.055);  //5 to 10%
 326
 327   // Cuts efficiency (from cuts variation)
 328   fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
 329   for(Int_t i=1;i<=20;i++) fCutsEff->SetBinContent(i,0.10); // 10%
 330
 331   // PID efficiency (from PID/noPID)
 332   fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
 333   for(Int_t i=1;i<=20;i++) fPIDEff->SetBinContent(i,0.05); // 5%
 334   fPIDEff->SetBinContent(3,0.13); // 13%
 335
 336
 337   // MC dN/dpt  (copied from D0 : will update later)
 338   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
 339   for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,(Float_t)i*0.006);
 340
 341
 342   // particle-antiparticle
 343   /*
 344   fPartAntipart = new TH1F("fPartAntipart","fPartAntipart",20,0,20);
 345   fPartAntipart->SetBinContent(1,1);
 346   fPartAntipart->SetBinContent(2,1);
 347   fPartAntipart->SetBinContent(3,0.12);
 348   for(Int_t i=4;i<=20;i++) fPartAntipart->SetBinContent(i,0.05);   //5 to 12%
 349   */
 350   return;
 351 }
 352 //--------------------------------------------------------------------------
 353 void AliHFSystErr::InitDplustoKpipi2010PbPb020() {
 354   //
 355   // D+->Kpipi syst errors. Responsible: ??
 356   //  2010 PbPb sample, 0-20 CC
 357   //
 358
 359
 360  // Normalization
 361   fNorm = new TH1F("fNorm","fNorm",20,0,20);
 362   for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
 363
 364   // Branching ratio
 365   fBR = new TH1F("fBR","fBR",20,0,20);
 366   for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.04); // 4% PDG2010
 367
 368   // Tracking efficiency
 369   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
 370   for(Int_t i=1;i<=20;i++) fTrackingEff->SetBinContent(i,0.03); // 3% (1% per track)
 371
 372
 373   // Raw yield extraction
 374   fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
 375   for(Int_t i=1;i<=20;i++) fRawYield->SetBinContent(i,1);  //5 to 10%
 376   for(Int_t i=9;i<=12;i++) fRawYield->SetBinContent(i,31./156.);
 377
 378   // Cuts efficiency (from cuts variation)
 379   fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
 380   for(Int_t i=1;i<=20;i++) fCutsEff->SetBinContent(i,0.10); // 10%
 381
 382   // PID efficiency (from PID/noPID)
 383   fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
 384   for(Int_t i=1;i<=20;i++) fPIDEff->SetBinContent(i,0.05); // 5%
 385   fPIDEff->SetBinContent(3,0.13); // 13%
 386
 387
 388   // MC dN/dpt  (copied from D0 : will update later)
 389   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
 390   for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,(Float_t)i*0.006);
 391
 392
 393   // particle-antiparticle
 394   /*
 395   fPartAntipart = new TH1F("fPartAntipart","fPartAntipart",20,0,20);
 396   fPartAntipart->SetBinContent(1,1);
 397   fPartAntipart->SetBinContent(2,1);
 398   fPartAntipart->SetBinContent(3,0.12);
 399   for(Int_t i=4;i<=20;i++) fPartAntipart->SetBinContent(i,0.05);   //5 to 12%
 400   */
 401   return;
 402 }
 403
 404 //--------------------------------------------------------------------------
 405 void AliHFSystErr::InitDplustoKpipi2010PbPb4080() {
 406   //
 407   // D+->Kpipi syst errors. Responsible: ??
 408   //  2010 PbPb sample, 40-80 CC
 409   //
 410
 411
 412  // Normalization
 413   fNorm = new TH1F("fNorm","fNorm",20,0,20);
 414   for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
 415
 416   // Branching ratio
 417   fBR = new TH1F("fBR","fBR",20,0,20);
 418   for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.04); // 4% PDG2010
 419
 420   // Tracking efficiency
 421   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
 422   for(Int_t i=1;i<=20;i++) fTrackingEff->SetBinContent(i,0.03); // 3% (1% per track)
 423
 424
 425   // Raw yield extraction
 426   fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
 427   fRawYield->SetBinContent(1,1);
 428   fRawYield->SetBinContent(2,1);
 429   fRawYield->SetBinContent(3,1);
 430   fRawYield->SetBinContent(4,1.);
 431   fRawYield->SetBinContent(5,4./23.);
 432   fRawYield->SetBinContent(6,4./23.);
 433   fRawYield->SetBinContent(7,5./20.);
 434   fRawYield->SetBinContent(8,5./20.);
 435   for(Int_t i=9;i<=12;i++) fRawYield->SetBinContent(i,1./16.);
 436   for(Int_t i=13;i<=20;i++) fRawYield->SetBinContent(i,1);  //5 to 10%
 437
 438   // Cuts efficiency (from cuts variation)
 439   fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
 440   for(Int_t i=1;i<=20;i++) fCutsEff->SetBinContent(i,0.10); // 10%
 441
 442   // PID efficiency (from PID/noPID)
 443   fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
 444   for(Int_t i=1;i<=20;i++) fPIDEff->SetBinContent(i,0.05); // 5%
 445   fPIDEff->SetBinContent(3,0.13); // 13%
 446
 447
 448   // MC dN/dpt  (copied from D0 : will update later)
 449   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
 450   for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,(Float_t)i*0.006);
 451
 452
 453   // particle-antiparticle
 454   /*
 455   fPartAntipart = new TH1F("fPartAntipart","fPartAntipart",20,0,20);
 456   fPartAntipart->SetBinContent(1,1);
 457   fPartAntipart->SetBinContent(2,1);
 458   fPartAntipart->SetBinContent(3,0.12);
 459   for(Int_t i=4;i<=20;i++) fPartAntipart->SetBinContent(i,0.05);   //5 to 12%
 460   */
 461   return;
 462 }
 463
 464 //--------------------------------------------------------------------------
 465 void AliHFSystErr::InitDstartoD0pi2010pp() {
 466   //
 467   // D*+->D0pi syst errors. Responsible: A. Grelli, Y. Wang
 468   //  2010 pp sample
 469   //
 470
 471  // Normalization
 472   fNorm = new TH1F("fNorm","fNorm",20,0,20);
 473   for(Int_t i=1;i<=20;i++) fNorm->SetBinContent(i,0.10); // 10% error on sigmaV0and
 474
 475   // Branching ratio
 476   fBR = new TH1F("fBR","fBR",20,0,20);
 477   for(Int_t i=1;i<=20;i++) fBR->SetBinContent(i,0.015); // 1.5% PDG2010
 478
 479   // Tracking efficiency
 480   fTrackingEff = new TH1F("fTrackingEff","fTrackingEff",20,0,20);
 481   fTrackingEff->SetBinContent(1,0.12);
 482   fTrackingEff->SetBinContent(2,0.08);
 483   fTrackingEff->SetBinContent(3,0.05);
 484   for(Int_t i=4;i<=20;i++) fTrackingEff->SetBinContent(i,0.03); // 3% (1% per track)
 485
 486
 487   // Raw yield extraction
 488   fRawYield = new TH1F("fRawYield","fRawYield",20,0,20);
 489   fRawYield->SetBinContent(1,1);
 490   fRawYield->SetBinContent(2,0.12);
 491   fRawYield->SetBinContent(3,0.09);
 492   fRawYield->SetBinContent(4,0.08);
 493   fRawYield->SetBinContent(5,0.06);
 494   for(Int_t i=5;i<=20;i++) fRawYield->SetBinContent(i,0.04);  //4%
 495
 496   // Cuts efficiency (from cuts variation)
 497   fCutsEff = new TH1F("fCutsEff","fCutsEff",20,0,20);
 498   for(Int_t i=1;i<=20;i++) fCutsEff->SetBinContent(i,0.10); // 10%
 499
 500   // PID efficiency (from PID/noPID)
 501   fPIDEff = new TH1F("fPIDEff","fPIDEff",20,0,20);
 502   for(Int_t i=1;i<=20;i++) fPIDEff->SetBinContent(i,0.04); // 3%
 503   fPIDEff->SetBinContent(1,1); // 100%
 504   fPIDEff->SetBinContent(2,1); // 100%
 505   fPIDEff->SetBinContent(3,0.05); // 5%
 506   fPIDEff->SetBinContent(4,0.05); // 5%
 507   fPIDEff->SetBinContent(5,0.05); // 5%
 508
 509
 510   // MC dN/dpt  (copied from D0 : will update later)
 511   fMCPtShape = new TH1F("fMCPtShape","fMCPtShape",20,0,20);
 512   for(Int_t i=1;i<=20;i++) fMCPtShape->SetBinContent(i,(Float_t)i*0.006);
 513
 514   return;
 515 }
 516 //--------------------------------------------------------------------------
 517 void AliHFSystErr::InitDstartoD0pi2010PbPb020() {
 518   //
 519   // D*+->D0pi syst errors. Responsible: ??
 520   //  2010 PbPb sample, 0-20 CC
 521   //
 522
 523   AliInfo("Not yet implemented");
 524   return;
 525 }
 526 //--------------------------------------------------------------------------
 527 void AliHFSystErr::InitDstartoD0pi2010PbPb4080() {
 528   //
 529   // D*+->D0pi syst errors. Responsible: ??
 530   //  2010 PbPb sample, 40-80 CC
 531   //
 532
 533   AliInfo("Not yet implemented");
 534   return;
 535 }
 536 //--------------------------------------------------------------------------
 537 Double_t AliHFSystErr::GetCutsEffErr(Double_t pt) const {
 538   //
 539   // Get error
 540   //
 541
 542   Int_t bin=fCutsEff->FindBin(pt);
 543
 544   return fCutsEff->GetBinContent(bin);
 545 }
 546 //--------------------------------------------------------------------------
 547 Double_t AliHFSystErr::GetMCPtShapeErr(Double_t pt) const {
 548   //
 549   // Get error
 550   //
 551
 552   Int_t bin=fMCPtShape->FindBin(pt);
 553
 554   return fMCPtShape->GetBinContent(bin);
 555 }
 556 //--------------------------------------------------------------------------
 557 Double_t AliHFSystErr::GetSeleEffErr(Double_t pt) const {
 558   //
 559   // Get error
 560   //
 561
 562   Double_t err=GetCutsEffErr(pt)*GetCutsEffErr(pt)+GetMCPtShapeErr(pt)*GetMCPtShapeErr(pt);
 563
 564   return TMath::Sqrt(err);
 565 }
 566 //--------------------------------------------------------------------------
 567 Double_t AliHFSystErr::GetPIDEffErr(Double_t pt) const {
 568   //
 569   // Get error
 570   //
 571
 572   Int_t bin=fPIDEff->FindBin(pt);
 573
 574   return fPIDEff->GetBinContent(bin);
 575 }
 576 //--------------------------------------------------------------------------
 577 Double_t AliHFSystErr::GetTrackingEffErr(Double_t pt) const {
 578   //
 579   // Get error
 580   //
 581
 582   Int_t bin=fTrackingEff->FindBin(pt);
 583
 584   return fTrackingEff->GetBinContent(bin);
 585 }
 586 //--------------------------------------------------------------------------
 587 Double_t AliHFSystErr::GetRawYieldErr(Double_t pt) const {
 588   //
 589   // Get error
 590   //
 591
 592   Int_t bin=fRawYield->FindBin(pt);
 593
 594   return fRawYield->GetBinContent(bin);
 595 }
 596 //--------------------------------------------------------------------------
 597 Double_t AliHFSystErr::GetPartAntipartErr(Double_t pt) const {
 598   //
 599   // Get error
 600   //
 601
 602   Int_t bin=fPartAntipart->FindBin(pt);
 603
 604   return fPartAntipart->GetBinContent(bin);
 605 }
 606 //--------------------------------------------------------------------------
 607 Double_t AliHFSystErr::GetTotalSystErr(Double_t pt,Double_t feeddownErr) const {
 608   //
 609   // Get total syst error (except norm. error)
 610   //
 611
 612   Double_t err=0.;
 613
 614   if(fRawYield) err += GetRawYieldErr(pt)*GetRawYieldErr(pt);
 615   if(fTrackingEff) err += GetTrackingEffErr(pt)*GetTrackingEffErr(pt);
 616   if(fBR) err += GetBRErr()*GetBRErr();
 617   if(fCutsEff) err += GetCutsEffErr(pt)*GetCutsEffErr(pt);
 618   if(fPIDEff) err += GetPIDEffErr(pt)*GetPIDEffErr(pt);
 619   if(fMCPtShape) err += GetMCPtShapeErr(pt)*GetMCPtShapeErr(pt);
 620   if(fPartAntipart) err += GetPartAntipartErr(pt)*GetPartAntipartErr(pt);
 621
 622   err += feeddownErr*feeddownErr;
 623
 624   return TMath::Sqrt(err);
 625 }
 626 //---------------------------------------------------------------------------
 627 void AliHFSystErr::DrawErrors(TGraphAsymmErrors *grErrFeeddown) const {
 628   //
 629   // Draw errors
 630   //
 631   gStyle->SetOptStat(0);
 632
 633   TCanvas *cSystErr = new TCanvas("cSystErr","Systematic Errors",300,80,640,500);
 634   cSystErr->Range(0.20,-0.5,18.4,0.34);
 635   cSystErr->SetRightMargin(0.318);
 636   cSystErr->SetFillColor(0);
 637
 638   TH2F *hFrame = new TH2F("hFrame","Systematic errors; p_{t} [GeV/c]; Relative Error",20,0,20,100,-1,+1);
 639   hFrame->SetAxisRange(2.,11.9,"X");
 640   hFrame->SetAxisRange(-0.5,0.5,"Y");
 641   hFrame->Draw();
 642
 643   TLegend *leg = new TLegend(0.69,0.44,0.98,0.86,NULL,"brNDC");
 644   leg->SetTextSize(0.03601695);
 645   leg->SetFillStyle(0);
 646   leg->SetBorderSize(0);
 647
 648   TH1F *hTotErr=new TH1F("hTotErr","",20,0,20);
 649   Int_t nbins = fNorm->GetNbinsX();
 650   TGraphAsymmErrors *gTotErr = new TGraphAsymmErrors(nbins);
 651   for(Int_t i=1;i<=20;i++) {
 652     Double_t pt = hTotErr->GetBinCenter(i);
 653     Double_t ptwidth = hTotErr->GetBinWidth(i);
 654
 655     if(grErrFeeddown) {
 656       Double_t x=0., y=0., errxl=0., errxh=0., erryl=0., erryh=0.;
 657       Double_t toterryl=0., toterryh=0.;
 658       for(Int_t j=0; j<grErrFeeddown->GetN(); j++) {
 659         grErrFeeddown->GetPoint(j,x,y);
 660         errxh = grErrFeeddown->GetErrorXhigh(j);
 661         errxl = grErrFeeddown->GetErrorXlow(j);
 662         if ( ( (x-errxl) <= pt) && ( (x+errxl) >= pt) ) {
 663           erryh = grErrFeeddown->GetErrorYhigh(j);
 664           erryl = grErrFeeddown->GetErrorYlow(j);
 665         }
 666       }
 667       if (erryl>=1e-3) toterryl = GetTotalSystErr(pt,erryl);
 668       else toterryl = GetTotalSystErr(pt);
 669       if (erryh>=1e-3) toterryh = GetTotalSystErr(pt,erryh);
 670       else toterryh = GetTotalSystErr(pt);
 671
 672       hTotErr->SetBinContent(i,toterryh);
 673       gTotErr->SetPoint(i,pt,0.);
 674       gTotErr->SetPointError(i,ptwidth/2.,ptwidth/2.,toterryl,toterryh); // i, exl, exh, eyl, eyh
 675     }
 676     else {
 677       hTotErr->SetBinContent(i,GetTotalSystErr(pt));
 678       gTotErr->SetPoint(i,pt,0.);
 679       gTotErr->SetPointError(i,ptwidth/2.,ptwidth/2.,GetTotalSystErr(pt),GetTotalSystErr(pt)); // i, exl, exh, eyl, eyh
 680     }
 681
 682   }
 683   gTotErr->SetLineColor(kBlack);
 684   gTotErr->SetFillColor(kRed);
 685   gTotErr->SetFillStyle(3002);
 686   gTotErr->Draw("2");
 687   leg->AddEntry(gTotErr,"Total (excl. norm.)","f");
 688 //   hTotErr->SetLineColor(1);
 689 //   hTotErr->SetLineWidth(3);
 690 //   hTotErr->Draw("same");
 691 //   leg->AddEntry(hTotErr,"Total (excl. norm.)","l");
 692
 693
 694   fNorm->SetFillColor(1);
 695   fNorm->SetFillStyle(3002);
 696   //fNorm->Draw("same");
 697   //TH1F *hNormRefl = ReflectHisto(fNorm);
 698   //hNormRefl->Draw("same");
 699   leg->AddEntry(fNorm,"Normalization (10%)","");
 700
 701   if(grErrFeeddown) {
 702     grErrFeeddown->SetFillColor(kTeal-8);
 703     grErrFeeddown->SetFillStyle(3001);
 704     grErrFeeddown->Draw("2");
 705     leg->AddEntry(grErrFeeddown,"Feed-down from B","f");
 706   }
 707   if(fTrackingEff) {
 708     fTrackingEff->SetFillColor(4);
 709     fTrackingEff->SetFillStyle(3006);
 710     fTrackingEff->Draw("same");
 711     TH1F *hTrackingEffRefl = ReflectHisto(fTrackingEff);
 712     hTrackingEffRefl->Draw("same");
 713     leg->AddEntry(fTrackingEff,"Tracking efficiency","f");
 714   }
 715   if(fBR) {
 716     fBR->SetFillColor(6);
 717     fBR->SetFillStyle(3005);
 718     //fBR->SetFillStyle(3020);
 719     fBR->Draw("same");
 720     TH1F *hBRRefl = ReflectHisto(fBR);
 721     hBRRefl->Draw("same");
 722     leg->AddEntry(fBR,"Branching ratio","f");
 723   }
 724   if(fRawYield) {
 725     Int_t ci;   // for color index setting
 726     ci = TColor::GetColor("#00cc00");
 727     fRawYield->SetLineColor(ci);
 728     //    fRawYield->SetLineColor(3);
 729     fRawYield->SetLineWidth(3);
 730     fRawYield->Draw("same");
 731     TH1F *hRawYieldRefl = ReflectHisto(fRawYield);
 732     hRawYieldRefl->Draw("same");
 733     leg->AddEntry(fRawYield,"Yield extraction","l");
 734   }
 735   if(fCutsEff) {
 736     fCutsEff->SetLineColor(4);
 737     fCutsEff->SetLineWidth(3);
 738     fCutsEff->Draw("same");
 739     TH1F *hCutsEffRefl = ReflectHisto(fCutsEff);
 740     hCutsEffRefl->Draw("same");
 741     leg->AddEntry(fCutsEff,"Cuts efficiency","l");
 742   }
 743   if(fPIDEff) {
 744     fPIDEff->SetLineColor(7);
 745     fPIDEff->SetLineWidth(3);
 746     fPIDEff->Draw("same");
 747     TH1F *hPIDEffRefl = ReflectHisto(fPIDEff);
 748     hPIDEffRefl->Draw("same");
 749     leg->AddEntry(fPIDEff,"PID efficiency","l");
 750   }
 751   if(fMCPtShape) {
 752     Int_t ci = TColor::GetColor("#9933ff");
 753     fMCPtShape->SetLineColor(ci);
 754     //    fMCPtShape->SetLineColor(8);
 755     fMCPtShape->SetLineWidth(3);
 756     fMCPtShape->Draw("same");
 757     TH1F *hMCPtShapeRefl = ReflectHisto(fMCPtShape);
 758     hMCPtShapeRefl->Draw("same");
 759     leg->AddEntry(fMCPtShape,"MC p_{t} shape","l");
 760   }
 761   if(fPartAntipart) {
 762     Int_t ci = TColor::GetColor("#ff6600");
 763     fPartAntipart->SetLineColor(ci);
 764     //    fPartAntipart->SetLineColor(9);
 765     fPartAntipart->SetLineWidth(3);
 766     fPartAntipart->Draw("same");
 767     TH1F *hPartAntipartRefl = ReflectHisto(fPartAntipart);
 768     hPartAntipartRefl->Draw("same");
 769     leg->AddEntry(fPartAntipart,"D = #bar{D}","l");
 770   }
 771
 772
 773   leg->Draw();
 774
 775   cSystErr->SaveAs("RelativeSystematics.eps");
 776
 777   return;
 778 }
 779 //-------------------------------------------------------------------------
 780 TH1F* AliHFSystErr::ReflectHisto(TH1F *hin) const {
 781   //
 782   // Clones and reflects histogram
 783   //
 784   TH1F *hout=(TH1F*)hin->Clone("hout");
 785   hout->Scale(-1.);
 786
 787   return hout;
 788 }