CMS 3D CMS Logo

TIBRing.cc
Go to the documentation of this file.
1 #include "TIBRing.h"
2 
4 
12 
13 #include "LayerCrossingSide.h"
14 #include "DetGroupMerger.h"
16 
17 using namespace std;
18 
20 
21 TIBRing::TIBRing(vector<const GeomDet*>& theGeomDets):
23  theDets(theGeomDets.begin(),theGeomDets.end())
24 {
25  //checkRadius( first, last);
26  //sort( theDets.begin(), theDets.end(), DetLessPhi());
27  //checkPeriodicity( theDets.begin(), theDets.end());
28 
29  theBinFinder = BinFinderType( theDets.front()->surface().position().phi(),
30  theDets.size());
31 
33 
35 
36 
37  LogDebug("TkDetLayers") << "==== DEBUG TIBRing =====" ;
38  LogDebug("TkDetLayers") << "radius, thickness, lenght: "
39  << theCylinder->radius() << " , "
40  << theCylinder->bounds().thickness() << " , "
41  << theCylinder->bounds().length() ;
42 
43  for (vector<const GeomDet*>::const_iterator i=theDets.begin();
44  i != theDets.end(); i++){
45  LogDebug("TkDetLayers") << "Ring's Det pos z,perp,eta,phi: "
46  << (**i).position().z() << " , "
47  << (**i).position().perp() << " , "
48  << (**i).position().eta() << " , "
49  << (**i).position().phi() ;
50  }
51  LogDebug("TkDetLayers") << "==== end DEBUG TIBRing =====" ;
52 
53 
54 }
55 
56 
57 const vector<const GeometricSearchDet*>&
59 {
60  throw DetLayerException("TIBRing doesn't have GeometricSearchDet components");
61 }
62 
63 void TIBRing::checkRadius(vector<const GeomDet*>::const_iterator first,
64  vector<const GeomDet*>::const_iterator last)
65 {
66  // check radius range
67  float rMin = 10000.;
68  float rMax = 0.;
69  for (vector<const GeomDet*>::const_iterator i=first; i!=last; i++) {
70  float r = (**i).surface().position().perp();
71  if (r < rMin) rMin = r;
72  if (r > rMax) rMax = r;
73  }
74  if (rMax - rMin > 0.1) throw DetLayerException(
75  "TIBRing construction failed: detectors not at constant radius");
76 }
77 
78 
79 void TIBRing::checkPeriodicity(vector<const GeomDet*>::const_iterator first,
80  vector<const GeomDet*>::const_iterator last)
81 {
82  vector<double> adj_diff(last-first-1);
83  for (int i=0; i < static_cast<int>(adj_diff.size()); i++) {
84  vector<const GeomDet*>::const_iterator curent = first + i;
85  adj_diff[i] = (**(curent+1)).surface().position().phi() -
86  (**curent).surface().position().phi();
87  }
88  double step = stat_mean( adj_diff);
89  double phi_step = 2.*Geom::pi()/(last-first);
90 
91  if ( std::abs(step-phi_step)/phi_step > 0.01) {
92  int ndets = last-first;
93  edm::LogError("TkDetLayers") << "TIBRing Warning: not periodic. ndets=" << ndets ;
94  for (int j=0; j<ndets; j++) {
95  edm::LogError("TkDetLayers") << "Dets(r,phi): (" << theDets[j]->surface().position().perp()
96  << "," << theDets[j]->surface().position().phi() << ") " ;
97  }
98  throw DetLayerException( "Error: TIBRing is not periodic");
99  }
100 }
101 
103 
104  const GeomDet& det = *theDets.front();
105  GlobalVector radial = det.surface().position() - GlobalPoint(0,0,0);
106  GlobalVector normal = det.surface().toGlobal( LocalVector(0,0,1));
107  if(normal.dot(radial)<=0)normal*=-1;
108 // edm::LogInfo(TkDetLayers) << "BarrelDetRing::computeHelicity: phi(normal) " << normal.phi()
109 // << " phi(radial) " << radial.phi() ;
110  if (PhiLess()( normal.phi(), radial.phi())) {
111  theHelicity = 1; // smaller phi angles mean "inner" group
112  }
113  else {
114  theHelicity = 0; // smaller phi angles mean "outer" group
115  }
116 }
117 
118 
120 
121 }
122 
123 
124 pair<bool, TrajectoryStateOnSurface>
126  const MeasurementEstimator&) const{
127  edm::LogError("TkDetLayers") << "temporary dummy implementation of TIBRing::compatible()!!" ;
128  return pair<bool,TrajectoryStateOnSurface>();
129 }
130 
131 
132 void
134  const Propagator& prop,
135  const MeasurementEstimator& est,
136  vector<DetGroup> & result) const
137 {
138  vector<DetGroup> closestResult;
139  SubRingCrossings crossings;
140  crossings = computeCrossings( tsos, prop.propagationDirection());
141  if(! crossings.isValid_) return;
142 
143  typedef CompatibleDetToGroupAdder Adder;
145  tsos, prop, est, closestResult);
146 
147  if(closestResult.empty()){
149  tsos, prop, est, result);
150  return;
151  }
152 
153  DetGroupElement closestGel( closestResult.front().front());
154  float window = computeWindowSize( closestGel.det(), closestGel.trajectoryState(), est);
155 
156  float detWidth = closestGel.det()->surface().bounds().width();
157  if (crossings.nextDistance < detWidth + window) {
158  vector<DetGroup> nextResult;
159  if (Adder::add( *theDets[theBinFinder.binIndex(crossings.nextIndex)],
160  tsos, prop, est, nextResult)) {
161  int crossingSide = LayerCrossingSide::barrelSide( tsos, prop);
162  if (crossings.closestIndex < crossings.nextIndex) {
163  DetGroupMerger::orderAndMergeTwoLevels( std::move(closestResult), std::move(nextResult),
164  result,
165  theHelicity, crossingSide);
166  }
167  else {
168  DetGroupMerger::orderAndMergeTwoLevels( std::move(nextResult), std::move(closestResult),
169  result,
170  theHelicity, crossingSide);
171  }
172  }
173  else {
174  result.swap(closestResult);
175  }
176  }else{
177  result.swap(closestResult);
178  }
179 
180  // only loop over neighbors (other than closest and next) if window is BIG
181  if (window > 0.5f*detWidth) {
182  searchNeighbors( tsos, prop, est, crossings, window, result);
183  }
184 }
185 
187  const Propagator& prop,
188  const MeasurementEstimator& est,
189  const SubRingCrossings& crossings,
190  float window,
191  vector<DetGroup>& result) const
192 {
193  typedef CompatibleDetToGroupAdder Adder;
194  int crossingSide = LayerCrossingSide::barrelSide( tsos, prop);
195  typedef DetGroupMerger Merger;
196 
197  int negStart = std::min( crossings.closestIndex, crossings.nextIndex) - 1;
198  int posStart = std::max( crossings.closestIndex, crossings.nextIndex) + 1;
199 
200  int quarter = theDets.size()/4;
201  for (int idet=negStart; idet >= negStart - quarter+1; idet--) {
202  const GeomDet* neighbor = theDets[theBinFinder.binIndex(idet)];
203  // if (!overlap( gCrossingPos, *neighbor, window)) break; // mybe not needed?
204  // maybe also add shallow crossing angle test here???
205  vector<DetGroup> tmp1;
206  if (!Adder::add( *neighbor, tsos, prop, est, tmp1)) break;
207  vector<DetGroup> tmp2; tmp2.swap(result);
208  vector<DetGroup> newResult;
209  Merger::orderAndMergeTwoLevels(std::move(tmp1), std::move(tmp2), newResult, theHelicity, crossingSide);
210  result.swap(newResult);
211  }
212  for (int idet=posStart; idet < posStart + quarter-1; idet++) {
213  const GeomDet* neighbor = theDets[theBinFinder.binIndex(idet)];
214  // if (!overlap( gCrossingPos, *neighbor, window)) break; // mybe not needed?
215  // maybe also add shallow crossing angle test here???
216  vector<DetGroup> tmp1;
217  if (!Adder::add( *neighbor, tsos, prop, est, tmp1)) break;
218  vector<DetGroup> tmp2; tmp2.swap(result);
219  vector<DetGroup> newResult;
220  Merger::orderAndMergeTwoLevels(std::move(tmp2), std::move(tmp1), newResult, theHelicity, crossingSide);
221  result.swap(newResult);
222  }
223 }
224 
225 
228  PropagationDirection propDir) const
229 {
230  typedef HelixBarrelPlaneCrossing2OrderLocal Crossing;
232 
233  GlobalPoint startPos( startingState.globalPosition());
234  GlobalVector startDir( startingState.globalMomentum());
235  auto rho = startingState.transverseCurvature();
236 
237  HelixBarrelCylinderCrossing cylCrossing( startPos, startDir, rho,
238  propDir,specificSurface(),
240 
241  if (!cylCrossing.hasSolution()) return SubRingCrossings();
242 
243  GlobalPoint cylPoint( cylCrossing.position());
244  GlobalVector cylDir( cylCrossing.direction());
245  int closestIndex = theBinFinder.binIndex(cylPoint.barePhi());
246 
247  const Plane& closestPlane( theDets[closestIndex]->surface());
248 
249  LocalPoint closestPos = Crossing::positionOnly( cylPoint, cylDir, rho, closestPlane);
250  float closestDist = closestPos.x(); // use fact that local X perp to global Z
251 
252  //int next = cylPoint.phi() - closestPlane.position().phi() > 0 ? closest+1 : closest-1;
253  int nextIndex = PhiLess()( closestPlane.position().barePhi(), cylPoint.barePhi()) ?
254  closestIndex+1 : closestIndex-1;
255 
256  const Plane& nextPlane( theDets[ theBinFinder.binIndex(nextIndex)]->surface());
257  LocalPoint nextPos = Crossing::positionOnly( cylPoint, cylDir, rho, nextPlane);
258  float nextDist = nextPos.x();
259 
260  if (std::abs(closestDist) < std::abs(nextDist)) {
261  return SubRingCrossings( closestIndex, nextIndex, nextDist);
262  }
263  else {
264  return SubRingCrossings( nextIndex, closestIndex, closestDist);
265  }
266 }
267 
269  const TrajectoryStateOnSurface& tsos,
270  const MeasurementEstimator& est) const
271 {
272  return est.maximalLocalDisplacement(tsos, det->surface()).x();
273 }
274 
275 
276 
#define LogDebug(id)
GlobalPoint toGlobal(const Point2DBase< Scalar, LocalTag > lp) const
Definition: Surface.h:106
Common base class.
Local3DVector LocalVector
Definition: LocalVector.h:12
~TIBRing() __attribute__((cold))
Definition: TIBRing.cc:119
static int barrelSide(const TrajectoryStateOnSurface &startingState, const Propagator &prop)
returns 0 if barrel layer crossed from inside, 1 if from outside
BinFinderType theBinFinder
Definition: TIBRing.h:87
virtual const BoundSurface & surface() const
The surface of the GeometricSearchDet.
Definition: TIBRing.h:19
TIBRing(std::vector< const GeomDet * > &theGeomDets) __attribute__((cold))
Definition: TIBRing.cc:21
void computeHelicity() __attribute__((cold))
Definition: TIBRing.cc:102
virtual std::pair< bool, TrajectoryStateOnSurface > compatible(const TrajectoryStateOnSurface &ts, const Propagator &, const MeasurementEstimator &) const __attribute__((cold))
Definition: TIBRing.cc:125
GeometricSearchDet::DetWithState DetWithState
Definition: TIBRing.cc:19
Geom::Phi< T > phi() const
Definition: PV3DBase.h:69
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
std::vector< const GeomDet * > theDets
Definition: TIBRing.h:89
GlobalPoint globalPosition() const
PreciseFloatType< T, U >::Type dot(const Vector3DBase< U, FrameTag > &v) const
Definition: Vector3DBase.h:107
Vector2DBase< float, LocalTag > Local2DVector
PropagationDirection
bool neighbor(int endcap, int sector, int SectIndex, int id, int sub, int station)
virtual const std::vector< const GeometricSearchDet * > & components() const __attribute__((cold))
Returns basic components, if any.
Definition: TIBRing.cc:58
double stat_mean(const CONT &cont)
Definition: simple_stat.h:9
const Plane & surface() const
The nominal surface of the GeomDet.
Definition: GeomDet.h:42
Definition: Plane.h:17
void checkRadius(std::vector< const GeomDet * >::const_iterator first, std::vector< const GeomDet * >::const_iterator last) __attribute__((cold))
Definition: TIBRing.cc:63
virtual int binIndex(T phi) const
returns an index in the valid range for the bin that contains phi
virtual void groupedCompatibleDetsV(const TrajectoryStateOnSurface &startingState, const Propagator &prop, const MeasurementEstimator &est, std::vector< DetGroup > &result) const __attribute__((hot))
Definition: TIBRing.cc:133
ReferenceCountingPointer< BoundCylinder > theCylinder
Definition: TIBRing.h:90
T barePhi() const
Definition: PV3DBase.h:68
int theHelicity
Definition: TIBRing.h:91
virtual PropagationDirection propagationDirection() const final
Definition: Propagator.h:151
float computeWindowSize(const GeomDet *det, const TrajectoryStateOnSurface &tsos, const MeasurementEstimator &est) const __attribute__((hot))
Definition: TIBRing.cc:268
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
double f[11][100]
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
#define end
Definition: vmac.h:37
T min(T a, T b)
Definition: MathUtil.h:58
void searchNeighbors(const TrajectoryStateOnSurface &tsos, const Propagator &prop, const MeasurementEstimator &est, const SubRingCrossings &crossings, float window, std::vector< DetGroup > &result) const __attribute__((hot))
Definition: TIBRing.cc:186
void add(std::map< std::string, TH1 * > &h, TH1 *hist)
PeriodicBinFinderInPhi< float > BinFinderType
Definition: TIBRing.h:86
virtual const BoundCylinder & specificSurface() const
Return the ring surface as a.
Definition: TIBRing.h:41
SubRingCrossings computeCrossings(const TrajectoryStateOnSurface &startingState, PropagationDirection propDir) const __attribute__((hot))
Definition: TIBRing.cc:227
#define begin
Definition: vmac.h:30
GlobalVector globalMomentum() const
virtual Local2DVector maximalLocalDisplacement(const TrajectoryStateOnSurface &ts, const Plane &plane) const =0
void checkPeriodicity(std::vector< const GeomDet * >::const_iterator first, std::vector< const GeomDet * >::const_iterator last) __attribute__((cold))
Definition: TIBRing.cc:79
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)
step
constexpr double pi()
Definition: Pi.h:31
T x() const
Definition: PV3DBase.h:62
const PositionType & position() const
def move(src, dest)
Definition: eostools.py:510