EsIdealFecTool.cc

Go to the documentation of this file.

#include "EsIdealFecTool.h"

#include "Event/ElecFecCrate.h"
#include "Event/ElecPulseCollection.h"
#include "Event/ElecRpcPulse.h"
#include "GaudiKernel/SystemOfUnits.h"
#include "GaudiKernel/IRndmGenSvc.h"

using namespace std;
using namespace Gaudi;
#define PulseThreshold 0.5

EsIdealFecTool::EsIdealFecTool(const std::string& type,
                               const std::string& name,
                               const IInterface* parent)
    : GaudiTool(type,name,parent)
{
    declareInterface< IEsFrontEndTool >(this) ;

    declareProperty("CableSvcName",m_cableSvcName="CableSvc",
          "Name of service to map between detector, hardware, and electronic IDs");
    declareProperty("SimDataSvcName",m_simDataSvcName="SimDataSvc",
          "Name of service to provide FEC channel properties for simulation");

    declareProperty("SimFrequency", m_simFrequency = DayaBay::BaseFrequencyHz * Units::hertz,
                    "Simulation Frequency for RPCs. Defaults to DayaBay::BaseFrequencyHz");
    //declareProperty("SimFrequency", m_simFrequency = 40 * Units::megahertz,
                    //"Simulation Frequency for RPCs. Defaults to DayaBay::TdcFrequencyHz");
    declareProperty("IdealNoise", m_idealNoise = 0, "Noise frequency per RPC strip");
    declareProperty("IdealEff", m_idealEff = 1, "RPC Efficiency, common for all RPC strips");
    declareProperty("ElecDelay", m_elecDelay = 0., "Signal delay inside the electronics");
    declareProperty("ElecStretch", m_elecStretch = 3, "Signal stretch after threshold crossing. In clock cycles.");
}

EsIdealFecTool::~EsIdealFecTool(){}

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

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

    // Get random generators
    IRndmGenSvc *rgs = 0;
    if (service("RndmGenSvc",rgs,true).isFailure()) {
        fatal() << "Failed to get random service" << endreq;
        return StatusCode::FAILURE;
    }
    if (m_uni.initialize(rgs, Rndm::Flat(0,1)).isFailure()) {
        fatal() << "Failed to initialize uniform random numbers" << endreq;
        return StatusCode::FAILURE;
    }

  //StatusCode sc = loadResponse();
  //if(sc.isFailure()) return sc;
  debug()<< "m_simFrequency[ns-1] = "<<m_simFrequency<<endreq;
    //<<", DayaBay::TdcFrequencyHz = "<<DayaBay::TdcFrequencyHz<<endreq;

  return StatusCode::SUCCESS;
}

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

