MDetector/Channel.cc

Go to the documentation of this file.

#include <mdet/Channel.h>

#include <mdet/MHierarchyInfo.h>
//
#include <utl/AugerUnits.h>
//
#include <fwk/RandomEngineRegistry.h>
//
#include <algorithm>
#include <vector>
#include <sstream>
#include <iostream>
#include <string>
//
#include <CLHEP/Random/RandGauss.h>

#include <boost/math/tools/roots.hpp>
#include <boost/math/tools/tuple.hpp>
using boost::math::tools::eps_tolerance; // Binary functor for specified number of bits.
using boost::math::tools::toms748_solve;

namespace {
  
  struct RootComp {
    /*
     * In order to keep things strictly standard, the comparisson
     * functor needs to be homogeneous: so then a dummy key needs to 
     * be built.
     *
     * http://www.sgi.com/tech/stl/lower_bound.html
     *
     * http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2001/n1313.html
     */ 
    template<class C>
    bool operator()(const C& r1, const C& r2) {
      // We go to a template, tough we're kinda assuming our Root type specificaly.
      return r1.fTime < r2.fTime;
    }
  };
  
  bool IsEqual(double a, double b, double eps) 
  { 
    return std::fabs(a - b) <  eps * std::max( std::fabs(a), std::fabs(b));
  };

  class EqualPredicate {
  public:
    EqualPredicate(double e) : fEps(e) {
    }
    bool operator()(double a, double b) {
      return IsEqual(a, b, fEps);
    }
  private:
    double fEps;
  };


  template <class T, class F>
  class MyFunctionObject {
  public:
    MyFunctionObject(T const& thr, mdet::Channel::Discriminator::Proxy<F> const& func) : fThr(thr), fFctor( func ) {} 
  
    T operator()(T const& t ) {
      return fFctor(t) - fThr;
    }

  private:
    T  fThr;
    mdet::Channel::Discriminator::Proxy<F> fFctor;
  };



  template<class BackInsertionSequence, class Thr,class Fctor>
  void RootFinder( MyFunctionObject<Thr,Fctor>& func, BackInsertionSequence& roots, double s, double e, double f_min, double f_max, int nest, double relTol) { 
  //void RootFinder( mdet::Channel::Discriminator::Proxy<Fctor> & func, BackInsertionSequence& roots, double s, double e, double f_min, double f_max, int nest, double relTol) { 

    using namespace std;  
    using namespace boost::math;
    //double guess  = (s+e)/2;
  
    int digits = std::numeric_limits<double>::digits; // Maximum possible binary digits accuracy for type double.
    // digits used to control how accurate to try to make the result.
    int get_digits = (digits * 3) /4; // (3/4) of maximum accuracy.
  
    const boost::uintmax_t maxit = 100;
    boost::uintmax_t it = maxit; // Initally our chosen max iterations, but updated with actual.
  
    eps_tolerance<double> tol(get_digits);
  
    //Solve 
    std::pair<double, double> bracket_root = toms748_solve(func, s, e, f_min, f_max, tol, it);
  
    double r = bracket_root.first + (bracket_root.second - bracket_root.first)/2;  // Midway between brackets.
    
    if (it < maxit ) {  
      //cout << "Converged after " << it << " iterations " << endl;
      //cout << "(" << nest << ") start=" << s << " end=" << e << " guess=" << guess << " root=" << r << endl;
     /*
      * Extended check including boundaries; and using IsEqual
      */ 
      if ( (s <= r || IsEqual(r, s, relTol)) && (r <= e || IsEqual(r, e, relTol)) ) {
        //cout << "ADDING ONE ROOT\n";
        roots.push_back( r ) ;
        if ( !IsEqual(r, s, relTol) && !IsEqual(r, e, relTol ) ) {
           RootFinder( func, roots, s, r, func(s), func(r), nest+1, relTol);
           RootFinder( func, roots, r, e, func(r), func(r), nest+1, relTol);
        }//end boundary check
      }//end if in interval
  
    }//end if it < maxit
    
    return;
  } 
}

namespace mdet {

  const char* const 
  Channel::kComponentName = MHierarchyInfo::kComponentsNames[7];

  const char* const  
  Channel::kComponentId   = MHierarchyInfo::kComponentsIds[7];

