]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGDQ/dielectron/AliDielectronSignalFunc.cxx
Update to FindFASTJET.cmake; now accepts also for -DFASTJET= option
[u/mrichter/AliRoot.git] / PWGDQ / dielectron / AliDielectronSignalFunc.cxx
index 5c68176cf7d9efcd0e6abddc11933ba9bd028b89..4f7c648c2d07ed3b9f09c818f02e0e880a9e3e85 100644 (file)
 //                                                                       //
 //                                                                       //
 /*
-Dielectron signal extraction class using functions as input.
 
-A function to describe the signal as well as one to describe the background
-has to be deployed by the user. Alternatively on of the default implementaions
-can be used.
+  Class used for extracting the signal from an invariant mass spectrum.
+  It implements the AliDielectronSignalBase and -Ext classes and it uses user provided
+  functions to fit the spectrum with a combined signa+background fit.
+  Used invariant mass spectra are provided via an array of histograms. There are serveral method
+  to estimate the background and to extract the raw yield from the background subtracted spectra.
+
+  Example usage:
+
+  AliDielectronSignalFunc *sig = new AliDielectronSignalFunc();
+
+
+  1) invariant mass input spectra
+
+  1.1) Assuming a AliDielectronCF container as data format (check class for more details)
+  AliDielectronCFdraw *cf = new AliDielectronCFdraw("path/to/the/output/file.root");
+  TObjArray *arrHists = cf->CollectMinvProj(cf->FindStep("Config"));
+
+  1.2) Assuming a AliDielectronHF grid as data format (check class for more details)
+  AliDielectronHFhelper *hf = new AliDielectronHFhelper("path/to/the/output/file.root", "ConfigName");
+  TObjArray *arrHists = hf->CollectHistos(AliDielectronVarManager::kM);
+
+  1.3) Assuming a single histograms
+  TObjArray *histoArray = new TObjArray();
+  arrHists->Add(signalPP);            // add the spectrum histograms to the array
+  arrHists->Add(signalPM);            // the order is important !!!
+  arrHists->Add(signalMM);
+
+
+  2) background estimation
+
+  2.1) set the method for the background estimation (methods can be found in AliDielectronSignalBase)
+  sig->SetMethod(AliDielectronSignalBase::kFitted);
+  2.2) rebin the spectras if needed
+  //  sig->SetRebin(2);
+  2.3) add any background function you like
+  TF1 *fB = new TF1("fitBgrd","pol3",minFit,maxFit);
+
+
+  3) configure the signal extraction
+
+  3.1) chose one of the signal functions (MCshape, CrystalBall, Gauss)
+  TF1 *fS = new TF1("fitSign",AliDielectronSignalFunc::PeakFunCB,minFit,maxFit,5); // has 5 parameters
+  //  TF1 *fS = new TF1("fitSign",AliDielectronSignalFunc::PeakFunGaus,minFit,maxFit,3); // has 3 parameters
+  //  sig->SetMCSignalShape(hMCsign);
+  //  TF1 *fS = new TF1("fitSign",AliDielectronSignalFunc::PeakFunMC,minFit,maxFit,1); // requires a MC shape
+  3.2) set the method for the signal extraction (methods can be found in AliDielectronSignalBase)
+  depending on the method serveral inputs are needed (e.g. MC shape, PDG code of the particle of interest)
+  //  sig->SetParticleOfInterest(443); //default is jpsi
+  //  sig->SetMCSignalShape(signalMC);
+  //  sig->SetIntegralRange(minInt, maxInt);
+  sig->SetExtractionMethod(AliDielectronSignal::BinCounting); // this is the default
+
+
+  4) combined fit of bgrd+signal
+
+  4.1) combine the two functions
+  sig->CombineFunc(fS,fB);
+  4.2) apply fitting ranges and the fit options
+  sig->SetFitRange(minFit, maxFit);
+  sig->SetFitOption("NR");
+
+
+  5) start the processing
+
+  sig->Process(arrHists);
+  sig->Print(""); // print values and errors extracted
+
+
+  6) access the spectra and values created
+
+  6.1) standard spectra as provided filled in AliDielectronSignalExt
+  TH1F *hsign = (TH1F*) sig->GetUnlikeSignHistogram();  // same as the input (rebinned)
+  TH1F *hbgrd = (TH1F*) sig->GetBackgroundHistogram();  // filled histogram with fitBgrd
+  TH1F *hextr = (TH1F*) sig->GetSignalHistogram();      // after backgound extraction (rebinned)
+  TObject *oPeak = (TObject*) sig->GetPeakShape();      // can be a TF1 or TH1 depending on the method
+  6.2) flow spectra
+  TF1 *fFitSign  = sig->GetCombinedFunction();                // combined fit function
+  TF1 *fFitExtr  = sig->GetSignalFunction();                  // signal function
+  TF1 *fFitBgrd  = sig->GetBackgroundFunction();              // background function
+  6.3) access the extracted values and errors
+  sig->GetValues();     or GetErrors();                 // yield extraction
 
 */
 //                                                                       //
