ROsFadcReadoutTool.cc

Go to the documentation of this file.

#include "ROsFadcReadoutTool.h"

#include "Event/ElecFeeCrate.h"

#include "Event/ReadoutPmtCrate.h"

#include "Conventions/Electronics.h"

#include <vector>
#include <map>


ROsFadcReadoutTool::ROsFadcReadoutTool(const std::string& type,
                                     const std::string& name, 
                                     const IInterface* parent)
  : GaudiTool(type,name,parent)
{
  declareInterface< IROsFadcReadoutTool >(this) ;
  declareProperty("SimFrequency", m_simFrequency = DayaBay::TdcFrequencyHz,"Simulation Frequency");
  declareProperty("ESumFrequency", m_eSumFrequency = DayaBay::EsumFrequencyHz,"Esum Frequency");
  declareProperty("FadcFrequency", m_fadcFrequency = DayaBay::FadcFrequencyHz,"Fadc Frequency");
  declareProperty("ADCOffset", m_adcOffset = 5,"Fadc Frequency");
  declareProperty("ADCRange", m_adcRange = 2.2e-6,"Fadc Range");
  declareProperty("ADCBits", m_adcBits = 8,"Fadc bits");
}

ROsFadcReadoutTool::~ROsFadcReadoutTool(){}

StatusCode ROsFadcReadoutTool::initialize()
{
  return StatusCode::SUCCESS;
}

StatusCode ROsFadcReadoutTool::finalize()
{
  return StatusCode::SUCCESS;
}

StatusCode ROsFadcReadoutTool::readoutFADC(const DayaBay::ElecFeeCrate *elecCrate,
                                           DayaBay::ReadoutPmtCrate *roCrate,
                                           unsigned int startCycle, 
                                           unsigned int stopCycle)
{  
  info() << "Doing FADC readout" << endreq;
  const DayaBay::ElecFeeCrate::AnalogMap& boardMap = elecCrate->esum();
  DayaBay::ReadoutPmtCrate::FadcReadouts& fadcMap = roCrate->fadcReadout();
  ROsFadcReadoutTool::analogFADCmap analogFADCMap;
 
  debug() << "Filling fadc map with esum map of size " << boardMap.size() << endreq;

  if ( fillFadcMap(boardMap, analogFADCMap).isFailure() )
    {
      fatal() << "fillFadcMap failed" << endreq;
      return StatusCode::FAILURE;
    }

  debug() << "Size of filled FADC map " << analogFADCMap.size() << endreq;
  
  // loop over FADC vectors and upsample them
  ROsFadcReadoutTool::analogFADCmap::iterator boardIt = analogFADCMap.begin();
  for( ; boardIt != analogFADCMap.end(); ++boardIt)
    {
      DayaBay::FadcChannelId *outputId = new DayaBay::FadcChannelId(boardIt->first);
      //DayaBay::FadcChannelId outputId(boardIt->first);  

      verbose() << "Processing: " << outputId->connector() << " w/ id = "
                << *outputId << endreq;
   
      DayaBay::AnalogSignal& fadcVector = boardIt->second;
      DayaBay::AnalogSignal interpFadcVec;
      
      DayaBay::DigitalSignal& digFadcVec = fadcMap[*outputId];
      
      debug() << "Upsampling vec of size " << fadcVector.size() << endreq;

      if ( upsample(fadcVector, interpFadcVec, m_simFrequency, m_fadcFrequency, startCycle, stopCycle).isFailure() )
        {
          fatal() << "upsample failed" << endreq;
          return StatusCode::FAILURE;
        }

      debug() << "Upsampled vec has size " << interpFadcVec.size() << endreq;      
      debug() << "Digitizing vec of size " << interpFadcVec.size() << endreq;
      verbose() << "interpVector " << interpFadcVec << endreq;

      if ( digitize(interpFadcVec, fadcMap[*outputId], m_adcRange, m_adcBits, m_adcOffset).isFailure() )
        {
          fatal() << "digitize failed" << endreq;
          return StatusCode::FAILURE;
        }
      //delete outputId;
      debug() << "Digitized vec has size " << digFadcVec.size() << endreq;

    }

    return StatusCode::SUCCESS;
}

// Helpers for FADC Readout

StatusCode ROsFadcReadoutTool::fillFadcMap(const DayaBay::ElecFeeCrate::AnalogMap& input,
                                           analogFADCmap& output)
{
  // Fills a map with the FADC channel IDs and ESum vectors
    debug() << "Filling FADC map " << endreq;
    DayaBay::ElecFeeCrate::AnalogMap::const_iterator boardIt = input.begin();
    for (;boardIt != input.end(); ++boardIt)
    {
        DayaBay::FeeChannelId boardId = boardIt->first;
        const DayaBay::AnalogSignal& boardEnergy = boardIt->second;
        //int boardNum = boardId.board();
        unsigned int boardIdx = boardEnergy.size();
         
        verbose() << "boardIdx = " <<  boardEnergy.size() << endreq;
 
        // Get the FADC input and output channel IDs
        DayaBay::FadcChannelId fadcChId(getFadcInputChId(boardId));
        DayaBay::FadcChannelId outputChannel(fadcChId.outputId());
        
        verbose() << " FeeBoard -> FADCInputChannel -> FADCOutputChannel : " 
                  << boardId.board()
                  << " -> " << fadcChId.connector() 
                  << " -> " << outputChannel.connector() << endreq;
        
        DayaBay::AnalogSignal& fadcVector = output[outputChannel.fullPackedData()];

        // Resize the FADC vectors
        if(fadcVector.size() != boardEnergy.size())
        {
            fadcVector.resize(boardEnergy.size());
        }

        for( unsigned int i = 0; i < boardIdx; i++)
        {
            fadcVector[i] += boardEnergy[i];
        }
            
        verbose() << "FADC vector size " << fadcVector.size() << endreq;
    }
    return StatusCode::SUCCESS;
}

