JHVGainGraph.hh

Go to the documentation of this file.

#ifndef __JCALIBRATE_JHVGAINGRAPH__
#define __JCALIBRATE_JHVGAINGRAPH__

#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::JValueOutOfRange;
  
  static const std::string HVGAINGRAPH_SUFFIX = ".HVxG";
  
  static const char*       HVGAINGRAPH_DATA   = "HVGAINGRAPH_DATA";
  static const char*       HVGAINGRAPH_RESULT = "HVGAINGRAPH_RESULT";
  
  /**
   * Auxiliary data structure to store high-voltage versus gain data and interpolate the nominal high-voltage.
   */
  struct JHVGainGraph :
    public TMultiGraph
  {
    /**
     * Constructor.
     *
     * \param  name               name
     * \param  logScale           set logarithmic scaling
     * \param  base               base for logarithmic scaling
     */
    JHVGainGraph(const char*  name     = "HVGainGraph",
                 const bool   logScale = false,
                 const double base     = 10.0) :
      loglog(logScale), base(base)
    {
      this->SetName(name);

      TGraphErrors* g0 = getResult();
      TGraphErrors* g1 = getData();

      g0->SetName(HVGAINGRAPH_RESULT);
      g0->SetLineColor  (kRed);
      g0->SetMarkerColor(kRed);
      g0->SetMarkerStyle(kFullDotSmall);

      g1->SetName(HVGAINGRAPH_DATA);
      g1->SetMarkerStyle(kFullDotSmall);
    }


    /**
     * Copy constructor.
     *
     * \param  object             TMultiGraph object
     * \param  logScale           set logarithmic scale
     * \param  base               base for logarithmic scaling
     */
    JHVGainGraph(const TMultiGraph&   object,
                 const bool           logScale,
                 const double         base = 10.0) :
      loglog(logScale), base(base)
    {
      TList* graphs = object.GetListOfGraphs();

      if (graphs != NULL) {

        TGraphErrors* g0 = (TGraphErrors*)graphs->FindObject(HVGAINGRAPH_RESULT);
        TGraphErrors* g1 = (TGraphErrors*)graphs->FindObject(HVGAINGRAPH_DATA);

        if (g0 != NULL) {
          this->Add((TGraphErrors*)g0->Clone());
        }

        if (g1 != NULL) {
          this->Add((TGraphErrors*)g1->Clone());
        }
      }
    }
    
    
    /**
     * 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();
      
      if (loglog) {

        g1->SetPoint     (N, log(fabs(HV)) / log(base), log(gain) / log(base));
        g1->SetPointError(N, 0.0,                       gainError / log(base) / gain);
        
      } else {

        g1->SetPoint     (N, 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();

      if (loglog) {
        
        g1->SetPoint     (i, log(fabs(HV)) / log(base), log(gain) / log(base));
        g1->SetPointError(i, 0.0,                       gainError / log(base) / gain);
        
      } else {

        g1->SetPoint     (i, 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)
    {
      TGraphErrors* g1 = getData();
      
      if ( (g1 != NULL) && ((i >= 0 && i < g1->GetN()) &&
                            (j >= 0 && i < g1->GetN())) ) {

        if (loglog) {
          
          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 ((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)
    {
      TGraphErrors* g1 = getData();
      
      if ((i >= 0 && i < g1->GetN()) &&
          (j >= 0 && i < g1->GetN())) {
        
        double HV_i, HV_j, gain_i, gain_j;
        
        if (loglog) {

          HV_i   = -pow(base, g1->GetPointX(i));
          HV_j   = -pow(base, g1->GetPointX(j));
          gain_i =  pow(base, g1->GetPointY(i));
          gain_j =  pow(base, g1->GetPointY(j));

        } else {

          HV_i   = g1->GetPointX(i);
          HV_j   = g1->GetPointX(j);
          gain_i = g1->GetPointY(i);
          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 interpolate(const double gainTarget)
    {
      TGraphErrors* g1 = getData();
      
      if (g1->GetN() < 2 || !checkGain(gainTarget)) {
        return false;
      }
      
      const bool isLogLog = loglog;
      
      if (!isLogLog) {
        setLogLog();
      }
      
      // Search for valid inter-/extrapolation points
      
      Int_t i = 0;
      Int_t j = 1;
      
      const double logGainTarget = log(gainTarget) / log(base);

      for (Int_t k = 1; k < g1->GetN(); ++k) {

        const double dlogG_i = fabs(g1->GetPointY(i) - logGainTarget);
        const double dlogG_j = fabs(g1->GetPointY(j) - logGainTarget);  
        const double dlogG_k = fabs(g1->GetPointY(k) - logGainTarget);
        
        if         (dlogG_k < dlogG_i) {
          
          j = (areValid(i, k) ? i : j);
          i = k;
          
        } else if ((dlogG_k < dlogG_j || !areValid(i, j)) &&
                   areValid(i, k)) {
          j = k;
        }
      }

      // Inter-/Extrapolate high-voltage corresponding to given gain

      bool valid = areValid(i, j);
      
      if (valid) {
        
        const double logHV0    = g1->GetPointX(i);
        const double logHV1    = g1->GetPointX(j);
        
        const double logG0     = g1->GetPointY(i);
        const double logG1     = g1->GetPointY(j);
        const double elogG0    = g1->GetErrorY(i);
        const double elogG1    = g1->GetErrorY(j);
        
        const double dlogG0    = logGainTarget - logG0;
        const double dlogG1    = logGainTarget - logG1;
        
        const double slope     = (logG1 - logG0) / (logHV1 - logHV0);
        
        const double logHVnom  = dlogG0 / slope + logHV0;
        const double elogHVnom = sqrt(dlogG1 * dlogG1 * elogG0 * elogG0 +
                                      dlogG0 * dlogG0 * elogG1 * elogG1) / fabs(slope * (logG1 - logG0));
        
        const double distance  = fabs((logHVnom - logHV0) / (logHV0 - logHV1));

        static const double maxDist = 2.0;
        
        valid = (logHVnom > 0.0 && checkHV(-pow(base, logHVnom)) && distance < maxDist);
        
        if (valid) {
          
          TGraphErrors* g0 = getResult();
          
          g0->SetPoint     (0, logHVnom,  logGainTarget);
          g0->SetPointError(0, elogHVnom, 0.0);
        }
      }
      
      if (!isLogLog) {
        unsetLogLog();
      }
      
      return valid;
    }


    /**
     * Get high-voltage corresponding to given target gain value.
     *
     * \param  gainTarget         target gain value for inter-/extrapolation
     * \return                    inter-/extrapolated high-voltage\n
     *                            corresponding to target gain value [V]
     */
    double getTargetHV(const double gainTarget)
    {
      static const double precision = 1e-3;

      const double target = (loglog ? log(gainTarget) / log(base) : gainTarget);
      
      TGraphErrors* g0 = getResult();
      
      if ((g0->GetN() > 0 && fabs(g0->GetPointY(0) - target < precision)) || interpolate(gainTarget)) {
        
        return (loglog ? -pow(base, g0->GetPointX(0)) : g0->GetPointX(0));
        
      } else {

        THROW(JValueOutOfRange, "JHVGainGraph::getTargetError(): No valid inter-/extrapolation candidate found for given target gain.");
      }
    }
    
    
    /**
     * Get error estimate on high-voltage corresponding to given target gain value.
     *
     * \param  gainTarget         target gain value for inter-/extrapolation
     * \return                    error on inter-/extrapolated high-voltage\n
     *                            corresponding to the target gain value [V]
     */
    double getTargetHVError(const double gainTarget)
    {
      static const double precision = 1e-3;

      const double target = (loglog ? log(gainTarget) / log(base) : gainTarget);
      
      TGraphErrors* g0 = getResult();
      
      if ((g0->GetN() > 0 && fabs(g0->GetPointY(0) - target) < precision) || interpolate(gainTarget)) {

        return (loglog ? g0->GetErrorX(0) * pow(base, g0->GetPointX(0)) * log(base) : g0->GetErrorX(0));
        
      } else {

        THROW(JValueOutOfRange, "JHVGainGraph::getTargetError(): No valid inter-/extrapolation candidate found for given target gain.");        
      }
    }

    
    /**
     * Set log-log scale.
     */
    void setLogLog()
    {
      if (!loglog) {
        
        TIter next(this->GetListOfGraphs());
      
        while (TGraphErrors* graph = (TGraphErrors*)next()) {

          if (graph != NULL) {
            
            for (Int_t i = 0; i < graph->GetN(); ++i) {
       
              Double_t absLogHV      = log(fabs(graph->GetPointX(i))) / log(base);
              Double_t logGain       = log(graph->GetPointY(i))       / log(base);

              Double_t absLogHVerror = graph->GetErrorX(i) / fabs(graph->GetPointX(i)) / log(base);
              Double_t logGainError  = graph->GetErrorY(i) / fabs(graph->GetPointY(i)) / log(base);
              
              graph->SetPoint     (i, absLogHV,      logGain);
              graph->SetPointError(i, absLogHVerror, logGainError);
            }
          }
        }

        loglog = true;
      }
    }


    /**
     * Unset log-log scale.
     */
    void unsetLogLog()
    {
      if (loglog) {

        TIter next(this->GetListOfGraphs());
        
        while (TGraphErrors* graph = (TGraphErrors*)next()) {

          if (graph != NULL) {
            
            for (Int_t i = 0; i < graph->GetN(); ++i) {
       
              Double_t HV            = -pow(base, graph->GetPointX(i));
              Double_t gain          =  pow(base, graph->GetPointY(i));

              Double_t HVerror       = graph->GetErrorX(i) * -HV  * log(base);
              Double_t gainError     = graph->GetErrorY(i) * gain * log(base);
              
              graph->SetPoint     (i, HV,      gain);
              graph->SetPointError(i, HVerror, gainError);
            }
          }
        }

        loglog = false;
      } 
    }


    /**
     * Set logarithmic base for log-log scaling.
     *
     * \param  base               logarithmic base for scaling
     */
    void setLogBase(const double base)
    {
      if (loglog) {
        
        TIter next(this->GetListOfGraphs());
      
        while (TGraphErrors* graph = (TGraphErrors*)next()) {

          if (graph != NULL) {
            
            for (Int_t i = 0; i < graph->GetN(); ++i) {
       
              Double_t absLogHV      = graph->GetPointX(i) * log(this->base) / log(base);
              Double_t logGain       = graph->GetPointY(i) * log(this->base) / log(base);
            
              Double_t abslogHVerror = graph->GetErrorX(i) * log(this->base) / log(base);
              Double_t logGainError  = graph->GetErrorY(i) * log(this->base) / log(base);
            
              graph->SetPoint     (i, absLogHV,      logGain);
              graph->SetPointError(i, abslogHVerror, logGainError);
            }
          }
        }
      }
      
      this->base = base;
    }


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

      if (list != NULL && list->FindObject(HVGAINGRAPH_DATA) != NULL) {
        
        g1 = (TGraphErrors*)list->FindObject(HVGAINGRAPH_DATA);
        
      } else {
        
        g1 = new TGraphErrors();
        
        g1->SetName(HVGAINGRAPH_DATA);
        
        g1->Set(0);
        
        this->Add(g1);
      }

      return g1;
    }


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

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

        g0 = new TGraphErrors();

        g0->SetName(HVGAINGRAPH_RESULT);

        g0->Set(0);
        
        this->Add(g0);
      }

      return g0;
    }


    /**
     * Clone this object.
     *
     * \param  newname            new name
     */
    JHVGainGraph* Clone(const char* newname = "") const
    {
      JHVGainGraph* clone = new JHVGainGraph(*this, loglog, base);
      
      clone->SetName(newname);
      
      return clone;
    }
    
    
    /**
     * 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:

    bool   loglog;                        //!< Log-log scale toggle
    double base;                          //!< Base for logarithmic scaling
    
    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         JHVGainGraph::dHVmin    = 2 * 3.14;
  JRange<double> JHVGainGraph::hvRange   = JRange<double>(-1500, -80);
  JRange<double> JHVGainGraph::gainRange = JRange<double>(FITTOT_GAIN_MIN,
                                                          FITTOT_GAIN_MAX);
}

#endif