<commit_before>// @(#)root/geom:$Name:  $:$Id: TGeoTube.cxx,v 1.71 2007/01/16 09:04:50 brun Exp $
// Author: <NAME>   24/10/01
// TGeoTube::Contains() and DistFromInside/In() implemented by <NAME>

/*************************************************************************
 * Copyright (C) 1995-2000, <NAME> and <NAME>.               *
 * 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.
//
//_____________________________________________________________________________
//Begin_Html
/*
<img src="gif/t_tube.gif">
*/
//End_Html
//Begin_Html
/*
<img src="gif/t_tubedivR.gif">
*/
//End_Html
//Begin_Html
/*
<img src="gif/t_tubedivstepR.gif">
*/
//End_Html
//Begin_Html
/*
<img src="gif/t_tubedivPHI.gif">
*/
//End_Html
//Begin_Html
/*
<img src="gif/t_tubedivstepPHI.gif">
*/
//End_Html
//Begin_Html
/*
<img src="gif/t_tubedivZ.gif">
*/
//End_Html
//Begin_Html
/*
<img src="gif/t_tubedivstepZ.gif">
*/
//End_Html
//_____________________________________________________________________________
// TGeoTubeSeg - a phi segment of a tube. Has 5 parameters :
//            - the same 3 as a tube;
//            - first phi limit (in degrees)
//            - second phi limit
//
//_____________________________________________________________________________
//Begin_Html
/*
<img src="gif/t_tubseg.gif">
*/
//End_Html
//Begin_Html
/*
<img src="gif/t_tubsegdivstepR.gif">
*/
//End_Html
//Begin_Html
/*
<img src="gif/t_tubsegdivPHI.gif">
*/
//End_Html
//Begin_Html
/*
<img src="gif/t_tubsegdivZ.gif">
*/
//End_Html
//_____________________________________________________________________________
// 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);
//
//_____________________________________________________________________________
//Begin_Html
/*
<img src="gif/t_ctub.gif">
*/
//End_Html

#include "Riostream.h"

#include "TGeoManager.h"
#include "TGeoVolume.h"
#include "TVirtualGeoPainter.h"
#include "TGeoTube.h"
#include "TVirtualPad.h"
#include "TBuffer3D.h"
#include "TBuffer3DTypes.h"
#include "TMath.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
}

//_____________________________________________________________________________
Double_t TGeoTube::Capacity() const
{
// Computes capacity of the shape in [length^3]
   return TGeoTube::Capacity(fRmin,fRmax, fDz);
}   

//_____________________________________________________________________________
Double_t TGeoTube::Capacity(Double_t rmin, Double_t rmax, Double_t dz)
{
// Computes capacity of the shape in [length^3]
   Double_t capacity = 2.*TMath::Pi()*(rmax*rmax-rmin*rmin)*dz;
   return capacity;
}   

//_____________________________________________________________________________
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;<commit_msg>From <NAME> and <NAME>: TGeoTube::MakeBuffer3D() assigns wrong number of segments if tube does not have inner radius.<commit_after>// @(#)root/geom:$Name:  $:$Id: TGeoTube.cxx,v 1.72 2007/04/23 08:58:53 brun Exp $
// Author: <NAME>   24/10/01
// TGeoTube::Contains() and DistFromInside/In() implemented by Mihaela Gheata

/*************************************************************************
 * Copyright (C) 1995-2000, <NAME> and <NAME>.               *
 * 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.
//
//_____________________________________________________________________________
//Begin_Html
/*
<img src="gif/t_tube.gif">
*/
//End_Html
//Begin_Html
/*
<img src="gif/t_tubedivR.gif">
*/
//End_Html
//Begin_Html
/*
<img src="gif/t_tubedivstepR.gif">
*/
//End_Html
//Begin_Html
/*
<img src="gif/t_tubedivPHI.gif">
*/
//End_Html
//Begin_Html
/*
<img src="gif/t_tubedivstepPHI.gif">
*/
//End_Html
//Begin_Html
/*
<img src="gif/t_tubedivZ.gif">
*/
//End_Html
//Begin_Html
/*
<img src="gif/t_tubedivstepZ.gif">
*/
//End_Html
//_____________________________________________________________________________
// TGeoTubeSeg - a phi segment of a tube. Has 5 parameters :
//            - the same 3 as a tube;
//            - first phi limit (in degrees)
//            - second phi limit
//
//_____________________________________________________________________________
//Begin_Html
/*
<img src="gif/t_tubseg.gif">
*/
//End_Html
//Begin_Html
/*
<img src="gif/t_tubsegdivstepR.gif">
*/
//End_Html
//Begin_Html
/*
<img src="gif/t_tubsegdivPHI.gif">
*/
//End_Html
//Begin_Html
/*
<img src="gif/t_tubsegdivZ.gif">
*/
//End_Html
//_____________________________________________________________________________
// 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);
//
//_____________________________________________________________________________
//Begin_Html
/*
<img src="gif/t_ctub.gif">
*/
//End_Html

#include "Riostream.h"

#include "TGeoManager.h"
#include "TGeoVolume.h"
#include "TVirtualGeoPainter.h"
#include "TGeoTube.h"
#include "TVirtualPad.h"
#include "TBuffer3D.h"
#include "TBuffer3DTypes.h"
#include "TMath.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
}

//_____________________________________________________________________________
Double_t TGeoTube::Capacity() const
{
// Computes capacity of the shape in [length^3]
   return TGeoTube::Capacity(fRmin,fRmax, fDz);
}   

//_____________________________________________________________________________
Double_t TGeoTube::Capacity(Double_t rmin, Double_t rmax, Double_t dz)
{
// Computes capacity of the shape in [length^3]
   Double_t capacity = 2.*TMath::Pi()*(rmax*rmax-rmin*rmin)*dz;
   return capacity;
}   

//_____________________________________________________________________________
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;
   if (!HasRmin()) {
      nbPnts = 2*(n+1);
      nbSegs = 5*n;
      nbPols = 3*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;
{'fraction_non_alphanumeric': 0.14763504401343697, 'fraction_numerical': 0.016747174067626217, 'mean_word_length': 3.039219570169357, 'pattern_counts': {'":': 0, '<': 124, '<?xml version=': 0, '>': 172, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4187749', 'n_tokens_mistral': 14258, 'n_tokens_neox': 13791, 'n_words': 3160}