JHVInterpolator.hh

Go to the documentation of this file.

#ifndef __JCALIBRATE_JHVINTERPOLATOR__
#define __JCALIBRATE_JHVINTERPOLATOR__

#include "JLang/JException.hh"

#include "JTools/JRange.hh"

#include "JCalibrate/JFitToT.hh"

#include "TKey.h"
#include "TGraph.h"
#include "TGraphErrors.h"
#include "TMultiGraph.h"
#include "TAttMarker.h"
#include "TColor.h"


/**
 * \author bjung
 */

namespace JCALIBRATE {

  using JTOOLS::JRange;
  using JLANG::JNoValue;
  using JLANG::JValueOutOfRange;
  using JLANG::JNullPointerException;
  
  static const std::string HVINTERPOLATOR_SUFFIX = ".HVxG";
  
  static const char*       HVINTERPOLATOR_DATA   = "HVINTERPOLATOR_DATA";
  static const char*       HVINTERPOLATOR_RESULT = "HVINTERPOLATOR_RESULT";
  
  /**
   * Auxiliary data structure to store high-voltage versus gain data and interpolate the nominal high-voltage.
   */
  struct JHVInterpolator
  {
    /**
     * Constructor.
     *
     * \param  object             TMultiGraph object
     */
    JHVInterpolator(TMultiGraph& object) :
      data(&object)
    {
      TGraphErrors* g0 = new TGraphErrors();

      g0->SetName(HVINTERPOLATOR_RESULT);
      g0->SetLineColor  (kRed);
      g0->SetMarkerColor(kRed);
      g0->SetMarkerStyle(kFullDotSmall);
      g0->Set(0);

      data->Add(g0);
      
      TGraphErrors* g1 = new TGraphErrors();
      
      g1->SetName(HVINTERPOLATOR_DATA);
      g1->SetMarkerStyle(kFullDotSmall);
      g1->Set(0);
      
      data->Add(g1);
    }
    
    
    /**
     * Add point to diagram.
     *
     * \param  HV                 high-voltage        [V]
     * \param  gain               gain value
     * \param  gainError          error on gain value
     */
    void AddPoint(Double_t HV, Double_t gain, Double_t gainError)
    {      
      TGraphErrors* g1 = getData();
      
      const Int_t   N  = g1->GetN();
      
      g1->SetPoint     (N, fabs(HV), gain);
      g1->SetPointError(N, 0.0, gainError);
    }
    
    
    /**
     * Set point with index i.
     *
     * \param  i                  index of point
     * \param  HV                 high-voltage        [V]
     * \param  gain               gain value
     * \param  gainError          error on gain value
     */
    void SetPoint(Int_t i, Double_t HV, Double_t gain, Double_t gainError)
    {
      TGraphErrors* g1 = getData();

      g1->SetPoint     (i, fabs(HV), gain);
      g1->SetPointError(i, 0.0, gainError);
    }


