EFChargedTrack.cxx

Go to the documentation of this file.

#include "GaudiKernel/MsgStream.h"
// #include "GaudiKernel/AlgFactory.h"
// #include "GaudiKernel/ISvcLocator.h"
// #include "GaudiKernel/IDataProviderSvc.h"
#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/SmartDataPtr.h"
// #include "GaudiKernel/PropertyMgr.h"
#include "EventModel/Event.h"
#include "EventModel/EventHeader.h"
#include "HltAlgorithms/EFChargedTrack.h"
#include "MdcRecEvent/RecMdcTrack.h"
 
using namespace Event;
 
EFChargedTrack::EFChargedTrack( const std::string& name, ISvcLocator* pSvcLocator )
    : IEFAlgorithm( name, pSvcLocator ) {
  int output = ( m_output % 100 ) / 10;
  // declareProperty("OutputLevel",m_output = MSG::NIL);
  MsgStream log( msgSvc(), name );
  msgSvc()->setOutputLevel( name, output );
  m_ntrk  = new CriteriaItemValue;
  m_acol  = new CriteriaItemValue;
  m_mbal  = new CriteriaItemValue;
  m_pmax1 = new CriteriaItemValue;
  m_pmax2 = new CriteriaItemValue;
  m_cost1 = new CriteriaItemValue;
  m_cost2 = new CriteriaItemValue;
  m_vr    = new CriteriaItemValue;
  m_vz    = new CriteriaItemValue;
}
 
EFChargedTrack::~EFChargedTrack() {
  delete m_ntrk;
  delete m_acol;
  delete m_mbal;
  delete m_pmax1;
  delete m_pmax2;
  delete m_cost1;
  delete m_cost2;
  delete m_vr;
  delete m_vz;
}
 
StatusCode EFChargedTrack::initialize() {
 
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "in initialize()" << endmsg;
 
  IEFAlgorithm::initialize();
 
  StatusCode sc;
  sc = m_HltStoreSvc->put( "ntrk", m_ntrk );
  if ( sc.isFailure() )
  {
    log << MSG::ERROR << "m_HltStoreSvc->put(ntrk) wrong" << endmsg;
    return sc;
  }
  sc = m_HltStoreSvc->put( "acol", m_acol );
  if ( sc.isFailure() )
  {
    log << MSG::ERROR << "m_HltStoreSvc->put(acol) wrong" << endmsg;
    return sc;
  }
  sc = m_HltStoreSvc->put( "mbal", m_mbal );
  if ( sc.isFailure() )
  {
    log << MSG::ERROR << "m_HltStoreSvc->put(mbal) wrong" << endmsg;
    return sc;
  }
  sc = m_HltStoreSvc->put( "pmax1", m_pmax1 );
  if ( sc.isFailure() )
  {
    log << MSG::ERROR << "m_HltStoreSvc->put(pmax1) wrong" << endmsg;
    return sc;
  }
  sc = m_HltStoreSvc->put( "pmax2", m_pmax2 );
  if ( sc.isFailure() )
  {
    log << MSG::ERROR << "m_HltStoreSvc->put(pmax2) wrong" << endmsg;
    return sc;
  }
  sc = m_HltStoreSvc->put( "cost1", m_cost1 );
  if ( sc.isFailure() )
  {
    log << MSG::ERROR << "m_HltStoreSvc->put(cost1) wrong" << endmsg;
    return sc;
  }
  sc = m_HltStoreSvc->put( "cost2", m_cost2 );
  if ( sc.isFailure() )
  {
    log << MSG::ERROR << "m_HltStoreSvc->put(cost2) wrong" << endmsg;
    return sc;
  }
  sc = m_HltStoreSvc->put( "vr", m_vr );
  if ( sc.isFailure() )
  {
    log << MSG::ERROR << "m_HltStoreSvc->put(vr) wrong" << endmsg;
    return sc;
  }
  sc = m_HltStoreSvc->put( "vz", m_vz );
  if ( sc.isFailure() )
  {
    log << MSG::ERROR << "m_HltStoreSvc->put(vz) wrong" << endmsg;
    return sc;
  }
 
  log << MSG::DEBUG << "finish initialize()" << endmsg;
 
  return StatusCode::SUCCESS;
}
 
