[u/mrichter/AliRoot.git] / PWGDQ / 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);
  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);
  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
45b2b1b8	252	// TFitResult *ppFitResult = 0x0;
45b2b1b8	253	// TFitResult *mmFitResult = 0x0;
bc75eeb5	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);
45b2b1b8	258	fHistDataPP->Fit(funcClonePP, fFitOpt.Data(), "", fFitMin, fFitMax);
	259	// TFitResultPtr ppFitPtr = fHistDataPP->Fit(funcClonePP, fFitOpt.Data(), "", fFitMin, fFitMax);
	260	// ppFitResult = ppFitPtr.Get();
	261	fHistDataMM->Fit(funcCloneMM, fFitOpt.Data(), "", fFitMin, fFitMax);
bc75eeb5	262	fHistDataMM->Fit(funcCloneMM, fFitOpt.Data(), "", fFitMin, fFitMax);
45b2b1b8	263	// TFitResultPtr mmFitPtr = fHistDataMM->Fit(funcCloneMM, fFitOpt.Data(), "", fFitMin, fFitMax);
45b2b1b8	264	// mmFitResult = mmFitPtr.Get();
3505bfad	265
bc75eeb5	266	for(Int_t iBin=1; iBin<=fHistDataPM->GetXaxis()->GetNbins(); iBin++) {
	267	Double_t m = fHistDataPM->GetBinCenter(iBin);
	268	Double_t pm = fHistDataPM->GetBinContent(iBin);
	269	Double_t pp = funcClonePP->Eval(m);
	270	Double_t mm = funcCloneMM->Eval(m);
	271	Double_t epm = fHistDataPM->GetBinError(iBin);
	272	Double_t epp = 0;
	273	for(Int_t iPar=0; iPar<funcClonePP->GetNpar(); iPar++) {
cafdf10b	274	/* TF1Helper problem on alien compilation
bc75eeb5	275	for(Int_t jPar=iPar; jPar<funcClonePP->GetNpar(); jPar++) {
bc75eeb5	276	TF1 gradientIpar("gradientIpar",
3505bfad	277	ROOT::TF1Helper::TGradientParFunction(iPar,funcClonePP),0,0,0);
bc75eeb5	278	TF1 gradientJpar("gradientJpar",
3505bfad	279	ROOT::TF1Helper::TGradientParFunction(jPar,funcClonePP),0,0,0);
	280	epp += ppFitResult->CovMatrix(iPar,jPar)*
	281	gradientIpar.Eval(m)gradientJpar.Eval(m)
	282	(iPar==jPar ? 1.0 : 2.0);
bc75eeb5	283	}
cafdf10b	284	*/
bc75eeb5	285	}
	286	Double_t emm = 0;
	287	for(Int_t iPar=0; iPar<funcCloneMM->GetNpar(); iPar++) {
cafdf10b	288	/* TF1Helper problem on alien compilation
bc75eeb5	289	for(Int_t jPar=iPar; jPar<funcCloneMM->GetNpar(); jPar++) {
bc75eeb5	290	TF1 gradientIpar("gradientIpar",
3505bfad	291	ROOT::TF1Helper::TGradientParFunction(iPar,funcCloneMM),0,0,0);
bc75eeb5	292	TF1 gradientJpar("gradientJpar",
3505bfad	293	ROOT::TF1Helper::TGradientParFunction(jPar,funcCloneMM),0,0,0);
	294	emm += mmFitResult->CovMatrix(iPar,jPar)*
	295	gradientIpar.Eval(m)gradientJpar.Eval(m)
	296	(iPar==jPar ? 1.0 : 2.0);
bc75eeb5	297	}
cafdf10b	298	*/
bc75eeb5	299	}
3505bfad	300
bc75eeb5	301	Double_t signal = pm-2.0TMath::Sqrt(ppmm);
	302	Double_t background = 2.0TMath::Sqrt(ppmm);
	303	// error propagation on the signal calculation above
	304	Double_t esignal = TMath::Sqrt(epmepm+(mm/pp)epp+(pp/mm)*emm);
	305	Double_t ebackground = TMath::Sqrt((mm/pp)epp+(pp/mm)emm);
	306	fHistSignal->SetBinContent(iBin, signal);
	307	fHistSignal->SetBinError(iBin, esignal);
	308	fHistBackground->SetBinContent(iBin, background);
	309	fHistBackground->SetBinError(iBin, ebackground);
	310	}
