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