]> 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
o update dielectron package
[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/*
21Dielectron signal extraction class using functions as input.
22
23A function to describe the signal as well as one to describe the background
24has to be deployed by the user. Alternatively on of the default implementaions
25can 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>
37#include <TList.h>
38#include <TFitResult.h>
39//#include <../hist/hist/src/TF1Helper.h>
40
41#include <AliLog.h>
42
43#include "AliDielectronSignalFunc.h"
44
45ClassImp(AliDielectronSignalFunc)
46
47TH1F* AliDielectronSignalFunc::fgHistSimPM=0x0;
48
49AliDielectronSignalFunc::AliDielectronSignalFunc() :
50AliDielectronSignalBase(),
51fFuncSignal(0x0),
52fFuncBackground(0x0),
53fFuncSigBack(0x0),
54fParMass(1),
55fParMassWidth(2),
56fFitOpt("SMNQE"),
57fUseIntegral(kFALSE),
58fPolDeg(0),
59fChi2Dof(0.0)
60{
61 //
62 // Default Constructor
63 //
64
65}
66
67//______________________________________________
68AliDielectronSignalFunc::AliDielectronSignalFunc(const char* name, const char* title) :
69AliDielectronSignalBase(name, title),
70fFuncSignal(0x0),
71fFuncBackground(0x0),
72fFuncSigBack(0x0),
73fParMass(1),
74fParMassWidth(2),
75fFitOpt("SMNQE"),
76fUseIntegral(kFALSE),
77fPolDeg(0),
78fChi2Dof(0.0)
79{
80 //
81 // Named Constructor
82 //
83}
84
85//______________________________________________
86AliDielectronSignalFunc::~AliDielectronSignalFunc()
87{
88 //
89 // Default Destructor
90 //
91 if(fFuncSignal) delete fFuncSignal;
92 if(fFuncBackground) delete fFuncBackground;
93 if(fFuncSigBack) delete fFuncSigBack;
94}
95
96
97//______________________________________________
98void AliDielectronSignalFunc::Process(TObjArray * const arrhist)
99{
100 //
101 // Fit the invariant mass histograms and retrieve the signal and background
102 //
103 switch(fMethod) {
104 case kFittedMC :
105 ProcessFitIKF(arrhist);
106 break;
107
108 case kFitted :
109 ProcessFit(arrhist);
110 break;
111
112 case kLikeSign :
113 ProcessLS(arrhist);
114 break;
115
116 case kEventMixing :
117 ProcessEM(arrhist);
118 break;
119
120 default :
121 AliError("Background substraction method not known!");
122 }
123}
124
125//______________________________________________
126void AliDielectronSignalFunc::ProcessFitIKF(TObjArray * const arrhist) {
127 //
128 // Fit the +- invariant mass distribution only
129 //
130 //
131
132 const Double_t bigNumber = 100000.;
133 Double_t chi2ndfm1 = bigNumber;
134 Double_t ratiom1 = bigNumber;
135 Double_t chi2ndf = bigNumber;
136 Int_t nDP =0;
137
138 Int_t maxPolDeg = 8;
139
140 fHistDataPM = (TH1F*)(arrhist->At(1))->Clone("histPM"); // +- SE
141 if(fRebin>1) fHistDataPM->Rebin(fRebin);
142
143 fgHistSimPM = (TH1F*)(arrhist->At(3))->Clone("histPMsim"); // +- mc shape
144 if (!fgHistSimPM) {
145 AliFatal("No mc peak shape found at idx 3.");
146 return;
147 }
148 if(fRebin>1) fgHistSimPM->Rebin(fRebin);
149
150 // try out the polynomial degrees
151 for (Int_t iPD=0; iPD<=maxPolDeg; iPD++) {
152 TH1F *hf1 = (TH1F *) fHistDataPM->Clone(Form("hf1_PD%d",iPD));
153
154 FitOneMinv(hf1, fgHistSimPM, iPD);
155 if (fChi2Dof > 0) chi2ndf = fChi2Dof;
156 AliInfo(Form("nDP: %d, iPD: %d, chi2ndf: %f", nDP, iPD, chi2ndf));
157
158 ratiom1 = TMath::Abs(fChi2Dof - 1);
159 if (chi2ndfm1 > ratiom1) { // search for the closest to 1.
160 chi2ndfm1 = ratiom1;
161 nDP = iPD;
162 }
163 }
164
165
166 // fit again with the best polynomial degree
167 TH1F *h2 = (TH1F *) fHistDataPM->Clone(Form("h2_PD%d",nDP));
168
169 FitOneMinv(h2, fgHistSimPM, nDP);
170 AliInfo(Form("Best Fit: PD %d, chi^2/ndf %.3f, S/B %.3f",nDP,fChi2Dof,fValues(3)));
171
172}
173//______________________________________________
174void AliDielectronSignalFunc::FitOneMinv(TH1F *hMinv, TH1F *hSim, Int_t pod) {
175 //
176 // main function to fit an inv mass spectrum
177 //
178
179 TObjArray *arrResults = new TObjArray;
180 arrResults->SetOwner();
181 arrResults->AddAt(hMinv,0);
182
183 // Degree of polynomial
184 fPolDeg = pod;
185
186 // inclusion and exclusion areas (values)
187 const Double_t kJPMass = 3.096916;
188 // inclusion and exclusion areas (bin numbers)
189 const Int_t binIntLo = hMinv->FindBin(fIntMin);
190 const Int_t binIntHi = hMinv->FindBin(fIntMax);
191 // for error calculation
192 Double_t intAreaEdgeLo = hMinv->GetBinLowEdge(binIntLo);
193 Double_t intAreaEdgeHi = hMinv->GetBinLowEdge(binIntHi)+hMinv->GetBinWidth(binIntHi);
194 Double_t norm = (binIntHi-binIntLo)/(fIntMax - fIntMin);
195
196 TH1F *hBfit = (TH1F *) hMinv->Clone(); // for bg only fit (excluding peak region)
197 TH1F *hSigF = (TH1F *) hMinv->Clone(); // signal with subtracted bg
198 TH1F *hBgrd = (TH1F *) hMinv->Clone(); // bg histogram
199
200 hBfit->Reset();
201 hSigF->Reset();
202 hBgrd->Reset();
203
204 // extract start parameter for the MC signal fit
205 Double_t bgvalAv = (hMinv->Integral(1,hMinv->GetNbinsX()+1) - hMinv->Integral(binIntLo,binIntHi)) / (hMinv->GetNbinsX()+1 - (binIntHi-binIntLo));
206 Double_t pkval = hMinv->GetBinContent(hMinv->FindBin(kJPMass)) - bgvalAv;
207 Double_t heightMC = hSim->GetBinContent(hSim->FindBin(kJPMass));
208 Double_t peakScale = (heightMC > 0. ? pkval/heightMC : 0.0);
209
210 Int_t nBgnd = 2 + fPolDeg; // degree + c1st oefficient + higher coefficients
211 Int_t nMinv = nBgnd + 1; // bgrd + peakscale
212
213 // Create the spectra without peak region
214 for (Int_t iBin = 0; iBin <= hMinv->GetNbinsX(); iBin++) {
215 if ((iBin < binIntLo) || (iBin > binIntHi)) {
216 hBfit->SetBinContent(iBin,hMinv->GetBinContent(iBin));
217 hBfit->SetBinError(iBin,hMinv->GetBinError(iBin));
218 }
219 }
220
221
222 // =======
223 // 1.
224 // =======
225 // Do the fit to the background spectrum
226 TF1 *fBo = new TF1("bgrd_fit",BgndFun,fFitMin,fFitMax,nBgnd);
227 for (Int_t iPar=0; iPar<nBgnd; iPar++) {
228 if (iPar == 0) fBo->FixParameter(0, fPolDeg);
229 if (iPar == 1) fBo->SetParameter(iPar, bgvalAv);
230 if (iPar >= 2) fBo->SetParameter(iPar, 0.);
231 }
232 hBfit->Fit(fBo,"0qR");
233 //hBfit->SetNameTitle("bgrd_fit");
234 arrResults->AddAt(fBo,1);
235
236
237 // =======
238 // 2.
239 // =======
240 // Fit the whole spectrum with peak and background
241 TF1 *fSB = new TF1("bgrd_peak_fit",MinvFun,fFitMin,fFitMax,nMinv);
242 fSB->FixParameter(0, fPolDeg);
243 fSB->SetParameter(1, peakScale);
244 // copy the polynomial parameters
245 for (Int_t iPar=0; iPar<=fPolDeg; iPar++)
246 fSB->SetParameter(2+iPar, fBo->GetParameter(iPar+1));
247 hMinv->Fit(fSB,"0qR");
248 arrResults->AddAt(fSB,2);
249
250
251 // =======
252 // 3.
253 // =======
254 // Create the background function
255 TF1 *fB = new TF1("bgrdOnly_fkt",BgndFun,fFitMin,fFitMax,nBgnd);
256 fB->FixParameter(0,fPolDeg);
257 for (Int_t iDeg=0; iDeg<=fPolDeg; iDeg++) {
258 fB->SetParameter(1+iDeg,fSB->GetParameter(2+iDeg));
259 fB->SetParError(1+iDeg,fSB->GetParError(2+iDeg));
260 }
261 // create background histogram from background function
262 hBgrd->Eval(fB);
263 hBgrd->Fit(fB,"0qR");
264 // calculate the integral and integral error from fit function
265 Double_t intc = fB->Integral(intAreaEdgeLo, intAreaEdgeHi) * norm;
266 Double_t inte = fB->IntegralError(intAreaEdgeLo, intAreaEdgeHi) * norm;
267 arrResults->AddAt(fB,3);
268
269 // Fill the background spectrum erros. Use the error from the fit function for the background fB
270 for (Int_t iBin = 0; iBin <= hBgrd->GetNbinsX(); iBin++) {
271 Double_t x = hBgrd->GetBinCenter(iBin);
272 if ((x >= fFitMin) && (x <= fFitMax)) {
273 Double_t binte = inte / TMath::Sqrt((binIntHi-binIntLo)+1);
274 hBgrd->SetBinError(iBin,binte);
275 }
276 }
277 arrResults->AddAt(hBgrd,4);
278
279 // =======
280 // 4.
281 // =======
282 // Subtract the background
283 hSigF->Add(hMinv,hBgrd,1.0,-1.0);
284 for (Int_t iBin = 0; iBin <= hSigF->GetNbinsX(); iBin++) {
285 Double_t x = hSigF->GetBinCenter(iBin);
286 if ((x < fFitMin) || (x > fFitMax)) {
287 hSigF->SetBinContent(iBin,0.0);
288 hSigF->SetBinError(iBin,0.0);
289 }
290 }
291 hSigF->SetNameTitle("peak_only","");
292 arrResults->AddAt(hSigF,5);
293
294 // =======
295 // 5.
296 // =======
297 // Fit the background-subtracted spectrum
298 TF1 *fS = new TF1("peakOnly_fit",PeakFunCB,fFitMin,fFitMax,5);
299 fS->SetParameters(-.05,1,kJPMass,.003,700);
300 fS->SetParNames("alpha","n","meanx","sigma","N");
301 hSigF->Fit(fS,"0qR");
302 arrResults->AddAt(fS,6);
303
304
305 // connect data members
306 fFuncSignal = (TF1*) arrResults->At(6)->Clone();
307 fFuncBackground = (TF1*) arrResults->At(3)->Clone();
308 fFuncSigBack = (TF1*) arrResults->At(2)->Clone();
309 fHistSignal = (TH1F*)arrResults->At(5)->Clone();
310 fHistBackground = (TH1F*)arrResults->At(4)->Clone();
311
312 // fit results
313 Double_t chi2 = fSB->GetChisquare();
314 Int_t ndf = fSB->GetNDF();
315 fChi2Dof = chi2/ndf;
316
317 // signal + signal error
318 fValues(0) = hSigF->IntegralAndError(binIntLo, binIntHi, fErrors(0));
319 fValues(1) = intc; // background
320 fErrors(1) = inte; // background error
321 // S/B (2) and significance (3)
322 SetSignificanceAndSOB();
323 fValues(4) = fS->GetParameter(2); // mass
324 fErrors(4) = fS->GetParError(2); // mass error
325 fValues(5) = fS->GetParameter(3); // mass wdth
326 fErrors(5) = fS->GetParError(3); // mass wdth error
327
328
329 delete arrResults;
330
331}
332//______________________________________________________________________________
333Double_t AliDielectronSignalFunc::BgndFun(const Double_t *x, const Double_t *par) {
334 // parameters
335 // [0]: degree of polynomial
336 // [1]: constant polynomial coefficient
337 // [2]..: higher polynomial coefficients
338
339 Int_t deg = ((Int_t) par[0]);
340
341 Double_t f = 0.0;
342 Double_t yy = 1.0;
343 Double_t xx = x[0];
344
345 for (Int_t i = 0; i <= deg; i++) {
346 f += par[i+1] * yy;
347 yy *= xx;
348 }
349
350
351 return f;
352}
353//______________________________________________________________________________
354Double_t AliDielectronSignalFunc::PeakFun(const Double_t *x, const Double_t *par) {
355 // Fit MC signal shape
356 // parameters
357 // [0]: scale for simpeak
358
359 Double_t xx = x[0];
360
361 TH1F *hPeak = fgHistSimPM;
362 if (!hPeak) {
363 printf("F-AliDielectronSignalFunc::PeakFun: No histogram for peak fit defined!\n");
364 }
365
366 Int_t idx = hPeak->FindBin(xx);
367 if ((idx <= 1) ||
368 (idx >= hPeak->GetNbinsX())) {
369 return 0.0;
370 }
371
372 return (par[0] * hPeak->GetBinContent(idx));
373
374}
375
376//______________________________________________________________________________
377Double_t AliDielectronSignalFunc::MinvFun(const Double_t *x, const Double_t *par) {
378 // parameters
379 // [0]: degree of polynomial -> [0] for BgndFun
380 // [1]: scale for simpeak -> [0] for PeakFun
381 // [2]: constant polynomial coefficient -> [1] for BgndFun
382 // [3]..: higher polynomial coefficients -> [2].. for BgndFun
383
384 Int_t deg = ((Int_t) par[0]);
385 Double_t parPK[25], parBG[25];
386
387 parBG[0] = par[0]; // degree of polynomial
388
389 parPK[0] = par[1]; // MC minv scale
390 for (Int_t i = 0; i <= deg; i++) parBG[i+1] = par[i+2]; // polynomial coefficients
391
392 Double_t peak = PeakFun(x,parPK);
393 Double_t bgnd = BgndFun(x,parBG);
394 Double_t f = peak + bgnd;
395
396 return f;
397}
398
399//______________________________________________________________________________
400Double_t AliDielectronSignalFunc::PeakFunCB(const Double_t *x, const Double_t *par) {
401 // Crystal Ball function fit
402
403 Double_t alpha = par[0];
404 Double_t n = par[1];
405 Double_t meanx = par[2];
406 Double_t sigma = par[3];
407 Double_t nn = par[4];
408
409 Double_t a = TMath::Power((n/TMath::Abs(alpha)), n) * TMath::Exp(-.5*alpha*alpha);
410 Double_t b = n/TMath::Abs(alpha) - TMath::Abs(alpha);
411
412 Double_t arg = (x[0] - meanx)/sigma;
413 Double_t fitval = 0;
414
415 if (arg > -1.*alpha) {
416 fitval = nn * TMath::Exp(-.5*arg*arg);
417 } else {
418 fitval = nn * a * TMath::Power((b-arg), (-1*n));
419 }
420
421 return fitval;
422}
423
424
425//______________________________________________
426void AliDielectronSignalFunc::ProcessFit(TObjArray * const arrhist) {
427 //
428 // Fit the +- invariant mass distribution only
429 // Here we assume that the combined fit function is a sum of the signal and background functions
430 // and that the signal function is always the first term of this sum
431 //
432
433 fHistDataPM = (TH1F*)(arrhist->At(1))->Clone("histPM"); // +- SE
434 fHistDataPM->Sumw2();
435 if(fRebin>1)
436 fHistDataPM->Rebin(fRebin);
437
438 fHistSignal = new TH1F("HistSignal", "Fit substracted signal",
439 fHistDataPM->GetXaxis()->GetNbins(),
440 fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
441 fHistBackground = new TH1F("HistBackground", "Fit contribution",
442 fHistDataPM->GetXaxis()->GetNbins(),
443 fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
444
445 // the starting parameters of the fit function and their limits can be tuned
446 // by the user in its macro
447 fHistDataPM->Fit(fFuncSigBack, fFitOpt.Data(), "", fFitMin, fFitMax);
448 TFitResultPtr pmFitPtr = fHistDataPM->Fit(fFuncSigBack, fFitOpt.Data(), "", fFitMin, fFitMax);
449 //TFitResult *pmFitResult = pmFitPtr.Get(); // used only with TF1Helper
450 fFuncSignal->SetParameters(fFuncSigBack->GetParameters());
451 fFuncBackground->SetParameters(fFuncSigBack->GetParameters()+fFuncSignal->GetNpar());
452
453 for(Int_t iBin=1; iBin<=fHistDataPM->GetXaxis()->GetNbins(); iBin++) {
454 Double_t m = fHistDataPM->GetBinCenter(iBin);
455 Double_t pm = fHistDataPM->GetBinContent(iBin);
456 Double_t epm = fHistDataPM->GetBinError(iBin);
457 Double_t bknd = fFuncBackground->Eval(m);
458 Double_t ebknd = 0;
459 for(Int_t iPar=fFuncSignal->GetNpar(); iPar<fFuncSigBack->GetNpar(); iPar++) {
460/* TF1Helper problem on alien compilation
461 for(Int_t jPar=iPar; jPar<fFuncSigBack->GetNpar(); jPar++) {
462 TF1 gradientIpar("gradientIpar",
463 ROOT::TF1Helper::TGradientParFunction(iPar-fFuncSignal->GetNpar(),fFuncBackground),0,0,0);
464 TF1 gradientJpar("gradientJpar",
465 ROOT::TF1Helper::TGradientParFunction(jPar-fFuncSignal->GetNpar(),fFuncBackground),0,0,0);
466 ebknd += pmFitResult->CovMatrix(iPar,jPar)*
467 gradientIpar.Eval(m)*gradientJpar.Eval(m)*
468 (iPar==jPar ? 1.0 : 2.0);
469 }
470*/
471 }
472 Double_t signal = pm-bknd;
473 Double_t error = TMath::Sqrt(epm*epm+ebknd);
474 fHistSignal->SetBinContent(iBin, signal);
475 fHistSignal->SetBinError(iBin, error);
476 fHistBackground->SetBinContent(iBin, bknd);
477 fHistBackground->SetBinError(iBin, TMath::Sqrt(ebknd));
478 }
479
480 if(fUseIntegral) {
481 // signal
482 fValues(0) = fFuncSignal->Integral(fIntMin, fIntMax)/fHistDataPM->GetBinWidth(1);
483 fErrors(0) = 0;
484 for(Int_t iPar=0; iPar<fFuncSignal->GetNpar(); iPar++) {
485/* TF1Helper problem on alien compilation
486 for(Int_t jPar=iPar; jPar<fFuncSignal->GetNpar(); jPar++) {
487 TF1 gradientIpar("gradientIpar",
488 ROOT::TF1Helper::TGradientParFunction(iPar,fFuncSignal),0,0,0);
489 TF1 gradientJpar("gradientJpar",
490 ROOT::TF1Helper::TGradientParFunction(jPar,fFuncSignal),0,0,0);
491 fErrors(0) += pmFitResult->CovMatrix(iPar,jPar)*
492 gradientIpar.Integral(fIntMin,fIntMax)*gradientJpar.Integral(fIntMin,fIntMax)*
493 (iPar==jPar ? 1.0 : 2.0);
494 }
495*/
496 }
497 // background
498 fValues(1) = fFuncBackground->Integral(fIntMin, fIntMax)/fHistDataPM->GetBinWidth(1);
499 fErrors(1) = 0;
500 for(Int_t iPar=fFuncSignal->GetNpar(); iPar<fFuncSigBack->GetNpar(); iPar++) {
501/* TF1Helper problem on alien compilation
502 for(Int_t jPar=iPar; jPar<fFuncSigBack->GetNpar(); jPar++) {
503 TF1 gradientIpar("gradientIpar",
504 ROOT::TF1Helper::TGradientParFunction(iPar-fFuncSignal->GetNpar(),fFuncBackground),0,0,0);
505 TF1 gradientJpar("gradientJpar",
506 ROOT::TF1Helper::TGradientParFunction(jPar-fFuncSignal->GetNpar(),fFuncBackground),0,0,0);
507 fErrors(1) += pmFitResult->CovMatrix(iPar,jPar)*
508 gradientIpar.Integral(fIntMin, fIntMax)*gradientJpar.Integral(fIntMin, fIntMax)*
509 (iPar==jPar ? 1.0 : 2.0);
510 }
511*/
512 }
513 }
514 else {
515 // signal
516 fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
517 fHistSignal->FindBin(fIntMax), fErrors(0));
518 // background
519 fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
520 fHistBackground->FindBin(fIntMax),
521 fErrors(1));
522 }
523 // S/B and significance
524 SetSignificanceAndSOB();
525 fValues(4) = fFuncSigBack->GetParameter(fParMass);
526 fErrors(4) = fFuncSigBack->GetParError(fParMass);
527 fValues(5) = fFuncSigBack->GetParameter(fParMassWidth);
528 fErrors(5) = fFuncSigBack->GetParError(fParMassWidth);
529
530 fProcessed = kTRUE;
531}
532
533//______________________________________________
534void AliDielectronSignalFunc::ProcessLS(TObjArray * const arrhist) {
535 //
536 // Substract background using the like-sign spectrum
537 //
538 fHistDataPP = (TH1F*)(arrhist->At(0))->Clone("histPP"); // ++ SE
539 fHistDataPM = (TH1F*)(arrhist->At(1))->Clone("histPM"); // +- SE
540 fHistDataMM = (TH1F*)(arrhist->At(2))->Clone("histMM"); // -- SE
541 if (fRebin>1) {
542 fHistDataPP->Rebin(fRebin);
543 fHistDataPM->Rebin(fRebin);
544 fHistDataMM->Rebin(fRebin);
545 }
546 fHistDataPP->Sumw2();
547 fHistDataPM->Sumw2();
548 fHistDataMM->Sumw2();
549
550 fHistSignal = new TH1F("HistSignal", "Like-Sign substracted signal",
551 fHistDataPM->GetXaxis()->GetNbins(),
552 fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
553 fHistBackground = new TH1F("HistBackground", "Like-sign contribution",
554 fHistDataPM->GetXaxis()->GetNbins(),
555 fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax());
556
557 // fit the +- mass distribution
558 fHistDataPM->Fit(fFuncSigBack, fFitOpt.Data(), "", fFitMin, fFitMax);
559 fHistDataPM->Fit(fFuncSigBack, fFitOpt.Data(), "", fFitMin, fFitMax);
560 // declare the variables where the like-sign fit results will be stored
561// TFitResult *ppFitResult = 0x0;
562// TFitResult *mmFitResult = 0x0;
563 // fit the like sign background
564 TF1 *funcClonePP = (TF1*)fFuncBackground->Clone("funcClonePP");
565 TF1 *funcCloneMM = (TF1*)fFuncBackground->Clone("funcCloneMM");
566 fHistDataPP->Fit(funcClonePP, fFitOpt.Data(), "", fFitMin, fFitMax);
567 fHistDataPP->Fit(funcClonePP, fFitOpt.Data(), "", fFitMin, fFitMax);
568// TFitResultPtr ppFitPtr = fHistDataPP->Fit(funcClonePP, fFitOpt.Data(), "", fFitMin, fFitMax);
569// ppFitResult = ppFitPtr.Get();
570 fHistDataMM->Fit(funcCloneMM, fFitOpt.Data(), "", fFitMin, fFitMax);
571 fHistDataMM->Fit(funcCloneMM, fFitOpt.Data(), "", fFitMin, fFitMax);
572// TFitResultPtr mmFitPtr = fHistDataMM->Fit(funcCloneMM, fFitOpt.Data(), "", fFitMin, fFitMax);
573// mmFitResult = mmFitPtr.Get();
574
575 for(Int_t iBin=1; iBin<=fHistDataPM->GetXaxis()->GetNbins(); iBin++) {
576 Double_t m = fHistDataPM->GetBinCenter(iBin);
577 Double_t pm = fHistDataPM->GetBinContent(iBin);
578 Double_t pp = funcClonePP->Eval(m);
579 Double_t mm = funcCloneMM->Eval(m);
580 Double_t epm = fHistDataPM->GetBinError(iBin);
581 Double_t epp = 0;
582 for(Int_t iPar=0; iPar<funcClonePP->GetNpar(); iPar++) {
583/* TF1Helper problem on alien compilation
584 for(Int_t jPar=iPar; jPar<funcClonePP->GetNpar(); jPar++) {
585 TF1 gradientIpar("gradientIpar",
586 ROOT::TF1Helper::TGradientParFunction(iPar,funcClonePP),0,0,0);
587 TF1 gradientJpar("gradientJpar",
588 ROOT::TF1Helper::TGradientParFunction(jPar,funcClonePP),0,0,0);
589 epp += ppFitResult->CovMatrix(iPar,jPar)*
590 gradientIpar.Eval(m)*gradientJpar.Eval(m)*
591 (iPar==jPar ? 1.0 : 2.0);
592 }
593*/
594 }
595 Double_t emm = 0;
596 for(Int_t iPar=0; iPar<funcCloneMM->GetNpar(); iPar++) {
597/* TF1Helper problem on alien compilation
598 for(Int_t jPar=iPar; jPar<funcCloneMM->GetNpar(); jPar++) {
599 TF1 gradientIpar("gradientIpar",
600 ROOT::TF1Helper::TGradientParFunction(iPar,funcCloneMM),0,0,0);
601 TF1 gradientJpar("gradientJpar",
602 ROOT::TF1Helper::TGradientParFunction(jPar,funcCloneMM),0,0,0);
603 emm += mmFitResult->CovMatrix(iPar,jPar)*
604 gradientIpar.Eval(m)*gradientJpar.Eval(m)*
605 (iPar==jPar ? 1.0 : 2.0);
606 }
607*/
608 }
609
610 Double_t signal = pm-2.0*TMath::Sqrt(pp*mm);
611 Double_t background = 2.0*TMath::Sqrt(pp*mm);
612 // error propagation on the signal calculation above
613 Double_t esignal = TMath::Sqrt(epm*epm+(mm/pp)*epp+(pp/mm)*emm);
614 Double_t ebackground = TMath::Sqrt((mm/pp)*epp+(pp/mm)*emm);
615 fHistSignal->SetBinContent(iBin, signal);
616 fHistSignal->SetBinError(iBin, esignal);
617 fHistBackground->SetBinContent(iBin, background);
618 fHistBackground->SetBinError(iBin, ebackground);
619 }
620
621 // signal
622 fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin),
623 fHistSignal->FindBin(fIntMax), fErrors(0));
624 // background
625 fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin),
626 fHistBackground->FindBin(fIntMax),
627 fErrors(1));
628 // S/B and significance
629 SetSignificanceAndSOB();
630 fValues(4) = fFuncSigBack->GetParameter(fParMass);
631 fErrors(4) = fFuncSigBack->GetParError(fParMass);
632 fValues(5) = fFuncSigBack->GetParameter(fParMassWidth);
633 fErrors(5) = fFuncSigBack->GetParError(fParMassWidth);
634
635 fProcessed = kTRUE;
636}
637
638//______________________________________________
639void AliDielectronSignalFunc::ProcessEM(TObjArray * const arrhist) {
640 //
641 // Substract background with the event mixing technique
642 //
643 arrhist->GetEntries(); // just to avoid the unused parameter warning
644 AliError("Event mixing for background substraction method not implemented!");
645}
646
647//______________________________________________
648void AliDielectronSignalFunc::SetFunctions(TF1 * const combined, TF1 * const sig, TF1 * const back,
649 Int_t parM, Int_t parMres)
650{
651 //
652 // Set the signal, background functions and combined fit function
653 // Note: The process method assumes that the first n parameters in the
654 // combined fit function correspond to the n parameters of the signal function
655 // and the n+1 to n+m parameters to the m parameters of the background function!!!
656
657 if (!sig||!back||!combined) {
658 AliError("Both, signal and background function need to be set!");
659 return;
660 }
661 fFuncSignal=sig;
662 fFuncBackground=back;
663 fFuncSigBack=combined;
664 fParMass=parM;
665 fParMassWidth=parMres;
666}
667
668//______________________________________________
669void AliDielectronSignalFunc::SetDefaults(Int_t type)
670{
671 //
672 // Setup some default functions:
673 // type = 0: gaus signal + linear background in 2.5 - 4 GeV inv. mass
674 // type = 1: gaus signal + exponential background in 2.5 - 4 GeV inv. mass
675 // type = 2: half gaussian, half exponential signal function
676 // type = 3: Crystal-Ball function
677 // type = 4: Crystal-Ball signal + exponential background
678 //
679
680 if (type==0){
681 fFuncSignal=new TF1("DieleSignal","gaus",2.5,4);
682 fFuncBackground=new TF1("DieleBackground","pol1",2.5,4);
683 fFuncSigBack=new TF1("DieleCombined","gaus+pol1(3)",2.5,4);
684
685 fFuncSigBack->SetParameters(1,3.1,.05,2.5,1);
686 fFuncSigBack->SetParLimits(0,0,10000000);
687 fFuncSigBack->SetParLimits(1,3.05,3.15);
688 fFuncSigBack->SetParLimits(2,.02,.1);
689 }
690 else if (type==1){
691 fFuncSignal=new TF1("DieleSignal","gaus",2.5,4);
692 fFuncBackground=new TF1("DieleBackground","[0]*exp(-(x-[1])/[2])",2.5,4);
693 fFuncSigBack=new TF1("DieleCombined","gaus+[3]*exp(-(x-[4])/[5])",2.5,4);
694
695 fFuncSigBack->SetParameters(1,3.1,.05,1,2.5,1);
696 fFuncSigBack->SetParLimits(0,0,10000000);
697 fFuncSigBack->SetParLimits(1,3.05,3.15);
698 fFuncSigBack->SetParLimits(2,.02,.1);
699 }
700 else if (type==2){
701 // half gaussian, half exponential signal function
702 // exponential background
703 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);
704 fFuncBackground = new TF1("DieleBackground","[0]*exp(-(x-[1])/[2])+[3]",2.5,4);
705 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);
706 fFuncSigBack->SetParameters(1.,3.1,.05,.1,1,2.5,1,0);
707
708 fFuncSigBack->SetParLimits(0,0,10000000);
709 fFuncSigBack->SetParLimits(1,3.05,3.15);
710 fFuncSigBack->SetParLimits(2,.02,.1);
711 fFuncSigBack->FixParameter(6,2.5);
712 fFuncSigBack->FixParameter(7,0);
713 }
714}
715
716
717//______________________________________________
718void AliDielectronSignalFunc::Draw(const Option_t* option)
719{
720 //
721 // Draw the fitted function
722 //
723
724 TString drawOpt(option);
725 drawOpt.ToLower();
726
727 Bool_t optStat=drawOpt.Contains("stat");
728
729 fFuncSigBack->SetNpx(200);
730 fFuncSigBack->SetRange(fIntMin,fIntMax);
731 fFuncBackground->SetNpx(200);
732 fFuncBackground->SetRange(fIntMin,fIntMax);
733
734 TGraph *grSig=new TGraph(fFuncSigBack);
735 grSig->SetFillColor(kGreen);
736 grSig->SetFillStyle(3001);
737
738 TGraph *grBack=new TGraph(fFuncBackground);
739 grBack->SetFillColor(kRed);
740 grBack->SetFillStyle(3001);
741
742 grSig->SetPoint(0,grBack->GetX()[0],grBack->GetY()[0]);
743 grSig->SetPoint(grSig->GetN()-1,grBack->GetX()[grBack->GetN()-1],grBack->GetY()[grBack->GetN()-1]);
744
745 grBack->SetPoint(0,grBack->GetX()[0],0.);
746 grBack->SetPoint(grBack->GetN()-1,grBack->GetX()[grBack->GetN()-1],0.);
747
748 fFuncSigBack->SetRange(fFitMin,fFitMax);
749 fFuncBackground->SetRange(fFitMin,fFitMax);
750
751 if (!drawOpt.Contains("same")){
752 if (fHistDataPM){
753 fHistDataPM->Draw();
754 grSig->Draw("f");
755 } else {
756 grSig->Draw("af");
757 }
758 } else {
759 grSig->Draw("f");
760 }
761 if(fMethod==kFitted || fMethod==kFittedMC) grBack->Draw("f");
762 fFuncSigBack->Draw("same");
763 fFuncSigBack->SetLineWidth(2);
764 if(fMethod==kLikeSign) {
765 fHistDataPP->SetLineWidth(2);
766 fHistDataPP->SetLineColor(6);
767 fHistDataPP->Draw("same");
768 fHistDataMM->SetLineWidth(2);
769 fHistDataMM->SetLineColor(8);
770 fHistDataMM->Draw("same");
771 }
772
773 if(fMethod==kFitted || fMethod==kFittedMC)
774 fFuncBackground->Draw("same");
775
776 if (optStat) DrawStats();
777
778}