o updates for PbPb analysis of B->J/psi (Fiorella)
[u/mrichter/AliRoot.git] / PWGDQ / dielectron / AliDielectronBtoJPSItoEleCDFfitHandler.cxx
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>
16 #include <TFitter.h>
17 #include <TMinuit.h>
18 #include <TStopwatch.h>
19 #include <TCanvas.h>
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
32 void CDFFunction(Int_t &npar, Double_t *gin, Double_t &f, Double_t *x, Int_t iflag);
33
34 ClassImp(AliDielectronBtoJPSItoEleCDFfitHandler)
35
36         //______________________________________________________________________________
37 void 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 //_________________________________________________________________________________________________
48 AliDielectronBtoJPSItoEleCDFfitHandler::AliDielectronBtoJPSItoEleCDFfitHandler():
49         fIsParamFixed(49),
50         fPrintStatus(kFALSE),
51         fUp(0),
52         fX(0x0),
53         fM(0x0),
54         fPt(0x0),
55         fType(0x0),
56         fLikely(0x0),
57         fNcand(0),
58         fContPlot1(0x0),
59         fContPlot2(0x0),
60         fContPlot3(0x0),
61         fitter(0)
62 {
63         //
64         // default constructor
65         //
66 }
67 //_________________________________________________________________________________________________
68 AliDielectronBtoJPSItoEleCDFfitHandler::AliDielectronBtoJPSItoEleCDFfitHandler(Double_t* decaytime, 
69                 Double_t* invariantmass, Double_t *pt,Int_t *type, Int_t ncand) :
70         fIsParamFixed(49),
71         fPrintStatus(kFALSE),
72         fUp(0),
73         fX(decaytime),
74         fM(invariantmass),
75         fPt(pt),
76         fType(type),
77         fLikely(0x0),
78         fNcand(ncand),
79         fContPlot1(0x0),
80         fContPlot2(0x0),
81         fContPlot3(0x0),
82         fitter(0)
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");
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)");
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)");
139
140         AliInfo(Form("\n+++\n+++ Number of candidates ---> %d\n+++\n ", ncand));
141 }
142 //___________________________________________________________________________
143 AliDielectronBtoJPSItoEleCDFfitHandler& 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;
154                 fPt           = c.fPt;
155                 fType         = c.fType;
156                 fLikely       = c.fLikely;
157                 fNcand        = c.fNcand;
158                 fContPlot1    = c.fContPlot1;
159                 fContPlot2    = c.fContPlot2;
160                 fContPlot3    = c.fContPlot3;
161         }
162         return *this;
163 }
164
165 //_______________________________________________________________________________________
166 AliDielectronBtoJPSItoEleCDFfitHandler::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),
173         fPt(c.fPt),
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)
181 {
182         //
183         // Copy Constructor
184         //
185 }
186 //_______________________________________________________________________________________
187 AliDielectronBtoJPSItoEleCDFfitHandler::~AliDielectronBtoJPSItoEleCDFfitHandler()
188 {
189         //
190         //destructor
191         //
192         delete fLikely;
193 }
194 //_______________________________________________________________________________________
195 Int_t AliDielectronBtoJPSItoEleCDFfitHandler::DoMinimization(Int_t step)
196 {
197         //
198         // performs the minimization
199         //
200         if(step == 0 || !fitter){ 
201                 fitter = (TFitter*)TVirtualFitter::Fitter(this,49);
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);
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);
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);
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); 
253                 }
254
255         for(UInt_t indexparam = 0; indexparam < 49; indexparam++){
256                 if(IsParamFixed(indexparam)) fitter->FixParameter((Int_t)indexparam); 
257                 else fitter->ReleaseParameter((Int_t)indexparam);
258         }
259         Double_t arglist[2]={10000,0.1}; Int_t iret = 0;
260         if(step == 2) {Int_t  iret1 = fitter->ExecuteCommand("MINOS", arglist ,1); return iret1;}
261         if(step == 0) { fitter->SetParameter(8,"fFsig",fParamStartValues[8], 1.e-10, 0., 1.); 
262                          iret=fitter->ExecuteCommand("MIGRAD", arglist ,2); }
263         fitter->PrintResults(4,0);
264
265         if(step == 3) {
266
267                 TMinuit* minuitPoint = fitter->GetMinuit(); 
268
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                 
301                 c2->Draw();
302                 c2->SaveAs("contourPlot.root"); 
303         }
304
305         AliInfo("Minimization procedure finished\n");
306         return iret;
307 }
308 //_______________________________________________________________________________________
309 void 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
322         f = fLikely->EvaluateLikelihood(fX,fM,fPt, fType,fNcand);
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 //_______________________________________________________________________________________
332 void AliDielectronBtoJPSItoEleCDFfitHandler::SetParamStartValues(Double_t inputparamvalues[49])
333 {
334         for(Int_t index=0; index < 49; index++) fParamStartValues[index] = inputparamvalues[index];
335 }
336 //_______________________________________________________________________________________
337 void AliDielectronBtoJPSItoEleCDFfitHandler::SetResolutionConstants(Double_t* resolutionConst, Int_t type)
338 {
339         //
340         // Sets constants for the resolution function
341         //
342         fLikely->SetResolutionConstants(resolutionConst,type);
343
344 }
345 //_______________________________________________________________________________________
346 void 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 }
354
355 //_______________________________________________________________________________________
356 void 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
366 //_______________________________________________________________________________________
367 void 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 //_______________________________________________________________________________________
375 void 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