GeomBuilder/src/buildTruncTubs.icc

0001 /*
0002  *  Compute parameters for a truncated cylinder section.
0003  *  In the current geometry TruncTubs are described starting from startPhi = 0,
0004  *  so that the volume at ~pi/2 has a nonzero rotation, contrary to any other
0005  *  volume type in the same geometry. This is not a problem as long as proper
0006  *  local to global transformation are used.
0007  *
0008  *
0009  *  \author N. Amapane - INFN Torino
0010  */
0011
0012 inline void volumeHandle::buildTruncTubs(double zhalf,
0013                                          double rIn,
0014                                          double rOut,
0015                                          double startPhi,
0016                                          double deltaPhi,
0017                                          double cutAtStart,
0018                                          double cutAtDelta,
0019                                          bool cutInside) {
0020   LogTrace("MagGeoBuilder") << "Building TruncTubs surfaces...: ";
0021   LogTrace("MagGeoBuilder") << "zhalf      " << zhalf << newln << "rIn        " << rIn << newln << "rOut       " << rOut
0022                             << newln << "startPhi   " << startPhi << newln << "deltaPhi   " << deltaPhi << newln
0023                             << "cutAtStart " << cutAtStart << newln << "cutAtDelta " << cutAtDelta << newln
0024                             << "cutInside  " << cutInside;
0025
0026   //
0027   //  GlobalVector planeXAxis = refPlane->toGlobal(LocalVector( 1, 0, 0));
0028   //  GlobalVector planeYAxis = refPlane->toGlobal(LocalVector( 0, 1, 0));
0029   GlobalVector planeZAxis = refPlane->toGlobal(LocalVector(0, 0, 1));
0030
0031   Sides cyl_side;
0032   Sides plane_side;
0033   double rCyl = 0;
0034   double rCentr = 0;
0035   if (cutInside) {
0036     cyl_side = outer;
0037     rCyl = rOut;
0038     plane_side = inner;
0039     rCentr = (max(max(rIn, cutAtStart), cutAtDelta) + rOut) / 2.;
0040   } else {
0041     cyl_side = inner;
0042     rCyl = rIn;
0043     plane_side = outer;
0044     rCentr = (rIn + min(min(rOut, cutAtStart), cutAtDelta)) / 2.;
0045   }
0046
0047   // Recalculate center: (for a DDTruncTubs, DDD gives 0,0,Z)
0048   // The R of center is in the middle of the arc of cylinder fully contained
0049   // in the trunctubs.
0050   Geom::Phi<double> phiCenter(startPhi + deltaPhi / 2.);
0051   center_ = refPlane->toGlobal(LocalPoint(rCentr * cos(phiCenter), rCentr * sin(phiCenter), 0.));
0052
0053   // FIXME!! Actually should recompute RN from pos_Rplane; should not
0054   // matter anyhow
0055   theRN = rCentr;
0056
0057   // For simplicity, the position of the cut plane is taken at one of the
0058   // phi edges, where R is known
0059   // FIXME! move to center of plane
0060   // FIXME: compute. in double prec (not float)
0061   GlobalPoint pos_Rplane(refPlane->toGlobal(LocalPoint(LocalPoint::Cylindrical(cutAtStart, startPhi, 0.))));
0062
0063   // Compute angle of the cut plane
0064   // Got any easier formula?
0065   // FIXME: check that it still holds for cutAtStart > cutAtDelta
0066   double c = sqrt(cutAtDelta * cutAtDelta + cutAtStart * cutAtStart - 2 * cutAtDelta * cutAtStart * cos(deltaPhi));
0067   double alpha = startPhi - asin(sin(deltaPhi) * cutAtDelta / c);
0068   GlobalVector x_Rplane = refPlane->toGlobal(LocalVector(cos(alpha), sin(alpha), 0));
0069   Surface::RotationType rot_R(planeZAxis, x_Rplane);
0070
0071   // FIXME: use builder
0072   surfaces[plane_side] = new Plane(pos_Rplane, rot_R);
0073   surfaces[cyl_side] = new Cylinder(rCyl, Surface::PositionType(0, 0, center_.z()), Surface::RotationType());
0074
0075   // Build lateral surfaces.
0076   // Note that with the choice of center descrived above, the
0077   // plane position (origin of r.f.) of the smallest phi face
0078   // will be its center, while this is not true for the largest phi face.
0079   buildPhiZSurf(startPhi, deltaPhi, zhalf, rCentr);
0080
0081   if (debug) {
0082     LogTrace("MagGeoBuilder") << "pos_Rplane    " << pos_Rplane << " " << pos_Rplane.perp() << " " << pos_Rplane.phi()
0083                               << newln << "rot_R         " << surfaces[plane_side]->toGlobal(LocalVector(0., 0., 1.))
0084                               << " phi " << surfaces[plane_side]->toGlobal(LocalVector(0., 0., 1.)).phi() << newln
0085                               << "cyl radius     " << rCyl;
0086
0087     //   // Check ordering.
0088     if ((pos_Rplane.perp() < rCyl) != cutInside) {
0089       LogTrace("MagGeoBuilder") << "*** WARNING: pos_outer < pos_inner ";
0090     }
0091   }
0092   // Save volume boundaries
0093   theRMin = rIn;
0094   theRMax = rOut;
0095   thePhiMin = surfaces[phiminus]->position().phi();
0096 }