src/plugins/DDEcalEndcapAlgo.cpp

0001 #include "DD4hep/DetFactoryHelper.h"
0002 #include "DDCMS/DDCMSPlugins.h"
0003 #include <Math/AxisAngle.h>
0004 #include <Math/Rotation3D.h>
0005 #include <Math/Vector3D.h>
0006
0007 #include <CLHEP/Geometry/Transform3D.h>
0008 #include <CLHEP/Units/PhysicalConstants.h>
0009 #include <CLHEP/Units/SystemOfUnits.h>
0010 #include <CLHEP/Geometry/Point3D.h>
0011 #include <CLHEP/Geometry/Plane3D.h>
0012 #include <CLHEP/Geometry/Vector3D.h>
0013 #include <CLHEP/Geometry/Transform3D.h>
0014 #include <CLHEP/Vector/EulerAngles.h>
0015
0016 #include <array>
0017 #include <string>
0018 #include <vector>
0019
0020 using namespace std;
0021 using namespace dd4hep;
0022
0023 using DD3Vector = ROOT::Math::DisplacementVector3D<ROOT::Math::Cartesian3D<double>>;
0024 using DDTranslation = ROOT::Math::DisplacementVector3D<ROOT::Math::Cartesian3D<double> >;
0025 using DDRotation = ROOT::Math::Rotation3D;
0026 using DDRotationMatrix = ROOT::Math::Rotation3D;
0027 using DDAxisAngle = ROOT::Math::AxisAngle;
0028
0029 namespace {
0030
0031   constexpr long double piRadians(M_PI);
0032   constexpr long double degPerRad = 180. / piRadians;  // Degrees per radian
0033   constexpr double operator""_mm(long double length) { return length * 0.1; }
0034   constexpr long double operator""_deg(long double value) { return value / degPerRad; }
0035
0036   // Define Endcap Supercrystal class
0037   class DDEcalEndcapTrap {
0038   public:
0039     DDEcalEndcapTrap(const int hand, const double front, const double rear, const double length);
0040     DDEcalEndcapTrap() = delete;
0041
0042     void rotate(const DDRotationMatrix& rot);
0043     //void rotate(const DDTranslation& frontCentre, const DDTranslation& rearCentre);
0044     void translate(const DDTranslation& trans);
0045
0046     void rotateX(const double angle);
0047     void rotateY(const double angle);
0048     void translate();
0049     void moveto(const DDTranslation& frontCentre, const DDTranslation& rearCentre);
0050     double elevationAngle(const DDTranslation& trans);
0051     double polarAngle(const DDTranslation& trans);
0052     double elevationAngle();
0053     double polarAngle();
0054     DDTranslation cornerPos(const int icorner);
0055     void cornerPos(const int icorner, const DDTranslation& cc);
0056     DDTranslation centrePos();
0057     DDTranslation fcentrePos();
0058     DDTranslation rcentrePos();
0059     void calculateCorners();
0060     void calculateCentres();
0061     DDRotationMatrix rotation() { return m_rotation; }
0062     //void print();
0063
0064   private:
0065     DDRotationMatrix m_rotation;
0066     DDTranslation m_translation;
0067
0068     double m_centre[4];
0069     double m_fcentre[4];
0070     double m_rcentre[4];
0071     double m_corners[25];
0072     double m_front;
0073     double m_rear;
0074     double m_length;
0075
0076     int m_hand;
0077     int m_update;
0078   };
0079   // Implementation of DDEcalEndcapTrap class
0080
0081   DDEcalEndcapTrap::DDEcalEndcapTrap(const int hand, const double front, const double rear, const double length) {
0082     //
0083     //  Initialise corners of supercrystal.
0084
0085     // Start out with bottom surface on (x,z) plane, front face in (x,y) plane.
0086
0087     double xsign;
0088
0089     if (hand == 2) {
0090       xsign = -1.;
0091     } else {
0092       xsign = 1.;
0093     }
0094
0095     m_hand = hand;
0096     m_front = front;
0097     m_rear = rear;
0098     m_length = length;
0099
0100     int icorner;
0101     icorner = 1;
0102     m_corners[3 * icorner - 3] = xsign * front;
0103     m_corners[3 * icorner - 2] = front;
0104     m_corners[3 * icorner - 1] = 0.;
0105     icorner = 2;
0106     m_corners[3 * icorner - 3] = xsign * front;
0107     m_corners[3 * icorner - 2] = 0.;
0108     m_corners[3 * icorner - 1] = 0.;
0109     icorner = 3;
0110     m_corners[3 * icorner - 3] = 0.;
0111     m_corners[3 * icorner - 2] = 0.;
0112     m_corners[3 * icorner - 1] = 0.;
0113     icorner = 4;
0114     m_corners[3 * icorner - 3] = 0.;
0115     m_corners[3 * icorner - 2] = front;
0116     m_corners[3 * icorner - 1] = 0.;
0117
0118     icorner = 5;
0119     m_corners[3 * icorner - 3] = xsign * rear;
0120     m_corners[3 * icorner - 2] = rear;
0121     m_corners[3 * icorner - 1] = length;
0122     icorner = 6;
0123     m_corners[3 * icorner - 3] = xsign * rear;
0124     m_corners[3 * icorner - 2] = 0.;
0125     m_corners[3 * icorner - 1] = length;
0126     icorner = 7;
0127     m_corners[3 * icorner - 3] = 0.;
0128     m_corners[3 * icorner - 2] = 0.;
0129     m_corners[3 * icorner - 1] = length;
0130     icorner = 8;
0131     m_corners[3 * icorner - 3] = 0.;
0132     m_corners[3 * icorner - 2] = rear;
0133     m_corners[3 * icorner - 1] = length;
0134
0135     calculateCentres();
0136
0137     // Move centre of SC to (0,0,0)
0138
0139     translate();
0140
0141     // Rotate into standard position (face centres on z axis)
0142
0143     //  this->rotate();
0144
0145     calculateCentres();
0146   }
0147 #if 0
0148   void DDEcalEndcapTrap::rotate(const DDTranslation& /* frontCentre */, const DDTranslation& /* rearCentre */) {
0149     //
0150     //  Rotate supercrystal to bring front and rear face centres to specified points
0151     //
0152   }
0153 #endif
0154   void DDEcalEndcapTrap::rotate(const DDRotationMatrix& rot) {
0155     //
0156     //  Rotate supercrystal by specified rotation about (0,0,0)
0157     //
0158
0159     int icorner;
0160     DDTranslation cc;
0161     for (icorner = 1; icorner <= 8; icorner++) {
0162       cc = cornerPos(icorner);
0163       cc = rot * cc;
0164       cornerPos(icorner, cc);
0165     }
0166     m_rotation = rot * m_rotation;
0167     calculateCentres();
0168   }
0169
0170   void DDEcalEndcapTrap::translate() {
0171     translate(-1. * centrePos());
0172   }
0173
0174   void DDEcalEndcapTrap::translate(const DDTranslation& trans) {
0175     //
0176     //  Translate supercrystal by specified amount
0177     //
0178
0179     DDTranslation tcorner;
0180     for (int icorner = 1; icorner <= 8; icorner++) {
0181       tcorner = cornerPos(icorner) + trans;
0182       cornerPos(icorner, tcorner);
0183     }
0184     calculateCentres();
0185     m_translation = trans + m_translation;
0186   }
0187
0188   void DDEcalEndcapTrap::moveto(const DDTranslation& frontCentre, const DDTranslation& rearCentre) {
0189     //
0190     //  Rotate (about X then about Y) and translate supercrystal to bring axis joining front and rear face centres parallel to line connecting specified points
0191     //
0192
0193     //  Get azimuthal and polar angles of current axis and target axis
0194     double currentTheta = elevationAngle();
0195     double currentPhi = polarAngle();
0196     double targetTheta = elevationAngle(frontCentre - rearCentre);
0197     double targetPhi = polarAngle(frontCentre - rearCentre);
0198
0199     //  Rotate to correct angle (X then Y)
0200     rotateX(targetTheta - currentTheta);
0201     rotateY(targetPhi - currentPhi);
0202
0203     //  Translate SC to final position
0204     DDTranslation targetCentre = 0.5 * (frontCentre + rearCentre);
0205     translate(targetCentre - centrePos());
0206   }
0207
0208   void DDEcalEndcapTrap::rotateX(const double angle) {
0209     //
0210     //  Rotate SC through given angle about X axis
0211     //
0212
0213     const CLHEP::HepRotation tmp(CLHEP::Hep3Vector(1., 0., 0.), angle);
0214
0215     rotate(DDRotationMatrix(tmp.xx(), tmp.xy(), tmp.xz(), tmp.yx(), tmp.yy(), tmp.yz(), tmp.zx(), tmp.zy(), tmp.zz()));
0216   }
0217
0218   void DDEcalEndcapTrap::rotateY(const double angle) {
0219     //
0220     //  Rotate SC through given angle about Y axis
0221     //
0222     const CLHEP::HepRotation tmp(CLHEP::Hep3Vector(0., 1., 0.), angle);
0223
0224     rotate(DDRotationMatrix(tmp.xx(), tmp.xy(), tmp.xz(), tmp.yx(), tmp.yy(), tmp.yz(), tmp.zx(), tmp.zy(), tmp.zz()));
0225   }
0226
0227   void DDEcalEndcapTrap::calculateCentres() {
0228     //
0229     //  Calculate crystal centre and front & rear face centres
0230     //
0231
0232     int ixyz, icorner;
0233
0234     for (ixyz = 0; ixyz < 3; ixyz++) {
0235       m_centre[ixyz] = 0;
0236       m_fcentre[ixyz] = 0;
0237       m_rcentre[ixyz] = 0;
0238     }
0239
0240     for (icorner = 1; icorner <= 4; icorner++) {
0241       for (ixyz = 0; ixyz < 3; ixyz++) {
0242         m_centre[ixyz] = m_centre[ixyz] + 0.125 * m_corners[3 * icorner - 3 + ixyz];
0243         m_fcentre[ixyz] = m_fcentre[ixyz] + 0.25 * m_corners[3 * icorner - 3 + ixyz];
0244       }
0245     }
0246     for (icorner = 5; icorner <= 8; icorner++) {
0247       for (ixyz = 0; ixyz < 3; ixyz++) {
0248         m_centre[ixyz] = m_centre[ixyz] + 0.125 * m_corners[3 * icorner - 3 + ixyz];
0249         m_rcentre[ixyz] = m_rcentre[ixyz] + 0.25 * m_corners[3 * icorner - 3 + ixyz];
0250       }
0251     }
0252   }
0253
0254   DDTranslation DDEcalEndcapTrap::cornerPos(const int icorner) {
0255     //
0256     //  Return specified corner as a DDTranslation
0257     //
0258     return DDTranslation(m_corners[3 * icorner - 3], m_corners[3 * icorner - 2], m_corners[3 * icorner - 1]);
0259   }
0260
0261   void DDEcalEndcapTrap::cornerPos(const int icorner, const DDTranslation& cornerxyz) {
0262     //
0263     //  Save position of specified corner.
0264     //
0265     for (int ixyz = 0; ixyz < 3; ixyz++) {
0266       m_corners[3 * icorner - 3 + ixyz] = (0 == ixyz ? cornerxyz.x() : (1 == ixyz ? cornerxyz.y() : cornerxyz.z()));
0267       ;
0268     }
0269   }
0270
0271   DDTranslation DDEcalEndcapTrap::centrePos() {
0272     //
0273     //  Return SC centre as a DDTranslation
0274     //
0275     return DDTranslation(m_centre[0], m_centre[1], m_centre[2]);
0276   }
0277
0278   DDTranslation DDEcalEndcapTrap::fcentrePos() {
0279     //
0280     //  Return SC front face centre as a DDTranslation
0281     //
0282     return DDTranslation(m_fcentre[0], m_fcentre[1], m_fcentre[2]);
0283   }
0284
0285   DDTranslation DDEcalEndcapTrap::rcentrePos() {
0286     //
0287     //  Return SC rear face centre as a DDTranslation
0288     //
0289     return DDTranslation(m_rcentre[0], m_rcentre[1], m_rcentre[2]);
0290   }
0291
0292   double DDEcalEndcapTrap::elevationAngle(const DDTranslation& trans) {
0293     //
0294     //  Return elevation angle (out of x-z plane) of a given translation (seen as a vector from the origin).
0295     //
0296     double sintheta = trans.y() / trans.r();
0297     return asin(sintheta);
0298   }
0299
0300   double DDEcalEndcapTrap::elevationAngle() {
0301     //
0302     //  Return elevation angle (out of x-z plane) of SC in current position.
0303     //
0304     DDTranslation current = fcentrePos() - rcentrePos();
0305     return elevationAngle(current);
0306   }
0307
0308   double DDEcalEndcapTrap::polarAngle(const DDTranslation& trans) {
0309     //
0310     //  Return polar angle (from x to z) of a given translation (seen as a vector from the origin).
0311     //
0312     double tanphi = trans.x() / trans.z();
0313     return atan(tanphi);
0314   }
0315
0316   double DDEcalEndcapTrap::polarAngle() {
0317     //
0318     //  Return elevation angle (out of x-z plane) of SC in current position.
0319     //
0320     DDTranslation current = fcentrePos() - rcentrePos();
0321     return polarAngle(current);
0322   }
0323 #if 0
0324   void DDEcalEndcapTrap::print() {
0325     //
0326     //  Print SC coordinates for debugging
0327     //
0328     for (int ic = 1; ic <= 8; ic++) {
0329       /* DDTranslation cc = */  cornerPos(ic);
0330     }
0331   }
0332 #endif
0333   namespace {
0334     struct Endcap {
0335       string mat;
0336       double zOff;
0337
0338       string quaName;
0339       string quaMat;
0340
0341       string crysMat;
0342       string wallMat;
0343
0344       double crysLength;
0345       double crysRear;
0346       double crysFront;
0347       double sCELength;
0348       double sCERear;
0349       double sCEFront;
0350       double sCALength;
0351       double sCARear;
0352       double sCAFront;
0353       double sCAWall;
0354       double sCHLength;
0355       double sCHSide;
0356
0357       double nSCTypes;
0358       vector<double> vecEESCProf;
0359       double nColumns;
0360       vector<double> vecEEShape;
0361       double nSCCutaway;
0362       vector<double> vecEESCCutaway;
0363       double nSCquad;
0364       vector<double> vecEESCCtrs;
0365       double nCRSC;
0366       vector<double> vecEECRCtrs;
0367
0368       array<double, 3> cutParms;
0369       string cutBoxName;
0370
0371       string envName;
0372       string alvName;
0373       string intName;
0374       string cryName;
0375
0376       DDTranslation cryFCtr[5][5];
0377       DDTranslation cryRCtr[5][5];
0378       DDTranslation scrFCtr[10][10];
0379       DDTranslation scrRCtr[10][10];
0380
0381       double pFHalf;
0382       double pFFifth;
0383       double pF45;
0384
0385       vector<double> vecEESCLims;
0386
0387       double iLength;
0388       double iXYOff;
0389       double cryZOff;
0390       double zFront;
0391     };
0392
0393     const Rotation3D& myrot(dd4hep::cms::Namespace& ns, const string& nam, const Rotation3D& r) {
0394       ns.addRotation(nam, r);
0395       return ns.rotation(ns.prepend(nam));
0396     }
0397
0398     string_view mynamespace(string_view input) {
0399       string_view v = input;
0400       auto trim_pos = v.find(':');
0401       if (trim_pos != v.npos)
0402         v.remove_suffix(v.size() - (trim_pos + 1));
0403       return v;
0404     }
0405   }  // namespace
0406
0407   static long algorithm(dd4hep::Detector& /* description */, dd4hep::cms::ParsingContext& ctxt, xml_h e,
0408                         SensitiveDetector& /* sens */) {
0409     dd4hep::cms::Namespace ns(ctxt, e, true);
0410     dd4hep::cms::AlgoArguments args(ctxt, e);
0411
0412     // TRICK!
0413     string myns{mynamespace(args.parentName()).data(), mynamespace(args.parentName()).size()};
0414
0415     Endcap ee;
0416     ee.mat = args.str("EEMat");
0417     ee.zOff = args.dble("EEzOff");
0418
0419     ee.quaName = args.str("EEQuaName");
0420     ee.quaMat = args.str("EEQuaMat");
0421     ee.crysMat = args.str("EECrysMat");
0422     ee.wallMat = args.str("EEWallMat");
0423     ee.crysLength = args.dble("EECrysLength");
0424     ee.crysRear = args.dble("EECrysRear");
0425     ee.crysFront = args.dble("EECrysFront");
0426     ee.sCELength = args.dble("EESCELength");
0427     ee.sCERear = args.dble("EESCERear");
0428     ee.sCEFront = args.dble("EESCEFront");
0429     ee.sCALength = args.dble("EESCALength");
0430     ee.sCARear = args.dble("EESCARear");
0431     ee.sCAFront = args.dble("EESCAFront");
0432     ee.sCAWall = args.dble("EESCAWall");
0433     ee.sCHLength = args.dble("EESCHLength");
0434     ee.sCHSide = args.dble("EESCHSide");
0435     ee.nSCTypes = args.dble("EEnSCTypes");
0436     ee.nColumns = args.dble("EEnColumns");
0437     ee.nSCCutaway = args.dble("EEnSCCutaway");
0438     ee.nSCquad = args.dble("EEnSCquad");
0439     ee.nCRSC = args.dble("EEnCRSC");
0440     ee.vecEESCProf = args.vecDble("EESCProf");
0441     ee.vecEEShape = args.vecDble("EEShape");
0442     ee.vecEESCCutaway = args.vecDble("EESCCutaway");
0443     ee.vecEESCCtrs = args.vecDble("EESCCtrs");
0444     ee.vecEECRCtrs = args.vecDble("EECRCtrs");
0445
0446     ee.cutBoxName = args.str("EECutBoxName");
0447
0448     ee.envName = args.str("EEEnvName");
0449     ee.alvName = args.str("EEAlvName");
0450     ee.intName = args.str("EEIntName");
0451     ee.cryName = args.str("EECryName");
0452
0453     ee.pFHalf = args.dble("EEPFHalf");
0454     ee.pFFifth = args.dble("EEPFFifth");
0455     ee.pF45 = args.dble("EEPF45");
0456
0457     ee.vecEESCLims = args.vecDble("EESCLims");
0458     ee.iLength = args.dble("EEiLength");
0459     ee.iXYOff = args.dble("EEiXYOff");
0460     ee.cryZOff = args.dble("EECryZOff");
0461     ee.zFront = args.dble("EEzFront");
0462
0463     //  Position supercrystals in EE Quadrant
0464
0465     //********************************* cutbox for trimming edge SCs
0466     const double cutWid(ee.sCERear / sqrt(2.));
0467     ee.cutParms[0] = cutWid;
0468     ee.cutParms[1] = cutWid;
0469     ee.cutParms[2] = ee.sCELength / sqrt(2.);
0470     Solid eeCutBox = Box(ee.cutBoxName, ee.cutParms[0], ee.cutParms[1], ee.cutParms[2]);
0471     //**************************************************************
0472
0473     const double zFix(ee.zFront - 3172.0_mm);  // fix for changing z offset
0474
0475     //** fill supercrystal front and rear center positions from xml input
0476     for (unsigned int iC(0); iC != (unsigned int)ee.nSCquad; ++iC) {
0477       const unsigned int iOff(8 * iC);
0478       const unsigned int ix((unsigned int)ee.vecEESCCtrs[iOff + 0]);
0479       const unsigned int iy((unsigned int)ee.vecEESCCtrs[iOff + 1]);
0480
0481       assert(ix > 0 && ix < 11 && iy > 0 && iy < 11);
0482
0483       ee.scrFCtr[ix - 1][iy - 1] =
0484         DDTranslation(ee.vecEESCCtrs[iOff + 2], ee.vecEESCCtrs[iOff + 4], ee.vecEESCCtrs[iOff + 6] + zFix);
0485
0486       ee.scrRCtr[ix - 1][iy - 1] =
0487         DDTranslation(ee.vecEESCCtrs[iOff + 3], ee.vecEESCCtrs[iOff + 5], ee.vecEESCCtrs[iOff + 7] + zFix);
0488     }
0489
0490     //** fill crystal front and rear center positions from xml input
0491     for (unsigned int iC(0); iC != 25; ++iC) {
0492       const unsigned int iOff(8 * iC);
0493       const unsigned int ix((unsigned int)ee.vecEECRCtrs[iOff + 0]);
0494       const unsigned int iy((unsigned int)ee.vecEECRCtrs[iOff + 1]);
0495
0496       assert(ix > 0 && ix < 6 && iy > 0 && iy < 6);
0497
0498       ee.cryFCtr[ix - 1][iy - 1] =
0499         DDTranslation(ee.vecEECRCtrs[iOff + 2], ee.vecEECRCtrs[iOff + 4], ee.vecEECRCtrs[iOff + 6]);
0500
0501       ee.cryRCtr[ix - 1][iy - 1] =
0502         DDTranslation(ee.vecEECRCtrs[iOff + 3], ee.vecEECRCtrs[iOff + 5], ee.vecEECRCtrs[iOff + 7]);
0503     }
0504
0505     Solid eeCRSolid = Trap(ee.cryName,
0506                            0.5 * ee.crysLength,
0507                            atan((ee.crysRear - ee.crysFront) / (sqrt(2.) * ee.crysLength)),
0508                            45._deg,
0509                            0.5 * ee.crysFront,
0510                            0.5 * ee.crysFront,
0511                            0.5 * ee.crysFront,
0512                            0._deg,
0513                            0.5 * ee.crysRear,
0514                            0.5 * ee.crysRear,
0515                            0.5 * ee.crysRear,
0516                            0._deg);
0517     Volume eeCRLog = Volume(ee.cryName, eeCRSolid, ns.material(ee.crysMat));
0518
0519     for (unsigned int isc(0); isc < ee.nSCTypes; ++isc) {
0520       unsigned int iSCType = isc + 1;
0521       const string anum(std::to_string(iSCType));
0522       const double eFront(0.5 * ee.sCEFront);
0523       const double eRear(0.5 * ee.sCERear);
0524       const double eAng(atan((ee.sCERear - ee.sCEFront) / (sqrt(2.) * ee.sCELength)));
0525       const double ffived(45._deg);
0526       const double zerod(0._deg);
0527       string eeSCEnvName(1 == iSCType ? ee.envName + std::to_string(iSCType)
0528                          : (ee.envName + std::to_string(iSCType) + "Tmp"));
0529       Solid eeSCEnv = ns.addSolidNS(
0530                                     eeSCEnvName,
0531                                     Trap(eeSCEnvName, 0.5 * ee.sCELength, eAng, ffived, eFront, eFront, eFront, zerod, eRear, eRear, eRear, zerod));
0532
0533       const double aFront(0.5 * ee.sCAFront);
0534       const double aRear(0.5 * ee.sCARear);
0535       const double aAng(atan((ee.sCARear - ee.sCAFront) / (sqrt(2.) * ee.sCALength)));
0536       string eeSCAlvName(
0537                          (1 == iSCType ? ee.alvName + std::to_string(iSCType) : (ee.alvName + std::to_string(iSCType) + "Tmp")));
0538       Solid eeSCAlv = ns.addSolidNS(
0539                                     eeSCAlvName,
0540                                     Trap(eeSCAlvName, 0.5 * ee.sCALength, aAng, ffived, aFront, aFront, aFront, zerod, aRear, aRear, aRear, zerod));
0541
0542       const double dwall(ee.sCAWall);
0543       const double iFront(aFront - dwall);
0544       const double iRear(iFront);
0545       const double iLen(ee.iLength);
0546       string eeSCIntName(1 == iSCType ? ee.intName + std::to_string(iSCType)
0547                          : (ee.intName + std::to_string(iSCType) + "Tmp"));
0548       Solid eeSCInt = ns.addSolidNS(eeSCIntName,
0549                                     Trap(eeSCIntName,
0550                                          iLen / 2.,
0551                                          atan((ee.sCARear - ee.sCAFront) / (sqrt(2.) * ee.sCALength)),
0552                                          ffived,
0553                                          iFront,
0554                                          iFront,
0555                                          iFront,
0556                                          zerod,
0557                                          iRear,
0558                                          iRear,
0559                                          iRear,
0560                                          zerod));
0561
0562       const double dz(-0.5 * (ee.sCELength - ee.sCALength));
0563       const double dxy(0.5 * dz * (ee.sCERear - ee.sCEFront) / ee.sCELength);
0564       const double zIOff(-(ee.sCALength - iLen) / 2.);
0565       const double xyIOff(ee.iXYOff);
0566
0567       Volume eeSCELog;
0568       Volume eeSCALog;
0569       Volume eeSCILog;
0570
0571       if (1 == iSCType) {  // standard SC in this block
0572         eeSCELog = ns.addVolumeNS(Volume(myns + ee.envName + std::to_string(iSCType), eeSCEnv, ns.material(ee.mat)));
0573         eeSCALog = Volume(ee.alvName + std::to_string(iSCType), eeSCAlv, ns.material(ee.wallMat));
0574         eeSCILog = Volume(ee.intName + std::to_string(iSCType), eeSCInt, ns.material(ee.mat));
0575       } else {  // partial SCs this block: create subtraction volumes as appropriate
0576         const double half(ee.cutParms[0] - ee.pFHalf * ee.crysRear);
0577         const double fifth(ee.cutParms[0] + ee.pFFifth * ee.crysRear);
0578         const double fac(ee.pF45);
0579
0580         const double zmm(0.0_mm);
0581
0582         DDTranslation cutTra(
0583                              2 == iSCType ? DDTranslation(zmm, half, zmm)
0584                              : (3 == iSCType ? DDTranslation(half, zmm, zmm)
0585                                 : (4 == iSCType ? DDTranslation(zmm, -fifth, zmm)
0586                                    : (5 == iSCType ? DDTranslation(-half * fac, -half * fac, zmm)
0587                                       : DDTranslation(-fifth, zmm, zmm)))));
0588
0589         const CLHEP::HepRotationZ cutm(ffived);
0590
0591         Rotation3D cutRot(5 != iSCType ? Rotation3D()
0592                           : myrot(ns,
0593                                   "EECry5Rot",
0594                                   Rotation3D(cutm.xx(),
0595                                              cutm.xy(),
0596                                              cutm.xz(),
0597                                              cutm.yx(),
0598                                              cutm.yy(),
0599                                              cutm.yz(),
0600                                              cutm.zx(),
0601                                              cutm.zy(),
0602                                              cutm.zz())));
0603
0604         Solid eeCutEnv = SubtractionSolid(ee.envName + std::to_string(iSCType),
0605                                           ns.solid(ee.envName + std::to_string(iSCType) + "Tmp"),
0606                                           eeCutBox,
0607                                           Transform3D(cutRot, cutTra));
0608
0609         const DDTranslation extra(dxy, dxy, dz);
0610
0611         Solid eeCutAlv = SubtractionSolid(ee.alvName + std::to_string(iSCType),
0612                                           ns.solid(ee.alvName + std::to_string(iSCType) + "Tmp"),
0613                                           eeCutBox,
0614                                           Transform3D(cutRot, cutTra - extra));
0615
0616         const double mySign(iSCType < 4 ? +1. : -1.);
0617
0618         const DDTranslation extraI(xyIOff + mySign * 2.0_mm, xyIOff + mySign * 2.0_mm, zIOff);
0619
0620         Solid eeCutInt = SubtractionSolid(ee.intName + std::to_string(iSCType),
0621                                           ns.solid(ee.intName + std::to_string(iSCType) + "Tmp"),
0622                                           eeCutBox,
0623                                           Transform3D(cutRot, cutTra - extraI));
0624
0625         eeSCELog = ns.addVolumeNS(Volume(myns + ee.envName + std::to_string(iSCType), eeCutEnv, ns.material(ee.mat)));
0626         eeSCALog = Volume(ee.alvName + std::to_string(iSCType), eeCutAlv, ns.material(ee.wallMat));
0627         eeSCILog = Volume(ee.intName + std::to_string(iSCType), eeCutInt, ns.material(ee.mat));
0628       }
0629       eeSCELog.placeVolume(eeSCALog, iSCType * 100 + 1, Position(dxy, dxy, dz));
0630       eeSCALog.placeVolume(eeSCILog, iSCType * 100 + 1, Position(xyIOff, xyIOff, zIOff));
0631
0632       // Position crystals within parent supercrystal interior volume
0633       static const unsigned int ncol(5);
0634
0635       if (iSCType > 0 && iSCType <= ee.nSCTypes) {
0636         const unsigned int icoffset((iSCType - 1) * ncol - 1);
0637
0638         // Loop over columns of SC
0639         for (unsigned int icol(1); icol <= ncol; ++icol) {
0640           // Get column limits for this SC type from xml input
0641           const int ncrcol((int)ee.vecEESCProf[icoffset + icol]);
0642
0643           const int imin(0 < ncrcol ? 1 : (0 > ncrcol ? ncol + ncrcol + 1 : 0));
0644           const int imax(0 < ncrcol ? ncrcol : (0 > ncrcol ? ncol : 0));
0645
0646           if (imax > 0) {
0647             // Loop over crystals in this row
0648             for (int irow(imin); irow <= imax; ++irow) {
0649               // Create crystal as a DDEcalEndcapTrap object and calculate rotation and
0650               // translation required to position it in the SC.
0651               DDEcalEndcapTrap crystal(1, ee.crysFront, ee.crysRear, ee.crysLength);
0652
0653               crystal.moveto(ee.cryFCtr[icol - 1][irow - 1], ee.cryRCtr[icol - 1][irow - 1]);
0654
0655               string rname("EECrRoC" + std::to_string(icol) + "R" + std::to_string(irow));
0656
0657               eeSCALog.placeVolume(
0658                                    eeCRLog,
0659                                    100 * iSCType + 10 * (icol - 1) + (irow - 1),
0660                                    Transform3D(
0661                                                myrot(ns, rname, crystal.rotation()),
0662                                                Position(crystal.centrePos().x(), crystal.centrePos().y(), crystal.centrePos().z() - ee.cryZOff)));
0663             }
0664           }
0665         }
0666       }
0667     }
0668
0669     //** Loop over endcap columns
0670     for (int icol = 1; icol <= int(ee.nColumns); icol++) {
0671       //**  Loop over SCs in column, using limits from xml input
0672       for (int irow = int(ee.vecEEShape[2 * icol - 2]); irow <= int(ee.vecEEShape[2 * icol - 1]); ++irow) {
0673         if (ee.vecEESCLims[0] <= icol && ee.vecEESCLims[1] >= icol && ee.vecEESCLims[2] <= irow &&
0674             ee.vecEESCLims[3] >= irow) {
0675           // Find SC type (complete or partial) for this location
0676           unsigned int isctype = 1;
0677
0678           for (unsigned int ii = 0; ii < (unsigned int)(ee.nSCCutaway); ++ii) {
0679             if ((ee.vecEESCCutaway[3 * ii] == icol) && (ee.vecEESCCutaway[3 * ii + 1] == irow)) {
0680               isctype = int(ee.vecEESCCutaway[3 * ii + 2]);
0681             }
0682           }
0683
0684           // Create SC as a DDEcalEndcapTrap object and calculate rotation and
0685           // translation required to position it in the endcap.
0686           DDEcalEndcapTrap scrys(1, ee.sCEFront, ee.sCERear, ee.sCELength);
0687           scrys.moveto(ee.scrFCtr[icol - 1][irow - 1], ee.scrRCtr[icol - 1][irow - 1]);
0688           scrys.translate(DDTranslation(0., 0., -ee.zOff));
0689
0690           string rname(ee.envName + std::to_string(isctype) + std::to_string(icol) + "R" + std::to_string(irow));
0691           // Position SC in endcap
0692           Volume quaLog = ns.volume(ee.quaName);
0693           Volume childEnvLog = ns.volume(myns + ee.envName + std::to_string(isctype));
0694           quaLog.placeVolume(childEnvLog,
0695                              100 * isctype + 10 * (icol - 1) + (irow - 1),
0696                              Transform3D(scrys.rotation(), scrys.centrePos()));
0697         }
0698       }
0699     }
0700
0701     return 1;
0702   }
0703 }
0704
0705 DECLARE_DDCMS_DETELEMENT(DDCMS_ecal_DDEcalEndcapAlgo, algorithm)