From 728661538134afa86ba54426a57b09daeb74ee6a Mon Sep 17 00:00:00 2001
From: Viktor Klochkov <v.klochkov@gsi.de>
Date: Fri, 19 Jul 2019 11:16:57 +0200
Subject: [PATCH] fix alignment & add openmp
---
CMakeLists.txt | 2 +-
glauber/Fitter.cpp | 538 ++++++++++++++++++++---------------------
glauber/Fitter.h | 174 ++++++-------
glauber/FitterHelper.h | 146 +++++------
glauber/main.cpp | 119 ++++-----
5 files changed, 494 insertions(+), 485 deletions(-)
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 359f4c8..04f3815 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -21,8 +21,8 @@ if (CMAKE_BUILD_TYPE MATCHES DEBUG)
endif ()
set(CMAKE_CXX_FLAGS_DEBUG "-O0 -ggdb -DDEBUG -D__DEBUG -Wall")
+set(CMAKE_CXX_FLAGS_RELEASE "-O3 -march=native -ftree-vectorize -ffast-math -DNODEBUG -fopenmp")
set(CMAKE_CXX_FLAGS_RELWITHDEBINFO "${CMAKE_CXX_FLAGS_RELEASE} -ggdb")
-set(CMAKE_CXX_FLAGS_RELEASE "-O3 -march=native -ftree-vectorize -ffast-math -DNODEBUG")
message(STATUS "Using CXX flags for ${CMAKE_BUILD_TYPE}: ${CMAKE_CXX_FLAGS_${CMAKE_BUILD_TYPE}}")
set(SOURCES
diff --git a/glauber/Fitter.cpp b/glauber/Fitter.cpp
index f8c54e9..91e86e2 100644
--- a/glauber/Fitter.cpp
+++ b/glauber/Fitter.cpp
@@ -13,78 +13,78 @@ ClassImp(Glauber::Fitter)
// ----- Default constructor -------------------------------------------
Glauber::Fitter::Fitter(std::unique_ptr<TTree> tree)
{
- fSimTree = std::move(tree);
- std::cout << fSimTree->GetEntries() << std::endl;
-
- if (!fSimTree) {
- std::cout << "SetSimHistos: *** Error - " << std::endl;
- exit(EXIT_FAILURE);
- }
-
- fSimTree->SetBranchAddress("Npart", &fNpart);
- fSimTree->SetBranchAddress("Ncoll", &fNcoll);
+ fSimTree = std::move(tree);
+ std::cout << fSimTree->GetEntries() << std::endl;
+
+ if (!fSimTree) {
+ std::cout << "SetSimHistos: *** Error - " << std::endl;
+ exit(EXIT_FAILURE);
+ }
+
+ fSimTree->SetBranchAddress("Npart", &fNpart);
+ fSimTree->SetBranchAddress("Ncoll", &fNcoll);
}
void Glauber::Fitter::Init(int nEntries)
{
- if ( nEntries < 0 || nEntries > fSimTree->GetEntries() ){
- std::cout << "Init: *** ERROR - number of entries < 0 or less that number of entries in input tree" << std::endl;
- std::cout << "Init: *** number of entries in input tree = " << fSimTree->GetEntries() << std::endl;
- exit(EXIT_FAILURE);
- }
-
- const int NpartMax = int (fSimTree->GetMaximum("Npart") );
- const int NcollMax = int (fSimTree->GetMaximum("Ncoll") );
-
- fNpartHisto = TH1F ("fNpartHisto", "Npart", NpartMax/fBinSize, 0, NpartMax );
- fNcollHisto = TH1F ("fNcollHisto", "Ncoll", NcollMax/fBinSize, 0, NcollMax );
-
- for (int i=0; i<nEntries; i++)
- {
- fSimTree->GetEntry(i);
- fNcollHisto.Fill(fNcoll);
- fNpartHisto.Fill(fNpart);
- }
- std::cout << fSimTree->GetEntries() << std::endl;
-
- fNbins = fDataHisto.GetNbinsX();
-
- while ( fDataHisto.GetBinContent(fNbins-1) == 0)
- fNbins--;
-
- fNbins++;
-
- const float min = fDataHisto.GetXaxis()->GetXmin();
- const float max = fDataHisto.GetXaxis()->GetXmax();
-
- fMaxValue = min + (max - min)*fNbins/fDataHisto.GetNbinsX() ;
-
- std::cout << "fNbins = " << fNbins << std::endl;
- std::cout << "fMaxValue = " << fMaxValue << std::endl;
-
+ if ( nEntries < 0 || nEntries > fSimTree->GetEntries() ){
+ std::cout << "Init: *** ERROR - number of entries < 0 or less that number of entries in input tree" << std::endl;
+ std::cout << "Init: *** number of entries in input tree = " << fSimTree->GetEntries() << std::endl;
+ exit(EXIT_FAILURE);
+ }
+
+ const int NpartMax = int (fSimTree->GetMaximum("Npart") );
+ const int NcollMax = int (fSimTree->GetMaximum("Ncoll") );
+
+ fNpartHisto = TH1F ("fNpartHisto", "Npart", NpartMax/fBinSize, 0, NpartMax );
+ fNcollHisto = TH1F ("fNcollHisto", "Ncoll", NcollMax/fBinSize, 0, NcollMax );
+
+ for (int i=0; i<nEntries; i++)
+ {
+ fSimTree->GetEntry(i);
+ fNcollHisto.Fill(fNcoll);
+ fNpartHisto.Fill(fNpart);
+ }
+ std::cout << fSimTree->GetEntries() << std::endl;
+
+ fNbins = fDataHisto.GetNbinsX();
+
+ while ( fDataHisto.GetBinContent(fNbins-1) == 0)
+ fNbins--;
+
+ fNbins++;
+
+ const float min = fDataHisto.GetXaxis()->GetXmin();
+ const float max = fDataHisto.GetXaxis()->GetXmax();
+
+ fMaxValue = min + (max - min)*fNbins/fDataHisto.GetNbinsX() ;
+
+ std::cout << "fNbins = " << fNbins << std::endl;
+ std::cout << "fMaxValue = " << fMaxValue << std::endl;
+
}
float Glauber::Fitter::Nancestors(float f) const
{
- if (fMode == "Default") return f*fNpart + (1-f)*fNcoll;
- else if (fMode == "PSD") return f - fNpart;
- else if (fMode == "Npart") return TMath::Power(fNpart, f);
- else if (fMode == "Ncoll") return TMath::Power(fNcoll, f);
-
- return -1.;
+ if (fMode == "Default") return f*fNpart + (1-f)*fNcoll;
+ else if (fMode == "PSD") return f - fNpart;
+ else if (fMode == "Npart") return TMath::Power(fNpart, f);
+ else if (fMode == "Ncoll") return TMath::Power(fNcoll, f);
+
+ return -1.;
}
float Glauber::Fitter::NancestorsMax(float f) const
{
- const int NpartMax = fNpartHisto.GetXaxis()->GetXmax() ; // some magic
- const int NcollMax = fNcollHisto.GetXaxis()->GetXmax() ; //TODO
-
- if (fMode == "Default") return f*NpartMax + (1-f)*NcollMax;
- else if (fMode == "PSD") return NpartMax;
- else if (fMode == "Npart") return TMath::Power(NpartMax, f);
- else if (fMode == "Ncoll") return TMath::Power(NcollMax, f);
-
- return -1.;
+ const int NpartMax = fNpartHisto.GetXaxis()->GetXmax() ; // some magic
+ const int NcollMax = fNcollHisto.GetXaxis()->GetXmax() ; //TODO
+
+ if (fMode == "Default") return f*NpartMax + (1-f)*NcollMax;
+ else if (fMode == "PSD") return NpartMax;
+ else if (fMode == "Npart") return TMath::Power(NpartMax, f);
+ else if (fMode == "Ncoll") return TMath::Power(NcollMax, f);
+
+ return -1.;
}
/*
@@ -93,52 +93,52 @@ float Glauber::Fitter::NancestorsMax(float f) const
*/
void Glauber::Fitter::SetGlauberFitHisto (float f, float mu, float k, int n, Bool_t Norm2Data)
{
- fGlauberFitHisto = TH1F("glaub", "", fNbins*1.3, 0, 1.3*fMaxValue);
- fGlauberFitHisto.SetName(Form("glaub_%4.2f_%6.4f_%4.2f_%d", f, mu, k, n ));
-
- SetNBDhist(mu, k);
-
- std::unique_ptr<TH1F> htemp {(TH1F*)fNbdHisto.Clone("htemp")}; // WTF??? Not working without pointer
- for (int i=0; i<n; i++)
- {
- fSimTree->GetEntry(i);
- const int Na = int(Nancestors(f));
-
- float nHits {0.};
- for (int j=0; j<Na; j++){
- nHits += int(htemp->GetRandom());
- }
- fGlauberFitHisto.Fill(nHits);
+ fGlauberFitHisto = TH1F("glaub", "", fNbins*1.3, 0, 1.3*fMaxValue);
+ fGlauberFitHisto.SetName(Form("glaub_%4.2f_%6.4f_%4.2f_%d", f, mu, k, n ));
+
+ SetNBDhist(mu, k);
+
+ std::unique_ptr<TH1F> htemp {(TH1F*)fNbdHisto.Clone("htemp")}; // WTF??? Not working without pointer
+ for (int i=0; i<n; i++)
+ {
+ fSimTree->GetEntry(i);
+ const int Na = int(Nancestors(f));
+
+ float nHits {0.};
+ for (int j=0; j<Na; j++){
+ nHits += int(htemp->GetRandom());
}
-
- if (Norm2Data)
- NormalizeGlauberFit();
+ fGlauberFitHisto.Fill(nHits);
+ }
+
+ if (Norm2Data)
+ NormalizeGlauberFit();
}
void Glauber::Fitter::NormalizeGlauberFit ()
{
-
- int fGlauberFitHistoInt {0};
- int fDataHistoInt {0};
-
- const int lowchibin = fFitMinBin;
- const int highchibin = fFitMaxBin<fNbins ? fFitMaxBin : fNbins;
-
- for (int i=lowchibin; i<highchibin; i++)
- {
- fGlauberFitHistoInt += fGlauberFitHisto.GetBinContent(i+1);
- fDataHistoInt += fDataHisto.GetBinContent(i+1);
- }
-
- const float ScaleFactor = (float)fDataHistoInt/fGlauberFitHistoInt;
-
-// std::cout << "Scale = " << Scale << std::endl;
- fGlauberFitHisto.Scale(ScaleFactor);
+
+ int fGlauberFitHistoInt {0};
+ int fDataHistoInt {0};
+
+ const int lowchibin = fFitMinBin;
+ const int highchibin = fFitMaxBin<fNbins ? fFitMaxBin : fNbins;
+
+ for (int i=lowchibin; i<highchibin; i++)
+ {
+ fGlauberFitHistoInt += fGlauberFitHisto.GetBinContent(i+1);
+ fDataHistoInt += fDataHisto.GetBinContent(i+1);
+ }
+
+ const float ScaleFactor = (float)fDataHistoInt/fGlauberFitHistoInt;
+
+ // std::cout << "Scale = " << Scale << std::endl;
+ fGlauberFitHisto.Scale(ScaleFactor);
}
/**
- *
+ *
* @param mu mean value of negative binominal distribution (we are looking for it)
* @param chi2 return value (indicates good match)
* @param mu_min lower search edge for mean value NBD
@@ -150,48 +150,48 @@ void Glauber::Fitter::NormalizeGlauberFit ()
*/
void Glauber::Fitter::FindMuGoldenSection (float *mu, float *chi2, float mu_min, float mu_max, float f, float k, int nEvents, int nIter)
{
- const float phi {(1+TMath::Sqrt(5))/2};
-
- /* left */
- float mu_1 = mu_max - (mu_max-mu_min)/phi;
-
- /* right */
- float mu_2 = mu_min + (mu_max-mu_min)/phi;
-
- SetGlauberFitHisto (f, mu_1, k, nEvents);
- float chi2_mu1 = GetChi2 ();
-
- SetGlauberFitHisto (f, mu_2, k, nEvents);
- float chi2_mu2 = GetChi2 ();
-
- for (int j=0; j<nIter; j++)
- {
- if (chi2_mu1 > chi2_mu2)
- {
- mu_min = mu_1;
- mu_1 = mu_2;
- mu_2 = mu_min + (mu_max-mu_min)/phi;
- chi2_mu1 = chi2_mu2;
- SetGlauberFitHisto (f, mu_2, k, nEvents);
- chi2_mu2 = GetChi2 ();
- }
- else
- {
- mu_max = mu_2;
- mu_2 = mu_1;
- mu_1 = mu_max - (mu_max-mu_min)/phi;
- chi2_mu2 = chi2_mu1;
- SetGlauberFitHisto (f, mu_1, k, nEvents);
- chi2_mu1 = GetChi2 ();
- }
-
- std::cout << "mu1 = " << mu_1 << " mu2 = " << mu_2 << " chi2_mu1 = " << chi2_mu1 << " chi2_mu2 = " << chi2_mu2 << std::endl;
+ const float phi {(1+TMath::Sqrt(5))/2};
+
+ /* left */
+ float mu_1 = mu_max - (mu_max-mu_min)/phi;
+
+ /* right */
+ float mu_2 = mu_min + (mu_max-mu_min)/phi;
+
+ SetGlauberFitHisto (f, mu_1, k, nEvents);
+ float chi2_mu1 = GetChi2 ();
+
+ SetGlauberFitHisto (f, mu_2, k, nEvents);
+ float chi2_mu2 = GetChi2 ();
+
+ for (int j=0; j<nIter; j++)
+ {
+ if (chi2_mu1 > chi2_mu2)
+ {
+ mu_min = mu_1;
+ mu_1 = mu_2;
+ mu_2 = mu_min + (mu_max-mu_min)/phi;
+ chi2_mu1 = chi2_mu2;
+ SetGlauberFitHisto (f, mu_2, k, nEvents);
+ chi2_mu2 = GetChi2 ();
}
-
- /* take min(mu) */
- *mu = (chi2_mu1 < chi2_mu2) ? mu_1 : mu_2;
- /* take min(chi2) */
- *chi2 = (chi2_mu1 < chi2_mu2) ? chi2_mu1 : chi2_mu2;
+ else
+ {
+ mu_max = mu_2;
+ mu_2 = mu_1;
+ mu_1 = mu_max - (mu_max-mu_min)/phi;
+ chi2_mu2 = chi2_mu1;
+ SetGlauberFitHisto (f, mu_1, k, nEvents);
+ chi2_mu1 = GetChi2 ();
+ }
+
+ std::cout << "mu1 = " << mu_1 << " mu2 = " << mu_2 << " chi2_mu1 = " << chi2_mu1 << " chi2_mu2 = " << chi2_mu2 << std::endl;
+ }
+
+ /* take min(mu) */
+ *mu = (chi2_mu1 < chi2_mu2) ? mu_1 : mu_2;
+ /* take min(chi2) */
+ *chi2 = (chi2_mu1 < chi2_mu2) ? chi2_mu1 : chi2_mu2;
}
/**
@@ -205,63 +205,63 @@ void Glauber::Fitter::FindMuGoldenSection (float *mu, float *chi2, float mu_min,
*/
float Glauber::Fitter::FitGlauber (float *par, Float_t f0, Int_t k0, Int_t k1, Int_t nEvents)
{
- float f_fit{-1};
- float mu_fit{-1};
- float k_fit{-1};
- float Chi2Min {1e10};
-
- const TString filename = Form ( "%s/fit_%4.2f_%d_%d_%d.root", fOutDirName.Data(), f0, k0, k1, fFitMinBin );
+ float f_fit{-1};
+ float mu_fit{-1};
+ float k_fit{-1};
+ float Chi2Min {1e10};
+
+ const TString filename = Form ( "%s/fit_%4.2f_%d_%d_%d.root", fOutDirName.Data(), f0, k0, k1, fFitMinBin );
+
+ // std::unique_ptr<TFile> file {TFile::Open(filename, "recreate")};
+ // std::unique_ptr<TTree> tree {new TTree("test_tree", "tree" )};
+
+ TFile* file {TFile::Open(filename, "recreate")};
+ TTree* tree {new TTree("test_tree", "tree" )};
+
+ TH1F h1("h1", "", fNbins, 0, fMaxValue);
+
+ float f, mu, k, chi2, sigma;
+
+ tree->Branch("histo", "TH1F", &h1);
+ tree->Branch("f", &f, "f/F");
+ tree->Branch("mu", &mu, "mu/F");
+ tree->Branch("k", &k, "k/F");
+ tree->Branch("chi2", &chi2, "chi2/F");
+ tree->Branch("sigma",&sigma,"sigma/F");
+
+ f = f0;
+ for (int j=k0; j<=k1; j++)
+ {
+ mu = fMaxValue / NancestorsMax(f) ;
+ k = j;
+ const float mu_min = 0.7*mu;
+ const float mu_max = 1.0*mu;
-// std::unique_ptr<TFile> file {TFile::Open(filename, "recreate")};
-// std::unique_ptr<TTree> tree {new TTree("test_tree", "tree" )};
-
- TFile* file {TFile::Open(filename, "recreate")};
- TTree* tree {new TTree("test_tree", "tree" )};
+ FindMuGoldenSection (&mu, &chi2, mu_min, mu_max, f, k, nEvents, 10);
+ sigma = ( mu/k + 1 ) * mu;
+ h1 = fGlauberFitHisto;
- TH1F h1("h1", "", fNbins, 0, fMaxValue);
-
- float f, mu, k, chi2, sigma;
-
- tree->Branch("histo", "TH1F", &h1);
- tree->Branch("f", &f, "f/F");
- tree->Branch("mu", &mu, "mu/F");
- tree->Branch("k", &k, "k/F");
- tree->Branch("chi2", &chi2, "chi2/F");
- tree->Branch("sigma",&sigma,"sigma/F");
-
- f = f0;
- for (int j=k0; j<=k1; j++)
- {
- mu = fMaxValue / NancestorsMax(f) ;
- k = j;
- const float mu_min = 0.7*mu;
- const float mu_max = 1.0*mu;
-
- FindMuGoldenSection (&mu, &chi2, mu_min, mu_max, f, k, nEvents, 10);
- sigma = ( mu/k + 1 ) * mu;
- h1 = fGlauberFitHisto;
-
- tree->Fill();
-
- if (chi2 < Chi2Min)
- {
- f_fit = f;
- mu_fit = mu;
- k_fit = k;
- Chi2Min = chi2;
- fBestFitHisto = fGlauberFitHisto;
- }
-
- }
- tree->Write();
- file->Write();
- file->Close();
-
- par[0] = f_fit;
- par[1] = mu_fit;
- par[2] = k_fit;
+ tree->Fill();
- return Chi2Min;
+ if (chi2 < Chi2Min)
+ {
+ f_fit = f;
+ mu_fit = mu;
+ k_fit = k;
+ Chi2Min = chi2;
+ fBestFitHisto = fGlauberFitHisto;
+ }
+
+ }
+ tree->Write();
+ file->Write();
+ file->Close();
+
+ par[0] = f_fit;
+ par[1] = mu_fit;
+ par[2] = k_fit;
+
+ return Chi2Min;
}
/**
@@ -270,22 +270,22 @@ float Glauber::Fitter::FitGlauber (float *par, Float_t f0, Int_t k0, Int_t k1, I
*/
float Glauber::Fitter::GetChi2 () const
{
- float chi2 {0.0};
-
- const int lowchibin = fFitMinBin;
- const int highchibin = fFitMaxBin<fNbins ? fFitMaxBin : fNbins;
-
- for (int i=lowchibin; i<=highchibin; ++i)
- {
- if (fDataHisto.GetBinContent(i) < 1.0) continue;
- const float error2 = TMath::Power(fDataHisto.GetBinError(i), 2) + TMath::Power(fGlauberFitHisto.GetBinError(i), 2);
- const float diff2 = TMath::Power((fGlauberFitHisto.GetBinContent(i) - fDataHisto.GetBinContent(i)), 2) / error2;
-
- chi2 += diff2;
- }
-
- chi2 = chi2 / (highchibin - lowchibin + 1);
- return chi2;
+ float chi2 {0.0};
+
+ const int lowchibin = fFitMinBin;
+ const int highchibin = fFitMaxBin<fNbins ? fFitMaxBin : fNbins;
+
+ for (int i=lowchibin; i<=highchibin; ++i)
+ {
+ if (fDataHisto.GetBinContent(i) < 1.0) continue;
+ const float error2 = TMath::Power(fDataHisto.GetBinError(i), 2) + TMath::Power(fGlauberFitHisto.GetBinError(i), 2);
+ const float diff2 = TMath::Power((fGlauberFitHisto.GetBinContent(i) - fDataHisto.GetBinContent(i)), 2) / error2;
+
+ chi2 += diff2;
+ }
+
+ chi2 = chi2 / (highchibin - lowchibin + 1);
+ return chi2;
}
/**
@@ -295,18 +295,18 @@ float Glauber::Fitter::GetChi2 () const
*/
void Glauber::Fitter::SetNBDhist(float mu, float k)
{
- // Interface for TH1F.
- const int nBins = (mu+1.)*3 < 10 ? 10 : (mu+1.)*3;
-
- fNbdHisto = TH1F ("fNbdHisto", "", nBins, 0, nBins);
- fNbdHisto.SetName(Form("nbd_%f_%f",mu,k));
-
- for (int i=0; i<nBins; ++i)
- {
- const float val = NBD(i, mu, k);
- if (val>1e-10) fNbdHisto.SetBinContent(i+1, val);
-// std::cout << "val " << val << std::endl;
- }
+ // Interface for TH1F.
+ const int nBins = (mu+1.)*3 < 10 ? 10 : (mu+1.)*3;
+
+ fNbdHisto = TH1F ("fNbdHisto", "", nBins, 0, nBins);
+ fNbdHisto.SetName(Form("nbd_%f_%f",mu,k));
+
+ for (int i=0; i<nBins; ++i)
+ {
+ const float val = NBD(i, mu, k);
+ if (val>1e-10) fNbdHisto.SetBinContent(i+1, val);
+ // std::cout << "val " << val << std::endl;
+ }
}
/**
@@ -318,27 +318,27 @@ void Glauber::Fitter::SetNBDhist(float mu, float k)
*/
float Glauber::Fitter::NBD(float n, float mu, float k) const
{
- // Compute NBD.
- float F;
- float 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 = 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;
+ // Compute NBD.
+ float F;
+ float 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 = 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;
}
/**
* Creates histo with a given model parameter distribution
@@ -350,37 +350,37 @@ float Glauber::Fitter::NBD(float n, float mu, float k) const
*/
std::unique_ptr<TH1F> Glauber::Fitter::GetModelHisto (const float range[2], TString name, const float par[3], int nEvents)
{
- const float f = par[0];
- const float mu = par[1];
- const float k = par[2];
-
- float modelpar{-999.};
- fSimTree->SetBranchAddress(name, &modelpar);
-
- SetNBDhist(mu, k);
+ const float f = par[0];
+ const float mu = par[1];
+ const float k = par[2];
-// TRandom random;
-// random.SetSeed(mu*k);
-
- std::unique_ptr<TH1F> hModel(new TH1F ("hModel", "name", 100, fSimTree->GetMinimum(name), fSimTree->GetMaximum(name)) );
-
- for (int i=0; i<nEvents; i++)
- {
- fSimTree->GetEntry(i);
- const int Na = int(Nancestors(f));
- float nHits{0.};
- for (int j=0; j<Na; ++j){
- nHits += (int)fNbdHisto.GetRandom();
- }
-
- if ( nHits > range[0] && nHits < range[1] ){
- hModel->Fill(modelpar);
- }
-
+ float modelpar{-999.};
+ fSimTree->SetBranchAddress(name, &modelpar);
+
+ SetNBDhist(mu, k);
+
+ // TRandom random;
+ // random.SetSeed(mu*k);
+
+ std::unique_ptr<TH1F> hModel(new TH1F ("hModel", "name", 100, fSimTree->GetMinimum(name), fSimTree->GetMaximum(name)) );
+
+ #pragma omp parallel for
+ for (int i=0; i<nEvents; i++)
+ {
+ fSimTree->GetEntry(i);
+ const int Na = int(Nancestors(f));
+ float nHits{0.};
+ for (int j=0; j<Na; ++j){
+ nHits += (int)fNbdHisto.GetRandom();
}
- return std::move(hModel);
-
+ if ( nHits > range[0] && nHits < range[1] ){
+ hModel->Fill(modelpar);
+ }
+ }
+
+ return std::move(hModel);
+
}
diff --git a/glauber/Fitter.h b/glauber/Fitter.h
index 8c73ccb..b0a96ed 100644
--- a/glauber/Fitter.h
+++ b/glauber/Fitter.h
@@ -1,9 +1,9 @@
/** @file Fitter.h
- @class Glauber::Fitter
- @author Viktor Klochkov (klochkov44@gmail.com)
- @author Ilya Selyuzhenkov (ilya.selyuzhenkov@gmail.com)
- @brief Class to fit histo with Glauber based function
-*/
+ * @class Glauber::Fitter
+ * @author Viktor Klochkov (klochkov44@gmail.com)
+ * @author Ilya Selyuzhenkov (ilya.selyuzhenkov@gmail.com)
+ * @brief Class to fit histo with Glauber based function
+ */
#ifndef GlauberFitter_H
#define GlauberFitter_H 1
@@ -18,88 +18,88 @@
namespace Glauber
{
- class Fitter
- {
-
- public:
-
- /** Default constructor **/
- Fitter() {};
- Fitter(std::unique_ptr<TTree> tree) ;
- /** Destructor **/
- virtual ~Fitter(){};
-
- void Init(int nEntries);
- void SetGlauberFitHisto (Float_t f, Float_t mu, Float_t k, Int_t n = 10000, Bool_t Norm2Data = true);
- void NormalizeGlauberFit ();
- void DrawHistos (Bool_t isSim = true, Bool_t isData = true, Bool_t isGlauber = false, Bool_t isNBD = false);
-
- float FitGlauber (float *par, Float_t f0, Int_t k0, Int_t k1, Int_t nEvents);
- void FindMuGoldenSection (Float_t *mu, Float_t *chi2, Float_t mu_min, Float_t mu_max, Float_t f, Float_t k, Int_t nEvents = 10000, Int_t nIter = 5);
-
- Float_t GetChi2 (void) const;
-
- Float_t NBD(Float_t n, Float_t mu, Float_t k) const;
- void SetNBDhist(Float_t mu, Float_t k);
- float Nancestors(float f) const;
- float NancestorsMax(float f) const;
-
- std::unique_ptr<TH1F> GetModelHisto (const Float_t range[2], TString name, const Float_t par[3], Int_t nEvents);
-
-//
-// Setters
-//
- void SetInputHisto (const TH1F &h) { fDataHisto = h; }
- void SetFitMinBin (Int_t min) { fFitMinBin = min; }
- void SetFitMaxBin (Int_t min) { fFitMaxBin = min; }
- void SetNormMinBin (Int_t min) { fNormMinBin = min; }
- void SetBinSize (Int_t size) { fBinSize = size; }
- void SetOutDirName (TString name) { fOutDirName = name; }
- void SetMode (const TString mode) { fMode = mode; }
-
-//
-// Getters
-//
- TH1F GetGlauberFitHisto () const { return fGlauberFitHisto; }
- TH1F GetDataHisto () const { return fDataHisto; }
- TH1F GetNBDHisto () const { return fNbdHisto; }
- TH1F GetNpartHisto () const { return fNpartHisto; }
- TH1F GetNcollHisto () const { return fNcollHisto; }
- TH1F GetBestFiHisto () const { return fBestFitHisto; }
-
-
- private:
-
- /** Data members **/
- TH1F fNpartHisto;
- TH1F fNcollHisto;
- TH1F fDataHisto;
- TH1F fNbdHisto;
- TH1F fGlauberFitHisto;
- TH1F fBestFitHisto;
-
- /* MC data */
- std::unique_ptr<TTree> fSimTree{nullptr};
-
- Float_t fNpart{-1.};
- Float_t fNcoll{-1.};
-
- Float_t fMaxValue{-1.};
-
- Int_t fNbins{-1};
- Int_t fBinSize{-1};
-
- Int_t fFitMinBin{-1};
- Int_t fFitMaxBin{-1};
-
- Int_t fNormMinBin{-1};
-
- TString fMode{"Default"};
-
- TString fOutDirName{""};
- ClassDef(Fitter, 2);
-
- };
+ class Fitter
+ {
+
+ public:
+
+ /** Default constructor **/
+ Fitter() {};
+ Fitter(std::unique_ptr<TTree> tree) ;
+ /** Destructor **/
+ virtual ~Fitter(){};
+
+ void Init(int nEntries);
+ void SetGlauberFitHisto (Float_t f, Float_t mu, Float_t k, Int_t n = 10000, Bool_t Norm2Data = true);
+ void NormalizeGlauberFit ();
+ void DrawHistos (Bool_t isSim = true, Bool_t isData = true, Bool_t isGlauber = false, Bool_t isNBD = false);
+
+ float FitGlauber (float *par, Float_t f0, Int_t k0, Int_t k1, Int_t nEvents);
+ void FindMuGoldenSection (Float_t *mu, Float_t *chi2, Float_t mu_min, Float_t mu_max, Float_t f, Float_t k, Int_t nEvents = 10000, Int_t nIter = 5);
+
+ Float_t GetChi2 (void) const;
+
+ Float_t NBD(Float_t n, Float_t mu, Float_t k) const;
+ void SetNBDhist(Float_t mu, Float_t k);
+ float Nancestors(float f) const;
+ float NancestorsMax(float f) const;
+
+ std::unique_ptr<TH1F> GetModelHisto (const Float_t range[2], TString name, const Float_t par[3], Int_t nEvents);
+
+ //
+ // Setters
+ //
+ void SetInputHisto (const TH1F &h) { fDataHisto = h; }
+ void SetFitMinBin (Int_t min) { fFitMinBin = min; }
+ void SetFitMaxBin (Int_t min) { fFitMaxBin = min; }
+ void SetNormMinBin (Int_t min) { fNormMinBin = min; }
+ void SetBinSize (Int_t size) { fBinSize = size; }
+ void SetOutDirName (TString name) { fOutDirName = name; }
+ void SetMode (const TString mode) { fMode = mode; }
+
+ //
+ // Getters
+ //
+ TH1F GetGlauberFitHisto () const { return fGlauberFitHisto; }
+ TH1F GetDataHisto () const { return fDataHisto; }
+ TH1F GetNBDHisto () const { return fNbdHisto; }
+ TH1F GetNpartHisto () const { return fNpartHisto; }
+ TH1F GetNcollHisto () const { return fNcollHisto; }
+ TH1F GetBestFiHisto () const { return fBestFitHisto; }
+
+
+ private:
+
+ /** Data members **/
+ TH1F fNpartHisto;
+ TH1F fNcollHisto;
+ TH1F fDataHisto;
+ TH1F fNbdHisto;
+ TH1F fGlauberFitHisto;
+ TH1F fBestFitHisto;
+
+ /* MC data */
+ std::unique_ptr<TTree> fSimTree{nullptr};
+
+ Float_t fNpart{-1.};
+ Float_t fNcoll{-1.};
+
+ Float_t fMaxValue{-1.};
+
+ Int_t fNbins{-1};
+ Int_t fBinSize{-1};
+
+ Int_t fFitMinBin{-1};
+ Int_t fFitMaxBin{-1};
+
+ Int_t fNormMinBin{-1};
+
+ TString fMode{"Default"};
+
+ TString fOutDirName{""};
+ ClassDef(Fitter, 2);
+
+ };
}
#endif
diff --git a/glauber/FitterHelper.h b/glauber/FitterHelper.h
index daa705c..4c6c4ad 100644
--- a/glauber/FitterHelper.h
+++ b/glauber/FitterHelper.h
@@ -1,8 +1,8 @@
/** @file FitterHelper.h
- @author Viktor Klochkov (klochkov44@gmail.com)
- @author Ilya Selyuzhenkov (ilya.selyuzhenkov@gmail.com)
- @brief Methods for fit QA
-*/
+ * @author Viktor Klochkov (klochkov44@gmail.com)
+ * @author Ilya Selyuzhenkov (ilya.selyuzhenkov@gmail.com)
+ * @brief Methods for fit QA
+ */
#ifndef GlauberFitterHelper_H
#define GlauberFitterHelper_H 1
@@ -18,79 +18,79 @@
namespace Glauber
{
- inline void DrawHistos (const Fitter& fit, Bool_t isSim, Bool_t isData, Bool_t isGlauber, Bool_t isNBD )
- {
- std::unique_ptr <TCanvas> c1 {new TCanvas("c1", "canvas", 1500, 900)};
-
- c1->Divide(2,2);
-
- std::unique_ptr <TPad> c1_1 {(TPad*) c1->GetListOfPrimitives()->FindObject("c1_1")};
- std::unique_ptr <TPad> c1_2 {(TPad*) c1->GetListOfPrimitives()->FindObject("c1_2")};
- std::unique_ptr <TPad> c1_4 {(TPad*) c1->GetListOfPrimitives()->FindObject("c1_4")};
-
- c1_1->SetLogy(1);
- c1_2->SetLogy(1);
- c1_4->SetLogy(1);
-
- /*const*/ TH1F hGlaub = fit.GetGlauberFitHisto();
- /*const*/ TH1F hData = fit.GetDataHisto();
- /*const*/ TH1F hNBD = fit.GetNBDHisto();
- /*const*/ TH1F hNcoll = fit.GetNcollHisto();
- /*const*/ TH1F hNpart = fit.GetNpartHisto();
- /*const*/ TH1F hBestFit = fit.GetBestFiHisto();
-
- std::unique_ptr <TFile> fOut{TFile::Open("glauber_qa.root", "recreate")};
-
- if (isSim){
- c1->cd(1);
- hNcoll.SetLineColor(2);
-
- hNcoll.Draw();
- hNpart.Draw("same");
-
- std::unique_ptr <TLegend> legSim { new TLegend(0.6,0.75,0.75,0.83) };
- legSim->AddEntry(&hNpart ,"Npart", "l");
- legSim->AddEntry(&hNcoll ,"hNcoll", "l");
- legSim->Draw("same");
-
- hNcoll.Write();
- hNpart.Write();
- }
-
- if (isData){
- c1->cd(2);
- hData.Draw();
- hData.Write();
- if (isGlauber){
- hBestFit.SetLineColor (kRed);
- hBestFit.Draw("same");
-
- std::unique_ptr <TLegend> legData { new TLegend(0.6,0.75,0.75,0.83) };
- legData->AddEntry(&hBestFit ,"Fit", "l");
- legData->AddEntry(&hData ,"Data", "l");
- legData->Draw("same");
- hBestFit.Write();
- }
- }
-
- if (isNBD){
- c1->cd(3);
- hNBD.Draw();
- hNBD.Write();
-
- }
+ inline void DrawHistos (const Fitter& fit, Bool_t isSim, Bool_t isData, Bool_t isGlauber, Bool_t isNBD )
+ {
+ std::unique_ptr <TCanvas> c1 {new TCanvas("c1", "canvas", 1500, 900)};
+
+ c1->Divide(2,2);
+
+ std::unique_ptr <TPad> c1_1 {(TPad*) c1->GetListOfPrimitives()->FindObject("c1_1")};
+ std::unique_ptr <TPad> c1_2 {(TPad*) c1->GetListOfPrimitives()->FindObject("c1_2")};
+ std::unique_ptr <TPad> c1_4 {(TPad*) c1->GetListOfPrimitives()->FindObject("c1_4")};
+
+ c1_1->SetLogy(1);
+ c1_2->SetLogy(1);
+ c1_4->SetLogy(1);
+
+ /*const*/ TH1F hGlaub = fit.GetGlauberFitHisto();
+ /*const*/ TH1F hData = fit.GetDataHisto();
+ /*const*/ TH1F hNBD = fit.GetNBDHisto();
+ /*const*/ TH1F hNcoll = fit.GetNcollHisto();
+ /*const*/ TH1F hNpart = fit.GetNpartHisto();
+ /*const*/ TH1F hBestFit = fit.GetBestFiHisto();
+
+ std::unique_ptr <TFile> fOut{TFile::Open("glauber_qa.root", "recreate")};
+
+ if (isSim){
+ c1->cd(1);
+ hNcoll.SetLineColor(2);
+
+ hNcoll.Draw();
+ hNpart.Draw("same");
+
+ std::unique_ptr <TLegend> legSim { new TLegend(0.6,0.75,0.75,0.83) };
+ legSim->AddEntry(&hNpart ,"Npart", "l");
+ legSim->AddEntry(&hNcoll ,"hNcoll", "l");
+ legSim->Draw("same");
+
+ hNcoll.Write();
+ hNpart.Write();
+ }
+
+ if (isData){
+ c1->cd(2);
+ hData.Draw();
+ hData.Write();
+ if (isGlauber){
+ hBestFit.SetLineColor (kRed);
+ hBestFit.Draw("same");
- if (isGlauber){
- c1->cd(4);
- hBestFit.Draw();
- }
-
- c1->Write();
- c1->SaveAs("glauber.pdf");
- fOut->Close();
+ std::unique_ptr <TLegend> legData { new TLegend(0.6,0.75,0.75,0.83) };
+ legData->AddEntry(&hBestFit ,"Fit", "l");
+ legData->AddEntry(&hData ,"Data", "l");
+ legData->Draw("same");
+ hBestFit.Write();
+ }
}
+ if (isNBD){
+ c1->cd(3);
+ hNBD.Draw();
+ hNBD.Write();
+
+ }
+
+ if (isGlauber){
+ c1->cd(4);
+ hBestFit.Draw();
+ }
+ c1->Write();
+ c1->SaveAs("glauber.pdf");
+ fOut->Close();
+ }
+
+
}
diff --git a/glauber/main.cpp b/glauber/main.cpp
index 7a57337..7f6477f 100644
--- a/glauber/main.cpp
+++ b/glauber/main.cpp
@@ -1,4 +1,5 @@
#include <iostream>
+#include <chrono>
#include "Fitter.h"
#include "FitterHelper.h"
@@ -13,69 +14,77 @@ using namespace Glauber;
int main(int argc, char *argv[])
{
- if (argc < 2)
- {
- std::cout << "Wrong number of parameters! Executable usage:" << std::endl;
- std::cout << " ./glauber f0 k0" << std::endl;
- return -1;
- }
- const Float_t f0 = atof( argv[1]);
- const Int_t k0 = atoi( argv[2] );
- const Int_t k1 = atoi( argv[2] );
-
- // *****************************************
- // Modify this part according to your needs
- // *****************************************
- /// | mode | function for Na |
- /// | Default | Npart + (1-f)*Ncoll |
- /// | PSD | f - Npart |
- /// | Npart | Npart^f |
- /// | Ncoll | Ncoll^f |
- const TString mode = "Default";
-
- const TString glauber_filename = "../input/glauber_auau_sigma_30_100k.root"; // input files
- const TString glauber_treename = "nt_Au_Au";
- const TString in_filename = "../input/test_input.root";
- const TString histoname = "hMreco";
+ if (argc < 2)
+ {
+ std::cout << "Wrong number of parameters! Executable usage:" << std::endl;
+ std::cout << " ./glauber f0 k0" << std::endl;
+ return -1;
+ }
+ const Float_t f0 = atof( argv[1]);
+ const Int_t k0 = atoi( argv[2] );
+ const Int_t k1 = atoi( argv[2] );
+
+ // *****************************************
+ // Modify this part according to your needs
+ // *****************************************
- const Int_t min_bin = 50; // not fitting low multiplicity region due to trigger bias, etc
- const Int_t max_bin = 10000; // very large number to fit the whole histo
- const Int_t nevents = 100000;
+ /// | mode | function for Na |
+ /// | Default | Npart + (1-f)*Ncoll |
+ /// | PSD | f - Npart |
+ /// | Npart | Npart^f |
+ /// | Ncoll | Ncoll^f |
+ const TString mode = "Default";
+
+ const TString glauber_filename = "../input/glauber_auau_sigma_30_100k.root"; // input files
+ const TString glauber_treename = "nt_Au_Au";
+ const TString in_filename = "../input/test_input.root";
+ const TString histoname = "hMreco";
- const TString outdir = ".";
- // *****************************************
- // *****************************************
+ const Int_t min_bin = 50; // not fitting low multiplicity region due to trigger bias, etc
+ const Int_t max_bin = 10000; // very large number to fit the whole histo
+ const Int_t nevents = 100000;
- std::unique_ptr<TFile> glauber_file{ TFile::Open(glauber_filename, "read") };
- std::unique_ptr<TTree> glauber_tree{ (TTree*) glauber_file->Get(glauber_treename) };
+ const TString outdir = ".";
+ // *****************************************
+ // *****************************************
- std::unique_ptr<TFile> f{TFile::Open(in_filename)};
- TH1F *hdata = (TH1F*)f->Get(histoname);
-
- Fitter fitter ( std::move(glauber_tree) );
+ std::unique_ptr<TFile> glauber_file{ TFile::Open(glauber_filename, "read") };
+ std::unique_ptr<TTree> glauber_tree{ (TTree*) glauber_file->Get(glauber_treename) };
- fitter.SetMode(mode);
- fitter.SetInputHisto(*hdata);
- fitter.SetBinSize(1);
- fitter.Init(nevents);
-
- fitter.SetFitMinBin(min_bin);
- fitter.SetFitMaxBin(max_bin);
- fitter.SetOutDirName(outdir);
+ std::unique_ptr<TFile> f{TFile::Open(in_filename)};
+ TH1F *hdata = (TH1F*)f->Get(histoname);
+
+ Fitter fitter ( std::move(glauber_tree) );
- float par[3];
- const float chi2 = fitter.FitGlauber(par, f0, k0, k1, nevents);
+ fitter.SetMode(mode);
+ fitter.SetInputHisto(*hdata);
+ fitter.SetBinSize(1);
+ fitter.Init(nevents);
+
+ fitter.SetFitMinBin(min_bin);
+ fitter.SetFitMaxBin(max_bin);
+ fitter.SetOutDirName(outdir);
- std::cout << "f = " << par[0] << " mu = " << par[1] << " k = " << par[2] << " chi2 = " << chi2 << std::endl;
-
- Glauber::DrawHistos(fitter, true, true, true, true);
+ float par[3];
+
+ auto start = std::chrono::system_clock::now();
+
+ const float chi2 = fitter.FitGlauber(par, f0, k0, k1, nevents);
+
+ auto end = std::chrono::system_clock::now();
+ std::chrono::duration<double> elapsed_seconds = end - start;
+ std::cout << "elapsed time: " << elapsed_seconds.count() << " s\n";
- const float range[2] = {300, 350.};
- std::unique_ptr<TH1F> hB(fitter.GetModelHisto (range, "B", par, 100000));
- hB->SaveAs( "b_test.root" );
-
- std::cout << "END!" << std::endl;
+ std::cout << "f = " << par[0] << " mu = " << par[1] << " k = " << par[2] << " chi2 = " << chi2 << std::endl;
+
+ Glauber::DrawHistos(fitter, true, true, true, true);
- return 0;
+ const float range[2] = {300, 350.};
+ std::unique_ptr<TH1F> hB(fitter.GetModelHisto (range, "B", par, 100000));
+ hB->SaveAs( "b_test.root" );
+
+ std::cout << "END!" << std::endl;
+
+ return 0;
}
--
GitLab