  Channel::Channel
    (int cId, const det::VManager::IndexMap& parentMap, const FrontEnd& parent) :
      MDetectorComponent<Channel>::Type(cId, parentMap),
      fFrontEnd(parent)
  { 
    Register(fInvertingInputResistance);
    Register(fFeedbackResistance);
    Register(fDCGain);
    Register(fLowCutoffFrequency);
    Register(fHighCutoffFrequency);
    Register(fThreshold);
    Register(fResponseTime);
    Register(fDiscriminatorLowLevel);
    Register(fDiscriminatorHiLevel);
    Register(fSlewRate);
    Register(fAbsoluteError);
    Register(fInitialIntervalLength);
    Register(fIterationsNumber);
    Register(fMaxNumberOfErrors);
    Register(fSignalShiftMean);
    Register(fSignalShiftStdDev);
  }

  double 
  Channel::GetInvertingInputResistance()
    const
  {
    return GetData(fInvertingInputResistance, "invertingInputResistance");
  }

  double 
  Channel::GetFeedbackResistance()
    const 
  {
    return GetData(fFeedbackResistance, "feedbackResistance");
  }
  
  double 
  Channel::GetDCGain()
    const
  {
    return GetData(fDCGain, "DCGain");

  }
  
  double 
  Channel::GetLowCutoffFrequency()
    const
  {
    return GetData(fLowCutoffFrequency, "lowCutoffFrequency");
  }

  double 
  Channel::GetHighCutoffFrequency()
    const
  {
    return GetData(fHighCutoffFrequency, "highCutoffFrequency");
  }

  double 
  Channel::GetThreshold()
    const
  {
    return GetData(fThreshold, "threshold");
  }

  double 
  Channel::GetDiscriminatorLowLevel() 
    const
  {
    return GetData(fDiscriminatorLowLevel, "discriminatorLowLevel");
  }
  
  double 
  Channel::GetDiscriminatorHiLevel() 
    const 
  {
    return GetData(fDiscriminatorHiLevel, "discriminatorHiLevel");
  }
    
  double 
  Channel::GetSlewRate() 
    const
  {
    return GetData(fSlewRate, "slewRate");
  }

  double 
  Channel::GetResponseTime()
    const
  {
    return GetData(fResponseTime, "responseTime");
  }

  double 
  Channel::GetAbsoluteError()
    const
  {
    return GetData(fAbsoluteError, "absoluteError");
  }

  double
  Channel::GetInitialIntervalLength()
    const
  {
    return GetData(fInitialIntervalLength, "initialIntervalLength");
  }

  unsigned int
  Channel::GetIterationsNumber()
    const
  {
    return GetData(fIterationsNumber, "iterationsNumber");
  }

  unsigned int
  Channel::GetMaxNumberOfErrors()
    const
  {
    return GetData(fMaxNumberOfErrors, "maxNumberOfErrors");
  }


  double 
  Channel::GetSignalShiftMean() 
    const
  {
    return  GetData(fSignalShiftMean, "signalShiftMean");
  }
    
  double 
  Channel::GetSignalShiftStdDev() 
    const
  {
    return  GetData(fSignalShiftStdDev, "signalShiftStdDev");

  }

  double
  Channel::ComputeSignalShift()
    const
  {
    utl::RandomEngine& rnd = fwk::RandomEngineRegistry::GetInstance().Get(fwk::RandomEngineRegistry::eDetector);
    /*
     * 1) Initially, as in other cases, the return value was forced to be positive (with a loop until that), but
     * for this cases (a shift) it doesn't matter if the number (likely or not) is negative.
     * Granted, a negative shift will mean something like that the output came before the signal
     * entering the device. Despite that being impossible physically, it's something quite different than
     * having a negative amplitude in the pulses.
     *
     * 2) From CLHEP's RandGauss.icc file:
     *
     * inline double RandGauss::shoot(HepRandomEngine* anEngine, double mean, double stdDev) {
     *   return shoot(anEngine)*stdDev + mean;
     * }
     *
     * so we can freely set the standard deviation to zero and generate safely a degenerated
     * spike distribution, with only the mean as possible outcome.
     */  
    return CLHEP::RandGauss::shoot(& rnd.GetEngine(), GetSignalShiftMean(), GetSignalShiftStdDev());
  }

