NavigationTestAlg.cxx

Go to the documentation of this file.

#include "EventNavigator/NavigationTestAlg.h"
 
#include "EventNavigator/EventNavigator.h"
 
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/SmartDataPtr.h"
#include <cmath>
// #include "PartPropSvc/PartPropSvc.h"
 
// Monte-Carlo data
#include "McTruth/McParticle.h"
#include "McTruth/MdcMcHit.h"
 
// MDC reconstructed data
#include "MdcRecEvent/RecMdcHit.h"
#include "MdcRecEvent/RecMdcTrack.h"
 
// EMC reconstructed data
#include "EmcRecEventModel/RecEmcShower.h"
 
#include "TFile.h"
#include "TH1F.h"
#include "TH2F.h"
 
using namespace std;
using namespace Event;
 
DECLARE_COMPONENT( NavigationTestAlg )
NavigationTestAlg::NavigationTestAlg( const std::string& name, ISvcLocator* pSvcLocator )
    : Algorithm( name, pSvcLocator ) {
  declareProperty( "OutputFileName", m_fout = "histo.root" );
}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
 
StatusCode NavigationTestAlg::initialize() {
  MsgStream log( msgSvc(), name() );
 
  log << MSG::INFO << " Creating histograms" << endmsg;
 
  // create histograms
  m_histo.push_back( new TH1F( "deltap", "#delta p, true momentum - rec momentum, MeV/c", 600,
                               -100, 500 ) ); // 0
 
  m_histo.push_back(
      new TH1F( "n_parents", "Number of parents for reconstructed track", 20, 0, 20 ) ); // 1
  m_histo.push_back(
      new TH1F( "n_tracks", "Number of tracks for one MC particle", 20, 0, 20 ) ); // 2
 
  m_histo.push_back( new TH1F( "hit_parent", "Number of parent MC hits for reconstructed hit",
                               20, -10, 10 ) ); // 3
  m_histo.push_back(
      new TH1F( "hit_track", "Number of reconstructed hits for MC hit", 20, -10, 10 ) ); // 4
 
  m_histo.push_back( new TH1F( "true_p", "true momentum, MeV/c", 1500, 0, 1500 ) ); // 5
  m_histo.push_back(
      new TH1F( "rec_p", "reconstructed momentum, MeV/c", 1500, 0, 1500 ) ); // 6
 
  m_histo.push_back(
      new TH1F( "n_emcparents", "# of parents for reconstructed shower", 20, 0, 20 ) ); // 7
  m_histo.push_back(
      new TH1F( "n_showers", "# of showers for one true particle", 20, 0, 20 ) ); // 8
 
  m_histo.push_back( new TH1F( "pdg_miss", "Pdg of particles with # trk==0 and #showers > 0",
                               6000, -3000, 3000 ) ); // 9
 
  m_histo.push_back(
      new TH1F( "dE_true", "True energy - shower energy", 40, -2000, 2000 ) ); // 10
  m_histo.push_back(
      new TH1F( "dE_rec", "shower energy - sqrt(Prec**2+m_true**2)", 40, -2000, 2000 ) ); // 11
 
  m_histo.push_back( new TH1F( "E_true", "True energy, MeV", 100, 0, 1000 ) ); // 12
  m_histo.push_back( new TH1F( "E_rec", "Rec energy, MeV", 100, 0, 1000 ) );   // 13
 
  m_histo2.push_back(
      new TH2F( "deltap_p", "#delta p/p vs. p, MeV", 100, 0, 1000, 20, 0, 2 ) );
 
  return StatusCode::SUCCESS;
}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
 
