PWGDQ/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>
  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(); // used only with TF1Helper
 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     for(Int_t iPar=fFuncSignal->GetNpar(); iPar<fFuncSigBack->GetNpar(); iPar++) {
 151 /* TF1Helper problem on alien compilation
 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     }
 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     for(Int_t iPar=0; iPar<fFuncSignal->GetNpar(); iPar++) {
 176 /* TF1Helper problem on alien compilation
 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     }
 188     // background
 189     fValues(1) = fFuncBackground->Integral(fIntMin, fIntMax)/fHistDataPM->GetBinWidth(1);
 190     fErrors(1) = 0;
 191     for(Int_t iPar=fFuncSignal->GetNpar(); iPar<fFuncSigBack->GetNpar(); iPar++) {
 192 /* TF1Helper problem on alien compilation
 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     }
 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   fHistDataPP->Fit(funcClonePP, fFitOpt.Data(), "", fFitMin, fFitMax);
 259 //   TFitResultPtr ppFitPtr = fHistDataPP->Fit(funcClonePP, fFitOpt.Data(), "", fFitMin, fFitMax);
 260 //   ppFitResult = ppFitPtr.Get();
 261   fHistDataMM->Fit(funcCloneMM, fFitOpt.Data(), "", fFitMin, fFitMax);
 262   fHistDataMM->Fit(funcCloneMM, fFitOpt.Data(), "", fFitMin, fFitMax);
 263 //   TFitResultPtr mmFitPtr = fHistDataMM->Fit(funcCloneMM, fFitOpt.Data(), "", fFitMin, fFitMax);
 264 //   mmFitResult = mmFitPtr.Get();
 265
 266   for(Int_t iBin=1; iBin<=fHistDataPM->GetXaxis()->GetNbins(); iBin++) {
 267     Double_t m = fHistDataPM->GetBinCenter(iBin);
 268     Double_t pm = fHistDataPM->GetBinContent(iBin);
 269     Double_t pp = funcClonePP->Eval(m);
 270     Double_t mm = funcCloneMM->Eval(m);
 271     Double_t epm = fHistDataPM->GetBinError(iBin);
 272     Double_t epp = 0;
 273     for(Int_t iPar=0; iPar<funcClonePP->GetNpar(); iPar++) {
 274 /* TF1Helper problem on alien compilation
 275       for(Int_t jPar=iPar; jPar<funcClonePP->GetNpar(); jPar++) {
 276         TF1 gradientIpar("gradientIpar",
 277                          ROOT::TF1Helper::TGradientParFunction(iPar,funcClonePP),0,0,0);
 278         TF1 gradientJpar("gradientJpar",
 279                          ROOT::TF1Helper::TGradientParFunction(jPar,funcClonePP),0,0,0);
 280         epp += ppFitResult->CovMatrix(iPar,jPar)*
 281           gradientIpar.Eval(m)*gradientJpar.Eval(m)*
 282           (iPar==jPar ? 1.0 : 2.0);
 283       }
 284 */
 285     }
 286     Double_t emm = 0;
 287     for(Int_t iPar=0; iPar<funcCloneMM->GetNpar(); iPar++) {
 288 /* TF1Helper problem on alien compilation
 289       for(Int_t jPar=iPar; jPar<funcCloneMM->GetNpar(); jPar++) {
 290         TF1 gradientIpar("gradientIpar",
 291                          ROOT::TF1Helper::TGradientParFunction(iPar,funcCloneMM),0,0,0);
 292         TF1 gradientJpar("gradientJpar",
 293                          ROOT::TF1Helper::TGradientParFunction(jPar,funcCloneMM),0,0,0);
 294         emm += mmFitResult->CovMatrix(iPar,jPar)*
 295           gradientIpar.Eval(m)*gradientJpar.Eval(m)*
 296           (iPar==jPar ? 1.0 : 2.0);
 297       }
 298 */
 299     }
 300
 301     Double_t signal = pm-2.0*TMath::Sqrt(pp*mm);
 302     Double_t background = 2.0*TMath::Sqrt(pp*mm);
 303     // error propagation on the signal calculation above
 304     Double_t esignal = TMath::Sqrt(epm*epm+(mm/pp)*epp+(pp/mm)*emm);
 305     Double_t ebackground = TMath::Sqrt((mm/pp)*epp+(pp/mm)*emm);
 306     fHistSignal->SetBinContent(iBin, signal);
 307     fHistSignal->SetBinError(iBin, esignal);
 308     fHistBackground->SetBinContent(iBin, background);
 309     fHistBackground->SetBinError(iBin, ebackground);
 310   }
 311
 312   // signal
 313   fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
 314                                              fHistSignal->FindBin(fIntMax), fErrors(0));
 315   // background
 316   fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
 317                                                  fHistBackground->FindBin(fIntMax),
 318                                                  fErrors(1));
 319   // S/B and significance
 320   SetSignificanceAndSOB();
 321   fValues(4) = fFuncSigBack->GetParameter(fParMass);
 322   fErrors(4) = fFuncSigBack->GetParError(fParMass);
 323   fValues(5) = fFuncSigBack->GetParameter(fParMassWidth);
 324   fErrors(5) = fFuncSigBack->GetParError(fParMassWidth);
 325
 326   fProcessed = kTRUE;
 327 }
 328
 329 //______________________________________________
 330 void AliDielectronSignalFunc::ProcessEM(TObjArray * const arrhist) {
 331   //
 332   // Substract background with the event mixing technique
 333   //
 334   arrhist->GetEntries();   // just to avoid the unused parameter warning
 335   AliError("Event mixing for background substraction method not implemented!");
 336 }
 337
 338 //______________________________________________
 339 void AliDielectronSignalFunc::SetFunctions(TF1 * const combined, TF1 * const sig, TF1 * const back,
 340                                            Int_t parM, Int_t parMres)
 341 {
 342   //
 343   // Set the signal, background functions and combined fit function
 344   // Note: The process method assumes that the first n parameters in the
 345   //       combined fit function correspond to the n parameters of the signal function
 346   //       and the n+1 to n+m parameters to the m parameters of the background function!!!
 347
 348   if (!sig||!back||!combined) {
 349     AliError("Both, signal and background function need to be set!");
 350     return;
 351   }
 352   fFuncSignal=sig;
 353   fFuncBackground=back;
 354   fFuncSigBack=combined;
 355   fParMass=parM;
 356   fParMassWidth=parMres;
 357 }
 358
 359 //______________________________________________
 360 void AliDielectronSignalFunc::SetDefaults(Int_t type)
 361 {
 362   //
 363   // Setup some default functions:
 364   // type = 0: gaus signal + linear background in 2.5 - 4 GeV inv. mass
 365   // type = 1: gaus signal + exponential background in 2.5 - 4 GeV inv. mass
 366   // type = 2: half gaussian, half exponential signal function
 367   // type = 3: Crystal-Ball function
 368   // type = 4: Crystal-Ball signal + exponential background
 369   //
 370
 371   if (type==0){
 372     fFuncSignal=new TF1("DieleSignal","gaus",2.5,4);
 373     fFuncBackground=new TF1("DieleBackground","pol1",2.5,4);
 374     fFuncSigBack=new TF1("DieleCombined","gaus+pol1(3)",2.5,4);
 375
 376     fFuncSigBack->SetParameters(1,3.1,.05,2.5,1);
 377     fFuncSigBack->SetParLimits(0,0,10000000);
 378     fFuncSigBack->SetParLimits(1,3.05,3.15);
 379     fFuncSigBack->SetParLimits(2,.02,.1);
 380   }
 381   else if (type==1){
 382     fFuncSignal=new TF1("DieleSignal","gaus",2.5,4);
 383     fFuncBackground=new TF1("DieleBackground","[0]*exp(-(x-[1])/[2])",2.5,4);
 384     fFuncSigBack=new TF1("DieleCombined","gaus+[3]*exp(-(x-[4])/[5])",2.5,4);
 385
 386     fFuncSigBack->SetParameters(1,3.1,.05,1,2.5,1);
 387     fFuncSigBack->SetParLimits(0,0,10000000);
 388     fFuncSigBack->SetParLimits(1,3.05,3.15);
 389     fFuncSigBack->SetParLimits(2,.02,.1);
 390   }
 391   else if (type==2){
 392     // half gaussian, half exponential signal function
 393     // exponential background
 394     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);
 395     fFuncBackground = new TF1("DieleBackground","[0]*exp(-(x-[1])/[2])+[3]",2.5,4);
 396     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);
 397     fFuncSigBack->SetParameters(1.,3.1,.05,.1,1,2.5,1,0);
 398
 399     fFuncSigBack->SetParLimits(0,0,10000000);
 400     fFuncSigBack->SetParLimits(1,3.05,3.15);
 401     fFuncSigBack->SetParLimits(2,.02,.1);
 402     fFuncSigBack->FixParameter(6,2.5);
 403     fFuncSigBack->FixParameter(7,0);
 404   }
 405 }
 406
 407
 408 //______________________________________________
 409 void AliDielectronSignalFunc::Draw(const Option_t* option)
 410 {
 411   //
 412   // Draw the fitted function
 413   //
 414
 415   TString drawOpt(option);
 416   drawOpt.ToLower();
 417
 418   Bool_t optStat=drawOpt.Contains("stat");
 419
 420   fFuncSigBack->SetNpx(200);
 421   fFuncSigBack->SetRange(fIntMin,fIntMax);
 422   fFuncBackground->SetNpx(200);
 423   fFuncBackground->SetRange(fIntMin,fIntMax);
 424
 425   TGraph *grSig=new TGraph(fFuncSigBack);
 426   grSig->SetFillColor(kGreen);
 427   grSig->SetFillStyle(3001);
 428
 429   TGraph *grBack=new TGraph(fFuncBackground);
 430   grBack->SetFillColor(kRed);
 431   grBack->SetFillStyle(3001);
 432
 433   grSig->SetPoint(0,grBack->GetX()[0],grBack->GetY()[0]);
 434   grSig->SetPoint(grSig->GetN()-1,grBack->GetX()[grBack->GetN()-1],grBack->GetY()[grBack->GetN()-1]);
 435
 436   grBack->SetPoint(0,grBack->GetX()[0],0.);
 437   grBack->SetPoint(grBack->GetN()-1,grBack->GetX()[grBack->GetN()-1],0.);
 438
 439   fFuncSigBack->SetRange(fFitMin,fFitMax);
 440   fFuncBackground->SetRange(fFitMin,fFitMax);
 441
 442   if (!drawOpt.Contains("same")){
 443     if (fHistDataPM){
 444       fHistDataPM->Draw();
 445       grSig->Draw("f");
 446     } else {
 447       grSig->Draw("af");
 448     }
 449   } else {
 450     grSig->Draw("f");
 451   }
 452   if(fMethod==kFitted) grBack->Draw("f");
 453   fFuncSigBack->Draw("same");
 454   fFuncSigBack->SetLineWidth(2);
 455   if(fMethod==kLikeSign) {
 456     fHistDataPP->SetLineWidth(2);
 457     fHistDataPP->SetLineColor(6);
 458     fHistDataPP->Draw("same");
 459     fHistDataMM->SetLineWidth(2);
 460     fHistDataMM->SetLineColor(8);
 461     fHistDataMM->Draw("same");
 462   }
 463
 464   if(fMethod==kFitted)
 465     fFuncBackground->Draw("same");
 466
 467   if (optStat) DrawStats();
 468
 469 }