In This Package:

ElecSimProc.cc

Go to the documentation of this file.

#include "ElecSimProc.h"

#include "HepMC/GenEvent.h"
#include "Event/GenHeader.h"
#include "Event/SimHeader.h"
#include "Event/SimHitCollection.h"
#include "Event/ElecHeader.h"
#include "Event/ElecFeeCrate.h"
#include "Event/ElecFecCrate.h"
#include "DataSvc/IRunDataSvc.h"
#include "Event/RunData.h"
#include "ElecSim/IEsPulseTool.h"
#include "ElecSim/IEsFrontEndTool.h"
#include "Conventions/Electronics.h"
#include "CLHEP/Units/SystemOfUnits.h"
#include "MuonProphet/MpMuonFate.h"

#include <map>
using namespace std;

ElecSimProc::ElecSimProc(const std::string& name, ISvcLocator* pSvcLocator)
    : StageProcessor<DayaBay::ElecHeader>(name,pSvcLocator)
    , m_pmtTool(0)
    , m_rpcTool(0)
    , m_feeTool(0)
    , m_fecTool(0)
{

    m_pullMode=true;

    // ElecSim
    declareProperty("Detectors",m_detectorNames,"List of active detectors");
    declareProperty("PmtTool",m_pmtToolName="EsPmtEffectPulseTool",
                    "Name of the PMT simulation tool");
    declareProperty("RpcTool",m_rpcToolName="EsIdealPulseTool",
                    "Name of the RPC simulation tool");
    declareProperty("FeeTool",m_feeToolName="EsIdealFeeTool",
                    "Name of the PMT Front-end electronics simulation tool");
    declareProperty("FecTool",m_fecToolName="EsIdealFecTool",
                    "Name of the RPC Front-end electronics simulation tool");
    // Default configuration is DayaBay near site, ADs only.
    m_detectorNames.push_back("DayaBayAD1");
    m_detectorNames.push_back("DayaBayAD2");
    m_detectorNames.push_back("DayaBayIWS");
    m_detectorNames.push_back("DayaBayOWS");
    m_detectorNames.push_back("DayaBayRPC");
    m_detectorNames.push_back("LingAoAD1");
    m_detectorNames.push_back("LingAoAD2");
    m_detectorNames.push_back("LingAoIWS");
    m_detectorNames.push_back("LingAoOWS");
    m_detectorNames.push_back("LingAoRPC");
    m_detectorNames.push_back("FarAD1");
    m_detectorNames.push_back("FarAD2");
    m_detectorNames.push_back("FarAD3");
    m_detectorNames.push_back("FarAD4");
    m_detectorNames.push_back("FarIWS");
    m_detectorNames.push_back("FarOWS");
    m_detectorNames.push_back("FarRPC");

    // special for this package
    declareProperty("MinGapTime", m_minTimeGap =
                    DayaBay::preTimeTolerance + DayaBay::postTimeTolerance,
                    "Minimum width for a gap that defines a chunk of data.");
    declareProperty("LongestJump", m_longest = 10*365*24*60*60, 
                    "A protection for longest jump in time stamp. Default is set to 10 year.");

    m_hitPipeline.clear();
    m_giantHits.clear();
    m_giantHitsInTime.clear();
    m_giantHitSheader.clear();

    m_currentTime = TimeStamp::GetBOT();
    m_realElecHeaderTime = TimeStamp::GetBOT();
    m_emptyElecHeaderTime = TimeStamp::GetBOT();
}

ElecSimProc::~ElecSimProc()
{
}

