    // @(#)root/geom:$Name:  $:$Id: TGeoChecker.cxx,v 1.33 2005/04/25 07:53:27 brun Exp $
    // Author: Andrei Gheata   01/11/01
    // CheckGeometry(), CheckOverlaps() 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.             *
     *************************************************************************/

    ////////////////////////////////////////////////////////////////////////////////
    //   Geometry checking package.
    //
    // TGeoChecker class provides several geometry checking methods. There are two
    // types of tests that can be performed. One is based on random sampling or
    // ray-tracing and provides a visual check on how navigation methods work for
    // a given geometry. The second actually checks the validity of the geometry
    // definition in terms of overlapping/extruding objects. Both types of checks
    // can be done for a given branch (starting with a given volume) as well as for
    // the geometry as a whole.
    //
    // 1. TGeoChecker::CheckPoint(Double_t x, Double_t y, Double_t z)
    //
    // This method can be called direcly from the TGeoManager class and print a
    // report on how a given point is classified by the modeller: which is the
    // full path to the deepest node containing it, and the (under)estimation
    // of the distance to the closest boundary (safety).
    //
    // 2. TGeoChecker::RandomPoints(Int_t npoints)
    //
    // Can be called from TGeoVolume class. It first draws the volume and its
    // content with the current visualization settings. Then randomly samples points
    // in its bounding box, plotting in the geometry display only the points
    // classified as belonging to visible volumes.
    //
    // 3. TGeoChecker::RandomRays(Int_t nrays, Double_t startx, starty, startz)
    //
    // Can be called and acts in the same way as the previous, but instead of points,
    // rays having random isotropic directions are generated from the given point.
    // A raytracing algorithm propagates all rays untill they exit geometry, plotting
    // all segments crossing visible nodes in the same color as these.
    //
    // 4. TGeoChecker::Test(Int_t npoints)
    //
    // Implementation of TGeoManager::Test(). Computes the time for the modeller
    // to find out "Where am I?" for a given number of random points.
    //
    // 5. TGeoChecker::LegoPlot(ntheta, themin, themax, nphi, phimin, phimax,...)
    //
    // Implementation of TGeoVolume::LegoPlot(). Draws a spherical radiation length
    // lego plot for a given volume, in a given theta/phi range.
    //
    // 6. TGeoChecker::Weigth(Double_t precision)
    //
    // Implementation of TGeoVolume::Weigth(). Estimates the total weigth of a given
    // volume by matrial sampling. Accepts as input the desired precision.
    //
    // Overlap checker
    //-----------------
    //
    //  -> add it from the doc.
    /*
    <img src="gif/t_checker.jpg">
    */
    //End_Html

    #include "TVirtualPad.h"
    #include "TNtuple.h"
    #include "TH2.h"
    #include "TRandom3.h"
    #include "TPolyMarker3D.h"
    #include "TPolyLine3D.h"
    #include "TStopwatch.h"

    #include "TGeoBBox.h"
    #include "TGeoManager.h"
    #include "TGeoOverlap.h"
    #include "TGeoPainter.h"
    #include "TGeoChecker.h"

    #include "TBuffer3D.h"
    #include "TBuffer3DTypes.h"

    // statics and globals

    ClassImp(TGeoChecker)

    //-----------------------------------------------------------------------------
     TGeoChecker::TGeoChecker()
    {
    // Default constructor
       fGeoManager = 0;
       fVsafe      = 0;
    }
    //-----------------------------------------------------------------------------
     TGeoChecker::TGeoChecker(TGeoManager *geom)
    {
    // Constructor for a given geometry
       fGeoManager = geom;
       fVsafe = 0;
    }
    //-----------------------------------------------------------------------------
     TGeoChecker::TGeoChecker(const char * /*treename*/, const char * /*filename*/)
    {
    // constructor
       fGeoManager = gGeoManager;
       fVsafe = 0;
    }
    //-----------------------------------------------------------------------------
     TGeoChecker::~TGeoChecker()
    {
    // Destructor
    }
    //-----------------------------------------------------------------------------
     void TGeoChecker::CheckGeometry(Int_t nrays, Double_t startx, Double_t starty, Double_t startz) const
    {
    // Shoot nrays with random directions from starting point (startx, starty, startz)
    // in the reference frame of this volume. Track each ray until exiting geometry, then
    // shoot backwards from exiting point and compare boundary crossing points.
       Int_t i, j;
       Double_t start[3], end[3];
       Double_t dir[3];
       Double_t dummy[3];
       Double_t eps = 0.;
       Double_t *array1 = new Double_t[3*1000];
       Double_t *array2 = new Double_t[3*1000];
       TObjArray *pma = new TObjArray();
       TPolyMarker3D *pm;
       pm = new TPolyMarker3D();
       pm->SetMarkerColor(2); // error > eps
       pm->SetMarkerStyle(8);
       pm->SetMarkerSize(0.4);
       pma->AddAt(pm, 0);
       pm = new TPolyMarker3D();
       pm->SetMarkerColor(4); // point not found
       pm->SetMarkerStyle(8);
       pm->SetMarkerSize(0.4);
       pma->AddAt(pm, 1);
       pm = new TPolyMarker3D();
       pm->SetMarkerColor(6); // extra point back
       pm->SetMarkerStyle(8);
       pm->SetMarkerSize(0.4);
       pma->AddAt(pm, 2);
       Int_t nelem1, nelem2;
       Int_t dim1=1000, dim2=1000;
       if ((startx==0) && (starty==0) && (startz==0)) eps=1E-3;
       start[0] = startx+eps;
       start[1] = starty+eps;
       start[2] = startz+eps;
       Int_t n10=nrays/10;
       Double_t theta,phi;
       Double_t dw, dwmin, dx, dy, dz;
       Int_t ist1, ist2, ifound;
       for (i=0; i<nrays; i++) {
          if (n10) {
             if ((i%n10) == 0) printf("%i percent\n", Int_t(100*i/nrays));
          }
          phi = 2*TMath::Pi()*gRandom->Rndm();
          theta= TMath::ACos(1.-2.*gRandom->Rndm());
          dir[0]=TMath::Sin(theta)*TMath::Cos(phi);
          dir[1]=TMath::Sin(theta)*TMath::Sin(phi);
          dir[2]=TMath::Cos(theta);
          // shoot direct ray
          nelem1=nelem2=0;
    //      printf("DIRECT %i\n", i);
          ShootRay(&start[0], dir[0], dir[1], dir[2], array1, nelem1, dim1);
          if (!nelem1) continue;
    //      for (j=0; j<nelem1; j++) printf("%i : %f %f %f\n", j, array1[3*j], array1[3*j+1], array1[3*j+2]);
          memcpy(&end[0], &array1[3*(nelem1-1)], 3*sizeof(Double_t));
          // shoot ray backwards
    //      printf("BACK %i\n", i);
          ShootRay(&end[0], -dir[0], -dir[1], -dir[2], array2, nelem2, dim2, &start[0]);
          if (!nelem2) {
             printf("#### NOTHING BACK ###########################\n");
             for (j=0; j<nelem1; j++) {
                pm = (TPolyMarker3D*)pma->At(0);
                pm->SetNextPoint(array1[3*j], array1[3*j+1], array1[3*j+2]);
    	 }
    	 continue;
          }
    //      printf("BACKWARDS\n");
          Int_t k=nelem2>>1;
          for (j=0; j<k; j++) {
             memcpy(&dummy[0], &array2[3*j], 3*sizeof(Double_t));
             memcpy(&array2[3*j], &array2[3*(nelem2-1-j)], 3*sizeof(Double_t));
             memcpy(&array2[3*(nelem2-1-j)], &dummy[0], 3*sizeof(Double_t));
          }
    //      for (j=0; j<nelem2; j++) printf("%i : %f ,%f ,%f   \n", j, array2[3*j], array2[3*j+1], array2[3*j+2]);
          if (nelem1!=nelem2) printf("### DIFFERENT SIZES : nelem1=%i nelem2=%i ##########\n", nelem1, nelem2);
          ist1 = ist2 = 0;
          // check first match

          dx = array1[3*ist1]-array2[3*ist2];
          dy = array1[3*ist1+1]-array2[3*ist2+1];
          dz = array1[3*ist1+2]-array2[3*ist2+2];
          dw = dx*dir[0]+dy*dir[1]+dz*dir[2];
          fGeoManager->SetCurrentPoint(&array1[3*ist1]);
          fGeoManager->FindNode();
    //      printf("%i : %s (%g, %g, %g)\n", ist1, fGeoManager->GetPath(),
    //             array1[3*ist1], array1[3*ist1+1], array1[3*ist1+2]);
          if (TMath::Abs(dw)<1E-4) {
    //         printf("   matching %i (%g, %g, %g)\n", ist2, array2[3*ist2], array2[3*ist2+1], array2[3*ist2+2]);
             ist2++;
          } else {
             printf("### NOT MATCHING %i f:(%f, %f, %f) b:(%f %f %f) DCLOSE=%f\n", ist2, array1[3*ist1], array1[3*ist1+1], array1[3*ist1+2], array2[3*ist2], array2[3*ist2+1], array2[3*ist2+2],dw);
             pm = (TPolyMarker3D*)pma->At(0);
             pm->SetNextPoint(array2[3*ist2], array2[3*ist2+1], array2[3*ist2+2]);
             if (dw<0) {
                // first boundary missed on way back
             } else {
                // first boundary different on way back
                ist2++;
             }
          }

          while ((ist1<nelem1-1) && (ist2<nelem2)) {
             fGeoManager->SetCurrentPoint(&array1[3*ist1+3]);
             fGeoManager->FindNode();
    //         printf("%i : %s (%g, %g, %g)\n", ist1+1, fGeoManager->GetPath(),
    //                array1[3*ist1+3], array1[3*ist1+4], array1[3*ist1+5]);

             dx = array1[3*ist1+3]-array1[3*ist1];
             dy = array1[3*ist1+4]-array1[3*ist1+1];
             dz = array1[3*ist1+5]-array1[3*ist1+2];
             // distance to next point
             dwmin = dx+dir[0]+dy*dir[1]+dz*dir[2];
             while (ist2<nelem2) {
                ifound = 0;
                dx = array2[3*ist2]-array1[3*ist1];
                dy = array2[3*ist2+1]-array1[3*ist1+1];
                dz = array2[3*ist2+2]-array1[3*ist1+2];
                dw = dx+dir[0]+dy*dir[1]+dz*dir[2];
    	    if (TMath::Abs(dw-dwmin)<1E-4) {
    	       ist1++;
    	       ist2++;
    	       break;
    	    }
                if (dw<dwmin) {
                // point found on way back. Check if close enough to ist1+1
                   ifound++;
                   dw = dwmin-dw;
                   if (dw<1E-4) {
                   // point is matching ist1+1
    //                  printf("   matching %i (%g, %g, %g) DCLOSE=%g\n", ist2, array2[3*ist2], array2[3*ist2+1], array2[3*ist2+2], dw);
                      ist2++;
                      ist1++;
                      break;
                   } else {
                   // extra boundary found on way back
                      fGeoManager->SetCurrentPoint(&array2[3*ist2]);
                      fGeoManager->FindNode();
         	          pm = (TPolyMarker3D*)pma->At(2);
    	          pm->SetNextPoint(array2[3*ist2], array2[3*ist2+1], array2[3*ist2+2]);
                      printf("### EXTRA BOUNDARY %i :  %s found at DCLOSE=%f\n", ist2, fGeoManager->GetPath(), dw);
                      ist2++;
                      continue;
                   }
                } else {
                   if (!ifound) {
                      // point ist1+1 not found on way back
                      fGeoManager->SetCurrentPoint(&array1[3*ist1+3]);
                      fGeoManager->FindNode();
    	          pm = (TPolyMarker3D*)pma->At(1);
    	          pm->SetNextPoint(array2[3*ist1+3], array2[3*ist1+4], array2[3*ist1+5]);
                      printf("### BOUNDARY MISSED BACK #########################\n");
                      ist1++;
                      break;
                   } else {
                      ist1++;
                      break;
                   }
                }
             }
          }
       }
       pm = (TPolyMarker3D*)pma->At(0);
       pm->Draw("SAME");
       pm = (TPolyMarker3D*)pma->At(1);
       pm->Draw("SAME");
       pm = (TPolyMarker3D*)pma->At(2);
       pm->Draw("SAME");
       if (gPad) {
          gPad->Modified();
          gPad->Update();
       }
       delete [] array1;
       delete [] array2;
    }

    //-----------------------------------------------------------------------------
     void TGeoChecker::CheckOverlaps(const TGeoVolume *vol, Double_t ovlp, Option_t * /*option*/) const
    {
    // Check illegal overlaps for volume VOL within a limit OVLP.
       if (vol->GetFinder()) return;
       UInt_t nd = vol->GetNdaughters();
       if (!nd) return;
       // first, test if daughters extrude their container
       TGeoShape *shapem = vol->GetShape();
       TGeoShape *shaped;
       TGeoNode * node;
       TGeoMatrix *matrix;
    //   X3DPoints *buff, *buffm;
       TBuffer3D *buff, *buffm;
       buff = new TBuffer3D(TBuffer3DTypes::kGeneric,500,3,0,0,0,0);
       buffm = new TBuffer3D(TBuffer3DTypes::kGeneric,500,3,0,0,0,0);
       Int_t numPoints;
       Bool_t extrude, isextrusion, isoverlapping;
       TGeoOverlap *nodeovlp = 0;
       Bool_t ismany;
       Double_t *points, *pointsm;
       Double_t local[3];
       Double_t point[3];
       Double_t safety = TGeoShape::Big();
       UInt_t id, ip;
       // first, test if any of container vertices is inside some daughter
       // first, test if daughters extrude their container
       numPoints = shapem->GetNmeshVertices();
       buffm->SetRawSizes(numPoints, 3*numPoints, 0, 0, 0, 0);
       pointsm = buffm->fPnts;
       shapem->SetPoints(pointsm);

       for (id=0; id<nd; id++) {
          node = vol->GetNode(id);
          shaped = node->GetVolume()->GetShape();
          numPoints = shaped->GetNmeshVertices();
          buff->SetRawSizes(numPoints, 3*numPoints, 0, 0, 0, 0);
          points = buff->fPnts;
          shaped->SetPoints(points);
          matrix = node->GetMatrix();
          ismany = node->IsOverlapping();
          isextrusion=kFALSE;
          // loop all points of the daughter
          for (ip=0; ip<buff->NbPnts(); ip++) {
             memcpy(local, &points[3*ip], 3*sizeof(Double_t));
    	      matrix->LocalToMaster(local, point);
    	      extrude = !shapem->Contains(point);
    	      if (extrude) {
    	         safety = shapem->Safety(point, kFALSE);
    	         if (safety<ovlp) extrude=kFALSE;
    	      }
    	      if (extrude) {
                if (!isextrusion) {
                   isextrusion = kTRUE;
                   char *name = new char[20];
                   sprintf(name,"%s_x_%i", vol->GetName(),id);
                   nodeovlp = new TGeoExtrusion(name, (TGeoVolume*)vol,id,safety);
                   nodeovlp->SetNextPoint(point[0],point[1],point[2]);
                   fGeoManager->AddOverlap(nodeovlp);
                } else {
                   if (safety>nodeovlp->GetOverlap()) nodeovlp->SetOverlap(safety);
                   nodeovlp->SetNextPoint(point[0],point[1],point[2]);
                }
    	      }
          }
          // loop all points of the mother
          for (ip=0; ip<buffm->NbPnts(); ip++) {
             memcpy(point, &pointsm[3*ip], 3*sizeof(Double_t));
    	      matrix->MasterToLocal(point, local);
             extrude = shaped->Contains(local);
             if (extrude) {
                // skip points on mother mesh that have no neghbourhood ouside mother
                safety = shapem->Safety(point,kTRUE);
                if (safety>1E-4) {
                   extrude = kFALSE;
                } else {
                   safety = shaped->Safety(local,kTRUE);
    	            if (safety<ovlp) extrude=kFALSE;
                }
             }
             if (extrude) {
                if (!isextrusion) {
                   isextrusion = kTRUE;
                   char *name = new char[20];
                   sprintf(name,"%s_mo_%i", vol->GetName(),id);
                   nodeovlp = new TGeoExtrusion(name, (TGeoVolume*)vol,id,safety);
                   nodeovlp->SetNextPoint(point[0],point[1],point[2]);
                   fGeoManager->AddOverlap(nodeovlp);
                } else {
                   if (safety>nodeovlp->GetOverlap()) nodeovlp->SetOverlap(safety);
                   nodeovlp->SetNextPoint(point[0],point[1],point[2]);
                }
             }
          }
       }
       // now check if the daughters overlap with each other
       if (nd<2) {
          delete buff;
          delete buffm;
          return;
       }
       TGeoVoxelFinder *vox = vol->GetVoxels();
       TGeoNode *node1;
       TGeoMatrix *matrix1;
       Bool_t ismany1, overlap;
       Int_t novlp;
       Int_t *ovlps;
       Int_t ko;
       UInt_t io;
       for (id=0; id<nd; id++) {  // loop all nodes (node, shapem, matrix, pointsm, buffm)
          node = vol->GetNode(id);
          ismany = node->IsOverlapping();
          if (ismany) continue;
          shapem = node->GetVolume()->GetShape();
          matrix = node->GetMatrix();
          if (!vox) continue;
          vox->FindOverlaps(id);
          ovlps = node->GetOverlaps(novlp);
          if (!ovlps) continue;
          numPoints = shapem->GetNmeshVertices();
          buffm->SetRawSizes(numPoints, 3*numPoints, 0, 0, 0, 0);    // node buffer
          pointsm = buffm->fPnts;
          shapem->SetPoints(pointsm);
          for (ko=0; ko<novlp; ko++) { // loop all possible overlapping candidates
             io = ovlps[ko];           // (node1, shaped, matrix1, points, buff)
             if (io<id) continue;
             node1 = vol->GetNode(io);
             ismany1 = node1->IsOverlapping();
             if (ismany1) continue;
             matrix1 = node1->GetMatrix();
             shaped = node1->GetVolume()->GetShape();
             numPoints = shaped->GetNmeshVertices(); // overlapping candidate buffer
             buff->SetRawSizes(numPoints, 3*numPoints, 0, 0, 0, 0);
             points = buff->fPnts;
             shaped->SetPoints(points);
             // loop all points of candidate
             overlap = kFALSE;
             isoverlapping = kFALSE;
             for (ip=0; ip<buff->NbPnts(); ip++) {
                memcpy(local, &points[3*ip], 3*sizeof(Double_t));
    	         matrix1->LocalToMaster(local, point);
                matrix->MasterToLocal(point, local); // now point in local reference of node
    	         overlap = shapem->Contains(local);
    	         if (overlap) {
    	            safety = shapem->Safety(local, kTRUE);
    	            if (safety<ovlp) overlap=kFALSE;
    	         }
    	         if (overlap) {
                   if (!isoverlapping) {
                      isoverlapping = kTRUE;
                      char *name = new char[20];
                      sprintf(name,"%s_o_%i_%i", vol->GetName(),id,io);
                      nodeovlp = new TGeoNodeOverlap(name, (TGeoVolume*)vol,id,io,safety);
                      nodeovlp->SetNextPoint(point[0],point[1],point[2]);
                      fGeoManager->AddOverlap(nodeovlp);
                   } else {
                      if (safety>nodeovlp->GetOverlap()) nodeovlp->SetOverlap(safety);
                      nodeovlp->SetNextPoint(point[0],point[1],point[2]);
                   }
    	         }
             }
             // loop all points of node
             for (ip=0; ip<buffm->NbPnts(); ip++) {
                memcpy(local, &pointsm[3*ip], 3*sizeof(Double_t));
    	         matrix->LocalToMaster(local, point);
                matrix1->MasterToLocal(point, local); // now point in local reference of node
    	         overlap = shaped->Contains(local);
    	         if (overlap) {
    	            safety = shaped->Safety(local, kTRUE);
    	            if (safety<ovlp) overlap=kFALSE;
    	         }
    	         if (overlap) {
                   if (!isoverlapping) {
                      isoverlapping = kTRUE;
                      char *name = new char[20];
                      sprintf(name,"%s_o_%i_%i", vol->GetName(),id,io);
                      nodeovlp = new TGeoNodeOverlap(name, (TGeoVolume*)vol,id,io,safety);
                      nodeovlp->SetNextPoint(point[0],point[1],point[2]);
                      fGeoManager->AddOverlap(nodeovlp);
                   } else {
                      if (safety>nodeovlp->GetOverlap()) nodeovlp->SetOverlap(safety);
                      nodeovlp->SetNextPoint(point[0],point[1],point[2]);
                   }
    	         }
             }
          }
          node->SetOverlaps(0,0);
       }
       delete buff;
       delete buffm;
    }

    //-----------------------------------------------------------------------------
     void TGeoChecker::PrintOverlaps() const
    {
    // Print the current list of overlaps held by the manager class.
       TIter next(fGeoManager->GetListOfOverlaps());
       TGeoOverlap *ov;
       printf("=== Overlaps for %s ===\n", fGeoManager->GetName());
       while ((ov=(TGeoOverlap*)next())) ov->PrintInfo();
    }

    //-----------------------------------------------------------------------------
     void TGeoChecker::CheckPoint(Double_t x, Double_t y, Double_t z, Option_t *)
    {
    //--- Draw point (x,y,z) over the picture of the daughers of the volume containing this point.
    //   Generates a report regarding the path to the node containing this point and the distance to
    //   the closest boundary.

       Double_t point[3];
       Double_t local[3];
       point[0] = x;
       point[1] = y;
       point[2] = z;
       TGeoVolume *vol = fGeoManager->GetTopVolume();
       if (fVsafe) {
          TGeoNode *old = fVsafe->GetNode("SAFETY_1");
          if (old) fVsafe->GetNodes()->RemoveAt(vol->GetNdaughters()-1);
       }
    //   if (vol != fGeoManager->GetMasterVolume()) fGeoManager->RestoreMasterVolume();
       TGeoNode *node = fGeoManager->FindNode(point[0], point[1], point[2]);
       fGeoManager->MasterToLocal(point, local);
       // get current node
       printf("===  Check current point : (%g, %g, %g) ===\n", point[0], point[1], point[2]);
       printf("  - path : %s\n", fGeoManager->GetPath());
       // get corresponding volume
       if (node) vol = node->GetVolume();
       // compute safety distance (distance to boundary ignored)
       Double_t close = fGeoManager->Safety();
       printf("Safety radius : %f\n", close);
       if (close>1E-4) {
          TGeoVolume *sph = fGeoManager->MakeSphere("SAFETY", vol->GetMedium(), 0, close, 0,180,0,360);
          sph->SetLineColor(2);
          sph->SetLineStyle(3);
          vol->AddNode(sph,1,new TGeoTranslation(local[0], local[1], local[2]));
          fVsafe = vol;
       }
       TPolyMarker3D *pm = new TPolyMarker3D();
       pm->SetMarkerColor(2);
       pm->SetMarkerStyle(8);
       pm->SetMarkerSize(0.5);
       pm->SetNextPoint(local[0], local[1], local[2]);
       if (vol->GetNdaughters()<2) fGeoManager->SetTopVisible();
       else fGeoManager->SetTopVisible(kFALSE);
       fGeoManager->SetVisLevel(1);
       if (!vol->IsVisible()) vol->SetVisibility(kTRUE);
       vol->Draw();
       pm->Draw("SAME");
       gPad->Modified();
       gPad->Update();
    }
    //-----------------------------------------------------------------------------
     TH2F *TGeoChecker::LegoPlot(Int_t ntheta, Double_t themin, Double_t themax,
                                Int_t nphi,   Double_t phimin, Double_t phimax,
                                Double_t /*rmin*/, Double_t /*rmax*/, Option_t *option)
    {
    // Generate a lego plot fot the top volume, according to option.
       TH2F *hist = new TH2F("lego", option, nphi, phimin, phimax, ntheta, themin, themax);

       Double_t degrad = TMath::Pi()/180.;
       Double_t theta, phi, step, matprop, x;
       Double_t start[3];
       Double_t dir[3];
       TGeoNode *startnode, *endnode;
       Int_t i;  // loop index for phi
       Int_t j;  // loop index for theta
       Int_t ntot = ntheta * nphi;
       Int_t n10 = ntot/10;
       Int_t igen = 0, iloop=0;
       printf("=== Lego plot sph. => nrays=%i\n", ntot);
       for (i=1; i<=nphi; i++) {
          for (j=1; j<=ntheta; j++) {
             igen++;
             if (n10) {
                if ((igen%n10) == 0) printf("%i percent\n", Int_t(100*igen/ntot));
             }
             x = 0;
             theta = hist->GetYaxis()->GetBinCenter(j);
             phi   = hist->GetXaxis()->GetBinCenter(i)+1E-3;
             start[0] = start[1] = start[2] = 1E-3;
             dir[0]=TMath::Sin(theta*degrad)*TMath::Cos(phi*degrad);
             dir[1]=TMath::Sin(theta*degrad)*TMath::Sin(phi*degrad);
             dir[2]=TMath::Cos(theta*degrad);
             fGeoManager->InitTrack(&start[0], &dir[0]);
             startnode = fGeoManager->GetCurrentNode();
             if (fGeoManager->IsOutside()) startnode=0;
             if (startnode) {
                matprop = startnode->GetVolume()->GetMaterial()->GetRadLen();
             } else {
                matprop = 0.;
             }
             fGeoManager->FindNextBoundary();
    //         fGeoManager->IsStepEntering();
             // find where we end-up
             endnode = fGeoManager->Step();
             step = fGeoManager->GetStep();
             while (step<1E10) {
                // now see if we can make an other step
                iloop=0;
                while (!fGeoManager->IsEntering()) {
                   iloop++;
                   fGeoManager->SetStep(1E-3);
                   step += 1E-3;
                   endnode = fGeoManager->Step();
                }
                if (iloop>1000) printf("%i steps\n", iloop);
                if (matprop>0) {
                   x += step/matprop;
                }
                if (endnode==0 && step>1E10) break;
                // generate an extra step to cross boundary
                startnode = endnode;
                if (startnode) {
                   matprop = startnode->GetVolume()->GetMaterial()->GetRadLen();
                } else {
                   matprop = 0.;
                }

                fGeoManager->FindNextBoundary();
                endnode = fGeoManager->Step();
                step = fGeoManager->GetStep();
             }
             hist->Fill(phi, theta, x);
          }
       }
       return hist;
    }
    //______________________________________________________________________________
     void TGeoChecker::RandomPoints(TGeoVolume *vol, Int_t npoints, Option_t *option)
    {
    // Draw random points in the bounding box of a volume.
       if (!vol) return;
       gRandom = new TRandom3();
       vol->VisibleDaughters(kTRUE);
       vol->Draw();
       TString opt = option;
       opt.ToLower();
       TObjArray *pm = new TObjArray(128);
       TPolyMarker3D *marker = 0;
       const TGeoShape *shape = vol->GetShape();
       TGeoBBox *box = (TGeoBBox *)shape;
       Double_t dx = box->GetDX();
       Double_t dy = box->GetDY();
       Double_t dz = box->GetDZ();
       Double_t ox = (box->GetOrigin())[0];
       Double_t oy = (box->GetOrigin())[1];
       Double_t oz = (box->GetOrigin())[2];
       Double_t *xyz = new Double_t[3];
       printf("Random box : %f, %f, %f\n", dx, dy, dz);
       TGeoNode *node = 0;
       printf("Start... %i points\n", npoints);
       Int_t i=0;
       Int_t igen=0;
       Int_t ic = 0;
       Int_t n10 = npoints/10;
       Double_t ratio=0;
       while (igen<npoints) {
          xyz[0] = ox-dx+2*dx*gRandom->Rndm();
          xyz[1] = oy-dy+2*dy*gRandom->Rndm();
          xyz[2] = oz-dz+2*dz*gRandom->Rndm();
          fGeoManager->SetCurrentPoint(xyz);
          igen++;
          if (n10) {
             if ((igen%n10) == 0) printf("%i percent\n", Int_t(100*igen/npoints));
          }
          node = fGeoManager->FindNode();
          if (!node) continue;
          if (!node->IsOnScreen()) continue;
          // draw only points in overlapping/non-overlapping volumes
          if (opt.Contains("many") && !node->IsOverlapping()) continue;
          if (opt.Contains("only") && node->IsOverlapping()) continue;
          ic = node->GetColour();
          if ((ic<0) || (ic>=128)) ic = 0;
          marker = (TPolyMarker3D*)pm->At(ic);
          if (!marker) {
             marker = new TPolyMarker3D();
             marker->SetMarkerColor(ic);
             marker->SetMarkerStyle(8);
             marker->SetMarkerSize(0.4);
             pm->AddAt(marker, ic);
          }
          marker->SetNextPoint(xyz[0], xyz[1], xyz[2]);
          i++;
       }
       printf("Number of visible points : %i\n", i);
       ratio = (Double_t)i/(Double_t)igen;
       printf("efficiency : %g\n", ratio);
       for (Int_t m=0; m<128; m++) {
          marker = (TPolyMarker3D*)pm->At(m);
          if (marker) marker->Draw("SAME");
       }
       fGeoManager->GetTopVolume()->VisibleDaughters(kFALSE);
       printf("---Daughters of %s made invisible.\n", fGeoManager->GetTopVolume()->GetName());
       printf("---Make them visible with : gGeoManager->GetTopVolume()->VisibleDaughters();\n");
       delete pm;
       delete xyz;
    }
    //-----------------------------------------------------------------------------
     void TGeoChecker::RandomRays(Int_t nrays, Double_t startx, Double_t starty, Double_t startz)
    {
    // Randomly shoot nrays from point (startx,starty,startz) and plot intersections
    // with surfaces for current top node.
       TObjArray *pm = new TObjArray(128);
       TPolyLine3D *line = 0;
       gRandom = new TRandom3();
       TGeoVolume *vol=fGeoManager->GetTopVolume();
       vol->VisibleDaughters(kTRUE);

       Double_t start[3];
       Double_t dir[3];
       Int_t istep= 0;
       Double_t *point = fGeoManager->GetCurrentPoint();
       vol->Draw();
       printf("Start... %i rays\n", nrays);
       TGeoNode *startnode, *endnode;
       Bool_t vis1,vis2;
       Int_t i=0;
       Int_t ipoint;
       Int_t itot=0;
       Int_t n10=nrays/10;
       Double_t theta,phi, step;
       while (itot<nrays) {
          itot++;
          ipoint = 0;
          if (n10) {
             if ((itot%n10) == 0) printf("%i percent\n", Int_t(100*itot/nrays));
          }
          start[0] = startx;
          start[1] = starty;
          start[2] = startz;
          phi = 2*TMath::Pi()*gRandom->Rndm();
          theta= TMath::ACos(1.-2.*gRandom->Rndm());
          dir[0]=TMath::Sin(theta)*TMath::Cos(phi);
          dir[1]=TMath::Sin(theta)*TMath::Sin(phi);
          dir[2]=TMath::Cos(theta);
          startnode = fGeoManager->InitTrack(start[0],start[1],start[2], dir[0],dir[1],dir[2]);
          line = 0;
          if (fGeoManager->IsOutside()) startnode=0;
          vis1 = (startnode)?(startnode->IsOnScreen()):kFALSE;
          if (vis1) {
             line = new TPolyLine3D(2);
             line->SetLineColor(startnode->GetVolume()->GetLineColor());
             line->SetPoint(ipoint++, startx, starty, startz);
             i++;
             pm->Add(line);
          }
          // find where we end-up
          fGeoManager->FindNextBoundary();
          step = fGeoManager->GetStep();
          endnode = fGeoManager->Step();
          vis2 = (endnode)?(endnode->IsOnScreen()):kFALSE;
          while (step<1E10) {
             istep = 0;
             while (!fGeoManager->IsEntering()) {
                istep++;
                if (istep>1E4) break;
                fGeoManager->SetStep(1E-3);
                endnode = fGeoManager->Step();
                step += 1E-3;
             }
             if (istep>1E4) break;
    //         if (istep) printf("ADDED : %f (%i steps)\n", istep*1E-3, istep);
             vis2 = (endnode)?(endnode->IsOnScreen()):kFALSE;
             if (ipoint>0) {
             // old visible node had an entry point -> finish segment
                line->SetPoint(ipoint, point[0], point[1], point[2]);
                ipoint = 0;
                line   = 0;
             }
             if (vis2) {
                // create new segment
                line = new TPolyLine3D(2);
                line->SetLineColor(endnode->GetVolume()->GetLineColor());
                line->SetPoint(ipoint++, point[0], point[1], point[2]);
                i++;
                pm->Add(line);
             }
             // now see if we can make an other step
             if (endnode==0 && step>1E10) break;
             // generate an extra step to cross boundary
             startnode = endnode;
             fGeoManager->FindNextBoundary();
             step = fGeoManager->GetStep();
             endnode = fGeoManager->Step();
          }
       }
       // draw all segments
       for (Int_t m=0; m<pm->GetEntriesFast(); m++) {
          line = (TPolyLine3D*)pm->At(m);
          if (line) line->Draw("SAME");
       }
       printf("number of segments : %i\n", i);
       fGeoManager->GetTopVolume()->VisibleDaughters(kFALSE);
       printf("---Daughters of %s made invisible.\n", fGeoManager->GetTopVolume()->GetName());
       printf("---Make them visible with : gGeoManager->GetTopVolume()->VisibleDaughters();\n");
       delete pm;
    }
    //-----------------------------------------------------------------------------
     TGeoNode *TGeoChecker::SamplePoints(Int_t npoints, Double_t &dist, Double_t epsil,
                                        const char* g3path)
    {
    // shoot npoints randomly in a box of 1E-5 arround current point.
    // return minimum distance to points outside
       // make sure that path to current node is updated
       // get the response of tgeo
       TGeoNode *node = fGeoManager->FindNode();
       TGeoNode *nodegeo = 0;
       TGeoNode *nodeg3 = 0;
       TGeoNode *solg3 = 0;
       if (!node) {dist=-1; return 0;}
       gRandom = new TRandom3();
       Bool_t hasg3 = kFALSE;
       if (strlen(g3path)) hasg3 = kTRUE;
       char geopath[200];
       sprintf(geopath, "%s\n", fGeoManager->GetPath());
       dist = 1E10;
       TString common = "";
       // cd to common path
       Double_t point[3];
       Double_t closest[3];
       TGeoNode *node1 = 0;
       TGeoNode *node_close = 0;
       dist = 1E10;
       Double_t dist1 = 0;
       // initialize size of random box to epsil
       Double_t eps[3];
       eps[0] = epsil; eps[1]=epsil; eps[2]=epsil;
       Double_t *pointg = fGeoManager->GetCurrentPoint();
       if (hasg3) {
          TString spath = geopath;
          TString name = "";
          Int_t index=0;
          while (index>=0) {
             index = spath.Index("/", index+1);
             if (index>0) {
                name = spath(0, index);
                if (strstr(g3path, name.Data())) {
                   common = name;
                   continue;
                } else break;
             }
          }
          // if g3 response was given, cd to common path
          if (strlen(common.Data())) {
             while (strcmp(fGeoManager->GetPath(), common.Data()) && fGeoManager->GetLevel()) {
                nodegeo = fGeoManager->GetCurrentNode();
                fGeoManager->CdUp();
             }
             fGeoManager->cd(g3path);
             solg3 = fGeoManager->GetCurrentNode();
             while (strcmp(fGeoManager->GetPath(), common.Data()) && fGeoManager->GetLevel()) {
                nodeg3 = fGeoManager->GetCurrentNode();
                fGeoManager->CdUp();
             }
             if (!nodegeo) return 0;
             if (!nodeg3) return 0;
             fGeoManager->cd(common.Data());
             fGeoManager->MasterToLocal(fGeoManager->GetCurrentPoint(), &point[0]);
             Double_t xyz[3], local[3];
             for (Int_t i=0; i<npoints; i++) {
                xyz[0] = point[0] - eps[0] + 2*eps[0]*gRandom->Rndm();
                xyz[1] = point[1] - eps[1] + 2*eps[1]*gRandom->Rndm();
                xyz[2] = point[2] - eps[2] + 2*eps[2]*gRandom->Rndm();
                nodeg3->MasterToLocal(&xyz[0], &local[0]);
                if (!nodeg3->GetVolume()->Contains(&local[0])) continue;
                dist1 = TMath::Sqrt((xyz[0]-point[0])*(xyz[0]-point[0])+
                       (xyz[1]-point[1])*(xyz[1]-point[1])+(xyz[2]-point[2])*(xyz[2]-point[2]));
                if (dist1<dist) {
                // save node and closest point
                   dist = dist1;
                   node_close = solg3;
                   // make the random box smaller
                   eps[0] = TMath::Abs(point[0]-pointg[0]);
                   eps[1] = TMath::Abs(point[1]-pointg[1]);
                   eps[2] = TMath::Abs(point[2]-pointg[2]);
                }
             }
          }
          if (!node_close) dist = -1;
          return node_close;
       }

    //   gRandom = new TRandom3();
       // save current point
       memcpy(&point[0], pointg, 3*sizeof(Double_t));
       for (Int_t i=0; i<npoints; i++) {
          // generate a random point in MARS
          pointg[0] = point[0] - eps[0] + 2*eps[0]*gRandom->Rndm();
          pointg[1] = point[1] - eps[1] + 2*eps[1]*gRandom->Rndm();
          pointg[2] = point[2] - eps[2] + 2*eps[2]*gRandom->Rndm();
          // check if new node is different from the old one
          if (node1!=node) {
             dist1 = TMath::Sqrt((point[0]-pointg[0])*(point[0]-pointg[0])+
                     (point[1]-pointg[1])*(point[1]-pointg[1])+(point[2]-pointg[2])*(point[2]-pointg[2]));
             if (dist1<dist) {
                dist = dist1;
                node_close = node1;
                memcpy(&closest[0], pointg, 3*sizeof(Double_t));
                // make the random box smaller
                eps[0] = TMath::Abs(point[0]-pointg[0]);
                eps[1] = TMath::Abs(point[1]-pointg[1]);
                eps[2] = TMath::Abs(point[2]-pointg[2]);
             }
          }
       }
       // restore the original point and path
       memcpy(pointg, &point[0], 3*sizeof(Double_t));
       fGeoManager->FindNode();  // really needed ?
       if (!node_close) dist=-1;
       return node_close;
    }
    //-----------------------------------------------------------------------------
     void TGeoChecker::ShootRay(Double_t *start, Double_t dirx, Double_t diry, Double_t dirz, Double_t *array, Int_t &nelem, Int_t &dim, Double_t *endpoint) const
    {
    // Shoot one ray from start point with direction (dirx,diry,dirz). Fills input array
    // with points just after boundary crossings.
    //   Int_t array_dimension = 3*dim;
       nelem = 0;
       Int_t istep = 0;
       if (!dim) {
          printf("empty input array\n");
          return;
       }
    //   fGeoManager->CdTop();
       Double_t *point = fGeoManager->GetCurrentPoint();
       TGeoNode *endnode;
       Bool_t is_entering;
       Double_t step, forward;
       Double_t dir[3];
       dir[0] = dirx;
       dir[1] = diry;
       dir[2] = dirz;
       fGeoManager->InitTrack(start, &dir[0]);
       fGeoManager->GetCurrentNode();
    //   printf("Start : (%f,%f,%f)\n", point[0], point[1], point[2]);
       fGeoManager->FindNextBoundary();
       step = fGeoManager->GetStep();
    //   printf("---next : at step=%f\n", step);
       if (step>1E10) return;
       endnode = fGeoManager->Step();
       is_entering = fGeoManager->IsEntering();
       while (step<1E10) {
          if (endpoint) {
             forward = dirx*(endpoint[0]-point[0])+diry*(endpoint[1]-point[1])+dirz*(endpoint[2]-point[2]);
             if (forward<1E-3) {
    //	    printf("exit : Passed start point. nelem=%i\n", nelem);
    	    return;
    	 }
          }
          if (is_entering) {
             if (nelem>=dim) {
                Double_t *temparray = new Double_t[3*(dim+20)];
                memcpy(temparray, array, 3*dim*sizeof(Double_t));
                delete [] array;
                array = temparray;
    	          dim += 20;
             }
             memcpy(&array[3*nelem], point, 3*sizeof(Double_t));
    //         printf("%i (%f, %f, %f) step=%f\n", nelem, point[0], point[1], point[2], step);
             nelem++;
          } else {
             if (endnode==0 && step>1E10) {
    //	    printf("exit : NULL endnode. nelem=%i\n", nelem);
    	    return;
    	 }
             if (!fGeoManager->IsEntering()) {
    //            if (startnode) printf("stepping %f from (%f, %f, %f) inside %s...\n", step,point[0], point[1], point[2], startnode->GetName());
    //            else printf("stepping %f from (%f, %f, %f) OUTSIDE...\n", step,point[0], point[1], point[2]);
                istep = 0;
             }
             while (!fGeoManager->IsEntering()) {
                istep++;
    	    if (istep>1E3) {
    //	       Error("ShootRay", "more than 1000 steps. Step was %f", step);
    	       nelem = 0;
    	       return;
    	    }
                fGeoManager->SetStep(1E-5);
                endnode = fGeoManager->Step();
             }
             if (istep>0) printf("%i steps\n", istep);
             if (nelem>=dim) {
                Double_t *temparray = new Double_t[3*(dim+20)];
                memcpy(temparray, array, 3*dim*sizeof(Double_t));
                delete [] array;
                array = temparray;
    	          dim += 20;
             }
             memcpy(&array[3*nelem], point, 3*sizeof(Double_t));
    //         if (endnode) printf("%i (%f, %f, %f) step=%f\n", nelem, point[0], point[1], point[2], step);
             nelem++;
             is_entering = kTRUE;
          }
          fGeoManager->FindNextBoundary();
          step = fGeoManager->GetStep();
    //      printf("---next at step=%f\n", step);
          endnode = fGeoManager->Step();
          is_entering = fGeoManager->IsEntering();
       }
    //   printf("exit : INFINITE step. nelem=%i\n", nelem);
    }
    //-----------------------------------------------------------------------------
     void TGeoChecker::Test(Int_t npoints, Option_t *option)
    {
       // Check time of finding "Where am I" for n points.
       gRandom= new TRandom3();
       Bool_t recheck = !strcmp(option, "RECHECK");
       if (recheck) printf("RECHECK\n");
       const TGeoShape *shape = fGeoManager->GetTopVolume()->GetShape();
       Double_t dx = ((TGeoBBox*)shape)->GetDX();
       Double_t dy = ((TGeoBBox*)shape)->GetDY();
       Double_t dz = ((TGeoBBox*)shape)->GetDZ();
       Double_t ox = (((TGeoBBox*)shape)->GetOrigin())[0];
       Double_t oy = (((TGeoBBox*)shape)->GetOrigin())[1];
       Double_t oz = (((TGeoBBox*)shape)->GetOrigin())[2];
       Double_t *xyz = new Double_t[3*npoints];
       TStopwatch *timer = new TStopwatch();
       printf("Random box : %f, %f, %f\n", dx, dy, dz);
       timer->Start(kFALSE);
       Int_t i;
       for (i=0; i<npoints; i++) {
          xyz[3*i] = ox-dx+2*dx*gRandom->Rndm();
          xyz[3*i+1] = oy-dy+2*dy*gRandom->Rndm();
          xyz[3*i+2] = oz-dz+2*dz*gRandom->Rndm();
       }
       timer->Stop();
       printf("Generation time :\n");
       timer->Print();
       timer->Reset();
       TGeoNode *node, *node1;
       printf("Start... %i points\n", npoints);
       timer->Start(kFALSE);
       for (i=0; i<npoints; i++) {
          fGeoManager->SetCurrentPoint(xyz+3*i);
          if (recheck) fGeoManager->CdTop();
          node = fGeoManager->FindNode();
          if (recheck) {
             node1 = fGeoManager->FindNode();
             if (node1 != node) {
                printf("Difference for x=%g y=%g z=%g\n", xyz[3*i], xyz[3*i+1], xyz[3*i+2]);
                printf(" from top : %s\n", node->GetName());
                printf(" redo     : %s\n", fGeoManager->GetPath());
             }
          }
       }
       timer->Stop();
       timer->Print();
       delete xyz;
       delete timer;
    }
    //-----------------------------------------------------------------------------
     void TGeoChecker::TestOverlaps(const char* path)
    {
    //--- Geometry overlap checker based on sampling.
       if (fGeoManager->GetTopVolume()!=fGeoManager->GetMasterVolume()) fGeoManager->RestoreMasterVolume();
       printf("Checking overlaps for path :\n");
       if (!fGeoManager->cd(path)) return;
       TGeoNode *checked = fGeoManager->GetCurrentNode();
       checked->InspectNode();
       // shoot 1E4 points in the shape of the current volume
       gRandom= new TRandom3();
       Int_t npoints = 1000000;
       Double_t big = 1E6;
       Double_t xmin = big;
       Double_t xmax = -big;
       Double_t ymin = big;
       Double_t ymax = -big;
       Double_t zmin = big;
       Double_t zmax = -big;
       TObjArray *pm = new TObjArray(128);
       TPolyMarker3D *marker = 0;
       TPolyMarker3D *markthis = new TPolyMarker3D();
       markthis->SetMarkerColor(5);
       TNtuple *ntpl = new TNtuple("ntpl","random points","x:y:z");
       TGeoShape *shape = fGeoManager->GetCurrentNode()->GetVolume()->GetShape();
       Double_t *point = new Double_t[3];
       Double_t dx = ((TGeoBBox*)shape)->GetDX();
       Double_t dy = ((TGeoBBox*)shape)->GetDY();
       Double_t dz = ((TGeoBBox*)shape)->GetDZ();
       Double_t ox = (((TGeoBBox*)shape)->GetOrigin())[0];
       Double_t oy = (((TGeoBBox*)shape)->GetOrigin())[1];
       Double_t oz = (((TGeoBBox*)shape)->GetOrigin())[2];
       Double_t *xyz = new Double_t[3*npoints];
       Int_t i=0;
       printf("Generating %i points inside %s\n", npoints, fGeoManager->GetPath());
       while (i<npoints) {
          point[0] = ox-dx+2*dx*gRandom->Rndm();
          point[1] = oy-dy+2*dy*gRandom->Rndm();
          point[2] = oz-dz+2*dz*gRandom->Rndm();
          if (!shape->Contains(point)) continue;
          // convert each point to MARS
    //      printf("local  %9.3f %9.3f %9.3f\n", point[0], point[1], point[2]);
          fGeoManager->LocalToMaster(point, &xyz[3*i]);
    //      printf("master %9.3f %9.3f %9.3f\n", xyz[3*i], xyz[3*i+1], xyz[3*i+2]);
          xmin = TMath::Min(xmin, xyz[3*i]);
          xmax = TMath::Max(xmax, xyz[3*i]);
          ymin = TMath::Min(ymin, xyz[3*i+1]);
          ymax = TMath::Max(ymax, xyz[3*i+1]);
          zmin = TMath::Min(zmin, xyz[3*i+2]);
          zmax = TMath::Max(zmax, xyz[3*i+2]);
          i++;
       }
       delete point;
       ntpl->Fill(xmin,ymin,zmin);
       ntpl->Fill(xmax,ymin,zmin);
       ntpl->Fill(xmin,ymax,zmin);
       ntpl->Fill(xmax,ymax,zmin);
       ntpl->Fill(xmin,ymin,zmax);
       ntpl->Fill(xmax,ymin,zmax);
       ntpl->Fill(xmin,ymax,zmax);
       ntpl->Fill(xmax,ymax,zmax);
       ntpl->Draw("z:y:x");

       // shoot the poins in the geometry
       TGeoNode *node;
       TString cpath;
       Int_t ic=0;
       TObjArray *overlaps = new TObjArray();
       printf("using FindNode...\n");
       for (Int_t j=0; j<npoints; j++) {
          // always start from top level (testing only)
          fGeoManager->CdTop();
          fGeoManager->SetCurrentPoint(&xyz[3*j]);
          node = fGeoManager->FindNode();
          cpath = fGeoManager->GetPath();
          if (cpath.Contains(path)) {
             markthis->SetNextPoint(xyz[3*j], xyz[3*j+1], xyz[3*j+2]);
             continue;
          }
          // current point is found in an overlapping node
          if (!node) ic=128;
          else ic = node->GetVolume()->GetLineColor();
          if (ic >= 128) ic = 0;
          marker = (TPolyMarker3D*)pm->At(ic);
          if (!marker) {
             marker = new TPolyMarker3D();
             marker->SetMarkerColor(ic);
             pm->AddAt(marker, ic);
          }
          // draw the overlapping point
          marker->SetNextPoint(xyz[3*j], xyz[3*j+1], xyz[3*j+2]);
          if (node) {
             if (overlaps->IndexOf(node) < 0) overlaps->Add(node);
          }
       }
       // draw all overlapping points
       for (Int_t m=0; m<128; m++) {
          marker = (TPolyMarker3D*)pm->At(m);
    //      if (marker) marker->Draw("SAME");
       }
       markthis->Draw("SAME");
       if (gPad) gPad->Update();
       // display overlaps
       if (overlaps->GetEntriesFast()) {
          printf("list of overlapping nodes :\n");
          for (i=0; i<overlaps->GetEntriesFast(); i++) {
             node = (TGeoNode*)overlaps->At(i);
             if (node->IsOverlapping()) printf("%s  MANY\n", node->GetName());
             else printf("%s  ONLY\n", node->GetName());
          }
       } else printf("No overlaps\n");
       delete ntpl;
       delete pm;
       delete xyz;
       delete overlaps;
    }

    //-----------------------------------------------------------------------------
     Double_t TGeoChecker::Weight(Double_t precision, Option_t *option)
    {
    // Estimate weight of top level volume with a precision SIGMA(W)/W
    // better than PRECISION. Option can be "v" - verbose (default).
       TList *matlist = fGeoManager->GetListOfMaterials();
       Int_t nmat = matlist->GetSize();
       if (!nmat) return 0;
       Int_t *nin = new Int_t[nmat];
       memset(nin, 0, nmat*sizeof(Int_t));
       gRandom = new TRandom3();
       TString opt = option;
       opt.ToLower();
       Bool_t isverbose = opt.Contains("v");
       TGeoBBox *box = (TGeoBBox *)fGeoManager->GetTopVolume()->GetShape();
       Double_t dx = box->GetDX();
       Double_t dy = box->GetDY();
       Double_t dz = box->GetDZ();
       Double_t ox = (box->GetOrigin())[0];
       Double_t oy = (box->GetOrigin())[1];
       Double_t oz = (box->GetOrigin())[2];
       Double_t x,y,z;
       TGeoNode *node;
       TGeoMaterial *mat;
       Double_t vbox = 0.000008*dx*dy*dz; // m3
       Bool_t end = kFALSE;
       Double_t weight=0, sigma, eps, dens;
       Double_t eps0=1.;
       Int_t indmat;
       Int_t igen=0;
       Int_t iin = 0;
       while (!end) {
          x = ox-dx+2*dx*gRandom->Rndm();
          y = oy-dy+2*dy*gRandom->Rndm();
          z = oz-dz+2*dz*gRandom->Rndm();
          node = fGeoManager->FindNode(x,y,z);
          igen++;
          if (!node) continue;
          mat = node->GetVolume()->GetMedium()->GetMaterial();
          indmat = mat->GetIndex();
          if (indmat<0) continue;
          nin[indmat]++;
          iin++;
          if ((iin%100000)==0 || igen>1E8) {
             weight = 0;
             sigma = 0;
             for (indmat=0; indmat<nmat; indmat++) {
                mat = (TGeoMaterial*)matlist->At(indmat);
                dens = mat->GetDensity(); //  [g/cm3]
                if (dens<1E-2) dens=0;
                dens *= 1000.;            // [kg/m3]
                weight += dens*Double_t(nin[indmat]);
                sigma  += dens*dens*nin[indmat];
             }
    	       sigma = TMath::Sqrt(sigma);
    	       eps = sigma/weight;
    	       weight *= vbox/Double_t(igen);
    	       sigma *= vbox/Double_t(igen);
    	       if (eps<precision || igen>1E8) {
    	          if (isverbose) {
    	             printf("=== Weight of %s : %g +/- %g [kg]\n",
    	                    fGeoManager->GetTopVolume()->GetName(), weight, sigma);
                }
    	          end = kTRUE;
    	       } else {
    	          if (isverbose && eps<0.5*eps0) {
    	             printf("%8dK: %14.7g kg  %g %%\n",
                           igen/1000, weight, eps*100);
                   eps0 = eps;
                }
    	       }
          }
       }
       delete [] nin;
       return weight;
    }
    //-----------------------------------------------------------------------------
     Double_t TGeoChecker::CheckVoxels(TGeoVolume *vol, TGeoVoxelFinder *voxels, Double_t *xyz, Int_t npoints)
    {
    // count voxel timing
       TStopwatch timer;
       Double_t time;
       TGeoShape *shape = vol->GetShape();
       TGeoNode *node;
       TGeoMatrix *matrix;
       Double_t *point;
       Double_t local[3];
       Int_t *checklist;
       Int_t ncheck;

       timer.Start();
       for (Int_t i=0; i<npoints; i++) {
          point = xyz + 3*i;
          if (!shape->Contains(point)) continue;
          checklist = voxels->GetCheckList(point, ncheck);
          if (!checklist) continue;
          if (!ncheck) continue;
          for (Int_t id=0; id<ncheck; id++) {
             node = vol->GetNode(checklist[id]);
             matrix = node->GetMatrix();
             matrix->MasterToLocal(point, &local[0]);
             if (node->GetVolume()->GetShape()->Contains(&local[0])) break;
          }
       }
       time = timer.CpuTime();
       return time;
    }
    //-----------------------------------------------------------------------------
     Bool_t TGeoChecker::TestVoxels(TGeoVolume *vol, Int_t npoints)
    {
    // Returns optimal voxelization type for volume vol.
    //   kFALSE - cartesian
    //   kTRUE  - cylindrical
       TGeoVoxelFinder *voxels = vol->GetVoxels();
       if (!voxels) return kFALSE;
       gRandom= new TRandom3();
       const TGeoShape *shape = vol->GetShape();
       Double_t dx = ((TGeoBBox*)shape)->GetDX();
       Double_t dy = ((TGeoBBox*)shape)->GetDY();
       Double_t dz = ((TGeoBBox*)shape)->GetDZ();
       Double_t ox = (((TGeoBBox*)shape)->GetOrigin())[0];
       Double_t oy = (((TGeoBBox*)shape)->GetOrigin())[1];
       Double_t oz = (((TGeoBBox*)shape)->GetOrigin())[2];
       Double_t *xyz = new Double_t[3*npoints];
       Int_t i;
       // generate npoints
       for (i=0; i<npoints; i++) {
          xyz[3*i] = ox-dx+2*dx*gRandom->Rndm();
          xyz[3*i+1] = oy-dy+2*dy*gRandom->Rndm();
          xyz[3*i+2] = oz-dz+2*dz*gRandom->Rndm();
       }
       Bool_t voxtype = vol->IsCylVoxels();
       Double_t time1, time2;
       TGeoVoxelFinder *vox1, *vox2;
       // build both voxelization types
       if (voxtype) {
          printf("   default voxelization was cylindrical.\n");
          vox2 = voxels;
          vox1 = new TGeoVoxelFinder(vol);
          vox1->Voxelize("");
          if (!vol->GetVoxels()) {
             vol->SetVoxelFinder(voxels);
             delete xyz;
             return voxtype;
          }
       } else {
          printf("   default voxelization was cartesian.\n");
          vox1 = voxels;
          vox2 = new TGeoCylVoxels(vol);
          vox2->Voxelize("");
          if (!vol->GetVoxels()) {
             vol->SetVoxelFinder(voxels);
             delete xyz;
             return voxtype;
          }
       }
       // count both voxelization timings
       time1 = CheckVoxels(vol, vox1, xyz, npoints);
       time2 = CheckVoxels(vol, vox2, xyz, npoints);

       printf("   --- time for XYZ : %g\n", time1);
       printf("   --- time for cyl : %g\n", time2);

       if (time1<time2) {
          printf("   best : XYZ\n");
          delete vox2;
          vol->SetVoxelFinder(vox1);
          vol->SetCylVoxels(kFALSE);
          delete xyz;
          return kFALSE;
       }
       printf("   best : cyl\n");
       delete vox1;
       vol->SetVoxelFinder(vox2);
       vol->SetCylVoxels(kTRUE);
       delete xyz;
       return kTRUE;
    }