    /**
     * Checks whether high-voltage is within range
     *
     * \param  HV                 high-voltage [V]
     * \return                    true if high-voltage is within range; else false
     */
    const bool checkHV(const double HV) const
    {
      return (HV > hvRange.getLowerLimit() &&
              HV < hvRange.getUpperLimit());
    }
    
    
    /**
     * Checks whether two high-voltage values are different.
     *
     * \param  HV1                first  high-voltage [V]
     * \param  HV2                second high-voltage [V]
     * \return                    true if absolute difference is greater than minimal required difference; else false
     */
    const bool checkHV(const double HV1, const double HV2) const
    {
      return (fabs(HV1 - HV2) > dHVmin);
    }
    
    
    /**
     * Checks if gain is within range.
     *
     * \param   gain              gain value
     * \return                    true if gain within allowed range; else false
     */
    const bool checkGain(const double gain) const
    {
      return (gain > gainRange.getLowerLimit() &&
              gain < gainRange.getUpperLimit());
    }
    
    
    /**
     * Checks whether the gains of two points are strictly increasing as function of their absolute high-voltage
     *
     * \param  i                  index of first  point
     * \param  j                  index of second point
     * \return                    true if gains of given points are strictly increasing\n 
     *                            as function of the absolute high-voltage; else false
     */
    const bool areIncreasing(const Int_t i, const Int_t j) const
    {
      TGraphErrors* g1 = getData();
      
      if ((i >= 0 && i < g1->GetN()) &&
          (j >= 0 && i < g1->GetN())) {
          
        return ((g1->GetPointX(i) > g1->GetPointX(j) && g1->GetPointY(i) > g1->GetPointY(j)) ||
                (g1->GetPointX(i) < g1->GetPointX(j) && g1->GetPointY(i) < g1->GetPointY(j)));
        
      } else {

        return false;
      }
    }
    
    
    /**
     * Checks whether two points are valid for inter-/extrapolation.
     *
     * \param  i                  index of first  point
     * \param  j                  index of second point
     * \return                    true if valid for inter-/extrapolation; else false
     */
    const bool areValid(const Int_t i, const Int_t j) const
    {
      TGraphErrors* g1 = getData();
      
      if ((i >= 0 && i < g1->GetN()) &&
          (j >= 0 && j < g1->GetN())) {

        const double HV_i   = -fabs(g1->GetPointX(i));
        const double HV_j   = -fabs(g1->GetPointX(j));
        const double gain_i =       g1->GetPointY(i);
        const double gain_j =       g1->GetPointY(j);
        
        return (areIncreasing(   i,    j) && checkHV(HV_i) && checkGain(gain_i) &&
                checkHV      (HV_i, HV_j) && checkHV(HV_j) && checkGain(gain_j));
        
      } else {
        
        return false;
      }
    }
    
    
    /**
     * Inter-/Extrapolate the high-voltage value corresponding to the target gain value.
     *
     * \param  gainTarget         target gain value for inter-/extrapolation
     * \return                    true if inter-/extrapolation successful; else false
     */
    bool interpolateHV(const double gainTarget)
    {
      TGraphErrors* g1 = getData();
      
      if (g1->GetN() < 2 || !checkGain(gainTarget)) {
        return false;
      }
      
      // Search for valid inter-/extrapolation points
      
      Int_t i = 0;
      Int_t j = 1;
      
      for (Int_t k = 2; k < g1->GetN(); ++k) {

        const double dGain_i = fabs(g1->GetPointY(i) - gainTarget);     
        const double dGain_j = fabs(g1->GetPointY(j) - gainTarget);
        const double dGain_k = fabs(g1->GetPointY(k) - gainTarget);

        if         (dGain_k < dGain_j) {
          i = (areValid(j,k) ? j : i);
          j = k;
        } else if ((dGain_k < dGain_i || !areValid(i,j)) && areValid(j,k)) {
          i = j;
          j = k;
        }
      }

      // Inter-/Extrapolate high-voltage corresponding to given gain
      
      if (areValid(i, j)) {
        
        const double logHV0 = log(fabs(g1->GetPointX(i)));
        const double logHV1 = log(fabs(g1->GetPointX(j)));
        
        const double logG0  = log(g1->GetPointY(i));
        const double logG1  = log(g1->GetPointY(j));
        const double elogG0 = g1->GetErrorY(i) / g1->GetPointY(i);
        const double elogG1 = g1->GetErrorY(j) / g1->GetPointY(j);
        
        const double dlogG0 = log(gainTarget) - logG0;
        const double dlogG1 = log(gainTarget) - logG1;
        
        const double slope  = (logG1 - logG0) / (logHV1 - logHV0);
        
        const double HVnom  = exp(dlogG0 / slope + logHV0);
        const double eHVnom = HVnom * sqrt(dlogG1 * dlogG1 * elogG0 * elogG0 +
                                           dlogG0 * dlogG0 * elogG1 * elogG1) / fabs(slope * (logG1 - logG0));
        
        const        double distance = fabs((log(HVnom) - logHV0) / (logHV0 - logHV1));
        static const double maxDist  = 2.0;
        
        if (checkHV(-HVnom) && distance < maxDist) {
          
          TGraphErrors* g0 = getResult();
          
          g0->SetPoint     (0,  HVnom, gainTarget);
          g0->SetPointError(0, eHVnom,        0.0);
          
          return true;
        }
      }
      
      return false;
    }
    
    
    /**
     * Get interpolated high-voltage.
     *
     * \return                    inter-/extrapolated high-voltage\n
     *                            corresponding to target gain value [V]
     */
    double getHV() const
    {
      TGraphErrors* g0 = getResult();
      
      if (g0->GetN() > 0) {
        
        return g0->GetPointX(0);
        
      } else {

        THROW(JNoValue, "JHVInterpolator::getHV(): Missing HV inter-/extrapolation point. Please call JHVInterpolator::interpolateHV() first.");
      }
    }
    
    
    /**
     * Get error estimate on interpolated high-voltage.
     *
     * \return                    error on inter-/extrapolated high-voltage\n
     *                            corresponding to the target gain value [V]
     */
    double getHVError() const
    {
      TGraphErrors* g0 = getResult();
      
      if (g0->GetN() > 0) {

        return g0->GetErrorX(0);
        
      } else {

        THROW(JNoValue, "JHVInterpolator::getHVError(): Missing HV inter-/extrpolation point. Please call JHVInterpolator::interpolateHV() first.");    
      }
    }


