KVGaussFitMassModifier.cpp

#include "KVGaussFitMassModifier.h"
#include "KVMultiGaussIsotopeFit.h"
 
 
 
KVGaussFitMassModifier::KVGaussFitMassModifier(KVIDZAFromZGrid* gr, int Z)
   : fGrid{gr}
{
   // Prepare mass modification object for the multi-gaussian fit for element \f$Z\f$ in the given grid
   fGrid->Initialize();
   fItvs = gr->GetIntervalSet(Z);
   fGfit = gr->GetMultiGaussFit(Z);
}
 
 
 
 
void KVGaussFitMassModifier::Modify(int deltaA, TH1* pid_dist, const KVNameValueList& mass_fit_parameters)
{
   // \param[in] deltaA mass modification \f$\Delta A\f$ to apply to all isotopes of element, can be positive or negative
   // \param[in] pid_dist a PID distribution corresponding to the grid being used
   // \param[in] mass_fit_parameters parameters to control fit; by default KVMultiGaussIsotopeFit::mass_fit_parameters
   //
   // \note For any KVIDZAGrid, a 1-dimensional PID distribution to be used with multi-gaussian fits
   // can be generated from a 2D identification map (e.g. \f$\Delta E-E\f$ matrix)
   // using KVIDZAFromZGrid::LinearizeHistoToPID()
 
   KVNumberList alist;
   std::vector<double> pidlist;
   TIter nxt_int(fItvs->GetIntervals());
   interval* intvl = nullptr;
   while ((intvl = (interval*)nxt_int())) {
      intvl->SetA(intvl->GetA()+deltaA);   // retrieve & modify the A of each interval set
      alist.Add(intvl->GetA());
      pidlist.push_back(intvl->GetPID());
   }
   // set PIDlist & modified Alist in the fit object
   fGfit->SetPIDlist(pidlist);
   fGfit->SetAlist(alist.GetArray());
   fGfit->UpdateGaussianCentroidParameters();
   fGfit->InitializeParameterLimitsForNewFit();
 
   bool fit_limited = false;
   // check user fit parameters
   if (mass_fit_parameters.GetBoolValue("Limit range of fit")) {
      // if the user's range is not valid for the current interval set, we ignore it
      if (mass_fit_parameters.GetDoubleValue("PID min for fit") >= fItvs->GetZ() - 0.5
            && mass_fit_parameters.GetDoubleValue("PID max for fit") <= fItvs->GetZ() + 0.5)
      {
         fGfit->SetFitRange(mass_fit_parameters.GetDoubleValue("PID min for fit"),
                             mass_fit_parameters.GetDoubleValue("PID max for fit"));
         fit_limited = true;
      }
   }
   fGfit->SetSigmaLimits(mass_fit_parameters.GetDoubleValue("Minimum #sigma"),
                          mass_fit_parameters.GetDoubleValue("Maximum #sigma"));
 
   // refit the PID-A relationship
   fGfit->FitCentroids();
 
   pid_dist->Fit(fGfit, "NR");
 
   // now release the centroids
   fGfit->ReleaseCentroids();
 
   pid_dist->Fit(fGfit, "NRME");
 
   // if fit was limited to a range, remove the limits
   if(fit_limited)
      fGfit->SetFitRange(fItvs->GetZ() - 0.5,fItvs->GetZ() + 0.5);
 
   // set interval limits according to regions of most probable mass
   int most_prob_A = 0;
   nxt_int.Reset();
   intvl = nullptr;
   // minimum probability for which isotopes are taken into account
   double min_proba = mass_fit_parameters.GetDoubleValue("Minimum probability [%]") / 100.;
   double delta_pid = 0.001;
   TList accepted_intervals;//any intervals not in this list at the end of the procedure will be removed
   for (double pid = fGfit->GetPIDmin() ; pid <= fGfit->GetPIDmax(); pid += delta_pid) {
      double proba;
      auto Amax = fGfit->GetMostProbableA(pid, proba);
      if(!Amax)
         continue;
      if (proba > min_proba) {
         if (most_prob_A) {
            if (*Amax > most_prob_A) {
               intvl->SetPIDmax(pid - delta_pid);
               most_prob_A = *Amax;
               intvl = (interval*)nxt_int();
               while (intvl->GetA() < most_prob_A) {
                  intvl = (interval*)nxt_int();
               }
               accepted_intervals.Add(intvl);
               intvl->SetPIDmin(pid);
            }
         }
         else {
            most_prob_A = *Amax;
            intvl = (interval*)nxt_int();
            while (intvl->GetA() < most_prob_A) {
               intvl = (interval*)nxt_int();
            }
            accepted_intervals.Add(intvl);
            intvl->SetPIDmin(pid);
         }
      }
      else if (most_prob_A) {
         intvl->SetPIDmax(pid - delta_pid);
         most_prob_A = 0;
      }
   }
   nxt_int.Reset();
   int ig(1);
   // update PID positions from fitted centroids
   nxt_int.Reset();
   KVNumberList remaining_gaussians, remaining_alist;
   auto vec_alist = alist.GetArray();
   while ((intvl = (interval*)nxt_int())) {
      // in case we removed some peaks
      while (vec_alist[ig - 1] < intvl->GetA()) {
         ++ig;
      }
      intvl->SetPID(fGfit->GetCentroid(ig));
      remaining_gaussians.Add(ig);
      remaining_alist.Add(intvl->GetA());
      ++ig;
   }
   // save results in grid parameters
   KVNumberList zlist;
   if (fGrid->GetParameters()->HasStringParameter("MASSFITS"))
      zlist.Set(fGrid->GetParameters()->GetStringValue("MASSFITS"));
   zlist.Add(fItvs->GetZ());
   fGrid->GetParameters()->SetValue("MASSFITS", zlist.AsString());
   TString massfit = Form("MASSFIT_%d", fItvs->GetZ());
   KVNameValueList fitparams;
   fitparams.SetValue("Ng", alist.GetNValues());
   fitparams.SetValue("Alist", alist.AsQuotedString());
   fitparams.SetValue("PIDmin", fGfit->GetPIDmin());
   fitparams.SetValue("PIDmax", fGfit->GetPIDmax());
   fitparams.SetValue("Bkg_cst", fGfit->GetBackgroundConstant());
   fitparams.SetValue("Bkg_slp", fGfit->GetBackgroundSlope());
   fitparams.SetValue("GausWid", fGfit->GetGaussianWidth(0));
   fitparams.SetValue("PIDvsA_a0", fGfit->GetPIDvsAfit_a0());
   fitparams.SetValue("PIDvsA_a1", fGfit->GetPIDvsAfit_a1());
   fitparams.SetValue("PIDvsA_a2", fGfit->GetPIDvsAfit_a2());
   for (ig = 1; ig <= alist.GetNValues(); ++ig) {
      fitparams.SetValue(Form("Norm_%d", ig), fGfit->GetGaussianNorm(ig));
   }
   auto sanitized = fitparams.Get().ReplaceAll("=", ":");
   fGrid->GetParameters()->SetValue(massfit, sanitized);
}
 
 
ClassImp(KVGaussFitMassModifier)