DDHGCalNoTaperEndcap

Line Code
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189
#include <algorithm>
#include <cmath>
#include <string>
#include <vector>

#include "DataFormats/Math/interface/angle_units.h"
#include "DetectorDescription/Core/interface/DDAlgorithm.h"
#include "DetectorDescription/Core/interface/DDAlgorithmFactory.h"
#include "DetectorDescription/Core/interface/DDCurrentNamespace.h"
#include "DetectorDescription/Core/interface/DDLogicalPart.h"
#include "DetectorDescription/Core/interface/DDMaterial.h"
#include "DetectorDescription/Core/interface/DDSolid.h"
#include "DetectorDescription/Core/interface/DDSplit.h"
#include "DetectorDescription/Core/interface/DDTypes.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/PluginManager/interface/PluginFactory.h"

//#define EDM_ML_DEBUG
using namespace angle_units::operators;

class DDHGCalNoTaperEndcap : public DDAlgorithm {
public:
  DDHGCalNoTaperEndcap(void);

  void initialize(const DDNumericArguments& nArgs,
                  const DDVectorArguments& vArgs,
                  const DDMapArguments& mArgs,
                  const DDStringArguments& sArgs,
                  const DDStringVectorArguments& vsArgs) override;

  void execute(DDCompactView& cpv) override;

private:
  int createQuarter(DDCompactView& cpv, int xQuadrant, int yQuadrant, int startCopyNo);

  double m_startAngle;        // Start angle
  double m_tiltAngle;         // Tilt  angle
  int m_invert;               // Inverted or forward
  double m_rMin;              // Inner radius
  double m_rMax;              // Outer radius
  double m_zoffset;           // Offset in z
  double m_xyoffset;          // Offset in x or y
  int m_n;                    // Mumber of copies
  int m_startCopyNo;          // Start copy Number
  int m_incrCopyNo;           // Increment copy Number
  std::string m_childName;    // Children name
  std::string m_idNameSpace;  // Namespace of this and ALL sub-parts
};

DDHGCalNoTaperEndcap::DDHGCalNoTaperEndcap() {
  edm::LogVerbatim("HGCalGeom") << "DDHGCalNoTaperEndcap test: Creating an instance";
}

void DDHGCalNoTaperEndcap::initialize(const DDNumericArguments& nArgs,
                                      const DDVectorArguments& vArgs,
                                      const DDMapArguments&,
                                      const DDStringArguments& sArgs,
                                      const DDStringVectorArguments&) {
  m_tiltAngle = nArgs["tiltAngle"];
  m_invert = int(nArgs["invert"]);
  m_rMin = int(nArgs["rMin"]);
  m_rMax = int(nArgs["rMax"]);
  m_zoffset = nArgs["zoffset"];
  m_xyoffset = nArgs["xyoffset"];
  m_n = int(nArgs["n"]);
  m_startCopyNo = int(nArgs["startCopyNo"]);
  m_incrCopyNo = int(nArgs["incrCopyNo"]);
  m_childName = sArgs["ChildName"];
#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HGCalGeom") << "Tilt Angle " << m_tiltAngle << " Invert " << m_invert << " R " << m_rMin << ":"
                                << m_rMax << " Offset " << m_zoffset << ":" << m_xyoffset << " N " << m_n << " Copy "
                                << m_startCopyNo << ":" << m_incrCopyNo << " Child " << m_childName;
#endif

  m_idNameSpace = DDCurrentNamespace::ns();
#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HGCalGeom") << "DDHGCalNoTaperEndcap: NameSpace " << m_idNameSpace << "\tParent "
                                << parent().name();
#endif
}

void DDHGCalNoTaperEndcap::execute(DDCompactView& cpv) {
  int lastCopyNo = m_startCopyNo;
#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HGCalGeom") << "Create quarter 1:1";
#endif
  lastCopyNo = createQuarter(cpv, 1, 1, lastCopyNo);
#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HGCalGeom") << "Create quarter -1:1";
#endif
  lastCopyNo = createQuarter(cpv, -1, 1, lastCopyNo);
#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HGCalGeom") << "Create quarter -1:-1";
#endif
  lastCopyNo = createQuarter(cpv, -1, -1, lastCopyNo);
#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HGCalGeom") << "Create quarter 1:-1";
#endif
  createQuarter(cpv, 1, -1, lastCopyNo);
}

