File indexing completed on 2024-04-06 12:14:58
0001
0002
0003
0004
0005
0006 #include "DataFormats/Math/interface/angle_units.h"
0007 #include "DetectorDescription/Core/interface/DDAlgorithm.h"
0008 #include "DetectorDescription/Core/interface/DDAlgorithmFactory.h"
0009 #include "DetectorDescription/Core/interface/DDCurrentNamespace.h"
0010 #include "DetectorDescription/Core/interface/DDLogicalPart.h"
0011 #include "DetectorDescription/Core/interface/DDMaterial.h"
0012 #include "DetectorDescription/Core/interface/DDSolid.h"
0013 #include "DetectorDescription/Core/interface/DDSplit.h"
0014 #include "DetectorDescription/Core/interface/DDTypes.h"
0015 #include "DetectorDescription/Core/interface/DDutils.h"
0016 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0017 #include "FWCore/PluginManager/interface/PluginFactory.h"
0018 #include "Geometry/HGCalCommonData/interface/HGCalGeomTools.h"
0019 #include "Geometry/HGCalCommonData/interface/HGCalParameters.h"
0020 #include "Geometry/HGCalCommonData/interface/HGCalTypes.h"
0021 #include "Geometry/HGCalCommonData/interface/HGCalWaferType.h"
0022
0023 #include <cmath>
0024 #include <memory>
0025 #include <string>
0026 #include <unordered_set>
0027 #include <vector>
0028
0029
0030 using namespace angle_units::operators;
0031
0032 class DDHGCalEEAlgo : public DDAlgorithm {
0033 public:
0034
0035 DDHGCalEEAlgo();
0036 ~DDHGCalEEAlgo() override;
0037
0038 void initialize(const DDNumericArguments& nArgs,
0039 const DDVectorArguments& vArgs,
0040 const DDMapArguments& mArgs,
0041 const DDStringArguments& sArgs,
0042 const DDStringVectorArguments& vsArgs) override;
0043 void execute(DDCompactView& cpv) override;
0044
0045 protected:
0046 void constructLayers(const DDLogicalPart&, DDCompactView& cpv);
0047 void positionSensitive(const DDLogicalPart& glog,
0048 double rin,
0049 double rout,
0050 double zpos,
0051 int layertype,
0052 int layercenter,
0053 DDCompactView& cpv);
0054
0055 private:
0056 HGCalGeomTools geomTools_;
0057 std::unique_ptr<HGCalWaferType> waferType_;
0058
0059 static constexpr double tol1_ = 0.01;
0060 static constexpr double tol2_ = 0.00001;
0061
0062 std::vector<std::string> wafers_;
0063 std::vector<std::string> materials_;
0064 std::vector<std::string> names_;
0065 std::vector<double> thick_;
0066 std::vector<int> copyNumber_;
0067 std::vector<int> layers_;
0068 std::vector<double> layerThick_;
0069 std::vector<int> layerType_;
0070 std::vector<int> layerSense_;
0071 std::vector<int> layerCenter_;
0072 int firstLayer_;
0073 int absorbMode_;
0074 int sensitiveMode_;
0075 double zMinBlock_;
0076 std::vector<double> rad100to200_;
0077 std::vector<double> rad200to300_;
0078 double zMinRadPar_;
0079 int choiceType_;
0080 int nCutRadPar_;
0081 double fracAreaMin_;
0082 double waferSize_;
0083 double waferSepar_;
0084 int sectors_;
0085 std::vector<double> slopeB_;
0086 std::vector<double> zFrontB_;
0087 std::vector<double> rMinFront_;
0088 std::vector<double> slopeT_;
0089 std::vector<double> zFrontT_;
0090 std::vector<double> rMaxFront_;
0091 std::string nameSpace_;
0092 std::unordered_set<int> copies_;
0093 double alpha_, cosAlpha_;
0094 };
0095
0096 DDHGCalEEAlgo::DDHGCalEEAlgo() {
0097 #ifdef EDM_ML_DEBUG
0098 edm::LogVerbatim("HGCalGeom") << "DDHGCalEEAlgo: Creating an instance";
0099 #endif
0100 }
0101
0102 DDHGCalEEAlgo::~DDHGCalEEAlgo() {}
0103
0104 void DDHGCalEEAlgo::initialize(const DDNumericArguments& nArgs,
0105 const DDVectorArguments& vArgs,
0106 const DDMapArguments&,
0107 const DDStringArguments& sArgs,
0108 const DDStringVectorArguments& vsArgs) {
0109 wafers_ = vsArgs["WaferNames"];
0110 #ifdef EDM_ML_DEBUG
0111 edm::LogVerbatim("HGCalGeom") << "DDHGCalEEAlgo: " << wafers_.size() << " wafers";
0112 for (unsigned int i = 0; i < wafers_.size(); ++i)
0113 edm::LogVerbatim("HGCalGeom") << "Wafer[" << i << "] " << wafers_[i];
0114 #endif
0115 materials_ = vsArgs["MaterialNames"];
0116 names_ = vsArgs["VolumeNames"];
0117 thick_ = vArgs["Thickness"];
0118 copyNumber_.resize(materials_.size(), 1);
0119 #ifdef EDM_ML_DEBUG
0120 edm::LogVerbatim("HGCalGeom") << "DDHGCalEEAlgo: " << materials_.size() << " types of volumes";
0121 for (unsigned int i = 0; i < names_.size(); ++i)
0122 edm::LogVerbatim("HGCalGeom") << "Volume [" << i << "] " << names_[i] << " of thickness " << thick_[i]
0123 << " filled with " << materials_[i] << " first copy number " << copyNumber_[i];
0124 #endif
0125 layers_ = dbl_to_int(vArgs["Layers"]);
0126 layerThick_ = vArgs["LayerThick"];
0127 #ifdef EDM_ML_DEBUG
0128 edm::LogVerbatim("HGCalGeom") << "There are " << layers_.size() << " blocks";
0129 for (unsigned int i = 0; i < layers_.size(); ++i)
0130 edm::LogVerbatim("HGCalGeom") << "Block [" << i << "] of thickness " << layerThick_[i] << " with " << layers_[i]
0131 << " layers";
0132 #endif
0133 layerType_ = dbl_to_int(vArgs["LayerType"]);
0134 layerSense_ = dbl_to_int(vArgs["LayerSense"]);
0135 firstLayer_ = (int)(nArgs["FirstLayer"]);
0136 absorbMode_ = (int)(nArgs["AbsorberMode"]);
0137 sensitiveMode_ = (int)(nArgs["SensitiveMode"]);
0138 #ifdef EDM_ML_DEBUG
0139 edm::LogVerbatim("HGCalGeom") << "First Layer " << firstLayer_ << " and "
0140 << "Absober:Sensitive mode " << absorbMode_ << ":" << sensitiveMode_;
0141 #endif
0142 layerCenter_ = dbl_to_int(vArgs["LayerCenter"]);
0143 #ifdef EDM_ML_DEBUG
0144 for (unsigned int i = 0; i < layerCenter_.size(); ++i)
0145 edm::LogVerbatim("HGCalGeom") << "LayerCenter [" << i << "] " << layerCenter_[i];
0146 #endif
0147 if (firstLayer_ > 0) {
0148 for (unsigned int i = 0; i < layerType_.size(); ++i) {
0149 if (layerSense_[i] > 0) {
0150 int ii = layerType_[i];
0151 copyNumber_[ii] = firstLayer_;
0152 #ifdef EDM_ML_DEBUG
0153 edm::LogVerbatim("HGCalGeom") << "First copy number for layer type " << i << ":" << ii << " with "
0154 << materials_[ii] << " changed to " << copyNumber_[ii];
0155 #endif
0156 break;
0157 }
0158 }
0159 } else {
0160 firstLayer_ = 1;
0161 }
0162 #ifdef EDM_ML_DEBUG
0163 edm::LogVerbatim("HGCalGeom") << "There are " << layerType_.size() << " layers";
0164 for (unsigned int i = 0; i < layerType_.size(); ++i)
0165 edm::LogVerbatim("HGCalGeom") << "Layer [" << i << "] with material type " << layerType_[i] << " sensitive class "
0166 << layerSense_[i];
0167 #endif
0168 zMinBlock_ = nArgs["zMinBlock"];
0169 rad100to200_ = vArgs["rad100to200"];
0170 rad200to300_ = vArgs["rad200to300"];
0171 zMinRadPar_ = nArgs["zMinForRadPar"];
0172 choiceType_ = (int)(nArgs["choiceType"]);
0173 nCutRadPar_ = (int)(nArgs["nCornerCut"]);
0174 fracAreaMin_ = nArgs["fracAreaMin"];
0175 waferSize_ = nArgs["waferSize"];
0176 waferSepar_ = nArgs["SensorSeparation"];
0177 sectors_ = (int)(nArgs["Sectors"]);
0178 alpha_ = (1._pi) / sectors_;
0179 cosAlpha_ = cos(alpha_);
0180 #ifdef EDM_ML_DEBUG
0181 edm::LogVerbatim("HGCalGeom") << "zStart " << zMinBlock_ << " radius for wafer type separation uses "
0182 << rad100to200_.size() << " parameters; zmin " << zMinRadPar_ << " cutoff "
0183 << choiceType_ << ":" << nCutRadPar_ << ":" << fracAreaMin_ << " wafer width "
0184 << waferSize_ << " separations " << waferSepar_ << " sectors " << sectors_ << ":"
0185 << convertRadToDeg(alpha_) << ":" << cosAlpha_;
0186 for (unsigned int k = 0; k < rad100to200_.size(); ++k)
0187 edm::LogVerbatim("HGCalGeom") << "[" << k << "] 100-200 " << rad100to200_[k] << " 200-300 " << rad200to300_[k];
0188 #endif
0189 slopeB_ = vArgs["SlopeBottom"];
0190 zFrontB_ = vArgs["ZFrontBottom"];
0191 rMinFront_ = vArgs["RMinFront"];
0192 slopeT_ = vArgs["SlopeTop"];
0193 zFrontT_ = vArgs["ZFrontTop"];
0194 rMaxFront_ = vArgs["RMaxFront"];
0195 #ifdef EDM_ML_DEBUG
0196 for (unsigned int i = 0; i < slopeB_.size(); ++i)
0197 edm::LogVerbatim("HGCalGeom") << "Block [" << i << "] Zmin " << zFrontB_[i] << " Rmin " << rMinFront_[i]
0198 << " Slope " << slopeB_[i];
0199 for (unsigned int i = 0; i < slopeT_.size(); ++i)
0200 edm::LogVerbatim("HGCalGeom") << "Block [" << i << "] Zmin " << zFrontT_[i] << " Rmax " << rMaxFront_[i]
0201 << " Slope " << slopeT_[i];
0202 #endif
0203 nameSpace_ = DDCurrentNamespace::ns();
0204 #ifdef EDM_ML_DEBUG
0205 edm::LogVerbatim("HGCalGeom") << "DDHGCalEEAlgo: NameSpace " << nameSpace_;
0206 #endif
0207
0208 waferType_ = std::make_unique<HGCalWaferType>(
0209 rad100to200_, rad200to300_, (waferSize_ + waferSepar_), zMinRadPar_, choiceType_, nCutRadPar_, fracAreaMin_);
0210 }
0211
0212
0213
0214
0215
0216 void DDHGCalEEAlgo::execute(DDCompactView& cpv) {
0217 #ifdef EDM_ML_DEBUG
0218 edm::LogVerbatim("HGCalGeom") << "==>> Constructing DDHGCalEEAlgo...";
0219 copies_.clear();
0220 #endif
0221 constructLayers(parent(), cpv);
0222 #ifdef EDM_ML_DEBUG
0223 edm::LogVerbatim("HGCalGeom") << "DDHGCalEEAlgo: " << copies_.size() << " different wafer copy numbers";
0224 int k(0);
0225 for (std::unordered_set<int>::const_iterator itr = copies_.begin(); itr != copies_.end(); ++itr, ++k) {
0226 edm::LogVerbatim("HGCalGeom") << "Copy [" << k << "] : " << (*itr);
0227 }
0228 copies_.clear();
0229 edm::LogVerbatim("HGCalGeom") << "<<== End of DDHGCalEEAlgo construction...";
0230 #endif
0231 }
0232
0233 void DDHGCalEEAlgo::constructLayers(const DDLogicalPart& module, DDCompactView& cpv) {
0234 #ifdef EDM_ML_DEBUG
0235 edm::LogVerbatim("HGCalGeom") << "DDHGCalEEAlgo: \t\tInside Layers";
0236 #endif
0237 double zi(zMinBlock_);
0238 int laymin(0);
0239 for (unsigned int i = 0; i < layers_.size(); i++) {
0240 double zo = zi + layerThick_[i];
0241 double routF = HGCalGeomTools::radius(zi, zFrontT_, rMaxFront_, slopeT_);
0242 int laymax = laymin + layers_[i];
0243 double zz = zi;
0244 double thickTot(0);
0245 for (int ly = laymin; ly < laymax; ++ly) {
0246 int ii = layerType_[ly];
0247 int copy = copyNumber_[ii];
0248 double hthick = 0.5 * thick_[ii];
0249 double rinB = HGCalGeomTools::radius(zo - tol1_, zFrontB_, rMinFront_, slopeB_);
0250 zz += hthick;
0251 thickTot += thick_[ii];
0252
0253 std::string name = names_[ii] + std::to_string(copy);
0254 #ifdef EDM_ML_DEBUG
0255 edm::LogVerbatim("HGCalGeom") << "DDHGCalEEAlgo: Layer " << ly << ":" << ii << " Front " << zi << ", " << routF
0256 << " Back " << zo << ", " << rinB << " superlayer thickness " << layerThick_[i];
0257 #endif
0258 DDName matName(DDSplit(materials_[ii]).first, DDSplit(materials_[ii]).second);
0259 DDMaterial matter(matName);
0260 DDLogicalPart glog;
0261 if (layerSense_[ly] < 1) {
0262 std::vector<double> pgonZ, pgonRin, pgonRout;
0263 if (layerSense_[ly] == 0 || absorbMode_ == 0) {
0264 double rmax = routF * cosAlpha_ - tol1_;
0265 pgonZ.emplace_back(-hthick);
0266 pgonZ.emplace_back(hthick);
0267 pgonRin.emplace_back(rinB);
0268 pgonRin.emplace_back(rinB);
0269 pgonRout.emplace_back(rmax);
0270 pgonRout.emplace_back(rmax);
0271 } else {
0272 HGCalGeomTools::radius(zz - hthick,
0273 zz + hthick,
0274 zFrontB_,
0275 rMinFront_,
0276 slopeB_,
0277 zFrontT_,
0278 rMaxFront_,
0279 slopeT_,
0280 -layerSense_[ly],
0281 pgonZ,
0282 pgonRin,
0283 pgonRout);
0284 #ifdef EDM_ML_DEBUG
0285 edm::LogVerbatim("HGCalGeom") << "DDHGCalEEAlgo: z " << (zz - hthick) << ":" << (zz + hthick) << " with "
0286 << pgonZ.size() << " palnes";
0287 for (unsigned int isec = 0; isec < pgonZ.size(); ++isec)
0288 edm::LogVerbatim("HGCalGeom")
0289 << "[" << isec << "] z " << pgonZ[isec] << " R " << pgonRin[isec] << ":" << pgonRout[isec];
0290 #endif
0291 for (unsigned int isec = 0; isec < pgonZ.size(); ++isec) {
0292 pgonZ[isec] -= zz;
0293 pgonRout[isec] = pgonRout[isec] * cosAlpha_ - tol1_;
0294 }
0295 }
0296 DDSolid solid =
0297 DDSolidFactory::polyhedra(DDName(name, nameSpace_), sectors_, -alpha_, 2._pi, pgonZ, pgonRin, pgonRout);
0298 glog = DDLogicalPart(solid.ddname(), matter, solid);
0299 #ifdef EDM_ML_DEBUG
0300 edm::LogVerbatim("HGCalGeom") << "DDHGCalEEAlgo: " << solid.name() << " polyhedra of " << sectors_
0301 << " sectors covering " << convertRadToDeg(-alpha_) << ":"
0302 << convertRadToDeg(-alpha_ + 2._pi) << " with " << pgonZ.size()
0303 << " sections and filled with " << matName << ":" << &matter;
0304 for (unsigned int k = 0; k < pgonZ.size(); ++k)
0305 edm::LogVerbatim("HGCalGeom") << "[" << k << "] z " << pgonZ[k] << " R " << pgonRin[k] << ":" << pgonRout[k];
0306 #endif
0307 } else {
0308 double rins =
0309 (sensitiveMode_ < 1) ? rinB : HGCalGeomTools::radius(zz + hthick - tol1_, zFrontB_, rMinFront_, slopeB_);
0310 double routs =
0311 (sensitiveMode_ < 1) ? routF : HGCalGeomTools::radius(zz - hthick, zFrontT_, rMaxFront_, slopeT_);
0312 DDSolid solid = DDSolidFactory::tubs(DDName(name, nameSpace_), hthick, rins, routs, 0.0, 2._pi);
0313 glog = DDLogicalPart(solid.ddname(), matter, solid);
0314 #ifdef EDM_ML_DEBUG
0315 edm::LogVerbatim("HGCalGeom") << "DDHGCalEEFileAlgo: " << solid.name() << " Tubs made of " << matName << ":"
0316 << &matter << " of dimensions " << rinB << ":" << rins << ", " << routF << ":"
0317 << routs << ", " << hthick << ", 0.0, 360.0 and position " << glog.name()
0318 << " number " << copy << ":" << layerCenter_[copy - firstLayer_];
0319 #endif
0320 positionSensitive(glog, rins, routs, zz, layerSense_[ly], layerCenter_[copy - firstLayer_], cpv);
0321 }
0322 DDTranslation r1(0, 0, zz);
0323 DDRotation rot;
0324 cpv.position(glog, module, copy, r1, rot);
0325 ++copyNumber_[ii];
0326 #ifdef EDM_ML_DEBUG
0327 edm::LogVerbatim("HGCalGeom") << "DDHGCalEEAlgo: " << glog.name() << " number " << copy << " positioned in "
0328 << module.name() << " at " << r1 << " with no rotation";
0329 #endif
0330 zz += hthick;
0331 }
0332 zi = zo;
0333 laymin = laymax;
0334 if (std::abs(thickTot - layerThick_[i]) >= tol2_) {
0335 if (thickTot > layerThick_[i]) {
0336 edm::LogError("HGCalGeom") << "Thickness of the partition " << layerThick_[i] << " is smaller than " << thickTot
0337 << ": thickness of all its components **** ERROR ****";
0338 } else {
0339 edm::LogWarning("HGCalGeom") << "Thickness of the partition " << layerThick_[i] << " does not match with "
0340 << thickTot << " of the components";
0341 }
0342 }
0343 }
0344 }
0345
0346 void DDHGCalEEAlgo::positionSensitive(const DDLogicalPart& glog,
0347 double rin,
0348 double rout,
0349 double zpos,
0350 int layertype,
0351 int layercenter,
0352 DDCompactView& cpv) {
0353 static const double sqrt3 = std::sqrt(3.0);
0354 #ifdef EDM_ML_DEBUG
0355 std::vector<double> zl = {10000., 10494., 10536., 10578., 10619., 10661., 10702., 10791., 10879.};
0356 auto il = std::lower_bound(zl.begin(), zl.end(), zpos);
0357 int layer = static_cast<int>(il - zl.begin());
0358 #endif
0359 double r = 0.5 * (waferSize_ + waferSepar_);
0360 double R = 2.0 * r / sqrt3;
0361 double dy = 0.75 * R;
0362 int N = (int)(0.5 * rout / r) + 2;
0363 const auto& xyoff = geomTools_.shiftXY(layercenter, (waferSize_ + waferSepar_));
0364 #ifdef EDM_ML_DEBUG
0365 int ium(0), ivm(0), iumAll(0), ivmAll(0), kount(0), ntot(0), nin(0);
0366 std::vector<int> ntype(6, 0);
0367 edm::LogVerbatim("HGCalGeom") << "DDHGCalEEAlgo: " << glog.ddname() << " rin:rout " << rin << ":" << rout << " zpos "
0368 << zpos << " N " << N << " for maximum u, v; r " << r << " R " << R << " dy " << dy
0369 << " Shift " << xyoff.first << ":" << xyoff.second << " WaferSize "
0370 << (waferSize_ + waferSepar_);
0371 #endif
0372 for (int u = -N; u <= N; ++u) {
0373 for (int v = -N; v <= N; ++v) {
0374 int nr = 2 * v;
0375 int nc = -2 * u + v;
0376 double xpos = xyoff.first + nc * r;
0377 double ypos = xyoff.second + nr * dy;
0378 const auto& corner = HGCalGeomTools::waferCorner(xpos, ypos, r, R, rin, rout, false);
0379 #ifdef EDM_ML_DEBUG
0380 int iu = std::abs(u);
0381 int iv = std::abs(v);
0382 ++ntot;
0383 if (((corner.first <= 0) && std::abs(u) < 5 && std::abs(v) < 5) || (std::abs(u) < 2 && std::abs(v) < 2)) {
0384 edm::LogVerbatim("HGCalGeom") << "DDHGCalEEAlgo: " << glog.ddname() << " R " << rin << ":" << rout << "\n Z "
0385 << zpos << " LayerType " << layertype << " u " << u << " v " << v << " with "
0386 << corner.first << " corners";
0387 }
0388 #endif
0389 if (corner.first > 0) {
0390 int type = waferType_->getType(xpos, ypos, zpos);
0391 int copy = HGCalTypes::packTypeUV(type, u, v);
0392 if (layertype > 1)
0393 type += 3;
0394 #ifdef EDM_ML_DEBUG
0395 if (iu > ium)
0396 ium = iu;
0397 if (iv > ivm)
0398 ivm = iv;
0399 kount++;
0400 if (copies_.count(copy) == 0)
0401 copies_.insert(copy);
0402 #endif
0403 if (corner.first == (int)(HGCalParameters::k_CornerSize)) {
0404 #ifdef EDM_ML_DEBUG
0405 if (iu > iumAll)
0406 iumAll = iu;
0407 if (iv > ivmAll)
0408 ivmAll = iv;
0409 ++nin;
0410 #endif
0411 DDTranslation tran(xpos, ypos, 0.0);
0412 DDRotation rotation;
0413 DDName name = DDName(DDSplit(wafers_[type]).first, DDSplit(wafers_[type]).second);
0414 cpv.position(name, glog.ddname(), copy, tran, rotation);
0415 #ifdef EDM_ML_DEBUG
0416 static const std::vector<std::string> thickstr = {"h120", "l200", "l300", "h200"};
0417 char posit[20];
0418 sprintf(posit, "%8.3f %8.3f", xpos, ypos);
0419 double phi = convertRadToDeg(std::atan2(ypos, -xpos));
0420 if (phi < 0)
0421 phi += 360.0;
0422 int cass = (layer < 7) ? static_cast<int>(phi / 60.0) : static_cast<int>(phi / 30.0);
0423 edm::LogVerbatim("HGCalGeomX") << std::setw(2) << layer << " " << std::setw(2) << HGCalTypes::WaferFull << " "
0424 << thickstr[type] << " " << posit << " " << std::setw(2)
0425 << HGCalTypes::WaferOrient0 << " " << std::setw(2) << u << " " << std::setw(2)
0426 << v << std::setw(2) << cass;
0427 edm::LogVerbatim("HGCalGeom") << " DDEHGCalEEAlgo " << name << " Number " << copy << " u|v " << u << ":" << v
0428 << " Poaition " << xpos << ":" << ypos << ":" << zpos << " Angle "
0429 << convertRadToDeg(std::atan2(ypos, -xpos));
0430 ++ntype[type];
0431 edm::LogVerbatim("HGCalGeom") << " DDHGCalEEAlgo: " << name << " number " << copy << " positioned in "
0432 << glog.ddname() << " at " << tran << " with no rotation";
0433 #endif
0434 }
0435 }
0436 }
0437 }
0438 #ifdef EDM_ML_DEBUG
0439 edm::LogVerbatim("HGCalGeom") << "DDHGCalEEAlgo: Maximum # of u " << ium << ":" << iumAll << " # of v " << ivm << ":"
0440 << ivmAll << " and " << nin << ":" << kount << ":" << ntot << " wafers (" << ntype[0]
0441 << ":" << ntype[1] << ":" << ntype[2] << ":" << ntype[3] << ":" << ntype[4] << ":"
0442 << ntype[5] << ") for " << glog.ddname() << " R " << rin << ":" << rout;
0443 #endif
0444 }
0445
0446 DEFINE_EDM_PLUGIN(DDAlgorithmFactory, DDHGCalEEAlgo, "hgcal:DDHGCalEEAlgo");