| 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 | } |
git://git.uio.no / u / mrichter / AliRoot.git / blame