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