File indexing completed on 2025-05-02 02:40:40
0001 #include "Geometry/HGCalCommonData/interface/HGCalGeomParameters.h"
0002
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 #include "FWCore/Utilities/interface/Exception.h"
0005
0006 #include "DataFormats/Math/interface/GeantUnits.h"
0007 #include "DataFormats/Math/interface/Point3D.h"
0008 #include "DetectorDescription/Core/interface/DDConstant.h"
0009 #include "DetectorDescription/Core/interface/DDFilter.h"
0010 #include "DetectorDescription/Core/interface/DDFilteredView.h"
0011 #include "DetectorDescription/Core/interface/DDSolid.h"
0012 #include "DetectorDescription/Core/interface/DDValue.h"
0013 #include "DetectorDescription/Core/interface/DDutils.h"
0014 #include "DetectorDescription/RegressionTest/interface/DDErrorDetection.h"
0015 #include "Geometry/HGCalCommonData/interface/HGCalProperty.h"
0016 #include "Geometry/HGCalCommonData/interface/HGCalTileIndex.h"
0017 #include "Geometry/HGCalCommonData/interface/HGCalTypes.h"
0018 #include "Geometry/HGCalCommonData/interface/HGCalWaferIndex.h"
0019 #include "Geometry/HGCalCommonData/interface/HGCalWaferMask.h"
0020 #include "Geometry/HGCalCommonData/interface/HGCalWaferType.h"
0021
0022 #include <algorithm>
0023 #include <sstream>
0024 #include <unordered_set>
0025
0026
0027 using namespace geant_units::operators;
0028
0029 const double tolerance = 0.001;
0030 const double tolmin = 1.e-20;
0031
0032 HGCalGeomParameters::HGCalGeomParameters() : sqrt3_(std::sqrt(3.0)) {
0033 #ifdef EDM_ML_DEBUG
0034 edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters::HGCalGeomParameters "
0035 << "constructor";
0036 #endif
0037 }
0038
0039 void HGCalGeomParameters::loadGeometryHexagon(const DDFilteredView& _fv,
0040 HGCalParameters& php,
0041 const std::string& sdTag1,
0042 const DDCompactView* cpv,
0043 const std::string& sdTag2,
0044 const std::string& sdTag3,
0045 HGCalGeometryMode::WaferMode mode) {
0046 DDFilteredView fv = _fv;
0047 bool dodet(true);
0048 std::map<int, HGCalGeomParameters::layerParameters> layers;
0049 std::vector<HGCalParameters::hgtrform> trforms;
0050 std::vector<bool> trformUse;
0051
0052 while (dodet) {
0053 const DDSolid& sol = fv.logicalPart().solid();
0054
0055 std::vector<int> copy = fv.copyNumbers();
0056 int nsiz = static_cast<int>(copy.size());
0057 int lay = (nsiz > 0) ? copy[nsiz - 1] : 0;
0058 int zp = (nsiz > 2) ? copy[nsiz - 3] : -1;
0059 if (zp != 1)
0060 zp = -1;
0061 if (lay == 0) {
0062 throw cms::Exception("DDException") << "Funny layer # " << lay << " zp " << zp << " in " << nsiz << " components";
0063 } else {
0064 if (std::find(php.layer_.begin(), php.layer_.end(), lay) == php.layer_.end())
0065 php.layer_.emplace_back(lay);
0066 auto itr = layers.find(lay);
0067 if (itr == layers.end()) {
0068 double rin(0), rout(0);
0069 double zz = HGCalParameters::k_ScaleFromDDD * fv.translation().Z();
0070 if ((sol.shape() == DDSolidShape::ddpolyhedra_rz) || (sol.shape() == DDSolidShape::ddpolyhedra_rrz)) {
0071 const DDPolyhedra& polyhedra = static_cast<DDPolyhedra>(sol);
0072 const std::vector<double>& rmin = polyhedra.rMinVec();
0073 const std::vector<double>& rmax = polyhedra.rMaxVec();
0074 rin = 0.5 * HGCalParameters::k_ScaleFromDDD * (rmin[0] + rmin[1]);
0075 rout = 0.5 * HGCalParameters::k_ScaleFromDDD * (rmax[0] + rmax[1]);
0076 } else if (sol.shape() == DDSolidShape::ddtubs) {
0077 const DDTubs& tube = static_cast<DDTubs>(sol);
0078 rin = HGCalParameters::k_ScaleFromDDD * tube.rIn();
0079 rout = HGCalParameters::k_ScaleFromDDD * tube.rOut();
0080 }
0081 HGCalGeomParameters::layerParameters laypar(rin, rout, zz);
0082 layers[lay] = laypar;
0083 }
0084 DD3Vector x, y, z;
0085 fv.rotation().GetComponents(x, y, z);
0086 const CLHEP::HepRep3x3 rotation(x.X(), y.X(), z.X(), x.Y(), y.Y(), z.Y(), x.Z(), y.Z(), z.Z());
0087 const CLHEP::HepRotation hr(rotation);
0088 double xx = HGCalParameters::k_ScaleFromDDD * fv.translation().X();
0089 if (std::abs(xx) < tolerance)
0090 xx = 0;
0091 double yy = HGCalParameters::k_ScaleFromDDD * fv.translation().Y();
0092 if (std::abs(yy) < tolerance)
0093 yy = 0;
0094 double zz = HGCalParameters::k_ScaleFromDDD * fv.translation().Z();
0095 const CLHEP::Hep3Vector h3v(xx, yy, zz);
0096 HGCalParameters::hgtrform mytrf;
0097 mytrf.zp = zp;
0098 mytrf.lay = lay;
0099 mytrf.sec = 0;
0100 mytrf.subsec = 0;
0101 mytrf.h3v = h3v;
0102 mytrf.hr = hr;
0103 trforms.emplace_back(mytrf);
0104 trformUse.emplace_back(false);
0105 }
0106 dodet = fv.next();
0107 }
0108
0109
0110
0111
0112
0113 std::unordered_map<int32_t, int32_t> copies;
0114 HGCalParameters::layer_map copiesInLayers(layers.size() + 1);
0115 std::vector<int32_t> wafer2copy;
0116 std::vector<HGCalGeomParameters::cellParameters> wafers;
0117 std::string attribute = "Volume";
0118 DDValue val1(attribute, sdTag2, 0.0);
0119 DDSpecificsMatchesValueFilter filter1{val1};
0120 DDFilteredView fv1(*cpv, filter1);
0121 bool ok = fv1.firstChild();
0122 if (!ok) {
0123 throw cms::Exception("DDException") << "Attribute " << val1 << " not found but needed.";
0124 } else {
0125 dodet = true;
0126 std::unordered_set<std::string> names;
0127 while (dodet) {
0128 const DDSolid& sol = fv1.logicalPart().solid();
0129 const std::string& name = fv1.logicalPart().name().name();
0130 std::vector<int> copy = fv1.copyNumbers();
0131 int nsiz = static_cast<int>(copy.size());
0132 int wafer = (nsiz > 0) ? copy[nsiz - 1] : 0;
0133 int layer = (nsiz > 1) ? copy[nsiz - 2] : 0;
0134 if (nsiz < 2) {
0135 throw cms::Exception("DDException") << "Funny wafer # " << wafer << " in " << nsiz << " components";
0136 } else if (layer > static_cast<int>(layers.size())) {
0137 edm::LogWarning("HGCalGeom") << "Funny wafer # " << wafer << " Layer " << layer << ":" << layers.size()
0138 << " among " << nsiz << " components";
0139 } else {
0140 auto itr = copies.find(wafer);
0141 auto cpy = copiesInLayers[layer].find(wafer);
0142 if (itr != copies.end() && cpy == copiesInLayers[layer].end()) {
0143 copiesInLayers[layer][wafer] = itr->second;
0144 }
0145 if (itr == copies.end()) {
0146 copies[wafer] = wafer2copy.size();
0147 copiesInLayers[layer][wafer] = wafer2copy.size();
0148 double xx = HGCalParameters::k_ScaleFromDDD * fv1.translation().X();
0149 if (std::abs(xx) < tolerance)
0150 xx = 0;
0151 double yy = HGCalParameters::k_ScaleFromDDD * fv1.translation().Y();
0152 if (std::abs(yy) < tolerance)
0153 yy = 0;
0154 wafer2copy.emplace_back(wafer);
0155 GlobalPoint p(xx, yy, HGCalParameters::k_ScaleFromDDD * fv1.translation().Z());
0156 HGCalGeomParameters::cellParameters cell(false, wafer, p);
0157 wafers.emplace_back(cell);
0158 if (names.count(name) == 0) {
0159 std::vector<double> zv, rv;
0160 if (mode == HGCalGeometryMode::Polyhedra) {
0161 const DDPolyhedra& polyhedra = static_cast<DDPolyhedra>(sol);
0162 zv = polyhedra.zVec();
0163 rv = polyhedra.rMaxVec();
0164 } else {
0165 const DDExtrudedPolygon& polygon = static_cast<DDExtrudedPolygon>(sol);
0166 zv = polygon.zVec();
0167 rv = polygon.xVec();
0168 }
0169 php.waferR_ = 2.0 * HGCalParameters::k_ScaleFromDDDToG4 * rv[0] * tan30deg_;
0170 php.waferSize_ = HGCalParameters::k_ScaleFromDDD * rv[0];
0171 double dz = 0.5 * HGCalParameters::k_ScaleFromDDDToG4 * (zv[1] - zv[0]);
0172 #ifdef EDM_ML_DEBUG
0173 edm::LogVerbatim("HGCalGeom")
0174 << "Mode " << mode << " R " << php.waferSize_ << ":" << php.waferR_ << " z " << dz;
0175 #endif
0176 HGCalParameters::hgtrap mytr;
0177 mytr.lay = 1;
0178 mytr.bl = php.waferR_;
0179 mytr.tl = php.waferR_;
0180 mytr.h = php.waferR_;
0181 mytr.dz = dz;
0182 mytr.alpha = 0.0;
0183 mytr.cellSize = waferSize_;
0184 php.fillModule(mytr, false);
0185 names.insert(name);
0186 }
0187 }
0188 }
0189 dodet = fv1.next();
0190 }
0191 }
0192
0193
0194 std::map<int, int> wafertype;
0195 std::map<int, HGCalGeomParameters::cellParameters> cellsf, cellsc;
0196 DDValue val2(attribute, sdTag3, 0.0);
0197 DDSpecificsMatchesValueFilter filter2{val2};
0198 DDFilteredView fv2(*cpv, filter2);
0199 ok = fv2.firstChild();
0200 if (!ok) {
0201 throw cms::Exception("DDException") << "Attribute " << val2 << " not found but needed.";
0202 } else {
0203 dodet = true;
0204 while (dodet) {
0205 const DDSolid& sol = fv2.logicalPart().solid();
0206 const std::string& name = sol.name().name();
0207 std::vector<int> copy = fv2.copyNumbers();
0208 int nsiz = static_cast<int>(copy.size());
0209 int cellx = (nsiz > 0) ? copy[nsiz - 1] : 0;
0210 int wafer = (nsiz > 1) ? copy[nsiz - 2] : 0;
0211 int cell = HGCalTypes::getUnpackedCell6(cellx);
0212 int type = HGCalTypes::getUnpackedCellType6(cellx);
0213 if (type != 1 && type != 2) {
0214 throw cms::Exception("DDException")
0215 << "Funny cell # " << cell << " type " << type << " in " << nsiz << " components";
0216 } else {
0217 auto ktr = wafertype.find(wafer);
0218 if (ktr == wafertype.end())
0219 wafertype[wafer] = type;
0220 bool newc(false);
0221 std::map<int, HGCalGeomParameters::cellParameters>::iterator itr;
0222 double cellsize = php.cellSize_[0];
0223 if (type == 1) {
0224 itr = cellsf.find(cell);
0225 newc = (itr == cellsf.end());
0226 } else {
0227 itr = cellsc.find(cell);
0228 newc = (itr == cellsc.end());
0229 cellsize = php.cellSize_[1];
0230 }
0231 if (newc) {
0232 bool half = (name.find("Half") != std::string::npos);
0233 double xx = HGCalParameters::k_ScaleFromDDD * fv2.translation().X();
0234 double yy = HGCalParameters::k_ScaleFromDDD * fv2.translation().Y();
0235 if (half) {
0236 math::XYZPointD p1(-2.0 * cellsize / 9.0, 0, 0);
0237 math::XYZPointD p2 = fv2.rotation()(p1);
0238 xx += (HGCalParameters::k_ScaleFromDDD * (p2.X()));
0239 yy += (HGCalParameters::k_ScaleFromDDD * (p2.Y()));
0240 #ifdef EDM_ML_DEBUG
0241 if (std::abs(p2.X()) < HGCalParameters::tol)
0242 p2.SetX(0.0);
0243 if (std::abs(p2.Z()) < HGCalParameters::tol)
0244 p2.SetZ(0.0);
0245 edm::LogVerbatim("HGCalGeom") << "Wafer " << wafer << " Type " << type << " Cell " << cellx << " local "
0246 << xx << ":" << yy << " new " << p1 << ":" << p2;
0247 #endif
0248 }
0249 HGCalGeomParameters::cellParameters cp(half, wafer, GlobalPoint(xx, yy, 0));
0250 if (type == 1) {
0251 cellsf[cell] = cp;
0252 } else {
0253 cellsc[cell] = cp;
0254 }
0255 }
0256 }
0257 dodet = fv2.next();
0258 }
0259 }
0260
0261 loadGeometryHexagon(
0262 layers, trforms, trformUse, copies, copiesInLayers, wafer2copy, wafers, wafertype, cellsf, cellsc, php);
0263 }
0264
0265 void HGCalGeomParameters::loadGeometryHexagon(const cms::DDCompactView* cpv,
0266 HGCalParameters& php,
0267 const std::string& sdTag1,
0268 const std::string& sdTag2,
0269 const std::string& sdTag3,
0270 HGCalGeometryMode::WaferMode mode) {
0271 const cms::DDFilter filter("Volume", sdTag1);
0272 cms::DDFilteredView fv((*cpv), filter);
0273 std::map<int, HGCalGeomParameters::layerParameters> layers;
0274 std::vector<HGCalParameters::hgtrform> trforms;
0275 std::vector<bool> trformUse;
0276 std::vector<std::pair<int, int> > trused;
0277
0278 while (fv.firstChild()) {
0279 const std::vector<double>& pars = fv.parameters();
0280
0281 std::vector<int> copy = fv.copyNos();
0282 int nsiz = static_cast<int>(copy.size());
0283 int lay = (nsiz > 0) ? copy[0] : 0;
0284 int zp = (nsiz > 2) ? copy[2] : -1;
0285 if (zp != 1)
0286 zp = -1;
0287 if (lay == 0) {
0288 throw cms::Exception("DDException") << "Funny layer # " << lay << " zp " << zp << " in " << nsiz << " components";
0289 } else {
0290 if (std::find(php.layer_.begin(), php.layer_.end(), lay) == php.layer_.end())
0291 php.layer_.emplace_back(lay);
0292 auto itr = layers.find(lay);
0293 double zz = HGCalParameters::k_ScaleFromDD4hep * fv.translation().Z();
0294 if (itr == layers.end()) {
0295 double rin(0), rout(0);
0296 if (dd4hep::isA<dd4hep::Polyhedra>(fv.solid())) {
0297 rin = 0.5 * HGCalParameters::k_ScaleFromDD4hep * (pars[5] + pars[8]);
0298 rout = 0.5 * HGCalParameters::k_ScaleFromDD4hep * (pars[6] + pars[9]);
0299 } else if (dd4hep::isA<dd4hep::Tube>(fv.solid())) {
0300 dd4hep::Tube tubeSeg(fv.solid());
0301 rin = HGCalParameters::k_ScaleFromDD4hep * tubeSeg.rMin();
0302 rout = HGCalParameters::k_ScaleFromDD4hep * tubeSeg.rMax();
0303 }
0304 HGCalGeomParameters::layerParameters laypar(rin, rout, zz);
0305 layers[lay] = laypar;
0306 }
0307 std::pair<int, int> layz(lay, zp);
0308 if (std::find(trused.begin(), trused.end(), layz) == trused.end()) {
0309 trused.emplace_back(layz);
0310 DD3Vector x, y, z;
0311 fv.rotation().GetComponents(x, y, z);
0312 const CLHEP::HepRep3x3 rotation(x.X(), y.X(), z.X(), x.Y(), y.Y(), z.Y(), x.Z(), y.Z(), z.Z());
0313 const CLHEP::HepRotation hr(rotation);
0314 double xx = HGCalParameters::k_ScaleFromDD4hep * fv.translation().X();
0315 if (std::abs(xx) < tolerance)
0316 xx = 0;
0317 double yy = HGCalParameters::k_ScaleFromDD4hep * fv.translation().Y();
0318 if (std::abs(yy) < tolerance)
0319 yy = 0;
0320 double zz = HGCalParameters::k_ScaleFromDD4hep * fv.translation().Z();
0321 const CLHEP::Hep3Vector h3v(xx, yy, zz);
0322 HGCalParameters::hgtrform mytrf;
0323 mytrf.zp = zp;
0324 mytrf.lay = lay;
0325 mytrf.sec = 0;
0326 mytrf.subsec = 0;
0327 mytrf.h3v = h3v;
0328 mytrf.hr = hr;
0329 trforms.emplace_back(mytrf);
0330 trformUse.emplace_back(false);
0331 }
0332 }
0333 }
0334
0335
0336
0337
0338
0339 std::unordered_map<int32_t, int32_t> copies;
0340 HGCalParameters::layer_map copiesInLayers(layers.size() + 1);
0341 std::vector<int32_t> wafer2copy;
0342 std::vector<HGCalGeomParameters::cellParameters> wafers;
0343 const cms::DDFilter filter1("Volume", sdTag2);
0344 cms::DDFilteredView fv1((*cpv), filter1);
0345 bool ok = fv1.firstChild();
0346 if (!ok) {
0347 throw cms::Exception("DDException") << "Attribute " << sdTag2 << " not found but needed.";
0348 } else {
0349 bool dodet = true;
0350 std::unordered_set<std::string> names;
0351 while (dodet) {
0352 const std::string name = static_cast<std::string>(fv1.name());
0353 std::vector<int> copy = fv1.copyNos();
0354 int nsiz = static_cast<int>(copy.size());
0355 int wafer = (nsiz > 0) ? copy[0] : 0;
0356 int layer = (nsiz > 1) ? copy[1] : 0;
0357 if (nsiz < 2) {
0358 throw cms::Exception("DDException") << "Funny wafer # " << wafer << " in " << nsiz << " components";
0359 } else if (layer > static_cast<int>(layers.size())) {
0360 edm::LogWarning("HGCalGeom") << "Funny wafer # " << wafer << " Layer " << layer << ":" << layers.size()
0361 << " among " << nsiz << " components";
0362 } else {
0363 auto itr = copies.find(wafer);
0364 auto cpy = copiesInLayers[layer].find(wafer);
0365 if (itr != copies.end() && cpy == copiesInLayers[layer].end()) {
0366 copiesInLayers[layer][wafer] = itr->second;
0367 }
0368 if (itr == copies.end()) {
0369 copies[wafer] = wafer2copy.size();
0370 copiesInLayers[layer][wafer] = wafer2copy.size();
0371 double xx = HGCalParameters::k_ScaleFromDD4hep * fv1.translation().X();
0372 if (std::abs(xx) < tolerance)
0373 xx = 0;
0374 double yy = HGCalParameters::k_ScaleFromDD4hep * fv1.translation().Y();
0375 if (std::abs(yy) < tolerance)
0376 yy = 0;
0377 wafer2copy.emplace_back(wafer);
0378 GlobalPoint p(xx, yy, HGCalParameters::k_ScaleFromDD4hep * fv1.translation().Z());
0379 HGCalGeomParameters::cellParameters cell(false, wafer, p);
0380 wafers.emplace_back(cell);
0381 if (names.count(name) == 0) {
0382 double zv[2], rv;
0383 const std::vector<double>& pars = fv1.parameters();
0384 if (mode == HGCalGeometryMode::Polyhedra) {
0385 zv[0] = pars[4];
0386 zv[1] = pars[7];
0387 rv = pars[6];
0388 } else {
0389 zv[0] = pars[3];
0390 zv[1] = pars[9];
0391 rv = pars[4];
0392 }
0393 php.waferR_ = 2.0 * HGCalParameters::k_ScaleFromDD4hepToG4 * rv * tan30deg_;
0394 php.waferSize_ = HGCalParameters::k_ScaleFromDD4hep * rv;
0395 double dz = 0.5 * HGCalParameters::k_ScaleFromDD4hepToG4 * (zv[1] - zv[0]);
0396 #ifdef EDM_ML_DEBUG
0397 edm::LogVerbatim("HGCalGeom")
0398 << "Mode " << mode << " R " << php.waferSize_ << ":" << php.waferR_ << " z " << dz;
0399 #endif
0400 HGCalParameters::hgtrap mytr;
0401 mytr.lay = 1;
0402 mytr.bl = php.waferR_;
0403 mytr.tl = php.waferR_;
0404 mytr.h = php.waferR_;
0405 mytr.dz = dz;
0406 mytr.alpha = 0.0;
0407 mytr.cellSize = waferSize_;
0408 php.fillModule(mytr, false);
0409 names.insert(name);
0410 }
0411 }
0412 }
0413 dodet = fv1.firstChild();
0414 }
0415 }
0416
0417
0418 std::map<int, int> wafertype;
0419 std::map<int, HGCalGeomParameters::cellParameters> cellsf, cellsc;
0420 const cms::DDFilter filter2("Volume", sdTag3);
0421 cms::DDFilteredView fv2((*cpv), filter2);
0422 ok = fv2.firstChild();
0423 if (!ok) {
0424 throw cms::Exception("DDException") << "Attribute " << sdTag3 << " not found but needed.";
0425 } else {
0426 bool dodet = true;
0427 while (dodet) {
0428 const std::string name = static_cast<std::string>(fv2.name());
0429 std::vector<int> copy = fv2.copyNos();
0430 int nsiz = static_cast<int>(copy.size());
0431 int cellx = (nsiz > 0) ? copy[0] : 0;
0432 int wafer = (nsiz > 1) ? copy[1] : 0;
0433 int cell = HGCalTypes::getUnpackedCell6(cellx);
0434 int type = HGCalTypes::getUnpackedCellType6(cellx);
0435 if (type != 1 && type != 2) {
0436 throw cms::Exception("DDException")
0437 << "Funny cell # " << cell << " type " << type << " in " << nsiz << " components";
0438 } else {
0439 auto ktr = wafertype.find(wafer);
0440 if (ktr == wafertype.end())
0441 wafertype[wafer] = type;
0442 bool newc(false);
0443 std::map<int, HGCalGeomParameters::cellParameters>::iterator itr;
0444 double cellsize = php.cellSize_[0];
0445 if (type == 1) {
0446 itr = cellsf.find(cell);
0447 newc = (itr == cellsf.end());
0448 } else {
0449 itr = cellsc.find(cell);
0450 newc = (itr == cellsc.end());
0451 cellsize = php.cellSize_[1];
0452 }
0453 if (newc) {
0454 bool half = (name.find("Half") != std::string::npos);
0455 double xx = HGCalParameters::k_ScaleFromDD4hep * fv2.translation().X();
0456 double yy = HGCalParameters::k_ScaleFromDD4hep * fv2.translation().Y();
0457 if (half) {
0458 math::XYZPointD p1(-2.0 * cellsize / 9.0, 0, 0);
0459 math::XYZPointD p2 = fv2.rotation()(p1);
0460 xx += (HGCalParameters::k_ScaleFromDDD * (p2.X()));
0461 yy += (HGCalParameters::k_ScaleFromDDD * (p2.Y()));
0462 #ifdef EDM_ML_DEBUG
0463 if (std::abs(p2.X()) < HGCalParameters::tol)
0464 p2.SetX(0.0);
0465 if (std::abs(p2.Z()) < HGCalParameters::tol)
0466 p2.SetZ(0.0);
0467 edm::LogVerbatim("HGCalGeom") << "Wafer " << wafer << " Type " << type << " Cell " << cellx << " local "
0468 << xx << ":" << yy << " new " << p1 << ":" << p2;
0469 #endif
0470 }
0471 HGCalGeomParameters::cellParameters cp(half, wafer, GlobalPoint(xx, yy, 0));
0472 if (type == 1) {
0473 cellsf[cell] = cp;
0474 } else {
0475 cellsc[cell] = cp;
0476 }
0477 }
0478 }
0479 dodet = fv2.firstChild();
0480 }
0481 }
0482
0483 loadGeometryHexagon(
0484 layers, trforms, trformUse, copies, copiesInLayers, wafer2copy, wafers, wafertype, cellsf, cellsc, php);
0485 }
0486
0487 void HGCalGeomParameters::loadGeometryHexagon(const std::map<int, HGCalGeomParameters::layerParameters>& layers,
0488 std::vector<HGCalParameters::hgtrform>& trforms,
0489 std::vector<bool>& trformUse,
0490 const std::unordered_map<int32_t, int32_t>& copies,
0491 const HGCalParameters::layer_map& copiesInLayers,
0492 const std::vector<int32_t>& wafer2copy,
0493 const std::vector<HGCalGeomParameters::cellParameters>& wafers,
0494 const std::map<int, int>& wafertype,
0495 const std::map<int, HGCalGeomParameters::cellParameters>& cellsf,
0496 const std::map<int, HGCalGeomParameters::cellParameters>& cellsc,
0497 HGCalParameters& php) {
0498 if (((cellsf.size() + cellsc.size()) == 0) || (wafers.empty()) || (layers.empty())) {
0499 throw cms::Exception("DDException") << "HGCalGeomParameters: mismatch between geometry and specpar: cells "
0500 << cellsf.size() << ":" << cellsc.size() << " wafers " << wafers.size()
0501 << " layers " << layers.size();
0502 }
0503
0504 for (unsigned int i = 0; i < layers.size(); ++i) {
0505 for (auto& layer : layers) {
0506 if (layer.first == static_cast<int>(i + php.firstLayer_)) {
0507 php.layerIndex_.emplace_back(i);
0508 php.rMinLayHex_.emplace_back(layer.second.rmin);
0509 php.rMaxLayHex_.emplace_back(layer.second.rmax);
0510 php.zLayerHex_.emplace_back(layer.second.zpos);
0511 break;
0512 }
0513 }
0514 }
0515
0516 for (unsigned int i = 0; i < php.layer_.size(); ++i) {
0517 for (unsigned int i1 = 0; i1 < trforms.size(); ++i1) {
0518 if (!trformUse[i1] && php.layerGroup_[trforms[i1].lay - 1] == static_cast<int>(i + 1)) {
0519 trforms[i1].h3v *= static_cast<double>(HGCalParameters::k_ScaleFromDDD);
0520 trforms[i1].lay = (i + 1);
0521 trformUse[i1] = true;
0522 php.fillTrForm(trforms[i1]);
0523 int nz(1);
0524 for (unsigned int i2 = i1 + 1; i2 < trforms.size(); ++i2) {
0525 if (!trformUse[i2] && trforms[i2].zp == trforms[i1].zp &&
0526 php.layerGroup_[trforms[i2].lay - 1] == static_cast<int>(i + 1)) {
0527 php.addTrForm(trforms[i2].h3v);
0528 nz++;
0529 trformUse[i2] = true;
0530 }
0531 }
0532 if (nz > 0) {
0533 php.scaleTrForm(double(1.0 / nz));
0534 }
0535 }
0536 }
0537 }
0538
0539 double rmin = HGCalParameters::k_ScaleFromDDD * php.waferR_;
0540 for (unsigned i = 0; i < wafer2copy.size(); ++i) {
0541 php.waferCopy_.emplace_back(wafer2copy[i]);
0542 php.waferPosX_.emplace_back(wafers[i].xyz.x());
0543 php.waferPosY_.emplace_back(wafers[i].xyz.y());
0544 auto ktr = wafertype.find(wafer2copy[i]);
0545 int typet = (ktr == wafertype.end()) ? 0 : (ktr->second);
0546 php.waferTypeT_.emplace_back(typet);
0547 double r = wafers[i].xyz.perp();
0548 int type(3);
0549 for (int k = 1; k < 4; ++k) {
0550 if ((r + rmin) <= php.boundR_[k]) {
0551 type = k;
0552 break;
0553 }
0554 }
0555 php.waferTypeL_.emplace_back(type);
0556 }
0557 php.copiesInLayers_ = copiesInLayers;
0558 php.nSectors_ = static_cast<int>(php.waferCopy_.size());
0559
0560 std::vector<HGCalGeomParameters::cellParameters>::const_iterator itrf = wafers.end();
0561 for (unsigned int i = 0; i < cellsf.size(); ++i) {
0562 auto itr = cellsf.find(i);
0563 if (itr == cellsf.end()) {
0564 throw cms::Exception("DDException") << "HGCalGeomParameters: missing info for fine cell number " << i;
0565 } else {
0566 double xx = (itr->second).xyz.x();
0567 double yy = (itr->second).xyz.y();
0568 int waf = (itr->second).wafer;
0569 std::pair<double, double> xy = cellPosition(wafers, itrf, waf, xx, yy);
0570 php.cellFineX_.emplace_back(xy.first);
0571 php.cellFineY_.emplace_back(xy.second);
0572 php.cellFineHalf_.emplace_back((itr->second).half);
0573 }
0574 }
0575 itrf = wafers.end();
0576 for (unsigned int i = 0; i < cellsc.size(); ++i) {
0577 auto itr = cellsc.find(i);
0578 if (itr == cellsc.end()) {
0579 throw cms::Exception("DDException") << "HGCalGeomParameters: missing info for coarse cell number " << i;
0580 } else {
0581 double xx = (itr->second).xyz.x();
0582 double yy = (itr->second).xyz.y();
0583 int waf = (itr->second).wafer;
0584 std::pair<double, double> xy = cellPosition(wafers, itrf, waf, xx, yy);
0585 php.cellCoarseX_.emplace_back(xy.first);
0586 php.cellCoarseY_.emplace_back(xy.second);
0587 php.cellCoarseHalf_.emplace_back((itr->second).half);
0588 }
0589 }
0590 int depth(0);
0591 for (unsigned int i = 0; i < php.layerGroup_.size(); ++i) {
0592 bool first(true);
0593 for (unsigned int k = 0; k < php.layerGroup_.size(); ++k) {
0594 if (php.layerGroup_[k] == static_cast<int>(i + 1)) {
0595 if (first) {
0596 php.depth_.emplace_back(i + 1);
0597 php.depthIndex_.emplace_back(depth);
0598 php.depthLayerF_.emplace_back(k);
0599 ++depth;
0600 first = false;
0601 }
0602 }
0603 }
0604 }
0605 HGCalParameters::hgtrap mytr = php.getModule(0, false);
0606 mytr.bl *= HGCalParameters::k_ScaleFromDDD;
0607 mytr.tl *= HGCalParameters::k_ScaleFromDDD;
0608 mytr.h *= HGCalParameters::k_ScaleFromDDD;
0609 mytr.dz *= HGCalParameters::k_ScaleFromDDD;
0610 mytr.cellSize *= HGCalParameters::k_ScaleFromDDD;
0611 double dz = mytr.dz;
0612 php.fillModule(mytr, true);
0613 mytr.dz = 2 * dz;
0614 php.fillModule(mytr, true);
0615 mytr.dz = 3 * dz;
0616 php.fillModule(mytr, true);
0617 #ifdef EDM_ML_DEBUG
0618 edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters finds " << php.zLayerHex_.size() << " layers";
0619 for (unsigned int i = 0; i < php.zLayerHex_.size(); ++i) {
0620 int k = php.layerIndex_[i];
0621 edm::LogVerbatim("HGCalGeom") << "Layer[" << i << ":" << k << ":" << php.layer_[k]
0622 << "] with r = " << php.rMinLayHex_[i] << ":" << php.rMaxLayHex_[i]
0623 << " at z = " << php.zLayerHex_[i];
0624 }
0625 edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters has " << php.depthIndex_.size() << " depths";
0626 for (unsigned int i = 0; i < php.depthIndex_.size(); ++i) {
0627 int k = php.depthIndex_[i];
0628 edm::LogVerbatim("HGCalGeom") << "Reco Layer[" << i << ":" << k << "] First Layer " << php.depthLayerF_[i]
0629 << " Depth " << php.depth_[k];
0630 }
0631 edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters finds " << php.nSectors_ << " wafers";
0632 for (unsigned int i = 0; i < php.waferCopy_.size(); ++i)
0633 edm::LogVerbatim("HGCalGeom") << "Wafer[" << i << ": " << php.waferCopy_[i] << "] type " << php.waferTypeL_[i]
0634 << ":" << php.waferTypeT_[i] << " at (" << php.waferPosX_[i] << ","
0635 << php.waferPosY_[i] << ",0)";
0636 edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: wafer radius " << php.waferR_ << " and dimensions of the "
0637 << "wafers:";
0638 edm::LogVerbatim("HGCalGeom") << "Sim[0] " << php.moduleLayS_[0] << " dx " << php.moduleBlS_[0] << ":"
0639 << php.moduleTlS_[0] << " dy " << php.moduleHS_[0] << " dz " << php.moduleDzS_[0]
0640 << " alpha " << php.moduleAlphaS_[0];
0641 for (unsigned int k = 0; k < php.moduleLayR_.size(); ++k)
0642 edm::LogVerbatim("HGCalGeom") << "Rec[" << k << "] " << php.moduleLayR_[k] << " dx " << php.moduleBlR_[k] << ":"
0643 << php.moduleTlR_[k] << " dy " << php.moduleHR_[k] << " dz " << php.moduleDzR_[k]
0644 << " alpha " << php.moduleAlphaR_[k];
0645 edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters finds " << php.cellFineX_.size() << " fine cells in a wafer";
0646 for (unsigned int i = 0; i < php.cellFineX_.size(); ++i)
0647 edm::LogVerbatim("HGCalGeom") << "Fine Cell[" << i << "] at (" << php.cellFineX_[i] << "," << php.cellFineY_[i]
0648 << ",0)";
0649 edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters finds " << php.cellCoarseX_.size()
0650 << " coarse cells in a wafer";
0651 for (unsigned int i = 0; i < php.cellCoarseX_.size(); ++i)
0652 edm::LogVerbatim("HGCalGeom") << "Coarse Cell[" << i << "] at (" << php.cellCoarseX_[i] << ","
0653 << php.cellCoarseY_[i] << ",0)";
0654 edm::LogVerbatim("HGCalGeom") << "Obtained " << php.trformIndex_.size() << " transformation matrices";
0655 for (unsigned int k = 0; k < php.trformIndex_.size(); ++k) {
0656 edm::LogVerbatim("HGCalGeom") << "Matrix[" << k << "] (" << std::hex << php.trformIndex_[k] << std::dec
0657 << ") Translation (" << php.trformTranX_[k] << ", " << php.trformTranY_[k] << ", "
0658 << php.trformTranZ_[k] << " Rotation (" << php.trformRotXX_[k] << ", "
0659 << php.trformRotYX_[k] << ", " << php.trformRotZX_[k] << ", " << php.trformRotXY_[k]
0660 << ", " << php.trformRotYY_[k] << ", " << php.trformRotZY_[k] << ", "
0661 << php.trformRotXZ_[k] << ", " << php.trformRotYZ_[k] << ", " << php.trformRotZZ_[k]
0662 << ")";
0663 }
0664 edm::LogVerbatim("HGCalGeom") << "Dump copiesInLayers for " << php.copiesInLayers_.size() << " layers";
0665 for (unsigned int k = 0; k < php.copiesInLayers_.size(); ++k) {
0666 const auto& theModules = php.copiesInLayers_[k];
0667 edm::LogVerbatim("HGCalGeom") << "Layer " << k << ":" << theModules.size();
0668 int k2(0);
0669 for (std::unordered_map<int, int>::const_iterator itr = theModules.begin(); itr != theModules.end(); ++itr, ++k2) {
0670 edm::LogVerbatim("HGCalGeom") << "[" << k2 << "] " << itr->first << ":" << itr->second;
0671 }
0672 }
0673 #endif
0674 }
0675
0676 void HGCalGeomParameters::loadGeometryHexagon8(const DDFilteredView& _fv, HGCalParameters& php, int firstLayer) {
0677 DDFilteredView fv = _fv;
0678 bool dodet(true);
0679 std::map<int, HGCalGeomParameters::layerParameters> layers;
0680 std::map<std::pair<int, int>, HGCalParameters::hgtrform> trforms;
0681 int levelTop = 3 + std::max(php.levelT_[0], php.levelT_[1]);
0682 #ifdef EDM_ML_DEBUG
0683 int ntot(0);
0684 #endif
0685 while (dodet) {
0686 #ifdef EDM_ML_DEBUG
0687 ++ntot;
0688 #endif
0689 std::vector<int> copy = fv.copyNumbers();
0690 int nsiz = static_cast<int>(copy.size());
0691 if (nsiz < levelTop) {
0692 int lay = copy[nsiz - 1];
0693 int zside = (nsiz > php.levelZSide_) ? copy[php.levelZSide_] : -1;
0694 if (zside != 1)
0695 zside = -1;
0696 const DDSolid& sol = fv.logicalPart().solid();
0697 #ifdef EDM_ML_DEBUG
0698 edm::LogVerbatim("HGCalGeom") << sol.name() << " shape " << sol.shape() << " size " << nsiz << ":" << levelTop
0699 << " lay " << lay << " z " << zside;
0700 #endif
0701 if (lay == 0) {
0702 throw cms::Exception("DDException")
0703 << "Funny layer # " << lay << " zp " << zside << " in " << nsiz << " components";
0704 } else if (sol.shape() == DDSolidShape::ddtubs) {
0705 if (std::find(php.layer_.begin(), php.layer_.end(), lay) == php.layer_.end())
0706 php.layer_.emplace_back(lay);
0707 const DDTubs& tube = static_cast<DDTubs>(sol);
0708 double rin = HGCalParameters::k_ScaleFromDDD * tube.rIn();
0709 double rout = HGCalParameters::k_ScaleFromDDD * tube.rOut();
0710 auto itr = layers.find(lay);
0711 if (itr == layers.end()) {
0712 double zp = HGCalParameters::k_ScaleFromDDD * fv.translation().Z();
0713 HGCalGeomParameters::layerParameters laypar(rin, rout, zp);
0714 layers[lay] = laypar;
0715 } else {
0716 (itr->second).rmin = std::min(rin, (itr->second).rmin);
0717 (itr->second).rmax = std::max(rout, (itr->second).rmax);
0718 }
0719 if (trforms.find(std::make_pair(lay, zside)) == trforms.end()) {
0720 DD3Vector x, y, z;
0721 fv.rotation().GetComponents(x, y, z);
0722 const CLHEP::HepRep3x3 rotation(x.X(), y.X(), z.X(), x.Y(), y.Y(), z.Y(), x.Z(), y.Z(), z.Z());
0723 const CLHEP::HepRotation hr(rotation);
0724 double xx =
0725 ((std::abs(fv.translation().X()) < tolerance) ? 0
0726 : HGCalParameters::k_ScaleFromDDD * fv.translation().X());
0727 double yy =
0728 ((std::abs(fv.translation().Y()) < tolerance) ? 0
0729 : HGCalParameters::k_ScaleFromDDD * fv.translation().Y());
0730 const CLHEP::Hep3Vector h3v(xx, yy, HGCalParameters::k_ScaleFromDDD * fv.translation().Z());
0731 HGCalParameters::hgtrform mytrf;
0732 mytrf.zp = zside;
0733 mytrf.lay = lay;
0734 mytrf.sec = 0;
0735 mytrf.subsec = 0;
0736 mytrf.h3v = h3v;
0737 mytrf.hr = hr;
0738 trforms[std::make_pair(lay, zside)] = mytrf;
0739 }
0740 }
0741 }
0742 dodet = fv.next();
0743 }
0744 #ifdef EDM_ML_DEBUG
0745 edm::LogVerbatim("HGCalGeom") << "Total # of views " << ntot;
0746 #endif
0747 loadGeometryHexagon8(layers, trforms, firstLayer, php);
0748 }
0749
0750 void HGCalGeomParameters::loadGeometryHexagon8(const cms::DDCompactView* cpv,
0751 HGCalParameters& php,
0752 const std::string& sdTag1,
0753 int firstLayer) {
0754 const cms::DDFilter filter("Volume", sdTag1);
0755 cms::DDFilteredView fv((*cpv), filter);
0756 std::map<int, HGCalGeomParameters::layerParameters> layers;
0757 std::map<std::pair<int, int>, HGCalParameters::hgtrform> trforms;
0758 int levelTop = 3 + std::max(php.levelT_[0], php.levelT_[1]);
0759 #ifdef EDM_ML_DEBUG
0760 int ntot(0);
0761 #endif
0762 while (fv.firstChild()) {
0763 #ifdef EDM_ML_DEBUG
0764 ++ntot;
0765 #endif
0766
0767 int nsiz = static_cast<int>(fv.level());
0768 if (nsiz < levelTop) {
0769 std::vector<int> copy = fv.copyNos();
0770 int lay = copy[0];
0771 int zside = (nsiz > php.levelZSide_) ? copy[nsiz - php.levelZSide_ - 1] : -1;
0772 if (zside != 1)
0773 zside = -1;
0774 #ifdef EDM_ML_DEBUG
0775 edm::LogVerbatim("HGCalGeom") << fv.name() << " shape " << cms::dd::name(cms::DDSolidShapeMap, fv.shape())
0776 << " size " << nsiz << ":" << levelTop << " lay " << lay << " z " << zside << ":"
0777 << php.levelZSide_;
0778 #endif
0779 if (lay == 0) {
0780 throw cms::Exception("DDException")
0781 << "Funny layer # " << lay << " zp " << zside << " in " << nsiz << " components";
0782 } else if (fv.shape() == cms::DDSolidShape::ddtubs) {
0783 if (std::find(php.layer_.begin(), php.layer_.end(), lay) == php.layer_.end())
0784 php.layer_.emplace_back(lay);
0785 const std::vector<double>& pars = fv.parameters();
0786 double rin = HGCalParameters::k_ScaleFromDD4hep * pars[0];
0787 double rout = HGCalParameters::k_ScaleFromDD4hep * pars[1];
0788 auto itr = layers.find(lay);
0789 if (itr == layers.end()) {
0790 double zp = HGCalParameters::k_ScaleFromDD4hep * fv.translation().Z();
0791 HGCalGeomParameters::layerParameters laypar(rin, rout, zp);
0792 layers[lay] = laypar;
0793 } else {
0794 (itr->second).rmin = std::min(rin, (itr->second).rmin);
0795 (itr->second).rmax = std::max(rout, (itr->second).rmax);
0796 }
0797 if (trforms.find(std::make_pair(lay, zside)) == trforms.end()) {
0798 DD3Vector x, y, z;
0799 fv.rotation().GetComponents(x, y, z);
0800 const CLHEP::HepRep3x3 rotation(x.X(), y.X(), z.X(), x.Y(), y.Y(), z.Y(), x.Z(), y.Z(), z.Z());
0801 const CLHEP::HepRotation hr(rotation);
0802 double xx = ((std::abs(fv.translation().X()) < tolerance)
0803 ? 0
0804 : HGCalParameters::k_ScaleFromDD4hep * fv.translation().X());
0805 double yy = ((std::abs(fv.translation().Y()) < tolerance)
0806 ? 0
0807 : HGCalParameters::k_ScaleFromDD4hep * fv.translation().Y());
0808 const CLHEP::Hep3Vector h3v(xx, yy, HGCalParameters::k_ScaleFromDD4hep * fv.translation().Z());
0809 HGCalParameters::hgtrform mytrf;
0810 mytrf.zp = zside;
0811 mytrf.lay = lay;
0812 mytrf.sec = 0;
0813 mytrf.subsec = 0;
0814 mytrf.h3v = h3v;
0815 mytrf.hr = hr;
0816 trforms[std::make_pair(lay, zside)] = mytrf;
0817 }
0818 }
0819 }
0820 }
0821 #ifdef EDM_ML_DEBUG
0822 edm::LogVerbatim("HGCalGeom") << "Total # of views " << ntot;
0823 #endif
0824 loadGeometryHexagon8(layers, trforms, firstLayer, php);
0825 }
0826
0827 void HGCalGeomParameters::loadGeometryHexagonModule(const DDCompactView* cpv,
0828 HGCalParameters& php,
0829 const std::string& sdTag1,
0830 const std::string& sdTag2,
0831 int firstLayer) {
0832 #ifdef EDM_ML_DEBUG
0833 edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters (DDD)::loadGeometryHexagonModule called with tags " << sdTag1
0834 << ":" << sdTag2 << " firstLayer " << firstLayer << ":" << php.firstMixedLayer_;
0835 int ntot1(0), ntot2(0);
0836 #endif
0837 std::map<int, HGCalGeomParameters::layerParameters> layers;
0838 std::map<std::pair<int, int>, double> zvals;
0839 std::map<std::pair<int, int>, HGCalParameters::hgtrform> trforms;
0840 int levelTop = php.levelT_[0];
0841
0842 std::string attribute = "Volume";
0843 DDValue val1(attribute, sdTag2, 0.0);
0844 DDSpecificsMatchesValueFilter filter1{val1};
0845 DDFilteredView fv1(*cpv, filter1);
0846 bool dodet = fv1.firstChild();
0847 while (dodet) {
0848 #ifdef EDM_ML_DEBUG
0849 ++ntot1;
0850 #endif
0851 std::vector<int> copy = fv1.copyNumbers();
0852 int nsiz = static_cast<int>(copy.size());
0853 if (levelTop < nsiz) {
0854 int lay = copy[levelTop];
0855 int zside = (nsiz > php.levelZSide_) ? copy[php.levelZSide_] : -1;
0856 if (zside != 1)
0857 zside = -1;
0858 if (lay == 0) {
0859 throw cms::Exception("DDException")
0860 << "Funny layer # " << lay << " zp " << zside << " in " << nsiz << " components";
0861 } else {
0862 if (zvals.find(std::make_pair(lay, zside)) == zvals.end()) {
0863 zvals[std::make_pair(lay, zside)] = HGCalParameters::k_ScaleFromDDD * fv1.translation().Z();
0864 #ifdef EDM_ML_DEBUG
0865 std::ostringstream st1;
0866 st1 << "Name0 " << fv1.name() << " LTop " << levelTop << ":" << lay << " ZSide " << zside << " # of levels "
0867 << nsiz;
0868 for (const auto& c : copy)
0869 st1 << ":" << c;
0870 st1 << " Z " << zvals[std::make_pair(lay, zside)];
0871 edm::LogVerbatim("HGCalGeom") << st1.str();
0872 #endif
0873 }
0874 }
0875 }
0876 dodet = fv1.next();
0877 }
0878
0879 DDValue val2(attribute, sdTag1, 0.0);
0880 DDSpecificsMatchesValueFilter filter2{val2};
0881 DDFilteredView fv2(*cpv, filter2);
0882 dodet = fv2.firstChild();
0883 while (dodet) {
0884 std::vector<int> copy = fv2.copyNumbers();
0885 int nsiz = static_cast<int>(copy.size());
0886 #ifdef EDM_ML_DEBUG
0887 ++ntot2;
0888 edm::LogVerbatim("HGCalGeom") << "loadGeometryHexagonModule:: nsiz " << nsiz << " Ltop " << levelTop;
0889 #endif
0890 if (levelTop < nsiz) {
0891 int lay = copy[levelTop];
0892 int zside = (nsiz > php.levelZSide_) ? copy[php.levelZSide_] : -1;
0893 if (zside != 1)
0894 zside = -1;
0895 const DDSolid& sol = fv2.logicalPart().solid();
0896 #ifdef EDM_ML_DEBUG
0897 std::ostringstream st2;
0898 st2 << "Name1 " << sol.name() << " shape " << sol.shape() << " LTop " << levelTop << ":" << lay << " ZSide "
0899 << zside << ":" << php.levelZSide_ << " # of levels " << nsiz;
0900 for (const auto& c : copy)
0901 st2 << ":" << c;
0902 edm::LogVerbatim("HGCalGeom") << st2.str();
0903 #endif
0904 if (lay == 0) {
0905 throw cms::Exception("DDException")
0906 << "Funny layer # " << lay << " zp " << zside << " in " << nsiz << " components";
0907 } else if ((sol.shape() == DDSolidShape::ddtubs) || (sol.shape() == DDSolidShape::ddbox)) {
0908 if (zvals.find(std::make_pair(lay, zside)) != zvals.end()) {
0909 if (std::find(php.layer_.begin(), php.layer_.end(), lay) == php.layer_.end())
0910 php.layer_.emplace_back(lay);
0911 auto itr = layers.find(lay);
0912 if (itr == layers.end()) {
0913 double rin(0), rout(0);
0914 if (sol.shape() == DDSolidShape::ddtubs) {
0915 const DDTubs& tube = static_cast<DDTubs>(sol);
0916 rin = HGCalParameters::k_ScaleFromDDD * tube.rIn();
0917 rout = (php.firstMixedLayer_ > 0 && lay >= php.firstMixedLayer_)
0918 ? php.radiusMixBoundary_[lay - php.firstMixedLayer_]
0919 : HGCalParameters::k_ScaleFromDDD * tube.rOut();
0920 } else {
0921 const DDBox& box = static_cast<DDBox>(sol);
0922 rout = HGCalParameters::k_ScaleFromDDD * box.halfX();
0923 }
0924 double zp = zvals[std::make_pair(lay, 1)];
0925 HGCalGeomParameters::layerParameters laypar(rin, rout, zp);
0926 layers[lay] = laypar;
0927 #ifdef EDM_ML_DEBUG
0928 std::ostringstream st3;
0929 st3 << "Name1 " << fv2.name() << " LTop " << levelTop << ":" << lay << " ZSide " << zside << " # of levels "
0930 << nsiz;
0931 for (const auto& c : copy)
0932 st3 << ":" << c;
0933 st3 << " R " << rin << ":" << rout;
0934 edm::LogVerbatim("HGCalGeom") << st3.str();
0935 #endif
0936 }
0937
0938 if (trforms.find(std::make_pair(lay, zside)) == trforms.end()) {
0939 DD3Vector x, y, z;
0940 fv2.rotation().GetComponents(x, y, z);
0941 const CLHEP::HepRep3x3 rotation(x.X(), y.X(), z.X(), x.Y(), y.Y(), z.Y(), x.Z(), y.Z(), z.Z());
0942 const CLHEP::HepRotation hr(rotation);
0943 double xx = ((std::abs(fv2.translation().X()) < tolerance)
0944 ? 0
0945 : HGCalParameters::k_ScaleFromDDD * fv2.translation().X());
0946 double yy = ((std::abs(fv2.translation().Y()) < tolerance)
0947 ? 0
0948 : HGCalParameters::k_ScaleFromDDD * fv2.translation().Y());
0949 const CLHEP::Hep3Vector h3v(xx, yy, zvals[std::make_pair(lay, zside)]);
0950 HGCalParameters::hgtrform mytrf;
0951 mytrf.zp = zside;
0952 mytrf.lay = lay;
0953 mytrf.sec = 0;
0954 mytrf.subsec = 0;
0955 mytrf.h3v = h3v;
0956 mytrf.hr = hr;
0957 trforms[std::make_pair(lay, zside)] = mytrf;
0958 #ifdef EDM_ML_DEBUG
0959 edm::LogVerbatim("HGCalGeom") << "Translation " << h3v;
0960 #endif
0961 }
0962 }
0963 }
0964 }
0965 dodet = fv2.next();
0966 }
0967 #ifdef EDM_ML_DEBUG
0968 edm::LogVerbatim("HGCalGeom") << "Total # of views " << ntot1 << ":" << ntot2;
0969 #endif
0970 loadGeometryHexagon8(layers, trforms, firstLayer, php);
0971 }
0972
0973 void HGCalGeomParameters::loadGeometryHexagonModule(const cms::DDCompactView* cpv,
0974 HGCalParameters& php,
0975 const std::string& sdTag1,
0976 const std::string& sdTag2,
0977 int firstLayer) {
0978 #ifdef EDM_ML_DEBUG
0979 edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters (DD4hep)::loadGeometryHexagonModule called with tags " << sdTag1
0980 << ":" << sdTag2 << " firstLayer " << firstLayer;
0981 int ntot1(0), ntot2(0);
0982 #endif
0983 std::map<int, HGCalGeomParameters::layerParameters> layers;
0984 std::map<std::pair<int, int>, HGCalParameters::hgtrform> trforms;
0985 std::map<std::pair<int, int>, double> zvals;
0986 int levelTop = php.levelT_[0];
0987
0988 const cms::DDFilter filter1("Volume", sdTag2);
0989 cms::DDFilteredView fv1((*cpv), filter1);
0990 while (fv1.firstChild()) {
0991 #ifdef EDM_ML_DEBUG
0992 ++ntot1;
0993 #endif
0994 int nsiz = static_cast<int>(fv1.level());
0995 if (nsiz > levelTop) {
0996 std::vector<int> copy = fv1.copyNos();
0997 int lay = copy[nsiz - levelTop - 1];
0998 int zside = (nsiz > php.levelZSide_) ? copy[nsiz - php.levelZSide_ - 1] : -1;
0999 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
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_]
1061 : HGCalParameters::k_ScaleFromDD4hep * pars[1];
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
1084 : HGCalParameters::k_ScaleFromDD4hep * fv2.translation().X());
1085 double yy = ((std::abs(fv2.translation().Y()) < tolerance)
1086 ? 0
1087 : HGCalParameters::k_ScaleFromDD4hep * fv2.translation().Y());
1088 const CLHEP::Hep3Vector h3v(xx, yy, zvals[std::make_pair(lay, zside)]);
1089 HGCalParameters::hgtrform mytrf;
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
1172 void HGCalGeomParameters::loadSpecParsHexagon(const DDFilteredView& fv,
1173 HGCalParameters& php,
1174 const DDCompactView* cpv,
1175 const std::string& sdTag1,
1176 const std::string& sdTag2) {
1177 DDsvalues_type sv(fv.mergedSpecifics());
1178 php.boundR_ = getDDDArray("RadiusBound", sv, 4);
1179 rescale(php.boundR_, HGCalParameters::k_ScaleFromDDD);
1180 php.rLimit_ = getDDDArray("RadiusLimits", sv, 2);
1181 rescale(php.rLimit_, HGCalParameters::k_ScaleFromDDD);
1182 php.levelT_ = dbl_to_int(getDDDArray("LevelTop", sv, 0));
1183
1184
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
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
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
1216 void HGCalGeomParameters::loadSpecParsHexagon(const cms::DDFilteredView& fv,
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");
1223 rescale(php.boundR_, HGCalParameters::k_ScaleFromDD4hep);
1224 php.rLimit_ = fv.get<std::vector<double> >(sdTag4, "RadiusLimits");
1225 rescale(php.rLimit_, HGCalParameters::k_ScaleFromDD4hep);
1226 php.levelT_ = dbl_to_int(fv.get<std::vector<double> >(sdTag4, "LevelTop"));
1227
1228
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
1237 const auto& dummy = fv.get<std::vector<double> >(sdTag2, "WaferSize");
1238 waferSize_ = dummy[0] * HGCalParameters::k_ScaleFromDD4hepToG4;
1239
1240
1241 php.cellSize_ = fv.get<std::vector<double> >(sdTag3, "CellSize");
1242 rescale(php.cellSize_, HGCalParameters::k_ScaleFromDD4hepToG4);
1243
1244
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
1255 void HGCalGeomParameters::loadSpecParsHexagon(const HGCalParameters& php) {
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
1276 void HGCalGeomParameters::loadSpecParsHexagon8(const DDFilteredView& fv, HGCalParameters& php) {
1277 DDsvalues_type sv(fv.mergedSpecifics());
1278 php.cellThickness_ = getDDDArray("CellThickness", sv, 3);
1279 rescale(php.cellThickness_, HGCalParameters::k_ScaleFromDDD);
1280 if ((php.mode_ == HGCalGeometryMode::Hexagon8Module) || (php.mode_ == HGCalGeometryMode::Hexagon8Cassette) ||
1281 (php.mode_ == HGCalGeometryMode::Hexagon8CalibCell)) {
1282 php.waferThickness_ = getDDDArray("WaferThickness", sv, 3);
1283 rescale(php.waferThickness_, HGCalParameters::k_ScaleFromDDD);
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];
1296 php.zMinForRad_ = HGCalParameters::k_ScaleFromDDD * dummy[3];
1297
1298 php.radiusMixBoundary_ = fv.vector("RadiusMixBoundary");
1299 rescale(php.radiusMixBoundary_, HGCalParameters::k_ScaleFromDDD);
1300
1301 php.slopeMin_ = getDDDArray("SlopeBottom", sv, 0);
1302 php.zFrontMin_ = getDDDArray("ZFrontBottom", sv, 0);
1303 rescale(php.zFrontMin_, HGCalParameters::k_ScaleFromDDD);
1304 php.rMinFront_ = getDDDArray("RMinFront", sv, 0);
1305 rescale(php.rMinFront_, HGCalParameters::k_ScaleFromDDD);
1306
1307 php.slopeTop_ = getDDDArray("SlopeTop", sv, 0);
1308 php.zFrontTop_ = getDDDArray("ZFrontTop", sv, 0);
1309 rescale(php.zFrontTop_, HGCalParameters::k_ScaleFromDDD);
1310 php.rMaxFront_ = getDDDArray("RMaxFront", sv, 0);
1311 rescale(php.rMaxFront_, HGCalParameters::k_ScaleFromDDD);
1312
1313 php.zRanges_ = fv.vector("ZRanges");
1314 rescale(php.zRanges_, HGCalParameters::k_ScaleFromDDD);
1315
1316 const auto& dummy2 = getDDDArray("LayerOffset", sv, 1);
1317 php.layerOffset_ = dummy2[0];
1318 php.layerCenter_ = dbl_to_int(fv.vector("LayerCenter"));
1319
1320 if (php.mode_ == HGCalGeometryMode::Hexagon8CalibCell) {
1321 const auto& dummy3 = fv.vector("CalibCellRadius");
1322 php.calibCellRHD_ = HGCalParameters::k_ScaleFromDDD * dummy3[0];
1323 php.calibCellFullHD_ = dbl_to_int(fv.vector("CalibCellFullHD"));
1324 php.calibCellPartHD_ = dbl_to_int(fv.vector("CalibCellPartHD"));
1325 php.calibCellRLD_ = HGCalParameters::k_ScaleFromDDD * dummy3[1];
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
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 }
1357 if ((php.mode_ == HGCalGeometryMode::Hexagon8Module) || (php.mode_ == HGCalGeometryMode::Hexagon8Cassette) ||
1358 (php.mode_ == HGCalGeometryMode::Hexagon8CalibCell)) {
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;
1367 rescale(php.cassetteShift_, HGCalParameters::k_ScaleFromDDD);
1368 loadSpecParsHexagon8(php, layerType, waferIndex, waferProperties);
1369 }
1370 }
1371
1372 void HGCalGeomParameters::loadSpecParsHexagon8(const cms::DDFilteredView& fv,
1373 const cms::DDVectorsMap& vmap,
1374 HGCalParameters& php,
1375 const std::string& sdTag1) {
1376 php.cellThickness_ = fv.get<std::vector<double> >(sdTag1, "CellThickness");
1377 rescale(php.cellThickness_, HGCalParameters::k_ScaleFromDD4hep);
1378 if ((php.mode_ == HGCalGeometryMode::Hexagon8Module) || (php.mode_ == HGCalGeometryMode::Hexagon8Cassette) ||
1379 (php.mode_ == HGCalGeometryMode::Hexagon8CalibCell)) {
1380 php.waferThickness_ = fv.get<std::vector<double> >(sdTag1, "WaferThickness");
1381 rescale(php.waferThickness_, HGCalParameters::k_ScaleFromDD4hep);
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];
1395 php.zMinForRad_ = HGCalParameters::k_ScaleFromDD4hep * dummy[3];
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");
1402 rescale(php.zFrontMin_, HGCalParameters::k_ScaleFromDD4hep);
1403 php.rMinFront_ = fv.get<std::vector<double> >(sdTag1, "RMinFront");
1404 rescale(php.rMinFront_, HGCalParameters::k_ScaleFromDD4hep);
1405
1406 php.slopeTop_ = fv.get<std::vector<double> >(sdTag1, "SlopeTop");
1407 php.zFrontTop_ = fv.get<std::vector<double> >(sdTag1, "ZFrontTop");
1408 rescale(php.zFrontTop_, HGCalParameters::k_ScaleFromDD4hep);
1409 php.rMaxFront_ = fv.get<std::vector<double> >(sdTag1, "RMaxFront");
1410 rescale(php.rMaxFront_, HGCalParameters::k_ScaleFromDD4hep);
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)
1430 php.radiusMixBoundary_.emplace_back(HGCalParameters::k_ScaleFromDD4hep * i);
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")) {
1438 php.calibCellRHD_ = HGCalParameters::k_ScaleFromDD4hep * it.second[0];
1439 php.calibCellRLD_ = HGCalParameters::k_ScaleFromDDD * it.second[1];
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
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 }
1498 if ((php.mode_ == HGCalGeometryMode::Hexagon8Module) || (php.mode_ == HGCalGeometryMode::Hexagon8Cassette) ||
1499 (php.mode_ == HGCalGeometryMode::Hexagon8CalibCell)) {
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;
1518 rescale(php.cassetteShift_, HGCalParameters::k_ScaleFromDD4hep);
1519 loadSpecParsHexagon8(php, layerType, waferIndex, waferProperties);
1520 }
1521 }
1522
1523 void HGCalGeomParameters::loadSpecParsHexagon8(HGCalParameters& php) {
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
1571 void HGCalGeomParameters::loadSpecParsHexagon8(HGCalParameters& php,
1572 const std::vector<int>& layerType,
1573 const std::vector<int>& waferIndex,
1574 const std::vector<int>& waferProperties) {
1575
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 =
1594 ((php.mode_ == HGCalGeometryMode::Hexagon8Cassette) || (php.mode_ == HGCalGeometryMode::Hexagon8CalibCell))
1595 ? HGCalProperty::waferOrient(waferProperties[k])
1596 : HGCalWaferMask::getRotation(php.waferZSide_, partial, HGCalProperty::waferOrient(waferProperties[k]));
1597 php.waferInfoMap_[waferIndex[k]] = HGCalParameters::waferInfo(HGCalProperty::waferThick(waferProperties[k]),
1598 partial,
1599 orient,
1600 HGCalProperty::waferCassette(waferProperties[k]));
1601 #ifdef EDM_ML_DEBUG
1602 edm::LogVerbatim("HGCalGeom") << "[" << k << ":" << waferIndex[k] << ":"
1603 << HGCalWaferIndex::waferLayer(waferIndex[k]) << ":"
1604 << HGCalWaferIndex::waferU(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
1612 void HGCalGeomParameters::loadSpecParsTrapezoid(const DDFilteredView& fv, HGCalParameters& php) {
1613 DDsvalues_type sv(fv.mergedSpecifics());
1614 php.radiusMixBoundary_ = fv.vector("RadiusMixBoundary");
1615 rescale(php.radiusMixBoundary_, HGCalParameters::k_ScaleFromDDD);
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);
1620 rescale(php.rMinLayerBH_, HGCalParameters::k_ScaleFromDDD);
1621 assert(php.nPhiBinBH_.size() > 1);
1622 php.nCellsFine_ = php.nPhiBinBH_[1];
1623 php.nCellsCoarse_ = php.nPhiBinBH_[0];
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);
1631 rescale(php.zFrontMin_, HGCalParameters::k_ScaleFromDDD);
1632 php.rMinFront_ = getDDDArray("RMinFront", sv, 0);
1633 rescale(php.rMinFront_, HGCalParameters::k_ScaleFromDDD);
1634
1635 php.slopeTop_ = getDDDArray("SlopeTop", sv, 0);
1636 php.zFrontTop_ = getDDDArray("ZFrontTop", sv, 0);
1637 rescale(php.zFrontTop_, HGCalParameters::k_ScaleFromDDD);
1638 php.rMaxFront_ = getDDDArray("RMaxFront", sv, 0);
1639 rescale(php.rMaxFront_, HGCalParameters::k_ScaleFromDDD);
1640
1641 php.zRanges_ = fv.vector("ZRanges");
1642 rescale(php.zRanges_, HGCalParameters::k_ScaleFromDDD);
1643
1644
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 << ":"
1652 << scintillatorCassette << ":" << scintillatorFineCell;
1653 #endif
1654
1655 if ((php.waferMaskMode_ == scintillatorFile) || (php.waferMaskMode_ == scintillatorCassette) ||
1656 (php.waferMaskMode_ == scintillatorFineCell)) {
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");
1671 rescale(tileRMin, HGCalParameters::k_ScaleFromDDD);
1672 rescale(tileRMax, HGCalParameters::k_ScaleFromDDD);
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");
1680 rescale(tileRMinFine, HGCalParameters::k_ScaleFromDDD);
1681 rescale(tileRMaxFine, HGCalParameters::k_ScaleFromDDD);
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");
1687 rescale(rectract, HGCalParameters::k_ScaleFromDDD);
1688 for (unsigned int k1 = 0; k1 < rectract.size(); ++k1) {
1689 php.cassetteRetractTile_.emplace_back(rectract[k1]);
1690 #ifdef EDM_ML_DEBUG
1691 edm::LogVerbatim("HGCalGeom") << "cassetteRetractTile_[" << k1 << "] " << rectract[k1];
1692 #endif
1693 }
1694 int n = 2 * php.cassettes_ * (php.firstLayer_ - 1);
1695 for (int k1 = 0; k1 < n; ++k1)
1696 cassetteShift.emplace_back(0.);
1697 #ifdef EDM_ML_DEBUG
1698 edm::LogVerbatim("HGCalGeom") << "First Layer " << php.firstLayer_ << " Rett size " << rectract.size() << " N "
1699 << n;
1700 #endif
1701 double dphi = M_PI / php.cassettes_;
1702 for (unsigned int k1 = 0; k1 < rectract.size(); ++k1) {
1703 for (int k2 = 0; k2 < php.cassettes_; ++k2) {
1704 double phi = (2 * k2 + 1) * dphi;
1705 cassetteShift.emplace_back(rectract[k1] * cos(phi));
1706 cassetteShift.emplace_back(rectract[k1] * sin(phi));
1707 }
1708 }
1709 } else if (php.waferMaskMode_ == scintillatorCassette) {
1710 if (php.cassettes_ > 0)
1711 php.nphiCassette_ = php.nCellsCoarse_ / php.cassettes_;
1712 cassetteShift = fv.vector("CassetteShiftHE");
1713 }
1714
1715 rescale(cassetteShift, HGCalParameters::k_ScaleFromDDD);
1716 if (php.waferMaskMode_ == scintillatorFineCell)
1717 php.cassetteShiftTile_ = cassetteShift;
1718 else
1719 php.cassetteShift_ = cassetteShift;
1720 loadSpecParsTrapezoid(php,
1721 tileIndx,
1722 tileProperty,
1723 tileHEX1,
1724 tileHEX2,
1725 tileHEX3,
1726 tileHEX4,
1727 tileHEX5,
1728 tileHEX6,
1729 tileRMin,
1730 tileRMax,
1731 tileRMinFine,
1732 tileRMaxFine,
1733 tileRingMin,
1734 tileRingMax,
1735 tileRingMinFine,
1736 tileRingMaxFine);
1737 }
1738 }
1739
1740 void HGCalGeomParameters::loadSpecParsTrapezoid(const cms::DDFilteredView& fv,
1741 const cms::DDVectorsMap& vmap,
1742 HGCalParameters& php,
1743 const std::string& sdTag1) {
1744 for (auto const& it : vmap) {
1745 if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "RadiusMixBoundary")) {
1746 for (const auto& i : it.second)
1747 php.radiusMixBoundary_.emplace_back(HGCalParameters::k_ScaleFromDD4hep * i);
1748 } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "ZRanges")) {
1749 for (const auto& i : it.second)
1750 php.zRanges_.emplace_back(HGCalParameters::k_ScaleFromDD4hep * i);
1751 } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "LayerCenter")) {
1752 for (const auto& i : it.second)
1753 php.layerCenter_.emplace_back(std::round(i));
1754 }
1755 }
1756
1757 php.nPhiBinBH_ = dbl_to_int(fv.get<std::vector<double> >(sdTag1, "NPhiBinBH"));
1758 php.layerFrontBH_ = dbl_to_int(fv.get<std::vector<double> >(sdTag1, "LayerFrontBH"));
1759 php.rMinLayerBH_ = fv.get<std::vector<double> >(sdTag1, "RMinLayerBH");
1760 rescale(php.rMinLayerBH_, HGCalParameters::k_ScaleFromDD4hep);
1761 assert(php.nPhiBinBH_.size() > 1);
1762 php.nCellsFine_ = php.nPhiBinBH_[1];
1763 php.nCellsCoarse_ = php.nPhiBinBH_[0];
1764 assert(0 != php.nCellsFine_);
1765 assert(0 != php.nCellsCoarse_);
1766 php.cellSize_.emplace_back(2.0 * M_PI / php.nCellsFine_);
1767 php.cellSize_.emplace_back(2.0 * M_PI / php.nCellsCoarse_);
1768
1769 php.slopeMin_ = fv.get<std::vector<double> >(sdTag1, "SlopeBottom");
1770 php.zFrontMin_ = fv.get<std::vector<double> >(sdTag1, "ZFrontBottom");
1771 rescale(php.zFrontMin_, HGCalParameters::k_ScaleFromDD4hep);
1772 php.rMinFront_ = fv.get<std::vector<double> >(sdTag1, "RMinFront");
1773 rescale(php.rMinFront_, HGCalParameters::k_ScaleFromDD4hep);
1774
1775 php.slopeTop_ = fv.get<std::vector<double> >(sdTag1, "SlopeTop");
1776 php.zFrontTop_ = fv.get<std::vector<double> >(sdTag1, "ZFrontTop");
1777 rescale(php.zFrontTop_, HGCalParameters::k_ScaleFromDD4hep);
1778 php.rMaxFront_ = fv.get<std::vector<double> >(sdTag1, "RMaxFront");
1779 rescale(php.rMaxFront_, HGCalParameters::k_ScaleFromDD4hep);
1780 unsigned int kmax = (php.zFrontTop_.size() - php.slopeTop_.size());
1781 for (unsigned int k = 0; k < kmax; ++k)
1782 php.slopeTop_.emplace_back(0);
1783
1784 const auto& dummy2 = fv.get<std::vector<double> >(sdTag1, "LayerOffset");
1785 php.layerOffset_ = dummy2[0];
1786
1787 loadSpecParsTrapezoid(php);
1788
1789
1790 if ((php.waferMaskMode_ == scintillatorFile) || (php.waferMaskMode_ == scintillatorCassette) ||
1791 (php.waferMaskMode_ == scintillatorFineCell)) {
1792 std::vector<int> tileIndx, tileProperty;
1793 std::vector<int> tileHEX1, tileHEX2, tileHEX3, tileHEX4, tileHEX5, tileHEX6;
1794 std::vector<double> tileRMin, tileRMax, tileRMinFine, tileRMaxFine;
1795 std::vector<int> tileRingMin, tileRingMax, tileRingMinFine, tileRingMaxFine;
1796 std::vector<double> cassetteShift;
1797 for (auto const& it : vmap) {
1798 if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileIndex")) {
1799 for (const auto& i : it.second)
1800 tileIndx.emplace_back(std::round(i));
1801 } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileProperty")) {
1802 for (const auto& i : it.second)
1803 tileProperty.emplace_back(std::round(i));
1804 } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "NPhiLayer")) {
1805 for (const auto& i : it.second)
1806 php.nPhiLayer_.emplace_back(std::round(i));
1807 } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileHEX1")) {
1808 for (const auto& i : it.second)
1809 tileHEX1.emplace_back(std::round(i));
1810 } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileHEX2")) {
1811 for (const auto& i : it.second)
1812 tileHEX2.emplace_back(std::round(i));
1813 } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileHEX3")) {
1814 for (const auto& i : it.second)
1815 tileHEX3.emplace_back(std::round(i));
1816 } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileHEX4")) {
1817 for (const auto& i : it.second)
1818 tileHEX4.emplace_back(std::round(i));
1819 } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileRMin")) {
1820 for (const auto& i : it.second)
1821 tileRMin.emplace_back(HGCalParameters::k_ScaleFromDD4hep * i);
1822 } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileRMax")) {
1823 for (const auto& i : it.second)
1824 tileRMax.emplace_back(HGCalParameters::k_ScaleFromDD4hep * i);
1825 } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileRingMin")) {
1826 for (const auto& i : it.second)
1827 tileRingMin.emplace_back(std::round(i));
1828 } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileRingMax")) {
1829 for (const auto& i : it.second)
1830 tileRingMax.emplace_back(std::round(i));
1831 }
1832 }
1833 if (php.waferMaskMode_ == scintillatorFineCell) {
1834 php.nphiCassette_ = php.nCellsCoarse_ / php.cassettes_;
1835 php.nphiFineCassette_ = php.nCellsFine_ / php.cassettes_;
1836 std::vector<double> rectract;
1837 for (auto const& it : vmap) {
1838 if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileHEX5")) {
1839 for (const auto& i : it.second)
1840 tileHEX5.emplace_back(std::round(i));
1841 } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileHEX6")) {
1842 for (const auto& i : it.second)
1843 tileHEX6.emplace_back(std::round(i));
1844 } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileRMin6")) {
1845 for (const auto& i : it.second)
1846 tileRMinFine.emplace_back(i);
1847 rescale(tileRMinFine, HGCalParameters::k_ScaleFromDDD);
1848 } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileRMax6")) {
1849 for (const auto& i : it.second)
1850 tileRMaxFine.emplace_back(i);
1851 rescale(tileRMaxFine, HGCalParameters::k_ScaleFromDDD);
1852 } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileRingMin6")) {
1853 for (const auto& i : it.second)
1854 tileRingMinFine.emplace_back(std::round(i));
1855 } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "TileRingMax6")) {
1856 for (const auto& i : it.second)
1857 tileRingMaxFine.emplace_back(std::round(i));
1858 } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "ScintRetract")) {
1859 for (const auto& i : it.second)
1860 rectract.emplace_back(i);
1861 rescale(rectract, HGCalParameters::k_ScaleFromDDD);
1862 for (unsigned int k1 = 0; k1 < rectract.size(); ++k1) {
1863 php.cassetteRetractTile_.emplace_back(rectract[k1]);
1864 #ifdef EDM_ML_DEBUG
1865 edm::LogVerbatim("HGCalGeom") << "cassetteRetractTile_[" << k1 << "] " << rectract[k1];
1866 #endif
1867 }
1868 int n = 2 * php.cassettes_ * (php.firstLayer_ - 1);
1869 for (int k1 = 0; k1 < n; ++k1)
1870 cassetteShift.emplace_back(0.);
1871 #ifdef EDM_ML_DEBUG
1872 edm::LogVerbatim("HGCalGeom") << "First Layer " << php.firstLayer_ << " Rett size " << rectract.size()
1873 << " N " << n;
1874 #endif
1875 double dphi = M_PI / php.cassettes_;
1876 for (unsigned int k1 = 0; k1 < rectract.size(); ++k1) {
1877 for (int k2 = 0; k2 < php.cassettes_; ++k2) {
1878 double phi = (2 * k2 + 1) * dphi;
1879 cassetteShift.emplace_back(rectract[k1] * cos(phi));
1880 cassetteShift.emplace_back(rectract[k1] * sin(phi));
1881 }
1882 }
1883 }
1884 }
1885 } else if (php.waferMaskMode_ == scintillatorCassette) {
1886 for (auto const& it : vmap) {
1887 if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "CassetteShiftHE")) {
1888 for (const auto& i : it.second)
1889 cassetteShift.emplace_back(i);
1890 }
1891 }
1892 }
1893
1894 rescale(cassetteShift, HGCalParameters::k_ScaleFromDD4hep);
1895 if (php.waferMaskMode_ == scintillatorFineCell)
1896 php.cassetteShiftTile_ = cassetteShift;
1897 else
1898 php.cassetteShift_ = cassetteShift;
1899 loadSpecParsTrapezoid(php,
1900 tileIndx,
1901 tileProperty,
1902 tileHEX1,
1903 tileHEX2,
1904 tileHEX3,
1905 tileHEX4,
1906 tileHEX5,
1907 tileHEX6,
1908 tileRMin,
1909 tileRMax,
1910 tileRMinFine,
1911 tileRMaxFine,
1912 tileRingMin,
1913 tileRingMax,
1914 tileRingMinFine,
1915 tileRingMaxFine);
1916 }
1917 }
1918
1919 void HGCalGeomParameters::loadSpecParsTrapezoid(HGCalParameters& php) {
1920 #ifdef EDM_ML_DEBUG
1921 edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters:nCells " << php.nCellsFine_ << ":" << php.nCellsCoarse_
1922 << " cellSize: " << php.cellSize_[0] << ":" << php.cellSize_[1];
1923 for (unsigned int k = 0; k < php.layerFrontBH_.size(); ++k)
1924 edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: Type[" << k << "] Front Layer = " << php.layerFrontBH_[k]
1925 << " rMin = " << php.rMinLayerBH_[k];
1926 for (unsigned int k = 0; k < php.radiusMixBoundary_.size(); ++k) {
1927 edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: Mix[" << k << "] R = " << php.radiusMixBoundary_[k]
1928 << " Nphi = " << php.scintCells(k + php.firstLayer_)
1929 << " dPhi = " << php.scintCellSize(k + php.firstLayer_);
1930 }
1931
1932 for (unsigned int k = 0; k < php.zFrontMin_.size(); ++k)
1933 edm::LogVerbatim("HGCalGeom") << "HGCalParameters: Boundary[" << k << "] Bottom Z = " << php.zFrontMin_[k]
1934 << " Slope = " << php.slopeMin_[k] << " rMax = " << php.rMinFront_[k];
1935
1936 for (unsigned int k = 0; k < php.zFrontTop_.size(); ++k)
1937 edm::LogVerbatim("HGCalGeom") << "HGCalParameters: Boundary[" << k << "] Top Z = " << php.zFrontTop_[k]
1938 << " Slope = " << php.slopeTop_[k] << " rMax = " << php.rMaxFront_[k];
1939
1940 std::ostringstream st1;
1941 for (unsigned int k = 0; k < php.zRanges_.size(); ++k)
1942 st1 << ":" << php.zRanges_[k];
1943 edm::LogVerbatim("HGCalGeom") << "HGCalParameters: Z-Boundary[" << php.zRanges_.size() << "] " << st1.str();
1944
1945 edm::LogVerbatim("HGCalGeom") << "HGCalParameters: LayerOffset " << php.layerOffset_ << " in array of size "
1946 << php.layerCenter_.size();
1947 for (unsigned int k = 0; k < php.layerCenter_.size(); ++k)
1948 edm::LogVerbatim("HGCalGeom") << "[" << k << "] " << php.layerCenter_[k];
1949 #endif
1950 }
1951
1952 void HGCalGeomParameters::loadSpecParsTrapezoid(HGCalParameters& php,
1953 const std::vector<int>& tileIndx,
1954 const std::vector<int>& tileProperty,
1955 const std::vector<int>& tileHEX1,
1956 const std::vector<int>& tileHEX2,
1957 const std::vector<int>& tileHEX3,
1958 const std::vector<int>& tileHEX4,
1959 const std::vector<int>& tileHEX5,
1960 const std::vector<int>& tileHEX6,
1961 const std::vector<double>& tileRMin,
1962 const std::vector<double>& tileRMax,
1963 const std::vector<double>& tileRMinFine,
1964 const std::vector<double>& tileRMaxFine,
1965 const std::vector<int>& tileRingMin,
1966 const std::vector<int>& tileRingMax,
1967 const std::vector<int>& tileRingMinFine,
1968 const std::vector<int>& tileRingMaxFine) {
1969
1970 for (unsigned int k = 0; k < tileIndx.size(); ++k) {
1971 int hex5 = (k < tileHEX5.size()) ? tileHEX5[k] : 0;
1972 int hex6 = (k < tileHEX6.size()) ? tileHEX6[k] : 0;
1973 php.tileInfoMap_[tileIndx[k]] = HGCalParameters::tileInfo(HGCalProperty::tileType(tileProperty[k]),
1974 HGCalProperty::tileSiPM(tileProperty[k]),
1975 tileHEX1[k],
1976 tileHEX2[k],
1977 tileHEX3[k],
1978 tileHEX4[k],
1979 hex5,
1980 hex6);
1981 #ifdef EDM_ML_DEBUG
1982 edm::LogVerbatim("HGCalGeom") << "Tile[" << k << ":" << tileIndx[k] << "] "
1983 << " Type " << HGCalProperty::tileType(tileProperty[k]) << " SiPM "
1984 << HGCalProperty::tileSiPM(tileProperty[k]) << " HEX " << std::hex << tileHEX1[k]
1985 << ":" << tileHEX2[k] << ":" << tileHEX3[k] << ":" << tileHEX4[k] << ":" << hex5
1986 << ":" << hex6 << std::dec;
1987 #endif
1988 }
1989
1990 for (unsigned int k = 0; k < tileRMinFine.size(); ++k) {
1991 php.tileRingFineR_.emplace_back(tileRMinFine[k], tileRMaxFine[k]);
1992 #ifdef EDM_ML_DEBUG
1993 edm::LogVerbatim("HGCalGeom") << "TileRingFineR[" << k << "] " << tileRMinFine[k] << ":" << tileRMaxFine[k];
1994 #endif
1995 }
1996 for (unsigned int k = 0; k < tileRMin.size(); ++k) {
1997 php.tileRingR_.emplace_back(tileRMin[k], tileRMax[k]);
1998 #ifdef EDM_ML_DEBUG
1999 edm::LogVerbatim("HGCalGeom") << "TileRingR[" << k << "] " << tileRMin[k] << ":" << tileRMax[k];
2000 #endif
2001 }
2002
2003 for (unsigned int k = 0; k < tileRingMinFine.size(); ++k) {
2004 php.tileRingFineRange_.emplace_back(tileRingMinFine[k], tileRingMaxFine[k]);
2005 #ifdef EDM_ML_DEBUG
2006 edm::LogVerbatim("HGCalGeom") << "TileRingFineRange[" << k << "] " << tileRingMinFine[k] << ":"
2007 << tileRingMaxFine[k];
2008 #endif
2009 }
2010 for (unsigned k = 0; k < tileRingMin.size(); ++k) {
2011 php.tileRingRange_.emplace_back(tileRingMin[k], tileRingMax[k]);
2012 #ifdef EDM_ML_DEBUG
2013 edm::LogVerbatim("HGCalGeom") << "TileRingRange[" << k << "] " << tileRingMin[k] << ":" << tileRingMax[k];
2014 #endif
2015 }
2016 #ifdef EDM_ML_DEBUG
2017 std::ostringstream st1, st2;
2018 st1 << "Scintillator CassetteShift for " << php.cassetteShiftTile_.size() << " tiles:";
2019 for (unsigned int k = 0; k < php.cassetteShiftTile_.size(); ++k)
2020 st1 << " " << php.cassetteShiftTile_[k];
2021 edm::LogVerbatim("HGCalGeom") << st1.str();
2022 st2 << "NPhi for scntillator " << php.nPhiLayer_.size() << " layers:";
2023 for (unsigned int k = 0; k < php.nPhiLayer_.size(); ++k)
2024 st2 << " " << php.nPhiLayer_[k];
2025 edm::LogVerbatim("HGCalGeom") << st2.str();
2026 #endif
2027 }
2028
2029 void HGCalGeomParameters::loadWaferHexagon(HGCalParameters& php) {
2030 double waferW(HGCalParameters::k_ScaleFromDDD * waferSize_), rmin(HGCalParameters::k_ScaleFromDDD * php.waferR_);
2031 double rin(php.rLimit_[0]), rout(php.rLimit_[1]), rMaxFine(php.boundR_[1]);
2032 #ifdef EDM_ML_DEBUG
2033 edm::LogVerbatim("HGCalGeom") << "Input waferWidth " << waferW << ":" << rmin << " R Limits: " << rin << ":" << rout
2034 << " Fine " << rMaxFine;
2035 #endif
2036
2037 php.waferCopy_.clear();
2038 php.waferTypeL_.clear();
2039 php.waferTypeT_.clear();
2040 php.waferPosX_.clear();
2041 php.waferPosY_.clear();
2042 double dx = 0.5 * waferW;
2043 double dy = 3.0 * dx * tan(30._deg);
2044 double rr = 2.0 * dx * tan(30._deg);
2045 int ncol = static_cast<int>(2.0 * rout / waferW) + 1;
2046 int nrow = static_cast<int>(rout / (waferW * tan(30._deg))) + 1;
2047 int ns2 = (2 * ncol + 1) * (2 * nrow + 1) * php.layer_.size();
2048 int incm(0), inrm(0);
2049 HGCalParameters::layer_map copiesInLayers(php.layer_.size() + 1);
2050 HGCalParameters::waferT_map waferTypes(ns2 + 1);
2051 #ifdef EDM_ML_DEBUG
2052 int kount(0), ntot(0);
2053 edm::LogVerbatim("HGCalGeom") << "Row " << nrow << " Column " << ncol;
2054 #endif
2055 for (int nr = -nrow; nr <= nrow; ++nr) {
2056 int inr = (nr >= 0) ? nr : -nr;
2057 for (int nc = -ncol; nc <= ncol; ++nc) {
2058 int inc = (nc >= 0) ? nc : -nc;
2059 if (inr % 2 == inc % 2) {
2060 double xpos = nc * dx;
2061 double ypos = nr * dy;
2062 std::pair<int, int> corner = HGCalGeomTools::waferCorner(xpos, ypos, dx, rr, rin, rout, true);
2063 double rpos = std::sqrt(xpos * xpos + ypos * ypos);
2064 int typet = (rpos < rMaxFine) ? 1 : 2;
2065 int typel(3);
2066 for (int k = 1; k < 4; ++k) {
2067 if ((rpos + rmin) <= php.boundR_[k]) {
2068 typel = k;
2069 break;
2070 }
2071 }
2072 #ifdef EDM_ML_DEBUG
2073 ++ntot;
2074 #endif
2075 if (corner.first > 0) {
2076 int copy = HGCalTypes::packTypeUV(typel, nc, nr);
2077 if (inc > incm)
2078 incm = inc;
2079 if (inr > inrm)
2080 inrm = inr;
2081 #ifdef EDM_ML_DEBUG
2082 kount++;
2083 edm::LogVerbatim("HGCalGeom") << kount << ":" << ntot << " Copy " << copy << " Type " << typel << ":" << typet
2084 << " Location " << corner.first << " Position " << xpos << ":" << ypos
2085 << " Layers " << php.layer_.size();
2086 #endif
2087 php.waferCopy_.emplace_back(copy);
2088 php.waferTypeL_.emplace_back(typel);
2089 php.waferTypeT_.emplace_back(typet);
2090 php.waferPosX_.emplace_back(xpos);
2091 php.waferPosY_.emplace_back(ypos);
2092 for (unsigned int il = 0; il < php.layer_.size(); ++il) {
2093 std::pair<int, int> corner =
2094 HGCalGeomTools::waferCorner(xpos, ypos, dx, rr, php.rMinLayHex_[il], php.rMaxLayHex_[il], true);
2095 if (corner.first > 0) {
2096 auto cpy = copiesInLayers[php.layer_[il]].find(copy);
2097 if (cpy == copiesInLayers[php.layer_[il]].end())
2098 copiesInLayers[php.layer_[il]][copy] =
2099 ((corner.first == static_cast<int>(HGCalParameters::k_CornerSize)) ? php.waferCopy_.size() : -1);
2100 }
2101 if ((corner.first > 0) && (corner.first < static_cast<int>(HGCalParameters::k_CornerSize))) {
2102 int wl = HGCalWaferIndex::waferIndex(php.layer_[il], copy, 0, true);
2103 waferTypes[wl] = corner;
2104 }
2105 }
2106 }
2107 }
2108 }
2109 }
2110 php.copiesInLayers_ = copiesInLayers;
2111 php.waferTypes_ = waferTypes;
2112 php.nSectors_ = static_cast<int>(php.waferCopy_.size());
2113 php.waferUVMax_ = 0;
2114 #ifdef EDM_ML_DEBUG
2115 edm::LogVerbatim("HGCalGeom") << "HGCalWaferHexagon: # of columns " << incm << " # of rows " << inrm << " and "
2116 << kount << ":" << ntot << " wafers; R " << rin << ":" << rout;
2117 edm::LogVerbatim("HGCalGeom") << "Dump copiesInLayers for " << php.copiesInLayers_.size() << " layers";
2118 for (unsigned int k = 0; k < copiesInLayers.size(); ++k) {
2119 const auto& theModules = copiesInLayers[k];
2120 edm::LogVerbatim("HGCalGeom") << "Layer " << k << ":" << theModules.size();
2121 int k2(0);
2122 for (std::unordered_map<int, int>::const_iterator itr = theModules.begin(); itr != theModules.end(); ++itr, ++k2) {
2123 edm::LogVerbatim("HGCalGeom") << "[" << k2 << "] " << itr->first << ":" << itr->second;
2124 }
2125 }
2126 #endif
2127 }
2128
2129 void HGCalGeomParameters::loadWaferHexagon8(HGCalParameters& php) {
2130 double waferW(php.waferSize_);
2131 double waferS(php.sensorSeparation_);
2132 auto wType = std::make_unique<HGCalWaferType>(php.radius100to200_,
2133 php.radius200to300_,
2134 HGCalParameters::k_ScaleToDDD * (waferW + waferS),
2135 HGCalParameters::k_ScaleToDDD * php.zMinForRad_,
2136 php.choiceType_,
2137 php.nCornerCut_,
2138 php.fracAreaMin_);
2139
2140 double rout(php.rLimit_[1]);
2141 #ifdef EDM_ML_DEBUG
2142 edm::LogVerbatim("HGCalGeom") << "Input waferWidth " << waferW << ":" << waferS << " R Max: " << rout;
2143 #endif
2144
2145 php.waferCopy_.clear();
2146 php.waferTypeL_.clear();
2147 php.waferTypeT_.clear();
2148 php.waferPosX_.clear();
2149 php.waferPosY_.clear();
2150 double r = 0.5 * (waferW + waferS);
2151 double R = 2.0 * r / sqrt3_;
2152 double dy = 0.75 * R;
2153 double r1 = 0.5 * waferW;
2154 double R1 = 2.0 * r1 / sqrt3_;
2155 int N = (r == 0) ? 3 : (static_cast<int>(0.5 * rout / r) + 3);
2156 int ns1 = (2 * N + 1) * (2 * N + 1);
2157 int ns2 = ns1 * php.zLayerHex_.size();
2158 #ifdef EDM_ML_DEBUG
2159 edm::LogVerbatim("HGCalGeom") << "wafer " << waferW << ":" << waferS << " r " << r << " dy " << dy << " N " << N
2160 << " sizes " << ns1 << ":" << ns2;
2161 std::vector<int> indtypes(ns1 + 1);
2162 indtypes.clear();
2163 #endif
2164 HGCalParameters::wafer_map wafersInLayers(ns1 + 1);
2165 HGCalParameters::wafer_map typesInLayers(ns2 + 1);
2166 HGCalParameters::waferT_map waferTypes(ns2 + 1);
2167 int ipos(0), lpos(0), uvmax(0), nwarn(0);
2168 std::vector<int> uvmx(php.zLayerHex_.size(), 0);
2169 for (int v = -N; v <= N; ++v) {
2170 for (int u = -N; u <= N; ++u) {
2171 int nr = 2 * v;
2172 int nc = -2 * u + v;
2173 double xpos = nc * r;
2174 double ypos = nr * dy;
2175 int indx = HGCalWaferIndex::waferIndex(0, u, v);
2176 php.waferCopy_.emplace_back(indx);
2177 php.waferPosX_.emplace_back(xpos);
2178 php.waferPosY_.emplace_back(ypos);
2179 wafersInLayers[indx] = ipos;
2180 ++ipos;
2181 std::pair<int, int> corner = HGCalGeomTools::waferCorner(xpos, ypos, r1, R1, 0, rout, false);
2182 if ((corner.first == static_cast<int>(HGCalParameters::k_CornerSize)) ||
2183 ((corner.first > 0) && php.defineFull_)) {
2184 uvmax = std::max(uvmax, std::max(std::abs(u), std::abs(v)));
2185 }
2186 for (unsigned int i = 0; i < php.zLayerHex_.size(); ++i) {
2187 int copy = i + php.layerOffset_;
2188 std::pair<double, double> xyoff = geomTools_.shiftXY(php.layerCenter_[copy], (waferW + waferS));
2189 int lay = php.layer_[php.layerIndex_[i]];
2190 double xpos0 = xpos + xyoff.first;
2191 double ypos0 = ypos + xyoff.second;
2192 double zpos = php.zLayerHex_[i];
2193 int kndx = HGCalWaferIndex::waferIndex(lay, u, v);
2194 int type(-1);
2195 if ((php.mode_ == HGCalGeometryMode::Hexagon8File) || (php.mode_ == HGCalGeometryMode::Hexagon8Module) ||
2196 (php.mode_ == HGCalGeometryMode::Hexagon8Cassette) || (php.mode_ == HGCalGeometryMode::Hexagon8CalibCell))
2197 type = wType->getType(kndx, php.waferInfoMap_);
2198 if (type < 0)
2199 type = wType->getType(HGCalParameters::k_ScaleToDDD * xpos0,
2200 HGCalParameters::k_ScaleToDDD * ypos0,
2201 HGCalParameters::k_ScaleToDDD * zpos);
2202 php.waferTypeL_.emplace_back(type);
2203 typesInLayers[kndx] = lpos;
2204 ++lpos;
2205 #ifdef EDM_ML_DEBUG
2206 indtypes.emplace_back(kndx);
2207 #endif
2208 std::pair<int, int> corner =
2209 HGCalGeomTools::waferCorner(xpos0, ypos0, r1, R1, php.rMinLayHex_[i], php.rMaxLayHex_[i], false);
2210 #ifdef EDM_ML_DEBUG
2211 if (((corner.first == 0) && std::abs(u) < 5 && std::abs(v) < 5) || (std::abs(u) < 2 && std::abs(v) < 2)) {
2212 edm::LogVerbatim("HGCalGeom") << "Layer " << lay << " R " << php.rMinLayHex_[i] << ":" << php.rMaxLayHex_[i]
2213 << " u " << u << " v " << v << " with " << corner.first << " corners";
2214 }
2215 #endif
2216 if ((corner.first == static_cast<int>(HGCalParameters::k_CornerSize)) ||
2217 ((corner.first > 0) && php.defineFull_)) {
2218 uvmx[i] = std::max(uvmx[i], std::max(std::abs(u), std::abs(v)));
2219 }
2220 if ((corner.first < static_cast<int>(HGCalParameters::k_CornerSize)) && (corner.first > 0)) {
2221 #ifdef EDM_ML_DEBUG
2222 edm::LogVerbatim("HGCalGeom") << "Layer " << lay << " u|v " << u << ":" << v << " with corner "
2223 << corner.first << ":" << corner.second;
2224 #endif
2225 int wl = HGCalWaferIndex::waferIndex(lay, u, v);
2226 if (php.waferMaskMode_ > 0) {
2227 bool v17OrLess = (php.mode_ < HGCalGeometryMode::Hexagon8CalibCell);
2228 std::pair<int, int> corner0 = HGCalWaferMask::getTypeMode(
2229 xpos0, ypos0, r1, R1, php.rMinLayHex_[i], php.rMaxLayHex_[i], type, php.waferMaskMode_, v17OrLess);
2230 if ((php.mode_ == HGCalGeometryMode::Hexagon8File) || (php.mode_ == HGCalGeometryMode::Hexagon8Module) ||
2231 (php.mode_ == HGCalGeometryMode::Hexagon8Cassette) ||
2232 (php.mode_ == HGCalGeometryMode::Hexagon8CalibCell)) {
2233 auto itr = php.waferInfoMap_.find(wl);
2234 if (itr != php.waferInfoMap_.end()) {
2235 int part = (itr->second).part;
2236 int orient = (itr->second).orient;
2237 bool ok = ((php.mode_ == HGCalGeometryMode::Hexagon8Cassette) ||
2238 (php.mode_ == HGCalGeometryMode::Hexagon8CalibCell))
2239 ? true
2240 : HGCalWaferMask::goodTypeMode(
2241 xpos0, ypos0, r1, R1, php.rMinLayHex_[i], php.rMaxLayHex_[i], part, orient, false);
2242 if (ok)
2243 corner0 = std::make_pair(part, (HGCalTypes::k_OffsetRotation + orient));
2244 #ifdef EDM_ML_DEBUG
2245 edm::LogVerbatim("HGCalGeom")
2246 << "Layer:u:v " << i << ":" << lay << ":" << u << ":" << v << " Part " << corner0.first << ":"
2247 << part << " Orient " << corner0.second << ":" << orient << " Position " << xpos0 << ":" << ypos0
2248 << " delta " << r1 << ":" << R1 << " Limit " << php.rMinLayHex_[i] << ":" << php.rMaxLayHex_[i]
2249 << " Compatibiliety Flag " << ok;
2250 #endif
2251 if (!ok)
2252 ++nwarn;
2253 }
2254 }
2255 waferTypes[wl] = corner0;
2256 #ifdef EDM_ML_DEBUG
2257 edm::LogVerbatim("HGCalGeom")
2258 << "Layer " << lay << " u|v " << u << ":" << v << " Index " << std::hex << wl << std::dec << " pos "
2259 << xpos0 << ":" << ypos0 << " R " << r1 << ":" << R1 << " Range " << php.rMinLayHex_[i] << ":"
2260 << php.rMaxLayHex_[i] << type << ":" << php.waferMaskMode_ << " corner " << corner.first << ":"
2261 << corner.second << " croner0 " << corner0.first << ":" << corner0.second;
2262 #endif
2263 } else {
2264 waferTypes[wl] = corner;
2265 #ifdef EDM_ML_DEBUG
2266 edm::LogVerbatim("HGCalGeom") << "Layer " << lay << " u|v " << u << ":" << v << " with corner "
2267 << corner.first << ":" << corner.second;
2268 #endif
2269 }
2270 }
2271 }
2272 }
2273 }
2274 if (nwarn > 0)
2275 edm::LogWarning("HGCalGeom") << "HGCalGeomParameters::loadWafer8: there are " << nwarn
2276 << " wafers with non-matching partial- orientation types";
2277 php.waferUVMax_ = uvmax;
2278 php.waferUVMaxLayer_ = uvmx;
2279 php.wafersInLayers_ = wafersInLayers;
2280 php.typesInLayers_ = typesInLayers;
2281 php.waferTypes_ = waferTypes;
2282 php.nSectors_ = static_cast<int>(php.waferCopy_.size());
2283 HGCalParameters::hgtrap mytr;
2284 mytr.lay = 1;
2285 mytr.bl = php.waferR_;
2286 mytr.tl = php.waferR_;
2287 mytr.h = php.waferR_;
2288 mytr.alpha = 0.0;
2289 mytr.cellSize = HGCalParameters::k_ScaleToDDD * php.waferSize_;
2290 for (auto const& dz : php.cellThickness_) {
2291 mytr.dz = 0.5 * HGCalParameters::k_ScaleToDDD * dz;
2292 php.fillModule(mytr, false);
2293 }
2294 for (unsigned k = 0; k < php.cellThickness_.size(); ++k) {
2295 HGCalParameters::hgtrap mytr = php.getModule(k, false);
2296 mytr.bl *= HGCalParameters::k_ScaleFromDDD;
2297 mytr.tl *= HGCalParameters::k_ScaleFromDDD;
2298 mytr.h *= HGCalParameters::k_ScaleFromDDD;
2299 mytr.dz *= HGCalParameters::k_ScaleFromDDD;
2300 mytr.cellSize *= HGCalParameters::k_ScaleFromDDD;
2301 php.fillModule(mytr, true);
2302 }
2303 #ifdef EDM_ML_DEBUG
2304 edm::LogVerbatim("HGCalGeom") << "HGCalGeomParameters: Total of " << php.waferCopy_.size() << " wafers";
2305 for (unsigned int k = 0; k < php.waferCopy_.size(); ++k) {
2306 int id = php.waferCopy_[k];
2307 edm::LogVerbatim("HGCalGeom") << "[" << k << "] " << std::hex << id << std::dec << ":"
2308 << HGCalWaferIndex::waferLayer(id) << ":" << HGCalWaferIndex::waferU(id) << ":"
2309 << HGCalWaferIndex::waferV(id) << " x " << php.waferPosX_[k] << " y "
2310 << php.waferPosY_[k] << " index " << php.wafersInLayers_[id];
2311 }
2312 edm::LogVerbatim("HGCalGeom") << "HGCalParameters: Total of " << php.waferTypeL_.size() << " wafer types";
2313 for (unsigned int k = 0; k < php.waferTypeL_.size(); ++k) {
2314 int id = indtypes[k];
2315 edm::LogVerbatim("HGCalGeom") << "[" << k << "] " << php.typesInLayers_[id] << ":" << php.waferTypeL_[k] << " ID "
2316 << std::hex << id << std::dec << ":" << HGCalWaferIndex::waferLayer(id) << ":"
2317 << HGCalWaferIndex::waferU(id) << ":" << HGCalWaferIndex::waferV(id);
2318 }
2319 #endif
2320
2321
2322 php.xLayerHex_.clear();
2323 php.yLayerHex_.clear();
2324 double waferSize = php.waferSize_ + php.sensorSeparation_;
2325 #ifdef EDM_ML_DEBUG
2326 edm::LogVerbatim("HGCalGeom") << "WaferSize " << waferSize;
2327 #endif
2328 for (unsigned int k = 0; k < php.zLayerHex_.size(); ++k) {
2329 int copy = k + php.layerOffset_;
2330 std::pair<double, double> xyoff = geomTools_.shiftXY(php.layerCenter_[copy], waferSize);
2331 php.xLayerHex_.emplace_back(xyoff.first);
2332 php.yLayerHex_.emplace_back(xyoff.second);
2333 #ifdef EDM_ML_DEBUG
2334 edm::LogVerbatim("HGCalGeom") << "Layer[" << k << "] Off " << copy << ":" << php.layerCenter_[copy] << " Shift "
2335 << xyoff.first << ":" << xyoff.second;
2336 #endif
2337 }
2338 }
2339
2340 void HGCalGeomParameters::loadCellParsHexagon(const DDCompactView* cpv, HGCalParameters& php) {
2341
2342 std::string attribute = "OnlyForHGCalNumbering";
2343 DDSpecificsHasNamedValueFilter filter1{attribute};
2344 DDFilteredView fv1(*cpv, filter1);
2345 bool ok = fv1.firstChild();
2346
2347 if (ok) {
2348 php.cellFine_ = dbl_to_int(cpv->vector("waferFine"));
2349 php.cellCoarse_ = dbl_to_int(cpv->vector("waferCoarse"));
2350 }
2351
2352 loadCellParsHexagon(php);
2353 }
2354
2355 void HGCalGeomParameters::loadCellParsHexagon(const cms::DDVectorsMap& vmap, HGCalParameters& php) {
2356 for (auto const& it : vmap) {
2357 if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "waferFine")) {
2358 for (const auto& i : it.second)
2359 php.cellFine_.emplace_back(std::round(i));
2360 } else if (dd4hep::dd::compareEqual(dd4hep::dd::noNamespace(it.first), "waferCoarse")) {
2361 for (const auto& i : it.second)
2362 php.cellCoarse_.emplace_back(std::round(i));
2363 }
2364 }
2365
2366 loadCellParsHexagon(php);
2367 }
2368
2369 void HGCalGeomParameters::loadCellParsHexagon(const HGCalParameters& php) {
2370 #ifdef EDM_ML_DEBUG
2371 edm::LogVerbatim("HGCalGeom") << "HGCalLoadCellPars: " << php.cellFine_.size() << " rows for fine cells";
2372 for (unsigned int k = 0; k < php.cellFine_.size(); ++k)
2373 edm::LogVerbatim("HGCalGeom") << "[" << k << "]: " << php.cellFine_[k];
2374 edm::LogVerbatim("HGCalGeom") << "HGCalLoadCellPars: " << php.cellCoarse_.size() << " rows for coarse cells";
2375 for (unsigned int k = 0; k < php.cellCoarse_.size(); ++k)
2376 edm::LogVerbatim("HGCalGeom") << "[" << k << "]: " << php.cellCoarse_[k];
2377 #endif
2378 }
2379
2380 void HGCalGeomParameters::loadCellTrapezoid(HGCalParameters& php) {
2381 php.xLayerHex_.resize(php.zLayerHex_.size(), 0);
2382 php.yLayerHex_.resize(php.zLayerHex_.size(), 0);
2383 #ifdef EDM_ML_DEBUG
2384 edm::LogVerbatim("HGCalGeom") << "HGCalParameters: x|y|zLayerHex in array of size " << php.zLayerHex_.size();
2385 for (unsigned int k = 0; k < php.zLayerHex_.size(); ++k)
2386 edm::LogVerbatim("HGCalGeom") << "Layer[" << k << "] Shift " << php.xLayerHex_[k] << ":" << php.yLayerHex_[k] << ":"
2387 << php.zLayerHex_[k];
2388 #endif
2389
2390 if ((php.mode_ == HGCalGeometryMode::TrapezoidFile) || (php.mode_ == HGCalGeometryMode::TrapezoidModule) ||
2391 (php.mode_ == HGCalGeometryMode::TrapezoidCassette) || (php.mode_ == HGCalGeometryMode::TrapezoidFineCell)) {
2392
2393 #ifdef EDM_ML_DEBUG
2394 edm::LogVerbatim("HGCalGeom") << "HGCalParameters: Mode " << php.mode_ << ":"
2395 << HGCalGeometryMode::TrapezoidFineCell << " Sizes " << php.tileRingFineR_.size()
2396 << ":" << php.tileRingR_.size();
2397 #endif
2398 for (unsigned int k = 0; k < 2; ++k) {
2399 bool fine = ((k == 1) && (php.mode_ == HGCalGeometryMode::TrapezoidFineCell));
2400 unsigned int sizeR = (fine) ? php.tileRingFineR_.size() : php.tileRingR_.size();
2401 for (unsigned int kk = 0; kk < sizeR; ++kk) {
2402 if (fine)
2403 php.radiusLayer_[k].emplace_back(php.tileRingFineR_[kk].first);
2404 else
2405 php.radiusLayer_[k].emplace_back(php.tileRingR_[kk].first);
2406 #ifdef EDM_ML_DEBUG
2407 double zv = ((k == 0) ? (php.zLayerHex_[php.layerFrontBH_[1] - php.firstLayer_])
2408 : (php.zLayerHex_[php.zLayerHex_.size() - 1]));
2409 double rv = php.radiusLayer_[k].back();
2410 double eta = -(std::log(std::tan(0.5 * std::atan(rv / zv))));
2411 edm::LogVerbatim("HGCalGeom") << "New [" << kk << "] new R = " << rv << " Eta = " << eta;
2412 #endif
2413 }
2414 if (fine)
2415 php.radiusLayer_[k].emplace_back(php.tileRingFineR_[php.tileRingFineR_.size() - 1].second);
2416 else
2417 php.radiusLayer_[k].emplace_back(php.tileRingR_[php.tileRingR_.size() - 1].second);
2418 }
2419
2420 #ifdef EDM_ML_DEBUG
2421 edm::LogVerbatim("HGCalGeom") << "Tile Ring Range size " << php.tileRingFineRange_.size() << ":"
2422 << php.tileRingRange_.size() << ":" << php.nPhiBinBH_.size() << ":"
2423 << php.zLayerHex_.size() << ":" << php.nPhiLayer_.size();
2424 #endif
2425 unsigned int k1(0), k2(0);
2426 for (unsigned int k = 0; k < php.zLayerHex_.size(); ++k) {
2427 if (!php.tileRingFineRange_.empty()) {
2428 unsigned int k0 = (k1 < (php.tileRingFineRange_.size() - 1)) ? k1 : (php.tileRingFineRange_.size() - 1);
2429 #ifdef EDM_ML_DEBUG
2430 edm::LogVerbatim("HGCalGeom") << "Layer " << k << " Fine " << k1 << ":" << k0 << ":"
2431 << php.tileRingFineRange_.size();
2432 #endif
2433 php.iradMinBHFine_.emplace_back(1 + php.tileRingFineRange_[k0].first);
2434 php.iradMaxBHFine_.emplace_back(1 + php.tileRingFineRange_[k0].second);
2435 }
2436 if (!php.tileRingRange_.empty()) {
2437 unsigned int k0 = (k2 < (php.tileRingRange_.size() - 1)) ? k2 : (php.tileRingRange_.size() - 1);
2438 #ifdef EDM_ML_DEBUG
2439 edm::LogVerbatim("HGCalGeom") << "Layer " << k << " Coarse " << k2 << ":" << k0 << ":"
2440 << php.tileRingRange_.size();
2441 #endif
2442 php.iradMinBH_.emplace_back(1 + php.tileRingRange_[k0].first);
2443 php.iradMaxBH_.emplace_back(1 + php.tileRingRange_[k0].second);
2444 }
2445 if (php.nPhiLayer_[k] > 288) {
2446 ++k1;
2447 } else {
2448 ++k2;
2449 }
2450 #ifdef EDM_ML_DEBUG
2451 int kk = php.scintType(php.firstLayer_ + static_cast<int>(k));
2452 if (php.nPhiLayer_[k] > 288)
2453 edm::LogVerbatim("HGCalGeom") << "New Layer " << k << " Type " << kk << " Low edge "
2454 << php.iradMinBHFine_.back() << " Top edge " << php.iradMaxBHFine_.back();
2455 else
2456 edm::LogVerbatim("HGCalGeom") << "New Layer " << k << " Type " << kk << " Low edge " << php.iradMinBH_.back()
2457 << " Top edge " << php.iradMaxBH_.back();
2458 #endif
2459 }
2460 } else {
2461
2462 for (unsigned int k = 0; k < 2; ++k) {
2463 double rmax = ((k == 0) ? (php.rMaxLayHex_[php.layerFrontBH_[1] - php.firstLayer_] - 1)
2464 : (php.rMaxLayHex_[php.rMaxLayHex_.size() - 1]));
2465 double rv = php.rMinLayerBH_[k];
2466 double zv = ((k == 0) ? (php.zLayerHex_[php.layerFrontBH_[1] - php.firstLayer_])
2467 : (php.zLayerHex_[php.zLayerHex_.size() - 1]));
2468 php.radiusLayer_[k].emplace_back(rv);
2469 #ifdef EDM_ML_DEBUG
2470 double eta = -(std::log(std::tan(0.5 * std::atan(rv / zv))));
2471 edm::LogVerbatim("HGCalGeom") << "Old [" << k << "] rmax " << rmax << " Z = " << zv
2472 << " dEta = " << php.cellSize_[k] << "\nOld[0] new R = " << rv << " Eta = " << eta;
2473 int kount(1);
2474 #endif
2475 while (rv < rmax) {
2476 double eta = -(php.cellSize_[k] + std::log(std::tan(0.5 * std::atan(rv / zv))));
2477 rv = zv * std::tan(2.0 * std::atan(std::exp(-eta)));
2478 php.radiusLayer_[k].emplace_back(rv);
2479 #ifdef EDM_ML_DEBUG
2480 edm::LogVerbatim("HGCalGeom") << "Old [" << kount << "] new R = " << rv << " Eta = " << eta;
2481 ++kount;
2482 #endif
2483 }
2484 }
2485
2486 for (unsigned int k = 0; k < php.zLayerHex_.size(); ++k) {
2487 int kk = php.scintType(php.firstLayer_ + static_cast<int>(k));
2488 std::vector<double>::iterator low, high;
2489 low = std::lower_bound(php.radiusLayer_[kk].begin(), php.radiusLayer_[kk].end(), php.rMinLayHex_[k]);
2490 #ifdef EDM_ML_DEBUG
2491 edm::LogVerbatim("HGCalGeom") << "Old [" << k << "] RLow = " << php.rMinLayHex_[k] << " pos "
2492 << static_cast<int>(low - php.radiusLayer_[kk].begin());
2493 #endif
2494 if (low == php.radiusLayer_[kk].begin())
2495 ++low;
2496 int irlow = static_cast<int>(low - php.radiusLayer_[kk].begin());
2497 double drlow = php.radiusLayer_[kk][irlow] - php.rMinLayHex_[k];
2498 #ifdef EDM_ML_DEBUG
2499 edm::LogVerbatim("HGCalGeom") << "irlow " << irlow << " dr " << drlow << " min " << php.minTileSize_;
2500 #endif
2501 if (drlow < php.minTileSize_) {
2502 ++irlow;
2503 #ifdef EDM_ML_DEBUG
2504 drlow = php.radiusLayer_[kk][irlow] - php.rMinLayHex_[k];
2505 edm::LogVerbatim("HGCalGeom") << "Modified irlow " << irlow << " dr " << drlow;
2506 #endif
2507 }
2508 high = std::lower_bound(php.radiusLayer_[kk].begin(), php.radiusLayer_[kk].end(), php.rMaxLayHex_[k]);
2509 #ifdef EDM_ML_DEBUG
2510 edm::LogVerbatim("HGCalGeom") << "Old [" << k << "] RHigh = " << php.rMaxLayHex_[k] << " pos "
2511 << static_cast<int>(high - php.radiusLayer_[kk].begin());
2512 #endif
2513 if (high == php.radiusLayer_[kk].end())
2514 --high;
2515 int irhigh = static_cast<int>(high - php.radiusLayer_[kk].begin());
2516 double drhigh = php.rMaxLayHex_[k] - php.radiusLayer_[kk][irhigh - 1];
2517 #ifdef EDM_ML_DEBUG
2518 edm::LogVerbatim("HGCalGeom") << "irhigh " << irhigh << " dr " << drhigh << " min " << php.minTileSize_;
2519 #endif
2520 if (drhigh < php.minTileSize_) {
2521 --irhigh;
2522 #ifdef EDM_ML_DEBUG
2523 drhigh = php.rMaxLayHex_[k] - php.radiusLayer_[kk][irhigh - 1];
2524 edm::LogVerbatim("HGCalGeom") << "Modified irhigh " << irhigh << " dr " << drhigh;
2525 #endif
2526 }
2527 php.iradMinBHFine_.emplace_back(irlow);
2528 php.iradMaxBHFine_.emplace_back(irhigh);
2529 php.iradMinBH_.emplace_back(irlow);
2530 php.iradMaxBH_.emplace_back(irhigh);
2531 #ifdef EDM_ML_DEBUG
2532 edm::LogVerbatim("HGCalGeom") << "Old Layer " << k << " Type " << kk << " Low edge " << irlow << ":" << drlow
2533 << " Top edge " << irhigh << ":" << drhigh;
2534 #endif
2535 }
2536 }
2537 #ifdef EDM_ML_DEBUG
2538 for (unsigned int k = 0; k < 2; ++k) {
2539 edm::LogVerbatim("HGCalGeom") << "Type " << k << " with " << php.radiusLayer_[k].size() << " radii";
2540 for (unsigned int kk = 0; kk < php.radiusLayer_[k].size(); ++kk)
2541 edm::LogVerbatim("HGCalGeom") << "Ring[" << k << "][" << kk << "] " << php.radiusLayer_[k][kk];
2542 }
2543 #endif
2544
2545
2546 int im(0);
2547 php.tileUVMax_ = php.tileUVMaxFine_ = 00;
2548 HGCalParameters::hgtrap mytr;
2549 mytr.alpha = 0.0;
2550 for (unsigned int k = 0; k < php.zLayerHex_.size(); ++k) {
2551 if (php.nPhiLayer_[k] > 288) {
2552 if (php.iradMaxBHFine_[k] > php.tileUVMaxFine_)
2553 php.tileUVMaxFine_ = php.iradMaxBHFine_[k];
2554 } else {
2555 if (php.iradMaxBH_[k] > php.tileUVMax_)
2556 php.tileUVMax_ = php.iradMaxBH_[k];
2557 }
2558 int kk = (php.nPhiLayer_[k] > 288) ? 1 : 0;
2559 int irm = php.radiusLayer_[kk].size() - 1;
2560 int irmin = (php.nPhiLayer_[k] > 288) ? php.iradMinBHFine_[k] : php.iradMinBH_[k];
2561 int irmax = (php.nPhiLayer_[k] > 288) ? php.iradMaxBHFine_[k] : php.iradMaxBH_[k];
2562 #ifdef EDM_ML_DEBUG
2563 double rmin = php.radiusLayer_[kk][std::max((irmin - 1), 0)];
2564 double rmax = php.radiusLayer_[kk][std::min(irmax, irm)];
2565 edm::LogVerbatim("HGCalGeom") << "Layer " << php.firstLayer_ + k << ":" << kk << " Radius range " << irmin << ":"
2566 << irmax << ":" << rmin << ":" << rmax << " Size " << php.radiusLayer_[kk].size();
2567 #endif
2568 mytr.lay = php.firstLayer_ + k;
2569 for (int irad = irmin; irad <= irmax; ++irad) {
2570 double rmin = php.radiusLayer_[kk][std::max((irad - 1), 0)];
2571 double rmax = php.radiusLayer_[kk][std::min(irad, irm)];
2572 #ifdef EDM_ML_DEBUG
2573 edm::LogVerbatim("HGCalGeom") << "irad " << irad << ":" << php.radiusLayer_[kk].size() << " R " << rmin << ":"
2574 << rmax;
2575 #endif
2576 mytr.bl = 0.5 * rmin * php.scintCellSize(mytr.lay);
2577 mytr.tl = 0.5 * rmax * php.scintCellSize(mytr.lay);
2578 mytr.h = 0.5 * (rmax - rmin);
2579 mytr.dz = 0.5 * php.waferThick_;
2580 mytr.cellSize = 0.5 * (rmax + rmin) * php.scintCellSize(mytr.lay);
2581 php.fillModule(mytr, true);
2582 mytr.bl *= HGCalParameters::k_ScaleToDDD;
2583 mytr.tl *= HGCalParameters::k_ScaleToDDD;
2584 mytr.h *= HGCalParameters::k_ScaleToDDD;
2585 mytr.dz *= HGCalParameters::k_ScaleToDDD;
2586 mytr.cellSize *= HGCalParameters::k_ScaleFromDDD;
2587 php.fillModule(mytr, false);
2588 if (irad == irmin)
2589 php.firstModule_.emplace_back(im);
2590 ++im;
2591 if (irad == irmax)
2592 php.lastModule_.emplace_back(im);
2593 }
2594 }
2595 php.nSectors_ = std::max(php.tileUVMax_, php.tileUVMaxFine_);
2596 #ifdef EDM_ML_DEBUG
2597 edm::LogVerbatim("HGCalGeom") << "Maximum radius index " << php.tileUVMax_ << ":" << php.tileUVMaxFine_;
2598 for (unsigned int k = 0; k < php.firstModule_.size(); ++k)
2599 edm::LogVerbatim("HGCalGeom") << "Layer " << k + php.firstLayer_ << " Modules " << php.firstModule_[k] << ":"
2600 << php.lastModule_[k];
2601
2602 #endif
2603 }
2604
2605 std::vector<double> HGCalGeomParameters::getDDDArray(const std::string& str, const DDsvalues_type& sv, const int nmin) {
2606 DDValue value(str);
2607 if (DDfetch(&sv, value)) {
2608 const std::vector<double>& fvec = value.doubles();
2609 int nval = fvec.size();
2610 if (nmin > 0) {
2611 if (nval < nmin) {
2612 throw cms::Exception("DDException")
2613 << "HGCalGeomParameters: # of " << str << " bins " << nval << " < " << nmin << " ==> illegal";
2614 }
2615 } else {
2616 if (nval < 1 && nmin == 0) {
2617 throw cms::Exception("DDException")
2618 << "HGCalGeomParameters: # of " << str << " bins " << nval << " < 1 ==> illegal"
2619 << " (nmin=" << nmin << ")";
2620 }
2621 }
2622 return fvec;
2623 } else {
2624 if (nmin >= 0) {
2625 throw cms::Exception("DDException") << "HGCalGeomParameters: cannot get array " << str;
2626 }
2627 std::vector<double> fvec;
2628 return fvec;
2629 }
2630 }
2631
2632 std::pair<double, double> HGCalGeomParameters::cellPosition(
2633 const std::vector<HGCalGeomParameters::cellParameters>& wafers,
2634 std::vector<HGCalGeomParameters::cellParameters>::const_iterator& itrf,
2635 int wafer,
2636 double xx,
2637 double yy) {
2638 if (itrf == wafers.end()) {
2639 for (std::vector<HGCalGeomParameters::cellParameters>::const_iterator itr = wafers.begin(); itr != wafers.end();
2640 ++itr) {
2641 if (itr->wafer == wafer) {
2642 itrf = itr;
2643 break;
2644 }
2645 }
2646 }
2647 double dx(0), dy(0);
2648 if (itrf != wafers.end()) {
2649 dx = (xx - itrf->xyz.x());
2650 if (std::abs(dx) < tolerance)
2651 dx = 0;
2652 dy = (yy - itrf->xyz.y());
2653 if (std::abs(dy) < tolerance)
2654 dy = 0;
2655 }
2656 return std::make_pair(dx, dy);
2657 }
2658
2659 void HGCalGeomParameters::rescale(std::vector<double>& v, const double s) {
2660 std::for_each(v.begin(), v.end(), [s](double& n) { n *= s; });
2661 }
2662
2663 void HGCalGeomParameters::resetZero(std::vector<double>& v) {
2664 for (auto& n : v) {
2665 if (std::abs(n) < tolmin)
2666 n = 0;
2667 }
2668 }