[u/mrichter/AliRoot.git] / PWG2 / EVCHAR / AliCentralityGlauberFit.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/*   Origin: Alberica Toia, CERN, Alberica.Toia@cern.ch                   */

///////////////////////////////////////////////////////////////////////////////
//                                                                           //
//  class to determine centrality percentiles from 1D distributions          // 
//                                                                           //
///////////////////////////////////////////////////////////////////////////////

#include <TNamed.h>
#include <TH1D.h>
#include <TString.h>
#include <TFile.h>
#include <TMath.h>
#include <TROOT.h>
#include <TH2F.h>
#include <TF1.h>
#include <TNtuple.h>
#include <TStyle.h>
#include <TGraphErrors.h>
#include <vector>
#include <TMinuit.h>
#include <TStopwatch.h>
#include <TRandom.h>
#include "AliCentralityGlauberFit.h"
#include "AliLog.h"

ClassImp(AliCentralityGlauberFit)  
 
//______________________________________________________________________________
AliCentralityGlauberFit::AliCentralityGlauberFit(const char *filename) : 
  fNmu(0),    
  fMulow(0),  
  fMuhigh(0), 
  fNk(0),     
  fKlow(0),   
  fKhigh(0),  
  fNalpha(0), 
  fAlphalow(0),   
  fAlphahigh(0),  
  fNeff(0),       
  fEfflow(0),     
  fEffhigh(0),    
  fRebinFactor(0),
  fMultmin(0),    
  fMultmax(0),    
  fGlauntuple(0), 
  fNpart(0),       
  fNcoll(0),       
  fB(0),           
  fTaa(0),         
  fEffi(0),        
  fhEffi(0),      
  fTempHist(0),   
  fGlauberHist(0), 
  fFastFit(0),      
  fNBD(0),         
  fUseChi2(kTRUE),      
  finrootfile(0),  
  foutrootfilename(0),
  fhistnames(), 
  fPercentXsec(0)
{
  // standard constructor
  // glauber file
  TFile *f = TFile::Open(filename);
  fGlauntuple = (TNtuple*) f->Get("nt_Pb_Pb");
  fGlauntuple->SetBranchAddress("Npart",&fNpart);
  fGlauntuple->SetBranchAddress("Ncoll",&fNcoll);
  fGlauntuple->SetBranchAddress("B",&fB);
  fGlauntuple->SetBranchAddress("tAA",&fTaa);
  
  fNBD = new TF1 ("fNBD", AliCentralityGlauberFit::NBDFunc, 0, 100,2);

}

//--------------------------------------------------------------------------------------------------
void AliCentralityGlauberFit::SetRangeToFit(Float_t fmultmin, Float_t fmultmax)
{
  // Set fit range.

  fMultmin=fmultmin;
  fMultmax=fmultmax;
}

//--------------------------------------------------------------------------------------------------
void AliCentralityGlauberFit::SetGlauberParam(
					      Int_t   Nmu,
					      Float_t mulow,
					      Float_t muhigh,
					      Int_t   Nk,
					      Float_t klow,
					      Float_t khigh,
					      Int_t   Nalpha,
					      Float_t alphalow,
					      Float_t alphahigh,
					      Int_t   Neff,
					      Float_t efflow,
					      Float_t effhigh)
{
  // Set Glauber parameters.

  fNmu=Nmu;
  fMulow=mulow;
  fMuhigh=muhigh;
  fNk=Nk;
  fKlow=klow;
  fKhigh=khigh;
  fNalpha=Nalpha;
  fAlphalow=alphalow;
  fAlphahigh=alphahigh;
  fNeff=Neff;
  fEfflow=efflow;
  fEffhigh=effhigh;
}

