CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
PixelHitMatcher.cc
Go to the documentation of this file.
1 
13 
14 #include <typeinfo>
15 
16 using namespace reco ;
17 using namespace std ;
18 
20  ( float phi1min, float phi1max,
21  float phi2minB, float phi2maxB, float phi2minF, float phi2maxF,
22  float z2minB, float z2maxB, float r2minF, float r2maxF,
23  float rMinI, float rMaxI, bool searchInTIDTEC)
24  : //zmin1 and zmax1 are dummy at this moment, set from beamspot later
25  meas1stBLayer(phi1min,phi1max,0.,0.), meas2ndBLayer(phi2minB,phi2maxB,z2minB,z2maxB),
26  meas1stFLayer(phi1min,phi1max,0.,0.), meas2ndFLayer(phi2minF,phi2maxF,r2minF,r2maxF),
27  startLayers(),
28  prop1stLayer(0), prop2ndLayer(0),theGeometricSearchTracker(0),theTrackerEvent(0),theTracker(0),vertex_(0.),
29  searchInTIDTEC_(searchInTIDTEC), useRecoVertex_(false)
30  {
31  meas1stFLayer.setRRangeI(rMinI,rMaxI) ;
32  meas2ndFLayer.setRRangeI(rMinI,rMaxI) ;
33  }
34 
36  {
37  delete prop1stLayer ;
38  delete prop2ndLayer ;
39  }
40 
41 void PixelHitMatcher::set1stLayer( float dummyphi1min, float dummyphi1max )
42  {
43  meas1stBLayer.setPhiRange(dummyphi1min,dummyphi1max) ;
44  meas1stFLayer.setPhiRange(dummyphi1min,dummyphi1max) ;
45  }
46 
47 void PixelHitMatcher::set1stLayerZRange( float zmin1, float zmax1 )
48  {
49  meas1stBLayer.setZRange(zmin1,zmax1) ;
50  meas1stFLayer.setRRange(zmin1,zmax1) ;
51  }
52 
53 void PixelHitMatcher::set2ndLayer( float dummyphi2minB, float dummyphi2maxB, float dummyphi2minF, float dummyphi2maxF )
54  {
55  meas2ndBLayer.setPhiRange(dummyphi2minB,dummyphi2maxB) ;
56  meas2ndFLayer.setPhiRange(dummyphi2minF,dummyphi2maxF) ;
57  }
58 
60  { useRecoVertex_ = val ; }
61 
63  {
64  theTrackerEvent = & trackerData;
65  theLayerMeasurements = LayerMeasurements(*theTracker,*theTrackerEvent);
66  }
68  ( const MagneticField * magField,
69  const MeasurementTracker * theMeasurementTracker,
70  const TrackerGeometry * trackerGeometry )
71  {
72  if (theMeasurementTracker)
73  {
74  theTracker = theMeasurementTracker;
75  theGeometricSearchTracker=theMeasurementTracker->geometricSearchTracker() ;
76  startLayers.setup(theGeometricSearchTracker) ;
77  }
78 
79  theMagField = magField ;
80  theTrackerGeometry = trackerGeometry ;
81  float mass=.000511 ; // electron propagation
82  if (prop1stLayer) delete prop1stLayer ;
83  prop1stLayer = new PropagatorWithMaterial(oppositeToMomentum,mass,theMagField) ;
84  if (prop2ndLayer) delete prop2ndLayer ;
85  prop2ndLayer = new PropagatorWithMaterial(alongMomentum,mass,theMagField) ;
86  }
87 
88 vector<CLHEP::Hep3Vector> PixelHitMatcher::predicted1Hits()
89  { return pred1Meas ; }
90 
91 vector<CLHEP::Hep3Vector> PixelHitMatcher::predicted2Hits()
92  { return pred2Meas ; }
93 
95  { return vertex_ ; }
96 
97 //CLHEP::Hep3Vector point_to_vector( const GlobalPoint & p )
98 // { return CLHEP::Hep3Vector(p.x(),p.y(),p.z()) ; }
99 
100 std::vector<SeedWithInfo>
102  ( TrajectorySeedCollection * seeds, const GlobalPoint & xmeas,
103  const GlobalPoint & vprim, float energy, float fcharge )
104  {
105  typedef std::unordered_map<const GeomDet *, TrajectoryStateOnSurface> DetTsosAssoc;
106  typedef std::unordered_map<std::pair<const GeomDet*,GlobalPoint>, TrajectoryStateOnSurface> PosTsosAssoc;
107  const int charge = int(fcharge) ;
108 
109  const double xmeas_phi = xmeas.phi();
110  FreeTrajectoryState fts = FTSFromVertexToPointFactory::get(*theMagField, xmeas, vprim, energy, charge);
112  TrajectoryStateOnSurface tsos(fts, *bpb(fts.position(), fts.momentum()));
113 
114  std::vector<SeedWithInfo> result ;
115 
116  mapTsos_fast_.clear();
117  mapTsos2_fast_.clear();
118  mapTsos_fast_.reserve(seeds->size()) ;
119  mapTsos2_fast_.reserve(seeds->size()) ;
120 
121  for(const auto& seed : *seeds) {
122  hit_gp_map_.clear();
123  if( seed.nHits() > 9 ) {
124  edm::LogWarning("GsfElectronAlgo|UnexpectedSeed") <<"We cannot deal with seeds having more than 9 hits." ;
125  continue;
126  }
127  const TrajectorySeed::range& hits = seed.recHits();
128  // cache the global points
129  for( auto it = hits.first; it != hits.second; ++it ) {
130  hit_gp_map_.emplace_back(it->globalPosition());
131  }
132  //iterate on the hits
133  for( auto it1 = hits.first; it1 != hits.second; ++it1 ) {
134  if( !it1->isValid() ) continue;
135  const unsigned idx1 = std::distance(hits.first,it1);
136  const DetId id1 = it1->geographicalId();
137  const GeomDet *geomdet1 = it1->det();
138  const GlobalPoint& hit1Pos = hit_gp_map_[idx1];
139 
140  const TrajectoryStateOnSurface* tsos1;
141  DetTsosAssoc::iterator tsos1_itr = mapTsos_fast_.find(geomdet1);
142  if( tsos1_itr != mapTsos_fast_.end() ) {
143  tsos1 = &(tsos1_itr->second);
144  } else {
145  auto empl_result =
146  mapTsos_fast_.emplace(geomdet1,prop1stLayer->propagate(tsos,geomdet1->surface()));
147  tsos1 = &(empl_result.first->second);
148  }
149  if( !tsos1->isValid() ) continue;
150  std::pair<bool, double> est = ( id1.subdetId() % 2 ?
151  meas1stBLayer.estimate(vprim, *tsos1, hit1Pos) :
152  meas1stFLayer.estimate(vprim, *tsos1, hit1Pos) );
153  if( !est.first ) continue;
154  if( std::abs(normalized_phi(hit1Pos.phi()-xmeas_phi))>2.5 ) continue;
155  EleRelPointPair pp1(hit1Pos,tsos1->globalParameters().position(),vprim);
156  const math::XYZPoint relHit1Pos(hit1Pos-vprim), relTSOSPos(tsos1->globalParameters().position() - vprim);
157  const int subDet1 = id1.subdetId();
158  const float dRz1 = ( id1.subdetId()%2 ? pp1.dZ() : pp1.dPerp() );
159  const float dPhi1 = pp1.dPhi();
160  // setup our vertex
161  double zVertex;
162  if (!useRecoVertex_) {
163  // we don't know the z vertex position, get it from linear extrapolation
164  // compute the z vertex from the cluster point and the found pixel hit
165  const double pxHit1z = hit1Pos.z();
166  const double pxHit1x = hit1Pos.x();
167  const double pxHit1y = hit1Pos.y();
168  const double r1diff = std::sqrt( (pxHit1x-vprim.x())*(pxHit1x-vprim.x()) +
169  (pxHit1y-vprim.y())*(pxHit1y-vprim.y()) );
170  const double r2diff = std::sqrt( (xmeas.x()-pxHit1x)*(xmeas.x()-pxHit1x) +
171  (xmeas.y()-pxHit1y)*(xmeas.y()-pxHit1y) );
172  zVertex = pxHit1z - r1diff*(xmeas.z()-pxHit1z)/r2diff;
173  } else {
174  // here use rather the reco vertex z position
175  zVertex = vprim.z();
176  }
177  GlobalPoint vertex(vprim.x(),vprim.y(),zVertex);
178  FreeTrajectoryState fts2 = FTSFromVertexToPointFactory::get(*theMagField, hit1Pos, vertex, energy, charge) ;
179  // now find the matching hit
180  for( auto it2 = it1+1; it2 != hits.second; ++it2 ) {
181  if( !it2->isValid() ) continue;
182  const unsigned idx2 = std::distance(hits.first,it2);
183  const DetId id2 = it2->geographicalId();
184  const GeomDet *geomdet2 = it2->det();
185  const std::pair<const GeomDet*,GlobalPoint> det_key(geomdet2,hit1Pos);
186  const TrajectoryStateOnSurface* tsos2;
187  PosTsosAssoc::iterator tsos2_itr = mapTsos2_fast_.find(det_key);
188  if( tsos2_itr != mapTsos2_fast_.end() ) {
189  tsos2 = &(tsos2_itr->second);
190  } else {
191  auto empl_result =
192  mapTsos2_fast_.emplace(det_key,prop2ndLayer->propagate(fts2,geomdet2->surface()));
193  tsos2 = &(empl_result.first->second);
194  }
195  if( !tsos2->isValid() ) continue;
196  const GlobalPoint& hit2Pos = hit_gp_map_[idx2];
197  std::pair<bool,double> est2 = ( id2.subdetId()%2 ?
198  meas2ndBLayer.estimate(vertex, *tsos2,hit2Pos) :
199  meas2ndFLayer.estimate(vertex, *tsos2,hit2Pos) );
200  if (est2.first) {
201  EleRelPointPair pp2(hit2Pos,tsos2->globalParameters().position(),vertex) ;
202  const int subDet2 = id2.subdetId();
203  const float dRz2 = (subDet2%2==1)?pp2.dZ():pp2.dPerp();
204  const float dPhi2 = pp2.dPhi();
205  const unsigned char hitsMask = (1<<idx1)|(1<<idx2);
206  result.push_back(SeedWithInfo(seed,hitsMask,subDet2,dRz2,dPhi2,subDet1,dRz1,dPhi1)) ;
207  }
208  }// inner loop on hits
209  }// outer loop on hits
210  }// loop on seeds
211 
212  mapTsos_fast_.clear() ;
213  mapTsos2_fast_.clear() ;
214 
215  return result ;
216  }
217 
218 //========================= OBSOLETE ? =========================
219 
220 vector< pair< RecHitWithDist, PixelHitMatcher::ConstRecHitPointer > >
222  ( const GlobalPoint & xmeas,
223  const GlobalPoint & vprim,
224  float energy, float fcharge,
225  const TrackerTopology *tTopo,
226  const NavigationSchool& navigationSchool)
227  {
228  float SCl_phi = xmeas.phi();
229 
230  int charge = int(fcharge);
231  // return all compatible RecHit pairs (vector< TSiPixelRecHit>)
232  vector<pair<RecHitWithDist, ConstRecHitPointer> > result;
233  LogDebug("") << "[PixelHitMatcher::compatibleHits] entering .. ";
234 
235  vector<TrajectoryMeasurement> validMeasurements;
236  vector<TrajectoryMeasurement> invalidMeasurements;
237 
238  typedef vector<TrajectoryMeasurement>::const_iterator aMeas;
239 
240  pred1Meas.clear();
241  pred2Meas.clear();
242 
243  typedef vector<const BarrelDetLayer*>::const_iterator BarrelLayerIterator;
244  BarrelLayerIterator firstLayer = startLayers.firstBLayer();
245 
246  FreeTrajectoryState fts = FTSFromVertexToPointFactory::get(*theMagField,xmeas, vprim, energy, charge);
247 
249  TrajectoryStateOnSurface tsos(fts, *bpb(fts.position(), fts.momentum()));
250 
251  if (tsos.isValid()) {
252  vector<TrajectoryMeasurement> pixelMeasurements =
253  theLayerMeasurements.measurements(**firstLayer,tsos,
254  *prop1stLayer, meas1stBLayer);
255 
256  LogDebug("") <<"[PixelHitMatcher::compatibleHits] nbr of hits compatible with extrapolation to first layer: " << pixelMeasurements.size();
257  for (aMeas m=pixelMeasurements.begin(); m!=pixelMeasurements.end(); m++){
258  if (m->recHit()->isValid()) {
259  float localDphi = normalized_phi(SCl_phi-m->forwardPredictedState().globalPosition().phi()) ;
260  if(std::abs(localDphi)>2.5)continue;
261  CLHEP::Hep3Vector prediction(m->forwardPredictedState().globalPosition().x(),
262  m->forwardPredictedState().globalPosition().y(),
263  m->forwardPredictedState().globalPosition().z());
264  LogDebug("") << "[PixelHitMatcher::compatibleHits] compatible hit position " << m->recHit()->globalPosition();
265  LogDebug("") << "[PixelHitMatcher::compatibleHits] predicted position " << m->forwardPredictedState().globalPosition();
266  pred1Meas.push_back( prediction);
267 
268  validMeasurements.push_back(*m);
269 
270  LogDebug("") <<"[PixelHitMatcher::compatibleHits] Found a rechit in layer ";
271  const BarrelDetLayer *bdetl = dynamic_cast<const BarrelDetLayer *>(*firstLayer);
272  if (bdetl) {
273  LogDebug("") <<" with radius "<<bdetl->specificSurface().radius();
274  }
275  else LogDebug("") <<"Could not downcast!!";
276  }
277  }
278 
279 
280  // check if there are compatible 1st hits in the second layer
281  firstLayer++;
282 
283  vector<TrajectoryMeasurement> pixel2Measurements =
284  theLayerMeasurements.measurements(**firstLayer,tsos,
285  *prop1stLayer, meas1stBLayer);
286 
287  for (aMeas m=pixel2Measurements.begin(); m!=pixel2Measurements.end(); m++){
288  if (m->recHit()->isValid()) {
289  float localDphi = normalized_phi(SCl_phi-m->forwardPredictedState().globalPosition().phi()) ;
290  if(std::abs(localDphi)>2.5)continue;
291  CLHEP::Hep3Vector prediction(m->forwardPredictedState().globalPosition().x(),
292  m->forwardPredictedState().globalPosition().y(),
293  m->forwardPredictedState().globalPosition().z());
294  pred1Meas.push_back( prediction);
295  LogDebug("") << "[PixelHitMatcher::compatibleHits] compatible hit position " << m->recHit()->globalPosition() << endl;
296  LogDebug("") << "[PixelHitMatcher::compatibleHits] predicted position " << m->forwardPredictedState().globalPosition() << endl;
297 
298  validMeasurements.push_back(*m);
299  LogDebug("") <<"[PixelHitMatcher::compatibleHits] Found a rechit in layer ";
300  const BarrelDetLayer *bdetl = dynamic_cast<const BarrelDetLayer *>(*firstLayer);
301  if (bdetl) {
302  LogDebug("") <<" with radius "<<bdetl->specificSurface().radius();
303  }
304  else LogDebug("") <<"Could not downcast!!";
305  }
306 
307  }
308  }
309 
310 
311  // check if there are compatible 1st hits the forward disks
312  typedef vector<const ForwardDetLayer*>::const_iterator ForwardLayerIterator;
313  ForwardLayerIterator flayer;
314 
315  TrajectoryStateOnSurface tsosfwd(fts, *bpb(fts.position(), fts.momentum()));
316  if (tsosfwd.isValid()) {
317 
318  for (int i=0; i<2; i++) {
319  i == 0 ? flayer = startLayers.pos1stFLayer() : flayer = startLayers.neg1stFLayer();
320 
321  if (i==0 && xmeas.z() < -100. ) continue;
322  if (i==1 && xmeas.z() > 100. ) continue;
323 
324  vector<TrajectoryMeasurement> pixelMeasurements =
325  theLayerMeasurements.measurements(**flayer, tsosfwd,
326  *prop1stLayer, meas1stFLayer);
327 
328  for (aMeas m=pixelMeasurements.begin(); m!=pixelMeasurements.end(); m++){
329  if (m->recHit()->isValid()) {
330  float localDphi = normalized_phi(SCl_phi-m->forwardPredictedState().globalPosition().phi());
331  if(std::abs(localDphi)>2.5)continue;
332  CLHEP::Hep3Vector prediction(m->forwardPredictedState().globalPosition().x(),
333  m->forwardPredictedState().globalPosition().y(),
334  m->forwardPredictedState().globalPosition().z());
335  pred1Meas.push_back( prediction);
336 
337  validMeasurements.push_back(*m);
338  }
339  }
340 
341  //check if there are compatible 1st hits the outer forward disks
342  if (searchInTIDTEC_) {
343  flayer++;
344 
345  vector<TrajectoryMeasurement> pixel2Measurements =
346  theLayerMeasurements.measurements(**flayer, tsosfwd,
347  *prop1stLayer, meas1stFLayer);
348 
349  for (aMeas m=pixel2Measurements.begin(); m!=pixel2Measurements.end(); m++){
350  if (m->recHit()->isValid()) {
351  float localDphi = normalized_phi(SCl_phi-m->forwardPredictedState().globalPosition().phi()) ;
352  if(std::abs(localDphi)>2.5)continue;
353  CLHEP::Hep3Vector prediction(m->forwardPredictedState().globalPosition().x(),
354  m->forwardPredictedState().globalPosition().y(),
355  m->forwardPredictedState().globalPosition().z());
356  pred1Meas.push_back( prediction);
357 
358  validMeasurements.push_back(*m);
359  }
360  // else{std::cout<<" hit non valid "<<std::endl; }
361  } //end 1st hit in outer f disk
362  }
363  }
364  }
365 
366  // now we have the vector of all valid measurements of the first point
367  for (unsigned i=0; i<validMeasurements.size(); i++){
368 
369  const DetLayer * newLayer = theGeometricSearchTracker->detLayer(validMeasurements[i].recHit()->det()->geographicalId());
370 
371  double zVertex ;
372  if (!useRecoVertex_)
373  {
374  // we don't know the z vertex position, get it from linear extrapolation
375  // compute the z vertex from the cluster point and the found pixel hit
376  double pxHit1z = validMeasurements[i].recHit()->det()->surface().toGlobal(
377  validMeasurements[i].recHit()->localPosition()).z();
378  double pxHit1x = validMeasurements[i].recHit()->det()->surface().toGlobal(
379  validMeasurements[i].recHit()->localPosition()).x();
380  double pxHit1y = validMeasurements[i].recHit()->det()->surface().toGlobal(
381  validMeasurements[i].recHit()->localPosition()).y();
382  double r1diff = (pxHit1x-vprim.x())*(pxHit1x-vprim.x()) + (pxHit1y-vprim.y())*(pxHit1y-vprim.y());
383  r1diff=sqrt(r1diff);
384  double r2diff = (xmeas.x()-pxHit1x)*(xmeas.x()-pxHit1x) + (xmeas.y()-pxHit1y)*(xmeas.y()-pxHit1y);
385  r2diff=sqrt(r2diff);
386  zVertex = pxHit1z - r1diff*(xmeas.z()-pxHit1z)/r2diff;
387  }
388  else
389  {
390  // here we use the reco vertex z position
391  zVertex = vprim.z();
392  }
393 
394  if (i==0)
395  { vertex_ = zVertex; }
396 
397  GlobalPoint vertexPred(vprim.x(),vprim.y(),zVertex) ;
398  GlobalPoint hitPos( validMeasurements[i].recHit()->det()->surface().toGlobal( validMeasurements[i].recHit()->localPosition() ) ) ;
399 
400  FreeTrajectoryState secondFTS = FTSFromVertexToPointFactory::get(*theMagField, hitPos, vertexPred, energy, charge);
401 
402  PixelMatchNextLayers secondHit(&theLayerMeasurements, newLayer, secondFTS,
403  prop2ndLayer, &meas2ndBLayer,&meas2ndFLayer,
404  tTopo,navigationSchool,searchInTIDTEC_);
405  vector<CLHEP::Hep3Vector> predictions = secondHit.predictionInNextLayers();
406 
407  for (unsigned it = 0; it < predictions.size(); it++) pred2Meas.push_back(predictions[it]);
408 
409  // we may get more than one valid second measurements here even for single electrons:
410  // two hits from the same layer/disk (detector overlap) or from the loop over the
411  // next layers in EPMatchLoopNextLayers. Take only the 1st hit.
412 
413  if(!secondHit.measurementsInNextLayers().empty()){
414  for(unsigned int shit=0; shit<secondHit.measurementsInNextLayers().size(); shit++)
415  {
416  float dphi = normalized_phi(pred1Meas[i].phi()-validMeasurements[i].recHit()->globalPosition().phi()) ;
417  if (std::abs(dphi)<2.5)
418  {
419  ConstRecHitPointer pxrh = validMeasurements[i].recHit();
420  RecHitWithDist rh(pxrh,dphi);
421 
422  // pxrh = secondHit.measurementsInNextLayers()[0].recHit();
423  pxrh = secondHit.measurementsInNextLayers()[shit].recHit();
424 
425  pair<RecHitWithDist,ConstRecHitPointer> compatiblePair = pair<RecHitWithDist,ConstRecHitPointer>(rh,pxrh) ;
426  result.push_back(compatiblePair);
427  break;
428  }
429  }
430  }
431  }
432  return result;
433 }
434 
435 
#define LogDebug(id)
GlobalPoint toGlobal(const Point2DBase< Scalar, LocalTag > lp) const
Definition: Surface.h:106
virtual const BoundSurface & surface() const =0
The surface of the GeometricSearchDet.
int i
Definition: DBlmapReader.cc:9
static FreeTrajectoryState get(MagneticField const &magField, GlobalPoint const &xmeas, GlobalPoint const &xvert, float momentum, TrackCharge charge)
void setEvent(const MeasurementTrackerEvent &event)
void setES(const MagneticField *, const MeasurementTracker *theMeasurementTracker, const TrackerGeometry *trackerGeometry)
Geom::Phi< T > phi() const
Definition: PV3DBase.h:69
T y() const
Definition: PV3DBase.h:63
std::vector< CLHEP::Hep3Vector > predictionInNextLayers() const
const Plane & surface() const
The nominal surface of the GeomDet.
Definition: GeomDet.h:40
TransientTrackingRecHit::ConstRecHitPointer ConstRecHitPointer
void set2ndLayer(float dummyphi2minB, float dummyphi2maxB, float dummyphi2minF, float dummyphi2maxF)
tuple result
Definition: mps_fire.py:95
void set1stLayerZRange(float zmin1, float zmax1)
PixelHitMatcher(float phi1min, float phi1max, float phi2minB, float phi2maxB, float phi2minF, float phi2maxF, float z2minB, float z2maxB, float r2minF, float r2maxF, float rMinI, float rMaxI, bool searchInTIDTEC)
std::vector< TrajectorySeed > TrajectorySeedCollection
RealType normalized_phi(RealType phi)
T sqrt(T t)
Definition: SSEVec.h:18
T z() const
Definition: PV3DBase.h:64
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
std::pair< const_iterator, const_iterator > range
std::vector< CLHEP::Hep3Vector > predicted1Hits()
int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:37
GlobalVector momentum() const
Definition: DetId.h:18
GlobalPoint position() const
std::vector< TrajectoryMeasurement > measurementsInNextLayers() const
const GlobalTrajectoryParameters & globalParameters() const
void set1stLayer(float dummyphi1min, float dummyphi1max)
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
virtual const BoundCylinder & specificSurface() const final
Extension of the interface.
Geom::Phi< T > phi() const
virtual ~PixelHitMatcher()
tuple MeasurementTracker
std::vector< SeedWithInfo > compatibleSeeds(TrajectorySeedCollection *seeds, const GlobalPoint &xmeas, const GlobalPoint &vprim, float energy, float charge)
std::vector< CLHEP::Hep3Vector > predicted2Hits()
T x() const
Definition: PV3DBase.h:62
std::vector< std::pair< RecHitWithDist, ConstRecHitPointer > > compatibleHits(const GlobalPoint &xmeas, const GlobalPoint &vprim, float energy, float charge, const TrackerTopology *tTopo, const NavigationSchool &navigationSchool)
void setUseRecoVertex(bool val)