1 /*************************************************************************
2 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 ///////////////////////////////////////////////////////////////////////////
17 // Dielectron SignalFunc //
21 Dielectron signal extraction class using functions as input.
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
29 ///////////////////////////////////////////////////////////////////////////
36 #include <TPaveText.h>
40 #include "AliDielectronSignalFunc.h"
42 ClassImp(AliDielectronSignalFunc)
44 AliDielectronSignalFunc::AliDielectronSignalFunc() :
45 AliDielectronSignalBase(),
59 // Default Constructor
64 //______________________________________________
65 AliDielectronSignalFunc::AliDielectronSignalFunc(const char* name, const char* title) :
66 AliDielectronSignalBase(name, title),
84 //______________________________________________
85 AliDielectronSignalFunc::~AliDielectronSignalFunc()
90 if (fSigBack) delete fSigBack;
94 //______________________________________________
95 void AliDielectronSignalFunc::Process(TH1 * const hist)
98 // Fit the mass histogram with the function and retrieve the parameters
101 //reset result arrays
106 AliError("Use 'SetFunctions(TF1*,TF1*)' or 'SetDefaults(Int_t)' to setup the signal functions first'!");
110 //set current signal hist pointer
113 //initial the fit function to its default parameters
114 if (fVInitParam.GetMatrixArray()) fSigBack->SetParameters(fVInitParam.GetMatrixArray());
116 //Perform fit and retrieve values
117 hist->Fit(fSigBack,fFitOpt.Data(),"",fFitMin,fFitMax);
119 //retrieve values and errors
120 //TODO: perhpas implement different methods to retrieve the valus
121 fSignal->SetParameters(fSigBack->GetParameters());
122 fBackground->SetParameters(fSigBack->GetParameters()+fSignal->GetNpar());
124 //TODO: proper error estimate?!?
125 // perhaps this is not possible in a general way?
127 // signal and background histograms
128 TH1 *hSignal=(TH1*)hist->Clone("DieleSignalHist");
129 TH1 *hBackground=(TH1*)hist->Clone("DieleBackgroundHist");
131 fSignal->SetRange(hist->GetXaxis()->GetXmin(),hist->GetXaxis()->GetXmax());
132 fBackground->SetRange(hist->GetXaxis()->GetXmin(),hist->GetXaxis()->GetXmax());
134 hSignal->Add(fBackground,-1);
135 hBackground->Add(fSignal,-1);
138 //get values and errors signal and background
139 Double_t signal=0, signal_err=0;
140 Double_t background=0, background_err=0;
141 Double_t sigPlusBack=0, sigPlusBack_err=0;
144 signal=hSignal->IntegralAndError(hSignal->FindBin(GetIntegralMin()),
145 hSignal->FindBin(GetIntegralMax()), signal_err);
146 background=hBackground->IntegralAndError(hBackground->FindBin(GetIntegralMin()),
147 hBackground->FindBin(GetIntegralMax()), background_err);
148 sigPlusBack=hist->IntegralAndError(hist->FindBin(GetIntegralMin()),
149 hist->FindBin(GetIntegralMax()), sigPlusBack_err);
151 signal=fSignal->Integral(GetIntegralMin(),GetIntegralMax())/hist->GetBinWidth(hist->FindBin((GetIntegralMax()-GetIntegralMin())/2.));
152 background=fBackground->Integral(GetIntegralMin(),GetIntegralMax())/hist->GetBinWidth(hist->FindBin((GetIntegralMax()-GetIntegralMin())/2.));
153 sigPlusBack=fSigBack->Integral(GetIntegralMin(),GetIntegralMax())/hist->GetBinWidth(hist->FindBin((GetIntegralMax()-GetIntegralMin())/2.));
154 //TODO: Error calculation
158 SetSignal(signal,signal_err);
159 SetBackground(background,background_err);
160 SetSignificanceAndSOB();
162 SetMass(fSigBack->GetParameter(fParM), fSigBack->GetParError(fParM));
163 SetMassWidth(fSigBack->GetParameter(fParMres), fSigBack->GetParError(fParMres));
170 //______________________________________________
171 void AliDielectronSignalFunc::Process(TObjArray * const arrhist)
174 // Fit the mass histogram with the function and retrieve the parameters
177 TH1 *hist=(TH1*)arrhist->At(1);
180 //______________________________________________
181 void AliDielectronSignalFunc::SetFunctions(TF1 * const sig, TF1 * const back, TF1 * const combined,
182 Int_t parM, Int_t parMres)
185 // Set the signal, background functions and combined fit function
186 // Note: The process method assumes that the first n parameters in the
187 // combined fit function correspond to the n parameters of the signal function
188 // and the n+1 to n+m parameters to the m parameters of the background function!!!
189 // if parM and parMres are set this information can be used in the drawing function
190 // to also show in the statistics box the mass and mass resolution
192 if (!sig||!back||!combined) {
193 AliError("Both, signal and background function need to be set!");
205 //______________________________________________
206 void AliDielectronSignalFunc::InitParams()
209 // initilise common fit parameters
211 fVInitParam.ResizeTo(fSigBack->GetNpar());
212 fVInitParam.SetElements(fSigBack->GetParameters());
215 //______________________________________________
216 void AliDielectronSignalFunc::SetDefaults(Int_t type)
219 // Setup some default functions:
220 // type = 0: gaus signal + linear background in 2.5 - 4 GeV inv. mass
221 // type = 1: gaus signal + exponential background in 2.5 - 4 GeV inv. mass
222 // type = 2: half gaussian, half exponential signal function
223 // type = 3: Crystal-Ball function
224 // type = 4: Crystal-Ball signal + exponential background
229 fSignal=new TF1("DieleSignal","gaus",2.5,4);
230 fBackground=new TF1("DieleBackground","pol1",2.5,4);
231 fSigBack=new TF1("DieleCombined","gaus+pol1(3)",2.5,4);
233 fSigBack->SetParameters(1,3.1,.05,2.5,1);
234 fSigBack->SetParLimits(0,0,10000000);
235 fSigBack->SetParLimits(1,3.05,3.15);
236 fSigBack->SetParLimits(2,.02,.1);
238 SetFunctions(fSignal,fBackground,fSigBack,1,2);
241 fSignal=new TF1("DieleSignal","gaus",2.5,4);
242 fBackground=new TF1("DieleBackground","[0]*exp(-(x-[1])/[2])",2.5,4);
243 fSigBack=new TF1("DieleCombined","gaus+[3]*exp(-(x-[4])/[5])",2.5,4);
245 fSigBack->SetParameters(1,3.1,.05,1,2.5,1);
246 fSigBack->SetParLimits(0,0,10000000);
247 fSigBack->SetParLimits(1,3.05,3.15);
248 fSigBack->SetParLimits(2,.02,.1);
250 SetFunctions(fSignal,fBackground,fSigBack,1,2);
253 // half gaussian, half exponential signal function
254 // exponential background
255 fSignal = 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);
256 fBackground=new TF1("DieleBackground","[0]*exp(-(x-[1])/[2])+[3]",2.5,4);
257 fSigBack=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);
258 fSigBack->SetParameters(1.,3.1,.05,.1,1,2.5,1,0);
260 fSigBack->SetParLimits(0,0,10000000);
261 fSigBack->SetParLimits(1,3.05,3.15);
262 fSigBack->SetParLimits(2,.02,.1);
263 fSigBack->FixParameter(6,2.5);
264 fSigBack->FixParameter(7,0);
265 SetFunctions(fSignal,fBackground,fSigBack,1,2);
272 //______________________________________________
273 void AliDielectronSignalFunc::Draw(const Option_t* option)
276 // Draw the fitted function
279 TString drawOpt(option);
282 Bool_t optStat=drawOpt.Contains("stat");
284 fSigBack->SetNpx(200);
285 fSigBack->SetRange(GetIntegralMin(),GetIntegralMax());
286 fBackground->SetNpx(200);
287 fBackground->SetRange(GetIntegralMin(),GetIntegralMax());
289 TGraph *grSig=new TGraph(fSigBack);
290 grSig->SetFillColor(kGreen);
291 grSig->SetFillStyle(3001);
293 TGraph *grBack=new TGraph(fBackground);
294 grBack->SetFillColor(kRed);
295 grBack->SetFillStyle(3001);
297 grSig->SetPoint(0,grBack->GetX()[0],grBack->GetY()[0]);
298 grSig->SetPoint(grSig->GetN()-1,grBack->GetX()[grBack->GetN()-1],grBack->GetY()[grBack->GetN()-1]);
300 grBack->SetPoint(0,grBack->GetX()[0],0.);
301 grBack->SetPoint(grBack->GetN()-1,grBack->GetX()[grBack->GetN()-1],0.);
303 fSigBack->SetRange(fFitMin,fFitMax);
304 fBackground->SetRange(fFitMin,fFitMax);
306 if (!drawOpt.Contains("same")){
317 fSigBack->Draw("same");
318 fBackground->Draw("same");
320 if (optStat) DrawStats();