BesNavigatorInit.cxx

Go to the documentation of this file.

#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/SmartDataPtr.h"
 
#include "EventNavigator/BesNavigatorInit.h"
#include "EventNavigator/EventNavigator.h"
 
// Monte-Carlo data
#include "McTruth/EmcMcHit.h"
#include "McTruth/McParticle.h"
#include "McTruth/MdcMcHit.h"
#include "McTruth/MucMcHit.h"
#include "McTruth/TofMcHit.h"
 
// MDC reconstructed data
#include "MdcRecEvent/RecMdcHit.h"
#include "MdcRecEvent/RecMdcKalTrack.h"
#include "MdcRecEvent/RecMdcTrack.h"
 
// EMC reconstructed data
#include "EmcRecEventModel/RecEmcShower.h"
 
// TOF reconstructed data
#include "TofRecEvent/RecBTofHit.h"
#include "TofRecEvent/RecETofHit.h"
#include "TofRecEvent/RecTofTrack.h"
 
// MUC reconstructed data
#include "MucRecEvent/MucRecHit.h"
#include "MucRecEvent/RecMucTrack.h"
 
// Digi information
#include "EmcRecEventModel/RecEmcDigit.h"
#include "MdcRawEvent/MdcDigi.h"
 
using namespace std;
using namespace EventModel;
using namespace Event;
 
/////////////////////////////////////////////////////////////////////////////
DECLARE_COMPONENT( BesNavigatorInit )
BesNavigatorInit::BesNavigatorInit( const std::string& name, ISvcLocator* pSvcLocator )
    : Algorithm( name, pSvcLocator ) {
  // Subdetector flags
  declareProperty( "FillMdcInfo", b_fillMdc = true );
  declareProperty( "FillTofInfo", b_fillTof = true );
  declareProperty( "FillEmcInfo", b_fillEmc = true );
  declareProperty( "FillMucInfo", b_fillMuc = true );
  declareProperty( "MdcCut", m_mdcCut = 7 );
  declareProperty( "PrepareHitMaps", b_fillHitMaps = false );
  state = SIMU;
}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
 
StatusCode BesNavigatorInit::initialize() { return StatusCode::SUCCESS; }
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
 
StatusCode BesNavigatorInit::execute() {
  MsgStream log( msgSvc(), name() );
 
  state = SIMU;
 
  SmartDataPtr<EventNavigator> nav( eventSvc(), "/Event/Navigator" );
  if ( !nav )
  {
    log << MSG::DEBUG << "Create EventNavigator" << endmsg;
    m_navigator = new EventNavigator;
    m_navigator->setMdcCut( m_mdcCut );
  }
  else
  {
    log << MSG::DEBUG << "Found EventNavigator (read from DST)" << endmsg;
    m_navigator = nav;
    state       = RECO;
    if ( log.level() < 3 ) m_navigator->Print();
  }
 
  SmartDataPtr<McParticleCol> mcParticles( eventSvc(), "/Event/MC/McParticleCol" );
  if ( !mcParticles )
  {
    log << MSG::ERROR << "Unable to retrieve McParticleCol" << endmsg;
    return StatusCode::FAILURE;
  }
  else
    log << MSG::DEBUG << "McParticleCol retrieved of size " << mcParticles->size() << endmsg;
 
  for ( McParticleCol::const_iterator it = mcParticles->begin(); it != mcParticles->end();
        it++ )
  { m_navigator->addIdLink( ( *it )->trackIndex(), *it ); }
 
  // store MDC relations
  if ( b_fillMdc ) fillMdcInfo();
 
  // store EMC relations
  if ( b_fillEmc ) fillEmcInfo();
 
  // store TOF relations
  if ( b_fillTof ) fillTofInfo();
 
  // store MUC relations
  if ( b_fillMuc ) fillMucInfo();
 
  // Register EventNavigator to Gaudi TDS
  if ( state == SIMU )
  {
    StatusCode sc = eventSvc()->registerObject( "/Event/Navigator", m_navigator );
    if ( sc != StatusCode::SUCCESS )
    {
      log << MSG::ERROR << "Could not register EventNavigator" << endmsg;
      return StatusCode::FAILURE;
    }
    else log << MSG::INFO << "EventNavigator registered successfully in TDS" << endmsg;
  }
 
  if ( log.level() < 3 )
  {
    log << MSG::DEBUG << "Write EventNavigator" << endmsg;
    m_navigator->Print();
  }
 
  return StatusCode::SUCCESS;
}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
 
StatusCode BesNavigatorInit::finalize() { return StatusCode::SUCCESS; }
 
void BesNavigatorInit::fillMdcInfo() {
  MsgStream log( msgSvc(), name() );
 
  //****** Process hits ******
  // If simulated hits are there, fill index map (read from DST if no hits in TDS)
  SmartDataPtr<MdcMcHitCol> mdcMcHits( eventSvc(), "/Event/MC/MdcMcHitCol" );
  if ( mdcMcHits )
  {
    log << MSG::DEBUG << "MdcMcHitsCol retrieved of size " << mdcMcHits->size() << endmsg;
 
    if ( mdcMcHits->size() != 0 ) m_navigator->getMcMdcMcHitsIdx().clear();
 
    for ( MdcMcHitCol::const_iterator it = mdcMcHits->begin(); it != mdcMcHits->end(); it++ )
    {
      // fill relation MdcHit id <-> McParticle
      m_navigator->getMcMdcMcHitsIdx().insert(
          pair<int, int>( ( *it )->identify().get_value(), ( *it )->getTrackIndex() ) );
 
      //     m_navigator->addIdLink( (*it)->identify().get_value(), *it );
    }
  }
  else { log << MSG::DEBUG << "Unable to retrieve MdcMcHitCol" << endmsg; }
 
  // If reconstructed hits are there, fill index map (read from DST if no hits in TDS)
  SmartDataPtr<RecMdcHitCol> mdcRecHits( eventSvc(), "/Event/Recon/RecMdcHitCol" );
  if ( mdcRecHits )
  {
    log << MSG::DEBUG << "MdcRecHitCol retrieved of size " << mdcRecHits->size() << endmsg;
 
    if ( mdcRecHits->size() != 0 ) m_navigator->getMcMdcTracksIdx().clear();
 
    IndexMap& mchits = m_navigator->getMcMdcMcHitsIdx();
    for ( RecMdcHitCol::iterator it = mdcRecHits->begin(); it != mdcRecHits->end(); it++ )
    {
      int mdcId = ( ( *it )->getMdcId() ).get_value();
      //      m_navigator->addIdLink( mdcId, *it);
 
      const pair<IndexMap::const_iterator, IndexMap::const_iterator> range =
          mchits.equal_range( mdcId );
      for ( IndexMap::const_iterator jt = range.first; jt != range.second; jt++ )
      {
        m_navigator->getMcMdcTracksIdx().insert(
            pair<int, int>( ( *jt ).second, ( *it )->getTrkId() ) );
      }
    }
  }
  else log << MSG::DEBUG << "Unable to retrieve RecMdcHitCol" << endmsg;
 
  // ****** Get MDC top level objects from TDS ******
  // MDC Reconstructed Tracks
  SmartDataPtr<RecMdcTrackCol> mdcTracks( eventSvc(), "/Event/Recon/RecMdcTrackCol" );
  if ( mdcTracks )
  {
    log << MSG::DEBUG << "MdcTrackCol retrieved of size " << mdcTracks->size() << endmsg;
    for ( RecMdcTrackCol::const_iterator it = mdcTracks->begin(); it != mdcTracks->end();
          it++ )
      m_navigator->addIdLink( ( *it )->trackId(), *it );
  }
  else { log << MSG::DEBUG << "Unable to retrieve RecMdcTrackCol" << endmsg; }
 
  // MDC Kalman Tracks
  SmartDataPtr<RecMdcKalTrackCol> mdcKalTracks( eventSvc(), "/Event/Recon/RecMdcKalTrackCol" );
  if ( mdcKalTracks )
  {
    log << MSG::DEBUG << "MdcKalTrackCol retrieved of size " << mdcKalTracks->size() << endmsg;
    for ( RecMdcKalTrackCol::const_iterator it = mdcKalTracks->begin();
          it != mdcKalTracks->end(); it++ )
      m_navigator->addIdLink( ( *it )->trackId(), *it );
  }
  else { log << MSG::DEBUG << "Unable to retrieve RecMdcTrackCol" << endmsg; }
 
  //****** Make main relations between top level objects ***
  IndexMap::const_iterator i;
  IndexMap&                index = m_navigator->getMcMdcTracksIdx();
  for ( i = index.begin(); i != index.end(); i++ )
  {
    // i.first = McParticle id; i.second = MdcTrack id
    const McParticle* mcParticle = m_navigator->getMcParticle( ( *i ).first );
 
    if ( mcParticle == 0 ) continue;
 
    const RecMdcTrack* mdcTrack = m_navigator->getMdcTrack( ( *i ).second );
 
    if ( mdcTrack != 0 )
    {
      m_navigator->addLink( mcParticle, mdcTrack );
      m_navigator->addLink( mdcTrack, mcParticle );
    }
 
    const RecMdcKalTrack* mdcKalTrack =
        m_navigator->getMdcKalTrack( ( *i ).second ); // same trkid
    if ( mdcKalTrack != 0 )
    {
      m_navigator->addLink( mcParticle, mdcKalTrack );
      m_navigator->addLink( mdcKalTrack, mcParticle );
    }
  }
}
 
