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