ROsFeeReadoutTool.cc

Go to the documentation of this file.

#include "ROsFeeReadoutTool.h"
#include "ROsFadcReadoutTool.h"

#include "Conventions/Detectors.h"
#include "Conventions/Electronics.h"
#include "Conventions/Trigger.h"
#include "Conventions/Site.h"

#include "Event/SimReadoutHeader.h"
#include "Event/ReadoutPmtCrate.h"
#include "Event/ReadoutTriggerDataPkg.h"

#include "Event/ElecCrateHeader.h"
#include "Event/ElecFeeCrate.h"
#include "Event/ElecFeeChannel.h"

#include "Event/SimTrigCommand.h"
#include "Event/SimTrigCommandHeader.h"
#include "Event/SimTrigCommandCollection.h"

#include "DataSvc/ICalibDataSvc.h"

#include <set>
#include <vector>

ROsFeeReadoutTool::ROsFeeReadoutTool(const std::string& type,
             const std::string& name, 
             const IInterface* parent)
  : GaudiTool(type,name,parent)
{
  declareInterface< IROsReadoutTool >(this) ;
  // for now add all but RPC's as default.
  m_detectorsToProcess.push_back("DayaBayAD1");
  m_detectorsToProcess.push_back("DayaBayAD2");
  m_detectorsToProcess.push_back("DayaBayIWS");
  m_detectorsToProcess.push_back("DayaBayOWS");
  m_detectorsToProcess.push_back("LingAoAD1");
  m_detectorsToProcess.push_back("LingAoAD2");
  m_detectorsToProcess.push_back("LingAoIWS");
  m_detectorsToProcess.push_back("LingAoOWS"); 
  m_detectorsToProcess.push_back("FarAD1");
  m_detectorsToProcess.push_back("FarAD2");
  m_detectorsToProcess.push_back("FarAD3");
  m_detectorsToProcess.push_back("FarAD4");
  m_detectorsToProcess.push_back("FarIWS");
  m_detectorsToProcess.push_back("FarOWS");
  
  declareProperty("SimDataSvcName",m_simDataSvcName="StaticSimDataSvc",
                  "Name of service to provide FEE channel properties for simulation");
  declareProperty("DetectorsToProcess",m_detectorsToProcess,
                  "List of detectors to process with this tool");
  declareProperty("EnablePeakReadout",m_enablePeakReadout=true,
                  "Switch on/off peak readout mode");
  declareProperty("EnableWaveformReadout",m_enableWaveformReadout=false,
                  "Switch on/off waveform readout mode");
  declareProperty("EnableFadcReadout", m_enableFadcReadout=false,
                  "Switch on/off FADC readout mode");
  declareProperty("WaveformTool", m_waveformToolName="ROsFeeAdcMultiTool", 
                  "Name of waveform readout tool");
  declareProperty("TdcTool", m_tdcToolName="ROsFeeTdcTool", 
                  "Name of tdc readout tool");
  declareProperty("FadcTool", m_fadcToolName="ROsFadcReadoutTool", 
                  "Name of fadc readout tool");      

  //Changed readout Length : (1200/1.5625 = 768)
  declareProperty("ReadoutLength",m_readoutLength=780,
                  "Length of readout window in 640Mhz clock cycles");
  //Changed trigger Latency to 980 (why ?)
  declareProperty("TriggerLatency",m_triggerLatency=966,
                  "Trigger latency in 640Mhz clock cycles");
  //Changed readout offset (200/1.5625 = 128)
  declareProperty("ReadoutOffset",m_readoutOffset=198,
                  "Offset of readout window from trigger tdc in 640Mhz clock cycles");


  declareProperty("PeakFindingLength",m_peakFindingLength=12,
                  "Length of peak finding window in 40Mhz clock cycles");
  declareProperty("PeakFindingOffset",m_peakFindingOffset=0,
                  "Offset of peak finding window from corresponding hit in 40Mhz clock cycles");                
  declareProperty("PeakFindingOverlap",m_peakFindingOverlap=true,
                  "Allow overlapping peakfinding for consecutive hits");
  declareProperty("PreAdcOffset",m_preAdcOffset=4,
                  "Offset of first PreADC cycle from start of ADC peak finding window in 40Mhz clock cycles");
  declareProperty("NhitCycles", m_nhitCycles = DayaBay::NhitCycles, "Nhit cycles");
  declareProperty("EsumCycles", m_eSumCycles = DayaBay::EsumCycles, "Esum cycles");
  declareProperty("FADCOffset", m_fadcOffset = 250, "FADC offset");
  declareProperty("FADCLength", m_fadcLength = 120, "FADC readout window legth");
  declareProperty("MaxFineAdc",m_maxFineAdc=3500,
                  "Maximum fine-range ADC value before using the coarse range");
}

