#include <TStyle.h>
#include <TGraphErrors.h>
#include <vector>
+#include <TMinuit.h>
+#include <TStopwatch.h>
+#include <TRandom.h>
#include "AliCentralityGlauberFit.h"
-#include "TMinuit.h"
#include "AliLog.h"
-#include "TStopwatch.h"
-#include "TRandom.h"
+
ClassImp(AliCentralityGlauberFit)
//______________________________________________________________________________
+AliCentralityGlauberFit::AliCentralityGlauberFit(const char *filename) :
+ fNmu(0),
+ fMulow(0),
+ fMuhigh(0),
+ fNk(0),
+ fKlow(0),
+ fKhigh(0),
+ fNalpha(0),
+ fAlphalow(0),
+ fAlphahigh(0),
+ fRebinFactor(0),
+ fScalemin(0),
+ fMultmin(0),
+ fMultmax(0),
+ fGlauntuple(0),
+ fNpart(0),
+ fNcoll(0),
+ fB(0),
+ fTaa(0),
+ fEffi(0),
+ fhEffi(0),
+ fTempHist(0),
+ fGlauberHist(0),
+ fFastFit(0),
+ fAncestor(2),
+ fNBD(0),
+ fUseChi2(kTRUE),
+ fUseAverage(kFALSE),
+ fhAncestor(0),
+ fNevents(100000),
+ fNtrials(100),
+ fInrootfilename(0),
+ fInntuplename(0),
+ fOutrootfilename(0),
+ fOutntuplename(0),
+ fAncfilename("ancestor_hists.root"),
+ fHistnames()
+{
+ // Standard constructor.
+ TFile *f = 0;
+ if (filename) { // glauber file
+ f = TFile::Open(filename);
+ fGlauntuple = (TNtuple*) f->Get("nt_Pb_Pb");
+ fGlauntuple->SetBranchAddress("Npart",&fNpart);
+ fGlauntuple->SetBranchAddress("Ncoll",&fNcoll);
+ fGlauntuple->SetBranchAddress("B",&fB);
+ fGlauntuple->SetBranchAddress("tAA",&fTaa);
+ }
-
-AliCentralityGlauberFit::AliCentralityGlauberFit(const char * filename) : fFastFit(0), fTempHist(0), fGlauberHist(0), fUseChi2(kTRUE) {
- // standard constructor
- // glauber file
- TFile *f = TFile::Open(filename);
- glau_ntuple = (TNtuple*) f->Get("nt_Pb_Pb");
- glau_ntuple->SetBranchAddress("Npart",&Npart);
- glau_ntuple->SetBranchAddress("Ncoll",&Ncoll);
- glau_ntuple->SetBranchAddress("B",&B);
- glau_ntuple->SetBranchAddress("tAA",&tAA);
-
fNBD = new TF1 ("fNBD", AliCentralityGlauberFit::NBDFunc, 0, 100,2);
-
-}
-
-//--------------------------------------------------------------------------------------------------
-
-AliCentralityGlauberFit::~AliCentralityGlauberFit() {
- // destructor
-}
-
-//--------------------------------------------------------------------------------------------------
-
-void AliCentralityGlauberFit::AddHisto(TString name) {
- histnames.push_back(name);
}
//--------------------------------------------------------------------------------------------------
+void AliCentralityGlauberFit::SetRangeToFit(Double_t fmultmin, Double_t fmultmax)
+{
+ // Set fit range.
-void AliCentralityGlauberFit::SetOutputFile(TString filename) {
- outrootfilename = filename;
+ fMultmin=fmultmin;
+ fMultmax=fmultmax;
}
//--------------------------------------------------------------------------------------------------
-
-void AliCentralityGlauberFit::SetRangeToFit(Float_t fmultmin, Float_t fmultmax)
+void AliCentralityGlauberFit::SetRangeToScale(Double_t fscalemin)
{
- multmin=fmultmin;
- multmax=fmultmax;
+ // Set range where to scale simulated histo to data
+
+ fScalemin=fscalemin;
}
//--------------------------------------------------------------------------------------------------
-
void AliCentralityGlauberFit::SetGlauberParam(
- Int_t fNmu,
- Float_t fmulow,
- Float_t fmuhigh,
- Int_t fNk,
- Float_t fklow,
- Float_t fkhigh,
- Int_t fNalpha,
- Float_t falphalow,
- Float_t falphahigh,
- Int_t fNeff,
- Float_t fefflow,
- Float_t feffhigh)
-
+ Int_t Nmu,
+ Double_t mulow,
+ Double_t muhigh,
+ Int_t Nk,
+ Double_t klow,
+ Double_t khigh,
+ Int_t Nalpha,
+ Double_t alphalow,
+ Double_t alphahigh)
{
- Nmu=fNmu;
- mulow=fmulow;
- muhigh=fmuhigh;
- Nk=fNk;
- klow=fklow;
- khigh=fkhigh;
- Nalpha=fNalpha;
- alphalow=falphalow;
- alphahigh=falphahigh;
- Neff=fNeff;
- efflow=fefflow;
- effhigh=feffhigh;
+ // Set Glauber parameters.
+
+ fNmu=Nmu;
+ fMulow=mulow;
+ fMuhigh=muhigh;
+ fNk=Nk;
+ fKlow=klow;
+ fKhigh=khigh;
+ fNalpha=Nalpha;
+ fAlphalow=alphalow;
+ fAlphahigh=alphahigh;
}
//--------------------------------------------------------------------------------------------------
+void AliCentralityGlauberFit::MakeFits()
+{
+ // Make fits.
-Float_t AliCentralityGlauberFit::GetPercentileCrossSection() {
- return percentXsec;
-}
-
-//--------------------------------------------------------------------------------------------------
-void AliCentralityGlauberFit::SetRebin(Int_t frebinFactor) {
- rebinFactor=frebinFactor;
-}
-//--------------------------------------------------------------------------------------------------
-
-void AliCentralityGlauberFit::MakeFits(TString infilename) {
- TH1D *hDATA;
- TH1D *thistGlau;
- FILE* fTxt = fopen ("parameters.txt","w");
- TFile *outrootfile;
+ TH1F *hDATA;
+ TH1F *thistGlau;
+ TFile *inrootfile;
+ TFile *outrootfile;
+ //FILE* fTxt = fopen ("parameters.txt","w");
// open inrootfile, outrootfile
- inrootfile = new TFile(infilename);
- outrootfile = new TFile(outrootfilename,"RECREATE");
+ std::cout << "input file " << fInrootfilename << std::endl;
+ std::cout << "output file " << fOutrootfilename << std::endl;
+ inrootfile = new TFile(fInrootfilename,"OPEN");
+ outrootfile = new TFile(fOutrootfilename,"RECREATE");
// loop over all distribution names
- vector<TString>::const_iterator hni;
- for(hni=histnames.begin(); hni!=histnames.end(); hni++) {
- hDATA = (TH1D*) (inrootfile->Get(*hni));
- //TList *list = (TList*) (inrootfile->Get("coutput1"));
- //hDATA = (TH1D*) (list->FindObject(*hni));
- hDATA->Rebin(rebinFactor);
- TH1D *hGLAU = new TH1D("hGLAU","hGLAU",hDATA->GetNbinsX(),0,hDATA->GetNbinsX()*hDATA->GetBinWidth(1));
- Float_t chi2min = 9999999.0;
- Float_t alpha_min=-1;
- Float_t mu_min=-1;
- Float_t k_min=-1;
- Float_t eff_min=-1;
- Float_t alpha, mu, k, eff, chi2;
-
- for (int nalpha=0;nalpha<Nalpha;nalpha++) {
- alpha = alphalow + ((float) nalpha ) * (alphahigh - alphalow ) / Nalpha;
+ std::vector<TString>::const_iterator hni;
+ for(hni=fHistnames.begin(); hni!=fHistnames.end(); hni++) {
+ hDATA = (TH1F*) (inrootfile->Get(*hni));
+ if (!hDATA) {
+ TList *list = (TList*) (inrootfile->Get("coutput1"));
+ hDATA = (TH1F*) (list->FindObject(*hni));
+ }
+ hDATA->Rebin(fRebinFactor);
+ //TH1F *hGLAU = new TH1F("hGLAU","hGLAU",hDATA->GetNbinsX(),0,hDATA->GetNbinsX()*hDATA->GetBinWidth(1));
+ Double_t chi2min = 9999999.0;
+ Double_t alpha_min=-1;
+ Double_t mu_min=-1;
+ Double_t k_min=-1;
+ Double_t alpha, mu, k, chi2;
+
+ for (Int_t nalpha=0;nalpha<fNalpha;nalpha++) {
+ alpha = fAlphalow + ((Double_t) nalpha ) * (fAlphahigh - fAlphalow ) / fNalpha;
+
mu=0.0;
- for (int nmu=0; nmu<Nmu; nmu++) {
- mu = (mulow*(1-nalpha*0.05) ) + ((float) nmu ) * (muhigh - mulow ) / Nmu;
- //mu = mulow + ((float) nmu ) * (muhigh - mulow ) / Nmu;
-
- for (int nk=0; nk<Nk; nk++) {
- k = klow + ((float) nk ) * (khigh - klow ) / Nk;
-
- for (int neff=0; neff<Neff; neff++) {
- eff = efflow + ((float) neff) * (effhigh - efflow) / Neff;
-
- thistGlau = GlauberHisto(mu,k,eff,alpha,hDATA,kFALSE);
- chi2 = CalculateChi2(hDATA,thistGlau,eff);
- fprintf(fTxt, "%3.3f %3.3f %3.3f %3.3f %3.3f\n",(float) eff, (float) mu, (float) k, (float) alpha, chi2);
-
- if ( chi2 < chi2min ) {
- chi2min = chi2;
- alpha_min=alpha;
- mu_min=mu;
- k_min=k;
- eff_min=eff;
- }
+ for (Int_t nmu=0; nmu<fNmu; nmu++) {
+ mu = fMulow + ((Double_t) nmu ) * (fMuhigh - fMulow ) / fNmu;
+
+ for (Int_t nk=0; nk<fNk; nk++) {
+ k = fKlow + ((Double_t) nk ) * (fKhigh - fKlow ) / fNk;
+
+ thistGlau = GlauberHisto(mu,k,alpha,hDATA,kFALSE);
+ chi2 = CalculateChi2(hDATA,thistGlau);
+ //fprintf(fTxt, "%3.3f %3.3f %3.3f %3.3f\n", (float) mu, (float) k, (float) alpha, chi2);
+ if ( chi2 < chi2min ) {
+ chi2min = chi2;
+ alpha_min=alpha;
+ mu_min=mu;
+ k_min=k;
}
}
}
}
-
- thistGlau = GlauberHisto(mu_min,k_min,eff_min,alpha_min,hDATA,kTRUE);
- hGLAU = (TH1D *) thistGlau->Clone("hGLAU");
- //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))));
+
+ thistGlau = GlauberHisto(mu_min,k_min,alpha_min,hDATA,kTRUE);
+
+ TH1F * hGLAU = 0x0;
+ hGLAU = (TH1F *) thistGlau->Clone("hGLAU");
hGLAU->SetName( ((TString)hDATA->GetName()).Append(Form("_GLAU")));
- hGLAU->SetTitle( ((TString)hDATA->GetName()).Append(Form("_GLAU_%.3f_%.3f_%.3f_%.3f",mu_min,k_min,alpha_min,eff_min)));
+ hGLAU->SetTitle( ((TString)hDATA->GetName()).Append(Form("_GLAU_%.3f_%.3f_%.3f",
+ mu_min,k_min,alpha_min)));
- float mcintegral = hGLAU->Integral(hDATA->FindBin(multmin),hGLAU->GetNbinsX());
- float scale = (hDATA->Integral(hDATA->FindBin(multmin),hDATA->GetNbinsX())/mcintegral) * ((float) eff_min);
+ Double_t mcintegral = hGLAU->Integral(hGLAU->FindBin(fScalemin),hGLAU->GetNbinsX());
+ Double_t scale = (hDATA->Integral(hDATA->FindBin(fScalemin),hDATA->GetNbinsX())/mcintegral);
hGLAU->Scale(scale);
- heffi = GetTriggerEfficiencyFunction(hDATA, hGLAU);
- SaveHisto(hDATA,hGLAU,heffi,outrootfile);
- fclose (fTxt);
-
- cout << "chi2 min is " << chi2min << endl;
- cout << "fitted " << hGLAU->Integral(hGLAU->FindBin(multmin),hGLAU->FindBin(multmax))/hGLAU->Integral() << " of the total cross section" << endl;
+ fhEffi = GetTriggerEfficiencyFunction(hDATA, hGLAU);
+ SaveHisto(hDATA,hGLAU,fhEffi,outrootfile);
+ //fclose (fTxt);
+ std::cout << "chi2 min is " << chi2min << std::endl;
+ std::cout << "fitted " << hGLAU->Integral(hGLAU->FindBin(fMultmin),
+ hGLAU->FindBin(fMultmax))/hGLAU->Integral()
+ << " of the total cross section" << std::endl;
fTempHist=hDATA;
fGlauberHist=hGLAU;
}
// close inrootfile, outrootfile
inrootfile->Close();
outrootfile->Close();
-
}
//--------------------------------------------------------------------------------------------------
-void AliCentralityGlauberFit::MakeFitsMinuitNBD(TString infilename) {
- TH1D *hDATA;
- TH1D *thistGlau;
+void AliCentralityGlauberFit::MakeFitsMinuitNBD(Double_t alpha, Double_t mu, Double_t k)
+{
+ // Make fits using Minuit.
+
+ if (alpha<0)
+ alpha = fAlphalow;
+ if (mu<0)
+ mu = fMulow;
+ if (k<0)
+ k = fKlow;
+ printf("Calling Minuit with starting values: %f %f %f\n", alpha, mu, k);
+
+ TH1F *hDATA;
+ TH1F *thistGlau;
+ TFile *inrootfile;
TFile *outrootfile;
// open inrootfile, outrootfile
- inrootfile = new TFile(infilename);
- outrootfile = new TFile(outrootfilename,"RECREATE");
-
+ std::cout << "input file " << fInrootfilename << std::endl;
+ std::cout << "output file " << fOutrootfilename << std::endl;
+ inrootfile = new TFile(fInrootfilename,"OPEN");
+ outrootfile = new TFile(fOutrootfilename,"RECREATE");
+
// loop over all distribution names
- vector<TString>::const_iterator hni;
- for(hni=histnames.begin(); hni!=histnames.end(); hni++) {
- fFastFit=0;
- fUseChi2 = 1;
- hDATA = (TH1D*) (inrootfile->Get(*hni));
- hDATA->Rebin(rebinFactor);
+ std::vector<TString>::const_iterator hni;
+ for(hni=fHistnames.begin(); hni!=fHistnames.end(); hni++) {
+ hDATA = (TH1F*) (inrootfile->Get(*hni));
+ hDATA->Rebin(fRebinFactor);
fTempHist=hDATA;
- TH1D *hGLAU = new TH1D("hGLAU","hGLAU",hDATA->GetNbinsX(),0,hDATA->GetNbinsX()*hDATA->GetBinWidth(1));
+ TH1F *hGLAU = new TH1F("hGLAU","hGLAU",hDATA->GetNbinsX(),0,hDATA->GetNbinsX()*hDATA->GetBinWidth(1));
hGLAU->Sumw2();
+
// Minimize here
if(gMinuit) delete gMinuit;
new TMinuit(3);
Double_t arglist[2]={0};
Int_t ierflg;
if (fUseChi2) arglist[0] = 1; // should be 1 if you want to minimize chi2
- else arglist[0] = 0.5; // should be 0.5 for ll
+ else arglist[0] = 0.5; // should be 0.5 for ll
gMinuit->mnexcm("SET ERR",arglist, 1, ierflg);
-
- // // Change verbosity
- // gMinuit->Command((TString("SET PRINTOUT ")+long(fVerbosity)).Data());
-
- // Set parameters
- // start from the middle of the allowed range, as a guess
- // gMinuit->mnparm(0,"alpha", (alphalow+alphahigh)/2, 1, alphalow, alphahigh, ierflg);
- // gMinuit->mnparm(1,"mu" , (mulow +muhigh )/2, 1, mulow , muhigh , ierflg);
- // gMinuit->mnparm(2,"k" , (klow +khigh )/2, 1, klow , khigh , ierflg);
- // gMinuit->mnparm(3,"eff" , (efflow +effhigh )/2, 1, efflow , effhigh , ierflg);
-
- // NBD with limits
- // gMinuit->mnparm(0,"alpha", 1.104, 0.1, alphalow, alphahigh, ierflg);
- // gMinuit->mnparm(1,"mu" , 65, 1, mulow , muhigh , ierflg);
- // gMinuit->mnparm(2,"k" , 0.41, 0.1, klow , khigh , ierflg);
- // // gMinuit->mnparm(3,"eff" , 1, 1, 0.9 , 1 , ierflg);
- // gMinuit->mnparm(3,"eff" , 0.95, 0.1, efflow , effhigh , ierflg);
-
- // // NBD, no limits
- // gMinuit->mnparm(0,"alpha", 1.104, 0.1, 0,0, ierflg);
- // gMinuit->mnparm(1,"mu" , 65, 1, 0,0, ierflg);
- // gMinuit->mnparm(2,"k" , 0.41, 0.1, 0,0, ierflg);
- // gMinuit->mnparm(3,"eff" , 0.95, 0.1, 0,0 , ierflg);
-
- if (!fFastFit) {
- // // NBD, Npart**alpha
- //gMinuit->mnparm(0,"alpha", 1.16, 0.1, 0,0, ierflg);
- //gMinuit->mnparm(1,"mu" , 27.55, 1, 0,0, ierflg);
- //gMinuit->mnparm(2,"k" , 0.64, 0.1, 0,0, ierflg);
- //gMinuit->mnparm(3,"eff" , 0.97, 0.1, 0,0 , ierflg);
- // V0 non rescaled
- // // NBD, alpha Npart * (1-alpha) Ncoll
- gMinuit->mnparm(0,"alpha", 0.832, 0.1, 0,0, ierflg);
- gMinuit->mnparm(1,"mu" , 29, 1, 0,0, ierflg);
- gMinuit->mnparm(2,"k" , 1.4, 0.1, 0,0, ierflg);
- gMinuit->mnparm(3,"eff" , 0.991, 0.1, 0,0 , ierflg);
- // SPD
- //gMinuit->mnparm(0,"alpha", 0.820, 0.1, 0.7,1, ierflg);
- //gMinuit->mnparm(1,"mu" , 8.267, 1, 1,100, ierflg);
- //gMinuit->mnparm(2,"k" , 0.530, 0.1, 0,0, ierflg);
- //gMinuit->mnparm(3,"eff" , 0.990, 0.1, 0,0 , ierflg);
- } else if (fFastFit == 2) {
- //gaussian
- //0.859136, 30.2057, 3.87565, 0.989747
- gMinuit->mnparm(0,"alpha", 0.59, 0.1, 0,1, ierflg);
- gMinuit->mnparm(1,"mu" , 30.2, 1, 0,0, ierflg);
- gMinuit->mnparm(2,"k" , 3.875, 0.1, 0,0, ierflg);
- // gMinuit->mnparm(3,"eff" , 1, 1, 0.9 , 1 , ierflg);
- gMinuit->mnparm(3,"eff" , 0.98947, 0.1, 0,0, ierflg);
- }
+ gMinuit->mnparm(0,"alpha", alpha, (fAlphahigh-fAlphalow)/fNalpha, fAlphalow, fAlphahigh, ierflg);
+ gMinuit->mnparm(1,"mu" , mu, (fMuhigh-fMulow)/fNmu, fMulow, fMuhigh, ierflg);
+ gMinuit->mnparm(2,"k" , k, (fKhigh-fKlow)/fNk, fKlow, fKhigh, ierflg);
// Call migrad
- // arglist[0] = 100000; // max calls
- // arglist[1] = 0.1; // tolerance
- // arglist[0] = 50; // max calls
- // arglist[1] = 1; // tolerance
- // gMinuit->mnexcm("MIGrad",arglist, 2, ierflg);
arglist[0] = 100; // max calls
arglist[1] = 0.1; // tolerance
gMinuit->mnexcm("SIMPLEX",arglist, 2, ierflg);
- // // fFastFit = 2;
arglist[0] = 1000; // max calls
- arglist[1] = 0.1; // tolerance
+ arglist[1] = 0.1; // tolerance
gMinuit->mnexcm("MIGrad",arglist, 2, ierflg);
- // gMinuit->mnexcm("IMPROVE",arglist, 0, ierflg);
+ //gMinuit->mnexcm("IMPROVE",arglist, 0, ierflg);
if (ierflg != 0) {
AliWarning("Abnormal termination of minimization.");
- // exit(0);
}
- // if(opt.Contains("M")) {
- // gMinuit->mnexcm("IMPROVE",arglist, 0, ierflg);
- // }
- // if(opt.Contains("E")) {
- // gMinuit->mnexcm("HESSE",arglist, 0, ierflg);
- // gMinuit->mnexcm("MINOS",arglist, 0, ierflg);
- // }
-
// ______________________ Get chi2 and Fit Status __________
Double_t amin,edm,errdef;
- int nvpar, nparx, icstat;
-
+ Int_t nvpar, nparx, icstat;
gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
Double_t chi2min = amin;
- cout << "Fit status " << icstat << endl;
+ std::cout << "Fit status " << icstat << std::endl;
- Double_t alpha_min, mu_min, k_min, eff_min;
- Double_t alpha_mine, mu_mine, k_mine, eff_mine;
+ Double_t alpha_min, mu_min, k_min;
+ Double_t alpha_mine, mu_mine, k_mine;
gMinuit->GetParameter(0, alpha_min , alpha_mine );
gMinuit->GetParameter(1, mu_min , mu_mine );
gMinuit->GetParameter(2, k_min , k_mine );
- gMinuit->GetParameter(3, eff_min , eff_mine );
- // Print some infos
- cout << "chi2 min is " << chi2min << ", " << alpha_min << ", "<< mu_min<< ", " << k_min << ", " << eff_min << endl;
+ // print some infos
+ std::cout << "chi2 min is " << chi2min << ", " << alpha_min << ", "<< mu_min<< ", "
+ << k_min << std::endl;
- thistGlau = GlauberHisto(mu_min,k_min,eff_min,alpha_min,hDATA,kTRUE);
- hGLAU = (TH1D *) thistGlau->Clone("hGLAU");
+ thistGlau = GlauberHisto(mu_min,k_min,alpha_min,hDATA,kTRUE);
+
+ hGLAU = (TH1F *) thistGlau->Clone("hGLAU");
hGLAU->SetName( ((TString)hDATA->GetName()).Append(Form("_GLAU")));
- hGLAU->SetTitle( ((TString)hDATA->GetName()).Append(Form("_GLAU_%.3f_%.3f_%.3f_%.3f",mu_min,k_min,alpha_min,eff_min)));
+ hGLAU->SetTitle( ((TString)hDATA->GetName()).Append(Form("_GLAU_%.3f_%.3f_%.3f",
+ mu_min,k_min,alpha_min)));
- cout << "fitted " << hGLAU->Integral(hGLAU->FindBin(multmin),hGLAU->FindBin(multmax))/hGLAU->Integral() << " of the total cross section" << endl;
- // eff_min=1;
-
- // Save histogram and ntuple
- // fFastFit=kFALSE; // Use a better aprox to get the final histo
+ std::cout << "fitted " << hGLAU->Integral(hGLAU->FindBin(fMultmin),
+ hGLAU->FindBin(fMultmax))/hGLAU->Integral()
+ << " of the total cross section" << std::endl;
- // float mcintegral = hGLAU->Integral(hDATA->FindBin(multmin),hGLAU->GetNbinsX());
- // float scale = (hDATA->Integral(hDATA->FindBin(multmin),hDATA->GetNbinsX())/mcintegral) * ((float) eff_min);
- // hGLAU->Scale(scale);
-
- cout << "Chi2 final" << CalculateChi2(hDATA,hGLAU,eff_min) << endl;
-
+ Double_t mcintegral = hGLAU->Integral(hGLAU->FindBin(fScalemin),hGLAU->GetNbinsX());
+ Double_t scale = (hDATA->Integral(hDATA->FindBin(fScalemin),hDATA->GetNbinsX())/mcintegral);
+ hGLAU->Scale(scale);
- heffi = GetTriggerEfficiencyFunction(hDATA, hGLAU);
- SaveHisto(hDATA,hGLAU,heffi,outrootfile);
+ std::cout << "Chi2 final " << CalculateChi2(hDATA,hGLAU) << std::endl;
+ fhEffi = GetTriggerEfficiencyFunction(hDATA, hGLAU);
+ SaveHisto(hDATA,hGLAU,fhEffi,outrootfile);
fGlauberHist=hGLAU;
}
-
-
// close inrootfile, outrootfile
inrootfile->Close();
outrootfile->Close();
-
}
-
//--------------------------------------------------------------------------------------------------
+TH1F *AliCentralityGlauberFit::GlauberHisto(Double_t mu, Double_t k, Double_t alpha,
+ TH1F *hDATA, Bool_t save)
+{
+ // Get Glauber histogram.
+ static TH1F *h1 = (TH1F*)hDATA->Clone();
+ h1->Reset();
+ h1->SetName(Form("fit_%.3f_%.3f_%.3f",mu,k,alpha));
+
+ if (fUseAverage) {
+ fhAncestor = MakeAncestor(alpha);
+ for (Int_t np=1; np<=fhAncestor->GetNbinsX(); ++np) {
+ Double_t nanc = fhAncestor->GetBinCenter(np);
+ Double_t weights = fhAncestor->GetBinContent(np);
+
+ if (weights <= 0) continue;
+ Int_t trials = (Int_t) (20 * nanc * (int) mu);
+ if (trials <=0) continue;
+ for (Int_t j=0; j<trials; j++) {
+ double nbdvalue = NBD(j, mu * nanc, k * nanc);
+ h1->Fill((double) j, nbdvalue * weights);
+ }
+ }
+ return h1;
+ }
-TH1D * AliCentralityGlauberFit::GlauberHisto(Float_t mu, Float_t k, Float_t eff, Float_t alpha, TH1D *hDATA, Bool_t save) {
-
- static TStopwatch sw;
+ TH1F *hSample = fFastFit != 2 ? NBDhist(mu,k) : 0;
- TH1D *hSample = fFastFit != 2 ? NBDhist(mu,k) : 0;
- // fNBD->SetParameters(mu,k);
- static TH1D *h1 = (TH1D*)hDATA->Clone();
- h1->Reset();
- // h1->SetName(Form("fit_%.3f_%.3f_%.3f_%.3f",mu,k,eff,alpha));
TFile *outFile = NULL;
TNtuple *ntuple = NULL;
if (save) {
- outFile = TFile::Open("pippo.root", "RECREATE");
+ outFile = new TFile(fOutntuplename,"RECREATE");
ntuple = new TNtuple("gnt", "Glauber ntuple", "Npart:Ncoll:B:tAA:ntot");
}
-
- Int_t nents = 100000;//glau_ntuple->GetEntries();
- //if(fFastFit > 0 && fFastFit < 2) nents = 10000;
- // Int_t nents = 100;//glau_ntuple->GetEntries();
- for (Int_t i=0;i<nents;++i) {
- glau_ntuple->GetEntry(i);
- //Int_t n = TMath::Nint(TMath::Power(Npart,alpha));
- //Int_t n = TMath::Nint(TMath::Power(Ncoll,alpha));
- Int_t n = alpha * Npart + (1-alpha) * Ncoll;
+ Int_t nents = 0;
+ if (fGlauntuple)
+ nents = fGlauntuple->GetEntries();
+
+ for (Int_t i=0;i<fNevents;++i) {
+ if (fGlauntuple)
+ fGlauntuple->GetEntry(i % nents);
+ else {
+ fNpart = 2;
+ fNcoll = 1;
+ }
+ Int_t n=0;
+ //if (fAncestor == 1) n = (Int_t)(TMath::Power(fNpart,alpha));
+ if (fAncestor == 1) n = (Int_t)(TMath::Power(fNcoll,alpha));
+ else if (fAncestor == 2) n = (Int_t)(alpha * fNpart + (1-alpha) * fNcoll);
+ else if (fAncestor == 3) n = (Int_t)((1-alpha) * fNpart/2 + alpha * fNcoll);
+
Int_t ntot=0;
- // ntot=n*fNBD->GetRandom();
if (fFastFit == 1) {
- // NBD
- ntot = n*hSample->GetRandom();
+ ntot = (Int_t)(n*hSample->GetRandom()); // NBD
}
else if (fFastFit == 2) {
- // Gaussian
- Double_t sigma = k*TMath::Sqrt(n*mu);
- ntot = gRandom->Gaus(n*mu,sigma);
+ Double_t sigma = k*TMath::Sqrt(n*mu);
+ ntot = (Int_t)(gRandom->Gaus(n*mu,sigma)); // Gaussian
}
else {
- for(Int_t j = 0; j<n; ++j)
- ntot+=hSample->GetRandom();
- // ntot+=fNBD->GetRandom();
+ for(Int_t j = 0; j<(Int_t)n; ++j)
+ ntot += (Int_t)hSample->GetRandom();
}
- h1->Fill(ntot);
-
+ h1->Fill(ntot);
+
if (save)
- ntuple->Fill(Npart,Ncoll,B,tAA,ntot);
-
+ ntuple->Fill(fNpart,fNcoll,fB,fTaa,ntot);
}
if (save) {
if (fFastFit != 2) delete hSample;
return h1;
-
}
//--------------------------------------------------------------------------------------------------
-
-Float_t AliCentralityGlauberFit::CalculateChi2(TH1D *hDATA, TH1D *thistGlau, Float_t eff) {
- // note that for different values of neff the mc distribution is identical, just re-scaled below
- // normalize the two histogram
- // scale MC up but leave data alone - that preserves error bars !!!!
+Double_t AliCentralityGlauberFit::CalculateChi2(TH1F *hDATA, TH1F *thistGlau)
+{
+ // Note that for different values of neff the mc distribution is identical, just re-scaled below
+ // normalize the two histogram, scale MC up but leave data alone - that preserves error bars !!!!
- int lowchibin = hDATA->FindBin(multmin);
- int highchibin = hDATA->FindBin(multmax);
+ Int_t lowchibin = hDATA->FindBin(fMultmin);
+ Int_t highchibin = hDATA->FindBin(fMultmax);
- float mcintegral = thistGlau->Integral(lowchibin,mcintegral);
- // float scale = (hDATA->Integral(lowchibin,highchibin)/mcintegral) * ((float) eff);
- float scale = (hDATA->Integral(lowchibin,highchibin)/mcintegral) * ((float) eff);
+ Double_t mcintegral = thistGlau->Integral(1,thistGlau->GetNbinsX());
+ Double_t scale = (hDATA->Integral(1,hDATA->GetNbinsX())/mcintegral);
thistGlau->Scale(scale);
- // FIXME: Kolmogorov
- //return hDATA->KolmogorovTest(thistGlau,"M");
-
- // // calculate the chi2 between MC and real data over some range ????
+ // calculate the chi2 between MC and real data over some range ????
if (fUseChi2) {
- float chi2 = 0.0;
- for (int i=lowchibin; i<=highchibin; i++) {
+ Double_t chi2 = 0.0;
+ for (Int_t i=lowchibin; i<=highchibin; i++) {
if (hDATA->GetBinContent(i) < 1.0) continue;
- float diff = pow( (thistGlau->GetBinContent(i) - hDATA->GetBinContent(i)) , 2);
- diff = diff / (pow(hDATA->GetBinError(i) , 2)+pow(thistGlau->GetBinError(i) , 2)); // FIXME squared distance if commented
+ Double_t diff = TMath::Power((thistGlau->GetBinContent(i) - hDATA->GetBinContent(i)),2);
+ diff = diff / (TMath::Power(hDATA->GetBinError(i),2) + TMath::Power(thistGlau->GetBinError(i),2));
chi2 += diff;
}
chi2 = chi2 / (highchibin - lowchibin + 1);
}
// "-2 log likelihood ratio(mu;n) = 2[(mu - n) + n ln(n/mu)]"
else {
- cout << "LL" << endl;
-
+ std::cout << "LL" << std::endl;
- float ll = 0.0;
- for (int i=lowchibin; i<=highchibin; i++) {
- // if (thistGlau->GetBinContent(i) < 1) continue;
- // if (hDATA->GetBinContent(i) < 1) continue;
- // cout << ll << " " << thistGlau->GetBinContent(i) <<" " << hDATA->GetBinContent(i) << endl;
-
- Float_t data = hDATA ->GetBinContent(i);
- Float_t mc = thistGlau->GetBinContent(i);
+ Double_t ll = 0.0;
+ for (Int_t i=lowchibin; i<=highchibin; i++) {
+ Double_t data = hDATA ->GetBinContent(i);
+ Double_t mc = thistGlau->GetBinContent(i);
Int_t idata = TMath::Nint(data);
if (mc < 1.e-9) mc = 1.e-9;
- Float_t fsub = - mc + idata * TMath::Log(mc);
- Float_t fobs = 0;
- Int_t imc = TMath::Nint(mc);
+ Double_t fsub = - mc + idata * TMath::Log(mc);
+ Double_t fobs = 0;
if (idata > 0) {
for(Int_t istep = 0; istep < idata; istep++){
if (istep > 1)
fobs += TMath::Log(istep);
}
}
- // cout << mc << " " << data << " " << fsub << " " << fobs << endl;
fsub -= fobs;
ll -= fsub ;
}
return 2*ll;
}
-
}
//--------------------------------------------------------------------------------------------------
+TH1F *AliCentralityGlauberFit::GetTriggerEfficiencyFunction(TH1F *hist1, TH1F *hist2)
+{
+ // Get efficiency.
-TH1D * AliCentralityGlauberFit::GetTriggerEfficiencyFunction(TH1D *hist1, TH1D *hist2) {
- heffi = (TH1D*)hist1->Clone("heffi");
- heffi->Divide(hist2);
- return heffi;
+ fhEffi = (TH1F*)hist1->Clone("heffi");
+ fhEffi->Divide(hist2);
+ return fhEffi;
}
//--------------------------------------------------------------------------------------------------
+Double_t AliCentralityGlauberFit::GetTriggerEfficiencyIntegral(TH1F *hist1, TH1F *hist2)
+{
+ // Get eff integral.
-Float_t AliCentralityGlauberFit::GetTriggerEfficiencyIntegral(TH1D *hist1, TH1D *hist2) {
- effi = hist1->Integral()/hist2->Integral();
- return effi;
+ fEffi = hist1->Integral()/hist2->Integral();
+ return fEffi;
}
//--------------------------------------------------------------------------------------------------
+void AliCentralityGlauberFit::SaveHisto(TH1F *hist1, TH1F *hist2, TH1F *heffi, TFile *outrootfile)
+{
+ // Save histograms.
-void AliCentralityGlauberFit::SaveHisto(TH1D *hist1, TH1D *hist2, TH1D *heffi, TFile *outrootfile) {
outrootfile->cd();
hist1->Write();
hist2->Write();
}
//--------------------------------------------------------------------------------------------------
-
-Double_t AliCentralityGlauberFit::NBD(Int_t n, Double_t mu, Double_t k)
+Double_t AliCentralityGlauberFit::NBD(Int_t n, Double_t mu, Double_t k) const
{
+ // Compute NBD.
+ double F;
+ double f;
+
+ if (n+k > 100.0) {
+ // log method for handling large numbers
+ F = TMath::LnGamma(n + k)- TMath::LnGamma(n + 1.)- TMath::LnGamma(k);
+ f = n * TMath::Log(mu/k) - (n + k) * TMath::Log(1.0 + mu/k);
+ F = F+f;
+ F = TMath::Exp(F);
+ } else {
+ F = TMath::Gamma(n + k) / ( TMath::Gamma(n + 1.) * TMath::Gamma(k) );
+ f = n * TMath::Log(mu/k) - (n + k) * TMath::Log(1.0 + mu/k);
+ f = TMath::Exp(f);
+ F *= f;
+ }
+
+ return F;
- Double_t mudk = mu/k;
- //Double_t ret = TMath::Gamma(n+k) / TMath::Gamma(k) / TMath::Factorial(n) *
- // TMath::Power(mudk,n) / TMath::Power(1+mudk,n+k);
- 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);
- return ret;
}
//--------------------------------------------------------------------------------------------------
-
-TH1D *AliCentralityGlauberFit::NBDhist(Double_t mu, Double_t k)
+TH1F *AliCentralityGlauberFit::NBDhist(Double_t mu, Double_t k)
{
- TH1D *h = new TH1D("htemp","",100,-0.5,299.5);
+ // Interface for TH1F.
+
+ TH1F *h = new TH1F("htemp","",100,-0.5,299.5);
h->SetName(Form("nbd_%f_%f",mu,k));
h->SetDirectory(0);
for (Int_t i=0;i<300;++i) {
return h;
}
+//--------------------------------------------------------------------------------------------------
+void AliCentralityGlauberFit::MinuitFcnNBD(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
+{
+ // FCN for minuit.
-void AliCentralityGlauberFit::MinuitFcnNBD(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag){
- // FCN for minuit
Double_t alpha = par[0];
Double_t mu = par[1];
Double_t k = par[2];
- Double_t eff = par[3];
- //Double_t eff = 1;//par[3];
+ if (0) { //avoid warning
+ gin=gin;
+ npar=npar;
+ iflag=iflag;
+ }
AliCentralityGlauberFit * obj = (AliCentralityGlauberFit *) gMinuit->GetObjectFit();
-
- // static TStopwatch sw;
- // sw.Start();
- TH1D * thistGlau = obj->GlauberHisto(mu,k,eff,alpha,obj->GetTempHist(),kFALSE);
- f = obj->CalculateChi2(obj->GetTempHist(),thistGlau,eff);
- // sw.Stop();
- // sw.Print();
-
- Printf("%f - %f - %f - %f - %f \n",f,alpha,mu,k,eff);
-
-
+ TH1F * thistGlau = obj->GlauberHisto(mu,k,alpha,obj->GetTempHist(),kFALSE);
+ f = obj->CalculateChi2(obj->GetTempHist(),thistGlau);
+ Printf("Minuit step: chi2=%f, alpha=%f, mu=%f, k=%f\n",f,alpha,mu,k);
}
-Double_t AliCentralityGlauberFit::NBDFunc(Double_t * x, Double_t *par) {
- // TF1 interface
+//--------------------------------------------------------------------------------------------------
+Double_t AliCentralityGlauberFit::NBDFunc(const Double_t *x, const Double_t *par)
+{
+ // TF1 interface.
Double_t mu = par[0];
Double_t k = par[1];
Double_t n = x[0];
- Double_t mudk = mu/k;
- //Double_t ret = TMath::Gamma(n+k) / TMath::Gamma(k) / TMath::Factorial(n) *
- // TMath::Power(mudk,n) / TMath::Power(1+mudk,n+k);
- 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);
-
+ Double_t ret = exp( TMath::LnGamma(n+k) - TMath::LnGamma(k) - TMath::LnGamma(n+1) ) *
+ TMath::Power(mu/(mu+k),n) * TMath::Power(1-mu/(mu+k),k);
return ret;
+}
+
+//--------------------------------------------------------------------------------------------------
+TH1F *AliCentralityGlauberFit::MakeAncestor(Double_t alpha)
+{
+ // Make ancestor histogram.
+ TString hname(Form("fhAncestor_%.3f",alpha));
+ if (fhAncestor) {
+ if (hname.CompareTo(fhAncestor->GetName())==0)
+ return fhAncestor;
+ }
+
+ delete fhAncestor;
+ fhAncestor = 0;
+ TFile *ancfile = TFile::Open(fAncfilename,"read");
+ if (ancfile && ancfile->IsOpen()) {
+ fhAncestor = dynamic_cast<TH1F*>(ancfile->Get(hname));
+ if (fhAncestor) {
+ fhAncestor->SetDirectory(0);
+ delete ancfile;
+ return fhAncestor;
+ }
+ }
+ delete ancfile;
+
+ fhAncestor = new TH1F(hname,hname,3000,0,3000);
+ fhAncestor->SetDirectory(0);
+ Int_t nents = fGlauntuple->GetEntries();
+ for (Int_t i=0;i<nents;++i) {
+ fGlauntuple->GetEntry(i % nents);
+ Int_t n=0;
+ //if (fAncestor == 1) n = (Int_t) (TMath::Power(fNpart,alpha));
+ if (fAncestor == 1) n = (Int_t) (TMath::Power(fNcoll,alpha));
+ else if (fAncestor == 2) n = (Int_t) (alpha * fNpart + (1-alpha) * fNcoll);
+ else if (fAncestor == 3) n = (Int_t) ((1-alpha) * fNpart/2 + alpha * fNcoll);
+ fhAncestor->Fill(n);
+ }
+
+ ancfile = TFile::Open(fAncfilename,"update");
+ if (ancfile && ancfile->IsOpen()) {
+ fhAncestor->Write();
+ }
+ delete ancfile;
+ return fhAncestor;
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff