CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
HGCalGeometry.cc
Go to the documentation of this file.
1 /* for High Granularity Calorimeter
2  * This geometry is essentially driven by topology,
3  * which is thus encapsulated in this class.
4  * This makes this geometry not suitable to be loaded
5  * by regular CaloGeometryLoader<T>
6  */
12 
13 #include <cmath>
14 
15 #include <Math/Transform3D.h>
16 #include <Math/EulerAngles.h>
17 
19 typedef std::vector<float> ParmVec;
20 
21 //#define EDM_ML_DEBUG
22 
24  : m_topology( topology_ ),
25  m_cellVec( topology_.totalGeomModules()),
26  m_validGeomIds( topology_.totalGeomModules()),
27  m_halfType( topology_.detectorType()),
28  m_subdet( topology_.subDetector()) {
29 
30  m_validIds.reserve( topology().totalModules());
31 #ifdef EDM_ML_DEBUG
32  std::cout << "Expected total # of Geometry Modules "
33  << topology().totalGeomModules() << std::endl;
34 #endif
35 }
36 
38 
40 
42 }
43 
45  const CCGFloat* pv,
46  unsigned int i,
47  Pt3D& ref) {
48  FlatTrd::localCorners( lc, pv, ref ) ;
49 }
50 
52  const GlobalPoint& f2 ,
53  const GlobalPoint& f3 ,
54  const CCGFloat* parm ,
55  const DetId& detId ) {
56 
57  DetId geomId;
58  int cells;
60  if (topology().dddConstants().geomMode() == HGCalGeometryMode::Square) {
61  geomId = (detId.subdetId() == HGCEE ?
62  (DetId)(HGCEEDetId(detId).geometryCell()) :
63  (DetId)(HGCHEDetId(detId).geometryCell()));
64  cells = topology().dddConstants().maxCells(id.iLay,true);
65  } else {
66  geomId = (DetId)(HGCalDetId(detId).geometryCell());
67  cells = topology().dddConstants().numberCellsHexagon(id.iSec);
68 #ifdef EDM_ML_DEBUG
69  std::cout << "NewCell " << HGCalDetId(detId) << " GEOM "
70  << HGCalDetId(geomId) << std::endl;
71 #endif
72  }
73  const uint32_t cellIndex (topology().detId2denseGeomId(detId));
74 
75  m_cellVec.at( cellIndex ) = FlatTrd( cornersMgr(), f1, f2, f3, parm ) ;
76  m_validGeomIds.at( cellIndex ) = geomId ;
77 
78 #ifdef EDM_ML_DEBUG
79  std::cout << "Store for DetId " << std::hex << detId.rawId() << " GeomId "
80  << geomId.rawId() << std::dec << " Index " << cellIndex
81  << " cells " << cells << std::endl;
82  unsigned int nOld = m_validIds.size();
83 #endif
84  for (int cell = 0; cell < cells; ++cell) {
85  id.iCell = cell;
86  DetId idc = topology().encode(id);
87  if (topology().valid(idc)) {
88  m_validIds.emplace_back(idc);
89  if ((topology().dddConstants().geomMode() == HGCalGeometryMode::Square) &&
90  (!m_halfType)) {
91  id.iSubSec = -id.iSubSec;
92  m_validIds.emplace_back( topology().encode(id));
93  id.iSubSec = -id.iSubSec;
94  }
95  }
96  }
97 
98 #ifdef EDM_ML_DEBUG
99  std::cout << "HGCalGeometry::newCell-> [" << cellIndex << "]"
100  << " front:" << f1.x() << '/' << f1.y() << '/' << f1.z()
101  << " back:" << f2.x() << '/' << f2.y() << '/' << f2.z()
102  << " eta|phi " << m_cellVec[cellIndex].etaPos() << ":"
103  << m_cellVec[cellIndex].phiPos() << " id:";
104  if (topology().dddConstants().geomMode() != HGCalGeometryMode::Square) {
105  std::cout << HGCalDetId(detId);
106  } else if (m_subdet == HGCEE) {
107  std::cout << HGCEEDetId(detId);
108  } else {
109  std::cout << HGCHEDetId(detId);
110  }
111  unsigned int nNew = m_validIds.size();
112  std::cout << " with valid DetId from " << nOld << " to " << nNew
113  << std::endl;
114  std::cout << "Cell[" << cellIndex << "] " << std::hex << geomId.rawId()
115  << ":" << m_validGeomIds[cellIndex].rawId() << std::dec
116  << std::endl;
117 #endif
118 }
119 
120 std::shared_ptr<const CaloCellGeometry> HGCalGeometry::getGeometry(const DetId& id) const {
121  if (id == DetId()) return nullptr; // nothing to get
122  DetId geoId;
123  if (topology().dddConstants().geomMode() == HGCalGeometryMode::Square) {
124  geoId = (id.subdetId() == HGCEE ?
125  (DetId)(HGCEEDetId(id).geometryCell()) :
126  (DetId)(HGCHEDetId(id).geometryCell()));
127  } else {
128  geoId = (DetId)(HGCalDetId(id).geometryCell());
129  }
130  const uint32_t cellIndex (topology().detId2denseGeomId(geoId));
131  const GlobalPoint pos = (id != geoId) ? getPosition(id) : GlobalPoint();
132  return cellGeomPtr (cellIndex, pos);
133 
134 }
135 
136 bool HGCalGeometry::present(const DetId& id) const {
137  if (id == DetId()) return false;
138  DetId geoId;
139  if (topology().dddConstants().geomMode() == HGCalGeometryMode::Square) {
140  geoId = (id.subdetId() == HGCEE ?
141  (DetId)(HGCEEDetId(id).geometryCell()) :
142  (DetId)(HGCHEDetId(id).geometryCell()));
143  } else {
144  geoId = (DetId)(HGCalDetId(id).geometryCell());
145  }
146  const uint32_t index (topology().detId2denseGeomId(geoId));
147  return (nullptr != getGeometryRawPtr(index)) ;
148 }
149 
151 
152  unsigned int cellIndex = indexFor(id);
153  GlobalPoint glob;
154  if (cellIndex < m_cellVec.size()) {
156  std::pair<float,float> xy;
157  if (topology().dddConstants().geomMode() == HGCalGeometryMode::Square) {
158  xy = topology().dddConstants().locateCell(id_.iCell,id_.iLay,id_.iSubSec,true);
159  } else {
160  xy = topology().dddConstants().locateCellHex(id_.iCell,id_.iSec,true);
161  }
162  const HepGeom::Point3D<float> lcoord(xy.first,xy.second,0);
163  glob = m_cellVec[cellIndex].getPosition(lcoord);
164 #ifdef EDM_ML_DEBUG
165  std::cout << "getPosition:: index " << cellIndex << " Local " << lcoord.x()
166  << ":" << lcoord.y() << " ID " << id_.iCell << ":" << id_.iLay
167  << " Global " << glob << std::endl;
168 #endif
169  }
170  return glob;
171 }
172 
174 
175  HGCalGeometry::CornersVec co (8, GlobalPoint(0,0,0));
176  unsigned int cellIndex = indexFor(id);
177  if (cellIndex < m_cellVec.size()) {
179  std::pair<float,float> xy;
180  if (topology().dddConstants().geomMode() == HGCalGeometryMode::Square) {
181  xy = topology().dddConstants().locateCell(id_.iCell,id_.iLay,id_.iSubSec,true);
182  } else {
183  xy = topology().dddConstants().locateCellHex(id_.iCell,id_.iSec,true);
184  }
185  float dz = m_cellVec[cellIndex].param()[0];
186  float dx = 0.5*m_cellVec[cellIndex].param()[11];
187  static const int signx[] = {-1,-1,1,1,-1,-1,1,1};
188  static const int signy[] = {-1,1,1,-1,-1,1,1,-1};
189  static const int signz[] = {-1,-1,-1,-1,1,1,1,1};
190  for (unsigned int i = 0; i != 8; ++i) {
191  const HepGeom::Point3D<float> lcoord(xy.first+signx[i]*dx,xy.second+signy[i]*dx,signz[i]*dz);
192  co[i] = m_cellVec[cellIndex].getPosition(lcoord);
193  }
194  }
195  return co;
196 }
197 
199  unsigned int cellIndex = getClosestCellIndex(r);
200  if (cellIndex < m_cellVec.size()) {
202  HepGeom::Point3D<float> local;
203  if (topology().dddConstants().geomMode() == HGCalGeometryMode::Square) {
204  local = m_cellVec[cellIndex].getLocal(r);
205  } else {
206  if (r.z() > 0) local = HepGeom::Point3D<float>(r.x(),r.y(),0);
207  else local = HepGeom::Point3D<float>(-r.x(),r.y(),0);
208  }
209  std::pair<int,int> kxy =
210  topology().dddConstants().assignCell(local.x(),local.y(),id_.iLay,
211  id_.iSubSec,true);
212  id_.iCell = kxy.second;
214  id_.iSubSec = kxy.first;
215  } else {
216  id_.iSec = kxy.first;
217  id_.iSubSec = topology().dddConstants().waferTypeT(kxy.first);
218  if (id_.iSubSec != 1) id_.iSubSec = -1;
219  }
220 #ifdef EDM_ML_DEBUG
221  std::cout << "getClosestCell: local " << local << " Id " << id_.zside
222  << ":" << id_.iLay << ":" << id_.iSec << ":" << id_.iSubSec
223  << ":" << id_.iCell << std::endl;
224 #endif
225 
226  //check if returned cell is valid
227  if (id_.iCell>=0) return topology().encode(id_);
228  }
229 
230  //if not valid or out of bounds return a null DetId
231  return DetId();
232 }
233 
236  return dss;
237 }
238 
240  if (m_subdet == HGCEE) return "HGCalEE";
241  else if (m_subdet == HGCHEF) return "HGCalHEFront";
242  else if (m_subdet == HGCHEB) return "HGCalHEBack";
243  else return "Unknown";
244 }
245 
246 unsigned int HGCalGeometry::indexFor(const DetId& id) const {
247  unsigned int cellIndex = m_cellVec.size();
248  if (id != DetId()) {
249  DetId geoId;
250  if (topology().dddConstants().geomMode() == HGCalGeometryMode::Square) {
251  geoId = (id.subdetId() == HGCEE ?
252  (DetId)(HGCEEDetId(id).geometryCell()) :
253  (DetId)(HGCHEDetId(id).geometryCell()));
254  } else {
255  geoId = (DetId)(HGCalDetId(id).geometryCell());
256  }
257  cellIndex = topology().detId2denseGeomId(geoId);
258 #ifdef EDM_ML_DEBUG
259  std::cout << "indexFor " << std::hex << id.rawId() << ":" << geoId.rawId()
260  << std::dec << " index " << cellIndex << std::endl;
261 #endif
262  }
263  return cellIndex;
264 }
265 
266 unsigned int HGCalGeometry::sizeForDenseIndex() const {
267  return topology().totalGeomModules();
268 }
269 
271  // Modify the RawPtr class
272  const CaloCellGeometry* cell(&m_cellVec[index]);
273  return (m_cellVec.size() < index ||
274  nullptr == cell->param() ? nullptr : cell);
275 }
276 
277 std::shared_ptr<const CaloCellGeometry> HGCalGeometry::cellGeomPtr(uint32_t index) const {
278  if ((index >= m_cellVec.size()) || (m_validGeomIds[index].rawId() == 0))
279  return nullptr;
280  static const auto do_not_delete = [](const void*){};
281  auto cell = std::shared_ptr<const CaloCellGeometry>(&m_cellVec[index],do_not_delete);
282  if (nullptr == cell->param()) return nullptr;
283  return cell;
284 }
285 
286 std::shared_ptr<const CaloCellGeometry> HGCalGeometry::cellGeomPtr(uint32_t index, const GlobalPoint& pos) const {
287  if ((index >= m_cellVec.size()) || (m_validGeomIds[index].rawId() == 0))
288  return nullptr;
289  if (pos == GlobalPoint()) return cellGeomPtr(index);
290  auto cell = std::make_shared<FlatTrd>(m_cellVec[index]);
291  cell->setPosition(pos);
292 #ifdef EDM_ML_DEBUG
293 //std::cout << "cellGeomPtr " << newcell << ":" << cell << std::endl;
294 #endif
295  if (nullptr == cell->param()) return nullptr;
296  return cell;
297 }
298 
300  edm::LogError("HGCalGeom") << "HGCalGeometry::addValidID is not implemented";
301 }
302 
303 
304 unsigned int HGCalGeometry::getClosestCellIndex (const GlobalPoint& r) const {
305 
306  float phip = r.phi();
307  float zp = r.z();
308  unsigned int cellIndex = m_cellVec.size();
309  float dzmin(9999), dphimin(9999), dphi10(0.175);
310  for (unsigned int k=0; k<m_cellVec.size(); ++k) {
311  float dphi = phip-m_cellVec[k].phiPos();
312  while (dphi > M_PI) dphi -= 2*M_PI;
313  while (dphi <= -M_PI) dphi += 2*M_PI;
314  if (fabs(dphi) < dphi10) {
315  float dz = fabs(zp - m_cellVec[k].getPosition().z());
316  if (dz < (dzmin+0.001)) {
317  dzmin = dz;
318  if (fabs(dphi) < (dphimin+0.01)) {
319  cellIndex = k;
320  dphimin = fabs(dphi);
321  } else {
322  if (cellIndex >= m_cellVec.size()) cellIndex = k;
323  }
324  }
325  }
326  }
327 #ifdef EDM_ML_DEBUG
328  std::cout << "getClosestCellIndex::Input " << zp << ":" << phip << " Index "
329  << cellIndex;
330  if (cellIndex < m_cellVec.size())
331  std::cout << " Cell z " << m_cellVec[cellIndex].getPosition().z() << ":"
332  << dzmin << " phi " << m_cellVec[cellIndex].phiPos() << ":"
333  << dphimin;
334  std::cout << std::endl;
335 
336 #endif
337  return cellIndex;
338 }
339 
340 // FIXME: Change sorting algorithm if needed
341 namespace {
342  struct rawIdSort {
343  bool operator()( const DetId& a, const DetId& b ) {
344  return( a.rawId() < b.rawId());
345  }
346  };
347 }
348 
350  m_validIds.shrink_to_fit();
351  std::sort( m_validIds.begin(), m_validIds.end(), rawIdSort());
352 }
353 
357  CaloSubdetectorGeometry::IVec& dinsVector ) const {
358 
359  unsigned int numberOfCells = topology().totalGeomModules(); // total Geom Modules both sides
362 
363  trVector.reserve( numberOfCells * numberOfTransformParms());
364  iVector.reserve( numberOfCells );
365  dimVector.reserve( numberOfShapes * numberOfParametersPerShape );
366  dinsVector.reserve( numberOfCells );
367 
368  if (topology().geomMode() == HGCalGeometryMode::Square) {
369  for (unsigned int k=0; k <topology().dddConstants().volumes(); ++k) {
372  params[0] = vol.dz;
373  params[1] = params[2] = 0;
374  params[3] = params[7] = vol.h;
375  params[4] = params[8] = vol.bl;
376  params[5] = params[9] = vol.tl;
377  params[6] = params[10]= vol.alpha;
378  params[11]= vol.cellSize;
379  dimVector.insert( dimVector.end(), params.begin(), params.end());
380  }
381  } else {
382  for (unsigned itr=0; itr<topology().dddConstants().getTrFormN(); ++itr) {
384  int layer = mytr.lay;
385  for (int wafer=0; wafer<topology().dddConstants().sectors(); ++wafer) {
386  if (topology().dddConstants().waferInLayer(wafer,layer,true)) {
387  HGCalParameters::hgtrap vol = topology().dddConstants().getModule(wafer, true, true);
389  params[0] = vol.dz;
390  params[1] = params[2] = 0;
391  params[3] = params[7] = vol.h;
392  params[4] = params[8] = vol.bl;
393  params[5] = params[9] = vol.tl;
394  params[6] = params[10]= vol.alpha;
395  params[11]= vol.cellSize;
396  dimVector.insert( dimVector.end(), params.begin(), params.end());
397  }
398  }
399  }
400  }
401 
402  for (unsigned int i( 0 ); i < numberOfCells; ++i) {
403  DetId detId = m_validGeomIds[i];
404  int layer = ((detId.subdetId() == ForwardSubdetector::HGCEE) ?
405  (HGCEEDetId(detId).layer()) :
406  (HGCHEDetId(detId).layer()));
407  dinsVector.emplace_back( topology().detId2denseGeomId( detId ));
408  iVector.emplace_back( layer );
409 
410  Tr3D tr;
411  auto ptr = cellGeomPtr( i );
412  if ( nullptr != ptr ) {
413  ptr->getTransform( tr, ( Pt3DVec* ) nullptr );
414 
415  if( Tr3D() == tr ) { // there is no rotation
416  const GlobalPoint& gp( ptr->getPosition());
417  tr = HepGeom::Translate3D( gp.x(), gp.y(), gp.z());
418  }
419 
420  const CLHEP::Hep3Vector tt( tr.getTranslation());
421  trVector.emplace_back( tt.x());
422  trVector.emplace_back( tt.y());
423  trVector.emplace_back( tt.z());
424  if (6 == numberOfTransformParms()) {
425  const CLHEP::HepRotation rr( tr.getRotation());
426  const ROOT::Math::Transform3D rtr( rr.xx(), rr.xy(), rr.xz(), tt.x(),
427  rr.yx(), rr.yy(), rr.yz(), tt.y(),
428  rr.zx(), rr.zy(), rr.zz(), tt.z());
430  rtr.GetRotation( ea );
431  trVector.emplace_back( ea.Phi());
432  trVector.emplace_back( ea.Theta());
433  trVector.emplace_back( ea.Psi());
434  }
435  }
436  }
437 }
438 
440 
std::vector< DetId > m_validGeomIds
A base class to handle the particular shape of HGCal volumes.
Definition: FlatTrd.h:19
std::vector< CCGFloat > DimVec
virtual unsigned int numberOfParametersPerShape() const
std::string cellElement() const
std::shared_ptr< const CaloCellGeometry > cellGeomPtr(uint32_t index) const override
unsigned int sizeForDenseIndex() const
HGCalParameters::hgtrform getTrForm(unsigned int k) const
unsigned int getClosestCellIndex(const GlobalPoint &r) const
~HGCalGeometry() override
Geom::Phi< T > phi() const
Definition: PV3DBase.h:69
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
std::vector< unsigned int > IVec
T y() const
Definition: PV3DBase.h:63
void newCell(const GlobalPoint &f1, const GlobalPoint &f2, const GlobalPoint &f3, const CCGFloat *parm, const DetId &detId) override
std::vector< float > ParmVec
std::vector< GlobalPoint > CornersVec
Definition: HGCalGeometry.h:37
std::vector< CCGFloat > TrVec
GlobalPoint getPosition(const DetId &id) const
int layer() const
get the layer #
Definition: HGCHEDetId.h:50
HepGeom::Transform3D Tr3D
std::pair< float, float > locateCell(int cell, int lay, int type, bool reco) const
CaloCellGeometry::CCGFloat CCGFloat
CornersVec getCorners(const DetId &id) const
Returns the corner points of this cell&#39;s volume.
CaloCellGeometry::Tr3D Tr3D
CaloCellGeometry::Pt3DVec Pt3DVec
Definition: HGCalGeometry.h:34
std::pair< float, float > locateCellHex(int cell, int wafer, bool reco) const
uint32_t rawId() const
get the raw id
Definition: DetId.h:43
unsigned int totalGeomModules() const
HGCalGeometry(const HGCalTopology &topology)
HGCalDetId geometryCell() const
Definition: HGCalDetId.h:36
void initializeParms() override
unsigned int getTrFormN() const
const CCGFloat * param() const
virtual uint32_t detId2denseGeomId(const DetId &id) const
int numberCellsHexagon(int wafer) const
virtual unsigned int numberOfShapes() const
HGCHEDetId geometryCell() const
Definition: HGCHEDetId.cc:46
std::vector< float > ParmVec
DetId encode(const DecodedDetId &id_) const
std::vector< DetId > m_validIds
T z() const
Definition: PV3DBase.h:64
const HGCalTopology & topology() const
Definition: HGCalGeometry.h:98
def pv(vc)
Definition: MetAnalyzer.py:6
std::set< DetId > DetIdSet
Definition: HGCalGeometry.h:36
HGCalGeometryMode::GeometryMode geomMode() const
int sectors() const
int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:37
static void localCorners(Pt3DVec &vec, const CCGFloat *pv, Pt3D &ref)
Definition: FlatTrd.cc:162
DetIdSet getCells(const GlobalPoint &r, double dR) const override
Get a list of all cells within a dR of the given cell.
void localCorners(Pt3DVec &lc, const CCGFloat *pv, unsigned int i, Pt3D &ref)
#define M_PI
int k[5][pyjets_maxn]
#define TYPELOOKUP_DATA_REG(_dataclass_)
Definition: typelookup.h:96
DecodedDetId decode(const DetId &id) const
Definition: DetId.h:18
AlgebraicVector EulerAngles
Definition: Definitions.h:36
int layer() const
get the layer #
Definition: HGCEEDetId.h:49
ForwardSubdetector m_subdet
const HGCalDDDConstants & dddConstants() const
HGCEEDetId geometryCell() const
Definition: HGCEEDetId.h:34
DetId getClosestCell(const GlobalPoint &r) const override
CaloCellGeometry::CornersMgr * cornersMgr()
void getSummary(CaloSubdetectorGeometry::TrVec &trVector, CaloSubdetectorGeometry::IVec &iVector, CaloSubdetectorGeometry::DimVec &dimVector, CaloSubdetectorGeometry::IVec &dinsVector) const override
double b
Definition: hdecay.h:120
def encode(args, files)
std::shared_ptr< const CaloCellGeometry > getGeometry(const DetId &id) const override
Get the cell geometry of a given detector id. Should return false if not found.
CaloCellGeometry::Tr3D Tr3D
CaloCellGeometry::Pt3D Pt3D
Definition: HGCalGeometry.h:33
double a
Definition: hdecay.h:121
virtual void fillNamedParams(DDFilteredView fv)
std::pair< int, int > assignCell(float x, float y, int lay, int subSec, bool reco) const
const CaloCellGeometry * getGeometryRawPtr(uint32_t index) const override
void addValidID(const DetId &id)
virtual unsigned int numberOfTransformParms() const
HGCalParameters::hgtrap getModule(unsigned int k, bool hexType, bool reco) const
int maxCells(bool reco) const
T x() const
Definition: PV3DBase.h:62
int waferTypeT(int wafer) const
unsigned int indexFor(const DetId &id) const override
bool present(const DetId &id) const override
is this detid present in the geometry?
CellVec m_cellVec
bool waferInLayer(int wafer, int lay, bool reco) const
unsigned int volumes() const