DigammaPreSelect.cxx

Go to the documentation of this file.

#include "GaudiKernel/AlgFactory.h"
#include "GaudiKernel/IDataProviderSvc.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/PropertyMgr.h"
#include "GaudiKernel/SmartDataPtr.h"
 
#include "EventModel/Event.h"
#include "EventModel/EventHeader.h"
#include "EventModel/EventModel.h"
#include "EvtRecEvent/EvtRecEvent.h"
#include "EvtRecEvent/EvtRecTrack.h"
 
#include "CLHEP/Vector/LorentzVector.h"
#include "CLHEP/Vector/ThreeVector.h"
#include "CLHEP/Vector/TwoVector.h"
#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/IHistogramSvc.h"
#include "GaudiKernel/INTupleSvc.h"
#include "GaudiKernel/NTuple.h"
#include "MdcRawEvent/MdcDigi.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 "DigammaPreSelect/DigammaPreSelect.h"
int selected = 0;
 
#include <vector>
 
typedef std::vector<int>              Vint;
typedef std::vector<HepLorentzVector> Vp4;
const double                          PI    = 3.1415927;
const double                          EBeam = 1.548; // temporary
/////////////////////////////////////////////////////////////////////////////
DECLARE_COMPONENT( DigammaPreSelect )
DigammaPreSelect::DigammaPreSelect( const std::string& name, ISvcLocator* pSvcLocator )
    : Algorithm( name, pSvcLocator ) {
 
  // Declare the properties
 
  declareProperty( "Vr0cut", m_vr0cut = 1.0 ); // suggest cut: |z0|<5cm && r0<1cm
  declareProperty( "Vz0cut", m_vz0cut = 5.0 );
 
  declareProperty( "lowEnergyShowerCut", m_lowEnergyShowerCut = 0.1 );
  declareProperty( "highEnergyShowerCut", m_highEnergyShowerCut = 0.5 );
  declareProperty( "matchThetaCut", m_matchThetaCut = 0.2 );
  declareProperty( "matchPhiCut", m_matchPhiCut = 0.2 );
 
  declareProperty( "highMomentumCut", m_highMomentumCut = 0.5 );
  declareProperty( "EoPMaxCut", m_EoPMaxCut = 1.3 );
  declareProperty( "EoPMinCut", m_EoPMinCut = 0.6 );
  declareProperty( "minAngShEnergyCut", m_minAngShEnergyCut = 0.2 );
  declareProperty( "minAngCut", m_minAngCut = 0.3 );
  declareProperty( "acolliCut", m_acolliCut = 0.03 );
  declareProperty( "eNormCut", m_eNormCut = 0.5 );
  declareProperty( "pNormCut", m_pNormCut = 0.5 );
  declareProperty( "BarrelOrEndcap", m_BarrelOrEndcap = 1 );
  declareProperty( "Output", m_output = false );
  declareProperty( "oneProngMomentumCut", m_oneProngMomentumCut = 1.2 );
}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
StatusCode DigammaPreSelect::initialize() {
  MsgStream log( msgSvc(), name() );
 
  log << MSG::INFO << "in initialize()" << endmsg;
 
  if ( m_output )
  {
    StatusCode status;
    NTuplePtr  nt1( ntupleSvc(), "FILE1/bhabha" );
    if ( nt1 ) m_tuple1 = nt1;
    else
    {
      m_tuple1 = ntupleSvc()->book( "FILE1/bhabha", CLID_ColumnWiseTuple, "N-Tuple example" );
      if ( m_tuple1 )
      {
        status = m_tuple1->addItem( "sh1_ene", m_sh1_ene );
        status = m_tuple1->addItem( "sh1_theta", m_sh1_theta );
        status = m_tuple1->addItem( "sh1_phi", m_sh1_phi );
        status = m_tuple1->addItem( "sh2_ene", m_sh2_ene );
        status = m_tuple1->addItem( "sh2_theta", m_sh2_theta );
        status = m_tuple1->addItem( "sh2_phi", m_sh2_phi );
        status = m_tuple1->addItem( "di_phi", m_di_phi );
        status = m_tuple1->addItem( "di_the", m_di_the );
        status = m_tuple1->addItem( "acolli", m_acolli );
        status = m_tuple1->addItem( "mdc_hit", m_mdc_hit );
        status = m_tuple1->addItem( "etot", m_etot );
      }
      else
      {
        log << MSG::ERROR << "    Cannot book N-tuple:" << long( m_tuple1 ) << endmsg;
        return StatusCode::FAILURE;
      }
    }
 
    NTuplePtr nt2( ntupleSvc(), "FILE1/bha1" );
    if ( nt2 ) m_tuple2 = nt2;
    else
    {
      m_tuple2 = ntupleSvc()->book( "FILE1/bha1", CLID_ColumnWiseTuple, "N-Tuple example" );
      if ( m_tuple2 )
      {
        status = m_tuple2->addItem( "sh_ene", m_sh_ene );
        status = m_tuple2->addItem( "sh_theta", m_sh_theta );
        status = m_tuple2->addItem( "sh_phi", m_sh_phi );
      }
      else
      {
        log << MSG::ERROR << "    Cannot book N-tuple:" << long( m_tuple2 ) << endmsg;
        return StatusCode::FAILURE;
      }
    }
  }
 
  //
  //--------end of book--------
  //
 
  m_rejected                    = 0;
  m_events                      = 0;
  m_oneProngsSelected           = 0;
  m_twoProngsMatchedSelected    = 0;
  m_twoProngsOneMatchedSelected = 0;
  m_selectedTrkID1              = 999;
  m_selectedTrkID2              = 999;
 
  log << MSG::INFO << "successfully return from initialize()" << endmsg;
  return StatusCode::SUCCESS;
}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
StatusCode DigammaPreSelect::execute() {
 
  setFilterPassed( false );
 
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "in execute()" << endmsg;
 
  m_events++;
 
  SmartDataPtr<Event::EventHeader> eventHeader( eventSvc(), "/Event/EventHeader" );
  if ( !eventHeader )
  {
    cout << "  eventHeader  " << endl;
    return StatusCode::FAILURE;
  }
 
  int run   = eventHeader->runNumber();
  int event = eventHeader->eventNumber();
  if ( event % 1000 == 0 ) cout << "  run,event:  " << run << "," << event << endl;
 
  SmartDataPtr<EvtRecEvent> evtRecEvent( eventSvc(), EventModel::EvtRec::EvtRecEvent );
  if ( !evtRecEvent )
  {
    cout << "  evtRecEvent  " << endl;
    return StatusCode::FAILURE;
  }
 
  log << MSG::DEBUG << "ncharg, nneu, tottks = " << evtRecEvent->totalCharged() << " , "
      << evtRecEvent->totalNeutral() << " , " << evtRecEvent->totalTracks() << endmsg;
  SmartDataPtr<EvtRecTrackCol> evtRecTrkCol( eventSvc(), EventModel::EvtRec::EvtRecTrackCol );
  if ( !evtRecTrkCol )
  {
    cout << "  evtRecTrkCol  " << endl;
    return StatusCode::FAILURE;
  }
 
  Vint iGood;
  iGood.clear();
 
  double   ene5x5, theta, phi, eseed;
  double   showerX, showerY, showerZ;
  long int thetaIndex, phiIndex;
 
  RecEmcID     showerId;
  unsigned int npart;
 
  Vint ipip, ipim;
  ipip.clear();
  ipim.clear();
  Vp4 ppip, ppim;
  ppip.clear();
  ppim.clear();
 
  vector<RecEmcShower*> GoodShowers;
  GoodShowers.clear();
  double etot = 0;
  for ( int i = 0; i < evtRecEvent->totalTracks(); i++ )
  {
    if ( i >= evtRecTrkCol->size() ) break;
    EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + i;
    if ( ( *itTrk )->isEmcShowerValid() )
    {
      RecEmcShower* theShower = ( *itTrk )->emcShower();
      etot += theShower->e5x5();
      showerId   = theShower->getShowerId();
      npart      = EmcID::barrel_ec( showerId );
      ene5x5     = theShower->e5x5();
      thetaIndex = EmcID::theta_module( showerId );
      phiIndex   = EmcID::phi_module( showerId );
      if ( ene5x5 > 0.4 && ene5x5 < 4 && npart == 1 && ( m_BarrelOrEndcap == 1 ) )
      {
        GoodShowers.push_back( theShower );
        iGood.push_back( ( *itTrk )->trackId() );
      }
      else if ( ene5x5 > 0.4 && ene5x5 < 4 && ( npart == 2 || npart == 0 ) &&
                ( m_BarrelOrEndcap == 2 ) )
      {
        GoodShowers.push_back( theShower );
        iGood.push_back( ( *itTrk )->trackId() );
      }
      else if ( ene5x5 > 0.4 && ene5x5 < 4 && ( m_BarrelOrEndcap == 0 ) )
      {
        GoodShowers.push_back( theShower );
        iGood.push_back( ( *itTrk )->trackId() );
      }
    }
    // good photon cut will be set here
  }
  // Finish Good Photon Selection
  //
 
  double MaxE( 0 ), MaxPhi, MaxThe;
  int    MaxId;
  double SecE( 0 ), SecPhi, SecThe;
  for ( int i = 0; i < GoodShowers.size(); i++ )
  {
    RecEmcShower* theShower = GoodShowers[i];
    double        eraw      = theShower->energy();
    if ( eraw > MaxE )
    {
      MaxId  = i;
      MaxE   = eraw;
      MaxPhi = theShower->phi();
      MaxThe = theShower->theta();
    }
  }
  for ( int i = 0; i < GoodShowers.size(); i++ )
  {
    RecEmcShower* theShower = GoodShowers[i];
    double        eraw      = theShower->energy();
    if ( i != MaxId && eraw > SecE )
    {
      SecE   = eraw;
      SecPhi = theShower->phi();
      SecThe = theShower->theta();
    }
  }
 
  double dphi = ( fabs( MaxPhi - SecPhi ) - PI ) * 180. / PI;
  double dthe = ( fabs( MaxThe + SecThe ) - PI ) * 180. / PI;
  if ( GoodShowers.size() >= 2 && SecE > 1.0 && dphi > -4 && dphi < 2 && abs( dthe ) < 3 &&
       etot > 2.7 )
  {
    setFilterPassed( true );
    selected++;
  }
 
  if ( m_output )
  {
    m_etot      = etot;
    m_sh1_ene   = MaxE;
    m_sh1_theta = MaxThe;
    m_sh1_phi   = MaxPhi;
    m_sh2_ene   = SecE;
    m_sh2_theta = SecThe;
    m_sh2_phi   = SecPhi;
    m_di_phi    = dphi;
    m_di_the    = dthe;
    m_tuple1->write();
  }
 
  return StatusCode::SUCCESS;
}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
StatusCode DigammaPreSelect::finalize() {
 
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "in finalize()" << endmsg;
  cout << "total events: " << m_events << endl;
  cout << "selected digamma: " << selected << endl;
 
  return StatusCode::SUCCESS;
}