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