/home/bes3soft/bes3soft/Boss/7.0.2/dist/7.0.2/Event/EventNavigator/EventNavigator-00-01-03/src/BesNavigatorInit.cxx

Go to the documentation of this file.

#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/SmartDataPtr.h"

#include "EventNavigator/EventNavigator.h"
#include "EventNavigator/BesNavigatorInit.h"

// Monte-Carlo data
#include "McTruth/McParticle.h"
#include "McTruth/EmcMcHit.h"
#include "McTruth/MdcMcHit.h"
#include "McTruth/MucMcHit.h"
#include "McTruth/TofMcHit.h"

// MDC reconstructed data
#include "MdcRecEvent/RecMdcTrack.h"
#include "MdcRecEvent/RecMdcKalTrack.h"
#include "MdcRecEvent/RecMdcHit.h"

// EMC reconstructed data
#include "EmcRecEventModel/RecEmcShower.h"

// TOF reconstructed data
#include "TofRecEvent/RecTofTrack.h"
#include "TofRecEvent/RecBTofHit.h"
#include "TofRecEvent/RecETofHit.h"

// MUC reconstructed data
#include "MucRecEvent/RecMucTrack.h"
#include "MucRecEvent/MucRecHit.h"

// Digi information
#include "EmcRecEventModel/RecEmcDigit.h"
#include "MdcRawEvent/MdcDigi.h"

using namespace std;
using namespace EventModel;
using namespace Event;


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" << endreq;      
      m_navigator = new EventNavigator;
      m_navigator->setMdcCut(m_mdcCut);
    }
  else 
    {
      log << MSG::DEBUG << "Found EventNavigator (read from DST)" << endreq;
      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" << endreq;
      return StatusCode::FAILURE;
    }
  else
    log << MSG::DEBUG << "McParticleCol retrieved of size "<< mcParticles->size() << endreq;
  
  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" << endreq;
          return StatusCode::FAILURE;
        }                
      else
        log << MSG::INFO << "EventNavigator registered successfully in TDS" << endreq;
    }

   if( log.level() < 3 )
   {  
       log << MSG::DEBUG << "Write EventNavigator" << endreq;
       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() << endreq;
      
      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" << endreq;
    }

  // 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() << endreq;

      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" << endreq;
  
  
  // ****** 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() << endreq;
      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" << endreq;
    }
  
  // MDC Kalman Tracks
  SmartDataPtr<RecMdcKalTrackCol> mdcKalTracks(eventSvc(),"/Event/Recon/RecMdcKalTrackCol");
  if( mdcKalTracks )
    {
      log << MSG::DEBUG << "MdcKalTrackCol retrieved of size "<< mdcKalTracks->size() << endreq;
      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" << endreq;
    }
  
  //****** 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() << endreq;

      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" << endreq;
    }

  // ****** 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() << endreq;
      
      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" << endreq;
    }
  

  //****** 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()
{
}

---