]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG2/EVCHAR/AliCentralityGlauberFit.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Coding conv
[u/mrichter/AliRoot.git] / PWG2 / EVCHAR / AliCentralityGlauberFit.cxx
CommitLineData
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>
36#include "AliCentralityGlauberFit.h"
d20baea6 37#include "TMinuit.h"
38#include "AliLog.h"
39#include "TStopwatch.h"
40#include "TRandom.h"
e6f3f2fe 41ClassImp(AliCentralityGlauberFit)
42
43//______________________________________________________________________________
44
45
d20baea6 46AliCentralityGlauberFit::AliCentralityGlauberFit(const char * filename) : fFastFit(0), fTempHist(0), fGlauberHist(0), fUseChi2(kTRUE) {
e6f3f2fe 47 // standard constructor
48 // glauber file
d20baea6 49 TFile *f = TFile::Open(filename);
50 glau_ntuple = (TNtuple*) f->Get("nt_Pb_Pb");
51 glau_ntuple->SetBranchAddress("Npart",&Npart);
52 glau_ntuple->SetBranchAddress("Ncoll",&Ncoll);
53 glau_ntuple->SetBranchAddress("B",&B);
54 glau_ntuple->SetBranchAddress("tAA",&tAA);
55
56 fNBD = new TF1 ("fNBD", AliCentralityGlauberFit::NBDFunc, 0, 100,2);
57
e6f3f2fe 58}
59
60//--------------------------------------------------------------------------------------------------
61
62AliCentralityGlauberFit::~AliCentralityGlauberFit() {
63 // destructor
64}
65
66//--------------------------------------------------------------------------------------------------
67
68void AliCentralityGlauberFit::AddHisto(TString name) {
69 histnames.push_back(name);
70}
71
72//--------------------------------------------------------------------------------------------------
73
74void AliCentralityGlauberFit::SetOutputFile(TString filename) {
75 outrootfilename = filename;
76}
77
78//--------------------------------------------------------------------------------------------------
79
d20baea6 80void AliCentralityGlauberFit::SetRangeToFit(Float_t fmultmin, Float_t fmultmax)
81{
82 multmin=fmultmin;
83 multmax=fmultmax;
84}
85
86//--------------------------------------------------------------------------------------------------
87
e6f3f2fe 88void AliCentralityGlauberFit::SetGlauberParam(
89 Int_t fNmu,
90 Float_t fmulow,
91 Float_t fmuhigh,
92 Int_t fNk,
93 Float_t fklow,
94 Float_t fkhigh,
95 Int_t fNalpha,
96 Float_t falphalow,
d20baea6 97 Float_t falphahigh,
98 Int_t fNeff,
99 Float_t fefflow,
100 Float_t feffhigh)
e6f3f2fe 101
102{
103 Nmu=fNmu;
104 mulow=fmulow;
105 muhigh=fmuhigh;
106 Nk=fNk;
107 klow=fklow;
108 khigh=fkhigh;
109 Nalpha=fNalpha;
110 alphalow=falphalow;
111 alphahigh=falphahigh;
d20baea6 112 Neff=fNeff;
113 efflow=fefflow;
114 effhigh=feffhigh;
e6f3f2fe 115}
116
117//--------------------------------------------------------------------------------------------------
118
119Float_t AliCentralityGlauberFit::GetPercentileCrossSection() {
120 return percentXsec;
121}
122
d20baea6 123//--------------------------------------------------------------------------------------------------
124void AliCentralityGlauberFit::SetRebin(Int_t frebinFactor) {
125 rebinFactor=frebinFactor;
126}
e6f3f2fe 127//--------------------------------------------------------------------------------------------------
128
129void AliCentralityGlauberFit::MakeFits(TString infilename) {
130 TH1D *hDATA;
d20baea6 131 TH1D *thistGlau;
132 FILE* fTxt = fopen ("parameters.txt","w");
e6f3f2fe 133 TFile *outrootfile;
134
135 // open inrootfile, outrootfile
136 inrootfile = new TFile(infilename);
137 outrootfile = new TFile(outrootfilename,"RECREATE");
138
139 // loop over all distribution names
140 vector<TString>::const_iterator hni;
141 for(hni=histnames.begin(); hni!=histnames.end(); hni++) {
142 hDATA = (TH1D*) (inrootfile->Get(*hni));
d20baea6 143 //TList *list = (TList*) (inrootfile->Get("coutput1"));
144 //hDATA = (TH1D*) (list->FindObject(*hni));
145 hDATA->Rebin(rebinFactor);
e6f3f2fe 146 TH1D *hGLAU = new TH1D("hGLAU","hGLAU",hDATA->GetNbinsX(),0,hDATA->GetNbinsX()*hDATA->GetBinWidth(1));
d20baea6 147 Float_t chi2min = 9999999.0;
e6f3f2fe 148 Float_t alpha_min=-1;
149 Float_t mu_min=-1;
150 Float_t k_min=-1;
151 Float_t eff_min=-1;
152 Float_t alpha, mu, k, eff, chi2;
153
154 for (int nalpha=0;nalpha<Nalpha;nalpha++) {
d20baea6 155 alpha = alphalow + ((float) nalpha ) * (alphahigh - alphalow ) / Nalpha;
e6f3f2fe 156 mu=0.0;
157 for (int nmu=0; nmu<Nmu; nmu++) {
d20baea6 158 mu = (mulow*(1-nalpha*0.05) ) + ((float) nmu ) * (muhigh - mulow ) / Nmu;
159 //mu = mulow + ((float) nmu ) * (muhigh - mulow ) / Nmu;
e6f3f2fe 160
161 for (int nk=0; nk<Nk; nk++) {
162 k = klow + ((float) nk ) * (khigh - klow ) / Nk;
163
d20baea6 164 for (int neff=0; neff<Neff; neff++) {
165 eff = efflow + ((float) neff) * (effhigh - efflow) / Neff;
e6f3f2fe 166
167 thistGlau = GlauberHisto(mu,k,eff,alpha,hDATA,kFALSE);
168 chi2 = CalculateChi2(hDATA,thistGlau,eff);
d20baea6 169 fprintf(fTxt, "%3.3f %3.3f %3.3f %3.3f %3.3f\n",(float) eff, (float) mu, (float) k, (float) alpha, chi2);
170
e6f3f2fe 171 if ( chi2 < chi2min ) {
172 chi2min = chi2;
173 alpha_min=alpha;
174 mu_min=mu;
175 k_min=k;
176 eff_min=eff;
177 }
178 }
179 }
180 }
181 }
182
183 thistGlau = GlauberHisto(mu_min,k_min,eff_min,alpha_min,hDATA,kTRUE);
184 hGLAU = (TH1D *) thistGlau->Clone("hGLAU");
d20baea6 185 //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))));
186 hGLAU->SetName( ((TString)hDATA->GetName()).Append(Form("_GLAU")));
187 hGLAU->SetTitle( ((TString)hDATA->GetName()).Append(Form("_GLAU_%.3f_%.3f_%.3f_%.3f",mu_min,k_min,alpha_min,eff_min)));
188
189 float mcintegral = hGLAU->Integral(hDATA->FindBin(multmin),hGLAU->GetNbinsX());
190 float scale = (hDATA->Integral(hDATA->FindBin(multmin),hDATA->GetNbinsX())/mcintegral) * ((float) eff_min);
191 hGLAU->Scale(scale);
e6f3f2fe 192
193 heffi = GetTriggerEfficiencyFunction(hDATA, hGLAU);
194 SaveHisto(hDATA,hGLAU,heffi,outrootfile);
d20baea6 195 fclose (fTxt);
e6f3f2fe 196
d20baea6 197 cout << "chi2 min is " << chi2min << endl;
198 cout << "fitted " << hGLAU->Integral(hGLAU->FindBin(multmin),hGLAU->FindBin(multmax))/hGLAU->Integral() << " of the total cross section" << endl;
199
200 fTempHist=hDATA;
201 fGlauberHist=hGLAU;
e6f3f2fe 202 }
d20baea6 203
e6f3f2fe 204 // close inrootfile, outrootfile
205 inrootfile->Close();
206 outrootfile->Close();
207
208}
209
210//--------------------------------------------------------------------------------------------------
d20baea6 211void AliCentralityGlauberFit::MakeFitsMinuitNBD(TString infilename) {
212 TH1D *hDATA;
213 TH1D *thistGlau;
214 TFile *outrootfile;
215
216 // open inrootfile, outrootfile
217 inrootfile = new TFile(infilename);
218 outrootfile = new TFile(outrootfilename,"RECREATE");
219
220 // loop over all distribution names
221 vector<TString>::const_iterator hni;
222 for(hni=histnames.begin(); hni!=histnames.end(); hni++) {
223 fFastFit=0;
224 fUseChi2 = 1;
225 hDATA = (TH1D*) (inrootfile->Get(*hni));
226 hDATA->Rebin(rebinFactor);
227 fTempHist=hDATA;
228 TH1D *hGLAU = new TH1D("hGLAU","hGLAU",hDATA->GetNbinsX(),0,hDATA->GetNbinsX()*hDATA->GetBinWidth(1));
229 hGLAU->Sumw2();
230 // Minimize here
231 if(gMinuit) delete gMinuit;
232 new TMinuit(3);
233 gMinuit->mncler();
234 gMinuit->SetFCN(AliCentralityGlauberFit::MinuitFcnNBD);
235 gMinuit->SetObjectFit(this);
236
237 Double_t arglist[2]={0};
238 Int_t ierflg;
239 if (fUseChi2) arglist[0] = 1; // should be 1 if you want to minimize chi2
240 else arglist[0] = 0.5; // should be 0.5 for ll
241 gMinuit->mnexcm("SET ERR",arglist, 1, ierflg);
242
243 // // Change verbosity
244 // gMinuit->Command((TString("SET PRINTOUT ")+long(fVerbosity)).Data());
245
246 // Set parameters
247 // start from the middle of the allowed range, as a guess
248 // gMinuit->mnparm(0,"alpha", (alphalow+alphahigh)/2, 1, alphalow, alphahigh, ierflg);
249 // gMinuit->mnparm(1,"mu" , (mulow +muhigh )/2, 1, mulow , muhigh , ierflg);
250 // gMinuit->mnparm(2,"k" , (klow +khigh )/2, 1, klow , khigh , ierflg);
251 // gMinuit->mnparm(3,"eff" , (efflow +effhigh )/2, 1, efflow , effhigh , ierflg);
252
253 // NBD with limits
254 // gMinuit->mnparm(0,"alpha", 1.104, 0.1, alphalow, alphahigh, ierflg);
255 // gMinuit->mnparm(1,"mu" , 65, 1, mulow , muhigh , ierflg);
256 // gMinuit->mnparm(2,"k" , 0.41, 0.1, klow , khigh , ierflg);
257 // // gMinuit->mnparm(3,"eff" , 1, 1, 0.9 , 1 , ierflg);
258 // gMinuit->mnparm(3,"eff" , 0.95, 0.1, efflow , effhigh , ierflg);
259
260 // // NBD, no limits
261 // gMinuit->mnparm(0,"alpha", 1.104, 0.1, 0,0, ierflg);
262 // gMinuit->mnparm(1,"mu" , 65, 1, 0,0, ierflg);
263 // gMinuit->mnparm(2,"k" , 0.41, 0.1, 0,0, ierflg);
264 // gMinuit->mnparm(3,"eff" , 0.95, 0.1, 0,0 , ierflg);
265
266 if (!fFastFit) {
267 // // NBD, Npart**alpha
268 //gMinuit->mnparm(0,"alpha", 1.16, 0.1, 0,0, ierflg);
269 //gMinuit->mnparm(1,"mu" , 27.55, 1, 0,0, ierflg);
270 //gMinuit->mnparm(2,"k" , 0.64, 0.1, 0,0, ierflg);
271 //gMinuit->mnparm(3,"eff" , 0.97, 0.1, 0,0 , ierflg);
272 // V0 non rescaled
273 // // NBD, alpha Npart * (1-alpha) Ncoll
274 gMinuit->mnparm(0,"alpha", 0.832, 0.1, 0,0, ierflg);
275 gMinuit->mnparm(1,"mu" , 29, 1, 0,0, ierflg);
276 gMinuit->mnparm(2,"k" , 1.4, 0.1, 0,0, ierflg);
277 gMinuit->mnparm(3,"eff" , 0.991, 0.1, 0,0 , ierflg);
278 // SPD
279 //gMinuit->mnparm(0,"alpha", 0.820, 0.1, 0.7,1, ierflg);
280 //gMinuit->mnparm(1,"mu" , 8.267, 1, 1,100, ierflg);
281 //gMinuit->mnparm(2,"k" , 0.530, 0.1, 0,0, ierflg);
282 //gMinuit->mnparm(3,"eff" , 0.990, 0.1, 0,0 , ierflg);
283 } else if (fFastFit == 2) {
284 //gaussian
285 //0.859136, 30.2057, 3.87565, 0.989747
286 gMinuit->mnparm(0,"alpha", 0.59, 0.1, 0,1, ierflg);
287 gMinuit->mnparm(1,"mu" , 30.2, 1, 0,0, ierflg);
288 gMinuit->mnparm(2,"k" , 3.875, 0.1, 0,0, ierflg);
289 // gMinuit->mnparm(3,"eff" , 1, 1, 0.9 , 1 , ierflg);
290 gMinuit->mnparm(3,"eff" , 0.98947, 0.1, 0,0, ierflg);
291 }
e6f3f2fe 292
d20baea6 293 // Call migrad
294 // arglist[0] = 100000; // max calls
295 // arglist[1] = 0.1; // tolerance
296 // arglist[0] = 50; // max calls
297 // arglist[1] = 1; // tolerance
298 // gMinuit->mnexcm("MIGrad",arglist, 2, ierflg);
299 arglist[0] = 100; // max calls
300 arglist[1] = 0.1; // tolerance
301 gMinuit->mnexcm("SIMPLEX",arglist, 2, ierflg);
302 // // fFastFit = 2;
303 arglist[0] = 1000; // max calls
304 arglist[1] = 0.1; // tolerance
305 gMinuit->mnexcm("MIGrad",arglist, 2, ierflg);
306 // gMinuit->mnexcm("IMPROVE",arglist, 0, ierflg);
307
308 if (ierflg != 0) {
309 AliWarning("Abnormal termination of minimization.");
310 // exit(0);
311 }
312
313 // if(opt.Contains("M")) {
314 // gMinuit->mnexcm("IMPROVE",arglist, 0, ierflg);
315 // }
316 // if(opt.Contains("E")) {
317 // gMinuit->mnexcm("HESSE",arglist, 0, ierflg);
318 // gMinuit->mnexcm("MINOS",arglist, 0, ierflg);
319 // }
320
321 // ______________________ Get chi2 and Fit Status __________
322 Double_t amin,edm,errdef;
323 int nvpar, nparx, icstat;
e6f3f2fe 324
d20baea6 325
326 gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
327 Double_t chi2min = amin;
328 cout << "Fit status " << icstat << endl;
329
330 Double_t alpha_min, mu_min, k_min, eff_min;
331 Double_t alpha_mine, mu_mine, k_mine, eff_mine;
332 gMinuit->GetParameter(0, alpha_min , alpha_mine );
333 gMinuit->GetParameter(1, mu_min , mu_mine );
334 gMinuit->GetParameter(2, k_min , k_mine );
335 gMinuit->GetParameter(3, eff_min , eff_mine );
336
337 // Print some infos
338 cout << "chi2 min is " << chi2min << ", " << alpha_min << ", "<< mu_min<< ", " << k_min << ", " << eff_min << endl;
339
340 thistGlau = GlauberHisto(mu_min,k_min,eff_min,alpha_min,hDATA,kTRUE);
341 hGLAU = (TH1D *) thistGlau->Clone("hGLAU");
342 hGLAU->SetName( ((TString)hDATA->GetName()).Append(Form("_GLAU")));
343 hGLAU->SetTitle( ((TString)hDATA->GetName()).Append(Form("_GLAU_%.3f_%.3f_%.3f_%.3f",mu_min,k_min,alpha_min,eff_min)));
344
345
346 cout << "fitted " << hGLAU->Integral(hGLAU->FindBin(multmin),hGLAU->FindBin(multmax))/hGLAU->Integral() << " of the total cross section" << endl;
347 // eff_min=1;
348
349 // Save histogram and ntuple
350 // fFastFit=kFALSE; // Use a better aprox to get the final histo
351
352 // float mcintegral = hGLAU->Integral(hDATA->FindBin(multmin),hGLAU->GetNbinsX());
353 // float scale = (hDATA->Integral(hDATA->FindBin(multmin),hDATA->GetNbinsX())/mcintegral) * ((float) eff_min);
354 // hGLAU->Scale(scale);
355
356 cout << "Chi2 final" << CalculateChi2(hDATA,hGLAU,eff_min) << endl;
357
358
359 heffi = GetTriggerEfficiencyFunction(hDATA, hGLAU);
360 SaveHisto(hDATA,hGLAU,heffi,outrootfile);
361
362 fGlauberHist=hGLAU;
363 }
364
365
366 // close inrootfile, outrootfile
367 inrootfile->Close();
368 outrootfile->Close();
369
370}
371
372
373//--------------------------------------------------------------------------------------------------
374
375TH1D * AliCentralityGlauberFit::GlauberHisto(Float_t mu, Float_t k, Float_t eff, Float_t alpha, TH1D *hDATA, Bool_t save) {
376
377 static TStopwatch sw;
378
379 TH1D *hSample = fFastFit != 2 ? NBDhist(mu,k) : 0;
380 // fNBD->SetParameters(mu,k);
381 static TH1D *h1 = (TH1D*)hDATA->Clone();
382 h1->Reset();
383 // h1->SetName(Form("fit_%.3f_%.3f_%.3f_%.3f",mu,k,eff,alpha));
e6f3f2fe 384 TFile *outFile = NULL;
385 TNtuple *ntuple = NULL;
386
387 if (save) {
388 outFile = TFile::Open("pippo.root", "RECREATE");
d20baea6 389 ntuple = new TNtuple("gnt", "Glauber ntuple", "Npart:Ncoll:B:tAA:ntot");
e6f3f2fe 390 }
391
d20baea6 392
393 Int_t nents = 100000;//glau_ntuple->GetEntries();
394 //if(fFastFit > 0 && fFastFit < 2) nents = 10000;
395 // Int_t nents = 100;//glau_ntuple->GetEntries();
e6f3f2fe 396 for (Int_t i=0;i<nents;++i) {
397 glau_ntuple->GetEntry(i);
d20baea6 398 //Int_t n = TMath::Nint(TMath::Power(Npart,alpha));
399 //Int_t n = TMath::Nint(TMath::Power(Ncoll,alpha));
400 Int_t n = alpha * Npart + (1-alpha) * Ncoll;
e6f3f2fe 401 Int_t ntot=0;
d20baea6 402 // ntot=n*fNBD->GetRandom();
403 if (fFastFit == 1) {
404 // NBD
405 ntot = n*hSample->GetRandom();
406 }
407 else if (fFastFit == 2) {
408 // Gaussian
409 Double_t sigma = k*TMath::Sqrt(n*mu);
410 ntot = gRandom->Gaus(n*mu,sigma);
411 }
412 else {
413 for(Int_t j = 0; j<n; ++j)
414 ntot+=hSample->GetRandom();
415 // ntot+=fNBD->GetRandom();
416 }
e6f3f2fe 417 h1->Fill(ntot);
418
419 if (save)
d20baea6 420 ntuple->Fill(Npart,Ncoll,B,tAA,ntot);
e6f3f2fe 421
422 }
d20baea6 423
e6f3f2fe 424 if (save) {
425 ntuple->Write();
426 outFile->Close();
427 }
428
d20baea6 429 if (fFastFit != 2) delete hSample;
e6f3f2fe 430 return h1;
431
432}
433
434//--------------------------------------------------------------------------------------------------
435
436Float_t AliCentralityGlauberFit::CalculateChi2(TH1D *hDATA, TH1D *thistGlau, Float_t eff) {
437 // note that for different values of neff the mc distribution is identical, just re-scaled below
438 // normalize the two histogram
439 // scale MC up but leave data alone - that preserves error bars !!!!
440
d20baea6 441 int lowchibin = hDATA->FindBin(multmin);
442 int highchibin = hDATA->FindBin(multmax);
443
444 float mcintegral = thistGlau->Integral(lowchibin,mcintegral);
445 // float scale = (hDATA->Integral(lowchibin,highchibin)/mcintegral) * ((float) eff);
446 float scale = (hDATA->Integral(lowchibin,highchibin)/mcintegral) * ((float) eff);
447 thistGlau->Scale(scale);
448
449 // FIXME: Kolmogorov
450 //return hDATA->KolmogorovTest(thistGlau,"M");
451
452 // // calculate the chi2 between MC and real data over some range ????
453 if (fUseChi2) {
454 float chi2 = 0.0;
455 for (int i=lowchibin; i<=highchibin; i++) {
456 if (hDATA->GetBinContent(i) < 1.0) continue;
457 float diff = pow( (thistGlau->GetBinContent(i) - hDATA->GetBinContent(i)) , 2);
458 diff = diff / (pow(hDATA->GetBinError(i) , 2)+pow(thistGlau->GetBinError(i) , 2)); // FIXME squared distance if commented
459 chi2 += diff;
460 }
461 chi2 = chi2 / (highchibin - lowchibin + 1);
462 return chi2;
e6f3f2fe 463 }
d20baea6 464 // "-2 log likelihood ratio(mu;n) = 2[(mu - n) + n ln(n/mu)]"
465 else {
466 cout << "LL" << endl;
467
468
469 float ll = 0.0;
470 for (int i=lowchibin; i<=highchibin; i++) {
471 // if (thistGlau->GetBinContent(i) < 1) continue;
472 // if (hDATA->GetBinContent(i) < 1) continue;
473 // cout << ll << " " << thistGlau->GetBinContent(i) <<" " << hDATA->GetBinContent(i) << endl;
474
475 Float_t data = hDATA ->GetBinContent(i);
476 Float_t mc = thistGlau->GetBinContent(i);
477 Int_t idata = TMath::Nint(data);
478 if (mc < 1.e-9) mc = 1.e-9;
479 Float_t fsub = - mc + idata * TMath::Log(mc);
480 Float_t fobs = 0;
481 Int_t imc = TMath::Nint(mc);
482 if (idata > 0) {
483 for(Int_t istep = 0; istep < idata; istep++){
484 if (istep > 1)
485 fobs += TMath::Log(istep);
486 }
487 }
488 // cout << mc << " " << data << " " << fsub << " " << fobs << endl;
489 fsub -= fobs;
490 ll -= fsub ;
491 }
492 return 2*ll;
e6f3f2fe 493 }
d20baea6 494
e6f3f2fe 495}
496
497//--------------------------------------------------------------------------------------------------
498
499TH1D * AliCentralityGlauberFit::GetTriggerEfficiencyFunction(TH1D *hist1, TH1D *hist2) {
500 heffi = (TH1D*)hist1->Clone("heffi");
501 heffi->Divide(hist2);
502 return heffi;
503}
504
505//--------------------------------------------------------------------------------------------------
506
507Float_t AliCentralityGlauberFit::GetTriggerEfficiencyIntegral(TH1D *hist1, TH1D *hist2) {
508 effi = hist1->Integral()/hist2->Integral();
509 return effi;
510}
511
512//--------------------------------------------------------------------------------------------------
513
514void AliCentralityGlauberFit::SaveHisto(TH1D *hist1, TH1D *hist2, TH1D *heffi, TFile *outrootfile) {
515 outrootfile->cd();
516 hist1->Write();
517 hist2->Write();
518 heffi->Write();
519}
520
521//--------------------------------------------------------------------------------------------------
522
523Double_t AliCentralityGlauberFit::NBD(Int_t n, Double_t mu, Double_t k)
524{
d20baea6 525
e6f3f2fe 526 Double_t mudk = mu/k;
d20baea6 527 //Double_t ret = TMath::Gamma(n+k) / TMath::Gamma(k) / TMath::Factorial(n) *
528 // TMath::Power(mudk,n) / TMath::Power(1+mudk,n+k);
529 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);
e6f3f2fe 530 return ret;
531}
532
533//--------------------------------------------------------------------------------------------------
534
535TH1D *AliCentralityGlauberFit::NBDhist(Double_t mu, Double_t k)
536{
d20baea6 537 TH1D *h = new TH1D("htemp","",100,-0.5,299.5);
e6f3f2fe 538 h->SetName(Form("nbd_%f_%f",mu,k));
539 h->SetDirectory(0);
d20baea6 540 for (Int_t i=0;i<300;++i) {
e6f3f2fe 541 Double_t val = NBD(i,mu,k);
d20baea6 542 if (val>1e-20) h->Fill(i,val);
e6f3f2fe 543 }
544 return h;
545}
546
547
d20baea6 548void AliCentralityGlauberFit::MinuitFcnNBD(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag){
549 // FCN for minuit
550 Double_t alpha = par[0];
551 Double_t mu = par[1];
552 Double_t k = par[2];
553 Double_t eff = par[3];
554 //Double_t eff = 1;//par[3];
555
556 AliCentralityGlauberFit * obj = (AliCentralityGlauberFit *) gMinuit->GetObjectFit();
557
558 // static TStopwatch sw;
559 // sw.Start();
560 TH1D * thistGlau = obj->GlauberHisto(mu,k,eff,alpha,obj->GetTempHist(),kFALSE);
561 f = obj->CalculateChi2(obj->GetTempHist(),thistGlau,eff);
562 // sw.Stop();
563 // sw.Print();
564
565 Printf("%f - %f - %f - %f - %f \n",f,alpha,mu,k,eff);
566
567
568}
569
570Double_t AliCentralityGlauberFit::NBDFunc(Double_t * x, Double_t *par) {
571 // TF1 interface
572
573 Double_t mu = par[0];
574 Double_t k = par[1];
575 Double_t n = x[0];
576 Double_t mudk = mu/k;
577 //Double_t ret = TMath::Gamma(n+k) / TMath::Gamma(k) / TMath::Factorial(n) *
578 // TMath::Power(mudk,n) / TMath::Power(1+mudk,n+k);
579 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);
580
581 return ret;
582
583}