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