StatusCode ElecSimProc::initialize()
{
    StatusCode sc = this->StageProcessor<DayaBay::ElecHeader>::initialize();
    if (sc.isFailure()) return sc;

    // Convert detector names to Detector IDs
    int nDetNames = m_detectorNames.size();
    for(int detIdx=0; detIdx<nDetNames; detIdx++){
        DayaBay::Detector det(m_detectorNames[detIdx]);
        if(det.site() == Site::kUnknown
           || det.detectorId() == DetectorId::kUnknown){
            error() << "Invalid detector name: " << m_detectorNames[detIdx] << endreq;
            return StatusCode::FAILURE;
        }
        m_detectors.push_back(det);
    }    

    // Get PMT Pulse tool
    try {
        m_pmtTool = tool<IEsPulseTool>(m_pmtToolName);
    }
    catch(const GaudiException& exg) {
        fatal() << "Failed to get pmt tool: \"" << m_pmtToolName << "\"" << endreq;
        return StatusCode::FAILURE;
    }
    info () << "Added tool " << m_pmtToolName << endreq;
  
    // Get RPC Pulse tool
    if(m_pmtToolName == m_rpcToolName){
        m_rpcTool = m_pmtTool;
    }else {
        try {
            m_rpcTool = tool<IEsPulseTool>(m_rpcToolName);
        }
        catch(const GaudiException& exg) {
            fatal() << "Failed to get rpc tool: \"" << m_rpcToolName << "\"" << endreq;
            return StatusCode::FAILURE;
        }
        info () << "Added tool " << m_rpcToolName << endreq;
    }

    // Get PMT Front-end Electronics tool
    try {
        m_feeTool = tool<IEsFrontEndTool>(m_feeToolName);
    }
    catch(const GaudiException& exg) {
        fatal() << "Failed to get fee tool: \"" << m_feeToolName << "\"" << endreq;
        return StatusCode::FAILURE;
    }
    info () << "Added tool " << m_feeToolName << endreq;

    // Get RPC Front-end Electronics tool
    try {
        m_fecTool = tool<IEsFrontEndTool>(m_fecToolName);
    }
    catch(const GaudiException& exg) {
        fatal() << "Failed to get fec tool: \"" << m_fecToolName << "\"" << endreq;
        return StatusCode::FAILURE;
    }
    info () << "Added tool " << m_fecToolName << endreq;
    
    // Set the 'partition' for this run in RunDataSvc
    IRunDataSvc* runDataSvc = svc<IRunDataSvc>("RunDataSvc",true);
    if(!runDataSvc){
        error() << "Failed to get RunDataSvc" << endreq;
        return StatusCode::FAILURE;
    }
    int task = 1;  // Tells RunDataSvc to use new simulation run data.
    Context emptyContext; // Need a dummy context for this service call
    ServiceMode svcMode(emptyContext, task);
    const DayaBay::RunData* runData = runDataSvc->runData(svcMode);
    if(!runData){
        error() << "Failed to get RunData for simulation" << endreq;
        return StatusCode::FAILURE;
    }
    DayaBay::RunData modifiedRunData(*runData);
    modifiedRunData.setDetectors(m_detectors);
    sc = runDataSvc->setRunData(modifiedRunData);
    if( !sc.isSuccess() ){
        error() << "Failed to set detectors in run" << endreq;
    }

    info()<<"minTimeGap "<<m_minTimeGap << " ns" <<endreq;

    return sc;
}

StatusCode ElecSimProc::execute()
{
    debug()<<"Executing ..."<<endreq;
    // Only do anything if the current stage time has advanced beyond
    // the last time we ran.
    static bool first_time = true;
    if (!first_time && m_currentTime > thisStage()->currentTime()) {
        return StatusCode::SUCCESS;
    }
    first_time = false;

    StatusCode sc;

    sc = this->findGap();
    if (sc.isFailure()) return sc;

    sc = this->fillHitHeader();
    if (sc.isFailure()) return sc;    

    sc = this->findGiantHitInTime();
    if (sc.isFailure()) return sc;
    
    // Number of normal hits is above zero
    if( m_count>0 )  {
        sc = this->preExecute();
        if (sc.isFailure()) return sc;
        sc = this->ES_execute();
        if (sc.isFailure()) return sc;
        sc = this->postExecute();
        if (sc.isFailure()) return sc;
    }
    
    // Number of giant hits is above zero
    if( m_giantHitsInTime.size()>0 )  {
        sc = this->fastElecSim();
        if (sc.isFailure()) return sc;
    }

    // update m_currentTime
    if( m_realElecHeaderTime > m_emptyElecHeaderTime )  {
        m_currentTime = m_realElecHeaderTime;
    } else {
        m_currentTime = m_emptyElecHeaderTime;
    }

    // blowChunks() must be the last one, because it will hold the last refcount of some
    // input headers. Before do this, it must be added to inputheaders in postExecute().
    sc = this->blowChunks();
    if (sc.isFailure()) return sc;    

    return StatusCode::SUCCESS;
}

