    // @(#)root/geom:$Name:  $:$Id: TGeoTube.cxx,v 1.63 2005/05/25 14:25:16 brun Exp $
    // Author: Andrei Gheata   24/10/01
    // TGeoTube::Contains() and DistFromInside/In() implemented by Mihaela Gheata

    /*************************************************************************
     * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers.               *
     * All rights reserved.                                                  *
     *                                                                       *
     * For the licensing terms see $ROOTSYS/LICENSE.                         *
     * For the list of contributors see $ROOTSYS/README/CREDITS.             *
     *************************************************************************/

    //_____________________________________________________________________________
    // TGeoTube - cylindrical tube class. It takes 3 parameters :
    //            inner radius, outer radius and half-length dz.
    //
    //_____________________________________________________________________________
    //

/* */ //

// /* */ //

// /* */ //

// /* */ //

// /* */ //

// /* */ //

// /* */ //

//_____________________________________________________________________________ // TGeoTubeSeg - a phi segment of a tube. Has 5 parameters : // - the same 3 as a tube; // - first phi limit (in degrees) // - second phi limit // //_____________________________________________________________________________ // /* */ //

// /* */ //

// /* */ //

// /* */ //

//_____________________________________________________________________________ // TGeoCtub - a tube segment cut with 2 planes. Has 11 parameters : // - the same 5 as a tube segment; // - x, y, z components of the normal to the -dZ cut plane in // point (0, 0, -dZ); // - x, y, z components of the normal to the +dZ cut plane in // point (0, 0, dZ); // //_____________________________________________________________________________ // /* */ //

