In This Package:

EsIdealFeeTool.cc

Go to the documentation of this file.

#include "EsIdealFeeTool.h"
#include "Event/ElecFeeCrate.h"
#include "Event/ElecPulseCollection.h"
#include "Event/ElecPmtPulse.h"
#include "GaudiKernel/SystemOfUnits.h"
#include "Conventions/Trigger.h"

#include <sstream>
#include <string.h>

EsIdealFeeTool::EsIdealFeeTool(const std::string& type,
                   const std::string& name, 
                   const IInterface* parent)
  : GaudiTool(type,name,parent),
    m_cableSvc(0)
{
    declareInterface< IEsFrontEndTool >(this);
    
    declareProperty("CableSvcName",m_cableSvcName="StaticCableSvc",
                    "Name of service to map between detector, hardware, and electronic IDs");
    declareProperty("SimDataSvcName",m_simDataSvcName="StaticSimDataSvc",
                    "Name of service to provide FEE channel properties for simulation");
    declareProperty("SimFrequency", m_simFrequency = DayaBay::TdcFrequencyHz,
                    "Simulation Frequency");
    declareProperty("ESumFrequency", m_eSumFrequency = DayaBay::EsumFrequencyHz,
                    "Esum Frequency");
    declareProperty("TriggerWindowCycles", 
                    m_triggerWindowCycles=DayaBay::TriggerWindowCycles,
                    "Number of trigger window ccyles");
    declareProperty("GainFactor", 
                    m_gainFactor=2.,
                    "Mean PMT gain (in units of 10^7)");
    declareProperty("NoiseAmp",m_noiseAmp=0.5e-3*Gaudi::Units::volt,
                    "Voltage amplitude of noise");
    declareProperty("SpeAmp",m_speAmp=0.0035,
                    "Spe pulse height in V at 1e7 gain");
    declareProperty("DiscThreshScale", m_discThreshScale=1,
                    "Factor to scale all channels' discriminator threshold");
    declareProperty("EnableNonlinearity", m_enableNonlinearity=true, 
                    "Turns on/off all nonlinear effects");
    declareProperty("EnableNoise",m_enableNoise=false,"Turn on/off noise individually");
    declareProperty("EnableDynamicWaveform", m_enableDynamicWaveform=true, 
                    "Turns on/off dynamic waveform individually");
    declareProperty("EnableRinging", m_enableRinging=true, 
                    "Turns on/off ringing individually");
    declareProperty("EnableOvershoot", m_enableOvershoot=true, 
                    "Turns on/off overshoot individually");
    declareProperty("EnableSaturation", m_enableSaturation=true, 
                    "Turns on/off saturation individually");
    declareProperty("EnableAcCoupling", m_enableAcCoupling=true, 
                    "Turns on/off high pass filter before discriminator");
    declareProperty("EnableESumTotal", m_enableESumTotal=true, 
                    "Turns on/off Total ESum simulation");
    declareProperty("EnableESumH", m_enableESumH=true, "Turns on/off upper ESum simulation");
    declareProperty("EnableESumL", m_enableESumL=true, "Turns on/off lower ESum simulation");
    declareProperty("EnableFastSimMode", m_enableFastSimMode=false, 
                    "Turns on/off fast simulation mode to save CPU time");
    declareProperty("PulseCountSlotWidth", m_pulseCountSlotWidth=5e-9*Gaudi::Units::second,
                    "pulse count time slot size");
    declareProperty("PulseCountWindow", m_pulseCountWindow=2, 
                    "Number of counted slots before and after a pulse");
    declareProperty("LinearityThreshold", m_linearityThreshold=20, 
                    "Threshold in P.E. for activating nonlinear pulse shape");
    declareProperty("EnablePretrigger", m_enablePretrigger=true, 
                    "Turns on/off pretrigger");
    declareProperty("TotalESumTriggerThreshold", 
                    m_eSumTriggerThreshold=DayaBay::Trigger::kADESumThreshold, 
                    "Total ESum trigger threshold as set in TrigSim");
    declareProperty("NHitTriggerThreshold", 
                    m_nHitTriggerThreshold=DayaBay::Trigger::kADthreshold, 
                    "nHit trigger threshold as set in TrigSim");
    declareProperty("ADCOffset", m_adcOffset = 500, "Esum Frequency");
    declareProperty("ADCRange", m_adcRange = 1, "Esum Range");
    declareProperty("ADCBits", m_adcBits = 10, "Esum Bits");
}

EsIdealFeeTool::~EsIdealFeeTool(){}

StatusCode EsIdealFeeTool::initialize()
{
  // Get Cable Service
  m_cableSvc = svc<ICableSvc>(m_cableSvcName,true);

  // Get PMT simulation input data service
  m_simDataSvc = svc<ISimDataSvc>(m_simDataSvcName,true);

  StatusCode sc = loadResponse();
  if(sc.isFailure()) return sc;

  // Random number service
  IRndmGenSvc *rgs = 0;
  if (service("RndmGenSvc",rgs,true).isFailure()) {
    fatal() << "Failed to get random service" << endreq;
    return StatusCode::FAILURE;        
  }
  
  sc = m_gauss.initialize(rgs, Rndm::Gauss(0,1));
  if (sc.isFailure()) {
    fatal() << "Failed to initialize gaussian random numbers" << endreq;
    return StatusCode::FAILURE;
  }
  // Set up pretrigger threshold
  int preThr = int(0.8 * m_nHitTriggerThreshold + 0.5);
  int esumPreThr = int (0.8 * 1500 * m_eSumTriggerThreshold/Gaudi::Units::volt + 0.5);
  if(preThr > esumPreThr) preThr = esumPreThr;
  m_pretriggerThreshold = preThr;
  if( m_enablePretrigger ){
    info() << "Pretrigger threshold is " << m_pretriggerThreshold << " pe " << endreq;
  }
  else {
    info() << "Pretrigger disabled" << endreq;
    info() << "Convolution threshold is " << m_pretriggerThreshold << " pe " << endreq;
  }
  
  if( m_enableFastSimMode){
    info() << "Fast electronics simulation mode enabled: Note that Upper and Lower ESum trigger do not yield meaningful results" << endreq;
    m_linearityThreshold = 100000;
    m_enableESumH = false;
    m_enableESumL = false;
  }
  if( !m_enableNonlinearity  ) {
    info() << "All nonlinear effects disabled." << endreq; 
    m_enableOvershoot = false;   
    m_enableSaturation = false;
    m_enableRinging = false;
    m_enableDynamicWaveform = false;
    m_enableAcCoupling = false;
  }
  else {
    info() << "Nonlinear model enabled for channels with more than " 
       << m_linearityThreshold << " hits " << endreq;
        
    if( !m_enableAcCoupling ) 
      info() << "AC coupling disabled" << endreq; 
    else 
      debug() << "AC coupling enabled" << endreq;
  
    if( !m_enableOvershoot ) 
      info() << "Overshoot disabled" << endreq; 
    else 
      debug() << "Overshoot enabled" << endreq;
      
    if( !m_enableRinging ) 
      info() << "Ringing disabled" << endreq; 
    else 
      debug() << "Ringing enabled" << endreq;
      
    if( !m_enableSaturation ) 
      info() << "Saturation disabled" << endreq; 
    else 
      debug() << "Saturation enabled" << endreq;
      
    if( !m_enableDynamicWaveform ) 
      info() << "Dynamic waveform disabled" << endreq; 
    else 
      debug() << "Dynamic waveform enabled" << endreq;
  }
  return StatusCode::SUCCESS;
}

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

