From 8ce8a616cec74b5b95d2f8aa17fef57dfe80ec25 Mon Sep 17 00:00:00 2001 From: Ilya Segal Date: Wed, 18 Sep 2019 17:32:55 +0200 Subject: [PATCH 01/10] Replace Fitter.cpp --- glauber/Fitter.cpp | 660 +++++++++++++++++++++++++++------------------ 1 file changed, 392 insertions(+), 268 deletions(-) diff --git a/glauber/Fitter.cpp b/glauber/Fitter.cpp index 854fce3..e350b4c 100644 --- a/glauber/Fitter.cpp +++ b/glauber/Fitter.cpp @@ -7,292 +7,415 @@ #include "TLegend.h" #include "TMath.h" #include "TRandom.h" +#include ClassImp(Glauber::Fitter) -// ----- Default constructor ------------------------------------------- -Glauber::Fitter::Fitter(std::unique_ptr tree) { - fSimTree = std::move(tree); - std::cout << fSimTree->GetEntries() << std::endl; +void Glauber::Fitter::Init(std::unique_ptr tree) +{ - if (!fSimTree) { - std::cout << "SetSimHistos: *** Error - " << std::endl; - exit(EXIT_FAILURE); - } + fSimTree = std::move(tree); + + if (!fSimTree) + { + std::cout << "SetSimHistos: *** Error - " << std::endl; + exit(EXIT_FAILURE); + } + + std::map ::iterator it0=Glauber_Parameters.begin(); + for (int i=0; iSetBranchAddress(it0->first, &(it0->second)); + it0++; + } - fSimTree->SetBranchAddress("Npart", &fNpart); - fSimTree->SetBranchAddress("Ncoll", &fNcoll); -} + if ( nEvents < 0 || nEvents > 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); + } -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); - } + it0=Glauber_Parameters.begin(); + for (int i=0; i ((it0->first), int (fSimTree->GetMaximum(it0->first))) ); + Glauber_ParametersMin.insert( std::pair ((it0->first), int (fSimTree->GetMinimum(it0->first))) ); + it0++; + } + + it0=Glauber_Parameters.begin(); + std::map ::iterator it1=Glauber_ParametersMin.begin(); + std::map ::iterator it2=Glauber_ParametersMax.begin(); + for (int j=0; jfirst) Glauber_Parameters_Histos.insert( std::pair ((it2->first), new TH1F (Form("%s_Histo", (gGlauberParameters[i].name).Data()), Form("%s;%s;counts", (gGlauberParameters[i].title).Data(), (gGlauberParameters[i].axis_title).Data()), 1.3*((it2->second)-(it1->second)+1)/gGlauberParameters[i].bin_value, 1.3*it1->second, 1.3*it2->second)) ); + it1++; + it2++; + } - const int NpartMax = int(fSimTree->GetMaximum("Npart")); - const int NcollMax = int(fSimTree->GetMaximum("Ncoll")); + it0=Glauber_Parameters.begin(); + it1=Glauber_ParametersMin.begin(); + it2=Glauber_ParametersMax.begin(); + std::map ::iterator it3=Glauber_Parameters_Histos.begin(); + for (int i=0; iGetEntry(i); + it0=Glauber_Parameters.begin(); + it3=Glauber_Parameters_Histos.begin(); + for (int j=0; jsecond)->Fill(it0->second); + it0++; + it3++; + } + } - fNpartHisto = TH1F("fNpartHisto", "Npart", NpartMax / fBinSize, 0, NpartMax); - fNcollHisto = TH1F("fNcollHisto", "Ncoll", NcollMax / fBinSize, 0, NcollMax); + std::cout << "Entries in GlauberTree:" << fSimTree->GetEntries() << std::endl; - 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(); - fNbins = fDataHisto.GetNbinsX(); + while ( fDataHisto.GetBinContent(fNbins-1) == 0) + fNbins--; - while (fDataHisto.GetBinContent(fNbins - 1) == 0) - fNbins--; + fNbins++; - fNbins++; + const float min = fDataHisto.GetXaxis()->GetXmin(); + const float max = fDataHisto.GetXaxis()->GetXmax(); - const auto min = fDataHisto.GetXaxis()->GetXmin(); - const auto max = fDataHisto.GetXaxis()->GetXmax(); + fMaxValue = min + (max - min)*fNbins/fDataHisto.GetNbinsX() ; - fMaxValue = min + (max - min) * fNbins / fDataHisto.GetNbinsX(); + it0=Glauber_Parameters.begin(); + it1=Glauber_ParametersMin.begin(); + it2=Glauber_ParametersMax.begin(); + it3=Glauber_Parameters_Histos.begin(); + for (int j=0; jfirst) Glauber_Parameters_VS_Multiplicity_Histos.insert( std::pair ((it2->first), new TH2F (Form("%s_VS_Multiplicity_Histo", (gGlauberParameters[i].name).Data()), Form("%s VS Multiplicity;Multiplicity;%s", (gGlauberParameters[i].title).Data(), (gGlauberParameters[i].axis_title).Data()), 1.3*fNbins, 0, 1.3*fMaxValue, 1.3*((it2->second)-(it1->second)+1)/gGlauberParameters[i].bin_value, 1.3*it1->second, 1.3*it2->second)) ); + it1++; + it2++; + } - std::cout << "fNbins = " << fNbins << std::endl; - std::cout << "fMaxValue = " << fMaxValue << std::endl; + it0=Glauber_Parameters.begin(); + it1=Glauber_ParametersMin.begin(); + it2=Glauber_ParametersMax.begin(); + it3=Glauber_Parameters_Histos.begin(); + for (int j=0; jfirst) Glauber_Parameters_VS_Multiplicity_BestHistos.insert( std::pair ((it2->first), new TH2F (Form("%s_VS_Multiplicity_BestHisto", (gGlauberParameters[i].name).Data()), Form("%s VS Multiplicity;Multiplicity;%s", (gGlauberParameters[i].title).Data(), (gGlauberParameters[i].axis_title).Data()), 1.3*fNbins, 0, 1.3*fMaxValue, 1.3*((it2->second)-(it1->second)+1)/gGlauberParameters[i].bin_value, 1.3*it1->second, 1.3*it2->second)) ); + it1++; + it2++; + } + std::cout << "fNbins = " << fNbins << std::endl; + std::cout << "fMaxValue = " << fMaxValue << std::endl; } -double Glauber::Fitter::Nancestors(double 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.; +float Glauber::Fitter::Nancestors(float f) const +{ + if (fAncestor_Mode == "Default") return f*Glauber_Parameters.at("Npart") + (1-f)*Glauber_Parameters.at("Ncoll"); + else if (fAncestor_Mode == "PSD") return f-Glauber_Parameters.at("Npart"); + else if (fAncestor_Mode == "Npart") return TMath::Power(Glauber_Parameters.at("Npart"), f); + else if (fAncestor_Mode == "Ncoll") return TMath::Power(Glauber_Parameters.at("Ncoll"), f); + + return -1.; } -double Glauber::Fitter::NancestorsMax(double 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.; +float Glauber::Fitter::NancestorsMax(float f) const +{ + const int NpartMax = Glauber_Parameters_Histos.at("Npart")->GetXaxis()->GetXmax() ; // some magic + const int NcollMax = Glauber_Parameters_Histos.at("Ncoll")->GetXaxis()->GetXmax() ; //TODO + + if (fAncestor_Mode == "Default") return f*NpartMax + (1-f)*NcollMax; + else if (fAncestor_Mode == "PSD") return f; + else if (fAncestor_Mode == "Npart") return TMath::Power(NpartMax, f); + else if (fAncestor_Mode == "Ncoll") return TMath::Power(NcollMax, f); + + return -1.; } /* * take Glauber MC data from fSimTree * Populate fGlauberFitHisto with NBD x Na */ -void Glauber::Fitter::SetGlauberFitHisto(double f, double mu, double 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 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)); +void Glauber::Fitter::SetGlauberFitHisto (float f, float mu, float k, Bool_t Norm2Data) +{ + fGlauberFitHisto = TH1F("glaub", "", fNbins*1.3, 0, 1.3*fMaxValue); + + fGlauberFitHisto.SetName("glaub_fit_histo"); + + SetNBDhist(mu, k); + + std::unique_ptr htemp {(TH1F*)fNbdHisto.Clone("htemp")}; // WTF??? Not working without pointer + for (int i=0; iGetEntry(i); + const int Na = int(Nancestors(f)); + + float nHits {0.}; + for (int j=0; jGetRandom()); + fGlauberFitHisto.Fill(nHits); + + std::map ::iterator it1=Glauber_Parameters.begin(); + std::map ::iterator it2=Glauber_Parameters_VS_Multiplicity_Histos.begin(); + for (int q=0; qsecond)->Fill(nHits, it1->second); + it1++; + it2++; + } - double nHits{0.}; - for (int j = 0; j < Na; j++) { - nHits += int(htemp->GetRandom()); } - fGlauberFitHisto.Fill(nHits); - } - - if (Norm2Data) - NormalizeGlauberFit(); + 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 double ScaleFactor = (double) fDataHistoInt / fGlauberFitHistoInt; +void Glauber::Fitter::NormalizeGlauberFit () +{ + + int fGlauberFitHistoInt {0}; + int fDataHistoInt {0}; + + const int lowchibin = fFitMinBin; + const int highchibin = fFitMaxBin file {TFile::Open(filename, "recreate")}; - // std::unique_ptr 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); - - double 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 double mu_min = 0.7 * mu; - const double 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; +float Glauber::Fitter::FitGlauber (float *par, Float_t f0, Float_t f1, Int_t k0, Int_t k1) +{ + float f_fit{-1}; + float mu_fit{-1}; + float k_fit{-1}; + float Chi2Min {1e10}; + float Chi2Min_error {0}; + + const TString filename = Form ( "%s/fit%s.root", fOutDirName.Data(), fOutFileIDName.Data() ); + +// std::unique_ptr file {TFile::Open(filename, "recreate")}; +// std::unique_ptr tree {new TTree("test_tree", "tree" )}; + + TFile* file {TFile::Open(filename, "recreate")}; + TTree* tree {new TTree("test_tree", "tree" )}; + + float f, mu, k, chi2, chi2_error, sigma; + + 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("chi2_error", &chi2_error, "chi2_error/F"); + tree->Branch("sigma",&sigma,"sigma/F"); + + int n=1; + for (float i=f0; i<=f1; i=i+0.1000000) + { + f=i; + 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, &chi2_error, mu_min, mu_max, f, k, n); + n=n+fnIter; + sigma = ( mu/k + 1 ) * mu; + + tree->Fill(); + + if (chi2 < Chi2Min) + { + f_fit = f; + mu_fit = mu; + k_fit = k; + Chi2Min = chi2; + Chi2Min_error = chi2_error; + fBestFitHisto = fGlauberFitHisto; + + std::map ::iterator it1=Glauber_Parameters_VS_Multiplicity_BestHistos.begin(); + std::map ::iterator it2=Glauber_Parameters_VS_Multiplicity_Histos.begin(); + for (int q=0; qsecond)=(it2->second); + it1++; + it2++; + } + } + + } } - } - tree->Write(); - file->Write(); - file->Close(); + tree->Write(); + file->Write(); + file->Close(); - par[0] = f_fit; - par[1] = mu_fit; - par[2] = k_fit; - - return Chi2Min; + par[0] = f_fit; + par[1] = mu_fit; + par[2] = k_fit; + par[3] = Chi2Min_error; + + return Chi2Min; } /** * Compare fGlauberFitHisto with fDataHisto * @return chi2 value */ -double Glauber::Fitter::GetChi2() const { - double 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 double error2 = TMath::Power(fDataHisto.GetBinError(i), 2) + TMath::Power(fGlauberFitHisto.GetBinError(i), 2); - const double diff2 = TMath::Power((fGlauberFitHisto.GetBinContent(i) - fDataHisto.GetBinContent(i)), 2) / error2; - - chi2 += diff2; - } +float Glauber::Fitter::GetChi2 () const +{ + float chi2 {0.0}; + + const int lowchibin = fFitMinBin; + const int highchibin = fFitMaxBin 1e-10) fNbdHisto.SetBinContent(i + 1, val); - // std::cout << "val " << val << std::endl; - } +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("nbd"); + + for (int i=0; i1e-10) fNbdHisto.SetBinContent(i+1, val); +// std::cout << "val " << val << std::endl; + } } /** @@ -302,65 +425,66 @@ void Glauber::Fitter::SetNBDhist(double mu, double k) { * @param k argument * @return NBD for a given parameters */ -double Glauber::Fitter::NBD(double n, double mu, double 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 = 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; +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; } /** * Creates histo with a given model parameter distribution * @param range observable range * @param name name of the MC-Glauber model parameter * @param par array with fit parameters - * @param Nevents * @return pointer to the histogram */ -std::unique_ptr Glauber::Fitter::GetModelHisto(const double range[2], - const TString &name, - const double par[3], - int nEvents) { - const double f = par[0]; - const double mu = par[1]; - const double k = par[2]; - - double modelpar{-999.}; - fSimTree->SetBranchAddress(name, &modelpar); - - SetNBDhist(mu, k); - std::unique_ptr 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)); - double nHits{0.}; - for (int j = 0; j < Na; ++j) { - nHits += (int) fNbdHisto.GetRandom(); - } - - if (nHits > range[0] && nHits < range[1]) { - hModel->Fill(modelpar); +std::unique_ptr Glauber::Fitter::GetModelHisto (const float range[2], TString name, const float par[3]) +{ + 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); + +// TRandom random; +// random.SetSeed(mu*k); + + std::unique_ptr hModel(new TH1F ("hModel", "name", 100, fSimTree->GetMinimum(name), fSimTree->GetMaximum(name)) ); + + for (int i=0; iGetEntry(i); + const int Na = int(Nancestors(f)); + float nHits{0.}; + for (int j=0; j range[0] && nHits < range[1] ){ + hModel->Fill(modelpar); + } + } - } - - return hModel; - + + return std::move(hModel); + } - - - From c4f34d00486ec85c3dd715adf6fd63d85c80eeb6 Mon Sep 17 00:00:00 2001 From: Ilya Segal Date: Wed, 18 Sep 2019 17:33:31 +0200 Subject: [PATCH 02/10] Replace Fitter.h --- glauber/Fitter.h | 321 +++++++++++++++++++++++++++++++++-------------- 1 file changed, 228 insertions(+), 93 deletions(-) diff --git a/glauber/Fitter.h b/glauber/Fitter.h index 07e02d2..ff9c859 100644 --- a/glauber/Fitter.h +++ b/glauber/Fitter.h @@ -1,109 +1,244 @@ /** @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 #include +#include #include "TString.h" #include "TNamed.h" #include "TH1F.h" +#include "TH2F.h" #include "TTree.h" // #include "TMinuit.h" -namespace Glauber { -class Fitter { - - public: - - /** Default constructor **/ - Fitter() = default;; - explicit Fitter(std::unique_ptr tree); - /** Destructor **/ - virtual ~Fitter() = default;; - - void Init(int nEntries); - void SetGlauberFitHisto(double f, double mu, double 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) {}; - - double FitGlauber(double *par, double f0, Int_t k0, Int_t k1, Int_t nEvents); - void FindMuGoldenSection(double *mu, - double *chi2, - double mu_min, - double mu_max, - double f, - double k, - Int_t nEvents = 10000, - Int_t nIter = 5); - - double GetChi2() const; - - double NBD(double n, double mu, double k) const; - void SetNBDhist(double mu, double k); - double Nancestors(double f) const; - double NancestorsMax(double f) const; - - std::unique_ptr GetModelHisto(const double range[2], const TString &name, const double 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(const 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 fSimTree{nullptr}; - - double fNpart{-1.}; - double fNcoll{-1.}; - - double 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); - -}; +namespace Glauber +{ + enum eGlauberParameters { + kNpart = 0, kNcoll, kNhard, kB, kBNN, kNcollpp, kNcollpn, kNcollnn, kVarX, kVarY, kVarXY, kNpartA, kNpartB, kNpart0, kAreaW, kPsi1, kEcc1, kPsi2, kEcc2, kPsi3, kEcc3, kPsi4, kEcc4, kPsi5, kEcc5, kAreaO, kAreaA, kX0, kY0, kPhi0, kLength, kMeanX, kMeanY, kMeanX2, kMeanY2, kMeanXY, kMeanXSystem, kMeanYSystem, kMeanXA, kMeanYA, kMeanXB, kMeanYB, kPhiA, kThetaA, kPhiB, kThetaB, kGP + }; + + const struct TGlauberParameters { + Int_t id; + TString name; + TString title; + TString axis_title; + Float_t bin_value; + } gGlauberParameters[kGP] = { + { .id = kNpart, .name = "Npart", .title = "Number of participating nucleons", .axis_title="N_{part}", .bin_value = 1 }, + { .id = kNcoll, .name = "Ncoll", .title = "Number of binary collisions", .axis_title="N_{coll}", .bin_value = 1 }, + { .id = kNhard, .name = "Nhard", .title = "Number of hard collisions", .axis_title="N_{hard}", .bin_value = 1 }, + { .id = kB, .name = "B", .title = "Impact parameter", .axis_title="#it{B}, fm", .bin_value = 0.01 }, + { .id = kBNN, .name = "BNN", .title = "Mean nucleon–nucleon impact parameter", .axis_title="#it{B_{NN}}, fm", .bin_value = 0.01 }, + { .id = kNcollpp, .name = "Ncollpp", .title = "Number of pp collisions", .axis_title="N^{pp}_{coll}", .bin_value = 1 }, + { .id = kNcollpn, .name = "Ncollpn", .title = "Number of pn collisions", .axis_title="N^{pn}_{coll}", .bin_value = 1 }, + { .id = kNcollnn, .name = "Ncollnn", .title = "Number of nn collisions", .axis_title="N^{nn}_{coll}", .bin_value = 1 }, + { .id = kVarX, .name = "VarX", .title = "Variance of x for wounded nucleons", .axis_title="#sigma^{2}_{x}, fm^{2}", .bin_value = 0.01 }, + { .id = kVarY, .name = "VarY", .title = "Variance of y for wounded nucleons", .axis_title="#sigma^{2}_{x}, fm^{2}", .bin_value = 0.01 }, + { .id = kVarXY, .name = "VarXY", .title = "Covariance of x and y for wounded nucleons", .axis_title="sigma_{xy}, fm^{2}", .bin_value = 0.01 }, + { .id = kNpartA, .name = "NpartA", .title = "Number of participating nucleons from nucleus A", .axis_title="N^{A}_{part}", .bin_value = 1 }, + { .id = kNpartB, .name = "NpartB", .title = "Number of participating nucleons from nucleus B", .axis_title="N^{B}_{part}", .bin_value = 1 }, + { .id = kNpart0, .name = "Npart0", .title = "Number of singly-wounded participating nucleons", .axis_title="N^{0}_{part}", .bin_value = 1 }, + { .id = kAreaW, .name = "AreaW", .title = "Area defined by width of participants", .axis_title="AreaW, fm^{2}", .bin_value = 0.01 }, + { .id = kPsi1, .name = "Psi1", .title = "Event plane angle of 1st harmonic", .axis_title="#psi_{1}, rad", .bin_value = 0.01 }, + { .id = kEcc1, .name = "Ecc1", .title = "Participant eccentricity for 1st harmonic", .axis_title="#epsilon_{1}", .bin_value = 0.01 }, + { .id = kPsi2, .name = "Psi2", .title = "Event plane angle of 2nd harmonic", .axis_title="#psi_{2}, rad", .bin_value = 0.01 }, + { .id = kEcc2, .name = "Ecc2", .title = "Participant eccentricity for 2nd harmonic", .axis_title="#epsilon_{2}", .bin_value = 0.01 }, + { .id = kPsi3, .name = "Psi3", .title = "Event plane angle of 3rd harmonic", .axis_title="#psi_{3}, rad", .bin_value = 0.01 }, + { .id = kEcc3, .name = "Ecc3", .title = "Participant eccentricity for 3rd harmonic", .axis_title="#epsilon_{3}", .bin_value = 0.01 }, + { .id = kPsi4, .name = "Psi4", .title = "Event plane angle of 4th harmonic", .axis_title="#psi_{4}, rad", .bin_value = 0.01 }, + { .id = kEcc4, .name = "Ecc4", .title = "Participant eccentricity for 4th harmonic", .axis_title="#epsilon_{4}", .bin_value = 0.01 }, + { .id = kPsi5, .name = "Psi5", .title = "Event plane angle of 5th harmonic", .axis_title="#psi_{5}, rad", .bin_value = 0.01 }, + { .id = kEcc5, .name = "Ecc5", .title = "Participant eccentricity for 5th harmonic", .axis_title="#epsilon_{5}", .bin_value = 0.01 }, + { .id = kAreaO, .name = "AreaO", .title = "Area by ”or” of participants in grid", .axis_title="Area OR", .bin_value = 0.01 }, + { .id = kAreaA, .name = "AreaA", .title = "Area by ”and” of participants in grid", .axis_title="Area AND", .bin_value = 0.01 }, + { .id = kX0, .name = "X0", .title = "Production point in x", .axis_title="x, fm", .bin_value = 0.01 }, + { .id = kY0, .name = "Y0", .title = "Production point in y", .axis_title="y, fm", .bin_value = 0.01 }, + { .id = kPhi0, .name = "Phi0", .title = "Direction in #varphi", .axis_title="#varphi, rad", .bin_value = 0.01 }, + { .id = kLength, .name = "Length", .title = "Length in #varphi", .axis_title="p_{T} (GeV/#it{c})", .bin_value = 0.01 }, + { .id = kMeanX, .name = "MeanX", .title = "Mean of x for wounded nucleons", .axis_title=", fm", .bin_value = 0.01 }, + { .id = kMeanY, .name = "MeanY", .title = "Mean of y for wounded nucleons", .axis_title=", fm", .bin_value = 0.01 }, + { .id = kMeanX2, .name = "MeanX2", .title = "Mean of x^{2} for wounded nucleons", .axis_title=", fm^2", .bin_value = 0.01 }, + { .id = kMeanXY, .name = "MeanY2", .title = "Mean of y^{2} for wounded nucleons", .axis_title=", fm^2", .bin_value = 0.01 }, + { .id = kMeanY2, .name = "MeanXY", .title = "Mean of xy for wounded nucleons", .axis_title=", fm^2", .bin_value = 0.01 }, + { .id = kMeanXSystem, .name = "MeanXSystem", .title = "Mean of x for all nucleons", .axis_title=", fm", .bin_value = 0.01 }, + { .id = kMeanYSystem, .name = "MeanYSystem", .title = "Mean of y for all nucleons", .axis_title=", fm", .bin_value = 0.01 }, + { .id = kMeanXA, .name = "MeanXA", .title = "Mean of x for nucleons in nucleus A", .axis_title=", fm", .bin_value = 0.01 }, + { .id = kMeanYA, .name = "MeanYA", .title = "Mean of y for nucleons in nucleus A", .axis_title=", fm", .bin_value = 0.01 }, + { .id = kMeanXB, .name = "MeanXB", .title = "Mean of x for nucleons in nucleus B", .axis_title=", fm", .bin_value = 0.01 }, + { .id = kMeanYB, .name = "MeanYB", .title = "Mean of y for nucleons in nucleus B", .axis_title=", fm", .bin_value = 0.01 }, + { .id = kPhiA, .name = "PhiA", .title = "Azimuthal angle for nucleus A", .axis_title="#varphi^{A}, rad", .bin_value = 0.01 }, + { .id = kThetaA, .name = "ThetaA", .title = "Polar angle for nucleus A", .axis_title="#theta^{A}, rad", .bin_value = 0.01 }, + { .id = kPhiB, .name = "PhiB", .title = "Azimuthal angle for nucleus B", .axis_title="#varphi^{B}, rad", .bin_value = 0.01 }, + { .id = kThetaB, .name = "ThetaB", .title = "Polar angle for nucleus B", .axis_title="#theta^{A}, rad", .bin_value = 0.01 } + }; + + class Fitter + { + + public: + + /** Default constructor **/ + Fitter() {}; + /** Destructor **/ + virtual ~Fitter(){}; + + void Init(std::unique_ptr tree); // Initialize all parameters of class Fitter + void SetGlauberFitHisto (Float_t f, Float_t mu, Float_t k, Bool_t Norm2Data = true); // Setting of model histogtam + void NormalizeGlauberFit (); // Function for normilize model histogram with data histogram + void DrawHistos (Bool_t isSim = true, Bool_t isData = true, Bool_t isGlauber = false, Bool_t isNBD = false); // Function for drawing output histograms and creating output files + + float FitGlauber (float *par, Float_t f0, Float_t f1, Int_t k0, Int_t k1); // Main function for fitting + void FindMuGoldenSection (Float_t *mu, Float_t *chi2, float*chi2_error, Float_t mu_min, Float_t mu_max, Float_t f, Float_t k, int n=0); // Function that find mu by iteration it to the minimum of chi2 + + Float_t GetChi2 (void) const; // Function for chi2 counting + Float_t GetChi2Error (void) const; // Function for chi2 error counting + + Float_t NBD(Float_t n, Float_t mu, Float_t k) const; // Function for NBD counting with parameters mu and k + void SetNBDhist(Float_t mu, Float_t k); // Setting of NBD distribution with parameters mu and k + + float Nancestors(float f) const; // Function for counting amount of ancestors in one event + float NancestorsMax(float f) const; // Function for counting maximum amount of ancestors in one event + + std::unique_ptr GetModelHisto (const Float_t range[2], TString name, const Float_t par[4]); // Returns model histogram with known parameters + +// +// Setters +// + void SetInputHisto (const TH1F &h) { fDataHisto = h; } // Setting of data histogram + void SetFitMinBin (Int_t min) { fFitMinBin = min; } // Setting of not fitting low multiplicity region due to trigger bias, etc + void SetFitMaxBin (Int_t max) { fFitMaxBin = max; } // Setting of very large number to fit the whole histo + void SetnIter (Int_t nIter) { fnIter = nIter; } // Setting of number of iterations of finding better mu (with same f and k) + void SetOutDirName (const TString name) { fOutDirName = name; } // Setting of the name of directory for output files + void SetOutFileIDName (const TString name) { fOutFileIDName = name; } // Setting of ID of output files (ending of all output files of one job will have same ID) + void SetAncestor_Mode (const TString Ancestor_Mode) { fAncestor_Mode = Ancestor_Mode; } // Setting of the mode for counting ancestors + void SetFitFunction_Mode (const TString FitFunction_Mode) { fFitFunction_Mode = FitFunction_Mode; } // Setting of the mode for fitting function (COMING SOON) + void SetGlauber_filename (const TString Glauber_filename) { fGlauber_filename = Glauber_filename; } // Setting of input file with Glauber Tree + void SetGlauber_treename (const TString Glauber_treename) { fGlauber_treename = Glauber_treename; } // Setting of the name of tree with GlauberMC model + void SetDataHisto_filename (const TString DataHisto_filename) { fDataHisto_filename = DataHisto_filename; } // Setting of input file with data histo + void SetDataHisto_name (const TString DataHisto_name) { fDataHisto_name = DataHisto_name; } // Setting of the name of data histo + void SetMassNumber (Float_t A) { fA = A; } // Setting of mass number of projectile nucleus + void SetNEvents (Int_t Events) { nEvents = Events; } // Setting of number of events that you need from GlauberMC + +// +// Getters +// + TString GetOutFileIDName () const { return fOutFileIDName; } // Returns ID of output files (its ending) + TString GetGlauber_filename () { return fGlauber_filename; } // Returns input file with Glauber Tree + TString GetGlauber_treename () { return fGlauber_treename; } // Returns the name of tree with + TString GetDataHisto_filename () { return fDataHisto_filename; } // Returns input file with data histo + TString GetDataHisto_name () { return fDataHisto_name; } // Returns the name of data histo + Int_t GetFitMinBin () const { return fFitMinBin; } // Returns not fitting low multiplicity region due to trigger bias, etc + Int_t GetFitMaxBin () const { return fFitMaxBin; } // Returns very large number to fit the whole histo + + TH1F GetGlauberFitHisto () const { return fGlauberFitHisto; } // Returns model histogtam + TH1F GetBestFitHisto () const { return fBestFitHisto; } // Returns model histogtam with minimal chi2 + TH1F GetDataHisto () const { return fDataHisto; } // Returns data histogram + TH1F GetNBDHisto () const { return fNbdHisto; } // Returns histogram of NBD distribution + + std::map GetMapOfGlauber_Parameters_Histos () const { return Glauber_Parameters_Histos; } // Returns map with historams of parameters of GlauberMC input file + std::map GetMapOfGlauber_Parameters_VS_Multiplicity_Histos () const { return Glauber_Parameters_VS_Multiplicity_Histos; } // Returns map with historams of parameters of GlauberMC input file versus modeled multiplicity + std::map GetMapOfGlauber_Parameters_VS_Multiplicity_BestHistos () const { return Glauber_Parameters_VS_Multiplicity_BestHistos; } // Returns map with historams of parameters of GlauberMC input file versus modeled multiplicity with minimal #chi^{2} + + private: + + /** Data members **/ + + TH1F fDataHisto; + TH1F fNbdHisto; + TH1F fGlauberFitHisto; + TH1F fBestFitHisto; + + std::map Glauber_Parameters_Histos; + std::map Glauber_Parameters_VS_Multiplicity_Histos; + std::map Glauber_Parameters_VS_Multiplicity_BestHistos; + + + /* MC data */ + std::unique_ptr fSimTree{nullptr}; + + Float_t fA{-1.}; //mass number + std::map Glauber_Parameters = { + {"Npart", -1.}, // Number of participating nucleons + {"Ncoll", -1.}, // Number of binary collisions + {"Nhard", -1.}, // Number of hard collisions (with hard particle production) + {"B", -1.}, // Generated impact parameter, fm + {"BNN", -1.}, // Mean nucleon–nucleon impact parameter, fm + {"Ncollpp", -1.}, // Number of pp collisions + {"Ncollpn", -1.}, // Number of pn collisions + {"Ncollnn", -1.}, // Number of nn collisions + {"VarX", -1.}, // Variance of x for wounded nucleons, σ^2_x, fm^2 + {"VarY", -1.}, // Variance of y for wounded nucleons, σ^2_y, fm^2 + {"VarXY", -1.}, // Covariance of x and y for wounded nucleons, σ_xy ≡ − , fm^2 + {"NpartA", -1.}, // Number of participating nucleons from nucleus A + {"NpartB", -1.}, // Number of participating nucleons from nucleus B + {"Npart0", -1.}, // Number of singly-wounded participating nucleons + {"AreaW", -1.}, // Area defined by width of participants, fm^2 + {"Psi1", -1.}, // Event plane angle of 1st harmonic, rad + {"Ecc1", -1.}, // Participant eccentricity for 1st harmonic + {"Psi2", -1.}, // Event plane angle of 2nd harmonic, rad + {"Ecc2", -1.}, // Participant eccentricity for 2nd harmonic + {"Psi3", -1.}, // Event plane angle of 3rd harmonic, rad + {"Ecc3", -1.}, // Participant eccentricity for 3rd harmonic + {"Psi4", -1.}, // Event plane angle of 4th harmonic, rad + {"Ecc4", -1.}, // Participant eccentricity for 4th harmonic + {"Psi5", -1.}, // Event plane angle of 5th harmonic, rad + {"Ecc5", -1.}, // Participant eccentricity for Nth harmonic + {"AreaO", -1.}, // Area by ”or” of participants in grid + {"AreaA", -1.}, // Area by ”and” of participants in grid + {"X0", -1.}, // Production point in x, fm + {"Y0", -1.}, // Production point in y, fm + {"Phi0", -1.}, // Direction in phi + {"Length", -1.}, // Length in phi0 + {"MeanX", -1.}, // Mean of x for wounded nucleons, , fm + {"MeanY", -1.}, // Mean of y for wounded nucleons, , fm + {"MeanX2", -1.}, // Mean of x^2 for wounded nucleons, , fm^2 + {"MeanY2", -1.}, // Mean of y^2 for wounded nucleons, , fm^2 + {"MeanXY", -1.}, // Mean of xy for wounded nucleons, , fm^2 + {"MeanXSystem", -1.}, // Mean of x for all nucleons, fm^2 + {"MeanYSystem", -1.}, // Mean of y for all nucleons, fm^2 + {"MeanXA", -1.}, // Mean of x for nucleons in nucleus A, fm^2 + {"MeanYA", -1.}, // Mean of y for nucleons in nucleus A, fm^2 + {"MeanXB", -1.}, // Mean of x for nucleons in nucleus B, fm^2 + {"MeanYB", -1.}, // Mean of y for nucleons in nucleus B, fm^2 + {"PhiA", -1.}, // Azimuthal angle for nucleus A, rad + {"ThetaA", -1.}, // Polar angle for nucleus A, rad + {"PhiB", -1.}, // Azimuthal angle for nucleus B, rad + {"ThetaB", -1.}, // Polar angle for nucleus B, rad + }; + + + std::map Glauber_ParametersMax; + std::map Glauber_ParametersMin; + + Float_t fMaxValue{-1.}; + + Int_t fNbins{-1 }; + + Int_t fFitMinBin{-1}; + Int_t fFitMaxBin{-1}; + Int_t nEvents; + + Int_t fnIter; + + TString fAncestor_Mode{"Default"}; + TString fFitFunction_Mode{"NBD"}; + TString fGlauber_filename{""}; + TString fGlauber_treename{""}; + TString fDataHisto_filename{""}; + TString fDataHisto_name{""}; + + TString fOutDirName{""}; + TString fOutFileIDName {""}; + ClassDef(Fitter, 2); + + }; } #endif From 986fa650d22864776a1c75b5f70b37d060323e82 Mon Sep 17 00:00:00 2001 From: Ilya Segal Date: Wed, 18 Sep 2019 17:33:54 +0200 Subject: [PATCH 03/10] Replace FitterHelper.h --- glauber/FitterHelper.h | 174 ++++++++++++++++++++++++----------------- 1 file changed, 101 insertions(+), 73 deletions(-) diff --git a/glauber/FitterHelper.h b/glauber/FitterHelper.h index 79addf2..ea8e5d3 100644 --- a/glauber/FitterHelper.h +++ b/glauber/FitterHelper.h @@ -1,92 +1,120 @@ /** @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 +#include #include "TCanvas.h" #include "TH1.h" +#include "TH2.h" #include "TPad.h" #include "TLegend.h" #include "TFile.h" -#include "glauber/Fitter.h" - -namespace Glauber { -inline void DrawHistos(const Fitter &fit, Bool_t isSim, Bool_t isData, Bool_t isGlauber, Bool_t isNBD) { - std::unique_ptr c1{new TCanvas("c1", "canvas", 1500, 900)}; - - c1->Divide(2, 2); - - std::unique_ptr c1_1{(TPad *) c1->GetListOfPrimitives()->FindObject("c1_1")}; - std::unique_ptr c1_2{(TPad *) c1->GetListOfPrimitives()->FindObject("c1_2")}; - std::unique_ptr 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 fOut{TFile::Open("glauber_qa.root", "recreate")}; - - if (isSim) { - c1->cd(1); - hNcoll.SetLineColor(2); - - hNcoll.Draw(); - hNpart.Draw("same"); - - std::unique_ptr 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 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(); +#include "Fitter.h" + + +namespace Glauber +{ + inline void DrawHistos (const Fitter& fit, float par[4], float chi2, Bool_t isSim, Bool_t isData, Bool_t isGlauber, Bool_t isNBD) + { + std::unique_ptr c1 {new TCanvas("c1", "canvas", 1500, 900)}; + + c1->Divide(2,2); + + std::unique_ptr c1_1 {(TPad*) c1->GetListOfPrimitives()->FindObject("c1_1")}; + std::unique_ptr c1_2 {(TPad*) c1->GetListOfPrimitives()->FindObject("c1_2")}; + std::unique_ptr 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.GetMapOfGlauber_Parameters_Histos()).at("Ncoll")); + /*const*/ TH1F hNpart = *((fit.GetMapOfGlauber_Parameters_Histos()).at("Npart")); + /*const*/ TH1F hBestFit = fit.GetBestFitHisto(); + + TH2F* BestHistos[kGP]; + for (int i=0; i fOut{TFile::Open(Form("glauber_qa%s.root", fit.GetOutFileIDName().Data()), "recreate")}; + + if (isSim){ + c1->cd(1); + hNcoll.SetLineColor(2); + + hNcoll.Draw(); + hNpart.Draw("same"); + + std::unique_ptr legSim { new TLegend(0.6,0.75,0.75,0.83) }; + legSim->AddEntry(&hNpart ,"Npart", "l"); + legSim->AddEntry(&hNcoll ,"Ncoll", "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 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(); + for (int i=0; iWrite(); + } + } + + if (isNBD){ + c1->cd(3); + hNBD.Draw(); + hNBD.Write(); + + } + + if (isGlauber){ + c1->cd(4); + hBestFit.Draw(); + } + + TTree *BestResult=new TTree("BestResult", "BestResult"); + Float_t mu, f, k, chi2_error; + BestResult -> Branch("mu", &mu); + BestResult -> Branch("f", &f); + BestResult -> Branch("k", &k); + BestResult -> Branch("chi2", &chi2); + BestResult -> Branch("chi2_error", &chi2_error); + + mu=par[1]; f=par[0]; k=par[2]; chi2_error=par[3]; + + BestResult -> Fill(); + BestResult -> Write(); + + c1->Write(); + c1->SaveAs(Form("glauber%s.pdf", fit.GetOutFileIDName().Data())); + fOut->Close(); } - } + + +} - if (isNBD) { - c1->cd(3); - hNBD.Draw(); - hNBD.Write(); - } - if (isGlauber) { - c1->cd(4); - hBestFit.Draw(); - } - c1->Write(); - c1->SaveAs("glauber.pdf"); - fOut->Close(); -} -} #endif From ff30de9884762a2cbf7a6cfaa235ce8fe4dbf5fa Mon Sep 17 00:00:00 2001 From: Ilya Segal Date: Wed, 18 Sep 2019 17:34:11 +0200 Subject: [PATCH 04/10] Replace main.cpp --- glauber/main.cpp | 149 ++++++++++++++++++++++++----------------------- 1 file changed, 77 insertions(+), 72 deletions(-) diff --git a/glauber/main.cpp b/glauber/main.cpp index 0229eb2..2ead99e 100644 --- a/glauber/main.cpp +++ b/glauber/main.cpp @@ -1,85 +1,90 @@ #include -#include #include "Fitter.h" #include "FitterHelper.h" #include "TH1.h" +#include "TH2.h" #include "TFile.h" #include "TLegend.h" #include "TH2.h" -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 double_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"; - - 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; - - const TString outdir = "."; - // ***************************************** - // ***************************************** - - std::unique_ptr glauber_file{TFile::Open(glauber_filename, "read")}; - std::unique_ptr glauber_tree{(TTree *) glauber_file->Get(glauber_treename)}; - - std::unique_ptr f{TFile::Open(in_filename)}; - TH1F *hdata = (TH1F *) f->Get(histoname); - Fitter fitter(std::move(glauber_tree)); +#include - fitter.SetMode(mode); - fitter.SetInputHisto(*hdata); - fitter.SetBinSize(1); - fitter.Init(nevents); - - fitter.SetFitMinBin(min_bin); - fitter.SetFitMaxBin(max_bin); - fitter.SetOutDirName(outdir); - - double par[3]; - - auto start = std::chrono::system_clock::now(); - - const double chi2 = fitter.FitGlauber(par, f0, k0, k1, nevents); - - auto end = std::chrono::system_clock::now(); - std::chrono::duration elapsed_seconds = end - start; - std::cout << "elapsed time: " << elapsed_seconds.count() << " s\n"; - - std::cout << "f = " << par[0] << " mu = " << par[1] << " k = " << par[2] << " chi2 = " << chi2 << std::endl; - - Glauber::DrawHistos(fitter, true, true, true, true); - - const double range[2] = {300, 350.}; - std::unique_ptr hB(fitter.GetModelHisto(range, "B", par, 100000)); - hB->SaveAs("b_test.root"); - - std::cout << "END!" << std::endl; +using namespace Glauber; - return 0; +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]); + std::cout << "f=" << f0 << std::endl; + const Float_t f1 = atof( argv[2]); + std::cout << "f=" << f1 << std::endl; + const Int_t k0 = atoi( argv[3] ); + std::cout << "k0=" << k0 << std::endl; + const Int_t k1 = atoi( argv[4] ); + std::cout << "k1=" << k1 << std::endl; + + Fitter fitter; + + // ***************************************** + // Modify this part according to your needs + // ***************************************** + + /// |Ancestor_Mode| function for Na | + /// | Default | f*Npart + (1-f)*Ncoll | + /// | PSD | f-Npart | + /// | Npart | Npart^f | + /// | Ncoll | Ncoll^f | + + + fitter.SetFitMinBin (60); // Set not fitting low multiplicity region due to trigger bias, etc + fitter.SetFitMaxBin (200); // Set very large number to fit the whole histo + fitter.SetGlauber_filename ("/home/vad/NA61/Centrality/glau_PbPb_13GeV.root"); // Set input file with Glauber Tree + fitter.SetGlauber_treename ("nt_Pb_Pb"); // Set name of tree with GlauberMC model + fitter.SetDataHisto_filename ("/home/vad/NA61/QA/qa_new.root"); // Set input file with data histo + fitter.SetDataHisto_name ("reco_info/hMreco"); // Set name of data histo + fitter.SetOutDirName ("."); // Set name of directory for output files + fitter.SetOutFileIDName ("_M13new_5"); // Set ID of output files (ending of all output files of one job will have same ID) + fitter.SetAncestor_Mode ("Default"); // Set mode for counting ancestors + fitter.SetFitFunction_Mode ("NBD"); // Set mode for fitting function (COMING SOON) + fitter.SetnIter (10); // Set number of iterations of finding better mu (with same f and k) + fitter.SetMassNumber (207); // Set mass number of projectile nucleus + fitter.SetNEvents (999999); // Set number of events that you need from GlauberMC (or you can uncomment the same setter a few lines down with default number of events: 10 times the number of events in the data histo in the fitting range) + + std::cout << "min_bin=" << fitter.GetFitMinBin() << " max_bin=" << fitter.GetFitMaxBin() << std::endl; + + // ***************************************** + // ***************************************** + + std::unique_ptr glauber_file{ TFile::Open(fitter.GetGlauber_filename(), "read") }; + std::unique_ptr glauber_tree{ (TTree*) glauber_file->Get(fitter.GetGlauber_treename()) }; + std::unique_ptr f{TFile::Open(fitter.GetDataHisto_filename())}; + TH1F *hdata = (TH1F*)f->Get(fitter.GetDataHisto_name()); + + fitter.SetNEvents(10*(int(hdata->Integral(fitter.GetFitMinBin(),fitter.GetFitMaxBin())))); + + fitter.SetInputHisto(*hdata); + + fitter.Init(std::move(glauber_tree)); + + float par[4]; + + float chi2=1e10; + chi2 = fitter.FitGlauber(par, f0, f1, k0, k1); + + std::cout << "f = " << par[0] << " mu = " << par[1] << " k = " << par[2] << " chi2 = " << chi2 << " chi2_error = " << par[3] << std::endl; + +// Glauber::DrawHistos(fitter, true, true, true, true); + + DrawHistos(fitter, par, chi2, true, true, true, true); + + std::cout << "END!" << std::endl; + + return 0; } From 7f2f84f093913e9533288dd1379c8c5fee595e09 Mon Sep 17 00:00:00 2001 From: Ilya Segal Date: Wed, 18 Sep 2019 18:06:43 +0200 Subject: [PATCH 05/10] Upload New File --- glauber/macro/Chi2.C | 108 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 108 insertions(+) create mode 100644 glauber/macro/Chi2.C diff --git a/glauber/macro/Chi2.C b/glauber/macro/Chi2.C new file mode 100644 index 0000000..e6452df --- /dev/null +++ b/glauber/macro/Chi2.C @@ -0,0 +1,108 @@ +void Chi2(TString InFileName) +{ + TFile *file = new TFile(InFileName); + TTree *tree = (TTree*) file->Get("test_tree"); + TTree *delta = new TTree("delta","delta"); + + Float_t f, mu, k; + Float_t chi2, chi2_error, chi2_min_error, CHI2, DELTA=1e10 ; + Float_t chi2_min=1e10; + Float_t f_min, mu_min, k_min; + Float_t f_max=0, mu_max=0, k_max=0; + Float_t f_delta=0, mu_delta=0, k_delta=0; + Float_t sigma; + + tree->SetBranchAddress("f", &f); + tree->SetBranchAddress("mu", &mu); + tree->SetBranchAddress("k", &k); + tree->SetBranchAddress("chi2", &chi2); + tree->SetBranchAddress("chi2_error", &chi2_error); + tree->SetBranchAddress("sigma", &sigma); + + delta->Branch("chi2", &chi2_min); + delta->Branch("delta_chi2", &chi2_min_error); + delta->Branch("mu", &mu_min); + delta->Branch("delta_mu", &mu_delta); + delta->Branch("k", &k_min); + delta->Branch("delta_k", &k_delta); + delta->Branch("f", &f_min); + delta->Branch("delta_f", &f_delta); + + Int_t n = tree->GetEntries(); + + std::cout << " GetEntries = " << n << std::endl; + + std::vector x; + std::vector y; + std::vector z; + + int X, Y; + for (Int_t i=0; iGetEntry(i); + + std::cout << " f = " << f << " mu = " << mu << " k = " << k << " sigma = " << sigma << " chi2 = " << chi2 << " chi2_error = " << chi2_error << std::endl; + + x.push_back(f); + y.push_back(k); + z.push_back(chi2); + + if (chi2SetName("#chi^{2} vs f, k"); + g->SetTitle("#chi^{2} vs f, k; f; k; #chi^{2}"); + + g->Draw("colz"); + + float F,K; + for (double xx=X; (xx<=X+0.2 && xx<=g->GetXmax()); xx+=0.001) { + for (double yy=Y; (yy<=Y+10.0 && yy<=g->GetYmax()); yy+=0.001) if (TMath::Abs(g->Interpolate(xx,yy)-CHI2)f_delta) f_delta=TMath::Abs(xx-f_min); + if (TMath::Abs(yy-k_min)>k_delta) k_delta=TMath::Abs(yy-k_min); + F=1e10; K=1e10; + for (Int_t ii=0; iiGetEntry(ii); if (TMath::Abs(f-xx)<=F && TMath::Abs(k-yy)<=K && (TMath::Abs(mu-mu_min))>mu_delta) mu_delta=TMath::Abs(mu-mu_min);} + DELTA=TMath::Abs(g->Interpolate(xx,yy)-CHI2); + } + for (double yy=Y; (yy>=Y-10.0 && yy>=g->GetYmin()); yy-=0.001) if (TMath::Abs(g->Interpolate(xx,yy)-CHI2)f_delta) f_delta=TMath::Abs(xx-f_min); + if (TMath::Abs(yy-k_min)>k_delta) k_delta=TMath::Abs(yy-k_min); + F=1e10; K=1e10; + for (Int_t ii=0; iiGetEntry(ii); if (TMath::Abs(f-xx)<=F && TMath::Abs(k-yy)<=K && (TMath::Abs(mu-mu_min))>mu_delta) mu_delta=TMath::Abs(mu-mu_min);} + DELTA=TMath::Abs(g->Interpolate(xx,yy)-CHI2); + } + } + + for (double xx=X; (xx>=X-0.2 && xx>=g->GetXmin()); xx-=0.001) { + for (double yy=Y; (yy<=Y+2.0 || yy<=g->GetYmax()); yy+=0.001) if (TMath::Abs(g->Interpolate(xx,yy)-CHI2)f_delta) f_delta=TMath::Abs(xx-f_min); + if (TMath::Abs(yy-k_min)>k_delta) k_delta=TMath::Abs(yy-k_min); + F=1e10; K=1e10; + for (Int_t ii=0; iiGetEntry(ii); if (TMath::Abs(f-xx)<=F && TMath::Abs(k-yy)<=K && (TMath::Abs(mu-mu_min))>mu_delta) mu_delta=TMath::Abs(mu-mu_min);} + DELTA=TMath::Abs(g->Interpolate(xx,yy)-CHI2); + } + for (double yy=Y; (yy>=Y-2.0 && yy>=g->GetYmin()); yy-=0.001) if (TMath::Abs(g->Interpolate(xx,yy)-CHI2)f_delta) f_delta=TMath::Abs(xx-f_min); + if (TMath::Abs(yy-k_min)>k_delta) k_delta=TMath::Abs(yy-k_min); + F=1e10; K=1e10; + for (Int_t ii=0; iiGetEntry(ii); if (TMath::Abs(f-xx)<=F && TMath::Abs(k-yy)<=K && (TMath::Abs(mu-mu_min))>mu_delta) mu_delta=TMath::Abs(mu-mu_min);} + DELTA=TMath::Abs(g->Interpolate(xx,yy)-CHI2); + } + } + + g->Draw("colz"); + + std::cout << " f = " << f_min << "+/-" << f_delta << " mu = " << mu_min << "+/-" << mu_delta << " k = " << k_min << "+/-" << k_delta << " chi2 = " << chi2_min << "+/-" << chi2_min_error << std::endl; + + TFile *f1 = new TFile("Fit_Errors_RPC.root", "recreate"); + delta->Write(); + g->Write(); + f1->Close(); + + +} From 2efe3f5c61e8e4fbcbeb5e9d2066995c24cd8e76 Mon Sep 17 00:00:00 2001 From: Ilya Segal Date: Wed, 18 Sep 2019 18:07:36 +0200 Subject: [PATCH 06/10] Update Chi2.C --- glauber/macro/Chi2.C | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/glauber/macro/Chi2.C b/glauber/macro/Chi2.C index e6452df..916f483 100644 --- a/glauber/macro/Chi2.C +++ b/glauber/macro/Chi2.C @@ -61,15 +61,15 @@ void Chi2(TString InFileName) g->Draw("colz"); float F,K; - for (double xx=X; (xx<=X+0.2 && xx<=g->GetXmax()); xx+=0.001) { - for (double yy=Y; (yy<=Y+10.0 && yy<=g->GetYmax()); yy+=0.001) if (TMath::Abs(g->Interpolate(xx,yy)-CHI2)GetXmax()); xx+=0.01) { + for (double yy=Y; (yy<=Y+10.0 && yy<=g->GetYmax()); yy+=0.01) if (TMath::Abs(g->Interpolate(xx,yy)-CHI2)f_delta) f_delta=TMath::Abs(xx-f_min); if (TMath::Abs(yy-k_min)>k_delta) k_delta=TMath::Abs(yy-k_min); F=1e10; K=1e10; for (Int_t ii=0; iiGetEntry(ii); if (TMath::Abs(f-xx)<=F && TMath::Abs(k-yy)<=K && (TMath::Abs(mu-mu_min))>mu_delta) mu_delta=TMath::Abs(mu-mu_min);} DELTA=TMath::Abs(g->Interpolate(xx,yy)-CHI2); } - for (double yy=Y; (yy>=Y-10.0 && yy>=g->GetYmin()); yy-=0.001) if (TMath::Abs(g->Interpolate(xx,yy)-CHI2)=Y-10.0 && yy>=g->GetYmin()); yy-=0.01) if (TMath::Abs(g->Interpolate(xx,yy)-CHI2)f_delta) f_delta=TMath::Abs(xx-f_min); if (TMath::Abs(yy-k_min)>k_delta) k_delta=TMath::Abs(yy-k_min); F=1e10; K=1e10; @@ -78,15 +78,15 @@ void Chi2(TString InFileName) } } - for (double xx=X; (xx>=X-0.2 && xx>=g->GetXmin()); xx-=0.001) { - for (double yy=Y; (yy<=Y+2.0 || yy<=g->GetYmax()); yy+=0.001) if (TMath::Abs(g->Interpolate(xx,yy)-CHI2)=X-0.2 && xx>=g->GetXmin()); xx-=0.01) { + for (double yy=Y; (yy<=Y+2.0 || yy<=g->GetYmax()); yy+=0.01) if (TMath::Abs(g->Interpolate(xx,yy)-CHI2)f_delta) f_delta=TMath::Abs(xx-f_min); if (TMath::Abs(yy-k_min)>k_delta) k_delta=TMath::Abs(yy-k_min); F=1e10; K=1e10; for (Int_t ii=0; iiGetEntry(ii); if (TMath::Abs(f-xx)<=F && TMath::Abs(k-yy)<=K && (TMath::Abs(mu-mu_min))>mu_delta) mu_delta=TMath::Abs(mu-mu_min);} DELTA=TMath::Abs(g->Interpolate(xx,yy)-CHI2); } - for (double yy=Y; (yy>=Y-2.0 && yy>=g->GetYmin()); yy-=0.001) if (TMath::Abs(g->Interpolate(xx,yy)-CHI2)=Y-2.0 && yy>=g->GetYmin()); yy-=0.01) if (TMath::Abs(g->Interpolate(xx,yy)-CHI2)f_delta) f_delta=TMath::Abs(xx-f_min); if (TMath::Abs(yy-k_min)>k_delta) k_delta=TMath::Abs(yy-k_min); F=1e10; K=1e10; From ea6273c254ed31854b38ee1befe5bb45a9a2923d Mon Sep 17 00:00:00 2001 From: Ilya Segal Date: Wed, 18 Sep 2019 18:18:50 +0200 Subject: [PATCH 07/10] Replace plot_chi2.C --- glauber/macro/plot_chi2.C | 164 +++++++++++++++++++------------------- 1 file changed, 83 insertions(+), 81 deletions(-) diff --git a/glauber/macro/plot_chi2.C b/glauber/macro/plot_chi2.C index 566e089..a98da14 100644 --- a/glauber/macro/plot_chi2.C +++ b/glauber/macro/plot_chi2.C @@ -1,84 +1,86 @@ - - -void plot_chi2(TString InFileName) { - gStyle->SetOptStat(0000); - - TString qa_filename = "/home/vklochkov/Data/na61/na61_30_qa.root"; - TString histoname = "reco_info/hMreco"; - - std::unique_ptr fqa{TFile::Open(qa_filename.Data())}; - TH1F *hData1 = (TH1F *) fqa->Get(histoname); - - TFile *file = new TFile(InFileName); - TTree *tree = (TTree *) file->Get("test_tree"); - - Float_t f, mu, k; - Float_t chi2; - Float_t sigma; - TH1F *h1 = new TH1F("h1", "", 500, 0, 500); - - tree->SetBranchAddress("f", &f); - tree->SetBranchAddress("mu", &mu); - tree->SetBranchAddress("k", &k); - tree->SetBranchAddress("chi2", &chi2); - tree->SetBranchAddress("sigma", &sigma); - tree->SetBranchAddress("histo", &h1); - - Int_t n = tree->GetEntries(); - - std::cout << " GetEntries = " << n << std::endl; - - std::vector x; - std::vector y; - std::vector z; - - for (Int_t i = 0; i < n; i++) { - tree->GetEntry(i); - - std::cout << " f = " << f << " mu = " << mu << " k = " << k << " sigma = " << sigma << " chi2 = " << chi2 - << std::endl; - - x.push_back(f); - y.push_back(k); - z.push_back(chi2); - - // if ( chi2 > 2.5 ) continue; - - if (false) { - TCanvas *c1 = new TCanvas("c1", "c1", 1200, 800); - - hData1->Draw(); - h1->Draw("same"); - h1->GetYaxis()->SetTitle("counts"); - h1->GetZaxis()->SetTitle("multiplicity"); - - gPad->SetLogy(); - gPad->Update(); - hData1->SetLineColor(kRed); - TLegend *leg1 = new TLegend(0.7, 0.75, 0.85, 0.89); - leg1->AddEntry(hData1, "M_{STS}", "l"); - leg1->AddEntry(h1, "MC-Glauber Fit", "l"); - leg1->Draw("same"); - gPad->Update(); - - c1->SaveAs("fit.root"); - c1->SaveAs("fit.C"); - + + +void plot_chi2(TString InFileName) +{ + gStyle->SetOptStat(0000); +/* + TString qa_filename = "/home/vklochkov/Data/na61/na61_30_qa.root"; + TString histoname = "reco_info/hMreco"; + + std::unique_ptr fqa{TFile::Open( qa_filename.Data() )}; + TH1F *hData1 = (TH1F*)fqa->Get( histoname ); +*/ + TFile *file = new TFile(InFileName); + TTree *tree = (TTree*) file->Get("test_tree"); + + Float_t f, mu, k; + Float_t chi2 ; + Float_t sigma; + TH1F *h1 = new TH1F ("h1", "", 500, 0, 500); + + tree->SetBranchAddress("f", &f); + tree->SetBranchAddress("mu", &mu); + tree->SetBranchAddress("k", &k); + tree->SetBranchAddress("chi2", &chi2); + tree->SetBranchAddress("sigma", &sigma); + tree->SetBranchAddress("histo", &h1); + + Int_t n = tree->GetEntries(); + + std::cout << " GetEntries = " << n << std::endl; + + std::vector x; + std::vector y; + std::vector z; + + for (Int_t i=0; iGetEntry(i); + + std::cout << " f = " << f << " mu = " << mu << " k = " << k << " sigma = " << sigma << " chi2 = " << chi2 << std::endl; + + x.push_back(f); + y.push_back(k); + z.push_back(chi2); + +// if ( chi2 > 2.5 ) continue; + + if ( false ) + { + TCanvas *c1 = new TCanvas ( "c1", "c1", 1200, 800 ); + +// hData1->Draw(); + h1->Draw("same"); + h1->GetYaxis()->SetTitle("counts"); + h1->GetZaxis()->SetTitle("multiplicity"); + + gPad->SetLogy(); + gPad->Update(); +// hData1->SetLineColor(kRed); + TLegend* leg1 = new TLegend(0.7, 0.75, 0.85, 0.89); +// leg1->AddEntry(hData1, "M_{STS}", "l"); + leg1->AddEntry(h1, "MC-Glauber Fit", "l"); + leg1->Draw("same"); + gPad->Update(); + + c1->SaveAs("fit.root"); + c1->SaveAs("fit.C"); + + } + +// Int_t jj; +// std::cin >> jj; +// if (jj == 0) break; } - - // Int_t jj; - // std::cin >> jj; - // if (jj == 0) break; - } - - TCanvas *c1 = new TCanvas("c1", "c1", 800, 800); - - TGraph2D *g = new TGraph2D(x.size(), &(x[0]), &(y[0]), &(z[0])); - g->SetName("#chi^{2} vs f, k"); - g->SetTitle("#chi^{2} vs f, k; f; k; #chi^{2}"); - - // gPad->SetLogz(); - g->Draw("colz"); - gPad->Update(); + + TCanvas *c1 = new TCanvas("c1", "c1", 800, 800); + + TGraph2D *g = new TGraph2D(x.size(), &(x[0]), &(y[0]), &(z[0])); + g->SetName("#chi^{2} vs f, k"); + g->SetTitle("#chi^{2} vs f, k; f; k; #chi^{2}"); + +// gPad->SetLogz(); + g->Draw("colz"); + gPad->Update(); } From 72734f49086c7996d538453f604eae5a7cd39884 Mon Sep 17 00:00:00 2001 From: Ilya Segal Date: Mon, 21 Oct 2019 12:17:41 +0200 Subject: [PATCH 08/10] Add Iteration method for finding mu and Tree eith all steps of fitting procedure --- glauber/Fitter.cpp | 116 +++++++++++++++++++++++++++++++++------------ 1 file changed, 87 insertions(+), 29 deletions(-) diff --git a/glauber/Fitter.cpp b/glauber/Fitter.cpp index e350b4c..a36d0ae 100644 --- a/glauber/Fitter.cpp +++ b/glauber/Fitter.cpp @@ -29,9 +29,9 @@ void Glauber::Fitter::Init(std::unique_ptr tree) it0++; } - if ( nEvents < 0 || nEvents > fSimTree->GetEntries() ) + if ( fnEvents < 0 || fnEvents > fSimTree->GetEntries() ) { - std::cout << "Init: *** ERROR - number of entries < 0 or less that number of entries in input tree" << std::endl; + std::cout << "Init: *** ERROR - number of entries < 0 or less than number of entries in input tree" << std::endl; std::cout << "Init: *** number of entries in input tree = " << fSimTree->GetEntries() << std::endl; exit(EXIT_FAILURE); } @@ -49,7 +49,7 @@ void Glauber::Fitter::Init(std::unique_ptr tree) std::map ::iterator it2=Glauber_ParametersMax.begin(); for (int j=0; jfirst) Glauber_Parameters_Histos.insert( std::pair ((it2->first), new TH1F (Form("%s_Histo", (gGlauberParameters[i].name).Data()), Form("%s;%s;counts", (gGlauberParameters[i].title).Data(), (gGlauberParameters[i].axis_title).Data()), 1.3*((it2->second)-(it1->second)+1)/gGlauberParameters[i].bin_value, 1.3*it1->second, 1.3*it2->second)) ); + for (int i=0; ifirst) Glauber_Parameters_Histos.insert( std::pair ((it2->first), new TH1F (Form("%s_Histo", (gGlauberParameters[i].name).Data()), Form("%s;%s;counts", (gGlauberParameters[i].title).Data(), (gGlauberParameters[i].axis_title).Data()), 1.3*((it2->second)-(it1->second))/gGlauberParameters[i].bin_value, 1.3*it1->second, 1.3*it2->second)) ); it1++; it2++; } @@ -58,7 +58,7 @@ void Glauber::Fitter::Init(std::unique_ptr tree) it1=Glauber_ParametersMin.begin(); it2=Glauber_ParametersMax.begin(); std::map ::iterator it3=Glauber_Parameters_Histos.begin(); - for (int i=0; iGetEntry(i); it0=Glauber_Parameters.begin(); @@ -91,7 +91,7 @@ void Glauber::Fitter::Init(std::unique_ptr tree) it3=Glauber_Parameters_Histos.begin(); for (int j=0; jfirst) Glauber_Parameters_VS_Multiplicity_Histos.insert( std::pair ((it2->first), new TH2F (Form("%s_VS_Multiplicity_Histo", (gGlauberParameters[i].name).Data()), Form("%s VS Multiplicity;Multiplicity;%s", (gGlauberParameters[i].title).Data(), (gGlauberParameters[i].axis_title).Data()), 1.3*fNbins, 0, 1.3*fMaxValue, 1.3*((it2->second)-(it1->second)+1)/gGlauberParameters[i].bin_value, 1.3*it1->second, 1.3*it2->second)) ); + for (int i=0; ifirst) Glauber_Parameters_VS_Multiplicity_Histos.insert( std::pair ((it2->first), new TH2F (Form("%s_VS_Multiplicity_Histo", (gGlauberParameters[i].name).Data()), Form("%s VS Multiplicity;Multiplicity;%s", (gGlauberParameters[i].title).Data(), (gGlauberParameters[i].axis_title).Data()), 1.3*fNbins, 0, 1.3*fMaxValue, 1.3*((it2->second)-(it1->second))/gGlauberParameters[i].bin_value, 1.3*it1->second, 1.3*it2->second)) ); it1++; it2++; } @@ -102,7 +102,7 @@ void Glauber::Fitter::Init(std::unique_ptr tree) it3=Glauber_Parameters_Histos.begin(); for (int j=0; jfirst) Glauber_Parameters_VS_Multiplicity_BestHistos.insert( std::pair ((it2->first), new TH2F (Form("%s_VS_Multiplicity_BestHisto", (gGlauberParameters[i].name).Data()), Form("%s VS Multiplicity;Multiplicity;%s", (gGlauberParameters[i].title).Data(), (gGlauberParameters[i].axis_title).Data()), 1.3*fNbins, 0, 1.3*fMaxValue, 1.3*((it2->second)-(it1->second)+1)/gGlauberParameters[i].bin_value, 1.3*it1->second, 1.3*it2->second)) ); + for (int i=0; ifirst) Glauber_Parameters_VS_Multiplicity_BestHistos.insert( std::pair ((it2->first), new TH2F (Form("%s_VS_Multiplicity_BestHisto", (gGlauberParameters[i].name).Data()), Form("%s VS Multiplicity;Multiplicity;%s", (gGlauberParameters[i].title).Data(), (gGlauberParameters[i].axis_title).Data()), 1.3*fNbins, 0, 1.3*fMaxValue, 1.3*((it2->second)-(it1->second))/gGlauberParameters[i].bin_value, 1.3*it1->second, 1.3*it2->second)) ); it1++; it2++; } @@ -148,7 +148,7 @@ void Glauber::Fitter::SetGlauberFitHisto (float f, float mu, float k, Bool_t Nor SetNBDhist(mu, k); std::unique_ptr htemp {(TH1F*)fNbdHisto.Clone("htemp")}; // WTF??? Not working without pointer - for (int i=0; iGetEntry(i); const int Na = int(Nancestors(f)); @@ -202,7 +202,7 @@ void Glauber::Fitter::NormalizeGlauberFit () * @param f parameter of Na * @param k NBD parameter */ -void Glauber::Fitter::FindMuGoldenSection (float *mu, float *chi2, float*chi2_error, float mu_min, float mu_max, float f, float k, int n) +void Glauber::Fitter::FindMuGoldenSection (TTree *tree, float *mu, float *chi2, float*chi2_error, int *n, float *sigma, float mu_min, float mu_max, float f, float k ) { const float phi {(1+TMath::Sqrt(5))/2}; @@ -216,12 +216,27 @@ void Glauber::Fitter::FindMuGoldenSection (float *mu, float *chi2, float*chi2_er SetGlauberFitHisto (f, mu_1, k); float chi2_mu1 = GetChi2 (); float chi2_mu1_error = GetChi2Error (); + + *mu = mu_1; + *chi2 = chi2_mu1; + *chi2_error = chi2_mu1; + *sigma = ( (*mu)/k + 1 ) * (*mu); + tree->Fill(); + std::cout << "n = " << (*n) << " f = " << f << " k = " << k << " mu = " << (*mu) << " chi2 = " << (*chi2) << " chi2_error = " << (*chi2_error) << std::endl; SetGlauberFitHisto (f, mu_2, k); float chi2_mu2 = GetChi2 (); float chi2_mu2_error = GetChi2Error (); - for (int j=0; jFill(); + std::cout << "n = " << (*n) << " f = " << f << " k = " << k << " mu = " << (*mu) << " chi2 = " << (*chi2) << " chi2_error = " << (*chi2_error) << std::endl; + + for (int j=0; j chi2_mu2) { @@ -232,6 +247,12 @@ void Glauber::Fitter::FindMuGoldenSection (float *mu, float *chi2, float*chi2_er SetGlauberFitHisto (f, mu_2, k); chi2_mu2 = GetChi2 (); chi2_mu2_error = GetChi2Error (); + *n = *n + 1; + *mu = mu_2; + *chi2 = chi2_mu2; + *chi2_error = chi2_mu2; + *sigma = ( (*mu)/k + 1 ) * (*mu); + tree->Fill(); } else { @@ -241,10 +262,16 @@ void Glauber::Fitter::FindMuGoldenSection (float *mu, float *chi2, float*chi2_er chi2_mu2 = chi2_mu1; SetGlauberFitHisto (f, mu_1, k); chi2_mu1 = GetChi2 (); - chi2_mu1_error = GetChi2Error (); + chi2_mu1_error = GetChi2Error (); + *n = *n + 1; + *mu = mu_2; + *chi2 = chi2_mu2; + *chi2_error = chi2_mu2; + *sigma = ( (*mu)/k + 1 ) * (*mu); + tree->Fill(); } - std::cout << "n = " << n+j << " f = " << f << " k = " << k << " mu1 = " << mu_1 << " mu2 = " << mu_2 << " chi2_mu1 = " << chi2_mu1 << " chi2_mu2 = " << chi2_mu2 << std::endl; + std::cout << "n = " << (*n) << " f = " << f << " k = " << k << " mu = " << (*mu) << " chi2 = " << (*chi2) << " chi2_error = " << (*chi2_error) << std::endl; } /* take min(mu) */ @@ -255,16 +282,44 @@ void Glauber::Fitter::FindMuGoldenSection (float *mu, float *chi2, float*chi2_er *chi2_error = (chi2_mu1 < chi2_mu2) ? chi2_mu1_error : chi2_mu2_error; } +/** + * + * @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 + * @param mu_max upper search edge for mean value NBD + * @param f parameter of Na + * @param k NBD parameter + */ +void Glauber::Fitter::FindMuIteration (TTree *tree, float mu, float *chi2, float *chi2_error, int *n, float *sigma, float f, float k) +{ + float chi2_best=1e10, chi2_error_best=1e10; + SetGlauberFitHisto (f, mu, k); + *sigma = ( mu/k + 1 ) * mu; + for (int fRange = fMinFitRange; fRange < (fMaxMultiplicity-fMinMultiplicity); fRange = fRange + fMultiplicityStep) + for (int fRangeCenter = fMinMultiplicity+fRange/2; fRangeCenter < (fMaxMultiplicity-fRange/2); fRangeCenter = fRangeCenter + fMultiplicityStep) { + fFitMinBin = fRangeCenter - fRange/2; + fFitMaxBin = fRangeCenter + fRange/2; + *chi2 = GetChi2 (); + *chi2_error = GetChi2Error (); + tree->Fill(); + std::cout << "n = " << (*n) << " f = " << f << " k = " << k << " mu = " << mu << " FitMinBin = " << fFitMinBin << " FitMaxBin = " << fFitMaxBin << " chi2 = " << (*chi2) << " chi2_error = " << (*chi2_error) << std::endl; + *n = *n + 1; + if ((*chi2) < chi2_best) { chi2_best=(*chi2); chi2_error_best=(*chi2_error); } + } + + /* take min(chi2) */ + *chi2 = chi2_best; + /* take min(chi2_error) */ + *chi2_error = chi2_error_best; +} + /** * Find the best match * * @param return value of best fit parameters - * @param f0 lower search edge for parameter of Na, for which chi2 will be calculated - * @param f1 upper search edge for parameter of Na, for which chi2 will be calculated - * @param k0 lower search edge for NBD parameter - * @param k1 upper search edge for NBD parameter */ -float Glauber::Fitter::FitGlauber (float *par, Float_t f0, Float_t f1, Int_t k0, Int_t k1) +float Glauber::Fitter::FitGlauber (float *par) { float f_fit{-1}; float mu_fit{-1}; @@ -281,7 +336,11 @@ float Glauber::Fitter::FitGlauber (float *par, Float_t f0, Float_t f1, Int_t k0, TTree* tree {new TTree("test_tree", "tree" )}; float f, mu, k, chi2, chi2_error, sigma; - + int n=1; + + tree->Branch("n", &n, "n/F"); + tree->Branch("MinBin", &fFitMinBin, "fFitMinBin/F"); + tree->Branch("MaxBin", &fFitMaxBin, "fFitMaxBin/F"); tree->Branch("f", &f, "f/F"); tree->Branch("mu", &mu, "mu/F"); tree->Branch("k", &k, "k/F"); @@ -289,22 +348,21 @@ float Glauber::Fitter::FitGlauber (float *par, Float_t f0, Float_t f1, Int_t k0, tree->Branch("chi2_error", &chi2_error, "chi2_error/F"); tree->Branch("sigma",&sigma,"sigma/F"); - int n=1; - for (float i=f0; i<=f1; i=i+0.1000000) + + + for (float i=f_fMin; i<=f_fMax; i=i+f_fStep) { - f=i; - for (int j=k0; j<=k1; j++) + f = i; + for (float j=f_kMin; j<=f_kMax; j=j+f_kStep) { mu = fMaxValue / NancestorsMax(f) ; k = j; - const float mu_min = 0.7*mu; - const float mu_max = 1.0*mu; + const float mu_min = f_MuMinPercentage*mu; + const float mu_max = f_MuMaxPercentage*mu; - FindMuGoldenSection (&mu, &chi2, &chi2_error, mu_min, mu_max, f, k, n); - n=n+fnIter; - sigma = ( mu/k + 1 ) * mu; - - tree->Fill(); + if (fFit_Mode == "GoldenSection") FindMuGoldenSection (tree, &mu, &chi2, &chi2_error, &n, &sigma, mu_min, mu_max, f, k); + else if (fFit_Mode == "Iteration") for (float mu=mu_min; mu<=mu_max; mu=mu+f_muStep) FindMuIteration (tree, mu, &chi2, &chi2_error, &n, &sigma, f, k); + else std::cout << "ERROR: ILLIGAL MU FINDING MODE" << std::endl; if (chi2 < Chi2Min) { @@ -472,7 +530,7 @@ std::unique_ptr Glauber::Fitter::GetModelHisto (const float range[2], TStr std::unique_ptr hModel(new TH1F ("hModel", "name", 100, fSimTree->GetMinimum(name), fSimTree->GetMaximum(name)) ); - for (int i=0; iGetEntry(i); const int Na = int(Nancestors(f)); From cd9e39e4ed6e3c91dfc979382641977df9ed91b2 Mon Sep 17 00:00:00 2001 From: Ilya Segal Date: Mon, 21 Oct 2019 12:18:40 +0200 Subject: [PATCH 09/10] Add Iteration method for finding mu and Tree with all steps of fitting procedure --- glauber/Fitter.h | 63 +++++++++++++++++++++++++++++++++++++++--------- 1 file changed, 52 insertions(+), 11 deletions(-) diff --git a/glauber/Fitter.h b/glauber/Fitter.h index ff9c859..da0d83d 100644 --- a/glauber/Fitter.h +++ b/glauber/Fitter.h @@ -4,7 +4,6 @@ @author Ilya Selyuzhenkov (ilya.selyuzhenkov@gmail.com) @brief Class to fit histo with Glauber based function */ - #ifndef GlauberFitter_H #define GlauberFitter_H 1 @@ -94,8 +93,9 @@ namespace Glauber void NormalizeGlauberFit (); // Function for normilize model histogram with data histogram void DrawHistos (Bool_t isSim = true, Bool_t isData = true, Bool_t isGlauber = false, Bool_t isNBD = false); // Function for drawing output histograms and creating output files - float FitGlauber (float *par, Float_t f0, Float_t f1, Int_t k0, Int_t k1); // Main function for fitting - void FindMuGoldenSection (Float_t *mu, Float_t *chi2, float*chi2_error, Float_t mu_min, Float_t mu_max, Float_t f, Float_t k, int n=0); // Function that find mu by iteration it to the minimum of chi2 + float FitGlauber (float *par); // Main function for fitting + void FindMuGoldenSection (TTree *tree, float *mu, float *chi2, float *chi2_error, int *n, float *sigma, float mu_min, float mu_max, float f, float k); // Function that find mu by iteration it to the minimum of chi2 + void FindMuIteration (TTree *tree, float mu, float *chi2, float *chi2_error, int *n, float *sigma, float f, float k); // Function that find mu by iteration it to the minimum of chi2 Float_t GetChi2 (void) const; // Function for chi2 counting Float_t GetChi2Error (void) const; // Function for chi2 error counting @@ -112,19 +112,33 @@ namespace Glauber // Setters // void SetInputHisto (const TH1F &h) { fDataHisto = h; } // Setting of data histogram - void SetFitMinBin (Int_t min) { fFitMinBin = min; } // Setting of not fitting low multiplicity region due to trigger bias, etc - void SetFitMaxBin (Int_t max) { fFitMaxBin = max; } // Setting of very large number to fit the whole histo - void SetnIter (Int_t nIter) { fnIter = nIter; } // Setting of number of iterations of finding better mu (with same f and k) + void SetFitMinBin (Int_t min) { fFitMinBin = min; } // Setting of minimal bin for fitting + void SetFitMaxBin (Int_t max) { fFitMaxBin = max; } // Setting of maximal bin for fitting + void SetnMuIter (Int_t nMuIter) { fnMuIter = nMuIter; } // Setting of number of iterations of finding better mu (with same f and k) void SetOutDirName (const TString name) { fOutDirName = name; } // Setting of the name of directory for output files void SetOutFileIDName (const TString name) { fOutFileIDName = name; } // Setting of ID of output files (ending of all output files of one job will have same ID) void SetAncestor_Mode (const TString Ancestor_Mode) { fAncestor_Mode = Ancestor_Mode; } // Setting of the mode for counting ancestors void SetFitFunction_Mode (const TString FitFunction_Mode) { fFitFunction_Mode = FitFunction_Mode; } // Setting of the mode for fitting function (COMING SOON) + void SetFit_Mode (const TString Fit_Mode) { fFit_Mode = Fit_Mode; } // Setting of the mode for fitting procedure ("GoldenSection" or "Iteration") void SetGlauber_filename (const TString Glauber_filename) { fGlauber_filename = Glauber_filename; } // Setting of input file with Glauber Tree void SetGlauber_treename (const TString Glauber_treename) { fGlauber_treename = Glauber_treename; } // Setting of the name of tree with GlauberMC model void SetDataHisto_filename (const TString DataHisto_filename) { fDataHisto_filename = DataHisto_filename; } // Setting of input file with data histo void SetDataHisto_name (const TString DataHisto_name) { fDataHisto_name = DataHisto_name; } // Setting of the name of data histo void SetMassNumber (Float_t A) { fA = A; } // Setting of mass number of projectile nucleus - void SetNEvents (Int_t Events) { nEvents = Events; } // Setting of number of events that you need from GlauberMC + void SetNEvents (Int_t Events) { fnEvents = Events; } // Setting of number of events that you need from GlauberMC + void SetMultiplicityStep ( Int_t MultiplicityStep ) { fMultiplicityStep = MultiplicityStep; }; // Setting of step of variating lower and upper multiplicity range + void SetMinMultiplicity ( Int_t MinMultiplicity ) { fMinMultiplicity = MinMultiplicity; }; // Setting of not fitting low multiplicity region due to trigger bias, etc + void SetMaxMultiplicity ( Int_t MaxMultiplicity ) { fMaxMultiplicity = MaxMultiplicity; }; // Setting of very large number to fit the whole histo + void SetMinFitRange ( Int_t MinFitRange ) { fMinFitRange = MinFitRange; }; // Setting of minimal multiplicity range for fitting + void Set_kMin ( Float_t kMin ) { f_kMin = kMin; }; // Setting of lower parameter k + void Set_kMax ( Float_t kMax ) { f_kMax = kMax; }; // Setting of upper parameter k + void Set_kStep ( Float_t kStep ) { f_kStep = kStep; }; // Setting of step of variating of parameter k + void Set_fMin ( Float_t fMin ) { f_fMin = fMin; }; // Setting of lower parameter f + void Set_fMax ( Float_t fMax ) { f_fMax = fMax; }; // Setting of upper parameter f + void Set_fStep ( Float_t fStep ) { f_fStep = fStep; }; // Setting of step of variating of parameter f + void Set_muStep ( Float_t muStep ) { f_muStep = muStep; }; // Setting of step of variating of parameter mu + void Set_MuMinPercentage ( Float_t MuMinPercentage ) { f_MuMinPercentage = MuMinPercentage; }; // Setting of lower mu in percent of mu=(MaxMultiplicity/MaxAncestors) for first iteration + void Set_MuMaxPercentage ( Float_t MuMaxPercentage ) { f_MuMaxPercentage = MuMaxPercentage; }; // Setting of upper mu in percent of mu=(MaxMultiplicity/MaxAncestors) for first iteration // // Getters @@ -134,8 +148,20 @@ namespace Glauber TString GetGlauber_treename () { return fGlauber_treename; } // Returns the name of tree with TString GetDataHisto_filename () { return fDataHisto_filename; } // Returns input file with data histo TString GetDataHisto_name () { return fDataHisto_name; } // Returns the name of data histo + Int_t GetFitMinBin () const { return fFitMinBin; } // Returns not fitting low multiplicity region due to trigger bias, etc Int_t GetFitMaxBin () const { return fFitMaxBin; } // Returns very large number to fit the whole histo + Int_t GetMultiplicityStep () const { return fMultiplicityStep; }; // Returns step of variating lower and upper multiplicity range + Int_t GetMinMultiplicity () const { return fMinMultiplicity; }; // Returns not fitting low multiplicity region due to trigger bias, etc + Int_t GetMaxMultiplicity () const { return fMaxMultiplicity; }; // Returns very large number to fit the whole histo + Int_t GetMinFitRange () const { return fMinFitRange; }; // Returns minimal multiplicity range for fitting + Float_t Get_kMin () const { return f_kMin; }; // Returns lower parameter k + Float_t Get_kMax () const { return f_kMax; }; // Returns upper parameter k + Float_t Get_kStep () const { return f_kStep; }; // Returns step of variating of parameter + Float_t Get_fMin () const { return f_fMin; }; // Returns lower parameter f + Float_t Get_fMax () const { return f_fMax; }; // Returns upper parameter f + Float_t Get_fStep () const { return f_fStep; }; // Returns step of variating of parameter f + Float_t Get_muStep () const { return f_muStep; }; // Returns step of variating of parameter mu TH1F GetGlauberFitHisto () const { return fGlauberFitHisto; } // Returns model histogtam TH1F GetBestFitHisto () const { return fBestFitHisto; } // Returns model histogtam with minimal chi2 @@ -221,14 +247,29 @@ namespace Glauber Int_t fNbins{-1 }; - Int_t fFitMinBin{-1}; - Int_t fFitMaxBin{-1}; - Int_t nEvents; + Int_t fnEvents; - Int_t fnIter; + Int_t fnMuIter{2}; + + Int_t fFitMinBin{-1}; + Int_t fFitMaxBin{-1}; + Int_t fMultiplicityStep{1}; + Int_t fMinMultiplicity{1}; + Int_t fMaxMultiplicity{100}; + Int_t fMinFitRange{4}; + Float_t f_kMin{1.0}; + Float_t f_kMax{1.0}; + Float_t f_kStep{0.2}; + Float_t f_fMin{0.0}; + Float_t f_fMax{1.0}; + Float_t f_fStep{0.1}; + Float_t f_muStep{0.01}; + Float_t f_MuMinPercentage{0.7}; + Float_t f_MuMaxPercentage{1.0}; TString fAncestor_Mode{"Default"}; TString fFitFunction_Mode{"NBD"}; + TString fFit_Mode{"GoldenSection"}; TString fGlauber_filename{""}; TString fGlauber_treename{""}; TString fDataHisto_filename{""}; From fadba71124cd02748d5abcb2f37be0125b62bc24 Mon Sep 17 00:00:00 2001 From: Ilya Segal Date: Mon, 21 Oct 2019 12:19:13 +0200 Subject: [PATCH 10/10] Add Iteration method for finding mu and Tree with all steps of fitting procedure --- glauber/main.cpp | 85 +++++++++++++++++++++++++----------------------- 1 file changed, 45 insertions(+), 40 deletions(-) diff --git a/glauber/main.cpp b/glauber/main.cpp index 2ead99e..dd8efb4 100644 --- a/glauber/main.cpp +++ b/glauber/main.cpp @@ -15,52 +15,61 @@ using namespace Glauber; int main(int argc, char *argv[]) { - if (argc < 2) + if (argc == 0) { std::cout << "Wrong number of parameters! Executable usage:" << std::endl; - std::cout << " ./glauber f0 k0" << std::endl; + std::cout << " ./glauber" << std::endl; return -1; } - const Float_t f0 = atof( argv[1]); - std::cout << "f=" << f0 << std::endl; - const Float_t f1 = atof( argv[2]); - std::cout << "f=" << f1 << std::endl; - const Int_t k0 = atoi( argv[3] ); - std::cout << "k0=" << k0 << std::endl; - const Int_t k1 = atoi( argv[4] ); - std::cout << "k1=" << k1 << std::endl; + Fitter fitter; - - // ***************************************** - // Modify this part according to your needs - // ***************************************** + + //**********************************************************************************************************************************************************************************************// + // Modify this part according to your needs // + //**********************************************************************************************************************************************************************************************// /// |Ancestor_Mode| function for Na | - /// | Default | f*Npart + (1-f)*Ncoll | + /// | Default | f*Npart + (1-f)*Ncoll | /// | PSD | f-Npart | /// | Npart | Npart^f | /// | Ncoll | Ncoll^f | - - - fitter.SetFitMinBin (60); // Set not fitting low multiplicity region due to trigger bias, etc - fitter.SetFitMaxBin (200); // Set very large number to fit the whole histo - fitter.SetGlauber_filename ("/home/vad/NA61/Centrality/glau_PbPb_13GeV.root"); // Set input file with Glauber Tree - fitter.SetGlauber_treename ("nt_Pb_Pb"); // Set name of tree with GlauberMC model - fitter.SetDataHisto_filename ("/home/vad/NA61/QA/qa_new.root"); // Set input file with data histo - fitter.SetDataHisto_name ("reco_info/hMreco"); // Set name of data histo - fitter.SetOutDirName ("."); // Set name of directory for output files - fitter.SetOutFileIDName ("_M13new_5"); // Set ID of output files (ending of all output files of one job will have same ID) - fitter.SetAncestor_Mode ("Default"); // Set mode for counting ancestors - fitter.SetFitFunction_Mode ("NBD"); // Set mode for fitting function (COMING SOON) - fitter.SetnIter (10); // Set number of iterations of finding better mu (with same f and k) - fitter.SetMassNumber (207); // Set mass number of projectile nucleus - fitter.SetNEvents (999999); // Set number of events that you need from GlauberMC (or you can uncomment the same setter a few lines down with default number of events: 10 times the number of events in the data histo in the fitting range) - - std::cout << "min_bin=" << fitter.GetFitMinBin() << " max_bin=" << fitter.GetFitMaxBin() << std::endl; - - // ***************************************** - // ***************************************** + + fitter.SetMultiplicityStep (1); // Set step of variating lower and upper multiplicity range + fitter.SetMinMultiplicity (5); // Set not fitting low multiplicity region due to trigger bias, etc + fitter.SetMaxMultiplicity (200); // Set very large number to fit the whole histo + fitter.SetMinFitRange (190); // Set minimal multiplicity range for fitting + fitter.Set_kMin (95.0); // Set lower parameter k + fitter.Set_kMax (100.0); // Set upper parameter k + fitter.Set_kStep (0.5); // Set step of variating of parameter k + fitter.Set_fMin (0.9); // Set lower parameter f + fitter.Set_fMax (1.0); // Set upper parameter f + fitter.Set_fStep (0.1); // Set of step of variating of parameter f + fitter.SetGlauber_filename ("/home/vad/NA61/Centrality/glau_PbPb_13GeV.root"); // Set input file with Glauber Tree + fitter.SetGlauber_treename ("nt_Pb_Pb"); // Set name of tree with GlauberMC model + fitter.SetDataHisto_filename ("/home/vad/NA61/QA/qa_new.root"); // Set input file with data histo + fitter.SetDataHisto_name ("reco_info/hMreco"); // Set name of data histo + fitter.SetOutDirName ("."); // Set name of directory for output files + fitter.SetOutFileIDName ("_MTEST"); // Set ID of output files (ending of all output files of one job will have same ID) + fitter.SetAncestor_Mode ("Default"); // Set mode for counting ancestors + fitter.SetFitFunction_Mode ("NBD"); // Set mode for fitting function (COMING SOON) + fitter.SetFit_Mode ("Iteration"); // Set mode for fitting procedure ("GoldenSection" or "Iteration") + fitter.SetnMuIter (10); // Set number of iterations of finding better mu (with same f and k) + fitter.SetMassNumber (207); // Set mass number of projectile nucleus + fitter.SetNEvents (999999); // Set number of events that you need from GlauberMC (or you can uncomment the same setter a few lines down with default number of events: 10 times the number of events in the data histo in the fitting range) + + //**********************************************************************************************************************************************************************************************// + //**********************************************************************************************************************************************************************************************// + + std::cout << "f0=" << fitter.Get_fMin() << std::endl; + std::cout << "f1=" << fitter.Get_fMax() << std::endl; + std::cout << "f_step=" << fitter.Get_fStep() << std::endl; + std::cout << "k0=" << fitter.Get_kMin() << std::endl; + std::cout << "k1=" << fitter.Get_kMax() << std::endl; + std::cout << "k_step=" << fitter.Get_kStep() << std::endl; + std::cout << "min_bin=" << fitter.GetMinMultiplicity() << std::endl; + std::cout << "max_bin=" << fitter.GetMaxMultiplicity() << std::endl; + std::cout << "bin_step=" << fitter.GetMultiplicityStep() << std::endl; std::unique_ptr glauber_file{ TFile::Open(fitter.GetGlauber_filename(), "read") }; std::unique_ptr glauber_tree{ (TTree*) glauber_file->Get(fitter.GetGlauber_treename()) }; @@ -76,11 +85,7 @@ int main(int argc, char *argv[]) float par[4]; float chi2=1e10; - chi2 = fitter.FitGlauber(par, f0, f1, k0, k1); - - std::cout << "f = " << par[0] << " mu = " << par[1] << " k = " << par[2] << " chi2 = " << chi2 << " chi2_error = " << par[3] << std::endl; - -// Glauber::DrawHistos(fitter, true, true, true, true); + chi2 = fitter.FitGlauber(par); DrawHistos(fitter, par, chi2, true, true, true, true);