PWGPP/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   fRebinFactor(0),
  56   fScalemin(0),
  57   fMultmin(0),
  58   fMultmax(0),
  59   fGlauntuple(0),
  60   fNpart(0),
  61   fNcoll(0),
  62   fB(0),
  63   fTaa(0),
  64   fEffi(0),
  65   fhEffi(0),
  66   fTempHist(0),
  67   fGlauberHist(0),
  68   fFastFit(0),
  69   fAncestor(2),
  70   fNBD(0),
  71   fUseChi2(kTRUE),
  72   fUseAverage(kFALSE),
  73   fhAncestor(0),
  74   fNevents(100000),
  75   fNtrials(100),
  76   fInrootfilename(0),
  77   fInntuplename(0),
  78   fOutrootfilename(0),
  79   fOutntuplename(0),
  80   fAncfilename("ancestor_hists.root"),
  81   fHistnames()
  82 {
  83   // Standard constructor.
  84   TFile *f = 0;
  85   if (filename) {  // glauber file
  86     f = TFile::Open(filename);
  87     fGlauntuple = (TNtuple*) f->Get("nt_Pb_Pb");
  88     fGlauntuple->SetBranchAddress("Npart",&fNpart);
  89     fGlauntuple->SetBranchAddress("Ncoll",&fNcoll);
  90     fGlauntuple->SetBranchAddress("B",&fB);
  91     fGlauntuple->SetBranchAddress("tAA",&fTaa);
  92   }
  93
  94   fNBD = new TF1 ("fNBD", AliCentralityGlauberFit::NBDFunc, 0, 100,2);
  95 }
  96
  97 //--------------------------------------------------------------------------------------------------
  98 void AliCentralityGlauberFit::SetRangeToFit(Double_t fmultmin, Double_t fmultmax)
  99 {
 100   // Set fit range.
 101
 102   fMultmin=fmultmin;
 103   fMultmax=fmultmax;
 104 }
 105
 106 //--------------------------------------------------------------------------------------------------
 107 void AliCentralityGlauberFit::SetRangeToScale(Double_t fscalemin)
 108 {
 109   // Set range where to scale simulated histo to data
 110
 111   fScalemin=fscalemin;
 112 }
 113
 114 //--------------------------------------------------------------------------------------------------
 115 void AliCentralityGlauberFit::SetGlauberParam(
 116                                               Int_t   Nmu,
 117                                               Double_t mulow,
 118                                               Double_t muhigh,
 119                                               Int_t   Nk,
 120                                               Double_t klow,
 121                                               Double_t khigh,
 122                                               Int_t   Nalpha,
 123                                               Double_t alphalow,
 124                                               Double_t alphahigh)
 125 {
 126   // Set Glauber parameters.
 127
 128   fNmu=Nmu;
 129   fMulow=mulow;
 130   fMuhigh=muhigh;
 131   fNk=Nk;
 132   fKlow=klow;
 133   fKhigh=khigh;
 134   fNalpha=Nalpha;
 135   fAlphalow=alphalow;
 136   fAlphahigh=alphahigh;
 137 }
 138
 139 //--------------------------------------------------------------------------------------------------
 140 void AliCentralityGlauberFit::MakeFits()
 141 {
 142   // Make fits.
 143
 144   TH1F *hDATA;
 145   TH1F *thistGlau;
 146   TFile *inrootfile;
 147   TFile *outrootfile;
 148   //FILE* fTxt = fopen ("parameters.txt","w");
 149
 150   // open inrootfile, outrootfile
 151   std::cout << "input file "  << fInrootfilename  << std::endl;
 152   std::cout << "output file " << fOutrootfilename << std::endl;
 153   inrootfile  = new TFile(fInrootfilename,"OPEN");
 154   outrootfile = new TFile(fOutrootfilename,"RECREATE");
 155
 156   // loop over all distribution names
 157   std::vector<TString>::const_iterator hni;
 158   for(hni=fHistnames.begin(); hni!=fHistnames.end(); hni++) {
 159     hDATA  = (TH1F*) (inrootfile->Get(*hni));
 160     if (!hDATA) {
 161       TList *list  = (TList*) (inrootfile->Get("CentralityStat"));
 162       //TList *list  = (TList*) (inrootfile->Get("VZEROEquaFactorStat"));
 163       hDATA  = (TH1F*) (list->FindObject(*hni));
 164     }
 165     hDATA->Rebin(fRebinFactor);
 166     //TH1F *hGLAU = new TH1F("hGLAU","hGLAU",hDATA->GetNbinsX(),0,hDATA->GetNbinsX()*hDATA->GetBinWidth(1));
 167     Double_t chi2min = 9999999.0;
 168     Double_t alpha_min=-1;
 169     Double_t mu_min=-1;
 170     Double_t k_min=-1;
 171     Double_t alpha, mu, k, chi2;
 172
 173     for (Int_t nalpha=0;nalpha<fNalpha;nalpha++) {
 174       alpha = fAlphalow   + ((Double_t) nalpha ) * (fAlphahigh  - fAlphalow ) / fNalpha;
 175
 176       mu=0.0;
 177       for (Int_t nmu=0; nmu<fNmu; nmu++) {
 178         mu = fMulow   + ((Double_t) nmu ) * (fMuhigh  - fMulow ) / fNmu;
 179
 180         for (Int_t nk=0; nk<fNk; nk++) {
 181           k = fKlow  + ((Double_t) nk ) * (fKhigh  - fKlow ) / fNk;
 182
 183           thistGlau = GlauberHisto(mu,k,alpha,hDATA,kFALSE);
 184           chi2 = CalculateChi2(hDATA,thistGlau);
 185           //fprintf(fTxt, "%3.3f %3.3f %3.3f %3.3f\n", (float) mu, (float) k, (float) alpha, chi2);
 186           if ( chi2 < chi2min ) {
 187             chi2min = chi2;
 188             alpha_min=alpha;
 189             mu_min=mu;
 190             k_min=k;
 191           }
 192         }
 193       }
 194     }
 195
 196     thistGlau = GlauberHisto(mu_min,k_min,alpha_min,hDATA,kTRUE);
 197
 198     TH1F * hGLAU = 0x0;
 199     hGLAU = (TH1F *) thistGlau->Clone("hGLAU");
 200     hGLAU->SetName( ((TString)hDATA->GetName()).Append(Form("_GLAU")));
 201     hGLAU->SetTitle( ((TString)hDATA->GetName()).Append(Form("_GLAU_%.3f_%.3f_%.3f",
 202                                                              mu_min,k_min,alpha_min)));
 203
 204     Double_t mcintegral = hGLAU->Integral(hGLAU->FindBin(fScalemin),hGLAU->GetNbinsX());
 205     Double_t scale = (hDATA->Integral(hDATA->FindBin(fScalemin),hDATA->GetNbinsX())/mcintegral);
 206     hGLAU->Scale(scale);
 207
 208     fhEffi = GetTriggerEfficiencyFunction(hDATA, hGLAU);
 209     SaveHisto(hDATA,hGLAU,fhEffi,outrootfile);
 210     //fclose (fTxt);
 211
 212     std::cout << "chi2 min is " << chi2min << std::endl;
 213     std::cout << "fitted " << hGLAU->Integral(hGLAU->FindBin(fMultmin),
 214                                               hGLAU->FindBin(fMultmax))/hGLAU->Integral()
 215               << " of the total cross section" << std::endl;
 216     fTempHist=hDATA;
 217     fGlauberHist=hGLAU;
 218   }
 219
 220   // close inrootfile, outrootfile
 221   inrootfile->Close();
 222   outrootfile->Close();
 223 }
 224
 225 //--------------------------------------------------------------------------------------------------
 226 void AliCentralityGlauberFit::MakeFitsMinuitNBD(Double_t alpha, Double_t mu, Double_t k)
 227 {
 228   // Make fits using Minuit.
 229
 230   if (alpha<0)
 231     alpha = fAlphalow;
 232   if (mu<0)
 233     mu = fMulow;
 234   if (k<0)
 235     k = fKlow;
 236   printf("Calling Minuit with starting values: %f %f %f\n", alpha, mu, k);
 237
 238   TH1F *hDATA;
 239   TH1F *thistGlau;
 240   TFile *inrootfile;
 241   TFile *outrootfile;
 242
 243   // open inrootfile, outrootfile
 244   std::cout << "input file "  << fInrootfilename  << std::endl;
 245   std::cout << "output file " << fOutrootfilename << std::endl;
 246   inrootfile  = new TFile(fInrootfilename,"OPEN");
 247   outrootfile = new TFile(fOutrootfilename,"RECREATE");
 248
 249   // loop over all distribution names
 250   std::vector<TString>::const_iterator hni;
 251   for(hni=fHistnames.begin(); hni!=fHistnames.end(); hni++) {
 252     hDATA  = (TH1F*) (inrootfile->Get(*hni));
 253     if (!hDATA) {
 254       TList *list  = (TList*) (inrootfile->Get("CentralityStat"));
 255       //TList *list  = (TList*) (inrootfile->Get("VZEROEquaFactorStat"));
 256       hDATA  = (TH1F*) (list->FindObject(*hni));
 257     }
 258     hDATA->Rebin(fRebinFactor);
 259
 260     fTempHist=hDATA;
 261     TH1F *hGLAU = new TH1F("hGLAU","hGLAU",hDATA->GetNbinsX(),0,hDATA->GetNbinsX()*hDATA->GetBinWidth(1));
 262     hGLAU->Sumw2();
 263
 264     // Minimize here
 265     if(gMinuit) delete gMinuit;
 266     new TMinuit(3);
 267     gMinuit->mncler();
 268     gMinuit->SetFCN(AliCentralityGlauberFit::MinuitFcnNBD);
 269     gMinuit->SetObjectFit(this);
 270
 271     Double_t arglist[2]={0};
 272     Int_t ierflg;
 273     if (fUseChi2) arglist[0] = 1; // should be 1 if you want to minimize chi2
 274     else arglist[0] = 0.5;        // should be 0.5 for ll
 275     gMinuit->mnexcm("SET ERR",arglist, 1, ierflg);
 276     gMinuit->mnparm(0,"alpha", alpha,  (fAlphahigh-fAlphalow)/fNalpha,  fAlphalow, fAlphahigh, ierflg);
 277     gMinuit->mnparm(1,"mu"   , mu,     (fMuhigh-fMulow)/fNmu, fMulow, fMuhigh, ierflg);
 278     gMinuit->mnparm(2,"k"    , k,      (fKhigh-fKlow)/fNk, fKlow, fKhigh, ierflg);
 279
 280     // Call migrad
 281     arglist[0] = 100; // max calls
 282     arglist[1] = 0.1; // tolerance
 283     gMinuit->mnexcm("SIMPLEX",arglist, 2, ierflg);
 284     arglist[0] = 1000; // max calls
 285     arglist[1] = 0.1;  // tolerance
 286     gMinuit->mnexcm("MIGrad",arglist, 2, ierflg);
 287     //gMinuit->mnexcm("IMPROVE",arglist, 0, ierflg);
 288
 289     if (ierflg != 0) {
 290       AliWarning("Abnormal termination of minimization.");
 291     }
 292
 293     // ______________________ Get chi2 and Fit Status __________
 294     Double_t amin,edm,errdef;
 295     Int_t    nvpar, nparx, icstat;
 296
 297     gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
 298     Double_t chi2min = amin;
 299     std::cout << "Fit status " << icstat << std::endl;
 300
 301     Double_t alpha_min, mu_min, k_min;
 302     Double_t alpha_mine, mu_mine, k_mine;
 303     gMinuit->GetParameter(0, alpha_min , alpha_mine );
 304     gMinuit->GetParameter(1, mu_min    , mu_mine    );
 305     gMinuit->GetParameter(2, k_min     , k_mine     );
 306
 307     // print some infos
 308     std::cout << "chi2 min is " << chi2min << ", " << alpha_min << ", "<< mu_min<< ", "
 309               << k_min << std::endl;
 310
 311     thistGlau = GlauberHisto(mu_min,k_min,alpha_min,hDATA,kTRUE);
 312
 313     hGLAU = (TH1F *) thistGlau->Clone("hGLAU");
 314     hGLAU->SetName( ((TString)hDATA->GetName()).Append(Form("_GLAU")));
 315     hGLAU->SetTitle( ((TString)hDATA->GetName()).Append(Form("_GLAU_%.3f_%.3f_%.3f",
 316                                                              mu_min,k_min,alpha_min)));
 317
 318
 319     std::cout << "fitted " << hGLAU->Integral(hGLAU->FindBin(fMultmin),
 320                                               hGLAU->FindBin(fMultmax))/hGLAU->Integral()
 321               << " of the total cross section" << std::endl;
 322
 323     Double_t mcintegral = hGLAU->Integral(hGLAU->FindBin(fScalemin),hGLAU->GetNbinsX());
 324     Double_t scale = (hDATA->Integral(hDATA->FindBin(fScalemin),hDATA->GetNbinsX())/mcintegral);
 325     hGLAU->Scale(scale);
 326
 327     std::cout << "Chi2 final " << CalculateChi2(hDATA,hGLAU) << std::endl;
 328
 329     fhEffi = GetTriggerEfficiencyFunction(hDATA, hGLAU);
 330     SaveHisto(hDATA,hGLAU,fhEffi,outrootfile);
 331     fGlauberHist=hGLAU;
 332   }
 333   // close inrootfile, outrootfile
 334   inrootfile->Close();
 335   outrootfile->Close();
 336 }
 337
 338 //--------------------------------------------------------------------------------------------------
 339 TH1F *AliCentralityGlauberFit::GlauberHisto(Double_t mu, Double_t k, Double_t alpha,
 340                                             TH1F *hDATA, Bool_t save)
 341 {
 342   // Get Glauber histogram.
 343   static TH1F *h1 = (TH1F*)hDATA->Clone();
 344   h1->Reset();
 345   h1->SetName(Form("fit_%.3f_%.3f_%.3f",mu,k,alpha));
 346
 347   if (fUseAverage) {
 348     fhAncestor = MakeAncestor(alpha);
 349     for (Int_t np=1; np<=fhAncestor->GetNbinsX(); ++np) {
 350       Double_t nanc = fhAncestor->GetBinCenter(np);
 351       Double_t weights = fhAncestor->GetBinContent(np);
 352
 353       if (weights <= 0) continue;
 354       Int_t trials = (Int_t) (20 * nanc * (Int_t) mu);
 355       if (trials <=0) continue;
 356       for (Int_t j=0; j<trials; j++) {
 357         Double_t nbdvalue = NBD(j, mu * nanc, k * nanc);
 358         h1->Fill((Double_t) j, nbdvalue * weights);
 359       }
 360     }
 361     return h1;
 362   }
 363
 364   TH1F *hSample = fFastFit != 2 ? NBDhist(mu,k) : 0;
 365
 366   TFile *outFile = NULL;
 367   TNtuple *ntuple = NULL;
 368
 369   if (save) {
 370     outFile = new TFile(fOutntuplename,"RECREATE");
 371     ntuple = new TNtuple("gnt", "Glauber ntuple", "Npart:Ncoll:B:tAA:ntot");
 372   }
 373
 374   Int_t nents = 0;
 375   if (fGlauntuple)
 376     nents = fGlauntuple->GetEntries();
 377
 378   for (Int_t i=0;i<fNevents;++i) {
 379     if (fGlauntuple)
 380       fGlauntuple->GetEntry(i % nents);
 381     else {
 382       fNpart = 2;
 383       fNcoll = 1;
 384     }
 385     Int_t n=0;
 386     if (fAncestor == 1)      n = (Int_t)(TMath::Power(fNpart,alpha));
 387     //if (fAncestor == 1)      n = (Int_t)(TMath::Power(fNcoll,alpha));
 388     else if (fAncestor == 2) n = (Int_t)(alpha * fNpart + (1-alpha) * fNcoll);
 389     else if (fAncestor == 3) n = (Int_t)((1-alpha) * fNpart/2 + alpha * fNcoll);
 390
 391     Int_t ntot=0;
 392     if (fFastFit == 1) {
 393       ntot = (Int_t)(n*hSample->GetRandom()); // NBD
 394     }
 395     else if (fFastFit == 2) {
 396       Double_t sigma = k*TMath::Sqrt(n*mu);
 397       ntot = (Int_t)(gRandom->Gaus(n*mu,sigma)); // Gaussian
 398     }
 399     else {
 400       for(Int_t j = 0; j<(Int_t)n; ++j)
 401         ntot += (Int_t)hSample->GetRandom();
 402     }
 403     h1->Fill(ntot);
 404
 405     if (save)
 406       ntuple->Fill(fNpart,fNcoll,fB,fTaa,ntot);
 407   }
 408
 409   if (save) {
 410     ntuple->Write();
 411     outFile->Close();
 412   }
 413
 414   if (fFastFit != 2) delete hSample;
 415   return h1;
 416 }
 417
 418 //--------------------------------------------------------------------------------------------------
 419 Double_t AliCentralityGlauberFit::CalculateChi2(TH1F *hDATA, TH1F *thistGlau)
 420 {
 421   // Note that for different values of neff the mc distribution is identical, just re-scaled below
 422   // normalize the two histogram, scale MC up but leave data alone - that preserves error bars !!!!
 423
 424   Int_t lowchibin =   hDATA->FindBin(fMultmin);
 425   Int_t highchibin =  hDATA->FindBin(fMultmax);
 426
 427   Double_t mcintegral = thistGlau->Integral(1,thistGlau->GetNbinsX());
 428   Double_t scale = (hDATA->Integral(1,hDATA->GetNbinsX())/mcintegral);
 429   thistGlau->Scale(scale);
 430
 431   // calculate the chi2 between MC and real data over some range ????
 432   if (fUseChi2) {
 433     Double_t chi2 = 0.0;
 434     for (Int_t i=lowchibin; i<=highchibin; i++) {
 435       if (hDATA->GetBinContent(i) < 1.0) continue;
 436       Double_t diff = TMath::Power((thistGlau->GetBinContent(i) - hDATA->GetBinContent(i)),2);
 437       diff = diff / (TMath::Power(hDATA->GetBinError(i),2) + TMath::Power(thistGlau->GetBinError(i),2));
 438       chi2 += diff;
 439     }
 440     chi2 = chi2 / (highchibin - lowchibin + 1);
 441     return chi2;
 442   }
 443   // "-2 log likelihood ratio(mu;n) = 2[(mu - n) + n ln(n/mu)]"
 444   else {
 445     std::cout << "LL" << std::endl;
 446
 447     Double_t ll = 0.0;
 448     for (Int_t i=lowchibin; i<=highchibin; i++) {
 449       Double_t data = hDATA    ->GetBinContent(i);
 450       Double_t mc   = thistGlau->GetBinContent(i);
 451       Int_t idata = TMath::Nint(data);
 452       if (mc < 1.e-9) mc = 1.e-9;
 453       Double_t fsub = - mc + idata * TMath::Log(mc);
 454       Double_t fobs = 0;
 455       if (idata > 0) {
 456         for(Int_t istep = 0; istep < idata; istep++){
 457           if (istep > 1)
 458             fobs += TMath::Log(istep);
 459         }
 460       }
 461       fsub -= fobs;
 462       ll -=  fsub ;
 463     }
 464     return 2*ll;
 465   }
 466 }
 467
 468 //--------------------------------------------------------------------------------------------------
 469 TH1F *AliCentralityGlauberFit::GetTriggerEfficiencyFunction(TH1F *hist1, TH1F *hist2)
 470 {
 471   // Get efficiency.
 472
 473   fhEffi = (TH1F*)hist1->Clone("heffi");
 474   fhEffi->Divide(hist2);
 475   return fhEffi;
 476 }
 477
 478 //--------------------------------------------------------------------------------------------------
 479 Double_t AliCentralityGlauberFit::GetTriggerEfficiencyIntegral(TH1F *hist1, TH1F *hist2)
 480 {
 481   // Get eff integral.
 482
 483   fEffi = hist1->Integral()/hist2->Integral();
 484   return fEffi;
 485 }
 486
 487 //--------------------------------------------------------------------------------------------------
 488 void AliCentralityGlauberFit::SaveHisto(TH1F *hist1, TH1F *hist2, TH1F *heffi, TFile *outrootfile)
 489 {
 490   // Save histograms.
 491
 492   outrootfile->cd();
 493   hist1->Write();
 494   hist2->Write();
 495   heffi->Write();
 496 }
 497
 498 //--------------------------------------------------------------------------------------------------
 499 Double_t AliCentralityGlauberFit::NBD(Int_t n, Double_t mu, Double_t k) const
 500 {
 501   // Compute NBD.
 502   Double_t F;
 503   Double_t f;
 504
 505   if (n+k > 100.0) {
 506     // log method for handling large numbers
 507     F  = TMath::LnGamma(n + k)- TMath::LnGamma(n + 1.)- TMath::LnGamma(k);
 508     f  = n * TMath::Log(mu/k) - (n + k) * TMath::Log(1.0 + mu/k);
 509     F = F+f;
 510     F = TMath::Exp(F);
 511   } else {
 512     F  = TMath::Gamma(n + k) / ( TMath::Gamma(n + 1.) * TMath::Gamma(k) );
 513     f  = n * TMath::Log(mu/k) - (n + k) * TMath::Log(1.0 + mu/k);
 514     f  = TMath::Exp(f);
 515     F *= f;
 516   }
 517
 518   return F;
 519
 520 }
 521
 522 //--------------------------------------------------------------------------------------------------
 523 TH1F *AliCentralityGlauberFit::NBDhist(Double_t mu, Double_t k)
 524 {
 525   // Interface for TH1F.
 526
 527   TH1F *h = new TH1F("htemp","",100,-0.5,299.5);
 528   h->SetName(Form("nbd_%f_%f",mu,k));
 529   h->SetDirectory(0);
 530   for (Int_t i=0;i<300;++i) {
 531     Double_t val = NBD(i,mu,k);
 532     if (val>1e-20) h->Fill(i,val);
 533   }
 534   return h;
 535 }
 536
 537 //--------------------------------------------------------------------------------------------------
 538 void AliCentralityGlauberFit::MinuitFcnNBD(Int_t &/*npar*/, Double_t */*gin*/, Double_t &f, Double_t *par, Int_t /*iflag*/)
 539 {
 540   // FCN for minuit.
 541
 542   Double_t alpha = par[0];
 543   Double_t mu    = par[1];
 544   Double_t k     = par[2];
 545
 546   AliCentralityGlauberFit * obj = (AliCentralityGlauberFit *) gMinuit->GetObjectFit();
 547   TH1F * thistGlau = obj->GlauberHisto(mu,k,alpha,obj->GetTempHist(),kFALSE);
 548   f = obj->CalculateChi2(obj->GetTempHist(),thistGlau);
 549   Printf("Minuit step: chi2=%f, alpha=%f,  mu=%f, k=%f\n",f,alpha,mu,k);
 550 }
 551
 552 //--------------------------------------------------------------------------------------------------
 553 Double_t AliCentralityGlauberFit::NBDFunc(const Double_t *x, const Double_t *par)
 554 {
 555   // TF1  interface.
 556
 557   Double_t mu = par[0];
 558   Double_t k  = par[1];
 559   Double_t n  = x[0];
 560   Double_t ret = exp( TMath::LnGamma(n+k) - TMath::LnGamma(k) - TMath::LnGamma(n+1) ) *
 561     TMath::Power(mu/(mu+k),n) * TMath::Power(1-mu/(mu+k),k);
 562   return ret;
 563 }
 564
 565 //--------------------------------------------------------------------------------------------------
 566 TH1F *AliCentralityGlauberFit::MakeAncestor(Double_t alpha)
 567 {
 568   // Make ancestor histogram.
 569
 570   TString hname(Form("fhAncestor_%.3f",alpha));
 571   if (fhAncestor) {
 572     if (hname.CompareTo(fhAncestor->GetName())==0)
 573       return fhAncestor;
 574   }
 575
 576   delete fhAncestor;
 577   fhAncestor = 0;
 578   TFile *ancfile = TFile::Open(fAncfilename,"read");
 579   if (ancfile && ancfile->IsOpen()) {
 580     fhAncestor = dynamic_cast<TH1F*>(ancfile->Get(hname));
 581     if (fhAncestor) {
 582       fhAncestor->SetDirectory(0);
 583       delete ancfile;
 584       return fhAncestor;
 585     }
 586   }
 587   delete ancfile;
 588
 589   fhAncestor = new TH1F(hname,hname,3000,0,3000);
 590   fhAncestor->SetDirectory(0);
 591   Int_t nents = fGlauntuple->GetEntries();
 592   for (Int_t i=0;i<nents;++i) {
 593     fGlauntuple->GetEntry(i % nents);
 594     Int_t n=0;
 595     if (fAncestor == 1)    n = (Int_t) (TMath::Power(fNpart,alpha));
 596     //if (fAncestor == 1)      n = (Int_t) (TMath::Power(fNcoll,alpha));
 597     else if (fAncestor == 2) n = (Int_t) (alpha * fNpart + (1-alpha) * fNcoll);
 598     else if (fAncestor == 3) n = (Int_t) ((1-alpha) * fNpart/2 + alpha * fNcoll);
 599     fhAncestor->Fill(n);
 600   }
 601
 602   ancfile = TFile::Open(fAncfilename,"update");
 603   if (ancfile && ancfile->IsOpen()) {
 604     fhAncestor->Write();
 605   }
 606   delete ancfile;
 607   return fhAncestor;
 608 }