  /*
   *********************************
   ***** Amplification stage. ******
   *********************************
   */ 

  std::complex<double> 
  Channel::ComputeFrequencyFactor(double freq, double freqLimit) 
    const
  {
    const std::complex<double> _1(1, 0); // real unit.
    const std::complex<double> i(0, 1);  // imaginary unit
    return _1 + i * (freq / freqLimit);
  }

  std::complex<double> 
  Channel::ComputeTransfer(double freq) 
    const
  {
    const std::complex<double> fDenominator = ComputeFrequencyFactor(freq, GetLowCutoffFrequency()) * 
                                              ComputeFrequencyFactor(freq, GetHighCutoffFrequency());
    const std::complex<double> openLoopGain = GetDCGain() / fDenominator;
    const std::complex<double> rDenominator = GetInvertingInputResistance() + 
                                              GetFeedbackResistance()       +
                                              openLoopGain * GetInvertingInputResistance();
    // Note the sign that defines the inverting characteristic.
    // This is related to the fact that the SPE pulses are negative peaks. 
    return - openLoopGain * GetInvertingInputResistance() * GetFeedbackResistance() / rDenominator;
  }

  /*
   *********************************
   ***** Discrimination stage. *****
   *********************************
   */ 

  Channel::Discriminator::Callable::~Callable()
  {
    // At first this destructor was made protected and non-virtual:
    //
    // Being an interface we protect (and non-virtualize)
    // the destructor.
    // http://www.gotw.ca/publications/mill18.htm
    //
    // One may get
    // ./mdet/Channel.h:59: warning: 
    // 'class mdet::Channel::Discriminator::Callable' has virtual functions but non-virtual destructor
    // but that's a phony alarm.
    // See for instance,
    // http://stackoverflow.com/questions/127426/gnu-compiler-warning-class-has-virtual-functions-but-non-virtual-destructor
    //
    // At last, that's an error, because we'll be deleting through it: so needs to be virtual, so
    // as to have the subclass' destructor be properly called.
    // Also, being protected you may get 
    // MDetector/Channel.cc: In destructor 'std::auto_ptr<_Tp>::~auto_ptr() [with _Tp = mdet::Channel::Discriminator::Callable]':
    // MDetector/Channel.cc:174:   instantiated from here
    // MDetector/Channel.cc:128: error: 'mdet::Channel::Discriminator::Callable::~Callable()' is protected
  }

  bool 
  Channel::Discriminator::OverThreshold(double t) 
    const
  {
    return ( (*fInputPulse)(t) - fChannel.GetThreshold() ) > 0;
  }

  double 
  Channel::Discriminator::ApplySaturation(double v) 
    const
  {
    if (v > fChannel.GetDiscriminatorHiLevel())
      return fChannel.GetDiscriminatorHiLevel();
    else if (v < fChannel.GetDiscriminatorLowLevel())
      return  fChannel.GetDiscriminatorLowLevel();
    else
      return v;
  }

  double
  Channel::Discriminator::ComputeOutput(double signReferenceTime, double deltaTime, double prevOutput)
    const
  {
    const bool over = OverThreshold(signReferenceTime);
    if (deltaTime > fTransitionTime) {
      if (over) 
        return fChannel.GetDiscriminatorHiLevel(); 
      else  
        return fChannel.GetDiscriminatorLowLevel(); 
    } else {
      double deltaVoltage = deltaTime * fChannel.GetSlewRate();
      double o = prevOutput;
      if (over)
        o += deltaVoltage;
      else
        o -= deltaVoltage;
      return ApplySaturation(o);
    }
  }
    
  Channel::Discriminator::AtThresholdConstIterator 
  Channel::Discriminator::AtThresholdBegin() 
    const
  {
    return fRawRoots.begin();
  }

  Channel::Discriminator::AtThresholdConstIterator 
  Channel::Discriminator::AtThresholdEnd() 
    const
  {
    return fRawRoots.end();
  }

