Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:08:11

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