const int ROsFadcReadoutTool::getFadcInputChId(const DayaBay::FeeChannelId& boardId){
    // convert FEE Board Id to Corresponding FADC input ID
    
    // Update to use Cable Mapping Svc...
    // for now create it on the fly.
    int boardNum=boardId.board();
    debug() << "getting input connector id for FEE board" << boardNum << endreq;
    
    DayaBay::FadcChannelId fadcChId(boardNum,boardId.site(),boardId.detectorId());
    return fadcChId.fullPackedData();
}

StatusCode ROsFadcReadoutTool::upsample(const DayaBay::AnalogSignal& signal, 
                                        DayaBay::AnalogSignal& output,
                                        int inputFrequency,
                                        int outputFrequency, int start, int stop)
{
      // Finds the points needed for the interpolation function
      // Make sure the output vector is properly sized to the size of the 
      // the readout window

      verbose() << "Upsampling FADC vector" << endreq;

      double timeStep = 1/((inputFrequency)*(1.e-9));
      double fadcStep = 1/((outputFrequency)*(1.e-9));
      // Factor to get the correct interpolation values
      double muFactor = fadcStep/timeStep;
      int maxIdx = signal.size();
      // Check to make sure the start and stop times are between 0 and the max index
      if (start < 0) start=0;
      if (stop < 0) stop=0;
      if (stop > maxIdx) stop=maxIdx;
      if (start > maxIdx) start=maxIdx;
      if (start > stop) return StatusCode::FAILURE;

      output.resize(stop - start);

      debug() << " start " << start << " stop " << stop << " muFactor " << muFactor << " timeStep " << timeStep <<endreq;
      // Only read out the vector in the readout window (between start and stop)

      for(int i = start; i < stop; i++)
        {  
          // Calculate the x and y values we want to interpolate between
          int ileftPoint = int (i/timeStep);
          int irightPoint = ileftPoint + 1;
          double y1 = signal[ileftPoint];
          double y2 = signal[irightPoint];
          double leftPoint = double( ileftPoint );
          double rightPoint = double( irightPoint );
          double mu = i*muFactor;

          // avoid rounding problems
          if ( abs(mu-leftPoint) < 1.e-6 ) mu = std::max(mu,leftPoint) ;
          if ( abs(mu-rightPoint) < 1.e-6 ) mu = std::min(mu,rightPoint) ; 

          debug() << "Try to interpolate. leftPoint " << leftPoint << " mu " << mu << " rightPoint " << rightPoint 
                  << " mu-leftPoint " << mu-leftPoint << " rightPoint-mu " << rightPoint-mu 
                  << " signal(leftPoint) " << y1 << " signal(rightPoint) " << y2 << endreq;

          // Only interpolate the values of mu that are between the left and right points
          if( leftPoint <= mu && mu <= rightPoint)
            {
              output[i - start] = linearInterpolate(y1, y2, mu, leftPoint, rightPoint);
            }

          else
            {
              fatal() << "Interpolation point out of bounds.leftPoint " << leftPoint << " mu " << mu 
                      << " rightPoint " << rightPoint << " mu-leftPoint " << mu-leftPoint << " rightPoint-mu " << rightPoint-mu << endreq;
              return StatusCode::FAILURE;
            }

          // add different options here if necessary
          debug() << "upsample " << output[i-start] << " <- "
                    << "linear Interp ( " << y1 << "," << y2 << "," << mu
                  << ")" << endreq;
      
        }

      return StatusCode::SUCCESS;
}

double ROsFadcReadoutTool::linearInterpolate(double y1, double y2, double mu, double leftPoint, double rightPoint)
{
  // simple linear interpolation function to upsample the signal to 1 GHz
  verbose() << "Using linear interpolation function" << endreq;
  if(leftPoint <= mu && mu <= rightPoint)
    {
      return( y1 + ((mu - leftPoint)/(rightPoint - leftPoint))*(y2 - y1) );
    }

  warning()<< "This should never happen! Interpolation point out of bounds. Returning value at left point." << endreq;
  return y1;
}

//
StatusCode ROsFadcReadoutTool::digitize(const DayaBay::AnalogSignal& signal,
                              DayaBay::DigitalSignal& output,
                              double range, double bits, double offset) 
{
  // Simple digitization algorithm
  verbose() << "digitizing vectors" << endreq;
  if(output.size() != signal.size()) output.resize(signal.size());
  int maxADC = int( pow(2, bits) );
  int signalN = signal.size();
  const double* inputStart = &signal[0];
  int* outputStart = &output[0];
  for (int i=0; i<signalN; i++) {
    int* curOutput = outputStart + i;
    *(curOutput) = int(((*(inputStart + i))/range)*maxADC + offset);
    //verbose() << "Trying to digitize " << *(inputStart + i) << endreq;
    if (*curOutput>maxADC)
      *curOutput=maxADC;
    if (*curOutput<0)
      *curOutput=0;
    /* Slower, but easier to read version of the same calculation
    outputStart[i] = int((signal[i]/range)*maxADC + offset);
    if (output[i]>maxADC)
      output[i]=maxADC;
    if (output[i]<0)
      output[i]=0;
    */
  }
  return StatusCode::SUCCESS;
}