  void 
  Channel::Discriminator::Init()
  {
    /*
     * Divide the interval and look up for root in each interval
     * individually.
     */ 
    const double absTol = fChannel.GetAbsoluteError();
    const double relTol = absTol / std::max( std::fabs(fBeginTime), std::fabs(fEndTime) );
    const double gap = fChannel.GetInitialIntervalLength();
    double currStart = fBeginTime;
    double currEnd   = currStart + gap;
    double shift     = 0.1*gap;//To overlap intervals

    Proxy<Callable> f( (*fInputPulse) );
    double thr = fChannel.GetThreshold();
    MyFunctionObject<double,Proxy<Callable> > functor(thr,f);

    //Solve in each interval
    for ( ; currEnd <= fEndTime; currStart += gap, currEnd = currStart+gap ) {

       double f_min = functor( currStart-shift ) ;
       double f_max = functor( currEnd  ) ;

       try {
         RootFinder( functor, fRawRoots, currStart-shift, currEnd, f_min, f_max, 0, relTol);
       } catch(const std::exception& e) {
         // Always useful to include try & catch blocks because default policies 
         // are to throw exceptions on arguments that cause errors like underflow, overflow. 
         // Lacking try & catch blocks, the program will abort without a message below,
         // which may give some helpful clues as to the cause of the exception.
         //if( verbose > 1 ) std::cout << "\nMessage from thrown exception was:\n   " << e.what() << std::endl;
       }//end try
    }//end for intervals

    /*
     * Sort: needed for the following duplicate removal and also because 
     * the operator() assumes the roots to be sorted.
     */
    std::sort(fRawRoots.begin(), fRawRoots.end());
    // Given that we may continue subdiving the intervals, some root can be found several times.
    RawRootsContainer::iterator last = std::unique_copy(
      fRawRoots.begin(), 
      fRawRoots.end(), 
      fRawRoots.begin(), 
      EqualPredicate(relTol) // create a functor for comparison
    );
    // Get rid of the remaining non-important elements.
    fRawRoots.erase(last, fRawRoots.end());
    // Preallocate:
    fRoots.reserve(fRawRoots.size());
    // Now fill the field with Root (time + output).
    double prevTime   = fBeginTime;
    // Take the low as the starting level (see operator()(double t)).
    double prevOutput = fChannel.GetDiscriminatorLowLevel();
    // (These 2 last values could be stored in our structure Root, 
    // but I feel that to be misleading).
    /*
     * Computing the roots is not the whole history, since the device has a finite-time
     * response we have to do some calculations so as to know if the roots are separated
     * enough (or not) for the device output signal change to the final desired value.
     */ 
    for (RawRootsContainer::const_iterator i = fRawRoots.begin(), e = fRawRoots.end(); i != e; ++i) {
      Root r;
      r.fTime               = *i;
      const double delta    = r.fTime - prevTime;
      // See the sign half way to the previous
      const double halfTime = ( prevTime + r.fTime ) / 2;
      r.fOutput = ComputeOutput(halfTime, delta, prevOutput);
      fRoots.push_back(r);
      // Update for the next loop.
      prevTime   = r.fTime;
      prevOutput = r.fOutput;
    }
  }

  Channel::Discriminator::Discriminator(const Discriminator& d) :
    fChannel(d.fChannel),
    fInputPulse(d.fInputPulse->Clone()), // This is the important one, the rest are just a copy-construction.
    fTransitionTime(d.fTransitionTime),
    fBeginTime(d.fBeginTime),
    fEndTime(d.fEndTime),
    fRoots(d.fRoots)
  {
  }

