Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:15:00

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