[u/mrichter/AliRoot.git] / PWG3 / dielectron / AliDielectronSignalFunc.cxx

/*************************************************************************
* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
*                                                                        *
* Author: The ALICE Off-line Project.                                    *
* Contributors are mentioned in the code where appropriate.              *
*                                                                        *
* Permission to use, copy, modify and distribute this software and its   *
* documentation strictly for non-commercial purposes is hereby granted   *
* without fee, provided that the above copyright notice appears in all   *
* copies and that both the copyright notice and this permission notice   *
* appear in the supporting documentation. The authors make no claims     *
* about the suitability of this software for any purpose. It is          *
* provided "as is" without express or implied warranty.                  *
**************************************************************************/

///////////////////////////////////////////////////////////////////////////
//                Dielectron SignalFunc                                  //
//                                                                       //
//                                                                       //
/*
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.

*/
//                                                                       //
///////////////////////////////////////////////////////////////////////////

#include <TF1.h>
#include <TH1.h>
#include <TGraph.h>
#include <TMath.h>
#include <TString.h>
#include <TPaveText.h>
#include <TList.h>
#include <TFitResult.h>
//#include <../hist/hist/src/TF1Helper.h>

#include <AliLog.h>

#include "AliDielectronSignalFunc.h"

ClassImp(AliDielectronSignalFunc)

AliDielectronSignalFunc::AliDielectronSignalFunc() :
AliDielectronSignalBase(),
fFuncSignal(0x0),
fFuncBackground(0x0),
fFuncSigBack(0x0),
fParMass(1),
fParMassWidth(2),
fFitOpt("SMNQE"),
fUseIntegral(kFALSE)
{
  //
  // Default Constructor
  //
  
}

//______________________________________________
AliDielectronSignalFunc::AliDielectronSignalFunc(const char* name, const char* title) :
AliDielectronSignalBase(name, title),
fFuncSignal(0x0),
fFuncBackground(0x0),
fFuncSigBack(0x0),
fParMass(1),
fParMassWidth(2),
fFitOpt("SMNQE"),
fUseIntegral(kFALSE)
{
  //
  // Named Constructor
  //
}

//______________________________________________
AliDielectronSignalFunc::~AliDielectronSignalFunc()
{
  //
  // Default Destructor
  //
  if(fFuncSignal) delete fFuncSignal;
  if(fFuncBackground) delete fFuncBackground;
  if(fFuncSigBack) delete fFuncSigBack;
}


//______________________________________________
void AliDielectronSignalFunc::Process(TObjArray * const arrhist)
{
  //
  // Fit the invariant mass histograms and retrieve the signal and background
  //
  switch(fMethod) {
  case kFitted :
    ProcessFit(arrhist);
    break;
    
  case kLikeSign :
    ProcessLS(arrhist);
    break;
    
  case kEventMixing :
    ProcessEM(arrhist);
    break;
    
  default :
    AliError("Background substraction method not known!");
  }
}


