00001 #include "RecoTracker/TkSeedGenerator/interface/SeedFromConsecutiveHitsCreator.h"
00002 #include "RecoTracker/TkTrackingRegions/interface/TrackingRegion.h"
00003
00004 #include "TrackingTools/KalmanUpdators/interface/KFUpdator.h"
00005 #include "RecoTracker/TkSeedGenerator/interface/FastHelix.h"
00006 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
00007 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
00008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00009 #include "MagneticField/Engine/interface/MagneticField.h"
00010 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
00011 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
00012 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
00013 #include "TrackingTools/MaterialEffects/interface/PropagatorWithMaterial.h"
00014 #include "FWCore/Framework/interface/ESHandle.h"
00015 #include "TrackingTools/GeomPropagators/interface/PropagationExceptions.h"
00016 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
00017 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
00018
00019 template <class T> T sqr( T t) {return t*t;}
00020
00021 const TrajectorySeed * SeedFromConsecutiveHitsCreator::trajectorySeed(
00022 TrajectorySeedCollection & seedCollection,
00023 const SeedingHitSet & hits,
00024 const TrackingRegion & region,
00025 const edm::EventSetup& es)
00026 {
00027 if ( hits.size() < 2) return 0;
00028
00029 GlobalTrajectoryParameters kine = initialKinematic(hits, region, es);
00030 float sinTheta = sin(kine.momentum().theta());
00031
00032 CurvilinearTrajectoryError error = initialError(region, sinTheta);
00033 FreeTrajectoryState fts(kine, error);
00034
00035 return buildSeed(seedCollection,hits,fts,es);
00036 }
00037
00038
00039 GlobalTrajectoryParameters SeedFromConsecutiveHitsCreator::initialKinematic(
00040 const SeedingHitSet & hits,
00041 const TrackingRegion & region,
00042 const edm::EventSetup& es) const
00043 {
00044 edm::ESHandle<MagneticField> bfield;
00045 es.get<IdealMagneticFieldRecord>().get(bfield);
00046
00047 GlobalTrajectoryParameters kine;
00048
00049 TransientTrackingRecHit::ConstRecHitPointer tth1 = hits[0];
00050 TransientTrackingRecHit::ConstRecHitPointer tth2 = hits[1];
00051 const GlobalPoint& vertexPos = region.origin();
00052
00053 FastHelix helix(tth2->globalPosition(), tth1->globalPosition(), vertexPos, es);
00054 if (helix.isValid()) {
00055 kine = helix.stateAtVertex().parameters();
00056 } else {
00057 GlobalVector initMomentum(tth2->globalPosition() - vertexPos);
00058 initMomentum *= (100./initMomentum.perp());
00059 kine = GlobalTrajectoryParameters(vertexPos, initMomentum, 1, &*bfield);
00060 }
00061
00062 bool isBOFF = ( std::abs(bfield->inTesla(GlobalPoint(0,0,0)).z()) < 1e-3 );
00063 if (isBOFF && (theBOFFMomentum > 0)) {
00064 kine = GlobalTrajectoryParameters(kine.position(),
00065 kine.momentum().unit() * theBOFFMomentum,
00066 kine.charge(),
00067 &*bfield);
00068 }
00069 return kine;
00070 }
00071
00072
00073
00074 CurvilinearTrajectoryError SeedFromConsecutiveHitsCreator::
00075 initialError( const TrackingRegion& region, float sinTheta) const
00076 {
00077
00078
00079 GlobalError vertexErr( sqr(region.originRBound()), 0, sqr(region.originRBound()),
00080 0, 0, sqr(region.originZBound()));
00081
00082 float ptMin = region.ptMin();
00083
00084
00085 AlgebraicSymMatrix C(5,1);
00086
00087
00088
00089
00090 float sin2th = sqr(sinTheta);
00091 float minC00 = 1.0;
00092 C[0][0] = std::max(sin2th/sqr(ptMin), minC00);
00093 float zErr = vertexErr.czz();
00094 float transverseErr = vertexErr.cxx();
00095 C[3][3] = transverseErr;
00096 C[4][4] = zErr*sin2th + transverseErr*(1-sin2th);
00097
00098 return CurvilinearTrajectoryError(C);
00099 }
00100
00101 const TrajectorySeed * SeedFromConsecutiveHitsCreator::buildSeed(
00102 TrajectorySeedCollection & seedCollection,
00103 const SeedingHitSet & hits,
00104 const FreeTrajectoryState & fts,
00105 const edm::EventSetup& es) const
00106 {
00107
00108 edm::ESHandle<TrackerGeometry> tracker;
00109 es.get<TrackerDigiGeometryRecord>().get(tracker);
00110
00111
00112 edm::ESHandle<Propagator> propagatorHandle;
00113 es.get<TrackingComponentsRecord>().get(thePropagatorLabel, propagatorHandle);
00114 const Propagator* propagator = &(*propagatorHandle);
00115
00116
00117 KFUpdator updator;
00118
00119
00120
00121 TrajectoryStateOnSurface updatedState;
00122 edm::OwnVector<TrackingRecHit> seedHits;
00123
00124 const TrackingRecHit* hit = 0;
00125 for ( unsigned int iHit = 0; iHit < hits.size(); iHit++) {
00126 hit = hits[iHit]->hit();
00127 TrajectoryStateOnSurface state = (iHit==0) ?
00128 propagator->propagate(fts,tracker->idToDet(hit->geographicalId())->surface())
00129 : propagator->propagate(updatedState, tracker->idToDet(hit->geographicalId())->surface());
00130 if (!state.isValid()) return 0;
00131
00132 TransientTrackingRecHit::ConstRecHitPointer tth = hits[iHit];
00133
00134 TransientTrackingRecHit::RecHitPointer newtth = refitHit( tth, state);
00135
00136 if (!checkHit(state,newtth,es)) return 0;
00137
00138 updatedState = updator.update(state, *newtth);
00139 if (!updatedState.isValid()) return 0;
00140
00141 seedHits.push_back(newtth->hit()->clone());
00142 }
00143
00144 TrajectoryStateTransform transformer;
00145 boost::shared_ptr<PTrajectoryStateOnDet> PTraj(
00146 transformer.persistentState(updatedState, hit->geographicalId().rawId()));
00147
00148 seedCollection.push_back( TrajectorySeed(*PTraj,seedHits,alongMomentum));
00149 return &seedCollection.back();
00150 }
00151
00152 TransientTrackingRecHit::RecHitPointer SeedFromConsecutiveHitsCreator::refitHit(
00153 const TransientTrackingRecHit::ConstRecHitPointer &hit,
00154 const TrajectoryStateOnSurface &state) const
00155 {
00156 return hit->clone(state);
00157 }