Back to home page

Project CMSSW displayed by LXR

 
 

    


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

0001 /*
0002  * DDHGCalModule.cc
0003  *
0004  *  Created on: 26-August-2019
0005  *      Author: Sunanda Banerjee
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 //#define EDM_ML_DEBUG
0019 #ifdef EDM_ML_DEBUG
0020 #include <unordered_set>
0021 #endif
0022 using namespace angle_units::operators;
0023 
0024 static long algorithm(dd4hep::Detector& /* description */, 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");    // Wafers
0031   auto materials = args.value<std::vector<std::string> >("MaterialNames");   // Materials
0032   const auto& names = args.value<std::vector<std::string> >("VolumeNames");  // Names
0033   const auto& thick = args.value<std::vector<double> >("Thickness");         // Thickness of the material
0034   std::vector<int> copyNumber;                                               // Initial copy numbers
0035   copyNumber.resize(materials.size(), 1);
0036   for (unsigned int i = 0; i < materials.size(); ++i) {
0037     if (materials[i] == "materials:M_NEMAFR4plate")
0038       materials[i] = "materials:M_NEMA FR4 plate";
0039   }
0040 #ifdef EDM_ML_DEBUG
0041   edm::LogVerbatim("HGCalGeom") << "DDHGCalModule: " << 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") << "DDHGCalModule: " << 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");             // Number of layers in a section
0051   const auto& layerThick = args.value<std::vector<double> >("LayerThick");  // Thickness of each section
0052   const auto& layerType = args.value<std::vector<int> >("LayerType");       // Type of the layer
0053   const auto& layerSense = args.value<std::vector<int> >("LayerSense");     // Content of a layer (sensitive?)
0054 #ifdef EDM_ML_DEBUG
0055   edm::LogVerbatim("HGCalGeom") << "DDHGCalModule: " << 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") << "DDHGCalModule: " << 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");  // Starting z-value of the block
0065   double rMaxFine = args.value<double>("rMaxFine");    // Maximum r-value for fine wafer
0066   double waferW = args.value<double>("waferW");        // Width of the wafer
0067   int sectors = args.value<int>("Sectors");            // Sectors
0068 #ifdef EDM_ML_DEBUG
0069   edm::LogVerbatim("HGCalGeom") << "DDHGCalModule: zStart " << cms::convert2mm(zMinBlock) << " rFineCoarse "
0070                                 << cms::convert2mm(rMaxFine) << " wafer width " << cms::convert2mm(waferW)
0071                                 << " sectors " << sectors;
0072 #endif
0073   const auto& slopeB = args.value<std::vector<double> >("SlopeBottom");   // Slope at the lower R
0074   const auto& slopeT = args.value<std::vector<double> >("SlopeTop");      // Slopes at the larger R
0075   const auto& zFront = args.value<std::vector<double> >("ZFront");        // Starting Z values for the slopes
0076   const auto& rMaxFront = args.value<std::vector<double> >("RMaxFront");  // Corresponding rMax's
0077 #ifdef EDM_ML_DEBUG
0078   edm::LogVerbatim("HGCalGeom") << "DDHGCalModule: Bottom slopes " << slopeB[0] << ":" << slopeB[1] << " and "
0079                                 << slopeT.size() << " slopes for top";
0080   for (unsigned int i = 0; i < slopeT.size(); ++i)
0081     edm::LogVerbatim("HGCalGeom") << "Block [" << i << "] Zmin " << cms::convert2mm(zFront[i]) << " Rmax "
0082                                   << cms::convert2mm(rMaxFront[i]) << " Slope " << slopeT[i];
0083 #endif
0084   std::string idNameSpace = static_cast<std::string>(ns.name());  // Namespace of this and ALL sub-parts
0085   const auto& idName = args.parentName();                         // Name of the "parent" volume.
0086 #ifdef EDM_ML_DEBUG
0087   std::unordered_set<int> copies;  // List of copy #'s
0088   edm::LogVerbatim("HGCalGeom") << "DDHGCalModule: NameSpace " << idNameSpace << " Mother " << idName;
0089 #endif
0090 
0091   // Mother module
0092   dd4hep::Volume module = ns.volume(idName);
0093 
0094   double zi(zMinBlock);
0095   int laymin(0);
0096   for (unsigned int i = 0; i < layers.size(); i++) {
0097     double zo = zi + layerThick[i];
0098     double routF = HGCalGeomTools::radius(zi, zFront, rMaxFront, slopeT);
0099     int laymax = laymin + layers[i];
0100     double zz = zi;
0101     double thickTot(0);
0102     for (int ly = laymin; ly < laymax; ++ly) {
0103       int ii = layerType[ly];
0104       int copy = copyNumber[ii];
0105       double rinB = (layerSense[ly] == 0) ? (zo * slopeB[0]) : (zo * slopeB[1]);
0106       zz += (0.5 * thick[ii]);
0107       thickTot += thick[ii];
0108 
0109       std::string name = "HGCal" + names[ii] + std::to_string(copy);
0110 #ifdef EDM_ML_DEBUG
0111       edm::LogVerbatim("HGCalGeom") << "DDHGCalModule: Layer " << ly << ":" << ii << " Front " << cms::convert2mm(zi)
0112                                     << ", " << cms::convert2mm(routF) << " Back " << cms::convert2mm(zo) << ", "
0113                                     << cms::convert2mm(rinB) << " superlayer thickness "
0114                                     << cms::convert2mm(layerThick[i]);
0115 #endif
0116       dd4hep::Material matter = ns.material(materials[ii]);
0117       dd4hep::Volume glog;
0118       if (layerSense[ly] == 0) {
0119         double alpha = 1._pi / sectors;
0120         double rmax = routF * cos(alpha) - tol;
0121         std::vector<double> pgonZ, pgonRin, pgonRout;
0122         pgonZ.emplace_back(-0.5 * thick[ii]);
0123         pgonZ.emplace_back(0.5 * thick[ii]);
0124         pgonRin.emplace_back(rinB);
0125         pgonRin.emplace_back(rinB);
0126         pgonRout.emplace_back(rmax);
0127         pgonRout.emplace_back(rmax);
0128         dd4hep::Solid solid = dd4hep::Polyhedra(sectors, -alpha, 2._pi, pgonZ, pgonRin, pgonRout);
0129         ns.addSolidNS(ns.prepend(name), solid);
0130         glog = dd4hep::Volume(solid.name(), solid, matter);
0131 #ifdef EDM_ML_DEBUG
0132         edm::LogVerbatim("HGCalGeom") << "DDHGCalModule: " << solid.name() << " polyhedra of " << sectors
0133                                       << " sectors covering " << convertRadToDeg(-alpha) << ":"
0134                                       << (360.0 + convertRadToDeg(-alpha)) << " with " << pgonZ.size() << " sections";
0135         for (unsigned int k = 0; k < pgonZ.size(); ++k)
0136           edm::LogVerbatim("HGCalGeom") << "[" << k << "] z " << cms::convert2mm(pgonZ[k]) << " R "
0137                                         << cms::convert2mm(pgonRin[k]) << ":" << cms::convert2mm(pgonRout[k]);
0138 #endif
0139       } else {
0140         dd4hep::Solid solid = dd4hep::Tube(0.5 * thick[ii], rinB, routF, 0.0, 2._pi);
0141         ns.addSolidNS(ns.prepend(name), solid);
0142         glog = dd4hep::Volume(solid.name(), solid, matter);
0143 #ifdef EDM_ML_DEBUG
0144         edm::LogVerbatim("HGCalGeom") << "DDHGCalModule: " << solid.name() << " Tubs made of " << materials[ii]
0145                                       << " of dimensions " << cms::convert2mm(rinB) << ", " << cms::convert2mm(routF)
0146                                       << ", " << cms::convert2mm(0.5 * thick[ii]) << ", 0.0, 360.0";
0147         edm::LogVerbatim("HGCalGeom") << "DDHGCalModule test position in: " << glog.name() << " number " << copy;
0148 #endif
0149         double dx = 0.5 * waferW;
0150         double dy = 3.0 * dx * tan(30._deg);
0151         double rr = 2.0 * dx * tan(30._deg);
0152         int ncol = static_cast<int>(2.0 * routF / waferW) + 1;
0153         int nrow = static_cast<int>(routF / (waferW * tan(30._deg))) + 1;
0154 #ifdef EDM_ML_DEBUG
0155         int incm(0), inrm(0), kount(0), ntot(0), nin(0), nfine(0), ncoarse(0);
0156         edm::LogVerbatim("HGCalGeom") << glog.name() << " rout " << cms::convert2mm(routF) << " Row " << nrow
0157                                       << " Column " << ncol;
0158 #endif
0159         for (int nr = -nrow; nr <= nrow; ++nr) {
0160           int inr = (nr >= 0) ? nr : -nr;
0161           for (int nc = -ncol; nc <= ncol; ++nc) {
0162             int inc = (nc >= 0) ? nc : -nc;
0163             if (inr % 2 == inc % 2) {
0164               double xpos = nc * dx;
0165               double ypos = nr * dy;
0166               std::pair<int, int> corner = HGCalGeomTools::waferCorner(xpos, ypos, dx, rr, rinB, routF, true);
0167 #ifdef EDM_ML_DEBUG
0168               ++ntot;
0169 #endif
0170               if (corner.first > 0) {
0171                 int copyL = HGCalTypes::packTypeUV(0, nc, nr);
0172 #ifdef EDM_ML_DEBUG
0173                 if (inc > incm)
0174                   incm = inc;
0175                 if (inr > inrm)
0176                   inrm = inr;
0177                 kount++;
0178                 copies.insert(copy);
0179 #endif
0180                 if (corner.first == (int)(HGCalParameters::k_CornerSize)) {
0181                   double rpos = std::sqrt(xpos * xpos + ypos * ypos);
0182                   dd4hep::Position tran(xpos, ypos, 0.0);
0183                   dd4hep::Rotation3D rotation;
0184                   dd4hep::Volume glog1 = (rpos < rMaxFine) ? ns.volume(wafer[0]) : ns.volume(wafer[1]);
0185                   glog.placeVolume(glog1, copyL, dd4hep::Transform3D(rotation, tran));
0186 #ifdef EDM_ML_DEBUG
0187                   ++nin;
0188                   if (rpos < rMaxFine)
0189                     ++nfine;
0190                   else
0191                     ++ncoarse;
0192                   edm::LogVerbatim("HGCalGeom")
0193                       << "DDHGCalModule: " << glog1.name() << " number " << copyL << " positioned in " << glog.name()
0194                       << " at (" << cms::convert2mm(xpos) << "," << cms::convert2mm(ypos) << ",0) with " << rotation;
0195 #endif
0196                 }
0197               }
0198             }
0199           }
0200         }
0201 #ifdef EDM_ML_DEBUG
0202         edm::LogVerbatim("HGCalGeom") << "DDHGCalModule: # of columns " << incm << " # of rows " << inrm << " and "
0203                                       << nin << ":" << kount << ":" << ntot << " wafers (" << nfine << ":" << ncoarse
0204                                       << ") for " << glog.name() << " R " << cms::convert2mm(rinB) << ":"
0205                                       << cms::convert2mm(routF);
0206 #endif
0207       }
0208       dd4hep::Position r1(0, 0, zz);
0209       dd4hep::Rotation3D rot;
0210       module.placeVolume(glog, copy, dd4hep::Transform3D(rot, r1));
0211       ++copyNumber[ii];
0212 #ifdef EDM_ML_DEBUG
0213       edm::LogVerbatim("HGCalGeom") << "DDHGCalModule: " << glog.name() << " number " << copy << " positioned in "
0214                                     << module.name() << " at (0,0," << cms::convert2mm(zz) << ") with no rotation";
0215 #endif
0216       zz += (0.5 * thick[ii]);
0217     }  // End of loop over layers in a block
0218     zi = zo;
0219     laymin = laymax;
0220     if (fabs(thickTot - layerThick[i]) > tol2) {
0221       if (thickTot > layerThick[i]) {
0222         edm::LogError("HGCalGeom") << "Thickness of the partition " << cms::convert2mm(layerThick[i])
0223                                    << " is smaller than thickness " << cms::convert2mm(thickTot)
0224                                    << " of all its components **** ERROR ****\n";
0225       } else {
0226         edm::LogWarning("HGCalGeom") << "Thickness of the partition " << cms::convert2mm(layerThick[i])
0227                                      << " does not match with " << cms::convert2mm(thickTot) << " of the components\n";
0228       }
0229     }
0230   }  // End of loop over blocks
0231 
0232 #ifdef EDM_ML_DEBUG
0233   edm::LogVerbatim("HGCalGeom") << copies.size() << " different wafer copy numbers";
0234   int k(0);
0235   for (std::unordered_set<int>::const_iterator itr = copies.begin(); itr != copies.end(); ++itr, ++k)
0236     edm::LogVerbatim("HGCalGeom") << "Copy[" << k << "] : " << (*itr);
0237   edm::LogVerbatim("HGCalGeom") << "<<== End of DDHGCalModule construction ...";
0238 #endif
0239 
0240   return cms::s_executed;
0241 }
0242 
0243 // first argument is the type from the xml file
0244 DECLARE_DDCMS_DETELEMENT(DDCMS_hgcal_DDHGCalModule, algorithm)