void BesNavigatorInit::fillEmcInfo() {
  MsgStream log( msgSvc(), name() );
 
  //****** Process hits ******
  // If simulated hits are there, fill index map (read from DST if no hits in TDS)
  SmartDataPtr<EmcMcHitCol> emcMcHits( eventSvc(), "/Event/MC/EmcMcHitCol" );
  if ( emcMcHits )
  {
    log << MSG::DEBUG << "EmcMcHitsCol retrieved of size " << emcMcHits->size() << endmsg;
 
    if ( emcMcHits->size() != 0 ) m_navigator->getMcEmcMcHitsIdx().clear();
 
    for ( EmcMcHitCol::const_iterator it = emcMcHits->begin(); it != emcMcHits->end(); it++ )
    {
      // first put Id of crystal, which particle hits initially
      int crystalId    = ( *it )->identify().get_value();
      int McParticleID = 0;
      if ( crystalId != 0 ) // not gamma converted in TOF
      {
        McParticleID = ( *it )->getTrackIndex();
        m_navigator->getMcEmcMcHitsIdx().insert( pair<int, int>( crystalId, McParticleID ) );
        //        m_navigator->addIdLink( crystalId, *it );
      }
 
      // then add all other crystals touched by shower
      const map<Identifier, double>&          hitmap = ( *it )->getHitMap();
      map<Identifier, double>::const_iterator jt;
      for ( jt = hitmap.begin(); jt != hitmap.end(); jt++ )
      {
        int crystalId = ( *jt ).first.get_value();
 
        if ( McParticleID != 0 ) // not gamma converted in TOF
        {
          m_navigator->getMcEmcMcHitsIdx().insert(
              pair<int, int>( crystalId, McParticleID ) ); // CrystalId, McParticleID
        }
        // m_navigator->addIdLink( (*jt).first.get_value(), *it );
      }
    }
  }
  else { log << MSG::DEBUG << "Unable to retrieve EmcMcHitCol" << endmsg; }
 
  // ****** Get EMC top level objects from TDS ******
  // EMC Reconstructed Showers
  SmartDataPtr<RecEmcShowerCol> emcRecShowers( eventSvc(), "/Event/Recon/RecEmcShowerCol" );
  if ( emcRecShowers )
  {
    log << MSG::DEBUG << "RecEmcShowerCol retrieved of size " << emcRecShowers->size()
        << endmsg;
 
    IndexMap& mchits = m_navigator->getMcEmcMcHitsIdx();
 
    for ( RecEmcShowerCol::const_iterator it = emcRecShowers->begin();
          it != emcRecShowers->end(); it++ )
    {
      int showerId = ( *it )->getShowerId().get_value();
      m_navigator->addIdLink( showerId, *it );
 
      // calculate relations only during reco and store them to DST.
      // EmcMcHit tables are removed then.
      if ( !mchits.empty() )
      {
        const RecEmcFractionMap& fractionMap = ( *it )->getFractionMap();
        for ( RecEmcFractionMap::const_iterator jt = fractionMap.begin();
              jt != fractionMap.end(); jt++ )
        {
          int crystalId = ( *jt ).first.get_value();
          const pair<IndexMap::const_iterator, IndexMap::const_iterator> range =
              mchits.equal_range( crystalId ); // McParticles which hit given crystal
 
          for ( IndexMap::const_iterator kt = range.first; kt != range.second; kt++ )
          {
            m_navigator->getMcEmcRecShowersIdx().insert( pair<int, int>(
                ( *kt ).second, showerId ) ); // McParticle id <-> RecShower ShowerId
          }
        }
      }
    }
  }
  else { log << MSG::DEBUG << "Unable to retrieve RecEmcShowerCol" << endmsg; }
 
  //****** Make main relations between top level objects ***
  IndexMap::const_iterator i;
  IndexMap&                index = m_navigator->getMcEmcRecShowersIdx();
  for ( i = index.begin(); i != index.end(); i++ )
  {
    // i.first = McParticle id; i.second = MdcTrack id
    const McParticle* mcParticle = m_navigator->getMcParticle( ( *i ).first );
 
    if ( mcParticle == 0 ) continue;
 
    const RecEmcShowerVector& emcShowers = m_navigator->getEmcRecShowers( ( *i ).second );
 
    if ( !emcShowers.empty() )
    {
      RecEmcShowerVector::const_iterator j;
      for ( j = emcShowers.begin(); j != emcShowers.end(); j++ )
      {
        m_navigator->addLink( mcParticle, *j );
        m_navigator->addLink( *j, mcParticle );
      }
    }
  }
}
 
void BesNavigatorInit::fillTofInfo() {}
 
void BesNavigatorInit::fillMucInfo() {}