//______________________________________________
void AliDielectronSignalFunc::ProcessFit(TObjArray * const arrhist) {
  //
  // Fit the +- invariant mass distribution only
  // Here we assume that the combined fit function is a sum of the signal and background functions
  //    and that the signal function is always the first term of this sum
  //
  
  fHistDataPM = (TH1F*)(arrhist->At(1))->Clone("histPM");  // +-    SE
  fHistDataPM->Sumw2();
  if(fRebin>1)
    fHistDataPM->Rebin(fRebin);
  
  fHistSignal = new TH1F("HistSignal", "Like-Sign substracted signal",
                         fHistDataPM->GetXaxis()->GetNbins(),
                         fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
  fHistBackground = new TH1F("HistBackground", "Like-sign contribution",
                             fHistDataPM->GetXaxis()->GetNbins(),
                             fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
  
  // the starting parameters of the fit function and their limits can be tuned
  // by the user in its macro
  fHistDataPM->Fit(fFuncSigBack, fFitOpt.Data(), "", fFitMin, fFitMax);
  TFitResultPtr pmFitPtr = fHistDataPM->Fit(fFuncSigBack, fFitOpt.Data(), "", fFitMin, fFitMax);
  //TFitResult *pmFitResult = pmFitPtr.Get(); // used only with TF1Helper
  fFuncSignal->SetParameters(fFuncSigBack->GetParameters());
  fFuncBackground->SetParameters(fFuncSigBack->GetParameters()+fFuncSignal->GetNpar());
  
  for(Int_t iBin=1; iBin<=fHistDataPM->GetXaxis()->GetNbins(); 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;
    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);
      }
*/
    }
    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) {
    // signal
    fValues(0) = fFuncSignal->Integral(fIntMin, fIntMax)/fHistDataPM->GetBinWidth(1);
    fErrors(0) = 0;
    for(Int_t iPar=0; iPar<fFuncSignal->GetNpar(); iPar++) {
/* TF1Helper problem on alien compilation
      for(Int_t jPar=iPar; jPar<fFuncSignal->GetNpar(); jPar++) {
        TF1 gradientIpar("gradientIpar",
                         ROOT::TF1Helper::TGradientParFunction(iPar,fFuncSignal),0,0,0);
        TF1 gradientJpar("gradientJpar",
                         ROOT::TF1Helper::TGradientParFunction(jPar,fFuncSignal),0,0,0);
        fErrors(0) += pmFitResult->CovMatrix(iPar,jPar)*
          gradientIpar.Integral(fIntMin,fIntMax)*gradientJpar.Integral(fIntMin,fIntMax)*
          (iPar==jPar ? 1.0 : 2.0);
      }
*/
    }
    // background
    fValues(1) = fFuncBackground->Integral(fIntMin, fIntMax)/fHistDataPM->GetBinWidth(1);
    fErrors(1) = 0;
    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);
        fErrors(1) += pmFitResult->CovMatrix(iPar,jPar)*
          gradientIpar.Integral(fIntMin, fIntMax)*gradientJpar.Integral(fIntMin, fIntMax)*
          (iPar==jPar ? 1.0 : 2.0);
      }
*/
    }
  }
  else {
    // signal
    fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
                                               fHistSignal->FindBin(fIntMax), fErrors(0));
    // background
    fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
      fHistBackground->FindBin(fIntMax),
      fErrors(1));
  }
  // S/B and significance
  SetSignificanceAndSOB();
  fValues(4) = fFuncSigBack->GetParameter(fParMass);
  fErrors(4) = fFuncSigBack->GetParError(fParMass);
  fValues(5) = fFuncSigBack->GetParameter(fParMassWidth);
  fErrors(5) = fFuncSigBack->GetParError(fParMassWidth);
  
  fProcessed = kTRUE;
}

