diff --git a/Common/Tools/Multiplicity/multGlauberNBDFitter.cxx b/Common/Tools/Multiplicity/multGlauberNBDFitter.cxx index ebabca167a1..18b8faeb4fa 100644 --- a/Common/Tools/Multiplicity/multGlauberNBDFitter.cxx +++ b/Common/Tools/Multiplicity/multGlauberNBDFitter.cxx @@ -8,7 +8,11 @@ // In applying this license CERN does not waive the privileges and immunities // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -// + +/// \file multGlauberNBDFitter.cxx +/// \brief Helper class to perform Glauber fits +/// \author David Dobrigkeit Chinellato + /********************************************** * * Class meant to do Glauber+NBD fits @@ -27,6 +31,8 @@ #include "multGlauberNBDFitter.h" +#include + #include #include #include @@ -34,8 +40,6 @@ #include #include #include -#include -#include #include #include #include @@ -43,73 +47,27 @@ #include #include -#include - -#include // FIXME -using namespace std; +#include ClassImp(multGlauberNBDFitter); -multGlauberNBDFitter::multGlauberNBDFitter() : TNamed(), - fNBD(0x0), - fhNanc(0x0), - fhNpNc(0x0), - fhV0M(0x0), - ffChanged(kTRUE), - fCurrentf(-1), - fAncestorMode(2), - fNpart(0x0), - fNcoll(0x0), - fContent(0x0), - fNNpNcPairs(-1), - fMaxNpNcPairs(1000000), - fMu(45), - fdMu(0.0), - fk(1.5), - ff(0.8), - fnorm(100), - fFitOptions("R0"), - fFitNpx(5000) -{ - // Constructor - fNpart = new Double_t[fMaxNpNcPairs]; - fNcoll = new Double_t[fMaxNpNcPairs]; - fContent = new Long_t[fMaxNpNcPairs]; - - // Ancestor histo - fhNanc = new TH1D("fhNanc", "", 1000, -0.5, 999.5); - - // NBD - fNBD = new TF1("fNBD", "ROOT::Math::negative_binomial_pdf(x,[0],[1])", 0, 45000); - fNBD->SetNpx(45000); - - // master function - fGlauberNBD = new TF1("fGlauberNBD", this, &multGlauberNBDFitter::ProbDistrib, - 0, 50000, 4, "multGlauberNBDFitter", "ProbDistrib"); - fGlauberNBD->SetParameter(0, fMu); - fGlauberNBD->SetParameter(1, fk); - fGlauberNBD->SetParameter(2, ff); - fGlauberNBD->SetParameter(3, fnorm); - - fGlauberNBD->SetParName(0, "mu"); - fGlauberNBD->SetParName(1, "k"); - fGlauberNBD->SetParName(2, "f"); - fGlauberNBD->SetParName(3, "norm"); - fGlauberNBD->SetParName(4, "dMu/dNanc"); -} - +multGlauberNBDFitter::multGlauberNBDFitter() : multGlauberNBDFitter("multGlauberNBDFitter") {} multGlauberNBDFitter::multGlauberNBDFitter(const char* name, const char* title) : TNamed(name, title), - fNBD(0x0), - fhNanc(0x0), - fhNpNc(0x0), - fhV0M(0x0), - ffChanged(kTRUE), + fNBD(new TF1("fNBD", "ROOT::Math::negative_binomial_pdf(x,[0],[1])", 0, 45000)), + fGlauberNBD(nullptr), + fTrentoNBD(nullptr), + fNBDFitterMode(NBDFitterMode::Glauber), + fhNanc(nullptr), + fhNpNc(nullptr), + fhNSources(nullptr), + fhV0M(nullptr), + ffChanged(true), fCurrentf(-1), - fAncestorMode(2), - fNpart(0x0), - fNcoll(0x0), - fContent(0x0), + fAncestorMode(AncestorMode::Continuous), + fNpart(nullptr), + fNcoll(nullptr), + fContent(nullptr), fNNpNcPairs(-1), fMaxNpNcPairs(1000000), fMu(45), @@ -121,134 +79,182 @@ multGlauberNBDFitter::multGlauberNBDFitter(const char* name, const char* title) fFitNpx(5000) { // Named constructor - fNpart = new Double_t[fMaxNpNcPairs]; - fNcoll = new Double_t[fMaxNpNcPairs]; - fContent = new Long_t[fMaxNpNcPairs]; - - // Ancestor histo - // fhNanc = new TH1D("fhNanc", "", fAncestorMode==2?10000:1000, -0.5, 999.5); + fNpart = new double[fMaxNpNcPairs]; + fNcoll = new double[fMaxNpNcPairs]; + fContent = new int64_t[fMaxNpNcPairs]; // NBD - fNBD = new TF1("fNBD", "ROOT::Math::negative_binomial_pdf(x,[0],[1])", 0, 45000); fNBD->SetNpx(45000); // master function - fGlauberNBD = new TF1("fGlauberNBD", this, &multGlauberNBDFitter::ProbDistrib, - 0, 50000, 5, "multGlauberNBDFitter", "ProbDistrib"); - fGlauberNBD->SetParameter(0, fMu); - fGlauberNBD->SetParameter(1, fk); - fGlauberNBD->SetParameter(2, ff); - fGlauberNBD->SetParameter(3, fnorm); - - fGlauberNBD->SetParName(0, "mu"); - fGlauberNBD->SetParName(1, "k"); - fGlauberNBD->SetParName(2, "f"); - fGlauberNBD->SetParName(3, "norm"); - fGlauberNBD->SetParName(4, "dMu/dNanc"); + InitGlauberNBD(fMu, fk, ff, fnorm); } + //________________________________________________________________ multGlauberNBDFitter::~multGlauberNBDFitter() { // Destructor - if (fNBD) { - delete fNBD; - fNBD = 0x0; - } - if (fhNanc) { - delete fhNanc; - fhNanc = 0x0; - } - if (fhNpNc) { - delete fhNpNc; - fhNpNc = 0x0; - } - if (fNpart) - delete[] fNpart; - if (fNcoll) - delete[] fNcoll; - if (fContent) - delete[] fContent; + delete fNBD; + delete fhNanc; + delete fhNpNc; + delete[] fNpart; + delete[] fNcoll; + delete[] fContent; +} + +void multGlauberNBDFitter::InitGlauberNBD(const float mu, const float k, const float f, const float norm) +{ + fNBDFitterMode = NBDFitterMode::Glauber; + fGlauberNBD = new TF1("fGlauberNBD", this, &multGlauberNBDFitter::GlauberProbDistrib, + 0, 50000, 5, "multGlauberNBDFitter", "GlauberProbDistrib"); + fGlauberNBD->SetParameter(Index(FitPar::mu), mu); + fGlauberNBD->SetParameter(Index(FitPar::k), k); + fGlauberNBD->SetParameter(Index(FitPar::f), f); + fGlauberNBD->SetParameter(Index(FitPar::norm), norm); + + fGlauberNBD->SetParName(Index(FitPar::mu), "mu"); + fGlauberNBD->SetParName(Index(FitPar::k), "k"); + fGlauberNBD->SetParName(Index(FitPar::f), "f"); + fGlauberNBD->SetParName(Index(FitPar::norm), "norm"); + fGlauberNBD->SetParName(Index(FitPar::dMu), "dMu/dNanc"); +} + +void multGlauberNBDFitter::InitTrentoNBD(const float mu, const float k, const float norm) +{ + fNBDFitterMode = NBDFitterMode::Trento; + fGlauberNBD = nullptr; + fTrentoNBD = new TF1("fTrentoNBD", this, &multGlauberNBDFitter::TrentoProbDistrib, + 0, 50000, 5, "multGlauberNBDFitter", "TrentoProbDistrib"); + fTrentoNBD->SetParameter(Index(FitPar::mu), mu); + fTrentoNBD->SetParameter(Index(FitPar::k), k); + fTrentoNBD->SetParameter(Index(FitPar::norm), norm); + fTrentoNBD->FixParameter(Index(FitPar::dMu), 0); + + fTrentoNBD->SetParName(Index(FitPar::mu), "mu"); + fTrentoNBD->SetParName(Index(FitPar::k), "k"); + fTrentoNBD->SetParName(Index(FitPar::f), "f"); + fTrentoNBD->SetParName(Index(FitPar::norm), "norm"); + fTrentoNBD->SetParName(Index(FitPar::dMu), "dMu/dNanc"); } //______________________________________________________ -Double_t multGlauberNBDFitter::ProbDistrib(Double_t* x, Double_t* par) +double multGlauberNBDFitter::GlauberProbDistrib(const double* x, const double* par) // Master fitter function { - Double_t lMultValue = x[0]; - Double_t lProbability = 0.0; - ffChanged = kTRUE; - const Double_t lAlmost0 = 1.e-13; + double lMultValue = x[0]; + double lProbability = 0.0; + ffChanged = true; + static constexpr double Almost0 = 1.e-13; // Comment this line in order to make the code evaluate Nancestor all the time - if (TMath::Abs(fCurrentf - par[2]) < lAlmost0) - ffChanged = kFALSE; + if (std::abs(fCurrentf - par[Index(FitPar::f)]) < Almost0) { + ffChanged = false; + } //______________________________________________________ // Recalculate the ancestor distribution in case f changed if (ffChanged) { - fCurrentf = par[2]; + fCurrentf = par[Index(FitPar::f)]; fhNanc->Reset(); for (int ibin = 0; ibin < fNNpNcPairs; ibin++) { - Double_t lOption0 = (Int_t)(fNpart[ibin] * par[2] + fNcoll[ibin] * (1.0 - par[2])); - Double_t lOption1 = TMath::Floor(fNpart[ibin] * par[2] + fNcoll[ibin] * (1.0 - par[2]) + 0.5); - Double_t lOption2 = (fNpart[ibin] * par[2] + fNcoll[ibin] * (1.0 - par[2])); - if (fAncestorMode == 0) + double lOption0 = static_cast(fNpart[ibin] * par[Index(FitPar::f)] + fNcoll[ibin] * (1.0 - par[Index(FitPar::f)])); + double lOption1 = std::floor(fNpart[ibin] * par[Index(FitPar::f)] + fNcoll[ibin] * (1.0 - par[Index(FitPar::f)]) + 0.5); + double lOption2 = (fNpart[ibin] * par[Index(FitPar::f)] + fNcoll[ibin] * (1.0 - par[Index(FitPar::f)])); + if (fAncestorMode == AncestorMode::Truncated) { fhNanc->Fill(lOption0, fContent[ibin]); - if (fAncestorMode == 1) + } + if (fAncestorMode == AncestorMode::Rounded) { fhNanc->Fill(lOption1, fContent[ibin]); - if (fAncestorMode == 2) + } + if (fAncestorMode == AncestorMode::Continuous) { fhNanc->Fill(lOption2, fContent[ibin]); + } } if (fhNanc->Integral() < 1) { - cout << "ERROR: ANCESTOR HISTOGRAM EMPTY" << endl; - cout << "Will not do anything. Call InitializeNpNc if you want to plot without fitting" << endl; + LOG(info) << "ERROR: ANCESTOR HISTOGRAM EMPTY"; + LOG(info) << "Will not do anything. Call InitNpNc if you want to plot without fitting"; return 0; } fhNanc->Scale(1. / fhNanc->Integral()); } //______________________________________________________ // Actually evaluate function - Int_t lStartBin = fhNanc->FindBin(0.0) + 1; - for (Long_t iNanc = lStartBin; iNanc < fhNanc->GetNbinsX() + 1; iNanc++) { - Double_t lNancestors = fhNanc->GetBinCenter(iNanc); - Double_t lNancestorCount = fhNanc->GetBinContent(iNanc); - // if(lNancestorCount<1e-12&&lNancestors>10) break; + int lStartBin = fhNanc->FindBin(0.0) + 1; + for (int64_t iNanc = lStartBin; iNanc < fhNanc->GetNbinsX() + 1; ++iNanc) { + double lNancestors = fhNanc->GetBinCenter(iNanc); + double lNancestorCount = fhNanc->GetBinContent(iNanc); // allow for variable mu in case requested - Double_t lThisMu = (((Double_t)lNancestors)) * (par[0] + par[4] * lNancestors); - Double_t lThisk = (((Double_t)lNancestors)) * par[1]; - Double_t lpval = TMath::Power(1.0 + lThisMu / lThisk, -1); - fNBD->SetParameter(1, lThisk); - fNBD->SetParameter(0, lpval); - Double_t lMult = 0.0; - if (lMultValue > 1e-6) - lMult = fAncestorMode != 2 ? fNBD->Eval(lMultValue) : ContinuousNBD(lMultValue, lThisMu, lThisk); + double lThisMu = lNancestors * (par[Index(FitPar::mu)] + par[Index(FitPar::dMu)] * lNancestors); + double lThisk = lNancestors * par[Index(FitPar::k)]; + double lpval = std::pow(1.0 + lThisMu / lThisk, -1); + fNBD->SetParameter(Index(NBDPar::k), lThisk); + fNBD->SetParameter(Index(NBDPar::p), lpval); + + double lMult = 0.0; + static constexpr double MultTolerance = 1e-6; + if (lMultValue > MultTolerance) { + lMult = fAncestorMode != AncestorMode::Continuous ? fNBD->Eval(lMultValue) : ContinuousNBD(lMultValue, lThisMu, lThisk); + } lProbability += lNancestorCount * lMult; } //______________________________________________________ - return par[3] * lProbability; + return par[Index(FitPar::norm)] * lProbability; +} + +double multGlauberNBDFitter::TrentoProbDistrib(const double* x, const double* par) +// Master fitter function +{ + double lMultValue = x[0]; + double lProbability = 0.0; + //______________________________________________________ + // Actually ealuate function + for (int64_t iNSrc = 1; iNSrc < fhNSources->GetNbinsX() + 1; ++iNSrc) { + double lNsources = fhNSources->GetBinCenter(iNSrc); + double lThisMu = lNsources * par[Index(FitPar::mu)]; + double lThisk = lNsources * par[Index(FitPar::k)]; + double lpval = std::pow(1 + lThisMu / lThisk, -1); + fNBD->SetParameter(Index(NBDPar::k), lThisk); + fNBD->SetParameter(Index(NBDPar::p), lpval); + double lMult = fNBD->Eval(lMultValue); + lProbability += fhNSources->GetBinContent(fhNSources->FindBin(iNSrc)) * lMult; + } + //______________________________________________________ + return par[Index(FitPar::norm)] * lProbability; } //________________________________________________________________ -Bool_t multGlauberNBDFitter::SetNpartNcollCorrelation(TH2* hNpNc) +bool multGlauberNBDFitter::SetNpartNcollCorrelation(TH2* hNpNc) { - Bool_t lReturnValue = kTRUE; + bool lReturnValue = true; if (hNpNc) { - fhNpNc = reinterpret_cast(hNpNc); + fhNpNc = hNpNc; } else { - lReturnValue = kFALSE; + lReturnValue = false; } return lReturnValue; } //________________________________________________________________ -Bool_t multGlauberNBDFitter::SetInputV0M(TH1* hV0M) +bool multGlauberNBDFitter::SetNSources(TH1* hNSources) { - Bool_t lReturnValue = kTRUE; + bool lReturnValue = true; + if (hNSources) { + fhNSources = hNSources; + } else { + lReturnValue = false; + } + return lReturnValue; +} + +//________________________________________________________________ +bool multGlauberNBDFitter::SetInputV0M(TH1* hV0M) +{ + bool lReturnValue = true; if (hV0M) { - fhV0M = reinterpret_cast(hV0M); + fhV0M = hV0M; } else { - lReturnValue = kFALSE; + lReturnValue = false; } return lReturnValue; } @@ -265,20 +271,30 @@ TF1* multGlauberNBDFitter::GetGlauberNBD() return fGlauberNBD; } +TF1* multGlauberNBDFitter::GetTrentoNBD() +{ + return fTrentoNBD; +} + //________________________________________________________________ -void multGlauberNBDFitter::SetFitRange(Double_t lMin, Double_t lMax) +void multGlauberNBDFitter::SetFitRange(const double lMin, const double lMax) { - fGlauberNBD->SetRange(lMin, lMax); + if (fGlauberNBD) { + fGlauberNBD->SetRange(lMin, lMax); + } + if (fTrentoNBD) { + fTrentoNBD->SetRange(lMin, lMax); + } } //________________________________________________________________ -void multGlauberNBDFitter::SetFitOptions(TString lOpt) +void multGlauberNBDFitter::SetFitOptions(const TString& lOpt) { fFitOptions = lOpt; } //________________________________________________________________ -void multGlauberNBDFitter::SetFitNpx(Long_t lNpx) +void multGlauberNBDFitter::SetFitNpx(const int64_t lNpx) { fFitNpx = lNpx; } @@ -286,79 +302,120 @@ void multGlauberNBDFitter::SetFitNpx(Long_t lNpx) //________________________________________________________________ void multGlauberNBDFitter::InitAncestor() { - if (!fhNanc) - fhNanc = new TH1D("fhNanc", "", fAncestorMode == 2 ? 10000 : 1000, -0.5, 999.5); + if (!fhNanc) { + fhNanc = new TH1D("fhNanc", "", fAncestorMode == AncestorMode::Continuous ? 10000 : 1000, -0.5, 999.5); + } } //________________________________________________________________ -Bool_t multGlauberNBDFitter::DoFit() +bool multGlauberNBDFitter::DoFit() { - InitAncestor(); // Try very hard, please TVirtualFitter::SetMaxIterations(5000000); + auto timer = new TStopwatch(); + timer->Start(true); + + bool lReturnValue = false; + switch (fNBDFitterMode) { + case NBDFitterMode::Glauber: + LOG(info) << "Doing Glauber fit!!"; + lReturnValue = DoGlauberFit(); + break; + case NBDFitterMode::Trento: + LOG(info) << "Doing Trento fit!!"; + lReturnValue = DoTrentoFit(); + break; + + default: + break; + } + + timer->Stop(); + double lTotalTime = timer->RealTime(); + if (lReturnValue) { + LOG(info) << "---> Fitting succeeded after " << lTotalTime << " seconds"; + } else { + LOG(info) << "---> Fitting failed after " << lTotalTime << " seconds"; + } + return lReturnValue; +} + +bool multGlauberNBDFitter::DoGlauberFit() +{ + InitAncestor(); if (!InitializeNpNc()) { - cout << "---> Initialization of Npart x Ncoll correlation info failed!" << endl; - return kFALSE; + LOG(info) << "---> Initialization of Npart x Ncoll correlation info failed!"; + return false; } - TStopwatch* timer = new TStopwatch(); - timer->Start(kTRUE); - if (fAncestorMode == 0) - cout << "---> Config: Nancestors will be truncated" << endl; - if (fAncestorMode == 1) - cout << "---> Config: Nancestors will be rounded" << endl; - if (fAncestorMode == 2) - cout << "---> Config: Nancestors will be taken as float" << endl; - cout << "---> Now fitting, please wait..." << endl; + if (fAncestorMode == AncestorMode::Truncated) { + LOG(info) << "---> Config: Nancestors will be truncated"; + } + if (fAncestorMode == AncestorMode::Rounded) { + LOG(info) << "---> Config: Nancestors will be rounded"; + } + if (fAncestorMode == AncestorMode::Continuous) { + LOG(info) << "---> Config: Nancestors will be taken as float"; + } + LOG(info) << "---> Now fitting, please wait..."; fGlauberNBD->SetNpx(fFitNpx); TFitResultPtr fitptr; fFitOptions.Append("S"); fitptr = fhV0M->Fit("fGlauberNBD", fFitOptions.Data()); - timer->Stop(); - Double_t lTotalTime = timer->RealTime(); - cout << "---> Fitting took " << lTotalTime << " seconds" << endl; - - fMu = fGlauberNBD->GetParameter(0); - fk = fGlauberNBD->GetParameter(1); - ff = fGlauberNBD->GetParameter(2); - fnorm = fGlauberNBD->GetParameter(3); - fdMu = fGlauberNBD->GetParameter(4); + fMu = fGlauberNBD->GetParameter(Index(FitPar::mu)); + fk = fGlauberNBD->GetParameter(Index(FitPar::k)); + ff = fGlauberNBD->GetParameter(Index(FitPar::f)); + fnorm = fGlauberNBD->GetParameter(Index(FitPar::norm)); + fdMu = fGlauberNBD->GetParameter(Index(FitPar::dMu)); + return fitptr.Get()->IsValid(); +} +bool multGlauberNBDFitter::DoTrentoFit() +{ + fTrentoNBD->SetNpx(fFitNpx); + TFitResultPtr fitptr; + fFitOptions.Append("S"); + fitptr = fhV0M->Fit("fTrentoNBD", fFitOptions.Data()); + fMu = fTrentoNBD->GetParameter(Index(FitPar::mu)); + fk = fTrentoNBD->GetParameter(Index(FitPar::k)); + fnorm = fTrentoNBD->GetParameter(Index(FitPar::norm)); return fitptr.Get()->IsValid(); } //________________________________________________________________ -Bool_t multGlauberNBDFitter::InitializeNpNc() +bool multGlauberNBDFitter::InitializeNpNc() { // This function initializes fhNpNc // Warning: X == Npart, Y == Ncoll - Bool_t lReturnValue = kFALSE; + bool lReturnValue = false; if (fhNpNc) { fNNpNcPairs = 0; // Sweep all allowed values of Npart, Ncoll; find counters - for (int xbin = 1; xbin < 500; xbin++) { - for (int ybin = 1; ybin < 3000; ybin++) { + static constexpr int NBinsX = 500; + static constexpr int NBinsY = 3000; + for (int xbin = 1; xbin < NBinsX; xbin++) { + for (int ybin = 1; ybin < NBinsY; ybin++) { if (fhNpNc->GetBinContent(fhNpNc->FindBin(xbin, ybin)) != 0) { fNpart[fNNpNcPairs] = xbin; fNcoll[fNNpNcPairs] = ybin; - fContent[fNNpNcPairs] = fhNpNc->GetBinContent(fhNpNc->FindBin(xbin, ybin)); + fContent[fNNpNcPairs] = static_cast(fhNpNc->GetBinContent(fhNpNc->FindBin(xbin, ybin))); fNNpNcPairs++; } } } - cout << "Initialized with number of (Npart, Ncoll) pairs: " << fNNpNcPairs << endl; - lReturnValue = kTRUE; + LOG(info) << "Initialized with number of (Npart, Ncoll) pairs: " << fNNpNcPairs; + lReturnValue = true; } else { - cout << "Failed to initialize! Please provide input histogram with (Npart, Ncoll) info!" << endl; - cout << "Please remember to call SetNpartNcollCorrelation before doing fit!" << endl; + LOG(info) << "Failed to initialize! Please provide input histogram with (Npart, Ncoll) info!"; + LOG(info) << "Please remember to call SetNpartNcollCorrelation before doing fit!"; } return lReturnValue; } //________________________________________________________________ -Double_t multGlauberNBDFitter::ContinuousNBD(Double_t n, Double_t mu, Double_t k) +double multGlauberNBDFitter::ContinuousNBD(const double n, const double mu, const double k) { // Adaptation of the negative binomial distribution // for non-integer arguments: analytical continuation @@ -367,43 +424,44 @@ Double_t multGlauberNBDFitter::ContinuousNBD(Double_t n, Double_t mu, Double_t k // in fact it is equivalent to that if 'n' is typecast as // an integer prior to use - Double_t F; - Double_t f; + double F{}; + double f{}; - if (n + k > 100.0) { + static constexpr double NumStabilityThreshold = 100.0; + if (n + k > NumStabilityThreshold) { // 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 = std::lgamma(n + k) - std::lgamma(n + 1.) - std::lgamma(k); + f = n * std::log(mu / k) - (n + k) * std::log1p(mu / k); F = F + f; - F = TMath::Exp(F); + F = std::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 = std::tgamma(n + k) / (std::tgamma(n + 1.) * std::tgamma(k)); + f = n * std::log(mu / k) - (n + k) * std::log1p(mu / k); + f = std::exp(f); F *= f; } return F; } -void multGlauberNBDFitter::CalculateAvNpNc(TProfile* lNPartProf, TProfile* lNCollProf, TH2F* lNPart2DPlot, TH2F* lNColl2DPlot, TH1F* hPercentileMap, Double_t lLoRange, Double_t lHiRange, TH3D* lNpNcEcc, TH2F* lEcc2DPlot, TH3D* lNpNcB, TH2F* lB2DPlot, TH2F* lNancestor2DPlot, Double_t fProbabilityCutoff) +void multGlauberNBDFitter::CalculateAvNpNc(TProfile* lNPartProf, TProfile* lNCollProf, TH2F* lNPart2DPlot, TH2F* lNColl2DPlot, TH1F* hPercentileMap, double lLoRange, double lHiRange, TH3D* lNpNcEcc, TH2F* lEcc2DPlot, TH3D* lNpNcB, TH2F* lB2DPlot, TH2F* lNancestor2DPlot, double fProbabilityCutoff) { - cout << "Calculating , in centrality bins..." << endl; - cout << "Range to calculate: " << lLoRange << " to " << lHiRange << endl; + LOG(info) << "Calculating , in centrality bins..."; + LOG(info) << "Range to calculate: " << lLoRange << " to " << lHiRange; - cout << "Acquiring values from the fit function..." << endl; + LOG(info) << "Acquiring values from the fit function..."; - fMu = fGlauberNBD->GetParameter(0); - fk = fGlauberNBD->GetParameter(1); - ff = fGlauberNBD->GetParameter(2); - fnorm = fGlauberNBD->GetParameter(3); - fdMu = fGlauberNBD->GetParameter(4); + fMu = fGlauberNBD->GetParameter(Index(FitPar::mu)); + fk = fGlauberNBD->GetParameter(Index(FitPar::k)); + ff = fGlauberNBD->GetParameter(Index(FitPar::f)); + fnorm = fGlauberNBD->GetParameter(Index(FitPar::norm)); + fdMu = fGlauberNBD->GetParameter(Index(FitPar::dMu)); - cout << "Please inspect now: " << endl; - cout << "Glauber NBD mu ............: " << fMu << endl; - cout << "Glauber NBD k .............: " << fk << endl; - cout << "Glauber NBD f .............: " << ff << endl; - cout << "Glauber NBD norm ..........: " << fnorm << endl; - cout << "Glauber NBD dmu/dNanc .....: " << fdMu << endl; + LOG(info) << "Please inspect now: "; + LOG(info) << "Glauber NBD mu ............: " << fMu; + LOG(info) << "Glauber NBD k .............: " << fk; + LOG(info) << "Glauber NBD f .............: " << ff; + LOG(info) << "Glauber NBD norm ..........: " << fnorm; + LOG(info) << "Glauber NBD dmu/dNanc .....: " << fdMu; // 2-fold nested loop: // + looping over all Nancestor combinations @@ -416,31 +474,36 @@ void multGlauberNBDFitter::CalculateAvNpNc(TProfile* lNPartProf, TProfile* lNCol } // bypass to zero for (int ibin = 0; ibin < fNNpNcPairs; ibin++) { - if (ibin % 200 == 0) - cout << "At NpNc pair #" << ibin << " of " << fNNpNcPairs << "..." << endl; - Double_t lNAncestors0 = (Int_t)(fNpart[ibin] * ff + fNcoll[ibin] * (1.0 - ff)); - Double_t lNAncestors1 = TMath::Floor(fNpart[ibin] * ff + fNcoll[ibin] * (1.0 - ff) + 0.5); - Double_t lNAncestors2 = (fNpart[ibin] * ff + fNcoll[ibin] * (1.0 - ff)); + static constexpr int LogInterval = 200; + if (ibin % LogInterval == 0) { + LOG(info) << "At NpNc pair #" << ibin << " of " << fNNpNcPairs << "..."; + } + double lNAncestors0 = static_cast(fNpart[ibin] * ff + fNcoll[ibin] * (1.0 - ff)); + double lNAncestors1 = std::floor(fNpart[ibin] * ff + fNcoll[ibin] * (1.0 - ff) + 0.5); + double lNAncestors2 = (fNpart[ibin] * ff + fNcoll[ibin] * (1.0 - ff)); // define ancestors officially - Double_t lNancestors = lNAncestors0; - if (fAncestorMode == 1) + double lNancestors = lNAncestors0; + if (fAncestorMode == AncestorMode::Rounded) { lNancestors = lNAncestors1; - if (fAncestorMode == 2) + } + if (fAncestorMode == AncestorMode::Continuous) { lNancestors = lNAncestors2; + } // eccentricity handling - TH1D* hEccentricity = 0x0; + TH1D* hEccentricity = nullptr; if (lNpNcEcc) { // locate the histogram that corresponds to the eccentricity distribution in this NpNc pair lNpNcEcc->GetXaxis()->SetRange(lNpNcEcc->GetXaxis()->FindBin(fNpart[ibin]), lNpNcEcc->GetXaxis()->FindBin(fNpart[ibin])); lNpNcEcc->GetYaxis()->SetRange(lNpNcEcc->GetYaxis()->FindBin(fNcoll[ibin]), lNpNcEcc->GetYaxis()->FindBin(fNcoll[ibin])); - hEccentricity = reinterpret_cast(lNpNcEcc->Project3D("z")); + hEccentricity = dynamic_cast(lNpNcEcc->Project3D("z")); hEccentricity->SetName(Form("hEccentricity_%i", ibin)); // normalize into unitary fractions - Double_t eccIntegral = hEccentricity->Integral(1, hEccentricity->GetNbinsX() + 1); - if (eccIntegral > 1e-6) { // no counts + double eccIntegral = hEccentricity->Integral(1, hEccentricity->GetNbinsX() + 1); + static constexpr double EccTolerance = 1e-6; + if (eccIntegral > EccTolerance) { // no counts hEccentricity->Scale(1. / eccIntegral); } else { hEccentricity->Scale(0.0); @@ -448,52 +511,58 @@ void multGlauberNBDFitter::CalculateAvNpNc(TProfile* lNPartProf, TProfile* lNCol } // impact parameter handling - TH1D* hImpactParameter = 0x0; + TH1D* hImpactParameter = nullptr; if (lNpNcB) { // locate the histogram that corresponds to the eccentricity distribution in this NpNc pair lNpNcB->GetXaxis()->SetRange(lNpNcB->GetXaxis()->FindBin(fNpart[ibin]), lNpNcB->GetXaxis()->FindBin(fNpart[ibin])); lNpNcB->GetYaxis()->SetRange(lNpNcB->GetYaxis()->FindBin(fNcoll[ibin]), lNpNcB->GetYaxis()->FindBin(fNcoll[ibin])); - hImpactParameter = reinterpret_cast(lNpNcB->Project3D("z")); + hImpactParameter = dynamic_cast(lNpNcB->Project3D("z")); hImpactParameter->SetName(Form("hImpactParameter_%i", ibin)); // normalize into unitary fractions - Double_t bIntegral = hImpactParameter->Integral(1, hImpactParameter->GetNbinsX() + 1); - if (bIntegral > 1e-6) { // no counts + double bIntegral = hImpactParameter->Integral(1, hImpactParameter->GetNbinsX() + 1); + static constexpr double ImpactParTolerance = 1e-6; + if (bIntegral > ImpactParTolerance) { // no counts hImpactParameter->Scale(1. / bIntegral); } else { hImpactParameter->Scale(0.0); } } - for (Long_t lMultValue = 1; lMultValue < lHiRange; lMultValue++) { - Double_t lNancestorCount = fContent[ibin]; - Double_t lThisMu = (((Double_t)lNancestors)) * fMu; - Double_t lThisk = (((Double_t)lNancestors)) * fk; - Double_t lpval = TMath::Power(1 + lThisMu / lThisk, -1); + for (int64_t lMultValue = 1; lMultValue < lHiRange; lMultValue++) { + double lNancestorCount = fContent[ibin]; + double lThisMu = lNancestors * fMu; + double lThisk = lNancestors * fk; + double lpval = std::pow(1 + lThisMu / lThisk, -1); fNBD->SetParameter(1, lThisk); fNBD->SetParameter(0, lpval); - Double_t lMult = 0.0; - if (lMultValue > 1e-6) - lMult = fAncestorMode != 2 ? fNBD->Eval(lMultValue) : ContinuousNBD(lMultValue, lThisMu, lThisk); - Double_t lProbability = lNancestorCount * lMult; + double lMult = 0.0; + static constexpr double MultTolerance = 1e-6; + if (lMultValue > MultTolerance) { + lMult = fAncestorMode != AncestorMode::Continuous ? fNBD->Eval(lMultValue) : ContinuousNBD(lMultValue, lThisMu, lThisk); + } + double lProbability = lNancestorCount * lMult; if (lProbability < fProbabilityCutoff) { continue; // skip if probability of contributing too small } - Double_t lMultValueToFill = lMultValue; - if (hPercentileMap) + double lMultValueToFill = lMultValue; + if (hPercentileMap) { lMultValueToFill = hPercentileMap->GetBinContent(hPercentileMap->FindBin(lMultValue)); + } lNPartProf->Fill(lMultValueToFill, fNpart[ibin], lProbability); lNCollProf->Fill(lMultValueToFill, fNcoll[ibin], lProbability); if (lNancestor2DPlot) { // fill cross-check histogram with lNancestorCount at lNancestors value lNancestor2DPlot->Fill(lMultValueToFill, lNancestors, lProbability); } - if (lNPart2DPlot) + if (lNPart2DPlot) { lNPart2DPlot->Fill(lMultValueToFill, fNpart[ibin], lProbability); - if (lNColl2DPlot) + } + if (lNColl2DPlot) { lNColl2DPlot->Fill(lMultValueToFill, fNcoll[ibin], lProbability); + } if (lNpNcEcc) { // collapse the entire eccentricity distribution for this combo for (int ib = 1; ib < hEccentricity->GetNbinsX() + 1; ib++) { diff --git a/Common/Tools/Multiplicity/multGlauberNBDFitter.h b/Common/Tools/Multiplicity/multGlauberNBDFitter.h index 598ef312ce9..6c006be3b07 100644 --- a/Common/Tools/Multiplicity/multGlauberNBDFitter.h +++ b/Common/Tools/Multiplicity/multGlauberNBDFitter.h @@ -8,10 +8,16 @@ // In applying this license CERN does not waive the privileges and immunities // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -// + +/// \file multGlauberNBDFitter.h +/// \brief Helper class to perform Glauber fits +/// \author David Dobrigkeit Chinellato + #ifndef COMMON_TOOLS_MULTIPLICITY_MULTGLAUBERNBDFITTER_H_ #define COMMON_TOOLS_MULTIPLICITY_MULTGLAUBERNBDFITTER_H_ +#include + #include #include #include @@ -21,67 +27,109 @@ #include #include -#include class multGlauberNBDFitter : public TNamed { - public: // basic functionality multGlauberNBDFitter(); explicit multGlauberNBDFitter(const char* name, const char* title = "Glauber+NBD fitter"); - ~multGlauberNBDFitter(); + ~multGlauberNBDFitter() override; + + enum class AncestorMode { + Truncated = 0, // Nancestors = (int) truncation + Rounded, // Nancestors = rounded to nearest integer + Continuous // Nancestors = kept as float, uses ContinuousNBD + }; + + enum class NBDFitterMode { + Glauber = 0, // Uses Npart & Ncoll from TGlauberMC + Trento // Uses trento entropy + }; + + enum class FitPar { + mu = 0, + k, + f, + norm, + dMu + }; + + enum class NBDPar { + p = 0, + k + }; + + template + static constexpr int Index(TEnum e) + { + return static_cast(e); + } // Master fitter function - Double_t ProbDistrib(Double_t* x, Double_t* par); + double GlauberProbDistrib(const double* x, const double* par); + double TrentoProbDistrib(const double* x, const double* par); void InitAncestor(); // Do Fit: where everything happens - Bool_t DoFit(); + bool DoFit(); + bool DoGlauberFit(); + bool DoTrentoFit(); // Set input characteristics: the 2D plot with Npart, Nanc - Bool_t SetNpartNcollCorrelation(TH2* hNpNc); + bool SetNpartNcollCorrelation(TH2* hNpNc); + // Set input characteristics: the 1D plot with trento entropy + bool SetNSources(TH1* hNSources); // Set main input to be fitted (the V0M distribution) - Bool_t SetInputV0M(TH1* hV0M); + bool SetInputV0M(TH1* hV0M); // Interface to get funtions if asked to TF1* GetNBD(); TF1* GetGlauberNBD(); + TF1* GetTrentoNBD(); // Helper - Bool_t InitializeNpNc(); + bool InitializeNpNc(); + void InitGlauberNBD(const float mu, const float k, const float f, const float norm); + void InitTrentoNBD(const float mu = 45, const float k = 1.5, const float norm = 100); // Interface for debug - void SetAncestorMode(Int_t lAncMode = 0) { fAncestorMode = lAncMode; } - Int_t GetAncestorMode() { return fAncestorMode; } + void SetAncestorMode(AncestorMode lAncMode = AncestorMode::Truncated) { fAncestorMode = lAncMode; } + void SetAncestorMode(int lAncMode = Index(AncestorMode::Truncated)) + { + if (lAncMode < Index(AncestorMode::Truncated) || lAncMode > Index(AncestorMode::Continuous)) { + LOG(fatal) << "Invalid ancestor mode: " << lAncMode; + return; + } + fAncestorMode = static_cast(lAncMode); + } + int GetAncestorMode() { return static_cast(fAncestorMode); } TH1D* GetAncestorHistogram() { return fhNanc; } // Interface to set vals - void SetMu(Double_t lVal) { fMu = lVal; } - void Setk(Double_t lVal) { fk = lVal; } - void Setf(Double_t lVal) { ff = lVal; } - void SetNorm(Double_t lVal) { fnorm = lVal; } + void SetMu(const double lVal) { fMu = lVal; } + void Setk(const double lVal) { fk = lVal; } + void Setf(const double lVal) { ff = lVal; } + void SetNorm(const double lVal) { fnorm = lVal; } // Interface to get vals - Double_t GetMu() { return fMu; } - Double_t Getk() { return fk; } - Double_t Getf() { return ff; } - Double_t GetNorm() { return fnorm; } + double GetMu() { return fMu; } + double Getk() { return fk; } + double Getf() { return ff; } + double GetNorm() { return fnorm; } - void SetFitRange(Double_t lMin, Double_t lMax); - void SetFitOptions(TString lOpt); - void SetFitNpx(Long_t lNpx); + void SetFitRange(const double lMin, const double lMax); + void SetFitOptions(const TString& lOpt); + void SetFitNpx(const int64_t lNpx); // For ancestor mode 2 - Double_t ContinuousNBD(Double_t n, Double_t mu, Double_t k); + double ContinuousNBD(const double n, const double mu, const double k); // For estimating Npart, Ncoll in multiplicity bins // also viable: eccentricity, impact parameter, ancestor cross-check plot - void CalculateAvNpNc(TProfile* lNPartProf, TProfile* lNCollProf, TH2F* lNPart2DPlot, TH2F* lNColl2DPlot, TH1F* hPercentileMap, Double_t lLoRange = -1, Double_t lHiRange = -1, TH3D* lNpNcEcc = 0x0, TH2F* lEcc2DPlot = 0x0, TH3D* lNpNcB = 0x0, TH2F* lB2DPlot = 0x0, TH2F* lNancestor2DPlot = 0x0, Double_t fProbabilityCutoff = -1); - - // void Print(Option_t *option="") const; + void CalculateAvNpNc(TProfile* lNPartProf, TProfile* lNCollProf, TH2F* lNPart2DPlot, TH2F* lNColl2DPlot, TH1F* hPercentileMap, double lLoRange = -1, double lHiRange = -1, TH3D* lNpNcEcc = nullptr, TH2F* lEcc2DPlot = nullptr, TH3D* lNpNcB = nullptr, TH2F* lB2DPlot = nullptr, TH2F* lNancestor2DPlot = nullptr, double fProbabilityCutoff = -1); private: // This function serves as the (analytical) NBD @@ -89,36 +137,41 @@ class multGlauberNBDFitter : public TNamed // This function is the key fitting function TF1* fGlauberNBD; + TF1* fTrentoNBD; + + // This is the fitting mode that we use + NBDFitterMode fNBDFitterMode; // Reference histo - TH1D* fhNanc; // basic ancestor distribution - TH2* fhNpNc; // correlation between Npart and Ncoll - TH1* fhV0M; // basic ancestor distribution + TH1D* fhNanc; // basic ancestor distribution + TH2* fhNpNc; // correlation between Npart and Ncoll + TH1* fhNSources; // Trento entropy + TH1* fhV0M; // basic ancestor distribution // Fitting utilities - Bool_t ffChanged; - Double_t fCurrentf; + bool ffChanged; + double fCurrentf; // 0: truncation, 1: rounding, 2: analytical continuation - Int_t fAncestorMode; + AncestorMode fAncestorMode; // Buffer for (Npart, Ncoll) pairs in memory - Double_t* fNpart; - Double_t* fNcoll; - Long_t* fContent; - Long_t fNNpNcPairs; // number of pairs to use - Long_t fMaxNpNcPairs; + double* fNpart; + double* fNcoll; + int64_t* fContent; + int64_t fNNpNcPairs; // number of pairs to use + int64_t fMaxNpNcPairs; // The actual output: mu, k, f, norm - Double_t fMu; - Double_t fdMu; // variable mu option - Double_t fk; - Double_t ff; - Double_t fnorm; + double fMu; + double fdMu; // variable mu option + double fk; + double ff; + double fnorm; TString fFitOptions; - Long_t fFitNpx; + int64_t fFitNpx; - ClassDef(multGlauberNBDFitter, 1); + ClassDefOverride(multGlauberNBDFitter, 1); }; #endif // COMMON_TOOLS_MULTIPLICITY_MULTGLAUBERNBDFITTER_H_