StatusCode EsIdealFeeTool::generateSignals(DayaBay::ElecPulseCollection* pulses,
                       DayaBay::ElecCrate* elecCrate)
{    
  // Ensure that crate is a FeeCrate
  DayaBay::ElecFeeCrate* crate = 0;
  crate = dynamic_cast<DayaBay::ElecFeeCrate*>(elecCrate);
  if(crate == NULL){
    error() << "generateSignals(): Crate is not a FEE Crate." << endreq;
    return StatusCode::FAILURE;
  }
  
  // Context, ServiceMode for this data
  Context context(pulses->detector().site(), SimFlag::kMC, 
          pulses->header()->header()->timeStamp(),
          pulses->detector().detectorId());
  int task = 0;
  ServiceMode svcMode(context, task);

  // Initialize crate if necessary
  if(crate->channelData().size() == 0){
    // First, get list of all connected FEE channels
    const std::vector<DayaBay::DetectorSensor>& pmtList 
      = m_cableSvc->sensors( svcMode );
    int nPmts = pmtList.size();
    for(int pmtIdx = 0; pmtIdx < nPmts; pmtIdx++){
      // Add channel to crate
      DayaBay::FeeChannelId channelId(m_cableSvc->elecChannelId(pmtList[pmtIdx], svcMode).fullPackedData());
      crate->addChannel(channelId);
      verbose() << "Adding channel to crate: " << channelId << endreq;
    }
  }

  // Prepare time range of the simulation
  TimeStamp earliest = pulses->header()->header()->earliest();
  TimeStamp latest = pulses->header()->header()->latest();
  // The time window for the crate signals
  double simTime = latest - earliest;
  verbose() << "Simulation time is " << simTime << endreq;
  // The number of samples in time given the simulation frequency
  int simSamples = int(simTime * m_simFrequency);
  verbose() << "Simulation has " << simSamples << " samples at frequency "
     << m_simFrequency << " Hz" << endreq;
  
  // Compute time windows that will possibly be readout
  int nTimeSlices = int(simTime / 150e-9 + 0.5);
  //std::vector<int> hitCountSlices(nTimeSlices,0);
  //int* hitCountSlices = new int[nTimeSlices];
  int hitCountSlices[nTimeSlices];
  for ( int i = 0; i < nTimeSlices; ++i ) hitCountSlices[i] = 0;

  // Organize Pulses by Channel
  std::map<DayaBay::ElecChannelId, std::vector<DayaBay::ElecPmtPulse*> > pulseMap;
  StatusCode sc;
  sc = mapPulsesByChannel( pulses->pulses(), pulseMap, hitCountSlices );      
  if(sc.isFailure()) return sc;
  debug() << "Mapped Pulses" << endreq;

  // Pretrigger:
  // Compute map of time windows that will possibly read out  
  m_adcConvolutionWindows.clear();
  m_esumConvolutionWindows.clear();
  int adcConvStart, adcConvEnd;
  int esumConvStart, esumConvEnd;
  std::map<int,int>::reverse_iterator sliceIt;
  for(int i = 1; i < nTimeSlices; i++){
    if( hitCountSlices[i] + hitCountSlices[i-1] < m_pretriggerThreshold ) continue;
    if (i - 5 < 0) adcConvStart = 0;
    else adcConvStart = (i - 5);
    if (i - 1 < 0) esumConvStart = 0;
    else esumConvStart = (i - 1);
    
    if (i + 5 > nTimeSlices) adcConvEnd = nTimeSlices;
    else adcConvEnd = (i + 5);
    if (i + 1 > nTimeSlices) esumConvEnd = nTimeSlices;
    else esumConvEnd = (i + 1);
    
    if( m_adcConvolutionWindows.empty() ) {
      m_adcConvolutionWindows[adcConvStart] = adcConvEnd;
      m_esumConvolutionWindows[esumConvStart] = esumConvEnd;
      continue;
    }
    sliceIt = m_adcConvolutionWindows.rbegin();
    if( (adcConvStart - sliceIt->second > 2) ) 
      m_adcConvolutionWindows[adcConvStart ] = adcConvEnd;
    else
      m_adcConvolutionWindows[sliceIt->first] = adcConvEnd;
      
    sliceIt = m_esumConvolutionWindows.rbegin();
    if( (esumConvStart - sliceIt->second > 2) ) 
      m_esumConvolutionWindows[esumConvStart ] = esumConvEnd;
    else                                         
      m_esumConvolutionWindows[sliceIt->first] = esumConvEnd;
  }
  debug() << "Number of time windows to be possibly read out: " << m_adcConvolutionWindows.size() << endreq;

  // If no time windows is considered of possible interest, skip the entire event
  if(m_adcConvolutionWindows.size()==0 && m_enablePretrigger) {
    debug() << "Event below pretrigger threshold, no waveform simulation" << endreq;
    return StatusCode::SUCCESS;
  }
  
  // Loop over channels, and fill with the simulated signals
  const DayaBay::ElecFeeCrate::ChannelData& channelData = crate->channelData();
  DayaBay::ElecFeeCrate::ChannelData::const_iterator channelIter, 
    done = channelData.end();  
  for(channelIter = channelData.begin(); channelIter != done; ++channelIter){
    DayaBay::FeeChannelId channelId = channelIter->first;
    DayaBay::ElecFeeChannel& channel = crate->channel(channelId);

    //Fill each channel with signals
    sc = generateOneChannel(channelId, pulseMap[channelId], channel, 
                simTime, svcMode);
    if( !sc.isSuccess() ) return sc;

  } // End loop over channels

  StatusCode status = updateBoardSignals(crate);
  if (status.isFailure()) return status;
  status = updateEsumSignals(crate);
  if (status.isFailure()) return status;  

  verbose() << "Crate\n" << *crate << endreq;

  return StatusCode::SUCCESS;
}