//--------------------------------------------------------------------------------------------------
void AliCentralityGlauberFit::MakeFits(TString infilename) 
{
  TH1D *hDATA;
  TH1D *thistGlau; 
  FILE* fTxt = fopen ("parameters.txt","w");
  TFile *outrootfile;
  
  // open inrootfile, outrootfile
  finrootfile  = new TFile(infilename);
  outrootfile = new TFile(foutrootfilename,"RECREATE");
  
  // loop over all distribution names  
  std::vector<TString>::const_iterator hni;
  for(hni=fhistnames.begin(); hni!=fhistnames.end(); hni++) {
    hDATA  = (TH1D*) (finrootfile->Get(*hni)); 
    //TList *list  = (TList*) (inrootfile->Get("coutput1")); 
    //hDATA  = (TH1D*) (list->FindObject(*hni)); 
    hDATA->Rebin(fRebinFactor);
    TH1D *hGLAU = new TH1D("hGLAU","hGLAU",hDATA->GetNbinsX(),0,hDATA->GetNbinsX()*hDATA->GetBinWidth(1));
    Float_t chi2min = 9999999.0;
    Float_t alpha_min=-1;
    Float_t mu_min=-1;
    Float_t k_min=-1;
    Float_t eff_min=-1;
    Float_t alpha, mu, k, eff, chi2;   

    for (int nalpha=0;nalpha<fNalpha;nalpha++) {
      alpha = fAlphalow   + ((float) nalpha ) * (fAlphahigh  - fAlphalow ) / fNalpha;
      mu=0.0;
      for (int nmu=0; nmu<fNmu; nmu++) {
	mu = (fMulow*(1-nalpha*0.05) )  + ((float) nmu ) * (fMuhigh  - fMulow ) / fNmu;
	//mu = mulow   + ((float) nmu ) * (muhigh  - mulow ) / Nmu;
	
	for (int nk=0; nk<fNk; nk++) {
	  k = fKlow  + ((float) nk ) * (fKhigh  - fKlow ) / fNk;
	  
	  for (int neff=0; neff<fNeff; neff++) {
	    eff = fEfflow + ((float) neff) * (fEffhigh - fEfflow) / fNeff;
	    
	    thistGlau = GlauberHisto(mu,k,eff,alpha,hDATA,kFALSE);
	    chi2 = CalculateChi2(hDATA,thistGlau,eff);
	    fprintf(fTxt, "%3.3f %3.3f %3.3f %3.3f %3.3f\n",(float) eff, (float) mu, (float) k, (float) alpha, chi2);

	    if ( chi2 < chi2min ) {
	      chi2min = chi2;
	      alpha_min=alpha;
	      mu_min=mu;
	      k_min=k;
	      eff_min=eff;
	    }
	  }
	}
      }
    }

    thistGlau = GlauberHisto(mu_min,k_min,eff_min,alpha_min,hDATA,kTRUE);
    hGLAU = (TH1D *) thistGlau->Clone("hGLAU");
    //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))));
    hGLAU->SetName( ((TString)hDATA->GetName()).Append(Form("_GLAU")));
    hGLAU->SetTitle( ((TString)hDATA->GetName()).Append(Form("_GLAU_%.3f_%.3f_%.3f_%.3f",mu_min,k_min,alpha_min,eff_min)));

    float mcintegral = hGLAU->Integral(1,hGLAU->GetNbinsX());
    float scale = (hDATA->Integral(1,hDATA->GetNbinsX())/mcintegral) * ((float) eff_min);
    hGLAU->Scale(scale);

    fhEffi = GetTriggerEfficiencyFunction(hDATA, hGLAU);
    SaveHisto(hDATA,hGLAU,fhEffi,outrootfile);
    fclose (fTxt);

    std::cout << "chi2 min is " << chi2min << std::endl;
    std::cout << "fitted " << hGLAU->Integral(hGLAU->FindBin(fMultmin),
                                              hGLAU->FindBin(fMultmax))/hGLAU->Integral() 
              << " of the total cross section" << std::endl; 
    fTempHist=hDATA;
    fGlauberHist=hGLAU;
  }

  // close inrootfile, outrootfile
  finrootfile->Close();
  outrootfile->Close();
}

