Cosmetics
[u/mrichter/AliRoot.git] / PWG2 / EVCHAR / AliCentralityGlauberFit.cxx
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>
36 #include <TMinuit.h>
37 #include <TStopwatch.h>
38 #include <TRandom.h>
39 #include "AliCentralityGlauberFit.h"
40 #include "AliLog.h"
41
42 ClassImp(AliCentralityGlauberFit)  
43  
44 //______________________________________________________________________________
45 AliCentralityGlauberFit::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 {
78   // standard constructor
79   // glauber file
80   TFile *f = TFile::Open(filename);
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);
86   
87   fNBD = new TF1 ("fNBD", AliCentralityGlauberFit::NBDFunc, 0, 100,2);
88
89 }
90
91 //--------------------------------------------------------------------------------------------------
92 void AliCentralityGlauberFit::SetRangeToFit(Float_t fmultmin, Float_t fmultmax)
93 {
94   // Set fit range.
95
96   fMultmin=fmultmin;
97   fMultmax=fmultmax;
98 }
99
100 //--------------------------------------------------------------------------------------------------
101 void AliCentralityGlauberFit::SetGlauberParam(
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)
114 {
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;
129 }
130
131 //--------------------------------------------------------------------------------------------------
132 void AliCentralityGlauberFit::MakeFits(TString infilename) 
133 {
134   TH1D *hDATA;
135   TH1D *thistGlau; 
136   FILE* fTxt = fopen ("parameters.txt","w");
137   TFile *outrootfile;
138   
139   // open inrootfile, outrootfile
140   finrootfile  = new TFile(infilename);
141   outrootfile = new TFile(foutrootfilename,"RECREATE");
142   
143   // loop over all distribution names  
144   std::vector<TString>::const_iterator hni;
145   for(hni=fhistnames.begin(); hni!=fhistnames.end(); hni++) {
146     hDATA  = (TH1D*) (finrootfile->Get(*hni)); 
147     //TList *list  = (TList*) (inrootfile->Get("coutput1")); 
148     //hDATA  = (TH1D*) (list->FindObject(*hni)); 
149     hDATA->Rebin(fRebinFactor);
150     TH1D *hGLAU = new TH1D("hGLAU","hGLAU",hDATA->GetNbinsX(),0,hDATA->GetNbinsX()*hDATA->GetBinWidth(1));
151     Float_t chi2min = 9999999.0;
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
158     for (int nalpha=0;nalpha<fNalpha;nalpha++) {
159       alpha = fAlphalow   + ((float) nalpha ) * (fAlphahigh  - fAlphalow ) / fNalpha;
160       mu=0.0;
161       for (int nmu=0; nmu<fNmu; nmu++) {
162         mu = (fMulow*(1-nalpha*0.05) )  + ((float) nmu ) * (fMuhigh  - fMulow ) / fNmu;
163         //mu = mulow   + ((float) nmu ) * (muhigh  - mulow ) / Nmu;
164         
165         for (int nk=0; nk<fNk; nk++) {
166           k = fKlow  + ((float) nk ) * (fKhigh  - fKlow ) / fNk;
167           
168           for (int neff=0; neff<fNeff; neff++) {
169             eff = fEfflow + ((float) neff) * (fEffhigh - fEfflow) / fNeff;
170             
171             thistGlau = GlauberHisto(mu,k,eff,alpha,hDATA,kFALSE);
172             chi2 = CalculateChi2(hDATA,thistGlau,eff);
173             fprintf(fTxt, "%3.3f %3.3f %3.3f %3.3f %3.3f\n",(float) eff, (float) mu, (float) k, (float) alpha, chi2);
174
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");
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
193     float mcintegral = hGLAU->Integral(hDATA->FindBin(fMultmin),hGLAU->GetNbinsX());
194     float scale = (hDATA->Integral(hDATA->FindBin(fMultmin),hDATA->GetNbinsX())/mcintegral) * ((float) eff_min);
195     hGLAU->Scale(scale);
196
197     fhEffi = GetTriggerEfficiencyFunction(hDATA, hGLAU);
198     SaveHisto(hDATA,hGLAU,fhEffi,outrootfile);
199     fclose (fTxt);
200
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; 
205     fTempHist=hDATA;
206     fGlauberHist=hGLAU;
207   }
208
209   // close inrootfile, outrootfile
210   finrootfile->Close();
211   outrootfile->Close();
212 }
213
214 //--------------------------------------------------------------------------------------------------
215 void AliCentralityGlauberFit::MakeFitsMinuitNBD(TString infilename) 
216 {
217   // Make fits using Minut.
218
219   TH1D *hDATA;
220   TH1D *thistGlau; 
221   TFile *outrootfile;
222   
223   // open inrootfile, outrootfile
224   finrootfile  = new TFile(infilename);
225   outrootfile = new TFile(foutrootfilename,"RECREATE");
226   
227   // loop over all distribution names  
228   std::vector<TString>::const_iterator hni;
229   for(hni=fhistnames.begin(); hni!=fhistnames.end(); hni++) {
230     fFastFit=0;
231     fUseChi2 = 1;
232     hDATA  = (TH1D*) (finrootfile->Get(*hni)); 
233     hDATA->Rebin(fRebinFactor);
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
250   // // Change verbosity
251   // gMinuit->Command((TString("SET PRINTOUT ")+long(fVerbosity)).Data());
252
253     // Set parameters
254     // start from the middle of the allowed range, as a guess
255     // gMinuit->mnparm(0,"alpha", (alphalow+alphahigh)/2, 1, alphalow, alphahigh, ierflg);
256     // gMinuit->mnparm(1,"mu"   , (mulow   +muhigh   )/2, 1, mulow   , muhigh   , ierflg);
257     // gMinuit->mnparm(2,"k"    , (klow    +khigh    )/2, 1, klow    , khigh    , ierflg);
258     // gMinuit->mnparm(3,"eff"  , (efflow  +effhigh  )/2, 1, efflow  , effhigh   , ierflg);
259
260     // NBD with limits
261     // gMinuit->mnparm(0,"alpha", 1.104, 0.1, alphalow, alphahigh, ierflg);
262     // gMinuit->mnparm(1,"mu"   , 65,  1,    mulow   , muhigh   , ierflg);
263     // gMinuit->mnparm(2,"k"    , 0.41, 0.1, klow    , khigh    , ierflg);
264     // //    gMinuit->mnparm(3,"eff"  , 1,   1, 0.9       , 1        , ierflg);
265     // gMinuit->mnparm(3,"eff"  , 0.95, 0.1, efflow  , effhigh   , ierflg);
266
267     // // NBD, no limits
268     // gMinuit->mnparm(0,"alpha", 1.104, 0.1, 0,0, ierflg);
269     // gMinuit->mnparm(1,"mu"   , 65,  1,     0,0, ierflg);
270     // gMinuit->mnparm(2,"k"    , 0.41, 0.1,  0,0, ierflg);
271     // gMinuit->mnparm(3,"eff"  , 0.95, 0.1,  0,0 , ierflg);
272
273     if (!fFastFit) {
274     // // NBD, Npart**alpha     
275       //gMinuit->mnparm(0,"alpha", 1.16, 0.1, 0,0, ierflg);
276       //gMinuit->mnparm(1,"mu"   , 27.55,   1,    0,0, ierflg);
277       //gMinuit->mnparm(2,"k"    , 0.64, 0.1,  0,0, ierflg);
278       //gMinuit->mnparm(3,"eff"  , 0.97, 0.1,  0,0 , ierflg);
279       // V0 non rescaled
280       // // NBD, alpha Npart * (1-alpha) Ncoll     
281       gMinuit->mnparm(0,"alpha", 0.832, 0.1, 0,0, ierflg);
282       gMinuit->mnparm(1,"mu"   , 29,   1,    0,0, ierflg);
283       gMinuit->mnparm(2,"k"    , 1.4, 0.1,  0,0, ierflg);
284       gMinuit->mnparm(3,"eff"  , 0.991, 0.1,  0,0 , ierflg);
285       // SPD
286       //gMinuit->mnparm(0,"alpha", 0.820, 0.1, 0.7,1, ierflg);
287       //gMinuit->mnparm(1,"mu"   , 8.267,   1,  1,100, ierflg);
288       //gMinuit->mnparm(2,"k"    , 0.530, 0.1,  0,0, ierflg);
289       //gMinuit->mnparm(3,"eff"  , 0.990, 0.1,  0,0 , ierflg);      
290     } else if (fFastFit == 2) {
291       //gaussian
292       //0.859136, 30.2057, 3.87565, 0.989747
293       gMinuit->mnparm(0,"alpha", 0.59, 0.1, 0,1, ierflg);
294       gMinuit->mnparm(1,"mu"   , 30.2,  1,     0,0, ierflg);
295       gMinuit->mnparm(2,"k"    , 3.875, 0.1,  0,0, ierflg);
296       //    gMinuit->mnparm(3,"eff"  , 1,   1, 0.9       , 1        , ierflg);
297       gMinuit->mnparm(3,"eff"  , 0.98947, 0.1,  0,0, ierflg);
298     }
299
300     // Call migrad
301     // arglist[0] = 100000; // max calls
302     // arglist[1] = 0.1; // tolerance    
303     // arglist[0] = 50; // max calls
304     // arglist[1] = 1; // tolerance    
305     // gMinuit->mnexcm("MIGrad",arglist, 2, ierflg);
306     arglist[0] = 100; // max calls
307     arglist[1] = 0.1; // tolerance    
308     gMinuit->mnexcm("SIMPLEX",arglist, 2, ierflg);
309     // //    fFastFit = 2;
310     arglist[0] = 1000; // max calls
311     arglist[1] = 0.1; // tolerance    
312     gMinuit->mnexcm("MIGrad",arglist, 2, ierflg);
313     //    gMinuit->mnexcm("IMPROVE",arglist, 0, ierflg);
314
315     if (ierflg != 0) {
316       AliWarning("Abnormal termination of minimization.");
317       //    exit(0);
318     }
319     
320     // if(opt.Contains("M")) {
321     //   gMinuit->mnexcm("IMPROVE",arglist, 0, ierflg);
322     // }
323     // if(opt.Contains("E")) {
324     // gMinuit->mnexcm("HESSE",arglist, 0, ierflg);
325     // gMinuit->mnexcm("MINOS",arglist, 0, ierflg);
326     // }
327
328     // ______________________ Get chi2 and Fit Status __________
329     Double_t amin,edm,errdef;
330     int      nvpar, nparx, icstat;
331     
332     
333     gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
334     Double_t chi2min = amin;
335     std::cout << "Fit status " << icstat << std::endl;
336
337     Double_t alpha_min, mu_min, k_min, eff_min;
338     Double_t alpha_mine, mu_mine, k_mine, eff_mine;
339     gMinuit->GetParameter(0, alpha_min , alpha_mine );
340     gMinuit->GetParameter(1, mu_min    , mu_mine    );
341     gMinuit->GetParameter(2, k_min     , k_mine     );
342     gMinuit->GetParameter(3, eff_min   , eff_mine   );
343
344     // Print some infos
345     std::cout << "chi2 min is " << chi2min << ", " << alpha_min << ", "<< mu_min<< ", "  << k_min << ", " <<  eff_min << std::endl;
346
347     thistGlau = GlauberHisto(mu_min,k_min,eff_min,alpha_min,hDATA,kTRUE);
348     hGLAU = (TH1D *) thistGlau->Clone("hGLAU");
349     hGLAU->SetName( ((TString)hDATA->GetName()).Append(Form("_GLAU")));
350     hGLAU->SetTitle( ((TString)hDATA->GetName()).Append(Form("_GLAU_%.3f_%.3f_%.3f_%.3f",mu_min,k_min,alpha_min,eff_min)));
351
352     
353     std::cout << "fitted " << hGLAU->Integral(hGLAU->FindBin(fMultmin),
354                                               hGLAU->FindBin(fMultmax))/hGLAU->Integral() 
355               << " of the total cross section" << std::endl; 
356     //    eff_min=1;
357
358     // Save histogram and ntuple
359     //    fFastFit=kFALSE; // Use a better aprox to get the final histo
360
361     // float mcintegral = hGLAU->Integral(hDATA->FindBin(multmin),hGLAU->GetNbinsX());
362     // float scale = (hDATA->Integral(hDATA->FindBin(multmin),hDATA->GetNbinsX())/mcintegral) * ((float) eff_min);
363     // hGLAU->Scale(scale);
364
365     std::cout << "Chi2 final" << CalculateChi2(hDATA,hGLAU,eff_min) << std::endl;
366
367     fhEffi = GetTriggerEfficiencyFunction(hDATA, hGLAU);
368     SaveHisto(hDATA,hGLAU,fhEffi,outrootfile);
369     fGlauberHist=hGLAU;
370   }
371   // close inrootfile, outrootfile
372   finrootfile->Close();
373   outrootfile->Close();
374 }
375
376 //--------------------------------------------------------------------------------------------------
377 TH1D *AliCentralityGlauberFit::GlauberHisto(Float_t mu, Float_t k, Float_t eff, Float_t alpha, 
378                                              TH1D *hDATA, Bool_t save) 
379 {
380   // Get Glauber histogram
381
382   eff=eff;
383   static TStopwatch sw;
384
385   TH1D *hSample = fFastFit != 2 ? NBDhist(mu,k) : 0;
386   //  fNBD->SetParameters(mu,k);
387   static TH1D *h1 = (TH1D*)hDATA->Clone();
388   h1->Reset();
389   //  h1->SetName(Form("fit_%.3f_%.3f_%.3f_%.3f",mu,k,eff,alpha));
390   TFile *outFile = NULL;
391   TNtuple *ntuple = NULL;
392  
393   if (save) {
394     outFile = TFile::Open("pippo.root", "RECREATE");
395     ntuple = new TNtuple("gnt", "Glauber ntuple", "Npart:Ncoll:B:tAA:ntot");
396   } 
397   
398
399   Int_t nents = 100000;//glau_ntuple->GetEntries();
400   //if(fFastFit > 0 && fFastFit < 2) nents = 10000;
401   //  Int_t nents = 100;//glau_ntuple->GetEntries();
402   for (Int_t i=0;i<nents;++i) {
403     fGlauntuple->GetEntry(i);
404     //Int_t n = TMath::Nint(TMath::Power(Npart,alpha));
405     //Int_t n = TMath::Nint(TMath::Power(Ncoll,alpha));
406     Int_t n = alpha * fNpart + (1-alpha) * fNcoll;
407     Int_t ntot=0;
408     //    ntot=n*fNBD->GetRandom();    
409     if (fFastFit == 1) {
410       // NBD
411       ntot = n*hSample->GetRandom();
412     }
413     else if (fFastFit == 2) {
414       // Gaussian
415       Double_t sigma = k*TMath::Sqrt(n*mu);  
416       ntot = gRandom->Gaus(n*mu,sigma);
417     }
418     else {
419       for(Int_t j = 0; j<n; ++j) 
420         ntot+=hSample->GetRandom();
421         //      ntot+=fNBD->GetRandom();
422     }
423     h1->Fill(ntot);
424     
425     if (save) 
426       ntuple->Fill(fNpart,fNcoll,fB,fTaa,ntot);
427    
428   }
429
430   if (save) {
431     ntuple->Write();
432     outFile->Close();
433   }
434   
435   if (fFastFit != 2) delete hSample;
436   return h1;
437   
438 }
439
440 //--------------------------------------------------------------------------------------------------
441 Float_t AliCentralityGlauberFit::CalculateChi2(TH1D *hDATA, TH1D *thistGlau, Float_t eff) 
442 {
443   // note that for different values of neff the mc distribution is identical, just re-scaled below
444   // normalize the two histogram
445   // scale MC up but leave data alone - that preserves error bars !!!!
446   
447   Int_t lowchibin =   hDATA->FindBin(fMultmin);
448   Int_t highchibin =  hDATA->FindBin(fMultmax);
449
450   Double_t mcintegral = thistGlau->Integral(lowchibin,highchibin);
451   Double_t scale = (hDATA->Integral(lowchibin,highchibin)/mcintegral) * ((Double_t)eff);
452   thistGlau->Scale(scale);
453
454   // FIXME: Kolmogorov
455   //return hDATA->KolmogorovTest(thistGlau,"M");
456
457   // calculate the chi2 between MC and real data over some range ????
458   if (fUseChi2) {
459     float chi2 = 0.0;
460     for (int i=lowchibin; i<=highchibin; i++) {
461       if (hDATA->GetBinContent(i) < 1.0) continue;
462       float diff = pow( (thistGlau->GetBinContent(i) - hDATA->GetBinContent(i)) , 2);
463       diff = diff / (pow(hDATA->GetBinError(i) , 2)+pow(thistGlau->GetBinError(i) , 2)); 
464       chi2 += diff;
465     }
466     chi2 = chi2 / (highchibin - lowchibin + 1);
467     return chi2;
468   } else {  // "-2 log likelihood ratio(mu;n) = 2[(mu - n) + n ln(n/mu)]"
469     std::cout << "LL" << std::endl;
470
471     float ll = 0.0;
472     for (int i=lowchibin; i<=highchibin; i++) {
473       Float_t data = hDATA    ->GetBinContent(i);
474       Float_t mc   = thistGlau->GetBinContent(i);
475       Int_t idata = TMath::Nint(data);
476       if (mc < 1.e-9) mc = 1.e-9;
477       Float_t fsub = - mc + idata * TMath::Log(mc);
478       Float_t fobs = 0;
479       //Int_t imc = TMath::Nint(mc);
480       if (idata > 0) {
481         for(Int_t istep = 0; istep < idata; istep++){
482           if (istep > 1) 
483             fobs += TMath::Log(istep);
484         }
485       }
486       fsub -= fobs;
487       ll -=  fsub ;
488     }
489     return 2*ll;
490   }
491 }
492
493 //--------------------------------------------------------------------------------------------------
494 TH1D * AliCentralityGlauberFit::GetTriggerEfficiencyFunction(TH1D *hist1, TH1D *hist2) 
495 {
496   // Get efficiency.
497
498   fhEffi = (TH1D*)hist1->Clone("heffi");
499   fhEffi->Divide(hist2);
500   return fhEffi;
501 }
502
503 //--------------------------------------------------------------------------------------------------
504
505 Float_t AliCentralityGlauberFit::GetTriggerEfficiencyIntegral(TH1D *hist1, TH1D *hist2) 
506 {
507   // Get eff integral.
508   fEffi = hist1->Integral()/hist2->Integral();
509   return fEffi;
510 }
511
512 //--------------------------------------------------------------------------------------------------
513
514 void AliCentralityGlauberFit::SaveHisto(TH1D *hist1, TH1D *hist2, TH1D *heffi, TFile *outrootfile) {
515   outrootfile->cd();
516   hist1->Write();
517   hist2->Write();
518   heffi->Write();
519 }
520
521 //--------------------------------------------------------------------------------------------------
522
523 Double_t AliCentralityGlauberFit::NBD(Int_t n, Double_t mu, Double_t k)
524 {
525   // Compute NBD.
526   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);
527   return ret;
528 }
529
530 //--------------------------------------------------------------------------------------------------
531
532 TH1D *AliCentralityGlauberFit::NBDhist(Double_t mu, Double_t k)
533 {
534   TH1D *h = new TH1D("htemp","",100,-0.5,299.5);
535   h->SetName(Form("nbd_%f_%f",mu,k));
536   h->SetDirectory(0);
537   for (Int_t i=0;i<300;++i) {
538     Double_t val = NBD(i,mu,k);
539     if (val>1e-20) h->Fill(i,val);
540   }
541   return h;
542 }
543
544 void AliCentralityGlauberFit::MinuitFcnNBD(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
545 {
546   // FCN for minuit
547   Double_t alpha = par[0];
548   Double_t mu    = par[1];
549   Double_t k     = par[2];
550   Double_t eff   = par[3];
551   //Double_t eff   = 1;//par[3];
552   if (0) { //avoid warning
553     gin=gin;
554     npar=npar;
555     iflag=iflag;
556   }
557   AliCentralityGlauberFit * obj = (AliCentralityGlauberFit *) gMinuit->GetObjectFit();
558
559   // static TStopwatch sw;
560   // sw.Start();
561   TH1D * thistGlau = obj->GlauberHisto(mu,k,eff,alpha,obj->GetTempHist(),kFALSE);
562   f = obj->CalculateChi2(obj->GetTempHist(),thistGlau,eff);
563   // sw.Stop();
564   // sw.Print();
565   Printf("%f - %f - %f - %f - %f \n",f,alpha,mu,k,eff);
566 }
567
568 Double_t AliCentralityGlauberFit::NBDFunc(Double_t * x, Double_t *par) 
569 {
570   // TF1  interface
571   Double_t mu = par[0];
572   Double_t k  = par[1];
573   Double_t n  = x[0];
574   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);
575   return ret;
576 }