StatusCode EsIdealFeeTool::generateOneChannel(
                 const DayaBay::FeeChannelId& channelId,
                 std::vector<DayaBay::ElecPmtPulse*>& channelPulses,
                 DayaBay::ElecFeeChannel& channel,
                 double simTime,
                 const ServiceMode& svcMode)
{
  // Generate the FEE signals caused by the PMT pulses for one channel

  // Get the simulation properties for this channel
  const DayaBay::FeeSimData* feeSimData = 
    m_simDataSvc->feeSimData(channelId, svcMode);
  if(!feeSimData){
    error() << "No Simulation input properties for FEE channel: " 
        << channelId << endreq;
    return StatusCode::FAILURE;
  }

  // The number of samples in time given the simulation frequency
  unsigned int simSamples = int(simTime * m_simFrequency);

  // Prepare Raw Signal
  if(m_rawSignal.size() != simSamples) m_rawSignal.resize( simSamples );
  if(m_discriminatorSignal.size() != simSamples) m_discriminatorSignal.resize( simSamples );
  if(m_shapedSignal.size() != simSamples) m_shapedSignal.resize( simSamples );
  double* rawStart = &m_rawSignal[0];
  double* discriminatorStart = &m_discriminatorSignal[0];
  double* shapedStart = &m_shapedSignal[0];
  for( unsigned int sigIdx = 0; sigIdx!=simSamples; sigIdx++){
    *(rawStart + sigIdx) = 0;
    *(discriminatorStart + sigIdx) = 0;
    *(shapedStart + sigIdx) = 0;
  } 
  
  // Add noise to raw signal if requested
  if( m_enableNoise ){
    for(unsigned int index=0; index < simSamples; index++){
      *(rawStart + index) += noise();
      *(discriminatorStart + index) += noise();
    }
  }

  // Pulses on this channel?
  if( channelPulses.size() > 0 ){
    verbose() << "Channel " << channelId << " has " << channelPulses.size() 
          << " pulses." << endreq;
    // Use linear or nonlinear pulse shapes?
    if( m_enableNonlinearity && (channelPulses.size() > m_linearityThreshold) ){
      // Use nonlinear pulse shape

      verbose() << "Applying nonlinearity:" << channelPulses.size() 
        << " pulses (threshold=" << m_linearityThreshold << ")"  
        << endreq;

      // SJ: Prepare time slots for pulse counting
      int numPulseTimeSlots = int(simTime*Gaudi::Units::second
                  /m_pulseCountSlotWidth) + 1;
      verbose() << "Number of time slots for afterpulse counting = " 
        << numPulseTimeSlots << std::endl; 
      std::vector<int> pulseTimeSlots(numPulseTimeSlots);

      // SJ: Fill pulse count time slots
      // Count the number of pulses for each time slot to model nonlinearity
      for (int i = 0; i < numPulseTimeSlots; i++ ) pulseTimeSlots[i] = 0;
      std::vector<DayaBay::ElecPmtPulse*>::const_iterator pulseIter, pulseDone 
    = channelPulses.end();
      for(pulseIter=channelPulses.begin(); pulseIter != pulseDone; ++pulseIter){
        DayaBay::ElecPmtPulse* pulse = *pulseIter;
        int timeSlot = int(pulse->time() * Gaudi::Units::nanosecond 
                   / m_pulseCountSlotWidth);
        int pulseNo = int(pulse->amplitude()+0.5);
        pulseTimeSlots[timeSlot]+=pulseNo;
        verbose() << "Added " << pulseNo << " pulses in time slot " << timeSlot << endreq;
      }

      // Loop over pulses, add each to raw signal
      for(pulseIter=channelPulses.begin(); pulseIter != pulseDone; ++pulseIter){
        DayaBay::ElecPmtPulse* pulse = *pulseIter;
        int tOffset = int(pulse->time()*1e-9 * m_simFrequency);
        float amplitude = pulse->amplitude();
    
        // SJ: Fill pulse time slots
        // Count the total number of pulses within a nearby time window
        // This number of pulses is used to determine the nonlinearity
        int nPulse=0;
        int timeSlot = int(pulse->time() * Gaudi::Units::nanosecond 
                   / m_pulseCountSlotWidth);
        for (int i = timeSlot - m_pulseCountWindow; 
             i <= timeSlot + m_pulseCountWindow; i++){
          if(i>=0 && i<numPulseTimeSlots) nPulse += pulseTimeSlots[i];
        }
        verbose() << "Number of pulses in time window: " << nPulse << endreq;   
    
        // Now add pulses to the raw signal
        int nPulseSamples = m_pmtPulse.size();
        for(int index = 0; index < nPulseSamples; index++){
          unsigned int sampleIdx = tOffset + index;
          double idxTime = index * (1. / m_simFrequency);
          if(sampleIdx>0 && sampleIdx<simSamples){
            // SJ: Add pulses to raw signal
            if( m_enableDynamicWaveform ){
              m_rawSignal[sampleIdx] += - amplitude * pmtPulse( idxTime, nPulse, m_enableOvershoot );
              m_discriminatorSignal[sampleIdx] += - amplitude * pmtPulse( idxTime, nPulse, false );
            }else{
              if(m_enableOvershoot) m_rawSignal[sampleIdx] += - amplitude * m_pmtPulse[index];
              else m_rawSignal[sampleIdx] += - amplitude * m_pmtPulse_noOvershoot[index];
              m_discriminatorSignal[sampleIdx] += - amplitude * m_pmtPulse_noOvershoot[index];
            }
          }
        }
      } // End loop over pulses

      // Add saturation and ringing to the raw signal
      if( (channelPulses.size() > 0) 
      && (m_enableSaturation || m_enableRinging) ){
        // Determine the raw signal maximum for use in ringing
        // At the same time, add saturation if needed 
        //double rawMax;
        double rawMax;
        int maxSample;
        rawMax = 0;
        for (unsigned int index=0; index<simSamples; index++) {
          // SJ: Simulate saturation of the raw signal
          if( m_rawSignal[index] > rawMax ) {
            rawMax = m_rawSignal[index];
            maxSample = index;
          }
          if(m_enableSaturation){
            m_rawSignal[index] *= satFactor(m_rawSignal[index]);
            m_discriminatorSignal[index] *= satFactor(m_discriminatorSignal[index]);
          }
        }
        
        // Add ringing if requested
        if( m_enableRinging && rawMax != 0 ){
          for (unsigned int index=maxSample; index<simSamples; index++) {
            double idxTime =  index * (1. / m_simFrequency);
            m_rawSignal[index] += - ringing ( idxTime, rawMax );
          }
        }
      } // End of saturation and ringing
      
      // Generate non-linear shaped signal
      shape(m_rawSignal,m_shapedSignal,m_shapingResponse,m_adcConvolutionWindows);
      m_shapedSignal.resize( simSamples );
      // SJ: Scale the resulting convolution to preserve the pulse area
      for(unsigned int index=0; index < simSamples; index++) {
        m_shapedSignal[index] /= m_shapingSum;
      }
      // End non-linear pulse addition

    }else{
      // Use linear pulse shape
      
      double* pmtPulseStart;
      if(m_enableOvershoot)  pmtPulseStart = &m_pmtPulse[0];
      else pmtPulseStart = &m_pmtPulse_noOvershoot[0];
      double* pmtDiscriminatorPulseStart = &m_pmtPulse_noOvershoot[0];
      double* shapedPmtPulseStart = &m_shapedPmtPulse[0];
      double* ringingStart = &m_ringing[0];
      double* shapedRingingStart = &m_shapedRinging[0];

      // Loop over pulses
      std::vector<DayaBay::ElecPmtPulse*>::const_iterator pulseIter, pulseDone 
    = channelPulses.end();
      for(pulseIter=channelPulses.begin(); pulseIter != pulseDone; ++pulseIter){
        DayaBay::ElecPmtPulse* pulse = *pulseIter;
        int tOffset = int(pulse->time()*1e-9 * m_simFrequency);
        float amplitude = pulse->amplitude();
    
        // Now add pulses to the raw and shaped signals
        int nPulseSamples = m_pmtPulse.size();
        for(int index = 0; index < nPulseSamples; index++){
          unsigned int sampleIdx = tOffset + index;
          
          if(sampleIdx>0 && sampleIdx<simSamples){
            *(rawStart + sampleIdx) -= amplitude * (*(pmtPulseStart + index));
            *(discriminatorStart + sampleIdx) -= amplitude * (*(pmtDiscriminatorPulseStart + index));
            *(shapedStart + sampleIdx) -= amplitude * (*(shapedPmtPulseStart+index));
          }
        }
      } // End loop over pulses

      // Add ringing to the raw signal
      if( (channelPulses.size() > 0)  && (m_enableSaturation || m_enableRinging) ){
        // Determine the raw signal maximum for use in ringing
        // At the same time, add saturation if needed 
        double rawMax;
        int maxSample;
        rawMax = 0;
        for (unsigned int index=0; index<simSamples; index++) {
          double* curValue = rawStart + index;
          if( *curValue > rawMax ) {
            rawMax = *curValue;
            maxSample = index;
          }
          if(m_enableFastSimMode && m_enableSaturation){
            *(rawStart + index) *= satFactor(*(rawStart + index));
            *(rawStart + index) *= satFactor(-10 * (*(rawStart + index)));
          }
        }
        // Add ringing and/or saturation if requested
        if( m_enableRinging && rawMax != 0 ){
          unsigned int ringEnd = m_ringing.size();
          double ringScaling = -rawMax / m_pulseMax;
          if( ringEnd + maxSample > simSamples ) 
            ringEnd = simSamples - maxSample;
          for (unsigned int ringIdx=0; ringIdx<ringEnd; ringIdx++) {
            int index = maxSample + ringIdx;
            if(m_enableFastSimMode && m_enableSaturation){
              *(discriminatorStart + index) *= satFactor(*(discriminatorStart + index));
              *(shapedStart + index) *= pow(satFactor(15.* (fabs(*(shapedStart + index))) ),0.75);
            }
            *(rawStart + index) += ringScaling * (*(ringingStart+ringIdx));
            *(shapedStart + index) += ringScaling * (*(shapedRingingStart+ringIdx));
          }
        }
      } // End saturation and ringing
      
    } // End addition of linear PMT pulse

  } // End addition of pulses to raw signal

  // Convert raw signals to TDC, hit, ADC, Esum values
  // TDC values
  if (m_enableAcCoupling) sample( m_discriminatorSignal, m_tdcSignal, 
                  m_simFrequency, DayaBay::TdcFrequencyHz );
  else sample( m_rawSignal, m_tdcSignal,
           m_simFrequency, DayaBay::TdcFrequencyHz );
  std::vector<int> tdcValues;
  discriminate( m_tdcSignal, feeSimData->m_channelThreshold, tdcValues, 
        DayaBay::Threshold::kAbove);
  channel.setTdc(tdcValues);
  verbose() << "Channel " << channelId << ": " << tdcValues.size() 
        << " TDC values." << endreq;
  
  // Assemble Hit signal (for trigger) based on TDC clock values
  // The hit signal is a binary 0 or 1 vs. time at the Nhit frequency.
  unsigned int hitSamples = convertClock( m_tdcSignal.size(),
                      DayaBay::TdcFrequencyHz,
                      DayaBay::NhitFrequencyHz ) + 1;
  // Resize and clear the buffer
  if( m_hitSignal.size() != hitSamples) m_hitSignal.resize( hitSamples );
  int* hitStart = &m_hitSignal[0];
  for( unsigned int hitIdx = 0; hitIdx != hitSamples; hitIdx++ ){
    *(hitStart+hitIdx) = 0;
  }
  if( tdcValues.size() > 0 ){
    std::vector<int>::const_iterator tdcIter, tdcDone = tdcValues.end();
    for(tdcIter = tdcValues.begin(); tdcIter != tdcDone; tdcIter++){
      int tdcClock = *tdcIter;
      int hitClock = convertClock( tdcClock, DayaBay::TdcFrequencyHz, 
                   DayaBay::NhitFrequencyHz );
      m_hitSignal[hitClock] = 1;
    }
  }
  channel.setHit(m_hitSignal);
  verbose() << "Hit signal size: " << m_hitSignal.size() << endreq;
    
  // ADC Values: Digitized signals vs. time for both the high and low gains
  // Reduce shaped pulse to ADC frequency
  sample( m_shapedSignal, m_adcAnalogSignal,
      m_simFrequency, DayaBay::AdcFrequencyHz );
  // 12-bit ADC after shaped signals. One for adcHigh and one for adcLow.
  int adcBits = 12;
  digitize( m_adcAnalogSignal, m_adcHighGain,
        feeSimData->m_adcRangeHigh, adcBits, 
        feeSimData->m_adcBaselineHigh );
  digitize( m_adcAnalogSignal, m_adcLowGain,
        feeSimData->m_adcRangeLow, adcBits, 
        feeSimData->m_adcBaselineLow );
  channel.setAdcHigh( m_adcHighGain );
  channel.setAdcLow( m_adcLowGain );
  verbose() << "ADC signal size: " << m_adcHighGain.size() << endreq;
  
  // Channel energy: contribution to the Esum signal 
  //sample( m_rawSignal, channel.energy(),
      //       m_simFrequency, m_eSumFrequency );
  //verbose() << "Energy signal size: " << m_energySignal.size() << endreq;
  channel.setEnergy(m_rawSignal);
  return StatusCode::SUCCESS;
}  // End generate One channel