//--------------------------------------------------------------------------------------------------
void AliCentralityGlauberFit::MakeFitsMinuitNBD(TString infilename) 
{
  // Make fits using Minut.

  TH1D *hDATA;
  TH1D *thistGlau; 
  TFile *outrootfile;
  
  // open inrootfile, outrootfile
  finrootfile  = new TFile(infilename);
  outrootfile = new TFile(foutrootfilename,"RECREATE");
  
  // loop over all distribution names  
  std::vector<TString>::const_iterator hni;
  for(hni=fhistnames.begin(); hni!=fhistnames.end(); hni++) {
    fFastFit=0;
    fUseChi2 = 1;
    hDATA  = (TH1D*) (finrootfile->Get(*hni)); 
    hDATA->Rebin(fRebinFactor);
    fTempHist=hDATA;
    TH1D *hGLAU = new TH1D("hGLAU","hGLAU",hDATA->GetNbinsX(),0,hDATA->GetNbinsX()*hDATA->GetBinWidth(1));
    hGLAU->Sumw2();
    // Minimize here
    if(gMinuit) delete gMinuit;
    new TMinuit(3);
    gMinuit->mncler();
    gMinuit->SetFCN(AliCentralityGlauberFit::MinuitFcnNBD);
    gMinuit->SetObjectFit(this); 

    Double_t arglist[2]={0};
    Int_t ierflg;
    if (fUseChi2) arglist[0] = 1; // should be 1 if you want to minimize chi2
    else arglist[0] = 0.5; // should be 0.5 for ll
    gMinuit->mnexcm("SET ERR",arglist, 1, ierflg);

    //Change verbosity
    if (0) {
      //gMinuit->Command((TString("SET PRINTOUT ")+long(fVerbosity)).Data());
    }

    if (!fFastFit) {
      gMinuit->mnparm(0,"alpha", fAlphalow,  0.1,  0.75,1,   ierflg);
      gMinuit->mnparm(1,"mu"   , fMulow,     0.2,  1,100,    ierflg);
      gMinuit->mnparm(2,"k"    , fKlow,      0.1,  0.5,2.5,  ierflg);
      gMinuit->mnparm(3,"eff"  , fEfflow,    0.1,  0.95,1.0, ierflg);      
    } else if (fFastFit == 2) {
      gMinuit->mnparm(0,"alpha", 0.59, 0.1, 0,1, ierflg);
      gMinuit->mnparm(1,"mu"   , 30.2,  1,     0,0, ierflg);
      gMinuit->mnparm(2,"k"    , 3.875, 0.1,  0,0, ierflg);
      gMinuit->mnparm(3,"eff"  , 0.98947, 0.1,  0,0, ierflg);
    }

    // Call migrad
    // arglist[0] = 100000; // max calls
    // arglist[1] = 0.1; // tolerance    
    // arglist[0] = 50; // max calls
    // arglist[1] = 1; // tolerance    
    // gMinuit->mnexcm("MIGrad",arglist, 2, ierflg);
    arglist[0] = 100; // max calls
    arglist[1] = 0.1; // tolerance    
    gMinuit->mnexcm("SIMPLEX",arglist, 2, ierflg);
    arglist[0] = 1000; // max calls
    arglist[1] = 0.1; // tolerance    
    gMinuit->mnexcm("MIGrad",arglist, 2, ierflg);
    //gMinuit->mnexcm("IMPROVE",arglist, 0, ierflg);

    if (ierflg != 0) {
      AliWarning("Abnormal termination of minimization.");
    }
    
    // if(opt.Contains("M")) {
    //   gMinuit->mnexcm("IMPROVE",arglist, 0, ierflg);
    // }
    // if(opt.Contains("E")) {
    // gMinuit->mnexcm("HESSE",arglist, 0, ierflg);
    // gMinuit->mnexcm("MINOS",arglist, 0, ierflg);
    // }

    // ______________________ Get chi2 and Fit Status __________
    Double_t amin,edm,errdef;
    int      nvpar, nparx, icstat;
    
    
    gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
    Double_t chi2min = amin;
    std::cout << "Fit status " << icstat << std::endl;

    Double_t alpha_min, mu_min, k_min, eff_min;
    Double_t alpha_mine, mu_mine, k_mine, eff_mine;
    gMinuit->GetParameter(0, alpha_min , alpha_mine );
    gMinuit->GetParameter(1, mu_min    , mu_mine    );
    gMinuit->GetParameter(2, k_min     , k_mine     );
    gMinuit->GetParameter(3, eff_min   , eff_mine   );

    // Print some infos
    std::cout << "chi2 min is " << chi2min << ", " << alpha_min << ", "<< mu_min<< ", "  << k_min << ", " <<  eff_min << std::endl;

    thistGlau = GlauberHisto(mu_min,k_min,eff_min,alpha_min,hDATA,kTRUE);
    hGLAU = (TH1D *) thistGlau->Clone("hGLAU");
    hGLAU->SetName( ((TString)hDATA->GetName()).Append(Form("_GLAU")));
    hGLAU->SetTitle( ((TString)hDATA->GetName()).Append(Form("_GLAU_%.3f_%.3f_%.3f_%.3f",mu_min,k_min,alpha_min,eff_min)));

    
    std::cout << "fitted " << hGLAU->Integral(hGLAU->FindBin(fMultmin),
                                              hGLAU->FindBin(fMultmax))/hGLAU->Integral() 
              << " of the total cross section" << std::endl; 

    float mcintegral = hGLAU->Integral(1,hGLAU->GetNbinsX());
    float scale = (hDATA->Integral(1,hDATA->GetNbinsX())/mcintegral) * ((float) eff_min);
    hGLAU->Scale(scale);

    std::cout << "Chi2 final " << CalculateChi2(hDATA,hGLAU,eff_min) << std::endl;

    fhEffi = GetTriggerEfficiencyFunction(hDATA, hGLAU);
    SaveHisto(hDATA,hGLAU,fhEffi,outrootfile);
    fGlauberHist=hGLAU;
  }
  // close inrootfile, outrootfile
  finrootfile->Close();
  outrootfile->Close();
}