#include "Riostream.h" #include "TROOT.h" #include "TGeoManager.h" #include "TGeoVolume.h" #include "TVirtualGeoPainter.h" #include "TGeoTube.h" #include "TVirtualPad.h" #include "TBuffer3D.h" #include "TBuffer3DTypes.h" ClassImp(TGeoTube) //_____________________________________________________________________________ TGeoTube::TGeoTube() { // Default constructor SetShapeBit(TGeoShape::kGeoTube); fRmin = 0.0; fRmax = 0.0; fDz = 0.0; } //_____________________________________________________________________________ TGeoTube::TGeoTube(Double_t rmin, Double_t rmax, Double_t dz) :TGeoBBox(0, 0, 0) { // Default constructor specifying minimum and maximum radius SetShapeBit(TGeoShape::kGeoTube); SetTubeDimensions(rmin, rmax, dz); if ((fDz<0) || (fRmin<0) || (fRmax<0)) { SetShapeBit(kGeoRunTimeShape); // if (fRmax<=fRmin) SetShapeBit(kGeoInvalidShape); // printf("tube : dz=%f rmin=%f rmax=%f\n", dz, rmin, rmax); } ComputeBBox(); } //_____________________________________________________________________________ TGeoTube::TGeoTube(const char *name, Double_t rmin, Double_t rmax, Double_t dz) :TGeoBBox(name, 0, 0, 0) { // Default constructor specifying minimum and maximum radius SetShapeBit(TGeoShape::kGeoTube); SetTubeDimensions(rmin, rmax, dz); if ((fDz<0) || (fRmin<0) || (fRmax<0)) { SetShapeBit(kGeoRunTimeShape); // if (fRmax<=fRmin) SetShapeBit(kGeoInvalidShape); // printf("tube : dz=%f rmin=%f rmax=%f\n", dz, rmin, rmax); } ComputeBBox(); } //_____________________________________________________________________________ TGeoTube::TGeoTube(Double_t *param) :TGeoBBox(0, 0, 0) { // Default constructor specifying minimum and maximum radius // param[0] = Rmin // param[1] = Rmax // param[2] = dz SetShapeBit(TGeoShape::kGeoTube); SetDimensions(param); if ((fDz<0) || (fRmin<0) || (fRmax<0)) SetShapeBit(kGeoRunTimeShape); ComputeBBox(); } //_____________________________________________________________________________ TGeoTube::~TGeoTube() { // destructor } //_____________________________________________________________________________ void TGeoTube::ComputeBBox() { // compute bounding box of the tube fDX = fDY = fRmax; fDZ = fDz; } //_____________________________________________________________________________ void TGeoTube::ComputeNormal(Double_t *point, Double_t *dir, Double_t *norm) { // Compute normal to closest surface from POINT. Double_t saf[3]; Double_t rsq = point[0]*point[0]+point[1]*point[1]; Double_t r = TMath::Sqrt(rsq); saf[0] = TMath::Abs(fDz-TMath::Abs(point[2])); saf[1] = (fRmin>1E-10)?TMath::Abs(r-fRmin):TGeoShape::Big(); saf[2] = TMath::Abs(fRmax-r); Int_t i = TMath::LocMin(3,saf); if (i==0) { norm[0] = norm[1] = 0.; norm[2] = TMath::Sign(1.,dir[2]); return; } norm[2] = 0; Double_t phi = TMath::ATan2(point[1], point[0]); norm[0] = TMath::Cos(phi); norm[1] = TMath::Sin(phi); if (norm[0]*dir[0]+norm[1]*dir[1]<0) { norm[0] = -norm[0]; norm[1] = -norm[1]; } } //_____________________________________________________________________________ void TGeoTube::ComputeNormalS(Double_t *point, Double_t *dir, Double_t *norm, Double_t /*rmin*/, Double_t /*rmax*/, Double_t /*dz*/) { // Compute normal to closest surface from POINT. norm[2] = 0; Double_t phi = TMath::ATan2(point[1], point[0]); norm[0] = TMath::Cos(phi); norm[1] = TMath::Sin(phi); if (norm[0]*dir[0]+norm[1]*dir[1]<0) { norm[0] = -norm[0]; norm[1] = -norm[1]; } } //_____________________________________________________________________________ Bool_t TGeoTube::Contains(Double_t *point) const { // test if point is inside this tube if (TMath::Abs(point[2]) > fDz) return kFALSE; Double_t r2 = point[0]*point[0]+point[1]*point[1]; if ((r2<fRmin*fRmin) || (r2>fRmax*fRmax)) return kFALSE; return kTRUE; } //_____________________________________________________________________________ Int_t TGeoTube::DistancetoPrimitive(Int_t px, Int_t py) { // compute closest distance from point px,py to each corner Int_t n = gGeoManager->GetNsegments(); Int_t numPoints = 4*n; if (!HasRmin()) numPoints = 2*(n+1); return ShapeDistancetoPrimitive(numPoints, px, py); } //_____________________________________________________________________________ Double_t TGeoTube::DistFromInsideS(Double_t *point, Double_t *dir, Double_t rmin, Double_t rmax, Double_t dz) { // Compute distance from inside point to surface of the tube (static) // Boundary safe algorithm. // compute distance to surface // Do Z Double_t sz = TGeoShape::Big(); if (dir[2]) { sz = (TMath::Sign(dz, dir[2])-point[2])/dir[2]; if (sz<=0) return 0.0; } // Do R Double_t nsq=dir[0]*dir[0]+dir[1]*dir[1]; if (TMath::Abs(nsq)<TGeoShape::Tolerance()) return sz; Double_t rsq=point[0]*point[0]+point[1]*point[1]; Double_t rdotn=point[0]*dir[0]+point[1]*dir[1]; Double_t b,d; Double_t sr = TGeoShape::Big(); // inner cylinder if (rmin>0) { // Protection in case point is actually outside the tube if (rsq <= rmin*rmin+TGeoShape::Tolerance()) { if (rdotn<0) return 0.0; } else { if (rdotn<0) { DistToTube(rsq,nsq,rdotn,rmin,b,d); if (d>0) { sr=-b-d; if (sr>0) return TMath::Min(sz,sr); } } } } // outer cylinder if (rsq >= rmax*rmax-TGeoShape::Tolerance()) { if (rdotn>=0) return 0.0; } DistToTube(rsq,nsq,rdotn,rmax,b,d); if (d>0) { sr=-b+d; if (sr>0) return TMath::Min(sz,sr); } return 0.; } //_____________________________________________________________________________ Double_t TGeoTube::DistFromInside(Double_t *point, Double_t *dir, Int_t iact, Double_t step, Double_t *safe) const { // Compute distance from inside point to surface of the tube // Boundary safe algorithm. if (iact<3 && safe) { *safe = Safety(point, kTRUE); if (iact==0) return TGeoShape::Big(); if ((iact==1) && (*safe>step)) return TGeoShape::Big(); } // compute distance to surface return DistFromInsideS(point, dir, fRmin, fRmax, fDz); } //_____________________________________________________________________________ Double_t TGeoTube::DistFromOutsideS(Double_t *point, Double_t *dir, Double_t rmin, Double_t rmax, Double_t dz) { // Static method to compute distance from outside point to a tube with given parameters // Boundary safe algorithm. // check Z planes Double_t xi,yi,zi; Double_t rmaxsq = rmax*rmax; Double_t rminsq = rmin*rmin; zi = dz - TMath::Abs(point[2]); Double_t s = TGeoShape::Big(); Bool_t in = kFALSE; Bool_t inz = (zi<0)?kFALSE:kTRUE; if (!inz) { if (point[2]*dir[2]>=0) return TGeoShape::Big(); s = -zi/TMath::Abs(dir[2]); xi = point[0]+s*dir[0]; yi = point[1]+s*dir[1]; Double_t r2=xi*xi+yi*yi; if ((rminsq<=r2) && (r2<=rmaxsq)) return s; } Double_t rsq = point[0]*point[0]+point[1]*point[1]; // check outer cyl. surface Double_t nsq=dir[0]*dir[0]+dir[1]*dir[1]; Double_t rdotn=point[0]*dir[0]+point[1]*dir[1]; Double_t b,d; Bool_t inrmax = kFALSE; Bool_t inrmin = kFALSE; if (rsq<=rmaxsq+TGeoShape::Tolerance()) inrmax = kTRUE; if (rsq>=rminsq-TGeoShape::Tolerance()) inrmin = kTRUE; in = inz & inrmin & inrmax; // If inside, we are most likely on a boundary within machine precision. if (in) { Bool_t checkout = kFALSE; Double_t r = TMath::Sqrt(rsq); if (zi<rmax-r) { if ((rmin==0) || (zi<r-rmin)) { if (point[2]*dir[2]<0) return 0.0; return TGeoShape::Big(); } } if ((rmaxsq-rsq) < (rsq-rminsq)) checkout = kTRUE; if (checkout) { if (rdotn>=0) return TGeoShape::Big(); return 0.0; } if (rmin==0) return 0.0; if (rdotn>=0) return 0.0; // Ray exiting rmin -> check (+) solution for inner tube if (TMath::Abs(nsq)<TGeoShape::Tolerance()) return TGeoShape::Big(); DistToTube(rsq, nsq, rdotn, rmin, b, d); if (d>0) { s=-b+d; if (s>0) { zi=point[2]+s*dir[2]; if (TMath::Abs(zi)<=dz) return s; } } return TGeoShape::Big(); } // Check outer cylinder (only r>rmax has to be considered) if (TMath::Abs(nsq)<TGeoShape::Tolerance()) return TGeoShape::Big(); if (!inrmax) { DistToTube(rsq, nsq, rdotn, rmax, b, d); if (d>0) { s=-b-d; if (s>0) { zi=point[2]+s*dir[2]; if (TMath::Abs(zi)<=dz) return s; } } } // check inner cylinder if (rmin>0) { DistToTube(rsq, nsq, rdotn, rmin, b, d); if (d>0) { s=-b+d; if (s>0) { zi=point[2]+s*dir[2]; if (TMath::Abs(zi)<=dz) return s; } } } return TGeoShape::Big(); } //_____________________________________________________________________________ Double_t TGeoTube::DistFromOutside(Double_t *point, Double_t *dir, Int_t iact, Double_t step, Double_t *safe) const { // Compute distance from outside point to surface of the tube and safe distance // Boundary safe algorithm. // fist localize point w.r.t tube if (iact<3 && safe) { *safe = Safety(point, kFALSE); if (iact==0) return TGeoShape::Big(); if ((iact==1) && (step<=*safe)) return TGeoShape::Big(); } // find distance to shape return DistFromOutsideS(point, dir, fRmin, fRmax, fDz); } //_____________________________________________________________________________ void TGeoTube::DistToTube(Double_t rsq, Double_t nsq, Double_t rdotn, Double_t radius, Double_t &b, Double_t &delta) { // Static method computing the distance to a tube with given radius, starting from // POINT along DIR director cosines. The distance is computed as : // RSQ = point[0]*point[0]+point[1]*point[1] // NSQ = dir[0]*dir[0]+dir[1]*dir[1] ---> should NOT be 0 !!! // RDOTN = point[0]*dir[0]+point[1]*dir[1] // The distance can be computed as : // D = -B +/- DELTA // where DELTA.GT.0 and D.GT.0 Double_t t1 = 1./nsq; Double_t t3=rsq-(radius*radius); b = t1*rdotn; Double_t c =t1*t3; delta = b*b-c; if (delta>0) { delta=TMath::Sqrt(delta); } else { delta = -1; } } //_____________________________________________________________________________ TGeoVolume *TGeoTube::Divide(TGeoVolume *voldiv, const char *divname, Int_t iaxis, Int_t ndiv, Double_t start, Double_t step) { //--- Divide this tube shape belonging to volume "voldiv" into ndiv volumes // called divname, from start position with the given step. Returns pointer // to created division cell volume in case of Z divisions. For radial division // creates all volumes with different shapes and returns pointer to volume that // was divided. In case a wrong division axis is supplied, returns pointer to // volume that was divided. TGeoShape *shape; //--- shape to be created TGeoVolume *vol; //--- division volume to be created TGeoVolumeMulti *vmulti; //--- generic divided volume TGeoPatternFinder *finder; //--- finder to be attached TString opt = ""; //--- option to be attached Int_t id; Double_t end = start+ndiv*step; switch (iaxis) { case 1: //--- R division finder = new TGeoPatternCylR(voldiv, ndiv, start, end); vmulti = gGeoManager->MakeVolumeMulti(divname, voldiv->GetMedium()); voldiv->SetFinder(finder); finder->SetDivIndex(voldiv->GetNdaughters()); for (id=0; id<ndiv; id++) { shape = new TGeoTube(start+id*step, start+(id+1)*step, fDz); vol = new TGeoVolume(divname, shape, voldiv->GetMedium()); vmulti->AddVolume(vol); opt = "R"; voldiv->AddNodeOffset(vol, id, 0, opt.Data()); ((TGeoNodeOffset*)voldiv->GetNodes()->At(voldiv->GetNdaughters()-1))->SetFinder(finder); } return vmulti; case 2: //--- Phi division finder = new TGeoPatternCylPhi(voldiv, ndiv, start, end); voldiv->SetFinder(finder); finder->SetDivIndex(voldiv->GetNdaughters()); shape = new TGeoTubeSeg(fRmin, fRmax, fDz, -step/2, step/2); vol = new TGeoVolume(divname, shape, voldiv->GetMedium()); vmulti = gGeoManager->MakeVolumeMulti(divname, voldiv->GetMedium()); vmulti->AddVolume(vol); opt = "Phi"; for (id=0; id<ndiv; id++) { voldiv->AddNodeOffset(vol, id, start+id*step+step/2, opt.Data()); ((TGeoNodeOffset*)voldiv->GetNodes()->At(voldiv->GetNdaughters()-1))->SetFinder(finder); } return vmulti; case 3: //--- Z division finder = new TGeoPatternZ(voldiv, ndiv, start, start+ndiv*step); voldiv->SetFinder(finder); finder->SetDivIndex(voldiv->GetNdaughters()); shape = new TGeoTube(fRmin, fRmax, step/2); vol = new TGeoVolume(divname, shape, voldiv->GetMedium()); vmulti = gGeoManager->MakeVolumeMulti(divname, voldiv->GetMedium()); vmulti->AddVolume(vol); opt = "Z"; for (id=0; id<ndiv; id++) { voldiv->AddNodeOffset(vol, id, start+step/2+id*step, opt.Data()); ((TGeoNodeOffset*)voldiv->GetNodes()->At(voldiv->GetNdaughters()-1))->SetFinder(finder); } return vmulti; default: Error("Divide", "In shape %s wrong axis type for division", GetName()); return 0; } } //_____________________________________________________________________________ const char *TGeoTube::GetAxisName(Int_t iaxis) const { // Returns name of axis IAXIS. switch (iaxis) { case 1: return "R"; case 2: return "PHI"; case 3: return "Z"; default: return "UNDEFINED"; } } //_____________________________________________________________________________ Double_t TGeoTube::GetAxisRange(Int_t iaxis, Double_t &xlo, Double_t &xhi) const { // Get range of shape for a given axis. xlo = 0; xhi = 0; Double_t dx = 0; switch (iaxis) { case 1: xlo = fRmin; xhi = fRmax; dx = xhi-xlo; return dx; case 2: xlo = 0; xhi = 360; dx = 360; return dx; case 3: xlo = -fDz; xhi = fDz; dx = xhi-xlo; return dx; } return dx; } //_____________________________________________________________________________ void TGeoTube::GetBoundingCylinder(Double_t *param) const { //--- Fill vector param[4] with the bounding cylinder parameters. The order // is the following : Rmin, Rmax, Phi1, Phi2, dZ param[0] = fRmin; // Rmin param[0] *= param[0]; param[1] = fRmax; // Rmax param[1] *= param[1]; param[2] = 0.; // Phi1 param[3] = 360.; // Phi1 } //_____________________________________________________________________________ TGeoShape *TGeoTube::GetMakeRuntimeShape(TGeoShape *mother, TGeoMatrix * /*mat*/) const { // in case shape has some negative parameters, these has to be computed // in order to fit the mother if (!TestShapeBit(kGeoRunTimeShape)) return 0; Double_t rmin, rmax, dz; Double_t xmin,xmax; rmin = fRmin; rmax = fRmax; dz = fDz; if (fDz<0) { mother->GetAxisRange(3,xmin,xmax); if (xmax<0) return 0; dz=xmax; } mother->GetAxisRange(1,xmin,xmax); if (fRmin<0) { if (xmin<0) return 0; rmin = xmin; } if (fRmax<0) { if (xmax<=0) return 0; rmax = xmax; } return (new TGeoTube(GetName(), rmin, rmax, dz)); } //_____________________________________________________________________________ void TGeoTube::InspectShape() const { // print shape parameters printf("*** Shape %s: TGeoTube ***\n", GetName()); printf(" Rmin = %11.5f\n", fRmin); printf(" Rmax = %11.5f\n", fRmax); printf(" dz = %11.5f\n", fDz); printf(" Bounding box:\n"); TGeoBBox::InspectShape(); } //_____________________________________________________________________________ TBuffer3D *TGeoTube::MakeBuffer3D() const { // Creates a TBuffer3D describing *this* shape. // Coordinates are in local reference frame. Int_t n = gGeoManager->GetNsegments(); Int_t NbPnts = 4*n; Int_t NbSegs = 8*n; Int_t NbPols = 4*n; TBuffer3D* buff = new TBuffer3D(TBuffer3DTypes::kGeneric, NbPnts, 3*NbPnts, NbSegs, 3*NbSegs, NbPols, 6*NbPols); if (buff) { SetPoints(buff->fPnts); SetSegsAndPols(*buff); } return buff; } //_____________________________________________________________________________ void TGeoTube::SetSegsAndPols(TBuffer3D &buffer) const { // Fill TBuffer3D structure for segments and polygons. Int_t i, j,indx; Int_t n = gGeoManager->GetNsegments(); Int_t c = (((buffer.fColor) %8) -1) * 4; if (c < 0) c = 0; if (HasRmin()) { // circle segments: // lower rmin circle: i=0, (0, n-1) // lower rmax circle: i=1, (n, 2n-1) // upper rmin circle: i=2, (2n, 3n-1) // upper rmax circle: i=1, (3n, 4n-1) for (i = 0; i < 4; i++) { for (j = 0; j < n; j++) { indx = 3*(i*n+j); buffer.fSegs[indx ] = c; buffer.fSegs[indx+1] = i*n+j; buffer.fSegs[indx+2] = i*n+(j+1)%n; } } // Z-parallel segments // inner: i=4, (4n, 5n-1) // outer: i=5, (5n, 6n-1) for (i = 4; i < 6; i++) { for (j = 0; j < n; j++) { indx = 3*(i*n+j); buffer.fSegs[indx ] = c+1; buffer.fSegs[indx+1] = (i-4)*n+j; buffer.fSegs[indx+2] = (i-2)*n+j; } } // Radial segments // lower: i=6, (6n, 7n-1) // upper: i=7, (7n, 8n-1) for (i = 6; i < 8; i++) { for (j = 0; j < n; j++) { indx = 3*(i*n+j); buffer.fSegs[indx ] = c; buffer.fSegs[indx+1] = 2*(i-6)*n+j; buffer.fSegs[indx+2] = (2*(i-6)+1)*n+j; } } // Polygons i=0; // Inner lateral (0, n-1) for (j = 0; j < n; j++) { indx = 6*(i*n+j); buffer.fPols[indx ] = c; buffer.fPols[indx+1] = 4; buffer.fPols[indx+2] = j; buffer.fPols[indx+3] = 4*n+(j+1)%n; buffer.fPols[indx+4] = 2*n+j; buffer.fPols[indx+5] = 4*n+j; } i=1; // Outer lateral (n,2n-1) for (j = 0; j < n; j++) { indx = 6*(i*n+j); buffer.fPols[indx ] = c+1; buffer.fPols[indx+1] = 4; buffer.fPols[indx+2] = n+j; buffer.fPols[indx+3] = 5*n+j; buffer.fPols[indx+4] = 3*n+j; buffer.fPols[indx+5] = 5*n+(j+1)%n; } i=2; // lower disc (2n, 3n-1) for (j = 0; j < n; j++) { indx = 6*(i*n+j); buffer.fPols[indx ] = c; buffer.fPols[indx+1] = 4; buffer.fPols[indx+2] = j; buffer.fPols[indx+3] = 6*n+j; buffer.fPols[indx+4] = n+j; buffer.fPols[indx+5] = 6*n+(j+1)%n; } i=3; // upper disc (3n, 4n-1) for (j = 0; j < n; j++) { indx = 6*(i*n+j); buffer.fPols[indx ] = c; buffer.fPols[indx+1] = 4; buffer.fPols[indx+2] = 2*n+j; buffer.fPols[indx+3] = 7*n+(j+1)%n; buffer.fPols[indx+4] = 3*n+j; buffer.fPols[indx+5] = 7*n+j; } return; } // Rmin=0 tubes // circle segments // lower rmax circle: i=0, (0, n-1) // upper rmax circle: i=1, (n, 2n-1) for (i = 0; i < 2; i++) { for (j = 0; j < n; j++) { indx = 3*(i*n+j); buffer.fSegs[indx ] = c; buffer.fSegs[indx+1] = 2+i*n+j; buffer.fSegs[indx+2] = 2+i*n+(j+1)%n; } } // Z-parallel segments (2n,3n-1) for (j = 0; j < n; j++) { indx = 3*(2*n+j); buffer.fSegs[indx ] = c+1; buffer.fSegs[indx+1] = 2+j; buffer.fSegs[indx+2] = 2+n+j; } // Radial segments // Lower circle: i=3, (3n,4n-1) // Upper circle: i=4, (4n,5n-1) for (i = 3; i < 5; i++) { for (j = 0; j < n; j++) { indx = 3*(i*n+j); buffer.fSegs[indx ] = c; buffer.fSegs[indx+1] = i-3; buffer.fSegs[indx+2] = 2+(i-3)*n+j; } } // Polygons // lateral (0,n-1) for (j = 0; j < n; j++) { indx = 6*j; buffer.fPols[indx ] = c+1; buffer.fPols[indx+1] = 4; buffer.fPols[indx+2] = j; buffer.fPols[indx+3] = 2*n+j; buffer.fPols[indx+4] = n+j; buffer.fPols[indx+5] = 2*n+(j+1)%n; } // bottom triangles (n,2n-1) for (j = 0; j < n; j++) { indx = 6*n + 5*j; buffer.fPols[indx ] = c; buffer.fPols[indx+1] = 3; buffer.fPols[indx+2] = j; buffer.fPols[indx+3] = 3*n+(j+1)%n; buffer.fPols[indx+4] = 3*n+j; } // top triangles (2n,3n-1) for (j = 0; j < n; j++) { indx = 6*n + 5*n + 5*j; buffer.fPols[indx ] = c; buffer.fPols[indx+1] = 3; buffer.fPols[indx+2] = n+j; buffer.fPols[indx+3] = 4*n+j; buffer.fPols[indx+4] = 4*n+(j+1)%n; } } //_____________________________________________________________________________ Double_t TGeoTube::Safety(Double_t *point, Bool_t in) const { // computes the closest distance from given point to this shape, according // to option. The matching point on the shape is stored in spoint. #ifndef NEVER Double_t r = TMath::Sqrt(point[0]*point[0]+point[1]*point[1]); Double_t safe, safrmin, safrmax; if (in) { safe = fDz-TMath::Abs(point[2]); // positive if inside if (fRmin>1E-10) { safrmin = r-fRmin; if (safrmin < safe) safe = safrmin; } safrmax = fRmax-r; if (safrmax < safe) safe = safrmax; } else { safe = -fDz+TMath::Abs(point[2]); if (fRmin>1E-10) { safrmin = -r+fRmin; if (safrmin > safe) safe = safrmin; } safrmax = -fRmax+r; if (safrmax > safe) safe = safrmax; } return safe; #else Double_t saf[3]; Double_t rsq = point[0]*point[0]+point[1]*point[1]; Double_t r = TMath::Sqrt(rsq); saf[0] = fDz-TMath::Abs(point[2]); // positive if inside saf[1] = (fRmin>1E-10)?(r-fRmin):TGeoShape::Big(); saf[2] = fRmax-r; if (in) return saf[TMath::LocMin(3,saf)]; for (Int_t i=0; i<3; i++) saf[i]=-saf[i]; return saf[TMath::LocMax(3,saf)]; #endif } //_____________________________________________________________________________ Double_t TGeoTube::SafetyS(Double_t *point, Bool_t in, Double_t rmin, Double_t rmax, Double_t dz, Int_t skipz) { // computes the closest distance from given point to this shape, according // to option. The matching point on the shape is stored in spoint. Double_t saf[3]; Double_t rsq = point[0]*point[0]+point[1]*point[1]; Double_t r = TMath::Sqrt(rsq); switch (skipz) { case 1: // skip lower Z plane saf[0] = dz - point[2]; break; case 2: // skip upper Z plane saf[0] = dz + point[2]; break; case 3: // skip both saf[0] = TGeoShape::Big(); break; default: saf[0] = dz-TMath::Abs(point[2]); } saf[1] = (rmin>1E-10)?(r-rmin):TGeoShape::Big(); saf[2] = rmax-r; // printf("saf0=%g saf1=%g saf2=%g in=%d skipz=%d\n", saf[0],saf[1],saf[2],in,skipz); if (in) return saf[TMath::LocMin(3,saf)]; for (Int_t i=0; i<3; i++) saf[i]=-saf[i]; return saf[TMath::LocMax(3,saf)]; } //_____________________________________________________________________________ void TGeoTube::SavePrimitive(ofstream &out, Option_t * /*option*/) { // Save a primitive as a C++ statement(s) on output stream "out". if (TObject::TestBit(kGeoSavePrimitive)) return; out << " // Shape: " << GetName() << " type: " << ClassName() << endl; out << " rmin = " << fRmin << ";" << endl; out << " rmax = " << fRmax << ";" << endl; out << " dz = " << fDz << ";" << endl; out << " TGeoShape *" << GetPointerName() << " = new TGeoTube(\"" << GetName() << "\",rmin,rmax,dz);" << endl; TObject::SetBit(TGeoShape::kGeoSavePrimitive); } //_____________________________________________________________________________ void TGeoTube::SetTubeDimensions(Double_t rmin, Double_t rmax, Double_t dz) { fRmin = rmin; fRmax = rmax; fDz = dz; if (fRmin>0 && fRmax>0 && fRmin>=fRmax) Error("SetTubeDimensions", "In shape %s wrong rmin=%g rmax=%g", GetName(), rmin,rmax); } //_____________________________________________________________________________ void TGeoTube::SetDimensions(Double_t *param) { Double_t rmin = param[0]; Double_t rmax = param[1]; Double_t dz = param[2]; SetTubeDimensions(rmin, rmax, dz); } //_____________________________________________________________________________ void TGeoTube::SetPoints(Double_t *points) const { // create tube mesh points Double_t dz; Int_t j, n; n = gGeoManager->GetNsegments(); Double_t dphi = 360./n; Double_t phi = 0; dz = fDz; Int_t indx = 0; if (points) { if (HasRmin()) { // 4*n points // (0,n-1) lower rmin circle // (2n, 3n-1) upper rmin circle for (j = 0; j < n; j++) { phi = j*dphi*TMath::DegToRad(); points[indx+6*n] = points[indx] = fRmin * TMath::Cos(phi); indx++; points[indx+6*n] = points[indx] = fRmin * TMath::Sin(phi); indx++; points[indx+6*n] = dz; points[indx] =-dz; indx++; } // (n, 2n-1) lower rmax circle // (3n, 4n-1) upper rmax circle for (j = 0; j < n; j++) { phi = j*dphi*TMath::DegToRad(); points[indx+6*n] = points[indx] = fRmax * TMath::Cos(phi); indx++; points[indx+6*n] = points[indx] = fRmax * TMath::Sin(phi); indx++; points[indx+6*n]= dz; points[indx] =-dz; indx++; } } else { // centers of lower/upper circles (0,1) points[indx++] = 0.; points[indx++] = 0.; points[indx++] = -dz; points[indx++] = 0.; points[indx++] = 0.; points[indx++] = dz; // lower rmax circle (2, 2+n-1) // upper rmax circle (2+n, 2+2n-1) for (j = 0; j < n; j++) { phi = j*dphi*TMath::DegToRad(); points[indx+3*n] = points[indx] = fRmax * TMath::Cos(phi); indx++; points[indx+3*n] = points[indx] = fRmax * TMath::Sin(phi); indx++; points[indx+3*n]= dz; points[indx] =-dz; indx++; } } } } //_____________________________________________________________________________ void TGeoTube::SetPoints(Float_t *points) const { // create tube mesh points Double_t dz; Int_t j, n; n = gGeoManager->GetNsegments(); Double_t dphi = 360./n; Double_t phi = 0; dz = fDz; Int_t indx = 0; if (points) { if (HasRmin()) { // 4*n points // (0,n-1) lower rmin circle // (2n, 3n-1) upper rmin circle for (j = 0; j < n; j++) { phi = j*dphi*TMath::DegToRad(); points[indx+6*n] = points[indx] = fRmin * TMath::Cos(phi); indx++; points[indx+6*n] = points[indx] = fRmin * TMath::Sin(phi); indx++; points[indx+6*n] = dz; points[indx] =-dz; indx++; } // (n, 2n-1) lower rmax circle // (3n, 4n-1) upper rmax circle for (j = 0; j < n; j++) { phi = j*dphi*TMath::DegToRad(); points[indx+6*n] = points[indx] = fRmax * TMath::Cos(phi); indx++; points[indx+6*n] = points[indx] = fRmax * TMath::Sin(phi); indx++; points[indx+6*n]= dz; points[indx] =-dz; indx++; } } else { // centers of lower/upper circles (0,1) points[indx++] = 0.; points[indx++] = 0.; points[indx++] = -dz; points[indx++] = 0.; points[indx++] = 0.; points[indx++] = dz; // lower rmax circle (2, 2+n-1) // upper rmax circle (2+n, 2+2n-1) for (j = 0; j < n; j++) { phi = j*dphi*TMath::DegToRad(); points[indx+3*n] = points[indx] = fRmax * TMath::Cos(phi); indx++; points[indx+3*n] = points[indx] = fRmax * TMath::Sin(phi); indx++; points[indx+3*n]= dz; points[indx] =-dz; indx++; } } } } //_____________________________________________________________________________ Int_t TGeoTube::GetNmeshVertices() const { // Return number of vertices of the mesh representation Int_t n = gGeoManager->GetNsegments(); Int_t numPoints = n*4; if (!HasRmin()) numPoints = 2*(n+1); return numPoints; } //_____________________________________________________________________________ void TGeoTube::Sizeof3D() const { ///// fill size of this 3-D object /// TVirtualGeoPainter *painter = gGeoManager->GetGeomPainter(); /// if (!painter) return; /// Int_t n = gGeoManager->GetNsegments(); /// Int_t numPoints = n*4; /// Int_t numSegs = n*8; /// Int_t numPolys = n*4; /// painter->AddSize3D(numPoints, numSegs, numPolys); } //_____________________________________________________________________________ const TBuffer3D & TGeoTube::GetBuffer3D(Int_t reqSections, Bool_t localFrame) const { static TBuffer3DTube buffer; TGeoBBox::FillBuffer3D(buffer, reqSections, localFrame); if (reqSections & TBuffer3D::kShapeSpecific) { buffer.fRadiusInner = fRmin; buffer.fRadiusOuter = fRmax; buffer.fHalfLength = fDz; buffer.SetSectionsValid(TBuffer3D::kShapeSpecific); } if (reqSections & TBuffer3D::kRawSizes) { Int_t n = gGeoManager->GetNsegments(); Int_t NbPnts = 4*n; Int_t NbSegs = 8*n; Int_t NbPols = 4*n; if (!HasRmin()) { NbPnts = 2*(n+1); NbSegs = 5*n; NbPols = 3*n; } if (buffer.SetRawSizes(NbPnts, 3*NbPnts, NbSegs, 3*NbSegs, NbPols, 6*NbPols)) { buffer.SetSectionsValid(TBuffer3D::kRawSizes); } } if ((reqSections & TBuffer3D::kRaw) && buffer.SectionsValid(TBuffer3D::kRawSizes)) { SetPoints(buffer.fPnts); if (!buffer.fLocalFrame) { TransformPoints(buffer.fPnts, buffer.NbPnts()); } SetSegsAndPols(buffer); buffer.SetSectionsValid(TBuffer3D::kRaw); } return buffer; } ClassImp(TGeoTubeSeg) //_____________________________________________________________________________ TGeoTubeSeg::TGeoTubeSeg() { // Default constructor SetShapeBit(TGeoShape::kGeoTubeSeg); fPhi1 = fPhi2 = 0.0; } //_____________________________________________________________________________ TGeoTubeSeg::TGeoTubeSeg(Double_t rmin, Double_t rmax, Double_t dz, Double_t phi1, Double_t phi2) :TGeoTube(rmin, rmax, dz) { // Default constructor specifying minimum and maximum radius SetShapeBit(TGeoShape::kGeoTubeSeg); SetTubsDimensions(rmin, rmax, dz, phi1, phi2); ComputeBBox(); } //_____________________________________________________________________________ TGeoTubeSeg::TGeoTubeSeg(const char *name, Double_t rmin, Double_t rmax, Double_t dz, Double_t phi1, Double_t phi2) :TGeoTube(name, rmin, rmax, dz) { // Default constructor specifying minimum and maximum radius SetShapeBit(TGeoShape::kGeoTubeSeg); SetTubsDimensions(rmin, rmax, dz, phi1, phi2); ComputeBBox(); } //_____________________________________________________________________________ TGeoTubeSeg::TGeoTubeSeg(Double_t *param) :TGeoTube(0, 0, 0) { // Default constructor specifying minimum and maximum radius // param[0] = Rmin // param[1] = Rmax // param[2] = dz // param[3] = phi1 // param[4] = phi2 SetShapeBit(TGeoShape::kGeoTubeSeg); SetDimensions(param); ComputeBBox(); } //_____________________________________________________________________________ TGeoTubeSeg::~TGeoTubeSeg() { // destructor } //_____________________________________________________________________________ void TGeoTubeSeg::ComputeBBox() { // compute bounding box of the tube segment Double_t xc[4]; Double_t yc[4]; xc[0] = fRmax*TMath::Cos(fPhi1*TMath::DegToRad()); yc[0] = fRmax*TMath::Sin(fPhi1*TMath::DegToRad()); xc[1] = fRmax*TMath::Cos(fPhi2*TMath::DegToRad()); yc[1] = fRmax*TMath::Sin(fPhi2*TMath::DegToRad()); xc[2] = fRmin*TMath::Cos(fPhi1*TMath::DegToRad()); yc[2] = fRmin*TMath::Sin(fPhi1*TMath::DegToRad()); xc[3] = fRmin*TMath::Cos(fPhi2*TMath::DegToRad()); yc[3] = fRmin*TMath::Sin(fPhi2*TMath::DegToRad()); Double_t xmin = xc[TMath::LocMin(4, &xc[0])]; Double_t xmax = xc[TMath::LocMax(4, &xc[0])]; Double_t ymin = yc[TMath::LocMin(4, &yc[0])]; Double_t ymax = yc[TMath::LocMax(4, &yc[0])]; Double_t dp = fPhi2-fPhi1; if (dp<0) dp+=360; Double_t ddp = -fPhi1; if (ddp<0) ddp+= 360; if (ddp>360) ddp-=360; if (ddp<=dp) xmax = fRmax; ddp = 90-fPhi1; if (ddp<0) ddp+= 360; if (ddp>360) ddp-=360; if (ddp<=dp) ymax = fRmax; ddp = 180-fPhi1; if (ddp<0) ddp+= 360; if (ddp>360) ddp-=360; if (ddp<=dp) xmin = -fRmax; ddp = 270-fPhi1; if (ddp<0) ddp+= 360; if (ddp>360) ddp-=360; if (ddp<=dp) ymin = -fRmax; fOrigin[0] = (xmax+xmin)/2; fOrigin[1] = (ymax+ymin)/2; fOrigin[2] = 0; fDX = (xmax-xmin)/2; fDY = (ymax-ymin)/2; fDZ = fDz; } //_____________________________________________________________________________ void TGeoTubeSeg::ComputeNormal(Double_t *point, Double_t *dir, Double_t *norm) { // Compute normal to closest surface from POINT. Double_t saf[3]; Double_t rsq = point[0]*point[0]+point[1]*point[1]; Double_t r = TMath::Sqrt(rsq); Double_t c1 = TMath::Cos(fPhi1*TMath::DegToRad()); Double_t s1 = TMath::Sin(fPhi1*TMath::DegToRad()); Double_t c2 = TMath::Cos(fPhi2*TMath::DegToRad()); Double_t s2 = TMath::Sin(fPhi2*TMath::DegToRad()); saf[0] = TMath::Abs(fDz-TMath::Abs(point[2])); saf[1] = (fRmin>1E-10)?TMath::Abs(r-fRmin):TGeoShape::Big(); saf[2] = TMath::Abs(fRmax-r); Int_t i = TMath::LocMin(3,saf); if (TGeoShape::IsCloseToPhi(saf[i], point,c1,s1,c2,s2)) { TGeoShape::NormalPhi(point,dir,norm,c1,s1,c2,s2); return; } if (i==0) { norm[0] = norm[1] = 0.; norm[2] = TMath::Sign(1.,dir[2]); return; } norm[2] = 0; Double_t phi = TMath::ATan2(point[1], point[0]); norm[0] = TMath::Cos(phi); norm[1] = TMath::Sin(phi); if (norm[0]*dir[0]+norm[1]*dir[1]<0) { norm[0] = -norm[0]; norm[1] = -norm[1]; } } //_____________________________________________________________________________ void TGeoTubeSeg::ComputeNormalS(Double_t *point, Double_t *dir, Double_t *norm, Double_t rmin, Double_t rmax, Double_t /*dz*/, Double_t c1, Double_t s1, Double_t c2, Double_t s2) { // Compute normal to closest surface from POINT. Double_t saf[2]; Double_t rsq = point[0]*point[0]+point[1]*point[1]; Double_t r = TMath::Sqrt(rsq); saf[0] = (rmin>1E-10)?TMath::Abs(r-rmin):TGeoShape::Big(); saf[1] = TMath::Abs(rmax-r); Int_t i = TMath::LocMin(2,saf); if (TGeoShape::IsCloseToPhi(saf[i], point,c1,s1,c2,s2)) { TGeoShape::NormalPhi(point,dir,norm,c1,s1,c2,s2); return; } norm[2] = 0; Double_t phi = TMath::ATan2(point[1], point[0]); norm[0] = TMath::Cos(phi); norm[1] = TMath::Sin(phi); if (norm[0]*dir[0]+norm[1]*dir[1]<0) { norm[0] = -norm[0]; norm[1] = -norm[1]; } } //_____________________________________________________________________________ Bool_t TGeoTubeSeg::Contains(Double_t *point) const { // test if point is inside this tube segment // first check if point is inside the tube if (!TGeoTube::Contains(point)) return kFALSE; return IsInPhiRange(point, fPhi1, fPhi2); } //_____________________________________________________________________________ Int_t TGeoTubeSeg::DistancetoPrimitive(Int_t px, Int_t py) { // compute closest distance from point px,py to each corner Int_t n = gGeoManager->GetNsegments()+1; const Int_t numPoints = 4*n; return ShapeDistancetoPrimitive(numPoints, px, py); } //_____________________________________________________________________________ Double_t TGeoTubeSeg::DistFromInsideS(Double_t *point, Double_t *dir, Double_t rmin, Double_t rmax, Double_t dz, Double_t c1, Double_t s1, Double_t c2, Double_t s2, Double_t cm, Double_t sm, Double_t cdfi) { // Compute distance from inside point to surface of the tube segment (static) // Boundary safe algorithm. // Do Z Double_t stube = TGeoTube::DistFromInsideS(point,dir,rmin,rmax,dz); if (stube<=0) return 0.0; Double_t sfmin = TGeoShape::Big(); Double_t rsq = point[0]*point[0]+point[1]*point[1]; Double_t r = TMath::Sqrt(rsq); Double_t cpsi=point[0]*cm+point[1]*sm; if (cpsi>r*cdfi+TGeoShape::Tolerance()) { sfmin = TGeoShape::DistToPhiMin(point, dir, s1, c1, s2, c2, sm, cm); return TMath::Min(stube,sfmin); } // Point on the phi boundary or outside // which one: phi1 or phi2 Double_t ddotn, xi, yi; if (TMath::Abs(point[1]-s1*r) < TMath::Abs(point[1]-s2*r)) { ddotn = s1*dir[0]-c1*dir[1]; if (ddotn>=0) return 0.0; ddotn = -s2*dir[0]+c2*dir[1]; if (ddotn<=0) return stube; sfmin = s2*point[0]-c2*point[1]; if (sfmin<=0) return stube; sfmin /= ddotn; if (sfmin >= stube) return stube; xi = point[0]+sfmin*dir[0]; yi = point[1]+sfmin*dir[1]; if (yi*cm-xi*sm<0) return stube; return sfmin; } ddotn = -s2*dir[0]+c2*dir[1]; if (ddotn>=0) return 0.0; ddotn = s1*dir[0]-c1*dir[1]; if (ddotn<=0) return stube; sfmin = -s1*point[0]+c1*point[1]; if (sfmin<=0) return stube; sfmin /= ddotn; if (sfmin >= stube) return stube; xi = point[0]+sfmin*dir[0]; yi = point[1]+sfmin*dir[1]; if (yi*cm-xi*sm>0) return stube; return sfmin; } //_____________________________________________________________________________ Double_t TGeoTubeSeg::DistFromInside(Double_t *point, Double_t *dir, Int_t iact, Double_t step, Double_t *safe) const { // Compute distance from inside point to surface of the tube segment // Boundary safe algorithm. if (iact<3 && safe) { *safe = SafetyS(point, kTRUE, fRmin, fRmax, fDz, fPhi1, fPhi2); if (iact==0) return TGeoShape::Big(); if ((iact==1) && (*safe>step)) return TGeoShape::Big(); } Double_t phi1 = fPhi1*TMath::DegToRad(); Double_t phi2 = fPhi2*TMath::DegToRad(); Double_t c1 = TMath::Cos(phi1); Double_t c2 = TMath::Cos(phi2); Double_t s1 = TMath::Sin(phi1); Double_t s2 = TMath::Sin(phi2); Double_t phim = 0.5*(phi1+phi2); Double_t cm = TMath::Cos(phim); Double_t sm = TMath::Sin(phim); Double_t dfi = 0.5*(phi2-phi1); Double_t cdfi = TMath::Cos(dfi); // compute distance to surface return TGeoTubeSeg::DistFromInsideS(point,dir,fRmin,fRmax,fDz,c1,s1,c2,s2,cm,sm,cdfi); } //_____________________________________________________________________________ Double_t TGeoTubeSeg::DistFromOutsideS(Double_t *point, Double_t *dir, Double_t rmin, Double_t rmax, Double_t dz, Double_t c1, Double_t s1, Double_t c2, Double_t s2, Double_t cm, Double_t sm, Double_t cdfi) { // Static method to compute distance to arbitrary tube segment from outside point // Boundary safe algorithm. Double_t r2, cpsi; // check Z planes Double_t xi, yi, zi; zi = dz - TMath::Abs(point[2]); Double_t rmaxsq = rmax*rmax; Double_t rminsq = rmin*rmin; Double_t s = TGeoShape::Big(); Double_t snxt=TGeoShape::Big(); Bool_t in = kFALSE; Bool_t inz = (zi<0)?kFALSE:kTRUE; if (!inz) { if (point[2]*dir[2]>=0) return TGeoShape::Big(); s = -zi/TMath::Abs(dir[2]); xi = point[0]+s*dir[0]; yi = point[1]+s*dir[1]; r2=xi*xi+yi*yi; if ((rminsq<=r2) && (r2<=rmaxsq)) { cpsi=(xi*cm+yi*sm)/TMath::Sqrt(r2); if (cpsi>=cdfi) return s; } } // check outer cyl. surface Double_t rsq = point[0]*point[0]+point[1]*point[1]; Double_t r = TMath::Sqrt(rsq); Double_t nsq=dir[0]*dir[0]+dir[1]*dir[1]; Double_t rdotn=point[0]*dir[0]+point[1]*dir[1]; Double_t b,d; Bool_t inrmax = kFALSE; Bool_t inrmin = kFALSE; Bool_t inphi = kFALSE; if (rsq<=rmaxsq+TGeoShape::Tolerance()) inrmax = kTRUE; if (rsq>=rminsq-TGeoShape::Tolerance()) inrmin = kTRUE; cpsi=point[0]*cm+point[1]*sm; if (cpsi>r*cdfi-TGeoShape::Tolerance()) inphi = kTRUE; in = inz & inrmin & inrmax & inphi; // If inside, we are most likely on a boundary within machine precision. if (in) { Bool_t checkout = kFALSE; Double_t safphi = (cpsi-r*cdfi)*TMath::Sqrt(1.-cdfi*cdfi); // Double_t sch, cch; // check if on Z boundaries if (zi<rmax-r) { if ((rmin==0) || (zi<r-rmin)) { if (zi<safphi) { if (point[2]*dir[2]<0) return 0.0; return TGeoShape::Big(); } } } if ((rmaxsq-rsq) < (rsq-rminsq)) checkout = kTRUE; // check if on Rmax boundary if (checkout && (rmax-r<safphi)) { if (rdotn>=0) return TGeoShape::Big(); return 0.0; } if (TMath::Abs(nsq)<TGeoShape::Tolerance()) return TGeoShape::Big(); // check if on phi boundary if ((rmin==0) || (safphi<r-rmin)) { // We may cross again a phi of rmin boundary // check first if we are on phi1 or phi2 Double_t un; if (TMath::Abs(point[1]-s1*r) < TMath::Abs(point[1]-s2*r)) { un = dir[0]*s1-dir[1]*c1; if (un < 0) return 0.0; if (cdfi>=0) return TGeoShape::Big(); un = -dir[0]*s2+dir[1]*c2; if (un<0) { s = -point[0]*s2+point[1]*c2; if (s>0) { s /= (-un); zi = point[2]+s*dir[2]; if (TMath::Abs(zi)<=dz) { xi = point[0]+s*dir[0]; yi = point[1]+s*dir[1]; r2=xi*xi+yi*yi; if ((rminsq<=r2) && (r2<=rmaxsq)) { if ((yi*cm-xi*sm)>0) return s; } } } } } else { un = -dir[0]*s2+dir[1]*c2; if (un < 0) return 0.0; if (cdfi>=0) return TGeoShape::Big(); un = dir[0]*s1-dir[1]*c1; if (un<0) { s = point[0]*s1-point[1]*c1; if (s>0) { s /= (-un); zi = point[2]+s*dir[2]; if (TMath::Abs(zi)<=dz) { xi = point[0]+s*dir[0]; yi = point[1]+s*dir[1]; r2=xi*xi+yi*yi; if ((rminsq<=r2) && (r2<=rmaxsq)) { if ((yi*cm-xi*sm)<0) return s; } } } } } // We may also cross rmin, (+) solution if (rdotn>=0) return TGeoShape::Big(); if (cdfi>=0) return TGeoShape::Big(); DistToTube(rsq, nsq, rdotn, rmin, b, d); if (d>0) { s=-b+d; if (s>0) { zi=point[2]+s*dir[2]; if (TMath::Abs(zi)<=dz) { xi=point[0]+s*dir[0]; yi=point[1]+s*dir[1]; if ((xi*cm+yi*sm) >= rmin*cdfi) return s; } } } return TGeoShape::Big(); }