ROsFeeReadoutTool::~ROsFeeReadoutTool(){}

StatusCode ROsFeeReadoutTool::initialize()
{
  std::vector<std::string>::iterator it;
  for(it=m_detectorsToProcess.begin();it!=m_detectorsToProcess.end();++it){
    short int detId = DayaBay::Detector::siteDetPackedFromString(*it);
    DayaBay::Detector det(detId);
    m_detectors.insert(det);
  }
  
  try {
      m_waveformTool = tool<IROsFeeWaveformTool>(m_waveformToolName) ;
    debug() << "Using \"" << m_waveformToolName 
            << "\" for waveform readout " << endreq; 
  }
  catch(const GaudiException& exg) {
      fatal() << "Failed to get Readout Tool: \"" 
              << m_waveformToolName << "\"" << endreq;
      return StatusCode::FAILURE;
  }
  try {
      m_tdcTool = tool<IROsFeeTdcTool>(m_tdcToolName) ;
      debug() << "Using \"" << m_tdcToolName 
            << "\" for TDC readout " << endreq;
  }
  catch(const GaudiException& exg) {
      fatal() << "Failed to get Readout Tool: \"" << m_tdcToolName << "\"" << endreq;
      return StatusCode::FAILURE;
  }
  try {
      m_fadcTool = tool<IROsFadcReadoutTool>(m_fadcToolName) ;
    debug() << "Using \"" << m_fadcToolName 
            << "\" for fadc readout " << endreq; 
  }
  catch(const GaudiException& exg) {
      fatal() << "Failed to get Readout Tool: \"" 
              << m_waveformToolName << "\"" << endreq;
      return StatusCode::FAILURE;
  }
  
  // Get PMT simulation input data service
  m_simDataSvc = svc<ISimDataSvc>(m_simDataSvcName,true);

  if( !m_enablePeakReadout && !m_enableWaveformReadout ) {
      error() << "All readout modes switched off" << endreq;
      return StatusCode::FAILURE;
  }
  if( m_enablePeakReadout )
      info() << "Peak readout mode enabled" << endreq;
  if( m_enableWaveformReadout )
      info() << "Waveform readout mode enabled" << endreq;
      
  return StatusCode::SUCCESS;
}

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

