]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - PWGDQ/dielectron/AliDielectronSignalFunc.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
filtering updates
[u/mrichter/AliRoot.git] / PWGDQ / dielectron / AliDielectronSignalFunc.cxx
... / ...
CommitLineData
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
22 Class used for extracting the signal from an invariant mass spectrum.
23 It implements the AliDielectronSignalBase and -Ext classes and it uses user provided
24 functions to fit the spectrum with a combined signa+background fit.
25 Used invariant mass spectra are provided via an array of histograms. There are serveral method
26 to estimate the background and to extract the raw yield from the background subtracted spectra.
27
28 Example usage:
29
30 AliDielectronSignalFunc *sig = new AliDielectronSignalFunc();
31
32
33 1) invariant mass input spectra
34
35 1.1) Assuming a AliDielectronCF container as data format (check class for more details)
36 AliDielectronCFdraw *cf = new AliDielectronCFdraw("path/to/the/output/file.root");
37 TObjArray *arrHists = cf->CollectMinvProj(cf->FindStep("Config"));
38
39 1.2) Assuming a AliDielectronHF grid as data format (check class for more details)
40 AliDielectronHFhelper *hf = new AliDielectronHFhelper("path/to/the/output/file.root", "ConfigName");
41 TObjArray *arrHists = hf->CollectHistos(AliDielectronVarManager::kM);
42
43 1.3) Assuming a single histograms
44 TObjArray *histoArray = new TObjArray();
45 arrHists->Add(signalPP); // add the spectrum histograms to the array
46 arrHists->Add(signalPM); // the order is important !!!
47 arrHists->Add(signalMM);
48
49
50 2) background estimation
51
52 2.1) set the method for the background estimation (methods can be found in AliDielectronSignalBase)
53 sig->SetMethod(AliDielectronSignalBase::kFitted);
54 2.2) rebin the spectras if needed
55 // sig->SetRebin(2);
56 2.3) add any background function you like
57 TF1 *fB = new TF1("fitBgrd","pol3",minFit,maxFit);
58
59
60 3) configure the signal extraction
61
62 3.1) chose one of the signal functions (MCshape, CrystalBall, Gauss)
63 TF1 *fS = new TF1("fitSign",AliDielectronSignalFunc::PeakFunCB,minFit,maxFit,5); // has 5 parameters
64 // TF1 *fS = new TF1("fitSign",AliDielectronSignalFunc::PeakFunGaus,minFit,maxFit,3); // has 3 parameters
65 // sig->SetMCSignalShape(hMCsign);
66 // TF1 *fS = new TF1("fitSign",AliDielectronSignalFunc::PeakFunMC,minFit,maxFit,1); // requires a MC shape
67 3.2) set the method for the signal extraction (methods can be found in AliDielectronSignalBase)
68 depending on the method serveral inputs are needed (e.g. MC shape, PDG code of the particle of interest)
69 // sig->SetParticleOfInterest(443); //default is jpsi
70 // sig->SetMCSignalShape(signalMC);
71 // sig->SetIntegralRange(minInt, maxInt);
72 sig->SetExtractionMethod(AliDielectronSignal::BinCounting); // this is the default
73
74
75 4) combined fit of bgrd+signal
76
77 4.1) combine the two functions
78 sig->CombineFunc(fS,fB);
79 4.2) apply fitting ranges and the fit options
80 sig->SetFitRange(minFit, maxFit);
81 sig->SetFitOption("NR");
82
83
84 5) start the processing
85
86 sig->Process(arrHists);
87 sig->Print(""); // print values and errors extracted
88
89
90 6) access the spectra and values created
91
92 6.1) standard spectra as provided filled in AliDielectronSignalExt
93 TH1F *hsign = (TH1F*) sig->GetUnlikeSignHistogram(); // same as the input (rebinned)
94 TH1F *hbgrd = (TH1F*) sig->GetBackgroundHistogram(); // filled histogram with fitBgrd
95 TH1F *hextr = (TH1F*) sig->GetSignalHistogram(); // after backgound extraction (rebinned)
96 TObject *oPeak = (TObject*) sig->GetPeakShape(); // can be a TF1 or TH1 depending on the method
97 6.2) flow spectra
98 TF1 *fFitSign = sig->GetCombinedFunction(); // combined fit function
99 TF1 *fFitExtr = sig->GetSignalFunction(); // signal function
100 TF1 *fFitBgrd = sig->GetBackgroundFunction(); // background function
101 6.3) access the extracted values and errors
102 sig->GetValues(); or GetErrors(); // yield extraction
103
104*/
105// //
106///////////////////////////////////////////////////////////////////////////
107
108#include <TF1.h>
109#include <TH1.h>
110#include <TGraph.h>
111#include <TMath.h>
112#include <TString.h>
113#include <TPaveText.h>
114#include <TList.h>
115#include <TFitResult.h>
116//#include <../hist/hist/src/TF1Helper.h>
117
118#include <AliLog.h>
119
120#include "AliDielectronSignalFunc.h"
121
122ClassImp(AliDielectronSignalFunc)
123
124//TH1F* AliDielectronSignalFunc::fgHistSimPM=0x0;
125TF1* AliDielectronSignalFunc::fFuncSignal=0x0;
126TF1* AliDielectronSignalFunc::fFuncBackground=0x0;
127Int_t AliDielectronSignalFunc::fNparPeak=0;
128Int_t AliDielectronSignalFunc::fNparBgnd=0;
129
130
131AliDielectronSignalFunc::AliDielectronSignalFunc() :
132AliDielectronSignalExt(),
133fFuncSigBack(0x0),
134fParMass(1),
135fParMassWidth(2),
136fFitOpt("SMNQE"),
137fUseIntegral(kFALSE),
138fDof(0),
139fChi2Dof(0.0)
140{
141 //
142 // Default Constructor
143 //
144
145}
146
147//______________________________________________
148AliDielectronSignalFunc::AliDielectronSignalFunc(const char* name, const char* title) :
149AliDielectronSignalExt(name, title),
150fFuncSigBack(0x0),
151fParMass(1),
152fParMassWidth(2),
153fFitOpt("SMNQE"),
154fUseIntegral(kFALSE),
155fDof(0),
156fChi2Dof(0.0)
157{
158 //
159 // Named Constructor
160 //
161}
162
163//______________________________________________
164AliDielectronSignalFunc::~AliDielectronSignalFunc()
165{
166 //
167 // Default Destructor
168 //
169 if(fFuncSigBack) delete fFuncSigBack;
170}
171
172
173//______________________________________________
174void AliDielectronSignalFunc::Process(TObjArray * const arrhist)
175{
176 //
177 // Fit the invariant mass histograms and retrieve the signal and background
178 //
179 switch(fMethod) {
180 case kFitted :
181 ProcessFit(arrhist);
182 break;
183
184 case kLikeSignFit :
185 ProcessFitLS(arrhist);
186 break;
187
188 case kEventMixingFit :
189 ProcessFitEM(arrhist);
190 break;
191
192 case kLikeSign :
193 case kLikeSignArithm :
194 case kLikeSignRcorr:
195 case kLikeSignArithmRcorr:
196 ProcessLS(arrhist); // process like-sign subtraction method
197 break;
198
199 case kEventMixing :
200 ProcessEM(arrhist); // process event mixing method
201 break;
202
203 case kRotation:
204 ProcessRotation(arrhist);
205 break;
206
207 default :
208 AliError("Background substraction method not known!");
209 }
210}
211//______________________________________________________________________________
212Double_t AliDielectronSignalFunc::PeakFunMC(const Double_t *x, const Double_t *par) {
213 // Fit MC signal shape
214 // parameters
215 // [0]: scale for simpeak
216
217 Double_t xx = x[0];
218
219 TH1F *hPeak = fgHistSimPM;
220 if (!hPeak) {
221 printf("E-AliDielectronSignalFunc::PeakFun: No histogram for peak fit defined!\n");
222 return 0.0;
223 }
224
225 Int_t idx = hPeak->FindBin(xx);
226 if ((idx <= 1) ||
227 (idx >= hPeak->GetNbinsX())) {
228 return 0.0;
229 }
230
231 return (par[0] * hPeak->GetBinContent(idx));
232
233}
234
235//______________________________________________________________________________
236Double_t AliDielectronSignalFunc::PeakFunCB(const Double_t *x, const Double_t *par) {
237 // Crystal Ball fit function
238
239 Double_t alpha = par[0];
240 Double_t n = par[1];
241 Double_t meanx = par[2];
242 Double_t sigma = par[3];
243 Double_t nn = par[4];
244
245 Double_t a = TMath::Power((n/TMath::Abs(alpha)), n) * TMath::Exp(-.5*alpha*alpha);
246 Double_t b = n/TMath::Abs(alpha) - TMath::Abs(alpha);
247
248 Double_t arg = (x[0] - meanx)/sigma;
249 Double_t fitval = 0;
250
251 if (arg > -1.*alpha) {
252 fitval = nn * TMath::Exp(-.5*arg*arg);
253 } else {
254 fitval = nn * a * TMath::Power((b-arg), (-1*n));
255 }
256
257 return fitval;
258}
259
260//______________________________________________________________________________
261Double_t AliDielectronSignalFunc::PeakFunGaus(const Double_t *x, const Double_t *par) {
262 // Gaussian fit function
263 //printf("fNparBgrd %d \n",fNparBgnd);
264 Double_t n = par[0];
265 Double_t mean = par[1];
266 Double_t sigma = par[2];
267 Double_t xx = x[0];
268
269 return ( n*TMath::Exp(-0.5*TMath::Power((xx-mean)/sigma,2)) );
270}
271
272//______________________________________________
273void AliDielectronSignalFunc::ProcessFit(TObjArray * const arrhist) {
274 //
275 // Fit the +- invariant mass distribution only
276 // Here we assume that the combined fit function is a sum of the signal and background functions
277 // and that the signal function is always the first term of this sum
278 //
279
280 fHistDataPM = (TH1F*)(arrhist->At(1))->Clone("histPM"); // +- SE
281 fHistDataPM->Sumw2();
282 if(fRebin>1)
283 fHistDataPM->Rebin(fRebin);
284
285 fHistSignal = new TH1F("HistSignal", "Fit substracted signal",
286 fHistDataPM->GetXaxis()->GetNbins(),
287 fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
288 fHistBackground = new TH1F("HistBackground", "Fit contribution",
289 fHistDataPM->GetXaxis()->GetNbins(),
290 fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
291
292 // the starting parameters of the fit function and their limits can be tuned
293 // by the user in its macro
294 fHistDataPM->Fit(fFuncSigBack, fFitOpt.Data(), "", fFitMin, fFitMax);
295 TFitResultPtr pmFitPtr = fHistDataPM->Fit(fFuncSigBack, fFitOpt.Data(), "", fFitMin, fFitMax);
296 //TFitResult *pmFitResult = pmFitPtr.Get(); // used only with TF1Helper
297 //fFuncBackground->SetParameters(fFuncSigBack->GetParameters());
298 fFuncSignal->SetParameters(fFuncSigBack->GetParameters());
299 fFuncBackground->SetParameters(fFuncSigBack->GetParameters()+fFuncSignal->GetNpar());
300
301 // fill the background spectrum
302 fHistBackground->Eval(fFuncBackground);
303 // set the error for the background histogram
304 fHistBackground->Fit(fFuncBackground,"0qR","",fFitMin,fFitMax);
305 Double_t inte = fFuncBackground->IntegralError(fIntMin, fIntMax)/fHistDataPM->GetBinWidth(1);
306 Double_t binte = inte / TMath::Sqrt((fHistDataPM->FindBin(fIntMax)-fHistDataPM->FindBin(fIntMin))+1);
307 for(Int_t iBin=fHistDataPM->FindBin(fIntMin); iBin<=fHistDataPM->FindBin(fIntMax); iBin++) {
308 fHistBackground->SetBinError(iBin, binte);
309 }
310
311 for(Int_t iBin=1; iBin<=fHistDataPM->GetXaxis()->GetNbins(); iBin++) {
312 // Double_t m = fHistDataPM->GetBinCenter(iBin);
313 Double_t pm = fHistDataPM->GetBinContent(iBin);
314 Double_t epm = fHistDataPM->GetBinError(iBin);
315 Double_t bknd = fHistBackground->GetBinContent(iBin);
316 Double_t ebknd = fHistBackground->GetBinError(iBin);
317 for(Int_t iPar=fFuncSignal->GetNpar(); iPar<fFuncSigBack->GetNpar(); iPar++) {
318 /* TF1Helper problem on alien compilation
319 for(Int_t jPar=iPar; jPar<fFuncSigBack->GetNpar(); jPar++) {
320 TF1 gradientIpar("gradientIpar",
321 ROOT::TF1Helper::TGradientParFunction(iPar-fFuncSignal->GetNpar(),fFuncBackground),0,0,0);
322 TF1 gradientJpar("gradientJpar",
323 ROOT::TF1Helper::TGradientParFunction(jPar-fFuncSignal->GetNpar(),fFuncBackground),0,0,0);
324 ebknd += pmFitResult->CovMatrix(iPar,jPar)*
325 gradientIpar.Eval(m)*gradientJpar.Eval(m)*
326 (iPar==jPar ? 1.0 : 2.0);
327 }
328 */
329 }
330 Double_t signal = pm-bknd;
331 Double_t error = TMath::Sqrt(epm*epm+ebknd);
332 fHistSignal->SetBinContent(iBin, signal);
333 fHistSignal->SetBinError(iBin, error);
334 }
335
336 if(fUseIntegral) {
337 // signal
338 fValues(0) = fFuncSignal->Integral(fIntMin, fIntMax)/fHistDataPM->GetBinWidth(1);
339 fErrors(0) = 0;
340 for(Int_t iPar=0; iPar<fFuncSignal->GetNpar(); iPar++) {
341/* TF1Helper problem on alien compilation
342 for(Int_t jPar=iPar; jPar<fFuncSignal->GetNpar(); jPar++) {
343 TF1 gradientIpar("gradientIpar",
344 ROOT::TF1Helper::TGradientParFunction(iPar,fFuncSignal),0,0,0);
345 TF1 gradientJpar("gradientJpar",
346 ROOT::TF1Helper::TGradientParFunction(jPar,fFuncSignal),0,0,0);
347 fErrors(0) += pmFitResult->CovMatrix(iPar,jPar)*
348 gradientIpar.Integral(fIntMin,fIntMax)*gradientJpar.Integral(fIntMin,fIntMax)*
349 (iPar==jPar ? 1.0 : 2.0);
350 }
351*/
352 }
353 // background
354 fValues(1) = fFuncBackground->Integral(fIntMin, fIntMax)/fHistDataPM->GetBinWidth(1);
355 fErrors(1) = fFuncBackground->IntegralError(fIntMin, fIntMax)/fHistDataPM->GetBinWidth(1);
356 for(Int_t iPar=fFuncSignal->GetNpar(); iPar<fFuncSigBack->GetNpar(); iPar++) {
357/* TF1Helper problem on alien compilation
358 for(Int_t jPar=iPar; jPar<fFuncSigBack->GetNpar(); jPar++) {
359 TF1 gradientIpar("gradientIpar",
360 ROOT::TF1Helper::TGradientParFunction(iPar-fFuncSignal->GetNpar(),fFuncBackground),0,0,0);
361 TF1 gradientJpar("gradientJpar",
362 ROOT::TF1Helper::TGradientParFunction(jPar-fFuncSignal->GetNpar(),fFuncBackground),0,0,0);
363 fErrors(1) += pmFitResult->CovMatrix(iPar,jPar)*
364 gradientIpar.Integral(fIntMin, fIntMax)*gradientJpar.Integral(fIntMin, fIntMax)*
365 (iPar==jPar ? 1.0 : 2.0);
366 }
367*/
368 }
369 }
370 else {
371 // signal
372 fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
373 fHistSignal->FindBin(fIntMax), fErrors(0));
374 // background
375 fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
376 fHistBackground->FindBin(fIntMax),
377 fErrors(1));
378 }
379 // S/B and significance
380 SetSignificanceAndSOB();
381 fValues(4) = fFuncSigBack->GetParameter(fParMass);
382 fErrors(4) = fFuncSigBack->GetParError(fParMass);
383 fValues(5) = fFuncSigBack->GetParameter(fParMassWidth);
384 fErrors(5) = fFuncSigBack->GetParError(fParMassWidth);
385
386 fProcessed = kTRUE;
387
388 fHistBackground->GetListOfFunctions()->Add(fFuncBackground);
389
390}
391
392//______________________________________________
393void AliDielectronSignalFunc::ProcessFitLS(TObjArray * const arrhist) {
394 //
395 // Substract background using the like-sign spectrum
396 //
397 fHistDataPP = (TH1F*)(arrhist->At(0))->Clone("histPP"); // ++ SE
398 fHistDataPM = (TH1F*)(arrhist->At(1))->Clone("histPM"); // +- SE
399 fHistDataMM = (TH1F*)(arrhist->At(2))->Clone("histMM"); // -- SE
400 if (fRebin>1) {
401 fHistDataPP->Rebin(fRebin);
402 fHistDataPM->Rebin(fRebin);
403 fHistDataMM->Rebin(fRebin);
404 }
405 fHistDataPP->Sumw2();
406 fHistDataPM->Sumw2();
407 fHistDataMM->Sumw2();
408
409 fHistSignal = new TH1F("HistSignal", "Like-Sign substracted signal",
410 fHistDataPM->GetXaxis()->GetNbins(),
411 fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
412 fHistBackground = new TH1F("HistBackground", "Like-sign contribution",
413 fHistDataPM->GetXaxis()->GetNbins(),
414 fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
415
416 // fit the +- mass distribution
417 fHistDataPM->Fit(fFuncSigBack, fFitOpt.Data(), "", fFitMin, fFitMax);
418 fHistDataPM->Fit(fFuncSigBack, fFitOpt.Data(), "", fFitMin, fFitMax);
419 // declare the variables where the like-sign fit results will be stored
420// TFitResult *ppFitResult = 0x0;
421// TFitResult *mmFitResult = 0x0;
422 // fit the like sign background
423 TF1 *funcClonePP = (TF1*)fFuncBackground->Clone("funcClonePP");
424 TF1 *funcCloneMM = (TF1*)fFuncBackground->Clone("funcCloneMM");
425 fHistDataPP->Fit(funcClonePP, fFitOpt.Data(), "", fFitMin, fFitMax);
426 fHistDataPP->Fit(funcClonePP, fFitOpt.Data(), "", fFitMin, fFitMax);
427// TFitResultPtr ppFitPtr = fHistDataPP->Fit(funcClonePP, fFitOpt.Data(), "", fFitMin, fFitMax);
428// ppFitResult = ppFitPtr.Get();
429 fHistDataMM->Fit(funcCloneMM, fFitOpt.Data(), "", fFitMin, fFitMax);
430 fHistDataMM->Fit(funcCloneMM, fFitOpt.Data(), "", fFitMin, fFitMax);
431// TFitResultPtr mmFitPtr = fHistDataMM->Fit(funcCloneMM, fFitOpt.Data(), "", fFitMin, fFitMax);
432// mmFitResult = mmFitPtr.Get();
433
434 for(Int_t iBin=1; iBin<=fHistDataPM->GetXaxis()->GetNbins(); iBin++) {
435 Double_t m = fHistDataPM->GetBinCenter(iBin);
436 Double_t pm = fHistDataPM->GetBinContent(iBin);
437 Double_t pp = funcClonePP->Eval(m);
438 Double_t mm = funcCloneMM->Eval(m);
439 Double_t epm = fHistDataPM->GetBinError(iBin);
440 Double_t epp = 0;
441 for(Int_t iPar=0; iPar<funcClonePP->GetNpar(); iPar++) {
442/* TF1Helper problem on alien compilation
443 for(Int_t jPar=iPar; jPar<funcClonePP->GetNpar(); jPar++) {
444 TF1 gradientIpar("gradientIpar",
445 ROOT::TF1Helper::TGradientParFunction(iPar,funcClonePP),0,0,0);
446 TF1 gradientJpar("gradientJpar",
447 ROOT::TF1Helper::TGradientParFunction(jPar,funcClonePP),0,0,0);
448 epp += ppFitResult->CovMatrix(iPar,jPar)*
449 gradientIpar.Eval(m)*gradientJpar.Eval(m)*
450 (iPar==jPar ? 1.0 : 2.0);
451 }
452*/
453 }
454 Double_t emm = 0;
455 for(Int_t iPar=0; iPar<funcCloneMM->GetNpar(); iPar++) {
456/* TF1Helper problem on alien compilation
457 for(Int_t jPar=iPar; jPar<funcCloneMM->GetNpar(); jPar++) {
458 TF1 gradientIpar("gradientIpar",
459 ROOT::TF1Helper::TGradientParFunction(iPar,funcCloneMM),0,0,0);
460 TF1 gradientJpar("gradientJpar",
461 ROOT::TF1Helper::TGradientParFunction(jPar,funcCloneMM),0,0,0);
462 emm += mmFitResult->CovMatrix(iPar,jPar)*
463 gradientIpar.Eval(m)*gradientJpar.Eval(m)*
464 (iPar==jPar ? 1.0 : 2.0);
465 }
466*/
467 }
468
469 Double_t signal = pm-2.0*TMath::Sqrt(pp*mm);
470 Double_t background = 2.0*TMath::Sqrt(pp*mm);
471 // error propagation on the signal calculation above
472 Double_t esignal = TMath::Sqrt(epm*epm+(mm/pp)*epp+(pp/mm)*emm);
473 Double_t ebackground = TMath::Sqrt((mm/pp)*epp+(pp/mm)*emm);
474 fHistSignal->SetBinContent(iBin, signal);
475 fHistSignal->SetBinError(iBin, esignal);
476 fHistBackground->SetBinContent(iBin, background);
477 fHistBackground->SetBinError(iBin, ebackground);
478 }
479
480 // signal
481 fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
482 fHistSignal->FindBin(fIntMax), fErrors(0));
483 // background
484 fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
485 fHistBackground->FindBin(fIntMax),
486 fErrors(1));
487 // S/B and significance
488 SetSignificanceAndSOB();
489 fValues(4) = fFuncSigBack->GetParameter(fParMass);
490 fErrors(4) = fFuncSigBack->GetParError(fParMass);
491 fValues(5) = fFuncSigBack->GetParameter(fParMassWidth);
492 fErrors(5) = fFuncSigBack->GetParError(fParMassWidth);
493
494 fProcessed = kTRUE;
495}
496
497//______________________________________________
498void AliDielectronSignalFunc::ProcessFitEM(TObjArray * const arrhist) {
499 //
500 // Substract background with the event mixing technique
501 //
502 arrhist->GetEntries(); // just to avoid the unused parameter warning
503 AliError("Event mixing for background substraction method not implemented!");
504}
505
506//______________________________________________
507void AliDielectronSignalFunc::SetFunctions(TF1 * const combined, TF1 * const sig, TF1 * const back,
508 Int_t parM, Int_t parMres)
509{
510 //
511 // Set the signal, background functions and combined fit function
512 // Note: The process method assumes that the first n parameters in the
513 // combined fit function correspond to the n parameters of the signal function
514 // and the n+1 to n+m parameters to the m parameters of the background function!!!
515
516 if (!sig||!back||!combined) {
517 AliError("Both, signal and background function need to be set!");
518 return;
519 }
520 fFuncSignal=sig;
521 fFuncBackground=back;
522 fFuncSigBack=combined;
523 fParMass=parM;
524 fParMassWidth=parMres;
525
526}
527
528//______________________________________________
529void AliDielectronSignalFunc::SetDefaults(Int_t type)
530{
531 //
532 // Setup some default functions:
533 // type = 0: gaus signal + linear background in 2.5 - 4 GeV inv. mass
534 // type = 1: gaus signal + exponential background in 2.5 - 4 GeV inv. mass
535 // type = 2: half gaussian, half exponential signal function
536 // type = 3: Crystal-Ball function
537 // type = 4: Crystal-Ball signal + exponential background
538 //
539
540 if (type==0){
541 fFuncSignal=new TF1("DieleSignal","gaus",2.5,4);
542 fFuncBackground=new TF1("DieleBackground","pol1",2.5,4);
543 fFuncSigBack=new TF1("DieleCombined","gaus+pol1(3)",2.5,4);
544
545 fFuncSigBack->SetParameters(1,3.1,.05,2.5,1);
546 fFuncSigBack->SetParLimits(0,0,10000000);
547 fFuncSigBack->SetParLimits(1,3.05,3.15);
548 fFuncSigBack->SetParLimits(2,.02,.1);
549 }
550 else if (type==1){
551 fFuncSignal=new TF1("DieleSignal","gaus",2.5,4);
552 fFuncBackground=new TF1("DieleBackground","[0]*exp(-(x-[1])/[2])",2.5,4);
553 fFuncSigBack=new TF1("DieleCombined","gaus+[3]*exp(-(x-[4])/[5])",2.5,4);
554
555 fFuncSigBack->SetParameters(1,3.1,.05,1,2.5,1);
556 fFuncSigBack->SetParLimits(0,0,10000000);
557 fFuncSigBack->SetParLimits(1,3.05,3.15);
558 fFuncSigBack->SetParLimits(2,.02,.1);
559 }
560 else if (type==2){
561 // half gaussian, half exponential signal function
562 // exponential background
563 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);
564 fFuncBackground = new TF1("DieleBackground","[0]*exp(-(x-[1])/[2])+[3]",2.5,4);
565 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);
566 fFuncSigBack->SetParameters(1.,3.1,.05,.1,1,2.5,1,0);
567
568 fFuncSigBack->SetParLimits(0,0,10000000);
569 fFuncSigBack->SetParLimits(1,3.05,3.15);
570 fFuncSigBack->SetParLimits(2,.02,.1);
571 fFuncSigBack->FixParameter(6,2.5);
572 fFuncSigBack->FixParameter(7,0);
573 }
574}
575
576
577//______________________________________________
578void AliDielectronSignalFunc::Draw(const Option_t* option)
579{
580 //
581 // Draw the fitted function
582 //
583
584 TString drawOpt(option);
585 drawOpt.ToLower();
586
587 Bool_t optStat=drawOpt.Contains("stat");
588
589 fFuncSigBack->SetNpx(200);
590 fFuncSigBack->SetRange(fIntMin,fIntMax);
591 fFuncBackground->SetNpx(200);
592 fFuncBackground->SetRange(fIntMin,fIntMax);
593
594 TGraph *grSig=new TGraph(fFuncSigBack);
595 grSig->SetFillColor(kGreen);
596 grSig->SetFillStyle(3001);
597
598 TGraph *grBack=new TGraph(fFuncBackground);
599 grBack->SetFillColor(kRed);
600 grBack->SetFillStyle(3001);
601
602 grSig->SetPoint(0,grBack->GetX()[0],grBack->GetY()[0]);
603 grSig->SetPoint(grSig->GetN()-1,grBack->GetX()[grBack->GetN()-1],grBack->GetY()[grBack->GetN()-1]);
604
605 grBack->SetPoint(0,grBack->GetX()[0],0.);
606 grBack->SetPoint(grBack->GetN()-1,grBack->GetX()[grBack->GetN()-1],0.);
607
608 fFuncSigBack->SetRange(fFitMin,fFitMax);
609 fFuncBackground->SetRange(fFitMin,fFitMax);
610
611 if (!drawOpt.Contains("same")){
612 if (fHistDataPM){
613 fHistDataPM->Draw();
614 grSig->Draw("f");
615 } else {
616 grSig->Draw("af");
617 }
618 } else {
619 grSig->Draw("f");
620 }
621 if(fMethod==kFitted || fMethod==kFittedMC) grBack->Draw("f");
622 fFuncSigBack->Draw("same");
623 fFuncSigBack->SetLineWidth(2);
624 if(fMethod==kLikeSign) {
625 fHistDataPP->SetLineWidth(2);
626 fHistDataPP->SetLineColor(6);
627 fHistDataPP->Draw("same");
628 fHistDataMM->SetLineWidth(2);
629 fHistDataMM->SetLineColor(8);
630 fHistDataMM->Draw("same");
631 }
632
633 if(fMethod==kFitted || fMethod==kFittedMC)
634 fFuncBackground->Draw("same");
635
636 if (optStat) DrawStats();
637
638}
639
640
641//______________________________________________________________________________
642void AliDielectronSignalFunc::CombineFunc(TF1 * const peak, TF1 * const bgnd) {
643 //
644 // combine the bgnd and the peak function
645 //
646
647 if (!peak||!bgnd) {
648 AliError("Both, signal and background function need to be set!");
649 return;
650 }
651 fFuncSignal=peak;
652 fFuncBackground=bgnd;
653
654 fNparPeak = fFuncSignal->GetNpar();
655 fNparBgnd = fFuncBackground->GetNpar();
656
657 fFuncSigBack = new TF1("BgndPeak",AliDielectronSignalFunc::PeakBgndFun, fFitMin,fFitMax, fNparPeak+fNparBgnd);
658 return;
659}
660
661//______________________________________________________________________________
662Double_t AliDielectronSignalFunc::PeakBgndFun(const Double_t *x, const Double_t *par) {
663 //
664 // merge peak and bgnd functions
665 //
666 return (fFuncSignal->EvalPar(x,par) + fFuncBackground->EvalPar(x,par+fNparPeak));
667}
668