StatusCode EFChargedTrack::execute() {
 
  reset();
 
  MsgStream log( msgSvc(), name() );
  // DataObject*  pObject;
 
  // Part 1: Get the event header, print out event and run number
  SmartDataPtr<Event::EventHeader> eventHeader( eventSvc(), "/Event/EventHeader" );
  if ( !eventHeader )
  {
    log << MSG::FATAL << "Could not find Event Header" << endmsg;
    return ( StatusCode::FAILURE );
  }
 
  // Part 2: Retrieve Track Collection
  SmartDataPtr<RecMdcTrackCol> mdcTrackCol( eventSvc(), "/Event/Recon/RecMdcTrackCol" );
  if ( !mdcTrackCol )
  {
    log << MSG::FATAL << "Could not find Mdc track collection!!" << endmsg;
    return ( StatusCode::FAILURE );
  }
  // Part 3:
  RecMdcTrackCol::iterator iterTrk = mdcTrackCol->begin();
  double                   max1    = 0.;
  double                   max2    = 0.;
  double                   cost1 = -999, cost2 = -999, phi1 = -999, phi2 = -999;
  int    hit1 = -999, hit2 = -999, shit1 = -999, shit2 = -999, charg1 = -999, charg2 = -999;
  double mdcbalance = -999;
 
  unsigned int ntrk  = mdcTrackCol->size();
  double       kappa = -999, tanl = -999, sint = -999;
  double       p = -999, theta = -999, phi = -999;
  double       vz1 = -999, vz2 = -999, vr1 = -999, vr2 = -999;
 
  // double theta1=0,theta2=0;
  for ( ; iterTrk != mdcTrackCol->end(); iterTrk++ )
  {
    // z0=(*iterTrk)->getZ0();
    // x0=(*iterTrk)->getX0();
    // y0=(*iterTrk)->getY0();
    phi   = ( *iterTrk )->helix( 1 );
    kappa = ( *iterTrk )->helix( 2 );
    tanl  = ( *iterTrk )->helix( 4 );
 
    theta = 0.5 * 3.1415926 - atan( tanl );
    sint  = sin( theta );
    // log << MSG::DEBUG << "FastTrk=>("<<x0<<","<<y0<<","<<z0")=>"
    //   <<(*iterTrk)->getId() << ":"<<(*iterTrk)->getPivot()<<":"
    //   <<kappa<<">"<<tanl<<endmsg;
    if ( abs( kappa ) > 0.001 && abs( sint ) > 0.01 ) { p = abs( 1. / kappa ) / sint; }
    else
    {
      p = 1000.;
      log << MSG::WARNING << "FastTrk=>"
          << " kappa=" << kappa << "; sint=" << sint << endmsg;
    }
    if ( p >= max1 )
    {
      max2   = max1;
      cost2  = cost1;
      phi2   = phi1;
      max1   = p;
      cost1  = cos( theta );
      phi1   = phi;
      vz2    = vz1;
      vz1    = ( *iterTrk )->helix( 3 );
      vr2    = vr1;
      vr1    = ( *iterTrk )->helix( 0 );
      hit2   = hit1;
      hit1   = ( *iterTrk )->getNhits();
      shit2  = shit1;
      shit1  = ( *iterTrk )->nster();
      charg2 = charg1;
      charg1 = ( *iterTrk )->charge();
    }
    else if ( p > max2 )
    {
      max2   = p;
      cost2  = cos( theta );
      phi2   = phi;
      vz2    = ( *iterTrk )->helix( 3 );
      vr2    = ( *iterTrk )->helix( 0 );
      hit2   = ( *iterTrk )->getNhits();
      shit2  = ( *iterTrk )->nster();
      charg2 = ( *iterTrk )->charge();
    }
    log << MSG::DEBUG << "p=" << p << ", "
        << "theta=" << theta << ", phi=" << phi << ", vz=" << ( *iterTrk )->helix( 3 )
        << ", vr=" << ( *iterTrk )->helix( 0 ) << endmsg;
    if ( cos( theta ) > 0 ) mdcbalance += 1.;
    else if ( cos( theta ) < 0 ) mdcbalance -= 1.;
  }
  if ( ntrk >= 2 ) mdcbalance /= ntrk;
  else mdcbalance = 1;
 
  double acol = 180.;
  if ( ntrk >= 2 )
  {
    acol = 180. -
           180. / 3.1415926 *
               acos( cos( phi1 ) * sin( acos( cost1 ) ) * cos( phi2 ) * sin( acos( cost2 ) ) +
                     sin( phi1 ) * sin( acos( cost1 ) ) * sin( phi2 ) * sin( acos( cost2 ) ) +
                     cost1 * cost2 );
  }
 
  log << MSG::INFO << "ntrk=" << ntrk << "; mdc balance=" << mdcbalance << "; pmax1=" << max1
      << "; pmax2=" << max2 << "; acol=" << acol << "; cost1=" << cost1 << "; cost2=" << cost2
      << endmsg;
 
  // Part 4: Put the criteria item(s) to HltStoreSvc here
  m_ntrk->setValue( ntrk );
  m_acol->setValue( acol );
  m_mbal->setValue( mdcbalance );
  m_pmax1->setValue( max1 );
  m_pmax2->setValue( max2 );
  m_cost1->setValue( cost1 );
  m_cost2->setValue( cost2 );
  m_vr->setValue( vr1 );
  m_vz->setValue( vz1 );
 
  m_ef->addToEFVec( ntrk, 6 );
  m_ef->appToEFVec( max1, 7 );
  m_ef->appToEFVec( cost1, 8 );
  m_ef->appToEFVec( vz1, 9 );
  m_ef->appToEFVec( phi1, 10 );
  m_ef->appToEFVec( vr1, 11 );
  m_ef->addToEFVec( hit1, 12 );
  m_ef->addToEFVec( shit1, 13 );
  m_ef->appToEFVec( max2, 14 );
  m_ef->appToEFVec( cost2, 15 );
  m_ef->appToEFVec( vz2, 16 );
  m_ef->appToEFVec( phi2, 17 );
  m_ef->appToEFVec( vr2, 18 );
  m_ef->addToEFVec( hit2, 19 );
  m_ef->addToEFVec( shit2, 20 );
  m_ef->appToEFVec( acol, 21 );
  m_ef->setVecBit( true, 0, 4 );
  if ( ntrk == 0 ) m_ef->addToEFVec( 1, 1 );
  else if ( ntrk == 1 ) m_ef->addToEFVec( 8, 1 );
  else if ( ntrk >= 2 ) m_ef->addToEFVec( 16, 1 );
  else m_ef->addToEFVec( 0, 1 );
 
  m_ef->addToEFVec( charg1, 23 );
  m_ef->addToEFVec( charg2, 24 );
 
  m_run = 1;
 
  return StatusCode::SUCCESS;
}
 
StatusCode EFChargedTrack::finalize() {
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "in finalize()" << endmsg;
  return StatusCode::SUCCESS;
}
 
void EFChargedTrack::reset() {
  // Reset
  if ( m_run )
  {
    m_ntrk->reset();
    m_acol->reset();
    m_mbal->reset();
    m_pmax1->reset();
    m_pmax2->reset();
    m_cost1->reset();
    m_cost2->reset();
    m_vr->reset();
    m_vz->reset();
    m_run = 0;
  }
  return;
}