// Private functions

StatusCode EsIdealFeeTool::loadResponse(){ 
  // SJ: Load the CR-(RC)^4 response data
  double dT_seconds = (1. / m_simFrequency);
  double dT_units = dT_seconds*CLHEP::second;
  int nShapingSamples = int(DayaBay::preTimeTolerance/dT_units);
  int nPulseSamples = 3*nShapingSamples;
  int nRingingSamples = int(DayaBay::postTimeTolerance/dT_units);
  m_shapingResponse.resize(nShapingSamples);
  m_esumResponse.resize(nPulseSamples);
  m_pmtPulse.resize(nPulseSamples);
  m_pmtPulse_noOvershoot.resize(nPulseSamples);
  m_ringing.resize(nRingingSamples);
  m_shapingSum = 0;
  m_esumShapingSum = 0;
  m_pulseMax = 0;

  for (int i=0; i<nShapingSamples; i++) {
    m_shapingResponse[i] = pulseShaping(i*dT_seconds);
    m_shapingSum += m_shapingResponse[i];
  }
  
  for (int i=0; i<nPulseSamples; i++){
    m_esumResponse[i] = esumShaping(i*dT_units/CLHEP::nanosecond);
    m_esumShapingSum += m_esumResponse[i];
    }

  for (int i=0; i<nPulseSamples; i++) {
    m_pmtPulse[i] = pmtPulse(i*dT_seconds,1,true);
    m_pmtPulse_noOvershoot[i] = pmtPulse(i*dT_seconds,1,false);
    // Keep track of 
    if( m_pmtPulse[i] < m_pulseMax ) m_pulseMax = m_pmtPulse[i];

  }
  // Invert to be a positive-voltage peak
  m_pulseMax *= -1;
 
  for (int i=0; i<nRingingSamples; i++) {
    m_ringing[i] = ringing(i*dT_seconds, m_pulseMax);
  }

  // Convolve the pmt pulse with the CR-(RC)^4 shaping response
  convolve( m_pmtPulse, m_shapingResponse, m_shapedPmtPulse );
  m_shapedPmtPulse.resize( m_pmtPulse.size() );
  convolve( m_pmtPulse_noOvershoot, m_shapingResponse, m_shapedPmtPulse_noOvershoot );
  m_shapedPmtPulse_noOvershoot.resize( m_pmtPulse.size() );
  // Scale the resulting convolution to preserve the pulse area
  for(unsigned int i=0; i < m_pmtPulse.size(); i++) {
    m_shapedPmtPulse[i] /= m_shapingSum;
    m_shapedPmtPulse_noOvershoot[i] /= m_shapingSum;
  }
  // Convolve the ringing with the CR-(RC)^4 shaping response
  convolve( m_ringing, m_shapingResponse, m_shapedRinging );
  m_shapedRinging.resize( m_ringing.size() );
  // Scale the resulting convolution to preserve normalization
  for(unsigned int i=0; i < m_shapedRinging.size(); i++) 
    m_shapedRinging[i] /= m_shapingSum;

  // Dump ideal pulse shapes
  {
    std::stringstream output;
    output << "pmtPulse = [ ";
    for(int i=0; i<nPulseSamples; i++){
      if(i!=0) output << ", ";
      output << m_pmtPulse[i];
    } 
    output << " ]";
    verbose() << output.str() << endreq;
  }
  
  {
    std::stringstream output;
    output << "Esum Shaping Pulse = [ ";
    for(int i=0; i<nPulseSamples; i++){
        if(i!=0) output << ", ";
    output << m_esumResponse[i];
    }
    output << " ]";
    verbose() << output.str() << endreq;
  }
  
  {
    std::stringstream output;
    output << "ringing = [ ";
    for(int i=0; i<nRingingSamples; i++){
      if(i!=0) output << ", ";
      output << m_ringing[i];
    }
    output << " ]";
    verbose() << output.str() << endreq;
  }
  {
    std::stringstream output;
    output << "shapedPmtPulse = [ ";
    for(int i=0; i<nPulseSamples; i++){
      if(i!=0) output << ", ";
      output << m_shapedPmtPulse[i];
    } 
    output << " ]";
    verbose() << output.str() << endreq;
  }
  {
    std::stringstream output;
    output << "shapedRinging = [ ";
    for(int i=0; i<nRingingSamples; i++){
      if(i!=0) output << ", ";
      output << m_shapedRinging[i];
    }
    output << " ]";
    verbose() << output.str() << endreq;
  }

  return StatusCode::SUCCESS;
}

