WaveformProcessor_service.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
/// \file WaveformProcessor_service.cc
///
///  see Calibrator/func/ADCShapeFitTable
///      Calibrator/art/Calibrator_service
///
/// \author ram2aq@virginia.edu
////////////////////////////////////////////////////////////////////////

#include <iostream>

#include "DDTCore/WaveformProcessor.h"

template<class T> inline T sqr(T x){return x*x;}

namespace novaddt
{

  // Number of points we expect to receive, to find the best peak match with
  const unsigned int kNumFineTimingADCPoints = 4; // Including the DCS origin

  const double kZeroOffsetSamples = -1;

  const double kSamplesPerOffset = 1./64;

  const unsigned int kNumTDCPerSample = 32;


  // this method is no longer used, but could come in handy later
  /*TDC OffsetToTDC(uint32_t const& t0, uint64_t const& high_time, 
                  uint8_t const& offset)
  {
    TDC retval;
    retval.val = high_time<<32 | (t0 + (uint64_t)(offset * 
                                 (kNumTDCPerSample * kSamplesPerOffset)));
    retval.frac = offset%2 * 50;
    return retval;
  }*/

  //......................................................................
  /// \brief Helper function for \a ADCShapeFit inner loop.
  ///
  /// Also used independently to calculate the best fit adc peak height (the
  /// return value).
  inline double GetExpectations(double t0,
                                double riseTime, double fallTime,
                                const int16_t* obs,
                                double* exps)
  {
    for(unsigned int t = 0; t < kNumFineTimingADCPoints; ++t){
      // This is the ideal ASIC formula from SimpleReadout
      exps[t] = (t > t0) ? exp(-(t-t0)/fallTime)*(1-exp(-(t-t0)/riseTime)) : 0;
    }
    // We can calculate the best fit baseline and normalization at this offset
    // analytically (just write the chisq expression, and take derivatives with
    // norm or base), like so:
    double U, V, W, X, Z;
    U = V = W = X = Z = 0;
    for(unsigned int i = 0; i < kNumFineTimingADCPoints; ++i){
      U += exps[i];
      V += 1;
      W += obs[i];
      X += exps[i]*exps[i];
      Z += obs[i]*exps[i];
    }
    const double norm = (W*U-V*Z)/(U*U-X*V);
    const double base = (W*X-Z*U)/(X*V-U*U);

    // Correct all the expectations for the best guess baseline and norm
    for(unsigned int n = 0; n < kNumFineTimingADCPoints; ++n){
      exps[n] = norm*exps[n]+base;
    }
    return norm;
  }

  /*
static float decodeVal8Bit(unsigned int input){
  
  // First get the exponent
  // unsigned int expo = input >> 5;
  
   // Next get the Mantisa
   // unsigned int mantisa = input & 0x1f;
   
   // The value that we return includes an averaging so that
   // we minimize the error (but this converts it to a floating point)
   // return (pow(2, (expo-1))*(mantisa + 32.5) - 0.5);
   
   // This line is the fully optimized version of this
 
   // Pass through values below our floor
   //if(input < 64){return input;}
 
   // This computes the value if its over the pass through floor
   // return (0x1 << ((input >> 5) -1)) *((input & 0x1f) + 32.5) - 0.5);
 
   // Final hand optimized version

  bool sign = input >> 8;
  input = input & 0xff;

  if(sign)
    return (input < 64) ? input : (0x1 << ((input >> 5) -1)) *((input & 0x1f) + 32.5) - 0.5;
  else{
    input = 256 - input;        
    return -((input < 64) ? input : (0x1 << ((input >> 5) -1)) *((input & 0x1f) + 32.5) - 0.5);
  }
}

*/

// Helper function that returns a mask with the Nth bit set
inline unsigned int bitMask(unsigned int n){
  return (n < 32 ) ? (0x1 << n) : 0x0 ;
}


// Helper function that reutrns a mask with the lower N bits set
inline unsigned int nbitMask(unsigned int n){
  unsigned int mask = 0xffffffff;
  mask = mask << (32 - n);
  mask = mask >> (32 - n);
  return mask;
}


// Helper function that returns the highest bit set in an unsigned int
int highestBit(unsigned int input){
  // Now determine the position of the highest non-zero bit
  int highBit_idx=31;
  while(highBit_idx > 0){
    if( bitMask(highBit_idx) & input ){
      // We have found a non-zero bit
      // so we break out
      return highBit_idx;
    }
    --highBit_idx;
  }
  return -1;
}