//--------------------------------------------------------------------------------------------------
TH1D *AliCentralityGlauberFit::GlauberHisto(Float_t mu, Float_t k, Float_t eff, Float_t alpha, 
                                            TH1D *hDATA, Bool_t save) 
{
  // Get Glauber histogram

  eff=eff;
  static TStopwatch sw;

  TH1D *hSample = fFastFit != 2 ? NBDhist(mu,k) : 0;
  //  fNBD->SetParameters(mu,k);
  static TH1D *h1 = (TH1D*)hDATA->Clone();
  h1->Reset();
  //  h1->SetName(Form("fit_%.3f_%.3f_%.3f_%.3f",mu,k,eff,alpha));
  TFile *outFile = NULL;
  TNtuple *ntuple = NULL;
 
  if (save) {
    outFile = TFile::Open("pippo.root", "RECREATE");
    ntuple = new TNtuple("gnt", "Glauber ntuple", "Npart:Ncoll:B:tAA:ntot");
  } 

  Int_t nents = 100000;//glau_ntuple->GetEntries();
  //if(fFastFit > 0 && fFastFit < 2) nents = 10000;
  //  Int_t nents = 100;//glau_ntuple->GetEntries();
  for (Int_t i=0;i<nents;++i) {
    fGlauntuple->GetEntry(i);
    //Int_t n = TMath::Nint(TMath::Power(Npart,alpha));
    //Int_t n = TMath::Nint(TMath::Power(Ncoll,alpha));
    Int_t n = alpha * fNpart + (1-alpha) * fNcoll;
    Int_t ntot=0;
    //    ntot=n*fNBD->GetRandom();    
    if (fFastFit == 1) {
      // NBD
      ntot = n*hSample->GetRandom();
    }
    else if (fFastFit == 2) {
      // Gaussian
      Double_t sigma = k*TMath::Sqrt(n*mu);  
      ntot = gRandom->Gaus(n*mu,sigma);
    }
    else {
      for(Int_t j = 0; j<n; ++j) 
	ntot+=hSample->GetRandom();
	//	ntot+=fNBD->GetRandom();
    }
    h1->Fill(ntot);
    
    if (save) 
      ntuple->Fill(fNpart,fNcoll,fB,fTaa,ntot);
   
  }

  if (save) {
    ntuple->Write();
    outFile->Close();
  }
  
  if (fFastFit != 2) delete hSample;
  return h1;
  
}

