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