2a14a7b1	311
bc75eeb5	312	// signal
bc75eeb5	313	fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
3505bfad	314	fHistSignal->FindBin(fIntMax), fErrors(0));
bc75eeb5	315	// background
bc75eeb5	316	fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
3505bfad	317	fHistBackground->FindBin(fIntMax),
3505bfad	318	fErrors(1));
bc75eeb5	319	// S/B and significance
	320	SetSignificanceAndSOB();
	321	fValues(4) = fFuncSigBack->GetParameter(fParMass);
	322	fErrors(4) = fFuncSigBack->GetParError(fParMass);
	323	fValues(5) = fFuncSigBack->GetParameter(fParMassWidth);
	324	fErrors(5) = fFuncSigBack->GetParError(fParMassWidth);
3505bfad	325
bc75eeb5	326	fProcessed = kTRUE;
572b0139	327	}
bc75eeb5	328
bc75eeb5	329	//______________________________________________
3505bfad	330	void AliDielectronSignalFunc::ProcessEM(TObjArray * const arrhist) {
bc75eeb5	331	//
	332	// Substract background with the event mixing technique
	333	//
3505bfad	334	arrhist->GetEntries(); // just to avoid the unused parameter warning
	335	AliError("Event mixing for background substraction method not implemented!");
	336	}
bc75eeb5	337
572b0139	338	//______________________________________________
bc75eeb5	339	void AliDielectronSignalFunc::SetFunctions(TF1 * const combined, TF1 * const sig, TF1 * const back,
572b0139	340	Int_t parM, Int_t parMres)
	341	{
	342	//
	343	// Set the signal, background functions and combined fit function
	344	// Note: The process method assumes that the first n parameters in the
	345	// combined fit function correspond to the n parameters of the signal function
	346	// and the n+1 to n+m parameters to the m parameters of the background function!!!
bc75eeb5	347
572b0139	348	if (!sig\|\|!back\|\|!combined) {
572b0139	349	AliError("Both, signal and background function need to be set!");
61d106d3	350	return;
572b0139	351	}
bc75eeb5	352	fFuncSignal=sig;
	353	fFuncBackground=back;
	354	fFuncSigBack=combined;
	355	fParMass=parM;
	356	fParMassWidth=parMres;
572b0139	357	}
	358
	359	//______________________________________________
	360	void AliDielectronSignalFunc::SetDefaults(Int_t type)
	361	{
	362	//
	363	// Setup some default functions:
	364	// type = 0: gaus signal + linear background in 2.5 - 4 GeV inv. mass
	365	// type = 1: gaus signal + exponential background in 2.5 - 4 GeV inv. mass
	366	// type = 2: half gaussian, half exponential signal function
	367	// type = 3: Crystal-Ball function
	368	// type = 4: Crystal-Ball signal + exponential background
	369	//
572b0139	370
572b0139	371	if (type==0){
bc75eeb5	372	fFuncSignal=new TF1("DieleSignal","gaus",2.5,4);
	373	fFuncBackground=new TF1("DieleBackground","pol1",2.5,4);
	374	fFuncSigBack=new TF1("DieleCombined","gaus+pol1(3)",2.5,4);
572b0139	375
bc75eeb5	376	fFuncSigBack->SetParameters(1,3.1,.05,2.5,1);
	377	fFuncSigBack->SetParLimits(0,0,10000000);
	378	fFuncSigBack->SetParLimits(1,3.05,3.15);
	379	fFuncSigBack->SetParLimits(2,.02,.1);
3505bfad	380	}
bc75eeb5	381	else if (type==1){
	382	fFuncSignal=new TF1("DieleSignal","gaus",2.5,4);
	383	fFuncBackground=new TF1("DieleBackground","[0]*exp(-(x-[1])/[2])",2.5,4);
	384	fFuncSigBack=new TF1("DieleCombined","gaus+[3]*exp(-(x-[4])/[5])",2.5,4);
572b0139	385
bc75eeb5	386	fFuncSigBack->SetParameters(1,3.1,.05,1,2.5,1);
	387	fFuncSigBack->SetParLimits(0,0,10000000);
	388	fFuncSigBack->SetParLimits(1,3.05,3.15);
	389	fFuncSigBack->SetParLimits(2,.02,.1);
3505bfad	390	}
bc75eeb5	391	else if (type==2){
572b0139	392	// half gaussian, half exponential signal function
572b0139	393	// exponential background
bc75eeb5	394	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);
	395	fFuncBackground = new TF1("DieleBackground","[0]*exp(-(x-[1])/[2])+[3]",2.5,4);
	396	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);
	397	fFuncSigBack->SetParameters(1.,3.1,.05,.1,1,2.5,1,0);
