BbEmc.cxx

Go to the documentation of this file.

#include "GaudiKernel/IDataProviderSvc.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/PropertyMgr.h"
#include "GaudiKernel/SmartDataPtr.h"
 
#include "EvTimeEvent/RecEsTime.h"
#include "EventModel/Event.h"
#include "EventModel/EventHeader.h"
#include "EventModel/EventModel.h"
#include "McTruth/McParticle.h"
#include "TrigEvent/TrigData.h"
#include "TrigEvent/TrigEvent.h"
 
#include "DstEvent/TofHitStatus.h"
#include "EvtRecEvent/EvtRecEvent.h"
#include "EvtRecEvent/EvtRecTrack.h"
 
#include "CLHEP/Vector/ThreeVector.h"
#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/IHistogramSvc.h"
#include "GaudiKernel/INTupleSvc.h"
#include "GaudiKernel/NTuple.h"
#include "TMath.h"
using CLHEP::Hep3Vector;
#include "CLHEP/Vector/LorentzVector.h"
using CLHEP::HepLorentzVector;
#include "CLHEP/Vector/TwoVector.h"
using CLHEP::Hep2Vector;
#include "CLHEP/Geometry/Point3D.h"
#ifndef ENABLE_BACKWARDS_COMPATIBILITY
typedef HepGeom::Point3D<double> HepPoint3D;
#endif
#include "HltEvent/DstHltInf.h"
#include "HltEvent/HltInf.h"
#include "ParticleID/ParticleID.h"
#include "VertexFit/KinematicFit.h"
#include "VertexFit/SecondVertexFit.h"
#include "VertexFit/VertexFit.h"
 
#include "BbEmc.h"
#include <vector>
typedef std::vector<int>              Vint;
typedef std::vector<HepLorentzVector> Vp4;
const double                          velc = 299.792458;
const double xmass[5] = { 0.000511, 0.105658, 0.139570, 0.493677, 0.938272 };
const double pai      = 3.1415926;
 
static long m_cout_all( 0 ), m_cout_tracks( 0 ), m_cout_good( 0 ), m_cout_dang( 0 );
////////////////////////////////////////////////////////////////////
DECLARE_COMPONENT( BbEmc )
BbEmc::BbEmc( const std::string& name, ISvcLocator* pSvcLocator )
    : Algorithm( name, pSvcLocator ) {
  // Declare the properties
  // control flag
  declareProperty( "hist", m_hist = false );
  declareProperty( "Trigger", m_trigger_flag = true );
  declareProperty( "Hlt", m_hlt_flag = true );
  declareProperty( "Estime", m_est_flag = true );
  // RecMdcKalTrack or RecMdcTrack
  declareProperty( "KalTrk", m_kalTrk_flag = true );
  // Good shower selection cut
  declareProperty( "EneCut", m_energy_cut = 1.2 );
  declareProperty( "MomCut", m_mom_cut = 0.04 );
  declareProperty( "DangCut", m_dangCut = 2.5 );
  // Mdc Track vertex cut
  declareProperty( "Vr0cut", m_vr0cut = 1.0 );
  declareProperty( "Vz0cut", m_vz0cut = 5.0 );
}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
StatusCode BbEmc::initialize() {
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "in initialize()" << endmsg;
 
  StatusCode status;
 
  if ( m_hist )
  {
    NTuplePtr nt1( ntupleSvc(), "FILEBbEmc/bbEmc" );
    if ( nt1 ) m_tuple1 = nt1;
    else
    {
      m_tuple1 = ntupleSvc()->book( "FILEBbEmc/bbEmc", CLID_ColumnWiseTuple,
                                    "BbEmc N-Tuple example" );
      if ( m_tuple1 )
      {
        status = m_tuple1->addItem( "run", m_run );
        status = m_tuple1->addItem( "event", m_event );
        status = m_tuple1->addItem( "numCtrk", m_num_Ctrk );
        status = m_tuple1->addItem( "numNtrk", m_num_Ntrk );
        status = m_tuple1->addItem( "numGtrk", m_num_Gtrk );
        // trigger and hlt
        status = m_tuple1->addItem( "trigindex", m_trig_index, 0, 48 );
        status = m_tuple1->addIndexedItem( "trigcond", m_trig_index, m_trig_cond );
        status = m_tuple1->addIndexedItem( "trigchan", m_trig_index, m_trig_chan );
        status = m_tuple1->addItem( "timewindow", m_trig_timewindow );
        status = m_tuple1->addItem( "timetype", m_trig_timetype );
        status = m_tuple1->addItem( "hlttype", m_hlt_type );
        // estime
        status = m_tuple1->addItem( "estime", m_est_start );
        status = m_tuple1->addItem( "status", m_est_status );
        status = m_tuple1->addItem( "quality", m_est_quality );
        // Emc shower
        status = m_tuple1->addItem( "dang", m_dang );
        status = m_tuple1->addItem( "index", m_index, 0, 2 );
        status = m_tuple1->addIndexedItem( "theta", m_index, m_theta );
        status = m_tuple1->addIndexedItem( "phi", m_index, m_phi );
        status = m_tuple1->addIndexedItem( "ene", m_index, m_ene );
      }
      else
      {
        log << MSG::ERROR << "Cannot book N-tuple theone" << endmsg;
        return StatusCode::FAILURE;
      }
    }
 
    log << MSG::INFO << "end initialize()" << endmsg;
    // finish book
  } // end of m_hist
 
  return StatusCode::SUCCESS;
}
 