StatusCode NavigationTestAlg::execute() {
  MsgStream log( msgSvc(), name() );
 
  // Get EventNavigator from the TDS
  SmartDataPtr<EventNavigator> navigator( eventSvc(), "/Event/Navigator" );
  if ( !navigator )
  {
    log << MSG::ERROR << " Unable to retrieve EventNavigator" << endmsg;
    return StatusCode::FAILURE;
  }
 
  log << MSG::INFO << "EventNavigator object" << endmsg;
  navigator->Print();
 
  // ======== Analyse Monte-Carlo information using navigator
  log << MSG::INFO
      << "=======================================================================" << endmsg;
  log << MSG::INFO
      << "MC Particles ==============================================================="
      << endmsg;
 
  // Get McParticle collection
  SmartDataPtr<McParticleCol> mcParticles( eventSvc(), "/Event/MC/McParticleCol" );
  if ( !mcParticles )
  {
    log << MSG::ERROR << " Unable to retrieve McParticleCol" << endmsg;
    return StatusCode::FAILURE;
  }
 
  // For each McParticle
  for ( McParticleCol::iterator it = mcParticles->begin(); it != mcParticles->end(); it++ )
  {
    // get charge of McParticles
    int pdg_code = ( *it )->particleProperty();
 
    // get true momentum value
    double true_mom = ( *it )->initialFourMomentum().vect().mag();
 
    // fill true momentum
    m_histo[5]->Fill( true_mom );
 
    log << MSG::INFO << "Retrieved McParticle # " << ( *it )->trackIndex() << " PDG "
        << pdg_code << " of true momentum " << true_mom << " GeV/c " << endmsg;
 
    // get list of reconstructed tracks, which correspond to this particle, using
    // EventNavigator
    RecMdcTrackVector    tracks  = navigator->getMdcTracks( *it );
    RecMdcKalTrackVector ktracks = navigator->getMdcKalTracks( *it );
 
    // fill the distribution of number of rec. tracks corresponding to single particle
    m_histo[2]->Fill( tracks.size() );
 
    log << MSG::INFO << "  Found " << tracks.size() << " tracks and " << ktracks.size()
        << " Kalman tracks" << endmsg;
 
    // for each track
    for ( unsigned int i = 0; i < ktracks.size(); i++ )
    {
      // Hep3Vector p = momentum(tracks[i]);
      double momentum = ktracks[i]->p();
 
      log << MSG::INFO << "\t Track # " << i << " id = " << ktracks[i]->trackId() << " Pid "
          << ktracks[i]->getPidType() << " Ptot " << momentum << " GeV/c" << endmsg;
 
      // fill difference of true momentum and rec momentum for the _same_ track
      m_histo[0]->Fill( true_mom - momentum );
 
      // fill delta_p/p wrt p for the _same_ track
      m_histo2[0]->Fill( true_mom, fabs( true_mom - momentum ) / true_mom );
    }
 
    // get list of reconstructed showers, which correspond to this particle, using
    // EventNavigator
    RecEmcShowerVector showers = navigator->getEmcRecShowers( *it );
 
    m_histo[8]->Fill( showers.size() );
 
    log << MSG::INFO << "  Found " << showers.size() << " showers" << endmsg;
 
    for ( unsigned int i = 0; i < showers.size(); i++ )
    {
      double true_energy = ( *it )->initialFourMomentum().e();
      double rec_energy  = showers[i]->energy() * 1000; // MeV
 
      log << MSG::INFO << "\t Shower # " << i
          << " Id = " << showers[i]->getShowerId().get_value()
          << " E = " << showers[i]->energy() * 1000 << " MeV " << endmsg;
 
      m_histo[12]->Fill( true_energy );
      m_histo[13]->Fill( rec_energy );
    }
  }
 
  log << MSG::INFO
      << "MDC Tracks ==============================================================" << endmsg;
  // ======== Analyse reconstructed track information using navigator
 
  SmartDataPtr<RecMdcKalTrackCol> mdcKalTracks( eventSvc(), "/Event/Recon/RecMdcKalTrackCol" );
  if ( !mdcKalTracks )
  {
    log << MSG::ERROR << " Unable to retrieve MdcKalTrackCol" << endmsg;
    return StatusCode::SUCCESS;
  }
 
  for ( RecMdcKalTrackCol::iterator it = mdcKalTracks->begin(); it != mdcKalTracks->end();
        it++ )
  {
    McParticleVector particles = navigator->getMcParticles( *it );
 
    log << MSG::INFO << "Retrieved " << particles.size()
        << " McParticles for for MdcKalTrack # " << ( *it )->trackId()
        << " of reconstructed momentum " << ( *it )->p()
        << " GeV/c (PID=" << ( *it )->getPidType() << ")" << endmsg;
 
    // fill rec momentum
    m_histo[6]->Fill( ( *it )->p() );
 
    // fill the distribution of number of parent particles corresponding to this track
    m_histo[1]->Fill( particles.size() );
 
    // for each parent particle
    for ( unsigned int i = 0; i < particles.size(); i++ )
    {
      int pdg_code = particles[i]->particleProperty();
 
      // get true momentum value
      double true_mom = particles[i]->initialFourMomentum().vect().mag();
 
      // get relevance
      int relevance = navigator->getMcParticleRelevance( *it, particles[i] );
 
      // dump particle names, momenta and their relevance
      log << MSG::INFO << "\t" << pdg_code << " momentum " << true_mom << " relevance "
          << relevance << endmsg;
    }
  }
 
  // ======== Analyse reconstructed shower information using navigator
  log << MSG::INFO
      << "EMC Showers =============================================================="
      << endmsg;
  SmartDataPtr<RecEmcShowerCol> emcShowers( eventSvc(), "/Event/Recon/RecEmcShowerCol" );
  if ( !emcShowers )
  {
    log << MSG::ERROR << " Unable to retrieve EmcRecShowerCol" << endmsg;
    return StatusCode::SUCCESS;
  }
 
  int ij = 0;
  for ( RecEmcShowerCol::iterator it = emcShowers->begin(); it != emcShowers->end(); it++ )
  {
    McParticleVector particles = navigator->getMcParticles( *it );
    log << MSG::INFO << "Retrieved McParticles for EmcShower # " << ij++ << " Id "
        << ( *it )->getShowerId().get_value() << " Energy = " << ( *it )->energy() << endmsg;
 
    for ( unsigned int i = 0; i < particles.size(); i++ )
    {
      int pdg_code = particles[i]->particleProperty();
 
      // get true energy
      double true_e = particles[i]->initialFourMomentum().e();
 
      // get relevance
      int relevance = navigator->getMcParticleRelevance( *it, particles[i] );
 
      log << "\t Particle " << i << " PDG " << pdg_code << " E=" << true_e << " Relevance "
          << relevance << endmsg;
    }
 
    // fill the distribution of number of parent particles corresponding to this track
    m_histo[7]->Fill( particles.size() );
  }
 
  log << MSG::INFO
      << "=======================================================================" << endmsg;
 
  return StatusCode::SUCCESS;
}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
 
StatusCode NavigationTestAlg::finalize() {
 
  // Save histograms to ROOT file
  TFile of( m_fout.c_str(), "RECREATE" );
  of.cd();
  for ( vector<TH1F*>::iterator it = m_histo.begin(); it != m_histo.end(); it++ )
    ( *it )->Write();
  for ( vector<TH2F*>::iterator it = m_histo2.begin(); it != m_histo2.end(); it++ )
    ( *it )->Write();
  of.Close();
 
  return StatusCode::SUCCESS;
}