StatusCode EsIdealFecTool::generateSignals(DayaBay::ElecPulseCollection* pulses,
                                           DayaBay::ElecCrate* elecCrate)
{
    if ( msgLevel(MSG::DEBUG) )
        debug() << "running generateSignals()" << endreq;


    //Ensure the crate is a FEC crate
    DayaBay::ElecFecCrate* crate = 0;
    crate = dynamic_cast<DayaBay::ElecFecCrate*>(elecCrate);
    if(crate == NULL){
        error() <<"generateSignals(): Crate is not a FEC 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->fecBoards().size() == 0) {
        DayaBay::ElecFecCrate::FecBoardMap brds;
        // First, get list of all connected FEC channels
        const std::vector<DayaBay::DetectorSensor>& rpcList
          = m_cableSvc->sensors( svcMode );
        std::vector<DayaBay::DetectorSensor>::const_iterator rpcIt, rpcDone = rpcList.end();
        for(rpcIt = rpcList.begin(); rpcIt != rpcDone; rpcIt++) {
          // Add channel to the corresponding FEC board
          DayaBay::FecChannelId channelId(m_cableSvc->elecChannelId(*rpcIt, svcMode).fullPackedData());
          DayaBay::FecChannelId boardId(channelId.boardId());
          if (brds.find(boardId) == brds.end()) // add this board since it does not exist yet
            brds[boardId].setBoardId(boardId);
          brds[boardId].addChannel(channelId); // this create a channel in the board's channel map
          if ( msgLevel(MSG::VERBOSE) )
            verbose() << "generateSignals(): Adding channel \""<<channelId
                <<"\" to the board \""<<channelId.boardId()<<"\""<< endreq;
        }
        // set the boards for the crate
        crate->setFecBoards(brds);
    }

    // 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) * Units::second;
    verbose() << "generateSignals(): Simulation time is " << simTime / Units::nanosecond <<" ns"<< endreq;
    // The number of samples in time given the simulation frequency
    int simSamples = int(simTime * m_simFrequency);
    verbose() << "generateSignals(): Simulation has " << simSamples << " samples at frequency "
            << m_simFrequency/Units::megahertz << " MHz" << endreq;

    // Calculate offset of the 0th clock tick.
    // tickTime - time in between ticks
    // residualTime - time from the last tick to the time of earliest hit (in nano seconds)
    // clockShift - shift of the 0th tick relative to the earliest hit
    // clockOffset - time from the ElecHeader's time and the 0th tick
    // I assume here the clock is synchronized to the beginning of a second.
    double tickTime = Units::nanosecond/m_simFrequency;
    int residualTime = int(earliest.GetNanoSec()) % int(tickTime/Units::nanosecond);
    double clockShift = 0.;
    if (residualTime > 0)
        clockShift = tickTime - residualTime * Units::nanosecond;
    double clockOffset = (context.GetTimeStamp() - earliest)*Units::second + clockShift;
    if ( msgLevel(MSG::DEBUG) ) {
        debug() << "genrateSignals(): clockShift = " << clockShift
            << ", clockOffset = " << clockOffset
            << ", tickTime = " << tickTime << endreq;
    }

  // Organize Pulses by board and Channel
  PulseMap pulseMap;
  StatusCode sc;
  sc = mapPulsesByChannel( pulses->pulses(), pulseMap );
  if(sc.isFailure()) {
    return sc;
  }
  debug() << "generateSignals():  Pulses Mapped" << endreq;

    // Fill template signal samples
    DayaBay::DigitalSignal tmpSignal(simSamples, 0);

    // noise to be simulated check
    if ( m_idealNoise > 0. ) {
        if ( msgLevel(MSG::DEBUG) ) {
            debug()<<"generateSignals(): Noise will be mixed to the signals"
                <<", noise rate is set to "<<m_idealNoise/Units::kilohertz<<" kHz."<<endreq;
        }
    }

    const DayaBay::ElecFecCrate::FecBoardMap& boards = crate->fecBoards();
    // loop over boards in the crate
    DayaBay::ElecFecCrate::FecBoardMap::const_iterator brdIt,
                    brdDone = boards.end();
    for(brdIt = boards.begin(); brdIt != brdDone; ++brdIt){
        DayaBay::FecChannelId boardId = brdIt->first;
        DayaBay::ElecFecBoard& board = crate->board(boardId);
        std::map<DayaBay::FecChannelId, std::vector<DayaBay::ElecRpcPulse*> >& pulseBrd = pulseMap[boardId];

        // Loop over channels, and fill with the simulated signals
        DayaBay::ElecFecBoard::ChannelMap::const_iterator chIt, chDone = board.channels().end();
        for (chIt = board.channels().begin(); chIt != chDone; chIt++) {
            const DayaBay::FecChannelId& channelId = chIt->first;
            DayaBay::ElecFecChannel& channel = board.channel(channelId);
            channel.setHits(tmpSignal); // initialize digital signal in the channel

            std::vector<DayaBay::ElecRpcPulse*>& pulses = pulseBrd[channelId];
            if(!pulses.empty()){  //Hits
                StatusCode sc = generateOneChannel(pulses, channel, clockOffset, svcMode);
                if (sc.isFailure()) {
                    error()<<"generateSignals(): Cannot generate pulses in the channel "<<channelId<<endreq;
                    return StatusCode::FAILURE;
                }
            }
            if ( m_idealNoise > 0. ) {
                // add noise to the signal
                if ( msgLevel(MSG::VERBOSE) ) {
                    verbose()<<"generateSignals(): Noise will be mixed to the signals."
                        <<"Noise rate is set to "<<m_idealNoise/Units::kilohertz
                        <<" kHz for channel "<<channelId<<endreq;
                }
                StatusCode sc = addChannelNoise(channelId, channel, simTime, svcMode);
                if (sc.isFailure()) {
                    error()<<"Failed to generate noise in the channel "<<channelId<<endreq;
                    return StatusCode::FAILURE;
                }
            } else {
                // no noise to be simulated
                if ( msgLevel(MSG::VERBOSE) ) {
                    verbose()<<"generateSignals(): No noise will be mixed to the signals in this channel, "
                        << channelId
                        <<", noise rate is set to "<<m_idealNoise<<endreq;
                }
            }
        }
    }
    debug() << "generateSignals(): Signals generated" << endreq;
    verbose() << "generateSignals(): Crate\n" << *crate << endreq;

    return StatusCode::SUCCESS;
}