  unsigned int NBitInt2Float(float _input, int inputBits, int outputBitsMantisa, int outputBitsExpo){
    // First construct the bit mask for the input
    unsigned int InputMask = 0;
    unsigned int nBitIn = 0;
    int highBit = 0;
    int mantisaShift = 0;
    unsigned int expo = 0;
    unsigned int fp_rep=0;
    
    int input = (int)_input;
    bool neg = false;
    if(input < 0){
      input = -1*input;
      neg = true;
    }

    if(input > 4095)
      input = 4095;

    // If the intput value is smaller than the smallest exponent, just return the 
    // value without further encoding
    if(input < (0x1 << (outputBitsMantisa+1) )){return (neg) ? 256-input: 256+input;}
    
    // Create an input mask to mask down the input (i.e. ensure that any random
    // high bits above the specified length are masked down
    InputMask = nbitMask(inputBits);
    // p(InputMask);
    
  // Now mask off to the correct bits
    nBitIn = input & InputMask;
    // p(nBitIn);
    
    // Find the highest non-zero bit in the input word
    highBit = highestBit(nBitIn);
    if(highBit < 0){return 0;}  // Check for the error code

  //printf("High Bit: %d\n",highBit);

  // Comput the shift that is required
  mantisaShift = highBit - outputBitsMantisa;

  // printf("M Shift: %d\n",mantisaShift);
  
  // Shift the input word down
  if(mantisaShift > 0){
    nBitIn = nBitIn >> mantisaShift;
  }
  // p(nBitIn);
  
  // Now mask off the top bit
  nBitIn &= nbitMask(outputBitsMantisa);

  // p(nBitIn);

  // Construct the exponent

  if(mantisaShift > 0){
    expo = mantisaShift + 1;
    expo &= nbitMask(outputBitsExpo);
  }
  //p(expo);

  // Shift the exponent into the correct postion
  expo = expo << outputBitsMantisa;

  //p(expo);

  // Finally combine the mantisa and exponent
  fp_rep = expo | nBitIn;

  //p(fp_rep);

  return (neg) ? 256 - fp_rep: 256 + fp_rep;
}



  //......................................................................
  uint8_t ADCShapeFit(int16_t adc1, int16_t adc2, int16_t adc3,
                      double riseTime, double fallTime,
                      bool& goodTime)
  {
    // The observations, in chronological order
    const int16_t obs[kNumFineTimingADCPoints] = {0, adc1, adc2, adc3};
    double bestchisq = 1e10; // infinity

    // Sane default
    uint8_t bestoffset = -kZeroOffsetSamples/kSamplesPerOffset;
    goodTime = false;

    // Scan through possible offsets of the peak, looking for the best match
    uint8_t offset = 0;
    do{
      // In TDC
      const double t0 = (double(offset)*kSamplesPerOffset+kZeroOffsetSamples);
      // Expectations
      double exps[kNumFineTimingADCPoints];
      const double norm = GetExpectations(t0, riseTime, fallTime, obs, exps);

      if(norm < 0) continue;
      // Assuming equal, uncorrelated, errors on all points. (Including
      // correlations was not found to improve resolution).

      double chisq = 0;
      for(unsigned int i = 0; i < kNumFineTimingADCPoints; ++i)
        chisq += sqr(obs[i]-exps[i]);
      assert(chisq >= 0);

      if(chisq < bestchisq){
        bestchisq = chisq;
        bestoffset = offset;
        goodTime = true;
      }
    } while(++offset != 0);

    if(bestoffset == 0 || bestoffset == 255){
      /*
      std::cerr << "Best time offset is at edge of range. "
                << "If this happens repeatedly there's a "
                << "timing problem. Details: "
                << adc0 << " " << adc1 << " " << adc2 << std::endl;
      */
      bestoffset = -kZeroOffsetSamples/kSamplesPerOffset;
      goodTime = false;
    }

    return bestoffset;
  }


  //---------------------------------------------------------------
  WaveformProcessor::WaveformProcessor(fhicl::ParameterSet const& pset,
                                   art::ActivityRegistry& reg)
     : fSampleTime  (pset.get< double > ("SampleTime"))
     , fRiseTime    (pset.get< double > ("RiseTime"))
     , fFallTime    (pset.get< double > ("FallTime"))
     , fUseMulti    (pset.get< bool >   ("UseMulti"))
     , fUseMulti4ADC(pset.get< bool >   ("UseMulti4ADC"))
     , fDoHighADCFits(pset.get< bool >  ("DoHighADCFits"))
  {
      this->reconfigure(pset);
      reg.sPostBeginRun.watch(this, &WaveformProcessor::postBeginRun);
  }
  //----------------------------------------------------------------
  WaveformProcessor::~WaveformProcessor()
  {
  }