// Get the next SimHeader
StatusCode ElecSimProc::getSimData(DayaBay::SimHeader*& output)
{
    // Currently we get it from the lower stage, but this function
    // will eventually handle getting it from AES when the SimHeaders
    // are read in from file.

    IStageData* stageData = 0;
    StatusCode sc = lowerStage()->nextElement(stageData);
    if (sc.isFailure()) {
        error() << "Failed to get data from detector simulation stage"
                << endreq;
        return sc;
    }
    
    DayaBay::SimHeader* shead = 0;
    
    SimData* input=0;
    input = dynamic_cast<SimData*>(stageData);
    
    if( input!=0 ) {
        shead = &input->header();
    } else {
        HeaderData* pHeaderData = dynamic_cast<HeaderData*>(stageData);
        DayaBay::HeaderObject* pHeaderObj = &pHeaderData->header();
        shead = dynamic_cast<DayaBay::SimHeader*>(pHeaderObj);
    }
    
    if( shead==0 ) {
        error() << "Failed to get SimHeader from lower stage" <<endreq;
    }

    debug() << "stage data's time "<< stageData->time()<<endreq;

    // we must temporarily hold on to this DataObject because this
    // input may have been created during previous exection cycle.
    shead->addRef();            // And because deleting the SimData
    delete stageData;           // drops the only other reference.  It
    stageData=0;                // will be released once it is no
                                // longer providing any hits to the
                                // pipeline

    debug() << "Current time of lower stage "<<lowerStage()->currentTime()<<endreq;
    debug() << "Got SimHeader from lower stage at " << shead->timeStamp() << endreq;

    output = shead;
    return StatusCode::SUCCESS;
}

// Fill hit pipeline and our SimHitHeader with given SimHeader's worth
// of hits
StatusCode ElecSimProc::fillPipeLine(DayaBay::SimHeader* shead)
{
    // pipeline
    TimeStamp realLatest = TimeStamp::GetBOT();

    // Get the SimHit Header
    const DayaBay::SimHitHeader* shitHeader = shead->hits();
    if (!shitHeader) {
        warning() << "Got NULL SimHitHeader with SimHeader @ "
                  << shead->timeStamp()
                  << " skipping..."
                  << endreq;
        return StatusCode::SUCCESS;
    }
    const DayaBay::SimHitHeader::hc_map& hcMap = shitHeader->hitCollection();

    const TimeStamp& absTime = shead->timeStamp();
    debug()<<"absTime = "<<absTime<<endreq;
    DayaBay::SimHitHeader::hc_map::const_iterator hcIter;
    for (hcIter=hcMap.begin(); hcIter != hcMap.end(); ++hcIter) {

        DayaBay::Detector det(hcIter->first);
        DayaBay::SimHitCollection* hits = hcIter->second;
        if(!hits) return StatusCode::FAILURE;
        
        const DayaBay::SimHitCollection::hit_container& hc = hits->collection();
        for (size_t ind=0; ind < hc.size(); ++ind) {
            DayaBay::SimHit* hit = hc[ind];

            DayaBay::AdPmtSensor pmtid(hit->sensDetId());
            if (pmtid.bogus()) {
                warning() << "SimHeaderCnvHits: warning got bogus PMTid: " << pmtid << endreq;
            }

            TimeStamp ts = absTime;
            if( hit->hitTime()/CLHEP::second < m_longest ) {
                ts.Add(hit->hitTime()/CLHEP::second);
                m_hitPipeline.insert(pair<TimeStamp,DayaBay::SimHit*>(ts,hit));
                //debug()<<hit->hitTime()/CLHEP::second<<" [sec]"<<endreq;
                // find the real latest time
                if( ts>realLatest ) realLatest = ts;
            } else {
                warning()<<"A hit in the far future "<<hit->hitTime()/CLHEP::second<<" [sec] is skipped."<<endreq;
            }
        }
    }

    const DayaBay::IHeader* iheader = (shead->inputHeaders()) [0]; 
    const DayaBay::GenHeader* genHeader = dynamic_cast<const DayaBay::GenHeader*>(iheader);
    const HepMC::GenEvent* genEvent = genHeader->event();

    if( genEvent )  {
        HepMC::GenEvent::particle_const_iterator pci, pci_end = genEvent->particles_end();
    
        for( pci = genEvent->particles_begin(); pci != pci_end; ++pci )  {
            if( (*pci)->pdg_id() == 13 || (*pci)->pdg_id() == -13 )  {   // muon
                if( (*pci)->status() != MpMuonFate::kNeedSim )   {   //  geant4-noble particle
                    debug()<<"Find one geant4-noble particle"<<endreq;
                    DayaBay::SimHit* aGiantHit = new DayaBay::SimHit;
                    double dt = 0 * CLHEP::nanosecond;
                    aGiantHit->setHitTime(dt);
                    aGiantHit->setWeight(0);
                    aGiantHit->setSensDetId( DetectorId::kUnknown );
                    debug()<<"absTime "<<absTime<<endreq;
                    TimeStamp ts = absTime;
                    ts.Add(dt/CLHEP::second);
                    debug()<<"ts      "<<ts<<endreq;
                    m_hitPipeline.insert(pair<TimeStamp,DayaBay::SimHit*>(ts,aGiantHit));
                    m_giantHits.insert(pair<TimeStamp,DayaBay::SimHit*>(ts,aGiantHit));
                    m_giantHitSheader[aGiantHit] = shead;
                    debug()<<"aGiantHit is added "<<aGiantHit<<" with SimHeader "<<shead<<" at time "<<ts<<endreq;
                    debug()<<"Now m_giantHitSheader has "<<m_giantHitSheader.size()<<" giant hits."<<endreq;
                    // find the real latest time
                    if( ts>realLatest ) realLatest = ts;
                    break; // only process one fast simulated muon
                }
            }
        }
    }
    
    debug()<<"Current m_giantHitSheader begin"<<endreq;
    map<const DayaBay::SimHit*, const DayaBay::SimHeader*>::iterator it,itend = m_giantHitSheader.end();
    for( it=m_giantHitSheader.begin(); it!=itend; it++ ) {
        debug()<<"aGiantHit "<<it->first<<" with SimHeader "<<it->second<<endreq;
    }

    // current SimHeaders
    m_currentHeaders.push_back(shead);
    debug()<<"m_currentHeaders has "<<m_currentHeaders.size()<<" SimHeaders."<<endreq;

    // record the real latest time
    m_realLatest.insert( pair<DayaBay::SimHeader*, TimeStamp>(shead,realLatest) );
    
    return StatusCode::SUCCESS;
}

