MdcCalRecTrk.cxx

Go to the documentation of this file.

#include "MdcCalibAlg/MdcCalRecTrk.h"
 
#include "CLHEP/Units/PhysicalConstants.h"
#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/IDataProviderSvc.h"
#include "GaudiKernel/IMessageSvc.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/SmartDataPtr.h"
#include "GaudiKernel/StatusCode.h"
#include "McTruth/MdcMcHit.h"
#include "MdcRawEvent/MdcDigi.h"
 
MdcCalRecTrk::MdcCalRecTrk( int pid ) { m_pid = pid; }
 
MdcCalRecTrk::~MdcCalRecTrk() {
  unsigned int i;
  for ( i = 0; i < m_rechit.size(); i++ ) { delete m_rechit[i]; }
  m_rechit.clear();
}
 
void MdcCalRecTrk::setRecTrk( RecMdcTrackCol::iterator it_trk ) {
  IMessageSvc* msgSvc;
  Gaudi::svcLocator()->service( "MessageSvc", msgSvc );
  MsgStream log( msgSvc, "MdcCalRecTrk" );
  log << MSG::DEBUG << "MdcCalRecTrk::setRecTrk()" << endmsg;
 
  m_dr           = ( *it_trk )->helix( 0 );
  m_phi0         = ( *it_trk )->helix( 1 );
  m_kappa        = ( *it_trk )->helix( 2 );
  m_dz           = ( *it_trk )->helix( 3 );
  m_tanl         = ( *it_trk )->helix( 4 );
  m_chisq        = ( *it_trk )->chi2();
  m_nhits        = ( *it_trk )->getNhits();
  m_helix        = ( *it_trk )->helix();
  m_helixerr     = ( *it_trk )->err();
  m_fgNoiseRatio = fgNoiseRatio( m_phi0 );
 
  double mass;
  if ( 0 == m_pid ) mass = 0.000511;
  else if ( 1 == m_pid ) mass = 0.105658;
  else if ( 2 == m_pid ) mass = 0.139570;
  else if ( 3 == m_pid ) mass = 0.493677;
  else if ( 4 == m_pid ) mass = 0.938272;
  else if ( 5 == m_pid ) mass = 0.139570; // cosmic-ray
  else log << MSG::FATAL << "wrong particle type" << endmsg;
  m_p4 = ( *it_trk )->p4( mass );
 
  m_dr *= 10.0; // cm -> mm
  m_dz *= 10.0; // cm -> mm
 
  m_pt = 1.0 / m_kappa;
  m_p  = m_pt * sqrt( m_tanl * m_tanl + 1.0 );
 
  HitRefVec           gothits = ( *it_trk )->getVecHits();
  HitRefVec::iterator it_hit  = gothits.begin();
  MdcCalRecHit*       rechit;
  for ( ; it_hit != gothits.end(); it_hit++ )
  {
    rechit = new MdcCalRecHit();
    rechit->setRecHit( it_hit );
    m_rechit.push_back( rechit );
  }
}
 
