GeomBuilder/src/buildTrap.icc

0001 /*
0002  *  Compute parameters for a trapezoid
0003  *
0004  *  \author N. Amapane - INFN Torino
0005  */
0006
0007 inline void volumeHandle::buildTrap(double x1,
0008                                     double x2,
0009                                     double x3,
0010                                     double x4,
0011                                     double y1,
0012                                     double y2,
0013                                     double theta,
0014                                     double phi,
0015                                     double halfZ,
0016                                     double alpha1,
0017                                     double alpha2) {
0018   LogTrace("MagGeoBuilder") << "Building trapezoid surfaces...: ";
0019   LogTrace("MagGeoBuilder") << "x1 " << x1 << newln << "x2 " << x2 << newln << "x3 " << x3 << newln << "x4 " << x4
0020                             << newln << "y1 " << y1 << newln << "y2 " << y2 << newln << "theta  " << theta << newln
0021                             << "phi    " << phi << newln << "halfZ  " << halfZ << newln << "alpha1 " << alpha1 << newln
0022                             << "alpha2 " << alpha2;
0023
0024   // Just check assumptions on parameters...
0025   if (debug) {
0026     const double epsilon = 1e-5;
0027     if (theta > epsilon || phi > epsilon || y1 - y2 > epsilon || x1 - x3 > epsilon || x2 - x4 > epsilon ||
0028         alpha1 - alpha2 > epsilon) {
0029       LogTrace("MagGeoBuilder") << "*** WARNING: unexpected trapezoid parameters.";
0030     }
0031   }
0032
0033   //  Used parameters halfZ, x1, x2, y1, alpha1
0034   GlobalVector planeXAxis = refPlane->toGlobal(LocalVector(1, 0, 0));
0035   GlobalVector planeYAxis = refPlane->toGlobal(LocalVector(0, 1, 0));
0036   GlobalVector planeZAxis = refPlane->toGlobal(LocalVector(0, 0, 1));
0037
0038   // For the correct definition of inner, outer, phiplus, phiminus,
0039   // Zplus, zminus the orientation of axes must be checked.
0040   // Normal (convention in version 85l_030919):
0041   //  planeXAxis in the direction of decreasing global phi;
0042   //  planeZAxis in the direction of global Z;
0043   // => Invert the sign of local X and Z if planeZAxis points to - global Z
0044
0045   // FIXME Assumption: it is assumed that local Y is always along global R
0046   // (true for version 1103l, not necessarily in the future)
0047   GlobalVector Rvol(refPlane->position().x(), refPlane->position().y(), refPlane->position().z());
0048   theRN = std::abs(planeYAxis.dot(Rvol));
0049
0050   if (planeZAxis.z() < 0.) {
0051     x1 = -x1;
0052     x2 = -x2;
0053     halfZ = -halfZ;
0054   }
0055
0056   double halfX((x1 + x2) / 2.);
0057
0058   GlobalPoint pos_outer(refPlane->toGlobal(LocalPoint(0., y1, 0.)));
0059   GlobalPoint pos_inner(refPlane->toGlobal(LocalPoint(0., -y1, 0.)));
0060   GlobalPoint pos_zplus(refPlane->toGlobal(LocalPoint(0., 0., halfZ)));
0061   GlobalPoint pos_zminus(refPlane->toGlobal(LocalPoint(0., 0., -halfZ)));
0062   GlobalPoint pos_phiplus(refPlane->toGlobal(LocalPoint(-halfX, 0., 0.)));
0063   GlobalPoint pos_phiminus(refPlane->toGlobal(LocalPoint(halfX, 0., 0.)));
0064
0065   LogTrace("MagGeoBuilder") << "RN            " << theRN << newln  // RN line needs extra space to line up
0066                             << "pos_outer    " << pos_outer << " " << pos_outer.perp() << " " << pos_outer.phi()
0067                             << newln << "pos_inner    " << pos_inner << " " << pos_inner.perp() << " "
0068                             << pos_inner.phi() << newln << "pos_zplus    " << pos_zplus << " " << pos_zplus.perp()
0069                             << " " << pos_zplus.phi() << newln << "pos_zminus   " << pos_zminus << " "
0070                             << pos_zminus.perp() << " " << pos_zminus.phi() << newln << "pos_phiplus  " << pos_phiplus
0071                             << " " << pos_phiplus.perp() << " " << pos_phiplus.phi() << newln << "pos_phiminus "
0072                             << pos_phiminus << " " << pos_phiminus.perp() << " " << pos_phiminus.phi();
0073
0074   // Check ordering.
0075   if (debug) {
0076     if (pos_outer.perp() < pos_inner.perp()) {
0077       LogTrace("MagGeoBuilder") << "*** WARNING: pos_outer < pos_inner for trapezoid";
0078     }
0079     if (pos_zplus.z() < pos_zminus.z()) {
0080       LogTrace("MagGeoBuilder") << "*** WARNING: pos_zplus < pos_zminus for trapezoid";
0081     }
0082     if (Geom::Phi<float>(pos_phiplus.phi() - pos_phiminus.phi()) < 0.) {
0083       LogTrace("MagGeoBuilder") << "*** WARNING: pos_phiplus < pos_phiminus for trapezoid";
0084     }
0085   }
0086
0087   // Local Y axis of the faces at +-phi. The local X is = global Z.
0088   GlobalVector y_phiplus = (refPlane->toGlobal(LocalVector((tan(alpha1) * y1 - (x2 - x1) / 2.), y1, 0.))).unit();
0089   GlobalVector y_phiminusV = (refPlane->toGlobal(LocalVector((tan(alpha1) * y1 + (x2 - x1) / 2.), y1, 0.)));
0090   GlobalVector y_phiminus = y_phiminusV.unit();
0091
0092   LogTrace("MagGeoBuilder") << "y_phiplus  " << y_phiplus << newln << "y_phiminus " << y_phiminus;
0093
0094   Surface::RotationType rot_R(planeZAxis, planeXAxis);
0095   Surface::RotationType rot_Z(planeXAxis, planeYAxis);
0096   Surface::RotationType rot_phiplus(planeZAxis, y_phiplus);
0097   Surface::RotationType rot_phiminus(planeZAxis, y_phiminus);
0098
0099   // FIXME: use builder
0100   surfaces[outer] = new Plane(pos_outer, rot_R);
0101   surfaces[inner] = new Plane(pos_inner, rot_R);
0102   surfaces[zplus] = new Plane(pos_zplus, rot_Z);
0103   surfaces[zminus] = new Plane(pos_zminus, rot_Z);
0104   surfaces[phiplus] = new Plane(pos_phiplus, rot_phiplus);
0105   surfaces[phiminus] = new Plane(pos_phiminus, rot_phiminus);
0106
0107   // Save volume boundaries
0108   theRMin = std::abs(surfaces[inner]->toLocal(GlobalPoint(0, 0, 0)).z());
0109   theRMax = std::abs(surfaces[outer]->toLocal(GlobalPoint(0, 0, 0)).z());
0110
0111   // Setting "absolute" phi boundary spots a problem in V85l
0112   // e.g. vol 139 (origin: small inconsistency in volumes.txt)
0113   // So let's keep phi at middle phi- plane...
0114   //  FIXME: Remove this workaround once support for v85l is dropped.
0115   if (name.substr(0, 3) == "V_Z") {
0116     thePhiMin = surfaces[phiminus]->position().phi();
0117   } else {
0118     thePhiMin = min((pos_phiminus + y_phiminusV).phi(), (pos_phiminus - y_phiminusV).phi());
0119   }
0120   LogTrace("MagGeoBuilder") << "rot_R    " << surfaces[outer]->toGlobal(LocalVector(0., 0., 1.)) << newln << "rot_Z    "
0121                             << surfaces[zplus]->toGlobal(LocalVector(0., 0., 1.)) << newln << "rot_phi+ "
0122                             << surfaces[phiplus]->toGlobal(LocalVector(0., 0., 1.)) << " phi "
0123                             << surfaces[phiplus]->toGlobal(LocalVector(0., 0., 1.)).phi() << newln << "rot_phi- "
0124                             << surfaces[phiminus]->toGlobal(LocalVector(0., 0., 1.)) << " phi "
0125                             << surfaces[phiminus]->toGlobal(LocalVector(0., 0., 1.)).phi();
0126   if (debug) {
0127     if (std::abs(surfaces[phiminus]->position().phi() - thePhiMin) > 1e-5) {
0128       LogTrace("MagGeoBuilder") << " Trapez. phiMin = " << thePhiMin;
0129     }
0130   }
0131 }