CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
TkGluedMeasurementDet.cc
Go to the documentation of this file.
10 #include "RecHitPropagator.h"
12 
14 #include <iostream>
15 
16 #include <typeinfo>
17 
18 // #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
20 
21 using namespace std;
22 
24  const SiStripRecHitMatcher* matcher,
25  const StripClusterParameterEstimator* cpe) :
26  MeasurementDet(gdet),
27  theMatcher(matcher), theCPE(cpe),
28  theMonoDet(nullptr), theStereoDet(nullptr)
29 {}
30 
32  const MeasurementDet* stereoDet) {
33  theMonoDet = dynamic_cast<const TkStripMeasurementDet *>(monoDet);
34  theStereoDet = dynamic_cast<const TkStripMeasurementDet *>(stereoDet);
35 
36  if ((theMonoDet == 0) || (theStereoDet == 0)) {
37  throw MeasurementDetException("TkGluedMeasurementDet ERROR: Trying to glue a det which is not a TkStripMeasurementDet");
38  }
39 }
40 
43 {
44 
46  HitCollectorForRecHits collector( &fastGeomDet(), theMatcher, theCPE, result );
47  collectRecHits(ts, collector);
48  return result;
49 }
50 
51 struct take_address { template<typename T> const T * operator()(const T &val) const { return &val; } };
52 
53 #ifdef DOUBLE_MATCH
54 template<typename Collector>
55 void
56 TkGluedMeasurementDet::collectRecHits( const TrajectoryStateOnSurface& ts, Collector & collector) const
57 {
58  doubleMatch(ts,collector);
59 }
60 #else
61 template<typename Collector>
62 void
63 TkGluedMeasurementDet::collectRecHits( const TrajectoryStateOnSurface& ts, Collector & collector) const
64 {
65  //------ WARNING: here ts is used as it is on the mono/stereo surface.
66  //----- A further propagation is necessary.
67  //----- To limit the problem, the SimpleCPE should be used
68  RecHitContainer monoHits = theMonoDet->recHits( ts);
69  GlobalVector glbDir = (ts.isValid() ? ts.globalParameters().momentum() : position()-GlobalPoint(0,0,0));
70 
71  //edm::LogWarning("TkGluedMeasurementDet::recHits") << "Query-for-detid-" << theGeomDet->geographicalId().rawId();
72 
73  //checkProjection(ts, monoHits, stereoHits);
74 
75  if (monoHits.empty()) {
76  // make stereo TTRHs and project them
77  projectOnGluedDet( collector, theStereoDet->recHits(ts), glbDir);
78  } else {
79  // collect simple stereo hits
80  static std::vector<SiStripRecHit2D> simpleSteroHitsByValue;
81  simpleSteroHitsByValue.clear();
82  theStereoDet->simpleRecHits(ts, simpleSteroHitsByValue);
83 
84  if (simpleSteroHitsByValue.empty()) {
85  projectOnGluedDet( collector, monoHits, glbDir);
86  } else {
87 
88  LocalVector tkDir = (ts.isValid() ? ts.localDirection() : surface().toLocal( position()-GlobalPoint(0,0,0)));
90  vsStereoHits.resize(simpleSteroHitsByValue.size());
91  std::transform(simpleSteroHitsByValue.begin(), simpleSteroHitsByValue.end(), vsStereoHits.begin(), take_address());
92 
93  // convert mono hits to type expected by matcher
94  for (RecHitContainer::const_iterator monoHit = monoHits.begin();
95  monoHit != monoHits.end(); ++monoHit) {
96  const TrackingRecHit* tkhit = (**monoHit).hit();
97  const SiStripRecHit2D* verySpecificMonoHit = reinterpret_cast<const SiStripRecHit2D*>(tkhit);
98  theMatcher->match( verySpecificMonoHit, vsStereoHits.begin(), vsStereoHits.end(),
99  collector.collector(), &specificGeomDet(), tkDir);
100 
101  if (collector.hasNewMatchedHits()) {
102  collector.clearNewMatchedHitsFlag();
103  } else {
104  collector.addProjected( **monoHit, glbDir );
105  }
106  } // loop on mono hit
107  }
108  //GIO// std::cerr << "TkGluedMeasurementDet hits " << monoHits.size() << "/" << stereoHits.size() << " => " << result.size() << std::endl;
109  }
110 }
111 #endif
112 
113 std::vector<TrajectoryMeasurement>
115  const TrajectoryStateOnSurface& startingState,
116  const Propagator&,
117  const MeasurementEstimator& est) const
118 {
119  std::vector<TrajectoryMeasurement> result;
120  if (theMonoDet->isActive() || theStereoDet->isActive()) {
121 
122  HitCollectorForFastMeasurements collector( &fastGeomDet(), theMatcher, theCPE, stateOnThisDet, est, result);
123  collectRecHits(stateOnThisDet, collector);
124 
125  if ( result.empty()) {
126  //LogDebug("TkStripMeasurementDet") << "No hit found on TkGlued. Testing strips... ";
127  const BoundPlane &gluedPlane = geomDet().surface();
128  if ( // sorry for the big IF, but I want to exploit short-circuiting of logic
129  stateOnThisDet.hasError() && ( /* do this only if the state has uncertainties, otherwise it will throw
130  (states without uncertainties are passed to this code from seeding */
131  (theMonoDet->isActive() &&
133  testStrips(stateOnThisDet,gluedPlane,*theMonoDet)
134  )
135  ) /*Mono OK*/ ||
136  (theStereoDet->isActive() &&
138  testStrips(stateOnThisDet,gluedPlane,*theStereoDet)
139  )
140  ) /*Stereo OK*/
141  ) /* State has errors */
142  ) {
143  result.push_back( TrajectoryMeasurement( stateOnThisDet,
145  } else {
146  result.push_back( TrajectoryMeasurement(stateOnThisDet,
148  }
149  } else {
150  // sort results according to estimator value
151  if ( result.size() > 1) {
152  sort( result.begin(), result.end(), TrajMeasLessEstim());
153  }
154  }
155  } else {
156  // LogDebug("TkStripMeasurementDet") << " DetID " << geomDet().geographicalId().rawId() << " (glued) fully inactive";
157  result.push_back( TrajectoryMeasurement( stateOnThisDet,
159  0.F));
160  }
161  return result;
162 
163 }
164 
167  const TrajectoryStateOnSurface& ts) const
168 {
169  if (hits.empty()) return hits;
172  for ( RecHitContainer::const_iterator ihit = hits.begin(); ihit!=hits.end(); ihit++) {
173  result.push_back( proj.project( **ihit, fastGeomDet(), ts));
174  }
175  return result;
176 }
177 
178 template<typename Collector>
179 void
181  const RecHitContainer& hits,
182  const GlobalVector & gdir ) const
183 {
184  for ( RecHitContainer::const_iterator ihit = hits.begin(); ihit!=hits.end(); ihit++) {
185  collector.addProjected( **ihit, gdir );
186  }
187 }
188 
190  const RecHitContainer& monoHits,
191  const RecHitContainer& stereoHits) const
192 {
193  for (RecHitContainer::const_iterator i=monoHits.begin(); i != monoHits.end(); ++i) {
194  checkHitProjection( **i, ts, fastGeomDet());
195  }
196  for (RecHitContainer::const_iterator i=stereoHits.begin(); i != stereoHits.end(); ++i) {
197  checkHitProjection( **i, ts, fastGeomDet());
198  }
199 }
200 
202  const TrajectoryStateOnSurface& ts,
203  const GeomDet& det) const
204 {
206  TransientTrackingRecHit::RecHitPointer projectedHit = proj.project( hit, det, ts);
207 
208  RecHitPropagator prop;
209  TrajectoryStateOnSurface propState = prop.propagate( hit, det.surface(), ts);
210 
211  if ((projectedHit->localPosition()-propState.localPosition()).mag() > 0.0001) {
212  cout << "PROBLEM: projected and propagated hit positions differ by "
213  << (projectedHit->localPosition()-propState.localPosition()).mag() << endl;
214  }
215 
216  LocalError le1 = projectedHit->localPositionError();
217  LocalError le2 = propState.localError().positionError();
218  double eps = 1.e-5;
219  double cutoff = 1.e-4; // if element below cutoff, use absolute instead of relative accuracy
220  double maxdiff = std::max( std::max( fabs(le1.xx() - le2.xx())/(cutoff+le1.xx()),
221  fabs(le1.xy() - le2.xy())/(cutoff+fabs(le1.xy()))),
222  fabs(le1.yy() - le2.yy())/(cutoff+le1.xx()));
223  if (maxdiff > eps) {
224  cout << "PROBLEM: projected and propagated hit errors differ by "
225  << maxdiff << endl;
226  }
227 
228 }
229 
230 bool
232  const BoundPlane &gluedPlane,
233  const TkStripMeasurementDet &mdet) const {
234  // from TrackingRecHitProjector
235  const GeomDet &det = mdet.fastGeomDet();
236  const BoundPlane &stripPlane = det.surface();
237 
238  //LocalPoint glp = tsos.localPosition();
239  LocalError err = tsos.localError().positionError();
240  /*LogDebug("TkStripMeasurementDet") <<
241  "Testing local pos glued: " << glp <<
242  " local err glued: " << tsos.localError().positionError() <<
243  " in? " << gluedPlane.bounds().inside(glp) <<
244  " in(3s)? " << gluedPlane.bounds().inside(glp, err, 3.0f);*/
245 
246  GlobalVector gdir = tsos.globalParameters().momentum();
247 
248  LocalPoint slp = stripPlane.toLocal(tsos.globalPosition());
249  LocalVector sld = stripPlane.toLocal(gdir);
250 
251  double delta = stripPlane.localZ( tsos.globalPosition());
252  LocalPoint pos = slp - sld * delta/sld.z();
253 
254 
255  // now the error
256  LocalVector hitXAxis = stripPlane.toLocal( gluedPlane.toGlobal( LocalVector(1,0,0)));
257  if (stripPlane.normalVector().dot( gluedPlane.normalVector()) < 0) {
258  // the two planes are inverted, and the correlation element must change sign
259  err = LocalError( err.xx(), -err.xy(), err.yy());
260  }
261  LocalError rotatedError = err.rotate( hitXAxis.x(), hitXAxis.y());
262 
263  /* // This is probably meaningless
264  LogDebug("TkStripMeasurementDet") <<
265  "Testing local pos on strip (SLP): " << slp <<
266  " in? :" << stripPlane.bounds().inside(slp) <<
267  " in(3s)? :" << stripPlane.bounds().inside(slp, rotatedError, 3.0f);
268  // but it helps to test bugs in the formula for POS */
269  /*LogDebug("TkStripMeasurementDet") <<
270  "Testing local pos strip: " << pos <<
271  " in? " << stripPlane.bounds().inside(pos) <<
272  " in(3s)? " << stripPlane.bounds().inside(pos, rotatedError, 3.0f);*/
273 
274  // now we need to convert to MeasurementFrame
275  const StripTopology &topo = mdet.specificGeomDet().specificTopology();
276  float utraj = topo.measurementPosition(pos).x();
277  float uerr = std::sqrt(topo.measurementError(pos,rotatedError).uu());
278  return mdet.testStrips(utraj, uerr);
279 }
280 
281 #include<boost/bind.hpp>
283  const SiStripRecHitMatcher * matcher, const StripClusterParameterEstimator* cpe,
285  geomDet_(geomDet), matcher_(matcher), cpe_(cpe),target_(target),
286  collector_(boost::bind(&HitCollectorForRecHits::add,boost::ref(*this),_1)),
287  hasNewHits_(false)
288 {
289 }
290 
291 void
293  const GlobalVector & gdir)
294 {
296  target_.push_back( proj.project( hit, *geomDet_, gdir));
297 }
298 
300  const SiStripRecHitMatcher * matcher, const StripClusterParameterEstimator* cpe,
301  const TrajectoryStateOnSurface& stateOnThisDet,
302  const MeasurementEstimator& est,
303  std::vector<TrajectoryMeasurement> & target) :
304  geomDet_(geomDet), matcher_(matcher), cpe_(cpe),stateOnThisDet_(stateOnThisDet), est_(est), target_(target),
305  collector_(boost::bind(&HitCollectorForFastMeasurements::add,boost::ref(*this),_1)),
306  hasNewHits_(false)
307 {
308 }
309 
310 void
312 {
313  static LocalCache<TSiStripMatchedRecHit> lcache; // in case of pool allocator it will be cleared centrally
314  std::auto_ptr<TSiStripMatchedRecHit> & cache = lcache.ptr;
315  TSiStripMatchedRecHit::buildInPlace( cache, geomDet_, &hit2d, matcher_, cpe_ );
316  std::pair<bool,double> diffEst = est_.estimate( stateOnThisDet_, *cache);
317  if ( diffEst.first) {
318  cache->clonePersistentHit(); // clone and take ownership of the persistent 2D hit
319  target_.push_back(
320  TrajectoryMeasurement( stateOnThisDet_,
321  RecHitPointer(cache.release()),
322  diffEst.second)
323  );
324  } else {
325  cache->clearPersistentHit(); // drop ownership
326  }
327  hasNewHits_ = true; //FIXME: see also what happens moving this within testAndPush
328 }
329 
330 void
332  const GlobalVector & gdir)
333 {
334  // here we're ok with some extra casual new's and delete's
336  RecHitPointer phit = proj.project( hit, *geomDet_, gdir );
337  std::pair<bool,double> diffEst = est_.estimate( stateOnThisDet_, *phit);
338  if ( diffEst.first) {
339  target_.push_back( TrajectoryMeasurement( stateOnThisDet_, phit, diffEst.second) );
340  }
341 
342 }
343 
344 
345 
346 #ifdef DOUBLE_MATCH
347 #include "doubleMatch.icc"
348 #endif
dbl * delta
Definition: mlp_gen.cc:36
GlobalPoint toGlobal(const Point2DBase< Scalar, LocalTag > lp) const
Definition: Surface.h:78
int i
Definition: DBlmapReader.cc:9
void collectRecHits(const TrajectoryStateOnSurface &, Collector &coll) const dso_internal
float xx() const
Definition: LocalError.h:24
const TkStripMeasurementDet * theMonoDet
HitCollectorForFastMeasurements(const GeomDet *geomDet, const SiStripRecHitMatcher *matcher, const StripClusterParameterEstimator *cpe, const TrajectoryStateOnSurface &stateOnThisDet, const MeasurementEstimator &est, std::vector< TrajectoryMeasurement > &target)
RecHitContainer projectOnGluedDet(const RecHitContainer &hits, const TrajectoryStateOnSurface &ts) const dso_internal
void addProjected(const TransientTrackingRecHit &hit, const GlobalVector &gdir)
void init(const MeasurementDet *monoDet, const MeasurementDet *stereoDet)
void checkHitProjection(const TransientTrackingRecHit &hit, const TrajectoryStateOnSurface &ts, const GeomDet &det) const dso_internal
TrajectoryStateOnSurface propagate(const TransientTrackingRecHit &hit, const Plane &plane, const TrajectoryStateOnSurface &ts) const
void checkProjection(const TrajectoryStateOnSurface &ts, const RecHitContainer &monoHits, const RecHitContainer &stereoHits) const
SiStripMatchedRecHit2D * match(const SiStripRecHit2D *monoRH, const SiStripRecHit2D *stereoRH, const GluedGeomDet *gluedDet, LocalVector trackdirection) const
T mag() const
The vector magnitude. Equivalent to sqrt(vec.mag2())
LocalVector localDirection() const
GlobalVector normalVector() const
Definition: Plane.h:47
TransientTrackingRecHit::ConstRecHitContainer RecHitContainer
#define nullptr
virtual const GeomDet & geomDet() const
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
T y() const
Definition: PV3DBase.h:62
float localZ(const GlobalPoint &gp) const
Fast access to distance from plane for a point.
Definition: Plane.h:52
GlobalPoint globalPosition() const
void add(const std::vector< const T * > &source, std::vector< const T * > &dest)
PreciseFloatType< T, U >::Type dot(const Vector3DBase< U, FrameTag > &v) const
Definition: Vector3DBase.h:107
const TkStripMeasurementDet * monoDet() const
void addProjected(const TransientTrackingRecHit &hit, const GlobalVector &gdir)
const BoundSurface & surface() const
virtual RecHitContainer recHits(const TrajectoryStateOnSurface &) const
virtual const StripTopology & specificTopology() const
Returns a reference to the strip proxy topology.
LocalError positionError() const
TkGluedMeasurementDet(const GluedGeomDet *gdet, const SiStripRecHitMatcher *matcher, const StripClusterParameterEstimator *cpe)
void doubleMatch(const TrajectoryStateOnSurface &ts, Collector &collector) const dso_internal
const GeomDet & fastGeomDet() const
RecHitPointer project(const TransientTrackingRecHit &hit, const GeomDet &det, const TrajectoryStateOnSurface &ts) const
float xy() const
Definition: LocalError.h:25
virtual MeasurementError measurementError(const LocalPoint &, const LocalError &) const =0
float yy() const
Definition: LocalError.h:26
const T & max(const T &a, const T &b)
T sqrt(T t)
Definition: SSEVec.h:46
LocalPoint toLocal(const GlobalPoint &gp) const
const GluedGeomDet & specificGeomDet() const
T z() const
Definition: PV3DBase.h:63
tuple result
Definition: query.py:137
float uu() const
bool testStrips(const TrajectoryStateOnSurface &tsos, const BoundPlane &gluedPlane, const TkStripMeasurementDet &mdet) const dso_internal
Test the strips on one of the two dets with projection.
virtual MeasurementPoint measurementPosition(const LocalPoint &) const =0
const LocalTrajectoryError & localError() const
bool hasAllGoodChannels() const
does this module have at least one bad strip, APV or channel?
const T * operator()(const T &val) const
const Surface::PositionType & position() const
virtual RecHitContainer recHits(const TrajectoryStateOnSurface &) const
const SiStripRecHitMatcher * theMatcher
const GlobalTrajectoryParameters & globalParameters() const
virtual std::vector< TrajectoryMeasurement > fastMeasurements(const TrajectoryStateOnSurface &stateOnThisDet, const TrajectoryStateOnSurface &startingState, const Propagator &, const MeasurementEstimator &) const
std::vector< const SiStripRecHit2D * > SimpleHitCollection
const TkStripMeasurementDet * theStereoDet
const StripGeomDetUnit & specificGeomDet() const
void simpleRecHits(const TrajectoryStateOnSurface &ts, std::vector< SiStripRecHit2D > &result) const
const BoundPlane & surface() const
The nominal surface of the GeomDet.
Definition: GeomDet.h:35
bool testStrips(float utraj, float uerr) const
return true if there are &#39;enough&#39; good strips in the utraj +/- 3 uerr range.
bool isActive() const
Is this module active in reconstruction? It must be both &#39;setActiveThisEvent&#39; and &#39;setActive&#39;...
tuple cout
Definition: gather_cfg.py:121
const TkStripMeasurementDet * stereoDet() const
static uInt32 F(BLOWFISH_CTX *ctx, uInt32 x)
Definition: blowfish.cc:281
LocalError rotate(float x, float y) const
Return a new LocalError, rotated by an angle defined by the direction (x,y)
Definition: LocalError.h:39
const StripClusterParameterEstimator * theCPE
T x() const
Definition: PV2DBase.h:44
long double T
T x() const
Definition: PV3DBase.h:61
HitCollectorForRecHits(const GeomDet *geomDet, const SiStripRecHitMatcher *matcher, const StripClusterParameterEstimator *cpe, RecHitContainer &target)
std::auto_ptr< T > ptr