plugins/dd4hep/DDAngular.cc

0001 #include "DD4hep/DetFactoryHelper.h"
0002 #include "DataFormats/Math/interface/CMSUnits.h"
0003 #include "DetectorDescription/DDCMS/interface/DDPlugins.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005 #include <Math/Rotation3D.h>
0006 #include <Math/AxisAngle.h>
0007 #include <Math/DisplacementVector3D.h>
0008
0009 using namespace std;
0010 using namespace dd4hep;
0011 using namespace cms;
0012 using namespace cms_units::operators;
0013
0014 using DD3Vector = ROOT::Math::DisplacementVector3D<ROOT::Math::Cartesian3D<double> >;
0015 using DDAxisAngle = ROOT::Math::AxisAngle;
0016
0017 namespace {
0018   DD3Vector fUnitVector(double theta, double phi) {
0019     return DD3Vector(cos(phi) * sin(theta), sin(phi) * sin(theta), cos(theta));
0020   }
0021 }  // namespace
0022
0023 static long algorithm(Detector& /* description */, cms::DDParsingContext& ctxt, xml_h e) {
0024   cms::DDNamespace ns(ctxt, e, true);
0025   DDAlgoArguments args(ctxt, e);
0026
0027   int n = args.value<int>("N");
0028   int startCopyNo = args.find("StartCopyNo") ? args.value<int>("StartCopyNo") : 1;
0029   int incrCopyNo = args.find("IncrCopyNo") ? args.value<int>("IncrCopyNo") : 1;
0030   double rangeAngle = args.value<double>("RangeAngle");
0031   double startAngle = args.value<double>("StartAngle");
0032   double radius = args.value<double>("Radius");
0033   vector<double> center = args.value<vector<double> >("Center");
0034   vector<double> rotateSolid = args.value<vector<double> >("RotateSolid");
0035   Volume mother = ns.volume(args.parentName());
0036   string childName = args.value<string>("ChildName");
0037   childName = ns.prepend(childName);
0038   Volume child = ns.volume(childName);
0039
0040   double delta = 0e0;
0041
0042   if (fabs(rangeAngle - 360.0_deg) < 0.001_deg) {
0043     delta = rangeAngle / double(n);
0044   } else if (n > 1) {
0045     delta = rangeAngle / double(n - 1);
0046   }
0047
0048   LogDebug("DDAlgorithm") << "debug: Parameters for positioning:: n " << n << " Start, Range, Delta "
0049                           << convertRadToDeg(startAngle) << " " << convertRadToDeg(rangeAngle) << " "
0050                           << convertRadToDeg(delta) << " Radius " << radius << " Centre " << center[0] << ", "
0051                           << center[1] << ", " << center[2] << ", Rotate solid " << rotateSolid[0] << ", "
0052                           << rotateSolid[1] << ", " << rotateSolid[2];
0053   LogDebug("DDAlgorithm") << "debug: Parent " << mother.name() << "\tChild " << child.name() << " NameSpace "
0054                           << ns.name();
0055
0056   Rotation3D solidRot = Rotation3D();  // Identity rotation
0057   auto sz = rotateSolid.size();
0058   if (sz % 3) {
0059     LogDebug("DDAlgorithm") << "\trotateSolid must occur 3*n times (defining n subsequent rotations)\n"
0060                             << "\t  currently it appears " << sz << " times!\n";
0061   }
0062   for (unsigned int i = 0; i < sz; i += 3) {
0063     const double thetaValue = rotateSolid[i];
0064     if ((thetaValue > 180._deg) || (thetaValue < 0._deg)) {
0065       LogDebug("DDAlgorithm") << "\trotateSolid \'theta\' must be in range [0,180*deg]\n"
0066                               << "\t  currently it is " << convertRadToDeg(thetaValue) << "*deg in rotateSolid["
0067                               << double(i) << "]!\n";
0068     }
0069     DDAxisAngle temp(fUnitVector(rotateSolid[i], rotateSolid[i + 1]), rotateSolid[i + 2]);
0070     LogDebug("DDAlgorithm") << "  rotsolid[" << i << "] axis=" << temp.Axis()
0071                             << " rot.angle=" << convertRadToDeg(temp.Angle());
0072     solidRot = temp * solidRot;
0073   }
0074
0075   double theta = 90._deg;
0076   int copy = startCopyNo;
0077   double phi = startAngle;
0078   for (int i = 0; i < n; ++i) {
0079     double phix = phi;
0080     double phiy = phix + 90._deg;
0081     double phideg = convertRadToDeg(phix);
0082
0083     Rotation3D rotation = makeRotation3D(theta, phix, theta, phiy, 0., 0.) * solidRot;
0084     string rotstr = ns.nsName(child.name()) + std::to_string(phideg * 10.);
0085     auto irot = ctxt.rotations.find(ns.prepend(rotstr));
0086     if (irot != ctxt.rotations.end()) {
0087       rotation = ns.rotation(ns.prepend(rotstr));
0088     }
0089
0090     double xpos = radius * cos(phi) + center[0];
0091     double ypos = radius * sin(phi) + center[1];
0092     double zpos = center[2];
0093     Position tran(xpos, ypos, zpos);
0094     mother.placeVolume(child, copy, Transform3D(rotation, tran));
0095     LogDebug("DDAlgorithm") << "test " << child.name() << " number " << copy << " positioned in " << mother.name()
0096                             << " at " << tran << " with " << rotstr << ": " << rotation;
0097     copy += incrCopyNo;
0098     phi += delta;
0099   }
0100   return 1;
0101 }
0102
0103 // first argument is the type from the xml file
0104 DECLARE_DDCMS_DETELEMENT(DDCMS_global_DDAngular, algorithm)