  //----------------------------------------------------------------
  void WaveformProcessor::reconfigure(const fhicl::ParameterSet& pset)
  {
    cet::search_path sp("FW_SEARCH_PATH");
    sp.find_file(pset.get< std::string >("ShapeTableFilename"), fShapeTableFilename);
    std::cout << pset.get< std::string > ("ShapeTableFilename") << std::endl;
    assert(!fShapeTableFilename.empty());
  }

  //---------------------------------------------------------------
  void WaveformProcessor::postBeginRun(art::Run const& run)
  {
    if (fUseMulti && !fTableLoaded) { //load the table if it doesn't already exist
      nTableLoads++;
      std::cout << "Loading table from file for " << nTableLoads << " time." << std::endl;
      TFile f(fShapeTableFilename.c_str());
      assert(!f.IsZombie());

      TTree* tr = (TTree*)f.Get("shape_lookup");
      assert(tr);

      const unsigned int N = tr->GetEntries();
      assert(N == kNumTableEntries);

      uint8_t offset;
      bool good;
      uint16_t adc;
      tr->SetBranchAddress("offset", &offset);
      tr->SetBranchAddress("good", &good);
      if (fUseMulti4ADC) 
        tr->SetBranchAddress("adc", &adc);

      for (unsigned int n=0; n<N; n++){
        tr->GetEntry(n);
        fTable[n] = offset;
        fGood[n] = good;
        if (fUseMulti4ADC)
          fADCTable[n] = adc;
      }

      fTableLoaded = true;
    } // end if (!fTableLoaded)
  }


  void WaveformProcessor::process(uint32_t const& t0, uint64_t const& high_time, 
                           int16_t const& adc1, int16_t const& adc2, 
                           int16_t const& adc3, ADC& adcpeak, 
                           bool& goodTime, TDC& retval)
  {
    uint8_t offset=0;


    int16_t _adc1 = NBitInt2Float(adc1, 12, 5, 3);
    int16_t _adc2 = NBitInt2Float(adc2, 12, 5, 3);
    int16_t _adc3 = NBitInt2Float(adc3, 12, 5, 3);

    if(_adc3 >= kADC3Min && _adc3 < kADC3Max &&
       _adc1 >= kADC1Min && _adc1 < _adc3+kADC1RelMax &&
       _adc2 >= kADC2Min && _adc2 < kADC2Max){
     
      // And this is the offset it'll be at, based on aforementioned geometry
      const int tableIdx = ((kADC2Max-kADC2Min)* // Wedge thickness
                            ((kADC3Min-kADC1Min)*(_adc3-kADC3Min+1) + // Rectangular part
                             kADC1RelMax*(_adc3-kADC3Min)+
                             ((_adc3-kADC3Min)*(_adc3-kADC3Min-1))/2+(_adc1-kADC3Min))+ // Triangular part
                            (_adc2-kADC2Min)); // Distance through wedge

      assert(tableIdx >= 0 && tableIdx < int(kNumTableEntries));
      offset = fTable[tableIdx];
      goodTime = fGood[tableIdx];

    }
    else{
      // It's not in our table, perform the fit if fDoHighADCFits is true
      if(fDoHighADCFits)
        offset = ADCShapeFit(adc1, adc2, adc3, fRiseTime, fFallTime, goodTime);
      else{
        offset = 0;
        goodTime = false;
      }
    }

    // We need to extract the best normalization, so repeat the relevant part
    // of the inner loop at that offset.
    // --- For speed the normalization has been saved in a different table,
    // --- the use of which is set as a fhicl parameter.
    if(goodTime){
      //Perform the normalization:
      // The observations, in chronological order
      //const int16_t obs[kNumFineTimingADCPoints] = {0, adc1, adc2, adc3};
      //const double t0 = (double(offset)*kSamplesPerOffset+kZeroOffsetSamples);
      //double junk[kNumFineTimingADCPoints]; // Expectations
      //const double norm = GetExpectations(t0, fRiseTime, fFallTime, obs, junk);

      // The value the peak of the curve would have before scaling by norm
      // Analytic calculation from the functional form in ReadoutSim
      //static const double peak = pow(fRiseTime/(fRiseTime+fFallTime), 
      //  fRiseTime/fFallTime)-pow(fRiseTime/(fRiseTime+fFallTime), 
      //  1+fRiseTime/fFallTime);
      //adcpeak = norm*peak;

      if (fUseMulti4ADC) { //Extract from lookup table for speed:
        // And this is the offset it'll be at, based on aforementioned geometry
        const int tableIdx = ((kADC2Max-kADC2Min)* // Wedge thickness
                              ((kADC3Min-kADC1Min)*(_adc3-kADC3Min+1) + // Rectangular part
                               kADC1RelMax*(_adc3-kADC3Min)+
                               ((_adc3-kADC3Min)*(_adc3-kADC3Min-1))/2+(_adc1-kADC3Min))+ // Triangular part
                              (_adc2-kADC2Min)); // Distance through wedge



        adcpeak = fADCTable[tableIdx];
      } // else case handled in initmulti()

    }
    else{
      // Unless the fit failed, in which case we're probably safest just doing
      // this.
      adcpeak = ADC(adc3);
    }

    retval.val = (high_time<<32 | (t0 + (uint64_t)(offset * 
                                 (kNumTDCPerSample * kSamplesPerOffset))));
    retval.frac = (uint8_t)(offset%2 * 50);
  }

