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>
  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();
 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       for(Int_t jPar=iPar; jPar<fFuncSigBack->GetNpar(); jPar++) {
 152         TF1 gradientIpar("gradientIpar",
 153                          ROOT::TF1Helper::TGradientParFunction(iPar-fFuncSignal->GetNpar(),fFuncBackground),0,0,0);
 154         TF1 gradientJpar("gradientJpar",
 155                          ROOT::TF1Helper::TGradientParFunction(jPar-fFuncSignal->GetNpar(),fFuncBackground),0,0,0);
 156         ebknd += pmFitResult->CovMatrix(iPar,jPar)*
 157           gradientIpar.Eval(m)*gradientJpar.Eval(m)*
 158           (iPar==jPar ? 1.0 : 2.0);
 159       }
 160     }
 161     Double_t signal = pm-bknd;
 162     Double_t error = TMath::Sqrt(epm*epm+ebknd);
 163     fHistSignal->SetBinContent(iBin, signal);
 164     fHistSignal->SetBinError(iBin, error);
 165     fHistBackground->SetBinContent(iBin, bknd);
 166     fHistBackground->SetBinError(iBin, TMath::Sqrt(ebknd));
 167   }
 168
 169   if(fUseIntegral) {
 170     // signal
 171     fValues(0) = fFuncSignal->Integral(fIntMin, fIntMax)/fHistDataPM->GetBinWidth(1);
 172     fErrors(0) = 0;
 173     for(Int_t iPar=0; iPar<fFuncSignal->GetNpar(); iPar++) {
 174       for(Int_t jPar=iPar; jPar<fFuncSignal->GetNpar(); jPar++) {
 175         TF1 gradientIpar("gradientIpar",
 176                          ROOT::TF1Helper::TGradientParFunction(iPar,fFuncSignal),0,0,0);
 177         TF1 gradientJpar("gradientJpar",
 178                          ROOT::TF1Helper::TGradientParFunction(jPar,fFuncSignal),0,0,0);
 179         fErrors(0) += pmFitResult->CovMatrix(iPar,jPar)*
 180           gradientIpar.Integral(fIntMin,fIntMax)*gradientJpar.Integral(fIntMin,fIntMax)*
 181           (iPar==jPar ? 1.0 : 2.0);
 182       }
 183     }
 184     // background
 185     fValues(1) = fFuncBackground->Integral(fIntMin, fIntMax)/fHistDataPM->GetBinWidth(1);
 186     fErrors(1) = 0;
 187     for(Int_t iPar=fFuncSignal->GetNpar(); iPar<fFuncSigBack->GetNpar(); iPar++) {
 188       for(Int_t jPar=iPar; jPar<fFuncSigBack->GetNpar(); jPar++) {
 189         TF1 gradientIpar("gradientIpar",
 190                          ROOT::TF1Helper::TGradientParFunction(iPar-fFuncSignal->GetNpar(),fFuncBackground),0,0,0);
 191         TF1 gradientJpar("gradientJpar",
 192                          ROOT::TF1Helper::TGradientParFunction(jPar-fFuncSignal->GetNpar(),fFuncBackground),0,0,0);
 193         fErrors(1) += pmFitResult->CovMatrix(iPar,jPar)*
 194           gradientIpar.Integral(fIntMin, fIntMax)*gradientJpar.Integral(fIntMin, fIntMax)*
 195           (iPar==jPar ? 1.0 : 2.0);
 196       }
 197     }
 198   }
 199   else {
 200     // signal
 201     fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
 202                                                fHistSignal->FindBin(fIntMax), fErrors(0));
 203     // background
 204     fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
 205       fHistBackground->FindBin(fIntMax),
 206       fErrors(1));
 207   }
 208   // S/B and significance
 209   SetSignificanceAndSOB();
 210   fValues(4) = fFuncSigBack->GetParameter(fParMass);
 211   fErrors(4) = fFuncSigBack->GetParError(fParMass);
 212   fValues(5) = fFuncSigBack->GetParameter(fParMassWidth);
 213   fErrors(5) = fFuncSigBack->GetParError(fParMassWidth);
 214
 215   fProcessed = kTRUE;
 216 }
 217
 218 //______________________________________________
 219 void AliDielectronSignalFunc::ProcessLS(TObjArray * const arrhist) {
 220   //
 221   // Substract background using the like-sign spectrum
 222   //
 223   fHistDataPP = (TH1F*)(arrhist->At(0))->Clone("histPP");  // ++    SE
 224   fHistDataPM = (TH1F*)(arrhist->At(1))->Clone("histPM");  // +-    SE
 225   fHistDataMM = (TH1F*)(arrhist->At(2))->Clone("histMM");  // --    SE
 226   if (fRebin>1) {
 227     fHistDataPP->Rebin(fRebin);
 228     fHistDataPM->Rebin(fRebin);
 229     fHistDataMM->Rebin(fRebin);
 230   }
 231   fHistDataPP->Sumw2();
 232   fHistDataPM->Sumw2();
 233   fHistDataMM->Sumw2();
 234
 235   fHistSignal = new TH1F("HistSignal", "Like-Sign substracted signal",
 236                          fHistDataPM->GetXaxis()->GetNbins(),
 237                          fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
 238   fHistBackground = new TH1F("HistBackground", "Like-sign contribution",
 239                              fHistDataPM->GetXaxis()->GetNbins(),
 240                              fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
 241
 242   // fit the +- mass distribution
 243   fHistDataPM->Fit(fFuncSigBack, fFitOpt.Data(), "", fFitMin, fFitMax);
 244   fHistDataPM->Fit(fFuncSigBack, fFitOpt.Data(), "", fFitMin, fFitMax);
 245   // declare the variables where the like-sign fit results will be stored
 246   TFitResult *ppFitResult = 0x0;
 247   TFitResult *mmFitResult = 0x0;
 248   // fit the like sign background
 249   TF1 *funcClonePP = (TF1*)fFuncBackground->Clone("funcClonePP");
 250   TF1 *funcCloneMM = (TF1*)fFuncBackground->Clone("funcCloneMM");
 251   fHistDataPP->Fit(funcClonePP, fFitOpt.Data(), "", fFitMin, fFitMax);
 252   TFitResultPtr ppFitPtr = fHistDataPP->Fit(funcClonePP, fFitOpt.Data(), "", fFitMin, fFitMax);
 253   ppFitResult = ppFitPtr.Get();
 254   fHistDataMM->Fit(funcCloneMM, fFitOpt.Data(), "", fFitMin, fFitMax);
 255   TFitResultPtr mmFitPtr = fHistDataMM->Fit(funcCloneMM, fFitOpt.Data(), "", fFitMin, fFitMax);
 256   mmFitResult = mmFitPtr.Get();
 257
 258   for(Int_t iBin=1; iBin<=fHistDataPM->GetXaxis()->GetNbins(); iBin++) {
 259     Double_t m = fHistDataPM->GetBinCenter(iBin);
 260     Double_t pm = fHistDataPM->GetBinContent(iBin);
 261     Double_t pp = funcClonePP->Eval(m);
 262     Double_t mm = funcCloneMM->Eval(m);
 263     Double_t epm = fHistDataPM->GetBinError(iBin);
 264     Double_t epp = 0;
 265     for(Int_t iPar=0; iPar<funcClonePP->GetNpar(); iPar++) {
 266       for(Int_t jPar=iPar; jPar<funcClonePP->GetNpar(); jPar++) {
 267         TF1 gradientIpar("gradientIpar",
 268                          ROOT::TF1Helper::TGradientParFunction(iPar,funcClonePP),0,0,0);
 269         TF1 gradientJpar("gradientJpar",
 270                          ROOT::TF1Helper::TGradientParFunction(jPar,funcClonePP),0,0,0);
 271         epp += ppFitResult->CovMatrix(iPar,jPar)*
 272           gradientIpar.Eval(m)*gradientJpar.Eval(m)*
 273           (iPar==jPar ? 1.0 : 2.0);
 274       }
 275     }
 276     Double_t emm = 0;
 277     for(Int_t iPar=0; iPar<funcCloneMM->GetNpar(); iPar++) {
 278       for(Int_t jPar=iPar; jPar<funcCloneMM->GetNpar(); jPar++) {
 279         TF1 gradientIpar("gradientIpar",
 280                          ROOT::TF1Helper::TGradientParFunction(iPar,funcCloneMM),0,0,0);
 281         TF1 gradientJpar("gradientJpar",
 282                          ROOT::TF1Helper::TGradientParFunction(jPar,funcCloneMM),0,0,0);
 283         emm += mmFitResult->CovMatrix(iPar,jPar)*
 284           gradientIpar.Eval(m)*gradientJpar.Eval(m)*
 285           (iPar==jPar ? 1.0 : 2.0);
 286       }
 287     }
 288
 289     Double_t signal = pm-2.0*TMath::Sqrt(pp*mm);
 290     Double_t background = 2.0*TMath::Sqrt(pp*mm);
 291     // error propagation on the signal calculation above
 292     Double_t esignal = TMath::Sqrt(epm*epm+(mm/pp)*epp+(pp/mm)*emm);
 293     Double_t ebackground = TMath::Sqrt((mm/pp)*epp+(pp/mm)*emm);
 294     fHistSignal->SetBinContent(iBin, signal);
 295     fHistSignal->SetBinError(iBin, esignal);
 296     fHistBackground->SetBinContent(iBin, background);
 297     fHistBackground->SetBinError(iBin, ebackground);
 298   }
 299
 300   // signal
 301   fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
 302                                              fHistSignal->FindBin(fIntMax), fErrors(0));
 303   // background
 304   fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
 305                                                  fHistBackground->FindBin(fIntMax),
 306                                                  fErrors(1));
 307   // S/B and significance
 308   SetSignificanceAndSOB();
 309   fValues(4) = fFuncSigBack->GetParameter(fParMass);
 310   fErrors(4) = fFuncSigBack->GetParError(fParMass);
 311   fValues(5) = fFuncSigBack->GetParameter(fParMassWidth);
 312   fErrors(5) = fFuncSigBack->GetParError(fParMassWidth);
 313
 314   fProcessed = kTRUE;
 315 }
 316
 317 //______________________________________________
 318 void AliDielectronSignalFunc::ProcessEM(TObjArray * const arrhist) {
 319   //
 320   // Substract background with the event mixing technique
 321   //
 322   arrhist->GetEntries();   // just to avoid the unused parameter warning
 323   AliError("Event mixing for background substraction method not implemented!");
 324 }
 325
 326 //______________________________________________
 327 void AliDielectronSignalFunc::SetFunctions(TF1 * const combined, TF1 * const sig, TF1 * const back,
 328                                            Int_t parM, Int_t parMres)
 329 {
 330   //
 331   // Set the signal, background functions and combined fit function
 332   // Note: The process method assumes that the first n parameters in the
 333   //       combined fit function correspond to the n parameters of the signal function
 334   //       and the n+1 to n+m parameters to the m parameters of the background function!!!
 335
 336   if (!sig||!back||!combined) {
 337     AliError("Both, signal and background function need to be set!");
 338     return;
 339   }
 340   fFuncSignal=sig;
 341   fFuncBackground=back;
 342   fFuncSigBack=combined;
 343   fParMass=parM;
 344   fParMassWidth=parMres;
 345 }
 346
 347 //______________________________________________
 348 void AliDielectronSignalFunc::SetDefaults(Int_t type)
 349 {
 350   //
 351   // Setup some default functions:
 352   // type = 0: gaus signal + linear background in 2.5 - 4 GeV inv. mass
 353   // type = 1: gaus signal + exponential background in 2.5 - 4 GeV inv. mass
 354   // type = 2: half gaussian, half exponential signal function
 355   // type = 3: Crystal-Ball function
 356   // type = 4: Crystal-Ball signal + exponential background
 357   //
 358
 359   if (type==0){
 360     fFuncSignal=new TF1("DieleSignal","gaus",2.5,4);
 361     fFuncBackground=new TF1("DieleBackground","pol1",2.5,4);
 362     fFuncSigBack=new TF1("DieleCombined","gaus+pol1(3)",2.5,4);
 363
 364     fFuncSigBack->SetParameters(1,3.1,.05,2.5,1);
 365     fFuncSigBack->SetParLimits(0,0,10000000);
 366     fFuncSigBack->SetParLimits(1,3.05,3.15);
 367     fFuncSigBack->SetParLimits(2,.02,.1);
 368   }
 369   else if (type==1){
 370     fFuncSignal=new TF1("DieleSignal","gaus",2.5,4);
 371     fFuncBackground=new TF1("DieleBackground","[0]*exp(-(x-[1])/[2])",2.5,4);
 372     fFuncSigBack=new TF1("DieleCombined","gaus+[3]*exp(-(x-[4])/[5])",2.5,4);
 373
 374     fFuncSigBack->SetParameters(1,3.1,.05,1,2.5,1);
 375     fFuncSigBack->SetParLimits(0,0,10000000);
 376     fFuncSigBack->SetParLimits(1,3.05,3.15);
 377     fFuncSigBack->SetParLimits(2,.02,.1);
 378   }
 379   else if (type==2){
 380     // half gaussian, half exponential signal function
 381     // exponential background
 382     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);
 383     fFuncBackground = new TF1("DieleBackground","[0]*exp(-(x-[1])/[2])+[3]",2.5,4);
 384     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);
 385     fFuncSigBack->SetParameters(1.,3.1,.05,.1,1,2.5,1,0);
 386
 387     fFuncSigBack->SetParLimits(0,0,10000000);
 388     fFuncSigBack->SetParLimits(1,3.05,3.15);
 389     fFuncSigBack->SetParLimits(2,.02,.1);
 390     fFuncSigBack->FixParameter(6,2.5);
 391     fFuncSigBack->FixParameter(7,0);
 392   }
 393 }
 394
 395
 396 //______________________________________________
 397 void AliDielectronSignalFunc::Draw(const Option_t* option)
 398 {
 399   //
 400   // Draw the fitted function
 401   //
 402
 403   TString drawOpt(option);
 404   drawOpt.ToLower();
 405
 406   Bool_t optStat=drawOpt.Contains("stat");
 407
 408   fFuncSigBack->SetNpx(200);
 409   fFuncSigBack->SetRange(fIntMin,fIntMax);
 410   fFuncBackground->SetNpx(200);
 411   fFuncBackground->SetRange(fIntMin,fIntMax);
 412
 413   TGraph *grSig=new TGraph(fFuncSigBack);
 414   grSig->SetFillColor(kGreen);
 415   grSig->SetFillStyle(3001);
 416
 417   TGraph *grBack=new TGraph(fFuncBackground);
 418   grBack->SetFillColor(kRed);
 419   grBack->SetFillStyle(3001);
 420
 421   grSig->SetPoint(0,grBack->GetX()[0],grBack->GetY()[0]);
 422   grSig->SetPoint(grSig->GetN()-1,grBack->GetX()[grBack->GetN()-1],grBack->GetY()[grBack->GetN()-1]);
 423
 424   grBack->SetPoint(0,grBack->GetX()[0],0.);
 425   grBack->SetPoint(grBack->GetN()-1,grBack->GetX()[grBack->GetN()-1],0.);
 426
 427   fFuncSigBack->SetRange(fFitMin,fFitMax);
 428   fFuncBackground->SetRange(fFitMin,fFitMax);
 429
 430   if (!drawOpt.Contains("same")){
 431     if (fHistDataPM){
 432       fHistDataPM->Draw();
 433       grSig->Draw("f");
 434     } else {
 435       grSig->Draw("af");
 436     }
 437   } else {
 438     grSig->Draw("f");
 439   }
 440   if(fMethod==kFitted) grBack->Draw("f");
 441   fFuncSigBack->Draw("same");
 442   fFuncSigBack->SetLineWidth(2);
 443   if(fMethod==kLikeSign) {
 444     fHistDataPP->SetLineWidth(2);
 445     fHistDataPP->SetLineColor(6);
 446     fHistDataPP->Draw("same");
 447     fHistDataMM->SetLineWidth(2);
 448     fHistDataMM->SetLineColor(8);
 449     fHistDataMM->Draw("same");
 450   }
 451
 452   if(fMethod==kFitted)
 453     fFuncBackground->Draw("same");
 454
 455   if (optStat) DrawStats();
 456
 457 }