void MdcCalRecTrk::setKalTrk( RecMdcKalTrackCol::iterator it_trk ) {
  IMessageSvc* msgSvc;
  Gaudi::svcLocator()->service( "MessageSvc", msgSvc );
  MsgStream log( msgSvc, "MdcCalRecTrk" );
  log << MSG::DEBUG << "MdcCalRecTrk::setKalTrk()" << endmsg;
 
  if ( 0 == m_pid ) RecMdcKalTrack::setPidType( RecMdcKalTrack::electron );
  else if ( 1 == m_pid ) RecMdcKalTrack::setPidType( RecMdcKalTrack::muon );
  else if ( 2 == m_pid ) RecMdcKalTrack::setPidType( RecMdcKalTrack::pion );
  else if ( 3 == m_pid ) RecMdcKalTrack::setPidType( RecMdcKalTrack::kaon );
  else if ( 4 == m_pid ) RecMdcKalTrack::setPidType( RecMdcKalTrack::proton );
  else if ( 5 == m_pid ) RecMdcKalTrack::setPidType( RecMdcKalTrack::muon ); // cosmic-ray
  else log << MSG::FATAL << "wrong particle type" << endmsg;
 
  m_dr           = ( *it_trk )->dr();
  m_phi0         = ( *it_trk )->fi0();
  m_kappa        = ( *it_trk )->kappa();
  m_dz           = ( *it_trk )->dz();
  m_tanl         = ( *it_trk )->tanl();
  m_chisq        = ( *it_trk )->chi2();
  m_fgNoiseRatio = fgNoiseRatio( m_phi0 );
 
  m_p4.setPx( -sin( m_phi0 ) / fabs( m_kappa ) );
  m_p4.setPy( cos( m_phi0 ) / fabs( m_kappa ) );
  m_p4.setPz( m_tanl / fabs( m_kappa ) );
  double p3 = m_p4.mag();
  double mass;
  if ( 0 == m_pid ) mass = 0.000511;
  else if ( 1 == m_pid ) mass = 0.105658;
  else if ( 2 == m_pid ) mass = 0.139570;
  else if ( 3 == m_pid ) mass = 0.493677;
  else if ( 4 == m_pid ) mass = 0.938272;
  else if ( 5 == m_pid ) mass = 0.139570; // cosmic-ray
  m_p4.setE( sqrt( p3 * p3 + mass * mass ) );
 
  m_dr *= 10.0; // cm -> mm
  m_dz *= 10.0; // cm -> mm
 
  m_pt = 1.0 / m_kappa;
  m_p  = m_pt * sqrt( m_tanl * m_tanl + 1.0 );
 
  //      cout << setw(15) << m_p << setw(15) << m_dr << setw(15) << m_phi0
  //      << setw(15) << m_kappa << setw(15) << m_dz << setw(15) << m_tanl << endl;
 
  HelixSegRefVec           gothelixsegs = ( *it_trk )->getVecHelixSegs();
  HelixSegRefVec::iterator it_hit       = gothelixsegs.begin();
  MdcCalRecHit*            rechit;
 
  int k = 0;
  for ( ; it_hit != gothelixsegs.end(); it_hit++ )
  {
    rechit = new MdcCalRecHit();
    rechit->setKalHit( it_hit );
    m_rechit.push_back( rechit );
 
    k++;
  }
  m_nhits = k;
  //      cout<<"hits "<<m_nhits <<endl;
}
 