StatusCode ElecSimProc::moar()
{
    debug() << "Getting MOAR!" << endreq;

    DayaBay::SimHeader* shead = 0;
    StatusCode sc = this->getSimData(shead);
    if (sc.isFailure()) return sc;
    
    sc = this->fillPipeLine(shead);
    if (sc.isFailure()) return sc;

    return StatusCode::SUCCESS;
}

StatusCode ElecSimProc::findGap()
{
    debug()<<"Entering findGap()"<<endreq;
    while ( m_hitPipeline.size() < 2 )  {  // at least two hits are needed to find a gap
        StatusCode sc = this->moar(); // need MOAR!
        if (sc.isFailure()) return sc;
    }

    HitPipeline::iterator current;
    HitPipeline::iterator previous;

    TimeStamp GapStop;

    bool GapFound=false;
    bool TimeTest=false;

    //
    // _________hhhhhh__hh_h____hhhhh_H ________________hhhh_H_hhh_Hhhhhhhhhhh___
    //          Chunk Start        Stop|GapStart    Stop
    //

    TimeStamp BeginTime;
    TimeStamp EndTime;
    // Gap must be found and there is no further hits from lower stage will get into this gap.
    do {
        debug()<<"Currently hit pipeline have SimHits #: "<<m_hitPipeline.size()<<endreq;
        BeginTime=m_hitPipeline.begin()->first;
        EndTime= m_hitPipeline.rbegin()->first;
        debug()<<"Pipeline begin time: "<<BeginTime<<endreq;
        debug()<<"Pipeline end time:   "<<EndTime<<endreq;

        current = m_hitPipeline.begin();
        // reset them to the beginning of hits pipeline
        m_chunkStart = current->first;
        m_chunkStop  = current->first;

        ++current; // Start from the second one
        for (;current != m_hitPipeline.end(); 
             ++current) {  // try all SimHits in one SimHeader
            previous = current;
            --previous;    // Always compare current one with previous one
            double dt_sec = (current->first - previous->first);  // in second
            double dt_units = dt_sec * CLHEP::second;            // with unit

            verbose() << "Previous hit at " << previous->first 
                      << " has deltaT " << dt_units/CLHEP::nanosecond << " ns" << endreq;
            verbose() << "Current hit at " << current->first << endreq;

            if (dt_units >= m_minTimeGap) {
                m_chunkStop =previous->first;             // the previous one's time, moving it forward
                GapStop=m_chunkStop;
                GapStop.Add(m_minTimeGap/CLHEP::second);  // where the gap stop

                GapFound=true;
                debug() << "Found gap of " << dt_units/CLHEP::nanosecond<< " ns" << endreq;
                if (lowerStage()->currentTime() >= GapStop) {
                    TimeTest=true;
                    debug() << "TimeTest true" <<endreq;
                    debug() << "chunk start " << m_chunkStart <<endreq;
                    debug() << "chunk stop  " << m_chunkStop <<endreq;
                    return StatusCode::SUCCESS;
                } else {
                    TimeTest=false;
                }               
                break;  // Leave the innermost for loop
            }
        } // All SimHits in pipeline
        // Try some hits even gap is found. By asking more hits from lower stage the time can be updated.
        if ( !GapFound || !TimeTest ) {
            StatusCode sc = this->moar(); // need MOAR!  
            if (sc.isFailure()) return sc;
            debug()<<"Right now hit pipeline have SimHits #: "<<m_hitPipeline.size()<<endreq;
            GapFound=false;
            TimeTest=false;
        }
    } while( !GapFound || !TimeTest ); //Both need to be satisfied.
     
    return StatusCode::FAILURE;  // It will never get here
}

