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