bool MdcCalRecTrk::fgNoiseRatio( double phi0 ) {
  IMessageSvc* msgSvc;
  Gaudi::svcLocator()->service( "MessageSvc", msgSvc );
  MsgStream log( msgSvc, "MdcCalib" );
 
  IDataProviderSvc* eventSvc = NULL;
  Gaudi::svcLocator()->service( "EventDataSvc", eventSvc );
  SmartDataPtr<MdcDigiCol> mdcDigiCol( eventSvc, "/Event/Digi/MdcDigiCol" );
  if ( !mdcDigiCol ) { log << MSG::FATAL << "Could not find event" << endmsg; }
 
  int  lay;
  int  cel;
  bool hitCel[43][288];
  for ( lay = 0; lay < 43; lay++ )
  {
    for ( cel = 0; cel < 288; cel++ ) { hitCel[lay][cel] = false; }
  }
 
  MdcDigiCol::iterator iter = mdcDigiCol->begin();
  unsigned             fgOverFlow;
  int                  digiId = 0;
  Identifier           id;
  for ( ; iter != mdcDigiCol->end(); iter++, digiId++ )
  {
    MdcDigi* aDigi = ( *iter );
    id             = ( aDigi )->identify();
    lay            = MdcID::layer( id );
    cel            = MdcID::wire( id );
    fgOverFlow     = ( aDigi )->getOverflow();
 
    bool goodTQ = true;
    if ( ( ( fgOverFlow & 3 ) != 0 ) || ( ( fgOverFlow & 12 ) != 0 ) ||
         ( aDigi->getTimeChannel() == 0x7FFFFFFF ) ||
         ( aDigi->getChargeChannel() == 0x7FFFFFFF ) )
      goodTQ = false;
 
    bool goodT = true;
    if ( ( ( fgOverFlow & 1 ) != 0 ) || ( aDigi->getTimeChannel() == 0x7FFFFFFF ) )
      goodT = false;
 
    //    if(!goodTQ) continue;
    if ( !goodT ) continue;
 
    hitCel[lay][cel] = true;
  }
 
  SmartDataPtr<RecMdcTrackCol> newtrkCol( eventSvc, "/Event/Recon/RecMdcTrackCol" );
  if ( !newtrkCol )
  {
    log << MSG::ERROR << "Could not find RecMdcTrackCol" << endmsg;
    return -1;
  }
 
  int                      nNoiRange = 4; // number of cells
  double                   dphi      = 1.0;
  RecMdcTrackCol::iterator it_trk    = newtrkCol->begin();
  for ( it_trk = newtrkCol->begin(); it_trk != newtrkCol->end(); it_trk++ )
  {
    HitRefVec           gothits  = ( *it_trk )->getVecHits();
    HitRefVec::iterator it_hit   = gothits.begin();
    HepVector           helix    = ( *it_trk )->helix();
    HepSymMatrix        helixErr = ( *it_trk )->err();
    double              phi0Rec  = ( *it_trk )->helix( 1 );
    double              delphi;
    if ( phi0Rec > phi0 ) delphi = phi0Rec - phi0;
    else delphi = phi0 - phi0Rec;
    if ( delphi > CLHEP::pi ) delphi -= CLHEP::pi;
    if ( delphi > ( CLHEP::pi * 0.17 ) ) continue;
 
    int nhitLay1 = 0; // lay0-3
    int nhitLay2 = 0; // lay4-7
    int nhitLay3 = 0; // lay8-19
    int nhitLay4 = 0; // lay20-42
    int nhitT1   = 0;
    int nhitT2   = 0;
    int nhitT3   = 0;
    int nhitT4   = 0;
    for ( lay = 0; lay < 8; lay++ )
    {
      double docamin = 0.7; // cm
      if ( lay > 7 ) docamin = 0.9;
      int celmin = -1;
      int ncel   = m_mdcGeomSvc->Layer( lay )->NCell();
      for ( cel = 0; cel < ncel; cel++ )
      {
        double            wphi;
        const MdcGeoWire* pWire = m_mdcGeomSvc->Wire( lay, cel );
        double            xx    = pWire->Backward().x();
        double            yy    = pWire->Backward().y();
        double            rr    = sqrt( ( xx * xx ) + ( yy * yy ) );
        if ( yy >= 0 ) wphi = acos( xx / rr );
        else wphi = CLHEP::twopi - acos( xx / rr );
 
        double phiTrk = phi0 + CLHEP::halfpi;
        if ( phiTrk > CLHEP::twopi ) phiTrk -= CLHEP::twopi;
        if ( !( ( fabs( wphi - phiTrk ) < dphi ) ||
                ( fabs( wphi + CLHEP::twopi - phiTrk ) < dphi ) ||
                ( fabs( wphi - CLHEP::twopi - phiTrk ) < dphi ) ) )
        { continue; }
 
        double doca = m_mdcUtilitySvc->doca( lay, cel, helix, helixErr );
        //          cout << setw(5) << lay << setw(5) << cel << setw(15) << doca <<
        // endl;
        if ( fabs( doca ) < fabs( docamin ) )
        {
          docamin = doca;
          celmin  = cel;
        }
      }
      if ( celmin > -1 )
      {
        if ( lay < 4 ) nhitLay1++;
        else if ( lay < 8 ) nhitLay2++;
        else if ( lay < 20 ) nhitLay3++;
        else nhitLay4++;
        for ( int ii = ( -nNoiRange ); ii <= nNoiRange; ii++ )
        {
          if ( 0 == ii ) continue;
          int icell = celmin + ii;
          if ( icell >= ncel ) icell -= ncel;
          if ( icell < 0 ) icell += ncel;
          //             cout << "hit " << setw(5) << lay << setw(5) << celmin <<
          // setw(5) << icell << setw(5) << hitCel[lay][icell]<<endl;
          if ( hitCel[lay][icell] )
          {
            if ( lay < 4 ) nhitT1++;
            else if ( lay < 8 ) nhitT2++;
            else if ( lay < 20 ) nhitT3++;
            else nhitT4++;
          }
        }
      }
    }
    if ( ( nhitLay1 <= 0 ) || ( nhitLay2 <= 0 ) || ( nhitLay3 <= 0 ) || ( nhitLay4 <= 0 ) )
      return false;
    double ratio1 = (double)nhitT1 / ( (double)nhitLay1 * (double)nNoiRange * 2.0 );
    double ratio2 = (double)nhitT2 / ( (double)nhitLay2 * (double)nNoiRange * 2.0 );
    double ratio3 = (double)nhitT3 / ( (double)nhitLay3 * (double)nNoiRange * 2.0 );
    double ratio4 = (double)nhitT4 / ( (double)nhitLay4 * (double)nNoiRange * 2.0 );
 
    if ( ( ratio1 > 0.08 ) || ( ratio2 > 0.08 ) || ( ratio3 > 0.03 ) || ( ratio4 > 0.03 ) )
      return false;
    else return true;
  }
  return false;
}