StatusCode ElecSimProc::fillHitHeader()
{
    // Take SimHitHeaders from m_currentHeaders and fill hit
    // collections of fake amalSimHitHeader with all hits that have times
    // before chunk stop.
    // A lot of bookkeeping
    debug()<<"Fill an amalgamated SimHitHeader"<<endreq;
    m_count=0;  // count for how many hits add to the fake SimHitHeader
    m_simHeaderCount.clear();

    m_realHitEarliest = TimeStamp::GetEOT();
    m_realHitLatest   = TimeStamp::GetBOT();

    // Clear out any previously filled hit collections
    DayaBay::SimHitHeader::hc_map& amal_hcmap = m_amalSimHitHeader.hitCollection(); // detectors collection
    DayaBay::SimHitHeader::hc_map::const_iterator it, done = amal_hcmap.end();
    for (it = amal_hcmap.begin(); it != done; ++it) {
        DayaBay::SimHitCollection* pHitCollection = it->second;
        pHitCollection->collection().clear();
        delete pHitCollection;
    }
    amal_hcmap.clear();

    // Loop over each previously saved header that contributed at
    // least one hit that came within trunk window to the pipeline.
    for (size_t ihead = 0; ihead < m_currentHeaders.size(); ++ihead) {
        m_simHeaderCount[m_currentHeaders[ihead]]=0;
        
        TimeStamp referenceTime = m_currentHeaders[ihead]->timeStamp();
        const DayaBay::SimHitHeader* shhead = m_currentHeaders[ihead]->hits();
        const DayaBay::SimHitHeader::hc_map& hcMap = shhead->hitCollection();
        DayaBay::SimHitHeader::hc_map::const_iterator hcit, hcdone = hcMap.end();

        // loop over each detectors hit collection
        for (hcit = hcMap.begin(); hcit != hcdone; ++hcit) {

            // Get source hit collection for this detector
            DayaBay::SimHitCollection* shcol = hcit->second;
            if (!shcol) { return StatusCode::FAILURE; }
            const DayaBay::SimHitCollection::hit_container& hc = shcol->collection();
            
            // copy requisite hits
            for (size_t ihit=0; ihit < hc.size(); ++ihit) {
                DayaBay::SimHit* hit = hc[ihit];
                TimeStamp ts = referenceTime;
                if( hit->hitTime()/CLHEP::second > m_longest ) continue;
                ts.Add(hit->hitTime()/CLHEP::second);
                
                //debug()<<"reference time "<< referenceTime<<endreq;
                //debug()<<"delta t "<<hit->hitTime()/CLHEP::second<<endreq;

                if (ts >= m_chunkStart && ts <= m_chunkStop) {
                    m_count++;
                    m_simHeaderCount[m_currentHeaders[ihead]]++;

                    // 
                    if( ts<m_realHitEarliest ) m_realHitEarliest = ts;
                    if( ts>m_realHitLatest )   m_realHitLatest = ts;

                    // Get/make our hit collection for this detector
                    DayaBay::SimHitHeader::hc_map::iterator it_det=amal_hcmap.find(hcit->first);
                    DayaBay::SimHitCollection* my_shcol;
                    if (it_det == amal_hcmap.end()) {
                        my_shcol = new DayaBay::SimHitCollection();
                        //....... what to set its SimHitHeader to? ......... 
                        my_shcol->setDetector(DayaBay::Detector(hcit->first));
                        amal_hcmap[hcit->first] = my_shcol;
                    } else {
                        my_shcol=it_det->second;
                    }
                    
                    if (!ihead) {
                        // Set SimHitHeader to be the one in the first SimHeader.
                        // Is this the right thing to do?
                        my_shcol->setHeader(const_cast<DayaBay::SimHitHeader*>(shhead));
                    }
                    
                    DayaBay::SimHitCollection::hit_container& my_hc = my_shcol->collection();
                    verbose()<<"hit with time "<<ts<<" added to amalgamated SimHitHeader"<<endreq;
                    my_hc.push_back(hit);
                }

            } // hits
        } // detectors
    } // headers

    debug()<<m_count<<" hits added to the amalgamated SimHitHeader"<<endreq;
    debug()<<"Size of m_amalSimHitHeader.hitCollection "<<m_amalSimHitHeader.hitCollection().size()<<endreq;
    debug()<<"Size of m_simHeaderCount "<<m_simHeaderCount.size()<<endreq;

    return StatusCode::SUCCESS;
}


