BhabhaPreSelect.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 "BhabhaPreSelect/BhabhaPreSelect.h"
#include "BhabhaPreSelect/BhabhaType.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
int BhabhaPreSelect::idmax[43] = { 40,  44,  48,  56,  64,  72,  80,  80,  76,  76,  88,
                                   88,  100, 100, 112, 112, 128, 128, 140, 140, 160, 160,
                                   160, 160, 176, 176, 176, 176, 208, 208, 208, 208, 240,
                                   240, 240, 240, 256, 256, 256, 256, 288, 288, 288 };
 
/////////////////////////////////////////////////////////////////////////////
DECLARE_COMPONENT( BhabhaPreSelect )
BhabhaPreSelect::BhabhaPreSelect( const std::string& name, ISvcLocator* pSvcLocator )
    : Algorithm( name, pSvcLocator ) {
 
  // Declare the properties
 
  declareProperty( "BarrelOrEndcap", m_BarrelOrEndcap = 1 );
  declareProperty( "Output", m_output = false );
  // cout<<"  BarrelOrEndcap  "<<m_BarrelOrEndcap<<endl;
}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
StatusCode BhabhaPreSelect::initialize() {
  MsgStream log( msgSvc(), name() );
 
  log << MSG::INFO << "in initialize()" << endmsg;
  cout << "  BarrelOrEndcap  " << m_BarrelOrEndcap << endl;
 
  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( "etot", m_etot );
        status = m_tuple1->addItem( "mdc_hit1", m_mdc_hit1 );
        status = m_tuple1->addItem( "mdc_hit2", m_mdc_hit2 );
      }
      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 BhabhaPreSelect::execute() {
 
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "in execute()" << endmsg;
 
  SmartDataPtr<Event::EventHeader> eventHeader( eventSvc(), "/Event/EventHeader" );
  // if(!eventHeader)
  //{
  //   cout<<"  eventHeader  "<<endl;
  //   return StatusCode::FAILURE;
  // }
 
  int run   = eventHeader->runNumber();
  int event = eventHeader->eventNumber();
  // cout<<"  event  "<<event<<endl;
  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;
  }
 
  m_events++;
  setFilterPassed( false );
 
  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;
 
  int total1 = 0; // mdc hit
  int total2 = 0; // mdc hit
  if ( GoodShowers.size() >= 2 && MaxE > 1.0 && abs( dthe ) < 3 &&
       ( ( dphi > -25 && dphi < -4 ) || ( dphi > 2 && dphi < 20 ) ) && etot > 2.7 )
  {
 
    double phi1 = MaxPhi < 0 ? MaxPhi + CLHEP::twopi : MaxPhi;
    double phi2 = SecPhi < 0 ? SecPhi + CLHEP::twopi : SecPhi;
 
    // Define sector (phi11,phi12) and (phi21,phi22)
    double phi11 = min( phi1, phi2 );
    double phi22 = max( phi1, phi2 );
    double phi12 = ( phi11 + phi22 - CLHEP::pi ) * 0.5;
    double phi21 = ( phi11 + phi22 + CLHEP::pi ) * 0.5;
    if ( phi12 < 0. ) phi12 += CLHEP::twopi;
    if ( phi21 > CLHEP::twopi ) phi21 -= CLHEP::twopi;
 
    SmartDataPtr<MdcDigiCol> mdcDigiCol( evtSvc(), "/Event/Digi/MdcDigiCol" );
    if ( !mdcDigiCol )
    {
      log << MSG::FATAL << "Could not find MdcDigiCol!" << endmsg;
      return StatusCode::FAILURE;
    }
    int hitnum = mdcDigiCol->size();
    for ( int i = 0; i < hitnum; i++ )
    {
      MdcDigi* digi   = dynamic_cast<MdcDigi*>( mdcDigiCol->containedObject( i ) );
      double   time   = digi->getTimeChannel();
      double   charge = digi->getChargeChannel();
      if ( time == 0x7FFFFFFF || charge == 0x7FFFFFFF ) continue;
      Identifier   id( digi->identify() );
      unsigned int iphi   = MdcID::wire( id );
      unsigned int ilayer = MdcID::layer( id );
      if ( ilayer >= 43 )
        log << MSG::ERROR << "MDC(" << ilayer << "," << iphi << ")" << endmsg;
      double phi = CLHEP::twopi * iphi / idmax[ilayer];
      if ( WhetherSector( phi, phi11, phi12 ) ) total1++;
      else if ( WhetherSector( phi, phi21, phi22 ) ) total2++;
      // cout<<"phi="<<phi<<",phi11="<<phi11<<",phi12="<<phi12
      //<<",phi21="<<phi21<<",phi22="<<phi22<<endl;
    }
    // cout<<"total1="<<total1<<"total2="<<total2<<endl;
    if ( ( m_BarrelOrEndcap == 1 && total1 > 15 && total2 > 15 ) ||
         ( m_BarrelOrEndcap != 1 && total1 > 5 && total2 > 5 ) )
    {
      setFilterPassed( true );
      selected++;
    }
  }
 
  if ( m_output )
  {
    m_etot      = etot;
    m_mdc_hit1  = total1;
    m_mdc_hit2  = total2;
    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 BhabhaPreSelect::finalize() {
 
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "in finalize()" << endmsg;
  cout << "m_events " << m_events << endl;
  cout << " selected " << selected << endl;
  cout << "m_rejected " << m_rejected << endl;
  cout << "m_oneProngsSelected " << m_oneProngsSelected << endl;
  cout << "m_twoProngsMatchedSelected " << m_twoProngsMatchedSelected << endl;
  cout << "m_twoProngsOneMatchedSelected " << m_twoProngsOneMatchedSelected << endl;
 
  return StatusCode::SUCCESS;
}
 
bool BhabhaPreSelect::WhetherSector( double ph, double ph1, double ph2 ) {
  double phi1  = min( ph1, ph2 );
  double phi2  = max( ph1, ph2 );
  double delta = 0.0610865; // 3.5*3.1415926/180.
  if ( ( phi2 - phi1 ) < CLHEP::pi )
  {
    phi1 -= delta;
    phi2 += delta;
    if ( phi1 < 0. ) phi1 += CLHEP::twopi;
    if ( phi2 > CLHEP::twopi ) phi2 -= CLHEP::twopi;
    double tmp1 = min( phi1, phi2 );
    double tmp2 = max( phi1, phi2 );
    phi1        = tmp1;
    phi2        = tmp2;
  }
  else
  {
    phi1 += delta;
    phi2 -= delta;
  }
  if ( ( phi2 - phi1 ) < CLHEP::pi )
  {
    if ( ph <= phi2 && ph >= phi1 ) return true;
    else return false;
  }
  else
  {
    if ( ph >= phi2 || ph <= phi1 ) return true;
    else return false;
  }
}