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 TH1F* AliDielectronSignalFunc::fgHistSimPM=0x0;
  48 TF1*  AliDielectronSignalFunc::fFuncSignal=0x0;
  49 TF1*  AliDielectronSignalFunc::fFuncBackground=0x0;
  50 Int_t AliDielectronSignalFunc::fNparPeak=0;
  51 Int_t AliDielectronSignalFunc::fNparBgnd=0;
  52
  53
  54 AliDielectronSignalFunc::AliDielectronSignalFunc() :
  55 AliDielectronSignalBase(),
  56 fFuncSigBack(0x0),
  57 fParMass(1),
  58 fParMassWidth(2),
  59 fFitOpt("SMNQE"),
  60 fUseIntegral(kFALSE),
  61 fPolDeg(0),
  62 fDof(0),
  63 fChi2Dof(0.0)
  64 {
  65   //
  66   // Default Constructor
  67   //
  68
  69 }
  70
  71 //______________________________________________
  72 AliDielectronSignalFunc::AliDielectronSignalFunc(const char* name, const char* title) :
  73 AliDielectronSignalBase(name, title),
  74 fFuncSigBack(0x0),
  75 fParMass(1),
  76 fParMassWidth(2),
  77 fFitOpt("SMNQE"),
  78 fUseIntegral(kFALSE),
  79 fPolDeg(0),
  80 fDof(0),
  81 fChi2Dof(0.0)
  82 {
  83   //
  84   // Named Constructor
  85   //
  86 }
  87
  88 //______________________________________________
  89 AliDielectronSignalFunc::~AliDielectronSignalFunc()
  90 {
  91   //
  92   // Default Destructor
  93   //
  94   if(fFuncSigBack) delete fFuncSigBack;
  95 }
  96
  97
  98 //______________________________________________
  99 void AliDielectronSignalFunc::Process(TObjArray * const arrhist)
 100 {
 101   //
 102   // Fit the invariant mass histograms and retrieve the signal and background
 103   //
 104   switch(fMethod) {
 105   case kFitted :
 106     ProcessFit(arrhist);
 107     break;
 108
 109   case kLikeSign :
 110     ProcessLS(arrhist);
 111     break;
 112
 113   case kEventMixing :
 114     ProcessEM(arrhist);
 115     break;
 116
 117   default :
 118     AliError("Background substraction method not known!");
 119   }
 120 }
 121 //______________________________________________________________________________
 122 Double_t AliDielectronSignalFunc::PeakFunMC(const Double_t *x, const Double_t *par) {
 123   // Fit MC signal shape
 124   // parameters
 125   // [0]:   scale for simpeak
 126
 127   Double_t xx  = x[0];
 128
 129   TH1F *hPeak = fgHistSimPM;
 130   if (!hPeak) {
 131     printf("E-AliDielectronSignalFunc::PeakFun: No histogram for peak fit defined!\n");
 132     return 0.0;
 133   }
 134
 135   Int_t idx = hPeak->FindBin(xx);
 136   if ((idx <= 1) ||
 137       (idx >= hPeak->GetNbinsX())) {
 138     return 0.0;
 139   }
 140
 141   return (par[0+fNparBgnd] * hPeak->GetBinContent(idx));
 142
 143 }
 144
 145 //______________________________________________________________________________
 146 Double_t AliDielectronSignalFunc::PeakFunCB(const Double_t *x, const Double_t *par) {
 147   // Crystal Ball fit function
 148
 149   Double_t alpha = par[0+fNparBgnd];
 150   Double_t     n = par[1+fNparBgnd];
 151   Double_t meanx = par[2+fNparBgnd];
 152   Double_t sigma = par[3+fNparBgnd];
 153   Double_t    nn = par[4+fNparBgnd];
 154
 155   Double_t a = TMath::Power((n/TMath::Abs(alpha)), n) * TMath::Exp(-.5*alpha*alpha);
 156   Double_t b = n/TMath::Abs(alpha) - TMath::Abs(alpha);
 157
 158   Double_t arg = (x[0] - meanx)/sigma;
 159   Double_t fitval = 0;
 160
 161   if (arg > -1.*alpha) {
 162     fitval = nn * TMath::Exp(-.5*arg*arg);
 163   } else {
 164     fitval = nn * a * TMath::Power((b-arg), (-1*n));
 165   }
 166
 167   return fitval;
 168 }
 169
 170 //______________________________________________________________________________
 171 Double_t AliDielectronSignalFunc::PeakFunGaus(const Double_t *x, const Double_t *par) {
 172   // Gaussian fit function
 173
 174   Double_t     n = par[0+fNparBgnd];
 175   Double_t  mean = par[1+fNparBgnd];
 176   Double_t sigma = par[2+fNparBgnd];
 177   Double_t    xx = x[0];
 178
 179   return ( n*exp(-0.5*TMath::Power((xx-mean)/sigma,2)) );
 180 }
 181
 182 //______________________________________________
 183 void AliDielectronSignalFunc::ProcessFit(TObjArray * const arrhist) {
 184   //
 185   // Fit the +- invariant mass distribution only
 186   // Here we assume that the combined fit function is a sum of the signal and background functions
 187   //    and that the signal function is always the first term of this sum
 188   //
 189
 190   fHistDataPM = (TH1F*)(arrhist->At(1))->Clone("histPM");  // +-    SE
 191   fHistDataPM->Sumw2();
 192   if(fRebin>1)
 193     fHistDataPM->Rebin(fRebin);
 194
 195   fHistSignal = new TH1F("HistSignal", "Fit substracted signal",
 196                          fHistDataPM->GetXaxis()->GetNbins(),
 197                          fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
 198   fHistBackground = new TH1F("HistBackground", "Fit contribution",
 199                              fHistDataPM->GetXaxis()->GetNbins(),
 200                              fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
 201
 202   // the starting parameters of the fit function and their limits can be tuned
 203   // by the user in its macro
 204   fHistDataPM->Fit(fFuncSigBack, fFitOpt.Data(), "", fFitMin, fFitMax);
 205   TFitResultPtr pmFitPtr = fHistDataPM->Fit(fFuncSigBack, fFitOpt.Data(), "", fFitMin, fFitMax);
 206   //TFitResult *pmFitResult = pmFitPtr.Get(); // used only with TF1Helper
 207   fFuncBackground->SetParameters(fFuncSigBack->GetParameters());
 208   fFuncSignal->SetParameters(fFuncSigBack->GetParameters()+fFuncBackground->GetNpar());
 209   //  fFuncBackground->SetParameters(fFuncSigBack->GetParameters()+fFuncSignal->GetNpar());
 210
 211   for(Int_t iBin=1; iBin<=fHistDataPM->GetXaxis()->GetNbins(); iBin++) {
 212     Double_t m = fHistDataPM->GetBinCenter(iBin);
 213     Double_t pm = fHistDataPM->GetBinContent(iBin);
 214     Double_t epm = fHistDataPM->GetBinError(iBin);
 215     Double_t bknd = fFuncBackground->Eval(m);
 216     Double_t ebknd = 0;
 217     for(Int_t iPar=fFuncSignal->GetNpar(); iPar<fFuncSigBack->GetNpar(); iPar++) {
 218 /* TF1Helper problem on alien compilation
 219       for(Int_t jPar=iPar; jPar<fFuncSigBack->GetNpar(); jPar++) {
 220         TF1 gradientIpar("gradientIpar",
 221                          ROOT::TF1Helper::TGradientParFunction(iPar-fFuncSignal->GetNpar(),fFuncBackground),0,0,0);
 222         TF1 gradientJpar("gradientJpar",
 223                          ROOT::TF1Helper::TGradientParFunction(jPar-fFuncSignal->GetNpar(),fFuncBackground),0,0,0);
 224         ebknd += pmFitResult->CovMatrix(iPar,jPar)*
 225           gradientIpar.Eval(m)*gradientJpar.Eval(m)*
 226           (iPar==jPar ? 1.0 : 2.0);
 227       }
 228 */
 229     }
 230     Double_t signal = pm-bknd;
 231     Double_t error = TMath::Sqrt(epm*epm+ebknd);
 232     fHistSignal->SetBinContent(iBin, signal);
 233     fHistSignal->SetBinError(iBin, error);
 234     fHistBackground->SetBinContent(iBin, bknd);
 235     fHistBackground->SetBinError(iBin, TMath::Sqrt(ebknd));
 236   }
 237
 238   if(fUseIntegral) {
 239     // signal
 240     fValues(0) = fFuncSignal->Integral(fIntMin, fIntMax)/fHistDataPM->GetBinWidth(1);
 241     fErrors(0) = 0;
 242     for(Int_t iPar=0; iPar<fFuncSignal->GetNpar(); iPar++) {
 243 /* TF1Helper problem on alien compilation
 244       for(Int_t jPar=iPar; jPar<fFuncSignal->GetNpar(); jPar++) {
 245         TF1 gradientIpar("gradientIpar",
 246                          ROOT::TF1Helper::TGradientParFunction(iPar,fFuncSignal),0,0,0);
 247         TF1 gradientJpar("gradientJpar",
 248                          ROOT::TF1Helper::TGradientParFunction(jPar,fFuncSignal),0,0,0);
 249         fErrors(0) += pmFitResult->CovMatrix(iPar,jPar)*
 250           gradientIpar.Integral(fIntMin,fIntMax)*gradientJpar.Integral(fIntMin,fIntMax)*
 251           (iPar==jPar ? 1.0 : 2.0);
 252       }
 253 */
 254     }
 255     // background
 256     fValues(1) = fFuncBackground->Integral(fIntMin, fIntMax)/fHistDataPM->GetBinWidth(1);
 257     fErrors(1) = 0;
 258     for(Int_t iPar=fFuncSignal->GetNpar(); iPar<fFuncSigBack->GetNpar(); iPar++) {
 259 /* TF1Helper problem on alien compilation
 260       for(Int_t jPar=iPar; jPar<fFuncSigBack->GetNpar(); jPar++) {
 261         TF1 gradientIpar("gradientIpar",
 262                          ROOT::TF1Helper::TGradientParFunction(iPar-fFuncSignal->GetNpar(),fFuncBackground),0,0,0);
 263         TF1 gradientJpar("gradientJpar",
 264                          ROOT::TF1Helper::TGradientParFunction(jPar-fFuncSignal->GetNpar(),fFuncBackground),0,0,0);
 265         fErrors(1) += pmFitResult->CovMatrix(iPar,jPar)*
 266           gradientIpar.Integral(fIntMin, fIntMax)*gradientJpar.Integral(fIntMin, fIntMax)*
 267           (iPar==jPar ? 1.0 : 2.0);
 268       }
 269 */
 270     }
 271   }
 272   else {
 273     // signal
 274     fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
 275                                                fHistSignal->FindBin(fIntMax), fErrors(0));
 276     // background
 277     fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
 278       fHistBackground->FindBin(fIntMax),
 279       fErrors(1));
 280   }
 281   // S/B and significance
 282   SetSignificanceAndSOB();
 283   fValues(4) = fFuncSigBack->GetParameter(fParMass);
 284   fErrors(4) = fFuncSigBack->GetParError(fParMass);
 285   fValues(5) = fFuncSigBack->GetParameter(fParMassWidth);
 286   fErrors(5) = fFuncSigBack->GetParError(fParMassWidth);
 287
 288   fProcessed = kTRUE;
 289
 290   fHistBackground->GetListOfFunctions()->Add(fFuncBackground);
 291
 292 }
 293
 294 //______________________________________________
 295 void AliDielectronSignalFunc::ProcessLS(TObjArray * const arrhist) {
 296   //
 297   // Substract background using the like-sign spectrum
 298   //
 299   fHistDataPP = (TH1F*)(arrhist->At(0))->Clone("histPP");  // ++    SE
 300   fHistDataPM = (TH1F*)(arrhist->At(1))->Clone("histPM");  // +-    SE
 301   fHistDataMM = (TH1F*)(arrhist->At(2))->Clone("histMM");  // --    SE
 302   if (fRebin>1) {
 303     fHistDataPP->Rebin(fRebin);
 304     fHistDataPM->Rebin(fRebin);
 305     fHistDataMM->Rebin(fRebin);
 306   }
 307   fHistDataPP->Sumw2();
 308   fHistDataPM->Sumw2();
 309   fHistDataMM->Sumw2();
 310
 311   fHistSignal = new TH1F("HistSignal", "Like-Sign substracted signal",
 312                          fHistDataPM->GetXaxis()->GetNbins(),
 313                          fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
 314   fHistBackground = new TH1F("HistBackground", "Like-sign contribution",
 315                              fHistDataPM->GetXaxis()->GetNbins(),
 316                              fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
 317
 318   // fit the +- mass distribution
 319   fHistDataPM->Fit(fFuncSigBack, fFitOpt.Data(), "", fFitMin, fFitMax);
 320   fHistDataPM->Fit(fFuncSigBack, fFitOpt.Data(), "", fFitMin, fFitMax);
 321   // declare the variables where the like-sign fit results will be stored
 322 //   TFitResult *ppFitResult = 0x0;
 323 //   TFitResult *mmFitResult = 0x0;
 324   // fit the like sign background
 325   TF1 *funcClonePP = (TF1*)fFuncBackground->Clone("funcClonePP");
 326   TF1 *funcCloneMM = (TF1*)fFuncBackground->Clone("funcCloneMM");
 327   fHistDataPP->Fit(funcClonePP, fFitOpt.Data(), "", fFitMin, fFitMax);
 328   fHistDataPP->Fit(funcClonePP, fFitOpt.Data(), "", fFitMin, fFitMax);
 329 //   TFitResultPtr ppFitPtr = fHistDataPP->Fit(funcClonePP, fFitOpt.Data(), "", fFitMin, fFitMax);
 330 //   ppFitResult = ppFitPtr.Get();
 331   fHistDataMM->Fit(funcCloneMM, fFitOpt.Data(), "", fFitMin, fFitMax);
 332   fHistDataMM->Fit(funcCloneMM, fFitOpt.Data(), "", fFitMin, fFitMax);
 333 //   TFitResultPtr mmFitPtr = fHistDataMM->Fit(funcCloneMM, fFitOpt.Data(), "", fFitMin, fFitMax);
 334 //   mmFitResult = mmFitPtr.Get();
 335
 336   for(Int_t iBin=1; iBin<=fHistDataPM->GetXaxis()->GetNbins(); iBin++) {
 337     Double_t m = fHistDataPM->GetBinCenter(iBin);
 338     Double_t pm = fHistDataPM->GetBinContent(iBin);
 339     Double_t pp = funcClonePP->Eval(m);
 340     Double_t mm = funcCloneMM->Eval(m);
 341     Double_t epm = fHistDataPM->GetBinError(iBin);
 342     Double_t epp = 0;
 343     for(Int_t iPar=0; iPar<funcClonePP->GetNpar(); iPar++) {
 344 /* TF1Helper problem on alien compilation
 345       for(Int_t jPar=iPar; jPar<funcClonePP->GetNpar(); jPar++) {
 346         TF1 gradientIpar("gradientIpar",
 347                          ROOT::TF1Helper::TGradientParFunction(iPar,funcClonePP),0,0,0);
 348         TF1 gradientJpar("gradientJpar",
 349                          ROOT::TF1Helper::TGradientParFunction(jPar,funcClonePP),0,0,0);
 350         epp += ppFitResult->CovMatrix(iPar,jPar)*
 351           gradientIpar.Eval(m)*gradientJpar.Eval(m)*
 352           (iPar==jPar ? 1.0 : 2.0);
 353       }
 354 */
 355     }
 356     Double_t emm = 0;
 357     for(Int_t iPar=0; iPar<funcCloneMM->GetNpar(); iPar++) {
 358 /* TF1Helper problem on alien compilation
 359       for(Int_t jPar=iPar; jPar<funcCloneMM->GetNpar(); jPar++) {
 360         TF1 gradientIpar("gradientIpar",
 361                          ROOT::TF1Helper::TGradientParFunction(iPar,funcCloneMM),0,0,0);
 362         TF1 gradientJpar("gradientJpar",
 363                          ROOT::TF1Helper::TGradientParFunction(jPar,funcCloneMM),0,0,0);
 364         emm += mmFitResult->CovMatrix(iPar,jPar)*
 365           gradientIpar.Eval(m)*gradientJpar.Eval(m)*
 366           (iPar==jPar ? 1.0 : 2.0);
 367       }
 368 */
 369     }
 370
 371     Double_t signal = pm-2.0*TMath::Sqrt(pp*mm);
 372     Double_t background = 2.0*TMath::Sqrt(pp*mm);
 373     // error propagation on the signal calculation above
 374     Double_t esignal = TMath::Sqrt(epm*epm+(mm/pp)*epp+(pp/mm)*emm);
 375     Double_t ebackground = TMath::Sqrt((mm/pp)*epp+(pp/mm)*emm);
 376     fHistSignal->SetBinContent(iBin, signal);
 377     fHistSignal->SetBinError(iBin, esignal);
 378     fHistBackground->SetBinContent(iBin, background);
 379     fHistBackground->SetBinError(iBin, ebackground);
 380   }
 381
 382   // signal
 383   fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
 384                                              fHistSignal->FindBin(fIntMax), fErrors(0));
 385   // background
 386   fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
 387                                                  fHistBackground->FindBin(fIntMax),
 388                                                  fErrors(1));
 389   // S/B and significance
 390   SetSignificanceAndSOB();
 391   fValues(4) = fFuncSigBack->GetParameter(fParMass);
 392   fErrors(4) = fFuncSigBack->GetParError(fParMass);
 393   fValues(5) = fFuncSigBack->GetParameter(fParMassWidth);
 394   fErrors(5) = fFuncSigBack->GetParError(fParMassWidth);
 395
 396   fProcessed = kTRUE;
 397 }
 398
 399 //______________________________________________
 400 void AliDielectronSignalFunc::ProcessEM(TObjArray * const arrhist) {
 401   //
 402   // Substract background with the event mixing technique
 403   //
 404   arrhist->GetEntries();   // just to avoid the unused parameter warning
 405   AliError("Event mixing for background substraction method not implemented!");
 406 }
 407
 408 //______________________________________________
 409 void AliDielectronSignalFunc::SetFunctions(TF1 * const combined, TF1 * const sig, TF1 * const back,
 410                                            Int_t parM, Int_t parMres)
 411 {
 412   //
 413   // Set the signal, background functions and combined fit function
 414   // Note: The process method assumes that the first n parameters in the
 415   //       combined fit function correspond to the n parameters of the signal function
 416   //       and the n+1 to n+m parameters to the m parameters of the background function!!!
 417
 418   if (!sig||!back||!combined) {
 419     AliError("Both, signal and background function need to be set!");
 420     return;
 421   }
 422   fFuncSignal=sig;
 423   fFuncBackground=back;
 424   fFuncSigBack=combined;
 425   fParMass=parM;
 426   fParMassWidth=parMres;
 427
 428 }
 429
 430 //______________________________________________
 431 void AliDielectronSignalFunc::SetDefaults(Int_t type)
 432 {
 433   //
 434   // Setup some default functions:
 435   // type = 0: gaus signal + linear background in 2.5 - 4 GeV inv. mass
 436   // type = 1: gaus signal + exponential background in 2.5 - 4 GeV inv. mass
 437   // type = 2: half gaussian, half exponential signal function
 438   // type = 3: Crystal-Ball function
 439   // type = 4: Crystal-Ball signal + exponential background
 440   //
 441
 442   if (type==0){
 443     fFuncSignal=new TF1("DieleSignal","gaus",2.5,4);
 444     fFuncBackground=new TF1("DieleBackground","pol1",2.5,4);
 445     fFuncSigBack=new TF1("DieleCombined","gaus+pol1(3)",2.5,4);
 446
 447     fFuncSigBack->SetParameters(1,3.1,.05,2.5,1);
 448     fFuncSigBack->SetParLimits(0,0,10000000);
 449     fFuncSigBack->SetParLimits(1,3.05,3.15);
 450     fFuncSigBack->SetParLimits(2,.02,.1);
 451   }
 452   else if (type==1){
 453     fFuncSignal=new TF1("DieleSignal","gaus",2.5,4);
 454     fFuncBackground=new TF1("DieleBackground","[0]*exp(-(x-[1])/[2])",2.5,4);
 455     fFuncSigBack=new TF1("DieleCombined","gaus+[3]*exp(-(x-[4])/[5])",2.5,4);
 456
 457     fFuncSigBack->SetParameters(1,3.1,.05,1,2.5,1);
 458     fFuncSigBack->SetParLimits(0,0,10000000);
 459     fFuncSigBack->SetParLimits(1,3.05,3.15);
 460     fFuncSigBack->SetParLimits(2,.02,.1);
 461   }
 462   else if (type==2){
 463     // half gaussian, half exponential signal function
 464     // exponential background
 465     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);
 466     fFuncBackground = new TF1("DieleBackground","[0]*exp(-(x-[1])/[2])+[3]",2.5,4);
 467     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);
 468     fFuncSigBack->SetParameters(1.,3.1,.05,.1,1,2.5,1,0);
 469
 470     fFuncSigBack->SetParLimits(0,0,10000000);
 471     fFuncSigBack->SetParLimits(1,3.05,3.15);
 472     fFuncSigBack->SetParLimits(2,.02,.1);
 473     fFuncSigBack->FixParameter(6,2.5);
 474     fFuncSigBack->FixParameter(7,0);
 475   }
 476 }
 477
 478
 479 //______________________________________________
 480 void AliDielectronSignalFunc::Draw(const Option_t* option)
 481 {
 482   //
 483   // Draw the fitted function
 484   //
 485
 486   TString drawOpt(option);
 487   drawOpt.ToLower();
 488
 489   Bool_t optStat=drawOpt.Contains("stat");
 490
 491   fFuncSigBack->SetNpx(200);
 492   fFuncSigBack->SetRange(fIntMin,fIntMax);
 493   fFuncBackground->SetNpx(200);
 494   fFuncBackground->SetRange(fIntMin,fIntMax);
 495
 496   TGraph *grSig=new TGraph(fFuncSigBack);
 497   grSig->SetFillColor(kGreen);
 498   grSig->SetFillStyle(3001);
 499
 500   TGraph *grBack=new TGraph(fFuncBackground);
 501   grBack->SetFillColor(kRed);
 502   grBack->SetFillStyle(3001);
 503
 504   grSig->SetPoint(0,grBack->GetX()[0],grBack->GetY()[0]);
 505   grSig->SetPoint(grSig->GetN()-1,grBack->GetX()[grBack->GetN()-1],grBack->GetY()[grBack->GetN()-1]);
 506
 507   grBack->SetPoint(0,grBack->GetX()[0],0.);
 508   grBack->SetPoint(grBack->GetN()-1,grBack->GetX()[grBack->GetN()-1],0.);
 509
 510   fFuncSigBack->SetRange(fFitMin,fFitMax);
 511   fFuncBackground->SetRange(fFitMin,fFitMax);
 512
 513   if (!drawOpt.Contains("same")){
 514     if (fHistDataPM){
 515       fHistDataPM->Draw();
 516       grSig->Draw("f");
 517     } else {
 518       grSig->Draw("af");
 519     }
 520   } else {
 521     grSig->Draw("f");
 522   }
 523   if(fMethod==kFitted || fMethod==kFittedMC) grBack->Draw("f");
 524   fFuncSigBack->Draw("same");
 525   fFuncSigBack->SetLineWidth(2);
 526   if(fMethod==kLikeSign) {
 527     fHistDataPP->SetLineWidth(2);
 528     fHistDataPP->SetLineColor(6);
 529     fHistDataPP->Draw("same");
 530     fHistDataMM->SetLineWidth(2);
 531     fHistDataMM->SetLineColor(8);
 532     fHistDataMM->Draw("same");
 533   }
 534
 535   if(fMethod==kFitted || fMethod==kFittedMC)
 536     fFuncBackground->Draw("same");
 537
 538   if (optStat) DrawStats();
 539
 540 }
 541
 542
 543 //______________________________________________________________________________
 544 void AliDielectronSignalFunc::CombineFunc(TF1 * const peak, TF1 * const bgnd) {
 545   //
 546   // combine the bgnd and the peak function
 547   //
 548
 549   if (!peak||!bgnd) {
 550     AliError("Both, signal and background function need to be set!");
 551     return;
 552   }
 553   fFuncSignal=peak;
 554   fFuncBackground=bgnd;
 555
 556   fNparPeak     = fFuncSignal->GetNpar();
 557   fNparBgnd     = fFuncBackground->GetNpar();
 558
 559   fFuncSigBack = new TF1("BgndPeak",AliDielectronSignalFunc::PeakBgndFun, 0.,5.,fNparPeak+fNparBgnd);
 560   return;
 561 }
 562
 563 //______________________________________________________________________________
 564 Double_t AliDielectronSignalFunc::PeakBgndFun(const Double_t *x, const Double_t *par) {
 565   //
 566   // merge peak and bgnd functions
 567   //
 568   return (fFuncBackground->EvalPar(x,par) + fFuncSignal->EvalPar(x,par));
 569 }
 570