LocalLambdaSelector.cxx

Go to the documentation of this file.

#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/IJobOptionsSvc.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/PropertyMgr.h"
 
#include "EvtRecEvent/EvtRecVeeVertex.h"
#include "DTagAlg/LocalLambdaSelector.h"
#include "VertexFit/VertexFit.h"
#include "VertexFitRefine/VertexFitRefine.h"
#include "VertexFit/SecondVertexFit.h"
#include "VertexFit/IVertexDbSvc.h"
#include "ParticleID/ParticleID.h"
 
 
LocalLambdaSelector::LocalLambdaSelector()
{
        IJobOptionsSvc* jobSvc;
        Gaudi::svcLocator()->service("JobOptionsSvc", jobSvc);
 
        PropertyMgr m_propMgr;
 
        m_propMgr.declareProperty("LambdaMinMassCut",  m_minMass = 1.10 );
        m_propMgr.declareProperty("LambdaMaxMassCut",  m_maxMass = 1.13 );
        m_propMgr.declareProperty("LambdaMaxChisq",    m_maxChisq = 100 );
        m_propMgr.declareProperty("UseProtonPID",      m_UseProtonPID = true );
        m_propMgr.declareProperty("ProtonPIDChiCut",   m_ChiCut  = 8);
 
        m_propMgr.declareProperty("DoSecondaryVFit",   m_doSecondaryVFit = false );
        m_propMgr.declareProperty("LambdaMaxVFitChisq",    m_maxVFitChisq = 100 );
        m_propMgr.declareProperty("LambdaMinFlightSig",    m_minFlightSig = 2.0 );
      m_propMgr.declareProperty("UseVFRefine",       m_useVFrefine = true);
 
      m_propMgr.declareProperty("UseBFieldCorr",     m_useBFC = true);
 
        jobSvc->setMyProperties("LocalLambdaSelector", &m_propMgr);
}
 
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;
        Gaudi::svcLocator()->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;
 
            
        }
 
 
 
 
 
 
}
 
LocalLambdaSelector lambdaSelector;