StatusCode ROsFeeReadoutTool::makeReadouts(DayaBay::SimReadoutHeader* roHeader,
                                           std::vector<DayaBay::ReadoutTriggerDataPkg*>& trigDataPkg,
                                           const DayaBay::ElecHeader& elecHeader)
{
  debug() << "running makeReadouts() in FeeReadoutTool" << endreq;
  
  DayaBay::SimReadoutHeader::SimReadoutContainer& outputReadouts 
    = roHeader->readouts();
  
  Context context = roHeader->context();
  //Get Crate Map
  const DayaBay::ElecCrateHeader *ech = elecHeader.crateHeader();
  const DayaBay::ElecCrateHeader::CrateMap crMap = ech->crates();
  DayaBay::ElecCrateHeader::CrateMap::const_iterator crIt;
    
  std::vector<DayaBay::ReadoutTriggerDataPkg*>::iterator pkgIt; 

  //Loop over All trigger data packages
  for(pkgIt = trigDataPkg.begin(); pkgIt != trigDataPkg.end(); ++pkgIt)
  {
    DayaBay::Detector det((*pkgIt)->detector());
    if(m_detectors.find(det) != m_detectors.end()){
      crIt = crMap.find(det);
      const DayaBay::ElecFeeCrate *crate; 
      if(crIt == crMap.end()){
        if((*pkgIt)->frames().size()==0){
          fatal() << "Trying to process a trigger for " << det.detName()
                  << " But Package Has No Frames" << endreq;
          return StatusCode::FAILURE;
        }else{
          fatal() << "Trying to process a trigger for " << det.detName()
                  << " But no crate exists -- skipping" << endreq;
          return StatusCode::FAILURE;
        }
        crate = 0;
      }else{
        crate = dynamic_cast<const DayaBay::ElecFeeCrate*>(crIt->second);
        debug() << "Processing a trigger for " << det.detName()
                  << " found crate at " << crate << endreq;
      }
      if(crate != 0){ // if the crate exists read it out.
        debug()<<"Reading out 1 trigger with: "
               <<((*pkgIt)->frames().size())-1
               <<" masked triggers."<<endreq;
        
        const DayaBay::ReadoutTriggerDataFrame *tdf = (*pkgIt)->frames().at(0);
        
        DayaBay::ReadoutPmtCrate *crateReadout = readoutCrate(tdf,crate,context);
       
        outputReadouts.push_back(new DayaBay::SimReadout(crateReadout,
                                                         roHeader)); 
        
        (*pkgIt)->setReadout(crateReadout);
        crateReadout->setTriggerDataPkg(*pkgIt);
        
        verbose() << "set readout:\n" << (**pkgIt) << endreq;
        
        if( m_enableFadcReadout )
        {
                debug() << " Trying FADC Readout " << endreq;
        
                unsigned int triggerCycle = tdf->cycle();        

                int start = int((triggerCycle*m_eSumCycles)/m_nhitCycles)+m_fadcOffset; 
        
                int stop = start+m_fadcLength;  

                m_fadcTool->readoutFADC(crate,crateReadout,start,stop);
                debug() << "readout map size: " << (crateReadout->fadcReadout()).size() << endreq;
        
        debug() << " Done reading out FADC" << endreq;
        }
      }
    }else{
      verbose() << "Found Trigger For " << det.detName()
                << "but this tool is not supposed "
                << "to process that detector type "
                << "check properties to configure."
                << endreq;
    }
  }
  return StatusCode::SUCCESS;                             
}

DayaBay::ReadoutPmtCrate* ROsFeeReadoutTool::readoutCrate(const DayaBay::ReadoutTriggerDataFrame *tdf, 
                                                          const DayaBay::ElecFeeCrate *cr, Context context)
{ 
  TimeStamp triggerTime(cr->header()->header()->earliest());
  triggerTime.Add( tdf->cycle() / double(DayaBay::NhitFrequencyHz) );
  DayaBay::ReadoutPmtCrate *out_readout 
    = new DayaBay::ReadoutPmtCrate(cr->detector(),0, triggerTime, tdf->triggerType());  
  
  DayaBay::ReadoutPmtCrate::PmtChannelReadouts ro_chMap;
  
  const DayaBay::ElecFeeCrate::ChannelData& chmap = cr->channelData();
  DayaBay::ElecFeeCrate::ChannelData::const_iterator chIt;
  
  for(chIt = chmap.begin(); chIt != chmap.end() ; ++chIt){
    DayaBay::ReadoutPmtChannel *ro_ch = readoutChannel(tdf,&(chIt->second),
                                                       chIt->first,context);
    if (ro_ch == 0) continue;         // test to see if channel has readout
    ro_ch->setChannelId(chIt->first); // set these here since an ElecFeeChannel
    ro_ch->setReadout(out_readout);   //     does not know about them
    ro_chMap[chIt->first] = *ro_ch;   // copy the value by it's address to the map
    ro_ch->setReadout(0);             // reset Readout, to avoid double deletion
    delete ro_ch;                     // delete unnecessary copy
  }
  
  out_readout->setChannelReadout(ro_chMap);
  return out_readout; 
}

