|
|
|
#include <GsfTrackProducerBase.h>
Public Member Functions | |
| GsfTrackProducerBase (bool trajectoryInEvent, bool split) | |
| Constructor. | |
| virtual void | putInEvt (edm::Event &, const Propagator *prop, const MeasurementTracker *measTk, std::auto_ptr< TrackingRecHitCollection > &, std::auto_ptr< reco::GsfTrackCollection > &, std::auto_ptr< reco::TrackExtraCollection > &, std::auto_ptr< reco::GsfTrackExtraCollection > &, std::auto_ptr< std::vector< Trajectory > > &, AlgoProductCollection &, const reco::BeamSpot &) |
| Put produced collections in the event. | |
Protected Member Functions | |
| void | fillMode (reco::GsfTrack &track, const TrajectoryStateOnSurface innertsos, const Propagator &gsfProp, const TransverseImpactPointExtrapolator &tipExtrapolator, TrajectoryStateClosestToBeamLineBuilder &tscblBuilder, const reco::BeamSpot &bs) const |
| void | fillStates (TrajectoryStateOnSurface tsos, std::vector< reco::GsfComponent5D > &states) const |
Private Member Functions | |
| bool | computeModeAtTM (const TrajectoryMeasurement &tm, reco::GsfTrackExtra::Point &position, reco::GsfTrackExtra::Vector &momentum, Measurement1D &deltaP) const |
| position, momentum and estimated deltaP at an intermediate measurement (true if successful) | |
| void | localParametersFromQpMode (const TrajectoryStateOnSurface tsos, AlgebraicVector5 ¶meters, AlgebraicSymMatrix55 &covariance) const |
| local parameters rescaled with q/p from mode | |
Private Attributes | |
| bool | useSplitting |
Produce Tracks from TrackCandidates
Definition at line 28 of file GsfTrackProducerBase.h.
| GsfTrackProducerBase::GsfTrackProducerBase | ( | bool | trajectoryInEvent, |
| bool | split | ||
| ) | [inline, explicit] |
Constructor.
Definition at line 32 of file GsfTrackProducerBase.h.
:
TrackProducerBase<reco::GsfTrack>(trajectoryInEvent),
useSplitting(split){}
| bool GsfTrackProducerBase::computeModeAtTM | ( | const TrajectoryMeasurement & | tm, |
| reco::GsfTrackExtra::Point & | position, | ||
| reco::GsfTrackExtra::Vector & | momentum, | ||
| Measurement1D & | deltaP | ||
| ) | const [private] |
position, momentum and estimated deltaP at an intermediate measurement (true if successful)
Definition at line 359 of file GsfTrackProducerBase.cc.
References TrajectoryMeasurement::backwardPredictedState(), CollinearFitAtTM2::deltaP(), TrajectoryMeasurement::forwardPredictedState(), TrajectoryStateOnSurface::globalPosition(), TrajectoryStateOnSurface::isValid(), TrajectoryStateOnSurface::localError(), TrajectoryStateOnSurface::localParameters(), localParametersFromQpMode(), LocalTrajectoryError::matrix(), pos, TrajectoryMeasurement::recHit(), query::result, TrajectoryMeasurement::updatedState(), LocalTrajectoryParameters::vector(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().
Referenced by putInEvt().
{
//
// states
//
TrajectoryStateOnSurface fwdState = tm.forwardPredictedState();
TrajectoryStateOnSurface bwdState = tm.backwardPredictedState();
TrajectoryStateOnSurface upState = tm.updatedState();
if ( !fwdState.isValid() || !bwdState.isValid() || !upState.isValid() ) {
return false;
}
//
// position from mean, momentum from mode (in cartesian coordinates)
// following PF code
//
GlobalPoint pos = upState.globalPosition();
position = reco::GsfTrackExtra::Point(pos.x(),pos.y(),pos.z());
MultiTrajectoryStateMode mts;
GlobalVector mom;
bool result = mts.momentumFromModeCartesian(upState,mom);
if ( !result ) {
// std::cout << "momentumFromModeCartesian failed" << std::endl;
return false;
}
momentum = reco::GsfTrackExtra::Vector(mom.x(),mom.y(),mom.z());
//
// calculation from deltaP from fit to forward & backward predictions
// (momentum from mode) and hit
//
// prepare input parameter vectors and covariance matrices
AlgebraicVector5 fwdPars = fwdState.localParameters().vector();
AlgebraicSymMatrix55 fwdCov = fwdState.localError().matrix();
localParametersFromQpMode(fwdState,fwdPars,fwdCov);
AlgebraicVector5 bwdPars = bwdState.localParameters().vector();
AlgebraicSymMatrix55 bwdCov = bwdState.localError().matrix();
localParametersFromQpMode(bwdState,bwdPars,bwdCov);
LocalPoint hitPos(0.,0.,0.);
LocalError hitErr(-1.,-1.,-1.);
if ( tm.recHit()->isValid() ) {
hitPos = tm.recHit()->localPosition();
hitErr = tm.recHit()->localPositionError();
}
CollinearFitAtTM2 collinearFit(fwdPars,fwdCov,bwdPars,bwdCov,hitPos,hitErr);
deltaP = collinearFit.deltaP();
return true;
}
| void GsfTrackProducerBase::fillMode | ( | reco::GsfTrack & | track, |
| const TrajectoryStateOnSurface | innertsos, | ||
| const Propagator & | gsfProp, | ||
| const TransverseImpactPointExtrapolator & | tipExtrapolator, | ||
| TrajectoryStateClosestToBeamLineBuilder & | tscblBuilder, | ||
| const reco::BeamSpot & | bs | ||
| ) | const [protected] |
Definition at line 262 of file GsfTrackProducerBase.cc.
References FreeTrajectoryState::charge(), FreeTrajectoryState::curvilinearError(), reco::GsfTrack::dimensionMode, reco::BeamSpot::dxdz(), reco::BeamSpot::dydz(), TransverseImpactPointExtrapolator::extrapolate(), CurvilinearTrajectoryError::matrix(), GaussianSumUtilities1D::mean(), SingleGaussianState1D::mean(), GaussianSumUtilities1D::mode(), GaussianSumUtilities1D::modeIsValid(), FreeTrajectoryState::momentum(), MultiGaussianStateTransform::multiState1D(), reco::BeamSpot::position(), Propagator::propagate(), reco::GsfTrack::setMode(), mathSSE::sqrt(), TrajectoryStateOnSurface::surface(), GaussianSumUtilities1D::variance(), SingleGaussianState1D::variance(), reco::TrackBase::vz(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().
Referenced by putInEvt().
{
// Get transverse impact parameter plane (from mean). This is a first approximation;
// the mode is then extrapolated to the
// final position closest to the beamline.
GlobalPoint bsPos(bs.position().x()+(track.vz()-bs.position().z())*bs.dxdz(),
bs.position().y()+(track.vz()-bs.position().z())*bs.dydz(),
track.vz());
TrajectoryStateOnSurface vtxTsos = tipExtrapolator.extrapolate(innertsos,bsPos);
if ( !vtxTsos.isValid() ) vtxTsos = innertsos;
// extrapolate mixture
vtxTsos = gsfProp.propagate(innertsos,vtxTsos.surface());
if ( !vtxTsos.isValid() ) return; // failed (GsfTrack keeps mode = mean)
// extract mode
// build perigee parameters (for covariance to be stored)
AlgebraicVector5 modeParameters;
AlgebraicSymMatrix55 modeCovariance;
// set parameters and variances for "mode" state (local parameters)
for ( unsigned int iv=0; iv<5; ++iv ) {
MultiGaussianState1D state1D = MultiGaussianStateTransform::multiState1D(vtxTsos,iv);
GaussianSumUtilities1D utils(state1D);
modeParameters(iv) = utils.mode().mean();
modeCovariance(iv,iv) = utils.mode().variance();
if ( !utils.modeIsValid() ) {
// if mode calculation fails: use mean
modeParameters(iv) = utils.mean();
modeCovariance(iv,iv) = utils.variance();
}
}
// complete covariance matrix
// approximation: use correlations from mean
const AlgebraicSymMatrix55& meanCovariance(vtxTsos.localError().matrix());
for ( unsigned int iv1=0; iv1<5; ++iv1 ) {
for ( unsigned int iv2=0; iv2<iv1; ++iv2 ) {
double cov12 = meanCovariance(iv1,iv2) *
sqrt(modeCovariance(iv1,iv1)/meanCovariance(iv1,iv1)*
modeCovariance(iv2,iv2)/meanCovariance(iv2,iv2));
modeCovariance(iv1,iv2) = modeCovariance(iv2,iv1) = cov12;
}
}
TrajectoryStateOnSurface modeTsos(LocalTrajectoryParameters(modeParameters,
vtxTsos.localParameters().pzSign()),
LocalTrajectoryError(modeCovariance),
vtxTsos.surface(),
vtxTsos.magneticField(),
vtxTsos.surfaceSide());
TrajectoryStateClosestToBeamLine tscbl = tscblBuilder(*modeTsos.freeState(),bs);
if ( !tscbl.isValid() ) return; // failed (GsfTrack keeps mode = mean)
//
// extract state at PCA and create momentum vector and covariance matrix
//
FreeTrajectoryState fts = tscbl.trackStateAtPCA();
GlobalVector tscblMom = fts.momentum();
reco::GsfTrack::Vector mom(tscblMom.x(),tscblMom.y(),tscblMom.z());
reco::GsfTrack::CovarianceMatrixMode cov;
const AlgebraicSymMatrix55& tscblCov = fts.curvilinearError().matrix();
for ( unsigned int iv1=0; iv1<reco::GsfTrack::dimensionMode; ++iv1 ) {
for ( unsigned int iv2=0; iv2<reco::GsfTrack::dimensionMode; ++iv2 ) {
cov(iv1,iv2) = tscblCov(iv1,iv2);
}
}
track.setMode(fts.charge(),mom,cov);
}
| void GsfTrackProducerBase::fillStates | ( | TrajectoryStateOnSurface | tsos, |
| std::vector< reco::GsfComponent5D > & | states | ||
| ) | const [protected] |
Definition at line 249 of file GsfTrackProducerBase.cc.
References TrajectoryStateOnSurface::components(), and i.
Referenced by putInEvt().
{
reco::GsfComponent5D::ParameterVector pLocS;
reco::GsfComponent5D::CovarianceMatrix cLocS;
std::vector<TrajectoryStateOnSurface> components(tsos.components());
for ( std::vector<TrajectoryStateOnSurface>::const_iterator i=components.begin();
i!=components.end(); ++i ) {
states.push_back(reco::GsfComponent5D(i->weight(),i->localParameters().vector(),i->localError().matrix()));
}
}
| void GsfTrackProducerBase::localParametersFromQpMode | ( | const TrajectoryStateOnSurface | tsos, |
| AlgebraicVector5 & | parameters, | ||
| AlgebraicSymMatrix55 & | covariance | ||
| ) | const [private] |
local parameters rescaled with q/p from mode
Definition at line 331 of file GsfTrackProducerBase.cc.
References i, TrajectoryStateOnSurface::localError(), TrajectoryStateOnSurface::localParameters(), LocalTrajectoryError::matrix(), SingleGaussianState1D::mean(), GaussianSumUtilities1D::mode(), GaussianSumUtilities1D::modeIsValid(), MultiGaussianStateTransform::multiState1D(), Parameters::parameters, mathSSE::sqrt(), SingleGaussianState1D::variance(), and LocalTrajectoryParameters::vector().
Referenced by computeModeAtTM().
{
//
// parameters and errors from combined state
//
parameters = tsos.localParameters().vector();
covariance = tsos.localError().matrix();
//
// mode for parameter 0 (q/p)
//
MultiGaussianState1D qpState(MultiGaussianStateTransform::multiState1D(tsos,0));
GaussianSumUtilities1D qpGS(qpState);
if ( !qpGS.modeIsValid() ) return;
double qp = qpGS.mode().mean();
double varQp = qpGS.mode().variance();
//
// replace q/p value and variance, rescale correlation terms
// (heuristic procedure - alternative would be mode in 5D ...)
//
double VarQpRatio = sqrt(varQp/covariance(0,0));
parameters(0) = qp;
covariance(0,0) = varQp;
for ( int i=1; i<5; ++i ) covariance(i,0) *= VarQpRatio;
}
| void GsfTrackProducerBase::putInEvt | ( | edm::Event & | evt, |
| const Propagator * | prop, | ||
| const MeasurementTracker * | measTk, | ||
| std::auto_ptr< TrackingRecHitCollection > & | selHits, | ||
| std::auto_ptr< reco::GsfTrackCollection > & | selTracks, | ||
| std::auto_ptr< reco::TrackExtraCollection > & | selTrackExtras, | ||
| std::auto_ptr< reco::GsfTrackExtraCollection > & | selGsfTrackExtras, | ||
| std::auto_ptr< std::vector< Trajectory > > & | selTrajectories, | ||
| AlgoProductCollection & | algoResults, | ||
| const reco::BeamSpot & | bs | ||
| ) | [virtual] |
Put produced collections in the event.
Definition at line 27 of file GsfTrackProducerBase.cc.
References reco::TrackExtraBase::add(), alongMomentum, AnalyticalPropagator_cfi::AnalyticalPropagator, anyDirection, TrajectoryStateOnSurface::components(), computeModeAtTM(), TrajectoryStateOnSurface::curvilinearError(), Trajectory::direction(), fillMode(), fillStates(), Trajectory::firstMeasurement(), edm::Event::getRefBeforePut(), TrajectoryStateOnSurface::globalParameters(), i, TrajectoryStateOnSurface::isValid(), edm::ESHandleBase::isValid(), j, Trajectory::lastMeasurement(), TrajectoryStateOnSurface::localParameters(), LogDebug, LogTrace, TrajectoryStateOnSurface::magneticField(), Trajectory::measurements(), GlobalTrajectoryParameters::momentum(), L1TEmulatorMonitor_cff::p, position, GlobalTrajectoryParameters::position(), edm::Event::put(), LocalTrajectoryParameters::pzSign(), TrajectoryMeasurement::recHit(), Trajectory::recHits(), TrackProducerBase< reco::GsfTrack >::rTracks_, Trajectory::seedRef(), reco::Track::setExtra(), reco::GsfTrack::setGsfExtra(), reco::TrackBase::setHitPattern(), TrackProducerBase< reco::GsfTrack >::setSecondHitPattern(), DetId::subdetId(), matplotRender::t, SiStripDetId::TEC, TrackProducerBase< reco::GsfTrack >::theSchool, SiStripDetId::TIB, SiStripDetId::TID, SiStripDetId::TOB, align::Tracker, TrackProducerBase< reco::GsfTrack >::trajectoryInEvent_, TrajectoryMeasurement::updatedState(), useSplitting, v, TrackValidation_HighPurity_cff::valid, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().
Referenced by GsfTrackProducer::produce(), and GsfTrackRefitter::produce().
{
TrackingRecHitRefProd rHits = evt.getRefBeforePut<TrackingRecHitCollection>();
reco::TrackExtraRefProd rTrackExtras = evt.getRefBeforePut<reco::TrackExtraCollection>();
reco::GsfTrackExtraRefProd rGsfTrackExtras = evt.getRefBeforePut<reco::GsfTrackExtraCollection>();
reco::GsfTrackRefProd rTracks = evt.getRefBeforePut<reco::GsfTrackCollection>();
edm::Ref<reco::TrackExtraCollection>::key_type idx = 0;
edm::Ref<reco::TrackExtraCollection>::key_type hidx = 0;
edm::Ref<reco::GsfTrackExtraCollection>::key_type idxGsf = 0;
edm::Ref<reco::GsfTrackCollection>::key_type iTkRef = 0;
edm::Ref< std::vector<Trajectory> >::key_type iTjRef = 0;
std::map<unsigned int, unsigned int> tjTkMap;
TSCBLBuilderNoMaterial tscblBuilder;
for(AlgoProductCollection::iterator i=algoResults.begin(); i!=algoResults.end();i++){
Trajectory * theTraj = (*i).first;
if(trajectoryInEvent_) {
selTrajectories->push_back(*theTraj);
iTjRef++;
}
// const TrajectoryFitter::RecHitContainer& transHits = theTraj->recHits(useSplitting); // NO: the return type in Trajectory is by VALUE
TrajectoryFitter::RecHitContainer transHits = theTraj->recHits(useSplitting);
reco::GsfTrack * theTrack = (*i).second.first;
PropagationDirection seedDir = (*i).second.second;
LogDebug("TrackProducer") << "In GsfTrackProducerBase::putInEvt - seedDir=" << seedDir;
reco::GsfTrack t = * theTrack;
selTracks->push_back( t );
iTkRef++;
// Store indices in local map (starts at 0)
if(trajectoryInEvent_) tjTkMap[iTjRef-1] = iTkRef-1;
//sets the outermost and innermost TSOSs
TrajectoryStateOnSurface outertsos;
TrajectoryStateOnSurface innertsos;
unsigned int innerId, outerId;
// --- NOTA BENE: the convention is to sort hits and measurements "along the momentum".
// This is consistent with innermost and outermost labels only for tracks from LHC collision
if (theTraj->direction() == alongMomentum) {
outertsos = theTraj->lastMeasurement().updatedState();
innertsos = theTraj->firstMeasurement().updatedState();
outerId = theTraj->lastMeasurement().recHit()->geographicalId().rawId();
innerId = theTraj->firstMeasurement().recHit()->geographicalId().rawId();
} else {
outertsos = theTraj->firstMeasurement().updatedState();
innertsos = theTraj->lastMeasurement().updatedState();
outerId = theTraj->firstMeasurement().recHit()->geographicalId().rawId();
innerId = theTraj->lastMeasurement().recHit()->geographicalId().rawId();
}
//build the TrackExtra
GlobalPoint v = outertsos.globalParameters().position();
GlobalVector p = outertsos.globalParameters().momentum();
math::XYZVector outmom( p.x(), p.y(), p.z() );
math::XYZPoint outpos( v.x(), v.y(), v.z() );
v = innertsos.globalParameters().position();
p = innertsos.globalParameters().momentum();
math::XYZVector inmom( p.x(), p.y(), p.z() );
math::XYZPoint inpos( v.x(), v.y(), v.z() );
reco::TrackExtraRef teref= reco::TrackExtraRef ( rTrackExtras, idx ++ );
reco::GsfTrack & track = selTracks->back();
track.setExtra( teref );
//======= I want to set the second hitPattern here =============
if (theSchool.isValid())
{
NavigationSetter setter( *theSchool );
setSecondHitPattern(theTraj,track,prop,measTk);
}
//==============================================================
selTrackExtras->push_back( reco::TrackExtra (outpos, outmom, true, inpos, inmom, true,
outertsos.curvilinearError(), outerId,
innertsos.curvilinearError(), innerId,
seedDir, theTraj->seedRef()));
reco::TrackExtra & tx = selTrackExtras->back();
size_t i = 0;
// --- NOTA BENE: the convention is to sort hits and measurements "along the momentum".
// This is consistent with innermost and outermost labels only for tracks from LHC collisions
if (theTraj->direction() == alongMomentum) {
for( TrajectoryFitter::RecHitContainer::const_iterator j = transHits.begin();
j != transHits.end(); j ++ ) {
if ((**j).hit()!=0){
TrackingRecHit * hit = (**j).hit()->clone();
track.setHitPattern( * hit, i ++ );
selHits->push_back( hit );
tx.add( TrackingRecHitRef( rHits, hidx ++ ) );
}
}
}else{
for( TrajectoryFitter::RecHitContainer::const_iterator j = transHits.end()-1;
j != transHits.begin()-1; --j ) {
if ((**j).hit()!=0){
TrackingRecHit * hit = (**j).hit()->clone();
track.setHitPattern( * hit, i ++ );
selHits->push_back( hit );
tx.add( TrackingRecHitRef( rHits, hidx ++ ) );
}
}
}
// ----
std::vector<reco::GsfTangent> tangents;
const Trajectory::DataContainer& measurements = theTraj->measurements();
if ( measurements.size()>2 ) {
tangents.reserve(measurements.size()-2);
Trajectory::DataContainer::const_iterator ibegin,iend;
int increment(0);
if (theTraj->direction() == alongMomentum) {
ibegin = measurements.begin() + 1;
iend = measurements.end() - 1;
increment = 1;
}
else {
ibegin = measurements.end() - 2;
iend = measurements.begin();
increment = -1;
}
math::XYZPoint position;
math::XYZVector momentum;
Measurement1D deltaP;
// only measurements on "mono" detectors
for ( Trajectory::DataContainer::const_iterator i=ibegin;
i!=iend; i+=increment ) {
if ( i->recHit().get() ) {
DetId detId(i->recHit()->geographicalId());
if ( detId.det()==DetId::Tracker ) {
int subdetId = detId.subdetId();
if ( subdetId==SiStripDetId::TIB || subdetId==SiStripDetId::TID ||
subdetId==SiStripDetId::TOB || subdetId==SiStripDetId::TEC ) {
if ( SiStripDetId(detId).stereo() ) continue;
}
}
}
bool valid = computeModeAtTM(*i,position,momentum,deltaP);
if ( valid ) {
tangents.push_back(reco::GsfTangent(position,momentum,deltaP));
}
}
}
//build the GsfTrackExtra
std::vector<reco::GsfComponent5D> outerStates;
outerStates.reserve(outertsos.components().size());
fillStates(outertsos,outerStates);
std::vector<reco::GsfComponent5D> innerStates;
innerStates.reserve(innertsos.components().size());
fillStates(innertsos,innerStates);
reco::GsfTrackExtraRef terefGsf = reco::GsfTrackExtraRef ( rGsfTrackExtras, idxGsf ++ );
track.setGsfExtra( terefGsf );
selGsfTrackExtras->push_back( reco::GsfTrackExtra (outerStates, outertsos.localParameters().pzSign(),
innerStates, innertsos.localParameters().pzSign(),
tangents));
if ( innertsos.isValid() ) {
GsfPropagatorAdapter gsfProp(AnalyticalPropagator(innertsos.magneticField(),anyDirection));
TransverseImpactPointExtrapolator tipExtrapolator(gsfProp);
fillMode(track,innertsos,gsfProp,tipExtrapolator,tscblBuilder,bs);
}
delete theTrack;
delete theTraj;
}
LogTrace("TrackingRegressionTest") << "========== TrackProducer Info ===================";
LogTrace("TrackingRegressionTest") << "number of finalGsfTracks: " << selTracks->size();
for (reco::GsfTrackCollection::const_iterator it = selTracks->begin(); it != selTracks->end(); it++) {
LogTrace("TrackingRegressionTest") << "track's n valid and invalid hit, chi2, pt : "
<< it->found() << " , "
<< it->lost() <<" , "
<< it->normalizedChi2() << " , "
<< it->pt() << " , "
<< it->eta() ;
}
LogTrace("TrackingRegressionTest") << "=================================================";
rTracks_ = evt.put( selTracks );
evt.put( selTrackExtras );
evt.put( selGsfTrackExtras );
evt.put( selHits );
if(trajectoryInEvent_) {
edm::OrphanHandle<std::vector<Trajectory> > rTrajs = evt.put(selTrajectories);
// Now Create traj<->tracks association map
std::auto_ptr<TrajGsfTrackAssociationCollection> trajTrackMap( new TrajGsfTrackAssociationCollection() );
for ( std::map<unsigned int, unsigned int>::iterator i = tjTkMap.begin();
i != tjTkMap.end(); i++ ) {
edm::Ref<std::vector<Trajectory> > trajRef( rTrajs, (*i).first );
edm::Ref<reco::GsfTrackCollection> tkRef( rTracks_, (*i).second );
trajTrackMap->insert( edm::Ref<std::vector<Trajectory> >( rTrajs, (*i).first ),
edm::Ref<reco::GsfTrackCollection>( rTracks_, (*i).second ) );
}
evt.put( trajTrackMap );
}
}
bool GsfTrackProducerBase::useSplitting [private] |
Definition at line 69 of file GsfTrackProducerBase.h.
Referenced by putInEvt().
1.7.3