    /**
     * Get graph with the input data for the interpolation.
     *
     * \return                   pointer to graph with input data
     */
    TGraphErrors* getData() const
    {
      TList* list = data->GetListOfGraphs();

      if (list != NULL && list->FindObject(HVINTERPOLATOR_DATA) != NULL) {
        
        return (TGraphErrors*)list->FindObject(HVINTERPOLATOR_DATA);
        
      } else if (list == NULL) {

        THROW(JNullPointerException, "JHVInterpolator()::getData(): Emtpy list of graphs.");

      } else {

        THROW(JNullPointerException, "JHVInterpolator()::getData(): No " << HVINTERPOLATOR_DATA << " graph.");
      }
    }


    /**
     * Get graph with the interpolation result.
     *
     * \return                   pointer to graph with the interpolation result
     */
    TGraphErrors* getResult() const
    {
      TList* list = data->GetListOfGraphs();

      if (list != NULL && list->FindObject(HVINTERPOLATOR_RESULT) != NULL) {
        
        return (TGraphErrors*)list->FindObject(HVINTERPOLATOR_RESULT);
        
      } else if (list == NULL) {

        THROW(JNullPointerException, "JHVInterpolator::getResult(): Empty list of graphs.");

      } else {
        
        THROW(JNullPointerException, "JHVInterpolator::getResult(): No " << HVINTERPOLATOR_RESULT << " graph.");
      }
    }

    
    /**
     * Set minimal separation distance for high-voltage.
     *
     * \param  minDist            minimal separation distance for high-voltage
     */
    static void setMinHVDistance(const double minDist)
    {
      dHVmin = minDist;
    }
    

    /**
     * Set valid gain range.
     *
     * \param  range              valid high-voltage range [V]
     */
    static void setHVRange(const JRange<double> range)
    {
      hvRange = range;
    }

    
    /**
     * Set valid gain range.
     *
     * \param  range              valid gain range
     */
    static void setGainRange(const JRange<double> range)
    {
      gainRange = range;
    }
    
    
  private:

    TMultiGraph*          data;           //!< HV-versus-gain data

    static double         dHVmin;         //!< Minimal high-voltage difference between two points [V]
    static JRange<double> hvRange;        //!< Allowed high-voltage range [V]
    static JRange<double> gainRange;      //!< Allowed gain range
  };


  /**
   * Default values.
   */
  double         JHVInterpolator::dHVmin    = 2 * 3.14;
  JRange<double> JHVInterpolator::hvRange   = JRange<double>(-1500, -80);
  JRange<double> JHVInterpolator::gainRange = JRange<double>(FITTOT_GAIN_MIN,
                                                             FITTOT_GAIN_MAX);
}

#endif