KVIsoscaling.cpp

#include "KVIsoscaling.h"
#include "KVString.h"
 
#include "TGraphErrors.h"
#include "TString.h"
#include "TObjArray.h"
#include "TObjString.h"
#include "TMath.h"
#include "TF1.h"
#include "TFile.h"
#include "TAxis.h"
 
#include <fstream>
 
ClassImp(KVIsoscaling)
 
 
// BEGIN_HTML <!--
/* -->
<h2>KVIsoscaling</h2>
<h4>KVClass</h4>
<!-- */
// --> END_HTML
 
//____________________________________________________________________________//
 
 
 
void KVIsoscaling::ReadYieldsFile(const std::string& system_name, const Char_t* file_path)
{
   // Read and save yields from file provided by the user and fill the corresponding structs and maps
   // Input ASCII file format is :
   // Z A   Yield   Yield_err
   //
   // \param[in] system name
   // \param[in] path for the input ASCII file
 
   Info("ReadYieldsFile", "===== Starting reading system %s =====\n", system_name.c_str());
 
   // --- Open the file for a new system ---
   std::ifstream my_file;
   my_file.open(file_path);
 
   if (my_file.is_open()) {
 
      // - Save file path -
      fSystemList_[system_name].file_path = std::string(file_path);
 
      // --- Read the file ---
      Int_t nline = 0;
      std::string line;
      while (getline(my_file, line)) {
         TString my_str(line);
 
         // --- Ignore comments ---
         if (my_str.BeginsWith("#")) continue;
 
         // --- Extract infos ---
         TObjArray* my_str_array = my_str.Tokenize(" \t,;");
 
         Int_t zz    = ((TObjString*) my_str_array->At(0))->String().Atoi();
         Int_t aa    = ((TObjString*) my_str_array->At(1))->String().Atoi();
         Float_t yy  = ((TObjString*) my_str_array->At(2))->String().Atof();
         Float_t yy_err = ((TObjString*) my_str_array->At(3))->String().Atof();
 
         // save new data
         fSystemList_[system_name].yields[zz][aa].yield     = yy;
         fSystemList_[system_name].yields[zz][aa].yield_err = yy_err;
 
         // --- Debug ---
         if (fdebug_) printf("(line=%d) (zz=%d, aa=%d, yy=%lf, yy_err=%lf)\n", nline, zz, aa, yy, yy_err);
 
         nline++;
      }//end getline
 
      // --- Close file ---
      my_file.close();
 
      //debug
      Info("ReadYieldsFile", "[system_name=%s] file %s closed, masses and yields saved\n", system_name.c_str(), file_path);
 
   }//end file.is_opened
 
   // --- Handle the errors ---
   else {
      Error("ReadYieldsFile", "!!! Can't open file '%s' !!!", file_path);
   }
 
   // --- Apply a gaussian fit to the provided yields ---
   BuildGaussianPlots(system_name);
}
 
 
 
 
void KVIsoscaling::BuildGaussianPlots(const std::string& system_name)
{
   // Method to build the plots used to check the gaussian approximation
   // Also compute the average masses for each charge \f$\langle A(Z)\rangle\f$
   //
   // \param[in] system name
 
   if (fdebug_) Info("BuildGaussianPlots", "====== STARTING BUILDING GAUSSIAN PLOTS FOR SYSTEM %s ======", system_name.c_str());
 
   // - Find the yields for the given system -
   YieldData_t my_yields = fSystemList_[system_name].yields;
 
   // - Iterate over the yields for the given element -
   for (auto it = my_yields.begin(); it != my_yields.end(); it++){
      Int_t zz = it->first;
      Element_t my_element = it->second;
 
      // - Create then the graph for the given element -
      TGraphErrors* gr = new TGraphErrors();
      gr->SetName(Form("gr_Yield_vs_N_Z%d_%s", zz, system_name.c_str()));
      gr->GetXaxis()->SetTitle("N");
      gr->GetYaxis()->SetTitle("Yield");
      gr->SetMarkerStyle(20);
 
      TGraphErrors* grln = new TGraphErrors();
      grln->SetName(Form("gr_LnYield_vs_N_Z%d_%s", zz, system_name.c_str()));
      grln->GetXaxis()->SetTitle("N");
      grln->GetYaxis()->SetTitle("Ln(Y)");
      grln->SetMarkerStyle(20);
 
      if (fdebug_) Info("BuildGaussianPlots", "Graph %s created", gr->GetName());
 
      // - Iterate over isotopes of the given element and fill the graph -
      Int_t npp = 0;
 
      for (auto itt = my_element.begin(); itt != my_element.end(); itt++){
         Int_t aa = itt->first;
         IsotopeYield_t my_isoyield = itt->second;
 
         Int_t nn         = aa-zz;
         Float_t yield    = my_isoyield.yield;
         Float_t yielderr = my_isoyield.yield_err;
 
         gr->SetPoint(npp, nn, yield);
         gr->SetPointError(npp, 0, yielderr);
 
         grln->SetPoint(npp, nn, TMath::Log(yield));
         grln->SetPointError(npp, 0, TMath::Abs(yielderr / yield));
 
         if (fdebug_) Info("BuildGaussianPlots", "(Z=%d, N=%d, A=%d) Yield=%lf +/- %lf", zz, nn, aa, yield, yielderr);
 
         npp++;
      }//end isotope iteration
 
      // --- Apply individual gaussian fit and save the results ---
      TF1* gaus = new TF1(Form("gausfit_yield_vs_N_Z%d_%s", zz, system_name.c_str()), "gaus(0)", 0, 60);
      gr->Fit(Form("gausfit_yield_vs_N_Z%d_%s", zz, system_name.c_str()), "Q");
 
      // --- Save results ---
      fSystemList_[system_name].graphs[zz]   = gr;
      fSystemList_[system_name].lngraphs[zz] = grln;
      fSystemList_[system_name].gauss[zz]    = gaus;
 
   }//end yields iteration
 
   if (fdebug_) Info("BuildGaussianPlots", "====== END BUILDING GAUSSIAN PLOTS FOR SYSTEM %s ======", system_name.c_str());
}
 
 
 
 
void KVIsoscaling::TestGaussianApprox(const std::string& system_name1, const std::string& system_name2, Int_t zz, Double_t tol)
{
   // Method to test the tolerance threshold of the gaussian approximation for a given charge (element) \f$Z\f$.
   // Assuming a tolerance \f$tol\f$ set by the user and an average mass \f$\langle A(Z)\rangle\f$ with a standard deviation \f$ rms(Z) \f$,
   // the region of experimental points will be limited to \f$ \langle A(Z)\rangle \pm tol \cdot rms(Z) \f$.
   // This method allows to test and plot the results of various tolerance thresholds.
   //
   // \param[in] system_name1 name of the n-deficient system
   // \param[in] system_name2 name of the n-rich system
   // \param[in] zz fixed charge for the isotopic yields to be tested
   // \param[in] tol threshold to be tested
 
   TGraphErrors* gr_sys1_all = fSystemList_[system_name1].graphs[zz];
   TGraphErrors* gr_sys2_all = fSystemList_[system_name2].graphs[zz];
   TF1* gaus1 = fSystemList_[system_name1].gauss[zz];
   TF1* gaus2 = fSystemList_[system_name2].gauss[zz];
 
   // --- Build the graphs with thresholds and ratios of the yields ---
   // Use only the A (N) shared between the 2 systems
   TGraphErrors* gr_sys1_tol = new TGraphErrors();
   gr_sys1_tol->SetName(Form("gr_yields_vs_N_Z%d_%s_tol", zz, system_name1.c_str()));
   TGraphErrors* gr_sys2_tol = new TGraphErrors();
   gr_sys2_tol->SetName(Form("gr_yields_vs_N_Z%d_%s_tol", zz, system_name1.c_str()));
 
   TGraphErrors* gr_ratio_all = new TGraphErrors();
   gr_ratio_all->SetName(Form("gr_ratio_all_Z%d_%s_%s", zz, system_name1.c_str(), system_name2.c_str()));
   TGraphErrors* gr_ratio_tol = new TGraphErrors();
   gr_ratio_tol->SetName(Form("gr_ratio_tol_Z%d_%s_%s", zz, system_name1.c_str(), system_name2.c_str()));
 
   // - Extract all isotopes shared between the 2 systems -
   KVNumberList nl_inter_a = GetSharedANumberList(system_name1.c_str(), system_name2.c_str(), zz);
 
   if (fdebug_) Info("TestGaussianApproxe", "=== Starting iteration over isotope list='%s' for Z=%d ===", nl_inter_a.GetList(), zz);
 
   // - Build graphs with tolerance applied -
   Int_t npoint_all = 0;
   Int_t np1 = 0;
   Int_t np2 = 0;
   Int_t npoint_tol = 0;
   nl_inter_a.Begin();
   while (!nl_inter_a.End()) {
      Int_t next_aa = nl_inter_a.Next();
      Int_t nn = next_aa - zz;
 
      Float_t yield1 = fSystemList_[system_name1].yields[zz][next_aa].yield;
      Float_t yield2 = fSystemList_[system_name2].yields[zz][next_aa].yield;
 
      Float_t yielderr1 = fSystemList_[system_name1].yields[zz][next_aa].yield_err;
      Float_t yielderr2 = fSystemList_[system_name2].yields[zz][next_aa].yield_err;;
 
      Float_t ratio     = yield2 / yield1;
      Float_t ratio_err = TMath::Sqrt(TMath::Power(yielderr2 / yield1, 2.) + TMath::Power(yielderr1 * yield2 / yield1 / yield1, 2.));
 
      // - Fill the graph of the ratios with all points -
      gr_ratio_all->SetPoint(npoint_all, nn, ratio);
      gr_ratio_all->SetPointError(npoint_all, 0, ratio_err);
      npoint_all++;
 
      // - Apply tolerance cut for the gaussian approximation -
      Double_t diff1 = TMath::Abs(gaus1->GetParameter(1) - nn);
      Double_t diff2 = TMath::Abs(gaus2->GetParameter(1) - nn);
 
      if (diff1 < tol * gaus1->GetParameter(2)) {
         gr_sys1_tol->SetPoint(np1, nn, yield1);
         gr_sys1_tol->SetPointError(np1, 0., yielderr1);
         np1++;
      }
 
      if (diff2 < tol * gaus2->GetParameter(2)) {
         gr_sys2_tol->SetPoint(np2, nn, yield2);
         gr_sys2_tol->SetPointError(np2, 0., yielderr2);
         np2++;
      }
 
      if ((diff1 < tol * gaus1->GetParameter(2)) && (diff2 < tol * gaus2->GetParameter(2))) {
         gr_ratio_tol->SetPoint(npoint_tol, nn, ratio);
         gr_ratio_tol->SetPointError(npoint_tol, 0, ratio_err);
         npoint_tol++;
      }
   }//end shared A iter
 
   // --- Build ratio of the gaussian fits ---
   TF1* gaus_ratio = new TF1(Form("gaus_ratio_Z%d_%ss_%s", zz, system_name1.c_str(), system_name2.c_str()), "gaus(0)/(gaus(3))", 0, 60);
   gaus_ratio->SetLineColor(kBlack);
   gaus_ratio->SetParameters(gaus2->GetParameter(0), gaus2->GetParameter(1), gaus2->GetParameter(2), gaus1->GetParameter(0), gaus1->GetParameter(1), gaus1->GetParameter(2));
 
   // --- Set some draw opt ---
   gr_sys1_all->SetMarkerColor(kRed);
   gr_sys1_all->SetMarkerStyle(24);
   gr_sys2_all->SetMarkerColor(kBlue);
   gr_sys2_all->SetMarkerStyle(26);
   gr_ratio_all->SetMarkerColor(kBlack);
   gr_ratio_all->SetMarkerStyle(25);
 
   gr_sys1_tol->SetMarkerColor(kRed);
   gr_sys1_tol->SetMarkerStyle(20);
   gr_sys2_tol->SetMarkerColor(kBlue);
   gr_sys2_tol->SetMarkerStyle(22);
   gr_ratio_tol->SetMarkerColor(kBlack);
   gr_ratio_tol->SetMarkerStyle(21);
 
   gaus1->SetLineColor(kRed);
   gaus2->SetLineColor(kBlue);
   gaus_ratio->SetLineColor(kBlack);
 
   // --- Draw ---
   gr_ratio_all->Draw("ap");
   //gPad->SetLogy();
   gr_sys1_all->Draw("p");
   gr_sys2_all->Draw("p");
   gaus1->Draw("same");
   gaus2->Draw("same");
   gaus_ratio->Draw("same");
   gr_sys1_tol->Draw("p");
   gr_sys2_tol->Draw("p");
   gr_ratio_tol->Draw("p");
}
 
 
 
 
void KVIsoscaling::BuildLnR21vsNPlots(const std::string& system_name1, const std::string& system_name2)
{
   // Method to compute the ratio \f$ln(R_{21})\f$ of the yields \f$Y_{i}(A,Z)\f$ of 2 asymetric reactions \f$(2)\f$ and \f$(1)\f$
   // (projectiles and targets of reactions \f$(2)\f$ and \f$(1)\f$ only differ in their respective number of neutrons).
   //
   // We are interested in extracting the alpha isoscaling parameter, thus we want to draw \f$ln(R_{21})\f$ as a function of \f$N\f$ for a given \f$Z\f$.
   // The tolerance parameter \f$tol\f$ is used to limit the number of points used to build the plots (see GetRMSTolerance() method).
   //
   // [in] system_name1 name of the n-deficient system
   // [in] system_name2 name of the n-rich system index
 
   if (fdebug_) Info("BuildLnR21vsNPlots", "====== STARTING BUILDING PLOTS FOR %s/%s ======", system_name1.c_str(), system_name2.c_str());
 
   // --- Check the systems ---
   if(system_name1!=system_name2){
      std::string my_str = system_name1+"_"+system_name2;
 
      // --- Create a vector of the ratios for each Z (element) available ---
      KVNumberList nl_inter = GetSharedZNumberList(system_name1, system_name2);
      if (fdebug_) Info("BuildLnR21vsNPlots", "%s/%s - list of shared Z ='%s'", system_name1.c_str(), system_name2.c_str(), nl_inter.GetList());
 
      // --- Find the corresponding masses, yields and yields errors ---
      Int_t ip = 0;
      nl_inter.Begin();
      while (!nl_inter.End()) {
         Int_t next_zz = nl_inter.Next();
 
         // - Create the associated graph for ratios -
         TGraphErrors* gr_all = new TGraphErrors();
         gr_all->SetName(Form("gr_lnR21_vs_N_Z%d_%s_%s_all", next_zz, system_name1.c_str(), system_name2.c_str()));
         gr_all->GetXaxis()->SetTitle("N");
         gr_all->GetYaxis()->SetTitle("ln(R21)");
         gr_all->SetMarkerStyle(20);
         gr_all->SetMarkerColor(next_zz % 9 + 1);
 
         TGraphErrors* gr = new TGraphErrors();
         gr->SetName(Form("gr_lnR21_vs_N_Z%d_%s_%s", next_zz, system_name1.c_str(), system_name2.c_str()));
         gr->GetXaxis()->SetTitle("N");
         gr->GetYaxis()->SetTitle("ln(R21)");
         gr->SetMarkerStyle(20);
         gr->SetMarkerColor(next_zz % 9 + 1);
 
         if(fdebug_) Info("BuildLnR21vsNPlots", "Graph %s created", gr->GetName());
 
         // - Find individual fits of the yields -
         TF1* gaus_np = fSystemList_[system_name1].gauss[next_zz];
         TF1* gaus_nr = fSystemList_[system_name2].gauss[next_zz];
         Double_t par0     = gaus_nr->GetParameter(0);
         Double_t par0_err = gaus_nr->GetParError(0);
         Double_t par1     = gaus_nr->GetParameter(1);
         Double_t par1_err = gaus_nr->GetParError(1);
         Double_t par2     = gaus_nr->GetParameter(2);
         Double_t par2_err = gaus_nr->GetParError(2);
         Double_t par3     = gaus_np->GetParameter(0);
         Double_t par3_err = gaus_np->GetParError(0);
         Double_t par4     = gaus_np->GetParameter(1);
         Double_t par4_err = gaus_np->GetParError(1);
         Double_t par5     = gaus_np->GetParameter(2);
         Double_t par5_err = gaus_np->GetParError(2);
 
         // - Compute gauss2/gauss1 function -
         TF1* gaus_ratio = new TF1(Form("gaus_ratio_Z%d_%s_%s", next_zz, system_name1.c_str(), system_name2.c_str()), "gaus(0)/(gaus(3))", 0, 60);
         gaus_ratio->SetLineColor(kBlack);
         gaus_ratio->SetParameters(par0, par1, par2, par3, par4, par5);
         fIsoscalingList_[my_str].gauss_ratios[next_zz] = gaus_ratio;
 
         TF1* gaus_ratio_ln = new TF1(Form("gaus_ratio_ln_Z%d_%s_%s", next_zz, system_name1.c_str(), system_name2.c_str()), "TMath::Log(gaus(0)/(gaus(3)))", 0, 60);
         gaus_ratio_ln->SetLineColor(kBlack);
         gaus_ratio_ln->SetParameters(par0, par1, par2, par3, par4, par5);
         fIsoscalingList_[my_str].gauss_ratiosln[next_zz] = gaus_ratio_ln;
 
         // - Compute ln(R21) -
         // - Extract all A shared between the 2 systems -
         KVNumberList nl_inter_a = GetSharedANumberList(system_name1.c_str(), system_name2.c_str(), next_zz);
 
         if (fdebug_) Info("BuildLnR21vsNPlots", "=== Starting iteration over shared A list='%s' for Z=%d ===", nl_inter_a.GetList(), next_zz);
 
         Int_t npp = 0;
         Int_t npp_all = 0;
         nl_inter_a.Begin();
         while (!nl_inter_a.End()) {
            Int_t next_aa = nl_inter_a.Next();
            Int_t nn      = next_aa - next_zz;
 
            Float_t yield1 = fSystemList_[system_name1].yields[next_zz][next_aa].yield;
            Float_t yield2 = fSystemList_[system_name2].yields[next_zz][next_aa].yield;
 
            Float_t lnr21  = TMath::Log(yield2 / yield1);
 
            Float_t yielderr1 = fSystemList_[system_name1].yields[next_zz][next_aa].yield_err;
            Float_t yielderr2 = fSystemList_[system_name1].yields[next_zz][next_aa].yield_err;
            Float_t err       = TMath::Sqrt(TMath::Power(yielderr2 / yield2, 2.) + TMath::Power(yielderr1 / yield1, 2.));
            //Float_t err       = TMath::Abs(yielderr2/yield2) + TMath::Abs(yielderr1/yield1);
 
            // --- Fill the graph with all points ---
            gr_all->SetPoint(npp_all, nn, lnr21);
            gr_all->SetPointError(npp_all, 0, err);
            npp_all++;
 
            // --- Apply tolerance cut for the gaussian approximation ---
            Double_t diff1 = TMath::Abs(par1 - nn);
            Double_t diff2 = TMath::Abs(par4 - nn);
            if (diff1 < ftol_ * par2 && diff2 < ftol_ * par5) {
               gr->SetPoint(npp, nn, lnr21);
               gr->SetPointError(npp, 0, err);
               npp++;
            }
            else if (fdebug_) {
               printf("[Z=%d, N=%d, A=%d] Point refused : (tol1=%lf, diff1=%lf), (tol2=%lf, diff2=%lf) - (mean1=%lf, sigma1=%lf), (mean2=%lf, sigma2=%lf)",
                      next_zz, next_aa - next_zz, next_aa, ftol_ * par2, diff1, ftol_ * par5, diff2, par1, par2, par4, par5);
            }
            if (fdebug_) Info("BuildLnR21vsNPlots", "[Z=%d, N=%d, A=%d] Yield1(%s)=%lf, Yield2(%s)=%lf", next_zz, next_aa - next_zz, next_aa, system_name1.c_str(), yield1, system_name2.c_str(), yield2);
         }
 
         ip++;
 
        // --- Save graph for the element ---
        fIsoscalingList_[my_str].graphs[next_zz] = gr;
 
      }//end z list iter
 
      if (fdebug_) Info("BuildLnR21vsNPlots", "====== END OF BUILDING PLOTS FOR %s/%s ======", system_name1.c_str(), system_name2.c_str());
 
   }//end of systems OK
   else {
      Error("BuildLnR21Plots", "!!! Something is wrong in the input systems !!!");
   }
 
   return;
}
 
 
 
 
void KVIsoscaling::FitLnR21vsNPlots(const std::string& system_name1, const std::string& system_name2, Option_t* option, Option_t* gooption)
{
   // This method applies a linear fit to the \f$ln(R_{21})\f$ ratios to extract \f$\alpha\f$ isoscaling parameters (for each element shared between the two systems).
   //
   // \param[in] system_name1 name of the n-deficient system
   // \param[in] system_name2 name of the n-rich system
   // \param[in] option fit options
   // \param[in] gooption graphic fit options
 
   if (fdebug_) Info("FitLnR21vsNPlots", "====== STARTING FITS FOR %s/%s ======", system_name1.c_str(), system_name2.c_str());
 
   // --- Init ---
   TGraphErrors* gr = nullptr;
   TF1* func        = nullptr;
 
   std::string my_str = system_name1+"_"+system_name2;
   IsoGraph_t my_isograph = fIsoscalingList_[my_str].graphs;
 
   // --- Iterate over all graphs and apply the fit ---
   for(auto it = my_isograph.begin(); it != my_isograph.end(); it++){
      Int_t zz = it->first;
      gr       = it->second;
 
      KVNumberList nl = GetSharedANumberList(system_name1, system_name2, zz);
 
      Int_t nmin = nl.First() - zz - 2;
      Int_t nmax = nl.Last() - zz + 2;
      func = new TF1(Form("fit_lnR21_vs_N_Z%d_%s_%s", zz, system_name1.c_str(), system_name2.c_str()), "pol1", float(nmin), float(nmax));
      func->SetLineColor(zz % 9 + 1);
 
      for (Int_t ii = 0; ii < 10; ii++) gr->Fit(func, option, gooption, float(nmin + 2), float(nmax - 2));
 
      // --- Save fit results ---
      fIsoscalingList_[my_str].fits[zz] = func;
 
   }//end graph iteration
 
   if (fdebug_) Info("BuildLnR21vsNPlots", "====== END OF FITS FOR %s/%s ======", system_name1.c_str(), system_name2.c_str());
}
 
 
 
 
void KVIsoscaling::DrawLnR21vsNFits(const std::string& system_name1, const std::string& system_name2)
{
   // This method draws the results of \f$ln(R_{21})\f$ vs \f$N\f$ linear fits for the given system pair.
   //
   // \param[in] system_name1 name of the n-deficient system
   // \param[in] system_name2 name of the n-rich system
 
   std::string my_str = system_name1+"_"+system_name2;
 
   TMultiGraph* mgr  = new TMultiGraph();
   TGraphErrors* gr  = nullptr;
   TF1*         func = nullptr;
 
   // --- Iterate over all graphs ---
   IsoGraph_t my_isograph = fIsoscalingList_[my_str].graphs;
   for(auto it = my_isograph.begin(); it != my_isograph.end(); it++){
      Int_t zz = it->first;
      gr       = it->second;
 
      mgr->Add(gr);
   }//end graph iter
 
   mgr->Draw("ap");
   mgr->GetXaxis()->SetTitle("N");
   mgr->GetYaxis()->SetTitle("ln(R21)");
 
   // --- Iterate over all fits ---
   IsoFit_t my_isofit = fIsoscalingList_[my_str].fits;
   for(auto itt = my_isofit.begin(); itt != my_isofit.end(); itt++){
      Int_t zz = itt->first;
      func     = itt->second;
 
      func->Draw("same");
   }//end fit iter
}
 
 
 
 
void KVIsoscaling::BuildIsoscalingPlots(const std::string& system_name1, const std::string& system_name2, Int_t mcolor, Int_t mstyle, Bool_t draw)
{
    // This method allows to create isoscaling results plots for a given system pair :
    //
    // \f$\alpha\f$ vs \f$Z\f$
    // \f$\Delta\f$ vs \f$Z\f$
    // \f$\alpha\f$ vs \f$\Delta\f$
    // \f$Csym/T\f$ vs \f$Z\f$
    //
    // To compute \f$C_{sym}/T\f$, we use the usual isoscaling formula (see for ex. Raduta & Gulminelli, PRC 75 044605 bis (2007), Eq.(16)) :
    //\f[
    // C_{sym}(\langle A(Z) \rangle)/T = \alpha(Z)/\Delta(Z)
    //\f]
    //
    // where
    //\f[
    //\Delta = (Z/\langle A_{1}(Z) \rangle)^{2}-(Z/\langle A_{2}(Z) \rangle)^{2}
    //\f]
    //
    // \param[in] system_name1 name of the n-deficient system
    // \param[in] system_name2 name of the n-rich system
    // \param[in] mcolor color to use for drawing
    // \param[in] mstyle style to use for drawing
    // \param[in] draw =kTRUE if we want to plot the result
 
    if (fdebug_) Info("BuildIsoscalingPlots", "====== STARTING BUILDING PLOTS FOR %s/%s ======", system_name1.c_str(), system_name2.c_str());
    std::string my_str = system_name1+"_"+system_name2;
 
    // --- Create the graphs ---
    TGraphErrors* gr_alpha = new TGraphErrors();
    gr_alpha->SetName(Form("Alpha_vs_Z_%s", my_str.c_str()));
    gr_alpha->SetMarkerStyle(20);
    gr_alpha->GetXaxis()->SetTitle("Z");
    gr_alpha->GetYaxis()->SetTitle("#alpha");
    gr_alpha->SetMarkerStyle(mstyle);
    gr_alpha->SetMarkerColor(mcolor);
 
    TGraphErrors* gr_delta = new TGraphErrors();
    gr_delta->SetName(Form("Delta_vs_Z_%s", my_str.c_str()));
    gr_delta->SetMarkerStyle(20);
    gr_delta->GetXaxis()->SetTitle("Z");
    gr_delta->GetYaxis()->SetTitle("#Delta");
    gr_delta->SetMarkerStyle(mstyle);
    gr_delta->SetMarkerColor(mcolor);
 
    TGraphErrors* gr_alpha_delta = new TGraphErrors();
    gr_alpha_delta->SetName(Form("Alpha_vs_Delta_%s", my_str.c_str()));
    gr_alpha_delta->SetMarkerStyle(20);
    gr_alpha_delta->GetYaxis()->SetTitle("#alpha");
    gr_alpha_delta->GetXaxis()->SetTitle("#Delta");
    gr_alpha_delta->SetMarkerStyle(mstyle);
    gr_alpha_delta->SetMarkerColor(mcolor);
 
    TGraphErrors* gr_csymT = new TGraphErrors();
    gr_csymT->SetName(Form("CsymOverT_vs_Z_%s", my_str.c_str()));
    gr_csymT->SetMarkerStyle(mstyle);
    gr_csymT->SetMarkerColor(mcolor);
    gr_csymT->GetXaxis()->SetTitle("Z");
    gr_csymT->GetYaxis()->SetTitle("#alpha/4#Delta");
 
    // --- Iterate over shared elements ---
    KVNumberList nl = GetSharedZNumberList(system_name1, system_name2);
 
    Int_t np = 0;
    nl.Begin();
    while (!nl.End()) {
       Int_t zz = nl.Next();
 
       Float_t alpha, alpha_err, delta, delta_err, csymT, csymT_err;
       Bool_t is_ok_alpha = GetAlpha(system_name1, system_name2, zz, alpha, alpha_err);
       Bool_t is_ok_delta = GetDeltaZA2(system_name1, system_name2, zz, delta, delta_err, fdebug_);
       Bool_t is_ok_csymT = GetCsymOverT(system_name1, system_name2, zz, csymT, csymT_err, fdebug_);
 
       if (is_ok_delta && is_ok_alpha && is_ok_csymT) {
 
          gr_alpha->SetPoint(np, zz, alpha);
          gr_alpha->SetPointError(np, 0., alpha_err);
 
          gr_delta->SetPoint(np, zz, delta);
          gr_delta->SetPointError(np, 0., delta_err);
 
          gr_alpha_delta->SetPoint(np, delta, alpha);
          gr_alpha_delta->SetPointError(np, delta_err, alpha_err);
 
          gr_csymT->SetPoint(np, zz, csymT);
          gr_csymT->SetPointError(np, 0., csymT_err);
 
          np++;
       }//end ok
    }//end element iteration
 
    // --- Save results ---
    fIsoscalingList_[my_str].gr_alpha = gr_alpha;
    fIsoscalingList_[my_str].gr_delta = gr_delta;
    fIsoscalingList_[my_str].gr_alpha_delta = gr_alpha_delta;
    fIsoscalingList_[my_str].gr_csymT = gr_csymT;
 
    // --- Draw ---
    //if (draw) gr->Draw("ap");
 
    if (fdebug_) Info("BuildIsoscalingPlots", "====== END BUILDING PLOTS FOR %s/%s ======", system_name1.c_str(), system_name2.c_str());
}
 
 
 
 
void KVIsoscaling::CreateCsymOverTMultiGraph(TMultiGraph* mgr)
{
   // This method creates a TMultiGraph containing all the \f$C_{sym}/T\f$ isoscaling results (for all combination of systems pairs already tested).
   //
   // \param[in] mgr pointer to the TMultiGraph provided by the user
   TGraphErrors* gr = nullptr;
 
   // --- Iterate over all systems ---
   for(auto it = fIsoscalingList_.begin(); it != fIsoscalingList_.end(); it++){
      std::string str_sys = it->first;
      Isoscaling_t my_iso = it->second;
 
      gr =  my_iso.gr_csymT;
      mgr->Add(gr);
 
      Info("CreateCsymOverTMultiGraph", "%s added to the MultiGraph", gr->GetName());
   }
}
 
 
 
 
void KVIsoscaling::SaveResultsROOT(const Char_t* file_name)
{
   // This method allows to save the yields and the isoscaling fits and graphs (\f$ln(R_{21})\f$, \f$\alpha\f$, \f$\Delta\f$, \f$C_{sym}/T\f$ graphs) in a *.root file.
   //
   // \param[in] file_name output file path
 
   TFile my_file(Form("%s", file_name), "RECREATE");
 
   // --- Save yields ---
   TGraphErrors *gr_yield = nullptr;
   TGraphErrors *gr_lnyield = nullptr;
   TF1 *gr_gauss = nullptr;
 
   for(auto it0 = fSystemList_.begin(); it0 != fSystemList_.end(); it0++){
      std::string str_sys = it0->first;
      System_t my_system  = it0->second;
 
      // - yields -
      YieldGraph_t my_graphs = my_system.graphs;
      for(auto itt0 = my_graphs.begin(); itt0 != my_graphs.end(); itt0++){
         Int_t zz = itt0->first;
         gr_yield = itt0->second;
 
         gr_yield->Write();
      }//end yields graph iteration
 
      // - ln(yield) -
      YieldGraph_t my_lngraphs = my_system.lngraphs;
      for(auto ittt0 = my_lngraphs.begin(); ittt0 != my_lngraphs.end(); ittt0++){
         Int_t zz   = ittt0->first;
         gr_lnyield = ittt0->second;
 
         gr_lnyield->Write();
      }//end Ln(yield) graphs iteration
 
      // - gaussian fits -
      YieldGauss_t my_gauss = my_system.gauss;
      for(auto itttt0 = my_gauss.begin(); itttt0 != my_gauss.end(); itttt0++){
         Int_t zz = itttt0->first;
          gr_gauss= itttt0->second;
 
         gr_gauss->Write();
      }//end fits iteration
   }//end yield iteration
 
   // --- Save isoscaling results ---
   TGraphErrors *gr_lnR21 = nullptr;
   TF1          *fit_lnR21= nullptr;
   TGraphErrors *gr_alpha = nullptr;
   TGraphErrors *gr_delta = nullptr;
   TGraphErrors *gr_alpha_delta = nullptr;
   TGraphErrors *gr_csymT = nullptr;
 
   for(auto it = fIsoscalingList_.begin(); it != fIsoscalingList_.end(); it++){
      std::string str_sys = it->first;
      Isoscaling_t my_iso = it->second;
 
      IsoGraph_t my_graphs = my_iso.graphs;
      for(auto itt = my_graphs.begin(); itt != my_graphs.end(); itt++){
         Int_t zz = itt->first;
         gr_lnR21 = itt->second;
 
         gr_lnR21->Write();
      }//end lnR21 graph iteration
 
      IsoFit_t my_fits = my_iso.fits;
      for(auto ittt = my_fits.begin(); ittt != my_fits.end(); ittt++){
         Int_t zz  = ittt->first;
         fit_lnR21 = ittt->second;
 
         fit_lnR21->Write();
      }//end lnR21 fits iteration
 
      gr_alpha = my_iso.gr_alpha;
      gr_delta = my_iso.gr_delta;
      gr_alpha_delta = my_iso.gr_alpha_delta;
      gr_csymT = my_iso.gr_csymT;
 
      gr_alpha->Write();
      gr_delta->Write();
      gr_csymT->Write();
      gr_alpha_delta->Write();
   }//end isoscaling list iter
 
   if (fdebug_) Info("SaveResultsROOT", "File '%s' saved", file_name);
}
 
 
 
 
void KVIsoscaling::SaveResultsASCII(const Char_t* file_name)
{
   // This method allows to save results in a *.txt file.
   //
   // \param[in] file_name output file path
   KVString fname(file_name);
 
   std::ofstream my_file;
   my_file.open(fname.Data(), std::ofstream::out);
   my_file << "### system_pair  Z   Alpha   d(Alpha)    <A_1>   d<A_1>  <A_2>   d<A_2>  denum   d(denum) Csym/T     d(Csym/T) ###" << std::endl;
 
   // --- Iterate over all isoscaling results ---
   for(auto it0 = fIsoscalingList_.begin(); it0 != fIsoscalingList_.end(); it0++){
      std::string str_sys = it0->first;
      Isoscaling_t my_iso = it0->second;
 
      // !!! only works if system_name does not contain _ !!!
      KVString mystr(str_sys);
      TObjArray *arr = mystr.Tokenize("_");
      std::string sys1(((TObjString*) arr->At(0))->GetString().Data());
      std::string sys2(((TObjString*) arr->At(1))->GetString().Data());
 
      my_file << "### " << sys1 << " " << sys2 << " ###" << std::endl;
 
      IsoFit_t my_isofit = my_iso.fits;
 
      // --- Iterate over all elements for the given isoscaling fit result ---
      for(auto it1 = my_isofit.begin(); it1 != my_isofit.end(); it1++){
          Int_t zz = it1->first;
          TF1 *fit = it1->second;
 
          // Alpha
          Float_t alpha = -666.;
          Float_t alpha_err = -666.;
          GetAlpha(sys1, sys2, zz, alpha, alpha_err);
 
          // Denum
          Float_t denum = -666.;
          Float_t denum_err = -666.;
          GetDeltaZA2(sys1, sys2, zz, denum, denum_err, kFALSE);
 
          //Csym/T
          Float_t csymT = -666.;
          Float_t csymT_err = -666.;
          GetCsymOverT(sys1, sys2, zz, csymT, csymT_err, kFALSE);
 
          //<A1> and <A2>
          Float_t meanA1 = -666.;
          Float_t meanA1_err = -666.;
          Float_t meanA2 =  -666.;
          Float_t meanA2_err = -666.;
          GetAMean(sys1, zz, meanA1, meanA1_err);
          GetAMean(sys2, zz, meanA2, meanA2_err);
 
 
          KVString line(Form("%s %s %d %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf", sys1.c_str(), sys2.c_str(), zz, alpha, alpha_err, meanA1, meanA1_err, meanA2, meanA2_err, denum, denum_err, csymT, csymT_err));
          my_file << line.Data() << std::endl;
      }//end element iteration
 
      if (fdebug_) Info("SaveResultsASCII", "[%s/%s] results saved\n", sys1.c_str(), sys2.c_str());
 
   }//end isoscaling results iteration
 
   my_file.close();
 
   if (fdebug_) Info("SaveResultsASCII", "File '%s' saved", fname.Data());
}
 
 
 
