KVMultiGaussIsotopeFit.cpp

#include "KVMultiGaussIsotopeFit.h"
#include "TGraph.h"
#include "TRandom.h"
#include <TCanvas.h>
 
 
 
double KVMultiGaussIsotopeFit::centroid_fit(double* x, double* p)
{
   // centroids of gaussians are expected to increase linearly with mass A but we allow a quadratic dependence
   //
   //~~~
   // x[0]=A
   // p[0]=offset
   // p[1]=slope
   // p[2]=quadratic term
   //~~~
 
   return p[0] + p[1] * x[0] + p[2] * x[0] * x[0];
}
 
 
 
 
double KVMultiGaussIsotopeFit::FitFunc(double* x, double* p)
{
   //~~~
   // x[0] = PID
   // p[0] = number of gaussians = Ng (fixed)
   // p[1],p[2] = background parameters: exp(p[1]+p[2]*x)
   // p[3]  = sigma for all gaussians
   // p[4],p[5],p[6] = offset & slope for centroid PID vs. A correlation
   // p[7]...p[6+Ng] = norm of each gaussian
   // p[6+Ng+1]...p[6+2*Ng] = A of each gaussian
   //~~~
   //
   // Total number of parameters is 7+2*Ng
 
   int Ng = p[0];
   double background = TMath::Exp(p[fit_param_index::bkg_cst] + p[fit_param_index::bkg_slp] * x[0]);
   for (int i = 1; i <= Ng; ++i) {
      background +=
     p[get_gauss_norm_index(i)] * TMath::Gaus(x[0], centroid_fit(&p[get_mass_index(i, Ng)],
           &p[fit_param_index::pidvsA_a0]), p[fit_param_index::gauss_wid], kTRUE);
   }
   return background;
}
 
   
 
 
std::optional<double> KVMultiGaussIsotopeFit::get_total_fit_for_pid(double PID) const
{
   // Calculates the total value (gaussians + background) of the fit for a given PID value.
   //
   // If this is less than 1, returns std::nullopt: in this case, the fit should not be used to
   // deduce an \f$A\f$ etc. for this PID.
   //
   // \note The identification functional is fitted to histogrammed data, therefore the resulting
   // complete function should always evaluate to some finite value \f$\geq 1\f$. If this is not the
   // case for certain values of PID, it suggests that the 'background' is not well fitted and the
   // functional cannot be used for identification (gives garbage results). In this case we do not
   // return a value.
   auto total = Eval(PID);
   if(total<1)
      return {};
   return total;
}
 
 
 
 
KVMultiGaussIsotopeFit::KVMultiGaussIsotopeFit(int z, int Ngauss, double PID_min, double PID_max, const KVNumberList& alist, std::vector<double> pidlist)
   : TF1(Form("MultiGaussIsotopeFit_Z=%d", z), this, &KVMultiGaussIsotopeFit::FitFunc, PID_min, PID_max, total_number_parameters(Ngauss)),
     Z{z},
     Niso{Ngauss},
     PIDmin{PID_min}, PIDmax{PID_max},
     Alist{alist.GetArray()},
     PIDlist{pidlist}
{
   // Constructor used to initialize and prepare a new fit of isotope PID spectrum
   //
   // \param[in] z atomic number \f$Z\f$ of element to fit
   // \param[in] Ngauss number of gaussians to use in the fit (corresponds to number of isotopes in PID spectrum)
   // \param[in] PID_min,PID_max limits of PID spectrum to which fit is performed
   // \param[in] alist ordered list of increasing \f$A\f$-values of isotopes in PID spectrum 
   // \param[in] pidlist corresponding initial values of PID for centroid of each isotope-gaussian in PID spectrum
    
   FixParameter(0, Niso);
   SetParLimits(fit_param_index::bkg_cst, -10., 25.);
   SetParameter(fit_param_index::bkg_cst, 4.);
   SetParName(fit_param_index::bkg_cst, "Norm");
   PositiveBkgSlope(false);
   SetParName(fit_param_index::bkg_slp, "Bkg. slope");
   SetParName(fit_param_index::gauss_wid, "Sigma");
   SetParLimits(fit_param_index::gauss_wid, min_sigma, max_sigma);
   SetParameter(fit_param_index::gauss_wid, 0.1);
 
   for (int ig = 1; ig <= Niso; ++ig) {
      SetParName(get_gauss_norm_index(ig), Form("Norm. A=%d", Alist[ig - 1]));
      SetParLimits(get_gauss_norm_index(ig), 1.e-3, 1.e+06);
      SetParameter(get_gauss_norm_index(ig), 1.);
      SetParName(get_mass_index(ig, Niso), Form("A_%d", ig));
      FixParameter(get_mass_index(ig, Niso), Alist[ig - 1]);
   }
 
   FitCentroids();
 
   SetLineColor(kBlack);
   SetLineWidth(2);
   SetNpx(500);
}
 
 
 
 
void KVMultiGaussIsotopeFit::InitializeParameterLimitsForNewFit()
{
   SetParLimits(fit_param_index::bkg_cst, -10., 25.);
   PositiveBkgSlope(GetParameter(fit_param_index::bkg_slp)>0);
   SetParLimits(fit_param_index::gauss_wid, min_sigma, max_sigma);
   for (int ig = 1; ig <= Niso; ++ig)
   {
      SetParLimits(get_gauss_norm_index(ig), 1.e-3, 1.e+06);
   }
}
 
 
 
 
KVMultiGaussIsotopeFit::KVMultiGaussIsotopeFit(int z, int Ngauss, double PID_min, double PID_max,
      const KVNumberList& alist, double bkg_cst, double bkg_slp,
      double gaus_wid, double pidvsa_a0, double pidvsa_a1, double pidvsa_a2)
   : TF1(Form("MultiGaussIsotopeFit_Z=%d", z), this, &KVMultiGaussIsotopeFit::FitFunc, PID_min, PID_max, total_number_parameters(Ngauss)),
     Z{z},
     Niso{Ngauss},
     PIDmin{PID_min}, PIDmax{PID_max},
     Alist{alist.GetArray()}
{
   // Constructor which can be used with existing fit results (not to perform new fits)
   //
   // Use SetGaussianNorm() to set the normalisation parameters for each gaussian
   //
   // \param[in] z atomic number \f$Z\f$ of element in fit
   // \param[in] Ngauss number of gaussians used in the fit (corresponds to number of isotopes in PID spectrum)
   // \param[in] PID_min,PID_max limits of PID spectrum to which fit was performed
   // \param[in] alist ordered list of increasing \f$A\f$-values of isotopes in PID spectrum 
   // \param[in] bkg_cst,bkg_slp constant and slope of fitted exponential background function
   // \param[in] gaus_wid width of gaussians used in fit
   // \param[in] pidvsa_a0,pidvsa_a1,pidvsa_a2 parameters of polynomial fit to PID-coordinate \f$A\f$-dependence
 
   FixParameter(0, Niso);
   SetParameter(fit_param_index::bkg_cst, bkg_cst);
   SetParameter(fit_param_index::bkg_slp, bkg_slp);
   SetParameter(fit_param_index::gauss_wid, gaus_wid);
   SetParameter(fit_param_index::pidvsA_a0, pidvsa_a0);
   SetParameter(fit_param_index::pidvsA_a1, pidvsa_a1);
   SetParameter(fit_param_index::pidvsA_a2, pidvsa_a2);
 
   SetParName(fit_param_index::bkg_cst, "Norm");
   SetParName(fit_param_index::bkg_slp, "Bkg. slope");
   SetParName(fit_param_index::gauss_wid, "Sigma");
 
   UpdateGaussianCentroidParameters();
 
   for (int i = 0; i < 3; ++i)
      SetParName(fit_param_index::pidvsA_a0 + i, Form("PIDvsA_a%d", i));
 
   SetLineColor(kBlack);
   SetLineWidth(2);
   SetNpx(500);
}
 
 
 
 
KVMultiGaussIsotopeFit::KVMultiGaussIsotopeFit(int Z, const KVNameValueList& fitparams)
   : KVMultiGaussIsotopeFit(Z, fitparams.GetIntValue("Ng"), fitparams.GetDoubleValue("PIDmin"),
                            fitparams.GetDoubleValue("PIDmax"), fitparams.GetStringValue("Alist"),
                            fitparams.GetDoubleValue("Bkg_cst"), fitparams.GetDoubleValue("Bkg_slp"),
                            fitparams.GetDoubleValue("GausWid"),
                            fitparams.GetDoubleValue("PIDvsA_a0"),
                            fitparams.GetDoubleValue("PIDvsA_a1"),
                            fitparams.GetDoubleValue("PIDvsA_a2"))
{
   // initialize from previous fit with parameters stored in KVNameValueList
   for (int ig = 1; ig <= fitparams.GetIntValue("Ng"); ++ig)
      SetGaussianNorm(ig, fitparams.GetDoubleValue(Form("Norm_%d", ig)));
}
 
 
 
 
void KVMultiGaussIsotopeFit::UpdateGaussianCentroidParameters()
{
   // Required in case the attribution of masses to the gaussians changes
   for (int ig = 1; ig <= Niso; ++ig) {
      SetParName(get_gauss_norm_index(ig), Form("Norm. A=%d", Alist[ig - 1]));
      SetParName(get_mass_index(ig, Niso), Form("A_%d", ig));
      FixParameter(get_mass_index(ig, Niso), Alist[ig - 1]);
   }
}
 
 
 
 
void KVMultiGaussIsotopeFit::FitCentroids()
{
   TGraph pid_vs_a;
   for (int ig = 1; ig <= Niso; ++ig) {
#if ROOT_VERSION_CODE >= ROOT_VERSION(6,24,0)
      pid_vs_a.AddPoint(Alist[ig - 1], PIDlist[ig - 1]);
#else
      pid_vs_a.SetPoint(pid_vs_a.GetN(), Alist[ig - 1], PIDlist[ig - 1]);
#endif
   }
   TF1 centroidFit("centroidFit", this, &KVMultiGaussIsotopeFit::centroid_fit, 0., 100., 3);
   centroidFit.SetParLimits(0, -50, 50);
   centroidFit.SetParameter(0, 0);
   centroidFit.SetParLimits(1, 1.e-2, 5.);
   centroidFit.SetParameter(1, 1.e-1);
   centroidFit.SetParLimits(2, -2.e-2, 1.);
   centroidFit.SetParameter(2, 1.e-3);
   pid_vs_a.Fit(&centroidFit, "N");
   for (int i = 0; i < 3; ++i) {
      FixParameter(fit_param_index::pidvsA_a0 + i, centroidFit.GetParameter(i));
      SetParName(fit_param_index::pidvsA_a0 + i, Form("PIDvsA_a%d", i));
   }
}
 
 
 
 
void KVMultiGaussIsotopeFit::ReleaseCentroids()
{
   // Release the constraint on the positions of the centroids
   SetParLimits(fit_param_index::pidvsA_a0, -50, 50);
   SetParLimits(fit_param_index::pidvsA_a1, 1.e-2, 5.);
   SetParLimits(fit_param_index::pidvsA_a2, -2.e-2, 5.);
}
 
 
 
 
void KVMultiGaussIsotopeFit::UnDraw(TVirtualPad* pad) const
{
   // Remove the graphical representation of this fit from the given pad
 
   auto old_fit = pad->FindObject(get_name_of_multifit(Z));
   if (old_fit) delete old_fit;
   for (auto a : Alist) {
      UnDrawGaussian(Z, a, pad);
   }
}
 
 
 
 
void KVMultiGaussIsotopeFit::UnDrawAnyGaussian(int z, TVirtualPad* pad)
{
   // Remove the graphical representation of any gaussian for this Z from the given pad
 
   TIter it(pad->GetListOfPrimitives());
   TObject* ob;
   KVList to_delete;
   while ((ob = it())) {
      TString obname = ob->GetName();
      if (obname.BeginsWith(get_root_name_of_isotope_gaussian(z))) to_delete.Add(ob);
   }
}
 
 
 
 
void KVMultiGaussIsotopeFit::DrawFitWithGaussians(Option_t* opt, const TString& fit_title) const
{
   // Draw the overall fit plus the individual gaussians for each isotope
   //
   // \param[in] opt drawing option, will be used for main fit (could be `"same"`)
   // \param[in] fit_title optional, if given will be used as title of fit function (displayed in TLegend and also above frame in TPad with option 'show histogram title')
 
   auto ff = DrawCopy(opt);
   ff->SetName(get_name_of_multifit(Z));
   if(!fit_title.IsNull())
      ff->SetTitle(fit_title);
   TF1 fgaus("fgaus", "gausn", PIDmin, PIDmax);
   // give a different colour to each gaussian
   auto cstep = TColor::GetPalette().GetSize() / (Niso + 1);
   int ig = 1;
   for (auto& a : Alist) {
      fgaus.SetParameters(GetGaussianNorm(ig), GetCentroid(ig), GetGaussianWidth(ig));
      fgaus.SetNpx(500);
      fgaus.SetLineColor(TColor::GetPalette()[cstep * ig]);
      fgaus.SetLineWidth(2);
      fgaus.SetLineStyle(9);
      auto gg = fgaus.DrawCopy("same");
      gg->SetName(get_name_of_isotope_gaussian(Z, a));
      gg->SetTitle(get_name_of_isotope_gaussian(Z, a));
      ++ig;
   }
}
 
 
 
 
std::optional<int> KVMultiGaussIsotopeFit::GetMostProbableA(double PID, double& P) const
{
   // For a given PID, calculate the most probable value of \f$A\f$, P is its probability.
   //
   // \returns most probable \f$A\f$
   //
   // \note The identification functional is fitted to histogrammed data, therefore the resulting
   // complete function should always evaluate to some finite value \f$\geq 1\f$. If this is not the
   // case for certain values of PID, it suggests that the 'background' is not well fitted and the
   // functional cannot be used for identification (gives garbage results). In this case we do not
   // return a value.
 
   auto total = get_total_fit_for_pid(PID);
   if(!total)
      return {};
   std::map<double, int> probabilities;
   int ig = 1;
   for (auto& a : Alist) {
      probabilities[evaluate_gaussian(ig, PID) / *total] = a;
      ++ig;
   }
   // the largest probability is now the last element in the map:
   // get a reverse iterator to the beginning of the reversed map
   auto it = probabilities.rbegin();
   P = it->first;
   return it->second;
}
 
 
 
 
std::optional<double> KVMultiGaussIsotopeFit::GetMeanA(double PID) const
{
   // for a given PID, calculate the mean value of \f$A\f$ from the weighted sum of all gaussians
   //
   // \returns mean value of \f$A\f$
   //
   // \note The identification functional is fitted to histogrammed data, therefore the resulting
   // complete function should always evaluate to some finite value \f$\geq 1\f$. If this is not the
   // case for certain values of PID, it suggests that the 'background' is not well fitted and the
   // functional cannot be used for identification (gives garbage results). In this case we do not
   // return a value.
 
   auto total = get_total_fit_for_pid(PID);
   if(!total)
      return {};
   int ig = 1;
   double amean(0), totprob(0);
   for (auto& a : Alist) {
      auto weight = evaluate_gaussian(ig, PID) / *total;
      amean += weight * a;
      totprob += weight;
      ++ig;
   }
   return totprob > 0 ? amean / totprob : -1.;
}
 
 
 
 
std::optional<std::map<int, double>> KVMultiGaussIsotopeFit::GetADistribution(double PID) const
{
   // For the given PID, the map is filled with all possible values of \f$A\f$
   //
   // \note The identification functional is fitted to histogrammed data, therefore the resulting
   // complete function should always evaluate to some finite value \f$\geq 1\f$. If this is not the
   // case for certain values of PID, it suggests that the 'background' is not well fitted and the
   // functional cannot be used for identification (gives garbage results). In this case we do not
   // return a value.
   //
   // \returns std::map containing probability distribution \f$P(A|PID)\f$
 
   auto total = get_total_fit_for_pid(PID);
   if(!total)
      return {};
   std::map<int, double> Adist;
   int ig = 1;
   for (auto& a : Alist) {
      Adist[a] = evaluate_gaussian(ig, PID) / *total;
      ++ig;
   }
   return Adist;
}
 
 
 
 
std::optional<int> KVMultiGaussIsotopeFit::GetA(double PID, double& P) const
{
   // Probabilistic method to determine \f$A\f$ from PID.
   //
   // The A returned will be drawn at random from the probability distribution given by the
   // sum of all gaussians (and the background) for the given PID.
   //
   // The result of the draw may be that this PID is part of the background noise:
   // in this case we do not return a value (returns std::nullopt : returned value will evaluate as false).
   //
   // P is the probability of the chosen result.
   //
   // \param[in] PID PID value from \f$Z\f$ identification
   // \param[out] P probablility the returned value of \f$A\f$ is correct
   // \returns randomly drawn \f$A\f$ for given PID
   //
   // \note The identification functional is fitted to histogrammed data, therefore the resulting
   // complete function should always evaluate to some finite value \f$\geq 1\f$. If this is not the
   // case for certain values of PID, it suggests that the 'background' is not well fitted and the
   // functional cannot be used for identification (gives garbage results). In this case we do not
   // return a value.
 
   auto total = get_total_fit_for_pid(PID);
   if(!total)
      return {};
   double p_tot = 0;
   int ig = 1;
   auto X = gRandom->Uniform();
   for (auto& a : Alist) {
      auto w = evaluate_gaussian(ig, PID) / *total;
      p_tot += w;
      if (X < w) {
         P = w;
         return a;
      }
      X -= w;
      ++ig;
   }
   P = 1. - p_tot;
   return {}; // background noise
}
 
 
 
 
double KVMultiGaussIsotopeFit::GetProbability(int A, double PID) const
{
   // \param[in] A isotope mass number
   // \param[in] PID value of PID associated with A
   // \return the probability that A is the correct mass number for a given PID value
   // \return zero if A is not associated with a Gaussian in the fit
 
   auto it = std::find(std::begin(Alist), std::end(Alist), A);
   if (it != std::end(Alist)) {
      int ig = std::distance(std::begin(Alist), it);
      return evaluate_gaussian(ig, PID) / Eval(PID);
   }
   return 0;
}
 
 
 
 
double KVMultiGaussIsotopeFit::GetInterpolatedA(double PID) const
{
   // \returns interpolated A from PID using the quadratic fit parameters
   // analytic inversion of \f$f(A_i) = a_0 + a_1 A_i + a_2 A_i^2\f$
 
   auto a = GetPIDvsAfit_a2();
   auto b = GetPIDvsAfit_a1();
   auto c = GetPIDvsAfit_a0() - PID;
   return (TMath::Sqrt(b * b - 4.*a * c) - b) / 2. / a;
}
 
 
 
 
double KVMultiGaussIsotopeFit::GetCentroid(int i) const
{
   // \returns the fitted centroid position of the ith gaussian (i=1,2,...,Niso) i.e. \f$f(A_i) = a_0 + a_1 A_i + a_2 A_i^2\f$
   assert(i > 0 && i <= Niso);
   return GetParameter(fit_param_index::pidvsA_a0)
         + (GetParameter(fit_param_index::pidvsA_a1)
            + GetParameter(fit_param_index::pidvsA_a2) * Alist[i - 1]) * Alist[i - 1];
}
 
 
 
 
double KVMultiGaussIsotopeFit::GetPIDFromInterpolatedA(double interpA)
{
   // The PID value stored in the KVIdentificationResult corresponding to this identification is the interpolated
   // value of \f$A\f$ calculated using GetInterpolatedA(). Calling this method with that value returns the
   // actual PID value which can be compared with the gaussians of the fit.
   return GetParameter(fit_param_index::pidvsA_a0)
         + (GetParameter(fit_param_index::pidvsA_a1)
            + GetParameter(fit_param_index::pidvsA_a2) * interpA) * interpA;
}
 
 
 
 
 
 
void KVMultiGaussIsotopeFit::SetFitRange(double min, double max)
{
   // Change range of fit
 
   SetRange(min, max);
   PIDmin = min;
   PIDmax = max;
}
 
 
 
 
void KVMultiGaussIsotopeFit::PositiveBkgSlope(bool yes)
{
   // if yes=true, only explore positive values for background slope
   //
   // if yes=false, only explore negative values
   if(yes)
   {
      SetParLimits(fit_param_index::bkg_slp, 0., 5.);
      SetParameter(fit_param_index::bkg_slp, 1.);
   }
   else
   {
      SetParLimits(fit_param_index::bkg_slp, -10, 0);
      SetParameter(fit_param_index::bkg_slp, -2.);
   }
}
 
 
ClassImp(KVMultiGaussIsotopeFit)