//*************************************************
StatusCode BbEmc::execute() {
 
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "in execute()" << endmsg;
  m_cout_all++;
 
  // save the events passed selection to a new file
  setFilterPassed( false );
 
  SmartDataPtr<Event::EventHeader> eventHeader( eventSvc(), "/Event/EventHeader" );
  int                              runNum   = eventHeader->runNumber();
  int                              eventNum = eventHeader->eventNumber();
  if ( m_cout_all % 1000 == 0 )
  { std::cout << name() << "::" << m_cout_all << " events executed" << std::endl; }
 
  SmartDataPtr<EvtRecEvent> evtRecEvent( eventSvc(), "/Event/EvtRec/EvtRecEvent" );
  log << MSG::DEBUG << "ncharg, nneu, tottks = " << evtRecEvent->totalCharged() << " , "
      << evtRecEvent->totalNeutral() << " , " << evtRecEvent->totalTracks() << endmsg;
 
  SmartDataPtr<EvtRecTrackCol> evtRecTrkCol( eventSvc(), "/Event/EvtRec/EvtRecTrackCol" );
  if ( ( evtRecEvent->totalTracks() < 2 ) || ( evtRecEvent->totalTracks() > 30 ) )
  { return StatusCode::SUCCESS; }
  m_cout_tracks++;
 
  // High level trigger
  int hlt_type = 0;
  if ( m_hlt_flag )
  {
    SmartDataPtr<HltInf> hlt( eventSvc(), "/Event/Hlt/HltInf" );
    DstHltInf*           aHlt;
    if ( !hlt )
    {
      log << MSG::WARNING << "Could not find Event HltInf, try DstHltInf now" << endmsg;
      SmartDataPtr<DstHltInf> hltDst( eventSvc(), "/Event/Hlt/DstHltInf" );
      if ( !hltDst )
      {
        log << MSG::ERROR << "Could not find Event DstHltInf too, please re-generate data"
            << endmsg;
        return StatusCode::FAILURE;
      }
      else { aHlt = hltDst; }
    }
    else { aHlt = hlt; }
    uint32_t temp_type( aHlt->getEventType() );
    int      mask( 1 );
    for ( int i = 0; i < 32; i++ )
    {
      if ( temp_type & mask )
      {
        hlt_type = i;
        break;
      }
      temp_type >>= 1;
    }
  }
 
  // Trigger
  SmartDataPtr<TrigData> trigData( eventSvc(), "/Event/Trig/TrigData" );
  if ( m_trigger_flag & runNum > 0 )
  {
    if ( !trigData )
    {
      log << MSG::FATAL << "Could not find TrigData in TDS" << endmsg;
      return StatusCode::FAILURE;
    }
    // Print trigger information once:
    log << MSG::DEBUG << "Trigger conditions: " << endmsg;
    for ( int trig_id = 0; trig_id < 48; trig_id++ )
    {
      log << MSG::DEBUG << "i:" << trig_id << "  name:" << trigData->getTrigCondName( trig_id )
          << "  cond:" << trigData->getTrigCondition( trig_id ) << endmsg;
    }
  }
 
  // Event start time
  SmartDataPtr<RecEsTimeCol> recEstimeCol( eventSvc(), "/Event/Recon/RecEsTimeCol" );
  if ( m_est_flag & runNum > 0 )
  {
    if ( !recEstimeCol )
    {
      log << MSG::WARNING << "Can not get RecEsTimeCol" << endmsg;
      return StatusCode::SUCCESS;
    }
    log << MSG::DEBUG << "size of EsTime: " << recEstimeCol->size() << endmsg;
  }
 
  // good shower selection
  Vp4              vGood;
  HepLorentzVector m_lv_ele;
  for ( int i = 0; i < evtRecEvent->totalTracks(); i++ )
  {
    if ( i >= evtRecTrkCol->size() ) break;
 
    EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + i;
    if ( !( *itTrk )->isEmcShowerValid() ) continue;
 
    RecEmcShower* emcTrk = ( *itTrk )->emcShower();
    if ( emcTrk->energy() < m_energy_cut ) continue;
 
    Hep3Vector emcpos( emcTrk->x(), emcTrk->y(), emcTrk->z() );
    m_lv_ele.setVect( emcpos );
    m_lv_ele.setE( emcTrk->energy() );
    vGood.push_back( m_lv_ele );
  }
 
  // num of good showers = 2
  if ( vGood.size() - 2 ) return StatusCode::SUCCESS;
  m_cout_good++;
 
  // angle between two showers
  double dang = vGood[0].vect().angle( vGood[1].vect() );
  if ( dang < m_dangCut ) return StatusCode::SUCCESS;
  m_cout_dang++;
 
  if ( vGood[0].e() > vGood[1].e() ) swap( vGood[0], vGood[1] );
 
  // check x0, y0, z0, r0
  // suggest cut: |z0|<5 && r0<1 (cm)
  double d0, z0, cosTheta, phi, mom;
  Vint   iGood;
  iGood.clear();
  for ( int i = 0; i < evtRecEvent->totalCharged(); i++ )
  {
    double              m_vz, m_vr;
    EvtRecTrackIterator itTrk;
    if ( m_kalTrk_flag )
    {
      itTrk = evtRecTrkCol->begin() + i;
      if ( !( *itTrk )->isMdcKalTrackValid() ) continue;
      RecMdcKalTrack* mdcTrk = ( *itTrk )->mdcKalTrack();
      m_vz                   = mdcTrk->z();
      m_vr                   = mdcTrk->r();
    }
    else
    {
      if ( !( *itTrk )->isMdcTrackValid() ) continue;
      RecMdcTrack* mdcTrk = ( *itTrk )->mdcTrack();
      m_vz                = mdcTrk->z();
      m_vr                = mdcTrk->r();
    }
    if ( fabs( m_vz ) >= m_vz0cut ) continue;
    if ( m_vr >= m_vr0cut ) continue;
    iGood.push_back( i );
  }
 
  if ( m_hist )
  {
    m_run   = runNum;
    m_event = eventNum;
 
    m_num_Ctrk = evtRecEvent->totalCharged();
    m_num_Ntrk = evtRecEvent->totalNeutral();
    m_num_Gtrk = iGood.size();
    if ( trigData )
    {
      for ( int trig_id = 0; trig_id < 48; trig_id++ )
      {
        m_trig_index = trig_id;
        if ( m_trig_index < 16 )
        { m_trig_chan[m_trig_index] = trigData->getTrigChannel( m_trig_index ); }
        m_trig_cond[m_trig_index] = trigData->getTrigCondition( m_trig_index );
      }
      m_trig_timewindow = trigData->getTimeWindow();
      m_trig_timetype   = trigData->getTimingType();
    }
 
    m_hlt_type = hlt_type;
 
    if ( recEstimeCol )
    {
      m_est_start   = ( *( recEstimeCol->begin() ) )->getTest();
      m_est_status  = ( *( recEstimeCol->begin() ) )->getStat();
      m_est_quality = ( *( recEstimeCol->begin() ) )->getQuality();
    }
 
    m_dang = dang;
    for ( int i = 0; i < 2; i++ )
    {
      m_index          = i;
      m_theta[m_index] = vGood[m_index].vect().theta();
      m_phi[m_index]   = vGood[m_index].vect().phi();
      m_ene[m_index]   = vGood[m_index].e();
    }
 
    m_tuple1->write();
  }
 
  setFilterPassed( true );
 
  return StatusCode::SUCCESS;
}
 
StatusCode BbEmc::finalize() {
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "in finalize()" << endmsg;
 
  std::cout << name() << " total event: " << m_cout_all << std::endl;
  std::cout << name() << " total tracks >= 2, <= 30: " << m_cout_tracks << std::endl;
  std::cout << name() << " good showers = 2: " << m_cout_good << std::endl;
  std::cout << name() << " angle between two showers: " << m_cout_dang << std::endl;
 
  return StatusCode::SUCCESS;
}