  double
  Channel::Discriminator::operator()(double t) 
    const
  {
    /*
     * Note that the real device may have a finite response time
     * (in the sense of detecting the signal over the threshold).
     * Given this, the real device won't start the transition
     * as soon as the signal crosses the threshold.
     * So far this hasn't been modeled. The changes needed may
     * need in these lines:
     * - Call this finite time 'response'
     * - While calculating the roots (in the init method) perform
     *   an analysis in the neighborhood: [root, root + response]
     * - If the time over / under the threshold (depending of the sign
     *   of the slope in the root) is not enough, the root may be
     *   discarded.
     * - If the time is enough, the root maybe taken as root + response;
     *   so as to start the transition when the required time is over.
     * - This shift seems to be enough to take into account the response
     *   time, and the algorithm may stay the same.
     */ 
    /*
     * In general, for this method we take that
     * initially (say in "turned-off" state)
     * the output is at the low level.
     */ 
    // Special case without roots
    if (fRoots.empty()) {
      double delta = t - fBeginTime;
      // 1st arg: Take the reference time for sign as the beginTime (any time will do).
      // 2nd arg: The delta, since there's no root, against the starting time.
      // 3rd arg: For the previous output take the low state (see comment at the beginning).
      return ComputeOutput(fBeginTime, delta, fChannel.GetDiscriminatorLowLevel());
    }
    RootComp comp; // Comparator.
    Root key;      // Dummy object to use for searching: only the time is relevant.
    key.fTime = t;
    /*
     * Find the lower bound for the given time
     * "The second version of lower_bound returns the furthermost iterator i in [first, last) 
     * such that, for every iterator j in [first, i), comp(*j, value) is true."
     * See http://www.sgi.com/tech/stl/lower_bound.html.
     *
     * So the iterator points to the first root that comes after the given time.
     */
    RootIterator i = std::lower_bound(fRoots.begin(), fRoots.end(), key, comp);
    // Difference time against the time of a root (or the start).
    double deltaTime;
    // Reference output taken as starting point to compute the value at the desired time.
    double prevOutput;
    if (i == fRoots.end()) {
      // The time comes after the last root.
      deltaTime = t - fRoots.back().fTime;
      prevOutput = fRoots.back().fOutput;
    } else if (i == fRoots.begin()) {
      // The time comes before the first root: so take difference against the starting time.
      deltaTime = t - fBeginTime;
      // Take as starting in the low level (see initial comment).
      prevOutput = fChannel.GetDiscriminatorLowLevel(); 
    } else {
      // See the previous root, that's the 
      // inmediately to the left of the given time.
      RootIterator prev = i - 1;
      deltaTime  = t - prev->fTime;
      prevOutput = prev->fOutput;
    }
    // The sign is checked for the given time (first argument).
    return ComputeOutput(t, deltaTime, prevOutput);
  }


  template<class Iterator>
  double 
  Channel::Discriminator::ComputeBorderTime(Iterator beg, Iterator end, double time) 
    const
  {
    /*
     * Note that strictly speaking the roots may only touch zero, ie they are an extrema.
     * So maybe none of the conditions hold, and that would mean that we wouldn't initialize
     * this variable in the following: so give it a value.
     * Important: note that in case we return this zero, both calls to this functions
     * (using in-order or reversed iterators) would return that zero. So there's
     * no problem, that's exactly what we need: if never the conditions holds that
     * would meant that all the roots are maxima, so actually the function never crosses
     * (tough it touches) zero. So, when calculating the difference between both calls
     * to this function we would obtain 0, which is correct.
     * Note that if all the roots are minima, then at least the first condition over
     * time (which isn't a root in both calls to this function) would be satisfied, and
     * so the return value initilized with that time.
     */
    double borderTime = 0;
    if (OverThreshold(time))
      borderTime = time;
    else {
      /*
       * Note that if there are no roots, then zero is returned: see comment at the beginning.
       * Also, if there are no roots, and the signal is over, then we would have entered
       * the 'then' condition of this if-then-else.
       */   
      double prevTime = time;
      for (Iterator i = beg; i != end; ++i) {
        double t    = i->fTime;
        // Check the sign halfway between these two representative times: 't' is always
        // a root and 'prevTime' is always a root but in the first loop.
        double test = (prevTime + t) / 2;
        if (OverThreshold(test)) {
          // If it's over the threshold, then it started doing so in the previous root.
          borderTime = prevTime;
          break;
        }
        // Keep the current root for the following loop.
        prevTime = t;
      }
    }
    return borderTime;
  }

  double
  Channel::Discriminator::GetOverThresholdTimeSpan() 
    const
  {
    double earliestTime = ComputeBorderTime(fRoots.begin(), fRoots.end(), fBeginTime);
    // Note that we're using here that we have a ReversibleContainer, which is the case.
    double latestTime   = ComputeBorderTime(fRoots.rbegin(), fRoots.rend(), fEndTime);
    return latestTime - earliestTime;
  }

  double 
  Channel::Discriminator::operator()(double *x, double * /*p*/)
    const
  {
    return (*fInputPulse)(x[0]) - fChannel.GetThreshold();
  }

}