//--------------------------------------------------------------------------------------------------
Float_t AliCentralityGlauberFit::CalculateChi2(TH1D *hDATA, TH1D *thistGlau, Float_t eff) 
{
  // note that for different values of neff the mc distribution is identical, just re-scaled below
  // normalize the two histogram
  // scale MC up but leave data alone - that preserves error bars !!!!
  
  int lowchibin =   hDATA->FindBin(fMultmin);
  int highchibin =  hDATA->FindBin(fMultmax);

  float mcintegral = thistGlau->Integral(1,thistGlau->GetNbinsX());
  float scale = (hDATA->Integral(1,hDATA->GetNbinsX())/mcintegral) * ((float) eff);
  thistGlau->Scale(scale);

  // FIXME: Kolmogorov
  //return hDATA->KolmogorovTest(thistGlau,"M");

  // // calculate the chi2 between MC and real data over some range ????
  if (fUseChi2) {
    float chi2 = 0.0;
    for (int i=lowchibin; i<=highchibin; i++) {
      if (hDATA->GetBinContent(i) < 1.0) continue;
      float diff = pow( (thistGlau->GetBinContent(i) - hDATA->GetBinContent(i)) , 2);
      diff = diff / (pow(hDATA->GetBinError(i) , 2)+pow(thistGlau->GetBinError(i) , 2)); // FIXME squared distance if commented
      chi2 += diff;
    }
    chi2 = chi2 / (highchibin - lowchibin + 1);
    return chi2;
  }
  // "-2 log likelihood ratio(mu;n) = 2[(mu - n) + n ln(n/mu)]"
  else {
    std::cout << "LL" << std::endl;
    

    float ll = 0.0;
    for (int i=lowchibin; i<=highchibin; i++) {
      //      if (thistGlau->GetBinContent(i) < 1) continue; 
      //      if (hDATA->GetBinContent(i) < 1) continue; 
      //    cout << ll << " " << thistGlau->GetBinContent(i) <<" " <<  hDATA->GetBinContent(i) << endl;
 
      Float_t data = hDATA    ->GetBinContent(i);
      Float_t mc   = thistGlau->GetBinContent(i);
      Int_t idata = TMath::Nint(data);
      if (mc < 1.e-9) mc = 1.e-9;
      Float_t fsub = - mc + idata * TMath::Log(mc);
      Float_t fobs = 0;
      //Int_t imc = TMath::Nint(mc);
      if (idata > 0) {
	for(Int_t istep = 0; istep < idata; istep++){
	  if (istep > 1) 
	    fobs += TMath::Log(istep);
	}
      }
      //      cout << mc << " " << data << " " << fsub << " " << fobs << endl;
      fsub -= fobs;
      ll -=  fsub ;
    }
    return 2*ll;
  }

}

//--------------------------------------------------------------------------------------------------
TH1D * AliCentralityGlauberFit::GetTriggerEfficiencyFunction(TH1D *hist1, TH1D *hist2) 
{
  // Get efficiency.

  fhEffi = (TH1D*)hist1->Clone("heffi");
  fhEffi->Divide(hist2);
  return fhEffi;
}

//--------------------------------------------------------------------------------------------------
Float_t AliCentralityGlauberFit::GetTriggerEfficiencyIntegral(TH1D *hist1, TH1D *hist2) 
{
  // Get eff integral.
  fEffi = hist1->Integral()/hist2->Integral();
  return fEffi;
}

