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(1,hGLAU->GetNbinsX());
 194     float scale = (hDATA->Integral(1,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     if (0) {
 252       //gMinuit->Command((TString("SET PRINTOUT ")+long(fVerbosity)).Data());
 253     }
 254
 255     if (!fFastFit) {
 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);
 260     } else if (fFastFit == 2) {
 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);
 264       gMinuit->mnparm(3,"eff"  , 0.98947, 0.1,  0,0, ierflg);
 265     }
 266
 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);
 276     arglist[0] = 1000; // max calls
 277     arglist[1] = 0.1; // tolerance
 278     gMinuit->mnexcm("MIGrad",arglist, 2, ierflg);
 279     //gMinuit->mnexcm("IMPROVE",arglist, 0, ierflg);
 280
 281     if (ierflg != 0) {
 282       AliWarning("Abnormal termination of minimization.");
 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;
 296
 297
 298     gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
 299     Double_t chi2min = amin;
 300     std::cout << "Fit status " << icstat << std::endl;
 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
 310     std::cout << "chi2 min is " << chi2min << ", " << alpha_min << ", "<< mu_min<< ", "  << k_min << ", " <<  eff_min << std::endl;
 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
 318     std::cout << "fitted " << hGLAU->Integral(hGLAU->FindBin(fMultmin),
 319                                               hGLAU->FindBin(fMultmax))/hGLAU->Integral()
 320               << " of the total cross section" << std::endl;
 321
 322     float mcintegral = hGLAU->Integral(1,hGLAU->GetNbinsX());
 323     float scale = (hDATA->Integral(1,hDATA->GetNbinsX())/mcintegral) * ((float) eff_min);
 324     hGLAU->Scale(scale);
 325
 326     std::cout << "Chi2 final " << CalculateChi2(hDATA,hGLAU,eff_min) << std::endl;
 327
 328     fhEffi = GetTriggerEfficiencyFunction(hDATA, hGLAU);
 329     SaveHisto(hDATA,hGLAU,fhEffi,outrootfile);
 330     fGlauberHist=hGLAU;
 331   }
 332   // close inrootfile, outrootfile
 333   finrootfile->Close();
 334   outrootfile->Close();
 335 }
 336
 337 //--------------------------------------------------------------------------------------------------
 338 TH1D *AliCentralityGlauberFit::GlauberHisto(Float_t mu, Float_t k, Float_t eff, Float_t alpha,
 339                                             TH1D *hDATA, Bool_t save)
 340 {
 341   // Get Glauber histogram
 342
 343   eff=eff;
 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));
 351   TFile *outFile = NULL;
 352   TNtuple *ntuple = NULL;
 353
 354   if (save) {
 355     outFile = TFile::Open("pippo.root", "RECREATE");
 356     ntuple = new TNtuple("gnt", "Glauber ntuple", "Npart:Ncoll:B:tAA:ntot");
 357   }
 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();
 362   for (Int_t i=0;i<nents;++i) {
 363     fGlauntuple->GetEntry(i);
 364     //Int_t n = TMath::Nint(TMath::Power(Npart,alpha));
 365     //Int_t n = TMath::Nint(TMath::Power(Ncoll,alpha));
 366     Int_t n = alpha * fNpart + (1-alpha) * fNcoll;
 367     Int_t ntot=0;
 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     }
 383     h1->Fill(ntot);
 384
 385     if (save)
 386       ntuple->Fill(fNpart,fNcoll,fB,fTaa,ntot);
 387
 388   }
 389
 390   if (save) {
 391     ntuple->Write();
 392     outFile->Close();
 393   }
 394
 395   if (fFastFit != 2) delete hSample;
 396   return h1;
 397
 398 }
 399
 400 //--------------------------------------------------------------------------------------------------
 401 Float_t AliCentralityGlauberFit::CalculateChi2(TH1D *hDATA, TH1D *thistGlau, Float_t eff)
 402 {
 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
 407   int lowchibin =   hDATA->FindBin(fMultmin);
 408   int highchibin =  hDATA->FindBin(fMultmax);
 409
 410   float mcintegral = thistGlau->Integral(1,thistGlau->GetNbinsX());
 411   float scale = (hDATA->Integral(1,hDATA->GetNbinsX())/mcintegral) * ((float) eff);
 412   thistGlau->Scale(scale);
 413
 414   // FIXME: Kolmogorov
 415   //return hDATA->KolmogorovTest(thistGlau,"M");
 416
 417   // // calculate the chi2 between MC and real data over some range ????
 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);
 423       diff = diff / (pow(hDATA->GetBinError(i) , 2)+pow(thistGlau->GetBinError(i) , 2)); // FIXME squared distance if commented
 424       chi2 += diff;
 425     }
 426     chi2 = chi2 / (highchibin - lowchibin + 1);
 427     return chi2;
 428   }
 429   // "-2 log likelihood ratio(mu;n) = 2[(mu - n) + n ln(n/mu)]"
 430   else {
 431     std::cout << "LL" << std::endl;
 432
 433
 434     float ll = 0.0;
 435     for (int i=lowchibin; i<=highchibin; i++) {
 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
 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;
 446       //Int_t imc = TMath::Nint(mc);
 447       if (idata > 0) {
 448         for(Int_t istep = 0; istep < idata; istep++){
 449           if (istep > 1)
 450             fobs += TMath::Log(istep);
 451         }
 452       }
 453       //      cout << mc << " " << data << " " << fsub << " " << fobs << endl;
 454       fsub -= fobs;
 455       ll -=  fsub ;
 456     }
 457     return 2*ll;
 458   }
 459
 460 }
 461
 462 //--------------------------------------------------------------------------------------------------
 463 TH1D * AliCentralityGlauberFit::GetTriggerEfficiencyFunction(TH1D *hist1, TH1D *hist2)
 464 {
 465   // Get efficiency.
 466
 467   fhEffi = (TH1D*)hist1->Clone("heffi");
 468   fhEffi->Divide(hist2);
 469   return fhEffi;
 470 }
 471
 472 //--------------------------------------------------------------------------------------------------
 473 Float_t AliCentralityGlauberFit::GetTriggerEfficiencyIntegral(TH1D *hist1, TH1D *hist2)
 474 {
 475   // Get eff integral.
 476   fEffi = hist1->Integral()/hist2->Integral();
 477   return fEffi;
 478 }
 479
 480 //--------------------------------------------------------------------------------------------------
 481 void AliCentralityGlauberFit::SaveHisto(TH1D *hist1, TH1D *hist2, TH1D *heffi, TFile *outrootfile)
 482 {
 483   outrootfile->cd();
 484   hist1->Write();
 485   hist2->Write();
 486   heffi->Write();
 487 }
 488
 489 //--------------------------------------------------------------------------------------------------
 490 Double_t AliCentralityGlauberFit::NBD(Int_t n, Double_t mu, Double_t k)
 491 {
 492   // Compute NBD.
 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);
 495   return ret;
 496 }
 497
 498 //--------------------------------------------------------------------------------------------------
 499 TH1D *AliCentralityGlauberFit::NBDhist(Double_t mu, Double_t k)
 500 {
 501   TH1D *h = new TH1D("htemp","",100,-0.5,299.5);
 502   h->SetName(Form("nbd_%f_%f",mu,k));
 503   h->SetDirectory(0);
 504   for (Int_t i=0;i<300;++i) {
 505     Double_t val = NBD(i,mu,k);
 506     if (val>1e-20) h->Fill(i,val);
 507   }
 508   return h;
 509 }
 510
 511 //--------------------------------------------------------------------------------------------------
 512 void AliCentralityGlauberFit::MinuitFcnNBD(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
 513 {
 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];
 520   if (0) { //avoid warning
 521     gin=gin;
 522     npar=npar;
 523     iflag=iflag;
 524   }
 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();
 533   Printf("%f - %f - %f - %f - %f \n",f,alpha,mu,k,eff);
 534 }
 535
 536 //--------------------------------------------------------------------------------------------------
 537 Double_t AliCentralityGlauberFit::NBDFunc(Double_t * x, Double_t *par)
 538 {
 539   // TF1  interface
 540   Double_t mu = par[0];
 541   Double_t k  = par[1];
 542   Double_t n  = x[0];
 543   Double_t ret = exp( TMath::LnGamma(n+k) - TMath::LnGamma(k) - TMath::LnGamma(n+1) ) * TMath::Power(mu/(mu+k),n) * TMath::Power(1-mu/(mu+k),k);
 544   return ret;
 545 }