]>
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 |
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) |
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) |
133 | { | |
e6f3f2fe | 134 | TH1D *hDATA; |
d20baea6 | 135 | TH1D *thistGlau; |
136 | FILE* fTxt = fopen ("parameters.txt","w"); | |
e6f3f2fe | 137 | TFile *outrootfile; |
138 | ||
139 | // open inrootfile, outrootfile | |
87dea6ac | 140 | finrootfile = new TFile(infilename); |
141 | outrootfile = new TFile(foutrootfilename,"RECREATE"); | |
e6f3f2fe | 142 | |
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")); |
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++) { |
159 | alpha = fAlphalow + ((float) nalpha ) * (fAlphahigh - fAlphalow ) / fNalpha; | |
e6f3f2fe | 160 | mu=0.0; |
87dea6ac | 161 | for (int nmu=0; nmu<fNmu; nmu++) { |
162 | mu = (fMulow*(1-nalpha*0.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++) { |
166 | k = fKlow + ((float) nk ) * (fKhigh - fKlow ) / fNk; | |
e6f3f2fe | 167 | |
87dea6ac | 168 | for (int neff=0; neff<fNeff; neff++) { |
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); |
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)(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 | ||
87dea6ac | 193 | float mcintegral = hGLAU->Integral(hDATA->FindBin(fMultmin),hGLAU->GetNbinsX()); |
194 | float scale = (hDATA->Integral(hDATA->FindBin(fMultmin),hDATA->GetNbinsX())/mcintegral) * ((float) eff_min); | |
d20baea6 | 195 | hGLAU->Scale(scale); |
e6f3f2fe | 196 | |
87dea6ac | 197 | fhEffi = GetTriggerEfficiencyFunction(hDATA, hGLAU); |
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; |
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); |
225 | outrootfile = new TFile(foutrootfilename,"RECREATE"); | |
d20baea6 | 226 | |
227 | // loop over all distribution names | |
87dea6ac | 228 | std::vector<TString>::const_iterator hni; |
229 | for(hni=fhistnames.begin(); hni!=fhistnames.end(); hni++) { | |
d20baea6 | 230 | fFastFit=0; |
231 | fUseChi2 = 1; | |
87dea6ac | 232 | hDATA = (TH1D*) (finrootfile->Get(*hni)); |
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 | ||
250 | // // Change verbosity | |
251 | // gMinuit->Command((TString("SET PRINTOUT ")+long(fVerbosity)).Data()); | |
252 | ||
253 | // Set parameters | |
254 | // start from the middle of the allowed range, as a guess | |
255 | // gMinuit->mnparm(0,"alpha", (alphalow+alphahigh)/2, 1, alphalow, alphahigh, ierflg); | |
256 | // gMinuit->mnparm(1,"mu" , (mulow +muhigh )/2, 1, mulow , muhigh , ierflg); | |
257 | // gMinuit->mnparm(2,"k" , (klow +khigh )/2, 1, klow , khigh , ierflg); | |
258 | // gMinuit->mnparm(3,"eff" , (efflow +effhigh )/2, 1, efflow , effhigh , ierflg); | |
259 | ||
260 | // NBD with limits | |
261 | // gMinuit->mnparm(0,"alpha", 1.104, 0.1, alphalow, alphahigh, ierflg); | |
262 | // gMinuit->mnparm(1,"mu" , 65, 1, mulow , muhigh , ierflg); | |
263 | // gMinuit->mnparm(2,"k" , 0.41, 0.1, klow , khigh , ierflg); | |
264 | // // gMinuit->mnparm(3,"eff" , 1, 1, 0.9 , 1 , ierflg); | |
265 | // gMinuit->mnparm(3,"eff" , 0.95, 0.1, efflow , effhigh , ierflg); | |
266 | ||
267 | // // NBD, no limits | |
268 | // gMinuit->mnparm(0,"alpha", 1.104, 0.1, 0,0, ierflg); | |
269 | // gMinuit->mnparm(1,"mu" , 65, 1, 0,0, ierflg); | |
270 | // gMinuit->mnparm(2,"k" , 0.41, 0.1, 0,0, ierflg); | |
271 | // gMinuit->mnparm(3,"eff" , 0.95, 0.1, 0,0 , ierflg); | |
272 | ||
273 | if (!fFastFit) { | |
274 | // // NBD, Npart**alpha | |
275 | //gMinuit->mnparm(0,"alpha", 1.16, 0.1, 0,0, ierflg); | |
276 | //gMinuit->mnparm(1,"mu" , 27.55, 1, 0,0, ierflg); | |
277 | //gMinuit->mnparm(2,"k" , 0.64, 0.1, 0,0, ierflg); | |
278 | //gMinuit->mnparm(3,"eff" , 0.97, 0.1, 0,0 , ierflg); | |
279 | // V0 non rescaled | |
280 | // // NBD, alpha Npart * (1-alpha) Ncoll | |
281 | gMinuit->mnparm(0,"alpha", 0.832, 0.1, 0,0, ierflg); | |
282 | gMinuit->mnparm(1,"mu" , 29, 1, 0,0, ierflg); | |
283 | gMinuit->mnparm(2,"k" , 1.4, 0.1, 0,0, ierflg); | |
284 | gMinuit->mnparm(3,"eff" , 0.991, 0.1, 0,0 , ierflg); | |
285 | // SPD | |
286 | //gMinuit->mnparm(0,"alpha", 0.820, 0.1, 0.7,1, ierflg); | |
287 | //gMinuit->mnparm(1,"mu" , 8.267, 1, 1,100, ierflg); | |
288 | //gMinuit->mnparm(2,"k" , 0.530, 0.1, 0,0, ierflg); | |
289 | //gMinuit->mnparm(3,"eff" , 0.990, 0.1, 0,0 , ierflg); | |
290 | } else if (fFastFit == 2) { | |
291 | //gaussian | |
292 | //0.859136, 30.2057, 3.87565, 0.989747 | |
293 | gMinuit->mnparm(0,"alpha", 0.59, 0.1, 0,1, ierflg); | |
294 | gMinuit->mnparm(1,"mu" , 30.2, 1, 0,0, ierflg); | |
295 | gMinuit->mnparm(2,"k" , 3.875, 0.1, 0,0, ierflg); | |
296 | // gMinuit->mnparm(3,"eff" , 1, 1, 0.9 , 1 , ierflg); | |
297 | gMinuit->mnparm(3,"eff" , 0.98947, 0.1, 0,0, ierflg); | |
298 | } | |
e6f3f2fe | 299 | |
d20baea6 | 300 | // Call migrad |
301 | // arglist[0] = 100000; // max calls | |
302 | // arglist[1] = 0.1; // tolerance | |
303 | // arglist[0] = 50; // max calls | |
304 | // arglist[1] = 1; // tolerance | |
305 | // gMinuit->mnexcm("MIGrad",arglist, 2, ierflg); | |
306 | arglist[0] = 100; // max calls | |
307 | arglist[1] = 0.1; // tolerance | |
308 | gMinuit->mnexcm("SIMPLEX",arglist, 2, ierflg); | |
309 | // // fFastFit = 2; | |
310 | arglist[0] = 1000; // max calls | |
311 | arglist[1] = 0.1; // tolerance | |
312 | gMinuit->mnexcm("MIGrad",arglist, 2, ierflg); | |
313 | // gMinuit->mnexcm("IMPROVE",arglist, 0, ierflg); | |
314 | ||
315 | if (ierflg != 0) { | |
316 | AliWarning("Abnormal termination of minimization."); | |
317 | // exit(0); | |
318 | } | |
319 | ||
320 | // if(opt.Contains("M")) { | |
321 | // gMinuit->mnexcm("IMPROVE",arglist, 0, ierflg); | |
322 | // } | |
323 | // if(opt.Contains("E")) { | |
324 | // gMinuit->mnexcm("HESSE",arglist, 0, ierflg); | |
325 | // gMinuit->mnexcm("MINOS",arglist, 0, ierflg); | |
326 | // } | |
327 | ||
328 | // ______________________ Get chi2 and Fit Status __________ | |
329 | Double_t amin,edm,errdef; | |
330 | int nvpar, nparx, icstat; | |
e6f3f2fe | 331 | |
d20baea6 | 332 | |
333 | gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat); | |
334 | Double_t chi2min = amin; | |
87dea6ac | 335 | std::cout << "Fit status " << icstat << std::endl; |
d20baea6 | 336 | |
337 | Double_t alpha_min, mu_min, k_min, eff_min; | |
338 | Double_t alpha_mine, mu_mine, k_mine, eff_mine; | |
339 | gMinuit->GetParameter(0, alpha_min , alpha_mine ); | |
340 | gMinuit->GetParameter(1, mu_min , mu_mine ); | |
341 | gMinuit->GetParameter(2, k_min , k_mine ); | |
342 | gMinuit->GetParameter(3, eff_min , eff_mine ); | |
343 | ||
344 | // Print some infos | |
87dea6ac | 345 | std::cout << "chi2 min is " << chi2min << ", " << alpha_min << ", "<< mu_min<< ", " << k_min << ", " << eff_min << std::endl; |
d20baea6 | 346 | |
347 | thistGlau = GlauberHisto(mu_min,k_min,eff_min,alpha_min,hDATA,kTRUE); | |
348 | hGLAU = (TH1D *) thistGlau->Clone("hGLAU"); | |
349 | hGLAU->SetName( ((TString)hDATA->GetName()).Append(Form("_GLAU"))); | |
350 | hGLAU->SetTitle( ((TString)hDATA->GetName()).Append(Form("_GLAU_%.3f_%.3f_%.3f_%.3f",mu_min,k_min,alpha_min,eff_min))); | |
351 | ||
352 | ||
87dea6ac | 353 | std::cout << "fitted " << hGLAU->Integral(hGLAU->FindBin(fMultmin), |
354 | hGLAU->FindBin(fMultmax))/hGLAU->Integral() | |
355 | << " of the total cross section" << std::endl; | |
d20baea6 | 356 | // eff_min=1; |
357 | ||
358 | // Save histogram and ntuple | |
359 | // fFastFit=kFALSE; // Use a better aprox to get the final histo | |
360 | ||
361 | // float mcintegral = hGLAU->Integral(hDATA->FindBin(multmin),hGLAU->GetNbinsX()); | |
362 | // float scale = (hDATA->Integral(hDATA->FindBin(multmin),hDATA->GetNbinsX())/mcintegral) * ((float) eff_min); | |
363 | // hGLAU->Scale(scale); | |
364 | ||
87dea6ac | 365 | std::cout << "Chi2 final" << CalculateChi2(hDATA,hGLAU,eff_min) << std::endl; |
d20baea6 | 366 | |
87dea6ac | 367 | fhEffi = GetTriggerEfficiencyFunction(hDATA, hGLAU); |
368 | SaveHisto(hDATA,hGLAU,fhEffi,outrootfile); | |
d20baea6 | 369 | fGlauberHist=hGLAU; |
370 | } | |
d20baea6 | 371 | // close inrootfile, outrootfile |
87dea6ac | 372 | finrootfile->Close(); |
d20baea6 | 373 | outrootfile->Close(); |
d20baea6 | 374 | } |
375 | ||
d20baea6 | 376 | //-------------------------------------------------------------------------------------------------- |
87dea6ac | 377 | TH1D *AliCentralityGlauberFit::GlauberHisto(Float_t mu, Float_t k, Float_t eff, Float_t alpha, |
378 | TH1D *hDATA, Bool_t save) | |
379 | { | |
380 | // Get Glauber histogram | |
d20baea6 | 381 | |
87dea6ac | 382 | eff=eff; |
d20baea6 | 383 | static TStopwatch sw; |
384 | ||
385 | TH1D *hSample = fFastFit != 2 ? NBDhist(mu,k) : 0; | |
386 | // fNBD->SetParameters(mu,k); | |
387 | static TH1D *h1 = (TH1D*)hDATA->Clone(); | |
388 | h1->Reset(); | |
389 | // h1->SetName(Form("fit_%.3f_%.3f_%.3f_%.3f",mu,k,eff,alpha)); | |
e6f3f2fe | 390 | TFile *outFile = NULL; |
391 | TNtuple *ntuple = NULL; | |
392 | ||
393 | if (save) { | |
394 | outFile = TFile::Open("pippo.root", "RECREATE"); | |
d20baea6 | 395 | ntuple = new TNtuple("gnt", "Glauber ntuple", "Npart:Ncoll:B:tAA:ntot"); |
e6f3f2fe | 396 | } |
397 | ||
d20baea6 | 398 | |
399 | Int_t nents = 100000;//glau_ntuple->GetEntries(); | |
400 | //if(fFastFit > 0 && fFastFit < 2) nents = 10000; | |
401 | // Int_t nents = 100;//glau_ntuple->GetEntries(); | |
e6f3f2fe | 402 | for (Int_t i=0;i<nents;++i) { |
87dea6ac | 403 | fGlauntuple->GetEntry(i); |
d20baea6 | 404 | //Int_t n = TMath::Nint(TMath::Power(Npart,alpha)); |
405 | //Int_t n = TMath::Nint(TMath::Power(Ncoll,alpha)); | |
87dea6ac | 406 | Int_t n = alpha * fNpart + (1-alpha) * fNcoll; |
e6f3f2fe | 407 | Int_t ntot=0; |
d20baea6 | 408 | // ntot=n*fNBD->GetRandom(); |
409 | if (fFastFit == 1) { | |
410 | // NBD | |
411 | ntot = n*hSample->GetRandom(); | |
412 | } | |
413 | else if (fFastFit == 2) { | |
414 | // Gaussian | |
415 | Double_t sigma = k*TMath::Sqrt(n*mu); | |
416 | ntot = gRandom->Gaus(n*mu,sigma); | |
417 | } | |
418 | else { | |
419 | for(Int_t j = 0; j<n; ++j) | |
420 | ntot+=hSample->GetRandom(); | |
421 | // ntot+=fNBD->GetRandom(); | |
422 | } | |
e6f3f2fe | 423 | h1->Fill(ntot); |
424 | ||
425 | if (save) | |
87dea6ac | 426 | ntuple->Fill(fNpart,fNcoll,fB,fTaa,ntot); |
e6f3f2fe | 427 | |
428 | } | |
d20baea6 | 429 | |
e6f3f2fe | 430 | if (save) { |
431 | ntuple->Write(); | |
432 | outFile->Close(); | |
433 | } | |
434 | ||
d20baea6 | 435 | if (fFastFit != 2) delete hSample; |
e6f3f2fe | 436 | return h1; |
437 | ||
438 | } | |
439 | ||
440 | //-------------------------------------------------------------------------------------------------- | |
441 | ||
442 | Float_t AliCentralityGlauberFit::CalculateChi2(TH1D *hDATA, TH1D *thistGlau, Float_t eff) { | |
443 | // note that for different values of neff the mc distribution is identical, just re-scaled below | |
444 | // normalize the two histogram | |
445 | // scale MC up but leave data alone - that preserves error bars !!!! | |
446 | ||
87dea6ac | 447 | int lowchibin = hDATA->FindBin(fMultmin); |
448 | int highchibin = hDATA->FindBin(fMultmax); | |
d20baea6 | 449 | |
450 | float mcintegral = thistGlau->Integral(lowchibin,mcintegral); | |
451 | // float scale = (hDATA->Integral(lowchibin,highchibin)/mcintegral) * ((float) eff); | |
452 | float scale = (hDATA->Integral(lowchibin,highchibin)/mcintegral) * ((float) eff); | |
453 | thistGlau->Scale(scale); | |
454 | ||
455 | // FIXME: Kolmogorov | |
456 | //return hDATA->KolmogorovTest(thistGlau,"M"); | |
457 | ||
458 | // // calculate the chi2 between MC and real data over some range ???? | |
459 | if (fUseChi2) { | |
460 | float chi2 = 0.0; | |
461 | for (int i=lowchibin; i<=highchibin; i++) { | |
462 | if (hDATA->GetBinContent(i) < 1.0) continue; | |
463 | float diff = pow( (thistGlau->GetBinContent(i) - hDATA->GetBinContent(i)) , 2); | |
464 | diff = diff / (pow(hDATA->GetBinError(i) , 2)+pow(thistGlau->GetBinError(i) , 2)); // FIXME squared distance if commented | |
465 | chi2 += diff; | |
466 | } | |
467 | chi2 = chi2 / (highchibin - lowchibin + 1); | |
468 | return chi2; | |
e6f3f2fe | 469 | } |
d20baea6 | 470 | // "-2 log likelihood ratio(mu;n) = 2[(mu - n) + n ln(n/mu)]" |
471 | else { | |
87dea6ac | 472 | std::cout << "LL" << std::endl; |
d20baea6 | 473 | |
474 | ||
475 | float ll = 0.0; | |
476 | for (int i=lowchibin; i<=highchibin; i++) { | |
477 | // if (thistGlau->GetBinContent(i) < 1) continue; | |
478 | // if (hDATA->GetBinContent(i) < 1) continue; | |
479 | // cout << ll << " " << thistGlau->GetBinContent(i) <<" " << hDATA->GetBinContent(i) << endl; | |
480 | ||
481 | Float_t data = hDATA ->GetBinContent(i); | |
482 | Float_t mc = thistGlau->GetBinContent(i); | |
483 | Int_t idata = TMath::Nint(data); | |
484 | if (mc < 1.e-9) mc = 1.e-9; | |
485 | Float_t fsub = - mc + idata * TMath::Log(mc); | |
486 | Float_t fobs = 0; | |
87dea6ac | 487 | //Int_t imc = TMath::Nint(mc); |
d20baea6 | 488 | if (idata > 0) { |
489 | for(Int_t istep = 0; istep < idata; istep++){ | |
490 | if (istep > 1) | |
491 | fobs += TMath::Log(istep); | |
492 | } | |
493 | } | |
494 | // cout << mc << " " << data << " " << fsub << " " << fobs << endl; | |
495 | fsub -= fobs; | |
496 | ll -= fsub ; | |
497 | } | |
498 | return 2*ll; | |
e6f3f2fe | 499 | } |
d20baea6 | 500 | |
e6f3f2fe | 501 | } |
502 | ||
503 | //-------------------------------------------------------------------------------------------------- | |
504 | ||
87dea6ac | 505 | TH1D * AliCentralityGlauberFit::GetTriggerEfficiencyFunction(TH1D *hist1, TH1D *hist2) |
506 | { | |
507 | // Get efficiency. | |
508 | ||
509 | fhEffi = (TH1D*)hist1->Clone("heffi"); | |
510 | fhEffi->Divide(hist2); | |
511 | return fhEffi; | |
e6f3f2fe | 512 | } |
513 | ||
514 | //-------------------------------------------------------------------------------------------------- | |
515 | ||
87dea6ac | 516 | Float_t AliCentralityGlauberFit::GetTriggerEfficiencyIntegral(TH1D *hist1, TH1D *hist2) |
517 | { | |
518 | // Get eff integral. | |
519 | fEffi = hist1->Integral()/hist2->Integral(); | |
520 | return fEffi; | |
e6f3f2fe | 521 | } |
522 | ||
523 | //-------------------------------------------------------------------------------------------------- | |
524 | ||
525 | void AliCentralityGlauberFit::SaveHisto(TH1D *hist1, TH1D *hist2, TH1D *heffi, TFile *outrootfile) { | |
526 | outrootfile->cd(); | |
527 | hist1->Write(); | |
528 | hist2->Write(); | |
529 | heffi->Write(); | |
530 | } | |
531 | ||
532 | //-------------------------------------------------------------------------------------------------- | |
533 | ||
534 | Double_t AliCentralityGlauberFit::NBD(Int_t n, Double_t mu, Double_t k) | |
535 | { | |
87dea6ac | 536 | // Compute NBD. |
d20baea6 | 537 | 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 | 538 | return ret; |
539 | } | |
540 | ||
541 | //-------------------------------------------------------------------------------------------------- | |
542 | ||
543 | TH1D *AliCentralityGlauberFit::NBDhist(Double_t mu, Double_t k) | |
544 | { | |
d20baea6 | 545 | TH1D *h = new TH1D("htemp","",100,-0.5,299.5); |
e6f3f2fe | 546 | h->SetName(Form("nbd_%f_%f",mu,k)); |
547 | h->SetDirectory(0); | |
d20baea6 | 548 | for (Int_t i=0;i<300;++i) { |
e6f3f2fe | 549 | Double_t val = NBD(i,mu,k); |
d20baea6 | 550 | if (val>1e-20) h->Fill(i,val); |
e6f3f2fe | 551 | } |
552 | return h; | |
553 | } | |
554 | ||
87dea6ac | 555 | void AliCentralityGlauberFit::MinuitFcnNBD(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) |
556 | { | |
d20baea6 | 557 | // FCN for minuit |
558 | Double_t alpha = par[0]; | |
559 | Double_t mu = par[1]; | |
560 | Double_t k = par[2]; | |
561 | Double_t eff = par[3]; | |
562 | //Double_t eff = 1;//par[3]; | |
563 | ||
564 | AliCentralityGlauberFit * obj = (AliCentralityGlauberFit *) gMinuit->GetObjectFit(); | |
565 | ||
566 | // static TStopwatch sw; | |
567 | // sw.Start(); | |
568 | TH1D * thistGlau = obj->GlauberHisto(mu,k,eff,alpha,obj->GetTempHist(),kFALSE); | |
569 | f = obj->CalculateChi2(obj->GetTempHist(),thistGlau,eff); | |
570 | // sw.Stop(); | |
571 | // sw.Print(); | |
87dea6ac | 572 | Printf("%f - %f - %f - %f - %f \n",f,alpha,mu,k,eff); |
d20baea6 | 573 | } |
574 | ||
87dea6ac | 575 | Double_t AliCentralityGlauberFit::NBDFunc(Double_t * x, Double_t *par) |
576 | { | |
d20baea6 | 577 | // TF1 interface |
d20baea6 | 578 | Double_t mu = par[0]; |
579 | Double_t k = par[1]; | |
580 | Double_t n = x[0]; | |
d20baea6 | 581 | 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 | 582 | return ret; |
d20baea6 | 583 | } |