StatusCode EsIdealFeeTool::mapPulsesByChannel(
            const DayaBay::ElecPulseCollection::PulseContainer& pulses,
            std::map<DayaBay::ElecChannelId, std::vector<DayaBay::ElecPmtPulse*> >& pulseMap,
            int* hitCountSlices)
{
  debug() << "Processing " << pulses.size() << " pmt pulses." << endreq;
  DayaBay::ElecPulseCollection::PulseContainer::const_iterator it, 
    pulseDone = pulses.end();
  int timeSliceNo;
  // loop over pulses in collection
  for (it=pulses.begin(); it != pulseDone; ++it) {
    DayaBay::ElecPmtPulse* pulse = dynamic_cast<DayaBay::ElecPmtPulse*>(*it);
    if(pulse == NULL){
      // Catch invalid pulses
      error() << "mapPulsesByChannel(): bad pulse." << endreq;
      return StatusCode::FAILURE;
    }
    (pulseMap[pulse->channelId()]).push_back(pulse);
    timeSliceNo = int( pulse->time()/150. - 0.5 );
    hitCountSlices[timeSliceNo]+=int(pulse->amplitude()+0.5);
  }
  return StatusCode::SUCCESS;
}

StatusCode EsIdealFeeTool::updateBoardSignals(DayaBay::ElecFeeCrate* crate)
{ 
  // Update the board Nhit and Esum info, according to the current
  // channel info
  
  DayaBay::ElecFeeCrate::DigitalMap& nhitMap = crate->nhit();
  DayaBay::ElecFeeCrate::AnalogMap& esumMap = crate->esum();

  // Loop over channels, and add data to the board-level signals
  const DayaBay::ElecFeeCrate::ChannelData& channelData = crate->channelData();
  DayaBay::ElecFeeCrate::ChannelData::const_iterator channelIter, 
    done = channelData.end();
  for(channelIter = channelData.begin(); channelIter != done; ++channelIter){
    DayaBay::FeeChannelId channelId = channelIter->first;
    const DayaBay::ElecFeeChannel& channel = channelIter->second;
    const DayaBay::DigitalSignal& hit = channel.hit();
    const DayaBay::AnalogSignal& energy = channel.energy();
    DayaBay::DigitalSignal& boardNhit = nhitMap[channelId.boardId()];
    DayaBay::AnalogSignal& boardEsum = esumMap[channelId.boardId()];
    
    unsigned int hitSamples = hit.size();

    // Create hitSync to mimic the syncronized discriminator output
    //DayaBay::DigitalSignal hitSync( hitSamples,0);
    if( m_hitSync.size() != hitSamples ) m_hitSync.resize( hitSamples );
    const int* hitStart = &hit[0];
    int* hitSyncStart = &m_hitSync[0];
    *hitSyncStart = 0;  // Set first sample to zero
    for (unsigned int i=1; i< hitSamples; i++){
      int* curHitSync = hitSyncStart + i;
      int* lastHitSync = curHitSync - 1;
      const int* curHit = hitStart + i;
      const int* lastHit = curHit - 1;
      if( *lastHitSync==0 && *lastHit==0 && *curHit==1 ) *curHitSync = 1;
      else *curHitSync = 0;
    }
    
    // The logic pulse from the channel level hit is extended for
    // multiple clock cycles (triggerWindowCycles) before going to the
    // trigger.
    //DayaBay::DigitalSignal hitHold( hitSamples , 0 );
    if( m_hitHold.size() != hitSamples ) m_hitHold.resize( hitSamples );
    int* hitHoldStart = &m_hitHold[0];
    for (unsigned int i=0; i< hitSamples; i++) *(hitHoldStart + i) = 0;
    
    if( boardNhit.size() != hitSamples ){
      boardNhit.resize( hitSamples );
      int* boardNhitStart = &boardNhit[0];
      for (unsigned int i=0; i< hitSamples; i++) *(boardNhitStart + i) = 0;
    }
      
    for( unsigned int idx = 0; idx < hitSamples; idx++ ){
      unsigned int tempIdx = 0;
      if( *(hitSyncStart+idx) != 0 ){
    for( int holdIdx = 0; holdIdx < m_triggerWindowCycles; holdIdx++ ){
      unsigned int currentIdx = idx + holdIdx;
      if( currentIdx < hitSamples ) *(hitHoldStart + currentIdx) = 1;
      tempIdx++;
    }
    if(tempIdx !=0) idx = idx +tempIdx-1 ;
      } // End hit
    }

    {
      int* boardNhitStart = &boardNhit[0];
      for(unsigned int idx = 0; idx < hitSamples; idx++)
    *(boardNhitStart + idx) += *(hitHoldStart + idx);
    }
    verbose() << "boardNhit " << boardNhit << endreq;
    
    // Add channel energy to board esum
    unsigned int esumSamples = energy.size();
    if( boardEsum.size() != esumSamples ){
      boardEsum.resize( esumSamples );
      double* boardEsumStart = &boardEsum[0];
      for (unsigned int i=0; i< esumSamples; i++) *(boardEsumStart + i) = 0;
    }
      const double* energyStart = &energy[0];
      double* boardEsumStart = &boardEsum[0];
      for( unsigned int idx = 0; idx < esumSamples; idx++ ){
        *(boardEsumStart + idx) += *(energyStart + idx);
      }
      
     // warning() << "board esum " << boardEsum << endreq;
     // warning() << "energy " << energy << endreq;

  } // End loop over channels

  return StatusCode::SUCCESS;
}

