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("F-AliDielectronSignalFunc::PeakFun: No histogram for peak fit defined!\n");
 364   }
 365
 366   Int_t idx = hPeak->FindBin(xx);
 367   if ((idx <= 1) ||
 368       (idx >= hPeak->GetNbinsX())) {
 369     return 0.0;
 370   }
 371
 372   return (par[0] * hPeak->GetBinContent(idx));
 373
 374 }
 375
 376 //______________________________________________________________________________
 377 Double_t AliDielectronSignalFunc::MinvFun(const Double_t *x, const Double_t *par) {
 378   // parameters
 379   // [0]:   degree of polynomial             -> [0]   for BgndFun
 380   // [1]:   scale for simpeak                -> [0]   for PeakFun
 381   // [2]:   constant polynomial coefficient  -> [1]   for BgndFun
 382   // [3]..: higher polynomial coefficients   -> [2].. for BgndFun
 383
 384   Int_t    deg = ((Int_t) par[0]);
 385   Double_t parPK[25], parBG[25];
 386
 387   parBG[0] = par[0]; // degree of polynomial
 388
 389   parPK[0] = par[1]; // MC minv scale
 390   for (Int_t i = 0; i <= deg; i++) parBG[i+1] = par[i+2]; // polynomial coefficients
 391
 392   Double_t peak = PeakFun(x,parPK);
 393   Double_t bgnd = BgndFun(x,parBG);
 394   Double_t f    = peak + bgnd;
 395
 396   return f;
 397 }
 398
 399 //______________________________________________________________________________
 400 Double_t AliDielectronSignalFunc::PeakFunCB(const Double_t *x, const Double_t *par) {
 401   // Crystal Ball function fit
 402
 403   Double_t alpha = par[0];
 404   Double_t     n = par[1];
 405   Double_t meanx = par[2];
 406   Double_t sigma = par[3];
 407   Double_t    nn = par[4];
 408
 409   Double_t a = TMath::Power((n/TMath::Abs(alpha)), n) * TMath::Exp(-.5*alpha*alpha);
 410   Double_t b = n/TMath::Abs(alpha) - TMath::Abs(alpha);
 411
 412   Double_t arg = (x[0] - meanx)/sigma;
 413   Double_t fitval = 0;
 414
 415   if (arg > -1.*alpha) {
 416     fitval = nn * TMath::Exp(-.5*arg*arg);
 417   } else {
 418     fitval = nn * a * TMath::Power((b-arg), (-1*n));
 419   }
 420
 421   return fitval;
 422 }
 423
 424
 425 //______________________________________________
 426 void AliDielectronSignalFunc::ProcessFit(TObjArray * const arrhist) {
 427   //
 428   // Fit the +- invariant mass distribution only
 429   // Here we assume that the combined fit function is a sum of the signal and background functions
 430   //    and that the signal function is always the first term of this sum
 431   //
 432
 433   fHistDataPM = (TH1F*)(arrhist->At(1))->Clone("histPM");  // +-    SE
 434   fHistDataPM->Sumw2();
 435   if(fRebin>1)
 436     fHistDataPM->Rebin(fRebin);
 437
 438   fHistSignal = new TH1F("HistSignal", "Fit substracted signal",
 439                          fHistDataPM->GetXaxis()->GetNbins(),
 440                          fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
 441   fHistBackground = new TH1F("HistBackground", "Fit contribution",
 442                              fHistDataPM->GetXaxis()->GetNbins(),
 443                              fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
 444
 445   // the starting parameters of the fit function and their limits can be tuned
 446   // by the user in its macro
 447   fHistDataPM->Fit(fFuncSigBack, fFitOpt.Data(), "", fFitMin, fFitMax);
 448   TFitResultPtr pmFitPtr = fHistDataPM->Fit(fFuncSigBack, fFitOpt.Data(), "", fFitMin, fFitMax);
 449   //TFitResult *pmFitResult = pmFitPtr.Get(); // used only with TF1Helper
 450   fFuncSignal->SetParameters(fFuncSigBack->GetParameters());
 451   fFuncBackground->SetParameters(fFuncSigBack->GetParameters()+fFuncSignal->GetNpar());
 452
 453   for(Int_t iBin=1; iBin<=fHistDataPM->GetXaxis()->GetNbins(); iBin++) {
 454     Double_t m = fHistDataPM->GetBinCenter(iBin);
 455     Double_t pm = fHistDataPM->GetBinContent(iBin);
 456     Double_t epm = fHistDataPM->GetBinError(iBin);
 457     Double_t bknd = fFuncBackground->Eval(m);
 458     Double_t ebknd = 0;
 459     for(Int_t iPar=fFuncSignal->GetNpar(); iPar<fFuncSigBack->GetNpar(); iPar++) {
 460 /* TF1Helper problem on alien compilation
 461       for(Int_t jPar=iPar; jPar<fFuncSigBack->GetNpar(); jPar++) {
 462         TF1 gradientIpar("gradientIpar",
 463                          ROOT::TF1Helper::TGradientParFunction(iPar-fFuncSignal->GetNpar(),fFuncBackground),0,0,0);
 464         TF1 gradientJpar("gradientJpar",
 465                          ROOT::TF1Helper::TGradientParFunction(jPar-fFuncSignal->GetNpar(),fFuncBackground),0,0,0);
 466         ebknd += pmFitResult->CovMatrix(iPar,jPar)*
 467           gradientIpar.Eval(m)*gradientJpar.Eval(m)*
 468           (iPar==jPar ? 1.0 : 2.0);
 469       }
 470 */
 471     }
 472     Double_t signal = pm-bknd;
 473     Double_t error = TMath::Sqrt(epm*epm+ebknd);
 474     fHistSignal->SetBinContent(iBin, signal);
 475     fHistSignal->SetBinError(iBin, error);
 476     fHistBackground->SetBinContent(iBin, bknd);
 477     fHistBackground->SetBinError(iBin, TMath::Sqrt(ebknd));
 478   }
 479
 480   if(fUseIntegral) {
 481     // signal
 482     fValues(0) = fFuncSignal->Integral(fIntMin, fIntMax)/fHistDataPM->GetBinWidth(1);
 483     fErrors(0) = 0;
 484     for(Int_t iPar=0; iPar<fFuncSignal->GetNpar(); iPar++) {
 485 /* TF1Helper problem on alien compilation
 486       for(Int_t jPar=iPar; jPar<fFuncSignal->GetNpar(); jPar++) {
 487         TF1 gradientIpar("gradientIpar",
 488                          ROOT::TF1Helper::TGradientParFunction(iPar,fFuncSignal),0,0,0);
 489         TF1 gradientJpar("gradientJpar",
 490                          ROOT::TF1Helper::TGradientParFunction(jPar,fFuncSignal),0,0,0);
 491         fErrors(0) += pmFitResult->CovMatrix(iPar,jPar)*
 492           gradientIpar.Integral(fIntMin,fIntMax)*gradientJpar.Integral(fIntMin,fIntMax)*
 493           (iPar==jPar ? 1.0 : 2.0);
 494       }
 495 */
 496     }
 497     // background
 498     fValues(1) = fFuncBackground->Integral(fIntMin, fIntMax)/fHistDataPM->GetBinWidth(1);
 499     fErrors(1) = 0;
 500     for(Int_t iPar=fFuncSignal->GetNpar(); iPar<fFuncSigBack->GetNpar(); iPar++) {
 501 /* TF1Helper problem on alien compilation
 502       for(Int_t jPar=iPar; jPar<fFuncSigBack->GetNpar(); jPar++) {
 503         TF1 gradientIpar("gradientIpar",
 504                          ROOT::TF1Helper::TGradientParFunction(iPar-fFuncSignal->GetNpar(),fFuncBackground),0,0,0);
 505         TF1 gradientJpar("gradientJpar",
 506                          ROOT::TF1Helper::TGradientParFunction(jPar-fFuncSignal->GetNpar(),fFuncBackground),0,0,0);
 507         fErrors(1) += pmFitResult->CovMatrix(iPar,jPar)*
 508           gradientIpar.Integral(fIntMin, fIntMax)*gradientJpar.Integral(fIntMin, fIntMax)*
 509           (iPar==jPar ? 1.0 : 2.0);
 510       }
 511 */
 512     }
 513   }
 514   else {
 515     // signal
 516     fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
 517                                                fHistSignal->FindBin(fIntMax), fErrors(0));
 518     // background
 519     fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
 520       fHistBackground->FindBin(fIntMax),
 521       fErrors(1));
 522   }
 523   // S/B and significance
 524   SetSignificanceAndSOB();
 525   fValues(4) = fFuncSigBack->GetParameter(fParMass);
 526   fErrors(4) = fFuncSigBack->GetParError(fParMass);
 527   fValues(5) = fFuncSigBack->GetParameter(fParMassWidth);
 528   fErrors(5) = fFuncSigBack->GetParError(fParMassWidth);
 529
 530   fProcessed = kTRUE;
 531 }
 532
 533 //______________________________________________
 534 void AliDielectronSignalFunc::ProcessLS(TObjArray * const arrhist) {
 535   //
 536   // Substract background using the like-sign spectrum
 537   //
 538   fHistDataPP = (TH1F*)(arrhist->At(0))->Clone("histPP");  // ++    SE
 539   fHistDataPM = (TH1F*)(arrhist->At(1))->Clone("histPM");  // +-    SE
 540   fHistDataMM = (TH1F*)(arrhist->At(2))->Clone("histMM");  // --    SE
 541   if (fRebin>1) {
 542     fHistDataPP->Rebin(fRebin);
 543     fHistDataPM->Rebin(fRebin);
 544     fHistDataMM->Rebin(fRebin);
 545   }
 546   fHistDataPP->Sumw2();
 547   fHistDataPM->Sumw2();
 548   fHistDataMM->Sumw2();
 549
 550   fHistSignal = new TH1F("HistSignal", "Like-Sign substracted signal",
 551                          fHistDataPM->GetXaxis()->GetNbins(),
 552                          fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
 553   fHistBackground = new TH1F("HistBackground", "Like-sign contribution",
 554                              fHistDataPM->GetXaxis()->GetNbins(),
 555                              fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
 556
 557   // fit the +- mass distribution
 558   fHistDataPM->Fit(fFuncSigBack, fFitOpt.Data(), "", fFitMin, fFitMax);
 559   fHistDataPM->Fit(fFuncSigBack, fFitOpt.Data(), "", fFitMin, fFitMax);
 560   // declare the variables where the like-sign fit results will be stored
 561 //   TFitResult *ppFitResult = 0x0;
 562 //   TFitResult *mmFitResult = 0x0;
 563   // fit the like sign background
 564   TF1 *funcClonePP = (TF1*)fFuncBackground->Clone("funcClonePP");
 565   TF1 *funcCloneMM = (TF1*)fFuncBackground->Clone("funcCloneMM");
 566   fHistDataPP->Fit(funcClonePP, fFitOpt.Data(), "", fFitMin, fFitMax);
 567   fHistDataPP->Fit(funcClonePP, fFitOpt.Data(), "", fFitMin, fFitMax);
 568 //   TFitResultPtr ppFitPtr = fHistDataPP->Fit(funcClonePP, fFitOpt.Data(), "", fFitMin, fFitMax);
 569 //   ppFitResult = ppFitPtr.Get();
 570   fHistDataMM->Fit(funcCloneMM, fFitOpt.Data(), "", fFitMin, fFitMax);
 571   fHistDataMM->Fit(funcCloneMM, fFitOpt.Data(), "", fFitMin, fFitMax);
 572 //   TFitResultPtr mmFitPtr = fHistDataMM->Fit(funcCloneMM, fFitOpt.Data(), "", fFitMin, fFitMax);
 573 //   mmFitResult = mmFitPtr.Get();
 574
 575   for(Int_t iBin=1; iBin<=fHistDataPM->GetXaxis()->GetNbins(); iBin++) {
 576     Double_t m = fHistDataPM->GetBinCenter(iBin);
 577     Double_t pm = fHistDataPM->GetBinContent(iBin);
 578     Double_t pp = funcClonePP->Eval(m);
 579     Double_t mm = funcCloneMM->Eval(m);
 580     Double_t epm = fHistDataPM->GetBinError(iBin);
 581     Double_t epp = 0;
 582     for(Int_t iPar=0; iPar<funcClonePP->GetNpar(); iPar++) {
 583 /* TF1Helper problem on alien compilation
 584       for(Int_t jPar=iPar; jPar<funcClonePP->GetNpar(); jPar++) {
 585         TF1 gradientIpar("gradientIpar",
 586                          ROOT::TF1Helper::TGradientParFunction(iPar,funcClonePP),0,0,0);
 587         TF1 gradientJpar("gradientJpar",
 588                          ROOT::TF1Helper::TGradientParFunction(jPar,funcClonePP),0,0,0);
 589         epp += ppFitResult->CovMatrix(iPar,jPar)*
 590           gradientIpar.Eval(m)*gradientJpar.Eval(m)*
 591           (iPar==jPar ? 1.0 : 2.0);
 592       }
 593 */
 594     }
 595     Double_t emm = 0;
 596     for(Int_t iPar=0; iPar<funcCloneMM->GetNpar(); iPar++) {
 597 /* TF1Helper problem on alien compilation
 598       for(Int_t jPar=iPar; jPar<funcCloneMM->GetNpar(); jPar++) {
 599         TF1 gradientIpar("gradientIpar",
 600                          ROOT::TF1Helper::TGradientParFunction(iPar,funcCloneMM),0,0,0);
 601         TF1 gradientJpar("gradientJpar",
 602                          ROOT::TF1Helper::TGradientParFunction(jPar,funcCloneMM),0,0,0);
 603         emm += mmFitResult->CovMatrix(iPar,jPar)*
 604           gradientIpar.Eval(m)*gradientJpar.Eval(m)*
 605           (iPar==jPar ? 1.0 : 2.0);
 606       }
 607 */
 608     }
 609
 610     Double_t signal = pm-2.0*TMath::Sqrt(pp*mm);
 611     Double_t background = 2.0*TMath::Sqrt(pp*mm);
 612     // error propagation on the signal calculation above
 613     Double_t esignal = TMath::Sqrt(epm*epm+(mm/pp)*epp+(pp/mm)*emm);
 614     Double_t ebackground = TMath::Sqrt((mm/pp)*epp+(pp/mm)*emm);
 615     fHistSignal->SetBinContent(iBin, signal);
 616     fHistSignal->SetBinError(iBin, esignal);
 617     fHistBackground->SetBinContent(iBin, background);
 618     fHistBackground->SetBinError(iBin, ebackground);
 619   }
 620
 621   // signal
 622   fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
 623                                              fHistSignal->FindBin(fIntMax), fErrors(0));
 624   // background
 625   fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
 626                                                  fHistBackground->FindBin(fIntMax),
 627                                                  fErrors(1));
 628   // S/B and significance
 629   SetSignificanceAndSOB();
 630   fValues(4) = fFuncSigBack->GetParameter(fParMass);
 631   fErrors(4) = fFuncSigBack->GetParError(fParMass);
 632   fValues(5) = fFuncSigBack->GetParameter(fParMassWidth);
 633   fErrors(5) = fFuncSigBack->GetParError(fParMassWidth);
 634
 635   fProcessed = kTRUE;
 636 }
 637
 638 //______________________________________________
 639 void AliDielectronSignalFunc::ProcessEM(TObjArray * const arrhist) {
 640   //
 641   // Substract background with the event mixing technique
 642   //
 643   arrhist->GetEntries();   // just to avoid the unused parameter warning
 644   AliError("Event mixing for background substraction method not implemented!");
 645 }
 646
 647 //______________________________________________
 648 void AliDielectronSignalFunc::SetFunctions(TF1 * const combined, TF1 * const sig, TF1 * const back,
 649                                            Int_t parM, Int_t parMres)
 650 {
 651   //
 652   // Set the signal, background functions and combined fit function
 653   // Note: The process method assumes that the first n parameters in the
 654   //       combined fit function correspond to the n parameters of the signal function
 655   //       and the n+1 to n+m parameters to the m parameters of the background function!!!
 656
 657   if (!sig||!back||!combined) {
 658     AliError("Both, signal and background function need to be set!");
 659     return;
 660   }
 661   fFuncSignal=sig;
 662   fFuncBackground=back;
 663   fFuncSigBack=combined;
 664   fParMass=parM;
 665   fParMassWidth=parMres;
 666 }
 667
 668 //______________________________________________
 669 void AliDielectronSignalFunc::SetDefaults(Int_t type)
 670 {
 671   //
 672   // Setup some default functions:
 673   // type = 0: gaus signal + linear background in 2.5 - 4 GeV inv. mass
 674   // type = 1: gaus signal + exponential background in 2.5 - 4 GeV inv. mass
 675   // type = 2: half gaussian, half exponential signal function
 676   // type = 3: Crystal-Ball function
 677   // type = 4: Crystal-Ball signal + exponential background
 678   //
 679
 680   if (type==0){
 681     fFuncSignal=new TF1("DieleSignal","gaus",2.5,4);
 682     fFuncBackground=new TF1("DieleBackground","pol1",2.5,4);
 683     fFuncSigBack=new TF1("DieleCombined","gaus+pol1(3)",2.5,4);
 684
 685     fFuncSigBack->SetParameters(1,3.1,.05,2.5,1);
 686     fFuncSigBack->SetParLimits(0,0,10000000);
 687     fFuncSigBack->SetParLimits(1,3.05,3.15);
 688     fFuncSigBack->SetParLimits(2,.02,.1);
 689   }
 690   else if (type==1){
 691     fFuncSignal=new TF1("DieleSignal","gaus",2.5,4);
 692     fFuncBackground=new TF1("DieleBackground","[0]*exp(-(x-[1])/[2])",2.5,4);
 693     fFuncSigBack=new TF1("DieleCombined","gaus+[3]*exp(-(x-[4])/[5])",2.5,4);
 694
 695     fFuncSigBack->SetParameters(1,3.1,.05,1,2.5,1);
 696     fFuncSigBack->SetParLimits(0,0,10000000);
 697     fFuncSigBack->SetParLimits(1,3.05,3.15);
 698     fFuncSigBack->SetParLimits(2,.02,.1);
 699   }
 700   else if (type==2){
 701     // half gaussian, half exponential signal function
 702     // exponential background
 703     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);
 704     fFuncBackground = new TF1("DieleBackground","[0]*exp(-(x-[1])/[2])+[3]",2.5,4);
 705     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);
 706     fFuncSigBack->SetParameters(1.,3.1,.05,.1,1,2.5,1,0);
 707
 708     fFuncSigBack->SetParLimits(0,0,10000000);
 709     fFuncSigBack->SetParLimits(1,3.05,3.15);
 710     fFuncSigBack->SetParLimits(2,.02,.1);
 711     fFuncSigBack->FixParameter(6,2.5);
 712     fFuncSigBack->FixParameter(7,0);
 713   }
 714 }
 715
 716
 717 //______________________________________________
 718 void AliDielectronSignalFunc::Draw(const Option_t* option)
 719 {
 720   //
 721   // Draw the fitted function
 722   //
 723
 724   TString drawOpt(option);
 725   drawOpt.ToLower();
 726
 727   Bool_t optStat=drawOpt.Contains("stat");
 728
 729   fFuncSigBack->SetNpx(200);
 730   fFuncSigBack->SetRange(fIntMin,fIntMax);
 731   fFuncBackground->SetNpx(200);
 732   fFuncBackground->SetRange(fIntMin,fIntMax);
 733
 734   TGraph *grSig=new TGraph(fFuncSigBack);
 735   grSig->SetFillColor(kGreen);
 736   grSig->SetFillStyle(3001);
 737
 738   TGraph *grBack=new TGraph(fFuncBackground);
 739   grBack->SetFillColor(kRed);
 740   grBack->SetFillStyle(3001);
 741
 742   grSig->SetPoint(0,grBack->GetX()[0],grBack->GetY()[0]);
 743   grSig->SetPoint(grSig->GetN()-1,grBack->GetX()[grBack->GetN()-1],grBack->GetY()[grBack->GetN()-1]);
 744
 745   grBack->SetPoint(0,grBack->GetX()[0],0.);
 746   grBack->SetPoint(grBack->GetN()-1,grBack->GetX()[grBack->GetN()-1],0.);
 747
 748   fFuncSigBack->SetRange(fFitMin,fFitMax);
 749   fFuncBackground->SetRange(fFitMin,fFitMax);
 750
 751   if (!drawOpt.Contains("same")){
 752     if (fHistDataPM){
 753       fHistDataPM->Draw();
 754       grSig->Draw("f");
 755     } else {
 756       grSig->Draw("af");
 757     }
 758   } else {
 759     grSig->Draw("f");
 760   }
 761   if(fMethod==kFitted || fMethod==kFittedMC) grBack->Draw("f");
 762   fFuncSigBack->Draw("same");
 763   fFuncSigBack->SetLineWidth(2);
 764   if(fMethod==kLikeSign) {
 765     fHistDataPP->SetLineWidth(2);
 766     fHistDataPP->SetLineColor(6);
 767     fHistDataPP->Draw("same");
 768     fHistDataMM->SetLineWidth(2);
 769     fHistDataMM->SetLineColor(8);
 770     fHistDataMM->Draw("same");
 771   }
 772
 773   if(fMethod==kFitted || fMethod==kFittedMC)
 774     fFuncBackground->Draw("same");
 775
 776   if (optStat) DrawStats();
 777
 778 }