// Get all geant4-noble particles in time window
StatusCode ElecSimProc::findGiantHitInTime()
{
    m_giantHitsInTime.clear();
    HitPipeline::iterator ci, ci_end = m_giantHits.end();
    
    for( ci=m_giantHits.begin(); ci!=ci_end; ++ci ) {
        
        TimeStamp ts = ci->first;
        DayaBay::SimHit* aHit = ci->second;
        
        if (ts >= m_chunkStart && ts <= m_chunkStop) {
            m_giantHitsInTime.insert( pair<TimeStamp,DayaBay::SimHit*>(ts,aHit) );
        }
    }
    return StatusCode::SUCCESS;
}


// It is time to run electronic simulation
StatusCode ElecSimProc::ES_execute()  
{
    DayaBay::ElecHeader* elecHeader = MakeHeaderObject();

    // Set input SimHeaders who has really contribute to the elecHeader
    // set context as the first simheader
    bool found;
    found=false;
    for (size_t ind=0; ind<m_currentHeaders.size(); ++ind) {
        if(m_simHeaderCount[m_currentHeaders[ind]]>0) {  // at least 1 hit is in use.
            this->AppendInputHeader(m_currentHeaders[ind]);
            if(!found) {
                elecHeader->setContext(m_currentHeaders[ind]->context());
                found=true;
            }
        }
    }

    // TimeStamps
    debug()<<"Before setContext "<<elecHeader->timeStamp()<<endreq;

    m_realHitEarliest.Add(-DayaBay::preTimeTolerance/CLHEP::second);  // earliest possible time
    m_realHitLatest.Add(DayaBay::postTimeTolerance/CLHEP::second);    // latest possible time

    elecHeader->setEarliest(m_realHitEarliest);
    elecHeader->setLatest(m_realHitLatest);
    elecHeader->setTimeStamp(m_realHitEarliest);

    m_realElecHeaderTime=m_realHitEarliest;               // They should be always be the same.
    debug()<<"After setContext  "<<elecHeader->timeStamp()<<endreq;

    // Add the pulse collection header
    DayaBay::ElecPulseHeader* pulseHeader = 
        new DayaBay::ElecPulseHeader(elecHeader);
    elecHeader->setPulseHeader(pulseHeader);

    // Add the crate header
    DayaBay::ElecCrateHeader* crateHeader = 
        new DayaBay::ElecCrateHeader(elecHeader);
    elecHeader->setCrateHeader(crateHeader);

    // Send each detector's hit collections through ElecSim tool chains
    // It is from m_amalSimHitHeader.
    const DayaBay::SimHitHeader::hc_map& hcMap = m_amalSimHitHeader.hitCollection();

    debug() << "Processing hit collections" <<endreq;

    DayaBay::SimHitHeader::hc_map::const_iterator hcIter;
    for (hcIter=hcMap.begin(); hcIter != hcMap.end(); ++hcIter) {

        DayaBay::Detector det(hcIter->first);
        debug() << "Checking " << det.detName() << " for hits." << endreq;
        DayaBay::SimHitCollection* hits = hcIter->second;
        if(!hits) return StatusCode::FAILURE;

        debug() << "Got hit collection from " << det.detName()<<" (id= "
                << det.siteDetPackedData() << ") " << " with "
                << hits->collection().size() << " hits." << endreq;

        StatusCode sc;
        // Check if this detector should be simulated.
        if( std::find(m_detectors.begin(),m_detectors.end(),det)
            == m_detectors.end()){
            debug() << "Detector " << det.detName() << " not simulated." << endreq;
            continue;
        }
        if( det.isAD() || det.isWaterShield() ){
            // Process PMT hits
            debug() << "Processing PMT hits" <<endreq;
            DayaBay::ElecPulseCollection* pulses =
                new DayaBay::ElecPulseCollection(pulseHeader, det);
            sc = m_pmtTool->generatePulses(hits, pulses);
            if( sc != StatusCode::SUCCESS ) return sc;
            pulseHeader->addPulseCollection(pulses);

            DayaBay::ElecFeeCrate* crate = new DayaBay::ElecFeeCrate(det,
                                                                     crateHeader);
            sc = m_feeTool->generateSignals(pulses, crate);
            if( sc != StatusCode::SUCCESS ) return sc;
            crateHeader->addCrate(crate);

        }else if(det.detectorId() == DetectorId::kRPC){
            // Process RPC hits
            debug() << "Processing RPC hits" <<endreq;
            DayaBay::ElecPulseCollection* pulses =
                new DayaBay::ElecPulseCollection(pulseHeader, det);
            sc = m_rpcTool->generatePulses(hits, pulses);
            if( sc != StatusCode::SUCCESS ) return sc;
            pulseHeader->addPulseCollection(pulses);

            DayaBay::ElecFecCrate* crate = new DayaBay::ElecFecCrate(det,
                                                                     crateHeader);
            sc = m_fecTool->generateSignals(pulses, crate);
            if( sc != StatusCode::SUCCESS ) return sc;
            crateHeader->addCrate(crate);
        }else{
            error() << "Unknown detector " << det << endreq;
            return StatusCode::FAILURE;
        }
    }

    // Send result to stage.
    ElecData* ed = new ElecData(*elecHeader);
    thisStage()->pushElement(ed);
    this->registerData(*ed);
    debug()<<"to grep: (Full Elec) new data pushed out at time "<<m_realHitEarliest<<endreq;
    
    return StatusCode::SUCCESS;
}