  void WaveformProcessor::initmulti() { // initialize TDC and ADC members 
                                     //using multipoint
    //check that multipoint can be used
    const uint32_t version = fNano->getHeader()->getVersion();

    static daqdataformats::NanoSliceVersionConvention conv;

    const unsigned int nSamples = conv.getNSamples(version);
    const unsigned int nPretrig = conv.getNPretriggeredSamples(version);

    if ( (nSamples < 3) || // too few samples
         (nSamples <= nPretrig) || //missing peak
         (version == 0) ) //non-multipoint
           { initdcs(); return; }


    TDC tempTDC;
    ADC tempADC;
    uint32_t t0 = fNano->getTimeStamp();
    uint64_t high_time = (uint64_t) fMicro->getHighWord();
    int16_t adc1 = (int16_t)fNano->getValue(1) - (int16_t)fNano->getValue(0);
    int16_t adc2 = (int16_t)fNano->getValue(2) - (int16_t)fNano->getValue(0);
    int16_t adc3 = (int16_t)fNano->getValue(3) - (int16_t)fNano->getValue(0);
    bool goodTime = true;
    process(t0, high_time, adc1, adc2, adc3, tempADC, goodTime, tempTDC);

    if(adc3 < 0) adc3 = 0;                 // Check if adc3 is negative, ADC is unsigned
    if(tempADC.val < 0) tempADC = (ADC)0;   // Check if adc3 is negative, ADC is unsigned
    if(!fUseMulti4ADC) tempADC = (ADC)adc3;

    fTDC = tempTDC;
    fADC = tempADC;
  }

  void WaveformProcessor::setup(uint32_t const& t0, uint64_t const& high_time, 
              int16_t const& adc1, int16_t const& adc2, int16_t const& adc3,
              int16_t const& dcs)
  {
    // For use by RawDigit2DAQHit in the offline.
    // dcs is identical to adc3 for 4-point readout
    if (!fUseMulti) { 
      fTDC.val = high_time<<32 | t0;
      fADC = ADC(dcs);
      return;
    }
    bool goodTime = true;
    TDC tempTDC;
    ADC tempADC;
    process(t0, high_time, adc1, adc2, adc3, tempADC, goodTime, tempTDC);
    if (!goodTime) { // something's gone wrong with the fit, use DCS
      fTDC.val = high_time<<32 | t0;
      fADC = ADC(adc3);
      return;
    }
    else {
      fTDC = tempTDC;
      if (fUseMulti4ADC)
        fADC = tempADC;
      else
        fADC = ADC(adc3);

    }
  }

  void WaveformProcessor::initdcs() { // initialize TDC and ADC members 
                                   //using dcs

    //auto trash = fMicro->getFloatingMicroSlice();
    fTDC.val = ((uint64_t)fMicro->getHighWord() << 32)
                | fNano->getTimeStamp();

    fTDC.frac = 0;
    int16_t adc3 = (int16_t)fNano->getValue(3) - (int16_t)fNano->getValue(0);
    if(adc3 > 0)
      fADC = (ADC)adc3;
    else
      fADC = (ADC)0;

  }

  DEFINE_ART_SERVICE(WaveformProcessor)

} // end namespace novaddt