StatusCode EsIdealFeeTool::updateEsumSignals(DayaBay::ElecFeeCrate* crate){
  // Loop over all boards in map
  DayaBay::ElecFeeCrate::AnalogMap& esumMap = crate->esum();
  DayaBay::ElecFeeCrate::AnalogMap::iterator boardIt = esumMap.begin();
  
  // Clear transient stores from previous Execution cycles
  m_esumH.clear();
  m_esumL.clear();
  m_esumTotal.clear();
  
  for (;boardIt != esumMap.end(); ++boardIt)
  {
      DayaBay::FeeChannelId boardId = boardIt->first;
      DayaBay::AnalogSignal& boardEnergy = boardIt->second;
      int boardNum = boardId.board();

      if( m_esumL.size() != boardEnergy.size() ){
          m_esumL.resize( boardEnergy.size() );
      }
      
      if( m_esumH.size() != boardEnergy.size() ){
        m_esumH.resize( boardEnergy.size() );
      }
      
      if( m_esumTotal.size() != boardEnergy.size() ){
        m_esumTotal.resize( boardEnergy.size() );
      }

      debug() << "board number " << boardNum << endreq;
      debug() << "Low Esum size " << m_esumL.size() << endreq;
      debug() << "board Energy size " << boardEnergy.size() << endreq;
      
      int boardIdx = boardEnergy.size();
      double* boardEsumStart = &boardEnergy[0];
      double* esumHStart = &m_esumH[0];
      double* esumLStart = &m_esumL[0];
      double* esumTtart  = &m_esumTotal[0];
      
      for( int i = 0; i < boardIdx; i++)
      {
          // update Upper Sum and total sum
          if (boardEsumType(boardNum) == DayaBay::ESumComp::kESumHigh){
              *(esumHStart+i) += *(boardEsumStart+i);
              *(esumTtart+i) += *(boardEsumStart+i);
              //m_esumH[i] += boardEnergy[i];
              //m_esumTotal[i] += boardEnergy[i];
          }
          
          // update Lower Sum and total sum
          if (boardEsumType(boardNum) == DayaBay::ESumComp::kESumLow){
              *(esumLStart+i) += *(boardEsumStart+i);
              *(esumTtart+i) += *(boardEsumStart+i);
              //m_esumL[i] += boardEnergy[i];
              //m_esumTotal[i] += boardEnergy[i];
          }
      }

  
  }// END LOOP over boards
    
  
  verbose() << "Lower ESum " << m_esumL << endreq;
  verbose() << "Total ESum " << m_esumTotal << endreq;
    
  // Shape Upper Sum
  if(m_enableESumH) {
    verbose() << "Upper ESum " << m_esumH << endreq;
    shape( m_esumH, m_shapedEsumH, m_esumResponse, m_esumConvolutionWindows );
    m_shapedEsumH.resize( m_esumH.size() );
  
    // NOTE: They should all be the same lenght so if we need to speed things up
    //        we could resacle all 3 in the same loop...
    
    // rescale after convolution
    double* esumHStart = &m_shapedEsumH[0];
    
    for(unsigned int i=0; i < m_shapedEsumH.size(); i++) 
    {  
      *(esumHStart+i) /= m_esumShapingSum;
      //m_shapedEsumH[i] /= m_esumShapingSum;  
    }
  
    // sample down to appropriate frequency and add to header 
    sample( m_shapedEsumH, crate->esumUpper(),
             m_simFrequency, m_eSumFrequency );
    
    verbose() << "Shaped upper sum " << crate->esumUpper() << endreq;
  }
  else {
    m_shapedEsumH.resize( m_esumH.size() );
    crate->setEsumUpper(m_shapedEsumH);
  }
  // Shape Lower Sum
  if(m_enableESumL) {
    shape( m_esumL, m_shapedEsumL, m_esumResponse, m_esumConvolutionWindows );
    m_shapedEsumL.resize( m_esumL.size() );
  
    // rescale after convolution
    double* esumLStart = &m_shapedEsumL[0];
    
    for(unsigned int i=0; i < m_shapedEsumL.size(); i++)
    {
        *(esumLStart+i) /= m_esumShapingSum;
        //m_shapedEsumL[i] /= m_esumShapingSum;
    }
    
    // sample down to appropriate frequency and add to header  
    sample( m_shapedEsumL, crate->esumLower(),
             m_simFrequency, m_eSumFrequency );
      
    verbose() << "Shaped lower sum " << crate->esumLower() << endreq;
  }
  else {
    m_shapedEsumL.resize( m_esumL.size() );
    crate->setEsumLower(m_shapedEsumL);
  }
  // Shape Total Sum 
  if(m_enableESumTotal) {
    shape( m_esumTotal, m_shapedEsumTotal, m_esumResponse, m_esumConvolutionWindows );
    m_shapedEsumTotal.resize( m_esumTotal.size() );
    
    // rescale after convolution
    double* esumTstart  = &m_shapedEsumTotal[0];
     
    for(unsigned int i=0; i < m_shapedEsumTotal.size(); i++) 
    {
       *(esumTstart+i) /= m_esumShapingSum;
       //m_shapedEsumTotal[i] /= m_esumShapingSum; 
    }  
  
    // sample down to appropriate frequency and add to header  
    sample( m_shapedEsumTotal, crate->esumTotal(),
             m_simFrequency, m_eSumFrequency );
        
    verbose() << "Shaped total sum " << crate->esumTotal() << endreq;
  }
  else {
    m_shapedEsumTotal.resize( m_esumTotal.size() );
    crate->setEsumTotal(m_shapedEsumTotal);
  }
  //Digitize the shaped total sum
  //FIXME: This needs to be fixed (edraeger@iit.edu)
  //         set as properties and/or get from service
  digitize(crate->esumTotal(), crate->esumADC(), m_adcRange, m_adcBits, m_adcOffset);
  
  verbose() << "Digitized signal " << crate->esumADC() << endreq;
  
  return StatusCode::SUCCESS;
}


