/search.css" rel="stylesheet" type="text/css"/> /search.js">
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

In This Package:

EsIdealFecTool.cc
Go to the documentation of this file.
00001 #include "EsIdealFecTool.h"
00002 
00003 #include "Event/ElecFecCrate.h"
00004 #include "Event/ElecPulseCollection.h"
00005 #include "Event/ElecRpcPulse.h"
00006 #include "GaudiKernel/SystemOfUnits.h"
00007 #include "GaudiKernel/IRndmGenSvc.h"
00008 
00009 using namespace std;
00010 using namespace Gaudi;
00011 #define PulseThreshold 0.5
00012 
00013 EsIdealFecTool::EsIdealFecTool(const std::string& type,
00014                                const std::string& name,
00015                                const IInterface* parent)
00016     : GaudiTool(type,name,parent)
00017 {
00018     declareInterface< IEsFrontEndTool >(this) ;
00019 
00020     declareProperty("CableSvcName",m_cableSvcName="CableSvc",
00021           "Name of service to map between detector, hardware, and electronic IDs");
00022     declareProperty("SimDataSvcName",m_simDataSvcName="SimDataSvc",
00023           "Name of service to provide FEC channel properties for simulation");
00024 
00025     declareProperty("SimFrequency", m_simFrequency = DayaBay::BaseFrequencyHz * Units::hertz,
00026                     "Simulation Frequency for RPCs. Defaults to DayaBay::BaseFrequencyHz");
00027     //declareProperty("SimFrequency", m_simFrequency = 40 * Units::megahertz,
00028                     //"Simulation Frequency for RPCs. Defaults to DayaBay::TdcFrequencyHz");
00029     declareProperty("IdealNoise", m_idealNoise = 0, "Noise frequency per RPC strip");
00030     declareProperty("IdealEff", m_idealEff = 1, "RPC Efficiency, common for all RPC strips");
00031     declareProperty("ElecDelay", m_elecDelay = 0., "Signal delay inside the electronics");
00032     declareProperty("ElecStretch", m_elecStretch = 3, "Signal stretch after threshold crossing. In clock cycles.");
00033 }
00034 
00035 EsIdealFecTool::~EsIdealFecTool(){}
00036 
00037 StatusCode EsIdealFecTool::initialize()
00038 {
00039   // Get Cable Service
00040   m_cableSvc = svc<ICableSvc>(m_cableSvcName,true);
00041 
00042   // Get RPC simulation input data service
00043   m_simDataSvc = svc<ISimDataSvc>(m_simDataSvcName,true);
00044 
00045     // Get random generators
00046     IRndmGenSvc *rgs = 0;
00047     if (service("RndmGenSvc",rgs,true).isFailure()) {
00048         fatal() << "Failed to get random service" << endreq;
00049         return StatusCode::FAILURE;
00050     }
00051     if (m_uni.initialize(rgs, Rndm::Flat(0,1)).isFailure()) {
00052         fatal() << "Failed to initialize uniform random numbers" << endreq;
00053         return StatusCode::FAILURE;
00054     }
00055 
00056   //StatusCode sc = loadResponse();
00057   //if(sc.isFailure()) return sc;
00058   debug()<< "m_simFrequency[ns-1] = "<<m_simFrequency<<endreq;
00059     //<<", DayaBay::TdcFrequencyHz = "<<DayaBay::TdcFrequencyHz<<endreq;
00060 
00061   return StatusCode::SUCCESS;
00062 }
00063 
00064 StatusCode EsIdealFecTool::finalize()
00065 {
00066   return StatusCode::SUCCESS;
00067 }
00068 
00069 StatusCode EsIdealFecTool::generateSignals(DayaBay::ElecPulseCollection* pulses,
00070                                            DayaBay::ElecCrate* elecCrate)
00071 {
00072     if ( msgLevel(MSG::DEBUG) )
00073         debug() << "running generateSignals()" << endreq;
00074 
00075 
00076     //Ensure the crate is a FEC crate
00077     DayaBay::ElecFecCrate* crate = 0;
00078     crate = dynamic_cast<DayaBay::ElecFecCrate*>(elecCrate);
00079     if(crate == NULL){
00080         error() <<"generateSignals(): Crate is not a FEC Crate."<< endreq;
00081         return StatusCode::FAILURE;
00082     }
00083 
00084     //Context, ServiceMode for this data
00085     Context context(pulses->detector().site(), SimFlag::kMC,
00086           pulses->header()->header()->timeStamp(),
00087           pulses->detector().detectorId());
00088     int task = 0;
00089     ServiceMode svcMode(context, task);
00090 
00091     // Initialize crate if necessary
00092     if(crate->fecBoards().size() == 0) {
00093         DayaBay::ElecFecCrate::FecBoardMap brds;
00094         // First, get list of all connected FEC channels
00095         const std::vector<DayaBay::DetectorSensor>& rpcList
00096           = m_cableSvc->sensors( svcMode );
00097         std::vector<DayaBay::DetectorSensor>::const_iterator rpcIt, rpcDone = rpcList.end();
00098         for(rpcIt = rpcList.begin(); rpcIt != rpcDone; rpcIt++) {
00099           // Add channel to the corresponding FEC board
00100           DayaBay::FecChannelId channelId(m_cableSvc->elecChannelId(*rpcIt, svcMode).fullPackedData());
00101           DayaBay::FecChannelId boardId(channelId.boardId());
00102           if (brds.find(boardId) == brds.end()) // add this board since it does not exist yet
00103             brds[boardId].setBoardId(boardId);
00104           brds[boardId].addChannel(channelId); // this create a channel in the board's channel map
00105           if ( msgLevel(MSG::VERBOSE) )
00106             verbose() << "generateSignals(): Adding channel \""<<channelId
00107                 <<"\" to the board \""<<channelId.boardId()<<"\""<< endreq;
00108         }
00109         // set the boards for the crate
00110         crate->setFecBoards(brds);
00111     }
00112 
00113     // Prepare time range of the simulation
00114     TimeStamp earliest = pulses->header()->header()->earliest();
00115     TimeStamp latest = pulses->header()->header()->latest();
00116     // The time window for the crate signals
00117     double simTime = (latest - earliest) * Units::second;
00118     verbose() << "generateSignals(): Simulation time is " << simTime / Units::nanosecond <<" ns"<< endreq;
00119     // The number of samples in time given the simulation frequency
00120     int simSamples = int(simTime * m_simFrequency);
00121     verbose() << "generateSignals(): Simulation has " << simSamples << " samples at frequency "
00122             << m_simFrequency/Units::megahertz << " MHz" << endreq;
00123 
00124     // Calculate offset of the 0th clock tick.
00125     // tickTime - time in between ticks
00126     // residualTime - time from the last tick to the time of earliest hit (in nano seconds)
00127     // clockShift - shift of the 0th tick relative to the earliest hit
00128     // clockOffset - time from the ElecHeader's time and the 0th tick
00129     // I assume here the clock is synchronized to the beginning of a second.
00130     double tickTime = Units::nanosecond/m_simFrequency;
00131     int residualTime = int(earliest.GetNanoSec()) % int(tickTime/Units::nanosecond);
00132     double clockShift = 0.;
00133     if (residualTime > 0)
00134         clockShift = tickTime - residualTime * Units::nanosecond;
00135     double clockOffset = (context.GetTimeStamp() - earliest)*Units::second + clockShift;
00136     if ( msgLevel(MSG::DEBUG) ) {
00137         debug() << "genrateSignals(): clockShift = " << clockShift
00138             << ", clockOffset = " << clockOffset
00139             << ", tickTime = " << tickTime << endreq;
00140     }
00141 
00142   // Organize Pulses by board and Channel
00143   PulseMap pulseMap;
00144   StatusCode sc;
00145   sc = mapPulsesByChannel( pulses->pulses(), pulseMap );
00146   if(sc.isFailure()) {
00147     return sc;
00148   }
00149   debug() << "generateSignals():  Pulses Mapped" << endreq;
00150 
00151     // Fill template signal samples
00152     DayaBay::DigitalSignal tmpSignal(simSamples, 0);
00153 
00154     // noise to be simulated check
00155     if ( m_idealNoise > 0. ) {
00156         if ( msgLevel(MSG::DEBUG) ) {
00157             debug()<<"generateSignals(): Noise will be mixed to the signals"
00158                 <<", noise rate is set to "<<m_idealNoise/Units::kilohertz<<" kHz."<<endreq;
00159         }
00160     }
00161 
00162     const DayaBay::ElecFecCrate::FecBoardMap& boards = crate->fecBoards();
00163     // loop over boards in the crate
00164     DayaBay::ElecFecCrate::FecBoardMap::const_iterator brdIt,
00165                     brdDone = boards.end();
00166     for(brdIt = boards.begin(); brdIt != brdDone; ++brdIt){
00167         DayaBay::FecChannelId boardId = brdIt->first;
00168         DayaBay::ElecFecBoard& board = crate->board(boardId);
00169         std::map<DayaBay::FecChannelId, std::vector<DayaBay::ElecRpcPulse*> >& pulseBrd = pulseMap[boardId];
00170 
00171         // Loop over channels, and fill with the simulated signals
00172         DayaBay::ElecFecBoard::ChannelMap::const_iterator chIt, chDone = board.channels().end();
00173         for (chIt = board.channels().begin(); chIt != chDone; chIt++) {
00174             const DayaBay::FecChannelId& channelId = chIt->first;
00175             DayaBay::ElecFecChannel& channel = board.channel(channelId);
00176             channel.setHits(tmpSignal); // initialize digital signal in the channel
00177 
00178             std::vector<DayaBay::ElecRpcPulse*>& pulses = pulseBrd[channelId];
00179             if(!pulses.empty()){  //Hits
00180                 StatusCode sc = generateOneChannel(pulses, channel, clockOffset, svcMode);
00181                 if (sc.isFailure()) {
00182                     error()<<"generateSignals(): Cannot generate pulses in the channel "<<channelId<<endreq;
00183                     return StatusCode::FAILURE;
00184                 }
00185             }
00186             if ( m_idealNoise > 0. ) {
00187                 // add noise to the signal
00188                 if ( msgLevel(MSG::VERBOSE) ) {
00189                     verbose()<<"generateSignals(): Noise will be mixed to the signals."
00190                         <<"Noise rate is set to "<<m_idealNoise/Units::kilohertz
00191                         <<" kHz for channel "<<channelId<<endreq;
00192                 }
00193                 StatusCode sc = addChannelNoise(channelId, channel, simTime, svcMode);
00194                 if (sc.isFailure()) {
00195                     error()<<"Failed to generate noise in the channel "<<channelId<<endreq;
00196                     return StatusCode::FAILURE;
00197                 }
00198             } else {
00199                 // no noise to be simulated
00200                 if ( msgLevel(MSG::VERBOSE) ) {
00201                     verbose()<<"generateSignals(): No noise will be mixed to the signals in this channel, "
00202                         << channelId
00203                         <<", noise rate is set to "<<m_idealNoise<<endreq;
00204                 }
00205             }
00206         }
00207     }
00208     debug() << "generateSignals(): Signals generated" << endreq;
00209     verbose() << "generateSignals(): Crate\n" << *crate << endreq;
00210 
00211     return StatusCode::SUCCESS;
00212 }
00213 
00214 StatusCode EsIdealFecTool::mapPulsesByChannel(
00215             const DayaBay::ElecPulseCollection::PulseContainer& pulses,
00216             PulseMap& pulseMap)
00217 {
00218   debug() << "mapPulsesByChannel(): Processing " << pulses.size() << " rpc pulses." << endreq;
00219   DayaBay::ElecPulseCollection::PulseContainer::const_iterator it,
00220     pulseDone = pulses.end();
00221   int timeSliceNo;
00222   // loop over pulses in collection
00223   for (it=pulses.begin(); it != pulseDone; ++it) {
00224     DayaBay::ElecRpcPulse* pulse = dynamic_cast<DayaBay::ElecRpcPulse*>(*it);
00225     if(pulse == NULL){
00226       // Catch invalid pulses
00227       error() << "mapPulsesByChannel(): bad pulse." << endreq;
00228       return StatusCode::FAILURE;
00229     }
00230     DayaBay::FecChannelId chId(pulse->channelId().fullPackedData());
00231     pulseMap[chId.boardId()][chId].push_back(pulse);
00232   }
00233   return StatusCode::SUCCESS;
00234 }
00235 
00236 StatusCode EsIdealFecTool::generateOneChannel(std::vector<DayaBay::ElecRpcPulse*>& channelPulses,
00237                     DayaBay::ElecFecChannel& channel,
00238                     double clockOffset,
00239                     const ServiceMode&)
00240 {
00241     DayaBay::DigitalSignal& digitSignal = channel.hits();
00242     std::vector<int>& indexes = channel.hitIndexes();
00243     // loop over all pulses in the channel
00244     int size = channelPulses.size();
00245     for (int i = 0; i < size; i++) {
00246         DayaBay::ElecRpcPulse* pulse = channelPulses[i];
00247         // time of the pulse relative to the 0th clock tick (earliest + clockShift)
00248         double offsetT = pulse->time() - clockOffset;
00249         // which clock tick corresponds to the time
00250         int clock = int ( offsetT * m_simFrequency + 1. );
00251         if ( clock >= digitSignal.size() ) {
00252             error()<<"generateOneChannel(): Calculated time of the pulse "<<i<<", "
00253                 <<offsetT/Units::nanosecond
00254                 <<" ns, is out of the simulation time window"<<endreq;
00255             return StatusCode::FAILURE;
00256         }
00257 
00258         if( (pulse->amplitude()) > 0.5 ) {
00259             // pulse crossed the threshold (virtual)
00260             double ran = m_uni();
00261             if(ran < m_idealEff){
00262                 for (int j = 0; j < m_elecStretch; j++) { // stretch the incoming signal
00263                     digitSignal[clock+j] = 1;
00264                     indexes.push_back(clock+j);
00265                 }
00266             }
00267         }
00268     }
00269     return StatusCode::SUCCESS;
00270 }
00271 
00272 StatusCode EsIdealFecTool::addChannelNoise(const DayaBay::FecChannelId& channelId,
00273                     DayaBay::ElecFecChannel& channel,
00274                     double simTime,
00275                     const ServiceMode&)
00276 {
00277     DayaBay::DigitalSignal& digitSignal = channel.hits();
00278     std::vector<int>& indexes = channel.hitIndexes();
00279 
00280     // Calculate mean number of hits in the time window
00281     double nhitsMean = simTime * m_idealNoise;
00282     // Maximum number of noise hits
00283     int nhitsMax = simTime * m_simFrequency;
00284 
00285     int nhits = getNoiseHits(nhitsMean, nhitsMax);
00286     if (nhits < 0) {
00287         error()<<"Cannot get number of noise hits for this channel: "
00288             << channelId << endreq;
00289         return StatusCode::FAILURE;
00290     }
00291 
00292     if ( msgLevel(MSG::VERBOSE) ) {
00293         verbose()<<"Generating "<<nhits<<" noise hits in channel "<<channelId<<endreq;
00294     }
00295 
00296     for (int i = 0; i < nhits; i++) {
00297         double ran = m_uni();
00298         int hit = ran * digitSignal.size();
00299         if (hit < digitSignal.size()) {
00300             for (int j = 0; j < m_elecStretch; j++) {
00301                 digitSignal[hit + j] = 1;
00302                 indexes.push_back(hit+j);
00303             }
00304         }
00305     }
00306 
00307     return StatusCode::SUCCESS;
00308 }
00309 
00310 int EsIdealFecTool::getNoiseHits(double mean, int max)
00311 {
00312     if (mean < 0.) {
00313         fatal() << "Trying to generate integer from poisson distribution with negative parameter: "
00314             << mean << endreq;
00315         return -1;
00316     }
00317     if (mean == 0.)
00318         return 0;
00319 
00320     double ran = m_uni() * exp(mean);
00321     int nhits = 0;
00322     double p = 1.;
00323     double F = 1.;
00324     while (F < ran && nhits < max) {
00325         nhits++;
00326         p *= mean/nhits;
00327         F += p;
00328     }
00329     if (nhits >= max) {
00330         info() << "Number of noise hits was about to exceed maximum allowed. "
00331             << "setting it to the maximum: " << max << endreq;
00332         return max;
00333     }
00334     return nhits;
00335 }
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

Generated on Fri May 16 2014 10:16:37 for ElecSim by doxygen 1.7.4