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