In This Package:

DigitizeAlg.cc

Go to the documentation of this file.

#include "DigitizeAlg/DigitizeAlg.h"

#include "Event/ElecHeader.h"
#include "Event/ElecFeeCrate.h"
#include "Event/ElecFecCrate.h"
#include "Event/RegistrationSequence.h"

#include "ElecSim/IEsPulseTool.h"
#include "ElecSim/IEsFrontEndTool.h"
#include "Conventions/Electronics.h"

#include "SimHitSplitSvc/SimHitSplitSvc.h"
#include "Event/SimReadoutHeader.h"

#include <map>
#include <algorithm>

using namespace DayaBay;

DigitizeAlg::DigitizeAlg(const string& name, ISvcLocator* pSvcLocator)
            :DybAlgorithm<SimReadoutHeader>(name,pSvcLocator)
    , m_pmtTool(0)
    , m_rpcTool(0)
    , m_feeTool(0)
    , m_fecTool(0)
    , m_triggerSorter(0)
    , m_roTrigPackTool(0)
    , m_simHitSplitSvc(0)
{
  //Earliest and Latest
  //TimeStamp initime(0.0);
  //m_hitEarliest = initime;
  //m_hitLatest   = initime;
  //Electronics
  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 of the Detectors
  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");

  //Trigger
  m_trigToolNames.push_back("TsMultTriggerTool");
  declareProperty("TrigTools",     m_trigToolNames,   "Tools to generate Triggers");
  declareProperty("TriggerFilter", m_trigFilterName = "TsSortTool",
                  "Tool to do final filtering/combining of triggers");
    
  //Readout
  m_roToolNames.push_back("ROsFeeReadoutTool");
  m_roToolNames.push_back("ROsFecReadoutTool");
  declareProperty("RoTools",          m_roToolNames,     "Tools to generate Triggers");
  declareProperty("TrigPackagerTool", m_roTrigPackName = "ROsTriggerDataPackerTool",
                  "Tool to convert trigger commands from TrigSim to TriggerDataPackage");
}

DigitizeAlg::~DigitizeAlg()
{
}

StatusCode DigitizeAlg::initialize()
{ 
  //ElecSim initialization
  //Convert detector names to Detector IDs
  for(vector<string>::size_type detIdx = 0; detIdx < m_detectorNames.size(); detIdx++)
  {
    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);
  }

  StatusCode sc = this->GaudiAlgorithm::initialize();
  if(sc.isFailure()) 
  { 
    error()<< "Error when initializing GaudiAlgorithm."<< endreq;
    return sc;
  }
   
  debug()<< "DigitizeAlg ininalize()"<< endreq;
 
  //Initialize the SimHitSplitSvc
  IService* isvc = 0;
  sc = serviceLocator()->service("SimHitSplitSvc", isvc, false);
  if(sc.isFailure()) 
  {
    error()<< "Service : SimHitSplitSvc could not be retrieved!"<< endreq;
    return sc;
  }
  isvc->addRef();
  sc = isvc->queryInterface(ISimHitSplitSvc::interfaceID(), (void**)&m_simHitSplitSvc);
  if(sc.isFailure())
  {
    error()<< "SimHitSplitSvc does not implement ISimHitSplitSvc"<< endreq;
    return sc;
  }

  //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;
  }
  debug()<< "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;
    }
    debug()<< "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;
  }
  debug()<< "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;
  }
  debug()<< "Added tool "<< m_fecToolName<< endreq;

  //Trigger
  for(vector<string>::size_type ind = 0; ind < m_trigToolNames.size(); ind++)
  {
    string tgr = m_trigToolNames[ind];
    try
    {
      m_trigTools.push_back(tool<ITsTriggerTool>(tgr));
    }
    catch(const GaudiException& exg) 
    {
      fatal()<< "Failed to get Trigger Tool: \""<< tgr<< "\""<< endreq;
      return StatusCode::FAILURE;
    }
    debug()<< "Added trigger tool "<< tgr<< endreq;
  }
  try
  {
    m_triggerSorter = tool<ITsTriggerTool>(m_trigFilterName);
  }
  catch(const GaudiException& exg)
  {
    fatal()<< "Failed to get TriggerFilterTool: \""<< m_trigFilterName<< "\""<< endreq;
    return StatusCode::FAILURE;
  }

  //Readout
  try
  {
    m_roTrigPackTool = tool<IROsTriggerDataPackerTool>(m_roTrigPackName);
  }
  catch(const GaudiException& exg)
  {
    fatal()<< "Failed to get Readout Trigger Package Tool: \"" 
           << m_roTrigPackName<< "\""<< endreq;
    return StatusCode::FAILURE;
  }
    
  for(vector<string>::size_type ind = 0; ind < m_roToolNames.size(); ind++)
  {
    string rotool = m_roToolNames[ind];
    try
    {
      m_roTools.push_back(tool<IROsReadoutTool>(rotool));
    }
    catch(const GaudiException& exg)
    {
      fatal()<< "Failed to get Readout Tool: \""<< rotool<< "\""<< endreq;
      return StatusCode::FAILURE;
    }
    debug()<< "Added readout tool "<< rotool<< endreq;
  }

  return sc;
}

