LumTau.cxx

Go to the documentation of this file.

#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/IDataProviderSvc.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/PropertyMgr.h"
#include "GaudiKernel/SmartDataPtr.h"
// #include "VertexFit/IVertexDbSvc.h"
 
#include "EventModel/Event.h"
#include "EventModel/EventHeader.h"
#include "EventModel/EventModel.h"
// #include "McTruth/McParticle.h"
 
// #include "CLHEP/Vector/LorentzVector.h"
// #include "CLHEP/Vector/ThreeVector.h"
// #include "CLHEP/Vector/TwoVector.h"
// #include "DstEvent/TofHitStatus.h"
#include "EvTimeEvent/RecEsTime.h"
#include "EventModel/EventHeader.h"
#include "EvtRecEvent/EvtRecEvent.h"
#include "EvtRecEvent/EvtRecTrack.h"
// #include "GaudiKernel/Bootstrap.h"
// #include "GaudiKernel/IHistogramSvc.h"
// #include "GaudiKernel/INTupleSvc.h"
// #include "GaudiKernel/NTuple.h"
// #include "TMath.h"
// using CLHEP::Hep2Vector;
// using CLHEP::Hep3Vector;
// using CLHEP::HepLorentzVector;
// #include "CLHEP/Geometry/Point3D.h"
// #ifndef ENABLE_BACKWARDS_COMPATIBILITY
// typedef HepGeom::Point3D<double> HepPoint3D;
// #endif
#include "LumTau.h"
// #include "ParticleID/ParticleID.h"
// #include "VertexFit/Helix.h"
// #include "VertexFit/KalmanKinematicFit.h"
// #include "VertexFit/KinematicFit.h"
// #include "VertexFit/SecondVertexFit.h"
// #include "VertexFit/VertexFit.h"
// #include <vector>
// typedef std::vector<int>              Vint;
// typedef std::vector<HepLorentzVector> Vp4;
// const double                          ecms = 3.097;
// const double                          velc = 299.792458;
// const double xmass[5] = { 0.000511, 0.105658, 0.139570, 0.493677, 0.938272 };
const double pai      = 3.1415926;
////////////////////////////////////////////////////////////////////
DECLARE_COMPONENT( LumTau )
LumTau::LumTau( const std::string& name, ISvcLocator* pSvcLocator )
    : Algorithm( name, pSvcLocator ) {}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