//===================================
//==        GETTERS YIELDS         ==
//===================================
 
 
Bool_t KVIsoscaling::GetAMean(const std::string& system_name, Int_t zz, Float_t& meanA, Float_t& meanA_err)
{
   // This method returns the average mass \f$ \langle A(Z) \rangle \f$ of the given element (charge \f$Z\f$) for the given system.
   // The average mass and associated uncertainty are obtained from the result of a gaussian fit applied to the isotopic distribution.
   //
   // \param[in] system_name name of the system
   // \param[in] zz charge of the element
   // \param[out] meanA average mass
   // \param[out] meanA_err error on the average mass
   Bool_t is_ok = kTRUE;
 
   // find the corresponding graph
   KVNumberList nl =  GetZNumberList(system_name);
 
   if(nl.Contains(zz)){
      TF1* my_gaus = fSystemList_[system_name].gauss[zz];
      meanA    = my_gaus->GetParameter(1);
      meanA_err= my_gaus->GetParError(1);
   }
   else{
      Error("GetAMean", "[sys=%s] Z=%d is out of range", system_name.c_str(), zz);
      is_ok = kFALSE;
   }
 
   return is_ok;
}
 
 
 
 
Int_t KVIsoscaling::FindZFromAmean(const std::string& system_name, Int_t aa)
{
   // Returns the charge \f$Z\f$ (element) associated to the provided mass \f$A\f$ such as \f$\langle A(Z) \rangle = A \f$.
   //
   // \param[in] system_name name of the system
   // \param[in] aa mass
 
   Int_t zz = -666;
   Double_t min_val = 10.;
 
   KVNumberList nl = GetZNumberList(system_name);
   nl.Begin();
   while (!nl.End()) {
      Int_t next_zz = nl.Next();
 
      // extract mean aa
      Float_t aa_real, aa_real_err;
      GetAMean(system_name, next_zz, aa_real, aa_real_err);
 
      Double_t delta_aa = TMath::Abs(aa - aa_real);
      if (delta_aa < min_val) {
         zz = next_zz;
         min_val = delta_aa;
 
         if (fdebug_) Info("FindZFromAmean", "[%s - A=%d] : Z=%d -> <A(Z)>=%lf +/- %lf [FOUND] delta=%lf, min_val=%lf", system_name.c_str(), aa, next_zz, aa_real, aa_real_err, delta_aa, min_val);
      }
      else if (fdebug_) Info("FindZFromAmean", "[%s - A=%d] : Z=%d -> <A(Z)>=%lf +/- %lf [NOT FOUND] delta=%lf, min_val=%lf", system_name.c_str(), aa, next_zz, aa_real, aa_real_err, delta_aa, min_val);
   }
 
   return zz;
}
 
 
 
