KVIDZAFromZGridMassHomogenizer.cpp

#include "KVGaussFitMassModifier.h"
#include "KVIDZAFromZGridMassHomogenizer.h"
#include "KVMultiGaussIsotopeFit.h"
#include "TCanvas.h"
 
 
 
void KVIDZAFromZGridMassHomogenizer::correct_masses_for_grid(const TString& grid_name, int Z, int deltaA,
                                                             const KVNameValueList& mass_fit_parameters)
{
   // Correct the masses for a given grid/telescope
   //
   // \param[in] Z the element atomic number to correct
   // \param[in] deltaA mass modification \f$\Delta A\f$ to apply to all isotopes of element, can be positive or negative
   // \param[in] mass_fit_parameters parameters to control fit; by default KVMultiGaussIsotopeFit::mass_fit_parameters
 
   auto gr = GetGrid(grid_name);
   if(gr)
   {
      auto h2d = dynamic_cast<TH2*>(theData.Get(grid_name));
      if(!h2d)
      {
         KVError::Warning(this, "correct_masses_for_grid",
                          "No data histo with name %s found for grid %s",
                          grid_name.Data(),gr->GetName());
         return;
      }
      auto h = gr->LinearizeHistoToPID(h2d,100);
      KVGaussFitMassModifier gfmm(gr, Z);
      gfmm.Modify(deltaA, h, mass_fit_parameters);
   }
}
 
 
 
 
KVIDZAFromZGridMassHomogenizer::KVIDZAFromZGridMassHomogenizer(const TString& gridfile, const TString& datafile)
   : theData{datafile}
{
   // Prepare mass homogeneity analysis
   //
   // \param[in] gridfile name/path to file containing grids to analyse
   // \param[in] datafile name/path to ROOT file containing 2D identification maps corresponding to the grids
   //
   // \note We assume that the histogram and grid for each telescope have the same name, which allows to associate them easily!
 
   KVIDGraph::SetAutoAdd(false);
   if(theGrids.ReadAsciiFile(gridfile))
   {
      theGrids.AddLastReadGrids();
      theGrids.SetName(gridfile);
      KVError::Info(this, "KVIDZAFromZGridMassHomogenizer",
                    "Read %d reference grids from file %s",
                    theGrids.GetGrids()->GetEntries(),
                    theGrids.GetName());
   }
   else
      KVError::Error(this, "KVIDZAFromZGridMassHomogenizer",
                    "Could not read file %s",
                     gridfile.Data());
}
 
 
 
 
void KVIDZAFromZGridMassHomogenizer::Analyse(KVIDZAFromZGrid* grid, double min_yield)
{
   // Analyse a grid according to the given minimum relative yield cut-off
 
   if(!grid)
      return;
   grid->Initialize();
   if(!grid->GetFits()->GetEntries())
      return;
   KVError::Info(this, "Analyse", "Grid = %s", grid->GetName());
   TIter nxt_fit(grid->GetFits());
   KVMultiGaussIsotopeFit* fit;
   while( (fit = (KVMultiGaussIsotopeFit*)nxt_fit()) )
   {
      KVNumberList alist;
      auto yields = fit->GetRankedYields();
      int i=0;
      for(auto a : fit->GetIsotopesRankedByYield())
      {
         if(yields[i]>=min_yield*yields[0])
            alist.Add(a);
         ++i;
         if(i==3)
            break;
      }
      all_mass_ranges[fit->GetZ()][alist.AsString()].push_back(grid->GetName());
   }
}
 
 
 
 
void KVIDZAFromZGridMassHomogenizer::Print() const
{
   // Print in the terminal the mass distribution statistics obtained from the analysis
   //
   // For each element \f$Z\f$ for which a fit exists in at least 1 grid, this prints out
   // each associated 'mass list' (list of \f$A\f$-values, KVNumberList format) and a (partial) list
   // of the telescopes/grids which are concerned
 
   for(auto& stats : all_mass_ranges)
   {
      std::cout << "Z = " << stats.first << " :\n";
      for(auto& stat : stats.second)
      {
         std::cout << "\t" << stat.first << " : " << stat.second.size() << "\t{";
         int i=0, imax=5;
         for(auto&s : stat.second)
         {
            if(i) std::cout << ", ";
            ++i;
            if(i>imax)
            {
               std::cout << "... ";
               break;
            }
            std::cout << s;
         }
         std::cout << "}" << std::endl;
      }
   }
}
 
 
 
 
void KVIDZAFromZGridMassHomogenizer::ShowFits(int Z, const std::string& alist)
{
   // For a given element and 'mass list' (see Print()), display the associated multi-gaussian fits for different grids/telescopes
   // (up to 20 can be displayed).
 
   if(all_mass_ranges[Z].find(alist)!=all_mass_ranges[Z].end())
   {
      int maxNtel = 20;
      auto Ntel = (int)all_mass_ranges[Z][alist].size();
      int nx=3,ny=2;
      auto tel_to_show = std::min(Ntel,maxNtel);
      if(tel_to_show<6)
      {
         if(tel_to_show<4)
         {
            ny=1;
            nx=tel_to_show;
         }
      }
      else if(tel_to_show>6)
      {
         if(tel_to_show<10)
         {
            ny=3;
         }
         else if(tel_to_show<13)
         {
            nx=4;
            ny=3;
         }
         else if(tel_to_show<16)
         {
            nx=5;
            ny=3;
         }
         else if(tel_to_show<21)
         {
            nx=5;
            ny=4;
         }
      }
      int can_dimx=1920, can_dimy=1080;
      if(tel_to_show<4)
      {
         can_dimy=600; can_dimx=tel_to_show*can_dimy;
      }
      auto c = new TCanvas(Form("fits_Z=%d_A=%s",Z,alist.c_str()),Form("fits_Z=%d_A=%s",Z,alist.c_str()),-can_dimx,can_dimy);
      c->Divide(nx,ny);
      int itel=1;
      for(auto s : all_mass_ranges[Z][alist])
      {
         auto gr = GetGrid(s.c_str());
         gr->Initialize();
         auto pad = c->cd(itel);
         gr->GetMultiGaussFit(Z)->DrawFitWithGaussians("",s.c_str());
         // legend height = .2/5*(ngauss+1)
         auto dy = .2*(gr->GetMultiGaussFit(Z)->GetNGaussians()+1)/5.;
         pad->BuildLegend(.1,.9-dy,.45,.9);
         pad->Modified();
         pad->Update();
         ++itel;
         if(itel>tel_to_show)
            break;
      }
   }
   else
   {
      KVError::Error(this,"ShowFits","No mass range '%s' for Z=%d",alist.c_str(),Z);
   }
}
 
 
 
 
void KVIDZAFromZGridMassHomogenizer::CorrectMasses(int Z, const std::string& alist, int deltaA,
                                                   const KVNameValueList& mass_fit_parameters)
{
   // Apply the given mass-shift to the isotopes for element \f$Z\f$ in all grids corresponding to the given 'mass list' (see Print()).
   //
   // \param[in] Z element atomic number
   // \param[in] alist 'mass list' attributed to grids to modify (see Print())
   // \param[in] deltaA mass modification \f$\Delta A\f$ to apply to all isotopes of element, can be positive or negative
   // \param[in] mass_fit_parameters parameters to control fit; by default KVMultiGaussIsotopeFit::mass_fit_parameters
 
   if(all_mass_ranges[Z].find(alist)!=all_mass_ranges[Z].end())
   {
      for(auto s : all_mass_ranges[Z][alist])
      {
         correct_masses_for_grid(s.c_str(),Z,deltaA,mass_fit_parameters);
      }
   }
}
 
 
 
 
void KVIDZAFromZGridMassHomogenizer::CorrectMassesForTelescope(int Z, const TString& idtel, int deltaA, const KVNameValueList& mass_fit_parameters)
{
   // Apply the given mass-shift to the isotopes for element \f$Z\f$ of the given telescope/grid
   //
   // \param[in] Z element atomic number
   // \param[in] idtel name of grid/telescope to modify
   // \param[in] deltaA mass modification \f$\Delta A\f$ to apply to all isotopes of element, can be positive or negative
   // \param[in] mass_fit_parameters parameters to control fit; by default KVMultiGaussIsotopeFit::mass_fit_parameters
 
   correct_masses_for_grid(idtel,Z,deltaA,mass_fit_parameters);
}
 
 
 
 
void KVIDZAFromZGridMassHomogenizer::SaveGrids(const TString& filename)
{
   // Save all grids: overwrites original file if no filename given (filename=""), otherwise writes a new file
   if(filename.IsNull())
      theGrids.WriteAsciiFile(theGrids.GetName());
   else
      theGrids.WriteAsciiFile(filename);
}
 
 
ClassImp(KVIDZAFromZGridMassHomogenizer)