File indexing completed on 2024-04-06 12:14:56
0001 #include "DataFormats/Math/interface/angle_units.h"
0002 #include "DD4hep/DetFactoryHelper.h"
0003 #include "DetectorDescription/DDCMS/interface/DDPlugins.h"
0004 #include "DetectorDescription/DDCMS/interface/DDutils.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "FWCore/Utilities/interface/Exception.h"
0007 #include "Geometry/HGCalCommonData/interface/HGCalGeomTools.h"
0008 #include "Geometry/HGCalCommonData/interface/HGCalParameters.h"
0009 #include "Geometry/HGCalCommonData/interface/HGCalTypes.h"
0010
0011
0012 #ifdef EDM_ML_DEBUG
0013 #include <unordered_set>
0014 #endif
0015 using namespace angle_units::operators;
0016
0017 static long algorithm(dd4hep::Detector& , cms::DDParsingContext& ctxt, xml_h e) {
0018 cms::DDNamespace ns(ctxt, e, true);
0019 cms::DDAlgoArguments args(ctxt, e);
0020 static constexpr double tol2 = 0.00001 * dd4hep::mm;
0021
0022 const auto& wafers = args.value<std::vector<std::string> >("WaferName");
0023 const auto& covers = args.value<std::vector<std::string> >("CoverName");
0024 #ifdef EDM_ML_DEBUG
0025 edm::LogVerbatim("HGCalGeom") << "DDHGCalTBModule: " << wafers.size() << " wafers";
0026 unsigned int i(0);
0027 for (auto wafer : wafers) {
0028 edm::LogVerbatim("HGCalGeom") << "Wafer[" << i << "] " << wafer;
0029 ++i;
0030 }
0031 edm::LogVerbatim("HGCalGeom") << "DDHGCalTBModule: " << covers.size() << " covers";
0032 i = 0;
0033 for (auto cover : covers) {
0034 edm::LogVerbatim("HGCalGeom") << "Cover[" << i << "] " << cover;
0035 ++i;
0036 }
0037 #endif
0038 const auto& materials = args.value<std::vector<std::string> >("MaterialNames");
0039 const auto& names = args.value<std::vector<std::string> >("VolumeNames");
0040 const auto& thick = args.value<std::vector<double> >("Thickness");
0041 std::vector<int> copyNumber;
0042 copyNumber.resize(materials.size(), 1);
0043 #ifdef EDM_ML_DEBUG
0044 edm::LogVerbatim("HGCalGeom") << "DDHGCalTBModule: " << materials.size() << " types of volumes";
0045 for (unsigned int i = 0; i < names.size(); ++i)
0046 edm::LogVerbatim("HGCalGeom") << "Volume [" << i << "] " << names[i] << " of thickness "
0047 << cms::convert2mm(thick[i]) << " filled with " << materials[i]
0048 << " first copy number " << copyNumber[i];
0049 #endif
0050 const auto& layers = args.value<std::vector<int> >("Layers");
0051 const auto& layerThick = args.value<std::vector<double> >("LayerThick");
0052 #ifdef EDM_ML_DEBUG
0053 edm::LogVerbatim("HGCalGeom") << "DDHGCalTBModule: " << layers.size() << " blocks";
0054 for (unsigned int i = 0; i < layers.size(); ++i)
0055 edm::LogVerbatim("HGCalGeom") << "Block [" << i << "] of thickness " << cms::convert2mm(layerThick[i]) << " with "
0056 << layers[i] << " layers";
0057 #endif
0058 const auto& layerType = args.value<std::vector<int> >("LayerType");
0059 const auto& layerSense = args.value<std::vector<int> >("LayerSense");
0060 #ifdef EDM_ML_DEBUG
0061 edm::LogVerbatim("HGCalGeom") << "DDHGCalTBModule: " << layerType.size() << " layers";
0062 for (unsigned int i = 0; i < layerType.size(); ++i)
0063 edm::LogVerbatim("HGCalGeom") << "Layer [" << i << "] with material type " << layerType[i] << " sensitive class "
0064 << layerSense[i];
0065 #endif
0066 const auto& zMinBlock = args.value<double>("zMinBlock");
0067 const auto& rMaxFine = args.value<double>("rMaxFine");
0068 const auto& waferW = args.value<double>("waferW");
0069 const auto& waferGap = args.value<double>("waferGap");
0070 const auto& absorbW = args.value<double>("absorberW");
0071 const auto& absorbH = args.value<double>("absorberH");
0072 #ifdef EDM_ML_DEBUG
0073 edm::LogVerbatim("HGCalGeom") << "DDHGCalTBModule: zStart " << cms::convert2mm(zMinBlock) << " rFineCoarse "
0074 << cms::convert2mm(rMaxFine) << " wafer width " << cms::convert2mm(waferW)
0075 << " gap among wafers " << cms::convert2mm(waferGap) << " absorber width "
0076 << cms::convert2mm(absorbW) << " absorber height " << cms::convert2mm(absorbH);
0077 #endif
0078 const auto& slopeB = args.value<std::vector<double> >("SlopeBottom");
0079 const auto& slopeT = args.value<std::vector<double> >("SlopeTop");
0080 const auto& zFront = args.value<std::vector<double> >("ZFront");
0081 const auto& rMaxFront = args.value<std::vector<double> >("RMaxFront");
0082 #ifdef EDM_ML_DEBUG
0083 edm::LogVerbatim("HGCalGeom") << "DDHGCalTBModule: Bottom slopes " << slopeB[0] << ":" << slopeB[1] << " and "
0084 << slopeT.size() << " slopes for top";
0085 for (unsigned int i = 0; i < slopeT.size(); ++i)
0086 edm::LogVerbatim("HGCalGeom") << "Block [" << i << "] Zmin " << cms::convert2mm(zFront[i]) << " Rmax "
0087 << cms::convert2mm(rMaxFront[i]) << " Slope " << slopeT[i];
0088 #endif
0089 std::string idNameSpace = static_cast<std::string>(ns.name());
0090 const auto& idName = args.parentName();
0091 #ifdef EDM_ML_DEBUG
0092 std::unordered_set<int> copies;
0093 edm::LogVerbatim("HGCalGeom") << "DDHGCalTBModule: NameSpace " << idNameSpace << " Mother " << idName;
0094 #endif
0095
0096
0097 dd4hep::Volume module = ns.volume(idName);
0098 #ifdef EDM_ML_DEBUG
0099 edm::LogVerbatim("HGCalGeom") << "==>> Constructing DDHGCalTBModule...";
0100 #endif
0101
0102 double zi(zMinBlock);
0103 double ww = (waferW + waferGap);
0104 double dx = 0.5 * ww;
0105 double dy = 3.0 * dx * tan(30._deg);
0106 double rr = 2.0 * dx * tan(30._deg);
0107 int laymin(0);
0108 for (unsigned int i = 0; i < layers.size(); i++) {
0109 double zo = zi + layerThick[i];
0110 double routF = HGCalGeomTools::radius(zi, zFront, rMaxFront, slopeT);
0111 int laymax = laymin + layers[i];
0112 double zz = zi;
0113 double thickTot(0);
0114 for (int ly = laymin; ly < laymax; ++ly) {
0115 int ii = layerType[ly];
0116 int copy = copyNumber[ii];
0117 double rinB = (layerSense[ly] == 0) ? (zo * slopeB[0]) : (zo * slopeB[1]);
0118 zz += (0.5 * thick[ii]);
0119 thickTot += thick[ii];
0120
0121 std::string name = "HGCal" + names[ii] + std::to_string(copy);
0122 #ifdef EDM_ML_DEBUG
0123 edm::LogVerbatim("HGCalGeom") << "DDHGCalTBModule: Layer " << ly << ":" << ii << " Front " << cms::convert2mm(zi)
0124 << ", " << cms::convert2mm(routF) << " Back " << cms::convert2mm(zo) << ", "
0125 << cms::convert2mm(rinB) << " superlayer thickness "
0126 << cms::convert2mm(layerThick[i]);
0127 #endif
0128
0129 dd4hep::Material matter = ns.material(materials[ii]);
0130 dd4hep::Volume glog;
0131 if (layerSense[ly] == 0) {
0132 dd4hep::Solid solid = dd4hep::Box(absorbW, absorbH, 0.5 * thick[ii]);
0133 ns.addSolidNS(ns.prepend(name), solid);
0134 glog = dd4hep::Volume(solid.name(), solid, matter);
0135 #ifdef EDM_ML_DEBUG
0136 edm::LogVerbatim("HGCalGeom") << "DDHGCalTBModule test: " << solid.name() << " box of dimension "
0137 << cms::convert2mm(absorbW) << ":" << cms::convert2mm(absorbH) << ":"
0138 << cms::convert2mm(0.5 * thick[ii]);
0139 #endif
0140 } else {
0141 dd4hep::Solid solid = dd4hep::Tube(rinB, routF, 0.5 * thick[ii], 0.0, 2._pi);
0142 ns.addSolidNS(ns.prepend(name), solid);
0143 glog = dd4hep::Volume(solid.name(), solid, matter);
0144 #ifdef EDM_ML_DEBUG
0145 edm::LogVerbatim("HGCalGeom") << "DDHGCalTBModule: " << solid.name() << " Tubs made of " << materials[ii]
0146 << " of dimensions " << cms::convert2mm(rinB) << ", " << cms::convert2mm(routF)
0147 << ", " << cms::convert2mm(0.5 * thick[ii]) << ", 0.0, 360.0";
0148 #endif
0149 int ncol = static_cast<int>(2.0 * routF / ww) + 1;
0150 int nrow = static_cast<int>(routF / (ww * tan(30._deg))) + 1;
0151 #ifdef EDM_ML_DEBUG
0152 int incm(0), inrm(0), kount(0), ntot(0), nin(0), nfine(0), ncoarse(0);
0153 edm::LogVerbatim("HGCalGeom") << glog.name() << " rout " << cms::convert2mm(routF) << " Row " << nrow
0154 << " Column " << ncol;
0155 #endif
0156 double xc[6], yc[6];
0157 for (int nr = -nrow; nr <= nrow; ++nr) {
0158 int inr = (nr >= 0) ? nr : -nr;
0159 for (int nc = -ncol; nc <= ncol; ++nc) {
0160 int inc = (nc >= 0) ? nc : -nc;
0161 if (inr % 2 == inc % 2) {
0162 double xpos = nc * dx;
0163 double ypos = nr * dy;
0164 xc[0] = xpos + dx;
0165 yc[0] = ypos - 0.5 * rr;
0166 xc[1] = xpos + dx;
0167 yc[1] = ypos + 0.5 * rr;
0168 xc[2] = xpos;
0169 yc[2] = ypos + rr;
0170 xc[3] = xpos - dx;
0171 yc[3] = ypos + 0.5 * rr;
0172 xc[4] = xpos + dx;
0173 yc[4] = ypos - 0.5 * rr;
0174 xc[5] = xpos;
0175 yc[5] = ypos - rr;
0176 bool cornerAll(true);
0177 for (int k = 0; k < 6; ++k) {
0178 double rpos = std::sqrt(xc[k] * xc[k] + yc[k] * yc[k]);
0179 if (rpos < rinB || rpos > routF)
0180 cornerAll = false;
0181 }
0182 #ifdef EDM_ML_DEBUG
0183 ++ntot;
0184 #endif
0185 if (cornerAll) {
0186 dd4hep::Volume glog1;
0187 if (layerSense[ly] == 1) {
0188 double rpos = std::sqrt(xpos * xpos + ypos * ypos);
0189 glog1 = (rpos < rMaxFine) ? ns.volume(wafers[0]) : ns.volume(wafers[1]);
0190 #ifdef EDM_ML_DEBUG
0191 ++nin;
0192 if (rpos < rMaxFine)
0193 ++nfine;
0194 else
0195 ++ncoarse;
0196 #endif
0197 } else {
0198 glog1 = ns.volume(covers[layerSense[ly] - 2]);
0199 }
0200 int copyL = HGCalTypes::packTypeUV(0, nc, nr);
0201 #ifdef EDM_ML_DEBUG
0202 if (inc > incm)
0203 incm = inc;
0204 if (inr > inrm)
0205 inrm = inr;
0206 kount++;
0207 copies.insert(copy);
0208 #endif
0209 dd4hep::Position tran(xpos, ypos, 0.0);
0210 glog.placeVolume(glog1, copyL, tran);
0211 #ifdef EDM_ML_DEBUG
0212 edm::LogVerbatim("HGCalGeom")
0213 << "DDHGCalModule: " << glog1.name() << " number " << copyL << " positioned in " << glog.name()
0214 << " at (" << cms::convert2mm(xpos) << "," << cms::convert2mm(ypos) << ",0)";
0215 #endif
0216 }
0217 }
0218 }
0219 }
0220 #ifdef EDM_ML_DEBUG
0221 edm::LogVerbatim("HGCalGeom") << "DDHGCalModule: # of columns " << incm << " # of rows " << inrm << " and "
0222 << nin << ":" << kount << ":" << ntot << " wafers (" << nfine << ":" << ncoarse
0223 << ") for " << glog.name() << " R " << cms::convert2mm(rinB) << ":"
0224 << cms::convert2mm(routF);
0225 #endif
0226 }
0227 dd4hep::Position r1(0, 0, zz);
0228 module.placeVolume(glog, copy, r1);
0229 ++copyNumber[ii];
0230 #ifdef EDM_ML_DEBUG
0231 edm::LogVerbatim("HGCalGeom") << "DDHGCalTBModule test: " << glog.name() << " number " << copy
0232 << " positioned in " << module.name() << " at (0,0," << cms::convert2mm(zz) << ")";
0233 #endif
0234 zz += (0.5 * thick[ii]);
0235 }
0236 zi = zo;
0237 laymin = laymax;
0238 if (fabs(thickTot - layerThick[i]) > tol2) {
0239 if (thickTot > layerThick[i]) {
0240 edm::LogError("HGCalGeom") << "Thickness of the partition " << cms::convert2mm(layerThick[i])
0241 << " is smaller than thickness " << cms::convert2mm(thickTot)
0242 << " of all its components **** ERROR ****\n";
0243 } else {
0244 edm::LogWarning("HGCalGeom") << "Thickness of the partition " << cms::convert2mm(layerThick[i])
0245 << " does not match with " << cms::convert2mm(thickTot) << " of the components\n";
0246 }
0247 }
0248 }
0249
0250 return cms::s_executed;
0251 }
0252
0253
0254 DECLARE_DDCMS_DETELEMENT(DDCMS_hgcal_DDHGCalTBModule, algorithm)