@@ -44,18 +121,21 @@ can be used.
 
 ClassImp(AliDielectronSignalFunc)
 
-TH1F* AliDielectronSignalFunc::fgHistSimPM=0x0;
+//TH1F* AliDielectronSignalFunc::fgHistSimPM=0x0;
+TF1*  AliDielectronSignalFunc::fFuncSignal=0x0;
+TF1*  AliDielectronSignalFunc::fFuncBackground=0x0;
+Int_t AliDielectronSignalFunc::fNparPeak=0;
+Int_t AliDielectronSignalFunc::fNparBgnd=0;
+
 
 AliDielectronSignalFunc::AliDielectronSignalFunc() :
-AliDielectronSignalBase(),
-fFuncSignal(0x0),
-fFuncBackground(0x0),
+AliDielectronSignalExt(),
 fFuncSigBack(0x0),
 fParMass(1),
 fParMassWidth(2),
 fFitOpt("SMNQE"),
 fUseIntegral(kFALSE),
-fPolDeg(0),
+fDof(0),
 fChi2Dof(0.0)
 {
   //
@@ -66,15 +146,13 @@ fChi2Dof(0.0)
 
 //______________________________________________
 AliDielectronSignalFunc::AliDielectronSignalFunc(const char* name, const char* title) :
-AliDielectronSignalBase(name, title),
-fFuncSignal(0x0),
-fFuncBackground(0x0),
+AliDielectronSignalExt(name, title),
 fFuncSigBack(0x0),
 fParMass(1),
 fParMassWidth(2),
 fFitOpt("SMNQE"),
 fUseIntegral(kFALSE),
-fPolDeg(0),
+fDof(0),
 fChi2Dof(0.0)
 {
   //
@@ -88,8 +166,6 @@ AliDielectronSignalFunc::~AliDielectronSignalFunc()
   //
   // Default Destructor
   //
-  if(fFuncSignal) delete fFuncSignal;
-  if(fFuncBackground) delete fFuncBackground;
   if(fFuncSigBack) delete fFuncSigBack;
 }
 
@@ -101,257 +177,39 @@ void AliDielectronSignalFunc::Process(TObjArray * const arrhist)
   // Fit the invariant mass histograms and retrieve the signal and background
   //
   switch(fMethod) {
-  case kFittedMC :
-       ProcessFitIKF(arrhist);
-       break;
-                 
   case kFitted :
     ProcessFit(arrhist);
     break;
-    
+
+  case kLikeSignFit :
+    ProcessFitLS(arrhist);
+    break;
+
+  case kEventMixingFit :
+    ProcessFitEM(arrhist);
+    break;
+
   case kLikeSign :
-    ProcessLS(arrhist);
+  case kLikeSignArithm :
+  case kLikeSignRcorr:
+  case kLikeSignArithmRcorr:
+    ProcessLS(arrhist);    // process like-sign subtraction method
     break;
-    
+
   case kEventMixing :
-    ProcessEM(arrhist);
+    ProcessEM(arrhist);    // process event mixing method
     break;
-    
+
+  case kRotation:
+    ProcessRotation(arrhist);
+    break;
+
   default :
     AliError("Background substraction method not known!");
   }
 }