int DDHGCalNoTaperEndcap::createQuarter(DDCompactView& cpv, int xQuadrant, int yQuadrant, int startCopyNo) {
  int copyNo = startCopyNo;
  double tiltAngle = m_tiltAngle;
  double xphi = xQuadrant * tiltAngle;
  double yphi = yQuadrant * tiltAngle;
  double theta = 90._deg;
  double phiX = 0.0;
  double phiY = theta;
  double phiZ = 3 * theta;
  double offsetZ = m_zoffset;
  double offsetXY = m_xyoffset;

  double offsetX = xQuadrant * 0.5 * offsetXY;
  double offsetY = yQuadrant * 0.5 * offsetXY;

#ifdef EDM_ML_DEBUG
  int rowmax(0), column(0);
#endif
  while (std::abs(offsetX) < m_rMax) {
#ifdef EDM_ML_DEBUG
    column++;
    int row(0);
#endif
    while (std::abs(offsetY) < m_rMax) {
#ifdef EDM_ML_DEBUG
      row++;
#endif
      double limit1 = sqrt((offsetX + 0.5 * xQuadrant * offsetXY) * (offsetX + 0.5 * xQuadrant * offsetXY) +
                           (offsetY + 0.5 * yQuadrant * offsetXY) * (offsetY + 0.5 * yQuadrant * offsetXY));
      double limit2 = sqrt((offsetX - 0.5 * xQuadrant * offsetXY) * (offsetX - 0.5 * xQuadrant * offsetXY) +
                           (offsetY - 0.5 * yQuadrant * offsetXY) * (offsetY - 0.5 * yQuadrant * offsetXY));
      // Make sure we do not add supermodules in rMin area
      if (limit2 > m_rMin && limit1 < m_rMax) {
#ifdef EDM_ML_DEBUG
        edm::LogVerbatim("HGCalGeom") << m_childName << " copyNo = " << copyNo << " (" << column << "," << row
                                      << "): offsetX,Y = " << offsetX << "," << offsetY << " limit=" << limit1 << ":"
                                      << limit2 << " rMin, rMax = " << m_rMin << "," << m_rMax;
#endif
        DDRotation rotation;
        std::string rotstr("NULL");

        // Check if we've already created the rotation matrix
        rotstr = "R";
        rotstr += std::to_string(copyNo);
        rotation = DDRotation(DDName(rotstr));
        if (!rotation) {
          rotation = DDrot(DDName(rotstr, m_idNameSpace),
                           std::make_unique<DDRotationMatrix>(
                               *DDcreateRotationMatrix(theta, phiX, theta + yphi, phiY, -yphi, phiZ) *
                               (*DDcreateRotationMatrix(theta + xphi, phiX, 90._deg, 90._deg, xphi, 0.0))));
        }

        DDTranslation tran(offsetX, offsetY, offsetZ);
#ifdef EDM_ML_DEBUG
        edm::LogVerbatim("HGCalGeom") << "Module " << copyNo << ": location = " << tran << " Rotation " << rotation;
#endif
        DDName parentName = parent().name();
        cpv.position(DDName(m_childName), parentName, copyNo, tran, rotation);

        copyNo += m_incrCopyNo;
      } else {
#ifdef EDM_ML_DEBUG
        edm::LogVerbatim("HGCalGeom") << " (" << column << "," << row << "): offsetX,Y = " << offsetX << "," << offsetY
                                      << " is out of limit=" << limit1 << ":" << limit2 << " rMin, rMax = " << m_rMin
                                      << "," << m_rMax;
#endif
      }

      yphi += yQuadrant * 2. * tiltAngle;
      offsetY += yQuadrant * offsetXY;
    }
#ifdef EDM_ML_DEBUG
    if (row > rowmax)
      rowmax = row;
#endif
    xphi += xQuadrant * 2. * tiltAngle;
    yphi = yQuadrant * tiltAngle;
    offsetY = yQuadrant * 0.5 * offsetXY;
    offsetX += xQuadrant * offsetXY;
  }
#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HGCalGeom") << rowmax << " rows and " << column << " columns in quadrant " << xQuadrant << ":"
                                << yQuadrant;
#endif
  return copyNo;
}

DEFINE_EDM_PLUGIN(DDAlgorithmFactory, DDHGCalNoTaperEndcap, "hgcal:DDHGCalNoTaperEndcap");