StatusCode LumTau::initialize() {
  MsgStream log( msgSvc(), name() );
 
  log << MSG::INFO << "in initialize()" << endmsg;
 
  StatusCode status;
 
  status = service( "InjSigIntervalSvc", m_intervalSvc );
  if ( status != StatusCode::SUCCESS )
  { log << MSG::FATAL << "can not use InjSigIntervalSvc" << endmsg; }
 
  NTuplePtr nt2( ntupleSvc(), "LumTau/event" );
  if ( nt2 ) m_tuple2 = nt2;
  else
  {
    m_tuple2 =
        ntupleSvc()->book( "LumTau/event", CLID_ColumnWiseTuple, "Bhabha N-Tuple signal" );
    if ( m_tuple2 )
    {
      status = m_tuple2->addItem( "topup", m_topup );
      status = m_tuple2->addItem( "run", m_run );
      status = m_tuple2->addItem( "rec", m_rec );
      status = m_tuple2->addItem( "time", m_time );
      status = m_tuple2->addItem( "etsT1", m_etsT1 );
      status = m_tuple2->addItem( "etot", m_etot );
      status = m_tuple2->addItem( "e1", m_e1 );
      status = m_tuple2->addItem( "e2", m_e2 );
      status = m_tuple2->addItem( "costht1", m_costht1 );
      status = m_tuple2->addItem( "costht2", m_costht2 );
      status = m_tuple2->addItem( "dltphi", m_dltphi );
      status = m_tuple2->addItem( "dlttht", m_dlttht );
      status = m_tuple2->addItem( "phi1", m_phi1 );
      status = m_tuple2->addItem( "phi2", m_phi2 );
    }
    else
    {
      log << MSG::ERROR << "Cannot book N-tuple2:" << long( m_tuple2 ) << endmsg;
      return StatusCode::FAILURE;
    }
  }
  return StatusCode::SUCCESS;
}
//*************************************************
StatusCode LumTau::execute() {
  StatusCode sc = StatusCode::SUCCESS;
 
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "in execute()" << endmsg;
 
  SmartDataPtr<Event::EventHeader> eventHeader( eventSvc(), "/Event/EventHeader" );
  SmartDataPtr<EvtRecEvent>        evtRecEvent( eventSvc(), EventModel::EvtRec::EvtRecEvent );
 
  log << MSG::DEBUG << "ncharg, nneu, tottks = " << evtRecEvent->totalCharged() << " , "
      << evtRecEvent->totalNeutral() << " , " << evtRecEvent->totalTracks() << endmsg;
  SmartDataPtr<EvtRecTrackCol> evtRecTrkCol( eventSvc(), EventModel::EvtRec::EvtRecTrackCol );
 
  log << MSG::DEBUG << "ncharg, nneu, tottks = " << evtRecEvent->totalCharged() << " , "
      << evtRecEvent->totalNeutral() << " , " << evtRecEvent->totalTracks() << endmsg;
  SmartDataPtr<RecEsTimeCol> evTimeCol( eventSvc(), "/Event/Recon/RecEsTimeCol" );
 
  int interval = 200;
  interval     = m_intervalSvc->getTInterval();
  log << MSG::DEBUG << "interval = " << interval << endmsg;
  m_topup = false;
  if ( interval <= 0 )
  {
    log << MSG::FATAL << "error in reading the interval!" << endmsg;
    return StatusCode::FAILURE;
  }
  else if ( interval < 100 )
  { // top-up mode
    m_topup = true;
  }
 
  double time = eventHeader->time();
  log << MSG::DEBUG << "time = " << time << endmsg;
  if ( time <= 0 ) return StatusCode::SUCCESS;
  m_time  = time;
  m_etsT1 = eventHeader->etsT1() / 2000000.; // convert to second
  log << MSG::DEBUG << "ets time = " << m_etsT1 << endmsg;
 
  m_run = eventHeader->runNumber();
  m_rec = eventHeader->eventNumber();
 
  if ( m_rec % 1000 == 0 )
    log << MSG::INFO << "Run   " << m_run << "     Event   " << m_rec << endmsg;
 
  double Emax1 = -1;
  double Emax2 = -1;
  int    Imax[2];
 
  if ( ( evtRecEvent->totalTracks() >= 2 ) )
  {
    double etot = 0.;
    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();
      double        Ener   = emcTrk->energy();
      if ( Ener > Emax2 )
      {
        Emax2   = Ener;
        Imax[1] = i;
      }
      if ( Ener > Emax1 )
      {
        Emax2   = Emax1;
        Imax[1] = Imax[0];
        Emax1   = Ener;
        Imax[0] = i;
      }
      etot += Ener;
    }
 
    m_etot = etot;
    m_e1   = Emax1;
    m_e2   = Emax2;
 
    log << MSG::INFO << "Emax1 = " << Emax1 << "Emax2= " << Emax2 << endmsg;
    if ( Emax1 > 0 && Emax2 > 0 )
    {
      double emcphi[2], emctht[2];
      for ( int i = 0; i < 2; i++ )
      {
        if ( i >= evtRecTrkCol->size() ) break;
        EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + Imax[i];
        if ( !( *itTrk )->isEmcShowerValid() ) continue;
        RecEmcShower* emcTrk = ( *itTrk )->emcShower();
        emcphi[i]            = emcTrk->phi();
        emctht[i]            = emcTrk->theta();
      }
 
      double dltphi = ( fabs( emcphi[0] - emcphi[1] ) - pai ) * 180. / pai;
      double dlttht = ( fabs( emctht[0] + emctht[1] ) - pai ) * 180. / pai;
      m_costht1     = cos( emctht[0] );
      m_costht2     = cos( emctht[1] );
      m_phi1        = emcphi[0];
      m_phi2        = emcphi[1];
      m_dlttht      = dlttht;
      m_dltphi      = dltphi;
    }
    else
    {
      m_etot    = -1;
      m_e1      = -1;
      m_e2      = -1;
      m_costht1 = -1;
      m_costht2 = -1;
      m_phi1    = -1;
      m_phi2    = -1;
      m_dlttht  = -1;
      m_dltphi  = -1;
    }
  }
  else
  {
    m_etot    = -1;
    m_e1      = -1;
    m_e2      = -1;
    m_costht1 = -1;
    m_costht2 = -1;
    m_phi1    = -1;
    m_phi2    = -1;
    m_dlttht  = -1;
    m_dltphi  = -1;
  }
 
  //    m_tuple2->write();
 
  StatusCode DiskWrite = m_tuple2->write();
  if ( DiskWrite != StatusCode::SUCCESS )
  {
    log << MSG::FATAL << "ERROR In LumTau DiskWrite!" << endmsg;
    exit( 1 );
  }
 
  return StatusCode::SUCCESS;
}
 
StatusCode LumTau::finalize() {
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "Event Finalize" << endmsg;
  return StatusCode::SUCCESS;
}