GEAR  1.9.0
 All Classes Namespaces Functions Variables Typedefs Enumerations Friends Pages
ZPlanarParametersImpl.cc
1 #include "gearimpl/ZPlanarParametersImpl.h"
2 #include <math.h>
3 #include <iostream>
4 
5 #define _EPSILON 0.0001
6 
7 namespace gear{
8 
9  // F. Gaede, 04.04.2008 : phi0 and offset have not been defined properly in the past
10  // what we want is: phi0: azimuthal angle of first ladder's normal and the offset
11  // going in the direction of increasing phi
12  // however what is used in the code is phi0 as the angle w/ the y-axis (and wrong orientation)
13  //
14  // -> we keep the internal representation of phi0 in the variable internalPhi0 and use this here
15 
16 
17  ZPlanarParametersImpl::ZPlanarParametersImpl
18  ( int type, double shellInnerRadius, double shellOuterRadius, double shellHalfLength, double shellGap, double shellRadLength ) :
19  _type( type ) ,
20  _shellInnerRadius( shellInnerRadius ) ,
21  _shellOuterRadius( shellOuterRadius ) ,
22  _shellHalfLength( shellHalfLength ) ,
23  _shellGap( shellGap ) ,
24  _shellRadLength( shellRadLength ) {}
25 
26 
27  bool ZPlanarParametersImpl::isPointInLadder( Vector3D p , bool sensitive , SensorID* sensorID ) const
28  {
29 
30 
31  // very first check for quick calculation
32  if( fabs( _shellHalfLength ) > 1e-3 ) { // shell parameters can be zero if no shell !
33  if( fabs( p[2] ) > ( _shellHalfLength + _shellGap/2 ) ) {
34  // Outside shell
35  return false ;
36  }
37  if( fabs( p[2] ) < ( _shellGap/2 ) ) {
38  // Inside gap
39  return false ;
40  }
41  }
42 
43  // get radius of point in projection 2D
44  double pRadius = sqrt( p[0]*p[0] + p[1]*p[1] ) ;
45 
46  // get phi of point in projection 2D
47  double pPhi = getPhiPoint( p ) ;
48 
49  for( int possibleLayer = 0 ; possibleLayer < _layer.getNLayers() ; possibleLayer++ ) {
50 
51  // first some quick checks to reduce cpu time
52  // check if radius is smaller then this layer
53  if( ( ( !sensitive ) && ( pRadius < ( _layer.getLadderDistance( possibleLayer ) - _EPSILON ) ) )
54  ||( ( sensitive ) && ( pRadius < ( _layer.getSensitiveDistance( possibleLayer ) - _EPSILON ) ) ) ) {
55 
56 // std::cout << "Layer " << possibleLayer << " skipped - pRadius " << pRadius << " < " << std::endl ;
57  continue ;
58  }
59 
60  // check if radius is bigger than max radius for this layer
61  if( pRadius > ( _layer.getMaxRadius( possibleLayer , sensitive ) + _EPSILON ) ) {
62 
63 // std::cout << "Layer " << possibleLayer << " skipped - pRadius " << pRadius << " > " << std::endl ;
64  continue ;
65  }
66 
67  // check z-plane
68  if( ( ( !sensitive ) && ( fabs( p[2] ) > ( _layer.getLadderLength( possibleLayer ) + _shellGap/2 + _EPSILON ) ) )
69  ||( ( sensitive ) && ( fabs( p[2] ) > (_layer.getSensitiveLength( possibleLayer ) + _shellGap/2 + _EPSILON ) ) ) ) {
70 
71 // std::cout << "Layer " << possibleLayer << " skipped - z " << p[2] << " ooR " << std::endl ;
72  continue ;
73  }
74 
75  // get extensions for layer
76  double zStart, zEnd, left, right, thick, d ;
77  if( !sensitive ) {
78  zEnd = _layer.getLadderLength( possibleLayer ) + _shellGap/2 ;
79  left = -_layer.getLadderWidth( possibleLayer ) / 2 - _layer.getLadderOffset( possibleLayer ) ; //-
80  right = _layer.getLadderWidth( possibleLayer ) / 2 - _layer.getLadderOffset( possibleLayer ) ;
81  thick = _layer.getLadderThickness( possibleLayer ) ;
82  d = _layer.getLadderDistance( possibleLayer ) ;
83  }
84  if (sensitive ) {
85  zEnd = _layer.getSensitiveLength( possibleLayer ) + _shellGap/2 ;
86  left = -_layer.getSensitiveWidth( possibleLayer ) / 2 - _layer.getSensitiveOffset( possibleLayer ) ;
87  right = _layer.getSensitiveWidth( possibleLayer ) / 2 - _layer.getSensitiveOffset( possibleLayer ) ;
88  thick = _layer.getSensitiveThickness( possibleLayer ) ;
89  d = _layer.getSensitiveDistance( possibleLayer ) ;
90  }
91 
92  zStart = _shellGap/2 ;
93 
94  // find out delta phi between each ladder
95  double deltaPhi = ( 2 * M_PI ) / _layer.getNLadders( possibleLayer ) ;
96  double startPhi = _layer.getStartInnerPhi( possibleLayer , sensitive ) ;
97  double endPhi = _layer.getEndInnerPhi( possibleLayer , sensitive ) ;
98 
99  // go through all ladders/sensitive
100  int nLadders = _layer.getNLadders( possibleLayer ) ;
101  for( int i = 0 ; i < nLadders ; i++ ) {
102 
103  double phi = _layer.getInternalPhi0( possibleLayer ) + i*deltaPhi ;
104  phi = correctPhiRange( phi ) ;
105 
106  // check start and end phi for this layer
107  double sPhi = correctPhiRange( startPhi + i*deltaPhi - _EPSILON ) ;
108  double ePhi = correctPhiRange( endPhi + i*deltaPhi + _EPSILON ) ;
109 
110  // check if point is in range of ladder
111  bool isInRange =
112  ( ( sPhi <= pPhi ) && ( pPhi <= ePhi ) )
113  ||( ( sPhi > ePhi ) && ( pPhi <= ePhi ) )
114  ||( ( sPhi > ePhi ) && ( pPhi >= sPhi ) ) ;
115 
116  if ( !isInRange ) {
117  // std::cout << " ladder/layer " << i << "/" << possibleLayer << " skipped ooR" << std::endl ;
118  continue ;
119  }
120 
121  // take normal vector of planes, Base, Side and z
122  Vector3D nBase, nSide, nZ ;
123 
124  nBase[0] = sin( phi ) ;
125  nBase[1] = cos( phi ) ;
126  nBase[2] = 0;
127 
128  nSide[0] = cos( phi ) ;
129  nSide[1] = -sin( phi ) ;
130  nSide[2] = 0 ;
131 
132  nZ[0] = 0 ;
133  nZ[1] = 0 ;
134  nZ[2] = 1 ;
135 
136  if( p[2] < 0 ) {
137  nZ[2] = - nZ[2] ;
138  }
139 
140  Vector3D r( d * nBase[0] , d * nBase[1] , d * nBase[2] ) ;
141 
142  Vector3D thisVector = distanceToPlane( p, r, nBase, nSide, nZ, left, right, zStart, zEnd ) ;
143 
144 
145  // calculate distance
146  double c = sqrt( thisVector[0] * thisVector[0] + thisVector[1] * thisVector[1] + thisVector[2] * thisVector[2] ) ;
147 
148  // debug
149 // std::cout << "**n[0] :" << nBase[0] << " n[1] " <<nBase[1] <<" n[2] " <<nBase[2] << std::endl ;
150 // std::cout << "**nv[0]:" << thisVector[0] << " nv[1]:" << thisVector[1] << "nv[2]:" << thisVector[2] << std::endl ;
151 // std::cout << "**c :" << c << " thick: " << thick << std::endl ; // debug
152 // std::cout << "**sPhi : " << sPhi*180/M_PI <<" pPhi " << pPhi*180/M_PI << " ePhi: " << ePhi*180/M_PI << std::endl ;
153 // std::cout << "**Phi : " << phi*180/M_PI << " atan() : " << atan( p[0]/p[1] ) << std::endl ;
154 // std::cout << "**isEqual[0] " << isEqual( thisVector[0] , -c*nBase[0] )
155 // << " , " << isEqual( thisVector[1] , -c*nBase[1] )
156 // << " , " << isEqual( thisVector[2] , -c*nBase[2] )
157 // << std::endl ;
158 
159  // check if vector is same direction as nBase and smaller than thick
160  // skip vectors if Basevector close to zero
161  if ( ( differsLess( nBase[0] , thisVector[0] ) || isEqual( thisVector[0] , -c*nBase[0] ) )
162  && ( differsLess( nBase[1] , thisVector[1] ) || isEqual( thisVector[1] , -c*nBase[1] ) )
163  && ( differsLess( nBase[2] , thisVector[2] ) || isEqual( thisVector[2] , -c*nBase[2] ) )
164  && ( ( c <= thick ) || isEqual( c , thick ) ) ) {
165 
166 
167  if( sensorID != 0 ){
168  // return an encoding of the actual sensor that contains the point
169 
170  sensorID->layer = possibleLayer ;
171  sensorID->side = ( p.z() >= 0. ? +1 : -1 ) ;
172 
173  //fg: it seems that this code counts clockwise (negative rotation!)...
174  sensorID->module = ( i ? nLadders - i : 0 ) ;
175 
176  sensorID->sensor = 0 ;
177  }
178 
179  return true ;
180  }
181 
182  } // all ladders
183 
184  } // all layers
185 
186  return false ;
187 
188  } // function isPointInLadder
189 
190 
191 
192 
193  Vector3D ZPlanarParametersImpl::distanceToNearestLadder( Vector3D p , bool sensitive ) const
194  {
195 
196  Vector3D origin( 0,0,0 ) ;
197 
198  // check if point is in already
199  if( isPointInLadder( p , sensitive ) ) {
200 
201  return origin ;
202  }
203 
204  // check if point is origin
205 
206  if( origin.x() == p.x() && origin.y() == p.y() && origin.z() == p.z() ) {
207 
208  if( _layer.getNLayers() > 0 ){
209 
210  // simply use first ladder - distance d , phi0 , gap/2 :
211  //fg: fixme: review definition of phi...
212  return Vector3D( _layer.getLadderDistance(0) , -_layer.getInternalPhi0(0) , _shellGap/2 , Vector3D::cylindrical ) ;
213  }
214 
215  return origin ;
216  }
217 
218  //fg: FIXME: this code needs to be cleaned up and checked
219  // it can also probably be made a bit faster ...
220 
221 
222  double pPhi = getPhiPoint( p ) ;
223 
224  // this Point is solution
225  Vector3D vPointPlane ;
226  double minDistance = 99999.99 ;
227 
228  // go through all layers
229  for ( int nearestLayer=0 ; nearestLayer < _layer.getNLayers() ; nearestLayer++ ) {
230 
231  // get deltaPhi between ladders
232  double deltaPhi = ( 2 * M_PI ) / _layer.getNLadders( nearestLayer ) ;
233 
234  // get extensions for layer
235  double dist, zStart, zEnd, left, right, thick ;
236  if( !sensitive ) {
237  zEnd = _layer.getLadderLength( nearestLayer ) + _shellGap/2 ;
238  left = -_layer.getLadderWidth( nearestLayer ) / 2 - _layer.getLadderOffset( nearestLayer ) ;
239  right = _layer.getLadderWidth( nearestLayer ) / 2 - _layer.getLadderOffset( nearestLayer ) ;
240  thick = _layer.getLadderThickness( nearestLayer ) ;
241  dist = _layer.getLadderDistance( nearestLayer ) ;
242  }
243  if (sensitive ) {
244  zEnd = _layer.getSensitiveLength( nearestLayer ) + _shellGap/2 ;
245  left = -_layer.getSensitiveWidth( nearestLayer ) / 2 - _layer.getSensitiveOffset( nearestLayer ) ;
246  right = _layer.getSensitiveWidth( nearestLayer ) / 2 - _layer.getSensitiveOffset( nearestLayer ) ;
247  thick = _layer.getSensitiveThickness( nearestLayer ) ;
248  dist = _layer.getSensitiveDistance( nearestLayer ) ;
249  }
250 
251  zStart = _shellGap/2 ;
252 
253  // go through all ladders/sensitive
254  for( int i = 0 ; i < _layer.getNLadders( nearestLayer ) ; i++ ) {
255 
256  double phi = _layer.getInternalPhi0( nearestLayer ) + i*deltaPhi ;
257  phi = correctPhiRange( phi ) ;
258 
259  // check start and end phi for this layer
260  double sPhi = correctPhiRange( phi - deltaPhi ) ;
261  double ePhi = correctPhiRange( phi + deltaPhi ) ;
262 
263  // check if point is in range of ladder
264  bool isInRange =
265  ( ( sPhi <= pPhi ) && ( pPhi <= ePhi ) )
266  ||( ( sPhi > ePhi ) && ( pPhi <= ePhi ) )
267  ||( ( sPhi > ePhi ) && ( pPhi >= sPhi ) ) ;
268 
269  if ( !isInRange ) {
270 
271  continue ;
272  }
273 
274  // take normal vector of planes, Base, Side, Z and spacePoint
275 //fg: ------- test vector3d ---------------------------
276 // Vector3D nBase, nSide, nZ ;
277 
278 // nBase[0] = sin( phi ) ;
279 // nBase[1] = cos( phi ) ;
280 // nBase[2] = 0;
281 
282 // nSide[0] = cos( phi ) ;
283 // nSide[1] = -sin( phi ) ;
284 // nSide[2] = 0 ;
285 
286 // nZ[0] = 0 ;
287 // nZ[1] = 0 ;
288 // nZ[2] = 1 ;
289 
290 // if( p[2] < 0 ) {
291 // nZ[2] = - nZ[2] ;
292 // }
293 
294 // Vector3D spacePoint ;
295 // spacePoint[0] = dist * nBase[0] ;
296 // spacePoint[1] = dist * nBase[1] ;
297 // spacePoint[2] = dist * nBase[2] ;
298 
299 
300  Vector3D nBase( sin( phi ), cos( phi ), 0.0 ) ;
301  Vector3D nSide( cos( phi ), -sin( phi ), 0.0 ) ;
302 
303  Vector3D nZ( 0. , 0. , 1. ) ;
304  if( p[2] < 0 ) {
305  nZ[2] = -nZ[2] ;
306  }
307  Vector3D spacePoint = dist * nBase ;
308 
309 //fg: ------- end test vector3d ---------------------------
310 
311 
312 
313  // upper and lower base
314  // inner (1) and outer (2) boundary
315  for ( int j = 1 ; j<=2 ; j ++ ) {
316 
317  double d = 0 ;
318  if( j == 2 ) {
319  d = thick ;
320  }
321 
322 //fg: ------- test vector3d ---------------------------
323 // Vector3D r ;
324 // r[0] = spacePoint[0] + d * nBase[0] ;
325 // r[1] = spacePoint[1] + d * nBase[1] ;
326 // r[2] = spacePoint[2] + d * nBase[2] ;
327  Vector3D r = spacePoint + d * nBase ;
328 //fg: ------- end test vector3d ---------------------------
329 
330 
331  Vector3D myVec = distanceToPlane( p, r, nBase, nSide, nZ, left, right, zStart, zEnd ) ;
332 
333  // debug
334 // Vector3D iP ( p[0]+myVec[0] , p[1]+myVec[1] , p[2]+myVec[2] ) ;
335 // bool isCorrect = isPointInLadder( iP ) ;
336 
337 // if ( !isCorrect ) {
338 // std::cout << "distanceToPlane returns sth wrong. j =" <<j << std::endl ;
339 // }
340 
341  double thisDistance = sqrt( myVec[0] * myVec[0] + myVec[1] * myVec[1] + myVec[2] * myVec[2] ) ;
342  if ( thisDistance <= minDistance ) {
343  minDistance = thisDistance ;
344  vPointPlane = myVec ;
345 
346 // // debug
347 // std::cout << "\nmyVec set " <<"base " << myVec[0] << " , " << myVec[1] << " , " << myVec[2] << std::endl ; //debug
348 // std::cout << " lay " << nearestLayer << " d " << d << " phi " << phi * 180/M_PI << " j " << j <<std::endl ;
349 // std::cout << " r " << r[0] << " , " << r[1] << " , " << r[2]
350 // << " nSide " << nSide[0] << " , " << nSide[1] << " , " << nSide[2]
351 // << " <--> " <
352 // << left << " - " << right << std::endl;
353  }
354 
355  } // for j -- lower & upper boundary
356 
357 
358  // left and right side
359  // 1 left, 2 right
360 
361  for ( int j = 1 ; j<=2 ; j ++ ) {
362 
363  double d = 0;
364  if ( (!sensitive) and ( j == 1 ) ) {
365  d = - _layer.getLadderWidth( nearestLayer ) / 2 - _layer.getLadderOffset( nearestLayer ) ;
366  }
367  if( (!sensitive) and ( j == 2 ) ) {
368  d = _layer.getLadderWidth( nearestLayer ) / 2 - _layer.getLadderOffset( nearestLayer ) ;
369  }
370  if( (sensitive) and ( j == 1 ) ) {
371  d = - _layer.getSensitiveWidth( nearestLayer ) / 2 - _layer.getSensitiveOffset( nearestLayer ) ;
372  }
373  if( (sensitive) and ( j == 2 ) ) {
374  d = _layer.getSensitiveWidth( nearestLayer ) / 2 - _layer.getSensitiveOffset( nearestLayer ) ;
375  }
376 
377  // lower left corner as r
378  Vector3D r ;
379  r[0] = spacePoint[0] + d * nSide[0] ;
380  r[1] = spacePoint[1] + d * nSide[1] ;
381  r[2] = spacePoint[2] + d * nSide[2] ;
382 
383  Vector3D myVec = distanceToPlane( p, r, nSide, nBase, nZ, 0, thick, zStart, zEnd ) ;
384 
385  double thisDistance = sqrt( myVec[0] * myVec[0] + myVec[1] * myVec[1] + myVec[2] * myVec[2] ) ;
386  if ( thisDistance < minDistance ) {
387  minDistance = thisDistance ;
388  vPointPlane = myVec ;
389 // std::cout << "myVec set " <<"left/right " << myVec[0] << " , " << myVec[1] << " , " << myVec[2] << std::endl ; //debug
390  }
391 
392  } // for j -- left and right corner
393 
394 
395 
396  // z Front and Back
397  for( int j = 1 ; j <= 2 ; j++ ) {
398 
399  double d = 0 ;
400 
401  if( j == 1 ) {
402  d = zEnd ;
403  }
404  else {
405  d = zStart ;
406  }
407 
408  // lower middle corner as r
409  Vector3D r;
410  r[0] = spacePoint[0] + d * nZ[0] ;
411  r[1] = spacePoint[1] + d * nZ[1] ;
412  r[2] = spacePoint[2] + d * nZ[2] ;
413 
414  Vector3D myVec = distanceToPlane( p, r, nZ, nBase, nSide, 0, thick, left, right ) ;
415 
416  double thisDistance = sqrt( myVec[0] * myVec[0] + myVec[1] * myVec[1] + myVec[2] * myVec[2] ) ;
417 
418  if ( thisDistance < minDistance ) {
419  minDistance = thisDistance ;
420  vPointPlane = myVec ;
421 // std::cout << "myVec set " <<" front/back " << myVec[0] << " , " << myVec[1] << " , " << myVec[2] << std::endl ; //debug
422 // std::cout << " thick " << thick << " left " << left << " right " << right << std::endl ;
423 // std::cout << " r " << r[0] << " , " << r[1] << " , " << r[2] << std::endl ;
424  }
425 
426  } // for j -- front and back
427 
428  } // for i -- all ladders in layer
429 
430  } // for nearestLayer -- all layers in detetor
431 
432  return vPointPlane ;
433 
434  } // function distanceToLadder
435 
436 
437 
438 
439 
440  Vector3D ZPlanarParametersImpl::intersectionLadder( Vector3D p , Vector3D v , bool sensitive ) const {
441 
442  // return values
443  double currentLength = 99999 ;
444  Vector3D currentPoint (0. , 0. , 0. ) ;
445 
446  for ( int takeLayer=0 ; takeLayer < _layer.getNLayers() ; takeLayer++ ) {
447 
448  // get deltaPhi between ladders
449  double deltaPhi = ( 2 * M_PI ) / _layer.getNLadders( takeLayer ) ;
450 
451  // get extensions for layer
452  double zStart, zEnd, left, right, thick ;
453  if( !sensitive ) {
454  zEnd = _layer.getLadderLength( takeLayer ) + _shellGap/2 ;
455  left = -_layer.getLadderWidth( takeLayer ) / 2 - _layer.getLadderOffset( takeLayer ) ;
456  right = _layer.getLadderWidth( takeLayer ) / 2 - _layer.getLadderOffset( takeLayer ) ;
457  thick = _layer.getLadderThickness( takeLayer ) ;
458  }
459  if (sensitive ) {
460  zEnd = _layer.getSensitiveLength( takeLayer ) + _shellGap/2 ;
461  left = -_layer.getSensitiveWidth( takeLayer ) / 2 - _layer.getSensitiveOffset( takeLayer ) ;
462  right = _layer.getSensitiveWidth( takeLayer ) / 2 - _layer.getSensitiveOffset( takeLayer ) ;
463  thick = _layer.getSensitiveThickness( takeLayer ) ;
464  }
465 
466  // go through both sides
467  for( int side = 0 ; side <2 ; side++ ) {
468  if( side == 0 ) {
469  zStart = _shellGap /2 ;
470  }
471  else {
472  zEnd = - zEnd ;
473  zStart = - _shellGap/2 ;
474  }
475 
476  // go through all ladders/sensitive
477  for( int i = 0 ; i < _layer.getNLadders( takeLayer ) ; i++ ) {
478 
479  double phi = _layer.getInternalPhi0( takeLayer ) + i*deltaPhi ;
480  phi = correctPhiRange( phi ) ;
481 
482  // take normal vector of planes, Base, Side, Z and spacePoint
483  Vector3D nBase, nSide, nZ, spacePoint;
484 
485  nBase[0] = sin( phi ) ;
486  nBase[1] = cos( phi ) ;
487  nBase[2] = 0;
488 
489  nSide[0] = cos( phi ) ;
490  nSide[1] = -sin( phi ) ;
491  nSide[2] = 0 ;
492 
493  nZ[0] = 0 ;
494  nZ[1] = 0 ;
495  nZ[2] = 1 ;
496 
497  // upper and lower base
498  // inner (1) and outer (2) boundary
499  for ( int j = 1 ; j<=2 ; j ++ ) {
500 
501  double d = 0 ;
502  if ( (!sensitive) and ( j == 1 ) ) {
503  d = _layer.getLadderDistance( takeLayer ) ;
504  }
505  if( (!sensitive) and ( j == 2 ) ) {
506  d = _layer.getLadderDistance( takeLayer ) + thick ;
507  }
508  if( (sensitive) and ( j == 1 ) ) {
509  d = _layer.getSensitiveDistance( takeLayer ) ;
510  }
511  if( (sensitive) and ( j == 2 ) ) {
512  d = _layer.getSensitiveDistance( takeLayer ) + thick ;
513  }
514  Vector3D r( d * nBase[0] , d * nBase[1] , d * nBase[2] ) ;
515 
516  // this vector is equivalent to spacepoint -- needed later on
517  if ( j == 1 ) spacePoint = r ;
518 
519  // get intersection point plane/line
520  Vector3D interSection = planeLineIntersection( r, nBase, p, v ) ;
521 
522  // get vector from point to intersection point
523  Vector3D vip ;
524  vip[0] = interSection[0] - p[0] ;
525  vip[1] = interSection[1] - p[1] ;
526  vip[2] = interSection[2] - p[2] ;
527 
528  // get vector from spacepoint to intersection point
529  Vector3D vn ;
530  vn[0] = interSection[0] - spacePoint[0] ;
531  vn[1] = interSection[1] - spacePoint[1] ;
532  vn[2] = interSection[2] - spacePoint[2] ;
533 
534  // correct in Plane
535  Vector3D inBorder = correctToBorderPoint( vn, nSide, nZ, left, right, zStart, zEnd ) ;
536 
537  // check if point has not changed -> intersects with detector
538  if ( isEqual( inBorder , vn ) ) {
539 
540  // check if this is closest intersection
541  double thisLength = sqrt( vip[0]*vip[0] + vip[1]*vip[1] +vip[2]*vip[2] ) ;
542  if( thisLength < currentLength ) {
543  currentPoint = interSection ;
544  currentLength = thisLength ;
545  }
546  }
547 
548  } // for j -- lower & upper boundary
549 
550 
551  // left and right side
552  // 1 left, 2 right
553 
554  for ( int j = 1 ; j<=2 ; j ++ ) {
555 
556  double d = 0;
557  if ( (!sensitive) and ( j == 1 ) ) {
558  d = -( _layer.getLadderWidth( takeLayer ) / 2 + _layer.getLadderOffset( takeLayer ) );
559  }
560  if( (!sensitive) and ( j == 2 ) ) {
561  d = _layer.getLadderWidth( takeLayer ) / 2 - _layer.getLadderOffset( takeLayer ) ;
562  }
563  if( (sensitive) and ( j == 1 ) ) {
564  d = -( _layer.getSensitiveWidth( takeLayer ) / 2 + _layer.getSensitiveOffset( takeLayer ) ) ;
565  }
566  if( (sensitive) and ( j == 2 ) ) {
567  d = _layer.getSensitiveWidth( takeLayer ) / 2 - _layer.getSensitiveOffset( takeLayer ) ;
568  }
569 
570  // lower left corner as r
571  Vector3D r ;
572  r[0] = spacePoint[0] + d * nSide[0] ;
573  r[1] = spacePoint[1] + d * nSide[1] ;
574  r[2] = spacePoint[2] + d * nSide[2] ;
575 
576  // get intersection point plane/line
577  Vector3D interSection = planeLineIntersection( r, nSide, p, v ) ;
578 
579  // get vector from point to intersection point
580  Vector3D vip ;
581  vip[0] = interSection[0] - p[0] ;
582  vip[1] = interSection[1] - p[1] ;
583  vip[2] = interSection[2] - p[2] ;
584 
585  // get vector from spacepoint to intersection point
586  Vector3D vn ;
587  vn[0] = interSection[0] - spacePoint[0] ;
588  vn[1] = interSection[1] - spacePoint[1] ;
589  vn[2] = interSection[2] - spacePoint[2] ;
590 
591  // correct in Plane
592  Vector3D inBorder = correctToBorderPoint( vn, nBase, nZ, 0, thick, zStart, zEnd ) ;
593 
594  if ( isEqual( inBorder , vn )) {
595  // check if this is closest intersection
596  double thisLength = sqrt( vip[0]*vip[0] + vip[1]*vip[1] + vip[2]*vip[2] ) ;
597  if( thisLength < currentLength ) {
598  currentPoint = interSection ;
599  currentLength = thisLength ;
600  }
601 
602  }
603 
604 
605  } // for j -- left and right corner
606 
607 
608  // front and back side in z
609 
610  for( int j = 1 ; j <= 2 ; j++ ) {
611 
612  double d = 0 ;
613 
614  if( j == 1 ) {
615  d = zEnd ;
616  }
617  else {
618  d = zStart ;
619  }
620 
621  // lower left corner as r
622  Vector3D r;
623  r[0] = spacePoint[0] + d * nZ[0] ;
624  r[1] = spacePoint[1] + d * nZ[1] ;
625  r[2] = spacePoint[2] + d * nZ[2] ;
626 
627  // get intersection point plane/line
628  Vector3D interSection = planeLineIntersection( r, nZ, p, v ) ;
629 
630  // get vector from point to intersection point
631  Vector3D vip ;
632  vip[0] = interSection[0] - p[0] ;
633  vip[1] = interSection[1] - p[1] ;
634  vip[2] = interSection[2] - p[2] ;
635 
636  // get vector from spacepoint to intersection point
637  Vector3D vn ;
638  vn[0] = interSection[0] - spacePoint[0] ;
639  vn[1] = interSection[1] - spacePoint[1] ;
640  vn[2] = interSection[2] - spacePoint[2] ;
641 
642  // correct in Plane
643  Vector3D inBorder = correctToBorderPoint( vn, nBase, nSide, 0, thick, left, right ) ;
644 
645  if ( isEqual( inBorder , vn ) ) {
646  // check if this is closest intersection
647  double thisLength = sqrt( vip[0]*vip[0] + vip[1]*vip[1] + vip[2]*vip[2] ) ;
648  if( thisLength < currentLength ) {
649  currentPoint = interSection ;
650  currentLength = thisLength ;
651  }
652  }
653 
654  } // for j -- front and back
655 
656 
657  } // for i -- all ladders in layer
658 
659  } // for all zSides
660 
661  } // for nearestLayer -- all layers in detector
662 
663  return currentPoint ;
664 
665  } // function isIntersectionLadder
666 
667 
668 
669  Vector3D ZPlanarParametersImpl::distanceToPlane( Vector3D p, Vector3D r, Vector3D n, Vector3D u, Vector3D v, float minU, float maxU, float minV, float maxV ) const {
670 
671  // check zeros
672 
673  if ( n[0]==0 && n[1]==0 && n[2]==0 ) {
674 
675  std::cout << "\ndistanceToPlane - Fatal error!" << std::endl ;
676  std::cout << "normal vector (0,0,0)" << std::endl ;
677  return Vector3D (0,0,0) ;
678  }
679 
680  // check u and v
681  if ( !isEqual( n[0]*( u[0]+v[0] ) + n[1]*( u[1]+v[1] ) + n[2]*( u[2]+v[2] ) , 0 ) ){
682  std::cout << "\ndistanceToPlan - Fatal error! "
683  << "n not normal to u and v " << std::endl ;
684  return Vector3D (0,0,0) ;
685  }
686 
687  // The normal n to the plane is n(n0,n1,n2).
688  // The nearest point in the plane is pn(pn0, pn1, pn2)
689  //
690  // The line defined by (pn0-p0)/n0 = (pn1-p1)/n1 = (pn2-p2)/n2 = t
691  // goes through p(p0,p1,p2) and is parallel to n.
692  //
693  // Solving for the point pn we get
694  //
695  // pn0 = n0*t+p0
696  // pn1 = n1*t+p1
697  // pn2 = n2*t+p2
698  //
699  // Now place these values in the equation for the plane:
700  //
701  // n0 * (n0*t+p0) + n1*(n1*t+p1) + n2*(n2*t+p2) + r = 0
702  //
703  // and solve for t:
704  //
705  // t = (-n0*p0-n1*p1-n2*p2-d) / ( n0 * n0 + n1 * n1 + n2 * n2 )
706 
707 
708  // the distance to plane is Vector to spacepoint scalar with n
709 
710  double d = ( r[0] * n[0] + r[1] * n[1] + r[2] * n[2] ) ;
711 
712  double t = - ( n[0] * p[0] + n[1] * p[1] + n[2] * p[2] - d ) ;
713 
714  // nearest point in plane
715  Vector3D pn ;
716  pn[0] = p[0] + n[0] * t;
717  pn[1] = p[1] + n[1] * t;
718  pn[2] = p[2] + n[2] * t;
719 
720  // vector from spacespoint to nearest point in plane
721  Vector3D vn ;
722  vn[0] = pn[0] - r[0] ;
723  vn[1] = pn[1] - r[1] ;
724  vn[2] = pn[2] - r[2] ;
725 
726  // correct to borderPoint
727  Vector3D vCor = correctToBorderPoint( vn , u, v, minU, maxU, minV, maxV ) ;
728 
729  // final vector from point to plane
730  Vector3D final ;
731  final[0] = r[0] + vCor[0] - p[0] ;
732  final[1] = r[1] + vCor[1] - p[1] ;
733  final[2] = r[2] + vCor[2] - p[2] ;
734 
735  return final ;
736 
737  } // function distanceToPlane
738 
739 
740 
741  Vector3D ZPlanarParametersImpl::planeLineIntersection( Vector3D r, Vector3D n, Vector3D linePoint, Vector3D lineDir) const {
742 
743  // calculate distance of plane
744  double distance = ( r[0]*n[0] + r[1]*n[1] + r[2]*n[2] ) / ( sqrt( n[0]*n[0] + n[1]*n[1] + n[2]*n[2] ) ) ;
745 
746  // equations for a point defined by the line with distance to plane = 0
747  // n * ( linePoint + t * lineDir ) - distance = 0
748  // solve for t
749 
750  double numerator = distance - n[0]*linePoint[0] - n[1]*linePoint[1] - n[2]*linePoint[2] ;
751  double denominator = n[0]*lineDir[0] + n[1]*lineDir[1] + n[2]*lineDir[2] ;
752 
753  if ( denominator == 0 ) {
754  // line paralell to plane -> no intersection
755  return Vector3D ( 0. , 0. , 0. ) ;
756  }
757 
758  double t = numerator/denominator ;
759 
760  if( t < 0 ) {
761  // intersection not in direction of line
762  return Vector3D ( 0. , 0. , 0. ) ;
763  }
764 
765  Vector3D final ;
766  final[0] = linePoint[0] + t * lineDir[0] ;
767  final[1] = linePoint[1] + t * lineDir[1] ;
768  final[2] = linePoint[2] + t * lineDir[2] ;
769 
770  return final ;
771 
772  }// function planLineIntersectin
773 
774 
775 
776  Vector3D ZPlanarParametersImpl::correctToBorderPoint( Vector3D vPlane , Vector3D u, Vector3D v, float minU, float maxU, float minV, float maxV ) const {
777 
778  // check if the given vector is within the boundaries
779  double neededUs = vPlane[0] * u[0] + vPlane[1] * u[1] + vPlane[2] * u[2] ;
780  double neededVs = vPlane[0] * v[0] + vPlane[1] * v[1] + vPlane[2] * v[2] ;
781 
782 
783  double goUs = neededUs, goVs = neededVs ;
784  bool changed = false ;
785 
786  if( neededUs < minU ) {
787  goUs = minU ;
788  changed = true ;
789  }
790  if( neededUs > maxU ) {
791  goUs = maxU ;
792  changed = true ;
793  }
794  if( neededVs < minV ) {
795  goVs = minV ;
796  changed = true ;
797  }
798  if( neededVs > maxV ) {
799  goVs = maxV ;
800  changed = true ;
801  }
802 
803  if( !changed ) return vPlane ;
804 
805  Vector3D final ;
806  final[0] = goUs * u[0] + goVs * v[0] ;
807  final[1] = goUs * u[1] + goVs * v[1] ;
808  final[2] = goUs * u[2] + goVs * v[2] ;
809 
810  return final ;
811 
812  } // correctToBorderPoint
813 
814 
815  double ZPlanarParametersImpl::correctPhiRange( double Phi ) const {
816 
817  if( Phi > M_PI ) {
818  return ( Phi - 2 * M_PI ) ;
819  }
820  if( Phi <= -M_PI ) {
821  return ( Phi + 2 * M_PI ) ;
822  }
823 
824  return Phi ;
825 
826  } // function correctPhiRange
827 
828 
829  double ZPlanarParametersImpl::getPhiPoint( Vector3D p ) const {
830 
831  //FG: this is the the angle with the negative y-axis
832  // see comment at the top !
833 
834  // get phi of point in projection 2D
835  double pPhi = 0. ;
836  if( ( p[0] >= 0 ) && ( p[1] == 0 ) )
837  pPhi = -M_PI/2 ;
838 
839  if( ( p[0] < 0 ) && ( p[1] == 0 ) )
840  pPhi = M_PI/2 ;
841 
842  if( ( p[0] == 0 ) && ( p[1] < 0 ) )
843  pPhi = 0 ;
844 
845  if( ( p[0] != 0 ) && ( p[1] < 0 ) )
846  pPhi = atan( p[0] / p[1] ) + M_PI ;
847 
848  else
849  pPhi = atan( p[0] / p[1] ) ;
850 
851  pPhi = correctPhiRange( pPhi ) ;
852 
853  return pPhi ;
854 
855  } // function getPhiPoint
856 
857 
858 
859  bool ZPlanarParametersImpl::isEqual( double valueOne, double valueTwo ) const
860  {
861 
862  // save calculating time if equal
863  if ( valueOne == valueTwo ) return true ;
864 
865 
866  // if one value is 0.0 we use the absolute distance
867  if( valueOne * valueTwo == 0.0 )
868 
869  return fabs( valueOne - valueTwo ) < _EPSILON ;
870 
871 
872  // if both values are smaller than epsilon they are equal for our purpose
873  if ( fabs(valueOne) < _EPSILON && fabs(valueTwo) < _EPSILON )
874 
875  return true ;
876 
877 
878  // otherwise we require the relative distance of the two values to be smaller than epsilon
879  return fabs( 2 * (valueOne - valueTwo) ) / (fabs(valueOne) + fabs(valueTwo) ) < _EPSILON ;
880 
881  }
882 
883 
884  bool ZPlanarParametersImpl::isEqual( Vector3D p1 , Vector3D p2 ) const
885  {
886  for( int i=0 ; i<3 ; i++ ) {
887  if( !( p1[i] == p2[i] ) ) {
888  return false ;
889  }
890  }
891  return true ;
892  } // function isEqual
893 
894 
895  bool ZPlanarParametersImpl::differsLess( double valueOne, double valueTwo ) const
896  {
897 
898  return ( ( fabs( valueOne ) + fabs( valueTwo ) ) < _EPSILON ) ;
899 
900  } // function differsLess
901 
902 
903 } // namespace
904 
905 
906 
gear::ZPlanarLayerLayoutImpl::getNLadders
virtual int getNLadders(int layerIndex) const
The number of ladders in the layer layerIndex - layer indexing starts at 0 for the layer closest to I...
Definition: ZPlanarLayerLayoutImpl.h:51
gear::ZPlanarLayerLayoutImpl::getSensitiveThickness
virtual double getSensitiveThickness(int layerIndex) const
The thickness in mm of the sensitive area in ladders in layer layerIndex.
Definition: ZPlanarLayerLayoutImpl.h:107
gear::ZPlanarLayerLayoutImpl::getSensitiveOffset
virtual double getSensitiveOffset(int layerIndex) const
Same as getLadderOffset() except for the sensitive part of the ladder.
Definition: ZPlanarLayerLayoutImpl.h:112
gear::ZPlanarLayerLayoutImpl::getMaxRadius
virtual double getMaxRadius(int layerIndex, bool sensitive=false) const
returns maximum radius for a given layer
Definition: ZPlanarLayerLayout.cc:49
gear::ZPlanarLayerLayoutImpl::getLadderDistance
virtual double getLadderDistance(int layerIndex) const
The distance of ladders in layer layerIndex from the IP - layer indexing starts at 0 for the layer cl...
Definition: ZPlanarLayerLayoutImpl.h:67
gear::ZPlanarLayerLayoutImpl::getSensitiveWidth
virtual double getSensitiveWidth(int layerIndex) const
The width of the sensitive area in ladders in layer layerIndex in mm.
Definition: ZPlanarLayerLayoutImpl.h:116
gear::ZPlanarLayerLayoutImpl::getSensitiveDistance
virtual double getSensitiveDistance(int layerIndex) const
The distance of sensitive area in ladders in layer layerIndex from the IP.
Definition: ZPlanarLayerLayoutImpl.h:103
gear::Vector3D
Simple three dimensional vector providing the components for cartesian, cylindrical and spherical coo...
Definition: Vector3D.h:18
gear::ZPlanarLayerLayoutImpl::getLadderWidth
virtual double getLadderWidth(int layerIndex) const
The width of the ladder in layer in mm for ladders in layer layerIndex - layer indexing starting at 0...
Definition: ZPlanarLayerLayoutImpl.h:86
gear::Vector3D::z
double z()
Cartesian cartesian z coordinate.
Definition: Vector3D.h:62
gear::ZPlanarLayerLayoutImpl::getNLayers
virtual int getNLayers() const
The total number of layers.
Definition: ZPlanarLayerLayoutImpl.h:46
gear::ZPlanarLayerLayoutImpl::getLadderThickness
virtual double getLadderThickness(int layerIndex) const
The thickness in mm of the ladders in layerIndex - layer indexing starting at 0 for the layer closest...
Definition: ZPlanarLayerLayoutImpl.h:72
gear::ZPlanarLayerLayoutImpl::getLadderOffset
virtual double getLadderOffset(int layerIndex) const
The offset of the ladder in mm defines the shift of the ladder in the direction of increasing phi per...
Definition: ZPlanarLayerLayoutImpl.h:81
gear::ZPlanarLayerLayoutImpl::getSensitiveLength
virtual double getSensitiveLength(int layerIndex) const
The length of the sensitive area in ladders in z direction in mm for ladders in layer layerIndex...
Definition: ZPlanarLayerLayoutImpl.h:120
gear::ZPlanarParametersImpl::distanceToNearestLadder
virtual Vector3D distanceToNearestLadder(Vector3D p) const
returns vector from given point p to nearest ladder
Definition: ZPlanarParametersImpl.h:125
gear::ZPlanarParametersImpl::intersectionLadder
virtual Vector3D intersectionLadder(Vector3D p, Vector3D v) const
returns the first point where a given strainght line (parameters point p and direction v) crosses a l...
Definition: ZPlanarParametersImpl.h:138
gear::SensorID
Helper struct for decoding a sensor ID.
Definition: GEAR.h:133
gear::ZPlanarLayerLayoutImpl::getEndInnerPhi
virtual double getEndInnerPhi(int layerIndex, bool sensitive=false) const
returns ending phi for first ladder in layer layerIndex (on side facing IP)
Definition: ZPlanarLayerLayout.cc:81
gear::ZPlanarLayerLayoutImpl::getLadderLength
virtual double getLadderLength(int layerIndex) const
The (half) length of the ladder in z direction in mm for ladders in layer layerIndex - layer indexing...
Definition: ZPlanarLayerLayoutImpl.h:94
gear::ZPlanarParametersImpl::isPointInLadder
virtual bool isPointInLadder(Vector3D p) const
returns whether a point is inside a ladder
Definition: ZPlanarParametersImpl.h:113
gear::ZPlanarParametersImpl::ZPlanarParametersImpl
ZPlanarParametersImpl(int type, double shellInnerRadius, double shellOuterRadius, double shellHalfLength, double shellGap, double shellRadLength)
C&#39;tor.
Definition: ZPlanarParametersImpl.cc:18
gear::ZPlanarLayerLayoutImpl::getStartInnerPhi
virtual double getStartInnerPhi(int layerIndex, bool sensitive=false) const
returns starting phi for first ladder in layer layerIndex (on side facing IP)
Definition: ZPlanarLayerLayout.cc:68
Generated by   doxygen 1.8.5