-
-//______________________________________________
-void AliDielectronSignalFunc::ProcessFitIKF(TObjArray * const arrhist) {
-  //
-  // Fit the +- invariant mass distribution only
-  // 
-  //
-  
-  const Double_t bigNumber = 100000.;
-  Double_t chi2ndfm1 = bigNumber;
-  Double_t ratiom1   = bigNumber;
-  Double_t chi2ndf   = bigNumber;
-  Int_t nDP =0;
-  
-  Int_t maxPolDeg = 8;
-      
-  fHistDataPM = (TH1F*)(arrhist->At(1))->Clone("histPM");  // +-    SE
-    if(fRebin>1) fHistDataPM->Rebin(fRebin);
-  
-  fgHistSimPM = (TH1F*)(arrhist->At(3))->Clone("histPMsim");  // +- mc shape
-  if (!fgHistSimPM) {
-       AliFatal("No mc peak shape found at idx 3.");
-       return;
-  }
-  if(fRebin>1) fgHistSimPM->Rebin(fRebin);
-  
-  // try out the polynomial degrees
-  for (Int_t iPD=0; iPD<=maxPolDeg; iPD++) {
-    TH1F *hf1 = (TH1F *) fHistDataPM->Clone(Form("hf1_PD%d",iPD));
-       
-       FitOneMinv(hf1, fgHistSimPM, iPD);
-       if (fChi2Dof > 0) chi2ndf = fChi2Dof;
-    AliInfo(Form("nDP: %d, iPD: %d, chi2ndf: %f", nDP, iPD, chi2ndf));
-       
-    ratiom1 = TMath::Abs(fChi2Dof - 1);
-    if (chi2ndfm1 > ratiom1) { // search for the closest to 1.
-      chi2ndfm1 = ratiom1;
-      nDP       = iPD;
-    }
-  }
-    
-  
-  // fit again with the best polynomial degree
-  TH1F *h2 = (TH1F *) fHistDataPM->Clone(Form("h2_PD%d",nDP));
-  
-  FitOneMinv(h2, fgHistSimPM, nDP);
-  AliInfo(Form("Best Fit: PD %d, chi^2/ndf %.3f, S/B %.3f",nDP,fChi2Dof,fValues(3)));
-    
-}
-//______________________________________________
-void AliDielectronSignalFunc::FitOneMinv(TH1F *hMinv, TH1F *hSim, Int_t pod) {
-  //
-  // main function to fit an inv mass spectrum
-  //
-  
-  TObjArray *arrResults = new TObjArray;
-  arrResults->SetOwner();
-  arrResults->AddAt(hMinv,0);
-  
-  // Degree of polynomial
-  fPolDeg = pod;
-    
-  // inclusion and exclusion areas (values)
-  const Double_t kJPMass    = 3.096916;
-  // inclusion and exclusion areas (bin numbers)
-  const Int_t binIntLo = hMinv->FindBin(fIntMin);
-  const Int_t binIntHi = hMinv->FindBin(fIntMax);
-  // for error calculation
-  Double_t intAreaEdgeLo = hMinv->GetBinLowEdge(binIntLo);
-  Double_t intAreaEdgeHi = hMinv->GetBinLowEdge(binIntHi)+hMinv->GetBinWidth(binIntHi);
-  Double_t norm = (binIntHi-binIntLo)/(fIntMax - fIntMin); 
-  
-  TH1F  *hBfit = (TH1F *) hMinv->Clone(); // for bg only fit (excluding peak region)
-  TH1F  *hSigF = (TH1F *) hMinv->Clone(); // signal with subtracted bg
-  TH1F  *hBgrd = (TH1F *) hMinv->Clone(); // bg histogram
-  
-  hBfit->Reset();
-  hSigF->Reset();
-  hBgrd->Reset();
-  
-  // extract start parameter for the MC signal fit
-  Double_t bgvalAv = (hMinv->Integral(1,hMinv->GetNbinsX()+1) - hMinv->Integral(binIntLo,binIntHi)) / (hMinv->GetNbinsX()+1 - (binIntHi-binIntLo));
-  Double_t pkval     = hMinv->GetBinContent(hMinv->FindBin(kJPMass)) - bgvalAv;
-  Double_t heightMC  = hSim->GetBinContent(hSim->FindBin(kJPMass));
-  Double_t peakScale = (heightMC > 0. ? pkval/heightMC : 0.0);
-  
-  Int_t nBgnd = 2 + fPolDeg; // degree + c1st oefficient + higher coefficients
-  Int_t nMinv = nBgnd + 1;  // bgrd + peakscale
-  
-  // Create the spectra without peak region
-  for (Int_t iBin = 0; iBin <= hMinv->GetNbinsX(); iBin++) {
-    if ((iBin < binIntLo) || (iBin > binIntHi)) {
-      hBfit->SetBinContent(iBin,hMinv->GetBinContent(iBin));
-      hBfit->SetBinError(iBin,hMinv->GetBinError(iBin));
-       }
-  }
-    
-  
-  // =======
-  //   1.
-  // =======
-  // Do the fit to the background spectrum
-  TF1 *fBo = new TF1("bgrd_fit",BgndFun,fFitMin,fFitMax,nBgnd);
-  for (Int_t iPar=0; iPar<nBgnd; iPar++) {
-    if (iPar == 0) fBo->FixParameter(0, fPolDeg);
-       if (iPar == 1) fBo->SetParameter(iPar, bgvalAv);
-    if (iPar >= 2) fBo->SetParameter(iPar, 0.);
-  }
-  hBfit->Fit(fBo,"0qR");
-  //hBfit->SetNameTitle("bgrd_fit");
-  arrResults->AddAt(fBo,1);
-  
-  
-  // =======
-  //   2.
-  // =======
-  // Fit the whole spectrum with peak and background
-  TF1 *fSB = new TF1("bgrd_peak_fit",MinvFun,fFitMin,fFitMax,nMinv);
-  fSB->FixParameter(0, fPolDeg);
-  fSB->SetParameter(1, peakScale);
-  // copy the polynomial parameters
-  for (Int_t iPar=0; iPar<=fPolDeg; iPar++) 
-       fSB->SetParameter(2+iPar, fBo->GetParameter(iPar+1));
-  hMinv->Fit(fSB,"0qR");
-  arrResults->AddAt(fSB,2);
-  
-  
-  // =======
-  //   3.
-  // =======
-  // Create the background function
-  TF1 *fB = new TF1("bgrdOnly_fkt",BgndFun,fFitMin,fFitMax,nBgnd);
-  fB->FixParameter(0,fPolDeg);
-  for (Int_t iDeg=0; iDeg<=fPolDeg; iDeg++) {
-       fB->SetParameter(1+iDeg,fSB->GetParameter(2+iDeg));
-       fB->SetParError(1+iDeg,fSB->GetParError(2+iDeg));
-  }
-  // create background histogram from background function
-  hBgrd->Eval(fB);
-  hBgrd->Fit(fB,"0qR");
-  // calculate the integral and integral error from fit function
-  Double_t intc = fB->Integral(intAreaEdgeLo, intAreaEdgeHi) * norm;
-  Double_t inte = fB->IntegralError(intAreaEdgeLo, intAreaEdgeHi) * norm;
-  arrResults->AddAt(fB,3);
-
-  // Fill the background spectrum erros. Use the error from the fit function for the background fB
-  for (Int_t iBin = 0; iBin <= hBgrd->GetNbinsX(); iBin++) {
-    Double_t x = hBgrd->GetBinCenter(iBin);
-    if ((x >= fFitMin) && (x <= fFitMax)) {
-         Double_t binte = inte / TMath::Sqrt((binIntHi-binIntLo)+1);
-         hBgrd->SetBinError(iBin,binte); 
-    }
-  }
-  arrResults->AddAt(hBgrd,4);
-  
-  // =======
-  //   4.
-  // =======
-  // Subtract the background  
-  hSigF->Add(hMinv,hBgrd,1.0,-1.0);
-  for (Int_t iBin = 0; iBin <= hSigF->GetNbinsX(); iBin++) {
-    Double_t x = hSigF->GetBinCenter(iBin);
-    if ((x < fFitMin) || (x > fFitMax)) {
-      hSigF->SetBinContent(iBin,0.0);
-      hSigF->SetBinError(iBin,0.0);
-    }
-  }
-  hSigF->SetNameTitle("peak_only","");
-  arrResults->AddAt(hSigF,5);
-  
-  // =======
-  //   5.
-  // =======
-  // Fit the background-subtracted spectrum
-  TF1 *fS  = new TF1("peakOnly_fit",PeakFunCB,fFitMin,fFitMax,5);
-  fS->SetParameters(-.05,1,kJPMass,.003,700);
-  fS->SetParNames("alpha","n","meanx","sigma","N");
-  hSigF->Fit(fS,"0qR");
-  arrResults->AddAt(fS,6);
-  
-  
-  // connect data members 
-  fFuncSignal     = (TF1*) arrResults->At(6)->Clone();
-  fFuncBackground = (TF1*) arrResults->At(3)->Clone();
-  fFuncSigBack    = (TF1*) arrResults->At(2)->Clone();
-  fHistSignal     = (TH1F*)arrResults->At(5)->Clone();
-  fHistBackground = (TH1F*)arrResults->At(4)->Clone();
-  
-  // fit results
-  Double_t chi2 = fSB->GetChisquare();
-  Int_t    ndf  = fSB->GetNDF();
-  fChi2Dof = chi2/ndf;
-  
-  // signal + signal error
-  fValues(0) = hSigF->IntegralAndError(binIntLo, binIntHi, fErrors(0));
-  fValues(1) = intc; // background
-  fErrors(1) = inte; // background error
-  // S/B (2) and significance (3)
-  SetSignificanceAndSOB();
-  fValues(4) = fS->GetParameter(2); // mass
-  fErrors(4) = fS->GetParError(2);  // mass error
-  fValues(5) = fS->GetParameter(3); // mass wdth
-  fErrors(5) = fS->GetParError(3);  // mass wdth error
-  
-    
-  delete arrResults;
-  
-}
-//______________________________________________________________________________
-Double_t AliDielectronSignalFunc::BgndFun(const Double_t *x, const Double_t *par) {
-  // parameters
-  // [0]:   degree of polynomial
-  // [1]:   constant polynomial coefficient
-  // [2]..: higher polynomial coefficients
-  
-  Int_t    deg   = ((Int_t) par[0]);
-  
-  Double_t f     = 0.0;
-  Double_t yy    = 1.0;
-  Double_t xx    = x[0];
-
-  for (Int_t i = 0; i <= deg; i++) {
-    f  += par[i+1] * yy;
-    yy *= xx;
-  } 
-  
-    
-  return f;
-}
 //______________________________________________________________________________
