CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
HcalGeometry.cc
Go to the documentation of this file.
7 #include <algorithm>
8 
10  theTopology ( new HcalTopology ),
11  m_ownsTopology ( true )
12 {
13 }
14 
16  theTopology ( topology ) ,
17  m_ownsTopology ( false )
18 {
19 }
20 
21 
23 {
24  if( m_ownsTopology ) delete theTopology ;
25 }
26 
27 
28 void
30 {
31  const std::vector<DetId>& baseIds ( CaloSubdetectorGeometry::getValidDetIds() ) ;
32  for( unsigned int i ( 0 ) ; i != baseIds.size() ; ++i )
33  {
34  const DetId id ( baseIds[i] );
35  if( id.subdetId() == HcalBarrel )
36  {
37  m_hbIds.push_back( id ) ;
38  }
39  else
40  {
41  if( id.subdetId() == HcalEndcap )
42  {
43  m_heIds.push_back( id ) ;
44  }
45  else
46  {
47  if( id.subdetId() == HcalOuter )
48  {
49  m_hoIds.push_back( id ) ;
50  }
51  else
52  {
53  if( id.subdetId() == HcalForward )
54  {
55  m_hfIds.push_back( id ) ;
56  }
57  }
58  }
59  }
60  }
61  std::sort( m_hbIds.begin(), m_hbIds.end() ) ;
62  std::sort( m_heIds.begin(), m_heIds.end() ) ;
63  std::sort( m_hoIds.begin(), m_hoIds.end() ) ;
64  std::sort( m_hfIds.begin(), m_hfIds.end() ) ;
65  m_emptyIds.resize( 0 ) ;
66 }
67 
68 const std::vector<DetId>&
70  int subdet ) const
71 {
72  if( 0 != subdet &&
73  0 == m_hbIds.size() ) fillDetIds() ;
74  return ( 0 == subdet ? CaloSubdetectorGeometry::getValidDetIds() :
75  ( HcalBarrel == subdet ? m_hbIds :
76  ( HcalEndcap == subdet ? m_heIds :
77  ( HcalOuter == subdet ? m_hoIds :
78  ( HcalForward == subdet ? m_hfIds : m_emptyIds ) ) ) ) ) ;
79 }
80 
82 
83  // Now find the closest eta_bin, eta value of a bin i is average
84  // of eta[i] and eta[i-1]
85  double abseta = fabs(r.eta());
86 
87  // figure out subdetector, giving preference to HE in HE/HF overlap region
89  if (abseta <= theHBHEEtaBounds[theTopology->lastHBRing()] ) {
90  bc = HcalBarrel;
91  } else if (abseta <= theHBHEEtaBounds[theTopology->lastHERing()] ) {
92  bc = HcalEndcap;
93  } else {
94  bc = HcalForward;
95  }
96 
97  if (bc == HcalForward) {
98  static const double z_short=1137.0;
99  int etaring = etaRing(bc, abseta); // This is safer
100  /*
101  static const double z_long=1115.0;
102  // determine front-face eta
103  double radius=sqrt(pow(r.x(),2)+pow(r.y(),2));
104  double trueAeta=asinh(z_long/radius);
105  // find eta bin
106  int etaring = etaRing(bc, trueAeta);
107  */
108  if (etaring>theTopology->lastHFRing()) etaring=theTopology->lastHFRing();
109 
110  int phibin = phiBin(r.phi(), etaring);
111 
112  // add a sign to the etaring
113  int etabin = (r.z() > 0) ? etaring : -etaring;
114  // Next line is premature depth 1 and 2 can coexist for large z-extent
115  HcalDetId bestId(bc,etabin,phibin,((fabs(r.z())>=z_short)?(2):(1)));
116  return bestId;
117  } else {
118 
119  // find eta bin
120  int etaring = etaRing(bc, abseta);
121 
122  int phibin = phiBin(r.phi(), etaring);
123 
124  // add a sign to the etaring
125  int etabin = (r.z() > 0) ? etaring : -etaring;
126 
127  //Now do depth if required
128  int dbin = 1;
129  double pointrz=0, drz=99999.;
130  HcalDetId currentId(bc, etabin, phibin, dbin);
131  if (bc == HcalBarrel) pointrz = r.mag();
132  else pointrz = std::abs(r.z());
133  HcalDetId bestId;
134  for ( ; currentId != HcalDetId(); theTopology->incrementDepth(currentId)) {
135  const CaloCellGeometry * cell = getGeometry(currentId);
136  assert(cell != 0);
137  double rz;
138  if (bc == HcalEndcap) rz = std::abs(cell->getPosition().z());
139  else rz = cell->getPosition().mag();
140  if (std::abs(pointrz-rz)<drz) {
141  bestId = currentId;
142  drz = std::abs(pointrz-rz);
143  }
144  }
145 
146  return bestId;
147  }
148 }
149 
150 
151 int HcalGeometry::etaRing(HcalSubdetector bc, double abseta) const
152 {
153  int etaring;
154  if( bc == HcalForward ) {
155  for(etaring = theTopology->firstHFRing();
156  etaring <= theTopology->lastHFRing(); ++etaring)
157  {
158  if(theHFEtaBounds[etaring-theTopology->firstHFRing()+1] > abseta) break;
159  }
160  }
161  else
162  {
163  for(etaring = 1;
164  etaring <= theTopology->lastHERing(); ++etaring)
165  {
166  if(theHBHEEtaBounds[etaring] >= abseta) break;
167  }
168  }
169 
170  return etaring;
171 }
172 
173 
174 int HcalGeometry::phiBin(double phi, int etaring) const
175 {
176  static const double twopi = M_PI+M_PI;
177  //put phi in correct range (0->2pi)
178  if(phi<0.0) phi += twopi;
179  if(phi>twopi) phi -= twopi;
180  int nphibins = theTopology->nPhiBins(etaring);
181  int phibin= static_cast<int>(phi/twopi*nphibins)+1;
182  int iphi;
183 
184  // rings 40 and 41 are offset wrt the other phi numbering
185  // 1 1 1 2
186  // ------------------------------
187  // 72 36 36 1
188  if(etaring >= theTopology->firstHFQuadPhiRing())
189  {
190  phi+=(twopi/36); //shift by half tower.
191  phibin=static_cast<int>(phi/twopi*nphibins);
192  if (phibin==0) phibin=18;
193  iphi=phibin*4-1; // 71,3,5,
194  } else {
195  // convert to the convention of numbering 1,3,5, in 36 phi bins
196  iphi=(phibin-1)*(72/nphibins) + 1;
197  }
198 
199  return iphi;
200 }
201 
204  double dR ) const
205 {
206  CaloSubdetectorGeometry::DetIdSet dis; // this is the return object
207 
208  if( 0.000001 < dR )
209  {
210  if( dR > M_PI/2. ) // this version needs "small" dR
211  {
212  dis = CaloSubdetectorGeometry::getCells( r, dR ) ; // base class version
213  }
214  else
215  {
216  const double dR2 ( dR*dR ) ;
217  const double reta ( r.eta() ) ;
218  const double rphi ( r.phi() ) ;
219  const double lowEta ( reta - dR ) ;
220  const double highEta ( reta + dR ) ;
221  const double lowPhi ( rphi - dR ) ;
222  const double highPhi ( rphi + dR ) ;
223 
224  const double hfEtaHi ( theHFEtaBounds[ theTopology->lastHFRing() -
225  theTopology->firstHFRing() + 1 ] ) ;
226 
227  if( highEta > -hfEtaHi &&
228  lowEta < hfEtaHi ) // in hcal
229  {
231 
232  for( unsigned int is ( 0 ) ; is != 4 ; ++is )
233  {
234  const int sign ( reta>0 ? 1 : -1 ) ;
235  const int ieta_center ( sign*etaRing( hs[is], fabs( reta ) ) ) ;
236  const int ieta_lo ( ( 0 < lowEta*sign ? sign : -sign )*etaRing( hs[is], fabs( lowEta ) ) ) ;
237  const int ieta_hi ( ( 0 < highEta*sign ? sign : -sign )*etaRing( hs[is], fabs( highEta ) ) ) ;
238  const int iphi_lo ( phiBin( lowPhi , ieta_center ) ) ;
239  const int iphi_hi ( phiBin( highPhi, ieta_center ) ) ;
240  const int jphi_lo ( iphi_lo>iphi_hi ? iphi_lo - 72 : iphi_lo ) ;
241  const int jphi_hi ( iphi_hi ) ;
242 
243  const int idep_lo ( 1 == is ? 4 : 1 ) ;
244  const int idep_hi ( 1 == is ? 4 :
245  ( 2 == is ? 3 : 2 ) ) ;
246  for( int ieta ( ieta_lo ) ; ieta <= ieta_hi ; ++ieta ) // over eta limits
247  {
248  if( ieta != 0 )
249  {
250  for( int jphi ( jphi_lo ) ; jphi <= jphi_hi ; ++jphi ) // over phi limits
251  {
252  const int iphi ( 1 > jphi ? jphi+72 : jphi ) ;
253 
254  for( int idep ( idep_lo ) ; idep <= idep_hi ; ++idep )
255  {
256  if( HcalDetId::validDetId( hs[is], ieta, iphi, idep ) )
257  {
258  const HcalDetId did ( hs[is], ieta, iphi, idep ) ;
259  const CaloCellGeometry* cell ( getGeometry( did ) );
260  if( 0 != cell )
261  {
262  const GlobalPoint& p ( cell->getPosition() ) ;
263  const double eta ( p.eta() ) ;
264  const double phi ( p.phi() ) ;
265  if( reco::deltaR2( eta, phi, reta, rphi ) < dR2 ) dis.insert( did ) ;
266  }
267  }
268  }
269  }
270  }
271  }
272  }
273  }
274  }
275  }
276  return dis;
277 }
278 
279 
280 unsigned int
282 {
283  const CaloGenericDetId gid ( id ) ;
284 
285  assert( gid.isHcal() ) ;
286 
287 
288  const HcalDetId hid ( id ) ;
289 
290  const int jz ( ( hid.zside() + 1 )/2 ) ;
291 
292  const int zoff ( jz*numberOfAlignments()/2 ) ;
293 
294  const int detoff ( zoff +
295  ( gid.isHB() ? 0 :
296  ( gid.isHE() ? numberOfBarrelAlignments()/2 :
297  ( gid.isHF() ? ( numberOfBarrelAlignments() +
301  numberOfForwardAlignments() )/2 ) ) ) ) ;
302 
303  const int iphi ( hid.iphi() ) ;
304 
305  unsigned int index ( numberOfAlignments() ) ;
306  if( gid.isHO() )
307  {
308  const int ieta ( hid.ieta() ) ;
309  const int ring ( ieta < -10 ? 0 :
310  ( ieta < -4 ? 1 :
311  ( ieta < 5 ? 2 :
312  ( ieta < 11 ? 3 : 4 ) ) ) ) ;
313 
314  index = detoff + 12*ring + ( iphi - 1 )%6 ;
315  }
316  else
317  {
318  index = detoff + ( iphi - 1 )%4 ;
319  }
320 
321  assert( index < numberOfAlignments() ) ;
322  return index ;
323 }
324 
325 unsigned int
327 {
328  return (unsigned int)DetId::Hcal - 1 ;
329 }
330 
331 std::vector<HepGeom::Point3D<double> >
332 HcalGeometry::localCorners( const double* pv,
333  unsigned int i,
334  HepGeom::Point3D<double> & ref )
335 {
336  const HcalDetId hid ( HcalDetId::detIdFromDenseIndex( i ) ) ;
337  const CaloGenericDetId cgid ( hid ) ;
338  if( cgid.isHF() )
339  {
340  return ( calogeom::IdealZPrism::localCorners( pv, ref ) ) ;
341  }
342  else
343  {
344  return ( calogeom::IdealObliquePrism::localCorners( pv, ref ) ) ;
345  }
346 }
347 
350  const GlobalPoint& f2 ,
351  const GlobalPoint& f3 ,
353  const double* parm ,
354  const DetId& detId )
355 {
356  const CaloGenericDetId cgid ( detId ) ;
357 
358  assert( cgid.isHcal() ) ;
359 
360  if( cgid.isHF() )
361  {
362  return ( new calogeom::IdealZPrism( f1, mgr, parm ) ) ;
363  }
364  else
365  {
366  return ( new calogeom::IdealObliquePrism( f1, mgr, parm ) ) ;
367  }
368 }
static unsigned int numberOfBarrelAlignments()
Definition: HcalGeometry.h:50
int firstHFRing() const
Definition: HcalTopology.h:65
std::vector< DetId > m_heIds
Definition: HcalGeometry.h:92
int i
Definition: DBlmapReader.cc:9
virtual DetId getClosestCell(const GlobalPoint &r) const
Definition: HcalGeometry.cc:81
static unsigned int alignmentTransformIndexLocal(const DetId &id)
bool isHE() const
std::vector< DetId > m_emptyIds
Definition: HcalGeometry.h:95
Geom::Phi< T > phi() const
Definition: PV3DBase.h:63
static std::vector< HepGeom::Point3D< double > > localCorners(const double *pv, HepGeom::Point3D< double > &ref)
Definition: IdealZPrism.cc:34
int nPhiBins(int etaRing) const
how many phi segments in this ring
#define abs(x)
Definition: mlp_lapack.h:159
int lastHBRing() const
Definition: HcalTopology.h:62
const HcalTopology * theTopology
Definition: HcalGeometry.h:89
static bool validDetId(HcalSubdetector subdet, int tower_ieta, int tower_iphi, int depth)
Definition: HcalDetId.cc:70
T eta() const
virtual const std::vector< DetId > & getValidDetIds(DetId::Detector det=DetId::Detector(0), int subdet=0) const
Get a list of valid detector ids (for the given subdetector)
Definition: HcalGeometry.cc:69
virtual CaloSubdetectorGeometry::DetIdSet getCells(const GlobalPoint &r, double dR) const
Get a list of all cells within a dR of the given cell.
virtual const std::vector< DetId > & getValidDetIds(DetId::Detector det=DetId::Detector(0), int subdet=0) const
Get a list of valid detector ids (for the given subdetector)
virtual const CaloCellGeometry * getGeometry(const DetId &id) const
Get the cell geometry of a given detector id. Should return false if not found.
std::vector< DetId > m_hbIds
Definition: HcalGeometry.h:91
static std::vector< HepGeom::Point3D< double > > localCorners(const double *pv, HepGeom::Point3D< double > &ref)
std::vector< DetId > m_hoIds
Definition: HcalGeometry.h:93
virtual DetIdSet getCells(const GlobalPoint &r, double dR) const
Get a list of all cells within a dR of the given cell.
T mag() const
Definition: PV3DBase.h:61
static HcalDetId detIdFromDenseIndex(uint32_t di)
Definition: HcalDetId.cc:155
T z() const
Definition: PV3DBase.h:58
int lastHFRing() const
Definition: HcalTopology.h:66
static unsigned int alignmentTransformIndexGlobal(const DetId &id)
bool m_ownsTopology
Definition: HcalGeometry.h:96
HcalSubdetector
Definition: HcalAssistant.h:32
std::vector< DetId > m_hfIds
Definition: HcalGeometry.h:94
bool isHcal() const
static std::vector< HepGeom::Point3D< double > > localCorners(const double *pv, unsigned int i, HepGeom::Point3D< double > &ref)
int etaRing(HcalSubdetector bc, double abseta) const
helper methods for getClosestCell
double deltaR2(double eta1, double phi1, double eta2, double phi2)
Definition: deltaR.h:13
Definition: DetId.h:20
bool incrementDepth(HcalDetId &id) const
static unsigned int numberOfEndcapAlignments()
Definition: HcalGeometry.h:52
#define M_PI
Definition: BFit3D.cc:3
static const double theHFEtaBounds[]
bool isHF() const
virtual ~HcalGeometry()
The HcalGeometry will delete all its cell geometries at destruction time.
Definition: HcalGeometry.cc:22
Detector
Definition: DetId.h:26
Definition: EZMgrFL.h:8
int firstHFQuadPhiRing() const
Definition: HcalTopology.h:71
T eta() const
Definition: PV3DBase.h:70
void fillDetIds() const
Definition: HcalGeometry.cc:29
static unsigned int numberOfForwardAlignments()
Definition: HcalGeometry.h:56
bool isHO() const
bool isHB() const
static CaloCellGeometry * newCell(const GlobalPoint &f1, const GlobalPoint &f2, const GlobalPoint &f3, CaloCellGeometry::CornersMgr *mgr, const double *parm, const DetId &detId)
const GlobalPoint & getPosition() const
static const double theHBHEEtaBounds[]
int phiBin(double phi, int etaring) const
static unsigned int numberOfAlignments()
Definition: HcalGeometry.h:58
int lastHERing() const
Definition: HcalTopology.h:64
Definition: DDAxes.h:10