Remove lgamma and replace by TMath::LnGamma which should work across platforms.
[u/mrichter/AliRoot.git] / PWG2 / EVCHAR / AliCentralityGlauberFit.cxx
CommitLineData
e6f3f2fe 1/**************************************************************************
2 * Copyright(c) 1998-1999, 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/* Origin: Alberica Toia, CERN, Alberica.Toia@cern.ch */
17
18///////////////////////////////////////////////////////////////////////////////
19// //
20// class to determine centrality percentiles from 1D distributions //
21// //
22///////////////////////////////////////////////////////////////////////////////
23
24#include <TNamed.h>
25#include <TH1D.h>
26#include <TString.h>
27#include <TFile.h>
28#include <TMath.h>
29#include <TROOT.h>
30#include <TH2F.h>
31#include <TF1.h>
32#include <TNtuple.h>
33#include <TStyle.h>
34#include <TGraphErrors.h>
35#include <vector>
87dea6ac 36#include <TMinuit.h>
37#include <TStopwatch.h>
38#include <TRandom.h>
e6f3f2fe 39#include "AliCentralityGlauberFit.h"
d20baea6 40#include "AliLog.h"
87dea6ac 41
e6f3f2fe 42ClassImp(AliCentralityGlauberFit)
43
44//______________________________________________________________________________
87dea6ac 45AliCentralityGlauberFit::AliCentralityGlauberFit(const char *filename) :
46 fNmu(0),
47 fMulow(0),
48 fMuhigh(0),
49 fNk(0),
50 fKlow(0),
51 fKhigh(0),
52 fNalpha(0),
53 fAlphalow(0),
54 fAlphahigh(0),
55 fNeff(0),
56 fEfflow(0),
57 fEffhigh(0),
58 fRebinFactor(0),
59 fMultmin(0),
60 fMultmax(0),
61 fGlauntuple(0),
62 fNpart(0),
63 fNcoll(0),
64 fB(0),
65 fTaa(0),
66 fEffi(0),
67 fhEffi(0),
68 fTempHist(0),
69 fGlauberHist(0),
70 fFastFit(0),
71 fNBD(0),
72 fUseChi2(kTRUE),
73 finrootfile(0),
74 foutrootfilename(0),
75 fhistnames(),
76 fPercentXsec(0)
77{
e6f3f2fe 78 // standard constructor
79 // glauber file
d20baea6 80 TFile *f = TFile::Open(filename);
87dea6ac 81 fGlauntuple = (TNtuple*) f->Get("nt_Pb_Pb");
82 fGlauntuple->SetBranchAddress("Npart",&fNpart);
83 fGlauntuple->SetBranchAddress("Ncoll",&fNcoll);
84 fGlauntuple->SetBranchAddress("B",&fB);
85 fGlauntuple->SetBranchAddress("tAA",&fTaa);
d20baea6 86
87 fNBD = new TF1 ("fNBD", AliCentralityGlauberFit::NBDFunc, 0, 100,2);
88
e6f3f2fe 89}
90
91//--------------------------------------------------------------------------------------------------
d20baea6 92void AliCentralityGlauberFit::SetRangeToFit(Float_t fmultmin, Float_t fmultmax)
93{
87dea6ac 94 // Set fit range.
95
96 fMultmin=fmultmin;
97 fMultmax=fmultmax;
d20baea6 98}
99
100//--------------------------------------------------------------------------------------------------
e6f3f2fe 101void AliCentralityGlauberFit::SetGlauberParam(
87dea6ac 102 Int_t Nmu,
103 Float_t mulow,
104 Float_t muhigh,
105 Int_t Nk,
106 Float_t klow,
107 Float_t khigh,
108 Int_t Nalpha,
109 Float_t alphalow,
110 Float_t alphahigh,
111 Int_t Neff,
112 Float_t efflow,
113 Float_t effhigh)
e6f3f2fe 114{
87dea6ac 115 // Set Glauber parameters.
116
117 fNmu=Nmu;
118 fMulow=mulow;
119 fMuhigh=muhigh;
120 fNk=Nk;
121 fKlow=klow;
122 fKhigh=khigh;
123 fNalpha=Nalpha;
124 fAlphalow=alphalow;
125 fAlphahigh=alphahigh;
126 fNeff=Neff;
127 fEfflow=efflow;
128 fEffhigh=effhigh;
e6f3f2fe 129}
130
131//--------------------------------------------------------------------------------------------------
87dea6ac 132void AliCentralityGlauberFit::MakeFits(TString infilename)
133{
e6f3f2fe 134 TH1D *hDATA;
d20baea6 135 TH1D *thistGlau;
136 FILE* fTxt = fopen ("parameters.txt","w");
e6f3f2fe 137 TFile *outrootfile;
138
139 // open inrootfile, outrootfile
87dea6ac 140 finrootfile = new TFile(infilename);
141 outrootfile = new TFile(foutrootfilename,"RECREATE");
e6f3f2fe 142
143 // loop over all distribution names
87dea6ac 144 std::vector<TString>::const_iterator hni;
145 for(hni=fhistnames.begin(); hni!=fhistnames.end(); hni++) {
146 hDATA = (TH1D*) (finrootfile->Get(*hni));
d20baea6 147 //TList *list = (TList*) (inrootfile->Get("coutput1"));
148 //hDATA = (TH1D*) (list->FindObject(*hni));
87dea6ac 149 hDATA->Rebin(fRebinFactor);
e6f3f2fe 150 TH1D *hGLAU = new TH1D("hGLAU","hGLAU",hDATA->GetNbinsX(),0,hDATA->GetNbinsX()*hDATA->GetBinWidth(1));
d20baea6 151 Float_t chi2min = 9999999.0;
e6f3f2fe 152 Float_t alpha_min=-1;
153 Float_t mu_min=-1;
154 Float_t k_min=-1;
155 Float_t eff_min=-1;
156 Float_t alpha, mu, k, eff, chi2;
157
87dea6ac 158 for (int nalpha=0;nalpha<fNalpha;nalpha++) {
159 alpha = fAlphalow + ((float) nalpha ) * (fAlphahigh - fAlphalow ) / fNalpha;
e6f3f2fe 160 mu=0.0;
87dea6ac 161 for (int nmu=0; nmu<fNmu; nmu++) {
162 mu = (fMulow*(1-nalpha*0.05) ) + ((float) nmu ) * (fMuhigh - fMulow ) / fNmu;
d20baea6 163 //mu = mulow + ((float) nmu ) * (muhigh - mulow ) / Nmu;
e6f3f2fe 164
87dea6ac 165 for (int nk=0; nk<fNk; nk++) {
166 k = fKlow + ((float) nk ) * (fKhigh - fKlow ) / fNk;
e6f3f2fe 167
87dea6ac 168 for (int neff=0; neff<fNeff; neff++) {
169 eff = fEfflow + ((float) neff) * (fEffhigh - fEfflow) / fNeff;
e6f3f2fe 170
171 thistGlau = GlauberHisto(mu,k,eff,alpha,hDATA,kFALSE);
172 chi2 = CalculateChi2(hDATA,thistGlau,eff);
d20baea6 173 fprintf(fTxt, "%3.3f %3.3f %3.3f %3.3f %3.3f\n",(float) eff, (float) mu, (float) k, (float) alpha, chi2);
174
e6f3f2fe 175 if ( chi2 < chi2min ) {
176 chi2min = chi2;
177 alpha_min=alpha;
178 mu_min=mu;
179 k_min=k;
180 eff_min=eff;
181 }
182 }
183 }
184 }
185 }
186
187 thistGlau = GlauberHisto(mu_min,k_min,eff_min,alpha_min,hDATA,kTRUE);
188 hGLAU = (TH1D *) thistGlau->Clone("hGLAU");
d20baea6 189 //hGLAU->SetName( ((TString)hDATA->GetName()).Append(Form("_GLAU_%d_%d_%d_%d",(int)(100*mu_min),(int)(100*k_min),(int)(100*alpha_min),(int)(100*eff_min))));
190 hGLAU->SetName( ((TString)hDATA->GetName()).Append(Form("_GLAU")));
191 hGLAU->SetTitle( ((TString)hDATA->GetName()).Append(Form("_GLAU_%.3f_%.3f_%.3f_%.3f",mu_min,k_min,alpha_min,eff_min)));
192
572a3604 193 float mcintegral = hGLAU->Integral(1,hGLAU->GetNbinsX());
194 float scale = (hDATA->Integral(1,hDATA->GetNbinsX())/mcintegral) * ((float) eff_min);
d20baea6 195 hGLAU->Scale(scale);
e6f3f2fe 196
87dea6ac 197 fhEffi = GetTriggerEfficiencyFunction(hDATA, hGLAU);
198 SaveHisto(hDATA,hGLAU,fhEffi,outrootfile);
d20baea6 199 fclose (fTxt);
e6f3f2fe 200
87dea6ac 201 std::cout << "chi2 min is " << chi2min << std::endl;
202 std::cout << "fitted " << hGLAU->Integral(hGLAU->FindBin(fMultmin),
203 hGLAU->FindBin(fMultmax))/hGLAU->Integral()
204 << " of the total cross section" << std::endl;
d20baea6 205 fTempHist=hDATA;
206 fGlauberHist=hGLAU;
e6f3f2fe 207 }
d20baea6 208
e6f3f2fe 209 // close inrootfile, outrootfile
87dea6ac 210 finrootfile->Close();
e6f3f2fe 211 outrootfile->Close();
e6f3f2fe 212}
213
214//--------------------------------------------------------------------------------------------------
87dea6ac 215void AliCentralityGlauberFit::MakeFitsMinuitNBD(TString infilename)
216{
217 // Make fits using Minut.
218
d20baea6 219 TH1D *hDATA;
220 TH1D *thistGlau;
221 TFile *outrootfile;
222
223 // open inrootfile, outrootfile
87dea6ac 224 finrootfile = new TFile(infilename);
225 outrootfile = new TFile(foutrootfilename,"RECREATE");
d20baea6 226
227 // loop over all distribution names
87dea6ac 228 std::vector<TString>::const_iterator hni;
229 for(hni=fhistnames.begin(); hni!=fhistnames.end(); hni++) {
d20baea6 230 fFastFit=0;
231 fUseChi2 = 1;
87dea6ac 232 hDATA = (TH1D*) (finrootfile->Get(*hni));
233 hDATA->Rebin(fRebinFactor);
d20baea6 234 fTempHist=hDATA;
235 TH1D *hGLAU = new TH1D("hGLAU","hGLAU",hDATA->GetNbinsX(),0,hDATA->GetNbinsX()*hDATA->GetBinWidth(1));
236 hGLAU->Sumw2();
237 // Minimize here
238 if(gMinuit) delete gMinuit;
239 new TMinuit(3);
240 gMinuit->mncler();
241 gMinuit->SetFCN(AliCentralityGlauberFit::MinuitFcnNBD);
242 gMinuit->SetObjectFit(this);
243
244 Double_t arglist[2]={0};
245 Int_t ierflg;
246 if (fUseChi2) arglist[0] = 1; // should be 1 if you want to minimize chi2
247 else arglist[0] = 0.5; // should be 0.5 for ll
248 gMinuit->mnexcm("SET ERR",arglist, 1, ierflg);
249
572a3604 250 //Change verbosity
251 if (0) {
252 //gMinuit->Command((TString("SET PRINTOUT ")+long(fVerbosity)).Data());
253 }
d20baea6 254
255 if (!fFastFit) {
572a3604 256 gMinuit->mnparm(0,"alpha", fAlphalow, 0.1, 0.75,1, ierflg);
257 gMinuit->mnparm(1,"mu" , fMulow, 0.2, 1,100, ierflg);
258 gMinuit->mnparm(2,"k" , fKlow, 0.1, 0.5,2.5, ierflg);
259 gMinuit->mnparm(3,"eff" , fEfflow, 0.1, 0.95,1.0, ierflg);
d20baea6 260 } else if (fFastFit == 2) {
d20baea6 261 gMinuit->mnparm(0,"alpha", 0.59, 0.1, 0,1, ierflg);
262 gMinuit->mnparm(1,"mu" , 30.2, 1, 0,0, ierflg);
263 gMinuit->mnparm(2,"k" , 3.875, 0.1, 0,0, ierflg);
d20baea6 264 gMinuit->mnparm(3,"eff" , 0.98947, 0.1, 0,0, ierflg);
265 }
e6f3f2fe 266
d20baea6 267 // Call migrad
268 // arglist[0] = 100000; // max calls
269 // arglist[1] = 0.1; // tolerance
270 // arglist[0] = 50; // max calls
271 // arglist[1] = 1; // tolerance
272 // gMinuit->mnexcm("MIGrad",arglist, 2, ierflg);
273 arglist[0] = 100; // max calls
274 arglist[1] = 0.1; // tolerance
275 gMinuit->mnexcm("SIMPLEX",arglist, 2, ierflg);
d20baea6 276 arglist[0] = 1000; // max calls
277 arglist[1] = 0.1; // tolerance
278 gMinuit->mnexcm("MIGrad",arglist, 2, ierflg);
572a3604 279 //gMinuit->mnexcm("IMPROVE",arglist, 0, ierflg);
d20baea6 280
281 if (ierflg != 0) {
282 AliWarning("Abnormal termination of minimization.");
d20baea6 283 }
284
285 // if(opt.Contains("M")) {
286 // gMinuit->mnexcm("IMPROVE",arglist, 0, ierflg);
287 // }
288 // if(opt.Contains("E")) {
289 // gMinuit->mnexcm("HESSE",arglist, 0, ierflg);
290 // gMinuit->mnexcm("MINOS",arglist, 0, ierflg);
291 // }
292
293 // ______________________ Get chi2 and Fit Status __________
294 Double_t amin,edm,errdef;
295 int nvpar, nparx, icstat;
e6f3f2fe 296
d20baea6 297
298 gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
299 Double_t chi2min = amin;
87dea6ac 300 std::cout << "Fit status " << icstat << std::endl;
d20baea6 301
302 Double_t alpha_min, mu_min, k_min, eff_min;
303 Double_t alpha_mine, mu_mine, k_mine, eff_mine;
304 gMinuit->GetParameter(0, alpha_min , alpha_mine );
305 gMinuit->GetParameter(1, mu_min , mu_mine );
306 gMinuit->GetParameter(2, k_min , k_mine );
307 gMinuit->GetParameter(3, eff_min , eff_mine );
308
309 // Print some infos
87dea6ac 310 std::cout << "chi2 min is " << chi2min << ", " << alpha_min << ", "<< mu_min<< ", " << k_min << ", " << eff_min << std::endl;
d20baea6 311
312 thistGlau = GlauberHisto(mu_min,k_min,eff_min,alpha_min,hDATA,kTRUE);
313 hGLAU = (TH1D *) thistGlau->Clone("hGLAU");
314 hGLAU->SetName( ((TString)hDATA->GetName()).Append(Form("_GLAU")));
315 hGLAU->SetTitle( ((TString)hDATA->GetName()).Append(Form("_GLAU_%.3f_%.3f_%.3f_%.3f",mu_min,k_min,alpha_min,eff_min)));
316
317
87dea6ac 318 std::cout << "fitted " << hGLAU->Integral(hGLAU->FindBin(fMultmin),
319 hGLAU->FindBin(fMultmax))/hGLAU->Integral()
320 << " of the total cross section" << std::endl;
d20baea6 321
572a3604 322 float mcintegral = hGLAU->Integral(1,hGLAU->GetNbinsX());
323 float scale = (hDATA->Integral(1,hDATA->GetNbinsX())/mcintegral) * ((float) eff_min);
324 hGLAU->Scale(scale);
d20baea6 325
572a3604 326 std::cout << "Chi2 final " << CalculateChi2(hDATA,hGLAU,eff_min) << std::endl;
d20baea6 327
87dea6ac 328 fhEffi = GetTriggerEfficiencyFunction(hDATA, hGLAU);
329 SaveHisto(hDATA,hGLAU,fhEffi,outrootfile);
d20baea6 330 fGlauberHist=hGLAU;
331 }
d20baea6 332 // close inrootfile, outrootfile
87dea6ac 333 finrootfile->Close();
d20baea6 334 outrootfile->Close();
d20baea6 335}
336
d20baea6 337//--------------------------------------------------------------------------------------------------
87dea6ac 338TH1D *AliCentralityGlauberFit::GlauberHisto(Float_t mu, Float_t k, Float_t eff, Float_t alpha,
572a3604 339 TH1D *hDATA, Bool_t save)
87dea6ac 340{
341 // Get Glauber histogram
d20baea6 342
87dea6ac 343 eff=eff;
d20baea6 344 static TStopwatch sw;
345
346 TH1D *hSample = fFastFit != 2 ? NBDhist(mu,k) : 0;
347 // fNBD->SetParameters(mu,k);
348 static TH1D *h1 = (TH1D*)hDATA->Clone();
349 h1->Reset();
350 // h1->SetName(Form("fit_%.3f_%.3f_%.3f_%.3f",mu,k,eff,alpha));
e6f3f2fe 351 TFile *outFile = NULL;
352 TNtuple *ntuple = NULL;
353
354 if (save) {
355 outFile = TFile::Open("pippo.root", "RECREATE");
d20baea6 356 ntuple = new TNtuple("gnt", "Glauber ntuple", "Npart:Ncoll:B:tAA:ntot");
e6f3f2fe 357 }
d20baea6 358
359 Int_t nents = 100000;//glau_ntuple->GetEntries();
360 //if(fFastFit > 0 && fFastFit < 2) nents = 10000;
361 // Int_t nents = 100;//glau_ntuple->GetEntries();
e6f3f2fe 362 for (Int_t i=0;i<nents;++i) {
87dea6ac 363 fGlauntuple->GetEntry(i);
d20baea6 364 //Int_t n = TMath::Nint(TMath::Power(Npart,alpha));
365 //Int_t n = TMath::Nint(TMath::Power(Ncoll,alpha));
87dea6ac 366 Int_t n = alpha * fNpart + (1-alpha) * fNcoll;
e6f3f2fe 367 Int_t ntot=0;
d20baea6 368 // ntot=n*fNBD->GetRandom();
369 if (fFastFit == 1) {
370 // NBD
371 ntot = n*hSample->GetRandom();
372 }
373 else if (fFastFit == 2) {
374 // Gaussian
375 Double_t sigma = k*TMath::Sqrt(n*mu);
376 ntot = gRandom->Gaus(n*mu,sigma);
377 }
378 else {
379 for(Int_t j = 0; j<n; ++j)
380 ntot+=hSample->GetRandom();
381 // ntot+=fNBD->GetRandom();
382 }
e6f3f2fe 383 h1->Fill(ntot);
384
385 if (save)
87dea6ac 386 ntuple->Fill(fNpart,fNcoll,fB,fTaa,ntot);
e6f3f2fe 387
388 }
d20baea6 389
e6f3f2fe 390 if (save) {
391 ntuple->Write();
392 outFile->Close();
393 }
394
d20baea6 395 if (fFastFit != 2) delete hSample;
e6f3f2fe 396 return h1;
397
398}
399
400//--------------------------------------------------------------------------------------------------
8dc7e599 401Float_t AliCentralityGlauberFit::CalculateChi2(TH1D *hDATA, TH1D *thistGlau, Float_t eff)
402{
e6f3f2fe 403 // note that for different values of neff the mc distribution is identical, just re-scaled below
404 // normalize the two histogram
405 // scale MC up but leave data alone - that preserves error bars !!!!
406
572a3604 407 int lowchibin = hDATA->FindBin(fMultmin);
408 int highchibin = hDATA->FindBin(fMultmax);
d20baea6 409
572a3604 410 float mcintegral = thistGlau->Integral(1,thistGlau->GetNbinsX());
411 float scale = (hDATA->Integral(1,hDATA->GetNbinsX())/mcintegral) * ((float) eff);
d20baea6 412 thistGlau->Scale(scale);
413
414 // FIXME: Kolmogorov
415 //return hDATA->KolmogorovTest(thistGlau,"M");
416
572a3604 417 // // calculate the chi2 between MC and real data over some range ????
d20baea6 418 if (fUseChi2) {
419 float chi2 = 0.0;
420 for (int i=lowchibin; i<=highchibin; i++) {
421 if (hDATA->GetBinContent(i) < 1.0) continue;
422 float diff = pow( (thistGlau->GetBinContent(i) - hDATA->GetBinContent(i)) , 2);
572a3604 423 diff = diff / (pow(hDATA->GetBinError(i) , 2)+pow(thistGlau->GetBinError(i) , 2)); // FIXME squared distance if commented
d20baea6 424 chi2 += diff;
425 }
426 chi2 = chi2 / (highchibin - lowchibin + 1);
427 return chi2;
572a3604 428 }
429 // "-2 log likelihood ratio(mu;n) = 2[(mu - n) + n ln(n/mu)]"
430 else {
87dea6ac 431 std::cout << "LL" << std::endl;
572a3604 432
d20baea6 433
434 float ll = 0.0;
435 for (int i=lowchibin; i<=highchibin; i++) {
572a3604 436 // if (thistGlau->GetBinContent(i) < 1) continue;
437 // if (hDATA->GetBinContent(i) < 1) continue;
438 // cout << ll << " " << thistGlau->GetBinContent(i) <<" " << hDATA->GetBinContent(i) << endl;
439
d20baea6 440 Float_t data = hDATA ->GetBinContent(i);
441 Float_t mc = thistGlau->GetBinContent(i);
442 Int_t idata = TMath::Nint(data);
443 if (mc < 1.e-9) mc = 1.e-9;
444 Float_t fsub = - mc + idata * TMath::Log(mc);
445 Float_t fobs = 0;
87dea6ac 446 //Int_t imc = TMath::Nint(mc);
d20baea6 447 if (idata > 0) {
448 for(Int_t istep = 0; istep < idata; istep++){
449 if (istep > 1)
450 fobs += TMath::Log(istep);
451 }
452 }
572a3604 453 // cout << mc << " " << data << " " << fsub << " " << fobs << endl;
d20baea6 454 fsub -= fobs;
455 ll -= fsub ;
456 }
457 return 2*ll;
e6f3f2fe 458 }
572a3604 459
e6f3f2fe 460}
461
462//--------------------------------------------------------------------------------------------------
87dea6ac 463TH1D * AliCentralityGlauberFit::GetTriggerEfficiencyFunction(TH1D *hist1, TH1D *hist2)
464{
465 // Get efficiency.
466
467 fhEffi = (TH1D*)hist1->Clone("heffi");
468 fhEffi->Divide(hist2);
469 return fhEffi;
e6f3f2fe 470}
471
472//--------------------------------------------------------------------------------------------------
87dea6ac 473Float_t AliCentralityGlauberFit::GetTriggerEfficiencyIntegral(TH1D *hist1, TH1D *hist2)
474{
475 // Get eff integral.
476 fEffi = hist1->Integral()/hist2->Integral();
477 return fEffi;
e6f3f2fe 478}
479
480//--------------------------------------------------------------------------------------------------
572a3604 481void AliCentralityGlauberFit::SaveHisto(TH1D *hist1, TH1D *hist2, TH1D *heffi, TFile *outrootfile)
482{
e6f3f2fe 483 outrootfile->cd();
484 hist1->Write();
485 hist2->Write();
486 heffi->Write();
487}
488
489//--------------------------------------------------------------------------------------------------
e6f3f2fe 490Double_t AliCentralityGlauberFit::NBD(Int_t n, Double_t mu, Double_t k)
491{
87dea6ac 492 // Compute NBD.
4cce943f 493 Double_t ret = TMath::Exp( TMath::LnGamma(n+k) - TMath::LnGamma(k) - TMath::LnGamma(n+1) ) *
494 TMath::Power(mu/(mu+k),n) * TMath::Power(1-mu/(mu+k),k);
e6f3f2fe 495 return ret;
496}
497
498//--------------------------------------------------------------------------------------------------
e6f3f2fe 499TH1D *AliCentralityGlauberFit::NBDhist(Double_t mu, Double_t k)
500{
d20baea6 501 TH1D *h = new TH1D("htemp","",100,-0.5,299.5);
e6f3f2fe 502 h->SetName(Form("nbd_%f_%f",mu,k));
503 h->SetDirectory(0);
d20baea6 504 for (Int_t i=0;i<300;++i) {
e6f3f2fe 505 Double_t val = NBD(i,mu,k);
d20baea6 506 if (val>1e-20) h->Fill(i,val);
e6f3f2fe 507 }
508 return h;
509}
510
572a3604 511//--------------------------------------------------------------------------------------------------
87dea6ac 512void AliCentralityGlauberFit::MinuitFcnNBD(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
513{
d20baea6 514 // FCN for minuit
515 Double_t alpha = par[0];
516 Double_t mu = par[1];
517 Double_t k = par[2];
518 Double_t eff = par[3];
519 //Double_t eff = 1;//par[3];
8ee6586a 520 if (0) { //avoid warning
521 gin=gin;
522 npar=npar;
523 iflag=iflag;
524 }
d20baea6 525 AliCentralityGlauberFit * obj = (AliCentralityGlauberFit *) gMinuit->GetObjectFit();
526
527 // static TStopwatch sw;
528 // sw.Start();
529 TH1D * thistGlau = obj->GlauberHisto(mu,k,eff,alpha,obj->GetTempHist(),kFALSE);
530 f = obj->CalculateChi2(obj->GetTempHist(),thistGlau,eff);
531 // sw.Stop();
532 // sw.Print();
87dea6ac 533 Printf("%f - %f - %f - %f - %f \n",f,alpha,mu,k,eff);
d20baea6 534}
535
572a3604 536//--------------------------------------------------------------------------------------------------
87dea6ac 537Double_t AliCentralityGlauberFit::NBDFunc(Double_t * x, Double_t *par)
538{
d20baea6 539 // TF1 interface
d20baea6 540 Double_t mu = par[0];
541 Double_t k = par[1];
542 Double_t n = x[0];
d20baea6 543 Double_t ret = exp( lgamma(n+k) - lgamma(k) - lgamma(n+1) ) * TMath::Power(mu/(mu+k),n) * TMath::Power(1-mu/(mu+k),k);
d20baea6 544 return ret;
d20baea6 545}