KVLightEnergyCsI.cpp

#include "KVLightEnergyCsI.h"
#include "TMath.h"
#include "TH1.h"
#include "TFitResultPtr.h"
#include "TFitResult.h"
 
ClassImp(KVLightEnergyCsI)
 
 
//___________________________________________________________________________
 
 
 
Double_t KVLightEnergyCsI::CalculLumiere(Double_t* x, Double_t* par)
{
   //Calculate total light output for given Z, A and energy
   //
   //~~~~~~~~~~~~~~~~~~
   // x[0] = energy (MeV)
   // par[0] = a1 : gain factor
   // par[1] = a2 : nuclear & recombination quenching term
   // par[2] = a3 : threshold (MeV/u) for delta-ray production
   // par[3] = a4 : fractional energy loss removed by delta rays
   //~~~~~~~~~~~~~~~~~~
   //
 
   Double_t energie = x[0];
   Double_t c1 = par[0];
   Double_t c2 = Z * Z * A * par[1];
   Double_t c3 = A * par[2];
   Double_t c4 = par[3];
   Double_t T = 8 * A;
   Double_t c4_new = c4 / (1. + TMath::Exp((c3 - energie) / T));
   Double_t c0 = c4 / (1. + TMath::Exp(c3 / T));
 
   Double_t lumcalc = c1 * energie;
   if (c2 > 0.0) {
      Double_t xm = -c1 * c0 * c2 * TMath::Log(c2 / (c2 + c3));
      lumcalc = lumcalc - c1 * c2 * TMath::Log(1. + energie / c2) + c1 * c2 * c4_new * TMath::Log((energie + c2) / (c3 + c2)) + xm;
   }
 
   return lumcalc;
}
 
 
 
 
Double_t KVLightEnergyCsI::InvCalc(Double_t*x, Double_t*par)
{
   // Inverse function for fitting data of calculated CsI energy vs. measured light output
   //
   // x[0] = light
   //
   // Same parameters as CalculLumiere, we plug their values into the calibration function
   // then return its inverted value
 
   GetCalibFunction()->SetParameters(par);
   return do_inversion(*x);
}
 
 
 
 
KVLightEnergyCsI::KVLightEnergyCsI(Bool_t make_func): KVCalibrator()
{
   // \param[in] make_func used by child class constructors to inhibit creation of internal calibration function
 
   SetType("LightEnergyCsI");
   if (make_func) SetCalibFunction(new TF1("fLight_CsI", this, &KVLightEnergyCsI::CalculLumiere, 0., 10000., 4));
   SetUseInverseFunction();
}
 
 
 
 
Double_t KVLightEnergyCsI::Compute(Double_t light, const KVNameValueList& z_and_a) const
{
   // Calculate the calibrated energy (in MeV) for a given total light output.
   //
   // The Z and A of the particle should be given as extra parameters:
   //
   //~~~~~~~~~~~~~~~{.cpp}
   //   Compute(light, "Z=3,A=7");
   //~~~~~~~~~~~~~~~
   //
 
   if (!CheckParameters(z_and_a)) {
      KVError::Error(this, "Compute", "[%s] Cannot compute energy for : ", GetName());
      z_and_a.ls();
      return -1;
   }
   Z = z_and_a.GetIntValue("Z");
   A = z_and_a.GetIntValue("A");
   return KVCalibrator::Compute(light, z_and_a);
}
 
 
 
 
Double_t KVLightEnergyCsI::Invert(Double_t energy, const KVNameValueList& z_and_a) const
{
   // Given the calibrated (or simulated) energy in MeV,
   // calculate the corresponding total light output according to the
   // calibration parameters (useful for filtering simulations).
   //
   // The Z and A of the particle should be given as extra parameters:
   //
   //~~~~~~~~~~~~~~~{.cpp}
   //   Invert(light, "Z=3,A=7");
   //~~~~~~~~~~~~~~~
   //
 
   if (!CheckParameters(z_and_a)) {
      KVError::Error(this, "Compute", "Cannot compute without knowing Z and A of nucleus");
      return -1;
   }
   Z = z_and_a.GetIntValue("Z");
   A = z_and_a.GetIntValue("A");
   return KVCalibrator::Invert(energy, z_and_a);
}
 
 
 
 
Bool_t KVLightEnergyCsI::CheckParameters(const KVNameValueList &z_and_a) const
{
   // \return kFALSE if parameter list does not contain "Z=..." and "A=..."
  return z_and_a.HasIntParameter("Z") && z_and_a.HasIntParameter("A");
}
 
 
 
 
void KVLightEnergyCsI::Fit(TH1 *lite_vs_e, int z, int a, std::optional<KVValueRange<Double_t> > fit_e_range)
{
   // Allows to easily fit a TProfile of mean total light output versus CsI energy
   // for a given isotope \f$(Z,A)\f$.
   //
   // Reasonable ranges for \f$a_{1}\f$ and \f$a_{2}\f$ parameters are set.
   //
   // Default value of \f$a_{3}\f$ is fixed (6); see fit_with_a3_fixed() and fit_with_a3_range() for options.
   //
   // Default value of \f$a_{4}\f$ is fixed (0.385); see fit_with_a4_fixed() and fit_with_a4_range() for options.
   Z=z; A=a;
   GetCalibFunction()->SetParameters(10.,0.5,a3_fit,a4_fit);
   GetCalibFunction()->SetParLimits(0,0.,1000.);
   GetCalibFunction()->SetParLimits(1,0.,1.);
   if(a3_range)
      GetCalibFunction()->SetParLimits(2,a3_range->Min(),a3_range->Max());
   else
      GetCalibFunction()->FixParameter(2,a3_fit);
   if(a4_range)
      GetCalibFunction()->SetParLimits(3,a4_range->Min(),a4_range->Max());
   else
      GetCalibFunction()->FixParameter(3,a4_fit);
   TFitResultPtr fit_result;
   if(fit_e_range)
      fit_result = lite_vs_e->Fit(GetCalibFunction(),"EMGS","",fit_e_range->Min(),fit_e_range->Max());
   else
      fit_result = lite_vs_e->Fit(GetCalibFunction(),"EMGS");
   std::cout << "chi**2/ndf=" << fit_result->Chi2()/fit_result->Ndf()
             << " proba=" << fit_result->Prob()*100 << "%\n";
   std::cout << GetCalibFunction()->GetParameter(0) << ","
                 << GetCalibFunction()->GetParameter(1) << ","
                     << GetCalibFunction()->GetParameter(2) << ","
                     << GetCalibFunction()->GetParameter(3) << std::endl;
}
 
 
 
 
void KVLightEnergyCsI::invFit(TH1*e_vs_lite, int z, int a, std::optional<KVValueRange<Double_t>> fit_lite_range)
{
   // Allows to easily fit a TProfile of mean calculated CsI energy vs. measured total light output
   // for a given isotope \f$(Z,A)\f$ using the inverted light-energy relation.
   //
   // Reasonable ranges for \f$a_{1}\f$ and \f$a_{2}\f$ parameters are set.
   //
   // Default value of \f$a_{3}\f$ is fixed (6); see fit_with_a3_fixed() and fit_with_a3_range() for options.
   //
   // Default value of \f$a_{4}\f$ is fixed (0.385); see fit_with_a4_fixed() and fit_with_a4_range() for options.
 
   if(!inv_fit_func)
      inv_fit_func = std::make_unique<TF1>("fInv_Light_CsI", this, &KVLightEnergyCsI::InvCalc, 0., 100000., 4);
   Z=z; A=a;
   inv_fit_func->SetParameters(10.,0.5,a3_fit,a4_fit);
   inv_fit_func->SetParLimits(0,0.,1000.);
   inv_fit_func->SetParLimits(1,0.,1.);
   if(a3_range)
      inv_fit_func->SetParLimits(2,a3_range->Min(),a3_range->Max());
   else
      inv_fit_func->FixParameter(2,a3_fit);
   if(a4_range)
      inv_fit_func->SetParLimits(3,a4_range->Min(),a4_range->Max());
   else
      inv_fit_func->FixParameter(3,a4_fit);
   TFitResultPtr fit_result;
   if(fit_lite_range)
      fit_result = e_vs_lite->Fit(inv_fit_func.get(),"EMGS","",fit_lite_range->Min(),fit_lite_range->Max());
   else
      fit_result = e_vs_lite->Fit(inv_fit_func.get(),"EMGS");
   std::cout << "chi**2/ndf=" << fit_result->Chi2()/fit_result->Ndf()
             << " proba=" << fit_result->Prob()*100 << "%\n";
   std::cout << inv_fit_func->GetParameter(0) << ","
                 << inv_fit_func->GetParameter(1) << ","
                     << inv_fit_func->GetParameter(2) << ","
                     << inv_fit_func->GetParameter(3) << std::endl;
}