DayaBay::ReadoutPmtChannel* ROsFeeReadoutTool::readoutChannel(
                     const DayaBay::ReadoutTriggerDataFrame *tdf, 
                     const DayaBay::ElecFeeChannel *feeChannel,
                     DayaBay::FeeChannelId channelId,
                     Context context)
{  
  // Get the simulation properties for this channel
   int task = 0;
  ServiceMode svcMode(context, task);
  const DayaBay::FeeSimData* feeSimData = 
    m_simDataSvc->feeSimData(channelId, svcMode);
  if(!feeSimData){
    error() << "No Simulation input properties for FEE channel: " 
            << channelId << endreq;
  }

  // convert trigger clock cycle from NhitFrequencyHz to TdcFrequencyHz
  DayaBay::ReadoutPmtChannel *ro_ch = new DayaBay::ReadoutPmtChannel();
  int nHitTrigCycle = tdf->cycle();
  int tdcTrigCyc = (int)(double(nHitTrigCycle) * 
     double(DayaBay::TdcCycles)/double(DayaBay::NhitCycles) + 0.5);
  verbose() << "Conv TDC Clock Cycle: " << tdcTrigCyc << endreq;
  std::vector<int> hit = feeChannel->hit(); 
  int tdcSize = hit.size() *(1+int(DayaBay::TdcCycles/DayaBay::NhitCycles));

  // Find clock cycle of readout window start and end in TDC frequency
  int tdcFirst;    
  if( (int)tdcTrigCyc - m_readoutOffset < 0){
    verbose() << "Invalid Readout TDC Start Cycle " 
        << (int)tdcTrigCyc - m_readoutOffset
              << " Using 0" << endreq;
    tdcFirst = 0;
  } else tdcFirst = tdcTrigCyc - m_readoutOffset;
  verbose() << "TDC First in ns: " << tdcFirst *1.5625 << endreq;
  
  int tdcLast;   
  if( tdcTrigCyc - m_readoutOffset + m_readoutLength > tdcSize){
     verbose() << "Invalid Readout TDC Last Cycle " 
         << tdcTrigCyc - m_readoutOffset + m_readoutLength
         << " Using " << tdcSize << endreq;
     tdcLast = tdcSize;
  } else tdcLast = tdcTrigCyc - m_readoutOffset + m_readoutLength;
  verbose() << "TDC Last in ns: " << tdcLast *1.5625 << endreq;
  verbose() << "reading out tdc window [ " 
         << tdcFirst << "," << tdcLast << " )" << endreq;
  
  // Compute common stop clock cycle for TDC values
  int tdcStopCycle = tdcTrigCyc + m_triggerLatency;
  const DayaBay::DigitalSignal* adc;
  int adcFirst;
  int adcLast;
  int adcPedFirst;
      
  // Compute variables for peak readout mode
  if ( m_enablePeakReadout ){                                            
    // Read out clock numbers which represent the time when signal crosses threshold 
    const std::vector<int>& tdcOrg = feeChannel->tdc();
    std::vector<int> tdc_vec;
    StatusCode sc = m_tdcTool->readoutTdc(tdcOrg,tdcFirst,tdcLast,tdc_vec);
    if( !sc.isSuccess() || tdc_vec.empty() ) {
        delete ro_ch;
        return 0;
    }
    
    // Loop over TDC values to compute corresponding peak ADC values and peaking clock cycles
    int nTdc = tdc_vec.size();
    std::vector<int> tdc_out(nTdc); 
    std::vector<int> tdcHitCount_out(nTdc);
    std::vector<int> adc_out(nTdc); 
    std::vector<int> adcCycle_out(nTdc);
    std::vector<int> adcRange_out(nTdc);
    std::vector<int> preAdc_out(nTdc);

    DayaBay::FeeGain::FeeGain_t gain; 
    int currentAdc;
    
    int currentTdc = -1000;  
    
    for ( int i = 0; i < nTdc; i++ ){
      // Check if preceeding peak finding is finished
      if ( tdc_vec[i] - currentTdc  > m_peakFindingLength * 16 || m_peakFindingOverlap == true){
        
        // Assume the high gain ADC will be read out 
        adc = &(feeChannel->adcHigh());         
        gain = DayaBay::FeeGain::kHigh;
        verbose() << "Maximum available ADC window [ " << 0 << "," << adc->size() 
                 << " )" << endreq;
        
        // Search for peak within 300 ns window from corresponding TDC clock cycle
        currentTdc = tdc_vec[i];
        currentAdc = (int) (double(currentTdc) * double(DayaBay::AdcCycles)/double(DayaBay::TdcCycles) );
        
        if( currentAdc - m_peakFindingOffset < 0){
          verbose() << "Invalid ADC Peak Finding Start Cycle " 
                    <<  currentAdc - m_peakFindingOffset
                    << " Using 0" << endreq;
          adcFirst = 0;
        } 
        else adcFirst = currentAdc - m_peakFindingOffset;
        verbose() << "ADC First in ns: " << adcFirst * 25. << endreq;
  
        if( currentAdc - m_peakFindingOffset + m_peakFindingLength > adc->size() ){
          verbose() << "Invalid ATDC Peak Finding Last Cycle " 
                    << currentAdc - m_peakFindingOffset + m_peakFindingLength
                    << " Using " << adc->size() -1 << endreq;
          adcLast = adc->size() -1;
        } 
        else adcLast = currentAdc - m_peakFindingOffset + m_peakFindingLength;
        verbose() << "ADC Last in ns: " << adcLast * 25. << endreq;
       
        verbose() << "Reading out peak finding window [ " 
                  << adcFirst << "," << adcLast +1 << ")" << endreq;
                 
        // Check if high gain assumption was correct
        for( int j = adcFirst; j < adcLast; ++j ){
          if( adc->at(j) >= m_maxFineAdc ){
            gain = DayaBay::FeeGain::kLow;
            adc = &(feeChannel->adcLow());
            debug() << "Switching to Low gain ADC since High was saturated" 
                   << endreq;
            break;
          }
        } // end check high gain adc
       
        verbose() << "Ready to readout ADC's" << endreq;
       
        // Find Peak ADC value
        int adcPeak = -1;
        int adcPeakCycle = 0; 
        for(int adcCycle = (int)adcFirst; adcCycle < (int)adcLast; ++adcCycle){
          if( adc->at(adcCycle) > adcPeak ){
            // Peak clock cycle will be read out wrt start of the peak finding window
            adcPeak = adc->at(adcCycle);
            adcPeakCycle = adcCycle-adcFirst;
            verbose() << "Found New Max ADC at " << adcPeakCycle << endreq;
          }
        }

        // Compute PreADC value    
        int preAdc = 0;
        adcPedFirst = adcFirst - m_preAdcOffset;
        for(int adcCycle = adcPedFirst; adcCycle<adcPedFirst + 4; adcCycle++){
          if(adcCycle < 0) {
            verbose() << "Invalid PreADC Cycle " 
                      <<  adcCycle
                      << " Using standard baseline value" << endreq;
            if(gain == DayaBay::FeeGain::kHigh) preAdc += int(feeSimData->m_adcBaselineHigh);
            if(gain == DayaBay::FeeGain::kLow) preAdc += int(feeSimData->m_adcBaselineLow);
          }
          else preAdc += adc->at(adcCycle);
        }
        verbose() << "First PreADC cycle in ns: " << adcPedFirst * 25. << endreq;
        preAdc = preAdc & 0x0fff;
        verbose() << "Reading out adc " <<  adcPeak
                  << " at cycle " << adcPeakCycle << endreq;
        verbose() << "Reading out preADC " << preAdc << endreq;
        
        adc_out[i] = adcPeak; 
        adcCycle_out[i] = adcPeakCycle; 
        adcRange_out[i] = gain;
        preAdc_out[i] = preAdc;

      } //end if -- check if preceding peak find alg is finished
      else {
       verbose() << "No peak finding (previous not finished)" << endreq;
       // Fill adc vectors with zeros
       adc_out[i] = 0; 
       adcCycle_out[i] = 0;
       adcRange_out[i] = 0;
       preAdc_out[i] = 0;
      } // end else
      
      // TDC will be read out wrt common stop time
      tdc_out[i] = tdcStopCycle - tdc_vec[i];  
      tdcHitCount_out[i] = nTdc - i;
    }// end tdc forloop
    
    // finalize peaking mode variables for current channel
    ro_ch->setTdc(tdc_out);
    ro_ch->setTdcHitCount(tdcHitCount_out);
    ro_ch->setAdc(adc_out);
    ro_ch->setAdcCycle(adcCycle_out);
    ro_ch->setAdcRange(adcRange_out);
    ro_ch->setPedestal(preAdc_out);
  }
  
  // Compute variables for waveform readout mode
  if ( m_enableWaveformReadout ){
    adcFirst = (int) (double(tdcFirst) * double(DayaBay::AdcCycles)/double(DayaBay::TdcCycles) + 0.5);
    adcLast = (int) (double(tdcLast) * double(DayaBay::AdcCycles)/double(DayaBay::TdcCycles) + 0.5);
    int adcStopCycle = (int) (double(tdcStopCycle) * double(DayaBay::AdcCycles)/double(DayaBay::TdcCycles) + 0.5);

    std::vector<int> waveAdcLow_out; 
    std::vector<int> waveAdcHigh_out;
    std::vector<int> waveAdcCycle_out;
    
    // readout low gain waveform
    adc = &(feeChannel->adcLow());
    StatusCode sc = m_waveformTool->readoutWaveform(*adc, adcFirst, adcLast, adcStopCycle, waveAdcLow_out, waveAdcCycle_out);
    if( !sc.isSuccess() || waveAdcLow_out.empty() ) {
        delete ro_ch;
        return 0;
    }
    waveAdcCycle_out.clear();
    
    // readout high gain waveform
    adc = &(feeChannel->adcHigh());
    sc = m_waveformTool->readoutWaveform(*adc, adcFirst, adcLast, adcStopCycle, waveAdcHigh_out, waveAdcCycle_out);
    if( !sc.isSuccess() || waveAdcHigh_out.empty() ){
        delete ro_ch;
        return 0;
    }
    
    // finalize waveform mode variables for current channel
    ro_ch->setWaveAdcLow(waveAdcLow_out);
    ro_ch->setWaveAdcHigh(waveAdcHigh_out);
    ro_ch->setWaveAdcCycle(waveAdcCycle_out);
  }        
  return ro_ch;
}

// provide mutate to satisfy interface IROsReadoutTool.
StatusCode ROsFeeReadoutTool::mutate(DayaBay::SimReadoutHeader* roHeader,
                                     std::vector<DayaBay::ReadoutTriggerDataPkg*>& trigDataPkg,
                                     const DayaBay::ElecHeader& elecHeader)
{     
    verbose() << "calling ROsFeeReadoutTool.mutate()" << endreq;
    return makeReadouts(roHeader,trigDataPkg,elecHeader);
}