KVIDZAFromZGridMassCorrector.cpp

#include "KVIDZAFromZGridMassCorrector.h"
#include "KVGaussFitMassModifier.h"
#include "KVMultiGaussIsotopeFit.h"
 
int KVIDZAFromZGridMassCorrector::fOnlyShowZ = 0;
 
 
 
std::vector<KVIDZAFromZGridMassCorrector::mass_correction_t>
KVIDZAFromZGridMassCorrector::DetermineMassCorrections(double _min_yield)
{
   // Determine the mass-shifts to apply to each fitted element \f$Z\f$ by comparing the mass distributions with those of the reference fits.
   //
   // \param[in] _min_yield minimum yield (expressed as a fraction of the largest yield in the reference fit) above which to consider isotopes in mass distributions
   //
   // \note For \f$Z=6\f$ `_min_yield` is fixed at 0.5
   //
   // Up to 3 isotope masses will be compared for the 2 fits of each element, depending on the threshold minimum yield parameter.
   //
   // If the sequence of masses for the reference & other fit can be related with a single offset (i.e. masses "24,25,26" for one,
   // masses "23,24,25" for the other implies an offset of \f$\pm 1\f$), then this is taken to be the required mass correction,
   // unless such a correction would lead to assign a spurious value corresponding to an undetectable resonance etc.
   //
   // \note set static member variable KVIDZAFromZGridMassCorrector::fOnlyShowZ to a non-zero value in order to only
   // consider the given \f$Z\f$. Reset KVIDZAFromZGridMassCorrector::fOnlyShowZ=0 in order to treat all elements.
   //
   // \returns a vector containing the mass correction required for each \f$Z\f$
 
   refGrid->Initialize();
   modGrid->Initialize();
   if(!modGrid->GetFits() || modGrid->GetFits()->IsEmpty())
   {
      KVError::Warning(this,"DetermineMassCorrections", "Grid %s has no multi-gaussian mass fits", modGrid->GetName());
      return {};
   }
   if(!refGrid->GetFits() || refGrid->GetFits()->IsEmpty())
   {
      KVError::Warning(this,"DetermineMassCorrections", "Grid %s has no multi-gaussian mass fits", refGrid->GetName());
      return {};
   }
   TIter next_fit(refGrid->GetFits());
   KVMultiGaussIsotopeFit* ref_fit;
   std::vector<mass_correction_t> corrections;
   while( (ref_fit = (KVMultiGaussIsotopeFit*)next_fit()) )
   {
      auto mod_fit = modGrid->GetMultiGaussFit(ref_fit->GetZ());
      if(!mod_fit) continue;
 
      double min_yield=_min_yield;
      if(ref_fit->GetZ()==6)
      {
         //KVError::Info(this,"DetermineMassCorrections", "Setting threshold to 50% for carbon");
         min_yield=0.5;
      }
 
      KVNumberList ref_masses,mod_masses;
      int ngauss = std::min(std::min(ref_fit->GetNGaussians(),mod_fit->GetNGaussians()),3);
      auto ref_A = ref_fit->GetIsotopesRankedByYield();
      auto mod_A = mod_fit->GetIsotopesRankedByYield();
      auto ref_yield = ref_fit->GetRankedYields();
      auto mod_yield = mod_fit->GetRankedYields();
      for(int i=0; i<ngauss; ++i)
      {
         // only consider the most-produced isotopes (maybe 1, 2 or 3)
         if(ref_yield[i]>min_yield*ref_yield[0])
         {
            ref_masses.Add(ref_A[i]);
            mod_masses.Add(mod_A[i]);
         }
      }
 
      if(ref_masses != mod_masses)
      {
         if(!fOnlyShowZ || (fOnlyShowZ && ref_fit->GetZ()==fOnlyShowZ))
         {
            std::cout << "Z=" << ref_fit->GetZ() << " : " << refGrid->GetName() << " A=" << ref_masses.AsString();
            std::cout << "\t\t--  " << modGrid->GetName() << " A=" << mod_masses.AsString();
 
            auto offset = mod_masses.FindOffset(ref_masses);
            if(offset) {
               // check that applying the offset to the complete mass list does not include resonances
               // or very short-lived isotopes that can never be identified
               auto alist = mod_fit->GetAlist();
               KVNucleus test_nuc;
               bool ok = true;
               KVNumberList new_masses;
               for(auto &A : alist)
               {
                  A+=*offset;
                  new_masses.Add(A);
                  if(test_nuc.IsResonance(ref_fit->GetZ(),A))
                     ok = false;
                  else if(test_nuc.GetLifeTime(ref_fit->GetZ(),A).value_or(1.e+06)<1.e-3)
                     ok = false;
                  else if(not_good[ref_fit->GetZ()].Contains(A))
                     ok = false;
               }
               if(ok)
               {
                  std::cout << "\t\t==> apply DeltaA = " << (*offset>0 ? "+" : "") << *offset;
                  corrections.push_back({ref_fit->GetZ(),*offset});
               }
               else
               {
                  std::cout << "\t\t==> applying DeltaA = " << (*offset>0 ? "+" : "") << *offset << " would lead to masses: " << new_masses.AsString() << " !!!\n";
               }
            }
            else
               std::cout << "  ==> cannot find constant DeltaA to apply...???";
            std::cout << std::endl;
         }
      }
   }
   return corrections;
}
 
 
 
 
void KVIDZAFromZGridMassCorrector::ApplyCorrections(const std::vector<KVIDZAFromZGridMassCorrector::mass_correction_t>& corrections,
                                                    TH1* pid_dist,
                                                    const KVNameValueList& mass_fit_parameters)
{
   // Takes the mass corrections previously determined by DetermineMassCorrections() and applies them to each fit
   // using KVGaussFitMassModifier
   //
   // \param[in] corrections see DetermineMassCorrections()
   // \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()
 
   for(auto& corr : corrections)
   {
      if(pid_dist)
      {
         std::cout << "Applying DeltaA=" << corr.second << " to isotopes for Z=" << corr.first << std::endl;
         KVGaussFitMassModifier gfmm(modGrid, corr.first);
         gfmm.Modify(corr.second, pid_dist, mass_fit_parameters);
      }
   }
}
 
 
ClassImp(KVIDZAFromZGridMassCorrector)