]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGDQ/dielectron/AliDielectronBtoJPSItoEleCDFfitHandler.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
o update dielectron package
[u/mrichter/AliRoot.git] / PWGDQ / dielectron / AliDielectronBtoJPSItoEleCDFfitHandler.cxx
CommitLineData
ba15fdfb 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#include <TVirtualFitter.h>
5720c765 16#include <TFitter.h>
17#include <TMinuit.h>
ba15fdfb 18#include <TStopwatch.h>
5720c765 19#include <TCanvas.h>
ba15fdfb 20#include "AliLog.h"
21#include "AliDielectronBtoJPSItoEleCDFfitHandler.h"
22#include "AliDielectronBtoJPSItoEleCDFfitFCN.h"
23
24//-------------------------------------------------------------------------
25// Class AliDielectronBtoJPSItoEleCDFfitHandler
26// Class to perform unbinned log-likelihood fit
27//
28// Origin: C. Di Giglio
29// Contact: carmelo.digiglio@ba.infn.it , giuseppe.bruno@ba.infn.it
30//-------------------------------------------------------------------------
31
32void CDFFunction(Int_t &npar, Double_t *gin, Double_t &f, Double_t *x, Int_t iflag);
33
34ClassImp(AliDielectronBtoJPSItoEleCDFfitHandler)
35
36 //______________________________________________________________________________
37void CDFFunction(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
38{
39 // This function is called by minuit
40 // The corresponding member method is called
41 // using SetObjectFit/GetObjectFit methods of TMinuit
42 AliDielectronBtoJPSItoEleCDFfitHandler* dummy = (AliDielectronBtoJPSItoEleCDFfitHandler *)TVirtualFitter::GetFitter()->GetObjectFit();
43 dummy->CdfFCN(npar, gin, f, par, iflag);
44}
45
46
47//_________________________________________________________________________________________________
48AliDielectronBtoJPSItoEleCDFfitHandler::AliDielectronBtoJPSItoEleCDFfitHandler():
5720c765 49 fIsParamFixed(45),
ba15fdfb 50 fPrintStatus(kFALSE),
51 fUp(0),
52 fX(0x0),
53 fM(0x0),
5720c765 54 fType(0x0),
55 fLikely(0x0),
56 fNcand(0),
57 fContPlot1(0x0),
58 fContPlot2(0x0),
59 fContPlot3(0x0),
60 fitter(0)
ba15fdfb 61{
62 //
63 // default constructor
64 //
65}
66//_________________________________________________________________________________________________
67AliDielectronBtoJPSItoEleCDFfitHandler::AliDielectronBtoJPSItoEleCDFfitHandler(Double_t* decaytime,
5720c765 68 Double_t* invariantmass, Int_t *type, Int_t ncand) :
69 fIsParamFixed(45),
ba15fdfb 70 fPrintStatus(kFALSE),
71 fUp(0),
72 fX(decaytime),
73 fM(invariantmass),
5720c765 74 fType(type),
ba15fdfb 75 fLikely(0x0),
5720c765 76 fNcand(ncand),
77 fContPlot1(0x0),
78 fContPlot2(0x0),
79 fContPlot3(0x0),
80 fitter(0)
ba15fdfb 81{
82 //
83 // constructor
84 //
85 AliInfo("\n+++\n+++ Minimization object AliDielectronBtoJPSItoEleCDFfitHandler created\n+++\n");
86 fLikely = new AliDielectronBtoJPSItoEleCDFfitFCN();
87 AliInfo("\n+++\n+++ CDF fit function object AliDielectronBtoJPSItoEleCDFfitFCN created\n+++\n");
88 AliInfo("Parameter 0 ----> fWeightRes");
89 AliInfo("Parameter 1 ----> fPos");
90 AliInfo("Parameter 2 ----> fNeg");
91 AliInfo("Parameter 3 ----> fSym");
92 AliInfo("Parameter 4 ----> fOneOvLamPlus");
93 AliInfo("Parameter 5 ----> fOneOvLamMinus");
94 AliInfo("Parameter 6 ----> fOneOvLamSym");
95 AliInfo("Parameter 7 ----> fB");
96 AliInfo("Parameter 8 ----> fFsig");
97 AliInfo("Parameter 9 ----> fMmean");
98 AliInfo("Parameter 10 ----> fNexp");
99 AliInfo("Parameter 11 ----> fSigma");
100 AliInfo("Parameter 12 ----> fAlpha");
101 AliInfo("Parameter 13 ----> fNorm");
102 AliInfo("Parameter 14 ----> fBkgNorm");
103 AliInfo("Parameter 15 ----> fBkgMean");
104 AliInfo("Parameter 16 ----> fBkgSlope");
105 AliInfo("Parameter 17 ----> fBkgConst");
5720c765 106 AliInfo("Parameter 18 ----> fNormGaus1 (FF)");
107 AliInfo("Parameter 19 ----> fNormGaus2 (FF)");
108 AliInfo("Parameter 20 ----> fMean1Res (FF)");
109 AliInfo("Parameter 21 ----> fsigma1Res (FF)");
110 AliInfo("Parameter 22 ----> fMean2Res (FF)");
111 AliInfo("Parameter 23 ----> fsigma2Res (FF)");
112 AliInfo("Parameter 24 ----> fAlfaRes (FF)");
113 AliInfo("Parameter 25 ----> fLambdaRes (FF)");
114 AliInfo("Parameter 26 ----> fNormResExp (FF)");
115 AliInfo("Parameter 27 ----> fNormGaus1 (FS)");
116 AliInfo("Parameter 28 ----> fNormGaus2 (FS)");
117 AliInfo("Parameter 29 ----> fMean1Res (FS)");
118 AliInfo("Parameter 30 ----> fsigma1Res (FS)");
119 AliInfo("Parameter 31 ----> fMean2Res (FS)");
120 AliInfo("Parameter 32 ----> fsigma2Res (FS)");
121 AliInfo("Parameter 33 ----> fAlfaRes (FS)");
122 AliInfo("Parameter 34 ----> fLambdaRes (FS)");
123 AliInfo("Parameter 35 ----> fNormResExp (FS)");
124 AliInfo("Parameter 36 ----> fNormGaus1 (SS)");
125 AliInfo("Parameter 37 ----> fNormGaus2 (SS)");
126 AliInfo("Parameter 38 ----> fMean1Res (SS)");
127 AliInfo("Parameter 39 ----> fsigma1Res (SS)");
128 AliInfo("Parameter 40 ----> fMean2Res (SS)");
129 AliInfo("Parameter 41 ----> fsigma2Res (SS)");
130 AliInfo("Parameter 42 ----> fAlfaRes (SS)");
131 AliInfo("Parameter 43 ----> fLambdaRes (SS)");
132 AliInfo("Parameter 44 ----> fNormResExp (SS)");
ba15fdfb 133
134 AliInfo(Form("\n+++\n+++ Number of candidates ---> %d\n+++\n ", ncand));
135}
136//___________________________________________________________________________
137AliDielectronBtoJPSItoEleCDFfitHandler& AliDielectronBtoJPSItoEleCDFfitHandler::operator=(const AliDielectronBtoJPSItoEleCDFfitHandler& c)
138{
139 //
140 // Assignment operator
141 //
142 if (this!=&c) {
143 fIsParamFixed = c.fIsParamFixed;
144 fPrintStatus = c.fPrintStatus;
145 fUp = c.fUp;
146 fX = c.fX;
147 fM = c.fM;
5720c765 148 fType = c.fType;
149 fLikely = c.fLikely;
ba15fdfb 150 fNcand = c.fNcand;
5720c765 151 fContPlot1 = c.fContPlot1;
152 fContPlot2 = c.fContPlot2;
153 fContPlot3 = c.fContPlot3;
ba15fdfb 154 }
155 return *this;
156}
157
158//_______________________________________________________________________________________
159AliDielectronBtoJPSItoEleCDFfitHandler::AliDielectronBtoJPSItoEleCDFfitHandler(const AliDielectronBtoJPSItoEleCDFfitHandler& c) :
160 TNamed(c),
161 fIsParamFixed(c.fIsParamFixed),
162 fPrintStatus(c.fPrintStatus),
163 fUp(c.fUp),
164 fX(c.fX),
165 fM(c.fM),
5720c765 166 fType(c.fType),
167 fLikely(c.fLikely),
168 fNcand(c.fNcand),
169 fContPlot1(c.fContPlot1),
170 fContPlot2(c.fContPlot2),
171 fContPlot3(c.fContPlot3),
172 fitter(c.fitter)
ba15fdfb 173{
174 //
175 // Copy Constructor
176 //
177}
178//_______________________________________________________________________________________
179AliDielectronBtoJPSItoEleCDFfitHandler::~AliDielectronBtoJPSItoEleCDFfitHandler()
180{
181 //
182 //destructor
183 //
184 delete fLikely;
185}
186//_______________________________________________________________________________________
5720c765 187Int_t AliDielectronBtoJPSItoEleCDFfitHandler::DoMinimization(Int_t step)
ba15fdfb 188{
189 //
190 // performs the minimization
191 //
5720c765 192 if(step == 0){
193 //fitter = TVirtualFitter::Fitter(this,20);
194 fitter = (TFitter*)TVirtualFitter::Fitter(this,45);
195 fitter->SetFCN(CDFFunction);
196 fitter->SetParameter(0,"fWeightRes",fParamStartValues[0], 1.e-08, 0., 1.e+06);
197 fitter->SetParameter(1,"fPos",fParamStartValues[1], 1.e-08, 0.,1.e+06);
198 fitter->SetParameter(2,"fNeg",fParamStartValues[2], 1.e-08, 0.,1.e+06);
199 fitter->SetParameter(3,"fSym",fParamStartValues[3], 1.e-08, 0.,1.e+06);
200 fitter->SetParameter(4,"fOneOvLamPlus",fParamStartValues[4], 1.e-10, 0.0000001, 5.e+01);
201 fitter->SetParameter(5,"fOneOvLamMinus",fParamStartValues[5], 1.e-10, 0.00000001, 5.e+01);
202 fitter->SetParameter(6,"fOneOvLamSym",fParamStartValues[6], 1.e-10, 0.00000001, 5.e+01);
203 fitter->SetParameter(7,"fB",fParamStartValues[7], 1.e-10, 0., 1.);
204 fitter->SetParameter(8,"fFsig",fParamStartValues[8], 1.e-10, 0., 1.);
205 fitter->SetParameter(9,"fMmean",fParamStartValues[9], 1.e-08, 0., 1.e+04);
206 fitter->SetParameter(10,"fNexp",fParamStartValues[10], 1.e-08, 0., 1.e+02);
207 fitter->SetParameter(11,"fSigma",fParamStartValues[11], 1.e-08, 0., 1.e+04);
208 fitter->SetParameter(12,"fAlpha",fParamStartValues[12], 1.e-08, 0., 1.e+04);
209 fitter->SetParameter(13,"fNorm",fParamStartValues[13], 1.e-08, 0., 1.e+04);
210 fitter->SetParameter(14,"fBkgNorm",fParamStartValues[14], 1.e-08, 0., 1.e+04);
211 fitter->SetParameter(15,"fBkgMean",fParamStartValues[15], 1.e-08, 0., 1.e+04);
212 fitter->SetParameter(16,"fBkgSlope",fParamStartValues[16], 1.e-08, 0., 1.e+04);
213 fitter->SetParameter(17,"fBkgConst",fParamStartValues[17], 1.e-08, 0., 1.e+04);
214 fitter->SetParameter(18,"fNormGaus1FF",fParamStartValues[18], 1.e-08, 0., 1.e+05);
215 fitter->SetParameter(19,"fNormGaus2FF",fParamStartValues[19], 1.e-08, 0., 1.e+05);
216 fitter->SetParameter(20,"fMean1ResFF",fParamStartValues[20], 1.e-08, 0., 1.e+05);
217 fitter->SetParameter(21,"fSigma1ResFF",fParamStartValues[21], 1.e-08, 0., 1.e+05);
218 fitter->SetParameter(22,"fMean2ResFF",fParamStartValues[22], 1.e-08, 0., 1.e+05);
219 fitter->SetParameter(23,"fSigma2ResFF",fParamStartValues[23], 1.e-08, 0., 1.e+05);
220 fitter->SetParameter(24,"fAlfaResFF",fParamStartValues[24], 1.e-08, 0., 1.e+05);
221 fitter->SetParameter(25,"fLambdaResFF",fParamStartValues[25], 1.e-08, 0., 1.e+05);
222 fitter->SetParameter(26,"fResNormExpFF",fParamStartValues[26], 1.e-08, 0., 1.e+05);
223 fitter->SetParameter(27,"fNormGaus1FS",fParamStartValues[27], 1.e-08, 0., 1.e+05);
224 fitter->SetParameter(28,"fNormGaus2FS",fParamStartValues[28], 1.e-08, 0., 1.e+05);
225 fitter->SetParameter(29,"fMean1ResFS",fParamStartValues[29], 1.e-08, 0., 1.e+05);
226 fitter->SetParameter(30,"fSigma1ResFS",fParamStartValues[30], 1.e-08, 0., 1.e+05);
227 fitter->SetParameter(31,"fMean2ResFS",fParamStartValues[31], 1.e-08, 0., 1.e+05);
228 fitter->SetParameter(32,"fSigma2ResFS",fParamStartValues[32], 1.e-08, 0., 1.e+05);
229 fitter->SetParameter(33,"fAlfaResFS",fParamStartValues[33], 1.e-08, 0., 1.e+05);
230 fitter->SetParameter(34,"fLambdaResFS",fParamStartValues[34], 1.e-08, 0., 1.e+05);
231 fitter->SetParameter(35,"fResNormExpFS",fParamStartValues[35], 1.e-08, 0., 1.e+05);
232 fitter->SetParameter(36,"fNormGaus1SS",fParamStartValues[36], 1.e-08, 0., 1.e+05);
233 fitter->SetParameter(37,"fNormGaus2SS",fParamStartValues[37], 1.e-08, 0., 1.e+05);
234 fitter->SetParameter(38,"fMean1ResSS",fParamStartValues[38], 1.e-08, 0., 1.e+05);
235 fitter->SetParameter(39,"fSigma1ResSS",fParamStartValues[39], 1.e-08, 0., 1.e+05);
236 fitter->SetParameter(40,"fMean2ResSS",fParamStartValues[40], 1.e-08, 0., 1.e+05);
237 fitter->SetParameter(41,"fSigma2ResSS",fParamStartValues[41], 1.e-08, 0., 1.e+05);
238 fitter->SetParameter(42,"fAlfaResSS",fParamStartValues[42], 1.e-08, 0., 1.e+05);
239 fitter->SetParameter(43,"fLambdaResSS",fParamStartValues[43], 1.e-08, 0., 1.e+05);
240 fitter->SetParameter(44,"fResNormExpSS",fParamStartValues[44], 1.e-08, 0., 1.e+05);
241 }
ba15fdfb 242
5720c765 243 for(UInt_t indexparam = 0; indexparam < 45; indexparam++){
244 if(IsParamFixed(indexparam)) fitter->FixParameter((Int_t)indexparam);
245 else fitter->ReleaseParameter((Int_t)indexparam);
246 }
247 Double_t arglist[2]={10000,0.1};
248 if(step == 2) {Int_t iret1 = fitter->ExecuteCommand("MINOS", arglist ,1); return iret1;}
ba15fdfb 249 Int_t iret=fitter->ExecuteCommand("MIGRAD", arglist ,2);
250 fitter->PrintResults(4,0);
251
5720c765 252 if(step == 3) {
253
254 TMinuit* minuitPoint = fitter->GetMinuit();
ba15fdfb 255
5720c765 256 TCanvas *c2 = new TCanvas("c2","contours",800,800);
257
258 //68.27% (1 sigma) confidence level for 2 parameters (fSIG versus fB)
259 minuitPoint->SetErrorDef(1.15); // 2.3/2
260 fContPlot1 = (TGraph*)minuitPoint->Contour(100,7,8);
261 fContPlot1->GetXaxis()->SetRange(0,1);
262 fContPlot1->GetYaxis()->SetRange(0,1);
263 fContPlot1->SetLineColor(42);
264 fContPlot1->SetLineWidth(3);
265
266 //95% (2 sigma) confidence level for 2 parameters (fSIG versus fB)
267 minuitPoint->SetErrorDef(2.995); // 5.99/2
268 fContPlot2 = (TGraph*)minuitPoint->Contour(100,7,8);
269 fContPlot2->GetXaxis()->SetRange(0,1);
270 fContPlot2->GetYaxis()->SetRange(0,1);
271 fContPlot2->SetLineColor(38);
272 fContPlot2->SetLineWidth(3);
273
274 //99.73% (3 sigma) confidence level for 2 parameters (fSIG versus fB)
275 minuitPoint->SetErrorDef(5.915); // 11.83/2
276 fContPlot3 = (TGraph*)minuitPoint->Contour(100,7,8);
277 fContPlot3->GetXaxis()->SetRange(0,1);
278 fContPlot3->GetXaxis()->SetTitle("f_{B}");
279 fContPlot3->GetYaxis()->SetTitle("f_{Sig}[2.4-4]");
280 fContPlot3->GetYaxis()->SetRange(0,1);
281 fContPlot3->SetLineColor(34);
282 fContPlot3->SetLineWidth(3);
283
284 fContPlot3->Draw("al");
285 fContPlot2->Draw("l");
286 fContPlot1->Draw("l");
287
288 c2->Draw();
289 }
290
291 AliInfo("Minimization procedure finished\n");
ba15fdfb 292 return iret;
293}
294//_______________________________________________________________________________________
295void AliDielectronBtoJPSItoEleCDFfitHandler::CdfFCN(Int_t & /* npar */,
296 Double_t * /* gin */,Double_t &f,Double_t *par,Int_t /* iflag */)
297{
298 //
299 // Definition of the FCN to be used by minuit
300 //
301 fLikely->SetAllParameters(par);
302 fLikely->ComputeMassIntegral();
303 if(fPrintStatus)fLikely->PrintStatus();
304
305 TStopwatch t;
306 t.Start();
307
5720c765 308 f = fLikely->EvaluateLikelihood(fX,fM,fType,fNcand);
ba15fdfb 309
310 t.Stop();
311 AliDebug(2,Form("Real time spent to calculate function == %f \n", t.RealTime()));
312 AliDebug(2,Form("CPU time spent to calculate function == %f \n", t.CpuTime()));
313 AliDebug(2,Form("Actual value of the AliDielectronBtoJPSItoEleCDFfitFCN function == %f \n", f));
314
315 return;
316}
317//_______________________________________________________________________________________
5720c765 318void AliDielectronBtoJPSItoEleCDFfitHandler::SetParamStartValues(Double_t inputparamvalues[45])
ba15fdfb 319{
5720c765 320 for(Int_t index=0; index < 45; index++) fParamStartValues[index] = inputparamvalues[index];
ba15fdfb 321}
322//_______________________________________________________________________________________
5720c765 323void AliDielectronBtoJPSItoEleCDFfitHandler::SetResolutionConstants(Double_t* resolutionConst, Int_t type)
ba15fdfb 324{
325 //
326 // Sets constants for the resolution function
327 //
5720c765 328 fLikely->SetResolutionConstants(resolutionConst,type);
ba15fdfb 329
330}
331//_______________________________________________________________________________________
332void AliDielectronBtoJPSItoEleCDFfitHandler::SetCrystalBallFunction(Bool_t okCB)
333{
334 //
335 // Sets the CB as the parametrization for the signal invariant mass spectrum
336 // (otherwise Landau is chosen)
337 //
338 fLikely->SetCrystalBallFunction(okCB);
339}
340//_______________________________________________________________________________________
341void AliDielectronBtoJPSItoEleCDFfitHandler::SetMassWndHigh(Double_t limit)
342{
343 //
344 // Sets upper limit for the invariant mass window (under J/PSI mass peak)
345 //
346 fLikely->SetMassWndHigh(limit);
347}
348//_______________________________________________________________________________________
349void AliDielectronBtoJPSItoEleCDFfitHandler::SetMassWndLow(Double_t limit)
350{
351 //
352 // Sets lower limit for the invariant mass window (under J/PSI mass peak)
353 //
354 fLikely->SetMassWndLow(limit);
355}
356