CMS 3D CMS Logo

ValidHitPairFilter.cc
Go to the documentation of this file.
3 
7 
9 
11 
14 
17 
20 
23 
25 
28 
30 
34 
36 
37 using namespace std;
38 
39 enum { BPix1=1, BPix2=2, BPix3=3,
42 
43 /*****************************************************************************/
44 float spin(float ph)
45 {
46  if(ph < 0) ph += 2 * M_PI;
47  return ph;
48 }
49 
50 /*****************************************************************************/
52  //Retrieve tracker topology from geometry
54  es.get<TrackerTopologyRcd>().get(tTopoHand);
55  tTopo = tTopoHand.product();
56 
57 
58  // Get tracker
59  edm::ESHandle<TrackerGeometry> trackerHandle;
60  es.get<TrackerDigiGeometryRecord>().get(trackerHandle);
61  theTracker = trackerHandle.product();
62 
63  // Get geometric search tracker
64  edm::ESHandle<GeometricSearchTracker> geometricSearchTrackerHandle;
65  es.get<TrackerRecoGeometryRecord>().get(geometricSearchTrackerHandle);
66  theGSTracker = geometricSearchTrackerHandle.product();
67 
68  // Get magnetic field
69  edm::ESHandle<MagneticField> magneticFieldHandle;
70  es.get<IdealMagneticFieldRecord>().get(magneticFieldHandle);
71  theMagneticField = magneticFieldHandle.product();
72 
73  // Get propagator
74  edm::ESHandle<Propagator> thePropagatorHandle;
75  es.get<TrackingComponentsRecord>().get("AnalyticalPropagator",
76  thePropagatorHandle);
77  thePropagator = thePropagatorHandle.product();
78 
79  // Bounds, hardwired FIXME
80  rzBounds[0].resize(8); phBounds[0].resize(20);
81  rzBounds[1].resize(8); phBounds[1].resize(32);
82  rzBounds[2].resize(8); phBounds[2].resize(44);
83 
84  rzBounds[3].resize(7); phBounds[3].resize(24);
85  rzBounds[4].resize(7); phBounds[4].resize(24);
86  rzBounds[5].resize(7); phBounds[5].resize(24);
87  rzBounds[6].resize(7); phBounds[6].resize(24);
88 
89  LogTrace("MinBiasTracking")
90  << " [ValidHitPairFilter] initializing pixel barrel";
91 
92  for(TrackerGeometry::DetContainer::const_iterator
93  det = theTracker->detsPXB().begin();
94  det!= theTracker->detsPXB().end(); det++)
95  {
96 
97  DetId pid=(*det)->geographicalId();
98  int il = tTopo->pxbLayer(pid) - 1;
99  int irz = tTopo->pxbModule(pid) - 1;
100  int iph = tTopo->pxbLadder(pid) - 1;
101 
102  rzBounds[il][irz] = (*det)->position().z();
103  phBounds[il][iph] = spin((*det)->position().phi());
104  }
105 
106  LogTrace("MinBiasTracking")
107  << " [ValidHitPairFilter] initializing pixel endcap";
108 
109  for(TrackerGeometry::DetContainer::const_iterator
110  det = theTracker->detsPXF().begin();
111  det!= theTracker->detsPXF().end(); det++)
112  {
113 
114 
115  DetId pid=(*det)->geographicalId();
116  int il = BPix3 + ((tTopo->pxfSide(pid) -1) << 1) + (tTopo->pxfDisk(pid) -1);
117  int irz = ((tTopo->pxfModule(pid)-1) << 1) + (tTopo->pxfPanel(pid)-1);
118  int iph = (tTopo->pxfBlade(pid) -1);
119 
120  rzBounds[il][irz] = (*det)->position().perp();
121  phBounds[il][iph] = spin((*det)->position().phi());
122  }
123 
124  for(int i = 0; i < 7; i++)
125  {
126  sort(rzBounds[i].begin(), rzBounds[i].end());
127  sort(phBounds[i].begin(), phBounds[i].end());
128  }
129 }
130 
131 /*****************************************************************************/
133 {
134 }
135 
136 /*****************************************************************************/
138 {
139  DetId id(recHit.geographicalId());
140 
141  if(id.subdetId() == int(PixelSubdetector::PixelBarrel))
142  {
143 
144  return tTopo->pxbLayer(id);
145  }
146  else
147  {
148 
149  return BPix3 + ((tTopo->pxfSide(id)-1) << 1) + tTopo->pxfDisk(id);
150  }
151 }
152 
153 /*****************************************************************************/
154 vector<int> ValidHitPairFilter::getMissingLayers(int a, int b) const
155 {
156  vector<int> l;
157  pair<int,int> c(a,b);
158 
159  if(c == pair<int,int>(BPix1,BPix2)) { l.push_back(int(BPix3));
160  l.push_back(int(FPix1_pos));
161  l.push_back(int(FPix1_neg)); return l; }
162  if(c == pair<int,int>(BPix1,BPix3)) { l.push_back(int(BPix2)); return l; }
163  if(c == pair<int,int>(BPix2,BPix3)) { l.push_back(int(BPix1)); return l; }
164  if(c == pair<int,int>(BPix1,FPix1_pos)) { l.push_back(int(BPix2));
165  l.push_back(int(FPix2_pos)); return l; }
166  if(c == pair<int,int>(BPix1,FPix1_neg)) { l.push_back(int(BPix2));
167  l.push_back(int(FPix2_neg)); return l; }
168  if(c == pair<int,int>(BPix1,FPix2_pos)) { l.push_back(int(BPix2));
169  l.push_back(int(FPix1_pos)); return l; }
170  if(c == pair<int,int>(BPix1,FPix2_neg)) { l.push_back(int(BPix2));
171  l.push_back(int(FPix1_neg)); return l; }
172  if(c == pair<int,int>(BPix2,FPix1_pos)) { l.push_back(int(BPix1));
173  l.push_back(int(FPix2_pos)); return l; }
174  if(c == pair<int,int>(BPix2,FPix1_neg)) { l.push_back(int(BPix1));
175  l.push_back(int(FPix2_neg)); return l; }
176  if(c == pair<int,int>(BPix2,FPix2_pos)) { l.push_back(int(BPix1));
177  l.push_back(int(FPix1_pos)); return l; }
178  if(c == pair<int,int>(BPix2,FPix2_neg)) { l.push_back(int(BPix1));
179  l.push_back(int(FPix1_neg)); return l; }
180  if(c == pair<int,int>(FPix1_pos,FPix2_pos)) { l.push_back(int(BPix1)); return l; }
181  if(c == pair<int,int>(FPix1_neg,FPix2_neg)) { l.push_back(int(BPix1)); return l; }
182 
183  return l;
184 }
185 
186 /*****************************************************************************/
188  (const reco::Track & track) const
189 {
190  GlobalPoint position(track.vertex().x(),
191  track.vertex().y(),
192  track.vertex().z());
193 
194  GlobalVector momentum(track.momentum().x(),
195  track.momentum().y(),
196  track.momentum().z());
197 
198  GlobalTrajectoryParameters gtp(position, momentum,
199  track.charge(), theMagneticField);
200 
201  FreeTrajectoryState fts(gtp);
202 
203  return fts;
204 }
205 
206 /*****************************************************************************/
207 vector<const GeomDet *> ValidHitPairFilter::getCloseDets
208  (int il,
209  float rz, const vector<float>& rzB,
210  float ph, const vector<float>& phB) const
211 {
212  vector<int> rzVec, phVec;
213 
214  // Radius or z
215  vector<float>::const_iterator irz = lower_bound(rzB.begin(),rzB.end(), rz);
216  if(irz > rzB.begin()) rzVec.push_back(irz - rzB.begin() - 1);
217  if(irz < rzB.end() ) rzVec.push_back(irz - rzB.begin() );
218 
219  // Phi, circular
220  vector<float>::const_iterator iph = lower_bound(phB.begin(),phB.end(), ph);
221  if(iph > phB.begin()) phVec.push_back(iph - phB.begin() - 1);
222  else phVec.push_back(phB.end() - phB.begin() - 1);
223  if(iph < phB.end() ) phVec.push_back(iph - phB.begin() );
224  else phVec.push_back(phB.begin() - phB.begin() );
225 
226  // Detectors
227  vector<const GeomDet *> dets;
228 
229  for(vector<int>::const_iterator irz = rzVec.begin(); irz!= rzVec.end(); irz++)
230  for(vector<int>::const_iterator iph = phVec.begin(); iph!= phVec.end(); iph++)
231  {
232  if(il+1 <= BPix3)
233  {
234  int layer = il + 1;
235  int ladder = *iph + 1;
236  int module = *irz + 1;
237 
238  LogTrace("MinBiasTracking")
239  << " [ValidHitPairFilter] third ("<<layer<< "|"<<ladder<<"|"<<module<<")";
240 
241  DetId id=tTopo->pxbDetId(layer,ladder,module);
242  dets.push_back(theTracker->idToDet(id));
243  }
244  else
245  {
246  int side = (il - BPix3) / 2 +1;
247  int disk = (il - BPix3) % 2 + 1;
248  int blade = *iph + 1;
249  int panel = (*irz) % 2 + 1;
250  int module = (*irz) / 2 + 1;
251 
252  LogTrace("MinBiasTracking")
253  << " [ValidHitPairFilter] third ("<<side<<"|"<<disk<<"|"<<blade<<"|"<<panel<<"|"<<module<<")";
254 
255 
256  DetId id=tTopo->pxfDetId(side,disk,blade,panel,module);
257  dets.push_back(theTracker->idToDet(id));
258 
259  }
260  }
261 
262  return dets;
263 }
264 
265 /*****************************************************************************/
267  (const reco::Track * track, const vector<const TrackingRecHit *>& recHits) const
268 {
269  bool hasGap = true;
270 
271  if(recHits.size() == 2)
272  {
273  LogTrace("MinBiasTracking")
274  << " [ValidHitPairFilter] pair" << HitInfo::getInfo(*(recHits[0]),tTopo)
275  << HitInfo::getInfo(*(recHits[1]),tTopo);
276 
277  float tol = 0.1; // cm
278  float sc = -1.; // scale, allow 'tol' of edge to count as outside
279  LocalError le(tol*tol, tol*tol, tol*tol);
280 
281  // determine missing layers
282  vector<int> missingLayers = getMissingLayers(getLayer(*(recHits[0])),
283  getLayer(*(recHits[1])));
284 
285  for(vector<int>::const_iterator missingLayer = missingLayers.begin();
286  missingLayer!= missingLayers.end();
287  missingLayer++)
288  {
289  int il = *missingLayer - 1;
290 
291  // propagate to missing layer
292  FreeTrajectoryState fts = getTrajectory(*track);
294  float rz = 0.;
295 
296  if(il < BPix3)
297  { // barrel
298  const BarrelDetLayer * layer =
299  (theGSTracker->pixelBarrelLayers())[il];
300 
301  tsos = thePropagator->propagate(fts, layer->surface());
302 
303  if(tsos.isValid())
304  rz = tsos.globalPosition().z();
305  }
306  else
307  { // endcap
308  const ForwardDetLayer * layer;
309  if(il - BPix3 < 2)
310  layer = (theGSTracker->negPixelForwardLayers())[il - BPix3 ];
311  else
312  layer = (theGSTracker->posPixelForwardLayers())[il - BPix3 - 2];
313 
314  tsos = thePropagator->propagate(fts, layer->surface());
315 
316  if(tsos.isValid())
317  rz = tsos.globalPosition().perp();
318  }
319 
320  if(tsos.isValid())
321  {
322  float phi = spin(tsos.globalPosition().phi());
323 
324  // check close dets
325  vector<const GeomDet *> closeDets =
326  getCloseDets(il, rz ,rzBounds[il], phi,phBounds[il]);
327 
328  for(vector<const GeomDet *>::const_iterator det = closeDets.begin();
329  det!= closeDets.end();
330  det++)
331  {
333  thePropagator->propagate(fts, (*det)->surface());
334 
335  if(tsos.isValid())
336  if((*det)->surface().bounds().inside(tsos.localPosition(), le, sc))
337  hasGap = false;
338  }
339  }
340  }
341  }
342 
343 
344  if(hasGap)
345  LogTrace("MinBiasTracking") << " [ValidHitPairFilter] has gap --> good";
346  else
347  LogTrace("MinBiasTracking") << " [ValidHitPairFilter] no gap --> rejected";
348 
349  return hasGap;
350 }
351 
T perp() const
Definition: PV3DBase.h:72
int getLayer(const TrackingRecHit &recHit) const
ValidHitPairFilter(const edm::EventSetup &es)
Geom::Phi< T > phi() const
Definition: PV3DBase.h:69
GlobalPoint globalPosition() const
const Vector & momentum() const
track momentum vector
Definition: TrackBase.h:714
const Point & vertex() const
reference point on the track. This method is DEPRECATED, please use referencePoint() instead ...
Definition: TrackBase.h:738
T z() const
Definition: PV3DBase.h:64
#define end
Definition: vmac.h:39
std::vector< int > getMissingLayers(int a, int b) const
#define LogTrace(id)
#define M_PI
Definition: DetId.h:18
double b
Definition: hdecay.h:120
#define begin
Definition: vmac.h:32
const BoundSurface & surface() const final
The surface of the GeometricSearchDet.
double a
Definition: hdecay.h:121
static int position[264][3]
Definition: ReadPGInfo.cc:509
T get() const
Definition: EventSetup.h:71
int charge() const
track electric charge
Definition: TrackBase.h:606
DetId geographicalId() const
static std::string getInfo(const DetId &id, const TrackerTopology *tTopo)
Definition: HitInfo.cc:23
T const * product() const
Definition: ESHandle.h:86
Definition: vlib.h:208
const BoundSurface & surface() const final
GeometricSearchDet interface.
float spin(float ph)
FreeTrajectoryState getTrajectory(const reco::Track &track) const
std::vector< const GeomDet * > getCloseDets(int il, float rz, const std::vector< float > &rzB, float ph, const std::vector< float > &phB) const