PWG3/dielectron/AliDielectronSignalFunc.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-2009, 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 //                Dielectron SignalFunc                                  //
  18 //                                                                       //
  19 //                                                                       //
  20 /*
  21 Dielectron signal extraction class using functions as input.
  22
  23 A function to describe the signal as well as one to describe the background
  24 has to be deployed by the user. Alternatively on of the default implementaions
  25 can be used.
  26
  27 */
  28 //                                                                       //
  29 ///////////////////////////////////////////////////////////////////////////
  30
  31 #include <TF1.h>
  32 #include <TH1.h>
  33 #include <TGraph.h>
  34 #include <TMath.h>
  35 #include <TString.h>
  36 #include <TPaveText.h>
  37 #include <TList.h>
  38 #include <TFitResult.h>
  39 //#include <../hist/hist/src/TF1Helper.h> //not supposed to be used!
  40
  41 #include <AliLog.h>
  42
  43 #include "AliDielectronSignalFunc.h"
  44
  45 ClassImp(AliDielectronSignalFunc)
  46
  47 AliDielectronSignalFunc::AliDielectronSignalFunc() :
  48 AliDielectronSignalBase(),
  49 fFuncSignal(0x0),
  50 fFuncBackground(0x0),
  51 fFuncSigBack(0x0),
  52 fParMass(1),
  53 fParMassWidth(2),
  54 fFitOpt("SMNQE"),
  55 fUseIntegral(kFALSE)
  56 {
  57   //
  58   // Default Constructor
  59   //
  60
  61 }
  62
  63 //______________________________________________
  64 AliDielectronSignalFunc::AliDielectronSignalFunc(const char* name, const char* title) :
  65 AliDielectronSignalBase(name, title),
  66 fFuncSignal(0x0),
  67 fFuncBackground(0x0),
  68 fFuncSigBack(0x0),
  69 fParMass(1),
  70 fParMassWidth(2),
  71 fFitOpt("SMNQE"),
  72 fUseIntegral(kFALSE)
  73 {
  74   //
  75   // Named Constructor
  76   //
  77 }
  78
  79 //______________________________________________
  80 AliDielectronSignalFunc::~AliDielectronSignalFunc()
  81 {
  82   //
  83   // Default Destructor
  84   //
  85   if(fFuncSignal) delete fFuncSignal;
  86   if(fFuncBackground) delete fFuncBackground;
  87   if(fFuncSigBack) delete fFuncSigBack;
  88 }
  89
  90
  91 //______________________________________________
  92 void AliDielectronSignalFunc::Process(TObjArray * const arrhist)
  93 {
  94   //
  95   // Fit the invariant mass histograms and retrieve the signal and background
  96   //
  97   switch(fMethod) {
  98   case kFitted :
  99     ProcessFit(arrhist);
 100     break;
 101
 102   case kLikeSign :
 103     ProcessLS(arrhist);
 104     break;
 105
 106   case kEventMixing :
 107     ProcessEM(arrhist);
 108     break;
 109
 110   default :
 111     AliError("Background substraction method not known!");
 112   }
 113 }
 114
 115
 116 //______________________________________________
 117 void AliDielectronSignalFunc::ProcessFit(TObjArray * const arrhist) {
 118   //
 119   // Fit the +- invariant mass distribution only
 120   // Here we assume that the combined fit function is a sum of the signal and background functions
 121   //    and that the signal function is always the first term of this sum
 122   //
 123
 124   fHistDataPM = (TH1F*)(arrhist->At(1))->Clone("histPM");  // +-    SE
 125   fHistDataPM->Sumw2();
 126   if(fRebin>1)
 127     fHistDataPM->Rebin(fRebin);
 128
 129   fHistSignal = new TH1F("HistSignal", "Like-Sign substracted signal",
 130                          fHistDataPM->GetXaxis()->GetNbins(),
 131                          fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
 132   fHistBackground = new TH1F("HistBackground", "Like-sign contribution",
 133                              fHistDataPM->GetXaxis()->GetNbins(),
 134                              fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
 135
 136   // the starting parameters of the fit function and their limits can be tuned
 137   // by the user in its macro
 138   fHistDataPM->Fit(fFuncSigBack, fFitOpt.Data(), "", fFitMin, fFitMax);
 139   TFitResultPtr pmFitPtr = fHistDataPM->Fit(fFuncSigBack, fFitOpt.Data(), "", fFitMin, fFitMax);
 140   //TFitResult *pmFitResult = pmFitPtr.Get(); //not used when TF1Helper out
 141   fFuncSignal->SetParameters(fFuncSigBack->GetParameters());
 142   fFuncBackground->SetParameters(fFuncSigBack->GetParameters()+fFuncSignal->GetNpar());
 143
 144   for(Int_t iBin=1; iBin<=fHistDataPM->GetXaxis()->GetNbins(); iBin++) {
 145     Double_t m = fHistDataPM->GetBinCenter(iBin);
 146     Double_t pm = fHistDataPM->GetBinContent(iBin);
 147     Double_t epm = fHistDataPM->GetBinError(iBin);
 148     Double_t bknd = fFuncBackground->Eval(m);
 149     Double_t ebknd = 0;
 150     /* to be revised ... TF1Helper
 151     for(Int_t iPar=fFuncSignal->GetNpar(); iPar<fFuncSigBack->GetNpar(); iPar++) {
 152       for(Int_t jPar=iPar; jPar<fFuncSigBack->GetNpar(); jPar++) {
 153         TF1 gradientIpar("gradientIpar",
 154                          ROOT::TF1Helper::TGradientParFunction(iPar-fFuncSignal->GetNpar(),fFuncBackground),0,0,0);
 155         TF1 gradientJpar("gradientJpar",
 156                          ROOT::TF1Helper::TGradientParFunction(jPar-fFuncSignal->GetNpar(),fFuncBackground),0,0,0);
 157         ebknd += pmFitResult->CovMatrix(iPar,jPar)*
 158           gradientIpar.Eval(m)*gradientJpar.Eval(m)*
 159           (iPar==jPar ? 1.0 : 2.0);
 160       }
 161     }
 162     */ // TF1Helper
 163     Double_t signal = pm-bknd;
 164     Double_t error = TMath::Sqrt(epm*epm+ebknd);
 165     fHistSignal->SetBinContent(iBin, signal);
 166     fHistSignal->SetBinError(iBin, error);
 167     fHistBackground->SetBinContent(iBin, bknd);
 168     fHistBackground->SetBinError(iBin, TMath::Sqrt(ebknd));
 169   }
 170
 171   if(fUseIntegral) {
 172     // signal
 173     fValues(0) = fFuncSignal->Integral(fIntMin, fIntMax)/fHistDataPM->GetBinWidth(1);
 174     fErrors(0) = 0;
 175     /* to be revised ... TF1Helper
 176     for(Int_t iPar=0; iPar<fFuncSignal->GetNpar(); iPar++) {
 177       for(Int_t jPar=iPar; jPar<fFuncSignal->GetNpar(); jPar++) {
 178         TF1 gradientIpar("gradientIpar",
 179                          ROOT::TF1Helper::TGradientParFunction(iPar,fFuncSignal),0,0,0);
 180         TF1 gradientJpar("gradientJpar",
 181                          ROOT::TF1Helper::TGradientParFunction(jPar,fFuncSignal),0,0,0);
 182         fErrors(0) += pmFitResult->CovMatrix(iPar,jPar)*
 183           gradientIpar.Integral(fIntMin,fIntMax)*gradientJpar.Integral(fIntMin,fIntMax)*
 184           (iPar==jPar ? 1.0 : 2.0);
 185       }
 186     }
 187     */ //TF1Helper
 188     // background
 189     fValues(1) = fFuncBackground->Integral(fIntMin, fIntMax)/fHistDataPM->GetBinWidth(1);
 190     fErrors(1) = 0;
 191     /* to be revised... TF1Helper
 192     for(Int_t iPar=fFuncSignal->GetNpar(); iPar<fFuncSigBack->GetNpar(); iPar++) {
 193       for(Int_t jPar=iPar; jPar<fFuncSigBack->GetNpar(); jPar++) {
 194         TF1 gradientIpar("gradientIpar",
 195                          ROOT::TF1Helper::TGradientParFunction(iPar-fFuncSignal->GetNpar(),fFuncBackground),0,0,0);
 196         TF1 gradientJpar("gradientJpar",
 197                          ROOT::TF1Helper::TGradientParFunction(jPar-fFuncSignal->GetNpar(),fFuncBackground),0,0,0);
 198         fErrors(1) += pmFitResult->CovMatrix(iPar,jPar)*
 199           gradientIpar.Integral(fIntMin, fIntMax)*gradientJpar.Integral(fIntMin, fIntMax)*
 200           (iPar==jPar ? 1.0 : 2.0);
 201       }
 202     }
 203     */ // TF1Helper
 204   }
 205   else {
 206     // signal
 207     fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
 208                                                fHistSignal->FindBin(fIntMax), fErrors(0));
 209     // background
 210     fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
 211       fHistBackground->FindBin(fIntMax),
 212       fErrors(1));
 213   }
 214   // S/B and significance
 215   SetSignificanceAndSOB();
 216   fValues(4) = fFuncSigBack->GetParameter(fParMass);
 217   fErrors(4) = fFuncSigBack->GetParError(fParMass);
 218   fValues(5) = fFuncSigBack->GetParameter(fParMassWidth);
 219   fErrors(5) = fFuncSigBack->GetParError(fParMassWidth);
 220
 221   fProcessed = kTRUE;
 222 }
 223
 224 //______________________________________________
 225 void AliDielectronSignalFunc::ProcessLS(TObjArray * const arrhist) {
 226   //
 227   // Substract background using the like-sign spectrum
 228   //
 229   fHistDataPP = (TH1F*)(arrhist->At(0))->Clone("histPP");  // ++    SE
 230   fHistDataPM = (TH1F*)(arrhist->At(1))->Clone("histPM");  // +-    SE
 231   fHistDataMM = (TH1F*)(arrhist->At(2))->Clone("histMM");  // --    SE
 232   if (fRebin>1) {
 233     fHistDataPP->Rebin(fRebin);
 234     fHistDataPM->Rebin(fRebin);
 235     fHistDataMM->Rebin(fRebin);
 236   }
 237   fHistDataPP->Sumw2();
 238   fHistDataPM->Sumw2();
 239   fHistDataMM->Sumw2();
 240
 241   fHistSignal = new TH1F("HistSignal", "Like-Sign substracted signal",
 242                          fHistDataPM->GetXaxis()->GetNbins(),
 243                          fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
 244   fHistBackground = new TH1F("HistBackground", "Like-sign contribution",
 245                              fHistDataPM->GetXaxis()->GetNbins(),
 246                              fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
 247
 248   // fit the +- mass distribution
 249   fHistDataPM->Fit(fFuncSigBack, fFitOpt.Data(), "", fFitMin, fFitMax);
 250   fHistDataPM->Fit(fFuncSigBack, fFitOpt.Data(), "", fFitMin, fFitMax);
 251   // declare the variables where the like-sign fit results will be stored
 252   TFitResult *ppFitResult = 0x0;
 253   TFitResult *mmFitResult = 0x0;
 254   // fit the like sign background
 255   TF1 *funcClonePP = (TF1*)fFuncBackground->Clone("funcClonePP");
 256   TF1 *funcCloneMM = (TF1*)fFuncBackground->Clone("funcCloneMM");
 257   fHistDataPP->Fit(funcClonePP, fFitOpt.Data(), "", fFitMin, fFitMax);
 258   TFitResultPtr ppFitPtr = fHistDataPP->Fit(funcClonePP, fFitOpt.Data(), "", fFitMin, fFitMax);
 259   ppFitResult = ppFitPtr.Get();
 260   fHistDataMM->Fit(funcCloneMM, fFitOpt.Data(), "", fFitMin, fFitMax);
 261   TFitResultPtr mmFitPtr = fHistDataMM->Fit(funcCloneMM, fFitOpt.Data(), "", fFitMin, fFitMax);
 262   mmFitResult = mmFitPtr.Get();
 263
 264   /* to be revised ... TF1Helper
 265   for(Int_t iBin=1; iBin<=fHistDataPM->GetXaxis()->GetNbins(); iBin++) {
 266     Double_t m = fHistDataPM->GetBinCenter(iBin);
 267     Double_t pm = fHistDataPM->GetBinContent(iBin);
 268     Double_t pp = funcClonePP->Eval(m);
 269     Double_t mm = funcCloneMM->Eval(m);
 270     Double_t epm = fHistDataPM->GetBinError(iBin);
 271     Double_t epp = 0;
 272     for(Int_t iPar=0; iPar<funcClonePP->GetNpar(); iPar++) {
 273       for(Int_t jPar=iPar; jPar<funcClonePP->GetNpar(); jPar++) {
 274         TF1 gradientIpar("gradientIpar",
 275                          ROOT::TF1Helper::TGradientParFunction(iPar,funcClonePP),0,0,0);
 276         TF1 gradientJpar("gradientJpar",
 277                          ROOT::TF1Helper::TGradientParFunction(jPar,funcClonePP),0,0,0);
 278         epp += ppFitResult->CovMatrix(iPar,jPar)*
 279           gradientIpar.Eval(m)*gradientJpar.Eval(m)*
 280           (iPar==jPar ? 1.0 : 2.0);
 281       }
 282     }
 283
 284     Double_t emm = 0;
 285     for(Int_t iPar=0; iPar<funcCloneMM->GetNpar(); iPar++) {
 286       for(Int_t jPar=iPar; jPar<funcCloneMM->GetNpar(); jPar++) {
 287         TF1 gradientIpar("gradientIpar",
 288                          ROOT::TF1Helper::TGradientParFunction(iPar,funcCloneMM),0,0,0);
 289         TF1 gradientJpar("gradientJpar",
 290                          ROOT::TF1Helper::TGradientParFunction(jPar,funcCloneMM),0,0,0);
 291         emm += mmFitResult->CovMatrix(iPar,jPar)*
 292           gradientIpar.Eval(m)*gradientJpar.Eval(m)*
 293           (iPar==jPar ? 1.0 : 2.0);
 294       }
 295     }
 296
 297     Double_t signal = pm-2.0*TMath::Sqrt(pp*mm);
 298     Double_t background = 2.0*TMath::Sqrt(pp*mm);
 299     // error propagation on the signal calculation above
 300     Double_t esignal = TMath::Sqrt(epm*epm+(mm/pp)*epp+(pp/mm)*emm);
 301     Double_t ebackground = TMath::Sqrt((mm/pp)*epp+(pp/mm)*emm);
 302     fHistSignal->SetBinContent(iBin, signal);
 303     fHistSignal->SetBinError(iBin, esignal);
 304     fHistBackground->SetBinContent(iBin, background);
 305     fHistBackground->SetBinError(iBin, ebackground);
 306   }
 307   */  // TF1Helper
 308
 309   // signal
 310   fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
 311                                              fHistSignal->FindBin(fIntMax), fErrors(0));
 312   // background
 313   fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
 314                                                  fHistBackground->FindBin(fIntMax),
 315                                                  fErrors(1));
 316   // S/B and significance
 317   SetSignificanceAndSOB();
 318   fValues(4) = fFuncSigBack->GetParameter(fParMass);
 319   fErrors(4) = fFuncSigBack->GetParError(fParMass);
 320   fValues(5) = fFuncSigBack->GetParameter(fParMassWidth);
 321   fErrors(5) = fFuncSigBack->GetParError(fParMassWidth);
 322
 323   fProcessed = kTRUE;
 324 }
 325
 326 //______________________________________________
 327 void AliDielectronSignalFunc::ProcessEM(TObjArray * const arrhist) {
 328   //
 329   // Substract background with the event mixing technique
 330   //
 331   arrhist->GetEntries();   // just to avoid the unused parameter warning
 332   AliError("Event mixing for background substraction method not implemented!");
 333 }
 334
 335 //______________________________________________
 336 void AliDielectronSignalFunc::SetFunctions(TF1 * const combined, TF1 * const sig, TF1 * const back,
 337                                            Int_t parM, Int_t parMres)
 338 {
 339   //
 340   // Set the signal, background functions and combined fit function
 341   // Note: The process method assumes that the first n parameters in the
 342   //       combined fit function correspond to the n parameters of the signal function
 343   //       and the n+1 to n+m parameters to the m parameters of the background function!!!
 344
 345   if (!sig||!back||!combined) {
 346     AliError("Both, signal and background function need to be set!");
 347     return;
 348   }
 349   fFuncSignal=sig;
 350   fFuncBackground=back;
 351   fFuncSigBack=combined;
 352   fParMass=parM;
 353   fParMassWidth=parMres;
 354 }
 355
 356 //______________________________________________
 357 void AliDielectronSignalFunc::SetDefaults(Int_t type)
 358 {
 359   //
 360   // Setup some default functions:
 361   // type = 0: gaus signal + linear background in 2.5 - 4 GeV inv. mass
 362   // type = 1: gaus signal + exponential background in 2.5 - 4 GeV inv. mass
 363   // type = 2: half gaussian, half exponential signal function
 364   // type = 3: Crystal-Ball function
 365   // type = 4: Crystal-Ball signal + exponential background
 366   //
 367
 368   if (type==0){
 369     fFuncSignal=new TF1("DieleSignal","gaus",2.5,4);
 370     fFuncBackground=new TF1("DieleBackground","pol1",2.5,4);
 371     fFuncSigBack=new TF1("DieleCombined","gaus+pol1(3)",2.5,4);
 372
 373     fFuncSigBack->SetParameters(1,3.1,.05,2.5,1);
 374     fFuncSigBack->SetParLimits(0,0,10000000);
 375     fFuncSigBack->SetParLimits(1,3.05,3.15);
 376     fFuncSigBack->SetParLimits(2,.02,.1);
 377   }
 378   else if (type==1){
 379     fFuncSignal=new TF1("DieleSignal","gaus",2.5,4);
 380     fFuncBackground=new TF1("DieleBackground","[0]*exp(-(x-[1])/[2])",2.5,4);
 381     fFuncSigBack=new TF1("DieleCombined","gaus+[3]*exp(-(x-[4])/[5])",2.5,4);
 382
 383     fFuncSigBack->SetParameters(1,3.1,.05,1,2.5,1);
 384     fFuncSigBack->SetParLimits(0,0,10000000);
 385     fFuncSigBack->SetParLimits(1,3.05,3.15);
 386     fFuncSigBack->SetParLimits(2,.02,.1);
 387   }
 388   else if (type==2){
 389     // half gaussian, half exponential signal function
 390     // exponential background
 391     fFuncSignal = new TF1("DieleSignal","(x<[1])*([0]*(exp(-0.5*((x-[1])/[2])^2)+exp((x-[1])/[3])*(1-exp(-0.5*((x-[1])/[2])^2))))+(x>=[1])*([0]*exp(-0.5*((x-[1])/[2])^2))",2.5,4);
 392     fFuncBackground = new TF1("DieleBackground","[0]*exp(-(x-[1])/[2])+[3]",2.5,4);
 393     fFuncSigBack = new TF1("DieleCombined","(x<[1])*([0]*(exp(-0.5*((x-[1])/[2])^2)+exp((x-[1])/[3])*(1-exp(-0.5*((x-[1])/[2])^2))))+(x>=[1])*([0]*exp(-0.5*((x-[1])/[2])^2))+[4]*exp(-(x-[5])/[6])+[7]",2.5,4);
 394     fFuncSigBack->SetParameters(1.,3.1,.05,.1,1,2.5,1,0);
 395
 396     fFuncSigBack->SetParLimits(0,0,10000000);
 397     fFuncSigBack->SetParLimits(1,3.05,3.15);
 398     fFuncSigBack->SetParLimits(2,.02,.1);
 399     fFuncSigBack->FixParameter(6,2.5);
 400     fFuncSigBack->FixParameter(7,0);
 401   }
 402 }
 403
 404
 405 //______________________________________________
 406 void AliDielectronSignalFunc::Draw(const Option_t* option)
 407 {
 408   //
 409   // Draw the fitted function
 410   //
 411
 412   TString drawOpt(option);
 413   drawOpt.ToLower();
 414
 415   Bool_t optStat=drawOpt.Contains("stat");
 416
 417   fFuncSigBack->SetNpx(200);
 418   fFuncSigBack->SetRange(fIntMin,fIntMax);
 419   fFuncBackground->SetNpx(200);
 420   fFuncBackground->SetRange(fIntMin,fIntMax);
 421
 422   TGraph *grSig=new TGraph(fFuncSigBack);
 423   grSig->SetFillColor(kGreen);
 424   grSig->SetFillStyle(3001);
 425
 426   TGraph *grBack=new TGraph(fFuncBackground);
 427   grBack->SetFillColor(kRed);
 428   grBack->SetFillStyle(3001);
 429
 430   grSig->SetPoint(0,grBack->GetX()[0],grBack->GetY()[0]);
 431   grSig->SetPoint(grSig->GetN()-1,grBack->GetX()[grBack->GetN()-1],grBack->GetY()[grBack->GetN()-1]);
 432
 433   grBack->SetPoint(0,grBack->GetX()[0],0.);
 434   grBack->SetPoint(grBack->GetN()-1,grBack->GetX()[grBack->GetN()-1],0.);
 435
 436   fFuncSigBack->SetRange(fFitMin,fFitMax);
 437   fFuncBackground->SetRange(fFitMin,fFitMax);
 438
 439   if (!drawOpt.Contains("same")){
 440     if (fHistDataPM){
 441       fHistDataPM->Draw();
 442       grSig->Draw("f");
 443     } else {
 444       grSig->Draw("af");
 445     }
 446   } else {
 447     grSig->Draw("f");
 448   }
 449   if(fMethod==kFitted) grBack->Draw("f");
 450   fFuncSigBack->Draw("same");
 451   fFuncSigBack->SetLineWidth(2);
 452   if(fMethod==kLikeSign) {
 453     fHistDataPP->SetLineWidth(2);
 454     fHistDataPP->SetLineColor(6);
 455     fHistDataPP->Draw("same");
 456     fHistDataMM->SetLineWidth(2);
 457     fHistDataMM->SetLineColor(8);
 458     fHistDataMM->Draw("same");
 459   }
 460
 461   if(fMethod==kFitted)
 462     fFuncBackground->Draw("same");
 463
 464   if (optStat) DrawStats();
 465
 466 }