StatusCode ElecSimProc::blowChunks()
{
    // clean pipeline, remove all the hits before m_chunkStop
    HitPipeline::iterator start = m_hitPipeline.begin();
    HitPipeline::iterator end   = m_hitPipeline.end();
    HitPipeline::iterator stop  = m_hitPipeline.end();  // chunkStop
    HitPipeline::iterator it;
    for( it=start; it!=end; ++it )  {
        if(it->first>m_chunkStop) {
            stop=it;
            break;  // Leave the innermost for loop
        }
    }
    //m_hitPipeline.find(m_chunkStop); // only find the first element in line, not good
    m_hitPipeline.erase(start,stop);    // [start, stop) are erased and stop is not included.

    start = m_giantHits.begin();
    end   = m_giantHits.end();
    stop  = m_giantHits.end(); // chunkStop

    for( it=start; it!=end; ++it )  {
        
        DayaBay::SimHit* aHit = it->second;
        TimeStamp ts = it->first;
        
        if( ts>m_chunkStop )  {
            stop = it;
            break;
        }
        
        m_giantHitSheader.erase( aHit );
        delete aHit;
    }
    m_giantHits.erase(start,stop);  // [start,stop)


    // clean used SimHeader if it has no hit with time less than m_chunkStop
    std::vector<const DayaBay::SimHeader*>::iterator ihd;
    for(ihd=m_currentHeaders.begin();ihd!=m_currentHeaders.end();++ihd) {

        // get the real latest
        TimeStamp latest=TimeStamp::GetEOT();
        map<const DayaBay::SimHeader*, TimeStamp>::iterator it = m_realLatest.find(*ihd);
        if( it!=m_realLatest.end() ) {
            latest = it->second;
        } else {
            error()<<"Can't find the latest time of a SimHeader object"<<endreq;
        }
        // use the real latest to remove header objects
        if( latest<=m_chunkStop ) {
            // release the reference to the SimHeader
            DayaBay::SimHeader* nonconst = const_cast<DayaBay::SimHeader*>(*ihd);
            nonconst->release();
            // remove its local copy
            m_currentHeaders.erase(ihd);
            ihd--;
            m_realLatest.erase(it);
        }
    }

    {
        debug()<<"Current m_giantHitSheader end"<<endreq;
        map<const DayaBay::SimHit*, const DayaBay::SimHeader*>::iterator it,itend = m_giantHitSheader.end();
        for( it=m_giantHitSheader.begin(); it!=itend; it++ ) {
            debug()<<"aGiantHit "<<it->first<<" with SimHeader "<<it->second<<endreq;
        }
    }

    return StatusCode::SUCCESS;
}


