]>
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 | /* | |
572b0139 | 21 | |
37a6f270 | 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 | |
572b0139 | 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> | |
bc75eeb5 | 114 | #include <TList.h> |
115 | #include <TFitResult.h> | |
cafdf10b | 116 | //#include <../hist/hist/src/TF1Helper.h> |
572b0139 | 117 | |
118 | #include <AliLog.h> | |
119 | ||
120 | #include "AliDielectronSignalFunc.h" | |
121 | ||
122 | ClassImp(AliDielectronSignalFunc) | |
123 | ||
37a6f270 | 124 | //TH1F* AliDielectronSignalFunc::fgHistSimPM=0x0; |
e43f22fb | 125 | TF1* AliDielectronSignalFunc::fFuncSignal=0x0; |
126 | TF1* AliDielectronSignalFunc::fFuncBackground=0x0; | |
127 | Int_t AliDielectronSignalFunc::fNparPeak=0; | |
128 | Int_t AliDielectronSignalFunc::fNparBgnd=0; | |
129 | ||
5720c765 | 130 | |
572b0139 | 131 | AliDielectronSignalFunc::AliDielectronSignalFunc() : |
37a6f270 | 132 | AliDielectronSignalExt(), |
3505bfad | 133 | fFuncSigBack(0x0), |
134 | fParMass(1), | |
135 | fParMassWidth(2), | |
136 | fFitOpt("SMNQE"), | |
5720c765 | 137 | fUseIntegral(kFALSE), |
e43f22fb | 138 | fDof(0), |
5720c765 | 139 | fChi2Dof(0.0) |
572b0139 | 140 | { |
141 | // | |
142 | // Default Constructor | |
143 | // | |
144 | ||
145 | } | |
146 | ||
147 | //______________________________________________ | |
148 | AliDielectronSignalFunc::AliDielectronSignalFunc(const char* name, const char* title) : | |
37a6f270 | 149 | AliDielectronSignalExt(name, title), |
3505bfad | 150 | fFuncSigBack(0x0), |
151 | fParMass(1), | |
152 | fParMassWidth(2), | |
153 | fFitOpt("SMNQE"), | |
5720c765 | 154 | fUseIntegral(kFALSE), |
e43f22fb | 155 | fDof(0), |
5720c765 | 156 | fChi2Dof(0.0) |
572b0139 | 157 | { |
158 | // | |
159 | // Named Constructor | |
160 | // | |
161 | } | |
162 | ||
163 | //______________________________________________ | |
164 | AliDielectronSignalFunc::~AliDielectronSignalFunc() | |
165 | { | |
166 | // | |
167 | // Default Destructor | |
168 | // | |
bc75eeb5 | 169 | if(fFuncSigBack) delete fFuncSigBack; |
572b0139 | 170 | } |
171 | ||
172 | ||
173 | //______________________________________________ | |
bc75eeb5 | 174 | void AliDielectronSignalFunc::Process(TObjArray * const arrhist) |
572b0139 | 175 | { |
176 | // | |
bc75eeb5 | 177 | // Fit the invariant mass histograms and retrieve the signal and background |
572b0139 | 178 | // |
bc75eeb5 | 179 | switch(fMethod) { |
180 | case kFitted : | |
181 | ProcessFit(arrhist); | |
182 | break; | |
37a6f270 | 183 | |
184 | case kLikeSignFit : | |
185 | ProcessFitLS(arrhist); | |
186 | break; | |
187 | ||
188 | case kEventMixingFit : | |
189 | ProcessFitEM(arrhist); | |
190 | break; | |
191 | ||
bc75eeb5 | 192 | case kLikeSign : |
37a6f270 | 193 | case kLikeSignArithm : |
194 | case kLikeSignRcorr: | |
195 | case kLikeSignArithmRcorr: | |
196 | ProcessLS(arrhist); // process like-sign subtraction method | |
bc75eeb5 | 197 | break; |
37a6f270 | 198 | |
bc75eeb5 | 199 | case kEventMixing : |
37a6f270 | 200 | ProcessEM(arrhist); // process event mixing method |
bc75eeb5 | 201 | break; |
37a6f270 | 202 | |
203 | case kRotation: | |
204 | ProcessRotation(arrhist); | |
205 | break; | |
206 | ||
bc75eeb5 | 207 | default : |
208 | AliError("Background substraction method not known!"); | |
572b0139 | 209 | } |
bc75eeb5 | 210 | } |
5720c765 | 211 | //______________________________________________________________________________ |
e43f22fb | 212 | Double_t AliDielectronSignalFunc::PeakFunMC(const Double_t *x, const Double_t *par) { |
5720c765 | 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) { | |
ac390e40 | 221 | printf("E-AliDielectronSignalFunc::PeakFun: No histogram for peak fit defined!\n"); |
222 | return 0.0; | |
5720c765 | 223 | } |
224 | ||
225 | Int_t idx = hPeak->FindBin(xx); | |
226 | if ((idx <= 1) || | |
227 | (idx >= hPeak->GetNbinsX())) { | |
228 | return 0.0; | |
229 | } | |
230 | ||
37a6f270 | 231 | return (par[0] * hPeak->GetBinContent(idx)); |
5720c765 | 232 | |
5720c765 | 233 | } |
234 | ||
235 | //______________________________________________________________________________ | |
236 | Double_t AliDielectronSignalFunc::PeakFunCB(const Double_t *x, const Double_t *par) { | |
e43f22fb | 237 | // Crystal Ball fit function |
5720c765 | 238 | |
37a6f270 | 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]; | |
5720c765 | 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 | ||
e43f22fb | 260 | //______________________________________________________________________________ |
261 | Double_t AliDielectronSignalFunc::PeakFunGaus(const Double_t *x, const Double_t *par) { | |
262 | // Gaussian fit function | |
37a6f270 | 263 | //printf("fNparBgrd %d \n",fNparBgnd); |
264 | Double_t n = par[0]; | |
265 | Double_t mean = par[1]; | |
266 | Double_t sigma = par[2]; | |
e43f22fb | 267 | Double_t xx = x[0]; |
268 | ||
37a6f270 | 269 | return ( n*TMath::Exp(-0.5*TMath::Power((xx-mean)/sigma,2)) ); |
e43f22fb | 270 | } |
572b0139 | 271 | |
bc75eeb5 | 272 | //______________________________________________ |
273 | void 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 | // | |
3505bfad | 279 | |
bc75eeb5 | 280 | fHistDataPM = (TH1F*)(arrhist->At(1))->Clone("histPM"); // +- SE |
281 | fHistDataPM->Sumw2(); | |
282 | if(fRebin>1) | |
283 | fHistDataPM->Rebin(fRebin); | |
572b0139 | 284 | |
5720c765 | 285 | fHistSignal = new TH1F("HistSignal", "Fit substracted signal", |
3505bfad | 286 | fHistDataPM->GetXaxis()->GetNbins(), |
287 | fHistDataPM->GetXaxis()->GetXmin(), fHistDataPM->GetXaxis()->GetXmax()); | |
5720c765 | 288 | fHistBackground = new TH1F("HistBackground", "Fit contribution", |
3505bfad | 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 | |
bc75eeb5 | 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); | |
cafdf10b | 296 | //TFitResult *pmFitResult = pmFitPtr.Get(); // used only with TF1Helper |
37a6f270 | 297 | //fFuncBackground->SetParameters(fFuncSigBack->GetParameters()); |
298 | fFuncSignal->SetParameters(fFuncSigBack->GetParameters()); | |
299 | fFuncBackground->SetParameters(fFuncSigBack->GetParameters()+fFuncSignal->GetNpar()); | |
300 | ||
37a6f270 | 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 | ||
bc75eeb5 | 311 | for(Int_t iBin=1; iBin<=fHistDataPM->GetXaxis()->GetNbins(); iBin++) { |
37a6f270 | 312 | // Double_t m = fHistDataPM->GetBinCenter(iBin); |
bc75eeb5 | 313 | Double_t pm = fHistDataPM->GetBinContent(iBin); |
314 | Double_t epm = fHistDataPM->GetBinError(iBin); | |
37a6f270 | 315 | Double_t bknd = fHistBackground->GetBinContent(iBin); |
316 | Double_t ebknd = fHistBackground->GetBinError(iBin); | |
bc75eeb5 | 317 | for(Int_t iPar=fFuncSignal->GetNpar(); iPar<fFuncSigBack->GetNpar(); iPar++) { |
37a6f270 | 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 | */ | |
bc75eeb5 | 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); | |
bc75eeb5 | 334 | } |
927480a1 | 335 | |
bc75eeb5 | 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++) { | |
cafdf10b | 341 | /* TF1Helper problem on alien compilation |
bc75eeb5 | 342 | for(Int_t jPar=iPar; jPar<fFuncSignal->GetNpar(); jPar++) { |
343 | TF1 gradientIpar("gradientIpar", | |
3505bfad | 344 | ROOT::TF1Helper::TGradientParFunction(iPar,fFuncSignal),0,0,0); |
bc75eeb5 | 345 | TF1 gradientJpar("gradientJpar", |
3505bfad | 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); | |
bc75eeb5 | 350 | } |
cafdf10b | 351 | */ |
3505bfad | 352 | } |
bc75eeb5 | 353 | // background |
354 | fValues(1) = fFuncBackground->Integral(fIntMin, fIntMax)/fHistDataPM->GetBinWidth(1); | |
37a6f270 | 355 | fErrors(1) = fFuncBackground->IntegralError(fIntMin, fIntMax)/fHistDataPM->GetBinWidth(1); |
bc75eeb5 | 356 | for(Int_t iPar=fFuncSignal->GetNpar(); iPar<fFuncSigBack->GetNpar(); iPar++) { |
cafdf10b | 357 | /* TF1Helper problem on alien compilation |
bc75eeb5 | 358 | for(Int_t jPar=iPar; jPar<fFuncSigBack->GetNpar(); jPar++) { |
359 | TF1 gradientIpar("gradientIpar", | |
3505bfad | 360 | ROOT::TF1Helper::TGradientParFunction(iPar-fFuncSignal->GetNpar(),fFuncBackground),0,0,0); |
bc75eeb5 | 361 | TF1 gradientJpar("gradientJpar", |
3505bfad | 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); | |
bc75eeb5 | 366 | } |
cafdf10b | 367 | */ |
bc75eeb5 | 368 | } |
572b0139 | 369 | } |
bc75eeb5 | 370 | else { |
371 | // signal | |
372 | fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin), | |
3505bfad | 373 | fHistSignal->FindBin(fIntMax), fErrors(0)); |
bc75eeb5 | 374 | // background |
375 | fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin), | |
37a6f270 | 376 | fHistBackground->FindBin(fIntMax), |
377 | fErrors(1)); | |
8df8e382 | 378 | } |
bc75eeb5 | 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); | |
3505bfad | 385 | |
bc75eeb5 | 386 | fProcessed = kTRUE; |
e43f22fb | 387 | |
388 | fHistBackground->GetListOfFunctions()->Add(fFuncBackground); | |
389 | ||
572b0139 | 390 | } |
391 | ||
392 | //______________________________________________ | |
37a6f270 | 393 | void AliDielectronSignalFunc::ProcessFitLS(TObjArray * const arrhist) { |
572b0139 | 394 | // |
bc75eeb5 | 395 | // Substract background using the like-sign spectrum |
572b0139 | 396 | // |
bc75eeb5 | 397 | fHistDataPP = (TH1F*)(arrhist->At(0))->Clone("histPP"); // ++ SE |
398 | fHistDataPM = (TH1F*)(arrhist->At(1))->Clone("histPM"); // +- SE | |
3505bfad | 399 | fHistDataMM = (TH1F*)(arrhist->At(2))->Clone("histMM"); // -- SE |
400 | if (fRebin>1) { | |
bc75eeb5 | 401 | fHistDataPP->Rebin(fRebin); |
402 | fHistDataPM->Rebin(fRebin); | |
403 | fHistDataMM->Rebin(fRebin); | |
3505bfad | 404 | } |
bc75eeb5 | 405 | fHistDataPP->Sumw2(); |
406 | fHistDataPM->Sumw2(); | |
407 | fHistDataMM->Sumw2(); | |
3505bfad | 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 | ||
bc75eeb5 | 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 | |
45b2b1b8 | 420 | // TFitResult *ppFitResult = 0x0; |
421 | // TFitResult *mmFitResult = 0x0; | |
bc75eeb5 | 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); | |
45b2b1b8 | 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); | |
bc75eeb5 | 430 | fHistDataMM->Fit(funcCloneMM, fFitOpt.Data(), "", fFitMin, fFitMax); |
45b2b1b8 | 431 | // TFitResultPtr mmFitPtr = fHistDataMM->Fit(funcCloneMM, fFitOpt.Data(), "", fFitMin, fFitMax); |
432 | // mmFitResult = mmFitPtr.Get(); | |
3505bfad | 433 | |
bc75eeb5 | 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++) { | |
cafdf10b | 442 | /* TF1Helper problem on alien compilation |
bc75eeb5 | 443 | for(Int_t jPar=iPar; jPar<funcClonePP->GetNpar(); jPar++) { |
444 | TF1 gradientIpar("gradientIpar", | |
3505bfad | 445 | ROOT::TF1Helper::TGradientParFunction(iPar,funcClonePP),0,0,0); |
bc75eeb5 | 446 | TF1 gradientJpar("gradientJpar", |
3505bfad | 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); | |
bc75eeb5 | 451 | } |
cafdf10b | 452 | */ |
bc75eeb5 | 453 | } |
454 | Double_t emm = 0; | |
455 | for(Int_t iPar=0; iPar<funcCloneMM->GetNpar(); iPar++) { | |
cafdf10b | 456 | /* TF1Helper problem on alien compilation |
bc75eeb5 | 457 | for(Int_t jPar=iPar; jPar<funcCloneMM->GetNpar(); jPar++) { |
458 | TF1 gradientIpar("gradientIpar", | |
3505bfad | 459 | ROOT::TF1Helper::TGradientParFunction(iPar,funcCloneMM),0,0,0); |
bc75eeb5 | 460 | TF1 gradientJpar("gradientJpar", |
3505bfad | 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); | |
bc75eeb5 | 465 | } |
cafdf10b | 466 | */ |
bc75eeb5 | 467 | } |
3505bfad | 468 | |
bc75eeb5 | 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 | } | |
2a14a7b1 | 479 | |
bc75eeb5 | 480 | // signal |
481 | fValues(0) = fHistSignal->IntegralAndError(fHistSignal->FindBin(fIntMin), | |
3505bfad | 482 | fHistSignal->FindBin(fIntMax), fErrors(0)); |
bc75eeb5 | 483 | // background |
484 | fValues(1) = fHistBackground->IntegralAndError(fHistBackground->FindBin(fIntMin), | |
3505bfad | 485 | fHistBackground->FindBin(fIntMax), |
486 | fErrors(1)); | |
bc75eeb5 | 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); | |
3505bfad | 493 | |
bc75eeb5 | 494 | fProcessed = kTRUE; |
572b0139 | 495 | } |
bc75eeb5 | 496 | |
497 | //______________________________________________ | |
37a6f270 | 498 | void AliDielectronSignalFunc::ProcessFitEM(TObjArray * const arrhist) { |
bc75eeb5 | 499 | // |
500 | // Substract background with the event mixing technique | |
501 | // | |
3505bfad | 502 | arrhist->GetEntries(); // just to avoid the unused parameter warning |
503 | AliError("Event mixing for background substraction method not implemented!"); | |
504 | } | |
bc75eeb5 | 505 | |
572b0139 | 506 | //______________________________________________ |
bc75eeb5 | 507 | void AliDielectronSignalFunc::SetFunctions(TF1 * const combined, TF1 * const sig, TF1 * const back, |
572b0139 | 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!!! | |
bc75eeb5 | 515 | |
572b0139 | 516 | if (!sig||!back||!combined) { |
517 | AliError("Both, signal and background function need to be set!"); | |
61d106d3 | 518 | return; |
572b0139 | 519 | } |
bc75eeb5 | 520 | fFuncSignal=sig; |
521 | fFuncBackground=back; | |
522 | fFuncSigBack=combined; | |
523 | fParMass=parM; | |
524 | fParMassWidth=parMres; | |
e43f22fb | 525 | |
572b0139 | 526 | } |
527 | ||
528 | //______________________________________________ | |
529 | void 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 | // | |
572b0139 | 539 | |
540 | if (type==0){ | |
bc75eeb5 | 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); | |
572b0139 | 544 | |
bc75eeb5 | 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); | |
3505bfad | 549 | } |
bc75eeb5 | 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); | |
572b0139 | 554 | |
bc75eeb5 | 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); | |
3505bfad | 559 | } |
bc75eeb5 | 560 | else if (type==2){ |
572b0139 | 561 | // half gaussian, half exponential signal function |
562 | // exponential background | |
bc75eeb5 | 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); | |
572b0139 | 567 | |
bc75eeb5 | 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); | |
572b0139 | 573 | } |
574 | } | |
575 | ||
576 | ||
577 | //______________________________________________ | |
578 | void AliDielectronSignalFunc::Draw(const Option_t* option) | |
579 | { | |
580 | // | |
581 | // Draw the fitted function | |
582 | // | |
bc75eeb5 | 583 | |
572b0139 | 584 | TString drawOpt(option); |
585 | drawOpt.ToLower(); | |
3505bfad | 586 | |
572b0139 | 587 | Bool_t optStat=drawOpt.Contains("stat"); |
588 | ||
bc75eeb5 | 589 | fFuncSigBack->SetNpx(200); |
590 | fFuncSigBack->SetRange(fIntMin,fIntMax); | |
591 | fFuncBackground->SetNpx(200); | |
592 | fFuncBackground->SetRange(fIntMin,fIntMax); | |
572b0139 | 593 | |
bc75eeb5 | 594 | TGraph *grSig=new TGraph(fFuncSigBack); |
572b0139 | 595 | grSig->SetFillColor(kGreen); |
596 | grSig->SetFillStyle(3001); | |
3505bfad | 597 | |
bc75eeb5 | 598 | TGraph *grBack=new TGraph(fFuncBackground); |
572b0139 | 599 | grBack->SetFillColor(kRed); |
600 | grBack->SetFillStyle(3001); | |
3505bfad | 601 | |
572b0139 | 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 | ||
bc75eeb5 | 608 | fFuncSigBack->SetRange(fFitMin,fFitMax); |
609 | fFuncBackground->SetRange(fFitMin,fFitMax); | |
572b0139 | 610 | |
611 | if (!drawOpt.Contains("same")){ | |
bc75eeb5 | 612 | if (fHistDataPM){ |
613 | fHistDataPM->Draw(); | |
8df8e382 | 614 | grSig->Draw("f"); |
615 | } else { | |
616 | grSig->Draw("af"); | |
617 | } | |
572b0139 | 618 | } else { |
619 | grSig->Draw("f"); | |
620 | } | |
5720c765 | 621 | if(fMethod==kFitted || fMethod==kFittedMC) grBack->Draw("f"); |
bc75eeb5 | 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); | |
3505bfad | 629 | fHistDataMM->SetLineColor(8); |
bc75eeb5 | 630 | fHistDataMM->Draw("same"); |
631 | } | |
3505bfad | 632 | |
5720c765 | 633 | if(fMethod==kFitted || fMethod==kFittedMC) |
bc75eeb5 | 634 | fFuncBackground->Draw("same"); |
3505bfad | 635 | |
8df8e382 | 636 | if (optStat) DrawStats(); |
bc75eeb5 | 637 | |
572b0139 | 638 | } |
e43f22fb | 639 | |
640 | ||
641 | //______________________________________________________________________________ | |
642 | void 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 | ||
37a6f270 | 657 | fFuncSigBack = new TF1("BgndPeak",AliDielectronSignalFunc::PeakBgndFun, fFitMin,fFitMax, fNparPeak+fNparBgnd); |
e43f22fb | 658 | return; |
659 | } | |
660 | ||
661 | //______________________________________________________________________________ | |
662 | Double_t AliDielectronSignalFunc::PeakBgndFun(const Double_t *x, const Double_t *par) { | |
663 | // | |
664 | // merge peak and bgnd functions | |
665 | // | |
37a6f270 | 666 | return (fFuncSignal->EvalPar(x,par) + fFuncBackground->EvalPar(x,par+fNparPeak)); |
e43f22fb | 667 | } |
668 |