//--------------------------------------------------------------------------------------------------
void AliCentralityGlauberFit::SaveHisto(TH1D *hist1, TH1D *hist2, TH1D *heffi, TFile *outrootfile) 
{
  outrootfile->cd();
  hist1->Write();
  hist2->Write();
  heffi->Write();
}

//--------------------------------------------------------------------------------------------------
Double_t AliCentralityGlauberFit::NBD(Int_t n, Double_t mu, Double_t k)
{
  // Compute NBD.
  Double_t ret = TMath::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);
  return ret;
}

//--------------------------------------------------------------------------------------------------
TH1D *AliCentralityGlauberFit::NBDhist(Double_t mu, Double_t k)
{
  TH1D *h = new TH1D("htemp","",100,-0.5,299.5);
  h->SetName(Form("nbd_%f_%f",mu,k));
  h->SetDirectory(0);
  for (Int_t i=0;i<300;++i) {
    Double_t val = NBD(i,mu,k);
    if (val>1e-20) h->Fill(i,val);
  }
  return h;
}

//--------------------------------------------------------------------------------------------------
void AliCentralityGlauberFit::MinuitFcnNBD(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
{
  // FCN for minuit
  Double_t alpha = par[0];
  Double_t mu    = par[1];
  Double_t k     = par[2];
  Double_t eff   = par[3];
  //Double_t eff   = 1;//par[3];
  if (0) { //avoid warning
    gin=gin;
    npar=npar;
    iflag=iflag;
  }
  AliCentralityGlauberFit * obj = (AliCentralityGlauberFit *) gMinuit->GetObjectFit();

  // static TStopwatch sw;
  // sw.Start();
  TH1D * thistGlau = obj->GlauberHisto(mu,k,eff,alpha,obj->GetTempHist(),kFALSE);
  f = obj->CalculateChi2(obj->GetTempHist(),thistGlau,eff);
  // sw.Stop();
  // sw.Print();
  Printf("%f - %f - %f - %f - %f \n",f,alpha,mu,k,eff);
}

//--------------------------------------------------------------------------------------------------
Double_t AliCentralityGlauberFit::NBDFunc(Double_t * x, Double_t *par) 
{
  // TF1  interface
  Double_t mu = par[0];
  Double_t k  = par[1];
  Double_t n  = x[0];
  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);
  return ret;
}
Commit	Line	Data
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	42	ClassImp(AliCentralityGlauberFit)
	43
	44	//______________________________________________________________________________
