Remove lgamma and replace by TMath::LnGamma which should work across platforms.
[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(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( lgamma(n+k) - lgamma(k) - lgamma(n+1) ) * TMath::Power(mu/(mu+k),n) * TMath::Power(1-mu/(mu+k),k);
544   return ret;
545 }