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
  49 AliDielectronSignalFunc::AliDielectronSignalFunc() :
  50 AliDielectronSignalBase(),
  51 fFuncSignal(0x0),
  52 fFuncBackground(0x0),
  53 fFuncSigBack(0x0),
  54 fParMass(1),
  55 fParMassWidth(2),
  56 fFitOpt("SMNQE"),
  57 fUseIntegral(kFALSE),
  58 fPolDeg(0),
  59 fChi2Dof(0.0)
  60 {
  61   //
  62   // Default Constructor
  63   //
  64
  65 }
  66
  67 //______________________________________________
  68 AliDielectronSignalFunc::AliDielectronSignalFunc(const char* name, const char* title) :
  69 AliDielectronSignalBase(name, title),
  70 fFuncSignal(0x0),
  71 fFuncBackground(0x0),
  72 fFuncSigBack(0x0),
  73 fParMass(1),
  74 fParMassWidth(2),
  75 fFitOpt("SMNQE"),
  76 fUseIntegral(kFALSE),
  77 fPolDeg(0),
  78 fChi2Dof(0.0)
  79 {
  80   //
  81   // Named Constructor
  82   //
  83 }
  84
  85 //______________________________________________
  86 AliDielectronSignalFunc::~AliDielectronSignalFunc()
  87 {
  88   //
  89   // Default Destructor
  90   //
  91   if(fFuncSignal) delete fFuncSignal;
  92   if(fFuncBackground) delete fFuncBackground;
  93   if(fFuncSigBack) delete fFuncSigBack;
  94 }
  95
  96
  97 //______________________________________________
  98 void AliDielectronSignalFunc::Process(TObjArray * const arrhist)
  99 {
 100   //
 101   // Fit the invariant mass histograms and retrieve the signal and background
 102   //
 103   switch(fMethod) {
 104   case kFittedMC :
 105         ProcessFitIKF(arrhist);
 106         break;
 107
 108   case kFitted :
 109     ProcessFit(arrhist);
 110     break;
 111
 112   case kLikeSign :
 113     ProcessLS(arrhist);
 114     break;
 115
 116   case kEventMixing :
 117     ProcessEM(arrhist);
 118     break;
 119
 120   default :
 121     AliError("Background substraction method not known!");
 122   }
 123 }
 124
 125 //______________________________________________
 126 void AliDielectronSignalFunc::ProcessFitIKF(TObjArray * const arrhist) {
 127   //
 128   // Fit the +- invariant mass distribution only
 129   //
 130   //
 131
 132   const Double_t bigNumber = 100000.;
 133   Double_t chi2ndfm1 = bigNumber;
 134   Double_t ratiom1   = bigNumber;
 135   Double_t chi2ndf   = bigNumber;
 136   Int_t nDP =0;
 137
 138   Int_t maxPolDeg = 8;
 139
 140   fHistDataPM = (TH1F*)(arrhist->At(1))->Clone("histPM");  // +-    SE
 141     if(fRebin>1) fHistDataPM->Rebin(fRebin);
 142
 143   fgHistSimPM = (TH1F*)(arrhist->At(3))->Clone("histPMsim");  // +- mc shape
 144   if (!fgHistSimPM) {
 145         AliFatal("No mc peak shape found at idx 3.");
 146         return;
 147   }
 148   if(fRebin>1) fgHistSimPM->Rebin(fRebin);
 149
 150   // try out the polynomial degrees
 151   for (Int_t iPD=0; iPD<=maxPolDeg; iPD++) {
 152     TH1F *hf1 = (TH1F *) fHistDataPM->Clone(Form("hf1_PD%d",iPD));
 153
 154         FitOneMinv(hf1, fgHistSimPM, iPD);
 155         if (fChi2Dof > 0) chi2ndf = fChi2Dof;
 156     AliInfo(Form("nDP: %d, iPD: %d, chi2ndf: %f", nDP, iPD, chi2ndf));
 157
 158     ratiom1 = TMath::Abs(fChi2Dof - 1);
 159     if (chi2ndfm1 > ratiom1) { // search for the closest to 1.
 160       chi2ndfm1 = ratiom1;
 161       nDP       = iPD;
 162     }
 163   }
 164
 165
 166   // fit again with the best polynomial degree
 167   TH1F *h2 = (TH1F *) fHistDataPM->Clone(Form("h2_PD%d",nDP));
 168
 169   FitOneMinv(h2, fgHistSimPM, nDP);
 170   AliInfo(Form("Best Fit: PD %d, chi^2/ndf %.3f, S/B %.3f",nDP,fChi2Dof,fValues(3)));
 171
 172 }
 173 //______________________________________________
 174 void AliDielectronSignalFunc::FitOneMinv(TH1F *hMinv, TH1F *hSim, Int_t pod) {
 175   //
 176   // main function to fit an inv mass spectrum
 177   //
 178
 179   TObjArray *arrResults = new TObjArray;
 180   arrResults->SetOwner();
 181   arrResults->AddAt(hMinv,0);
 182
 183   // Degree of polynomial
 184   fPolDeg = pod;
 185
 186   // inclusion and exclusion areas (values)
 187   const Double_t kJPMass    = 3.096916;
 188   // inclusion and exclusion areas (bin numbers)
 189   const Int_t binIntLo = hMinv->FindBin(fIntMin);
 190   const Int_t binIntHi = hMinv->FindBin(fIntMax);
 191   // for error calculation
 192   Double_t intAreaEdgeLo = hMinv->GetBinLowEdge(binIntLo);
 193   Double_t intAreaEdgeHi = hMinv->GetBinLowEdge(binIntHi)+hMinv->GetBinWidth(binIntHi);
 194   Double_t norm = (binIntHi-binIntLo)/(fIntMax - fIntMin);
 195
 196   TH1F  *hBfit = (TH1F *) hMinv->Clone(); // for bg only fit (excluding peak region)
 197   TH1F  *hSigF = (TH1F *) hMinv->Clone(); // signal with subtracted bg
 198   TH1F  *hBgrd = (TH1F *) hMinv->Clone(); // bg histogram
 199
 200   hBfit->Reset();
 201   hSigF->Reset();
 202   hBgrd->Reset();
 203
 204   // extract start parameter for the MC signal fit
 205   Double_t bgvalAv = (hMinv->Integral(1,hMinv->GetNbinsX()+1) - hMinv->Integral(binIntLo,binIntHi)) / (hMinv->GetNbinsX()+1 - (binIntHi-binIntLo));
 206   Double_t pkval     = hMinv->GetBinContent(hMinv->FindBin(kJPMass)) - bgvalAv;
 207   Double_t heightMC  = hSim->GetBinContent(hSim->FindBin(kJPMass));
 208   Double_t peakScale = (heightMC > 0. ? pkval/heightMC : 0.0);
 209
 210   Int_t nBgnd = 2 + fPolDeg; // degree + c1st oefficient + higher coefficients
 211   Int_t nMinv = nBgnd + 1;  // bgrd + peakscale
 212
 213   // Create the spectra without peak region
 214   for (Int_t iBin = 0; iBin <= hMinv->GetNbinsX(); iBin++) {
 215     if ((iBin < binIntLo) || (iBin > binIntHi)) {
 216       hBfit->SetBinContent(iBin,hMinv->GetBinContent(iBin));
 217       hBfit->SetBinError(iBin,hMinv->GetBinError(iBin));
 218         }
 219   }
 220
 221
 222   // =======
 223   //   1.
 224   // =======
 225   // Do the fit to the background spectrum
 226   TF1 *fBo = new TF1("bgrd_fit",BgndFun,fFitMin,fFitMax,nBgnd);
 227   for (Int_t iPar=0; iPar<nBgnd; iPar++) {
 228     if (iPar == 0) fBo->FixParameter(0, fPolDeg);
 229         if (iPar == 1) fBo->SetParameter(iPar, bgvalAv);
 230     if (iPar >= 2) fBo->SetParameter(iPar, 0.);
 231   }
 232   hBfit->Fit(fBo,"0qR");
 233   //hBfit->SetNameTitle("bgrd_fit");
 234   arrResults->AddAt(fBo,1);
 235
 236
 237   // =======
 238   //   2.
 239   // =======
 240   // Fit the whole spectrum with peak and background
 241   TF1 *fSB = new TF1("bgrd_peak_fit",MinvFun,fFitMin,fFitMax,nMinv);
 242   fSB->FixParameter(0, fPolDeg);
 243   fSB->SetParameter(1, peakScale);
 244   // copy the polynomial parameters
 245   for (Int_t iPar=0; iPar<=fPolDeg; iPar++)
 246         fSB->SetParameter(2+iPar, fBo->GetParameter(iPar+1));
 247   hMinv->Fit(fSB,"0qR");
 248   arrResults->AddAt(fSB,2);
 249
 250
 251   // =======
 252   //   3.
 253   // =======
 254   // Create the background function
 255   TF1 *fB = new TF1("bgrdOnly_fkt",BgndFun,fFitMin,fFitMax,nBgnd);
 256   fB->FixParameter(0,fPolDeg);
 257   for (Int_t iDeg=0; iDeg<=fPolDeg; iDeg++) {
 258         fB->SetParameter(1+iDeg,fSB->GetParameter(2+iDeg));
 259         fB->SetParError(1+iDeg,fSB->GetParError(2+iDeg));
 260   }
 261   // create background histogram from background function
 262   hBgrd->Eval(fB);
 263   hBgrd->Fit(fB,"0qR");
 264   // calculate the integral and integral error from fit function
 265   Double_t intc = fB->Integral(intAreaEdgeLo, intAreaEdgeHi) * norm;
 266   Double_t inte = fB->IntegralError(intAreaEdgeLo, intAreaEdgeHi) * norm;
 267   arrResults->AddAt(fB,3);
 268
 269   // Fill the background spectrum erros. Use the error from the fit function for the background fB
 270   for (Int_t iBin = 0; iBin <= hBgrd->GetNbinsX(); iBin++) {
 271     Double_t x = hBgrd->GetBinCenter(iBin);
 272     if ((x >= fFitMin) && (x <= fFitMax)) {
 273           Double_t binte = inte / TMath::Sqrt((binIntHi-binIntLo)+1);
 274           hBgrd->SetBinError(iBin,binte);
 275     }
 276   }
 277   arrResults->AddAt(hBgrd,4);
 278
 279   // =======
 280   //   4.
 281   // =======
 282   // Subtract the background
 283   hSigF->Add(hMinv,hBgrd,1.0,-1.0);
 284   for (Int_t iBin = 0; iBin <= hSigF->GetNbinsX(); iBin++) {
 285     Double_t x = hSigF->GetBinCenter(iBin);
 286     if ((x < fFitMin) || (x > fFitMax)) {
 287       hSigF->SetBinContent(iBin,0.0);
 288       hSigF->SetBinError(iBin,0.0);
 289     }
 290   }
 291   hSigF->SetNameTitle("peak_only","");
 292   arrResults->AddAt(hSigF,5);
 293
 294   // =======
 295   //   5.
 296   // =======
 297   // Fit the background-subtracted spectrum
 298   TF1 *fS  = new TF1("peakOnly_fit",PeakFunCB,fFitMin,fFitMax,5);
 299   fS->SetParameters(-.05,1,kJPMass,.003,700);
 300   fS->SetParNames("alpha","n","meanx","sigma","N");
 301   hSigF->Fit(fS,"0qR");
 302   arrResults->AddAt(fS,6);
 303
 304
 305   // connect data members
 306   fFuncSignal     = (TF1*) arrResults->At(6)->Clone();
 307   fFuncBackground = (TF1*) arrResults->At(3)->Clone();
 308   fFuncSigBack    = (TF1*) arrResults->At(2)->Clone();
 309   fHistSignal     = (TH1F*)arrResults->At(5)->Clone();
 310   fHistBackground = (TH1F*)arrResults->At(4)->Clone();
 311
 312   // fit results
 313   Double_t chi2 = fSB->GetChisquare();
 314   Int_t    ndf  = fSB->GetNDF();
 315   fChi2Dof = chi2/ndf;
 316
 317   // signal + signal error
 318   fValues(0) = hSigF->IntegralAndError(binIntLo, binIntHi, fErrors(0));
 319   fValues(1) = intc; // background
 320   fErrors(1) = inte; // background error
 321   // S/B (2) and significance (3)
 322   SetSignificanceAndSOB();
 323   fValues(4) = fS->GetParameter(2); // mass
 324   fErrors(4) = fS->GetParError(2);  // mass error
 325   fValues(5) = fS->GetParameter(3); // mass wdth
 326   fErrors(5) = fS->GetParError(3);  // mass wdth error
 327
 328
 329   delete arrResults;
 330
 331 }
 332 //______________________________________________________________________________
 333 Double_t AliDielectronSignalFunc::BgndFun(const Double_t *x, const Double_t *par) {
 334   // parameters
 335   // [0]:   degree of polynomial
 336   // [1]:   constant polynomial coefficient
 337   // [2]..: higher polynomial coefficients
 338
 339   Int_t    deg   = ((Int_t) par[0]);
 340
 341   Double_t f     = 0.0;
 342   Double_t yy    = 1.0;
 343   Double_t xx    = x[0];
 344
 345   for (Int_t i = 0; i <= deg; i++) {
 346     f  += par[i+1] * yy;
 347     yy *= xx;
 348   }
 349
 350
 351   return f;
 352 }
 353 //______________________________________________________________________________
 354 Double_t AliDielectronSignalFunc::PeakFun(const Double_t *x, const Double_t *par) {
 355   // Fit MC signal shape
 356   // parameters
 357   // [0]:   scale for simpeak
 358
 359   Double_t xx  = x[0];
 360
 361   TH1F *hPeak = fgHistSimPM;
 362   if (!hPeak) {
 363     printf("E-AliDielectronSignalFunc::PeakFun: No histogram for peak fit defined!\n");
 364     return 0.0;
 365   }
 366
 367   Int_t idx = hPeak->FindBin(xx);
 368   if ((idx <= 1) ||
 369       (idx >= hPeak->GetNbinsX())) {
 370     return 0.0;
 371   }
 372
 373   return (par[0] * hPeak->GetBinContent(idx));
 374
 375 }
 376
 377 //______________________________________________________________________________
 378 Double_t AliDielectronSignalFunc::MinvFun(const Double_t *x, const Double_t *par) {
 379   // parameters
 380   // [0]:   degree of polynomial             -> [0]   for BgndFun
 381   // [1]:   scale for simpeak                -> [0]   for PeakFun
 382   // [2]:   constant polynomial coefficient  -> [1]   for BgndFun
 383   // [3]..: higher polynomial coefficients   -> [2].. for BgndFun
 384
 385   Int_t    deg = ((Int_t) par[0]);
 386   Double_t parPK[25], parBG[25];
 387
 388   parBG[0] = par[0]; // degree of polynomial
 389
 390   parPK[0] = par[1]; // MC minv scale
 391   for (Int_t i = 0; i <= deg; i++) parBG[i+1] = par[i+2]; // polynomial coefficients
 392
 393   Double_t peak = PeakFun(x,parPK);
 394   Double_t bgnd = BgndFun(x,parBG);
 395   Double_t f    = peak + bgnd;
 396
 397   return f;
 398 }
 399
 400 //______________________________________________________________________________
 401 Double_t AliDielectronSignalFunc::PeakFunCB(const Double_t *x, const Double_t *par) {
 402   // Crystal Ball function fit
 403
 404   Double_t alpha = par[0];
 405   Double_t     n = par[1];
 406   Double_t meanx = par[2];
 407   Double_t sigma = par[3];
 408   Double_t    nn = par[4];
 409
 410   Double_t a = TMath::Power((n/TMath::Abs(alpha)), n) * TMath::Exp(-.5*alpha*alpha);
 411   Double_t b = n/TMath::Abs(alpha) - TMath::Abs(alpha);
 412
 413   Double_t arg = (x[0] - meanx)/sigma;
 414   Double_t fitval = 0;
 415
 416   if (arg > -1.*alpha) {
 417     fitval = nn * TMath::Exp(-.5*arg*arg);
 418   } else {
 419     fitval = nn * a * TMath::Power((b-arg), (-1*n));
 420   }
 421
 422   return fitval;
 423 }
 424
 425
 426 //______________________________________________
 427 void AliDielectronSignalFunc::ProcessFit(TObjArray * const arrhist) {
 428   //
 429   // Fit the +- invariant mass distribution only
 430   // Here we assume that the combined fit function is a sum of the signal and background functions
 431   //    and that the signal function is always the first term of this sum
 432   //
 433
 434   fHistDataPM = (TH1F*)(arrhist->At(1))->Clone("histPM");  // +-    SE
 435   fHistDataPM->Sumw2();
 436   if(fRebin>1)
 437     fHistDataPM->Rebin(fRebin);
 438
 439   fHistSignal = new TH1F("HistSignal", "Fit substracted signal",
 440                          fHistDataPM->GetXaxis()->GetNbins(),
 441                          fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
 442   fHistBackground = new TH1F("HistBackground", "Fit contribution",
 443                              fHistDataPM->GetXaxis()->GetNbins(),
 444                              fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
 445
 446   // the starting parameters of the fit function and their limits can be tuned
 447   // by the user in its macro
 448   fHistDataPM->Fit(fFuncSigBack, fFitOpt.Data(), "", fFitMin, fFitMax);
 449   TFitResultPtr pmFitPtr = fHistDataPM->Fit(fFuncSigBack, fFitOpt.Data(), "", fFitMin, fFitMax);
 450   //TFitResult *pmFitResult = pmFitPtr.Get(); // used only with TF1Helper
 451   fFuncSignal->SetParameters(fFuncSigBack->GetParameters());
 452   fFuncBackground->SetParameters(fFuncSigBack->GetParameters()+fFuncSignal->GetNpar());
 453
 454   for(Int_t iBin=1; iBin<=fHistDataPM->GetXaxis()->GetNbins(); iBin++) {
 455     Double_t m = fHistDataPM->GetBinCenter(iBin);
 456     Double_t pm = fHistDataPM->GetBinContent(iBin);
 457     Double_t epm = fHistDataPM->GetBinError(iBin);
 458     Double_t bknd = fFuncBackground->Eval(m);
 459     Double_t ebknd = 0;
 460     for(Int_t iPar=fFuncSignal->GetNpar(); iPar<fFuncSigBack->GetNpar(); iPar++) {
 461 /* TF1Helper problem on alien compilation
 462       for(Int_t jPar=iPar; jPar<fFuncSigBack->GetNpar(); jPar++) {
 463         TF1 gradientIpar("gradientIpar",
 464                          ROOT::TF1Helper::TGradientParFunction(iPar-fFuncSignal->GetNpar(),fFuncBackground),0,0,0);
 465         TF1 gradientJpar("gradientJpar",
 466                          ROOT::TF1Helper::TGradientParFunction(jPar-fFuncSignal->GetNpar(),fFuncBackground),0,0,0);
 467         ebknd += pmFitResult->CovMatrix(iPar,jPar)*
 468           gradientIpar.Eval(m)*gradientJpar.Eval(m)*
 469           (iPar==jPar ? 1.0 : 2.0);
 470       }
 471 */
 472     }
 473     Double_t signal = pm-bknd;
 474     Double_t error = TMath::Sqrt(epm*epm+ebknd);
 475     fHistSignal->SetBinContent(iBin, signal);
 476     fHistSignal->SetBinError(iBin, error);
 477     fHistBackground->SetBinContent(iBin, bknd);
 478     fHistBackground->SetBinError(iBin, TMath::Sqrt(ebknd));
 479   }
 480
 481   if(fUseIntegral) {
 482     // signal
 483     fValues(0) = fFuncSignal->Integral(fIntMin, fIntMax)/fHistDataPM->GetBinWidth(1);
 484     fErrors(0) = 0;
 485     for(Int_t iPar=0; iPar<fFuncSignal->GetNpar(); iPar++) {
 486 /* TF1Helper problem on alien compilation
 487       for(Int_t jPar=iPar; jPar<fFuncSignal->GetNpar(); jPar++) {
 488         TF1 gradientIpar("gradientIpar",
 489                          ROOT::TF1Helper::TGradientParFunction(iPar,fFuncSignal),0,0,0);
 490         TF1 gradientJpar("gradientJpar",
 491                          ROOT::TF1Helper::TGradientParFunction(jPar,fFuncSignal),0,0,0);
 492         fErrors(0) += pmFitResult->CovMatrix(iPar,jPar)*
 493           gradientIpar.Integral(fIntMin,fIntMax)*gradientJpar.Integral(fIntMin,fIntMax)*
 494           (iPar==jPar ? 1.0 : 2.0);
 495       }
 496 */
 497     }
 498     // background
 499     fValues(1) = fFuncBackground->Integral(fIntMin, fIntMax)/fHistDataPM->GetBinWidth(1);
 500     fErrors(1) = 0;
 501     for(Int_t iPar=fFuncSignal->GetNpar(); iPar<fFuncSigBack->GetNpar(); iPar++) {
 502 /* TF1Helper problem on alien compilation
 503       for(Int_t jPar=iPar; jPar<fFuncSigBack->GetNpar(); jPar++) {
 504         TF1 gradientIpar("gradientIpar",
 505                          ROOT::TF1Helper::TGradientParFunction(iPar-fFuncSignal->GetNpar(),fFuncBackground),0,0,0);
 506         TF1 gradientJpar("gradientJpar",
 507                          ROOT::TF1Helper::TGradientParFunction(jPar-fFuncSignal->GetNpar(),fFuncBackground),0,0,0);
 508         fErrors(1) += pmFitResult->CovMatrix(iPar,jPar)*
 509           gradientIpar.Integral(fIntMin, fIntMax)*gradientJpar.Integral(fIntMin, fIntMax)*
 510           (iPar==jPar ? 1.0 : 2.0);
 511       }
 512 */
 513     }
 514   }
 515   else {
 516     // signal
 517     fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
 518                                                fHistSignal->FindBin(fIntMax), fErrors(0));
 519     // background
 520     fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
 521       fHistBackground->FindBin(fIntMax),
 522       fErrors(1));
 523   }
 524   // S/B and significance
 525   SetSignificanceAndSOB();
 526   fValues(4) = fFuncSigBack->GetParameter(fParMass);
 527   fErrors(4) = fFuncSigBack->GetParError(fParMass);
 528   fValues(5) = fFuncSigBack->GetParameter(fParMassWidth);
 529   fErrors(5) = fFuncSigBack->GetParError(fParMassWidth);
 530
 531   fProcessed = kTRUE;
 532 }
 533
 534 //______________________________________________
 535 void AliDielectronSignalFunc::ProcessLS(TObjArray * const arrhist) {
 536   //
 537   // Substract background using the like-sign spectrum
 538   //
 539   fHistDataPP = (TH1F*)(arrhist->At(0))->Clone("histPP");  // ++    SE
 540   fHistDataPM = (TH1F*)(arrhist->At(1))->Clone("histPM");  // +-    SE
 541   fHistDataMM = (TH1F*)(arrhist->At(2))->Clone("histMM");  // --    SE
 542   if (fRebin>1) {
 543     fHistDataPP->Rebin(fRebin);
 544     fHistDataPM->Rebin(fRebin);
 545     fHistDataMM->Rebin(fRebin);
 546   }
 547   fHistDataPP->Sumw2();
 548   fHistDataPM->Sumw2();
 549   fHistDataMM->Sumw2();
 550
 551   fHistSignal = new TH1F("HistSignal", "Like-Sign substracted signal",
 552                          fHistDataPM->GetXaxis()->GetNbins(),
 553                          fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
 554   fHistBackground = new TH1F("HistBackground", "Like-sign contribution",
 555                              fHistDataPM->GetXaxis()->GetNbins(),
 556                              fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
 557
 558   // fit the +- mass distribution
 559   fHistDataPM->Fit(fFuncSigBack, fFitOpt.Data(), "", fFitMin, fFitMax);
 560   fHistDataPM->Fit(fFuncSigBack, fFitOpt.Data(), "", fFitMin, fFitMax);
 561   // declare the variables where the like-sign fit results will be stored
 562 //   TFitResult *ppFitResult = 0x0;
 563 //   TFitResult *mmFitResult = 0x0;
 564   // fit the like sign background
 565   TF1 *funcClonePP = (TF1*)fFuncBackground->Clone("funcClonePP");
 566   TF1 *funcCloneMM = (TF1*)fFuncBackground->Clone("funcCloneMM");
 567   fHistDataPP->Fit(funcClonePP, fFitOpt.Data(), "", fFitMin, fFitMax);
 568   fHistDataPP->Fit(funcClonePP, fFitOpt.Data(), "", fFitMin, fFitMax);
 569 //   TFitResultPtr ppFitPtr = fHistDataPP->Fit(funcClonePP, fFitOpt.Data(), "", fFitMin, fFitMax);
 570 //   ppFitResult = ppFitPtr.Get();
 571   fHistDataMM->Fit(funcCloneMM, fFitOpt.Data(), "", fFitMin, fFitMax);
 572   fHistDataMM->Fit(funcCloneMM, fFitOpt.Data(), "", fFitMin, fFitMax);
 573 //   TFitResultPtr mmFitPtr = fHistDataMM->Fit(funcCloneMM, fFitOpt.Data(), "", fFitMin, fFitMax);
 574 //   mmFitResult = mmFitPtr.Get();
 575
 576   for(Int_t iBin=1; iBin<=fHistDataPM->GetXaxis()->GetNbins(); iBin++) {
 577     Double_t m = fHistDataPM->GetBinCenter(iBin);
 578     Double_t pm = fHistDataPM->GetBinContent(iBin);
 579     Double_t pp = funcClonePP->Eval(m);
 580     Double_t mm = funcCloneMM->Eval(m);
 581     Double_t epm = fHistDataPM->GetBinError(iBin);
 582     Double_t epp = 0;
 583     for(Int_t iPar=0; iPar<funcClonePP->GetNpar(); iPar++) {
 584 /* TF1Helper problem on alien compilation
 585       for(Int_t jPar=iPar; jPar<funcClonePP->GetNpar(); jPar++) {
 586         TF1 gradientIpar("gradientIpar",
 587                          ROOT::TF1Helper::TGradientParFunction(iPar,funcClonePP),0,0,0);
 588         TF1 gradientJpar("gradientJpar",
 589                          ROOT::TF1Helper::TGradientParFunction(jPar,funcClonePP),0,0,0);
 590         epp += ppFitResult->CovMatrix(iPar,jPar)*
 591           gradientIpar.Eval(m)*gradientJpar.Eval(m)*
 592           (iPar==jPar ? 1.0 : 2.0);
 593       }
 594 */
 595     }
 596     Double_t emm = 0;
 597     for(Int_t iPar=0; iPar<funcCloneMM->GetNpar(); iPar++) {
 598 /* TF1Helper problem on alien compilation
 599       for(Int_t jPar=iPar; jPar<funcCloneMM->GetNpar(); jPar++) {
 600         TF1 gradientIpar("gradientIpar",
 601                          ROOT::TF1Helper::TGradientParFunction(iPar,funcCloneMM),0,0,0);
 602         TF1 gradientJpar("gradientJpar",
 603                          ROOT::TF1Helper::TGradientParFunction(jPar,funcCloneMM),0,0,0);
 604         emm += mmFitResult->CovMatrix(iPar,jPar)*
 605           gradientIpar.Eval(m)*gradientJpar.Eval(m)*
 606           (iPar==jPar ? 1.0 : 2.0);
 607       }
 608 */
 609     }
 610
 611     Double_t signal = pm-2.0*TMath::Sqrt(pp*mm);
 612     Double_t background = 2.0*TMath::Sqrt(pp*mm);
 613     // error propagation on the signal calculation above
 614     Double_t esignal = TMath::Sqrt(epm*epm+(mm/pp)*epp+(pp/mm)*emm);
 615     Double_t ebackground = TMath::Sqrt((mm/pp)*epp+(pp/mm)*emm);
 616     fHistSignal->SetBinContent(iBin, signal);
 617     fHistSignal->SetBinError(iBin, esignal);
 618     fHistBackground->SetBinContent(iBin, background);
 619     fHistBackground->SetBinError(iBin, ebackground);
 620   }
 621
 622   // signal
 623   fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
 624                                              fHistSignal->FindBin(fIntMax), fErrors(0));
 625   // background
 626   fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
 627                                                  fHistBackground->FindBin(fIntMax),
 628                                                  fErrors(1));
 629   // S/B and significance
 630   SetSignificanceAndSOB();
 631   fValues(4) = fFuncSigBack->GetParameter(fParMass);
 632   fErrors(4) = fFuncSigBack->GetParError(fParMass);
 633   fValues(5) = fFuncSigBack->GetParameter(fParMassWidth);
 634   fErrors(5) = fFuncSigBack->GetParError(fParMassWidth);
 635
 636   fProcessed = kTRUE;
 637 }
 638
 639 //______________________________________________
 640 void AliDielectronSignalFunc::ProcessEM(TObjArray * const arrhist) {
 641   //
 642   // Substract background with the event mixing technique
 643   //
 644   arrhist->GetEntries();   // just to avoid the unused parameter warning
 645   AliError("Event mixing for background substraction method not implemented!");
 646 }
 647
 648 //______________________________________________
 649 void AliDielectronSignalFunc::SetFunctions(TF1 * const combined, TF1 * const sig, TF1 * const back,
 650                                            Int_t parM, Int_t parMres)
 651 {
 652   //
 653   // Set the signal, background functions and combined fit function
 654   // Note: The process method assumes that the first n parameters in the
 655   //       combined fit function correspond to the n parameters of the signal function
 656   //       and the n+1 to n+m parameters to the m parameters of the background function!!!
 657
 658   if (!sig||!back||!combined) {
 659     AliError("Both, signal and background function need to be set!");
 660     return;
 661   }
 662   fFuncSignal=sig;
 663   fFuncBackground=back;
 664   fFuncSigBack=combined;
 665   fParMass=parM;
 666   fParMassWidth=parMres;
 667 }
 668
 669 //______________________________________________
 670 void AliDielectronSignalFunc::SetDefaults(Int_t type)
 671 {
 672   //
 673   // Setup some default functions:
 674   // type = 0: gaus signal + linear background in 2.5 - 4 GeV inv. mass
 675   // type = 1: gaus signal + exponential background in 2.5 - 4 GeV inv. mass
 676   // type = 2: half gaussian, half exponential signal function
 677   // type = 3: Crystal-Ball function
 678   // type = 4: Crystal-Ball signal + exponential background
 679   //
 680
 681   if (type==0){
 682     fFuncSignal=new TF1("DieleSignal","gaus",2.5,4);
 683     fFuncBackground=new TF1("DieleBackground","pol1",2.5,4);
 684     fFuncSigBack=new TF1("DieleCombined","gaus+pol1(3)",2.5,4);
 685
 686     fFuncSigBack->SetParameters(1,3.1,.05,2.5,1);
 687     fFuncSigBack->SetParLimits(0,0,10000000);
 688     fFuncSigBack->SetParLimits(1,3.05,3.15);
 689     fFuncSigBack->SetParLimits(2,.02,.1);
 690   }
 691   else if (type==1){
 692     fFuncSignal=new TF1("DieleSignal","gaus",2.5,4);
 693     fFuncBackground=new TF1("DieleBackground","[0]*exp(-(x-[1])/[2])",2.5,4);
 694     fFuncSigBack=new TF1("DieleCombined","gaus+[3]*exp(-(x-[4])/[5])",2.5,4);
 695
 696     fFuncSigBack->SetParameters(1,3.1,.05,1,2.5,1);
 697     fFuncSigBack->SetParLimits(0,0,10000000);
 698     fFuncSigBack->SetParLimits(1,3.05,3.15);
 699     fFuncSigBack->SetParLimits(2,.02,.1);
 700   }
 701   else if (type==2){
 702     // half gaussian, half exponential signal function
 703     // exponential background
 704     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);
 705     fFuncBackground = new TF1("DieleBackground","[0]*exp(-(x-[1])/[2])+[3]",2.5,4);
 706     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);
 707     fFuncSigBack->SetParameters(1.,3.1,.05,.1,1,2.5,1,0);
 708
 709     fFuncSigBack->SetParLimits(0,0,10000000);
 710     fFuncSigBack->SetParLimits(1,3.05,3.15);
 711     fFuncSigBack->SetParLimits(2,.02,.1);
 712     fFuncSigBack->FixParameter(6,2.5);
 713     fFuncSigBack->FixParameter(7,0);
 714   }
 715 }
 716
 717
 718 //______________________________________________
 719 void AliDielectronSignalFunc::Draw(const Option_t* option)
 720 {
 721   //
 722   // Draw the fitted function
 723   //
 724
 725   TString drawOpt(option);
 726   drawOpt.ToLower();
 727
 728   Bool_t optStat=drawOpt.Contains("stat");
 729
 730   fFuncSigBack->SetNpx(200);
 731   fFuncSigBack->SetRange(fIntMin,fIntMax);
 732   fFuncBackground->SetNpx(200);
 733   fFuncBackground->SetRange(fIntMin,fIntMax);
 734
 735   TGraph *grSig=new TGraph(fFuncSigBack);
 736   grSig->SetFillColor(kGreen);
 737   grSig->SetFillStyle(3001);
 738
 739   TGraph *grBack=new TGraph(fFuncBackground);
 740   grBack->SetFillColor(kRed);
 741   grBack->SetFillStyle(3001);
 742
 743   grSig->SetPoint(0,grBack->GetX()[0],grBack->GetY()[0]);
 744   grSig->SetPoint(grSig->GetN()-1,grBack->GetX()[grBack->GetN()-1],grBack->GetY()[grBack->GetN()-1]);
 745
 746   grBack->SetPoint(0,grBack->GetX()[0],0.);
 747   grBack->SetPoint(grBack->GetN()-1,grBack->GetX()[grBack->GetN()-1],0.);
 748
 749   fFuncSigBack->SetRange(fFitMin,fFitMax);
 750   fFuncBackground->SetRange(fFitMin,fFitMax);
 751
 752   if (!drawOpt.Contains("same")){
 753     if (fHistDataPM){
 754       fHistDataPM->Draw();
 755       grSig->Draw("f");
 756     } else {
 757       grSig->Draw("af");
 758     }
 759   } else {
 760     grSig->Draw("f");
 761   }
 762   if(fMethod==kFitted || fMethod==kFittedMC) grBack->Draw("f");
 763   fFuncSigBack->Draw("same");
 764   fFuncSigBack->SetLineWidth(2);
 765   if(fMethod==kLikeSign) {
 766     fHistDataPP->SetLineWidth(2);
 767     fHistDataPP->SetLineColor(6);
 768     fHistDataPP->Draw("same");
 769     fHistDataMM->SetLineWidth(2);
 770     fHistDataMM->SetLineColor(8);
 771     fHistDataMM->Draw("same");
 772   }
 773
 774   if(fMethod==kFitted || fMethod==kFittedMC)
 775     fFuncBackground->Draw("same");
 776
 777   if (optStat) DrawStats();
 778
 779 }