StatusCode EsIdealFecTool::mapPulsesByChannel(
            const DayaBay::ElecPulseCollection::PulseContainer& pulses,
            PulseMap& pulseMap)
{
  debug() << "mapPulsesByChannel(): Processing " << pulses.size() << " rpc 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::ElecRpcPulse* pulse = dynamic_cast<DayaBay::ElecRpcPulse*>(*it);
    if(pulse == NULL){
      // Catch invalid pulses
      error() << "mapPulsesByChannel(): bad pulse." << endreq;
      return StatusCode::FAILURE;
    }
    DayaBay::FecChannelId chId(pulse->channelId().fullPackedData());
    pulseMap[chId.boardId()][chId].push_back(pulse);
  }
  return StatusCode::SUCCESS;
}

StatusCode EsIdealFecTool::generateOneChannel(std::vector<DayaBay::ElecRpcPulse*>& channelPulses,
                    DayaBay::ElecFecChannel& channel,
                    double clockOffset,
                    const ServiceMode&)
{
    DayaBay::DigitalSignal& digitSignal = channel.hits();
    std::vector<int>& indexes = channel.hitIndexes();
    // loop over all pulses in the channel
    int size = channelPulses.size();
    for (int i = 0; i < size; i++) {
        DayaBay::ElecRpcPulse* pulse = channelPulses[i];
        // time of the pulse relative to the 0th clock tick (earliest + clockShift)
        double offsetT = pulse->time() - clockOffset;
        // which clock tick corresponds to the time
        int clock = int ( offsetT * m_simFrequency + 1. );
        if ( clock >= digitSignal.size() ) {
            error()<<"generateOneChannel(): Calculated time of the pulse "<<i<<", "
                <<offsetT/Units::nanosecond
                <<" ns, is out of the simulation time window"<<endreq;
            return StatusCode::FAILURE;
        }

        if( (pulse->amplitude()) > 0.5 ) {
            // pulse crossed the threshold (virtual)
            double ran = m_uni();
            if(ran < m_idealEff){
                for (int j = 0; j < m_elecStretch; j++) { // stretch the incoming signal
                    digitSignal[clock+j] = 1;
                    indexes.push_back(clock+j);
                }
            }
        }
    }
    return StatusCode::SUCCESS;
}

StatusCode EsIdealFecTool::addChannelNoise(const DayaBay::FecChannelId& channelId,
                    DayaBay::ElecFecChannel& channel,
                    double simTime,
                    const ServiceMode&)
{
    DayaBay::DigitalSignal& digitSignal = channel.hits();
    std::vector<int>& indexes = channel.hitIndexes();

    // Calculate mean number of hits in the time window
    double nhitsMean = simTime * m_idealNoise;
    // Maximum number of noise hits
    int nhitsMax = simTime * m_simFrequency;

    int nhits = getNoiseHits(nhitsMean, nhitsMax);
    if (nhits < 0) {
        error()<<"Cannot get number of noise hits for this channel: "
            << channelId << endreq;
        return StatusCode::FAILURE;
    }

    if ( msgLevel(MSG::VERBOSE) ) {
        verbose()<<"Generating "<<nhits<<" noise hits in channel "<<channelId<<endreq;
    }

    for (int i = 0; i < nhits; i++) {
        double ran = m_uni();
        int hit = ran * digitSignal.size();
        if (hit < digitSignal.size()) {
            for (int j = 0; j < m_elecStretch; j++) {
                digitSignal[hit + j] = 1;
                indexes.push_back(hit+j);
            }
        }
    }

    return StatusCode::SUCCESS;
}

int EsIdealFecTool::getNoiseHits(double mean, int max)
{
    if (mean < 0.) {
        fatal() << "Trying to generate integer from poisson distribution with negative parameter: "
            << mean << endreq;
        return -1;
    }
    if (mean == 0.)
        return 0;

    double ran = m_uni() * exp(mean);
    int nhits = 0;
    double p = 1.;
    double F = 1.;
    while (F < ran && nhits < max) {
        nhits++;
        p *= mean/nhits;
        F += p;
    }
    if (nhits >= max) {
        info() << "Number of noise hits was about to exceed maximum allowed. "
            << "setting it to the maximum: " << max << endreq;
        return max;
    }
    return nhits;
}