CMS 3D CMS Logo

HGCalGeomParameters.cc
Go to the documentation of this file.
2 
5 
21 
22 #include <algorithm>
23 #include <sstream>
24 #include <unordered_set>
25 
26 //#define EDM_ML_DEBUG
27 using namespace geant_units::operators;
28 
29 const double tolerance = 0.001;
30 const double tolmin = 1.e-20;
31 
33 #ifdef EDM_ML_DEBUG
34  edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters::HGCalGeomParameters "
35  << "constructor";
36 #endif
37 }
38 
40  HGCalParameters& php,
41  const std::string& sdTag1,
42  const DDCompactView* cpv,
43  const std::string& sdTag2,
44  const std::string& sdTag3,
46  DDFilteredView fv = _fv;
47  bool dodet(true);
48  std::map<int, HGCalGeomParameters::layerParameters> layers;
49  std::vector<HGCalParameters::hgtrform> trforms;
50  std::vector<bool> trformUse;
51 
52  while (dodet) {
53  const DDSolid& sol = fv.logicalPart().solid();
54  // Layers first
55  std::vector<int> copy = fv.copyNumbers();
56  int nsiz = static_cast<int>(copy.size());
57  int lay = (nsiz > 0) ? copy[nsiz - 1] : 0;
58  int zp = (nsiz > 2) ? copy[nsiz - 3] : -1;
59  if (zp != 1)
60  zp = -1;
61  if (lay == 0) {
62  throw cms::Exception("DDException") << "Funny layer # " << lay << " zp " << zp << " in " << nsiz << " components";
63  } else {
64  if (std::find(php.layer_.begin(), php.layer_.end(), lay) == php.layer_.end())
65  php.layer_.emplace_back(lay);
66  auto itr = layers.find(lay);
67  if (itr == layers.end()) {
68  double rin(0), rout(0);
70  if ((sol.shape() == DDSolidShape::ddpolyhedra_rz) || (sol.shape() == DDSolidShape::ddpolyhedra_rrz)) {
71  const DDPolyhedra& polyhedra = static_cast<DDPolyhedra>(sol);
72  const std::vector<double>& rmin = polyhedra.rMinVec();
73  const std::vector<double>& rmax = polyhedra.rMaxVec();
74  rin = 0.5 * HGCalParameters::k_ScaleFromDDD * (rmin[0] + rmin[1]);
75  rout = 0.5 * HGCalParameters::k_ScaleFromDDD * (rmax[0] + rmax[1]);
76  } else if (sol.shape() == DDSolidShape::ddtubs) {
77  const DDTubs& tube = static_cast<DDTubs>(sol);
78  rin = HGCalParameters::k_ScaleFromDDD * tube.rIn();
79  rout = HGCalParameters::k_ScaleFromDDD * tube.rOut();
80  }
81  HGCalGeomParameters::layerParameters laypar(rin, rout, zz);
82  layers[lay] = laypar;
83  }
84  DD3Vector x, y, z;
85  fv.rotation().GetComponents(x, y, z);
86  const CLHEP::HepRep3x3 rotation(x.X(), y.X(), z.X(), x.Y(), y.Y(), z.Y(), x.Z(), y.Z(), z.Z());
87  const CLHEP::HepRotation hr(rotation);
89  if (std::abs(xx) < tolerance)
90  xx = 0;
92  if (std::abs(yy) < tolerance)
93  yy = 0;
95  const CLHEP::Hep3Vector h3v(xx, yy, zz);
97  mytrf.zp = zp;
98  mytrf.lay = lay;
99  mytrf.sec = 0;
100  mytrf.subsec = 0;
101  mytrf.h3v = h3v;
102  mytrf.hr = hr;
103  trforms.emplace_back(mytrf);
104  trformUse.emplace_back(false);
105  }
106  dodet = fv.next();
107  }
108 
109  // Then wafers
110  // This assumes layers are build starting from 1 (which on 25 Jan 2016, they
111  // were) to ensure that new copy numbers are always added to the end of the
112  // list.
113  std::unordered_map<int32_t, int32_t> copies;
114  HGCalParameters::layer_map copiesInLayers(layers.size() + 1);
115  std::vector<int32_t> wafer2copy;
116  std::vector<HGCalGeomParameters::cellParameters> wafers;
117  std::string attribute = "Volume";
118  DDValue val1(attribute, sdTag2, 0.0);
119  DDSpecificsMatchesValueFilter filter1{val1};
120  DDFilteredView fv1(*cpv, filter1);
121  bool ok = fv1.firstChild();
122  if (!ok) {
123  throw cms::Exception("DDException") << "Attribute " << val1 << " not found but needed.";
124  } else {
125  dodet = true;
126  std::unordered_set<std::string> names;
127  while (dodet) {
128  const DDSolid& sol = fv1.logicalPart().solid();
129  const std::string& name = fv1.logicalPart().name().name();
130  std::vector<int> copy = fv1.copyNumbers();
131  int nsiz = static_cast<int>(copy.size());
132  int wafer = (nsiz > 0) ? copy[nsiz - 1] : 0;
133  int layer = (nsiz > 1) ? copy[nsiz - 2] : 0;
134  if (nsiz < 2) {
135  throw cms::Exception("DDException") << "Funny wafer # " << wafer << " in " << nsiz << " components";
136  } else if (layer > static_cast<int>(layers.size())) {
137  edm::LogWarning("HGCalGeom") << "Funny wafer # " << wafer << " Layer " << layer << ":" << layers.size()
138  << " among " << nsiz << " components";
139  } else {
140  auto itr = copies.find(wafer);
141  auto cpy = copiesInLayers[layer].find(wafer);
142  if (itr != copies.end() && cpy == copiesInLayers[layer].end()) {
143  copiesInLayers[layer][wafer] = itr->second;
144  }
145  if (itr == copies.end()) {
146  copies[wafer] = wafer2copy.size();
147  copiesInLayers[layer][wafer] = wafer2copy.size();
148  double xx = HGCalParameters::k_ScaleFromDDD * fv1.translation().X();
149  if (std::abs(xx) < tolerance)
150  xx = 0;
151  double yy = HGCalParameters::k_ScaleFromDDD * fv1.translation().Y();
152  if (std::abs(yy) < tolerance)
153  yy = 0;
154  wafer2copy.emplace_back(wafer);
156  HGCalGeomParameters::cellParameters cell(false, wafer, p);
157  wafers.emplace_back(cell);
158  if (names.count(name) == 0) {
159  std::vector<double> zv, rv;
161  const DDPolyhedra& polyhedra = static_cast<DDPolyhedra>(sol);
162  zv = polyhedra.zVec();
163  rv = polyhedra.rMaxVec();
164  } else {
165  const DDExtrudedPolygon& polygon = static_cast<DDExtrudedPolygon>(sol);
166  zv = polygon.zVec();
167  rv = polygon.xVec();
168  }
171  double dz = 0.5 * HGCalParameters::k_ScaleFromDDDToG4 * (zv[1] - zv[0]);
172 #ifdef EDM_ML_DEBUG
173  edm::LogVerbatim("HGCalGeom")
174  << "Mode " << mode << " R " << php.waferSize_ << ":" << php.waferR_ << " z " << dz;
175 #endif
177  mytr.lay = 1;
178  mytr.bl = php.waferR_;
179  mytr.tl = php.waferR_;
180  mytr.h = php.waferR_;
181  mytr.dz = dz;
182  mytr.alpha = 0.0;
183  mytr.cellSize = waferSize_;
184  php.fillModule(mytr, false);
185  names.insert(name);
186  }
187  }
188  }
189  dodet = fv1.next();
190  }
191  }
192 
193  // Finally the cells
194  std::map<int, int> wafertype;
195  std::map<int, HGCalGeomParameters::cellParameters> cellsf, cellsc;
196  DDValue val2(attribute, sdTag3, 0.0);
197  DDSpecificsMatchesValueFilter filter2{val2};
198  DDFilteredView fv2(*cpv, filter2);
199  ok = fv2.firstChild();
200  if (!ok) {
201  throw cms::Exception("DDException") << "Attribute " << val2 << " not found but needed.";
202  } else {
203  dodet = true;
204  while (dodet) {
205  const DDSolid& sol = fv2.logicalPart().solid();
206  const std::string& name = sol.name().name();
207  std::vector<int> copy = fv2.copyNumbers();
208  int nsiz = static_cast<int>(copy.size());
209  int cellx = (nsiz > 0) ? copy[nsiz - 1] : 0;
210  int wafer = (nsiz > 1) ? copy[nsiz - 2] : 0;
211  int cell = HGCalTypes::getUnpackedCell6(cellx);
213  if (type != 1 && type != 2) {
214  throw cms::Exception("DDException")
215  << "Funny cell # " << cell << " type " << type << " in " << nsiz << " components";
216  } else {
217  auto ktr = wafertype.find(wafer);
218  if (ktr == wafertype.end())
219  wafertype[wafer] = type;
220  bool newc(false);
221  std::map<int, HGCalGeomParameters::cellParameters>::iterator itr;
222  double cellsize = php.cellSize_[0];
223  if (type == 1) {
224  itr = cellsf.find(cell);
225  newc = (itr == cellsf.end());
226  } else {
227  itr = cellsc.find(cell);
228  newc = (itr == cellsc.end());
229  cellsize = php.cellSize_[1];
230  }
231  if (newc) {
232  bool half = (name.find("Half") != std::string::npos);
233  double xx = HGCalParameters::k_ScaleFromDDD * fv2.translation().X();
234  double yy = HGCalParameters::k_ScaleFromDDD * fv2.translation().Y();
235  if (half) {
236  math::XYZPointD p1(-2.0 * cellsize / 9.0, 0, 0);
237  math::XYZPointD p2 = fv2.rotation()(p1);
240 #ifdef EDM_ML_DEBUG
241  if (std::abs(p2.X()) < HGCalParameters::tol)
242  p2.SetX(0.0);
243  if (std::abs(p2.Z()) < HGCalParameters::tol)
244  p2.SetZ(0.0);
245  edm::LogVerbatim("HGCalGeom") << "Wafer " << wafer << " Type " << type << " Cell " << cellx << " local "
246  << xx << ":" << yy << " new " << p1 << ":" << p2;
247 #endif
248  }
250  if (type == 1) {
251  cellsf[cell] = cp;
252  } else {
253  cellsc[cell] = cp;
254  }
255  }
256  }
257  dodet = fv2.next();
258  }
259  }
260 
262  layers, trforms, trformUse, copies, copiesInLayers, wafer2copy, wafers, wafertype, cellsf, cellsc, php);
263 }
264 
266  HGCalParameters& php,
267  const std::string& sdTag1,
268  const std::string& sdTag2,
269  const std::string& sdTag3,
271  const cms::DDFilter filter("Volume", sdTag1);
272  cms::DDFilteredView fv((*cpv), filter);
273  std::map<int, HGCalGeomParameters::layerParameters> layers;
274  std::vector<HGCalParameters::hgtrform> trforms;
275  std::vector<bool> trformUse;
276  std::vector<std::pair<int, int> > trused;
277 
278  while (fv.firstChild()) {
279  const std::vector<double>& pars = fv.parameters();
280  // Layers first
281  std::vector<int> copy = fv.copyNos();
282  int nsiz = static_cast<int>(copy.size());
283  int lay = (nsiz > 0) ? copy[0] : 0;
284  int zp = (nsiz > 2) ? copy[2] : -1;
285  if (zp != 1)
286  zp = -1;
287  if (lay == 0) {
288  throw cms::Exception("DDException") << "Funny layer # " << lay << " zp " << zp << " in " << nsiz << " components";
289  } else {
290  if (std::find(php.layer_.begin(), php.layer_.end(), lay) == php.layer_.end())
291  php.layer_.emplace_back(lay);
292  auto itr = layers.find(lay);
294  if (itr == layers.end()) {
295  double rin(0), rout(0);
296  if (dd4hep::isA<dd4hep::Polyhedra>(fv.solid())) {
297  rin = 0.5 * HGCalParameters::k_ScaleFromDD4hep * (pars[5] + pars[8]);
298  rout = 0.5 * HGCalParameters::k_ScaleFromDD4hep * (pars[6] + pars[9]);
299  } else if (dd4hep::isA<dd4hep::Tube>(fv.solid())) {
300  dd4hep::Tube tubeSeg(fv.solid());
301  rin = HGCalParameters::k_ScaleFromDD4hep * tubeSeg.rMin();
302  rout = HGCalParameters::k_ScaleFromDD4hep * tubeSeg.rMax();
303  }
304  HGCalGeomParameters::layerParameters laypar(rin, rout, zz);
305  layers[lay] = laypar;
306  }
307  std::pair<int, int> layz(lay, zp);
308  if (std::find(trused.begin(), trused.end(), layz) == trused.end()) {
309  trused.emplace_back(layz);
310  DD3Vector x, y, z;
311  fv.rotation().GetComponents(x, y, z);
312  const CLHEP::HepRep3x3 rotation(x.X(), y.X(), z.X(), x.Y(), y.Y(), z.Y(), x.Z(), y.Z(), z.Z());
313  const CLHEP::HepRotation hr(rotation);
315  if (std::abs(xx) < tolerance)
316  xx = 0;
318  if (std::abs(yy) < tolerance)
319  yy = 0;
321  const CLHEP::Hep3Vector h3v(xx, yy, zz);
323  mytrf.zp = zp;
324  mytrf.lay = lay;
325  mytrf.sec = 0;
326  mytrf.subsec = 0;
327  mytrf.h3v = h3v;
328  mytrf.hr = hr;
329  trforms.emplace_back(mytrf);
330  trformUse.emplace_back(false);
331  }
332  }
333  }
334 
335  // Then wafers
336  // This assumes layers are build starting from 1 (which on 25 Jan 2016, they
337  // were) to ensure that new copy numbers are always added to the end of the
338  // list.
339  std::unordered_map<int32_t, int32_t> copies;
340  HGCalParameters::layer_map copiesInLayers(layers.size() + 1);
341  std::vector<int32_t> wafer2copy;
342  std::vector<HGCalGeomParameters::cellParameters> wafers;
343  const cms::DDFilter filter1("Volume", sdTag2);
344  cms::DDFilteredView fv1((*cpv), filter1);
345  bool ok = fv1.firstChild();
346  if (!ok) {
347  throw cms::Exception("DDException") << "Attribute " << sdTag2 << " not found but needed.";
348  } else {
349  bool dodet = true;
350  std::unordered_set<std::string> names;
351  while (dodet) {
352  const std::string name = static_cast<std::string>(fv1.name());
353  std::vector<int> copy = fv1.copyNos();
354  int nsiz = static_cast<int>(copy.size());
355  int wafer = (nsiz > 0) ? copy[0] : 0;
356  int layer = (nsiz > 1) ? copy[1] : 0;
357  if (nsiz < 2) {
358  throw cms::Exception("DDException") << "Funny wafer # " << wafer << " in " << nsiz << " components";
359  } else if (layer > static_cast<int>(layers.size())) {
360  edm::LogWarning("HGCalGeom") << "Funny wafer # " << wafer << " Layer " << layer << ":" << layers.size()
361  << " among " << nsiz << " components";
362  } else {
363  auto itr = copies.find(wafer);
364  auto cpy = copiesInLayers[layer].find(wafer);
365  if (itr != copies.end() && cpy == copiesInLayers[layer].end()) {
366  copiesInLayers[layer][wafer] = itr->second;
367  }
368  if (itr == copies.end()) {
369  copies[wafer] = wafer2copy.size();
370  copiesInLayers[layer][wafer] = wafer2copy.size();
372  if (std::abs(xx) < tolerance)
373  xx = 0;
375  if (std::abs(yy) < tolerance)
376  yy = 0;
377  wafer2copy.emplace_back(wafer);
379  HGCalGeomParameters::cellParameters cell(false, wafer, p);
380  wafers.emplace_back(cell);
381  if (names.count(name) == 0) {
382  double zv[2], rv;
383  const std::vector<double>& pars = fv1.parameters();
385  zv[0] = pars[4];
386  zv[1] = pars[7];
387  rv = pars[6];
388  } else {
389  zv[0] = pars[3];
390  zv[1] = pars[9];
391  rv = pars[4];
392  }
395  double dz = 0.5 * HGCalParameters::k_ScaleFromDD4hepToG4 * (zv[1] - zv[0]);
396 #ifdef EDM_ML_DEBUG
397  edm::LogVerbatim("HGCalGeom")
398  << "Mode " << mode << " R " << php.waferSize_ << ":" << php.waferR_ << " z " << dz;
399 #endif
401  mytr.lay = 1;
402  mytr.bl = php.waferR_;
403  mytr.tl = php.waferR_;
404  mytr.h = php.waferR_;
405  mytr.dz = dz;
406  mytr.alpha = 0.0;
407  mytr.cellSize = waferSize_;
408  php.fillModule(mytr, false);
409  names.insert(name);
410  }
411  }
412  }
413  dodet = fv1.firstChild();
414  }
415  }
416 
417  // Finally the cells
418  std::map<int, int> wafertype;
419  std::map<int, HGCalGeomParameters::cellParameters> cellsf, cellsc;
420  const cms::DDFilter filter2("Volume", sdTag3);
421  cms::DDFilteredView fv2((*cpv), filter2);
422  ok = fv2.firstChild();
423  if (!ok) {
424  throw cms::Exception("DDException") << "Attribute " << sdTag3 << " not found but needed.";
425  } else {
426  bool dodet = true;
427  while (dodet) {
428  const std::string name = static_cast<std::string>(fv2.name());
429  std::vector<int> copy = fv2.copyNos();
430  int nsiz = static_cast<int>(copy.size());
431  int cellx = (nsiz > 0) ? copy[0] : 0;
432  int wafer = (nsiz > 1) ? copy[1] : 0;
433  int cell = HGCalTypes::getUnpackedCell6(cellx);
435  if (type != 1 && type != 2) {
436  throw cms::Exception("DDException")
437  << "Funny cell # " << cell << " type " << type << " in " << nsiz << " components";
438  } else {
439  auto ktr = wafertype.find(wafer);
440  if (ktr == wafertype.end())
441  wafertype[wafer] = type;
442  bool newc(false);
443  std::map<int, HGCalGeomParameters::cellParameters>::iterator itr;
444  double cellsize = php.cellSize_[0];
445  if (type == 1) {
446  itr = cellsf.find(cell);
447  newc = (itr == cellsf.end());
448  } else {
449  itr = cellsc.find(cell);
450  newc = (itr == cellsc.end());
451  cellsize = php.cellSize_[1];
452  }
453  if (newc) {
454  bool half = (name.find("Half") != std::string::npos);
457  if (half) {
458  math::XYZPointD p1(-2.0 * cellsize / 9.0, 0, 0);
459  math::XYZPointD p2 = fv2.rotation()(p1);
462 #ifdef EDM_ML_DEBUG
463  if (std::abs(p2.X()) < HGCalParameters::tol)
464  p2.SetX(0.0);
465  if (std::abs(p2.Z()) < HGCalParameters::tol)
466  p2.SetZ(0.0);
467  edm::LogVerbatim("HGCalGeom") << "Wafer " << wafer << " Type " << type << " Cell " << cellx << " local "
468  << xx << ":" << yy << " new " << p1 << ":" << p2;
469 #endif
470  }
472  if (type == 1) {
473  cellsf[cell] = cp;
474  } else {
475  cellsc[cell] = cp;
476  }
477  }
478  }
479  dodet = fv2.firstChild();
480  }
481  }
482 
484  layers, trforms, trformUse, copies, copiesInLayers, wafer2copy, wafers, wafertype, cellsf, cellsc, php);
485 }
486 
487 void HGCalGeomParameters::loadGeometryHexagon(const std::map<int, HGCalGeomParameters::layerParameters>& layers,
488  std::vector<HGCalParameters::hgtrform>& trforms,
489  std::vector<bool>& trformUse,
490  const std::unordered_map<int32_t, int32_t>& copies,
491  const HGCalParameters::layer_map& copiesInLayers,
492  const std::vector<int32_t>& wafer2copy,
493  const std::vector<HGCalGeomParameters::cellParameters>& wafers,
494  const std::map<int, int>& wafertype,
495  const std::map<int, HGCalGeomParameters::cellParameters>& cellsf,
496  const std::map<int, HGCalGeomParameters::cellParameters>& cellsc,
497  HGCalParameters& php) {
498  if (((cellsf.size() + cellsc.size()) == 0) || (wafers.empty()) || (layers.empty())) {
499  throw cms::Exception("DDException") << "HGCalGeomParameters: mismatch between geometry and specpar: cells "
500  << cellsf.size() << ":" << cellsc.size() << " wafers " << wafers.size()
501  << " layers " << layers.size();
502  }
503 
504  for (unsigned int i = 0; i < layers.size(); ++i) {
505  for (auto& layer : layers) {
506  if (layer.first == static_cast<int>(i + php.firstLayer_)) {
507  php.layerIndex_.emplace_back(i);
508  php.rMinLayHex_.emplace_back(layer.second.rmin);
509  php.rMaxLayHex_.emplace_back(layer.second.rmax);
510  php.zLayerHex_.emplace_back(layer.second.zpos);
511  break;
512  }
513  }
514  }
515 
516  for (unsigned int i = 0; i < php.layer_.size(); ++i) {
517  for (unsigned int i1 = 0; i1 < trforms.size(); ++i1) {
518  if (!trformUse[i1] && php.layerGroup_[trforms[i1].lay - 1] == static_cast<int>(i + 1)) {
519  trforms[i1].h3v *= static_cast<double>(HGCalParameters::k_ScaleFromDDD);
520  trforms[i1].lay = (i + 1);
521  trformUse[i1] = true;
522  php.fillTrForm(trforms[i1]);
523  int nz(1);
524  for (unsigned int i2 = i1 + 1; i2 < trforms.size(); ++i2) {
525  if (!trformUse[i2] && trforms[i2].zp == trforms[i1].zp &&
526  php.layerGroup_[trforms[i2].lay - 1] == static_cast<int>(i + 1)) {
527  php.addTrForm(trforms[i2].h3v);
528  nz++;
529  trformUse[i2] = true;
530  }
531  }
532  if (nz > 0) {
533  php.scaleTrForm(double(1.0 / nz));
534  }
535  }
536  }
537  }
538 
539  double rmin = HGCalParameters::k_ScaleFromDDD * php.waferR_;
540  for (unsigned i = 0; i < wafer2copy.size(); ++i) {
541  php.waferCopy_.emplace_back(wafer2copy[i]);
542  php.waferPosX_.emplace_back(wafers[i].xyz.x());
543  php.waferPosY_.emplace_back(wafers[i].xyz.y());
544  auto ktr = wafertype.find(wafer2copy[i]);
545  int typet = (ktr == wafertype.end()) ? 0 : (ktr->second);
546  php.waferTypeT_.emplace_back(typet);
547  double r = wafers[i].xyz.perp();
548  int type(3);
549  for (int k = 1; k < 4; ++k) {
550  if ((r + rmin) <= php.boundR_[k]) {
551  type = k;
552  break;
553  }
554  }
555  php.waferTypeL_.emplace_back(type);
556  }
557  php.copiesInLayers_ = copiesInLayers;
558  php.nSectors_ = static_cast<int>(php.waferCopy_.size());
559 
560  std::vector<HGCalGeomParameters::cellParameters>::const_iterator itrf = wafers.end();
561  for (unsigned int i = 0; i < cellsf.size(); ++i) {
562  auto itr = cellsf.find(i);
563  if (itr == cellsf.end()) {
564  throw cms::Exception("DDException") << "HGCalGeomParameters: missing info for fine cell number " << i;
565  } else {
566  double xx = (itr->second).xyz.x();
567  double yy = (itr->second).xyz.y();
568  int waf = (itr->second).wafer;
569  std::pair<double, double> xy = cellPosition(wafers, itrf, waf, xx, yy);
570  php.cellFineX_.emplace_back(xy.first);
571  php.cellFineY_.emplace_back(xy.second);
572  php.cellFineHalf_.emplace_back((itr->second).half);
573  }
574  }
575  itrf = wafers.end();
576  for (unsigned int i = 0; i < cellsc.size(); ++i) {
577  auto itr = cellsc.find(i);
578  if (itr == cellsc.end()) {
579  throw cms::Exception("DDException") << "HGCalGeomParameters: missing info for coarse cell number " << i;
580  } else {
581  double xx = (itr->second).xyz.x();
582  double yy = (itr->second).xyz.y();
583  int waf = (itr->second).wafer;
584  std::pair<double, double> xy = cellPosition(wafers, itrf, waf, xx, yy);
585  php.cellCoarseX_.emplace_back(xy.first);
586  php.cellCoarseY_.emplace_back(xy.second);
587  php.cellCoarseHalf_.emplace_back((itr->second).half);
588  }
589  }
590  int depth(0);
591  for (unsigned int i = 0; i < php.layerGroup_.size(); ++i) {
592  bool first(true);
593  for (unsigned int k = 0; k < php.layerGroup_.size(); ++k) {
594  if (php.layerGroup_[k] == static_cast<int>(i + 1)) {
595  if (first) {
596  php.depth_.emplace_back(i + 1);
597  php.depthIndex_.emplace_back(depth);
598  php.depthLayerF_.emplace_back(k);
599  ++depth;
600  first = false;
601  }
602  }
603  }
604  }
605  HGCalParameters::hgtrap mytr = php.getModule(0, false);
611  double dz = mytr.dz;
612  php.fillModule(mytr, true);
613  mytr.dz = 2 * dz;
614  php.fillModule(mytr, true);
615  mytr.dz = 3 * dz;
616  php.fillModule(mytr, true);
617 #ifdef EDM_ML_DEBUG
618  edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters finds " << php.zLayerHex_.size() << " layers";
619  for (unsigned int i = 0; i < php.zLayerHex_.size(); ++i) {
620  int k = php.layerIndex_[i];
621  edm::LogVerbatim("HGCalGeom") << "Layer[" << i << ":" << k << ":" << php.layer_[k]
622  << "] with r = " << php.rMinLayHex_[i] << ":" << php.rMaxLayHex_[i]
623  << " at z = " << php.zLayerHex_[i];
624  }
625  edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters has " << php.depthIndex_.size() << " depths";
626  for (unsigned int i = 0; i < php.depthIndex_.size(); ++i) {
627  int k = php.depthIndex_[i];
628  edm::LogVerbatim("HGCalGeom") << "Reco Layer[" << i << ":" << k << "] First Layer " << php.depthLayerF_[i]
629  << " Depth " << php.depth_[k];
630  }
631  edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters finds " << php.nSectors_ << " wafers";
632  for (unsigned int i = 0; i < php.waferCopy_.size(); ++i)
633  edm::LogVerbatim("HGCalGeom") << "Wafer[" << i << ": " << php.waferCopy_[i] << "] type " << php.waferTypeL_[i]
634  << ":" << php.waferTypeT_[i] << " at (" << php.waferPosX_[i] << ","
635  << php.waferPosY_[i] << ",0)";
636  edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: wafer radius " << php.waferR_ << " and dimensions of the "
637  << "wafers:";
638  edm::LogVerbatim("HGCalGeom") << "Sim[0] " << php.moduleLayS_[0] << " dx " << php.moduleBlS_[0] << ":"
639  << php.moduleTlS_[0] << " dy " << php.moduleHS_[0] << " dz " << php.moduleDzS_[0]
640  << " alpha " << php.moduleAlphaS_[0];
641  for (unsigned int k = 0; k < php.moduleLayR_.size(); ++k)
642  edm::LogVerbatim("HGCalGeom") << "Rec[" << k << "] " << php.moduleLayR_[k] << " dx " << php.moduleBlR_[k] << ":"
643  << php.moduleTlR_[k] << " dy " << php.moduleHR_[k] << " dz " << php.moduleDzR_[k]
644  << " alpha " << php.moduleAlphaR_[k];
645  edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters finds " << php.cellFineX_.size() << " fine cells in a wafer";
646  for (unsigned int i = 0; i < php.cellFineX_.size(); ++i)
647  edm::LogVerbatim("HGCalGeom") << "Fine Cell[" << i << "] at (" << php.cellFineX_[i] << "," << php.cellFineY_[i]
648  << ",0)";
649  edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters finds " << php.cellCoarseX_.size()
650  << " coarse cells in a wafer";
651  for (unsigned int i = 0; i < php.cellCoarseX_.size(); ++i)
652  edm::LogVerbatim("HGCalGeom") << "Coarse Cell[" << i << "] at (" << php.cellCoarseX_[i] << ","
653  << php.cellCoarseY_[i] << ",0)";
654  edm::LogVerbatim("HGCalGeom") << "Obtained " << php.trformIndex_.size() << " transformation matrices";
655  for (unsigned int k = 0; k < php.trformIndex_.size(); ++k) {
656  edm::LogVerbatim("HGCalGeom") << "Matrix[" << k << "] (" << std::hex << php.trformIndex_[k] << std::dec
657  << ") Translation (" << php.trformTranX_[k] << ", " << php.trformTranY_[k] << ", "
658  << php.trformTranZ_[k] << " Rotation (" << php.trformRotXX_[k] << ", "
659  << php.trformRotYX_[k] << ", " << php.trformRotZX_[k] << ", " << php.trformRotXY_[k]
660  << ", " << php.trformRotYY_[k] << ", " << php.trformRotZY_[k] << ", "
661  << php.trformRotXZ_[k] << ", " << php.trformRotYZ_[k] << ", " << php.trformRotZZ_[k]
662  << ")";
663  }
664  edm::LogVerbatim("HGCalGeom") << "Dump copiesInLayers for " << php.copiesInLayers_.size() << " layers";
665  for (unsigned int k = 0; k < php.copiesInLayers_.size(); ++k) {
666  const auto& theModules = php.copiesInLayers_[k];
667  edm::LogVerbatim("HGCalGeom") << "Layer " << k << ":" << theModules.size();
668  int k2(0);
669  for (std::unordered_map<int, int>::const_iterator itr = theModules.begin(); itr != theModules.end(); ++itr, ++k2) {
670  edm::LogVerbatim("HGCalGeom") << "[" << k2 << "] " << itr->first << ":" << itr->second;
671  }
672  }
673 #endif
674 }
675 
677  DDFilteredView fv = _fv;
678  bool dodet(true);
679  std::map<int, HGCalGeomParameters::layerParameters> layers;
680  std::map<std::pair<int, int>, HGCalParameters::hgtrform> trforms;
681  int levelTop = 3 + std::max(php.levelT_[0], php.levelT_[1]);
682 #ifdef EDM_ML_DEBUG
683  int ntot(0);
684 #endif
685  while (dodet) {
686 #ifdef EDM_ML_DEBUG
687  ++ntot;
688 #endif
689  std::vector<int> copy = fv.copyNumbers();
690  int nsiz = static_cast<int>(copy.size());
691  if (nsiz < levelTop) {
692  int lay = copy[nsiz - 1];
693  int zside = (nsiz > php.levelZSide_) ? copy[php.levelZSide_] : -1;
694  if (zside != 1)
695  zside = -1;
696  const DDSolid& sol = fv.logicalPart().solid();
697 #ifdef EDM_ML_DEBUG
698  edm::LogVerbatim("HGCalGeom") << sol.name() << " shape " << sol.shape() << " size " << nsiz << ":" << levelTop
699  << " lay " << lay << " z " << zside;
700 #endif
701  if (lay == 0) {
702  throw cms::Exception("DDException")
703  << "Funny layer # " << lay << " zp " << zside << " in " << nsiz << " components";
704  } else if (sol.shape() == DDSolidShape::ddtubs) {
705  if (std::find(php.layer_.begin(), php.layer_.end(), lay) == php.layer_.end())
706  php.layer_.emplace_back(lay);
707  const DDTubs& tube = static_cast<DDTubs>(sol);
708  double rin = HGCalParameters::k_ScaleFromDDD * tube.rIn();
709  double rout = HGCalParameters::k_ScaleFromDDD * tube.rOut();
710  auto itr = layers.find(lay);
711  if (itr == layers.end()) {
712  double zp = HGCalParameters::k_ScaleFromDDD * fv.translation().Z();
713  HGCalGeomParameters::layerParameters laypar(rin, rout, zp);
714  layers[lay] = laypar;
715  } else {
716  (itr->second).rmin = std::min(rin, (itr->second).rmin);
717  (itr->second).rmax = std::max(rout, (itr->second).rmax);
718  }
719  if (trforms.find(std::make_pair(lay, zside)) == trforms.end()) {
720  DD3Vector x, y, z;
721  fv.rotation().GetComponents(x, y, z);
722  const CLHEP::HepRep3x3 rotation(x.X(), y.X(), z.X(), x.Y(), y.Y(), z.Y(), x.Z(), y.Z(), z.Z());
723  const CLHEP::HepRotation hr(rotation);
724  double xx =
725  ((std::abs(fv.translation().X()) < tolerance) ? 0
727  double yy =
728  ((std::abs(fv.translation().Y()) < tolerance) ? 0
730  const CLHEP::Hep3Vector h3v(xx, yy, HGCalParameters::k_ScaleFromDDD * fv.translation().Z());
732  mytrf.zp = zside;
733  mytrf.lay = lay;
734  mytrf.sec = 0;
735  mytrf.subsec = 0;
736  mytrf.h3v = h3v;
737  mytrf.hr = hr;
738  trforms[std::make_pair(lay, zside)] = mytrf;
739  }
740  }
741  }
742  dodet = fv.next();
743  }
744 #ifdef EDM_ML_DEBUG
745  edm::LogVerbatim("HGCalGeom") << "Total # of views " << ntot;
746 #endif
747  loadGeometryHexagon8(layers, trforms, firstLayer, php);
748 }
749 
751  HGCalParameters& php,
752  const std::string& sdTag1,
753  int firstLayer) {
754  const cms::DDFilter filter("Volume", sdTag1);
755  cms::DDFilteredView fv((*cpv), filter);
756  std::map<int, HGCalGeomParameters::layerParameters> layers;
757  std::map<std::pair<int, int>, HGCalParameters::hgtrform> trforms;
758  int levelTop = 3 + std::max(php.levelT_[0], php.levelT_[1]);
759 #ifdef EDM_ML_DEBUG
760  int ntot(0);
761 #endif
762  while (fv.firstChild()) {
763 #ifdef EDM_ML_DEBUG
764  ++ntot;
765 #endif
766  // Layers first
767  int nsiz = static_cast<int>(fv.level());
768  if (nsiz < levelTop) {
769  std::vector<int> copy = fv.copyNos();
770  int lay = copy[0];
771  int zside = (nsiz > php.levelZSide_) ? copy[nsiz - php.levelZSide_ - 1] : -1;
772  if (zside != 1)
773  zside = -1;
774 #ifdef EDM_ML_DEBUG
775  edm::LogVerbatim("HGCalGeom") << fv.name() << " shape " << cms::dd::name(cms::DDSolidShapeMap, fv.shape())
776  << " size " << nsiz << ":" << levelTop << " lay " << lay << " z " << zside << ":"
777  << php.levelZSide_;
778 #endif
779  if (lay == 0) {
780  throw cms::Exception("DDException")
781  << "Funny layer # " << lay << " zp " << zside << " in " << nsiz << " components";
782  } else if (fv.shape() == cms::DDSolidShape::ddtubs) {
783  if (std::find(php.layer_.begin(), php.layer_.end(), lay) == php.layer_.end())
784  php.layer_.emplace_back(lay);
785  const std::vector<double>& pars = fv.parameters();
786  double rin = HGCalParameters::k_ScaleFromDD4hep * pars[0];
787  double rout = HGCalParameters::k_ScaleFromDD4hep * pars[1];
788  auto itr = layers.find(lay);
789  if (itr == layers.end()) {
790  double zp = HGCalParameters::k_ScaleFromDD4hep * fv.translation().Z();
791  HGCalGeomParameters::layerParameters laypar(rin, rout, zp);
792  layers[lay] = laypar;
793  } else {
794  (itr->second).rmin = std::min(rin, (itr->second).rmin);
795  (itr->second).rmax = std::max(rout, (itr->second).rmax);
796  }
797  if (trforms.find(std::make_pair(lay, zside)) == trforms.end()) {
798  DD3Vector x, y, z;
799  fv.rotation().GetComponents(x, y, z);
800  const CLHEP::HepRep3x3 rotation(x.X(), y.X(), z.X(), x.Y(), y.Y(), z.Y(), x.Z(), y.Z(), z.Z());
801  const CLHEP::HepRotation hr(rotation);
802  double xx = ((std::abs(fv.translation().X()) < tolerance)
803  ? 0
805  double yy = ((std::abs(fv.translation().Y()) < tolerance)
806  ? 0
808  const CLHEP::Hep3Vector h3v(xx, yy, HGCalParameters::k_ScaleFromDD4hep * fv.translation().Z());
810  mytrf.zp = zside;
811  mytrf.lay = lay;
812  mytrf.sec = 0;
813  mytrf.subsec = 0;
814  mytrf.h3v = h3v;
815  mytrf.hr = hr;
816  trforms[std::make_pair(lay, zside)] = mytrf;
817  }
818  }
819  }
820  }
821 #ifdef EDM_ML_DEBUG
822  edm::LogVerbatim("HGCalGeom") << "Total # of views " << ntot;
823 #endif
824  loadGeometryHexagon8(layers, trforms, firstLayer, php);
825 }
826 
828  HGCalParameters& php,
829  const std::string& sdTag1,
830  const std::string& sdTag2,
831  int firstLayer) {
832 #ifdef EDM_ML_DEBUG
833  edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters (DDD)::loadGeometryHexagonModule called with tags " << sdTag1
834  << ":" << sdTag2 << " firstLayer " << firstLayer << ":" << php.firstMixedLayer_;
835  int ntot1(0), ntot2(0);
836 #endif
837  std::map<int, HGCalGeomParameters::layerParameters> layers;
838  std::map<std::pair<int, int>, double> zvals;
839  std::map<std::pair<int, int>, HGCalParameters::hgtrform> trforms;
840  int levelTop = php.levelT_[0];
841 
842  std::string attribute = "Volume";
843  DDValue val1(attribute, sdTag2, 0.0);
844  DDSpecificsMatchesValueFilter filter1{val1};
845  DDFilteredView fv1(*cpv, filter1);
846  bool dodet = fv1.firstChild();
847  while (dodet) {
848 #ifdef EDM_ML_DEBUG
849  ++ntot1;
850 #endif
851  std::vector<int> copy = fv1.copyNumbers();
852  int nsiz = static_cast<int>(copy.size());
853  if (levelTop < nsiz) {
854  int lay = copy[levelTop];
855  int zside = (nsiz > php.levelZSide_) ? copy[php.levelZSide_] : -1;
856  if (zside != 1)
857  zside = -1;
858  if (lay == 0) {
859  throw cms::Exception("DDException")
860  << "Funny layer # " << lay << " zp " << zside << " in " << nsiz << " components";
861  } else {
862  if (zvals.find(std::make_pair(lay, zside)) == zvals.end()) {
863  zvals[std::make_pair(lay, zside)] = HGCalParameters::k_ScaleFromDDD * fv1.translation().Z();
864 #ifdef EDM_ML_DEBUG
865  std::ostringstream st1;
866  st1 << "Name0 " << fv1.name() << " LTop " << levelTop << ":" << lay << " ZSide " << zside << " # of levels "
867  << nsiz;
868  for (const auto& c : copy)
869  st1 << ":" << c;
870  st1 << " Z " << zvals[std::make_pair(lay, zside)];
871  edm::LogVerbatim("HGCalGeom") << st1.str();
872 #endif
873  }
874  }
875  }
876  dodet = fv1.next();
877  }
878 
879  DDValue val2(attribute, sdTag1, 0.0);
880  DDSpecificsMatchesValueFilter filter2{val2};
881  DDFilteredView fv2(*cpv, filter2);
882  dodet = fv2.firstChild();
883  while (dodet) {
884  std::vector<int> copy = fv2.copyNumbers();
885  int nsiz = static_cast<int>(copy.size());
886 #ifdef EDM_ML_DEBUG
887  ++ntot2;
888  edm::LogVerbatim("HGCalGeom") << "loadGeometryHexagonModule:: nsiz " << nsiz << " Ltop " << levelTop;
889 #endif
890  if (levelTop < nsiz) {
891  int lay = copy[levelTop];
892  int zside = (nsiz > php.levelZSide_) ? copy[php.levelZSide_] : -1;
893  if (zside != 1)
894  zside = -1;
895  const DDSolid& sol = fv2.logicalPart().solid();
896 #ifdef EDM_ML_DEBUG
897  std::ostringstream st2;
898  st2 << "Name1 " << sol.name() << " shape " << sol.shape() << " LTop " << levelTop << ":" << lay << " ZSide "
899  << zside << ":" << php.levelZSide_ << " # of levels " << nsiz;
900  for (const auto& c : copy)
901  st2 << ":" << c;
902  edm::LogVerbatim("HGCalGeom") << st2.str();
903 #endif
904  if (lay == 0) {
905  throw cms::Exception("DDException")
906  << "Funny layer # " << lay << " zp " << zside << " in " << nsiz << " components";
907  } else if ((sol.shape() == DDSolidShape::ddtubs) || (sol.shape() == DDSolidShape::ddbox)) {
908  if (zvals.find(std::make_pair(lay, zside)) != zvals.end()) {
909  if (std::find(php.layer_.begin(), php.layer_.end(), lay) == php.layer_.end())
910  php.layer_.emplace_back(lay);
911  auto itr = layers.find(lay);
912  if (itr == layers.end()) {
913  double rin(0), rout(0);
914  if (sol.shape() == DDSolidShape::ddtubs) {
915  const DDTubs& tube = static_cast<DDTubs>(sol);
916  rin = HGCalParameters::k_ScaleFromDDD * tube.rIn();
917  rout = (php.firstMixedLayer_ > 0 && lay >= php.firstMixedLayer_)
918  ? php.radiusMixBoundary_[lay - php.firstMixedLayer_]
920  } else {
921  const DDBox& box = static_cast<DDBox>(sol);
922  rout = HGCalParameters::k_ScaleFromDDD * box.halfX();
923  }
924  double zp = zvals[std::make_pair(lay, 1)];
925  HGCalGeomParameters::layerParameters laypar(rin, rout, zp);
926  layers[lay] = laypar;
927 #ifdef EDM_ML_DEBUG
928  std::ostringstream st3;
929  st3 << "Name1 " << fv2.name() << " LTop " << levelTop << ":" << lay << " ZSide " << zside << " # of levels "
930  << nsiz;
931  for (const auto& c : copy)
932  st3 << ":" << c;
933  st3 << " R " << rin << ":" << rout;
934  edm::LogVerbatim("HGCalGeom") << st3.str();
935 #endif
936  }
937 
938  if (trforms.find(std::make_pair(lay, zside)) == trforms.end()) {
939  DD3Vector x, y, z;
940  fv2.rotation().GetComponents(x, y, z);
941  const CLHEP::HepRep3x3 rotation(x.X(), y.X(), z.X(), x.Y(), y.Y(), z.Y(), x.Z(), y.Z(), z.Z());
942  const CLHEP::HepRotation hr(rotation);
943  double xx = ((std::abs(fv2.translation().X()) < tolerance)
944  ? 0
946  double yy = ((std::abs(fv2.translation().Y()) < tolerance)
947  ? 0
949  const CLHEP::Hep3Vector h3v(xx, yy, zvals[std::make_pair(lay, zside)]);
951  mytrf.zp = zside;
952  mytrf.lay = lay;
953  mytrf.sec = 0;
954  mytrf.subsec = 0;
955  mytrf.h3v = h3v;
956  mytrf.hr = hr;
957  trforms[std::make_pair(lay, zside)] = mytrf;
958 #ifdef EDM_ML_DEBUG
959  edm::LogVerbatim("HGCalGeom") << "Translation " << h3v;
960 #endif
961  }
962  }
963  }
964  }
965  dodet = fv2.next();
966  }
967 #ifdef EDM_ML_DEBUG
968  edm::LogVerbatim("HGCalGeom") << "Total # of views " << ntot1 << ":" << ntot2;
969 #endif
970  loadGeometryHexagon8(layers, trforms, firstLayer, php);
971 }
972 
974  HGCalParameters& php,
975  const std::string& sdTag1,
976  const std::string& sdTag2,
977  int firstLayer) {
978 #ifdef EDM_ML_DEBUG
979  edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters (DD4hep)::loadGeometryHexagonModule called with tags " << sdTag1
980  << ":" << sdTag2 << " firstLayer " << firstLayer;
981  int ntot1(0), ntot2(0);
982 #endif
983  std::map<int, HGCalGeomParameters::layerParameters> layers;
984  std::map<std::pair<int, int>, HGCalParameters::hgtrform> trforms;
985  std::map<std::pair<int, int>, double> zvals;
986  int levelTop = php.levelT_[0];
987 
988  const cms::DDFilter filter1("Volume", sdTag2);
989  cms::DDFilteredView fv1((*cpv), filter1);
990  while (fv1.firstChild()) {
991 #ifdef EDM_ML_DEBUG
992  ++ntot1;
993 #endif
994  int nsiz = static_cast<int>(fv1.level());
995  if (nsiz > levelTop) {
996  std::vector<int> copy = fv1.copyNos();
997  int lay = copy[nsiz - levelTop - 1];
998  int zside = (nsiz > php.levelZSide_) ? copy[nsiz - php.levelZSide_ - 1] : -1;
999  if (zside != 1)
1000  zside = -1;
1001  if (lay == 0) {
1002  throw cms::Exception("DDException")
1003  << "Funny layer # " << lay << " zp " << zside << " in " << nsiz << " components";
1004  } else {
1005  if (zvals.find(std::make_pair(lay, zside)) == zvals.end()) {
1006  zvals[std::make_pair(lay, zside)] = HGCalParameters::k_ScaleFromDD4hep * fv1.translation().Z();
1007 #ifdef EDM_ML_DEBUG
1008  std::ostringstream st1;
1009  st1 << "Name0 " << fv1.name() << " LTop " << levelTop << ":" << lay << " ZSide " << zside << " # of levels "
1010  << nsiz;
1011  for (const auto& c : copy)
1012  st1 << ":" << c;
1013  st1 << " Z " << zvals[std::make_pair(lay, zside)];
1014  edm::LogVerbatim("HGCalGeom") << st1.str();
1015 #endif
1016  }
1017  }
1018  }
1019  }
1020 
1021  const cms::DDFilter filter2("Volume", sdTag1);
1022  cms::DDFilteredView fv2((*cpv), filter2);
1023  while (fv2.firstChild()) {
1024  // Layers first
1025  int nsiz = static_cast<int>(fv2.level());
1026 #ifdef EDM_ML_DEBUG
1027  ++ntot2;
1028  edm::LogVerbatim("HGCalGeom") << "loadGeometryHexagonModule:: nsiz " << nsiz << " Ltop " << levelTop;
1029 #endif
1030  if (nsiz > levelTop) {
1031  std::vector<int> copy = fv2.copyNos();
1032  int lay = copy[nsiz - levelTop - 1];
1033  int zside = (nsiz > php.levelZSide_) ? copy[nsiz - php.levelZSide_ - 1] : -1;
1034  if (zside != 1)
1035  zside = -1;
1036 #ifdef EDM_ML_DEBUG
1037  std::ostringstream st2;
1038  st2 << "Name1 " << fv2.name() << "Shape " << cms::dd::name(cms::DDSolidShapeMap, fv2.shape()) << " LTop "
1039  << levelTop << ":" << lay << " ZSide " << zside << ":" << php.levelZSide_ << " # of levels " << nsiz;
1040  for (const auto& c : copy)
1041  st2 << ":" << c;
1042  edm::LogVerbatim("HGCalGeom") << st2.str();
1043 #endif
1044  if (lay == 0) {
1045  throw cms::Exception("DDException")
1046  << "Funny layer # " << lay << " zp " << zside << " in " << nsiz << " components";
1047  } else {
1048  if (zvals.find(std::make_pair(lay, zside)) != zvals.end()) {
1049  if (std::find(php.layer_.begin(), php.layer_.end(), lay) == php.layer_.end())
1050  php.layer_.emplace_back(lay);
1051  auto itr = layers.find(lay);
1052  if (itr == layers.end()) {
1053  const std::vector<double>& pars = fv2.parameters();
1054  double rin(0), rout(0);
1055  if (dd4hep::isA<dd4hep::Box>(fv2.solid())) {
1056  rout = HGCalParameters::k_ScaleFromDD4hep * pars[0];
1057  } else {
1058  rin = HGCalParameters::k_ScaleFromDD4hep * pars[0];
1059  rout = (php.firstMixedLayer_ > 0 && lay >= php.firstMixedLayer_)
1060  ? php.radiusMixBoundary_[lay - php.firstMixedLayer_]
1062  }
1063  double zp = zvals[std::make_pair(lay, 1)];
1064  HGCalGeomParameters::layerParameters laypar(rin, rout, zp);
1065  layers[lay] = laypar;
1066 #ifdef EDM_ML_DEBUG
1067  std::ostringstream st3;
1068  st3 << "Name2 " << fv2.name() << " LTop " << levelTop << ":" << lay << " ZSide " << zside << " # of levels "
1069  << nsiz;
1070  for (const auto& c : copy)
1071  st3 << ":" << c;
1072  st3 << " R " << rin << ":" << rout;
1073  edm::LogVerbatim("HGCalGeom") << st3.str();
1074 #endif
1075  }
1076 
1077  if (trforms.find(std::make_pair(lay, zside)) == trforms.end()) {
1078  DD3Vector x, y, z;
1079  fv2.rotation().GetComponents(x, y, z);
1080  const CLHEP::HepRep3x3 rotation(x.X(), y.X(), z.X(), x.Y(), y.Y(), z.Y(), x.Z(), y.Z(), z.Z());
1081  const CLHEP::HepRotation hr(rotation);
1082  double xx = ((std::abs(fv2.translation().X()) < tolerance)
1083  ? 0
1085  double yy = ((std::abs(fv2.translation().Y()) < tolerance)
1086  ? 0
1088  const CLHEP::Hep3Vector h3v(xx, yy, zvals[std::make_pair(lay, zside)]);
1090  mytrf.zp = zside;
1091  mytrf.lay = lay;
1092  mytrf.sec = 0;
1093  mytrf.subsec = 0;
1094  mytrf.h3v = h3v;
1095  mytrf.hr = hr;
1096  trforms[std::make_pair(lay, zside)] = mytrf;
1097 #ifdef EDM_ML_DEBUG
1098  edm::LogVerbatim("HGCalGeom") << "Translation " << h3v;
1099 #endif
1100  }
1101  }
1102  }
1103  }
1104  }
1105 #ifdef EDM_ML_DEBUG
1106  edm::LogVerbatim("HGCalGeom") << "Total # of views " << ntot1 << ":" << ntot2;
1107 #endif
1108  loadGeometryHexagon8(layers, trforms, firstLayer, php);
1109 }
1110 
1111 void HGCalGeomParameters::loadGeometryHexagon8(const std::map<int, HGCalGeomParameters::layerParameters>& layers,
1112  std::map<std::pair<int, int>, HGCalParameters::hgtrform>& trforms,
1113  const int& firstLayer,
1114  HGCalParameters& php) {
1115  double rmin(0), rmax(0);
1116  for (unsigned int i = 0; i < layers.size(); ++i) {
1117  for (auto& layer : layers) {
1118  if (layer.first == static_cast<int>(i + firstLayer)) {
1119  php.layerIndex_.emplace_back(i);
1120  php.rMinLayHex_.emplace_back(layer.second.rmin);
1121  php.rMaxLayHex_.emplace_back(layer.second.rmax);
1122  php.zLayerHex_.emplace_back(layer.second.zpos);
1123  if (i == 0) {
1124  rmin = layer.second.rmin;
1125  rmax = layer.second.rmax;
1126  } else {
1127  if (rmin > layer.second.rmin)
1128  rmin = layer.second.rmin;
1129  if (rmax < layer.second.rmax)
1130  rmax = layer.second.rmax;
1131  }
1132  break;
1133  }
1134  }
1135  }
1136  php.rLimit_.emplace_back(rmin);
1137  php.rLimit_.emplace_back(rmax);
1138  php.depth_ = php.layer_;
1139  php.depthIndex_ = php.layerIndex_;
1140  php.depthLayerF_ = php.layerIndex_;
1141 
1142  for (unsigned int i = 0; i < php.layer_.size(); ++i) {
1143  for (auto& trform : trforms) {
1144  if (trform.first.first == static_cast<int>(i + firstLayer)) {
1145  php.fillTrForm(trform.second);
1146  }
1147  }
1148  }
1149 
1150 #ifdef EDM_ML_DEBUG
1151  edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: Minimum/maximum R " << php.rLimit_[0] << ":" << php.rLimit_[1];
1152  edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters finds " << php.zLayerHex_.size() << " layers";
1153  for (unsigned int i = 0; i < php.zLayerHex_.size(); ++i) {
1154  int k = php.layerIndex_[i];
1155  edm::LogVerbatim("HGCalGeom") << "Layer[" << i << ":" << k << ":" << php.layer_[k]
1156  << "] with r = " << php.rMinLayHex_[i] << ":" << php.rMaxLayHex_[i]
1157  << " at z = " << php.zLayerHex_[i];
1158  }
1159  edm::LogVerbatim("HGCalGeom") << "Obtained " << php.trformIndex_.size() << " transformation matrices";
1160  for (unsigned int k = 0; k < php.trformIndex_.size(); ++k) {
1161  edm::LogVerbatim("HGCalGeom") << "Matrix[" << k << "] (" << std::hex << php.trformIndex_[k] << std::dec
1162  << ") Translation (" << php.trformTranX_[k] << ", " << php.trformTranY_[k] << ", "
1163  << php.trformTranZ_[k] << " Rotation (" << php.trformRotXX_[k] << ", "
1164  << php.trformRotYX_[k] << ", " << php.trformRotZX_[k] << ", " << php.trformRotXY_[k]
1165  << ", " << php.trformRotYY_[k] << ", " << php.trformRotZY_[k] << ", "
1166  << php.trformRotXZ_[k] << ", " << php.trformRotYZ_[k] << ", " << php.trformRotZZ_[k]
1167  << ")";
1168  }
1169 #endif
1170 }
1171 
1173  HGCalParameters& php,
1174  const DDCompactView* cpv,
1175  const std::string& sdTag1,
1176  const std::string& sdTag2) {
1178  php.boundR_ = getDDDArray("RadiusBound", sv, 4);
1180  php.rLimit_ = getDDDArray("RadiusLimits", sv, 2);
1182  php.levelT_ = dbl_to_int(getDDDArray("LevelTop", sv, 0));
1183 
1184  // Grouping of layers
1185  php.layerGroup_ = dbl_to_int(getDDDArray("GroupingZFine", sv, 0));
1186  php.layerGroupM_ = dbl_to_int(getDDDArray("GroupingZMid", sv, 0));
1187  php.layerGroupO_ = dbl_to_int(getDDDArray("GroupingZOut", sv, 0));
1188  php.slopeMin_ = getDDDArray("Slope", sv, 1);
1189  const auto& dummy2 = getDDDArray("LayerOffset", sv, 0);
1190  if (!dummy2.empty())
1191  php.layerOffset_ = dummy2[0];
1192  else
1193  php.layerOffset_ = 0;
1194 
1195  // Wafer size
1196  std::string attribute = "Volume";
1197  DDSpecificsMatchesValueFilter filter1{DDValue(attribute, sdTag1, 0.0)};
1198  DDFilteredView fv1(*cpv, filter1);
1199  if (fv1.firstChild()) {
1200  DDsvalues_type sv(fv1.mergedSpecifics());
1201  const auto& dummy = getDDDArray("WaferSize", sv, 0);
1202  waferSize_ = dummy[0];
1203  }
1204 
1205  // Cell size
1206  DDSpecificsMatchesValueFilter filter2{DDValue(attribute, sdTag2, 0.0)};
1207  DDFilteredView fv2(*cpv, filter2);
1208  if (fv2.firstChild()) {
1209  DDsvalues_type sv(fv2.mergedSpecifics());
1210  php.cellSize_ = getDDDArray("CellSize", sv, 0);
1211  }
1212 
1213  loadSpecParsHexagon(php);
1214 }
1215 
1217  HGCalParameters& php,
1218  const std::string& sdTag1,
1219  const std::string& sdTag2,
1220  const std::string& sdTag3,
1221  const std::string& sdTag4) {
1222  php.boundR_ = fv.get<std::vector<double> >(sdTag4, "RadiusBound");
1224  php.rLimit_ = fv.get<std::vector<double> >(sdTag4, "RadiusLimits");
1226  php.levelT_ = dbl_to_int(fv.get<std::vector<double> >(sdTag4, "LevelTop"));
1227 
1228  // Grouping of layers
1229  php.layerGroup_ = dbl_to_int(fv.get<std::vector<double> >(sdTag1, "GroupingZFine"));
1230  php.layerGroupM_ = dbl_to_int(fv.get<std::vector<double> >(sdTag1, "GroupingZMid"));
1231  php.layerGroupO_ = dbl_to_int(fv.get<std::vector<double> >(sdTag1, "GroupingZOut"));
1232  php.slopeMin_ = fv.get<std::vector<double> >(sdTag4, "Slope");
1233  if (php.slopeMin_.empty())
1234  php.slopeMin_.emplace_back(0);
1235 
1236  // Wafer size
1237  const auto& dummy = fv.get<std::vector<double> >(sdTag2, "WaferSize");
1239 
1240  // Cell size
1241  php.cellSize_ = fv.get<std::vector<double> >(sdTag3, "CellSize");
1243 
1244  // Layer Offset
1245  const auto& dummy2 = fv.get<std::vector<double> >(sdTag1, "LayerOffset");
1246  if (!dummy2.empty()) {
1247  php.layerOffset_ = dummy2[0];
1248  } else {
1249  php.layerOffset_ = 0;
1250  }
1251 
1252  loadSpecParsHexagon(php);
1253 }
1254 
1256 #ifdef EDM_ML_DEBUG
1257  edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: wafer radius ranges"
1258  << " for cell grouping " << php.boundR_[0] << ":" << php.boundR_[1] << ":"
1259  << php.boundR_[2] << ":" << php.boundR_[3];
1260  edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: Minimum/maximum R " << php.rLimit_[0] << ":" << php.rLimit_[1];
1261  edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: LevelTop " << php.levelT_[0];
1262  edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: minimum slope " << php.slopeMin_[0] << " and layer groupings "
1263  << "for the 3 ranges:";
1264  for (unsigned int k = 0; k < php.layerGroup_.size(); ++k)
1265  edm::LogVerbatim("HGCalGeom") << "[" << k << "] " << php.layerGroup_[k] << ":" << php.layerGroupM_[k] << ":"
1266  << php.layerGroupO_[k];
1267  edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: Wafer Size: " << waferSize_;
1268  edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: " << php.cellSize_.size() << " cells of sizes:";
1269  for (unsigned int k = 0; k < php.cellSize_.size(); ++k)
1270  edm::LogVerbatim("HGCalGeom") << " [" << k << "] " << php.cellSize_[k];
1271  edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: First Layer " << php.firstLayer_ << " and layer offset "
1272  << php.layerOffset_;
1273 #endif
1274 }
1275 
1278  php.cellThickness_ = getDDDArray("CellThickness", sv, 3);
1282  php.waferThickness_ = getDDDArray("WaferThickness", sv, 3);
1284  } else {
1285  for (unsigned int k = 0; k < php.cellThickness_.size(); ++k)
1286  php.waferThickness_.emplace_back(php.waferThick_);
1287  }
1288 
1289  php.radius100to200_ = getDDDArray("Radius100to200", sv, 5);
1290  php.radius200to300_ = getDDDArray("Radius200to300", sv, 5);
1291 
1292  const auto& dummy = getDDDArray("RadiusCuts", sv, 4);
1293  php.choiceType_ = static_cast<int>(dummy[0]);
1294  php.nCornerCut_ = static_cast<int>(dummy[1]);
1295  php.fracAreaMin_ = dummy[2];
1297 
1298  php.radiusMixBoundary_ = fv.vector("RadiusMixBoundary");
1300 
1301  php.slopeMin_ = getDDDArray("SlopeBottom", sv, 0);
1302  php.zFrontMin_ = getDDDArray("ZFrontBottom", sv, 0);
1304  php.rMinFront_ = getDDDArray("RMinFront", sv, 0);
1306 
1307  php.slopeTop_ = getDDDArray("SlopeTop", sv, 0);
1308  php.zFrontTop_ = getDDDArray("ZFrontTop", sv, 0);
1310  php.rMaxFront_ = getDDDArray("RMaxFront", sv, 0);
1312 
1313  php.zRanges_ = fv.vector("ZRanges");
1315 
1316  const auto& dummy2 = getDDDArray("LayerOffset", sv, 1);
1317  php.layerOffset_ = dummy2[0];
1318  php.layerCenter_ = dbl_to_int(fv.vector("LayerCenter"));
1319 
1321  const auto& dummy3 = fv.vector("CalibCellRadius");
1323  php.calibCellFullHD_ = dbl_to_int(fv.vector("CalibCellFullHD"));
1324  php.calibCellPartHD_ = dbl_to_int(fv.vector("CalibCellPartHD"));
1326  php.calibCellFullLD_ = dbl_to_int(fv.vector("CalibCellFullLD"));
1327  php.calibCellPartLD_ = dbl_to_int(fv.vector("CalibCellPartLD"));
1328  } else {
1329  php.calibCellRHD_ = php.calibCellRLD_ = 0;
1330  php.calibCellFullHD_.clear();
1331  php.calibCellPartHD_.clear();
1332  php.calibCellFullLD_.clear();
1333  php.calibCellPartLD_.clear();
1334  }
1335 
1336  loadSpecParsHexagon8(php);
1337 
1338  // Read in parameters from Philip's file
1339  if (php.waferMaskMode_ > 1) {
1340  std::vector<int> layerType, waferIndex, waferProperties;
1341  std::vector<double> cassetteShift;
1342  if (php.waferMaskMode_ == siliconFileEE) {
1343  waferIndex = dbl_to_int(fv.vector("WaferIndexEE"));
1344  waferProperties = dbl_to_int(fv.vector("WaferPropertiesEE"));
1345  } else if (php.waferMaskMode_ == siliconCassetteEE) {
1346  waferIndex = dbl_to_int(fv.vector("WaferIndexEE"));
1347  waferProperties = dbl_to_int(fv.vector("WaferPropertiesEE"));
1348  cassetteShift = fv.vector("CassetteShiftEE");
1349  } else if (php.waferMaskMode_ == siliconFileHE) {
1350  waferIndex = dbl_to_int(fv.vector("WaferIndexHE"));
1351  waferProperties = dbl_to_int(fv.vector("WaferPropertiesHE"));
1352  } else if (php.waferMaskMode_ == siliconCassetteHE) {
1353  waferIndex = dbl_to_int(fv.vector("WaferIndexHE"));
1354  waferProperties = dbl_to_int(fv.vector("WaferPropertiesHE"));
1355  cassetteShift = fv.vector("CassetteShiftHE");
1356  }
1359  if ((php.waferMaskMode_ == siliconFileEE) || (php.waferMaskMode_ == siliconCassetteEE)) {
1360  layerType = dbl_to_int(fv.vector("LayerTypesEE"));
1361  } else if ((php.waferMaskMode_ == siliconFileHE) || (php.waferMaskMode_ == siliconCassetteHE)) {
1362  layerType = dbl_to_int(fv.vector("LayerTypesHE"));
1363  }
1364  }
1365 
1366  php.cassetteShift_ = cassetteShift;
1368  loadSpecParsHexagon8(php, layerType, waferIndex, waferProperties);
1369  }
1370 }
1371 
1373  const cms::DDVectorsMap& vmap,
1374  HGCalParameters& php,
1375  const std::string& sdTag1) {
1376  php.cellThickness_ = fv.get<std::vector<double> >(sdTag1, "CellThickness");
1380  php.waferThickness_ = fv.get<std::vector<double> >(sdTag1, "WaferThickness");
1382  } else {
1383  for (unsigned int k = 0; k < php.cellThickness_.size(); ++k)
1384  php.waferThickness_.emplace_back(php.waferThick_);
1385  }
1386 
1387  php.radius100to200_ = fv.get<std::vector<double> >(sdTag1, "Radius100to200");
1388  php.radius200to300_ = fv.get<std::vector<double> >(sdTag1, "Radius200to300");
1389 
1390  const auto& dummy = fv.get<std::vector<double> >(sdTag1, "RadiusCuts");
1391  if (dummy.size() > 3) {
1392  php.choiceType_ = static_cast<int>(dummy[0]);
1393  php.nCornerCut_ = static_cast<int>(dummy[1]);
1394  php.fracAreaMin_ = dummy[2];
1396  } else {
1397  php.choiceType_ = php.nCornerCut_ = php.fracAreaMin_ = php.zMinForRad_ = 0;
1398  }
1399 
1400  php.slopeMin_ = fv.get<std::vector<double> >(sdTag1, "SlopeBottom");
1401  php.zFrontMin_ = fv.get<std::vector<double> >(sdTag1, "ZFrontBottom");
1403  php.rMinFront_ = fv.get<std::vector<double> >(sdTag1, "RMinFront");
1405 
1406  php.slopeTop_ = fv.get<std::vector<double> >(sdTag1, "SlopeTop");
1407  php.zFrontTop_ = fv.get<std::vector<double> >(sdTag1, "ZFrontTop");
1409  php.rMaxFront_ = fv.get<std::vector<double> >(sdTag1, "RMaxFront");
1411  unsigned int kmax = (php.zFrontTop_.size() - php.slopeTop_.size());
1412  for (unsigned int k = 0; k < kmax; ++k)
1413  php.slopeTop_.emplace_back(0);
1414 
1415  const auto& dummy2 = fv.get<std::vector<double> >(sdTag1, "LayerOffset");
1416  if (!dummy2.empty()) {
1417  php.layerOffset_ = dummy2[0];
1418  } else {
1419  php.layerOffset_ = 0;
1420  }
1421 
1422  php.calibCellRHD_ = php.calibCellRLD_ = 0;
1423  php.calibCellFullHD_.clear();
1424  php.calibCellPartHD_.clear();
1425  php.calibCellFullLD_.clear();
1426  php.calibCellPartLD_.clear();
1427  for (auto const& it : vmap) {
1428  if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "RadiusMixBoundary")) {
1429  for (const auto& i : it.second)
1431  } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "ZRanges")) {
1432  for (const auto& i : it.second)
1433  php.zRanges_.emplace_back(HGCalParameters::k_ScaleFromDD4hep * i);
1434  } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "LayerCenter")) {
1435  for (const auto& i : it.second)
1436  php.layerCenter_.emplace_back(std::round(i));
1437  } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "CalibCellRadius")) {
1440  } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "CalibCellFullHD")) {
1441  for (const auto& i : it.second)
1442  php.calibCellFullHD_.emplace_back(std::round(i));
1443  } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "CalibCellPartHD")) {
1444  for (const auto& i : it.second)
1445  php.calibCellPartHD_.emplace_back(std::round(i));
1446  } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "CalibCellFullLD")) {
1447  for (const auto& i : it.second)
1448  php.calibCellFullLD_.emplace_back(std::round(i));
1449  } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "CalibCellPartLD")) {
1450  for (const auto& i : it.second)
1451  php.calibCellPartLD_.emplace_back(std::round(i));
1452  }
1453  }
1454 
1455  loadSpecParsHexagon8(php);
1456 
1457  // Read in parameters from Philip's file
1458  if (php.waferMaskMode_ > 1) {
1459  std::vector<int> layerType, waferIndex, waferProperties;
1460  std::vector<double> cassetteShift;
1461  if ((php.waferMaskMode_ == siliconFileEE) || (php.waferMaskMode_ == siliconCassetteEE)) {
1462  for (auto const& it : vmap) {
1463  if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "WaferIndexEE")) {
1464  for (const auto& i : it.second)
1465  waferIndex.emplace_back(std::round(i));
1466  } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "WaferPropertiesEE")) {
1467  for (const auto& i : it.second)
1468  waferProperties.emplace_back(std::round(i));
1469  }
1470  }
1471  if (php.waferMaskMode_ == siliconCassetteEE) {
1472  for (auto const& it : vmap) {
1473  if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "CassetteShiftEE")) {
1474  for (const auto& i : it.second)
1475  cassetteShift.emplace_back(i);
1476  }
1477  }
1478  }
1479  } else if ((php.waferMaskMode_ == siliconFileHE) || (php.waferMaskMode_ == siliconCassetteHE)) {
1480  for (auto const& it : vmap) {
1481  if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "WaferIndexHE")) {
1482  for (const auto& i : it.second)
1483  waferIndex.emplace_back(std::round(i));
1484  } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "WaferPropertiesHE")) {
1485  for (const auto& i : it.second)
1486  waferProperties.emplace_back(std::round(i));
1487  }
1488  }
1489  if (php.waferMaskMode_ == siliconCassetteHE) {
1490  for (auto const& it : vmap) {
1491  if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "CassetteShiftHE")) {
1492  for (const auto& i : it.second)
1493  cassetteShift.emplace_back(i);
1494  }
1495  }
1496  }
1497  }
1500  if ((php.waferMaskMode_ == siliconFileEE) || (php.waferMaskMode_ == siliconCassetteEE)) {
1501  for (auto const& it : vmap) {
1502  if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "LayerTypesEE")) {
1503  for (const auto& i : it.second)
1504  layerType.emplace_back(std::round(i));
1505  }
1506  }
1507  } else if ((php.waferMaskMode_ == siliconFileHE) || (php.waferMaskMode_ == siliconCassetteHE)) {
1508  for (auto const& it : vmap) {
1509  if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "LayerTypesHE")) {
1510  for (const auto& i : it.second)
1511  layerType.emplace_back(std::round(i));
1512  }
1513  }
1514  }
1515  }
1516 
1517  php.cassetteShift_ = cassetteShift;
1519  loadSpecParsHexagon8(php, layerType, waferIndex, waferProperties);
1520  }
1521 }
1522 
1524 #ifdef EDM_ML_DEBUG
1525  for (unsigned int k = 0; k < php.waferThickness_.size(); ++k)
1526  edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: wafer[" << k << "] Thickness " << php.waferThickness_[k];
1527  for (unsigned int k = 0; k < php.cellThickness_.size(); ++k)
1528  edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: cell[" << k << "] Thickness " << php.cellThickness_[k];
1529  edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: Polynomial "
1530  << "parameters for 120 to 200 micron "
1531  << "transition with " << php.radius100to200_.size() << " elements";
1532  for (unsigned int k = 0; k < php.radius100to200_.size(); ++k)
1533  edm::LogVerbatim("HGCalGeom") << "Element [" << k << "] " << php.radius100to200_[k];
1534  edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: Polynomial "
1535  << "parameters for 200 to 300 micron "
1536  << "transition with " << php.radius200to300_.size() << " elements";
1537  for (unsigned int k = 0; k < php.radius200to300_.size(); ++k)
1538  edm::LogVerbatim("HGCalGeom") << "Element [" << k << "] " << php.radius200to300_[k];
1539  edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: Parameters for the"
1540  << " transition " << php.choiceType_ << ":" << php.nCornerCut_ << ":"
1541  << php.fracAreaMin_ << ":" << php.zMinForRad_;
1542  for (unsigned int k = 0; k < php.radiusMixBoundary_.size(); ++k)
1543  edm::LogVerbatim("HGCalGeom") << "HGCalParameters: Mix[" << k << "] R = " << php.radiusMixBoundary_[k];
1544  for (unsigned int k = 0; k < php.zFrontMin_.size(); ++k)
1545  edm::LogVerbatim("HGCalGeom") << "HGCalParameters: Boundary[" << k << "] Bottom Z = " << php.zFrontMin_[k]
1546  << " Slope = " << php.slopeMin_[k] << " rMax = " << php.rMinFront_[k];
1547  for (unsigned int k = 0; k < php.zFrontTop_.size(); ++k)
1548  edm::LogVerbatim("HGCalGeom") << "HGCalParameters: Boundary[" << k << "] Top Z = " << php.zFrontTop_[k]
1549  << " Slope = " << php.slopeTop_[k] << " rMax = " << php.rMaxFront_[k];
1550  std::ostringstream st1;
1551  for (unsigned int k = 0; k < php.zRanges_.size(); ++k)
1552  st1 << ":" << php.zRanges_[k];
1553  edm::LogVerbatim("HGCalGeom") << "HGCalParameters: Z-Boundary[" << php.zRanges_.size() << "] " << st1.str();
1554  edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: LayerOffset " << php.layerOffset_ << " in array of size "
1555  << php.layerCenter_.size();
1556  for (unsigned int k = 0; k < php.layerCenter_.size(); ++k)
1557  edm::LogVerbatim("HGCalGeom") << "[" << k << "] " << php.layerCenter_[k];
1558  edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: Calibration cells in HD having radius " << php.calibCellRHD_
1559  << " in arrays of size " << php.calibCellFullHD_.size() << ":"
1560  << php.calibCellPartHD_.size();
1561  for (unsigned int k = 0; k < php.calibCellFullHD_.size(); ++k)
1562  edm::LogVerbatim("HGCalGeom") << "[" << k << "] " << php.calibCellFullHD_[k] << ":" << php.calibCellPartHD_[k];
1563  edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: Calibration cells in LD having radius " << php.calibCellRLD_
1564  << " in arrays of size " << php.calibCellFullLD_.size() << ":"
1565  << php.calibCellPartLD_.size();
1566  for (unsigned int k = 0; k < php.calibCellFullLD_.size(); ++k)
1567  edm::LogVerbatim("HGCalGeom") << "[" << k << "] " << php.calibCellFullLD_[k] << ":" << php.calibCellPartLD_[k];
1568 #endif
1569 }
1570 
1572  const std::vector<int>& layerType,
1573  const std::vector<int>& waferIndex,
1574  const std::vector<int>& waferProperties) {
1575  // Store parameters from Philip's file
1576  for (unsigned int k = 0; k < layerType.size(); ++k) {
1577  php.layerType_.emplace_back(HGCalTypes::layerType(layerType[k]));
1578 #ifdef EDM_ML_DEBUG
1579  edm::LogVerbatim("HGCalGeom") << "Layer[" << k << "] Type " << layerType[k] << ":" << php.layerType_.back();
1580 #endif
1581  }
1582  for (unsigned int k = 0; k < php.layerType_.size(); ++k) {
1583  double cth = (php.layerType_[k] == HGCalTypes::WaferCenterR) ? cos(php.layerRotation_) : 1.0;
1584  double sth = (php.layerType_[k] == HGCalTypes::WaferCenterR) ? sin(php.layerRotation_) : 0.0;
1585  php.layerRotV_.emplace_back(std::make_pair(cth, sth));
1586 #ifdef EDM_ML_DEBUG
1587  edm::LogVerbatim("HGCalGeom") << "Layer[" << k << "] Type " << php.layerType_[k] << " cos|sin(Theta) "
1588  << php.layerRotV_.back().first << ":" << php.layerRotV_.back().second;
1589 #endif
1590  }
1591  for (unsigned int k = 0; k < waferIndex.size(); ++k) {
1592  int partial = HGCalProperty::waferPartial(waferProperties[k]);
1593  int orient =
1595  ? HGCalProperty::waferOrient(waferProperties[k])
1596  : HGCalWaferMask::getRotation(php.waferZSide_, partial, HGCalProperty::waferOrient(waferProperties[k]));
1598  partial,
1599  orient,
1600  HGCalProperty::waferCassette(waferProperties[k]));
1601 #ifdef EDM_ML_DEBUG
1602  edm::LogVerbatim("HGCalGeom") << "[" << k << ":" << waferIndex[k] << ":"
1605  << HGCalWaferIndex::waferV(waferIndex[k]) << "] Thickness type "
1606  << HGCalProperty::waferThick(waferProperties[k]) << " Partial type " << partial
1607  << " Orientation " << HGCalProperty::waferOrient(waferProperties[k]) << ":" << orient;
1608 #endif
1609  }
1610 }
1611 
1614  php.radiusMixBoundary_ = fv.vector("RadiusMixBoundary");
1616 
1617  php.nPhiBinBH_ = dbl_to_int(getDDDArray("NPhiBinBH", sv, 0));
1618  php.layerFrontBH_ = dbl_to_int(getDDDArray("LayerFrontBH", sv, 0));
1619  php.rMinLayerBH_ = getDDDArray("RMinLayerBH", sv, 0);
1621  assert(php.nPhiBinBH_.size() > 1);
1622  php.nCellsFine_ = php.nPhiBinBH_[0];
1623  php.nCellsCoarse_ = php.nPhiBinBH_[1];
1624  assert(0 != php.nCellsFine_);
1625  assert(0 != php.nCellsCoarse_);
1626  php.cellSize_.emplace_back(2.0 * M_PI / php.nCellsFine_);
1627  php.cellSize_.emplace_back(2.0 * M_PI / php.nCellsCoarse_);
1628 
1629  php.slopeMin_ = getDDDArray("SlopeBottom", sv, 0);
1630  php.zFrontMin_ = getDDDArray("ZFrontBottom", sv, 0);
1632  php.rMinFront_ = getDDDArray("RMinFront", sv, 0);
1634 
1635  php.slopeTop_ = getDDDArray("SlopeTop", sv, 0);
1636  php.zFrontTop_ = getDDDArray("ZFrontTop", sv, 0);
1638  php.rMaxFront_ = getDDDArray("RMaxFront", sv, 0);
1640 
1641  php.zRanges_ = fv.vector("ZRanges");
1643 
1644  // Offsets
1645  const auto& dummy2 = getDDDArray("LayerOffset", sv, 1);
1646  php.layerOffset_ = dummy2[0];
1647  php.layerCenter_ = dbl_to_int(fv.vector("LayerCenter"));
1648 
1649  loadSpecParsTrapezoid(php);
1650 #ifdef EDM_ML_DEBUG
1651  edm::LogVerbatim("HGCalGeom") << "WaferMaskMode " << php.waferMaskMode_ << " Compare " << scintillatorFile << ":"
1653 #endif
1654  // tile parameters from Katja's file
1657  std::vector<int> tileIndx, tileProperty;
1658  std::vector<int> tileHEX1, tileHEX2, tileHEX3, tileHEX4, tileHEX5, tileHEX6;
1659  std::vector<double> tileRMin, tileRMax, tileRMinFine, tileRMaxFine;
1660  std::vector<int> tileRingMin, tileRingMax, tileRingMinFine, tileRingMaxFine;
1661  std::vector<double> cassetteShift;
1662  php.nPhiLayer_ = dbl_to_int(fv.vector("NPhiLayer"));
1663  tileIndx = dbl_to_int(fv.vector("TileIndex"));
1664  tileProperty = dbl_to_int(fv.vector("TileProperty"));
1665  tileHEX1 = dbl_to_int(fv.vector("TileHEX1"));
1666  tileHEX2 = dbl_to_int(fv.vector("TileHEX2"));
1667  tileHEX3 = dbl_to_int(fv.vector("TileHEX3"));
1668  tileHEX4 = dbl_to_int(fv.vector("TileHEX4"));
1669  tileRMin = fv.vector("TileRMin");
1670  tileRMax = fv.vector("TileRMax");
1673  tileRingMin = dbl_to_int(fv.vector("TileRingMin"));
1674  tileRingMax = dbl_to_int(fv.vector("TileRingMax"));
1675  if (php.waferMaskMode_ == scintillatorFineCell) {
1676  tileHEX5 = dbl_to_int(fv.vector("TileHEX5"));
1677  tileHEX6 = dbl_to_int(fv.vector("TileHEX6"));
1678  tileRMinFine = fv.vector("TileRMin6");
1679  tileRMaxFine = fv.vector("TileRMax6");
1682  tileRingMinFine = dbl_to_int(fv.vector("TileRingMin6"));
1683  tileRingMaxFine = dbl_to_int(fv.vector("TileRingMax6"));
1684  php.nphiCassette_ = php.nCellsCoarse_ / php.cassettes_;
1685  php.nphiFineCassette_ = php.nCellsFine_ / php.cassettes_;
1686  std::vector<double> rectract = fv.vector("ScintRetract");
1688  double dphi = M_PI / php.cassettes_;
1689  for (int k = 0; k < php.cassettes_; ++k) {
1690  double phi = (2 * k + 1) * dphi;
1691  cassetteShift.emplace_back(rectract[k] * cos(phi));
1692  cassetteShift.emplace_back(rectract[k] * sin(phi));
1693  }
1694  } else if (php.waferMaskMode_ == scintillatorCassette) {
1695  if (php.cassettes_ > 0)
1696  php.nphiCassette_ = php.nCellsCoarse_ / php.cassettes_;
1697  cassetteShift = fv.vector("CassetteShiftHE");
1698  }
1699 
1700  rescale(cassetteShift, HGCalParameters::k_ScaleFromDDD);
1702  php.cassetteShiftTile_ = cassetteShift;
1703  else
1704  php.cassetteShift_ = cassetteShift;
1706  tileIndx,
1707  tileProperty,
1708  tileHEX1,
1709  tileHEX2,
1710  tileHEX3,
1711  tileHEX4,
1712  tileHEX5,
1713  tileHEX6,
1714  tileRMin,
1715  tileRMax,
1716  tileRMinFine,
1717  tileRMaxFine,
1718  tileRingMin,
1719  tileRingMax,
1720  tileRingMinFine,
1721  tileRingMaxFine);
1722  }
1723 }
1724 
1726  const cms::DDVectorsMap& vmap,
1727  HGCalParameters& php,
1728  const std::string& sdTag1) {
1729  for (auto const& it : vmap) {
1730  if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "RadiusMixBoundary")) {
1731  for (const auto& i : it.second)
1733  } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "ZRanges")) {
1734  for (const auto& i : it.second)
1735  php.zRanges_.emplace_back(HGCalParameters::k_ScaleFromDD4hep * i);
1736  } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "LayerCenter")) {
1737  for (const auto& i : it.second)
1738  php.layerCenter_.emplace_back(std::round(i));
1739  }
1740  }
1741 
1742  php.nPhiBinBH_ = dbl_to_int(fv.get<std::vector<double> >(sdTag1, "NPhiBinBH"));
1743  php.layerFrontBH_ = dbl_to_int(fv.get<std::vector<double> >(sdTag1, "LayerFrontBH"));
1744  php.rMinLayerBH_ = fv.get<std::vector<double> >(sdTag1, "RMinLayerBH");
1746  assert(php.nPhiBinBH_.size() > 1);
1747  php.nCellsFine_ = php.nPhiBinBH_[0];
1748  php.nCellsCoarse_ = php.nPhiBinBH_[1];
1749  assert(0 != php.nCellsFine_);
1750  assert(0 != php.nCellsCoarse_);
1751  php.cellSize_.emplace_back(2.0 * M_PI / php.nCellsFine_);
1752  php.cellSize_.emplace_back(2.0 * M_PI / php.nCellsCoarse_);
1753 
1754  php.slopeMin_ = fv.get<std::vector<double> >(sdTag1, "SlopeBottom");
1755  php.zFrontMin_ = fv.get<std::vector<double> >(sdTag1, "ZFrontBottom");
1757  php.rMinFront_ = fv.get<std::vector<double> >(sdTag1, "RMinFront");
1759 
1760  php.slopeTop_ = fv.get<std::vector<double> >(sdTag1, "SlopeTop");
1761  php.zFrontTop_ = fv.get<std::vector<double> >(sdTag1, "ZFrontTop");
1763  php.rMaxFront_ = fv.get<std::vector<double> >(sdTag1, "RMaxFront");
1765  unsigned int kmax = (php.zFrontTop_.size() - php.slopeTop_.size());
1766  for (unsigned int k = 0; k < kmax; ++k)
1767  php.slopeTop_.emplace_back(0);
1768 
1769  const auto& dummy2 = fv.get<std::vector<double> >(sdTag1, "LayerOffset");
1770  php.layerOffset_ = dummy2[0];
1771 
1772  loadSpecParsTrapezoid(php);
1773 
1774  // tile parameters from Katja's file
1777  std::vector<int> tileIndx, tileProperty;
1778  std::vector<int> tileHEX1, tileHEX2, tileHEX3, tileHEX4, tileHEX5, tileHEX6;
1779  std::vector<double> tileRMin, tileRMax, tileRMinFine, tileRMaxFine;
1780  std::vector<int> tileRingMin, tileRingMax, tileRingMinFine, tileRingMaxFine;
1781  std::vector<double> cassetteShift;
1782  for (auto const& it : vmap) {
1783  if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileIndex")) {
1784  for (const auto& i : it.second)
1785  tileIndx.emplace_back(std::round(i));
1786  } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileProperty")) {
1787  for (const auto& i : it.second)
1788  tileProperty.emplace_back(std::round(i));
1789  } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "NPhiLayer")) {
1790  for (const auto& i : it.second)
1791  php.nPhiLayer_.emplace_back(std::round(i));
1792  } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileHEX1")) {
1793  for (const auto& i : it.second)
1794  tileHEX1.emplace_back(std::round(i));
1795  } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileHEX2")) {
1796  for (const auto& i : it.second)
1797  tileHEX2.emplace_back(std::round(i));
1798  } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileHEX3")) {
1799  for (const auto& i : it.second)
1800  tileHEX3.emplace_back(std::round(i));
1801  } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileHEX4")) {
1802  for (const auto& i : it.second)
1803  tileHEX4.emplace_back(std::round(i));
1804  } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileRMin")) {
1805  for (const auto& i : it.second)
1806  tileRMin.emplace_back(HGCalParameters::k_ScaleFromDD4hep * i);
1807  } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileRMax")) {
1808  for (const auto& i : it.second)
1809  tileRMax.emplace_back(HGCalParameters::k_ScaleFromDD4hep * i);
1810  } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileRingMin")) {
1811  for (const auto& i : it.second)
1812  tileRingMin.emplace_back(std::round(i));
1813  } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileRingMax")) {
1814  for (const auto& i : it.second)
1815  tileRingMax.emplace_back(std::round(i));
1816  }
1817  }
1818  if (php.waferMaskMode_ == scintillatorFineCell) {
1819  php.nphiCassette_ = php.nCellsCoarse_ / php.cassettes_;
1820  php.nphiFineCassette_ = php.nCellsFine_ / php.cassettes_;
1821  std::vector<double> rectract;
1822  for (auto const& it : vmap) {
1823  if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileHEX5")) {
1824  for (const auto& i : it.second)
1825  tileHEX5.emplace_back(std::round(i));
1826  } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileHEX6")) {
1827  for (const auto& i : it.second)
1828  tileHEX6.emplace_back(std::round(i));
1829  } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileRMin6")) {
1830  for (const auto& i : it.second)
1831  tileRMinFine.emplace_back(i);
1833  } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileRMax6")) {
1834  for (const auto& i : it.second)
1835  tileRMaxFine.emplace_back(i);
1837  } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileRingMin6")) {
1838  for (const auto& i : it.second)
1839  tileRingMinFine.emplace_back(std::round(i));
1840  } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileRingMax6")) {
1841  for (const auto& i : it.second)
1842  tileRingMaxFine.emplace_back(std::round(i));
1843  } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "ScintRetract")) {
1844  for (const auto& i : it.second)
1845  rectract.emplace_back(i);
1846  double dphi = M_PI / php.cassettes_;
1847  for (int k = 0; k < php.cassettes_; ++k) {
1848  double phi = (2 * k + 1) * dphi;
1849  cassetteShift.emplace_back(rectract[k] * cos(phi));
1850  cassetteShift.emplace_back(rectract[k] * sin(phi));
1851  }
1852  }
1853  }
1854  } else if (php.waferMaskMode_ == scintillatorCassette) {
1855  for (auto const& it : vmap) {
1856  if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "CassetteShiftHE")) {
1857  for (const auto& i : it.second)
1858  cassetteShift.emplace_back(i);
1859  }
1860  }
1861  }
1862 
1865  php.cassetteShiftTile_ = cassetteShift;
1866  else
1867  php.cassetteShift_ = cassetteShift;
1869  tileIndx,
1870  tileProperty,
1871  tileHEX1,
1872  tileHEX2,
1873  tileHEX3,
1874  tileHEX4,
1875  tileHEX5,
1876  tileHEX6,
1877  tileRMin,
1878  tileRMax,
1879  tileRMinFine,
1880  tileRMaxFine,
1881  tileRingMin,
1882  tileRingMax,
1883  tileRingMinFine,
1884  tileRingMaxFine);
1885  }
1886 }
1887 
1889 #ifdef EDM_ML_DEBUG
1890  edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters:nCells " << php.nCellsFine_ << ":" << php.nCellsCoarse_
1891  << " cellSize: " << php.cellSize_[0] << ":" << php.cellSize_[1];
1892  for (unsigned int k = 0; k < php.layerFrontBH_.size(); ++k)
1893  edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: Type[" << k << "] Front Layer = " << php.layerFrontBH_[k]
1894  << " rMin = " << php.rMinLayerBH_[k];
1895  for (unsigned int k = 0; k < php.radiusMixBoundary_.size(); ++k) {
1896  edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: Mix[" << k << "] R = " << php.radiusMixBoundary_[k]
1897  << " Nphi = " << php.scintCells(k + php.firstLayer_)
1898  << " dPhi = " << php.scintCellSize(k + php.firstLayer_);
1899  }
1900 
1901  for (unsigned int k = 0; k < php.zFrontMin_.size(); ++k)
1902  edm::LogVerbatim("HGCalGeom") << "HGCalParameters: Boundary[" << k << "] Bottom Z = " << php.zFrontMin_[k]
1903  << " Slope = " << php.slopeMin_[k] << " rMax = " << php.rMinFront_[k];
1904 
1905  for (unsigned int k = 0; k < php.zFrontTop_.size(); ++k)
1906  edm::LogVerbatim("HGCalGeom") << "HGCalParameters: Boundary[" << k << "] Top Z = " << php.zFrontTop_[k]
1907  << " Slope = " << php.slopeTop_[k] << " rMax = " << php.rMaxFront_[k];
1908 
1909  std::ostringstream st1;
1910  for (unsigned int k = 0; k < php.zRanges_.size(); ++k)
1911  st1 << ":" << php.zRanges_[k];
1912  edm::LogVerbatim("HGCalGeom") << "HGCalParameters: Z-Boundary[" << php.zRanges_.size() << "] " << st1.str();
1913 
1914  edm::LogVerbatim("HGCalGeom") << "HGCalParameters: LayerOffset " << php.layerOffset_ << " in array of size "
1915  << php.layerCenter_.size();
1916  for (unsigned int k = 0; k < php.layerCenter_.size(); ++k)
1917  edm::LogVerbatim("HGCalGeom") << "[" << k << "] " << php.layerCenter_[k];
1918 #endif
1919 }
1920 
1922  const std::vector<int>& tileIndx,
1923  const std::vector<int>& tileProperty,
1924  const std::vector<int>& tileHEX1,
1925  const std::vector<int>& tileHEX2,
1926  const std::vector<int>& tileHEX3,
1927  const std::vector<int>& tileHEX4,
1928  const std::vector<int>& tileHEX5,
1929  const std::vector<int>& tileHEX6,
1930  const std::vector<double>& tileRMin,
1931  const std::vector<double>& tileRMax,
1932  const std::vector<double>& tileRMinFine,
1933  const std::vector<double>& tileRMaxFine,
1934  const std::vector<int>& tileRingMin,
1935  const std::vector<int>& tileRingMax,
1936  const std::vector<int>& tileRingMinFine,
1937  const std::vector<int>& tileRingMaxFine) {
1938  // tile parameters from Katja's file
1939  for (unsigned int k = 0; k < tileIndx.size(); ++k) {
1940  int hex5 = (k < tileHEX5.size()) ? tileHEX5[k] : 0;
1941  int hex6 = (k < tileHEX6.size()) ? tileHEX6[k] : 0;
1944  tileHEX1[k],
1945  tileHEX2[k],
1946  tileHEX3[k],
1947  tileHEX4[k],
1948  hex5,
1949  hex6);
1950 #ifdef EDM_ML_DEBUG
1951  edm::LogVerbatim("HGCalGeom") << "Tile[" << k << ":" << tileIndx[k] << "] "
1952  << " Type " << HGCalProperty::tileType(tileProperty[k]) << " SiPM "
1953  << HGCalProperty::tileSiPM(tileProperty[k]) << " HEX " << std::hex << tileHEX1[k]
1954  << ":" << tileHEX2[k] << ":" << tileHEX3[k] << ":" << tileHEX4[k] << ":" << hex5
1955  << ":" << hex6 << std::dec;
1956 #endif
1957  }
1958 
1959  for (unsigned int k = 0; k < tileRMinFine.size(); ++k) {
1960  php.tileRingFineR_.emplace_back(tileRMinFine[k], tileRMaxFine[k]);
1961 #ifdef EDM_ML_DEBUG
1962  edm::LogVerbatim("HGCalGeom") << "TileRingFineR[" << k << "] " << tileRMinFine[k] << ":" << tileRMaxFine[k];
1963 #endif
1964  }
1965  for (unsigned int k = 0; k < tileRMin.size(); ++k) {
1966  php.tileRingR_.emplace_back(tileRMin[k], tileRMax[k]);
1967 #ifdef EDM_ML_DEBUG
1968  edm::LogVerbatim("HGCalGeom") << "TileRingR[" << k << "] " << tileRMin[k] << ":" << tileRMax[k];
1969 #endif
1970  }
1971 
1972  for (unsigned int k = 0; k < tileRingMinFine.size(); ++k) {
1973  php.tileRingFineRange_.emplace_back(tileRingMinFine[k], tileRingMaxFine[k]);
1974 #ifdef EDM_ML_DEBUG
1975  edm::LogVerbatim("HGCalGeom") << "TileRingFineRange[" << k << "] " << tileRingMinFine[k] << ":"
1976  << tileRingMaxFine[k];
1977 #endif
1978  }
1979  for (unsigned k = 0; k < tileRingMin.size(); ++k) {
1980  php.tileRingRange_.emplace_back(tileRingMin[k], tileRingMax[k]);
1981 #ifdef EDM_ML_DEBUG
1982  edm::LogVerbatim("HGCalGeom") << "TileRingRange[" << k << "] " << tileRingMin[k] << ":" << tileRingMax[k];
1983 #endif
1984  }
1985 #ifdef EDM_ML_DEBUG
1986  std::ostringstream st1, st2;
1987  st1 << "Scintillator CassetteShift for " << php.cassetteShiftTile_.size() << " tiles:";
1988  for (unsigned int k = 0; k < php.cassetteShiftTile_.size(); ++k)
1989  st1 << " " << php.cassetteShiftTile_[k];
1990  edm::LogVerbatim("HGCalGeom") << st1.str();
1991  st2 << "NPhi for scntillator " << php.nPhiLayer_.size() << " layers:";
1992  for (unsigned int k = 0; k < php.nPhiLayer_.size(); ++k)
1993  st2 << " " << php.nPhiLayer_[k];
1994  edm::LogVerbatim("HGCalGeom") << st2.str();
1995 #endif
1996 }
1997 
2000  double rin(php.rLimit_[0]), rout(php.rLimit_[1]), rMaxFine(php.boundR_[1]);
2001 #ifdef EDM_ML_DEBUG
2002  edm::LogVerbatim("HGCalGeom") << "Input waferWidth " << waferW << ":" << rmin << " R Limits: " << rin << ":" << rout
2003  << " Fine " << rMaxFine;
2004 #endif
2005  // Clear the vectors
2006  php.waferCopy_.clear();
2007  php.waferTypeL_.clear();
2008  php.waferTypeT_.clear();
2009  php.waferPosX_.clear();
2010  php.waferPosY_.clear();
2011  double dx = 0.5 * waferW;
2012  double dy = 3.0 * dx * tan(30._deg);
2013  double rr = 2.0 * dx * tan(30._deg);
2014  int ncol = static_cast<int>(2.0 * rout / waferW) + 1;
2015  int nrow = static_cast<int>(rout / (waferW * tan(30._deg))) + 1;
2016  int ns2 = (2 * ncol + 1) * (2 * nrow + 1) * php.layer_.size();
2017  int incm(0), inrm(0);
2018  HGCalParameters::layer_map copiesInLayers(php.layer_.size() + 1);
2019  HGCalParameters::waferT_map waferTypes(ns2 + 1);
2020 #ifdef EDM_ML_DEBUG
2021  int kount(0), ntot(0);
2022  edm::LogVerbatim("HGCalGeom") << "Row " << nrow << " Column " << ncol;
2023 #endif
2024  for (int nr = -nrow; nr <= nrow; ++nr) {
2025  int inr = (nr >= 0) ? nr : -nr;
2026  for (int nc = -ncol; nc <= ncol; ++nc) {
2027  int inc = (nc >= 0) ? nc : -nc;
2028  if (inr % 2 == inc % 2) {
2029  double xpos = nc * dx;
2030  double ypos = nr * dy;
2031  std::pair<int, int> corner = HGCalGeomTools::waferCorner(xpos, ypos, dx, rr, rin, rout, true);
2032  double rpos = std::sqrt(xpos * xpos + ypos * ypos);
2033  int typet = (rpos < rMaxFine) ? 1 : 2;
2034  int typel(3);
2035  for (int k = 1; k < 4; ++k) {
2036  if ((rpos + rmin) <= php.boundR_[k]) {
2037  typel = k;
2038  break;
2039  }
2040  }
2041 #ifdef EDM_ML_DEBUG
2042  ++ntot;
2043 #endif
2044  if (corner.first > 0) {
2045  int copy = HGCalTypes::packTypeUV(typel, nc, nr);
2046  if (inc > incm)
2047  incm = inc;
2048  if (inr > inrm)
2049  inrm = inr;
2050 #ifdef EDM_ML_DEBUG
2051  kount++;
2052  edm::LogVerbatim("HGCalGeom") << kount << ":" << ntot << " Copy " << copy << " Type " << typel << ":" << typet
2053  << " Location " << corner.first << " Position " << xpos << ":" << ypos
2054  << " Layers " << php.layer_.size();
2055 #endif
2056  php.waferCopy_.emplace_back(copy);
2057  php.waferTypeL_.emplace_back(typel);
2058  php.waferTypeT_.emplace_back(typet);
2059  php.waferPosX_.emplace_back(xpos);
2060  php.waferPosY_.emplace_back(ypos);
2061  for (unsigned int il = 0; il < php.layer_.size(); ++il) {
2062  std::pair<int, int> corner =
2063  HGCalGeomTools::waferCorner(xpos, ypos, dx, rr, php.rMinLayHex_[il], php.rMaxLayHex_[il], true);
2064  if (corner.first > 0) {
2065  auto cpy = copiesInLayers[php.layer_[il]].find(copy);
2066  if (cpy == copiesInLayers[php.layer_[il]].end())
2067  copiesInLayers[php.layer_[il]][copy] =
2068  ((corner.first == static_cast<int>(HGCalParameters::k_CornerSize)) ? php.waferCopy_.size() : -1);
2069  }
2070  if ((corner.first > 0) && (corner.first < static_cast<int>(HGCalParameters::k_CornerSize))) {
2071  int wl = HGCalWaferIndex::waferIndex(php.layer_[il], copy, 0, true);
2072  waferTypes[wl] = corner;
2073  }
2074  }
2075  }
2076  }
2077  }
2078  }
2079  php.copiesInLayers_ = copiesInLayers;
2080  php.waferTypes_ = waferTypes;
2081  php.nSectors_ = static_cast<int>(php.waferCopy_.size());
2082  php.waferUVMax_ = 0;
2083 #ifdef EDM_ML_DEBUG
2084  edm::LogVerbatim("HGCalGeom") << "HGCalWaferHexagon: # of columns " << incm << " # of rows " << inrm << " and "
2085  << kount << ":" << ntot << " wafers; R " << rin << ":" << rout;
2086  edm::LogVerbatim("HGCalGeom") << "Dump copiesInLayers for " << php.copiesInLayers_.size() << " layers";
2087  for (unsigned int k = 0; k < copiesInLayers.size(); ++k) {
2088  const auto& theModules = copiesInLayers[k];
2089  edm::LogVerbatim("HGCalGeom") << "Layer " << k << ":" << theModules.size();
2090  int k2(0);
2091  for (std::unordered_map<int, int>::const_iterator itr = theModules.begin(); itr != theModules.end(); ++itr, ++k2) {
2092  edm::LogVerbatim("HGCalGeom") << "[" << k2 << "] " << itr->first << ":" << itr->second;
2093  }
2094  }
2095 #endif
2096 }
2097 
2099  double waferW(php.waferSize_);
2100  double waferS(php.sensorSeparation_);
2101  auto wType = std::make_unique<HGCalWaferType>(php.radius100to200_,
2102  php.radius200to300_,
2103  HGCalParameters::k_ScaleToDDD * (waferW + waferS),
2105  php.choiceType_,
2106  php.nCornerCut_,
2107  php.fracAreaMin_);
2108 
2109  double rout(php.rLimit_[1]);
2110 #ifdef EDM_ML_DEBUG
2111  edm::LogVerbatim("HGCalGeom") << "Input waferWidth " << waferW << ":" << waferS << " R Max: " << rout;
2112 #endif
2113  // Clear the vectors
2114  php.waferCopy_.clear();
2115  php.waferTypeL_.clear();
2116  php.waferTypeT_.clear();
2117  php.waferPosX_.clear();
2118  php.waferPosY_.clear();
2119  double r = 0.5 * (waferW + waferS);
2120  double R = 2.0 * r / sqrt3_;
2121  double dy = 0.75 * R;
2122  double r1 = 0.5 * waferW;
2123  double R1 = 2.0 * r1 / sqrt3_;
2124  int N = (r == 0) ? 3 : (static_cast<int>(0.5 * rout / r) + 3);
2125  int ns1 = (2 * N + 1) * (2 * N + 1);
2126  int ns2 = ns1 * php.zLayerHex_.size();
2127 #ifdef EDM_ML_DEBUG
2128  edm::LogVerbatim("HGCalGeom") << "wafer " << waferW << ":" << waferS << " r " << r << " dy " << dy << " N " << N
2129  << " sizes " << ns1 << ":" << ns2;
2130  std::vector<int> indtypes(ns1 + 1);
2131  indtypes.clear();
2132 #endif
2133  HGCalParameters::wafer_map wafersInLayers(ns1 + 1);
2134  HGCalParameters::wafer_map typesInLayers(ns2 + 1);
2135  HGCalParameters::waferT_map waferTypes(ns2 + 1);
2136  int ipos(0), lpos(0), uvmax(0), nwarn(0);
2137  std::vector<int> uvmx(php.zLayerHex_.size(), 0);
2138  for (int v = -N; v <= N; ++v) {
2139  for (int u = -N; u <= N; ++u) {
2140  int nr = 2 * v;
2141  int nc = -2 * u + v;
2142  double xpos = nc * r;
2143  double ypos = nr * dy;
2144  int indx = HGCalWaferIndex::waferIndex(0, u, v);
2145  php.waferCopy_.emplace_back(indx);
2146  php.waferPosX_.emplace_back(xpos);
2147  php.waferPosY_.emplace_back(ypos);
2148  wafersInLayers[indx] = ipos;
2149  ++ipos;
2150  std::pair<int, int> corner = HGCalGeomTools::waferCorner(xpos, ypos, r1, R1, 0, rout, false);
2151  if ((corner.first == static_cast<int>(HGCalParameters::k_CornerSize)) ||
2152  ((corner.first > 0) && php.defineFull_)) {
2153  uvmax = std::max(uvmax, std::max(std::abs(u), std::abs(v)));
2154  }
2155  for (unsigned int i = 0; i < php.zLayerHex_.size(); ++i) {
2156  int copy = i + php.layerOffset_;
2157  std::pair<double, double> xyoff = geomTools_.shiftXY(php.layerCenter_[copy], (waferW + waferS));
2158  int lay = php.layer_[php.layerIndex_[i]];
2159  double xpos0 = xpos + xyoff.first;
2160  double ypos0 = ypos + xyoff.second;
2161  double zpos = php.zLayerHex_[i];
2162  int kndx = HGCalWaferIndex::waferIndex(lay, u, v);
2163  int type(-1);
2166  type = wType->getType(kndx, php.waferInfoMap_);
2167  if (type < 0)
2168  type = wType->getType(HGCalParameters::k_ScaleToDDD * xpos0,
2171  php.waferTypeL_.emplace_back(type);
2172  typesInLayers[kndx] = lpos;
2173  ++lpos;
2174 #ifdef EDM_ML_DEBUG
2175  indtypes.emplace_back(kndx);
2176 #endif
2177  std::pair<int, int> corner =
2178  HGCalGeomTools::waferCorner(xpos0, ypos0, r1, R1, php.rMinLayHex_[i], php.rMaxLayHex_[i], false);
2179 #ifdef EDM_ML_DEBUG
2180  if (((corner.first == 0) && std::abs(u) < 5 && std::abs(v) < 5) || (std::abs(u) < 2 && std::abs(v) < 2)) {
2181  edm::LogVerbatim("HGCalGeom") << "Layer " << lay << " R " << php.rMinLayHex_[i] << ":" << php.rMaxLayHex_[i]
2182  << " u " << u << " v " << v << " with " << corner.first << " corners";
2183  }
2184 #endif
2185  if ((corner.first == static_cast<int>(HGCalParameters::k_CornerSize)) ||
2186  ((corner.first > 0) && php.defineFull_)) {
2187  uvmx[i] = std::max(uvmx[i], std::max(std::abs(u), std::abs(v)));
2188  }
2189  if ((corner.first < static_cast<int>(HGCalParameters::k_CornerSize)) && (corner.first > 0)) {
2190 #ifdef EDM_ML_DEBUG
2191  edm::LogVerbatim("HGCalGeom") << "Layer " << lay << " u|v " << u << ":" << v << " with corner "
2192  << corner.first << ":" << corner.second;
2193 #endif
2194  int wl = HGCalWaferIndex::waferIndex(lay, u, v);
2195  if (php.waferMaskMode_ > 0) {
2196  bool v17OrLess = (php.mode_ < HGCalGeometryMode::Hexagon8CalibCell);
2197  std::pair<int, int> corner0 = HGCalWaferMask::getTypeMode(
2198  xpos0, ypos0, r1, R1, php.rMinLayHex_[i], php.rMaxLayHex_[i], type, php.waferMaskMode_, v17OrLess);
2202  auto itr = php.waferInfoMap_.find(wl);
2203  if (itr != php.waferInfoMap_.end()) {
2204  int part = (itr->second).part;
2205  int orient = (itr->second).orient;
2206  bool ok = ((php.mode_ == HGCalGeometryMode::Hexagon8Cassette) ||
2208  ? true
2210  xpos0, ypos0, r1, R1, php.rMinLayHex_[i], php.rMaxLayHex_[i], part, orient, false);
2211  if (ok)
2212  corner0 = std::make_pair(part, (HGCalTypes::k_OffsetRotation + orient));
2213 #ifdef EDM_ML_DEBUG
2214  edm::LogVerbatim("HGCalGeom")
2215  << "Layer:u:v " << i << ":" << lay << ":" << u << ":" << v << " Part " << corner0.first << ":"
2216  << part << " Orient " << corner0.second << ":" << orient << " Position " << xpos0 << ":" << ypos0
2217  << " delta " << r1 << ":" << R1 << " Limit " << php.rMinLayHex_[i] << ":" << php.rMaxLayHex_[i]
2218  << " Compatibiliety Flag " << ok;
2219 #endif
2220  if (!ok)
2221  ++nwarn;
2222  }
2223  }
2224  waferTypes[wl] = corner0;
2225 #ifdef EDM_ML_DEBUG
2226  edm::LogVerbatim("HGCalGeom")
2227  << "Layer " << lay << " u|v " << u << ":" << v << " Index " << std::hex << wl << std::dec << " pos "
2228  << xpos0 << ":" << ypos0 << " R " << r1 << ":" << R1 << " Range " << php.rMinLayHex_[i] << ":"
2229  << php.rMaxLayHex_[i] << type << ":" << php.waferMaskMode_ << " corner " << corner.first << ":"
2230  << corner.second << " croner0 " << corner0.first << ":" << corner0.second;
2231 #endif
2232  } else {
2233  waferTypes[wl] = corner;
2234 #ifdef EDM_ML_DEBUG
2235  edm::LogVerbatim("HGCalGeom") << "Layer " << lay << " u|v " << u << ":" << v << " with corner "
2236  << corner.first << ":" << corner.second;
2237 #endif
2238  }
2239  }
2240  }
2241  }
2242  }
2243  if (nwarn > 0)
2244  edm::LogWarning("HGCalGeom") << "HGCalGeomParameters::loadWafer8: there are " << nwarn
2245  << " wafers with non-matching partial- orientation types";
2246  php.waferUVMax_ = uvmax;
2247  php.waferUVMaxLayer_ = uvmx;
2248  php.wafersInLayers_ = wafersInLayers;
2249  php.typesInLayers_ = typesInLayers;
2250  php.waferTypes_ = waferTypes;
2251  php.nSectors_ = static_cast<int>(php.waferCopy_.size());
2253  mytr.lay = 1;
2254  mytr.bl = php.waferR_;
2255  mytr.tl = php.waferR_;
2256  mytr.h = php.waferR_;
2257  mytr.alpha = 0.0;
2259  for (auto const& dz : php.cellThickness_) {
2260  mytr.dz = 0.5 * HGCalParameters::k_ScaleToDDD * dz;
2261  php.fillModule(mytr, false);
2262  }
2263  for (unsigned k = 0; k < php.cellThickness_.size(); ++k) {
2264  HGCalParameters::hgtrap mytr = php.getModule(k, false);
2270  php.fillModule(mytr, true);
2271  }
2272 #ifdef EDM_ML_DEBUG
2273  edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: Total of " << php.waferCopy_.size() << " wafers";
2274  for (unsigned int k = 0; k < php.waferCopy_.size(); ++k) {
2275  int id = php.waferCopy_[k];
2276  edm::LogVerbatim("HGCalGeom") << "[" << k << "] " << std::hex << id << std::dec << ":"
2277  << HGCalWaferIndex::waferLayer(id) << ":" << HGCalWaferIndex::waferU(id) << ":"
2278  << HGCalWaferIndex::waferV(id) << " x " << php.waferPosX_[k] << " y "
2279  << php.waferPosY_[k] << " index " << php.wafersInLayers_[id];
2280  }
2281  edm::LogVerbatim("HGCalGeom") << "HGCalParameters: Total of " << php.waferTypeL_.size() << " wafer types";
2282  for (unsigned int k = 0; k < php.waferTypeL_.size(); ++k) {
2283  int id = indtypes[k];
2284  edm::LogVerbatim("HGCalGeom") << "[" << k << "] " << php.typesInLayers_[id] << ":" << php.waferTypeL_[k] << " ID "
2285  << std::hex << id << std::dec << ":" << HGCalWaferIndex::waferLayer(id) << ":"
2286  << HGCalWaferIndex::waferU(id) << ":" << HGCalWaferIndex::waferV(id);
2287  }
2288 #endif
2289 
2290  //Wafer offset
2291  php.xLayerHex_.clear();
2292  php.yLayerHex_.clear();
2293  double waferSize = php.waferSize_ + php.sensorSeparation_;
2294 #ifdef EDM_ML_DEBUG
2295  edm::LogVerbatim("HGCalGeom") << "WaferSize " << waferSize;
2296 #endif
2297  for (unsigned int k = 0; k < php.zLayerHex_.size(); ++k) {
2298  int copy = k + php.layerOffset_;
2299  std::pair<double, double> xyoff = geomTools_.shiftXY(php.layerCenter_[copy], waferSize);
2300  php.xLayerHex_.emplace_back(xyoff.first);
2301  php.yLayerHex_.emplace_back(xyoff.second);
2302 #ifdef EDM_ML_DEBUG
2303  edm::LogVerbatim("HGCalGeom") << "Layer[" << k << "] Off " << copy << ":" << php.layerCenter_[copy] << " Shift "
2304  << xyoff.first << ":" << xyoff.second;
2305 #endif
2306  }
2307 }
2308 
2310  // Special parameters for cell parameters
2311  std::string attribute = "OnlyForHGCalNumbering";
2312  DDSpecificsHasNamedValueFilter filter1{attribute};
2313  DDFilteredView fv1(*cpv, filter1);
2314  bool ok = fv1.firstChild();
2315 
2316  if (ok) {
2317  php.cellFine_ = dbl_to_int(cpv->vector("waferFine"));
2318  php.cellCoarse_ = dbl_to_int(cpv->vector("waferCoarse"));
2319  }
2320 
2321  loadCellParsHexagon(php);
2322 }
2323 
2325  for (auto const& it : vmap) {
2326  if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "waferFine")) {
2327  for (const auto& i : it.second)
2328  php.cellFine_.emplace_back(std::round(i));
2329  } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "waferCoarse")) {
2330  for (const auto& i : it.second)
2331  php.cellCoarse_.emplace_back(std::round(i));
2332  }
2333  }
2334 
2335  loadCellParsHexagon(php);
2336 }
2337 
2339 #ifdef EDM_ML_DEBUG
2340  edm::LogVerbatim("HGCalGeom") << "HGCalLoadCellPars: " << php.cellFine_.size() << " rows for fine cells";
2341  for (unsigned int k = 0; k < php.cellFine_.size(); ++k)
2342  edm::LogVerbatim("HGCalGeom") << "[" << k << "]: " << php.cellFine_[k];
2343  edm::LogVerbatim("HGCalGeom") << "HGCalLoadCellPars: " << php.cellCoarse_.size() << " rows for coarse cells";
2344  for (unsigned int k = 0; k < php.cellCoarse_.size(); ++k)
2345  edm::LogVerbatim("HGCalGeom") << "[" << k << "]: " << php.cellCoarse_[k];
2346 #endif
2347 }
2348 
2350  php.xLayerHex_.resize(php.zLayerHex_.size(), 0);
2351  php.yLayerHex_.resize(php.zLayerHex_.size(), 0);
2352 #ifdef EDM_ML_DEBUG
2353  edm::LogVerbatim("HGCalGeom") << "HGCalParameters: x|y|zLayerHex in array of size " << php.zLayerHex_.size();
2354  for (unsigned int k = 0; k < php.zLayerHex_.size(); ++k)
2355  edm::LogVerbatim("HGCalGeom") << "Layer[" << k << "] Shift " << php.xLayerHex_[k] << ":" << php.yLayerHex_[k] << ":"
2356  << php.zLayerHex_[k];
2357 #endif
2358  // Find the radius of each eta-partitions
2361  //Ring radii for each partition
2362 #ifdef EDM_ML_DEBUG
2363  edm::LogVerbatim("HGCalGeom") << "HGCalParameters: Mode " << php.mode_ << ":"
2364  << HGCalGeometryMode::TrapezoidFineCell << " Sizes " << php.tileRingFineR_.size()
2365  << ":" << php.tileRingR_.size();
2366 #endif
2367  for (unsigned int k = 0; k < 2; ++k) {
2368  bool fine = ((k == 0) && (php.mode_ == HGCalGeometryMode::TrapezoidFineCell));
2369  unsigned int sizeR = (fine) ? php.tileRingFineR_.size() : php.tileRingR_.size();
2370  for (unsigned int kk = 0; kk < sizeR; ++kk) {
2371  if (fine)
2372  php.radiusLayer_[k].emplace_back(php.tileRingFineR_[kk].first);
2373  else
2374  php.radiusLayer_[k].emplace_back(php.tileRingR_[kk].first);
2375 #ifdef EDM_ML_DEBUG
2376  double zv = ((k == 0) ? (php.zLayerHex_[php.layerFrontBH_[1] - php.firstLayer_])
2377  : (php.zLayerHex_[php.zLayerHex_.size() - 1]));
2378  double rv = php.radiusLayer_[k].back();
2379  double eta = -(std::log(std::tan(0.5 * std::atan(rv / zv))));
2380  edm::LogVerbatim("HGCalGeom") << "New [" << kk << "] new R = " << rv << " Eta = " << eta;
2381 #endif
2382  }
2383  if (fine)
2384  php.radiusLayer_[k].emplace_back(php.tileRingFineR_[php.tileRingFineR_.size() - 1].second);
2385  else
2386  php.radiusLayer_[k].emplace_back(php.tileRingR_[php.tileRingR_.size() - 1].second);
2387  }
2388  // Minimum and maximum radius index for each layer
2389 #ifdef EDM_ML_DEBUG
2390  edm::LogVerbatim("HGCalGeom") << "Till Ring Range size " << php.tileRingFineRange_.size() << ":"
2391  << php.tileRingRange_.size() << ":" << php.nPhiBinBH_.size() << ":"
2392  << php.zLayerHex_.size() << ":" << php.nPhiLayer_.size();
2393 #endif
2394  unsigned int k1(0), k2(0);
2395  for (unsigned int k = 0; k < php.zLayerHex_.size(); ++k) {
2396  if (php.nPhiLayer_[k] > 288) {
2397  php.iradMinBH_.emplace_back(1 + php.tileRingFineRange_[k1].first);
2398  php.iradMaxBH_.emplace_back(1 + php.tileRingFineRange_[k1].second);
2399  ++k1;
2400  } else {
2401  php.iradMinBH_.emplace_back(1 + php.tileRingRange_[k2].first);
2402  php.iradMaxBH_.emplace_back(1 + php.tileRingRange_[k2].second);
2403  ++k2;
2404  }
2405 #ifdef EDM_ML_DEBUG
2406  int kk = php.scintType(php.firstLayer_ + static_cast<int>(k));
2407  edm::LogVerbatim("HGCalGeom") << "New Layer " << k << " Type " << kk << " Low edge " << php.iradMinBH_.back()
2408  << " Top edge " << php.iradMaxBH_.back();
2409 #endif
2410  }
2411  } else {
2412  //Ring radii for each partition
2413  for (unsigned int k = 0; k < 2; ++k) {
2414  double rmax = ((k == 0) ? (php.rMaxLayHex_[php.layerFrontBH_[1] - php.firstLayer_] - 1)
2415  : (php.rMaxLayHex_[php.rMaxLayHex_.size() - 1]));
2416  double rv = php.rMinLayerBH_[k];
2417  double zv = ((k == 0) ? (php.zLayerHex_[php.layerFrontBH_[1] - php.firstLayer_])
2418  : (php.zLayerHex_[php.zLayerHex_.size() - 1]));
2419  php.radiusLayer_[k].emplace_back(rv);
2420 #ifdef EDM_ML_DEBUG
2421  double eta = -(std::log(std::tan(0.5 * std::atan(rv / zv))));
2422  edm::LogVerbatim("HGCalGeom") << "Old [" << k << "] rmax " << rmax << " Z = " << zv
2423  << " dEta = " << php.cellSize_[k] << "\nOld[0] new R = " << rv << " Eta = " << eta;
2424  int kount(1);
2425 #endif
2426  while (rv < rmax) {
2427  double eta = -(php.cellSize_[k] + std::log(std::tan(0.5 * std::atan(rv / zv))));
2428  rv = zv * std::tan(2.0 * std::atan(std::exp(-eta)));
2429  php.radiusLayer_[k].emplace_back(rv);
2430 #ifdef EDM_ML_DEBUG
2431  edm::LogVerbatim("HGCalGeom") << "Old [" << kount << "] new R = " << rv << " Eta = " << eta;
2432  ++kount;
2433 #endif
2434  }
2435  }
2436  // Find minimum and maximum radius index for each layer
2437  for (unsigned int k = 0; k < php.zLayerHex_.size(); ++k) {
2438  int kk = php.scintType(php.firstLayer_ + static_cast<int>(k));
2439  std::vector<double>::iterator low, high;
2440  low = std::lower_bound(php.radiusLayer_[kk].begin(), php.radiusLayer_[kk].end(), php.rMinLayHex_[k]);
2441 #ifdef EDM_ML_DEBUG
2442  edm::LogVerbatim("HGCalGeom") << "Old [" << k << "] RLow = " << php.rMinLayHex_[k] << " pos "
2443  << static_cast<int>(low - php.radiusLayer_[kk].begin());
2444 #endif
2445  if (low == php.radiusLayer_[kk].begin())
2446  ++low;
2447  int irlow = static_cast<int>(low - php.radiusLayer_[kk].begin());
2448  double drlow = php.radiusLayer_[kk][irlow] - php.rMinLayHex_[k];
2449 #ifdef EDM_ML_DEBUG
2450  edm::LogVerbatim("HGCalGeom") << "irlow " << irlow << " dr " << drlow << " min " << php.minTileSize_;
2451 #endif
2452  if (drlow < php.minTileSize_) {
2453  ++irlow;
2454 #ifdef EDM_ML_DEBUG
2455  drlow = php.radiusLayer_[kk][irlow] - php.rMinLayHex_[k];
2456  edm::LogVerbatim("HGCalGeom") << "Modified irlow " << irlow << " dr " << drlow;
2457 #endif
2458  }
2459  high = std::lower_bound(php.radiusLayer_[kk].begin(), php.radiusLayer_[kk].end(), php.rMaxLayHex_[k]);
2460 #ifdef EDM_ML_DEBUG
2461  edm::LogVerbatim("HGCalGeom") << "Old [" << k << "] RHigh = " << php.rMaxLayHex_[k] << " pos "
2462  << static_cast<int>(high - php.radiusLayer_[kk].begin());
2463 #endif
2464  if (high == php.radiusLayer_[kk].end())
2465  --high;
2466  int irhigh = static_cast<int>(high - php.radiusLayer_[kk].begin());
2467  double drhigh = php.rMaxLayHex_[k] - php.radiusLayer_[kk][irhigh - 1];
2468 #ifdef EDM_ML_DEBUG
2469  edm::LogVerbatim("HGCalGeom") << "irhigh " << irhigh << " dr " << drhigh << " min " << php.minTileSize_;
2470 #endif
2471  if (drhigh < php.minTileSize_) {
2472  --irhigh;
2473 #ifdef EDM_ML_DEBUG
2474  drhigh = php.rMaxLayHex_[k] - php.radiusLayer_[kk][irhigh - 1];
2475  edm::LogVerbatim("HGCalGeom") << "Modified irhigh " << irhigh << " dr " << drhigh;
2476 #endif
2477  }
2478  php.iradMinBH_.emplace_back(irlow);
2479  php.iradMaxBH_.emplace_back(irhigh);
2480 #ifdef EDM_ML_DEBUG
2481  edm::LogVerbatim("HGCalGeom") << "Old Layer " << k << " Type " << kk << " Low edge " << irlow << ":" << drlow
2482  << " Top edge " << irhigh << ":" << drhigh;
2483 #endif
2484  }
2485  }
2486 #ifdef EDM_ML_DEBUG
2487  for (unsigned int k = 0; k < 2; ++k) {
2488  edm::LogVerbatim("HGCalGeom") << "Type " << k << " with " << php.radiusLayer_[k].size() << " radii";
2489  for (unsigned int kk = 0; kk < php.radiusLayer_[k].size(); ++kk)
2490  edm::LogVerbatim("HGCalGeom") << "Ring[" << kk << "] " << php.radiusLayer_[k][kk];
2491  }
2492 #endif
2493 
2494  // Now define the volumes
2495  int im(0);
2496  php.waferUVMax_ = 0;
2498  mytr.alpha = 0.0;
2499  for (unsigned int k = 0; k < php.zLayerHex_.size(); ++k) {
2500  if (php.iradMaxBH_[k] > php.waferUVMax_)
2501  php.waferUVMax_ = php.iradMaxBH_[k];
2502  int kk = ((php.firstLayer_ + static_cast<int>(k)) < php.layerFrontBH_[1]) ? 0 : 1;
2503  int irm = php.radiusLayer_[kk].size() - 1;
2504 #ifdef EDM_ML_DEBUG
2505  double rmin = php.radiusLayer_[kk][std::max((php.iradMinBH_[k] - 1), 0)];
2506  double rmax = php.radiusLayer_[kk][std::min(php.iradMaxBH_[k], irm)];
2507  edm::LogVerbatim("HGCalGeom") << "Layer " << php.firstLayer_ + k << ":" << kk << " Radius range "
2508  << php.iradMinBH_[k] << ":" << php.iradMaxBH_[k] << ":" << rmin << ":" << rmax;
2509 #endif
2510  mytr.lay = php.firstLayer_ + k;
2511  for (int irad = php.iradMinBH_[k]; irad <= php.iradMaxBH_[k]; ++irad) {
2512  double rmin = php.radiusLayer_[kk][std::max((irad - 1), 0)];
2513  double rmax = php.radiusLayer_[kk][std::min(irad, irm)];
2514  mytr.bl = 0.5 * rmin * php.scintCellSize(mytr.lay);
2515  mytr.tl = 0.5 * rmax * php.scintCellSize(mytr.lay);
2516  mytr.h = 0.5 * (rmax - rmin);
2517  mytr.dz = 0.5 * php.waferThick_;
2518  mytr.cellSize = 0.5 * (rmax + rmin) * php.scintCellSize(mytr.lay);
2519  php.fillModule(mytr, true);
2525  php.fillModule(mytr, false);
2526  if (irad == php.iradMinBH_[k])
2527  php.firstModule_.emplace_back(im);
2528  ++im;
2529  if (irad == php.iradMaxBH_[k] - 1)
2530  php.lastModule_.emplace_back(im);
2531  }
2532  }
2533  php.nSectors_ = php.waferUVMax_;
2534 #ifdef EDM_ML_DEBUG
2535  edm::LogVerbatim("HGCalGeom") << "Maximum radius index " << php.waferUVMax_;
2536  for (unsigned int k = 0; k < php.firstModule_.size(); ++k)
2537  edm::LogVerbatim("HGCalGeom") << "Layer " << k + php.firstLayer_ << " Modules " << php.firstModule_[k] << ":"
2538  << php.lastModule_[k];
2539 #endif
2540 }
2541 
2542 std::vector<double> HGCalGeomParameters::getDDDArray(const std::string& str, const DDsvalues_type& sv, const int nmin) {
2543  DDValue value(str);
2544  if (DDfetch(&sv, value)) {
2545  const std::vector<double>& fvec = value.doubles();
2546  int nval = fvec.size();
2547  if (nmin > 0) {
2548  if (nval < nmin) {
2549  throw cms::Exception("DDException")
2550  << "HGCalGeomParameters: # of " << str << " bins " << nval << " < " << nmin << " ==> illegal";
2551  }
2552  } else {
2553  if (nval < 1 && nmin == 0) {
2554  throw cms::Exception("DDException")
2555  << "HGCalGeomParameters: # of " << str << " bins " << nval << " < 1 ==> illegal"
2556  << " (nmin=" << nmin << ")";
2557  }
2558  }
2559  return fvec;
2560  } else {
2561  if (nmin >= 0) {
2562  throw cms::Exception("DDException") << "HGCalGeomParameters: cannot get array " << str;
2563  }
2564  std::vector<double> fvec;
2565  return fvec;
2566  }
2567 }
2568 
2569 std::pair<double, double> HGCalGeomParameters::cellPosition(
2570  const std::vector<HGCalGeomParameters::cellParameters>& wafers,
2571  std::vector<HGCalGeomParameters::cellParameters>::const_iterator& itrf,
2572  int wafer,
2573  double xx,
2574  double yy) {
2575  if (itrf == wafers.end()) {
2576  for (std::vector<HGCalGeomParameters::cellParameters>::const_iterator itr = wafers.begin(); itr != wafers.end();
2577  ++itr) {
2578  if (itr->wafer == wafer) {
2579  itrf = itr;
2580  break;
2581  }
2582  }
2583  }
2584  double dx(0), dy(0);
2585  if (itrf != wafers.end()) {
2586  dx = (xx - itrf->xyz.x());
2587  if (std::abs(dx) < tolerance)
2588  dx = 0;
2589  dy = (yy - itrf->xyz.y());
2590  if (std::abs(dy) < tolerance)
2591  dy = 0;
2592  }
2593  return std::make_pair(dx, dy);
2594 }
2595 
2596 void HGCalGeomParameters::rescale(std::vector<double>& v, const double s) {
2597  std::for_each(v.begin(), v.end(), [s](double& n) { n *= s; });
2598 }
2599 
2600 void HGCalGeomParameters::resetZero(std::vector<double>& v) {
2601  for (auto& n : v) {
2602  if (std::abs(n) < tolmin)
2603  n = 0;
2604  }
2605 }
std::vector< int > iradMaxBH_
std::vector< double > waferPosY_
Log< level::Info, true > LogVerbatim
std::vector< int > layer_
static constexpr int scintillatorCassette
std::vector< double > moduleDzR_
std::vector< int > depthLayerF_
std::vector< int > depth_
std::vector< double > zFrontMin_
int32_t tileType(const int32_t property)
nav_type copyNumbers() const
return the stack of copy numbers
hgtrap getModule(unsigned int k, bool reco) const
std::vector< double > moduleHR_
std::vector< double > rMinVec(void) const
Definition: DDSolid.cc:335
layer_map copiesInLayers_
static constexpr int k_OffsetRotation
Definition: HGCalTypes.h:93
std::vector< double > const & vector(std::string_view iKey) const
The DDVector information.
std::vector< std::pair< double, double > > layerRotV_
std::vector< bool > cellCoarseHalf_
int scintType(const int layer) const
std::vector< bool > cellFineHalf_
const std::vector< int > copyNos() const
The list of the volume copy numbers.
const cms::DDSolidShape shape() const
std::vector< int > moduleLayR_
double rIn(void) const
Definition: DDSolid.cc:456
static constexpr int siliconFileEE
int32_t waferU(const int32_t index)
std::vector< int > cellFine_
void rescale(std::vector< double > &, const double s)
for(int i=first, nt=offsets[nh];i< nt;i+=gridDim.x *blockDim.x)
int32_t waferLayer(const int32_t index)
const double tolerance
static int32_t getUnpackedCell6(int id)
Definition: HGCalTypes.cc:40
std::vector< double > moduleHS_
std::string name(Mapping a, V value)
Definition: DDSolidShapes.h:31
std::vector< double > xVec(void) const
Definition: DDSolid.cc:378
std::vector< double > const & vector(std::string_view iKey) const
returns an empty container if not found
std::vector< double > zVec(void) const
Definition: DDSolid.cc:328
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
std::vector< double > trformTranY_
void loadGeometryHexagonModule(const DDCompactView *cpv, HGCalParameters &php, const std::string &sdTag1, const std::string &sdTag2, int firstLayer)
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
std::vector< double > rMaxVec(void) const
Definition: DDSolid.cc:342
double rOut(void) const
Definition: DDSolid.cc:458
std::vector< double > cellFineY_
void loadSpecParsHexagon(const DDFilteredView &fv, HGCalParameters &php, const DDCompactView *cpv, const std::string &sdTag1, const std::string &sdTag2)
const std::string & name() const
The name of a logical-part of the current node in the filtered-view.
std::vector< double > trformRotZY_
int zside(DetId const &)
static int getRotation(int zside, int type, int rotn)
void resetZero(std::vector< double > &)
std::vector< uint32_t > trformIndex_
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:19
assert(be >=bs)
std::vector< int > layerGroupM_
Compact representation of the geometrical detector hierarchy.
Definition: DDCompactView.h:81
std::vector< std::pair< int, int > > tileRingFineRange_
bool DDfetch(const DDsvalues_type *, DDValue &)
helper for retrieving DDValues from DDsvalues_type *.
Definition: DDsvalues.cc:79
HGCalGeometryMode::GeometryMode mode_
const std::string names[nVars_]
int32_t waferOrient(const int32_t property)
std::pair< double, double > shiftXY(int waferPosition, double waferSize) const
std::vector< double > trformRotXX_
std::vector< int > nPhiBinBH_
std::unordered_map< std::string, std::vector< double > > DDVectorsMap
Definition: DDNamespace.h:20
A DDSolid represents the shape of a part.
Definition: DDSolid.h:39
static constexpr double tol
void fillTrForm(const hgtrform &mytr)
void loadSpecParsHexagon8(const DDFilteredView &fv, HGCalParameters &php)
static constexpr uint32_t k_CornerSize
std::vector< double > zVec(void) const
Definition: DDSolid.cc:387
wafer_map wafersInLayers_
static int32_t getUnpackedCellType6(int id)
Definition: HGCalTypes.cc:38
std::vector< double > rMinLayerBH_
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double > > DD3Vector
std::vector< double > trformRotZX_
double scintCellSize(const int layer) const
std::string_view name() const
std::vector< double > cellCoarseX_
static constexpr int scintillatorFile
const double tolmin
int32_t waferCassette(const int32_t property)
std::vector< std::pair< int, int > > tileRingRange_
static constexpr int siliconCassetteHE
static std::pair< int32_t, int32_t > waferCorner(double xpos, double ypos, double r, double R, double rMin, double rMax, bool oldBug=false)
std::vector< int > firstModule_
std::vector< int > cellCoarse_
std::vector< double > trformRotYZ_
std::vector< std::pair< unsigned int, DDValue > > DDsvalues_type
Definition: DDsvalues.h:12
void loadCellTrapezoid(HGCalParameters &php)
std::pair< double, double > cellPosition(const std::vector< cellParameters > &wafers, std::vector< cellParameters >::const_iterator &itrf, int wafer, double xx, double yy)
std::vector< double > boundR_
const std::string & name() const
Returns the name.
Definition: DDName.cc:41
std::vector< double > cellSize_
std::vector< int > waferUVMaxLayer_
std::vector< double > moduleDzS_
bool next()
set current node to the next node in the filtered tree
std::vector< int > layerIndex_
std::vector< double > moduleAlphaR_
static std::pair< int, int > getTypeMode(const double &xpos, const double &ypos, const double &delX, const double &delY, const double &rin, const double &rout, const int &waferType, const int &mode, const bool &v17, const bool &debug=false)
std::vector< double > yLayerHex_
void loadCellParsHexagon(const DDCompactView *cpv, HGCalParameters &php)
T sqrt(T t)
Definition: SSEVec.h:23
std::vector< double > trformRotXY_
std::vector< double > rMaxFront_
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
std::vector< double > trformRotYX_
static constexpr int siliconCassetteEE
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
std::vector< double > getDDDArray(const std::string &str, const DDsvalues_type &sv, const int nmin)
const int level() const
get Iterator level
int32_t tileSiPM(const int32_t property)
std::vector< double > waferThickness_
std::vector< double > slopeTop_
std::vector< int > calibCellPartHD_
std::vector< int > layerCenter_
static constexpr double k_ScaleFromDD4hep
std::vector< double > moduleBlR_
int32_t waferThick(const int32_t property)
const std::array< const cms::dd::NameValuePair< DDSolidShape >, 21 > DDSolidShapeMap
Definition: DDSolidShapes.h:99
Definition: value.py:1
std::vector< double > rMinLayHex_
int32_t waferIndex(int32_t layer, int32_t waferU, int32_t waferV, bool old=false)
void fillModule(const hgtrap &mytr, bool reco)
std::vector< double > moduleTlS_
bool firstChild()
set the current node to the first child
std::vector< int > calibCellFullHD_
std::vector< double > zLayerHex_
std::vector< int > layerType_
std::vector< std::pair< double, double > > tileRingFineR_
#define M_PI
std::vector< int > dbl_to_int(const std::vector< double > &vecdbl)
Converts a std::vector of doubles to a std::vector of int.
Definition: DDutils.h:7
static constexpr double tan30deg_
void loadWaferHexagon(HGCalParameters &php)
waferT_map waferTypes_
DDsvalues_type mergedSpecifics() const
const N & name() const
Definition: DDBase.h:58
std::vector< double > rMaxLayHex_
std::vector< double > trformTranX_
static constexpr double k_ScaleFromDD4hepToG4
std::unordered_map< int32_t, std::pair< int32_t, int32_t > > waferT_map
const DDLogicalPart & logicalPart() const
The logical-part of the current node in the filtered-view.
std::vector< double > zRanges_
std::vector< double > slopeMin_
static constexpr int32_t WaferCenterR
Definition: HGCalTypes.h:27
T get(const std::string &)
extract attribute value
std::vector< int > lastModule_
static constexpr double k_ScaleToDDD
static bool goodTypeMode(const double &xpos, const double &ypos, const double &delX, const double &delY, const double &rin, const double &rout, const int &part, const int &rotn, const bool &v17, const bool &debug=false)
std::vector< double > radiusMixBoundary_
constexpr float sol
Definition: Config.h:13
#define N
Definition: blowfish.cc:9
std::vector< double > cellThickness_
std::vector< double > trformRotXZ_
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double > > XYZPointD
point in space with cartesian internal representation
Definition: Point3D.h:8
void scaleTrForm(double)
std::vector< int > layerGroup_
std::unordered_map< int32_t, int32_t > wafer_map
Interface to a Box.
Definition: DDSolid.h:167
part
Definition: HCALResponse.h:20
std::vector< double > radius200to300_
std::vector< double > radius100to200_
const DDSolid & solid(void) const
Returns a reference object of the solid being the shape of this LogicalPart.
std::vector< double > rMinFront_
dd4hep::Solid solid() const
std::vector< int > iradMinBH_
std::vector< double > trformRotYY_
std::vector< double > cellFineX_
static constexpr double k_ScaleFromDDDToG4
wafer_map typesInLayers_
const DDRotationMatrix & rotation() const
The absolute rotation of the current node.
std::vector< double > trformRotZZ_
int32_t waferPartial(const int32_t property)
const RotationMatrix rotation() const
static constexpr double k_ScaleFromDDD
double halfX(void) const
Definition: DDSolid.cc:212
std::vector< double > moduleAlphaS_
std::vector< int > layerGroupO_
void loadGeometryHexagon(const DDFilteredView &_fv, HGCalParameters &php, const std::string &sdTag1, const DDCompactView *cpv, const std::string &sdTag2, const std::string &sdTag3, HGCalGeometryMode::WaferMode mode)
void loadSpecParsTrapezoid(const DDFilteredView &fv, HGCalParameters &php)
std::vector< double > moduleBlS_
int32_t waferV(const int32_t index)
bool firstChild()
set the current node to the first child ...
std::vector< int > waferCopy_
std::vector< double > cassetteShift_
std::vector< int > depthIndex_
std::vector< int > layerFrontBH_
std::vector< int > calibCellPartLD_
static int32_t packTypeUV(int type, int u, int v)
Definition: HGCalTypes.cc:3
std::vector< double > rLimit_
std::vector< double > zFrontTop_
std::vector< double > radiusLayer_[2]
std::vector< int > waferTypeT_
std::vector< int > levelT_
std::vector< double > cellCoarseY_
Log< level::Warning, false > LogWarning
std::vector< std::pair< double, double > > tileRingR_
std::vector< int > nPhiLayer_
waferInfo_map waferInfoMap_
static int32_t layerType(int type)
Definition: HGCalTypes.cc:42
std::vector< int > moduleLayS_
static constexpr int siliconFileHE
std::vector< double > trformTranZ_
const DDTranslation & translation() const
The absolute translation of the current node.
#define str(s)
CLHEP::HepRotation hr
std::vector< double > waferPosX_
int scintCells(const int layer) const
void addTrForm(const CLHEP::Hep3Vector &h3v)
tileInfo_map tileInfoMap_
const std::vector< double > parameters() const
extract shape parameters
if(threadIdxLocalY==0 &&threadIdxLocalX==0)
std::vector< double > moduleTlR_
std::vector< int > waferTypeL_
std::vector< double > cassetteShiftTile_
std::vector< double > xLayerHex_
int32_t tileProperty(const int32_t type, const int32_t sipm)
void loadWaferHexagon8(HGCalParameters &php)
std::vector< std::unordered_map< int32_t, int32_t > > layer_map
void loadGeometryHexagon8(const DDFilteredView &_fv, HGCalParameters &php, int firstLayer)
def cp(fromDir, toDir, listOfFiles, overwrite=False, smallList=False)
std::vector< int > calibCellFullLD_
const Translation translation() const
static constexpr int scintillatorFineCell