LocalLambdaSelector.cxx

Go to the documentation of this file.

#include "GaudiKernel/Bootstrap.h"
 
#include "EvtRecEvent/EvtRecVeeVertex.h"
#include "LocalLambdaSelector.h"
#include "ParticleID/ParticleID.h"
#include "VertexDbSvc/IVertexDbSvc.h"
#include "VertexFit/SecondVertexFit.h"
#include "VertexFit/VertexFit.h"
#include "VertexFitRefine/VertexFitRefine.h"
 
DECLARE_COMPONENT( LocalLambdaSelector )
 
LocalLambdaSelector::LocalLambdaSelector( const std::string& type, const std::string& name,
                                          const IInterface* parent )
    : base_class( type, name, parent ) {
  declareProperty( "LambdaMinMassCut", m_minMass = 1.10 );
  declareProperty( "LambdaMaxMassCut", m_maxMass = 1.13 );
  declareProperty( "LambdaMaxChisq", m_maxChisq = 100 );
  declareProperty( "UseProtonPID", m_UseProtonPID = true );
  declareProperty( "ProtonPIDChiCut", m_ChiCut = 8 );
 
  declareProperty( "DoSecondaryVFit", m_doSecondaryVFit = false );
  declareProperty( "LambdaMaxVFitChisq", m_maxVFitChisq = 100 );
  declareProperty( "LambdaMinFlightSig", m_minFlightSig = 2.0 );
  declareProperty( "UseVFRefine", m_useVFrefine = true );
 
  declareProperty( "UseBFieldCorr", m_useBFC = true );
}
 
