TrackerCommonData/plugins/DDCutTubsFromPoints.cc

0001 ///////////////////////////////////////////////////////////////////////////////
0002 // File: DDCutTubsFromPoints.cc
0003 // Description: Create a ring of CutTubs segments from points on the rim.
0004 ///////////////////////////////////////////////////////////////////////////////
0005
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 #include "DetectorDescription/Core/interface/DDCurrentNamespace.h"
0008 #include "DetectorDescription/Core/interface/DDSolid.h"
0009 #include "DetectorDescription/Core/interface/DDVector.h"
0010 #include "DetectorDescription/Core/interface/DDTypes.h"
0011 #include "DetectorDescription/Core/interface/DDAlgorithm.h"
0012 #include "DetectorDescription/Core/interface/DDAlgorithmFactory.h"
0013 #include "DetectorDescription/Core/interface/DDMaterial.h"
0014 #include <CLHEP/Units/GlobalPhysicalConstants.h>
0015 #include <CLHEP/Units/SystemOfUnits.h>
0016
0017 #include <cmath>
0018 #include <map>
0019 #include <string>
0020 #include <vector>
0021
0022 // This algorithm creates a ring made of CutTubs segments from the phi,z points
0023 // of the rings "corners".
0024 // The algorithm only defines and places two copies of a single Solid with the given name.
0025 class DDCutTubsFromPoints : public DDAlgorithm {
0026 public:
0027   //Constructor and Destructor
0028   DDCutTubsFromPoints();
0029   ~DDCutTubsFromPoints() override;
0030
0031   void initialize(const DDNumericArguments& nArgs,
0032                   const DDVectorArguments& vArgs,
0033                   const DDMapArguments& mArgs,
0034                   const DDStringArguments& sArgs,
0035                   const DDStringVectorArguments& vsArgs) override;
0036
0037   void execute(DDCompactView& cpv) override;
0038
0039 private:
0040   struct Section {
0041     double phi;  // phi position of this edge
0042     double z_l;  // -Z end (cuttubs l plane)
0043     double z_t;  // +Z end (cuttubs t plane)
0044     // radius is implicitly r_min
0045   };
0046
0047   double r_min;
0048   double r_max;
0049   double z_pos;
0050   DDMaterial material;
0051
0052   // a segment is produced between each two consecutive sections that have a
0053   // non-zero phi distance. Sections with zero phi distance can be used to
0054   // create sharp jumps.
0055   std::vector<Section> sections;
0056   // this solid will be defined.
0057   DDName solidOutput;
0058 };
0059
0060 DDCutTubsFromPoints::DDCutTubsFromPoints() {
0061   LogDebug("TrackerGeom") << "DDCutTubsFromPoints info: Creating an instance";
0062 }
0063
0064 DDCutTubsFromPoints::~DDCutTubsFromPoints() {}
0065
0066 void DDCutTubsFromPoints::initialize(const DDNumericArguments& nArgs,
0067                                      const DDVectorArguments& vArgs,
0068                                      const DDMapArguments&,
0069                                      const DDStringArguments& sArgs,
0070                                      const DDStringVectorArguments&) {
0071   r_min = nArgs["rMin"];
0072   r_max = nArgs["rMax"];
0073   z_pos = nArgs["zPos"];
0074   material = DDMaterial(sArgs["Material"]);
0075
0076   auto phis_name = DDName(sArgs["Phi"]);
0077   auto z_ls_name = DDName(sArgs["z_l"]);
0078   auto z_ts_name = DDName(sArgs["z_t"]);
0079   DDVector phis(phis_name);
0080   DDVector z_ls(z_ls_name);
0081   DDVector z_ts(z_ts_name);
0082
0083   assert(phis.size() == z_ls.size());
0084   assert(phis.size() == z_ts.size());
0085
0086   for (unsigned i = 0; i < phis.size(); i++) {
0087     Section s = {phis[i], z_ls[i], z_ts[i]};
0088     sections.emplace_back(s);
0089   }
0090   assert(!sections.empty());
0091
0092   solidOutput = DDName(sArgs["SolidName"]);
0093
0094   std::string idNameSpace = DDCurrentNamespace::ns();
0095   DDName parentName = parent().name();
0096   LogDebug("TrackerGeom") << "DDCutTubsFromPoints debug: Parent " << parentName << "\tSolid " << solidOutput
0097                           << " NameSpace " << idNameSpace << "\tnumber of sections " << sections.size();
0098 }
0099
0100 static double square(double x) { return x * x; }
0101
0102 void DDCutTubsFromPoints::execute(DDCompactView& cpv) {
0103   // radius for plane calculations
0104   // We use r_max here, since P3 later has a Z that is always more inside
0105   // than the extreme points. This means the cutting planes have outwards
0106   // slopes in r-Z, and the corner at r_max could stick out of the bounding
0107   // volume otherwise.
0108   double r = r_max;
0109
0110   // min and max z for the placement in the end
0111   double min_z = 1e9;
0112   double max_z = -1e9;
0113
0114   // counter of actually produced segments (excluding skipped ones)
0115   int segment = 0;
0116
0117   // the segments and their corresponding offset (absolute, as in the input)
0118   std::vector<DDSolid> segments;
0119   std::vector<double> offsets;
0120
0121   Section s1 = sections[0];
0122   for (Section s2 : sections) {
0123     if (s1.phi != s2.phi) {
0124       segment++;
0125       // produce segment s1-s2.
0126       DDName segname(solidOutput.name() + "_seg_" + std::to_string(segment), solidOutput.ns());
0127
0128       double phi1 = s1.phi;
0129       double phi2 = s2.phi;
0130
0131       // track the min/max to properly place&align later
0132       if (s2.z_l < min_z)
0133         min_z = s2.z_l;
0134       if (s2.z_t > max_z)
0135         max_z = s2.z_t;
0136
0137       double P1_z_l = s1.z_l;
0138       double P1_z_t = s1.z_t;
0139       double P2_z_l = s2.z_l;
0140       double P2_z_t = s2.z_t;
0141
0142       double P1_x_t = cos(phi1) * r;
0143       double P1_x_l = cos(phi1) * r;
0144       double P1_y_t = sin(phi1) * r;
0145       double P1_y_l = sin(phi1) * r;
0146
0147       double P2_x_t = cos(phi2) * r;
0148       double P2_x_l = cos(phi2) * r;
0149       double P2_y_t = sin(phi2) * r;
0150       double P2_y_l = sin(phi2) * r;
0151
0152       // each cutting plane is defined by P1-3. P1-2 are corners of the
0153       // segment, P3 is at r=0 with the "average" z to get a nice cut.
0154       double P3_z_l = (P1_z_l + P2_z_l) / 2;
0155       double P3_z_t = (P1_z_t + P2_z_t) / 2;
0156
0157       // we only have one dz to position both planes. The anchor is implicitly
0158       // between the P3's, we have to use an offset later to make the segments
0159       // line up correctly.
0160       double dz = (P3_z_t - P3_z_l) / 2;
0161       double offset = (P3_z_t + P3_z_l) / 2;
0162
0163       // the plane is defined by P1-P3 and P2-P3; since P3 is at r=0 we
0164       // only need the z.
0165       double D1_z_l = P1_z_l - P3_z_l;
0166       double D2_z_l = P2_z_l - P3_z_l;
0167
0168       // the normal is then the cross product...
0169       double n_x_l = (P1_y_l * D2_z_l) - (D1_z_l * P2_y_l);
0170       double n_y_l = (D1_z_l * P2_x_l) - (P1_x_l * D2_z_l);
0171       double n_z_l = (P1_x_l * P2_y_l) - (P1_y_l * P2_x_l);
0172
0173       // ... normalized.
0174       // flip the sign here (but not for t) since root wants it like that.
0175       double norm = -sqrt(square(n_x_l) + square(n_y_l) + square(n_z_l));
0176       n_x_l /= norm;
0177       n_y_l /= norm;
0178       n_z_l /= norm;
0179
0180       // same game for the t side.
0181       double D1_z_t = P1_z_t - P3_z_t;
0182       double D2_z_t = P2_z_t - P3_z_t;
0183
0184       double n_x_t = (P1_y_t * D2_z_t) - (D1_z_t * P2_y_t);
0185       double n_y_t = (D1_z_t * P2_x_t) - (P1_x_t * D2_z_t);
0186       double n_z_t = (P1_x_t * P2_y_t) - (P1_y_t * P2_x_t);
0187
0188       norm = sqrt(square(n_x_t) + square(n_y_t) + square(n_z_t));
0189       n_x_t /= norm;
0190       n_y_t /= norm;
0191       n_z_t /= norm;
0192
0193       // the cuttubs wants a delta phi
0194       double dphi = phi2 - phi1;
0195
0196       DDSolid seg =
0197           DDSolidFactory::cuttubs(segname, dz, r_min, r_max, phi1, dphi, n_x_l, n_y_l, n_z_l, n_x_t, n_y_t, n_z_t);
0198       segments.emplace_back(seg);
0199       offsets.emplace_back(offset);
0200     }
0201     s1 = s2;
0202   }
0203
0204   assert(segments.size() >= 2);  // less would be special cases
0205
0206   DDSolid solid = segments[0];
0207   // placment happens relative to the first member of the union
0208   double shift = offsets[0];
0209
0210   for (unsigned i = 1; i < segments.size() - 1; i++) {
0211     // each sub-union needs a name. Well.
0212     DDName unionname(solidOutput.name() + "_uni_" + std::to_string(i + 1), solidOutput.ns());
0213     solid = DDSolidFactory::unionSolid(
0214         unionname, solid, segments[i], DDTranslation(0, 0, offsets[i] - shift), DDRotation());
0215   }
0216
0217   solid = DDSolidFactory::unionSolid(solidOutput,
0218                                      solid,
0219                                      segments[segments.size() - 1],
0220                                      DDTranslation(0, 0, offsets[segments.size() - 1] - shift),
0221                                      DDRotation());
0222
0223   // remove the common offset from the input, to get sth. aligned at z=0.
0224   double offset = -shift + (min_z + (max_z - min_z) / 2);
0225
0226   DDLogicalPart logical(solidOutput, material, solid, DDEnums::support);
0227
0228   // This is not as generic as the name promises, but it is hard to make a
0229   // solid w/o placing it.
0230   DDRotation rot180("pixfwdCommon:Z180");
0231   cpv.position(logical, parent(), 1, DDTranslation(0, 0, z_pos - offset), DDRotation());
0232   cpv.position(logical, parent(), 2, DDTranslation(0, 0, z_pos - offset), rot180);
0233 }
0234
0235 DEFINE_EDM_PLUGIN(DDAlgorithmFactory, DDCutTubsFromPoints, "track:DDCutTubsFromPoints");