DayaBay::ESumComp::ESumComp_t EsIdealFeeTool::boardEsumType(int boardId){
  // FIXME: This should probably be moved into a service.
  //         Probably into the CableMappingSvc
  if (boardId >= 1 && boardId < 7)
    {
      return DayaBay::ESumComp::kESumLow;
    }

  if (boardId >= 7 && boardId < 13)
    {
      return DayaBay::ESumComp::kESumHigh;
    }

  if (boardId == 13)
    {
      return DayaBay::ESumComp::kESumNone;
    }
      
  return DayaBay::ESumComp::kUnknown;
}


void EsIdealFeeTool::convolve(const DayaBay::AnalogSignal& signal,
                              const DayaBay::AnalogSignal& response,
                  DayaBay::AnalogSignal& output) 
{
  // Primitive convolution of signal and response, to produce
  // output.  Assumes values outside of range are zero.  A few simple
  // optimizations are use to speed up the process.
  int sN = signal.size();
  int rN = response.size();
  output.resize(sN+rN);
  const double* sigD = &(signal[0]);
  const double* resD = &(response[0]);
  double* outD = &(output[0]);
  for (int i=0; i<(sN+rN); i++) {
    double sum = 0;
    for (int j=0; j<rN; j++) {
      if (resD[j]==0)
        continue;
      if (i-j>=0 && i-j<sN)
        sum+=sigD[i-j]*resD[j];
    }
    outD[i]=sum;
  }
}

// Only convolve time intervalls with minimum number of hits
void EsIdealFeeTool::shape(const DayaBay::AnalogSignal& rawSignal,
                DayaBay::AnalogSignal& shapedSignal, 
                const DayaBay::AnalogSignal& response, std::map<int,int> convolutionWindows){
  std::map<int,int>::iterator convIt = convolutionWindows.begin();
  int convStart, convEnd;
  shapedSignal.clear();
  int samples = rawSignal.size();
  shapedSignal.resize(samples, 0);
  DayaBay::AnalogSignal rawFrag;
  DayaBay::AnalogSignal shapedFrag;
  for(; convIt != convolutionWindows.end(); convIt++){
    convStart = convIt->first * 96;
    convEnd = convIt->second * 96;
    if (convEnd > samples) convEnd = samples;
    rawFrag.resize( (convEnd - convStart) );
    shapedFrag.resize( (convEnd - convStart) );
    std::copy( rawSignal.begin() + convStart, rawSignal.begin() + convEnd , rawFrag.begin() );
    convolve( rawFrag, response, shapedFrag );
    shapedFrag.resize( (convEnd - convStart) );
    std::copy( shapedFrag.begin(), shapedFrag.end(), shapedSignal.begin() + convStart );
    debug() << "Convolve from " << convStart *1.5625 << " ns to " << convEnd * 1.5625 << " ns" << endreq;
  }
}
void EsIdealFeeTool::sample(const DayaBay::AnalogSignal& signal,
                DayaBay::AnalogSignal& output,
                            int inputFrequency,
                            int outputFrequency)
{
  // Fill the output array using the signal values sampled at a
  // regular frequency
  if(inputFrequency == outputFrequency){
    if( output.size() != signal.size() ) output.resize( signal.size() );
    memcpy(&output[0],&signal[0],output.size()*sizeof(double));
    return;
  }
  unsigned int nSamples = convertClock(signal.size(), 
                       inputFrequency, 
                       outputFrequency);
  if(output.size() != nSamples){
    if(output.capacity() > nSamples){
      output.resize(nSamples);
    }else{
      // Faster memory allocation
      DayaBay::AnalogSignal* newSignal = new DayaBay::AnalogSignal(nSamples);
      output.swap(*newSignal);
      delete newSignal;
    }
  }

  if(inputFrequency < outputFrequency){
    warning() << "sample(): Input frequency " << inputFrequency 
          << " is less than the sampling frequency " << outputFrequency
          << endreq;
  }
  if( (inputFrequency % outputFrequency) == 0 ){
    // Check for integer relationship between frequencies
    int relativeFrequency = inputFrequency / outputFrequency;
    double* outputStart = &output[0]; 
    const double* inputStart = &signal[0]; 
    int inputIdx = 0;
    unsigned int outputIdx = 0;
    while(outputIdx != nSamples){
      *(outputStart + outputIdx) = *(inputStart + inputIdx);
      outputIdx++;
      inputIdx+=relativeFrequency;
    }
    return;
  }else{
    // Frequencies not coherent.  Use signal sample with nearest lower index.
    double* outputStart = &output[0];
    const double* inputStart = &signal[0];
    int inputIdx = 0;
    unsigned int outputIdx = 0;
    double freqScale = inputFrequency / outputFrequency;
    while( outputIdx != nSamples ){
      inputIdx = int(outputIdx * freqScale);
      *(outputStart + outputIdx) = *(inputStart + inputIdx);
      outputIdx++;
    }
  }
}

void EsIdealFeeTool::digitize(const DayaBay::AnalogSignal& signal,
                              DayaBay::DigitalSignal& output,
                              double range, int bits, double offset) 
{
  // Simple digitization algorithm
  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);
    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;
    */
  }
}

