Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:14:59

0001 ///////////////////////////////////////////////////////////////////////////////
0002 // File: DDHGCalModule.cc
0003 // Description: Geometry factory class for HGCal (EE and HESil)
0004 ///////////////////////////////////////////////////////////////////////////////
0005 
0006 #include <algorithm>
0007 #include <cmath>
0008 #include <map>
0009 #include <string>
0010 #include <unordered_set>
0011 #include <vector>
0012 
0013 #include "DataFormats/Math/interface/angle_units.h"
0014 #include "DetectorDescription/Core/interface/DDAlgorithm.h"
0015 #include "DetectorDescription/Core/interface/DDAlgorithmFactory.h"
0016 #include "DetectorDescription/Core/interface/DDCurrentNamespace.h"
0017 #include "DetectorDescription/Core/interface/DDLogicalPart.h"
0018 #include "DetectorDescription/Core/interface/DDMaterial.h"
0019 #include "DetectorDescription/Core/interface/DDSolid.h"
0020 #include "DetectorDescription/Core/interface/DDSplit.h"
0021 #include "DetectorDescription/Core/interface/DDTypes.h"
0022 #include "DetectorDescription/Core/interface/DDutils.h"
0023 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0024 #include "FWCore/PluginManager/interface/PluginFactory.h"
0025 #include "Geometry/HGCalCommonData/interface/HGCalGeomTools.h"
0026 #include "Geometry/HGCalCommonData/interface/HGCalParameters.h"
0027 #include "Geometry/HGCalCommonData/interface/HGCalTypes.h"
0028 
0029 //#define EDM_ML_DEBUG
0030 using namespace angle_units::operators;
0031 
0032 class DDHGCalModule : public DDAlgorithm {
0033 public:
0034   // Constructor and Destructor
0035   DDHGCalModule();  // const std::string & name);
0036 
0037   void initialize(const DDNumericArguments& nArgs,
0038                   const DDVectorArguments& vArgs,
0039                   const DDMapArguments& mArgs,
0040                   const DDStringArguments& sArgs,
0041                   const DDStringVectorArguments& vsArgs) override;
0042   void execute(DDCompactView& cpv) override;
0043 
0044 protected:
0045   void constructLayers(const DDLogicalPart&, DDCompactView& cpv);
0046   double rMax(double z);
0047   void positionSensitive(DDLogicalPart& glog, double rin, double rout, DDCompactView& cpv);
0048 
0049 private:
0050   static constexpr double tol_ = 0.00001;
0051 
0052   std::vector<std::string> wafer_;      // Wafers
0053   std::vector<std::string> materials_;  // Materials
0054   std::vector<std::string> names_;      // Names
0055   std::vector<double> thick_;           // Thickness of the material
0056   std::vector<int> copyNumber_;         // Initial copy numbers
0057   std::vector<int> layers_;             // Number of layers in a section
0058   std::vector<double> layerThick_;      // Thickness of each section
0059   std::vector<int> layerType_;          // Type of the layer
0060   std::vector<int> layerSense_;         // COntent of a layer (sensitive?)
0061   double zMinBlock_;                    // Starting z-value of the block
0062   double rMaxFine_;                     // Maximum r-value for fine wafer
0063   double waferW_;                       // Width of the wafer
0064   int sectors_;                         // Sectors
0065   std::vector<double> slopeB_;          // Slope at the lower R
0066   std::vector<double> slopeT_;          // Slopes at the larger R
0067   std::vector<double> zFront_;          // Starting Z values for the slopes
0068   std::vector<double> rMaxFront_;       // Corresponding rMax's
0069   std::string idNameSpace_;             // Namespace of this and ALL sub-parts
0070   std::unordered_set<int> copies_;      // List of copy #'s
0071 };
0072 
0073 DDHGCalModule::DDHGCalModule() {
0074 #ifdef EDM_ML_DEBUG
0075   edm::LogVerbatim("HGCalGeom") << "DDHGCalModule: Creating an instance";
0076 #endif
0077 }
0078 
0079 void DDHGCalModule::initialize(const DDNumericArguments& nArgs,
0080                                const DDVectorArguments& vArgs,
0081                                const DDMapArguments&,
0082                                const DDStringArguments& sArgs,
0083                                const DDStringVectorArguments& vsArgs) {
0084   wafer_ = vsArgs["WaferName"];
0085 #ifdef EDM_ML_DEBUG
0086   edm::LogVerbatim("HGCalGeom") << "DDHGCalModule: " << wafer_.size() << " wafers";
0087   for (unsigned int i = 0; i < wafer_.size(); ++i)
0088     edm::LogVerbatim("HGCalGeom") << "Wafer[" << i << "] " << wafer_[i];
0089 #endif
0090   materials_ = vsArgs["MaterialNames"];
0091   names_ = vsArgs["VolumeNames"];
0092   thick_ = vArgs["Thickness"];
0093   copyNumber_.resize(materials_.size(), 1);
0094 #ifdef EDM_ML_DEBUG
0095   edm::LogVerbatim("HGCalGeom") << "DDHGCalModule: " << materials_.size() << " types of volumes";
0096   for (unsigned int i = 0; i < names_.size(); ++i)
0097     edm::LogVerbatim("HGCalGeom") << "Volume [" << i << "] " << names_[i] << " of thickness " << thick_[i]
0098                                   << " filled with " << materials_[i] << " first copy number " << copyNumber_[i];
0099 #endif
0100   layers_ = dbl_to_int(vArgs["Layers"]);
0101   layerThick_ = vArgs["LayerThick"];
0102 #ifdef EDM_ML_DEBUG
0103   edm::LogVerbatim("HGCalGeom") << "DDHGCalModule: " << layers_.size() << " blocks";
0104   for (unsigned int i = 0; i < layers_.size(); ++i)
0105     edm::LogVerbatim("HGCalGeom") << "Block [" << i << "] of thickness " << layerThick_[i] << " with " << layers_[i]
0106                                   << " layers";
0107 #endif
0108   layerType_ = dbl_to_int(vArgs["LayerType"]);
0109   layerSense_ = dbl_to_int(vArgs["LayerSense"]);
0110 #ifdef EDM_ML_DEBUG
0111   edm::LogVerbatim("HGCalGeom") << "DDHGCalModule: " << layerType_.size() << " layers";
0112   for (unsigned int i = 0; i < layerType_.size(); ++i)
0113     edm::LogVerbatim("HGCalGeom") << "Layer [" << i << "] with material type " << layerType_[i] << " sensitive class "
0114                                   << layerSense_[i];
0115 #endif
0116   zMinBlock_ = nArgs["zMinBlock"];
0117   rMaxFine_ = nArgs["rMaxFine"];
0118   waferW_ = nArgs["waferW"];
0119   sectors_ = (int)(nArgs["Sectors"]);
0120 #ifdef EDM_ML_DEBUG
0121   edm::LogVerbatim("HGCalGeom") << "DDHGCalModule: zStart " << zMinBlock_ << " rFineCoarse " << rMaxFine_
0122                                 << " wafer width " << waferW_ << " sectors " << sectors_;
0123 #endif
0124   slopeB_ = vArgs["SlopeBottom"];
0125   slopeT_ = vArgs["SlopeTop"];
0126   zFront_ = vArgs["ZFront"];
0127   rMaxFront_ = vArgs["RMaxFront"];
0128 #ifdef EDM_ML_DEBUG
0129   edm::LogVerbatim("HGCalGeom") << "DDHGCalModule: Bottom slopes " << slopeB_[0] << ":" << slopeB_[1] << " and "
0130                                 << slopeT_.size() << " slopes for top";
0131   for (unsigned int i = 0; i < slopeT_.size(); ++i)
0132     edm::LogVerbatim("HGCalGeom") << "Block [" << i << "] Zmin " << zFront_[i] << " Rmax " << rMaxFront_[i] << " Slope "
0133                                   << slopeT_[i];
0134 #endif
0135   idNameSpace_ = DDCurrentNamespace::ns();
0136 #ifdef EDM_ML_DEBUG
0137   edm::LogVerbatim("HGCalGeom") << "DDHGCalModule: NameSpace " << idNameSpace_;
0138 #endif
0139 }
0140 
0141 ////////////////////////////////////////////////////////////////////
0142 // DDHGCalModule methods...
0143 ////////////////////////////////////////////////////////////////////
0144 
0145 void DDHGCalModule::execute(DDCompactView& cpv) {
0146 #ifdef EDM_ML_DEBUG
0147   edm::LogVerbatim("HGCalGeom") << "==>> Constructing DDHGCalModule...";
0148 #endif
0149   copies_.clear();
0150   constructLayers(parent(), cpv);
0151 #ifdef EDM_ML_DEBUG
0152   edm::LogVerbatim("HGCalGeom") << copies_.size() << " different wafer copy numbers";
0153   int k(0);
0154   for (std::unordered_set<int>::const_iterator itr = copies_.begin(); itr != copies_.end(); ++itr, ++k)
0155     edm::LogVerbatim("HGCalGeom") << "Copy[" << k << "] : " << (*itr);
0156 #endif
0157   copies_.clear();
0158 #ifdef EDM_ML_DEBUG
0159   edm::LogVerbatim("HGCalGeom") << "<<== End of DDHGCalModule construction ...";
0160 #endif
0161 }
0162 
0163 void DDHGCalModule::constructLayers(const DDLogicalPart& module, DDCompactView& cpv) {
0164 #ifdef EDM_ML_DEBUG
0165   edm::LogVerbatim("HGCalGeom") << "DDHGCalModule: \t\tInside Layers";
0166 #endif
0167   double zi(zMinBlock_);
0168   int laymin(0);
0169   const double tol(0.01);
0170   for (unsigned int i = 0; i < layers_.size(); i++) {
0171     double zo = zi + layerThick_[i];
0172     double routF = rMax(zi);
0173     int laymax = laymin + layers_[i];
0174     double zz = zi;
0175     double thickTot(0);
0176     for (int ly = laymin; ly < laymax; ++ly) {
0177       int ii = layerType_[ly];
0178       int copy = copyNumber_[ii];
0179       double rinB = (layerSense_[ly] == 0) ? (zo * slopeB_[0]) : (zo * slopeB_[1]);
0180       zz += (0.5 * thick_[ii]);
0181       thickTot += thick_[ii];
0182 
0183       std::string name = "HGCal" + names_[ii] + std::to_string(copy);
0184 #ifdef EDM_ML_DEBUG
0185       edm::LogVerbatim("HGCalGeom") << "DDHGCalModule: Layer " << ly << ":" << ii << " Front " << zi << ", " << routF
0186                                     << " Back " << zo << ", " << rinB << " superlayer thickness " << layerThick_[i];
0187 #endif
0188       DDName matName(DDSplit(materials_[ii]).first, DDSplit(materials_[ii]).second);
0189       DDMaterial matter(matName);
0190       DDLogicalPart glog;
0191       if (layerSense_[ly] == 0) {
0192         double alpha = 1._pi / sectors_;
0193         double rmax = routF * cos(alpha) - tol;
0194         std::vector<double> pgonZ, pgonRin, pgonRout;
0195         pgonZ.emplace_back(-0.5 * thick_[ii]);
0196         pgonZ.emplace_back(0.5 * thick_[ii]);
0197         pgonRin.emplace_back(rinB);
0198         pgonRin.emplace_back(rinB);
0199         pgonRout.emplace_back(rmax);
0200         pgonRout.emplace_back(rmax);
0201         DDSolid solid =
0202             DDSolidFactory::polyhedra(DDName(name, idNameSpace_), sectors_, -alpha, 2._pi, pgonZ, pgonRin, pgonRout);
0203         glog = DDLogicalPart(solid.ddname(), matter, solid);
0204 #ifdef EDM_ML_DEBUG
0205         edm::LogVerbatim("HGCalGeom") << "DDHGCalModule: " << solid.name() << " polyhedra of " << sectors_
0206                                       << " sectors covering " << convertRadToDeg(-alpha) << ":"
0207                                       << (360.0 + convertRadToDeg(-alpha)) << " with " << pgonZ.size() << " sections";
0208         for (unsigned int k = 0; k < pgonZ.size(); ++k)
0209           edm::LogVerbatim("HGCalGeom") << "[" << k << "] z " << pgonZ[k] << " R " << pgonRin[k] << ":" << pgonRout[k];
0210 #endif
0211       } else {
0212         DDSolid solid = DDSolidFactory::tubs(DDName(name, idNameSpace_), 0.5 * thick_[ii], rinB, routF, 0.0, 2._pi);
0213         glog = DDLogicalPart(solid.ddname(), matter, solid);
0214 #ifdef EDM_ML_DEBUG
0215         edm::LogVerbatim("HGCalGeom") << "DDHGCalModule: " << solid.name() << " Tubs made of " << matName
0216                                       << " of dimensions " << rinB << ", " << routF << ", " << 0.5 * thick_[ii]
0217                                       << ", 0.0, 360.0";
0218         edm::LogVerbatim("HGCalGeom") << "DDHGCalModule test position in: " << glog.name() << " number " << copy;
0219 #endif
0220         positionSensitive(glog, rinB, routF, cpv);
0221       }
0222       DDTranslation r1(0, 0, zz);
0223       DDRotation rot;
0224       cpv.position(glog, module, copy, r1, rot);
0225       ++copyNumber_[ii];
0226 #ifdef EDM_ML_DEBUG
0227       edm::LogVerbatim("HGCalGeom") << "DDHGCalModule: " << glog.name() << " number " << copy << " positioned in "
0228                                     << module.name() << " at " << r1 << " with " << rot;
0229 #endif
0230       zz += (0.5 * thick_[ii]);
0231     }  // End of loop over layers in a block
0232     zi = zo;
0233     laymin = laymax;
0234     if (fabs(thickTot - layerThick_[i]) > tol_) {
0235       if (thickTot > layerThick_[i]) {
0236         edm::LogError("HGCalGeom") << "Thickness of the partition " << layerThick_[i] << " is smaller than thickness "
0237                                    << thickTot << " of all its components **** ERROR ****\n";
0238       } else {
0239         edm::LogWarning("HGCalGeom") << "Thickness of the partition " << layerThick_[i] << " does not match with "
0240                                      << thickTot << " of the components\n";
0241       }
0242     }
0243   }  // End of loop over blocks
0244 }
0245 
0246 double DDHGCalModule::rMax(double z) {
0247   double r(0);
0248 #ifdef EDM_ML_DEBUG
0249   unsigned int ik(0);
0250 #endif
0251   for (unsigned int k = 0; k < slopeT_.size(); ++k) {
0252     if (z < zFront_[k])
0253       break;
0254     r = rMaxFront_[k] + (z - zFront_[k]) * slopeT_[k];
0255 #ifdef EDM_ML_DEBUG
0256     ik = k;
0257 #endif
0258   }
0259 #ifdef EDM_ML_DEBUG
0260   edm::LogVerbatim("HGCalGeom") << "rMax : " << z << ":" << ik << ":" << r;
0261 #endif
0262   return r;
0263 }
0264 
0265 void DDHGCalModule::positionSensitive(DDLogicalPart& glog, double rin, double rout, DDCompactView& cpv) {
0266   double dx = 0.5 * waferW_;
0267   double dy = 3.0 * dx * tan(30._deg);
0268   double rr = 2.0 * dx * tan(30._deg);
0269   int ncol = (int)(2.0 * rout / waferW_) + 1;
0270   int nrow = (int)(rout / (waferW_ * tan(30._deg))) + 1;
0271   int incm(0), inrm(0);
0272 #ifdef EDM_ML_DEBUG
0273   int kount(0), ntot(0), nin(0), nfine(0), ncoarse(0);
0274   edm::LogVerbatim("HGCalGeom") << glog.ddname() << " rout " << rout << " Row " << nrow << " Column " << ncol;
0275 #endif
0276   for (int nr = -nrow; nr <= nrow; ++nr) {
0277     int inr = (nr >= 0) ? nr : -nr;
0278     for (int nc = -ncol; nc <= ncol; ++nc) {
0279       int inc = (nc >= 0) ? nc : -nc;
0280       if (inr % 2 == inc % 2) {
0281         double xpos = nc * dx;
0282         double ypos = nr * dy;
0283         auto const& corner = HGCalGeomTools::waferCorner(xpos, ypos, dx, rr, rin, rout, true);
0284 #ifdef EDM_ML_DEBUG
0285         ++ntot;
0286 #endif
0287         if (corner.first > 0) {
0288           int copy = HGCalTypes::packTypeUV(0, nc, nr);
0289           if (inc > incm)
0290             incm = inc;
0291           if (inr > inrm)
0292             inrm = inr;
0293 #ifdef EDM_ML_DEBUG
0294           kount++;
0295 #endif
0296           if (copies_.count(copy) == 0)
0297             copies_.insert(copy);
0298           if (corner.first == (int)(HGCalParameters::k_CornerSize)) {
0299             double rpos = std::sqrt(xpos * xpos + ypos * ypos);
0300             DDTranslation tran(xpos, ypos, 0.0);
0301             DDRotation rotation;
0302 #ifdef EDM_ML_DEBUG
0303             ++nin;
0304 #endif
0305             DDName name = (rpos < rMaxFine_) ? DDName(DDSplit(wafer_[0]).first, DDSplit(wafer_[0]).second)
0306                                              : DDName(DDSplit(wafer_[1]).first, DDSplit(wafer_[1]).second);
0307             cpv.position(name, glog.ddname(), copy, tran, rotation);
0308 #ifdef EDM_ML_DEBUG
0309             if (rpos < rMaxFine_)
0310               ++nfine;
0311             else
0312               ++ncoarse;
0313             edm::LogVerbatim("HGCalGeom") << "DDHGCalModule: " << name << " number " << copy << " positioned in "
0314                                           << glog.ddname() << " at " << tran << " with " << rotation;
0315 #endif
0316           }
0317         }
0318       }
0319     }
0320   }
0321 #ifdef EDM_ML_DEBUG
0322   edm::LogVerbatim("HGCalGeom") << "DDHGCalModule: # of columns " << incm << " # of rows " << inrm << " and " << nin
0323                                 << ":" << kount << ":" << ntot << " wafers (" << nfine << ":" << ncoarse << ") for "
0324                                 << glog.ddname() << " R " << rin << ":" << rout;
0325 #endif
0326 }
0327 
0328 DEFINE_EDM_PLUGIN(DDAlgorithmFactory, DDHGCalModule, "hgcal:DDHGCalModule");