//______________________________________________
void AliDielectronSignalFunc::ProcessLS(TObjArray * const arrhist) {
  //
  // Substract background using the like-sign spectrum
  //
  fHistDataPP = (TH1F*)(arrhist->At(0))->Clone("histPP");  // ++    SE
  fHistDataPM = (TH1F*)(arrhist->At(1))->Clone("histPM");  // +-    SE
  fHistDataMM = (TH1F*)(arrhist->At(2))->Clone("histMM");  // --    SE
  if (fRebin>1) {
    fHistDataPP->Rebin(fRebin);
    fHistDataPM->Rebin(fRebin);
    fHistDataMM->Rebin(fRebin);
  }
  fHistDataPP->Sumw2();
  fHistDataPM->Sumw2();
  fHistDataMM->Sumw2();
  
  fHistSignal = new TH1F("HistSignal", "Like-Sign substracted signal",
                         fHistDataPM->GetXaxis()->GetNbins(),
                         fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
  fHistBackground = new TH1F("HistBackground", "Like-sign contribution",
                             fHistDataPM->GetXaxis()->GetNbins(),
                             fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
  
  // fit the +- mass distribution
  fHistDataPM->Fit(fFuncSigBack, fFitOpt.Data(), "", fFitMin, fFitMax);
  fHistDataPM->Fit(fFuncSigBack, fFitOpt.Data(), "", fFitMin, fFitMax);
  // declare the variables where the like-sign fit results will be stored
  TFitResult *ppFitResult = 0x0;
  TFitResult *mmFitResult = 0x0;
  // fit the like sign background
  TF1 *funcClonePP = (TF1*)fFuncBackground->Clone("funcClonePP");
  TF1 *funcCloneMM = (TF1*)fFuncBackground->Clone("funcCloneMM");
  fHistDataPP->Fit(funcClonePP, fFitOpt.Data(), "", fFitMin, fFitMax);
  TFitResultPtr ppFitPtr = fHistDataPP->Fit(funcClonePP, fFitOpt.Data(), "", fFitMin, fFitMax);
  ppFitResult = ppFitPtr.Get();
  fHistDataMM->Fit(funcCloneMM, fFitOpt.Data(), "", fFitMin, fFitMax);
  TFitResultPtr mmFitPtr = fHistDataMM->Fit(funcCloneMM, fFitOpt.Data(), "", fFitMin, fFitMax);
  mmFitResult = mmFitPtr.Get();
  
  for(Int_t iBin=1; iBin<=fHistDataPM->GetXaxis()->GetNbins(); iBin++) {
    Double_t m = fHistDataPM->GetBinCenter(iBin);
    Double_t pm = fHistDataPM->GetBinContent(iBin);
    Double_t pp = funcClonePP->Eval(m);
    Double_t mm = funcCloneMM->Eval(m);
    Double_t epm = fHistDataPM->GetBinError(iBin);
    Double_t epp = 0;
    for(Int_t iPar=0; iPar<funcClonePP->GetNpar(); iPar++) {
/* TF1Helper problem on alien compilation
      for(Int_t jPar=iPar; jPar<funcClonePP->GetNpar(); jPar++) {
        TF1 gradientIpar("gradientIpar",
                         ROOT::TF1Helper::TGradientParFunction(iPar,funcClonePP),0,0,0);
        TF1 gradientJpar("gradientJpar",
                         ROOT::TF1Helper::TGradientParFunction(jPar,funcClonePP),0,0,0);
        epp += ppFitResult->CovMatrix(iPar,jPar)*
          gradientIpar.Eval(m)*gradientJpar.Eval(m)*
          (iPar==jPar ? 1.0 : 2.0);
      }
*/
    }
    Double_t emm = 0;
    for(Int_t iPar=0; iPar<funcCloneMM->GetNpar(); iPar++) {
/* TF1Helper problem on alien compilation
      for(Int_t jPar=iPar; jPar<funcCloneMM->GetNpar(); jPar++) {
        TF1 gradientIpar("gradientIpar",
                         ROOT::TF1Helper::TGradientParFunction(iPar,funcCloneMM),0,0,0);
        TF1 gradientJpar("gradientJpar",
                         ROOT::TF1Helper::TGradientParFunction(jPar,funcCloneMM),0,0,0);
        emm += mmFitResult->CovMatrix(iPar,jPar)*
          gradientIpar.Eval(m)*gradientJpar.Eval(m)*
          (iPar==jPar ? 1.0 : 2.0);
      }
*/
    }
    
    Double_t signal = pm-2.0*TMath::Sqrt(pp*mm);
    Double_t background = 2.0*TMath::Sqrt(pp*mm);
    // error propagation on the signal calculation above
    Double_t esignal = TMath::Sqrt(epm*epm+(mm/pp)*epp+(pp/mm)*emm);
    Double_t ebackground = TMath::Sqrt((mm/pp)*epp+(pp/mm)*emm);
    fHistSignal->SetBinContent(iBin, signal);
    fHistSignal->SetBinError(iBin, esignal);
    fHistBackground->SetBinContent(iBin, background);
    fHistBackground->SetBinError(iBin, ebackground);
  }
  
  // signal
  fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
                                             fHistSignal->FindBin(fIntMax), fErrors(0));
  // background
  fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
                                                 fHistBackground->FindBin(fIntMax),
                                                 fErrors(1));
  // S/B and significance
  SetSignificanceAndSOB();
  fValues(4) = fFuncSigBack->GetParameter(fParMass);
  fErrors(4) = fFuncSigBack->GetParError(fParMass);
  fValues(5) = fFuncSigBack->GetParameter(fParMassWidth);
  fErrors(5) = fFuncSigBack->GetParError(fParMassWidth);
  
  fProcessed = kTRUE;
}

