Back to home page

Project CMSSW displayed by LXR

 
 

    


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

0001 ///////////////////////////////////////////////////////////////////////////////
0002 // File: DDHCalEndcapModuleAlgo.cc
0003 //   adapted from CCal(G4)HcalEndcap.cc
0004 // Description: Geometry factory class for Hcal Endcap
0005 ///////////////////////////////////////////////////////////////////////////////
0006 
0007 #include <cmath>
0008 #include <algorithm>
0009 #include <string>
0010 #include <vector>
0011 
0012 #include "DataFormats/Math/interface/angle_units.h"
0013 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0014 #include "DD4hep/DetFactoryHelper.h"
0015 #include "DetectorDescription/DDCMS/interface/DDPlugins.h"
0016 #include "DetectorDescription/DDCMS/interface/DDutils.h"
0017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0018 
0019 //#define EDM_ML_DEBUG
0020 using namespace angle_units::operators;
0021 
0022 struct HCalEndcapModuleAlgo {
0023   std::string genMaterial;             //General material
0024   std::string absorberMat;             //Absorber material
0025   std::string plasticMat;              //Plastic material cover
0026   std::string scintMat;                //Scintillator material
0027   std::string rotstr;                  //Rotation matrix to place in mother
0028   int sectors;                         //Number of potenital straight edges
0029   double zMinBlock;                    //Minimum z extent of the block
0030   double zMaxBlock;                    //Maximum
0031   double z1Beam;                       //Position of gap end   along z-axis
0032   double ziDip;                        //Starting Z of dipped part of body
0033   double dzStep;                       //Width in Z of a layer
0034   double moduleThick;                  //Thickness of a layer (air/absorber)
0035   double layerThick;                   //Thickness of a layer (plastic)
0036   double scintThick;                   //Thickness of scinitllator
0037   double rMaxBack;                     //Maximum R after  the dip
0038   double rMaxFront;                    //Maximum R before the dip
0039   double slopeBot;                     //Slope of the bottom edge
0040   double slopeTop;                     //Slope of the top edge
0041   double slopeTopF;                    //Slope of the top front edge
0042   double trimLeft;                     //Trim of the left edge
0043   double trimRight;                    //Trim of the right edge
0044   double tolAbs;                       //Tolerance for absorber
0045   int modType;                         //Type of module
0046   int modNumber;                       //Module number
0047   int layerType;                       //Layer type
0048   std::vector<int> layerNumber;        //layer numbers
0049   std::vector<std::string> phiName;    //Name of Phi sections
0050   std::vector<std::string> layerName;  //Layer Names
0051 
0052   std::string idName;   //Name of the "parent" volume.
0053   std::string modName;  //Module Name
0054   int idOffset;         // Geant4 ID's...    = 4000;
0055 
0056   struct HcalEndcapPar {
0057     double yh1, bl1, tl1, yh2, bl2, tl2, alp, theta, phi, xpos, ypos, zpos;
0058     HcalEndcapPar(double yh1v = 0,
0059                   double bl1v = 0,
0060                   double tl1v = 0,
0061                   double yh2v = 0,
0062                   double bl2v = 0,
0063                   double tl2v = 0,
0064                   double alpv = 0,
0065                   double thv = 0,
0066                   double fiv = 0,
0067                   double x = 0,
0068                   double y = 0,
0069                   double z = 0)
0070         : yh1(yh1v),
0071           bl1(bl1v),
0072           tl1(tl1v),
0073           yh2(yh2v),
0074           bl2(bl2v),
0075           tl2(tl2v),
0076           alp(alpv),
0077           theta(thv),
0078           phi(fiv),
0079           xpos(x),
0080           ypos(y),
0081           zpos(z) {}
0082   };
0083 
0084   HCalEndcapModuleAlgo() = delete;
0085 
0086   HCalEndcapModuleAlgo(cms::DDParsingContext& ctxt, xml_h e) {
0087     cms::DDNamespace ns(ctxt, e, true);
0088     cms::DDAlgoArguments args(ctxt, e);
0089 
0090     genMaterial = args.value<std::string>("MaterialName");
0091     absorberMat = args.value<std::string>("AbsorberMat");
0092     plasticMat = args.value<std::string>("PlasticMat");
0093     scintMat = args.value<std::string>("ScintMat");
0094     rotstr = args.value<std::string>("Rotation");
0095     sectors = args.value<int>("Sectors");
0096 #ifdef EDM_ML_DEBUG
0097     edm::LogVerbatim("HCalGeom") << "DDHCalEndcapModuleAlgo: General material " << genMaterial << "\tAbsorber "
0098                                  << absorberMat << "\tPlastic " << plasticMat << "\tScintillator " << scintMat
0099                                  << "\tRotation " << rotstr << "\tSectors " << sectors;
0100 #endif
0101     zMinBlock = args.value<double>("ZMinBlock");
0102     zMaxBlock = args.value<double>("ZMaxBlock");
0103     z1Beam = args.value<double>("Z1Beam");
0104     ziDip = args.value<double>("ZiDip");
0105     dzStep = args.value<double>("DzStep");
0106     moduleThick = args.value<double>("ModuleThick");
0107     layerThick = args.value<double>("LayerThick");
0108     scintThick = args.value<double>("ScintThick");
0109 #ifdef EDM_ML_DEBUG
0110     edm::LogVerbatim("HCalGeom") << "DDHCalEndcapModuleAlgo: Zmin " << cms::convert2mm(zMinBlock) << "\tZmax "
0111                                  << cms::convert2mm(zMaxBlock) << "\tZ1Beam " << cms::convert2mm(z1Beam) << "\tZiDip "
0112                                  << cms::convert2mm(ziDip) << "\tDzStep " << cms::convert2mm(dzStep) << "\tModuleThick "
0113                                  << cms::convert2mm(moduleThick) << "\tLayerThick " << cms::convert2mm(layerThick)
0114                                  << "\tScintThick " << cms::convert2mm(scintThick);
0115 #endif
0116     rMaxFront = args.value<double>("RMaxFront");
0117     rMaxBack = args.value<double>("RMaxBack");
0118     trimLeft = args.value<double>("TrimLeft");
0119     trimRight = args.value<double>("TrimRight");
0120     tolAbs = args.value<double>("TolAbs");
0121 #ifdef EDM_ML_DEBUG
0122     edm::LogVerbatim("HCalGeom") << "DDHCalEndcapModuleAlgo: RMaxFront " << cms::convert2mm(rMaxFront) << "\tRmaxBack "
0123                                  << cms::convert2mm(rMaxBack) << "\tTrims " << cms::convert2mm(trimLeft) << ":"
0124                                  << cms::convert2mm(trimRight) << "\tTolAbs " << cms::convert2mm(tolAbs);
0125 #endif
0126     slopeBot = args.value<double>("SlopeBottom");
0127     slopeTop = args.value<double>("SlopeTop");
0128     slopeTopF = args.value<double>("SlopeTopFront");
0129     modType = args.value<int>("ModType");
0130     modNumber = args.value<int>("ModNumber");
0131     layerType = args.value<int>("LayerType");
0132 #ifdef EDM_ML_DEBUG
0133     edm::LogVerbatim("HCalGeom") << "DDHCalEndcapModuleAlgo: slopeBot " << slopeBot << "\tslopeTop " << slopeTop
0134                                  << "\tslopeTopF " << slopeTopF << "\tmodType " << modType << "\tmodNumber "
0135                                  << modNumber << "\tlayerType " << layerType;
0136 #endif
0137     layerNumber = args.value<std::vector<int> >("LayerNumber");
0138 #ifdef EDM_ML_DEBUG
0139     edm::LogVerbatim("HCalGeom") << "DDHCalEndcapModuleAlgo: " << layerNumber.size() << " layer Numbers";
0140     for (unsigned int i = 0; i < layerNumber.size(); ++i)
0141       edm::LogVerbatim("HCalGeom") << "LayerNumber[" << i << "] = " << layerNumber[i];
0142 #endif
0143     phiName = args.value<std::vector<std::string> >("PhiName");
0144 #ifdef EDM_ML_DEBUG
0145     edm::LogVerbatim("HCalGeom") << "DDHCalEndcapModuleAlgo: " << phiName.size() << " phi sectors";
0146     for (unsigned int i = 0; i < phiName.size(); ++i)
0147       edm::LogVerbatim("HCalGeom") << "PhiName[" << i << "] = " << phiName[i];
0148 #endif
0149     layerName = args.value<std::vector<std::string> >("LayerName");
0150 #ifdef EDM_ML_DEBUG
0151     edm::LogVerbatim("HCalGeom") << "DDHCalEndcapModuleAlgo: " << layerName.size() << " layers";
0152     for (unsigned int i = 0; i < layerName.size(); ++i)
0153       edm::LogVerbatim("HCalGeom") << "LayerName[" << i << "] = " << layerName[i];
0154 #endif
0155     idName = args.value<std::string>("MotherName");
0156     idOffset = args.value<int>("IdOffset");
0157     modName = args.value<std::string>("ModName");
0158 #ifdef EDM_ML_DEBUG
0159     edm::LogVerbatim("HCalGeom") << "DDHCalEndcapModuleAlgo: Parent " << args.parentName() << "   " << modName
0160                                  << " idName " << idName << " NameSpace " << ns.name() << " Offset " << idOffset;
0161 #endif
0162 
0163 #ifdef EDM_ML_DEBUG
0164     edm::LogVerbatim("HCalGeom") << "==>> Constructing DDHCalEndcapModuleAlgo...";
0165 #endif
0166 
0167     dd4hep::Volume mother = ns.volume(args.parentName());
0168     if (modType == 0)
0169       constructInsideModule0(ctxt, e, mother);
0170     else
0171       constructInsideModule(ctxt, e, mother);
0172 #ifdef EDM_ML_DEBUG
0173     edm::LogVerbatim("HCalGeom") << "<<== End of DDHCalEndcapModuleAlgo construction ...";
0174 #endif
0175   }
0176 
0177   void constructInsideModule0(cms::DDParsingContext& ctxt, xml_h e, dd4hep::Volume& module) {
0178     cms::DDNamespace ns(ctxt, e, true);
0179 
0180 #ifdef EDM_ML_DEBUG
0181     edm::LogVerbatim("HCalGeom") << "DDHCalEndcapModuleAlgo: \t\tInside module0";
0182 #endif
0183     ///////////////////////////////////////////////////////////////
0184     //Pointers to the Materials
0185     dd4hep::Material matabsorbr = ns.material(absorberMat);
0186     dd4hep::Material matplastic = ns.material(plasticMat);
0187     dd4hep::Rotation3D rot = getRotation(rotstr, ns);
0188 
0189     int layer = layerNumber[0];
0190     int layer0 = layerNumber[1];
0191     std::string name;
0192     dd4hep::Solid solid;
0193     dd4hep::Volume glog;
0194     for (unsigned int iphi = 0; iphi < phiName.size(); iphi++) {
0195       HCalEndcapModuleAlgo::HcalEndcapPar parm = parameterLayer0(iphi);
0196       name = idName + modName + layerName[0] + phiName[iphi];
0197       solid = dd4hep::Trap(ns.prepend(name),
0198                            0.5 * layerThick,
0199                            0,
0200                            0,
0201                            parm.yh1,
0202                            parm.bl1,
0203                            parm.tl1,
0204                            parm.alp,
0205                            parm.yh2,
0206                            parm.bl1,
0207                            parm.tl2,
0208                            parm.alp);
0209 #ifdef EDM_ML_DEBUG
0210       edm::LogVerbatim("HCalGeom") << "DDHCalEndcapModuleAlgo: " << solid.name() << " Trap made of " << plasticMat
0211                                    << " of dimensions " << cms::convert2mm(0.5 * layerThick) << ", 0, 0, "
0212                                    << cms::convert2mm(parm.yh1) << ", " << cms::convert2mm(parm.bl1) << ", "
0213                                    << cms::convert2mm(parm.tl1) << ", " << convertRadToDeg(parm.alp) << ", "
0214                                    << cms::convert2mm(parm.yh2) << ", " << cms::convert2mm(parm.bl2) << ", "
0215                                    << cms::convert2mm(parm.tl2) << ", " << convertRadToDeg(parm.alp);
0216 #endif
0217       glog = dd4hep::Volume(solid.name(), solid, matplastic);
0218 
0219       dd4hep::Position r1(parm.xpos, parm.ypos, parm.zpos);
0220       module.placeVolume(glog, idOffset + layer + 1, dd4hep::Transform3D(rot, r1));
0221 #ifdef EDM_ML_DEBUG
0222       edm::LogVerbatim("HCalGeom") << "DDHCalEndcapModuleAlgo: " << glog.name() << " number " << (idOffset + layer + 1)
0223                                    << " positioned in " << module.name() << " at (" << cms::convert2mm(parm.xpos)
0224                                    << ", " << cms::convert2mm(parm.ypos) << ", " << cms::convert2mm(parm.zpos)
0225                                    << ") with rotation: " << rot;
0226 #endif
0227       //Now construct the layer of scintillator inside this
0228       int copyNo = layer0 * 10 + layerType;
0229       name = modName + layerName[0] + phiName[iphi];
0230       constructScintLayer(glog, scintThick, parm, name, copyNo, ns);
0231     }
0232 
0233     //Now the absorber layer
0234     double zi = zMinBlock + layerThick;
0235     double zo = zi + 0.5 * dzStep;
0236     double rinF, routF, rinB, routB;
0237     if (modNumber == 0) {
0238       rinF = zi * slopeTopF;
0239       routF = (zi - z1Beam) * slopeTop;
0240       rinB = zo * slopeTopF;
0241       routB = (zo - z1Beam) * slopeTop;
0242     } else if (modNumber > 0) {
0243       rinF = zi * slopeBot;
0244       routF = zi * slopeTopF;
0245       rinB = zo * slopeBot;
0246       routB = zo * slopeTopF;
0247     } else {
0248       rinF = zi * slopeBot;
0249       routF = (zi - z1Beam) * slopeTop;
0250       rinB = zo * slopeBot;
0251       routB = (zo - z1Beam) * slopeTop;
0252     }
0253 #ifdef EDM_ML_DEBUG
0254     edm::LogVerbatim("HCalGeom") << "DDHCalEndcapModuleAlgo: Front " << cms::convert2mm(zi) << ", "
0255                                  << cms::convert2mm(rinF) << ", " << cms::convert2mm(routF) << " Back "
0256                                  << cms::convert2mm(zo) << ", " << cms::convert2mm(rinB) << ", "
0257                                  << cms::convert2mm(routB);
0258 #endif
0259     HCalEndcapModuleAlgo::HcalEndcapPar parm = parameterLayer(0, rinF, routF, rinB, routB, zi, zo);
0260 #ifdef EDM_ML_DEBUG
0261     edm::LogVerbatim("HCalGeom") << "DDHCalEndcapModuleAlgo: Trim " << cms::convert2mm(tolAbs) << " Param "
0262                                  << cms::convert2mm(parm.yh1) << ", " << cms::convert2mm(parm.bl1) << ", "
0263                                  << cms::convert2mm(parm.tl1) << ", " << cms::convert2mm(parm.yh2) << ", "
0264                                  << cms::convert2mm(parm.bl2) << ", " << cms::convert2mm(parm.tl2);
0265 #endif
0266     parm.bl1 -= tolAbs;
0267     parm.tl1 -= tolAbs;
0268     parm.bl2 -= tolAbs;
0269     parm.tl2 -= tolAbs;
0270 
0271     name = idName + modName + layerName[0] + "Absorber";
0272     solid = dd4hep::Trap(ns.prepend(name),
0273                          0.5 * moduleThick,
0274                          parm.theta,
0275                          parm.phi,
0276                          parm.yh1,
0277                          parm.bl1,
0278                          parm.tl1,
0279                          parm.alp,
0280                          parm.yh2,
0281                          parm.bl2,
0282                          parm.tl2,
0283                          parm.alp);
0284 #ifdef EDM_ML_DEBUG
0285     edm::LogVerbatim("HCalGeom") << "DDHCalEndcapModuleAlgo: " << solid.name() << " Trap made of " << matabsorbr.name()
0286                                  << " of dimensions " << cms::convert2mm(0.5 * moduleThick) << ", "
0287                                  << convertRadToDeg(parm.theta) << ", " << convertRadToDeg(parm.phi) << ", "
0288                                  << cms::convert2mm(parm.yh1) << ", " << cms::convert2mm(parm.bl1) << ", "
0289                                  << cms::convert2mm(parm.tl1) << ", " << convertRadToDeg(parm.alp) << ", "
0290                                  << cms::convert2mm(parm.yh2) << ", " << cms::convert2mm(parm.bl2) << ", "
0291                                  << cms::convert2mm(parm.tl2) << ", " << convertRadToDeg(parm.alp);
0292 #endif
0293     glog = dd4hep::Volume(solid.name(), solid, matabsorbr);
0294 
0295     dd4hep::Position r2(parm.xpos, parm.ypos, parm.zpos);
0296     module.placeVolume(glog, idOffset + layer + 1, dd4hep::Transform3D(rot, r2));
0297 #ifdef EDM_ML_DEBUG
0298     edm::LogVerbatim("HCalGeom") << "DDHCalEndcapModuleAlgo: " << glog.name() << " number 1 positioned in "
0299                                  << module.name() << " at (" << cms::convert2mm(parm.xpos) << ", "
0300                                  << cms::convert2mm(parm.ypos) << ", " << cms::convert2mm(parm.zpos)
0301                                  << ") with rotation: " << rot;
0302 #endif
0303   }
0304 
0305   void constructInsideModule(cms::DDParsingContext& ctxt, xml_h e, dd4hep::Volume& module) {
0306     cms::DDNamespace ns(ctxt, e, true);
0307 #ifdef EDM_ML_DEBUG
0308     edm::LogVerbatim("HCalGeom") << "DDHCalEndcapModuleAlgo: \t\tInside module";
0309 #endif
0310     ///////////////////////////////////////////////////////////////
0311     //Pointers to the Rotation Matrices and to the Materials
0312     dd4hep::Material matter = ns.material(genMaterial);
0313     dd4hep::Material matplastic = ns.material(plasticMat);
0314     dd4hep::Rotation3D rot = getRotation(rotstr, ns);
0315 
0316     double alpha = (1._pi) / sectors;
0317     double zi = zMinBlock;
0318 
0319     for (unsigned int i = 0; i < layerName.size(); i++) {
0320       std::string name;
0321       dd4hep::Solid solid;
0322       dd4hep::Volume glog, plog;
0323       int layer = layerNumber[i];
0324       double zo = zi + 0.5 * dzStep;
0325 
0326       for (unsigned int iphi = 0; iphi < phiName.size(); iphi++) {
0327         double ziAir = zo - moduleThick;
0328         double rinF, rinB;
0329         if (modNumber == 0) {
0330           rinF = ziAir * slopeTopF;
0331           rinB = zo * slopeTopF;
0332         } else {
0333           rinF = ziAir * slopeBot;
0334           rinB = zo * slopeBot;
0335         }
0336         double routF = getRout(ziAir);
0337         double routB = getRout(zo);
0338 #ifdef EDM_ML_DEBUG
0339         edm::LogVerbatim("HCalGeom") << "DDHCalEndcapModuleAlgo: Layer " << i << " Phi " << iphi << " Front "
0340                                      << cms::convert2mm(ziAir) << ", " << cms::convert2mm(rinF) << ", "
0341                                      << cms::convert2mm(routF) << " Back " << cms::convert2mm(zo) << ", "
0342                                      << cms::convert2mm(rinB) << ", " << cms::convert2mm(routB);
0343 #endif
0344         HCalEndcapModuleAlgo::HcalEndcapPar parm = parameterLayer(iphi, rinF, routF, rinB, routB, ziAir, zo);
0345 
0346         name = idName + modName + layerName[i] + phiName[iphi] + "Air";
0347         solid = dd4hep::Trap(ns.prepend(name),
0348                              0.5 * moduleThick,
0349                              parm.theta,
0350                              parm.phi,
0351                              parm.yh1,
0352                              parm.bl1,
0353                              parm.tl1,
0354                              parm.alp,
0355                              parm.yh2,
0356                              parm.bl2,
0357                              parm.tl2,
0358                              parm.alp);
0359 #ifdef EDM_ML_DEBUG
0360         edm::LogVerbatim("HCalGeom") << "DDHCalEndcapModuleAlgo: " << solid.name() << " Trap made of " << matter.name()
0361                                      << " of dimensions " << cms::convert2mm(0.5 * moduleThick) << ", "
0362                                      << convertRadToDeg(parm.theta) << ", " << convertRadToDeg(parm.phi) << ", "
0363                                      << cms::convert2mm(parm.yh1) << ", " << cms::convert2mm(parm.bl1) << ", "
0364                                      << cms::convert2mm(parm.tl1) << ", " << convertRadToDeg(parm.alp) << ", "
0365                                      << cms::convert2mm(parm.yh2) << ", " << cms::convert2mm(parm.bl2) << ", "
0366                                      << cms::convert2mm(parm.tl2) << ", " << convertRadToDeg(parm.alp);
0367 #endif
0368         glog = dd4hep::Volume(solid.name(), solid, matter);
0369 
0370         dd4hep::Position r1(parm.xpos, parm.ypos, parm.zpos);
0371         module.placeVolume(glog, layer + 1, dd4hep::Transform3D(rot, r1));
0372 #ifdef EDM_ML_DEBUG
0373         edm::LogVerbatim("HCalGeom") << "DDHCalEndcapModuleAlgo: " << glog.name() << " number " << (layer + 1)
0374                                      << " positioned in " << module.name() << " at (" << cms::convert2mm(parm.xpos)
0375                                      << ", " << cms::convert2mm(parm.ypos) << ", " << cms::convert2mm(parm.zpos)
0376                                      << ") with rotation: " << rot;
0377 #endif
0378         //Now the plastic with scintillators
0379         parm.yh1 = 0.5 * (routF - rinB) - getTrim(iphi);
0380         parm.bl1 = 0.5 * rinB * tan(alpha) - getTrim(iphi);
0381         parm.tl1 = 0.5 * routF * tan(alpha) - getTrim(iphi);
0382         name = idName + modName + layerName[i] + phiName[iphi];
0383         solid = dd4hep::Trap(ns.prepend(name),
0384                              0.5 * layerThick,
0385                              0,
0386                              0,
0387                              parm.yh1,
0388                              parm.bl1,
0389                              parm.tl1,
0390                              parm.alp,
0391                              parm.yh1,
0392                              parm.bl1,
0393                              parm.tl1,
0394                              parm.alp);
0395 #ifdef EDM_ML_DEBUG
0396         edm::LogVerbatim("HCalGeom") << "DDHCalEndcapModuleAlgo: " << solid.name() << " Trap made of "
0397                                      << matplastic.name() << " of dimensions " << cms::convert2mm(0.5 * layerThick)
0398                                      << ", 0, 0, " << cms::convert2mm(parm.yh1) << ", " << cms::convert2mm(parm.bl1)
0399                                      << ", " << cms::convert2mm(parm.tl1) << ", " << convertRadToDeg(parm.alp) << ", "
0400                                      << cms::convert2mm(parm.yh1) << ", " << cms::convert2mm(parm.bl1) << ", "
0401                                      << cms::convert2mm(parm.tl1) << ", " << convertRadToDeg(parm.alp);
0402 #endif
0403         plog = dd4hep::Volume(solid.name(), solid, matplastic);
0404 
0405         double ypos = 0.5 * (routF + rinB) - parm.xpos;
0406         dd4hep::Position r2(0., ypos, 0.);
0407         glog.placeVolume(plog, idOffset + layer + 1, r2);
0408 #ifdef EDM_ML_DEBUG
0409         edm::LogVerbatim("HCalGeom") << "DDHCalEndcapModuleAlgo: " << plog.name() << " number "
0410                                      << (idOffset + layer + 1) << " positioned in " << glog.name() << " at (0, "
0411                                      << cms::convert2mm(ypos) << ", 0) with no rotation";
0412 #endif
0413         //Constructin the scintillators inside
0414         int copyNo = layer * 10 + layerType;
0415         name = modName + layerName[i] + phiName[iphi];
0416         constructScintLayer(plog, scintThick, parm, name, copyNo, ns);
0417         zo += 0.5 * dzStep;
0418       }  // End of loop over phi indices
0419       zi = zo - 0.5 * dzStep;
0420     }  // End of loop on layers
0421   }
0422 
0423   HcalEndcapPar parameterLayer0(unsigned int iphi) {
0424     HCalEndcapModuleAlgo::HcalEndcapPar parm;
0425     //Given module and layer number compute parameters of trapezoid
0426     //and positioning parameters
0427     double alpha = (1._pi) / sectors;
0428 #ifdef EDM_ML_DEBUG
0429     edm::LogVerbatim("HCalGeom") << "Input " << iphi << " Alpha " << convertRadToDeg(alpha);
0430 #endif
0431     double zi, zo;
0432     if (iphi == 0) {
0433       zi = zMinBlock;
0434       zo = zi + layerThick;
0435     } else {
0436       zo = zMaxBlock;
0437       zi = zo - layerThick;
0438     }
0439     double rin, rout;
0440     if (modNumber == 0) {
0441       rin = zo * slopeTopF;
0442       rout = (zi - z1Beam) * slopeTop;
0443     } else if (modNumber > 0) {
0444       rin = zo * slopeBot;
0445       rout = zi * slopeTopF;
0446     } else {
0447       rin = zo * slopeBot;
0448       rout = (zi - z1Beam) * slopeTop;
0449     }
0450 #ifdef EDM_ML_DEBUG
0451     edm::LogVerbatim("HCalGeom") << "ModNumber " << modNumber << " " << cms::convert2mm(zi) << " "
0452                                  << cms::convert2mm(zo) << " " << slopeTopF << " " << slopeTop << " " << slopeBot << " "
0453                                  << cms::convert2mm(rin) << " " << cms::convert2mm(rout) << " "
0454                                  << cms::convert2mm(getTrim(iphi));
0455 #endif
0456     double yh = 0.5 * (rout - rin);
0457     double bl = 0.5 * rin * tan(alpha);
0458     double tl = 0.5 * rout * tan(alpha);
0459     parm.xpos = 0.5 * (rin + rout);
0460     parm.ypos = 0.5 * (bl + tl);
0461     parm.zpos = 0.5 * (zi + zo);
0462     parm.yh1 = parm.yh2 = yh - getTrim(iphi);
0463     parm.bl1 = parm.bl2 = bl - getTrim(iphi);
0464     parm.tl1 = parm.tl2 = tl - getTrim(iphi);
0465     parm.alp = atan(0.5 * tan(alpha));
0466     if (iphi == 0) {
0467       parm.ypos = -parm.ypos;
0468     } else {
0469       parm.alp = -parm.alp;
0470     }
0471 #ifdef EDM_ML_DEBUG
0472     edm::LogVerbatim("HCalGeom") << "Output Dimensions " << cms::convert2mm(parm.yh1) << " "
0473                                  << cms::convert2mm(parm.bl1) << " " << cms::convert2mm(parm.tl1) << " "
0474                                  << convertRadToDeg(parm.alp) << " Position " << cms::convert2mm(parm.xpos) << " "
0475                                  << cms::convert2mm(parm.ypos) << " " << cms::convert2mm(parm.zpos);
0476 #endif
0477     return parm;
0478   }
0479 
0480   HcalEndcapPar parameterLayer(
0481       unsigned int iphi, double rinF, double routF, double rinB, double routB, double zi, double zo) {
0482     HCalEndcapModuleAlgo::HcalEndcapPar parm;
0483     //Given rin, rout compute parameters of the trapezoid and
0484     //position of the trapezoid for a standrd layer
0485     double alpha = (1._pi) / sectors;
0486 #ifdef EDM_ML_DEBUG
0487     edm::LogVerbatim("HCalGeom") << "Input " << iphi << " Front " << cms::convert2mm(rinF) << " "
0488                                  << cms::convert2mm(routF) << " " << cms::convert2mm(zi) << " Back "
0489                                  << cms::convert2mm(rinB) << " " << cms::convert2mm(routB) << " " << cms::convert2mm(zo)
0490                                  << " Alpha " << convertRadToDeg(alpha);
0491 #endif
0492     parm.yh1 = 0.5 * (routF - rinB);
0493     parm.bl1 = 0.5 * rinB * tan(alpha);
0494     parm.tl1 = 0.5 * routF * tan(alpha);
0495     parm.yh2 = 0.5 * (routF - rinB);
0496     parm.bl2 = 0.5 * rinB * tan(alpha);
0497     parm.tl2 = 0.5 * routF * tan(alpha);
0498     double dx = 0.25 * (parm.bl2 + parm.tl2 - parm.bl1 - parm.tl1);
0499     double dy = 0.5 * (rinB + routF - rinB - routF);
0500     parm.xpos = 0.25 * (rinB + routF + rinB + routF);
0501     parm.ypos = 0.25 * (parm.bl2 + parm.tl2 + parm.bl1 + parm.tl1);
0502     parm.zpos = 0.5 * (zi + zo);
0503     parm.alp = atan(0.5 * tan(alpha));
0504     if (iphi == 0) {
0505       parm.ypos = -parm.ypos;
0506     } else {
0507       parm.alp = -parm.alp;
0508       dx = -dx;
0509     }
0510     double r = sqrt(dx * dx + dy * dy);
0511 #ifdef EDM_ML_DEBUG
0512     edm::LogVerbatim("HCalGeom") << "dx|dy|r " << cms::convert2mm(dx) << ":" << cms::convert2mm(dy) << ":"
0513                                  << cms::convert2mm(r);
0514 #endif
0515     if (r > 1.0e-8) {
0516       parm.theta = atan(r / (zo - zi));
0517       parm.phi = atan2(dy, dx);
0518     } else {
0519       parm.theta = parm.phi = 0;
0520     }
0521 #ifdef EDM_ML_DEBUG
0522     edm::LogVerbatim("HCalGeom") << "Output Dimensions " << cms::convert2mm(parm.yh1) << " "
0523                                  << cms::convert2mm(parm.bl1) << " " << cms::convert2mm(parm.tl1) << " "
0524                                  << cms::convert2mm(parm.yh2) << " " << cms::convert2mm(parm.bl2) << " "
0525                                  << cms::convert2mm(parm.tl2) << " " << convertRadToDeg(parm.alp) << " "
0526                                  << convertRadToDeg(parm.theta) << " " << convertRadToDeg(parm.phi) << " Position "
0527                                  << cms::convert2mm(parm.xpos) << " " << cms::convert2mm(parm.ypos) << " "
0528                                  << cms::convert2mm(parm.zpos);
0529 #endif
0530     return parm;
0531   }
0532 
0533   void constructScintLayer(dd4hep::Volume& detector,
0534                            double dz,
0535                            HCalEndcapModuleAlgo::HcalEndcapPar parm,
0536                            const std::string& nm,
0537                            int id,
0538                            cms::DDNamespace& ns) {
0539     dd4hep::Material matter = ns.material(scintMat);
0540     std::string name = idName + "Scintillator" + nm;
0541 
0542     dd4hep::Solid solid = dd4hep::Trap(
0543         ns.prepend(name), 0.5 * dz, 0, 0, parm.yh1, parm.bl1, parm.tl1, parm.alp, parm.yh1, parm.bl1, parm.tl1, parm.alp);
0544 #ifdef EDM_ML_DEBUG
0545     edm::LogVerbatim("HCalGeom") << "DDHCalEndcapModuleAlgo: " << solid.name() << " Trap made of " << scintMat
0546                                  << " of dimensions " << cms::convert2mm(0.5 * dz) << ", 0, 0, "
0547                                  << cms::convert2mm(parm.yh1) << ", " << cms::convert2mm(parm.bl1) << ", "
0548                                  << cms::convert2mm(parm.tl1) << ", " << convertRadToDeg(parm.alp) << ", "
0549                                  << cms::convert2mm(parm.yh1) << ", " << cms::convert2mm(parm.bl1) << ", "
0550                                  << cms::convert2mm(parm.tl1) << ", " << convertRadToDeg(parm.alp);
0551 #endif
0552     dd4hep::Volume glog(solid.name(), solid, matter);
0553 
0554     detector.placeVolume(glog, id);
0555 #ifdef EDM_ML_DEBUG
0556     edm::LogVerbatim("HCalGeom") << "DDHCalEndcapModuleAlgo: " << glog.name() << " number " << id << " positioned in "
0557                                  << detector.name() << " at (0,0,0) with no rotation";
0558 #endif
0559   }
0560 
0561   double getTrim(unsigned int j) const {
0562     if (j == 0)
0563       return trimLeft;
0564     else
0565       return trimRight;
0566   }
0567 
0568   double getRout(double z) const {
0569     double r = (modNumber >= 0) ? ((z - z1Beam) * slopeTop) : z * slopeTopF;
0570     if (z > ziDip) {
0571       if (r > rMaxBack)
0572         r = rMaxBack;
0573     } else {
0574       if (r > rMaxFront)
0575         r = rMaxFront;
0576     }
0577     return r;
0578   }
0579 
0580   dd4hep::Rotation3D getRotation(const std::string& rotstr, cms::DDNamespace& ns) {
0581     std::string rot = (strchr(rotstr.c_str(), NAMESPACE_SEP) == nullptr) ? ("rotations:" + rotstr) : rotstr;
0582 #ifdef EDM_ML_DEBUG
0583     edm::LogVerbatim("HCalGeom") << "getRotation: " << rotstr << ":" << rot << ":" << ns.rotation(rot);
0584 #endif
0585     return ns.rotation(rot);
0586   }
0587 };
0588 
0589 static long algorithm(dd4hep::Detector& /* description */, cms::DDParsingContext& ctxt, xml_h e) {
0590   HCalEndcapModuleAlgo hcalendcapalgo(ctxt, e);
0591   return cms::s_executed;
0592 }
0593 
0594 // first argument is the type from the xml file
0595 DECLARE_DDCMS_DETELEMENT(DDCMS_hcal_DDHCalEndcapModuleAlgo, algorithm)