87dea6ac	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	{
e6f3f2fe	78	// standard constructor
e6f3f2fe	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	92	void AliCentralityGlauberFit::SetRangeToFit(Float_t fmultmin, Float_t fmultmax)
d20baea6	93	{
87dea6ac	94	// Set fit range.
	95
	96	fMultmin=fmultmin;
	97	fMultmax=fmultmax;
d20baea6	98	}
	99
	100	//--------------------------------------------------------------------------------------------------
e6f3f2fe	101	void 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	132	void AliCentralityGlauberFit::MakeFits(TString infilename)
87dea6ac	133	{
e6f3f2fe	134	TH1D *hDATA;
d20baea6	135	TH1D *thistGlau;
d20baea6	136	FILE* fTxt = fopen ("parameters.txt","w");
e6f3f2fe	137	TFile *outrootfile;
	138
	139	// open inrootfile, outrootfile
87dea6ac	140	finrootfile = new TFile(infilename);
87dea6ac	141	outrootfile = new TFile(foutrootfilename,"RECREATE");
e6f3f2fe	142
e6f3f2fe	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"));
d20baea6	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++) {
87dea6ac	159	alpha = fAlphalow + ((float) nalpha ) * (fAlphahigh - fAlphalow ) / fNalpha;
e6f3f2fe	160	mu=0.0;
87dea6ac	161	for (int nmu=0; nmu<fNmu; nmu++) {
87dea6ac	162	mu = (fMulow(1-nalpha0.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++) {
87dea6ac	166	k = fKlow + ((float) nk ) * (fKhigh - fKlow ) / fNk;
e6f3f2fe	167
87dea6ac	168	for (int neff=0; neff<fNeff; neff++) {
87dea6ac	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);
d20baea6	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)(100mu_min),(int)(100k_min),(int)(100alpha_min),(int)(100eff_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());
572a3604	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);
87dea6ac	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;
d20baea6	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	215	void 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);
87dea6ac	225	outrootfile = new TFile(foutrootfilename,"RECREATE");
d20baea6	226
d20baea6	227	// loop over all distribution names
87dea6ac	228	std::vector<TString>::const_iterator hni;
87dea6ac	229	for(hni=fhistnames.begin(); hni!=fhistnames.end(); hni++) {
d20baea6	230	fFastFit=0;
d20baea6	231	fUseChi2 = 1;
87dea6ac	232	hDATA = (TH1D) (finrootfile->Get(hni));
87dea6ac	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
d20baea6	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);
d20baea6	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);
87dea6ac	329	SaveHisto(hDATA,hGLAU,fhEffi,outrootfile);
d20baea6	330	fGlauberHist=hGLAU;
d20baea6	331	}
d20baea6	332	// close inrootfile, outrootfile
87dea6ac	333	finrootfile->Close();
d20baea6	334	outrootfile->Close();
d20baea6	335	}
d20baea6	336
d20baea6	337	//--------------------------------------------------------------------------------------------------
87dea6ac	338	TH1D *AliCentralityGlauberFit::GlauberHisto(Float_t mu, Float_t k, Float_t eff, Float_t alpha,
572a3604	339	TH1D *hDATA, Bool_t save)
87dea6ac	340	{
87dea6ac	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));
d20baea6	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 = kTMath::Sqrt(nmu);
	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
e6f3f2fe	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	401	Float_t AliCentralityGlauberFit::CalculateChi2(TH1D hDATA, TH1D thistGlau, Float_t eff)
8dc7e599	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);
572a3604	408	int highchibin = hDATA->FindBin(fMultmax);
d20baea6	409
572a3604	410	float mcintegral = thistGlau->Integral(1,thistGlau->GetNbinsX());
572a3604	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	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;
e6f3f2fe	470	}
	471
	472	//--------------------------------------------------------------------------------------------------
87dea6ac	473	Float_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	481	void AliCentralityGlauberFit::SaveHisto(TH1D hist1, TH1D hist2, TH1D heffi, TFile outrootfile)
572a3604	482	{
e6f3f2fe	483	outrootfile->cd();
	484	hist1->Write();
	485	hist2->Write();
	486	heffi->Write();
	487	}
	488
	489	//--------------------------------------------------------------------------------------------------
e6f3f2fe	490	Double_t AliCentralityGlauberFit::NBD(Int_t n, Double_t mu, Double_t k)
e6f3f2fe	491	{
87dea6ac	492	// Compute NBD.
4cce943f	493	Double_t ret = TMath::Exp( TMath::LnGamma(n+k) - TMath::LnGamma(k) - TMath::LnGamma(n+1) ) *
4cce943f	494	TMath::Power(mu/(mu+k),n) * TMath::Power(1-mu/(mu+k),k);
e6f3f2fe	495	return ret;
	496	}
	497
	498	//--------------------------------------------------------------------------------------------------
e6f3f2fe	499	TH1D *AliCentralityGlauberFit::NBDhist(Double_t mu, Double_t k)
e6f3f2fe	500	{
d20baea6	501	TH1D *h = new TH1D("htemp","",100,-0.5,299.5);
e6f3f2fe	502	h->SetName(Form("nbd_%f_%f",mu,k));
e6f3f2fe	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	512	void AliCentralityGlauberFit::MinuitFcnNBD(Int_t &npar, Double_t gin, Double_t &f, Double_t par, Int_t iflag)
87dea6ac	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	}
d20baea6	535
572a3604	536	//--------------------------------------------------------------------------------------------------
87dea6ac	537	Double_t AliCentralityGlauberFit::NBDFunc(Double_t * x, Double_t *par)
87dea6ac	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	}