GeomBuilder/src/buildPseudoTrap.icc

0001 /*
0002  *  Compute parameters for a pseudotrapezoid .
0003  *  NOTE: The pseudotrapezoid is no longer supported in DD4hep.
0004  *  It is handled here for legacy geonetries that are still in use.
0005  *  For pseudotraps the refPlane is parallel to beam line; transformation:
0006  *
0007  *  Global(for vol at pi/2)    Local
0008  *  +R (+Y)                    +Z
0009  *  +phi(-X)                   +X
0010  *  +Z                         +Y
0011  *
0012  *  \author N. Amapane - INFN Torino
0013  */
0014
0015 inline void volumeHandle::buildPseudoTrap(
0016     double x1, double x2, double y1, double y2, double halfZ, double radius, bool atMinusZ) {
0017   LogTrace("MagGeoBuilder") << "Building PseudoTrap surfaces...: ";
0018   LogTrace("MagGeoBuilder") << "halfZ     " << halfZ << newln << "x1        " << x1 << newln << "x2        " << x2
0019                             << newln << "y1        " << y1 << newln << "y2        " << y2 << newln << "radius    "
0020                             << radius << newln << "atMinusZ  " << atMinusZ;
0021
0022   // Just check assumptions on parameters...
0023   const double epsilon = 1e-5;
0024   if (debug) {
0025     if (y1 - y2 > epsilon) {
0026       LogTrace("MagGeoBuilder") << "*** WARNING: unexpected pseudotrapezoid parameters.";
0027     }
0028
0029     // Check that volume is convex (no concave volume in current geometry...)
0030     if (radius * (atMinusZ ? -1. : +1) < 0.) {
0031       LogTrace("MagGeoBuilder") << "*** WARNING: pseudotrapezoid is concave";
0032     }
0033   }
0034   // CAVEAT: pseudotraps are rotated in a different way as traps,
0035   // since they have local Z along global R...
0036   GlobalVector planeXAxis = refPlane->toGlobal(LocalVector(1, 0, 0));
0037   GlobalVector planeYAxis = refPlane->toGlobal(LocalVector(0, 1, 0));
0038   GlobalVector planeZAxis = refPlane->toGlobal(LocalVector(0, 0, 1));
0039
0040   //FIXME assumption: here we do use the assumption on the orientation /
0041   //of local Z (see above)
0042   GlobalVector Rvol(refPlane->position().x(), refPlane->position().y(), refPlane->position().z());
0043   theRN = std::abs(planeZAxis.dot(Rvol));
0044
0045   double fR = std::abs(radius);
0046   Sides cyl_side;
0047   Sides plane_side;
0048   if (atMinusZ) {
0049     cyl_side = inner;
0050     plane_side = outer;
0051   } else {
0052     cyl_side = outer;
0053     plane_side = inner;
0054   }
0055   GlobalPoint pos_Rplane(refPlane->toGlobal(LocalPoint(0., 0., (atMinusZ ? +halfZ : -halfZ))));
0056   GlobalPoint pos_zplus(refPlane->toGlobal(LocalPoint(0., y1, 0.)));
0057   GlobalPoint pos_zminus(refPlane->toGlobal(LocalPoint(0., -y1, 0.)));
0058   double halfX((x1 + x2) / 2.);
0059   GlobalPoint pos_phiplus(refPlane->toGlobal(LocalPoint(+halfX, 0., 0.)));
0060   GlobalPoint pos_phiminus(refPlane->toGlobal(LocalPoint(-halfX, 0., 0.)));
0061
0062   //Check that cylinder is centered on beam axis...
0063   float rcheck;
0064   if (atMinusZ)
0065     rcheck = refPlane->toGlobal(LocalPoint(x1, 0., -halfZ)).perp();
0066   else
0067     rcheck = refPlane->toGlobal(LocalPoint(x2, 0., +halfZ)).perp();
0068
0069   if (debug) {
0070     if (std::abs(rcheck - fR) > 100. * epsilon) {  //FIXME!
0071       std::stringstream formatOutput;
0072       formatOutput << setprecision(10) << rcheck << " " << fR << " " << std::abs(rcheck - fR);
0073       LogTrace("MagGeoBuilder") << "*** WARNING: Cylinder surface not centered on beam axis " << formatOutput.str();
0074     }
0075   }
0076
0077   if (debug) {
0078     LogTrace("MagGeoBuilder") << "RN           " << theRN << newln << "pos_Rplane   " << pos_Rplane << " "
0079                               << pos_Rplane.perp() << " " << pos_Rplane.phi() << newln << "pos_zplus    " << pos_zplus
0080                               << " " << pos_zplus.perp() << " " << pos_zplus.phi() << newln << "pos_zminus   "
0081                               << pos_zminus << " " << pos_zminus.perp() << " " << pos_zminus.phi() << newln
0082                               << "pos_phiplus  " << pos_phiplus << " " << pos_phiplus.perp() << " " << pos_phiplus.phi()
0083                               << newln << "pos_phiminus " << pos_phiminus << " " << pos_phiminus.perp() << " "
0084                               << pos_phiminus.phi();
0085   }
0086
0087   // Check ordering.
0088   if ((pos_Rplane.perp() < radius) == atMinusZ) {
0089     LogTrace("MagGeoBuilder") << "*** WARNING: pos_outer < pos_inner for pseudotrapezoid";
0090   }
0091   if (pos_zplus.z() < pos_zminus.z()) {
0092     LogTrace("MagGeoBuilder") << "*** WARNING: pos_zplus < pos_zminus for pseudotrapezoid";
0093   }
0094   if (Geom::Phi<float>(pos_phiplus.phi() - pos_phiminus.phi()) < 0.) {
0095     LogTrace("MagGeoBuilder") << "*** WARNING: pos_phiplus < pos_phiminus for pseudotrapezoid";
0096   }
0097
0098   GlobalVector z_phiplus = (refPlane->toGlobal(LocalVector((x2 - x1) / 2., 0., halfZ))).unit();
0099   GlobalVector z_phiminus = (refPlane->toGlobal(LocalVector(-(x2 - x1) / 2., 0., halfZ))).unit();
0100
0101   LogTrace("MagGeoBuilder") << "z_phiplus  " << z_phiplus << " " << z_phiplus.phi() << newln << "z_phiminus "
0102                             << z_phiminus << " " << z_phiminus.phi();
0103
0104   Surface::RotationType rot_R(planeYAxis, planeXAxis);
0105   Surface::RotationType rot_Z(planeXAxis, planeZAxis);
0106   Surface::RotationType rot_phiplus(planeYAxis, z_phiplus);
0107   Surface::RotationType rot_phiminus(planeYAxis, z_phiminus);
0108
0109   // FIXME: use builder
0110   surfaces[plane_side] = new Plane(pos_Rplane, rot_R);
0111   surfaces[cyl_side] = new Cylinder(fR, Surface::PositionType(0, 0, center_.z()), Surface::RotationType());
0112
0113   surfaces[zplus] = new Plane(pos_zplus, rot_Z);
0114   surfaces[zminus] = new Plane(pos_zminus, rot_Z);
0115   surfaces[phiplus] = new Plane(pos_phiplus, rot_phiplus);
0116   surfaces[phiminus] = new Plane(pos_phiminus, rot_phiminus);
0117
0118   LogTrace("MagGeoBuilder") << "rot_R " << surfaces[plane_side]->toGlobal(LocalVector(0., 0., 1.)) << newln << "rot_Z "
0119                             << surfaces[zplus]->toGlobal(LocalVector(0., 0., 1.)) << newln << "rot_phi+ "
0120                             << surfaces[phiplus]->toGlobal(LocalVector(0., 0., 1.)) << " phi "
0121                             << surfaces[phiplus]->toGlobal(LocalVector(0., 0., 1.)).phi() << newln << "rot_phi- "
0122                             << surfaces[phiminus]->toGlobal(LocalVector(0., 0., 1.)) << " phi "
0123                             << surfaces[phiminus]->toGlobal(LocalVector(0., 0., 1.)).phi();
0124
0125   // Save volume boundaries
0126   double R1 = std::abs(surfaces[plane_side]->toLocal(GlobalPoint(0, 0, 0)).z());
0127   theRMin = min(fR, R1);
0128   theRMax = max(fR, R1);
0129   // FIXME: use phi of middle plane of phiminus surface. Is not the absolute phimin!
0130   thePhiMin = surfaces[phiminus]->position().phi();
0131 }