572b0139	398
bc75eeb5	399	fFuncSigBack->SetParLimits(0,0,10000000);
	400	fFuncSigBack->SetParLimits(1,3.05,3.15);
	401	fFuncSigBack->SetParLimits(2,.02,.1);
	402	fFuncSigBack->FixParameter(6,2.5);
	403	fFuncSigBack->FixParameter(7,0);
572b0139	404	}
	405	}
	406
	407
	408	//______________________________________________
	409	void AliDielectronSignalFunc::Draw(const Option_t* option)
	410	{
	411	//
	412	// Draw the fitted function
	413	//
bc75eeb5	414
572b0139	415	TString drawOpt(option);
572b0139	416	drawOpt.ToLower();
3505bfad	417
572b0139	418	Bool_t optStat=drawOpt.Contains("stat");
572b0139	419
bc75eeb5	420	fFuncSigBack->SetNpx(200);
	421	fFuncSigBack->SetRange(fIntMin,fIntMax);
	422	fFuncBackground->SetNpx(200);
	423	fFuncBackground->SetRange(fIntMin,fIntMax);
572b0139	424
bc75eeb5	425	TGraph *grSig=new TGraph(fFuncSigBack);
572b0139	426	grSig->SetFillColor(kGreen);
572b0139	427	grSig->SetFillStyle(3001);
3505bfad	428
bc75eeb5	429	TGraph *grBack=new TGraph(fFuncBackground);
572b0139	430	grBack->SetFillColor(kRed);
572b0139	431	grBack->SetFillStyle(3001);
3505bfad	432
572b0139	433	grSig->SetPoint(0,grBack->GetX()[0],grBack->GetY()[0]);
	434	grSig->SetPoint(grSig->GetN()-1,grBack->GetX()[grBack->GetN()-1],grBack->GetY()[grBack->GetN()-1]);
	435
	436	grBack->SetPoint(0,grBack->GetX()[0],0.);
	437	grBack->SetPoint(grBack->GetN()-1,grBack->GetX()[grBack->GetN()-1],0.);
	438
bc75eeb5	439	fFuncSigBack->SetRange(fFitMin,fFitMax);
bc75eeb5	440	fFuncBackground->SetRange(fFitMin,fFitMax);
572b0139	441
572b0139	442	if (!drawOpt.Contains("same")){
bc75eeb5	443	if (fHistDataPM){
bc75eeb5	444	fHistDataPM->Draw();
8df8e382	445	grSig->Draw("f");
	446	} else {
	447	grSig->Draw("af");
	448	}
572b0139	449	} else {
	450	grSig->Draw("f");
	451	}
bc75eeb5	452	if(fMethod==kFitted) grBack->Draw("f");
	453	fFuncSigBack->Draw("same");
	454	fFuncSigBack->SetLineWidth(2);
	455	if(fMethod==kLikeSign) {
	456	fHistDataPP->SetLineWidth(2);
	457	fHistDataPP->SetLineColor(6);
	458	fHistDataPP->Draw("same");
	459	fHistDataMM->SetLineWidth(2);
3505bfad	460	fHistDataMM->SetLineColor(8);
bc75eeb5	461	fHistDataMM->Draw("same");
bc75eeb5	462	}
3505bfad	463
3505bfad	464	if(fMethod==kFitted)
bc75eeb5	465	fFuncBackground->Draw("same");
3505bfad	466
8df8e382	467	if (optStat) DrawStats();
bc75eeb5	468
572b0139	469	}