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
  22   Class used for extracting the signal from an invariant mass spectrum.
  23   It implements the AliDielectronSignalBase and -Ext classes and it uses user provided
  24   functions to fit the spectrum with a combined signa+background fit.
  25   Used invariant mass spectra are provided via an array of histograms. There are serveral method
  26   to estimate the background and to extract the raw yield from the background subtracted spectra.
  27
  28   Example usage:
  29
  30   AliDielectronSignalFunc *sig = new AliDielectronSignalFunc();
  31
  32
  33   1) invariant mass input spectra
  34
  35   1.1) Assuming a AliDielectronCF container as data format (check class for more details)
  36   AliDielectronCFdraw *cf = new AliDielectronCFdraw("path/to/the/output/file.root");
  37   TObjArray *arrHists = cf->CollectMinvProj(cf->FindStep("Config"));
  38
  39   1.2) Assuming a AliDielectronHF grid as data format (check class for more details)
  40   AliDielectronHFhelper *hf = new AliDielectronHFhelper("path/to/the/output/file.root", "ConfigName");
  41   TObjArray *arrHists = hf->CollectHistos(AliDielectronVarManager::kM);
  42
  43   1.3) Assuming a single histograms
  44   TObjArray *histoArray = new TObjArray();
  45   arrHists->Add(signalPP);            // add the spectrum histograms to the array
  46   arrHists->Add(signalPM);            // the order is important !!!
  47   arrHists->Add(signalMM);
  48
  49
  50   2) background estimation
  51
  52   2.1) set the method for the background estimation (methods can be found in AliDielectronSignalBase)
  53   sig->SetMethod(AliDielectronSignalBase::kFitted);
  54   2.2) rebin the spectras if needed
  55   //  sig->SetRebin(2);
  56   2.3) add any background function you like
  57   TF1 *fB = new TF1("fitBgrd","pol3",minFit,maxFit);
  58
  59
  60   3) configure the signal extraction
  61
  62   3.1) chose one of the signal functions (MCshape, CrystalBall, Gauss)
  63   TF1 *fS = new TF1("fitSign",AliDielectronSignalFunc::PeakFunCB,minFit,maxFit,5); // has 5 parameters
  64   //  TF1 *fS = new TF1("fitSign",AliDielectronSignalFunc::PeakFunGaus,minFit,maxFit,3); // has 3 parameters
  65   //  sig->SetMCSignalShape(hMCsign);
  66   //  TF1 *fS = new TF1("fitSign",AliDielectronSignalFunc::PeakFunMC,minFit,maxFit,1); // requires a MC shape
  67   3.2) set the method for the signal extraction (methods can be found in AliDielectronSignalBase)
  68   depending on the method serveral inputs are needed (e.g. MC shape, PDG code of the particle of interest)
  69   //  sig->SetParticleOfInterest(443); //default is jpsi
  70   //  sig->SetMCSignalShape(signalMC);
  71   //  sig->SetIntegralRange(minInt, maxInt);
  72   sig->SetExtractionMethod(AliDielectronSignal::BinCounting); // this is the default
  73
  74
  75   4) combined fit of bgrd+signal
  76
  77   4.1) combine the two functions
  78   sig->CombineFunc(fS,fB);
  79   4.2) apply fitting ranges and the fit options
  80   sig->SetFitRange(minFit, maxFit);
  81   sig->SetFitOption("NR");
  82
  83
  84   5) start the processing
  85
  86   sig->Process(arrHists);
  87   sig->Print(""); // print values and errors extracted
  88
  89
  90   6) access the spectra and values created
  91
  92   6.1) standard spectra as provided filled in AliDielectronSignalExt
  93   TH1F *hsign = (TH1F*) sig->GetUnlikeSignHistogram();  // same as the input (rebinned)
  94   TH1F *hbgrd = (TH1F*) sig->GetBackgroundHistogram();  // filled histogram with fitBgrd
  95   TH1F *hextr = (TH1F*) sig->GetSignalHistogram();      // after backgound extraction (rebinned)
  96   TObject *oPeak = (TObject*) sig->GetPeakShape();      // can be a TF1 or TH1 depending on the method
  97   6.2) flow spectra
  98   TF1 *fFitSign  = sig->GetCombinedFunction();                // combined fit function
  99   TF1 *fFitExtr  = sig->GetSignalFunction();                  // signal function
 100   TF1 *fFitBgrd  = sig->GetBackgroundFunction();              // background function
 101   6.3) access the extracted values and errors
 102   sig->GetValues();     or GetErrors();                 // yield extraction
 103
 104 */
 105 //                                                                       //
 106 ///////////////////////////////////////////////////////////////////////////
 107
 108 #include <TF1.h>
 109 #include <TH1.h>
 110 #include <TGraph.h>
 111 #include <TMath.h>
 112 #include <TString.h>
 113 #include <TPaveText.h>
 114 #include <TList.h>
 115 #include <TFitResult.h>
 116 //#include <../hist/hist/src/TF1Helper.h>
 117
 118 #include <AliLog.h>
 119
 120 #include "AliDielectronSignalFunc.h"
 121
 122 ClassImp(AliDielectronSignalFunc)
 123
 124 //TH1F* AliDielectronSignalFunc::fgHistSimPM=0x0;
 125 TF1*  AliDielectronSignalFunc::fFuncSignal=0x0;
 126 TF1*  AliDielectronSignalFunc::fFuncBackground=0x0;
 127 Int_t AliDielectronSignalFunc::fNparPeak=0;
 128 Int_t AliDielectronSignalFunc::fNparBgnd=0;
 129
 130
 131 AliDielectronSignalFunc::AliDielectronSignalFunc() :
 132 AliDielectronSignalExt(),
 133 fFuncSigBack(0x0),
 134 fParMass(1),
 135 fParMassWidth(2),
 136 fFitOpt("SMNQE"),
 137 fUseIntegral(kFALSE),
 138 fDof(0),
 139 fChi2Dof(0.0)
 140 {
 141   //
 142   // Default Constructor
 143   //
 144
 145 }
 146
 147 //______________________________________________
 148 AliDielectronSignalFunc::AliDielectronSignalFunc(const char* name, const char* title) :
 149 AliDielectronSignalExt(name, title),
 150 fFuncSigBack(0x0),
 151 fParMass(1),
 152 fParMassWidth(2),
 153 fFitOpt("SMNQE"),
 154 fUseIntegral(kFALSE),
 155 fDof(0),
 156 fChi2Dof(0.0)
 157 {
 158   //
 159   // Named Constructor
 160   //
 161 }
 162
 163 //______________________________________________
 164 AliDielectronSignalFunc::~AliDielectronSignalFunc()
 165 {
 166   //
 167   // Default Destructor
 168   //
 169   if(fFuncSigBack) delete fFuncSigBack;
 170 }
 171
 172
 173 //______________________________________________
 174 void AliDielectronSignalFunc::Process(TObjArray * const arrhist)
 175 {
 176   //
 177   // Fit the invariant mass histograms and retrieve the signal and background
 178   //
 179   switch(fMethod) {
 180   case kFitted :
 181     ProcessFit(arrhist);
 182     break;
 183
 184   case kLikeSignFit :
 185     ProcessFitLS(arrhist);
 186     break;
 187
 188   case kEventMixingFit :
 189     ProcessFitEM(arrhist);
 190     break;
 191
 192   case kLikeSign :
 193   case kLikeSignArithm :
 194   case kLikeSignRcorr:
 195   case kLikeSignArithmRcorr:
 196     ProcessLS(arrhist);    // process like-sign subtraction method
 197     break;
 198
 199   case kEventMixing :
 200     ProcessEM(arrhist);    // process event mixing method
 201     break;
 202
 203   case kRotation:
 204     ProcessRotation(arrhist);
 205     break;
 206
 207   default :
 208     AliError("Background substraction method not known!");
 209   }
 210 }
 211 //______________________________________________________________________________
 212 Double_t AliDielectronSignalFunc::PeakFunMC(const Double_t *x, const Double_t *par) {
 213   // Fit MC signal shape
 214   // parameters
 215   // [0]:   scale for simpeak
 216
 217   Double_t xx  = x[0];
 218
 219   TH1F *hPeak = fgHistSimPM;
 220   if (!hPeak) {
 221     printf("E-AliDielectronSignalFunc::PeakFun: No histogram for peak fit defined!\n");
 222     return 0.0;
 223   }
 224
 225   Int_t idx = hPeak->FindBin(xx);
 226   if ((idx <= 1) ||
 227       (idx >= hPeak->GetNbinsX())) {
 228     return 0.0;
 229   }
 230
 231   return (par[0] * hPeak->GetBinContent(idx));
 232
 233 }
 234
 235 //______________________________________________________________________________
 236 Double_t AliDielectronSignalFunc::PeakFunCB(const Double_t *x, const Double_t *par) {
 237   // Crystal Ball fit function
 238
 239   Double_t alpha = par[0];
 240   Double_t     n = par[1];
 241   Double_t meanx = par[2];
 242   Double_t sigma = par[3];
 243   Double_t    nn = par[4];
 244
 245   Double_t a = TMath::Power((n/TMath::Abs(alpha)), n) * TMath::Exp(-.5*alpha*alpha);
 246   Double_t b = n/TMath::Abs(alpha) - TMath::Abs(alpha);
 247
 248   Double_t arg = (x[0] - meanx)/sigma;
 249   Double_t fitval = 0;
 250
 251   if (arg > -1.*alpha) {
 252     fitval = nn * TMath::Exp(-.5*arg*arg);
 253   } else {
 254     fitval = nn * a * TMath::Power((b-arg), (-1*n));
 255   }
 256
 257   return fitval;
 258 }
 259
 260 //______________________________________________________________________________
 261 Double_t AliDielectronSignalFunc::PeakFunGaus(const Double_t *x, const Double_t *par) {
 262   // Gaussian fit function
 263   //printf("fNparBgrd %d \n",fNparBgnd);
 264   Double_t     n = par[0];
 265   Double_t  mean = par[1];
 266   Double_t sigma = par[2];
 267   Double_t    xx = x[0];
 268
 269   return ( n*TMath::Exp(-0.5*TMath::Power((xx-mean)/sigma,2)) );
 270 }
 271
 272 //______________________________________________
 273 void AliDielectronSignalFunc::ProcessFit(TObjArray * const arrhist) {
 274   //
 275   // Fit the +- invariant mass distribution only
 276   // Here we assume that the combined fit function is a sum of the signal and background functions
 277   //    and that the signal function is always the first term of this sum
 278   //
 279
 280   fHistDataPM = (TH1F*)(arrhist->At(1))->Clone("histPM");  // +-    SE
 281   fHistDataPM->Sumw2();
 282   if(fRebin>1)
 283     fHistDataPM->Rebin(fRebin);
 284
 285   fHistSignal = new TH1F("HistSignal", "Fit substracted signal",
 286                          fHistDataPM->GetXaxis()->GetNbins(),
 287                          fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
 288   fHistBackground = new TH1F("HistBackground", "Fit contribution",
 289                              fHistDataPM->GetXaxis()->GetNbins(),
 290                              fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
 291
 292   // the starting parameters of the fit function and their limits can be tuned
 293   // by the user in its macro
 294   fHistDataPM->Fit(fFuncSigBack, fFitOpt.Data(), "", fFitMin, fFitMax);
 295   TFitResultPtr pmFitPtr = fHistDataPM->Fit(fFuncSigBack, fFitOpt.Data(), "", fFitMin, fFitMax);
 296   //TFitResult *pmFitResult = pmFitPtr.Get(); // used only with TF1Helper
 297   //fFuncBackground->SetParameters(fFuncSigBack->GetParameters());
 298   fFuncSignal->SetParameters(fFuncSigBack->GetParameters());
 299   fFuncBackground->SetParameters(fFuncSigBack->GetParameters()+fFuncSignal->GetNpar());
 300
 301   // fill the background spectrum
 302   fHistBackground->Eval(fFuncBackground);
 303   // set the error for the background histogram
 304   fHistBackground->Fit(fFuncBackground,"0qR","",fFitMin,fFitMax);
 305   Double_t inte  = fFuncBackground->IntegralError(fIntMin, fIntMax)/fHistDataPM->GetBinWidth(1);
 306   Double_t binte = inte / TMath::Sqrt((fHistDataPM->FindBin(fIntMax)-fHistDataPM->FindBin(fIntMin))+1);
 307   for(Int_t iBin=fHistDataPM->FindBin(fIntMin); iBin<=fHistDataPM->FindBin(fIntMax); iBin++) {
 308     fHistBackground->SetBinError(iBin, binte);
 309   }
 310
 311   for(Int_t iBin=1; iBin<=fHistDataPM->GetXaxis()->GetNbins(); iBin++) {
 312     //    Double_t m = fHistDataPM->GetBinCenter(iBin);
 313     Double_t pm = fHistDataPM->GetBinContent(iBin);
 314     Double_t epm = fHistDataPM->GetBinError(iBin);
 315     Double_t bknd = fHistBackground->GetBinContent(iBin);
 316     Double_t ebknd = fHistBackground->GetBinError(iBin);
 317     for(Int_t iPar=fFuncSignal->GetNpar(); iPar<fFuncSigBack->GetNpar(); iPar++) {
 318       /* TF1Helper problem on alien compilation
 319          for(Int_t jPar=iPar; jPar<fFuncSigBack->GetNpar(); jPar++) {
 320          TF1 gradientIpar("gradientIpar",
 321          ROOT::TF1Helper::TGradientParFunction(iPar-fFuncSignal->GetNpar(),fFuncBackground),0,0,0);
 322          TF1 gradientJpar("gradientJpar",
 323          ROOT::TF1Helper::TGradientParFunction(jPar-fFuncSignal->GetNpar(),fFuncBackground),0,0,0);
 324          ebknd += pmFitResult->CovMatrix(iPar,jPar)*
 325          gradientIpar.Eval(m)*gradientJpar.Eval(m)*
 326          (iPar==jPar ? 1.0 : 2.0);
 327          }
 328       */
 329     }
 330     Double_t signal = pm-bknd;
 331     Double_t error = TMath::Sqrt(epm*epm+ebknd);
 332     fHistSignal->SetBinContent(iBin, signal);
 333     fHistSignal->SetBinError(iBin, error);
 334   }
 335
 336   if(fUseIntegral) {
 337     // signal
 338     fValues(0) = fFuncSignal->Integral(fIntMin, fIntMax)/fHistDataPM->GetBinWidth(1);
 339     fErrors(0) = 0;
 340     for(Int_t iPar=0; iPar<fFuncSignal->GetNpar(); iPar++) {
 341 /* TF1Helper problem on alien compilation
 342       for(Int_t jPar=iPar; jPar<fFuncSignal->GetNpar(); jPar++) {
 343         TF1 gradientIpar("gradientIpar",
 344                          ROOT::TF1Helper::TGradientParFunction(iPar,fFuncSignal),0,0,0);
 345         TF1 gradientJpar("gradientJpar",
 346                          ROOT::TF1Helper::TGradientParFunction(jPar,fFuncSignal),0,0,0);
 347         fErrors(0) += pmFitResult->CovMatrix(iPar,jPar)*
 348           gradientIpar.Integral(fIntMin,fIntMax)*gradientJpar.Integral(fIntMin,fIntMax)*
 349           (iPar==jPar ? 1.0 : 2.0);
 350       }
 351 */
 352     }
 353     // background
 354     fValues(1) = fFuncBackground->Integral(fIntMin, fIntMax)/fHistDataPM->GetBinWidth(1);
 355     fErrors(1) = fFuncBackground->IntegralError(fIntMin, fIntMax)/fHistDataPM->GetBinWidth(1);
 356     for(Int_t iPar=fFuncSignal->GetNpar(); iPar<fFuncSigBack->GetNpar(); iPar++) {
 357 /* TF1Helper problem on alien compilation
 358       for(Int_t jPar=iPar; jPar<fFuncSigBack->GetNpar(); jPar++) {
 359         TF1 gradientIpar("gradientIpar",
 360                          ROOT::TF1Helper::TGradientParFunction(iPar-fFuncSignal->GetNpar(),fFuncBackground),0,0,0);
 361         TF1 gradientJpar("gradientJpar",
 362                          ROOT::TF1Helper::TGradientParFunction(jPar-fFuncSignal->GetNpar(),fFuncBackground),0,0,0);
 363         fErrors(1) += pmFitResult->CovMatrix(iPar,jPar)*
 364           gradientIpar.Integral(fIntMin, fIntMax)*gradientJpar.Integral(fIntMin, fIntMax)*
 365           (iPar==jPar ? 1.0 : 2.0);
 366       }
 367 */
 368     }
 369   }
 370   else {
 371     // signal
 372     fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
 373                                                fHistSignal->FindBin(fIntMax), fErrors(0));
 374     // background
 375     fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
 376                                                    fHistBackground->FindBin(fIntMax),
 377                                                    fErrors(1));
 378   }
 379   // S/B and significance
 380   SetSignificanceAndSOB();
 381   fValues(4) = fFuncSigBack->GetParameter(fParMass);
 382   fErrors(4) = fFuncSigBack->GetParError(fParMass);
 383   fValues(5) = fFuncSigBack->GetParameter(fParMassWidth);
 384   fErrors(5) = fFuncSigBack->GetParError(fParMassWidth);
 385
 386   fProcessed = kTRUE;
 387
 388   fHistBackground->GetListOfFunctions()->Add(fFuncBackground);
 389
 390 }
 391
 392 //______________________________________________
 393 void AliDielectronSignalFunc::ProcessFitLS(TObjArray * const arrhist) {
 394   //
 395   // Substract background using the like-sign spectrum
 396   //
 397   fHistDataPP = (TH1F*)(arrhist->At(0))->Clone("histPP");  // ++    SE
 398   fHistDataPM = (TH1F*)(arrhist->At(1))->Clone("histPM");  // +-    SE
 399   fHistDataMM = (TH1F*)(arrhist->At(2))->Clone("histMM");  // --    SE
 400   if (fRebin>1) {
 401     fHistDataPP->Rebin(fRebin);
 402     fHistDataPM->Rebin(fRebin);
 403     fHistDataMM->Rebin(fRebin);
 404   }
 405   fHistDataPP->Sumw2();
 406   fHistDataPM->Sumw2();
 407   fHistDataMM->Sumw2();
 408
 409   fHistSignal = new TH1F("HistSignal", "Like-Sign substracted signal",
 410                          fHistDataPM->GetXaxis()->GetNbins(),
 411                          fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
 412   fHistBackground = new TH1F("HistBackground", "Like-sign contribution",
 413                              fHistDataPM->GetXaxis()->GetNbins(),
 414                              fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
 415
 416   // fit the +- mass distribution
 417   fHistDataPM->Fit(fFuncSigBack, fFitOpt.Data(), "", fFitMin, fFitMax);
 418   fHistDataPM->Fit(fFuncSigBack, fFitOpt.Data(), "", fFitMin, fFitMax);
 419   // declare the variables where the like-sign fit results will be stored
 420 //   TFitResult *ppFitResult = 0x0;
 421 //   TFitResult *mmFitResult = 0x0;
 422   // fit the like sign background
 423   TF1 *funcClonePP = (TF1*)fFuncBackground->Clone("funcClonePP");
 424   TF1 *funcCloneMM = (TF1*)fFuncBackground->Clone("funcCloneMM");
 425   fHistDataPP->Fit(funcClonePP, fFitOpt.Data(), "", fFitMin, fFitMax);
 426   fHistDataPP->Fit(funcClonePP, fFitOpt.Data(), "", fFitMin, fFitMax);
 427 //   TFitResultPtr ppFitPtr = fHistDataPP->Fit(funcClonePP, fFitOpt.Data(), "", fFitMin, fFitMax);
 428 //   ppFitResult = ppFitPtr.Get();
 429   fHistDataMM->Fit(funcCloneMM, fFitOpt.Data(), "", fFitMin, fFitMax);
 430   fHistDataMM->Fit(funcCloneMM, fFitOpt.Data(), "", fFitMin, fFitMax);
 431 //   TFitResultPtr mmFitPtr = fHistDataMM->Fit(funcCloneMM, fFitOpt.Data(), "", fFitMin, fFitMax);
 432 //   mmFitResult = mmFitPtr.Get();
 433
 434   for(Int_t iBin=1; iBin<=fHistDataPM->GetXaxis()->GetNbins(); iBin++) {
 435     Double_t m = fHistDataPM->GetBinCenter(iBin);
 436     Double_t pm = fHistDataPM->GetBinContent(iBin);
 437     Double_t pp = funcClonePP->Eval(m);
 438     Double_t mm = funcCloneMM->Eval(m);
 439     Double_t epm = fHistDataPM->GetBinError(iBin);
 440     Double_t epp = 0;
 441     for(Int_t iPar=0; iPar<funcClonePP->GetNpar(); iPar++) {
 442 /* TF1Helper problem on alien compilation
 443       for(Int_t jPar=iPar; jPar<funcClonePP->GetNpar(); jPar++) {
 444         TF1 gradientIpar("gradientIpar",
 445                          ROOT::TF1Helper::TGradientParFunction(iPar,funcClonePP),0,0,0);
 446         TF1 gradientJpar("gradientJpar",
 447                          ROOT::TF1Helper::TGradientParFunction(jPar,funcClonePP),0,0,0);
 448         epp += ppFitResult->CovMatrix(iPar,jPar)*
 449           gradientIpar.Eval(m)*gradientJpar.Eval(m)*
 450           (iPar==jPar ? 1.0 : 2.0);
 451       }
 452 */
 453     }
 454     Double_t emm = 0;
 455     for(Int_t iPar=0; iPar<funcCloneMM->GetNpar(); iPar++) {
 456 /* TF1Helper problem on alien compilation
 457       for(Int_t jPar=iPar; jPar<funcCloneMM->GetNpar(); jPar++) {
 458         TF1 gradientIpar("gradientIpar",
 459                          ROOT::TF1Helper::TGradientParFunction(iPar,funcCloneMM),0,0,0);
 460         TF1 gradientJpar("gradientJpar",
 461                          ROOT::TF1Helper::TGradientParFunction(jPar,funcCloneMM),0,0,0);
 462         emm += mmFitResult->CovMatrix(iPar,jPar)*
 463           gradientIpar.Eval(m)*gradientJpar.Eval(m)*
 464           (iPar==jPar ? 1.0 : 2.0);
 465       }
 466 */
 467     }
 468
 469     Double_t signal = pm-2.0*TMath::Sqrt(pp*mm);
 470     Double_t background = 2.0*TMath::Sqrt(pp*mm);
 471     // error propagation on the signal calculation above
 472     Double_t esignal = TMath::Sqrt(epm*epm+(mm/pp)*epp+(pp/mm)*emm);
 473     Double_t ebackground = TMath::Sqrt((mm/pp)*epp+(pp/mm)*emm);
 474     fHistSignal->SetBinContent(iBin, signal);
 475     fHistSignal->SetBinError(iBin, esignal);
 476     fHistBackground->SetBinContent(iBin, background);
 477     fHistBackground->SetBinError(iBin, ebackground);
 478   }
 479
 480   // signal
 481   fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
 482                                              fHistSignal->FindBin(fIntMax), fErrors(0));
 483   // background
 484   fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
 485                                                  fHistBackground->FindBin(fIntMax),
 486                                                  fErrors(1));
 487   // S/B and significance
 488   SetSignificanceAndSOB();
 489   fValues(4) = fFuncSigBack->GetParameter(fParMass);
 490   fErrors(4) = fFuncSigBack->GetParError(fParMass);
 491   fValues(5) = fFuncSigBack->GetParameter(fParMassWidth);
 492   fErrors(5) = fFuncSigBack->GetParError(fParMassWidth);
 493
 494   fProcessed = kTRUE;
 495 }
 496
 497 //______________________________________________
 498 void AliDielectronSignalFunc::ProcessFitEM(TObjArray * const arrhist) {
 499   //
 500   // Substract background with the event mixing technique
 501   //
 502   arrhist->GetEntries();   // just to avoid the unused parameter warning
 503   AliError("Event mixing for background substraction method not implemented!");
 504 }
 505
 506 //______________________________________________
 507 void AliDielectronSignalFunc::SetFunctions(TF1 * const combined, TF1 * const sig, TF1 * const back,
 508                                            Int_t parM, Int_t parMres)
 509 {
 510   //
 511   // Set the signal, background functions and combined fit function
 512   // Note: The process method assumes that the first n parameters in the
 513   //       combined fit function correspond to the n parameters of the signal function
 514   //       and the n+1 to n+m parameters to the m parameters of the background function!!!
 515
 516   if (!sig||!back||!combined) {
 517     AliError("Both, signal and background function need to be set!");
 518     return;
 519   }
 520   fFuncSignal=sig;
 521   fFuncBackground=back;
 522   fFuncSigBack=combined;
 523   fParMass=parM;
 524   fParMassWidth=parMres;
 525
 526 }
 527
 528 //______________________________________________
 529 void AliDielectronSignalFunc::SetDefaults(Int_t type)
 530 {
 531   //
 532   // Setup some default functions:
 533   // type = 0: gaus signal + linear background in 2.5 - 4 GeV inv. mass
 534   // type = 1: gaus signal + exponential background in 2.5 - 4 GeV inv. mass
 535   // type = 2: half gaussian, half exponential signal function
 536   // type = 3: Crystal-Ball function
 537   // type = 4: Crystal-Ball signal + exponential background
 538   //
 539
 540   if (type==0){
 541     fFuncSignal=new TF1("DieleSignal","gaus",2.5,4);
 542     fFuncBackground=new TF1("DieleBackground","pol1",2.5,4);
 543     fFuncSigBack=new TF1("DieleCombined","gaus+pol1(3)",2.5,4);
 544
 545     fFuncSigBack->SetParameters(1,3.1,.05,2.5,1);
 546     fFuncSigBack->SetParLimits(0,0,10000000);
 547     fFuncSigBack->SetParLimits(1,3.05,3.15);
 548     fFuncSigBack->SetParLimits(2,.02,.1);
 549   }
 550   else if (type==1){
 551     fFuncSignal=new TF1("DieleSignal","gaus",2.5,4);
 552     fFuncBackground=new TF1("DieleBackground","[0]*exp(-(x-[1])/[2])",2.5,4);
 553     fFuncSigBack=new TF1("DieleCombined","gaus+[3]*exp(-(x-[4])/[5])",2.5,4);
 554
 555     fFuncSigBack->SetParameters(1,3.1,.05,1,2.5,1);
 556     fFuncSigBack->SetParLimits(0,0,10000000);
 557     fFuncSigBack->SetParLimits(1,3.05,3.15);
 558     fFuncSigBack->SetParLimits(2,.02,.1);
 559   }
 560   else if (type==2){
 561     // half gaussian, half exponential signal function
 562     // exponential background
 563     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);
 564     fFuncBackground = new TF1("DieleBackground","[0]*exp(-(x-[1])/[2])+[3]",2.5,4);
 565     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);
 566     fFuncSigBack->SetParameters(1.,3.1,.05,.1,1,2.5,1,0);
 567
 568     fFuncSigBack->SetParLimits(0,0,10000000);
 569     fFuncSigBack->SetParLimits(1,3.05,3.15);
 570     fFuncSigBack->SetParLimits(2,.02,.1);
 571     fFuncSigBack->FixParameter(6,2.5);
 572     fFuncSigBack->FixParameter(7,0);
 573   }
 574 }
 575
 576
 577 //______________________________________________
 578 void AliDielectronSignalFunc::Draw(const Option_t* option)
 579 {
 580   //
 581   // Draw the fitted function
 582   //
 583
 584   TString drawOpt(option);
 585   drawOpt.ToLower();
 586
 587   Bool_t optStat=drawOpt.Contains("stat");
 588
 589   fFuncSigBack->SetNpx(200);
 590   fFuncSigBack->SetRange(fIntMin,fIntMax);
 591   fFuncBackground->SetNpx(200);
 592   fFuncBackground->SetRange(fIntMin,fIntMax);
 593
 594   TGraph *grSig=new TGraph(fFuncSigBack);
 595   grSig->SetFillColor(kGreen);
 596   grSig->SetFillStyle(3001);
 597
 598   TGraph *grBack=new TGraph(fFuncBackground);
 599   grBack->SetFillColor(kRed);
 600   grBack->SetFillStyle(3001);
 601
 602   grSig->SetPoint(0,grBack->GetX()[0],grBack->GetY()[0]);
 603   grSig->SetPoint(grSig->GetN()-1,grBack->GetX()[grBack->GetN()-1],grBack->GetY()[grBack->GetN()-1]);
 604
 605   grBack->SetPoint(0,grBack->GetX()[0],0.);
 606   grBack->SetPoint(grBack->GetN()-1,grBack->GetX()[grBack->GetN()-1],0.);
 607
 608   fFuncSigBack->SetRange(fFitMin,fFitMax);
 609   fFuncBackground->SetRange(fFitMin,fFitMax);
 610
 611   if (!drawOpt.Contains("same")){
 612     if (fHistDataPM){
 613       fHistDataPM->Draw();
 614       grSig->Draw("f");
 615     } else {
 616       grSig->Draw("af");
 617     }
 618   } else {
 619     grSig->Draw("f");
 620   }
 621   if(fMethod==kFitted || fMethod==kFittedMC) grBack->Draw("f");
 622   fFuncSigBack->Draw("same");
 623   fFuncSigBack->SetLineWidth(2);
 624   if(fMethod==kLikeSign) {
 625     fHistDataPP->SetLineWidth(2);
 626     fHistDataPP->SetLineColor(6);
 627     fHistDataPP->Draw("same");
 628     fHistDataMM->SetLineWidth(2);
 629     fHistDataMM->SetLineColor(8);
 630     fHistDataMM->Draw("same");
 631   }
 632
 633   if(fMethod==kFitted || fMethod==kFittedMC)
 634     fFuncBackground->Draw("same");
 635
 636   if (optStat) DrawStats();
 637
 638 }
 639
 640
 641 //______________________________________________________________________________
 642 void AliDielectronSignalFunc::CombineFunc(TF1 * const peak, TF1 * const bgnd) {
 643   //
 644   // combine the bgnd and the peak function
 645   //
 646
 647   if (!peak||!bgnd) {
 648     AliError("Both, signal and background function need to be set!");
 649     return;
 650   }
 651   fFuncSignal=peak;
 652   fFuncBackground=bgnd;
 653
 654   fNparPeak     = fFuncSignal->GetNpar();
 655   fNparBgnd     = fFuncBackground->GetNpar();
 656
 657   fFuncSigBack = new TF1("BgndPeak",AliDielectronSignalFunc::PeakBgndFun, fFitMin,fFitMax, fNparPeak+fNparBgnd);
 658   return;
 659 }
 660
 661 //______________________________________________________________________________
 662 Double_t AliDielectronSignalFunc::PeakBgndFun(const Double_t *x, const Double_t *par) {
 663   //
 664   // merge peak and bgnd functions
 665   //
 666   return (fFuncSignal->EvalPar(x,par) + fFuncBackground->EvalPar(x,par+fNparPeak));
 667 }
 668