CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
EcalEndcapGeometry.cc
Go to the documentation of this file.
6 #include <CLHEP/Geometry/Point3D.h>
7 #include <CLHEP/Geometry/Plane3D.h>
8 
12 typedef HepGeom::Plane3D<CCGFloat> Pl3D ;
13 
14 EcalEndcapGeometry::EcalEndcapGeometry( void )
15  : _nnmods( 316 ),
16  _nncrys( 25 ),
17  zeP( 0. ),
18  zeN( 0. ),
19  m_wref( 0. ),
20  m_del( 0. ),
21  m_nref( 0 ),
22  m_borderMgr( 0 ),
23  m_borderPtrVec( 0 ),
24  m_avgZ( -1 ),
25  m_cellVec( k_NumberOfCellsForCorners )
26 {
27  m_xlo[0] = 999.;
28  m_xlo[1] = 999.;
29  m_xhi[0] = -999.;
30  m_xhi[1] = -999.;
31  m_ylo[0] = 999.;
32  m_ylo[1] = 999.;
33  m_yhi[0] = -999.;
34  m_yhi[1] = -999.;
35  m_xoff[0] = 0.;
36  m_xoff[1] = 0.;
37  m_yoff[0] = 0.;
38  m_yoff[1] = 0.;
39 }
40 
41 EcalEndcapGeometry::~EcalEndcapGeometry()
42 {
43  delete m_borderPtrVec ;
44  delete m_borderMgr ;
45 }
46 
47 unsigned int
48 EcalEndcapGeometry::alignmentTransformIndexLocal( const DetId& id )
49 {
50  const CaloGenericDetId gid ( id ) ;
51 
52  assert( gid.isEE() ) ;
53  unsigned int index ( EEDetId(id).ix()/51 + ( EEDetId(id).zside()<0 ? 0 : 2 ) ) ;
54 
55  return index ;
56 }
57 
58 DetId
59 EcalEndcapGeometry::detIdFromLocalAlignmentIndex( unsigned int iLoc )
60 {
61  return EEDetId( 20 + 50*( iLoc%2 ), 50, 2*( iLoc/2 ) - 1 ) ;
62 }
63 
64 unsigned int
65 EcalEndcapGeometry::alignmentTransformIndexGlobal( const DetId& /*id*/ )
66 {
67  return (unsigned int)DetId::Ecal - 1 ;
68 }
69 
70 void
71 EcalEndcapGeometry::initializeParms()
72 {
73  zeP=0.;
74  zeN=0.;
75  unsigned nP=0;
76  unsigned nN=0;
77  m_nref = 0 ;
78 
79  for( uint32_t i ( 0 ) ; i != m_cellVec.size() ; ++i )
80  {
81  const CaloCellGeometry* cell ( cellGeomPtr(i) ) ;
82  if( 0 != cell )
83  {
84  const CCGFloat z ( cell->getPosition().z() ) ;
85  if(z>0.)
86  {
87  zeP+=z;
88  ++nP;
89  }
90  else
91  {
92  zeN+=z;
93  ++nN;
94  }
95  const EEDetId myId ( EEDetId::detIdFromDenseIndex(i) ) ;
96  const unsigned int ix ( myId.ix() ) ;
97  const unsigned int iy ( myId.iy() ) ;
98  if( ix > m_nref ) m_nref = ix ;
99  if( iy > m_nref ) m_nref = iy ;
100  }
101  }
102  if( 0 < nP ) zeP/=(CCGFloat)nP;
103  if( 0 < nN ) zeN/=(CCGFloat)nN;
104 
105  m_xlo[0] = 999 ;
106  m_xhi[0] = -999 ;
107  m_ylo[0] = 999 ;
108  m_yhi[0] = -999 ;
109  m_xlo[1] = 999 ;
110  m_xhi[1] = -999 ;
111  m_ylo[1] = 999 ;
112  m_yhi[1] = -999 ;
113  for( uint32_t i ( 0 ) ; i != m_cellVec.size() ; ++i )
114  {
115  const CaloCellGeometry* cell ( cellGeomPtr(i) ) ;
116  if( 0 != cell )
117  {
118  const GlobalPoint p ( cell->getPosition() ) ;
119  const CCGFloat z ( p.z() ) ;
120  const CCGFloat zz ( 0 > z ? zeN : zeP ) ;
121  const CCGFloat x ( p.x()*zz/z ) ;
122  const CCGFloat y ( p.y()*zz/z ) ;
123 
124  if( 0 > z && x < m_xlo[0] ) m_xlo[0] = x ;
125  if( 0 < z && x < m_xlo[1] ) m_xlo[1] = x ;
126  if( 0 > z && y < m_ylo[0] ) m_ylo[0] = y ;
127  if( 0 < z && y < m_ylo[1] ) m_ylo[1] = y ;
128 
129  if( 0 > z && x > m_xhi[0] ) m_xhi[0] = x ;
130  if( 0 < z && x > m_xhi[1] ) m_xhi[1] = x ;
131  if( 0 > z && y > m_yhi[0] ) m_yhi[0] = y ;
132  if( 0 < z && y > m_yhi[1] ) m_yhi[1] = y ;
133  }
134  }
135 
136  m_xoff[0] = ( m_xhi[0] + m_xlo[0] )/2. ;
137  m_xoff[1] = ( m_xhi[1] + m_xlo[1] )/2. ;
138  m_yoff[0] = ( m_yhi[0] + m_ylo[0] )/2. ;
139  m_yoff[1] = ( m_yhi[1] + m_ylo[1] )/2. ;
140 
141  m_del = ( m_xhi[0] - m_xlo[0] + m_xhi[1] - m_xlo[1] +
142  m_yhi[0] - m_ylo[0] + m_yhi[1] - m_ylo[1] ) ;
143 
144  if( 1 != m_nref ) m_wref = m_del/(4.*(m_nref-1)) ;
145 
146  m_xlo[0] -= m_wref/2 ;
147  m_xlo[1] -= m_wref/2 ;
148  m_xhi[0] += m_wref/2 ;
149  m_xhi[1] += m_wref/2 ;
150 
151  m_ylo[0] -= m_wref/2 ;
152  m_ylo[1] -= m_wref/2 ;
153  m_yhi[0] += m_wref/2 ;
154  m_yhi[1] += m_wref/2 ;
155 
156  m_del += m_wref ;
157 /*
158  std::cout<<"zeP="<<zeP<<", zeN="<<zeN<<", nP="<<nP<<", nN="<<nN<<std::endl ;
159 
160  std::cout<<"xlo[0]="<<m_xlo[0]<<", xlo[1]="<<m_xlo[1]<<", xhi[0]="<<m_xhi[0]<<", xhi[1]="<<m_xhi[1]
161  <<"\nylo[0]="<<m_ylo[0]<<", ylo[1]="<<m_ylo[1]<<", yhi[0]="<<m_yhi[0]<<", yhi[1]="<<m_yhi[1]<<std::endl ;
162 
163  std::cout<<"xoff[0]="<<m_xoff[0]<<", xoff[1]"<<m_xoff[1]<<", yoff[0]="<<m_yoff[0]<<", yoff[1]"<<m_yoff[1]<<std::endl ;
164 
165  std::cout<<"nref="<<m_nref<<", m_wref="<<m_wref<<std::endl ;
166 */
167 }
168 
169 
170 unsigned int
171 EcalEndcapGeometry::xindex( CCGFloat x,
172  CCGFloat z ) const
173 {
174  const CCGFloat xlo ( 0 > z ? m_xlo[0] : m_xlo[1] ) ;
175  const int i ( 1 + int( ( x - xlo )/m_wref ) ) ;
176 
177  return ( 1 > i ? 1 :
178  ( m_nref < (unsigned int) i ? m_nref : (unsigned int) i ) ) ;
179 
180 }
181 
182 unsigned int
183 EcalEndcapGeometry::yindex( CCGFloat y,
184  CCGFloat z ) const
185 {
186  const CCGFloat ylo ( 0 > z ? m_ylo[0] : m_ylo[1] ) ;
187  const int i ( 1 + int( ( y - ylo )/m_wref ) ) ;
188 
189  return ( 1 > i ? 1 :
190  ( m_nref < (unsigned int) i ? m_nref : (unsigned int) i ) ) ;
191 }
192 
193 EEDetId
194 EcalEndcapGeometry::gId( float x,
195  float y,
196  float z ) const
197 {
198  const CCGFloat fac ( fabs( ( 0 > z ? zeN : zeP )/z ) ) ;
199  const unsigned int ix ( xindex( x*fac, z ) ) ;
200  const unsigned int iy ( yindex( y*fac, z ) ) ;
201  const unsigned int iz ( z>0 ? 1 : -1 ) ;
202 
203  if( EEDetId::validDetId( ix, iy, iz ) )
204  {
205  return EEDetId( ix, iy, iz ) ; // first try is on target
206  }
207  else // try nearby coordinates, spiraling out from center
208  {
209  for( unsigned int i ( 1 ) ; i != 6 ; ++i )
210  {
211  for( unsigned int k ( 0 ) ; k != 8 ; ++k )
212  {
213  const int jx ( 0 == k || 4 == k || 5 == k ? +i :
214  ( 1 == k || 5 < k ? -i : 0 ) ) ;
215  const int jy ( 2 == k || 4 == k || 6 == k ? +i :
216  ( 3 == k || 5 == k || 7 == k ? -i : 0 ) ) ;
217  if( EEDetId::validDetId( ix + jx, iy + jy, iz ) )
218  {
219  return EEDetId( ix + jx, iy + jy, iz ) ;
220  }
221  }
222  }
223  }
224  return EEDetId() ; // nowhere near any crystal
225 }
226 
227 
228 // Get closest cell, etc...
229 DetId
230 EcalEndcapGeometry::getClosestCell( const GlobalPoint& r ) const
231 {
232  try
233  {
234  EEDetId mycellID ( gId( r.x(), r.y(), r.z() ) ) ; // educated guess
235 
236  if( EEDetId::validDetId( mycellID.ix(),
237  mycellID.iy(),
238  mycellID.zside() ) )
239  {
240  // now get points in convenient ordering
241 
242  Pt3D A;
243  Pt3D B;
244  Pt3D C;
245  Pt3D point(r.x(),r.y(),r.z());
246  // D.K. : equation of plane : AA*x+BB*y+CC*z+DD=0;
247  // finding equation for each edge
248 
249  // ================================================================
250  CCGFloat x,y,z;
251  unsigned offset=0;
252  int zsign=1;
253  //================================================================
254 
255  // compute the distance of the point with respect of the 4 crystal lateral planes
256 
257 
258  if( 0 != getGeometry(mycellID) )
259  {
260  const GlobalPoint& myPosition=getGeometry(mycellID)->getPosition();
261 
262  x=myPosition.x();
263  y=myPosition.y();
264  z=myPosition.z();
265 
266  offset=0;
267  // This will disappear when Andre has applied his fix
268  zsign=1;
269 
270  if(z>0)
271  {
272  if(x>0&&y>0)
273  offset=1;
274  else if(x<0&&y>0)
275  offset=2;
276  else if(x>0&&y<0)
277  offset=0;
278  else if (x<0&&y<0)
279  offset=3;
280  zsign=1;
281  }
282  else
283  {
284  if(x>0&&y>0)
285  offset=3;
286  else if(x<0&&y>0)
287  offset=2;
288  else if(x>0&&y<0)
289  offset=0;
290  else if(x<0&&y<0)
291  offset=1;
292  zsign=-1;
293  }
294  GlobalPoint corners[8];
295  for(unsigned ic=0;ic<4;++ic)
296  {
297  corners[ic]=getGeometry(mycellID)->getCorners()[(unsigned)((zsign*ic+offset)%4)];
298  corners[4+ic]=getGeometry(mycellID)->getCorners()[(unsigned)(4+(zsign*ic+offset)%4)];
299  }
300 
301  CCGFloat SS[4];
302  for (short i=0; i < 4 ; ++i)
303  {
304  A = Pt3D(corners[i%4].x(),corners[i%4].y(),corners[i%4].z());
305  B = Pt3D(corners[(i+1)%4].x(),corners[(i+1)%4].y(),corners[(i+1)%4].z());
306  C = Pt3D(corners[4+(i+1)%4].x(),corners[4+(i+1)%4].y(),corners[4+(i+1)%4].z());
307  Pl3D plane(A,B,C);
308  plane.normalize();
309  CCGFloat distance = plane.distance(point);
310  if (corners[0].z()<0.) distance=-distance;
311  SS[i] = distance;
312  }
313 
314  // Only one move in necessary direction
315 
316  const bool yout ( 0 > SS[0]*SS[2] ) ;
317  const bool xout ( 0 > SS[1]*SS[3] ) ;
318 
319  if( yout || xout )
320  {
321  const int ydel ( !yout ? 0 : ( 0 < SS[0] ? -1 : 1 ) ) ;
322  const int xdel ( !xout ? 0 : ( 0 < SS[1] ? -1 : 1 ) ) ;
323  const unsigned int ix ( mycellID.ix() + xdel ) ;
324  const unsigned int iy ( mycellID.iy() + ydel ) ;
325  const unsigned int iz ( mycellID.zside() ) ;
326  if( EEDetId::validDetId( ix, iy, iz ) )
327  mycellID = EEDetId( ix, iy, iz ) ;
328  }
329 
330  return mycellID;
331  }
332  return DetId(0);
333  }
334  }
335  catch ( cms::Exception &e )
336  {
337  return DetId(0);
338  }
339  return DetId(0);
340 }
341 
343 EcalEndcapGeometry::getCells( const GlobalPoint& r,
344  double dR ) const
345 {
346  CaloSubdetectorGeometry::DetIdSet dis ; // return object
347  if( 0.000001 < dR )
348  {
349  if( dR > M_PI/2. ) // this code assumes small dR
350  {
351  dis = CaloSubdetectorGeometry::getCells( r, dR ) ; // base class version
352  }
353  else
354  {
355  const float dR2 ( dR*dR ) ;
356  const float reta ( r.eta() ) ;
357  const float rphi ( r.phi() ) ;
358  const float rx ( r.x() ) ;
359  const float ry ( r.y() ) ;
360  const float rz ( r.z() ) ;
361  const float fac ( std::abs( zeP/rz ) ) ;
362  const float xx ( rx*fac ) ; // xyz at endcap z
363  const float yy ( ry*fac ) ;
364  const float zz ( rz*fac ) ;
365 
366  const float xang ( std::atan( xx/zz ) ) ;
367  const float lowX ( zz>0 ? zz*std::tan( xang - dR ) : zz*std::tan( xang + dR ) ) ;
368  const float highX ( zz>0 ? zz*std::tan( xang + dR ) : zz*std::tan( xang - dR ) ) ;
369  const float yang ( std::atan( yy/zz ) ) ;
370  const float lowY ( zz>0 ? zz*std::tan( yang - dR ) : zz*std::tan( yang + dR ) ) ;
371  const float highY ( zz>0 ? zz*std::tan( yang + dR ) : zz*std::tan( yang - dR ) ) ;
372 
373  const float refxlo ( 0 > rz ? m_xlo[0] : m_xlo[1] ) ;
374  const float refxhi ( 0 > rz ? m_xhi[0] : m_xhi[1] ) ;
375  const float refylo ( 0 > rz ? m_ylo[0] : m_ylo[1] ) ;
376  const float refyhi ( 0 > rz ? m_yhi[0] : m_yhi[1] ) ;
377 
378  if( lowX < refxhi && // proceed if any possible overlap with the endcap
379  lowY < refyhi &&
380  highX > refxlo &&
381  highY > refylo )
382  {
383  const int ix_ctr ( xindex( xx, rz ) ) ;
384  const int iy_ctr ( yindex( yy, rz ) ) ;
385  const int iz ( rz>0 ? 1 : -1 ) ;
386 
387  const int ix_hi ( ix_ctr + int( ( highX - xx )/m_wref ) + 2 ) ;
388  const int ix_lo ( ix_ctr - int( ( xx - lowX )/m_wref ) - 2 ) ;
389 
390  const int iy_hi ( iy_ctr + int( ( highY - yy )/m_wref ) + 2 ) ;
391  const int iy_lo ( iy_ctr - int( ( yy - lowY )/m_wref ) - 2 ) ;
392 
393  for( int kx ( ix_lo ) ; kx <= ix_hi ; ++kx )
394  {
395  if( kx > 0 &&
396  kx <= (int) m_nref )
397  {
398  for( int ky ( iy_lo ) ; ky <= iy_hi ; ++ky )
399  {
400  if( ky > 0 &&
401  ky <= (int) m_nref )
402  {
403  if( EEDetId::validDetId( kx, ky, iz ) ) // reject invalid ids
404  {
405  const EEDetId id ( kx, ky, iz ) ;
406  const CaloCellGeometry* cell = &m_cellVec[ id.denseIndex()];
407  const float eta (cell->etaPos() ) ;
408  const float phi (cell->phiPos() ) ;
409  if( reco::deltaR2( eta, phi, reta, rphi ) < dR2 ) dis.insert( id ) ;
410  }
411  }
412  }
413  }
414  }
415  }
416  }
417  }
418  return dis;
419 }
420 
421 const EcalEndcapGeometry::OrderedListOfEBDetId*
422 EcalEndcapGeometry::getClosestBarrelCells( EEDetId id ) const
423 {
424  OrderedListOfEBDetId* ptr ( 0 ) ;
425  if( 0 != id.rawId() &&
426  0 != getGeometry( id ) )
427  {
428  const float phi ( 370. +
429  getGeometry( id )->getPosition().phi().degrees() );
430  const int iPhi ( 1 + int(phi)%360 ) ;
431  const int iz ( id.zside() ) ;
432  if( 0 == m_borderMgr )
433  {
434  m_borderMgr = new EZMgrFL<EBDetId>( 720*9, 9 ) ;
435  }
436  if( 0 == m_borderPtrVec )
437  {
438  m_borderPtrVec = new VecOrdListEBDetIdPtr() ;
439  m_borderPtrVec->reserve( 720 ) ;
440  for( unsigned int i ( 0 ) ; i != 720 ; ++i )
441  {
442  const int kz ( 360>i ? -1 : 1 ) ;
443  const int iEta ( kz*85 ) ;
444  const int iEtam1 ( kz*84 ) ;
445  const int iEtam2 ( kz*83 ) ;
446  const int jPhi ( i%360 + 1 ) ;
447  OrderedListOfEBDetId& olist ( *new OrderedListOfEBDetId( m_borderMgr ) );
448  olist[0]=EBDetId( iEta , jPhi ) ;
449  olist[1]=EBDetId( iEta , myPhi( jPhi+1 ) ) ;
450  olist[2]=EBDetId( iEta , myPhi( jPhi-1 ) ) ;
451  olist[3]=EBDetId( iEtam1, jPhi ) ;
452  olist[4]=EBDetId( iEtam1, myPhi( jPhi+1 ) ) ;
453  olist[5]=EBDetId( iEtam1, myPhi( jPhi-1 ) ) ;
454  olist[6]=EBDetId( iEta , myPhi( jPhi+2 ) ) ;
455  olist[7]=EBDetId( iEta , myPhi( jPhi-2 ) ) ;
456  olist[8]=EBDetId( iEtam2, jPhi ) ;
457  m_borderPtrVec->push_back( &olist ) ;
458  }
459  }
460  ptr = (*m_borderPtrVec)[ ( iPhi - 1 ) + ( 0>iz ? 0 : 360 ) ] ;
461  }
462  return ptr ;
463 }
464 
465 void
466 EcalEndcapGeometry::localCorners( Pt3DVec& lc ,
467  const CCGFloat* pv ,
468  unsigned int /*i*/,
469  Pt3D& ref )
470 {
471  TruncatedPyramid::localCorners( lc, pv, ref ) ;
472 }
473 
474 void
475 EcalEndcapGeometry::newCell( const GlobalPoint& f1 ,
476  const GlobalPoint& f2 ,
477  const GlobalPoint& f3 ,
478  const CCGFloat* parm ,
479  const DetId& detId )
480 {
481  const unsigned int cellIndex ( EEDetId( detId ).denseIndex() ) ;
482  m_cellVec[ cellIndex ] =
483  TruncatedPyramid( cornersMgr(), f1, f2, f3, parm ) ;
484  m_validIds.push_back( detId ) ;
485 }
486 
487 
488 CCGFloat
489 EcalEndcapGeometry::avgAbsZFrontFaceCenter() const
490 {
491  if( 0 > m_avgZ )
492  {
493  CCGFloat sum ( 0 ) ;
494  for( unsigned int i ( 0 ) ; i != m_cellVec.size() ; ++i )
495  {
496  const CaloCellGeometry* cell ( cellGeomPtr(i) ) ;
497  if( 0 != cell )
498  {
499  sum += fabs( cell->getPosition().z() ) ;
500  }
501  }
502  m_avgZ = sum/m_cellVec.size() ;
503  }
504  return m_avgZ ;
505 }
506 
507 const CaloCellGeometry*
508 EcalEndcapGeometry::cellGeomPtr( uint32_t index ) const
509 {
510  const CaloCellGeometry* cell ( &m_cellVec[ index ] ) ;
511  return ( m_cellVec.size() < index ||
512  0 == cell->param() ? 0 : cell ) ;
513 }
static EEDetId detIdFromDenseIndex(uint32_t din)
Definition: EEDetId.h:221
double degrees(double radiants)
def degrees
int i
Definition: DBlmapReader.cc:9
double_binary B
Definition: DDStreamer.cc:234
float phiPos() const
Geom::Phi< T > phi() const
Definition: PV3DBase.h:69
T y() const
Definition: PV3DBase.h:63
#define abs(x)
Definition: mlp_lapack.h:159
std::vector< Pt3D > Pt3DVec
T eta() const
float float float z
HepGeom::Plane3D< CCGFloat > Pl3D
virtual DetIdSet getCells(const GlobalPoint &r, double dR) const
Get a list of all cells within a dR of the given cell.
T z() const
Definition: PV3DBase.h:64
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
CaloCellGeometry::CCGFloat CCGFloat
unsigned int offset(bool)
int k[5][pyjets_maxn]
CaloCellGeometry::Pt3D Pt3D
Definition: DetId.h:20
#define M_PI
Definition: BFit3D.cc:3
static bool validDetId(int crystal_ix, int crystal_iy, int iz)
Definition: EEDetId.h:249
CaloCellGeometry::Pt3DVec Pt3DVec
HepGeom::Point3D< CCGFloat > Pt3D
T1 deltaR2(T1 eta1, T2 phi1, T3 eta2, T4 phi2)
Definition: deltaR.h:58
A base class to handle the particular shape of Ecal Xtals. Taken from ORCA Calorimetry Code...
T eta() const
Definition: PV3DBase.h:76
float etaPos() const
Definition: DDAxes.h:10
T x() const
Definition: PV3DBase.h:62
*vegas h *****************************************************used in the default bin number in original ***version of VEGAS is ***a higher bin number might help to derive a more precise ***grade subtle point
Definition: invegas.h:5
Definition: DDAxes.h:10