StatusCode DigitizeAlg::execute()
{
  debug()<< "Digitalization execute()...."<< endreq;

  // Using some functions in RegistrationSequence
  //RegistrationSequence* seq = getTES<RegistrationSequence>(RegistrationSequenceLocation::Default);
  //debug()<< "ATTENTION! RegistrationSequence Size = "<< seq->size()<< endreq;

  //Not storing SimHeader
  /*IRegistrationSequence::Registrations& regs = const_cast<IRegistrationSequence::Registrations&>(seq->registrations());
  IRegistrationSequence::Registrations::iterator regit;
  for(regit = regs.begin(); regit != regs.end(); regit++)
  {
    debug()<< regit->path()<< endreq;
    if(regit->path() == "/Event/Sim/SimHeader")
    {
      //DataObject* child = const_cast<DataObject*>(regit->object());   
      //SimHeader*  sh    = dynamic_cast<SimHeader*>(child);
      regit->setStore(false);
      debug()<< "Not Storing: "<< regit->path()<< endreq;
    }
    }*/

  //Get one SimHitHeader and a vector of SimHeaders
  debug()<< "Now get SimHitHeaders!"<< endreq;

  DayaBay::SimHeader* sh = getTES<SimHeader>(SimHeaderLocation::Default);
  vector<SimHitHeader*> shhs;

  //From SimHeader to SimReadoutHeader
  SimReadoutHeader* simReadoutHeader = MakeHeaderObject();
  simReadoutHeader->setContext(sh->context());
  SimReadoutHeader::SimReadoutContainer &readouts = simReadoutHeader->readouts();

  //Empty SimHeader
  TimeStamp earliest(sh->earliest());
  TimeStamp latest(sh->latest());
  TimeStamp timeStamp(sh->timeStamp());

  StatusCode sc = m_simHitSplitSvc->getSimHitHeader(sh, shhs, m_hitEarliest, m_hitLatest);
  if(sc.isFailure()) return sc;

  debug()<< "Number of SimHits segments: "<< shhs.size()<< endreq;
  //None-empty SimHeader
  for(vector<SimHitHeader*>::size_type i = 0; i < shhs.size(); i++)
  {
    //Get one ElecHeader from ElecSim
    debug()<< "Now get one ElecHeader!"<< endreq;
    ElecHeader* eh = new ElecHeader();
    debug()<< "ElecHeader at "<< eh<< "\n"
          << "Electronics simulation...."<< endreq;
    StatusCode sc = this->ElecSim(shhs[i], i, eh);
    if(sc.isFailure()) 
    {
      error()<< "Error in ElecSim():"<< endreq;
      return sc;
    }
    //delete shhs[i];

    //Define output SimReadoutHeader
    //Do trigger and readout simulation
    debug()<< "Now get one SimReadout....."<< endreq;
    SimReadoutHeader* srh = new SimReadoutHeader();
    debug()<< "SimReadoutHeader at "<< srh<< "\n"
          << "SimReadout simulation...."<< endreq;
    sc = this->TrigExecute(eh, srh);
    if(sc.isFailure()) 
    {
      error()<< "Error in TrigExecute():"<< endreq;
      return sc;
    }
    delete eh;

    if(0 == i)
    {
      earliest  = srh->earliest();
      timeStamp = srh->timeStamp();
    }
    if(i == (shhs.size() - 1)) latest = srh->latest();

    SimReadoutHeader::SimReadoutContainer &src = srh->readouts();
    for(unsigned j = 0; j < src.size(); j++)
    {
      src[j]->setHeader(simReadoutHeader);
      readouts.push_back(src[j]);
    }
    
    srh->readouts().clear();
    delete srh;
  }

  if(shhs.size() != 1 || shhs[0] != sh->hits())
  {
    for(vector<SimHitHeader*>::size_type i = 0; i < shhs.size(); i++)
    {
      SimHitHeader::hc_map& maphc = shhs[i]->hitCollection();
      for(SimHitHeader::hc_map::iterator it = maphc.begin(); it != maphc.end(); ++it)
      {
        it->second->collection().clear();
      }
      delete shhs[i];
    }
  }

  simReadoutHeader->setEarliest(earliest);
  simReadoutHeader->setLatest(latest);
  simReadoutHeader->setTimeStamp(timeStamp);

  m_hitEarliest.clear();
  m_hitLatest.clear();

  return StatusCode::SUCCESS;
}