void EsIdealFeeTool::discriminate(const DayaBay::AnalogSignal& signal,
                  double threshold,
                  std::vector<int>& output,
                  DayaBay::Threshold::Threshold_t mode)
{
  // Take input signal array and discriminate using threshold.  Return
  // vector of indices which satisfy the discriminator.
  //  Current Modes:
  //    - Threshold crossed
  //    - Above threshold
  threshold *= m_speAmp;
  threshold *= m_gainFactor;
  threshold *= m_discThreshScale;
  
  output.clear();
  const double* signalStart = &signal[0];
  unsigned int nSamples = signal.size();
  if (DayaBay::Threshold::kCrossing == mode) {
    bool aboveThreshold = false;
    for (unsigned int i=0;
         i < nSamples;
         ++i) {
      if (!aboveThreshold && (*(signalStart+i))>=threshold) {
        output.push_back(i);
        aboveThreshold = true;
      } else if (aboveThreshold && (*(signalStart+i))<threshold)
        aboveThreshold=false;
    }
  } else if (DayaBay::Threshold::kAbove == mode) {
    for (unsigned int i=0;
         i < nSamples;
         ++i)
      if ((*(signalStart+i))>=threshold)
        output.push_back(i);
  } else {
    error() << "discriminate(): Unknown Discrimination Mode" << endreq;
  }
}

double EsIdealFeeTool::pmtPulse(double deltaT, int nPulse, bool addOvershoot) {
  // SJ: New pmt pulse parametrization, cf DocDB 3779
  double shift = -2.5e-9; // time shift parameter
  if (deltaT-shift<0) return 0.;
  double v0 = -m_speAmp; // single p.e. pulse amplitude in V
  double width; // pulse width parameter
  double mu; // parameter determining the degree of asymmetry
  if ( nPulse != 0 ){
    width = getPulseWidth( nPulse); 
    mu    = getPulseForm( nPulse);
  }
  else {
    width = 7.5e-9;
    mu    = 0.5;
  }
  if(deltaT > 100 * width) return 0;
  double amplitude = 0;
  if(deltaT < 10 * width){ 
    // Add primary pulse 
    amplitude += m_gainFactor* v0 * exp( -pow( log( (deltaT-shift)/width),2) 
               / ( 2*pow(mu,2) ) ) * Gaudi::Units::volt;
  }
  if( addOvershoot ){
    // Add overshoot 
    amplitude += overshoot( deltaT, nPulse );
  }
  return amplitude;
}

double EsIdealFeeTool::pulseShaping(double deltaT) {
  // Gaussian CR-(RC)^4 pulse shaping
  // See Knoll, Radiation Detection and Measurement
  // Returns the pulse amplitude at time deltaT after a positive step
  // function of amplitude 1
  if (deltaT<0)
    return 0.;
  double t0 = 25.0e-9; // RC time constant in seconds
  double r = deltaT/t0;
  return pow( r, 4 ) * exp( -r );
}

double EsIdealFeeTool::esumShaping(double deltaT) {
  //ESum shaping function
  if (deltaT<0)
    return 0.;
  double t0 = 56.0; // RC time constant in seconds
  double r = deltaT/t0;
  return exp( -r )*(1/t0);
}

double EsIdealFeeTool::noise() {
  // FIXME: do discrete FT for range of frequencies
  // Simple noise model
  // Approximate Gaussian-distributed noise
  return m_gauss() * m_noiseAmp; // no offset
}

int EsIdealFeeTool::convertClock(int clock, int inputFrequency, int outputFrequency){
  // Convert a clock cycle from one frequency to another.  Round down
  // for incoherent frequencies.
  return int( ( clock / double(inputFrequency) ) * outputFrequency );
}

double EsIdealFeeTool::satFactor(double amp){
//SJ: Returns saturation factor as a function of the raw amplitude
  if (amp<0) return 1.;
  double y = 1.e6 * amp;
  if(y < 1) 
    return 1. - 0.5311049 * y + 0.1818968 * pow(y,2);
  if(y < 9.5) 
    return 0.778264 - 0.1400074 * y + 0.01303484 * pow(y,2) - 0.0004701175 * pow(y,3) + 0.000002450386 * pow(y,4);
  return 0.4212711 - 2.560691e-02 * y + 8.003103e-04 * pow(y,2) - 1.226207e-05 * pow(y,3) + 7.282761e-08 * pow(y,4);
}

double EsIdealFeeTool::overshoot(double deltaT, int nPulse) {
  // SJ: Simple overshoot parametrization
  if (deltaT < 0) return 0.;
  double v0  = getOvershootAmp( nPulse ); // Overshoot amplitude in V 
  double tau = 155.e-9; // Overshoot time constant in s
  double expo = v0*exp(-deltaT/tau);
  double t0   = 50e-9; // Onset parameters
  double t1   = 10e-9;
  double fermi = 1. / (exp( (t0 - deltaT) / t1) + 1.);
  return fermi * expo * Gaudi::Units::volt;
}

double EsIdealFeeTool::ringing(double deltaT, double pulseMax) {
  // SJ: Simple ringing parametrization
  if (deltaT<0) return 0.;
  double v0 = 0.00154 * pulseMax; // Ringing amplitude in V
  double period =  1.63e-6; //Ringing period in s
  double shift= 0.99e-6; // Phase shift
  double expo = exp(-deltaT/12.e-6);
  double t0 = 0.6e-6; // Onset parameters
  double t1 = 0.1e-6;
  double fermi = 1./(exp((t0-deltaT)/t1)+1);
  return v0 * fermi * expo * cos(6.28*(deltaT-shift)/period);
}

double EsIdealFeeTool::getPulseWidth(int nPulse){
  // SJ: Single pulse width as a function of the number of coincidental pulses
  if( nPulse > 2000) return 20.7842e-9;
  double x  = double(nPulse);
  double t0 = 1200.;
  double t1 = 200.;
  double f  = 1. / (exp( (t0-x)/t1 ) + 1.);
  return (7.+12.*f+x/1000.)*1.e-9;
}

double EsIdealFeeTool::getPulseForm(int nPulse){
  // SJ: Single pulse width as a function of the number of coincidental pulses
  if( nPulse > 3000) return 0.2;
  double x = double(nPulse);
  double t0 = 900.;
  double t1 = 200.;
  double f  = 0.45 - 0.25 * 1. / (exp( (t0-x)/t1) + 1.);
  return f;
}

double EsIdealFeeTool::getOvershootAmp(int nPulse){
  // SJ: Single pulse overshoot amplitude as a function of the number of coincidental pulses
  if( nPulse > 2000) return 2.4033e-4;
  double x  = double(nPulse);
  double t0 = -0.0148227;
  double t1 = -1186.36;
  double t2 =  287.821;
  double t3 =  0.0150628;
  double f  = t0 / (exp( (t1-x)/t2) + 1.) + t3;
  return f;
}

---