//______________________________________________
void AliDielectronSignalFunc::ProcessEM(TObjArray * const arrhist) {
  //
  // Substract background with the event mixing technique
  //
  arrhist->GetEntries();   // just to avoid the unused parameter warning
  AliError("Event mixing for background substraction method not implemented!");
}

//______________________________________________
void AliDielectronSignalFunc::SetFunctions(TF1 * const combined, TF1 * const sig, TF1 * const back,
                                           Int_t parM, Int_t parMres)
{
  //
  // Set the signal, background functions and combined fit function
  // Note: The process method assumes that the first n parameters in the
  //       combined fit function correspond to the n parameters of the signal function
  //       and the n+1 to n+m parameters to the m parameters of the background function!!!
  
  if (!sig||!back||!combined) {
    AliError("Both, signal and background function need to be set!");
    return;
  }
  fFuncSignal=sig;
  fFuncBackground=back;
  fFuncSigBack=combined;
  fParMass=parM;
  fParMassWidth=parMres;
}

//______________________________________________
void AliDielectronSignalFunc::SetDefaults(Int_t type)
{
  //
  // Setup some default functions:
  // type = 0: gaus signal + linear background in 2.5 - 4 GeV inv. mass
  // type = 1: gaus signal + exponential background in 2.5 - 4 GeV inv. mass
  // type = 2: half gaussian, half exponential signal function
  // type = 3: Crystal-Ball function
  // type = 4: Crystal-Ball signal + exponential background
  //
  
  if (type==0){
    fFuncSignal=new TF1("DieleSignal","gaus",2.5,4);
    fFuncBackground=new TF1("DieleBackground","pol1",2.5,4);
    fFuncSigBack=new TF1("DieleCombined","gaus+pol1(3)",2.5,4);
    
    fFuncSigBack->SetParameters(1,3.1,.05,2.5,1);
    fFuncSigBack->SetParLimits(0,0,10000000);
    fFuncSigBack->SetParLimits(1,3.05,3.15);
    fFuncSigBack->SetParLimits(2,.02,.1);
  }
  else if (type==1){
    fFuncSignal=new TF1("DieleSignal","gaus",2.5,4);
    fFuncBackground=new TF1("DieleBackground","[0]*exp(-(x-[1])/[2])",2.5,4);
    fFuncSigBack=new TF1("DieleCombined","gaus+[3]*exp(-(x-[4])/[5])",2.5,4);
    
    fFuncSigBack->SetParameters(1,3.1,.05,1,2.5,1);
    fFuncSigBack->SetParLimits(0,0,10000000);
    fFuncSigBack->SetParLimits(1,3.05,3.15);
    fFuncSigBack->SetParLimits(2,.02,.1);
  }
  else if (type==2){
    // half gaussian, half exponential signal function
    // exponential background
    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);
    fFuncBackground = new TF1("DieleBackground","[0]*exp(-(x-[1])/[2])+[3]",2.5,4);
    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);
    fFuncSigBack->SetParameters(1.,3.1,.05,.1,1,2.5,1,0);
    
    fFuncSigBack->SetParLimits(0,0,10000000);
    fFuncSigBack->SetParLimits(1,3.05,3.15);
    fFuncSigBack->SetParLimits(2,.02,.1);
    fFuncSigBack->FixParameter(6,2.5);
    fFuncSigBack->FixParameter(7,0);
  }
}


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