//=======================================
//==        GETTERS ISOSCALING         ==
//=======================================
 
 
Bool_t KVIsoscaling::GetAlpha(const std::string& system_name1, const std::string& system_name2, Int_t zz, Float_t& alpha, Float_t& alpha_err)
{
   // This method extracts the isoscaling \f$\alpha\f$ parameter.
   //
   // \param[in] system_name1 name of the n-deficient system
   // \param[in] system_name2 name of the n-rich system
   // \param[in] zz charge for the element the isoscaling fit was applied
   // \param[out] alpha isoscaling fit parameter \f$\alpha\f$
   // \param[out] alpha_err error on the isoscaling fit parameter \f$\alpha\f$
 
   Bool_t is_ok = kTRUE;
 
   // --- Find the corresponding alpha parameter in fit function list ---
   std::string my_str = system_name1+"_"+system_name2;
   TF1* func = fIsoscalingList_[my_str].fits[zz];
 
   if(func!=nullptr){
      alpha     = func->GetParameter(1);
      alpha_err = func->GetParError(1);
   }
   else{
       Error("GetAlpha", "!!! (%s / %s) [Z=%d] : Fit function not found !!!", system_name1.c_str(), system_name2.c_str(), zz);
       is_ok=kFALSE;
   }
 
   return is_ok;
}
 
 
 
 
Bool_t KVIsoscaling::GetDeltaZA2(const std::string& system_name1, const std::string& system_name2, Int_t zz, Float_t& denum, Float_t& denum_err, Bool_t debug)
{
   // This method returns the value of the asymmetry parameter \f$ \Delta \f$ (denumerator for \f$ C_{sym}/T \f$ estimation), such as
   //\f[
   //\Delta = (Z/\langle A_{1} \rangle)^{2}-(Z/ \langle A_{2} \rangle)^{2}
   //\f]
   //
   // \param[in] system_name1 name of the n-deficient system
   // \param[in] system_name2 name of the n-rich system
   // \param[in] zz charge for the element the isoscaling fit was applied
   // \param[out] denum asymmetry parameter \f$\Delta\f$
   // \param[out] denum_err error on the asymmetry parameter \f$\Delta\f$
 
   Bool_t is_ok = kFALSE;
 
   Float_t meanA1     = -666.;
   Float_t meanA2     = -666.;
   Float_t meanA1_err = -666.;
   Float_t meanA2_err = -666.;
 
   Bool_t ok1 = GetAMean(system_name1, zz, meanA1, meanA1_err);
   Bool_t ok2 = GetAMean(system_name2, zz, meanA2, meanA2_err);
 
   if (ok1 && ok2) {
      is_ok = kTRUE;
 
      denum     = TMath::Power(1.*zz / meanA1, 2.) - TMath::Power(1.*zz / meanA2, 2.);
      denum_err = TMath::Sqrt(TMath::Power(zz * zz * 2.*meanA1_err / meanA1 / meanA1 / meanA1, 2.) + TMath::Power(zz * zz * 2.*meanA2_err / meanA2 / meanA2 / meanA2, 2.));
 
      if (debug) Info("GetDeltaZA2", "(%s / %s) [Z=%d] : denum=%lf +/- %lf", system_name1.c_str(), system_name2.c_str(), zz, denum, denum_err);
   }
   else {
      if (!ok1) Warning("GetDeltaZA2", "<A1> was not found for Z=%d of system %s/%s", zz, system_name1.c_str(), system_name2.c_str());
      if (!ok2) Warning("GetDeltaZA2", "<A2> was not found for Z=%d of system %s/%s", zz, system_name1.c_str(), system_name2.c_str());
   }
 
   return is_ok;
}
 
 
 
 
Bool_t KVIsoscaling::GetCsymOverT(const std::string& system_name1, const std::string& system_name2, Int_t zz, Float_t& csymT, Float_t& csymT_err, Bool_t debug)
{
   // This method returns the value of \f$C_{sym}/T\f$ from isoscaling method such as (see for ex. Raduta & Gulminelli, PRC 75 044605 (2007) Eq.(16))
   //\f[
   // C_{sym}/T = \alpha/4\Delta
   //\f]
   //
   // \param[in] system_name1 name of the n-deficient system
   // \param[in] system_name2 name of the n-rich system
   // \param[in] zz charge for the element the isoscaling fit was applied
   // \param[out] csymT value of the ratio \f$C_{sym}/T = \alpha/4\Delta\f$
   // \param[out] csymT_err error on the ratio \f$C_{sym}/T = \alpha/4\Delta\f$
 
   Bool_t is_ok = kFALSE;
   Bool_t is_oka = kFALSE;
   Bool_t is_okd = kFALSE;
 
   // --- Get corresponding alpha isoscaling parameter ---
   Float_t alpha       = -666.;
   Float_t alpha_err = -666.;
   is_oka = GetAlpha(system_name1, system_name2, zz, alpha, alpha_err);
 
   // --- Get the denum ---
   Float_t denum       = -666.;
   Float_t denum_err = -666.;
   is_okd = GetDeltaZA2(system_name1, system_name2, zz, denum, denum_err, debug);
 
   // --- Compute Csym/T ---
   if (is_oka && is_okd) {
      is_ok = kTRUE;
      csymT     = alpha / 4. / denum;
      csymT_err = TMath::Sqrt(TMath::Power(alpha_err / 4. / denum, 2.) + TMath::Power(alpha * denum_err / 4. / denum / denum, 2.));
 
      if (debug) {
         Info("GetCsymOverT", "(%s / %s) [Z=%d] : alpha=%lf +/- %lf, denum=%lf +/- %lf, Csym/T=%lf +/- %lf", system_name1.c_str(), system_name2.c_str(), zz, alpha, alpha_err, denum, denum_err, csymT, csymT_err);
      }
   }
 
   return is_ok;
}
 
 
 
 
KVNumberList KVIsoscaling::GetZNumberList(const std::string& system_name)
{
    // This method returns the list of elements of a given system
    //
    // \param[in] system_name name of the system
    //
    // \returns KVNumberList of the all the charges \f$Z\f$
 
    KVNumberList nl;
    YieldData_t my_yields = fSystemList_[system_name].yields;
 
    for(auto it = my_yields.begin(); it != my_yields.end(); it++){
       Int_t zz = it->first;
       nl.Add(zz);
    }
 
    return nl;
}
 
 
 
 
KVNumberList KVIsoscaling::GetSharedZNumberList(const std::string& system_name1, const std::string& system_name2)
{
   // This method returns the list of all elements shared between the two provided systems.
   //
   // \param[in] system_name1 name of the first system
   // \param[in] system_name2 name of the second system
   //
   // \returns KVNumberList of the shared charges \f$Z\f$
 
   KVNumberList nl_z1 = GetZNumberList(system_name1);
   KVNumberList nl_z2 = GetZNumberList(system_name2);
   KVNumberList nl_inter;
   nl_z1.Copy(nl_inter);
   nl_inter.Inter(nl_z2);
 
   if (fdebug_) {
      Info("GetSharedZNumberList", "Zlist[%s]= %s, Zlist[%s]= %s, shared Zlist= %s",
           system_name1.c_str(), nl_z1.GetList(), system_name2.c_str(), nl_z2.GetList(), nl_inter.GetList()) ;
   }
 
   return nl_inter;
}
 
 
 
 
KVNumberList KVIsoscaling::GetANumberList(const std::string& system_name, Int_t zz)
{
    // This method returns the list of all isotopes for the given element.
    //
    // \param[in] system_name name of the system
    // \param[in] zz charge of the element
    //
    // \returns KVNumberList of the masses \f$A\f$
 
    KVNumberList nl;
    Element_t my_element = fSystemList_[system_name].yields[zz];
 
    for(auto it = my_element.begin(); it != my_element.end(); it++){
       Int_t aa = it->first;
       nl.Add(aa);
    }
 
    return nl;
}
 
 
 
 
KVNumberList KVIsoscaling::GetSharedANumberList(const std::string& system_name1, const std::string& system_name2, Int_t zz)
{
   // This method returns the list of all isotopes shared between the two provided systems for the given element.
   //
   // \param[in] system_name1 name of the first system
   // \param[in] system_name2 name of the second system
   // \param[in] zz charge of the element
   //
   // \returns KVNumberList of the shared masses \f$A\f$
 
   KVNumberList nl_a1 = GetANumberList(system_name1, zz);
   KVNumberList nl_a2 = GetANumberList(system_name2, zz);
   KVNumberList nl_inter_a;
   nl_a1.Copy(nl_inter_a);
   nl_inter_a.Inter(nl_a2);
 
   if (fdebug_){
      Info("GetSharedANumberList", "Z=%d : Alist[%s]= %s, Alist[%s]= %s, shared Alist= %s",
           zz, system_name1.c_str(), nl_a1.GetList(), system_name2.c_str(), nl_a2.GetList(), nl_inter_a.GetList()) ;
   }
 
   return nl_inter_a;
}
 
 
//============================
//==        PRINTERS        ==
//============================
 
 
void KVIsoscaling::PrintYieldsList()
{
   Info("PrintYields", "Printing the list of yields : ");
 
   // --- Iterate over systems ---
   for(auto it0 = fSystemList_.begin(); it0 != fSystemList_.end(); it0++){
 
      std::string sys_name = it0->first;
      System_t my_system   = it0->second;
 
      std::string my_path   = my_system.file_path;
      YieldData_t my_yields = my_system.yields;
 
      printf("\n====== %s =====\n", sys_name.c_str());
      printf("extracted from path = '%s'\n", my_path.c_str());
 
      // --- Iterate over elements for the given system ---
      for(auto it1 = my_yields.begin(); it1 != my_yields.end(); it1++){
         Int_t zz         = it1->first;
         Element_t my_ele = it1->second;
 
         // --- Iterate over isostopes of the given element ---
         for(auto it2 = my_ele.begin(); it2 != my_ele.end(); it2++){
             Int_t aa             = it2->first;
             IsotopeYield_t yield = it2->second;
 
             printf("   Y(Z=%d, A=%d) = %lf +/- %lf\n", zz, aa, yield.yield, yield.yield_err);
         }//end isostope iteration
      }//end element iteration
   }//end system iteration
 
   Info("PrintYields", "That's all folks !");
}
 
 
 
 
void KVIsoscaling::PrintSystemsList()
{
   Info("PrintSystemsList", "Printing the list of available systems : ");
 
   // --- Iterate over systems ---
   for(auto it0 = fSystemList_.begin(); it0 != fSystemList_.end(); it0++){
 
      std::string sys_name = it0->first;
      System_t my_system   = it0->second;
 
      std::string my_path   = my_system.file_path;
      YieldData_t my_yields = my_system.yields;
 
      printf("\n====== %s =====\n", sys_name.c_str());
      printf("extracted from path = '%s'\n", my_path.c_str());
   }//end system iteration
}