bool LocalLambdaSelector::operator()( CDLambda& aLambda ) {
 
  aLambda.setUserTag( 1 );
  EvtRecVeeVertex* lambda = const_cast<EvtRecVeeVertex*>( aLambda.navLambda() );
 
  if ( fabs( lambda->vertexId() ) != 3122 ) return false;
 
  if ( !m_useBFC )
  {
    double mass = lambda->mass();
    if ( ( mass <= m_minMass ) || ( mass >= m_maxMass ) ) return false;
    if ( lambda->chi2() >= m_maxChisq ) return false;
  }
 
  // PID
  int index[2] = { 0, 1 };
  if ( lambda->vertexId() < 0 )
  {
    index[0] = 1;
    index[1] = 0;
  }
 
  EvtRecTrack* recTrk_proton = lambda->daughter( index[0] );
  EvtRecTrack* recTrk_pion   = lambda->daughter( index[1] );
 
  if ( m_UseProtonPID )
  {
    ParticleID* pid = ParticleID::instance();
    pid->init();
    pid->setMethod( pid->methodProbability() );
    pid->setChiMinCut( m_ChiCut );
    pid->setRecTrack( recTrk_proton );
    // pid->usePidSys(pid->useDedx() | pid->useTof1() | pid->useTof2() | pid->useTof());
    pid->usePidSys( pid->useDedx() | pid->useTofCorr() );
    pid->identify( pid->onlyPion() | pid->onlyKaon() | pid->onlyProton() );
    pid->calculate();
    if ( !( pid->probProton() > 0. && pid->probProton() > pid->probPion() &&
            pid->probProton() > pid->probKaon() ) )
      return false;
  }
 
  if ( !m_useBFC )
  {
    if ( !m_doSecondaryVFit ) return true;
  }
  // --------------------------------------------------
  // Do a secondary vertex fit and check the flight significance
  // --------------------------------------------------
 
  RecMdcKalTrack* mdcKalTrk_proton = recTrk_proton->mdcKalTrack();
  mdcKalTrk_proton->setPidType( RecMdcKalTrack::proton );
  WTrackParameter WTrk_proton;
  WTrk_proton = WTrackParameter( 0.9314940054, mdcKalTrk_proton->getZHelixP(),
                                 mdcKalTrk_proton->getZErrorP() );
 
  RecMdcKalTrack* mdcKalTrk_pion = recTrk_pion->mdcKalTrack();
  mdcKalTrk_pion->setPidType( RecMdcKalTrack::pion );
  WTrackParameter WTrk_pion =
      WTrackParameter( 0.13957018, mdcKalTrk_pion->getZHelix(), mdcKalTrk_pion->getZError() );
 
  // VertexParameter wideVertex;
  HepPoint3D   vWideVertex( 0., 0., 0. );
  HepSymMatrix evWideVertex( 3, 0 );
 
  evWideVertex[0][0] = 1.0e12;
  evWideVertex[1][1] = 1.0e12;
  evWideVertex[2][2] = 1.0e12;
 
  // add the primary vertex
  HepPoint3D   vBeamSpot( 0., 0., 0. );
  HepSymMatrix evBeamSpot( 3, 0 );
 
  IVertexDbSvc* vtxsvc;
  StatusCode    sc = serviceLocator()->service( "VertexDbSvc", vtxsvc );
  if ( vtxsvc->isVertexValid() )
  {
    double* dbv = vtxsvc->PrimaryVertex();
    double* vv  = vtxsvc->SigmaPrimaryVertex();
    for ( unsigned int ivx = 0; ivx < 3; ivx++ )
    {
      vBeamSpot[ivx]       = dbv[ivx];
      evBeamSpot[ivx][ivx] = vv[ivx] * vv[ivx];
    }
  }
  else
  {
    cout << "LambdaSELECTOR ERROR:  Could not find a vertex from VertexDbSvc" << endl;
    return false;
  }
 
  if ( m_useVFrefine )
  {
 
    VertexParameter wideVertex;
 
    wideVertex.setVx( vWideVertex );
    wideVertex.setEvx( evWideVertex );
 
    // First, perform a vertex fit refine
    VertexFitRefine* vtxfitr = VertexFitRefine::instance();
    vtxfitr->init();
 
    vtxfitr->AddTrack( 0, mdcKalTrk_proton, RecMdcKalTrack::proton );
    vtxfitr->AddTrack( 1, mdcKalTrk_pion, RecMdcKalTrack::pion );
    vtxfitr->AddVertex( 0, wideVertex, 0, 1 );
 
    bool fitok = vtxfitr->Fit();
 
    HepLorentzVector ppr     = vtxfitr->pfit( 0 );
    HepLorentzVector ppi     = vtxfitr->pfit( 1 );
    HepLorentzVector plambda = ppr + ppi;
 
    double mass     = plambda.m();
    double chisqvtx = vtxfitr->chisq( 0 );
 
    if ( m_useBFC )
    {
      if ( ( mass <= m_minMass ) || ( mass >= m_maxMass ) ) return false;
      // if ( lambda->chi2() >= m_maxChisq ) return false;
      if ( chisqvtx >= m_maxChisq ) return false;
 
      if ( !m_doSecondaryVFit ) return true;
    }
 
    // Now perform the secondary vertex fit
    SecondVertexFit* svtxfit = SecondVertexFit::instance();
    svtxfit->init();
 
    // add the primary vertex
    VertexParameter beamSpot;
 
    beamSpot.setVx( vBeamSpot );
    beamSpot.setEvx( evBeamSpot );
 
    VertexParameter primaryVertex( beamSpot );
    svtxfit->setPrimaryVertex( primaryVertex );
 
    // add the secondary vertex
    svtxfit->setVpar( vtxfitr->vpar( 0 ) );
 
    // add the Ks or lambda
    svtxfit->AddTrack( 0, vtxfitr->wVirtualTrack( 0 ) );
 
    // do the second vertex fit
    // if the fit fails, the default values will not be changed
    if ( !svtxfit->Fit() ) return false;
 
    // save the new ks parameters
    double vfitlength = svtxfit->decayLength();
    double vfiterror  = svtxfit->decayLengthError();
    double vfitchi2   = svtxfit->chisq();
    double flightsig  = 0;
    if ( vfiterror != 0 ) flightsig = vfitlength / vfiterror;
 
    // if the ks does not meet the criteria, the information from
    // the secondary vertex fit will not be filled (default = -999)
    if ( vfitchi2 > m_maxVFitChisq ) return false;
    if ( flightsig < m_minFlightSig ) return false;
 
    return true;
  }
  else
  {
    VertexParameter wideVertex;
 
    wideVertex.setVx( vWideVertex );
    wideVertex.setEvx( evWideVertex );
 
    // First, perform a vertex fit
    VertexFit* vtxfit = VertexFit::instance();
    vtxfit->init();
 
    // add the daughters
    vtxfit->AddTrack( 0, WTrk_proton );
    vtxfit->AddTrack( 1, WTrk_pion );
    vtxfit->AddVertex( 0, wideVertex, 0, 1 );
 
    // do the fit
    vtxfit->Fit( 0 );
    vtxfit->Swim( 0 );
    vtxfit->BuildVirtualParticle( 0 );
 
    WTrackParameter wlambda = vtxfit->wVirtualTrack( 0 );
 
    HepLorentzVector plambda = wlambda.p();
 
    double mass     = plambda.m();
    double chisqvtx = vtxfit->chisq( 0 );
 
    if ( m_useBFC )
    {
      if ( ( mass <= m_minMass ) || ( mass >= m_maxMass ) ) return false;
      // if ( lambda->chi2() >= m_maxChisq ) return false;
      if ( chisqvtx >= m_maxChisq ) return false;
 
      if ( !m_doSecondaryVFit ) return true;
    }
 
    // Now perform the secondary vertex fit
    SecondVertexFit* svtxfit = SecondVertexFit::instance();
    svtxfit->init();
 
    // add the primary vertex
    VertexParameter beamSpot;
    beamSpot.setVx( vBeamSpot );
    beamSpot.setEvx( evBeamSpot );
 
    VertexParameter primaryVertex( beamSpot );
    svtxfit->setPrimaryVertex( primaryVertex );
 
    // add the secondary vertex
    svtxfit->setVpar( vtxfit->vpar( 0 ) );
 
    // add the Ks or lambda
    svtxfit->AddTrack( 0, vtxfit->wVirtualTrack( 0 ) );
 
    // do the second vertex fit
    // if the fit fails, the default values will not be changed
    if ( !svtxfit->Fit() ) return false;
 
    // save the new ks parameters
    double vfitlength = svtxfit->decayLength();
    double vfiterror  = svtxfit->decayLengthError();
    double vfitchi2   = svtxfit->chisq();
    double flightsig  = 0;
    if ( vfiterror != 0 ) flightsig = vfitlength / vfiterror;
 
    // if the ks does not meet the criteria, the information from
    // the secondary vertex fit will not be filled (default = -999)
    if ( vfitchi2 > m_maxVFitChisq ) return false;
    if ( flightsig < m_minFlightSig ) return false;
 
    return true;
  }
}