//ElecSim
StatusCode DigitizeAlg::ElecSim(const SimHitHeader* shh, int index, ElecHeader*& eh)
{
  //setTimestamp must be after setContext.
  debug()<< "Before setContext(): "<< eh->timeStamp().GetSec()
         << "--"<< eh->timeStamp().GetNanoSec()<< endreq;

  m_hitEarliest[index].Add(-preTimeTolerance/CLHEP::second);
  m_hitLatest[index].Add(postTimeTolerance/CLHEP::second);

  eh->setEarliest(m_hitEarliest[index]);
  eh->setLatest(m_hitLatest[index]);
  eh->setTimeStamp(m_hitEarliest[index]);

  debug()<< "After setContext(): TimeStamp: "<< eh->timeStamp().GetSec()
        << "--"<< eh->timeStamp().GetNanoSec()<< endreq;
  debug()<< "After setContext(): Earliest: "<< eh->earliest().GetSec()
        << "--"<< eh->earliest().GetNanoSec()<< endreq;
  debug()<< "After setContext(): Latest: "<< eh->latest().GetSec()
        << "--"<< eh->latest().GetNanoSec()<< endreq;

  //Add the pulse header
  ElecPulseHeader* eph = new ElecPulseHeader(eh);
  eh->setPulseHeader(eph);

  //Add the crate header
  ElecCrateHeader* ech = new ElecCrateHeader(eh);
  eh->setCrateHeader(ech);

  //Send each detector's hit collections through ElecSim tool chains
  //From the new SimHitHeader.
  map<short int, SimHitCollection*> mshc = shh->hitCollection();

  debug()<< "Processing SimHit collections: "<< mshc.size()<< endreq;

  //Loop over SimHitCollections in all detectors
  map<short int, SimHitCollection*>::iterator mit;
  for(mit = mshc.begin(); mit != mshc.end(); mit++)
  {
    Detector det(mit->first);
    debug()<< "Checking "<< det.detName()<< " for hits."<< endreq;
    SimHitCollection* hits = mit->second;
    if(!hits) return StatusCode::FAILURE;

    if(det.isAD())
    {
      const vector<SimHit*> simhits = hits->collection();
      debug()<< "There are "<< simhits.size()<< " ADPMTSimHits"<< endreq;
      /*for(vector<SimHit*>::const_iterator vit = simhits.begin(); vit != simhits.end(); vit++)
      {
        AdPmtSensor sensDetId((*vit)->sensDetId());
        debug()<< "SimHit properties for PMT: ring: "<< sensDetId.ring()<< " column: "<< sensDetId.column() 
              << " "<< sensDetId.sensorId()<< endreq;
              }*/
    }
    else if(det.isWaterShield())
    {
      const vector<SimHit*> simhits = hits->collection();
      debug()<< "There are "<< simhits.size()<< " WSPMTSimHits"<< endreq;      
      /*for(vector<SimHit*>::const_iterator vit = simhits.begin(); vit != simhits.end(); vit++)
      {
        PoolPmtSensor sensDetId((*vit)->sensDetId());
        debug()<< "SimHit properties for PMT: wallNum: "<< sensDetId.wallNumber()
              << " wallSpot: "<< sensDetId.wallSpot() << " "<< sensDetId.sensorId()<< endreq;
              }*/
    }

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

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

      ElecFeeCrate* crate = new ElecFeeCrate(det, ech);
      debug()<< "Before generateSignals: Crate: "<< crate<< endreq;
      sc = m_feeTool->generateSignals(pulses, crate);
      
      if(sc != StatusCode::SUCCESS) return sc;

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

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

  debug()<< "To grep: (Full Elec) new data pushed out at time "<< m_hitEarliest<< endreq;
    
  return StatusCode::SUCCESS;
}

StatusCode DigitizeAlg::TrigSim(const ElecHeader& eh, SimTrigHeader*& th)
{
  th->setCommandHeader(new SimTrigCommandHeader(th));
        
  //Set triggers to have same validity range as the crate it belongs to.
  th->setContext(eh.context());
  th->setTimeStamp(eh.timeStamp());
  th->setEarliest(eh.earliest());
  th->setLatest(eh.latest());

  debug()<< "TrigSim time: "<< endreq;
  debug()<< "Earliest: "<< th->earliest().GetSec()<< "--"<< th->earliest().GetNanoSec()<< " "
        << "Latest: "<< th->latest().GetSec()<< "--"<< th->latest().GetNanoSec()<< endreq;
    
  //Let each tool do its thing on the event
  for(vector<ITsTriggerTool*>::size_type ind = 0; ind < m_trigTools.size(); ind++)
  {
    debug()<< "Running trigger tool #"<< ind<< " "<< m_trigToolNames[ind]<< endreq;

    StatusCode sc = m_trigTools[ind]->mutate(th, eh);
    if(sc.isFailure())
    {
      fatal()<< "Trigger Tool "<< m_trigToolNames[ind]<< " failed"<< endreq;
      delete th; 
      th = 0;
      return StatusCode::FAILURE;
    }
  }
  if(m_triggerSorter->mutate(th, eh).isFailure())
  {
    fatal()<< "Trigger Tool "<< m_trigFilterName<< " "<< m_triggerSorter<< " failed"<< endreq;
    return StatusCode::FAILURE;
  }
    
  return StatusCode::SUCCESS;
}

StatusCode DigitizeAlg::ReadoutSim(const ElecHeader& eh, const SimTrigHeader& th, SimReadoutHeader*& srh)
{
  srh->setContext(eh.context());
  srh->setTimeStamp(eh.timeStamp());
  srh->setEarliest(eh.earliest());
  srh->setLatest(eh.latest());

  debug()<< "ReadoutSim time: "<< endreq;
  debug()<< srh->earliest().GetSec()<< "--"<< srh->earliest().GetNanoSec()<< " "
        << srh->latest().GetSec()<< "--"<< srh->latest().GetNanoSec()<< endreq;
    
  //Fill ReadoutHeader With TriggerDataPackage
  //With this way of doing it the individual readout tools theoretically
  //No longer need to know about the trigger header.
  debug()<< "Creating Trigger Data Packages."<< endreq;
  //Empty trigger data package object (use tool to fill it)
  vector<ReadoutTriggerDataPkg*> trigPkgs;
  //Fill trigger data packages 
  StatusCode sc = m_roTrigPackTool->fillDataPakages(trigPkgs, th);
  if(sc.isFailure())
  {
    fatal()<< "Readout Tool "<< m_roTrigPackName<< " "<< m_roTrigPackTool<< " failed"<< endreq;
    return StatusCode::FAILURE;
  }
  
  //Let each tool do its thing on the event
  for(vector<IROsReadoutTool*>::size_type ind = 0; ind < m_roTools.size(); ++ind)
  {
    debug()<< "Running readout tool #"<< ind<< " "<< m_roToolNames[ind]<< endreq;

    sc = m_roTools[ind]->mutate(srh, trigPkgs, eh);
    if(sc.isFailure())
    {
      fatal()<< "Readout Tool "<< m_roToolNames[ind]<< " failed"<< endreq;
      delete srh; 
      srh = 0;
      return StatusCode::FAILURE;
    }
  }

  return StatusCode::SUCCESS;
}

StatusCode DigitizeAlg::TrigExecute(const ElecHeader* eh, SimReadoutHeader*& srh)
{
  debug()<< "Running Trigger and Readout Simulation!"<< endreq;

  const ElecPulseHeader::PulseCollectionMap& pcmap = eh->pulseHeader()->pulseCollection();
  debug()<< "TrigExecute() get "<< pcmap.size()<< " pulse collections:"<< endreq;
  ElecPulseHeader::PulseCollectionMap::const_iterator it, done = pcmap.end();
  for(it = pcmap.begin(); it != done; ++it)
  {
    debug()<< "detector: "<< it->first<< " has "<< it->second->pulses().size()<< endreq;
  }

  SimTrigHeader* th = new SimTrigHeader();
  StatusCode sc = this->TrigSim(*eh, th);
  if(sc.isFailure() || !th)
  {
    error()<< "TrigExecute(): Failed to run TrigSim!"<< endreq;
    return StatusCode::FAILURE;
  }
  
  //put(th, th->defaultLocation());
  //set input headers
  //vector<const IHeader*> iElecHeaders;
  //iElecHeaders.push_back(eh);
  //th->setInputHeaders(iElecHeaders);
  //set execution number
  th->setExecNumber(srh->execNumber());

  debug()<< "exec num "<< srh->execNumber()<< endreq;

  sc = this->ReadoutSim(*eh, *th, srh);
  if(sc.isFailure() || !srh)
  {
    error()<< "TrigExecute(): Failed to run ReadoutSim"<< endreq;
    return StatusCode::FAILURE;
  }
  // set input headers
  //this->AppendInputHeader(th);
  delete th;
  
  return StatusCode::SUCCESS;
}

StatusCode DigitizeAlg::finalize()
{
  debug()<< "DigitizeAlg::finalize()...."<< endreq;
  this->DybAlgorithm<SimReadoutHeader>::finalize();
  return StatusCode::SUCCESS;
}

---