Change Mult binning scheme
[u/mrichter/AliRoot.git] / PWGDQ / dielectron / AliDielectronBtoJPSItoEleCDFfitFCNfitter.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 "AliDielectronBtoJPSItoEleCDFfitFCNfitter.h"
16 #include "TArrayD.h"
17 #include "TFormula.h"
18 #include "TF1.h"
19 #include "TF2.h"
20 #include "TH1F.h"
21 #include "TH2F.h"
22 #include "TMath.h"
23
24 AliDielectronBtoJPSItoEleCDFfitFCNfitter::AliDielectronBtoJPSItoEleCDFfitFCNfitter():
25  fInvMass(0x0),
26  fX(0x0),
27  fX2D(0x0),
28  fFCN(0x0),
29  fFitOpt("")
30 {
31  //
32  // ctor
33  //
34  fFCN = new AliDielectronBtoJPSItoEleCDFfitFCN();
35  fFCN->SetCrystalBallFunction(kTRUE);
36  for(Int_t i=0; i<kInvMassTotal; i++) fParameterInvMassToFix[i]=kFALSE;
37  for(Int_t i=0; i<kPseudo; i++)  fParameterXToFix[i]=kFALSE;
38  for(Int_t i=0; i<kPseudoBkg; i++)  fParameterXbkgToFix[i]=kFALSE;
39 }
40 //__________________________________________________________________________________________
41 AliDielectronBtoJPSItoEleCDFfitFCNfitter::~AliDielectronBtoJPSItoEleCDFfitFCNfitter()
42 {
43  //
44  // dtor
45  //
46  if(fInvMass) delete fInvMass;
47  if(fX) delete fX;
48  if(fX2D) delete fX2D;
49  if(fFCN) delete fFCN;
50 }
51 //__________________________________________________________________________________________
52 void AliDielectronBtoJPSItoEleCDFfitFCNfitter::SetInvMassSignalParameters(const Double_t massPar[5])
53 {
54  //
55  // Setter of the inv mass signal distribution parameters
56  //
57 fFCN->SetCrystalBallMmean(massPar[0]); //9
58 fFCN->SetCrystalBallNexp(massPar[1]); //10
59 fFCN->SetCrystalBallSigma(massPar[2]); //11
60 fFCN->SetCrystalBallAlpha(massPar[3]); //12
61 fFCN->SetCrystalBallNorm(massPar[4]); //13
62  //for(Int_t i=0; i<kInvMassSignal; i++) printf(" par%i %f \n",i,massPar[i]);
63 }
64 //__________________________________________________________________________________________
65 void AliDielectronBtoJPSItoEleCDFfitFCNfitter::SetInvMassParameters(const Double_t massPar[9])
66 {
67  //
68  // Setter of the inv mass total distribution parameters
69  //
70  Double_t massParSig[5];
71  Double_t massParBkg[4];
72
73  for(Int_t iPar=0; iPar<9; iPar++){
74   if(iPar<5) massParSig[iPar]=massPar[iPar];
75   else massParBkg[iPar-5]=massPar[iPar];
76  }
77  SetInvMassSignalParameters(massParSig);
78  SetInvMassBkgParameters(massParBkg);
79 }
80 //__________________________________________________________________________________________
81 void AliDielectronBtoJPSItoEleCDFfitFCNfitter::SetInvMassBkgParameters(const Double_t massPar[4])
82 {
83  //
84  // Setter of the Inv Mass background distribution parameters
85  //
86 fFCN->SetBkgInvMassNorm(massPar[0]);//14
87 fFCN->SetBkgInvMassMean(massPar[1]);//15
88 fFCN->SetBkgInvMassSlope(massPar[2]);//16
89 fFCN->SetBkgInvMassConst(massPar[3]);//17
90 }
91 //__________________________________________________________________________________________
92 void AliDielectronBtoJPSItoEleCDFfitFCNfitter::GetInvMassParameters(TArrayD &massPar)
93 {
94  //
95  // Getter of the inv mass total distribution parameters
96  //
97  massPar.Reset();
98  massPar.Set(kInvMassTotal);
99  massPar.AddAt(fFCN->GetCrystalBallMmean(),0); //9
100  massPar.AddAt(fFCN->GetCrystalBallNexp(),1); //10
101  massPar.AddAt(fFCN->GetCrystalBallSigma(),2); //11
102  massPar.AddAt(fFCN->GetCrystalBallAlpha(),3); //12
103  massPar.AddAt(fFCN->GetCrystalBallNorm(),4); //13
104  massPar.AddAt(fFCN->GetBkgInvMassNorm(),5);//14
105  massPar.AddAt(fFCN->GetBkgInvMassMean(),6);//15
106  massPar.AddAt(fFCN->GetBkgInvMassSlope(),7); //16
107  massPar.AddAt(fFCN->GetBkgInvMassConst(),8); //17
108 }
109 //__________________________________________________________________________________________
110 void AliDielectronBtoJPSItoEleCDFfitFCNfitter::GetInvMassSignalParameters(TArrayD &massPar)
111 {
112  //
113  // Getter of the Inv Mass signal distribution parameters
114  //
115  massPar.Reset();
116  massPar.Set(kInvMassSignal);
117  massPar.AddAt(fFCN->GetCrystalBallMmean(),0); //9
118  massPar.AddAt(fFCN->GetCrystalBallNexp(),1); //10
119  massPar.AddAt(fFCN->GetCrystalBallSigma(),2); //11
120  massPar.AddAt(fFCN->GetCrystalBallAlpha(),3); //12
121  massPar.AddAt(fFCN->GetCrystalBallNorm(),4); //13
122
123 }
124 //__________________________________________________________________________________________
125 void AliDielectronBtoJPSItoEleCDFfitFCNfitter::GetInvMassBkgParameters(TArrayD &massPar)
126 {
127  //
128  // Getter of the Inv Mass background distribution parameters
129  //
130  massPar.Reset();
131  massPar.Set(kInvMassBkg);
132  massPar.AddAt(fFCN->GetBkgInvMassNorm(),0); //14
133  massPar.AddAt(fFCN->GetBkgInvMassMean(),1); //15
134  massPar.AddAt(fFCN->GetCrystalBallSigma(),2); //11
135  massPar.AddAt(fFCN->GetCrystalBallAlpha(),3); //12
136 }
137
138 //__________________________________________________________________________________________
139 void AliDielectronBtoJPSItoEleCDFfitFCNfitter::SetPseudoProperBkgParameters(Double_t xBkgPar[7])
140 {
141  //
142  // Setter of the psudoproper background distribution parameters
143  //
144  fFCN->SetResWeight(xBkgPar[0]);//0
145  fFCN->SetFPlus(xBkgPar[1]);//1
146  fFCN->SetFMinus(xBkgPar[2]);//2
147  fFCN->SetFSym(xBkgPar[3]);//3
148  fFCN->SetLamPlus(xBkgPar[4]); //4
149  fFCN->SetLamMinus(xBkgPar[5]); // 5    
150  fFCN->SetLamSym(xBkgPar[6]); // 6
151
152 }
153 //__________________________________________________________________________________________
154 void AliDielectronBtoJPSItoEleCDFfitFCNfitter::GetPseudoProperBkgParameters(TArrayD &xBkgPar)
155 {
156  //
157  // Getter of the psudoproper background distribution parameters
158  //
159  xBkgPar.Reset();
160  xBkgPar.Set(kPseudoBkg);
161  xBkgPar.AddAt(fFCN->GetResWeight(),0);
162  xBkgPar.AddAt(fFCN->GetFPlus(),1);
163  xBkgPar.AddAt(fFCN->GetFMinus(),2);
164  xBkgPar.AddAt(fFCN->GetFSym(),3);
165  xBkgPar.AddAt(fFCN->GetLamPlus(),4);
166  xBkgPar.AddAt(fFCN->GetLamMinus(),5);
167  xBkgPar.AddAt(fFCN->GetLamSym(),6);
168 }
169 //__________________________________________________________________________________________
170 TH1F* AliDielectronBtoJPSItoEleCDFfitFCNfitter::FitInvMassSignal(Double_t norm,Double_t mMin, Double_t mMax)
171 {
172  //
173  // Fit on the invariant mass signal distribution ( Crystal ball params ) 
174  //
175  if(!fInvMass){
176   printf("no histogram...exiting\n");
177   return 0x0;
178  }
179
180  TArrayD pars;
181  GetInvMassSignalParameters(pars);
182
183  TF1 * invMass = new TF1("invMassSignal",this,&AliDielectronBtoJPSItoEleCDFfitFCNfitter::CDFInvMassSignal,mMin,mMax,kInvMassSignal+1);
184  for(Int_t i=0; i<kInvMassSignal ; i++) {
185   if(fParameterInvMassToFix[i]) invMass->FixParameter(i,pars.At(i));
186   else invMass->SetParameter(i,pars.At(i));
187  }
188  invMass->SetParameter(kInvMassSignal,norm);
189  fInvMass->Fit("invMassSignal",Form("R%s",fFitOpt.Data()),"",mMin,mMax);
190  return fInvMass;
191 }
192 //__________________________________________________________________________________________
193 TH1F * AliDielectronBtoJPSItoEleCDFfitFCNfitter::FitInvMass(Double_t norm[2],Double_t mMin, Double_t mMax)
194 {
195  //
196  // Fit on the invariant mass total distribution
197  //
198  if(!fInvMass){
199   printf("no histogram...exiting\n");
200   return 0x0;
201  }
202
203  TArrayD pars;
204  GetInvMassParameters(pars);
205
206  TF1 * invMassTot = new TF1("invMassTot",this,&AliDielectronBtoJPSItoEleCDFfitFCNfitter::CDFInvMassTotal,mMin,mMax,11);
207
208  for(Int_t i=0; i<kInvMassTotal ; i++) {
209   if(fParameterInvMassToFix[i]) invMassTot->FixParameter(i,pars.At(i));
210   else invMassTot->SetParameter(i,pars.At(i));
211  }
212
213  if(!fParameterInvMassToFix[5]) invMassTot->SetParLimits(5,0,999999999);
214
215
216  invMassTot->SetParameter(9,norm[0]);
217  invMassTot->SetParLimits(9,0,999999999);
218  invMassTot->SetParameter(10,norm[1]);
219  invMassTot->SetParLimits(10,0,999999999);
220
221  fInvMass->Fit("invMassTot",Form("R%s",fFitOpt.Data()),"",mMin,mMax);
222
223  //printf("Inv Mass Params : ");
224  //for(Int_t i=0; i< 11; i++) printf("p%i %f ",i,invMassTot->GetParameter(i));
225  //printf("\n");
226
227  return fInvMass;
228 }
229 //__________________________________________________________________________________________
230 void AliDielectronBtoJPSItoEleCDFfitFCNfitter::GetPseudoProperParameters(TArrayD &xPar, Int_t type)
231 {
232  //
233  // Getter of the psudoproper background distribution parameters per dielectron type (FF, FS, SS)
234  //
235  xPar.Reset();
236  xPar.Set(kPseudo);
237  xPar.AddAt(fFCN->GetNormGaus1ResFunc(type),0); 
238  xPar.AddAt(fFCN->GetNormGaus2ResFunc(type),3); 
239  xPar.AddAt(fFCN->GetResMean1(type),1); 
240  xPar.AddAt(fFCN->GetResSigma1(type),2); 
241  xPar.AddAt(fFCN->GetResMean2(type),4); 
242  xPar.AddAt(fFCN->GetResSigma2(type),5); 
243  xPar.AddAt(fFCN->GetResAlfa(type),6); 
244  xPar.AddAt(fFCN->GetResLambda(type),7); 
245  xPar.AddAt(fFCN->GetResNormExp(type),8); 
246 }
247 //__________________________________________________________________________________________
248 TH1F * AliDielectronBtoJPSItoEleCDFfitFCNfitter::FitResolutionFunction(Double_t xMin, Double_t xMax, Int_t type, Double_t norm)
249 {
250  //
251  // Fit on the resolution function distribution per dielectron type
252  //
253  if(!fX){
254   printf("no Pseudoproper decayhistogram...exiting\n");
255   return 0x0;
256  }
257
258  TArrayD pars;
259  GetPseudoProperParameters(pars,type);
260
261  TF1 * xPrompt = new TF1("xPrompt",this,&AliDielectronBtoJPSItoEleCDFfitFCNfitter::CDFResolutionFunction,xMin,xMax,11);
262  for(Int_t i=0; i<kPseudo ; i++) {
263   if(fParameterXToFix[i]) xPrompt->FixParameter(i,pars.At(i));
264   else xPrompt->SetParameter(i,pars.At(i));
265   printf("%f ",pars.At(i));
266  }
267  printf("\n");
268  if(!fParameterXToFix[3]) xPrompt->SetParLimits(3,0,999999999);
269  xPrompt->FixParameter(9,type);
270  xPrompt->SetParameter(10,norm);
271  fX->Fit("xPrompt",Form("R%s",fFitOpt.Data()),"",xMin,xMax);
272
273 // printf("Res Func type %i : ",type);
274 // for(Int_t i=0; i< 11; i++) printf("p%i %6.3f ",i,xPrompt->GetParameter(i));
275 // printf("\n");
276  return fX;
277 }
278
279 //__________________________________________________________________________________________
280 TH2F * AliDielectronBtoJPSItoEleCDFfitFCNfitter::FitBkgPsudoProper(Double_t xMin, Double_t xMax, Double_t norm)
281 {
282  //
283  // Fit on the resolution function distribution of the background ( in 2 D, whose one D is the dielectron
284  // type)
285  //
286
287  if(!fX2D){
288   printf("no BKG Pseudoproper decay 2D histogram ...exiting\n");
289   return 0x0;
290  }
291
292  TArrayD bkgParam;
293  GetPseudoProperBkgParameters(bkgParam);
294
295  TF2 *psProperBkgnd = new TF2("bkgPseudoProper",this,
296    &AliDielectronBtoJPSItoEleCDFfitFCNfitter::PsProperBackFitFunc,xMin,xMax,-0.5,2.5,8,"AliDielectronBtoJPSItoEleCDFfitFCNfitter","PsProperBackFitFunc");
297
298  for(Int_t i=0; i< kPseudoBkg ; i++) {
299   if(fParameterXbkgToFix[i]) psProperBkgnd->FixParameter(i,bkgParam.At(i));
300   else psProperBkgnd->SetParameter(i,bkgParam.At(i));
301  }
302  psProperBkgnd->SetParameter(kPseudoBkg,norm);
303
304  fX2D->Fit("bkgPseudoProper",Form("R%s",fFitOpt.Data()),"LEGO2",xMin,xMax);
305
306  return fX2D;
307 }
308
309 //__________________________________________________________________________________________
310 Double_t AliDielectronBtoJPSItoEleCDFfitFCNfitter::CDFInvMassSignal(const Double_t *x,const Double_t *par)
311 {
312  //
313  // function
314  //
315  SetInvMassSignalParameters(par);
316  return par[5]*fFCN->EvaluateCDFInvMassSigDistr(x[0]);
317 }
318 //__________________________________________________________________________________________
319 Double_t AliDielectronBtoJPSItoEleCDFfitFCNfitter::CDFInvMassBkg(const Double_t *x,const Double_t *par)
320 {
321  //
322  // function
323  //
324  SetInvMassBkgParameters(par);
325  return fFCN->EvaluateCDFInvMassBkgDistr(x[0]);
326 }
327 //__________________________________________________________________________________________
328 Double_t AliDielectronBtoJPSItoEleCDFfitFCNfitter::CDFInvMassTotal(const Double_t *x,const Double_t *par)
329 {
330  //
331  // function (it doesn't use the function integral->pay attention to it)
332  //
333  SetInvMassParameters(par);
334  //return ((EvaluateCDFInvMassSigDistr(x[0])))*par[9] + (EvaluateCDFInvMassBkgDistr(x[0]))*par[10];
335  return (((fFCN->EvaluateCDFInvMassSigDistr(x[0])))*par[9] + (1-par[9])*(fFCN->EvaluateCDFInvMassBkgDistr(x[0])))*par[10];
336 }
337 //__________________________________________________________________________________________
338 Double_t AliDielectronBtoJPSItoEleCDFfitFCNfitter::CDFResolutionFunction(const Double_t *x,const Double_t *par)
339 {
340  //
341  // function 
342  //
343 //  Double_t resParam[kPseudo];
344 //  for(Int_t i=0; i<kPseudo; i++) {
345 //   resParam[i] = par[i];
346 //  }
347  fFCN->SetResolutionConstants(par,(Int_t)par[9]);
348  return par[10]*fFCN->ResolutionFunc(x[0],(Int_t)(par[9]));
349 }
350 //__________________________________________________________________________________________
351 Double_t AliDielectronBtoJPSItoEleCDFfitFCNfitter::PsProperBackFitFunc(const Double_t* x, const Double_t* par)
352 {       
353  //
354  // function 
355  //
356  Int_t type = TMath::Nint(x[1]);
357  Double_t params[7];
358  for(Int_t i=0; i<7; i++) params[i]=par[i]; 
359  SetPseudoProperBkgParameters(params);          
360  return fFCN->EvaluateCDFDecayTimeBkgDistr(x[0],type)*par[7]*fFCN->GetResWeight(type);
361 }