Back to home page

Project CMSSW displayed by LXR

 
 

    


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

0001 ///////////////////////////////////////////////////////////////////////////////
0002 // File: DDHCalEndcapAlgo.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/Core/interface/DDSplit.h"
0016 #include "DetectorDescription/DDCMS/interface/DDPlugins.h"
0017 #include "DetectorDescription/DDCMS/interface/DDutils.h"
0018 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0019 
0020 //#define EDM_ML_DEBUG
0021 using namespace angle_units::operators;
0022 
0023 struct HCalEndcapAlgo {
0024   std::string genMaterial;             //General material
0025   int nsectors;                        //Number of potenital straight edges
0026   int nsectortot;                      //Number of straight edges (actual)
0027   int nEndcap;                         //Number of endcaps
0028   std::vector<int> eModule;            //Modules to be present in part i (?)
0029   std::string rotHalf;                 //Rotation matrix for half
0030   std::string rotation;                //Rotation matrix to place in mother
0031   double zShift;                       //needed for TB setup (move HE)
0032   double zFront;                       //Z of the front section
0033   double zEnd;                         //Outer Z of the HE
0034   double ziNose;                       //Starting Z of the nose
0035   double ziL0Nose;                     //Starting Z of layer 0 at nose
0036   double ziBody;                       //Starting Z of the body
0037   double ziL0Body;                     //Starting Z of layer 0 at body
0038   double ziKink;                       //Position of the kink point
0039   double z0Beam;                       //Position of gap front along z-axis
0040   double z1Beam;                       //Position of gap end   along z-axis
0041   double ziDip;                        //Starting Z of dipped part of body
0042   double dzStep;                       //Width in Z of a layer
0043   double dzShift;                      //Shift in Z for HE
0044   double zShiftHac2;                   //needed for TB (remove part Hac2)
0045   double rout;                         //Outer R of the HE
0046   double riKink;                       //Inner radius at kink point
0047   double riDip;                        //Inner radius at the dip point
0048   double roDip;                        //Outer radius at the dip point
0049   double heboxDepth;                   //Depth of the HE box
0050   double drEnd;                        //Shift in R for the end absorber
0051   double angTop;                       //Angle of top end of HE
0052   double angBot;                       //Angle of the bottom end of HE
0053   double angGap;                       //Gap angle (in degrees)
0054   double slope;                        //Slope of the gap on HE side
0055   std::string absMat;                  //Absorber     material
0056   int modules;                         //Number of modules
0057   std::vector<std::string> modName;    //Name
0058   std::vector<std::string> modMat;     //Material
0059   std::vector<int> modType;            //Type (0/1 for front/standard)
0060   std::vector<int> sectionModule;      //Number of sections in a module
0061   std::vector<int> layerN;             //Number of layers
0062   std::vector<int> layerN0;            //Layer numbers in section 0
0063   std::vector<int> layerN1;            //Layer numbers in section 1
0064   std::vector<int> layerN2;            //Layer numbers in section 2
0065   std::vector<int> layerN3;            //Layer numbers in section 3
0066   std::vector<int> layerN4;            //Layer numbers in section 4
0067   std::vector<int> layerN5;            //Layer numbers in section 5
0068   std::vector<double> thick;           //Thickness of absorber/air
0069   std::vector<double> trimLeft;        //Trimming of left  layers in module
0070   std::vector<double> trimRight;       //Trimming of right layers in module
0071   std::vector<double> zminBlock;       //Minimum Z
0072   std::vector<double> zmaxBlock;       //Maximum Z
0073   std::vector<double> rinBlock1;       //Inner Radius
0074   std::vector<double> routBlock1;      //Outer Radius at zmin
0075   std::vector<double> rinBlock2;       //Inner Radius
0076   std::vector<double> routBlock2;      //Outer Radius at zmax
0077   int phiSections;                     //Number of phi sections
0078   std::vector<std::string> phiName;    //Name of Phi sections
0079   int layers;                          //Number of layers
0080   std::vector<std::string> layerName;  //Layer Names
0081   std::vector<int> layerType;          //Detector type in each layer
0082   std::vector<double> layerT;          //Layer thickness (plastic + scint.)
0083   std::vector<double> scintT;          //Scintillator thickness
0084   std::string plastMat;                //Plastic      material
0085   std::string scintMat;                //Scintillator material
0086   std::string rotmat;                  //Rotation matrix for positioning
0087   std::string idName;                  //Name of the "parent" volume.
0088   int idOffset;                        // Geant4 ID's...    = 4000;
0089   double tolPos, tolAbs;               //Tolerances
0090 
0091   HCalEndcapAlgo() = delete;
0092 
0093   HCalEndcapAlgo(cms::DDParsingContext& ctxt, xml_h e) {
0094     cms::DDNamespace ns(ctxt, e, true);
0095     cms::DDAlgoArguments args(ctxt, e);
0096 
0097     genMaterial = args.value<std::string>("MaterialName");
0098     rotation = args.value<std::string>("Rotation");
0099     nsectors = args.value<int>("Sector");
0100     nsectortot = args.value<int>("SectorTot");
0101     nEndcap = args.value<int>("Endcap");
0102     rotHalf = args.value<std::string>("RotHalf");
0103     zShift = args.value<double>("ZShift");
0104     zFront = args.value<double>("ZFront");
0105     zEnd = args.value<double>("ZEnd");
0106     ziNose = args.value<double>("ZiNose");
0107     ziL0Nose = args.value<double>("ZiL0Nose");
0108     ziBody = args.value<double>("ZiBody");
0109     ziL0Body = args.value<double>("ZiL0Body");
0110     z0Beam = args.value<double>("Z0Beam");
0111     ziDip = args.value<double>("ZiDip");
0112     dzStep = args.value<double>("DzStep");
0113     zShiftHac2 = args.value<double>("ZShiftHac2");
0114     double gap = args.value<double>("Gap");
0115     double z1 = args.value<double>("Z1");
0116     double r1 = args.value<double>("R1");
0117     rout = args.value<double>("Rout");
0118     heboxDepth = args.value<double>("HEboxDepth");
0119     drEnd = args.value<double>("DrEnd");
0120     double etamin = args.value<double>("Etamin");
0121     angBot = args.value<double>("AngBot");
0122     angGap = args.value<double>("AngGap");
0123 #ifdef EDM_ML_DEBUG
0124     edm::LogVerbatim("HCalGeom") << "DDHCalEndcapAlgo: General material " << genMaterial << "\tSectors " << nsectors
0125                                  << ",  " << nsectortot << "\tEndcaps " << nEndcap << "\tRotation matrix for half "
0126                                  << rotHalf << "\n\tzFront " << cms::convert2mm(zFront) << " zEnd "
0127                                  << cms::convert2mm(zEnd) << " ziNose " << cms::convert2mm(ziNose) << " ziL0Nose "
0128                                  << cms::convert2mm(ziL0Nose) << " ziBody " << cms::convert2mm(ziBody) << " ziL0Body "
0129                                  << cms::convert2mm(ziL0Body) << " z0Beam " << cms::convert2mm(z0Beam) << " ziDip "
0130                                  << cms::convert2mm(ziDip) << " dzStep " << cms::convert2mm(dzStep) << " Gap "
0131                                  << cms::convert2mm(gap) << " z1 " << cms::convert2mm(z1) << "\n\tr1 "
0132                                  << cms::convert2mm(r1) << " rout " << cms::convert2mm(rout) << " HeboxDepth "
0133                                  << cms::convert2mm(heboxDepth) << " drEnd " << cms::convert2mm(drEnd) << "\tetamin "
0134                                  << etamin << " Bottom angle " << angBot << " Gap angle " << angGap << " Z-Shift "
0135                                  << cms::convert2mm(zShift) << " " << cms::convert2mm(zShiftHac2);
0136 #endif
0137 
0138     //Derived quantities
0139     angTop = 2.0 * atan(exp(-etamin));
0140     slope = tan(angGap);
0141     z1Beam = z1 - r1 / slope;
0142     ziKink = z1Beam + rout / slope;
0143     riKink = ziKink * tan(angBot);
0144     riDip = ziDip * tan(angBot);
0145     roDip = rout - heboxDepth;
0146     dzShift = (z1Beam - z0Beam) - gap / sin(angGap);
0147 #ifdef EDM_ML_DEBUG
0148     edm::LogVerbatim("HCalGeom") << "DDHCalEndcapAlgo: angTop " << convertRadToDeg(angTop) << "\tSlope " << slope
0149                                  << "\tDzShift " << cms::convert2mm(dzShift) << "\n\tz1Beam " << cms::convert2mm(z1Beam)
0150                                  << "\tziKink" << cms::convert2mm(ziKink) << "\triKink " << cms::convert2mm(riKink)
0151                                  << "\triDip " << cms::convert2mm(riDip) << "\n\troDip " << cms::convert2mm(roDip)
0152                                  << "\tRotation " << rotation;
0153 #endif
0154 
0155     ///////////////////////////////////////////////////////////////
0156     //Modules
0157     absMat = args.value<std::string>("AbsMat");
0158     modules = args.value<int>("Modules");
0159 #ifdef EDM_ML_DEBUG
0160     edm::LogVerbatim("HCalGeom") << "DDHCalEndcapAlgo: Number of modules " << modules << " and absorber material "
0161                                  << absMat;
0162 #endif
0163 
0164     modName = args.value<std::vector<std::string> >("ModuleName");
0165     modMat = args.value<std::vector<std::string> >("ModuleMat");
0166     modType = args.value<std::vector<int> >("ModuleType");
0167     sectionModule = args.value<std::vector<int> >("SectionModule");
0168     thick = args.value<std::vector<double> >("ModuleThick");
0169     trimLeft = args.value<std::vector<double> >("TrimLeft");
0170     trimRight = args.value<std::vector<double> >("TrimRight");
0171     eModule = args.value<std::vector<int> >("EquipModule");
0172     layerN = args.value<std::vector<int> >("LayerN");
0173     layerN0 = args.value<std::vector<int> >("LayerN0");
0174     layerN1 = args.value<std::vector<int> >("LayerN1");
0175     layerN2 = args.value<std::vector<int> >("LayerN2");
0176     layerN3 = args.value<std::vector<int> >("LayerN3");
0177     layerN4 = args.value<std::vector<int> >("LayerN4");
0178     layerN5 = args.value<std::vector<int> >("LayerN5");
0179 #ifdef EDM_ML_DEBUG
0180     for (int i = 0; i < modules; i++) {
0181       edm::LogVerbatim("HCalGeom") << "DDHCalEndcapAlgo: " << modName[i] << " type " << modType[i] << " Sections "
0182                                    << sectionModule[i] << " thickness of absorber/air " << cms::convert2mm(thick[i])
0183                                    << " trim " << cms::convert2mm(trimLeft[i]) << ", " << cms::convert2mm(trimRight[i])
0184                                    << " equip module " << eModule[i] << " with " << layerN[i] << " layers";
0185       if (i == 0) {
0186         for (int j = 0; j < layerN[i]; j++) {
0187           edm::LogVerbatim("HCalGeom") << "\t " << layerN0[j] << "/" << layerN0[j + 1];
0188         }
0189       } else if (i == 1) {
0190         for (int j = 0; j < layerN[i]; j++) {
0191           edm::LogVerbatim("HCalGeom") << "\t " << layerN1[j] << "/" << layerN1[j + 1];
0192         }
0193       } else if (i == 2) {
0194         for (int j = 0; j < layerN[i]; j++) {
0195           edm::LogVerbatim("HCalGeom") << "\t " << layerN2[j];
0196         }
0197       } else if (i == 3) {
0198         for (int j = 0; j < layerN[i]; j++) {
0199           edm::LogVerbatim("HCalGeom") << "\t " << layerN3[j];
0200         }
0201       } else if (i == 4) {
0202         for (int j = 0; j < layerN[i]; j++) {
0203           edm::LogVerbatim("HCalGeom") << "\t " << layerN4[j];
0204         }
0205       } else if (i == 5) {
0206         for (int j = 0; j < layerN[i]; j++) {
0207           edm::LogVerbatim("HCalGeom") << "\t " << layerN5[j];
0208         }
0209       }
0210     }
0211 #endif
0212 
0213     ///////////////////////////////////////////////////////////////
0214     //Layers
0215     phiSections = args.value<int>("PhiSections");
0216     phiName = args.value<std::vector<std::string> >("PhiName");
0217     layers = args.value<int>("Layers");
0218     layerName = args.value<std::vector<std::string> >("LayerName");
0219     layerType = args.value<std::vector<int> >("LayerType");
0220     layerT = args.value<std::vector<double> >("LayerT");
0221     scintT = args.value<std::vector<double> >("ScintT");
0222     scintMat = args.value<std::string>("ScintMat");
0223     plastMat = args.value<std::string>("PlastMat");
0224     rotmat = args.value<std::string>("RotMat");
0225 #ifdef EDM_ML_DEBUG
0226     edm::LogVerbatim("HCalGeom") << "DDHCalEndcapAlgo: Phi Sections " << phiSections;
0227     for (int i = 0; i < phiSections; i++)
0228       edm::LogVerbatim("HCalGeom") << "\tName[" << i << "] : " << phiName[i];
0229     edm::LogVerbatim("HCalGeom") << "\tPlastic: " << plastMat << "\tScintillator: " << scintMat << "\tRotation matrix "
0230                                  << rotmat << "\n\tNumber of layers " << layers;
0231     for (int i = 0; i < layers; i++) {
0232       edm::LogVerbatim("HCalGeom") << "\t" << layerName[i] << "\tType " << layerType[i] << "\tThickness "
0233                                    << cms::convert2mm(layerT[i]) << "\tScint.Thick " << cms::convert2mm(scintT[i]);
0234     }
0235 #endif
0236 
0237     ///////////////////////////////////////////////////////////////
0238     // Derive bounding of the modules
0239     int module = 0;
0240     // Layer 0 (Nose)
0241     if (modules > 0) {
0242       zminBlock.emplace_back(ziL0Nose);
0243       zmaxBlock.emplace_back(zminBlock[module] + layerT[0] + 0.5 * dzStep);
0244       rinBlock1.emplace_back(zminBlock[module] * tan(angTop));
0245       rinBlock2.emplace_back(zmaxBlock[module] * tan(angTop));
0246       routBlock1.emplace_back((zminBlock[module] - z1Beam) * slope);
0247       routBlock2.emplace_back((zmaxBlock[module] - z1Beam) * slope);
0248       module++;
0249     }
0250     // Layer 0 (Body)
0251     if (modules > 1) {
0252       zminBlock.emplace_back(ziL0Body);
0253       zmaxBlock.emplace_back(zminBlock[module] + layerT[0] + 0.5 * dzStep);
0254       rinBlock1.emplace_back(zminBlock[module] * tan(angBot));
0255       rinBlock2.emplace_back(zmaxBlock[module] * tan(angBot));
0256       routBlock1.emplace_back(zminBlock[module] * tan(angTop));
0257       routBlock2.emplace_back(zmaxBlock[module] * tan(angTop));
0258       module++;
0259     }
0260     // Hac1
0261     if (modules > 2) {
0262       zminBlock.emplace_back(ziNose);
0263       zmaxBlock.emplace_back(ziBody);
0264       rinBlock1.emplace_back(zminBlock[module] * tan(angTop));
0265       rinBlock2.emplace_back(zmaxBlock[module] * tan(angTop));
0266       routBlock1.emplace_back((zminBlock[module] - z1Beam) * slope);
0267       routBlock2.emplace_back((zmaxBlock[module] - z1Beam) * slope);
0268       module++;
0269     }
0270     // Hac2
0271     if (modules > 3) {
0272       zminBlock.emplace_back(ziBody);
0273       zmaxBlock.emplace_back(zminBlock[module] + layerN[3] * dzStep);
0274       rinBlock1.emplace_back(zminBlock[module] * tan(angBot));
0275       rinBlock2.emplace_back(zmaxBlock[module] * tan(angBot));
0276       routBlock1.emplace_back((zmaxBlock[module - 1] - z1Beam) * slope);
0277       routBlock2.emplace_back(rout);
0278       module++;
0279     }
0280     // Hac3
0281     if (modules > 4) {
0282       zminBlock.emplace_back(zmaxBlock[module - 1]);
0283       zmaxBlock.emplace_back(zminBlock[module] + layerN[4] * dzStep);
0284       rinBlock1.emplace_back(zminBlock[module] * tan(angBot));
0285       rinBlock2.emplace_back(zmaxBlock[module] * tan(angBot));
0286       routBlock1.emplace_back(rout);
0287       routBlock2.emplace_back(rout);
0288       module++;
0289     }
0290     // Hac4
0291     if (modules > 5) {
0292       zminBlock.emplace_back(zmaxBlock[module - 1]);
0293       zmaxBlock.emplace_back(zminBlock[module] + layerN[5] * dzStep);
0294       rinBlock1.emplace_back(zminBlock[module] * tan(angBot));
0295       rinBlock2.emplace_back(zmaxBlock[module] * tan(angBot));
0296       routBlock1.emplace_back(rout);
0297       routBlock2.emplace_back(roDip);
0298       module++;
0299     }
0300 #ifdef EDM_ML_DEBUG
0301     for (int i = 0; i < module; i++)
0302       edm::LogVerbatim("HCalGeom") << "DDHCalEndcapAlgo: Module " << i << "\tZ/Rin/Rout "
0303                                    << cms::convert2mm(zminBlock[i]) << ", " << cms::convert2mm(zmaxBlock[i]) << "/ "
0304                                    << cms::convert2mm(rinBlock1[i]) << ", " << cms::convert2mm(rinBlock2[i]) << "/ "
0305                                    << cms::convert2mm(routBlock1[i]) << ", " << cms::convert2mm(routBlock2[i]);
0306 #endif
0307 
0308     idName = args.value<std::string>("MotherName");
0309     idOffset = args.value<int>("IdOffset");
0310 #ifdef EDM_ML_DEBUG
0311     edm::LogVerbatim("HCalGeom") << "DDHCalEndcapAlgo: Parent " << args.parentName() << " idName " << idName
0312                                  << " NameSpace " << ns.name() << " Offset " << idOffset;
0313 #endif
0314 
0315     tolPos = args.value<double>("TolPos");
0316     tolAbs = args.value<double>("TolAbs");
0317 #ifdef EDM_ML_DEBUG
0318     edm::LogVerbatim("HCalGeom") << "DDHCalEndcapAlgo: Tolerances - Positioning " << cms::convert2mm(tolPos)
0319                                  << " Absorber " << cms::convert2mm(tolAbs);
0320     edm::LogVerbatim("HCalGeom") << "==>> Constructing DDHCalEndcapAlgo...";
0321 #endif
0322 
0323     dd4hep::Volume parent = ns.volume(args.parentName());
0324     constructGeneralVolume(ns, parent);
0325 #ifdef EDM_ML_DEBUG
0326     edm::LogVerbatim("HCalGeom") << "<<== End of DDHCalEndcapAlgo construction ...";
0327 #endif
0328   }
0329 
0330   void constructGeneralVolume(cms::DDNamespace& ns, dd4hep::Volume& parent) {
0331 #ifdef EDM_ML_DEBUG
0332     edm::LogVerbatim("HCalGeom") << "DDHCalEndcapAlgo: General volume...";
0333 #endif
0334 
0335     bool proto = true;
0336     for (int i = 0; i < 3; i++)
0337       if (eModule[i] > 0)
0338         proto = false;
0339 
0340     dd4hep::Rotation3D rot = getRotation(rotation, ns);
0341 #ifdef EDM_ML_DEBUG
0342     edm::LogVerbatim("HCalGeom") << " Rotation matrix " << rotation << " Rotation " << rot;
0343 #endif
0344 
0345     dd4hep::Position r0(0, 0, zShift);
0346     double alpha = (1._pi) / nsectors;
0347     double dphi = nsectortot * (2._pi) / nsectors;
0348 
0349     //!!!!!!!!!!!!!!!!!Should be zero. And removed as soon as
0350     //vertical walls are allowed in SolidPolyhedra
0351     double delz = 0;
0352 
0353     std::vector<double> pgonZ, pgonRmin, pgonRmax;
0354     if (proto) {
0355       double zf = ziBody + zShiftHac2;
0356       pgonZ.emplace_back(zf - dzShift);
0357       pgonRmin.emplace_back(zf * tan(angBot));
0358       pgonRmax.emplace_back((zf - z1Beam) * slope);
0359     } else {
0360       pgonZ.emplace_back(zFront - dzShift);
0361       pgonRmin.emplace_back(zFront * tan(angTop));
0362       pgonRmax.emplace_back((zFront - z1Beam) * slope);
0363       pgonZ.emplace_back(ziL0Body - dzShift);
0364       pgonRmin.emplace_back(ziL0Body * tan(angTop));
0365       pgonRmax.emplace_back((ziL0Body - z1Beam) * slope);
0366       pgonZ.emplace_back(ziL0Body - dzShift);
0367       pgonRmin.emplace_back(ziL0Body * tan(angBot));
0368       pgonRmax.emplace_back((ziL0Body - z1Beam) * slope);
0369     }
0370     pgonZ.emplace_back(ziKink - dzShift);
0371     pgonRmin.emplace_back(riKink);
0372     pgonRmax.emplace_back(rout);
0373     pgonZ.emplace_back(ziDip - dzShift);
0374     pgonRmin.emplace_back(riDip);
0375     pgonRmax.emplace_back(rout);
0376     pgonZ.emplace_back(ziDip - dzShift + delz);
0377     pgonRmin.emplace_back(riDip);
0378     pgonRmax.emplace_back(roDip);
0379     pgonZ.emplace_back(zEnd - dzShift);
0380     pgonRmin.emplace_back(zEnd * tan(angBot));
0381     pgonRmax.emplace_back(roDip);
0382     pgonZ.emplace_back(zEnd);
0383     pgonRmin.emplace_back(zEnd * tan(angBot));
0384     pgonRmax.emplace_back(roDip);
0385 
0386     std::string name("Null");
0387     dd4hep::Solid solid = dd4hep::Polyhedra(ns.prepend(idName), nsectortot, -alpha, dphi, pgonZ, pgonRmin, pgonRmax);
0388 #ifdef EDM_ML_DEBUG
0389     edm::LogVerbatim("HCalGeom") << "DDHCalEndcapAlgo: " << solid.name() << " Polyhedra made of " << genMaterial
0390                                  << " with " << nsectortot << " sectors from " << convertRadToDeg(-alpha) << " to "
0391                                  << convertRadToDeg(-alpha + dphi) << " and with " << pgonZ.size() << " sections";
0392     for (unsigned int i = 0; i < pgonZ.size(); i++)
0393       edm::LogVerbatim("HCalGeom") << "\t\tZ = " << cms::convert2mm(pgonZ[i])
0394                                    << "\tRmin = " << cms::convert2mm(pgonRmin[i])
0395                                    << "\tRmax = " << cms::convert2mm(pgonRmax[i]);
0396 #endif
0397     dd4hep::Material matter = ns.material(genMaterial);
0398     dd4hep::Volume genlogic(solid.name(), solid, matter);
0399 
0400     parent.placeVolume(genlogic, 1, dd4hep::Transform3D(rot, r0));
0401 #ifdef EDM_ML_DEBUG
0402     edm::LogVerbatim("HCalGeom") << "DDHCalEndcapAlgo: " << genlogic.name() << " number 1 positioned in "
0403                                  << parent.name() << " at (0, 0, " << cms::convert2mm(zShift)
0404                                  << ") with rotation: " << rot;
0405 #endif
0406 
0407     if (nEndcap != 1) {
0408       rot = getRotation(rotHalf, ns);
0409       parent.placeVolume(genlogic, 2, dd4hep::Transform3D(rot, r0));
0410 #ifdef EDM_ML_DEBUG
0411       edm::LogVerbatim("HCalGeom") << "DDHCalEndcapAlgo: " << genlogic.name() << " number 2 "
0412                                    << "positioned in " << parent.name() << " at (0, 0, " << cms::convert2mm(zShift)
0413                                    << ") with rotation: " << rot;
0414 #endif
0415     }
0416 
0417     //Forward half
0418     name = idName + "Front";
0419     std::vector<double> pgonZMod, pgonRminMod, pgonRmaxMod;
0420     for (unsigned int i = 0; i < (pgonZ.size() - 1); i++) {
0421       pgonZMod.emplace_back(pgonZ[i] + dzShift);
0422       pgonRminMod.emplace_back(pgonRmin[i]);
0423       pgonRmaxMod.emplace_back(pgonRmax[i]);
0424     }
0425     solid = dd4hep::Polyhedra(ns.prepend(name), nsectortot, -alpha, dphi, pgonZMod, pgonRminMod, pgonRmaxMod);
0426 #ifdef EDM_ML_DEBUG
0427     edm::LogVerbatim("HCalGeom") << "DDHCalEndcapAlgo: " << solid.name() << " Polyhedra made of " << genMaterial
0428                                  << " with " << nsectortot << " sectors from " << convertRadToDeg(-alpha) << " to "
0429                                  << convertRadToDeg(-alpha + dphi) << " and with " << pgonZMod.size() << " sections ";
0430     for (unsigned int i = 0; i < pgonZMod.size(); i++)
0431       edm::LogVerbatim("HCalGeom") << "\t\tZ = " << cms::convert2mm(pgonZMod[i])
0432                                    << "\tRmin = " << cms::convert2mm(pgonRminMod[i])
0433                                    << "\tRmax = " << cms::convert2mm(pgonRmaxMod[i]);
0434 #endif
0435 
0436     dd4hep::Volume genlogich(solid.name(), solid, matter);
0437     ns.addVolumeNS(genlogich);
0438     genlogic.placeVolume(genlogich, 1, dd4hep::Position(0, 0, -dzShift));
0439 #ifdef EDM_ML_DEBUG
0440     edm::LogVerbatim("HCalGeom") << "DDHCalEndcapAlgo: " << genlogich.name() << " number 1 positioned in "
0441                                  << genlogic.name() << " at (0,0," << -cms::convert2mm(dzShift) << ") with no rotation";
0442 #endif
0443 
0444     //Construct sector (from -alpha to +alpha)
0445     name = idName + "Module";
0446     solid = dd4hep::Polyhedra(ns.prepend(name), 1, -alpha, 2 * alpha, pgonZMod, pgonRminMod, pgonRmaxMod);
0447 #ifdef EDM_ML_DEBUG
0448     edm::LogVerbatim("HCalGeom") << "DDHCalEndcapAlgo: " << solid.name() << " Polyhedra made of " << genMaterial
0449                                  << " with 1 sector from " << convertRadToDeg(-alpha) << " to "
0450                                  << convertRadToDeg(alpha) << " and with " << pgonZMod.size() << " sections";
0451     for (unsigned int i = 0; i < pgonZMod.size(); i++)
0452       edm::LogVerbatim("HCalGeom") << "\t\tZ = " << cms::convert2mm(pgonZMod[i])
0453                                    << "\tRmin = " << cms::convert2mm(pgonRminMod[i])
0454                                    << "\tRmax = " << cms::convert2mm(pgonRmaxMod[i]);
0455 #endif
0456 
0457     dd4hep::Volume seclogic(solid.name(), solid, matter);
0458     for (int ii = 0; ii < nsectortot; ii++) {
0459       double phi = ii * 2 * alpha;
0460       dd4hep::Rotation3D rot0;
0461       if (phi != 0) {
0462         rot0 = dd4hep::RotationZ(phi);
0463 #ifdef EDM_ML_DEBUG
0464         edm::LogVerbatim("HCalGeom") << "DDHCalEndcapAlgo: Creating a new rotation \t 90," << convertRadToDeg(phi)
0465                                      << ", 90," << convertRadToDeg(phi + 90._deg) << ", 0, 0";
0466 #endif
0467       }
0468       genlogich.placeVolume(seclogic, ii + 1, rot0);
0469 #ifdef EDM_ML_DEBUG
0470       edm::LogVerbatim("HCalGeom") << "DDHCalEndcapAlgo: " << seclogic.name() << " number " << (ii + 1)
0471                                    << " positioned in " << genlogich.name() << " at (0, 0, 0) with rotation: " << rot0;
0472 #endif
0473     }
0474 
0475     //Construct the things inside the sector
0476     constructInsideSector(ns, seclogic);
0477 
0478     //Backward half
0479     name = idName + "Back";
0480     std::vector<double> pgonZBack, pgonRminBack, pgonRmaxBack;
0481     pgonZBack.emplace_back(zEnd - dzShift);
0482     pgonRminBack.emplace_back(pgonZBack[0] * tan(angBot) + drEnd);
0483     pgonRmaxBack.emplace_back(roDip);
0484     pgonZBack.emplace_back(zEnd);
0485     pgonRminBack.emplace_back(pgonZBack[1] * tan(angBot) + drEnd);
0486     pgonRmaxBack.emplace_back(roDip);
0487     solid = dd4hep::Polyhedra(ns.prepend(name), nsectortot, -alpha, dphi, pgonZBack, pgonRminBack, pgonRmaxBack);
0488 #ifdef EDM_ML_DEBUG
0489     edm::LogVerbatim("HCalGeom") << "DDHCalEndcapAlgo: " << solid.name() << " Polyhedra made of " << absMat << " with "
0490                                  << nsectortot << " sectors from " << convertRadToDeg(-alpha) << " to "
0491                                  << convertRadToDeg(-alpha + dphi) << " and with " << pgonZBack.size() << " sections";
0492     for (unsigned int i = 0; i < pgonZBack.size(); i++)
0493       edm::LogVerbatim("HCalGeom") << "\t\tZ = " << cms::convert2mm(pgonZBack[i])
0494                                    << "\tRmin = " << cms::convert2mm(pgonRminBack[i])
0495                                    << "\tRmax = " << cms::convert2mm(pgonRmaxBack[i]);
0496 #endif
0497 
0498     dd4hep::Material absMatter = ns.material(absMat);
0499     dd4hep::Volume glog(solid.name(), solid, absMatter);
0500     genlogic.placeVolume(glog, 1);
0501 #ifdef EDM_ML_DEBUG
0502     edm::LogVerbatim("HCalGeom") << "DDHCalEndcapAlgo: " << glog.name() << " number 1 positioned in " << genlogic.name()
0503                                  << " at (0,0,0) with no rotation";
0504 #endif
0505   }
0506 
0507   void constructInsideSector(cms::DDNamespace& ns, dd4hep::Volume& sector) {
0508 #ifdef EDM_ML_DEBUG
0509     edm::LogVerbatim("HCalGeom") << "DDHCalEndcapAlgo: Modules (" << modules << ") ...";
0510 #endif
0511 
0512     double alpha = (1._pi) / nsectors;
0513     for (int i = 0; i < modules; i++) {
0514       std::string name = idName + modName[i];
0515       dd4hep::Material matter = ns.material(modMat[i]);
0516 
0517       if (eModule[i] > 0) {
0518         int nsec = sectionModule[i];
0519 
0520         //!!!!!!!!!!!!!!!!!Should be zero. And removed as soon as
0521         //vertical walls are allowed in SolidPolyhedra
0522         double deltaz = 0;
0523 
0524         std::vector<double> pgonZ, pgonRmin, pgonRmax;
0525         if (nsec == 3) {
0526           double zf = zminBlock[i] + zShiftHac2;
0527           pgonZ.emplace_back(zf);
0528           pgonRmin.emplace_back(zf * tan(angBot));
0529           pgonRmax.emplace_back((zf - z1Beam) * slope);
0530           pgonZ.emplace_back(ziKink);
0531           pgonRmin.emplace_back(riKink);
0532           pgonRmax.emplace_back(rout);
0533         } else {
0534           pgonZ.emplace_back(zminBlock[i]);
0535           pgonRmin.emplace_back(rinBlock1[i]);
0536           pgonRmax.emplace_back(routBlock1[i]);
0537         }
0538         if (nsec == 4) {
0539           pgonZ.emplace_back(ziDip);
0540           pgonRmin.emplace_back(riDip);
0541           pgonRmax.emplace_back(rout);
0542           pgonZ.emplace_back(pgonZ[1] + deltaz);
0543           pgonRmin.emplace_back(pgonRmin[1]);
0544           pgonRmax.emplace_back(roDip);
0545         }
0546         pgonZ.emplace_back(zmaxBlock[i]);
0547         pgonRmin.emplace_back(rinBlock2[i]);
0548         pgonRmax.emplace_back(routBlock2[i]);
0549 
0550         //Solid & volume
0551         dd4hep::Solid solid = dd4hep::Polyhedra(ns.prepend(name), 1, -alpha, 2 * alpha, pgonZ, pgonRmin, pgonRmax);
0552 #ifdef EDM_ML_DEBUG
0553         edm::LogVerbatim("HCalGeom") << "DDHCalEndcapAlgo: " << solid.name() << " Polyhedra made of " << modMat[i]
0554                                      << " with 1 sector from " << convertRadToDeg(-alpha) << " to "
0555                                      << convertRadToDeg(alpha) << " and with " << nsec << " sections";
0556         for (unsigned int k = 0; k < pgonZ.size(); k++)
0557           edm::LogVerbatim("HCalGeom") << "\t\tZ = " << cms::convert2mm(pgonZ[k])
0558                                        << "\tRmin = " << cms::convert2mm(pgonRmin[k])
0559                                        << "\tRmax = " << cms::convert2mm(pgonRmax[k]);
0560 #endif
0561 
0562         dd4hep::Volume glog(solid.name(), solid, matter);
0563 
0564         sector.placeVolume(glog, i + 1);
0565 #ifdef EDM_ML_DEBUG
0566         edm::LogVerbatim("HCalGeom") << "DDHCalEndcapAlgo: " << glog.name() << " number " << (i + 1)
0567                                      << " positioned in " << sector.name() << " at (0,0,0) with no rotation";
0568 #endif
0569 
0570         if (modType[i] == 0)
0571           constructInsideModule0(ns, glog, i);
0572         else
0573           constructInsideModule(ns, glog, i);
0574       }
0575     }
0576   }
0577 
0578   void constructInsideModule0(cms::DDNamespace& ns, dd4hep::Volume& module, int mod) {
0579 #ifdef EDM_ML_DEBUG
0580     edm::LogVerbatim("HCalGeom") << "DDHCalEndcapAlgo: \t\tInside module0 ..." << mod;
0581 #endif
0582 
0583     ///////////////////////////////////////////////////////////////
0584     //Pointers to the Rotation Matrices and to the Materials
0585     dd4hep::Rotation3D rot = getRotation(rotmat, ns);
0586     dd4hep::Material matabsorbr = ns.material(absMat);
0587     dd4hep::Material matplastic = ns.material(plastMat);
0588 
0589     int layer = getLayer(mod, 0);
0590     int layer0 = getLayer(mod, 1);
0591     std::string name;
0592     double xpos, ypos, zpos;
0593     dd4hep::Solid solid;
0594     dd4hep::Volume glog;
0595     for (int iphi = 0; iphi < phiSections; iphi++) {
0596       double yh, bl, tl, alp;
0597       parameterLayer0(mod, layer, iphi, yh, bl, tl, alp, xpos, ypos, zpos);
0598       name = DDSplit(module.name()).first + layerName[layer] + phiName[iphi];
0599       solid = dd4hep::Trap(ns.prepend(name), 0.5 * layerT[layer], 0, 0, yh, bl, tl, alp, yh, bl, tl, alp);
0600 #ifdef EDM_ML_DEBUG
0601       edm::LogVerbatim("HCalGeom") << "DDHCalEndcapAlgo: " << solid.name() << " Trap made of " << plastMat
0602                                    << " of dimensions " << cms::convert2mm(0.5 * layerT[layer]) << ", 0, 0, "
0603                                    << cms::convert2mm(yh) << ", " << cms::convert2mm(bl) << ", " << cms::convert2mm(tl)
0604                                    << ", " << convertRadToDeg(alp) << ", " << cms::convert2mm(yh) << ", "
0605                                    << cms::convert2mm(bl) << ", " << cms::convert2mm(tl) << ", "
0606                                    << convertRadToDeg(alp);
0607 #endif
0608 
0609       glog = dd4hep::Volume(solid.name(), solid, matplastic);
0610 
0611       dd4hep::Position r1(xpos, ypos, zpos);
0612       module.placeVolume(glog, idOffset + layer + 1, dd4hep::Transform3D(rot, r1));
0613 
0614 #ifdef EDM_ML_DEBUG
0615       edm::LogVerbatim("HCalGeom") << "DDHCalEndcapAlgo: " << glog.name() << " number " << (idOffset + layer + 1)
0616                                    << " positioned in " << module.name() << " at (" << cms::convert2mm(xpos) << ", "
0617                                    << cms::convert2mm(ypos) << ", " << cms::convert2mm(zpos)
0618                                    << " with rotation: " << rot;
0619 #endif
0620 
0621       //Now construct the layer of scintillator inside this
0622       int copyNo = layer0 * 10 + layerType[layer];
0623       name = modName[mod] + layerName[layer] + phiName[iphi];
0624       constructScintLayer(ns, glog, scintT[layer], yh, bl, tl, alp, name, copyNo);
0625     }
0626 
0627     //Now the absorber layer
0628     double zi = zminBlock[mod] + layerT[layer];
0629     double zo = zi + 0.5 * dzStep;
0630     double rinF, routF, rinB, routB;
0631     if (mod == 0) {
0632       rinF = zi * tan(angTop);
0633       routF = (zi - z1Beam) * slope;
0634       rinB = zo * tan(angTop);
0635       routB = (zo - z1Beam) * slope;
0636     } else {
0637       rinF = zi * tan(angBot);
0638       routF = zi * tan(angTop);
0639       rinB = zo * tan(angBot);
0640       routB = zo * tan(angTop);
0641     }
0642 
0643 #ifdef EDM_ML_DEBUG
0644     edm::LogVerbatim("HCalGeom") << "DDHCalEndcapAlgo: Module " << mod << " Front " << cms::convert2mm(zi) << ", "
0645                                  << cms::convert2mm(rinF) << ", " << cms::convert2mm(routF) << " Back "
0646                                  << cms::convert2mm(zo) << ", " << cms::convert2mm(rinB) << ", "
0647                                  << cms::convert2mm(routB);
0648 #endif
0649 
0650     double yh1, bl1, tl1, yh2, bl2, tl2, theta, phi, alp;
0651     parameterLayer(
0652         0, rinF, routF, rinB, routB, zi, zo, yh1, bl1, tl1, yh2, bl2, tl2, alp, theta, phi, xpos, ypos, zpos);
0653     double fact = tolAbs;
0654 
0655 #ifdef EDM_ML_DEBUG
0656     edm::LogVerbatim("HCalGeom") << "DDHCalEndcapAlgo: Trim " << cms::convert2mm(fact) << " Param "
0657                                  << cms::convert2mm(yh1) << ", " << cms::convert2mm(bl1) << ", " << cms::convert2mm(tl1)
0658                                  << ", " << cms::convert2mm(yh2) << ", " << cms::convert2mm(bl2) << ", "
0659                                  << cms::convert2mm(tl2);
0660 #endif
0661 
0662     bl1 -= fact;
0663     tl1 -= fact;
0664     bl2 -= fact;
0665     tl2 -= fact;
0666 
0667     name = DDSplit(module.name()).first + "Absorber";
0668     solid = dd4hep::Trap(ns.prepend(name), 0.5 * thick[mod], theta, phi, yh1, bl1, tl1, alp, yh2, bl2, tl2, alp);
0669 #ifdef EDM_ML_DEBUG
0670     edm::LogVerbatim("HCalGeom") << "DDHCalEndcapAlgo: " << solid.name() << " Trap made of " << absMat
0671                                  << " of dimensions " << cms::convert2mm(0.5 * thick[mod]) << ", "
0672                                  << convertRadToDeg(theta) << ", " << convertRadToDeg(phi) << ", "
0673                                  << cms::convert2mm(yh1) << ", " << cms::convert2mm(bl1) << ", " << cms::convert2mm(tl1)
0674                                  << ", " << convertRadToDeg(alp) << ", " << cms::convert2mm(yh2) << ", "
0675                                  << cms::convert2mm(bl2) << ", " << cms::convert2mm(tl2) << ", "
0676                                  << convertRadToDeg(alp);
0677 #endif
0678 
0679     glog = dd4hep::Volume(solid.name(), solid, matabsorbr);
0680 
0681     dd4hep::Position r2(xpos, ypos, zpos);
0682     module.placeVolume(glog, 1, dd4hep::Transform3D(rot, r2));
0683 
0684 #ifdef EDM_ML_DEBUG
0685     edm::LogVerbatim("HCalGeom") << "DDHCalEndcapAlgo: " << glog.name() << " number 1 positioned in " << module.name()
0686                                  << " at (" << cms::convert2mm(xpos) << ", " << cms::convert2mm(ypos) << ", "
0687                                  << cms::convert2mm(zpos) << ") with rotation: " << rot;
0688 #endif
0689   }
0690 
0691   void constructInsideModule(cms::DDNamespace& ns, dd4hep::Volume& module, int mod) {
0692 #ifdef EDM_ML_DEBUG
0693     edm::LogVerbatim("HCalGeom") << "DDHCalEndcapAlgo: \t\tInside module ..." << mod;
0694 #endif
0695 
0696     ///////////////////////////////////////////////////////////////
0697     //Pointers to the Rotation Matrices and to the Materials
0698     dd4hep::Rotation3D rot = getRotation(rotmat, ns);
0699     dd4hep::Material matter = ns.material(genMaterial);
0700     dd4hep::Material matplastic = ns.material(plastMat);
0701 
0702     double alpha = (1._pi) / nsectors;
0703     double zi = zminBlock[mod];
0704 
0705     for (int i = 0; i < layerN[mod]; i++) {
0706       std::string name;
0707       dd4hep::Solid solid;
0708       dd4hep::Volume glog, plog;
0709       int layer = getLayer(mod, i);
0710       double zo = zi + 0.5 * dzStep;
0711 
0712       for (int iphi = 0; iphi < phiSections; iphi++) {
0713         double ziAir = zo - thick[mod];
0714         double rinF, rinB;
0715         if (layer == 1) {
0716           rinF = ziAir * tan(angTop);
0717           rinB = zo * tan(angTop);
0718         } else {
0719           rinF = ziAir * tan(angBot);
0720           rinB = zo * tan(angBot);
0721         }
0722         double routF = (ziAir - z1Beam) * slope;
0723         double routB = (zo - z1Beam) * slope;
0724         if (routF > routBlock2[mod])
0725           routF = routBlock2[mod];
0726         if (routB > routBlock2[mod])
0727           routB = routBlock2[mod];
0728 
0729 #ifdef EDM_ML_DEBUG
0730         edm::LogVerbatim("HCalGeom") << "DDHCalEndcapAlgo: Layer " << i << " Phi " << iphi << " Front "
0731                                      << cms::convert2mm(ziAir) << ", " << cms::convert2mm(rinF) << ", "
0732                                      << cms::convert2mm(routF) << " Back " << cms::convert2mm(zo) << ", "
0733                                      << cms::convert2mm(rinB) << ", " << cms::convert2mm(routB);
0734 #endif
0735 
0736         double yh1, bl1, tl1, yh2, bl2, tl2, theta, phi, alp;
0737         double xpos, ypos, zpos;
0738         parameterLayer(
0739             iphi, rinF, routF, rinB, routB, ziAir, zo, yh1, bl1, tl1, yh2, bl2, tl2, alp, theta, phi, xpos, ypos, zpos);
0740 
0741         name = DDSplit(module.name()).first + layerName[layer] + phiName[iphi] + "Air";
0742         solid = dd4hep::Trap(ns.prepend(name), 0.5 * thick[mod], theta, phi, yh1, bl1, tl1, alp, yh2, bl2, tl2, alp);
0743 #ifdef EDM_ML_DEBUG
0744         edm::LogVerbatim("HCalGeom") << "DDHCalEndcapAlgo: " << solid.name() << " Trap made of " << matter.name()
0745                                      << " of dimensions " << cms::convert2mm(0.5 * thick[mod]) << ", "
0746                                      << convertRadToDeg(theta) << ", " << convertRadToDeg(phi) << ", "
0747                                      << cms::convert2mm(yh1) << ", " << cms::convert2mm(bl1) << ", "
0748                                      << cms::convert2mm(tl1) << ", " << convertRadToDeg(alp) << ", "
0749                                      << cms::convert2mm(yh2) << ", " << cms::convert2mm(bl2) << ", "
0750                                      << cms::convert2mm(tl2) << ", " << convertRadToDeg(alp);
0751 #endif
0752 
0753         glog = dd4hep::Volume(solid.name(), solid, matter);
0754         dd4hep::Position r1(xpos, ypos, zpos);
0755         module.placeVolume(glog, layer + 1, dd4hep::Transform3D(rot, r1));
0756 
0757 #ifdef EDM_ML_DEBUG
0758         edm::LogVerbatim("HCalGeom") << "DDHCalEndcapAlgo: " << glog.name() << " number " << (layer + 1)
0759                                      << " positioned in " << module.name() << " at (" << cms::convert2mm(xpos) << ", "
0760                                      << cms::convert2mm(ypos) << ", " << cms::convert2mm(zpos)
0761                                      << ") with rotation: " << rot;
0762 #endif
0763 
0764         //Now the plastic with scintillators
0765         double yh = 0.5 * (routF - rinB) - getTrim(mod, iphi);
0766         double bl = 0.5 * rinB * tan(alpha) - getTrim(mod, iphi);
0767         double tl = 0.5 * routF * tan(alpha) - getTrim(mod, iphi);
0768         name = DDSplit(module.name()).first + layerName[layer] + phiName[iphi];
0769         solid = dd4hep::Trap(ns.prepend(name), 0.5 * layerT[layer], 0, 0, yh, bl, tl, alp, yh, bl, tl, alp);
0770 #ifdef EDM_ML_DEBUG
0771         edm::LogVerbatim("HCalGeom") << "DDHCalEndcapAlgo: " << solid.name() << " Trap made of " << plastMat
0772                                      << " of dimensions " << cms::convert2mm(0.5 * layerT[layer]) << ", 0, 0, "
0773                                      << cms::convert2mm(yh) << ", " << cms::convert2mm(bl) << ", "
0774                                      << cms::convert2mm(tl) << ", " << convertRadToDeg(alp) << ", "
0775                                      << cms::convert2mm(yh) << ", " << cms::convert2mm(bl) << ", "
0776                                      << cms::convert2mm(tl) << ", " << convertRadToDeg(alp);
0777 #endif
0778 
0779         plog = dd4hep::Volume(solid.name(), solid, matplastic);
0780         ypos = 0.5 * (routF + rinB) - xpos;
0781         glog.placeVolume(plog, idOffset + layer + 1, dd4hep::Position(0., ypos, 0.));
0782 
0783 #ifdef EDM_ML_DEBUG
0784         edm::LogVerbatim("HCalGeom") << "DDHCalEndcapAlgo: " << plog.name() << " number " << (idOffset + layer + 1)
0785                                      << " positioned in " << glog.name() << " at (0, " << cms::convert2mm(ypos)
0786                                      << ", 0) with no rotation";
0787 #endif
0788 
0789         //Constructing the scintillators inside
0790         int copyNo = layer * 10 + layerType[layer];
0791         name = modName[mod] + layerName[layer] + phiName[iphi];
0792         constructScintLayer(ns, plog, scintT[layer], yh, bl, tl, alp, name, copyNo);
0793         zo += 0.5 * dzStep;
0794       }  // End of loop over phi indices
0795       zi = zo - 0.5 * dzStep;
0796     }  // End of loop on layers
0797   }
0798 
0799   void constructScintLayer(cms::DDNamespace& ns,
0800                            dd4hep::Volume& detector,
0801                            double dz,
0802                            double yh,
0803                            double bl,
0804                            double tl,
0805                            double alp,
0806                            const std::string& nm,
0807                            int id) {
0808     dd4hep::Material matter = ns.material(scintMat);
0809     std::string name = idName + "Scintillator" + nm;
0810 
0811     dd4hep::Solid solid = dd4hep::Trap(ns.prepend(name), 0.5 * dz, 0, 0, yh, bl, tl, alp, yh, bl, tl, alp);
0812 #ifdef EDM_ML_DEBUG
0813     edm::LogVerbatim("HCalGeom") << "DDHCalEndcapAlgo: " << solid.name() << " Trap made of " << scintMat
0814                                  << " of dimensions " << cms::convert2mm(0.5 * dz) << ", 0, 0, " << cms::convert2mm(yh)
0815                                  << ", " << cms::convert2mm(bl) << ", " << cms::convert2mm(tl) << ", "
0816                                  << convertRadToDeg(alp) << ", " << cms::convert2mm(yh) << ", " << cms::convert2mm(bl)
0817                                  << ", " << cms::convert2mm(tl) << ", " << convertRadToDeg(alp);
0818 #endif
0819 
0820     dd4hep::Volume glog(solid.name(), solid, matter);
0821     detector.placeVolume(glog, id);
0822 #ifdef EDM_ML_DEBUG
0823     edm::LogVerbatim("HCalGeom") << "DDHCalEndcapAlgo: " << glog.name() << " number " << id << " positioned in "
0824                                  << detector.name() << " at (0,0,0) with no rotation";
0825 #endif
0826   }
0827 
0828   int getLayer(unsigned int i, unsigned int j) const {
0829     switch (i) {
0830       case 0:
0831         return layerN0[j];
0832         break;
0833       case 1:
0834         return layerN1[j];
0835         break;
0836       case 2:
0837         return layerN2[j];
0838         break;
0839       case 3:
0840         return layerN3[j];
0841         break;
0842       case 4:
0843         return layerN4[j];
0844         break;
0845       case 5:
0846         return layerN5[j];
0847         break;
0848       default:
0849         return 0;
0850     }
0851   }
0852 
0853   double getTrim(unsigned int i, unsigned int j) const {
0854     if (j == 0)
0855       return trimLeft[i];
0856     else
0857       return trimRight[j];
0858   }
0859 
0860   void parameterLayer0(int mod,
0861                        int layer,
0862                        int iphi,
0863                        double& yh,
0864                        double& bl,
0865                        double& tl,
0866                        double& alp,
0867                        double& xpos,
0868                        double& ypos,
0869                        double& zpos) {
0870     //Given module and layer number compute parameters of trapezoid
0871     //and positioning parameters
0872     double alpha = (1._pi) / nsectors;
0873 #ifdef EDM_ML_DEBUG
0874     edm::LogVerbatim("HCalGeom") << "Input " << iphi << " " << layer << " " << iphi << " Alpha "
0875                                  << convertRadToDeg(alpha);
0876 #endif
0877 
0878     double zi, zo;
0879     if (iphi == 0) {
0880       zi = zminBlock[mod];
0881       zo = zi + layerT[layer];
0882     } else {
0883       zo = zmaxBlock[mod];
0884       zi = zo - layerT[layer];
0885     }
0886     double rin, rout;
0887     if (mod == 0) {
0888       rin = zo * tan(angTop);
0889       rout = (zi - z1Beam) * slope;
0890     } else {
0891       rin = zo * tan(angBot);
0892       rout = zi * tan(angTop);
0893     }
0894     yh = 0.5 * (rout - rin);
0895     bl = 0.5 * rin * tan(alpha);
0896     tl = 0.5 * rout * tan(alpha);
0897     xpos = 0.5 * (rin + rout);
0898     ypos = 0.5 * (bl + tl);
0899     zpos = 0.5 * (zi + zo);
0900     yh -= getTrim(mod, iphi);
0901     bl -= getTrim(mod, iphi);
0902     tl -= getTrim(mod, iphi);
0903     alp = atan(0.5 * tan(alpha));
0904     if (iphi == 0) {
0905       ypos = -ypos;
0906     } else {
0907       alp = -alp;
0908     }
0909 
0910 #ifdef EDM_ML_DEBUG
0911     edm::LogVerbatim("HCalGeom") << "Output Dimensions " << cms::convert2mm(yh) << " " << cms::convert2mm(bl) << " "
0912                                  << cms::convert2mm(tl) << " " << convertRadToDeg(alp) << " Position "
0913                                  << cms::convert2mm(xpos) << " " << cms::convert2mm(ypos) << " "
0914                                  << cms::convert2mm(zpos);
0915 #endif
0916   }
0917 
0918   void parameterLayer(int iphi,
0919                       double rinF,
0920                       double routF,
0921                       double rinB,
0922                       double routB,
0923                       double zi,
0924                       double zo,
0925                       double& yh1,
0926                       double& bl1,
0927                       double& tl1,
0928                       double& yh2,
0929                       double& bl2,
0930                       double& tl2,
0931                       double& alp,
0932                       double& theta,
0933                       double& phi,
0934                       double& xpos,
0935                       double& ypos,
0936                       double& zpos) {
0937     //Given rin, rout compute parameters of the trapezoid and
0938     //position of the trapezoid for a standrd layer
0939     double alpha = (1._pi) / nsectors;
0940 
0941 #ifdef EDM_ML_DEBUG
0942     edm::LogVerbatim("HCalGeom") << "Input " << iphi << " Front " << cms::convert2mm(rinF) << " "
0943                                  << cms::convert2mm(routF) << " " << cms::convert2mm(zi) << " Back "
0944                                  << cms::convert2mm(rinB) << " " << cms::convert2mm(routB) << " " << cms::convert2mm(zo)
0945                                  << " Alpha " << convertRadToDeg(alpha);
0946 #endif
0947 
0948     yh1 = 0.5 * (routF - rinB);
0949     bl1 = 0.5 * rinB * tan(alpha);
0950     tl1 = 0.5 * routF * tan(alpha);
0951     yh2 = 0.5 * (routF - rinB);
0952     bl2 = 0.5 * rinB * tan(alpha);
0953     tl2 = 0.5 * routF * tan(alpha);
0954     double dx = 0.25 * (bl2 + tl2 - bl1 - tl1);
0955     double dy = 0.5 * (rinB + routF - rinB - routF);
0956     xpos = 0.25 * (rinB + routF + rinB + routF);
0957     ypos = 0.25 * (bl2 + tl2 + bl1 + tl1);
0958     zpos = 0.5 * (zi + zo);
0959     alp = atan(0.5 * tan(alpha));
0960     //  ypos-= tolPos;
0961     if (iphi == 0) {
0962       ypos = -ypos;
0963     } else {
0964       alp = -alp;
0965       dx = -dx;
0966     }
0967     double r = sqrt(dx * dx + dy * dy);
0968     theta = atan(r / (zo - zi));
0969     phi = atan2(dy, dx);
0970 
0971 #ifdef EDM_ML_DEBUG
0972     edm::LogVerbatim("HCalGeom") << "Output Dimensions " << cms::convert2mm(yh1) << " " << cms::convert2mm(bl1) << " "
0973                                  << cms::convert2mm(tl1) << " " << cms::convert2mm(yh2) << " " << cms::convert2mm(bl2)
0974                                  << " " << cms::convert2mm(tl2) << " " << convertRadToDeg(alp) << " "
0975                                  << convertRadToDeg(theta) << " " << convertRadToDeg(phi) << " Position "
0976                                  << cms::convert2mm(xpos) << " " << cms::convert2mm(ypos) << " "
0977                                  << cms::convert2mm(zpos);
0978 #endif
0979   }
0980 
0981   dd4hep::Rotation3D getRotation(const std::string& rotation, cms::DDNamespace& ns) {
0982     std::string rot = (strchr(rotation.c_str(), NAMESPACE_SEP) == nullptr) ? ("rotations:" + rotation) : rotation;
0983 #ifdef EDM_ML_DEBUG
0984     edm::LogVerbatim("HCalGeom") << "getRotation: " << rotation << ":" << rot << ":" << ns.rotation(rot);
0985 #endif
0986     return ns.rotation(rot);
0987   }
0988 };
0989 
0990 static long algorithm(dd4hep::Detector& /* description */, cms::DDParsingContext& ctxt, xml_h e) {
0991   HCalEndcapAlgo hcalendcapalgo(ctxt, e);
0992   return cms::s_executed;
0993 }
0994 
0995 // first argument is the type from the xml file
0996 DECLARE_DDCMS_DETELEMENT(DDCMS_hcal_DDHCalEndcapAlgo, algorithm)