CMS 3D CMS Logo

Phase1PixelBlade.cc
Go to the documentation of this file.
1 #include "Phase1PixelBlade.h"
2 
4 
6 #include "LayerCrossingSide.h"
7 #include "DetGroupMerger.h"
10 
14 
15 using namespace std;
16 
18 
20 
21 Phase1PixelBlade::Phase1PixelBlade(vector<const GeomDet*>& frontDets,
22  vector<const GeomDet*>& backDets):
24  theFrontDets(frontDets),
25  theBackDets(backDets),
26  front_radius_range_(std::make_pair(0,0)),
27  back_radius_range_(std::make_pair(0,0))
28 {
29  theDets.assign(theFrontDets.begin(),theFrontDets.end());
30  theDets.insert(theDets.end(),theBackDets.begin(),theBackDets.end());
31 
35 
36  //--------- DEBUG INFO --------------
37  LogDebug("TkDetLayers") << "DEBUG INFO for Phase1PixelBlade" ;
38  LogDebug("TkDetLayers") << "Blade z, perp, innerRadius, outerR[disk, front, back]: "
39  << this->position().z() << " , "
40  << this->position().perp() << " , ("
41  << theDiskSector->innerRadius() << " , "
42  << theDiskSector->outerRadius() << "), ("
43  << theFrontDiskSector->innerRadius() << " , "
44  << theFrontDiskSector->outerRadius() << "), ("
45  << theBackDiskSector->innerRadius() << " , "
46  << theBackDiskSector->outerRadius() << ")" << std::endl;
47 
48  for(vector<const GeomDet*>::const_iterator it=theFrontDets.begin();
49  it!=theFrontDets.end(); it++){
50  LogDebug("TkDetLayers") << "frontDet phi,z,r: "
51  << (*it)->position().phi() << " , "
52  << (*it)->position().z() << " , "
53  << (*it)->position().perp() << " , "
54  << " rmin: " << (*it)->surface().rSpan().first
55  << " rmax: " << (*it)->surface().rSpan().second
56  << std::endl;
57  }
58 
59  for(vector<const GeomDet*>::const_iterator it=theBackDets.begin();
60  it!=theBackDets.end(); it++){
61  LogDebug("TkDetLayers") << "backDet phi,z,r: "
62  << (*it)->position().phi() << " , "
63  << (*it)->position().z() << " , "
64  << (*it)->position().perp() << " , "
65  << " rmin: " << (*it)->surface().rSpan().first
66  << " rmax: " << (*it)->surface().rSpan().second
67  << std::endl;
68  }
69 
70 }
71 
72 
73 const vector<const GeometricSearchDet*>&
75  throw DetLayerException("TOBRod doesn't have GeometricSearchDet components");
76 }
77 
78 pair<bool, TrajectoryStateOnSurface>
80  const MeasurementEstimator&) const{
81  edm::LogError("TkDetLayers") << "temporary dummy implementation of Phase1PixelBlade::compatible()!!" ;
82  return pair<bool,TrajectoryStateOnSurface>();
83 }
84 
85 
86 
87 void
89  const Propagator& prop,
90  const MeasurementEstimator& est,
91  std::vector<DetGroup> & result) const{
92 
93 
94  SubLayerCrossings crossings;
95 
96  crossings = computeCrossings( tsos, prop.propagationDirection());
97  if(! crossings.isValid()) return;
98 
99  vector<DetGroup> closestResult;
100  addClosest( tsos, prop, est, crossings.closest(), closestResult);
101 
102  if (closestResult.empty()) {
103  vector<DetGroup> nextResult;
104  addClosest( tsos, prop, est, crossings.other(), nextResult);
105  if(nextResult.empty()) return;
106 
107  //DetGroupElement nextGel( nextResult.front().front());
108  //int crossingSide = LayerCrossingSide().endcapSide( nextGel.trajectoryState(), prop);
109 
110  DetGroupMerger::orderAndMergeTwoLevels( std::move(closestResult), std::move(nextResult), result,
111  0,0);//fixme gc patched for SLHC - already correctly sorted for SLHC
112  //crossings.closestIndex(), crossingSide);
113  }
114  else {
115  DetGroupElement closestGel( closestResult.front().front());
116  float window = computeWindowSize( closestGel.det(), closestGel.trajectoryState(), est);
117 
118  searchNeighbors( tsos, prop, est, crossings.closest(), window,
119  closestResult, false);
120 
121  vector<DetGroup> nextResult;
122  searchNeighbors( tsos, prop, est, crossings.other(), window,
123  nextResult, true);
124 
125  //int crossingSide = LayerCrossingSide().endcapSide( closestGel.trajectoryState(), prop);
126  DetGroupMerger::orderAndMergeTwoLevels( std::move(closestResult), std::move(nextResult), result,
127  0,0);//fixme gc patched for SLHC - already correctly sorted for SLHC
128  //crossings.closestIndex(), crossingSide);
129  }
130 }
131 
134  PropagationDirection propDir) const
135 {
136  HelixPlaneCrossing::PositionType startPos( startingState.globalPosition());
137  HelixPlaneCrossing::DirectionType startDir( startingState.globalMomentum());
138  double rho( startingState.transverseCurvature());
139  HelixArbitraryPlaneCrossing crossing( startPos, startDir, rho, propDir);
140 
141  pair<bool,double> innerPath = crossing.pathLength( *theFrontDiskSector);
142  if (!innerPath.first) return SubLayerCrossings();
143 
144  GlobalPoint gInnerPoint( crossing.position(innerPath.second));
145  //Code for use of binfinder
146  //int innerIndex = theInnerBinFinder.binIndex(gInnerPoint.perp());
147  //float innerDist = std::abs( theInnerBinFinder.binPosition(innerIndex) - gInnerPoint.z());
148 
149  // int innerIndex = findBin(gInnerPoint.perp(),0);
150  int innerIndex = findBin2(gInnerPoint,0);
151 
152  //fixme gc patched for SLHC - force order here to be in z
153  //float innerDist = std::abs( findPosition(innerIndex,0).perp() - gInnerPoint.perp());
154  float innerDist = ( startingState.globalPosition() - gInnerPoint).mag();
155  SubLayerCrossing innerSLC( 0, innerIndex, gInnerPoint);
156 
157  pair<bool,double> outerPath = crossing.pathLength( *theBackDiskSector);
158  if (!outerPath.first) return SubLayerCrossings();
159 
160  GlobalPoint gOuterPoint( crossing.position(outerPath.second));
161  //Code for use of binfinder
162  //int outerIndex = theOuterBinFinder.binIndex(gOuterPoint.perp());
163  //float outerDist = std::abs( theOuterBinFinder.binPosition(outerIndex) - gOuterPoint.perp());
164  // int outerIndex = findBin(gOuterPoint.perp(),1);
165  int outerIndex = findBin2(gOuterPoint,1);
166 
167  //fixme gc patched for SLHC - force order here to be in z
168  //float outerDist = std::abs( findPosition(outerIndex,1).perp() - gOuterPoint.perp());
169  float outerDist = ( startingState.globalPosition() - gOuterPoint).mag();
170  SubLayerCrossing outerSLC( 1, outerIndex, gOuterPoint);
171 
172  if (innerDist < outerDist) {
173  return SubLayerCrossings( innerSLC, outerSLC, 0);
174  }
175  else {
176  return SubLayerCrossings( outerSLC, innerSLC, 1);
177  }
178 }
179 
180 
181 
182 bool
184  const Propagator& prop,
185  const MeasurementEstimator& est,
186  const SubLayerCrossing& crossing,
187  vector<DetGroup>& result) const
188 {
189 
190  const vector<const GeomDet*>& sBlade( subBlade( crossing.subLayerIndex()));
191 
192  return CompatibleDetToGroupAdder().add( *sBlade[crossing.closestDetIndex()],
193  tsos, prop, est, result);
194 }
195 
196 
198  const TrajectoryStateOnSurface& tsos,
199  const MeasurementEstimator& est) const
200 {
201  return
202  est.maximalLocalDisplacement(tsos, det->surface()).x();
203 }
204 
205 
206 
207 
209  const Propagator& prop,
210  const MeasurementEstimator& est,
211  const SubLayerCrossing& crossing,
212  float window,
213  vector<DetGroup>& result,
214  bool checkClosest) const
215 {
216  GlobalPoint gCrossingPos = crossing.position();
217  const vector<const GeomDet*>& sBlade( subBlade( crossing.subLayerIndex()));
218 
219  int closestIndex = crossing.closestDetIndex();
220  int negStartIndex = closestIndex-1;
221  int posStartIndex = closestIndex+1;
222 
223  if (checkClosest) { // must decide if the closest is on the neg or pos side
224  if (gCrossingPos.perp() < sBlade[closestIndex]->surface().position().perp()) {
225  posStartIndex = closestIndex;
226  }
227  else {
228  negStartIndex = closestIndex;
229  }
230  }
231 
232  typedef CompatibleDetToGroupAdder Adder;
233  for (int idet=negStartIndex; idet >= 0; idet--) {
234  if (!overlap( gCrossingPos, *sBlade[idet], window)) break;
235  if (!Adder::add( *sBlade[idet], tsos, prop, est, result)) break;
236  }
237  for (int idet=posStartIndex; idet < static_cast<int>(sBlade.size()); idet++) {
238  if (!overlap( gCrossingPos, *sBlade[idet], window)) break;
239  if (!Adder::add( *sBlade[idet], tsos, prop, est, result)) break;
240  }
241 }
242 
243 
244 
245 bool Phase1PixelBlade::overlap( const GlobalPoint& crossPoint, const GeomDet& det, float window) const
246 {
247  // check if the z window around TSOS overlaps with the detector theDet (with a 1% margin added)
248 
249  // const float tolerance = 0.1;
250  const float relativeMargin = 1.01;
251 
252  LocalPoint localCrossPoint( det.surface().toLocal(crossPoint));
253  // if (std::abs(localCrossPoint.z()) > tolerance) {
254  // edm::LogInfo(TkDetLayers) << "Phase1PixelBlade::overlap calculation assumes point on surface, but it is off by "
255  // << localCrossPoint.z() ;
256  // }
257 
258  float localX = localCrossPoint.x();
259  float detHalfLength = det.surface().bounds().length()/2.;
260 
261  // edm::LogInfo(TkDetLayers) << "Phase1PixelBlade::overlap: Det at " << det.position() << " hit at " << localY
262  // << " Window " << window << " halflength " << detHalfLength ;
263 
264  if ( ( std::abs(localX)-window) < relativeMargin*detHalfLength ) { // FIXME: margin hard-wired!
265  return true;
266  } else {
267  return false;
268  }
269 }
270 
271 int
272 Phase1PixelBlade::findBin( float R,int diskSectorIndex) const
273 {
274  vector<const GeomDet*> localDets = diskSectorIndex==0 ? theFrontDets : theBackDets;
275 
276 
277  int theBin = 0;
278  float rDiff = std::abs( R - localDets.front()->surface().position().perp());
279  for (vector<const GeomDet*>::const_iterator i=localDets.begin(); i !=localDets.end(); i++){
280  float testDiff = std::abs( R - (**i).surface().position().perp());
281  if ( testDiff < rDiff) {
282  rDiff = testDiff;
283  theBin = i - localDets.begin();
284  }
285  }
286  return theBin;
287 }
288 
289 int
290 Phase1PixelBlade::findBin2( GlobalPoint thispoint,int diskSectorIndex) const
291 {
292  const vector<const GeomDet*> & localDets = diskSectorIndex==0 ? theFrontDets : theBackDets;
293 
294 
295 
296  int theBin = 0;
297  float sDiff = (thispoint - localDets.front()->surface().position()).mag();
298 
299  for (vector<const GeomDet*>::const_iterator i=localDets.begin(); i !=localDets.end(); i++){
300  float testDiff = ( thispoint - (**i).surface().position()).mag();
301  if ( testDiff < sDiff) {
302  sDiff = testDiff;
303  theBin = i - localDets.begin();
304  }
305  }
306  return theBin;
307 }
308 
309 
310 
312 Phase1PixelBlade::findPosition(int index,int diskSectorType) const
313 {
314  vector<const GeomDet*> diskSector = diskSectorType == 0 ? theFrontDets : theBackDets;
315  return (diskSector[index])->surface().position();
316 }
317 
318 std::pair<float, float>
319 Phase1PixelBlade::computeRadiusRanges(const std::vector<const GeomDet*> &current_dets) {
320 
322  Vector posSum(0,0,0);
323  for (auto i : current_dets)
324  posSum += (*i).surface().position().basicVector();
325 
326  Surface::PositionType meanPos( 0.,0.,posSum.z()/float(current_dets.size()));
327 
328  // temporary plane - for the computation of bounds
329  const Plane& tmpplane = current_dets.front()->surface();
330 
334 
335  GlobalVector planeXAxis = tmpplane.toGlobal(LocalVector(1, 0, 0));
336  GlobalPoint planePosition = tmpplane.position();
337 
338  if (planePosition.x() * planeXAxis.x()
339  + planePosition.y() * planeXAxis.y() > 0.){
340  yAxis = planeXAxis;
341  } else {
342  yAxis = -planeXAxis;
343  }
344 
345  GlobalVector planeZAxis = tmpplane.toGlobal(LocalVector( 0, 0, 1));
346  if (planeZAxis.z() * planePosition.z() > 0.) {
347  zAxis = planeZAxis;
348  } else {
349  zAxis = -planeZAxis;
350  }
351 
352  xAxis = yAxis.cross(zAxis);
353 
355  Plane plane(meanPos, rotation);
356 
357  Surface::PositionType tmpPos = current_dets.front()->surface().position();
358 
359  float rmin(plane.toLocal(tmpPos).perp());
360  float rmax(plane.toLocal(tmpPos).perp());
361 
362  for (auto it : current_dets) {
363  vector<GlobalPoint> corners = BoundingBox().corners(it->specificSurface());
364  for (auto i : corners) {
365  float r = plane.toLocal(i).perp();
366  rmin = min(rmin, r);
367  rmax = max(rmax, r);
368  }
369  }
370 
371  return std::make_pair(rmin, rmax);
372 }
#define LogDebug(id)
float computeWindowSize(const GeomDet *det, const TrajectoryStateOnSurface &tsos, const MeasurementEstimator &est) const
GlobalPoint toGlobal(const Point2DBase< Scalar, LocalTag > lp) const
Definition: Surface.h:106
Common base class.
virtual float length() const =0
Local3DVector LocalVector
Definition: LocalVector.h:12
std::vector< const GeomDet * > theBackDets
std::pair< float, float > computeRadiusRanges(const std::vector< const GeomDet * > &)
T perp() const
Definition: PV3DBase.h:72
Phase1PixelBlade(std::vector< const GeomDet * > &frontDets, std::vector< const GeomDet * > &backDets) __attribute__((cold))
int closestDetIndex() const
ROOT::Math::Plane3D::Vector Vector
Definition: EcalHitMaker.cc:29
T mag() const
The vector magnitude. Equivalent to sqrt(vec.mag2())
T y() const
Definition: PV3DBase.h:63
GlobalPoint globalPosition() const
const Bounds & bounds() const
Definition: Surface.h:120
virtual const std::vector< const GeometricSearchDet * > & components() const __attribute__((cold))
Returns basic components, if any.
PropagationDirection
std::pair< bool, TrajectoryStateOnSurface > compatible(const TrajectoryStateOnSurface &ts, const Propagator &, const MeasurementEstimator &) const __attribute__((cold))
const Plane & surface() const
The nominal surface of the GeomDet.
Definition: GeomDet.h:42
Definition: Plane.h:17
const GlobalPoint & position() const
const std::vector< const GeomDet * > & subBlade(int ind) const
virtual PropagationDirection propagationDirection() const final
Definition: Propagator.h:151
bool addClosest(const TrajectoryStateOnSurface &tsos, const Propagator &prop, const MeasurementEstimator &est, const SubLayerCrossing &crossing, std::vector< DetGroup > &result) const __attribute__((hot))
LocalPoint toLocal(const GlobalPoint &gp) const
int subLayerIndex() const
Vector3DBase< typename PreciseFloatType< T, U >::Type, FrameTag > cross(const Vector3DBase< U, FrameTag > &v) const
Definition: Vector3DBase.h:119
T z() const
Definition: PV3DBase.h:64
GeometricSearchDet::DetWithState DetWithState
int findBin2(GlobalPoint thispoint, int layer) const
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
std::vector< const GeomDet * > theDets
def window(xmin, xmax, ymin, ymax, x=0, y=0, width=100, height=100, xlogbase=None, ylogbase=None, minusInfinity=-1000, flipx=False, flipy=True)
Definition: svgfig.py:642
SubLayerCrossings computeCrossings(const TrajectoryStateOnSurface &tsos, PropagationDirection propDir) const __attribute__((hot))
T min(T a, T b)
Definition: MathUtil.h:58
bool overlap(const GlobalPoint &gpos, const GeomDet &det, float phiWin) const
const SubLayerCrossing & other() const
static bool add(const GeometricSearchDet &det, const TrajectoryStateOnSurface &tsos, const Propagator &prop, const MeasurementEstimator &est, std::vector< DetGroup > &result) __attribute__((hot))
void add(std::map< std::string, TH1 * > &h, TH1 *hist)
~Phase1PixelBlade() __attribute__((cold))
static BoundDiskSector * build(const std::vector< const GeomDet * > &dets) __attribute__((cold))
virtual const Surface::PositionType & position() const
Returns position of the surface.
ReferenceCountingPointer< BoundDiskSector > theDiskSector
GlobalPoint findPosition(int index, int diskSectorIndex) const
static std::vector< GlobalPoint > corners(const Plane &)
Definition: BoundingBox.cc:24
ReferenceCountingPointer< BoundDiskSector > theFrontDiskSector
const SubLayerCrossing & closest() const
GlobalVector globalMomentum() const
int findBin(float R, int layer) const
virtual Local2DVector maximalLocalDisplacement(const TrajectoryStateOnSurface &ts, const Plane &plane) const =0
virtual void groupedCompatibleDetsV(const TrajectoryStateOnSurface &tsos, const Propagator &prop, const MeasurementEstimator &est, std::vector< DetGroup > &result) const __attribute__((hot))
TkRotation< float > RotationType
std::pair< const GeomDet *, TrajectoryStateOnSurface > DetWithState
static void orderAndMergeTwoLevels(std::vector< DetGroup > &&one, std::vector< DetGroup > &&two, std::vector< DetGroup > &result, int firstIndex, int firstCrossed)
void searchNeighbors(const TrajectoryStateOnSurface &tsos, const Propagator &prop, const MeasurementEstimator &est, const SubLayerCrossing &crossing, float window, std::vector< DetGroup > &result, bool checkClosest) const __attribute__((hot))
T x() const
Definition: PV3DBase.h:62
const PositionType & position() const
std::vector< const GeomDet * > theFrontDets
def move(src, dest)
Definition: eostools.py:510
ReferenceCountingPointer< BoundDiskSector > theBackDiskSector