-Double_t AliDielectronSignalFunc::PeakFun(const Double_t *x, const Double_t *par) {
+Double_t AliDielectronSignalFunc::PeakFunMC(const Double_t *x, const Double_t *par) {
   // Fit MC signal shape
   // parameters
   // [0]:   scale for simpeak
@@ -360,7 +218,8 @@ Double_t AliDielectronSignalFunc::PeakFun(const Double_t *x, const Double_t *par
   
   TH1F *hPeak = fgHistSimPM;
   if (!hPeak) {
-    printf("F-AliDielectronSignalFunc::PeakFun: No histogram for peak fit defined!\n");
+    printf("E-AliDielectronSignalFunc::PeakFun: No histogram for peak fit defined!\n");
+    return 0.0;
   }
   
   Int_t idx = hPeak->FindBin(xx);
@@ -373,32 +232,9 @@ Double_t AliDielectronSignalFunc::PeakFun(const Double_t *x, const Double_t *par
   
 }
 
-//______________________________________________________________________________
-Double_t AliDielectronSignalFunc::MinvFun(const Double_t *x, const Double_t *par) {
-  // parameters
-  // [0]:   degree of polynomial             -> [0]   for BgndFun
-  // [1]:   scale for simpeak                -> [0]   for PeakFun
-  // [2]:   constant polynomial coefficient  -> [1]   for BgndFun
-  // [3]..: higher polynomial coefficients   -> [2].. for BgndFun
-  
-  Int_t    deg = ((Int_t) par[0]);
-  Double_t parPK[25], parBG[25];
-  
-  parBG[0] = par[0]; // degree of polynomial
-  
-  parPK[0] = par[1]; // MC minv scale
-  for (Int_t i = 0; i <= deg; i++) parBG[i+1] = par[i+2]; // polynomial coefficients
-  
-  Double_t peak = PeakFun(x,parPK);
-  Double_t bgnd = BgndFun(x,parBG);
-  Double_t f    = peak + bgnd;
-  
-  return f;
-}
-
 //______________________________________________________________________________
 Double_t AliDielectronSignalFunc::PeakFunCB(const Double_t *x, const Double_t *par) {
-  // Crystal Ball function fit
+  // Crystal Ball fit function
   
   Double_t alpha = par[0];
   Double_t     n = par[1];
@@ -421,6 +257,17 @@ Double_t AliDielectronSignalFunc::PeakFunCB(const Double_t *x, const Double_t *p
   return fitval;
 }
 
+//______________________________________________________________________________
+Double_t AliDielectronSignalFunc::PeakFunGaus(const Double_t *x, const Double_t *par) {
+  // Gaussian fit function
+  //printf("fNparBgrd %d \n",fNparBgnd);
+  Double_t     n = par[0];
+  Double_t  mean = par[1];
+  Double_t sigma = par[2];
+  Double_t    xx = x[0];
+
+  return ( n*TMath::Exp(-0.5*TMath::Power((xx-mean)/sigma,2)) );
+}
 
 //______________________________________________
 void AliDielectronSignalFunc::ProcessFit(TObjArray * const arrhist) {
@@ -447,34 +294,43 @@ void AliDielectronSignalFunc::ProcessFit(TObjArray * const arrhist) {
   fHistDataPM->Fit(fFuncSigBack, fFitOpt.Data(), "", fFitMin, fFitMax);
   TFitResultPtr pmFitPtr = fHistDataPM->Fit(fFuncSigBack, fFitOpt.Data(), "", fFitMin, fFitMax);
   //TFitResult *pmFitResult = pmFitPtr.Get(); // used only with TF1Helper
+  //fFuncBackground->SetParameters(fFuncSigBack->GetParameters());
   fFuncSignal->SetParameters(fFuncSigBack->GetParameters());
   fFuncBackground->SetParameters(fFuncSigBack->GetParameters()+fFuncSignal->GetNpar());
-  
+
+  // fill the background spectrum
+  fHistBackground->Eval(fFuncBackground);
+  // set the error for the background histogram
+  fHistBackground->Fit(fFuncBackground,"0qR","",fFitMin,fFitMax);
+  Double_t inte  = fFuncBackground->IntegralError(fIntMin, fIntMax)/fHistDataPM->GetBinWidth(1);
+  Double_t binte = inte / TMath::Sqrt((fHistDataPM->FindBin(fIntMax)-fHistDataPM->FindBin(fIntMin))+1);
+  for(Int_t iBin=fHistDataPM->FindBin(fIntMin); iBin<=fHistDataPM->FindBin(fIntMax); iBin++) {
+    fHistBackground->SetBinError(iBin, binte);
+  }
+
   for(Int_t iBin=1; iBin<=fHistDataPM->GetXaxis()->GetNbins(); iBin++) {
-    Double_t m = fHistDataPM->GetBinCenter(iBin);
+    //    Double_t m = fHistDataPM->GetBinCenter(iBin);
     Double_t pm = fHistDataPM->GetBinContent(iBin);
     Double_t epm = fHistDataPM->GetBinError(iBin);
-    Double_t bknd = fFuncBackground->Eval(m);
-    Double_t ebknd = 0;
+    Double_t bknd = fHistBackground->GetBinContent(iBin);
+    Double_t ebknd = fHistBackground->GetBinError(iBin);
     for(Int_t iPar=fFuncSignal->GetNpar(); iPar<fFuncSigBack->GetNpar(); iPar++) {
-/* TF1Helper problem on alien compilation
-      for(Int_t jPar=iPar; jPar<fFuncSigBack->GetNpar(); jPar++) {
-        TF1 gradientIpar("gradientIpar",
-                         ROOT::TF1Helper::TGradientParFunction(iPar-fFuncSignal->GetNpar(),fFuncBackground),0,0,0);
-        TF1 gradientJpar("gradientJpar",
-                         ROOT::TF1Helper::TGradientParFunction(jPar-fFuncSignal->GetNpar(),fFuncBackground),0,0,0);
-        ebknd += pmFitResult->CovMatrix(iPar,jPar)*
-          gradientIpar.Eval(m)*gradientJpar.Eval(m)*
-          (iPar==jPar ? 1.0 : 2.0);
-      }
-*/
+      /* TF1Helper problem on alien compilation
+        for(Int_t jPar=iPar; jPar<fFuncSigBack->GetNpar(); jPar++) {
+        TF1 gradientIpar("gradientIpar",
+        ROOT::TF1Helper::TGradientParFunction(iPar-fFuncSignal->GetNpar(),fFuncBackground),0,0,0);
+        TF1 gradientJpar("gradientJpar",
+        ROOT::TF1Helper::TGradientParFunction(jPar-fFuncSignal->GetNpar(),fFuncBackground),0,0,0);
+        ebknd += pmFitResult->CovMatrix(iPar,jPar)*
+        gradientIpar.Eval(m)*gradientJpar.Eval(m)*
+        (iPar==jPar ? 1.0 : 2.0);
+        }
+      */
     }
     Double_t signal = pm-bknd;
     Double_t error = TMath::Sqrt(epm*epm+ebknd);
     fHistSignal->SetBinContent(iBin, signal);
     fHistSignal->SetBinError(iBin, error);
-    fHistBackground->SetBinContent(iBin, bknd);
-    fHistBackground->SetBinError(iBin, TMath::Sqrt(ebknd));
   }
   
   if(fUseIntegral) {
@@ -496,7 +352,7 @@ void AliDielectronSignalFunc::ProcessFit(TObjArray * const arrhist) {
     }
     // background
     fValues(1) = fFuncBackground->Integral(fIntMin, fIntMax)/fHistDataPM->GetBinWidth(1);
-    fErrors(1) = 0;
+    fErrors(1) = fFuncBackground->IntegralError(fIntMin, fIntMax)/fHistDataPM->GetBinWidth(1);
     for(Int_t iPar=fFuncSignal->GetNpar(); iPar<fFuncSigBack->GetNpar(); iPar++) {
 /* TF1Helper problem on alien compilation
       for(Int_t jPar=iPar; jPar<fFuncSigBack->GetNpar(); jPar++) {
@@ -517,8 +373,8 @@ void AliDielectronSignalFunc::ProcessFit(TObjArray * const arrhist) {
                                                fHistSignal->FindBin(fIntMax), fErrors(0));
     // background
     fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
-      fHistBackground->FindBin(fIntMax),
-      fErrors(1));
+                                                  fHistBackground->FindBin(fIntMax),
+                                                  fErrors(1));
   }
   // S/B and significance
   SetSignificanceAndSOB();
@@ -528,10 +384,13 @@ void AliDielectronSignalFunc::ProcessFit(TObjArray * const arrhist) {
   fErrors(5) = fFuncSigBack->GetParError(fParMassWidth);
   
   fProcessed = kTRUE;
+  
+  fHistBackground->GetListOfFunctions()->Add(fFuncBackground);
+
 }
 
 //______________________________________________
-void AliDielectronSignalFunc::ProcessLS(TObjArray * const arrhist) {
+void AliDielectronSignalFunc::ProcessFitLS(TObjArray * const arrhist) {
   //
   // Substract background using the like-sign spectrum
   //
@@ -636,7 +495,7 @@ void AliDielectronSignalFunc::ProcessLS(TObjArray * const arrhist) {
 }
 
 //______________________________________________
-void AliDielectronSignalFunc::ProcessEM(TObjArray * const arrhist) {
+void AliDielectronSignalFunc::ProcessFitEM(TObjArray * const arrhist) {
   //
   // Substract background with the event mixing technique
   //
@@ -663,6 +522,7 @@ void AliDielectronSignalFunc::SetFunctions(TF1 * const combined, TF1 * const sig
   fFuncSigBack=combined;
   fParMass=parM;
   fParMassWidth=parMres;
+
 }
 
 //______________________________________________
@@ -776,3 +636,33 @@ void AliDielectronSignalFunc::Draw(const Option_t* option)
   if (optStat) DrawStats();
   
 }
+
+
+//______________________________________________________________________________
+void AliDielectronSignalFunc::CombineFunc(TF1 * const peak, TF1 * const bgnd) {
+  //
+  // combine the bgnd and the peak function
+  //
+
+  if (!peak||!bgnd) {
+    AliError("Both, signal and background function need to be set!");
+    return;
+  }
+  fFuncSignal=peak;
+  fFuncBackground=bgnd;
+
+  fNparPeak     = fFuncSignal->GetNpar();
+  fNparBgnd     = fFuncBackground->GetNpar();
+
+  fFuncSigBack = new TF1("BgndPeak",AliDielectronSignalFunc::PeakBgndFun, fFitMin,fFitMax, fNparPeak+fNparBgnd);
+  return;
+}
+
+//______________________________________________________________________________
+Double_t AliDielectronSignalFunc::PeakBgndFun(const Double_t *x, const Double_t *par) {
+  //
+  // merge peak and bgnd functions
+  //
+  return (fFuncSignal->EvalPar(x,par) + fFuncBackground->EvalPar(x,par+fNparPeak));
+}
+