// Fast electronic simulation, produce empty ElecHeader for MuonProphet muon
StatusCode ElecSimProc::fastElecSim()
{
    StatusCode sc;

    TimeStamp TimeOfLast = TimeStamp::GetBOT();
    HitPipeline::iterator it, it_end = m_giantHitsInTime.end();

    for( it=m_giantHitsInTime.begin(); it!=it_end; ++it )  {
        
        sc = this->preExecute();
        if (sc.isFailure()) return sc;

        DayaBay::ElecHeader* elecHeader = MakeHeaderObject();
        DayaBay::SimHit* aHit = it->second;
        TimeStamp ts = it->first;

        TimeStamp earliest = ts;
        TimeStamp latest =   ts;
        earliest.Add(-DayaBay::preTimeTolerance/CLHEP::second);  // earliest possible time
        earliest.Add(-1e-9); // 1ns earlier than its real simulated products
                             // Then it can always be popped out and simulated first
        latest.Add(DayaBay::postTimeTolerance/CLHEP::second);
        
        //const DayaBay::SimHeader* sheader = m_giantHitSheader[aHit];
        map<const DayaBay::SimHit*, const DayaBay::SimHeader*>::iterator itmap = m_giantHitSheader.find(aHit);
        const DayaBay::SimHeader* sheader;
        if( itmap != m_giantHitSheader.end() ) {
            sheader = itmap->second;
        } else {
            error()<<"Can't find aGiantHit's parent"<<endreq;
        }

        debug()<<"earliest = "<<earliest<<endreq;
        debug()<<"aHit at "<<aHit<<"   SimHeader at "<<sheader<<endreq;

        elecHeader->setContext(sheader->context());
        elecHeader->setEarliest(earliest);
        elecHeader->setLatest(latest);
        elecHeader->setTimeStamp(earliest);
        
        if( earliest > TimeOfLast ) {
            TimeOfLast = earliest;
        }

        this->AppendInputHeader( sheader );

        // Send result to stage.
        thisStage()->pushElement(new ElecData(*elecHeader));
        debug()<<"to grep: (Fast Elec) new data pushed out at time "<<earliest<<endreq;
        debug()<<"elecHeader at "<<elecHeader<<endreq;
        // Add ElecPulseHeader and ElecCrateHeader
        DayaBay::ElecPulseHeader* pPulse = 0; // new DayaBay::ElecPulseHeader;
        DayaBay::ElecCrateHeader* pCrate = 0; // new DayaBay::ElecCrateHeader;
        elecHeader->setPulseHeader(pPulse);
        elecHeader->setCrateHeader(pCrate);

        sc = this->postExecute();
        if (sc.isFailure()) return sc;
    }

    if( m_emptyElecHeaderTime < TimeOfLast ) {
        m_emptyElecHeaderTime = TimeOfLast;
    } else {
        error()<<"New simulated data has earlier time. Why?"<<endreq;
    }
    
    return StatusCode::SUCCESS;
}

// Local Variables:
// c-basic-offset: 4
// End:

---