TGeoChecker.cxx

Go to the documentation of this file.

// @(#)root/geom:$Id$
// 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.             *
 *************************************************************************/
 
/** \class TGeoChecker
\ingroup Geometry_painter
 
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.
 
#### TGeoChecker::CheckPoint(Double_t x, Double_t y, Double_t z)
 
This method can be called directly 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).
 
#### 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.
 
#### 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 until they exit geometry, plotting
all segments crossing visible nodes in the same color as these.
 
#### 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.
 
#### 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.
 
#### TGeoChecker::Weight(Double_t precision)
 
Implementation of TGeoVolume::Weight(). Estimates the total weight of a given
volume by material sampling. Accepts as input the desired precision.
 
#### TGeoChecker::CheckOverlaps(Double_t ovlp, Option_t *option)
Checks the geometry definition for illegal overlaps and illegal extrusions
within a selected geometry hierarchy.
 
An **illegal overlap** is reported when two placed volumes (typically sibling
daughters of the same mother, or a daughter against the mother depending on
the traversal) occupy a common region of space larger than the allowed tolerance.
In practical terms, a set of test points generated on the surface of one
candidate is found **inside** the other candidate volume by more than the configured
tolerance.
 
An **illegal extrusion** is reported when a placed daughter volume has surface
points that lie outside its expected container region (e.g. outside its mother
or outside the allowed placement region as defined by the geometry navigation
context), again beyond the configured tolerance.
 
The overlap checking is performed in three stages:
 
1. Candidate search and filtering
The checker traverses the requested hierarchy and generates candidate pairs.
Candidate pairs are aggressively reduced using fast bounding checks (including
oriented bounding box filters) before any expensive navigation tests are performed.
 
2. Surface point generation and caching
For each distinct TGeoShape appearing in the candidate list, the checker generates
a set of surface sampling points and caches them. The cached sampling depends on
the current meshing policy (see tuning below).
 
3. Exact overlap/extrusion tests (navigation-based)
The final test uses TGeo navigation (Contains, Safety, frame transforms) to evaluate
whether sampled surface points violate the overlap/extrusion constraints. This final
stage can run in parallel when ROOT Implicit MT is enabled.
 
* Usage:
- Check the full geometry
  - Call from the manager:
 
~~~{.cpp}
gGeoManager->CheckOverlaps(ovlp);
~~~
 
- Check only a sub-hierarchy
  - Call from a volume:
 
~~~{.cpp}
gGeoManager->GetTopVolume()->CheckOverlaps(ovlp);
~~~
 
This checks overlaps/extrusions for the selected volume and its daughters
(recursively), without requiring a full-geometry check.
 
* Parameters:
~~~{.cpp}
ovlp
~~~
The allowed overlap tolerance. Violations larger than ovlp are reported
as illegal overlaps/extrusions. The default value is 1.e-1, but a thorough overlap
check may use 1e-6
 
(The option string is intentionally not documented here; legacy sampling modes are
deprecated in favor of dedicated sampling entry points.)
 
* Tuning the surface sampling
 
The quality and cost of the surface sampling used in overlap checking can be tuned
via TGeoManager:
 
~~~{.cpp}
TGeoManager::SetNsegments(Int_t nseg) (default 20)
~~~
Controls the segmentation used by shapes whose surface discretization depends on
angular subdivision (e.g. tubes, cones, polycones, polygonal shapes). Increasing
nseg increases the geometric fidelity of the discretization and usually increases
the time spent in the point-generation and checking stages.
 
~~~{.cpp}
TGeoManager::SetNmeshPoints(Int_t npoints) (default 1000)
~~~
Controls the number of additional surface points requested via GetPointsOnSegments.
Increasing npoints improves sampling density (potentially catching smaller overlaps)
at the cost of more computation in the point-generation and overlap checking stages.
 
Both nseg and npoints affect the cached mesh points. Changing them changes the
sampling policy and triggers regeneration of cached surface points at the next
CheckOverlaps() call.
 
* Parallel execution
 
When ROOT Implicit MT is enabled via:
~~~{.cpp}
ROOT::EnableImplicitMT(N);
~~~
the checker runs the final navigation-based overlap/extrusion evaluation in parallel.
Each worker thread initializes its own navigation state to avoid shared state in
TGeo navigation.
*/
 
#include "TCanvas.h"
#include "TStyle.h"
#include "TFile.h"
#include "TNtuple.h"
#include "TH2.h"
#include "TRandom3.h"
#include "TPolyMarker3D.h"
#include "TPolyLine3D.h"
#include "TStopwatch.h"
 
#include "TGeoVoxelFinder.h"
#include "TGeoBBox.h"
#include "TGeoPcon.h"
#include "TGeoTessellated.h"
#include "TGeoManager.h"
#include "TGeoOverlap.h"
#include "TGeoChecker.h"
 
#include "TBuffer3D.h"
#include "TBuffer3DTypes.h"
#include "TMath.h"
 
#include <cstdlib>
 
// statics and globals
 
////////////////////////////////////////////////////////////////////////////////
/// Default constructor
 
TGeoChecker::TGeoChecker()
   : TVirtualGeoChecker(),
     fGeoManager(nullptr),
     fVsafe(nullptr),
     fFullCheck(kFALSE),
     fVal1(nullptr),
     fVal2(nullptr),
     fFlags(nullptr),
     fTimer(nullptr),
     fSelectedNode(nullptr),
     fNchecks(0)
{
}
 
////////////////////////////////////////////////////////////////////////////////
/// Constructor for a given geometry
 
TGeoChecker::TGeoChecker(TGeoManager *geom)
   : TVirtualGeoChecker(),
     fGeoManager(geom),
     fVsafe(nullptr),
     fFullCheck(kFALSE),
     fVal1(nullptr),
     fVal2(nullptr),
     fFlags(nullptr),
     fTimer(nullptr),
     fSelectedNode(nullptr),
     fNchecks(0)
{
}
 
////////////////////////////////////////////////////////////////////////////////
/// Destructor
 
TGeoChecker::~TGeoChecker()
{
   if (fTimer)
      delete fTimer;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Print current operation progress.
 
void TGeoChecker::OpProgress(const char *opname, Long64_t current, Long64_t size, TStopwatch *watch, Bool_t last,
                             Bool_t refresh, const char *msg)
{
   static Long64_t icount = 0;
   static TString oname;
   static TString nname;
   static Long64_t ocurrent = 0;
   static Long64_t osize = 0;
   static Int_t oseconds = 0;
   static TStopwatch *owatch = nullptr;
   static Bool_t oneoftwo = kFALSE;
   static Int_t nrefresh = 0;
   const char symbol[4] = {'=', '\\', '|', '/'};
   char progress[11] = "          ";
   Int_t ichar = icount % 4;
   TString message(msg);
   message += "         ";
 
   if (!refresh) {
      nrefresh = 0;
      if (!size)
         return;
      owatch = watch;
      oname = opname;
      ocurrent = TMath::Abs(current);
      osize = TMath::Abs(size);
      if (ocurrent > osize)
         ocurrent = osize;
   } else {
      nrefresh++;
      if (!osize)
         return;
   }
   icount++;
   Double_t time = 0.;
   Int_t hours = 0;
   Int_t minutes = 0;
   Int_t seconds = 0;
   if (owatch && !last) {
      owatch->Stop();
      time = owatch->RealTime();
      hours = (Int_t)(time / 3600.);
      time -= 3600 * hours;
      minutes = (Int_t)(time / 60.);
      time -= 60 * minutes;
      seconds = (Int_t)time;
      if (refresh) {
         if (oseconds == seconds) {
            owatch->Continue();
            return;
         }
         oneoftwo = !oneoftwo;
      }
      oseconds = seconds;
   }
   if (refresh && oneoftwo) {
      nname = oname;
      if (fNchecks <= nrefresh)
         fNchecks = nrefresh + 1;
      Int_t pctdone = (Int_t)(100. * nrefresh / fNchecks);
      oname = TString::Format("     == %3d%% ==", pctdone);
   }
   Double_t percent = 100.0 * ocurrent / osize;
   Int_t nchar = Int_t(percent / 10);
   if (nchar > 10)
      nchar = 10;
   Int_t i;
   for (i = 0; i < nchar; i++)
      progress[i] = '=';
   progress[nchar] = symbol[ichar];
   for (i = nchar + 1; i < 10; i++)
      progress[i] = ' ';
   progress[10] = '\0';
   oname += "                    ";
   oname.Remove(20);
   if (size < 10000)
      fprintf(stderr, "%s [%10s] %4lld ", oname.Data(), progress, ocurrent);
   else if (size < 100000)
      fprintf(stderr, "%s [%10s] %5lld ", oname.Data(), progress, ocurrent);
   else
      fprintf(stderr, "%s [%10s] %7lld ", oname.Data(), progress, ocurrent);
   if (time > 0.)
      fprintf(stderr, "[%6.2f %%]   TIME %.2d:%.2d:%.2d  %s\r", percent, hours, minutes, seconds, message.Data());
   else
      fprintf(stderr, "[%6.2f %%]  %s\r", percent, message.Data());
   if (refresh && oneoftwo)
      oname = nname;
   if (owatch)
      owatch->Continue();
   if (last) {
      icount = 0;
      owatch = nullptr;
      ocurrent = 0;
      osize = 0;
      oseconds = 0;
      oneoftwo = kFALSE;
      nrefresh = 0;
      fprintf(stderr, "\n");
   }
}
 
////////////////////////////////////////////////////////////////////////////////
/// Check pushes and pulls needed to cross the next boundary with respect to the
/// position given by FindNextBoundary. If radius is not mentioned the full bounding
/// box will be sampled.
 
void TGeoChecker::CheckBoundaryErrors(Int_t ntracks, Double_t radius)
{
   TGeoVolume *tvol = fGeoManager->GetTopVolume();
   Info("CheckBoundaryErrors", "Top volume is %s", tvol->GetName());
   const TGeoShape *shape = tvol->GetShape();
   TGeoBBox *box = (TGeoBBox *)shape;
   Double_t dl[3];
   Double_t ori[3];
   Double_t xyz[3];
   Double_t nxyz[3];
   Double_t dir[3];
   Double_t relp;
   Char_t path[1024];
   Char_t cdir[10];
 
   // Tree part
   TFile *f = new TFile("geobugs.root", "recreate");
   TTree *bug = new TTree("bug", "Geometrical problems");
   bug->Branch("pos", xyz, "xyz[3]/D");
   bug->Branch("dir", dir, "dir[3]/D");
   bug->Branch("push", &relp, "push/D");
   bug->Branch("path", &path, "path/C");
   bug->Branch("cdir", &cdir, "cdir/C");
 
   dl[0] = box->GetDX();
   dl[1] = box->GetDY();
   dl[2] = box->GetDZ();
   ori[0] = (box->GetOrigin())[0];
   ori[1] = (box->GetOrigin())[1];
   ori[2] = (box->GetOrigin())[2];
   if (radius > 0)
      dl[0] = dl[1] = dl[2] = radius;
 
   TDirectory::TContext ctx{nullptr}; // Disable directory registration for histograms
   TH1F *hnew = new TH1F("hnew", "Precision pushing", 30, -20., 10.);
   TH1F *hold = new TH1F("hold", "Precision pulling", 30, -20., 10.);
   TH2F *hplotS = new TH2F("hplotS", "Problematic points", 100, -dl[0], dl[0], 100, -dl[1], dl[1]);
   gStyle->SetOptStat(111111);
 
   TGeoNode *node = nullptr;
   Long_t igen = 0;
   Long_t itry = 0;
   Long_t n100 = ntracks / 100;
   Double_t rad = TMath::Sqrt(dl[0] * dl[0] + dl[1] * dl[1]);
   printf("Random box : %f, %f, %f, %f, %f, %f\n", ori[0], ori[1], ori[2], dl[0], dl[1], dl[2]);
   printf("Start... %i points\n", ntracks);
   if (!fTimer)
      fTimer = new TStopwatch();
   fTimer->Reset();
   fTimer->Start();
   while (igen < ntracks) {
      Double_t phi1 = TMath::TwoPi() * gRandom->Rndm();
      Double_t r = rad * gRandom->Rndm();
      xyz[0] = ori[0] + r * TMath::Cos(phi1);
      xyz[1] = ori[1] + r * TMath::Sin(phi1);
      Double_t z = (1. - 2. * gRandom->Rndm());
      xyz[2] = ori[2] + dl[2] * z * TMath::Abs(z);
      ++itry;
      fGeoManager->SetCurrentPoint(xyz);
      node = fGeoManager->FindNode();
      if (!node || node == fGeoManager->GetTopNode())
         continue;
      ++igen;
      if (n100 && !(igen % n100))
         OpProgress("Sampling progress:", igen, ntracks, fTimer);
      Double_t cost = 1. - 2. * gRandom->Rndm();
      Double_t sint = TMath::Sqrt((1. + cost) * (1. - cost));
      Double_t phi = TMath::TwoPi() * gRandom->Rndm();
      dir[0] = sint * TMath::Cos(phi);
      dir[1] = sint * TMath::Sin(phi);
      dir[2] = cost;
      fGeoManager->SetCurrentDirection(dir);
      fGeoManager->FindNextBoundary();
      Double_t step = fGeoManager->GetStep();
 
      relp = 1.e-21;
      for (Int_t i = 0; i < 30; ++i) {
         relp *= 10.;
         for (Int_t j = 0; j < 3; ++j)
            nxyz[j] = xyz[j] + step * (1. + relp) * dir[j];
         if (!fGeoManager->IsSameLocation(nxyz[0], nxyz[1], nxyz[2])) {
            hnew->Fill(i - 20.);
            if (i > 15) {
               const Double_t *norm = fGeoManager->FindNormal();
               strncpy(path, fGeoManager->GetPath(), 1024);
               path[1023] = '\0';
               Double_t dotp = norm[0] * dir[0] + norm[1] * dir[1] + norm[2] * dir[2];
               printf("Forward error i=%d p=%5.4f %5.4f %5.4f s=%5.4f dot=%5.4f path=%s\n", i, xyz[0], xyz[1], xyz[2],
                      step, dotp, path);
               hplotS->Fill(xyz[0], xyz[1], (Double_t)i);
               strncpy(cdir, "Forward", 10);
               bug->Fill();
            }
            break;
         }
      }
 
      relp = -1.e-21;
      for (Int_t i = 0; i < 30; ++i) {
         relp *= 10.;
         for (Int_t j = 0; j < 3; ++j)
            nxyz[j] = xyz[j] + step * (1. + relp) * dir[j];
         if (fGeoManager->IsSameLocation(nxyz[0], nxyz[1], nxyz[2])) {
            hold->Fill(i - 20.);
            if (i > 15) {
               const Double_t *norm = fGeoManager->FindNormal();
               strncpy(path, fGeoManager->GetPath(), 1024);
               path[1023] = '\0';
               Double_t dotp = norm[0] * dir[0] + norm[1] * dir[1] + norm[2] * dir[2];
               printf("Backward error i=%d p=%5.4f %5.4f %5.4f s=%5.4f dot=%5.4f path=%s\n", i, xyz[0], xyz[1], xyz[2],
                      step, dotp, path);
               strncpy(cdir, "Backward", 10);
               bug->Fill();
            }
            break;
         }
      }
   }
   fTimer->Stop();
 
   if (itry)
      printf("CPU time/point = %5.2emus: Real time/point = %5.2emus\n", 1000000. * fTimer->CpuTime() / itry,
             1000000. * fTimer->RealTime() / itry);
   bug->Write();
   delete bug;
   bug = nullptr;
   delete f;
 
   CheckBoundaryReference();
 
   if (itry)
      printf("Effic = %3.1f%%\n", (100. * igen) / itry);
   TCanvas *c1 = new TCanvas("c1", "Results", 600, 800);
   c1->Divide(1, 2);
   c1->cd(1);
   gPad->SetLogy();
   hold->Draw();
   c1->cd(2);
   gPad->SetLogy();
   hnew->Draw();
   /*TCanvas *c3 = */ new TCanvas("c3", "Plot", 600, 600);
   hplotS->Draw("cont0");
}
 
////////////////////////////////////////////////////////////////////////////////
/// Check the boundary errors reference file created by CheckBoundaryErrors method.
/// The shape for which the crossing failed is drawn with the starting point in red
/// and the extrapolated point to boundary (+/- failing push/pull) in yellow.
 
void TGeoChecker::CheckBoundaryReference(Int_t icheck)
{
   Double_t xyz[3];
   Double_t nxyz[3];
   Double_t dir[3];
   Double_t lnext[3];
   Double_t push;
   Char_t path[1024];
   Char_t cdir[10];
   // Tree part
   TFile *f = new TFile("geobugs.root", "read");
   TTree *bug = (TTree *)f->Get("bug");
   bug->SetBranchAddress("pos", xyz);
   bug->SetBranchAddress("dir", dir);
   bug->SetBranchAddress("push", &push);
   bug->SetBranchAddress("path", &path);
   bug->SetBranchAddress("cdir", &cdir);
 
   Int_t nentries = (Int_t)bug->GetEntries();
   printf("nentries %d\n", nentries);
   if (icheck < 0) {
      for (Int_t i = 0; i < nentries; i++) {
         bug->GetEntry(i);
         printf("%-9s error push=%g p=%5.4f %5.4f %5.4f s=%5.4f dot=%5.4f path=%s\n", cdir, push, xyz[0], xyz[1],
                xyz[2], 1., 1., path);
      }
   } else {
      if (icheck >= nentries)
         return;
      Int_t idebug = TGeoManager::GetVerboseLevel();
      TGeoManager::SetVerboseLevel(5);
      bug->GetEntry(icheck);
      printf("%-9s error push=%g p=%5.4f %5.4f %5.4f s=%5.4f dot=%5.4f path=%s\n", cdir, push, xyz[0], xyz[1], xyz[2],
             1., 1., path);
      fGeoManager->SetCurrentPoint(xyz);
      fGeoManager->SetCurrentDirection(dir);
      fGeoManager->FindNode();
      TGeoNode *next = fGeoManager->FindNextBoundary();
      Double_t step = fGeoManager->GetStep();
      for (Int_t j = 0; j < 3; j++)
         nxyz[j] = xyz[j] + step * (1. + 0.1 * push) * dir[j];
      Bool_t change = !fGeoManager->IsSameLocation(nxyz[0], nxyz[1], nxyz[2]);
      printf("step=%g in: %s\n", step, fGeoManager->GetPath());
      printf("  -> next = %s push=%g  change=%d\n", next->GetName(), push, (UInt_t)change);
      next->GetVolume()->InspectShape();
      next->GetVolume()->DrawOnly();
      if (next != fGeoManager->GetCurrentNode()) {
         Int_t index1 = fGeoManager->GetCurrentVolume()->GetIndex(next);
         if (index1 >= 0)
            fGeoManager->CdDown(index1);
      }
      TPolyMarker3D *pm = new TPolyMarker3D();
      fGeoManager->MasterToLocal(xyz, lnext);
      pm->SetNextPoint(lnext[0], lnext[1], lnext[2]);
      pm->SetMarkerStyle(2);
      pm->SetMarkerSize(0.2);
      pm->SetMarkerColor(kRed);
      pm->Draw("SAME");
      TPolyMarker3D *pm1 = new TPolyMarker3D();
      for (Int_t j = 0; j < 3; j++)
         nxyz[j] = xyz[j] + step * dir[j];
      fGeoManager->MasterToLocal(nxyz, lnext);
      pm1->SetNextPoint(lnext[0], lnext[1], lnext[2]);
      pm1->SetMarkerStyle(2);
      pm1->SetMarkerSize(0.2);
      pm1->SetMarkerColor(kYellow);
      pm1->Draw("SAME");
      TGeoManager::SetVerboseLevel(idebug);
   }
   delete bug;
   delete f;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Geometry checking. Optional overlap checkings (by sampling and by mesh). Optional
/// boundary crossing check + timing per volume.
///
/// STAGE 1: extensive overlap checking by sampling per volume. Stdout need to be
///  checked by user to get report, then TGeoVolume::CheckOverlaps(0.01, "s") can
///  be called for the suspicious volumes.
///
/// STAGE 2: normal overlap checking using the shapes mesh - fills the list of
///  overlaps.
///
/// STAGE 3: shooting NRAYS rays from VERTEX and counting the total number of
///  crossings per volume (rays propagated from boundary to boundary until
///  geometry exit). Timing computed and results stored in a histo.
///
/// STAGE 4: shooting 1 mil. random rays inside EACH volume and calling
///  FindNextBoundary() + Safety() for each call. The timing is normalized by the
///  number of crossings computed at stage 2 and presented as percentage.
///  One can get a picture on which are the most "burned" volumes during
///  transportation from geometry point of view. Another plot of the timing per
///  volume vs. number of daughters is produced.
///
/// All histos are saved in the file statistics.root
 
void TGeoChecker::CheckGeometryFull(Bool_t checkoverlaps, Bool_t checkcrossings, Int_t ntracks, const Double_t *vertex)
{
   Int_t nuid = fGeoManager->GetListOfUVolumes()->GetEntries();
   if (!fTimer)
      fTimer = new TStopwatch();
   Int_t i;
   Double_t value;
   fFlags = new Bool_t[nuid];
   memset(fFlags, 0, nuid * sizeof(Bool_t));
   TGeoVolume *vol;
   TCanvas *c = new TCanvas("overlaps", "Overlaps by sampling", 800, 800);
 
   // STAGE 1: Overlap checking by sampling per volume
   if (checkoverlaps) {
      printf("====================================================================\n");
      printf("STAGE 1: Overlap checking by sampling within 10 microns\n");
      printf("====================================================================\n");
      fGeoManager->CheckOverlaps(0.001, "s");
 
      // STAGE 2: Global overlap/extrusion checking
      printf("====================================================================\n");
      printf("STAGE 2: Global overlap/extrusion checking within 10 microns\n");
      printf("====================================================================\n");
      fGeoManager->CheckOverlaps(0.001);
   }
 
   if (!checkcrossings) {
      delete[] fFlags;
      fFlags = nullptr;
      delete c;
      return;
   }
 
   fVal1 = new Double_t[nuid];
   fVal2 = new Double_t[nuid];
   memset(fVal1, 0, nuid * sizeof(Double_t));
   memset(fVal2, 0, nuid * sizeof(Double_t));
   // STAGE 3: How many crossings per volume in a realistic case ?
   // Ignore volumes with no daughters
 
   // Generate rays from vertex in phi=[0,2*pi] theta=[0,pi]
   //   Int_t ntracks = 1000000;
   printf("====================================================================\n");
   printf("STAGE 3: Propagating %i tracks starting from vertex\n and counting number of boundary crossings...\n",
          ntracks);
   printf("====================================================================\n");
   Int_t nbound = 0;
   Double_t theta, phi;
   Double_t point[3], dir[3];
   memset(point, 0, 3 * sizeof(Double_t));
   if (vertex)
      memcpy(point, vertex, 3 * sizeof(Double_t));
 
   fTimer->Reset();
   fTimer->Start();
   for (i = 0; i < ntracks; i++) {
      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);
      if ((i % 100) == 0)
         OpProgress("Transporting tracks", i, ntracks, fTimer);
      //      if ((i%1000)==0) printf("... remaining tracks %i\n", ntracks-i);
      nbound += PropagateInGeom(point, dir);
   }
   if (!nbound) {
      printf("No boundary crossed\n");
      return;
   }
   fTimer->Stop();
   Double_t time1 = fTimer->CpuTime() * 1.E6;
   Double_t time2 = time1 / ntracks;
   Double_t time3 = time1 / nbound;
   OpProgress("Transporting tracks", ntracks, ntracks, fTimer, kTRUE);
   printf("Time for crossing %i boundaries: %g [ms]\n", nbound, time1);
   printf("Time per track for full geometry traversal: %g [ms], per crossing: %g [ms]\n", time2, time3);
 
   // STAGE 4: How much time per volume:
 
   printf("====================================================================\n");
   printf("STAGE 4: How much navigation time per volume per next+safety call\n");
   printf("====================================================================\n");
   TGeoIterator next(fGeoManager->GetTopVolume());
   TGeoNode *current;
   TString path;
   vol = fGeoManager->GetTopVolume();
   memset(fFlags, 0, nuid * sizeof(Bool_t));
   TStopwatch timer;
   timer.Start();
   i = 0;
   char volname[30];
   strncpy(volname, vol->GetName(), 15);
   volname[15] = '\0';
   OpProgress(volname, i++, nuid, &timer);
   Score(vol, 1, TimingPerVolume(vol));
   while ((current = next())) {
      vol = current->GetVolume();
      Int_t uid = vol->GetNumber();
      if (fFlags[uid])
         continue;
      fFlags[uid] = kTRUE;
      next.GetPath(path);
      fGeoManager->cd(path.Data());
      strncpy(volname, vol->GetName(), 15);
      volname[15] = '\0';
      OpProgress(volname, i++, nuid, &timer);
      Score(vol, 1, TimingPerVolume(vol));
   }
   OpProgress("STAGE 4 completed", i, nuid, &timer, kTRUE);
   // Draw some histos
   Double_t time_tot_pertrack = 0.;
   TCanvas *c1 = new TCanvas("c2", "ncrossings", 10, 10, 900, 500);
   c1->SetGrid();
   c1->SetTopMargin(0.15);
   TFile *f = new TFile("statistics.root", "RECREATE");
   TH1F *h = new TH1F("h", "number of boundary crossings per volume", 3, 0, 3);
   h->SetStats(false);
   h->SetFillColor(38);
   h->SetCanExtend(TH1::kAllAxes);
 
   memset(fFlags, 0, nuid * sizeof(Bool_t));
   for (i = 0; i < nuid; i++) {
      vol = fGeoManager->GetVolume(i);
      if (!vol->GetNdaughters())
         continue;
      time_tot_pertrack += fVal1[i] * fVal2[i];
      h->Fill(vol->GetName(), (Int_t)fVal1[i]);
   }
   time_tot_pertrack /= ntracks;
   h->LabelsDeflate();
   h->LabelsOption(">", "X");
   h->Draw();
 
   TCanvas *c2 = new TCanvas("c3", "time spent per volume in navigation", 10, 10, 900, 500);
   c2->SetGrid();
   c2->SetTopMargin(0.15);
   TH2F *h2 = new TH2F("h2", "time per FNB call vs. ndaughters", 100, 0, 100, 100, 0, 15);
   h2->SetStats(false);
   h2->SetMarkerStyle(2);
   TH1F *h1 = new TH1F("h1", "percent of time spent per volume", 3, 0, 3);
   h1->SetStats(false);
   h1->SetFillColor(38);
   h1->SetCanExtend(TH1::kAllAxes);
   for (i = 0; i < nuid; i++) {
      vol = fGeoManager->GetVolume(i);
      if (!vol || !vol->GetNdaughters())
         continue;
      value = fVal1[i] * fVal2[i] / ntracks / time_tot_pertrack;
      h1->Fill(vol->GetName(), value);
      h2->Fill(vol->GetNdaughters(), fVal2[i]);
   }
   h1->LabelsDeflate();
   h1->LabelsOption(">", "X");
   h1->Draw();
   TCanvas *c3 = new TCanvas("c4", "timing vs. ndaughters", 10, 10, 900, 500);
   c3->SetGrid();
   c3->SetTopMargin(0.15);
   h2->Draw();
   f->Write();
   delete[] fFlags;
   fFlags = nullptr;
   delete[] fVal1;
   fVal1 = nullptr;
   delete[] fVal2;
   fVal2 = nullptr;
   delete fTimer;
   fTimer = nullptr;
   delete c;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Propagate from START along DIR from boundary to boundary until exiting
/// geometry. Fill array of hits.
 
Int_t TGeoChecker::PropagateInGeom(Double_t *start, Double_t *dir)
{
   fGeoManager->InitTrack(start, dir);
   TGeoNode *current = nullptr;
   Int_t nzero = 0;
   Int_t nhits = 0;
   while (!fGeoManager->IsOutside()) {
      if (nzero > 3) {
         // Problems in trying to cross a boundary
         printf("Error in trying to cross boundary of %s\n", current->GetName());
         return nhits;
      }
      current = fGeoManager->FindNextBoundaryAndStep(TGeoShape::Big(), kFALSE);
      if (!current || fGeoManager->IsOutside())
         return nhits;
      Double_t step = fGeoManager->GetStep();
      if (step < 2. * TGeoShape::Tolerance()) {
         nzero++;
         continue;
      } else
         nzero = 0;
      // Generate the hit
      nhits++;
      TGeoVolume *vol = current->GetVolume();
      Score(vol, 0, 1.);
      Int_t iup = 1;
      TGeoNode *mother = fGeoManager->GetMother(iup++);
      while (mother && mother->GetVolume()->IsAssembly()) {
         Score(mother->GetVolume(), 0, 1.);
         mother = fGeoManager->GetMother(iup++);
      }
   }
   return nhits;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Score a hit for VOL
 
void TGeoChecker::Score(TGeoVolume *vol, Int_t ifield, Double_t value)
{
   Int_t uid = vol->GetNumber();
   switch (ifield) {
   case 0: fVal1[uid] += value; break;
   case 1: fVal2[uid] += value;
   }
}
 
////////////////////////////////////////////////////////////////////////////////
/// Compute timing per "FindNextBoundary" + "Safety" call. Volume must be
/// in the current path.
 
Double_t TGeoChecker::TimingPerVolume(TGeoVolume *vol)
{
   fTimer->Reset();
   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 point[3], dir[3], lpt[3], ldir[3];
   Double_t pstep = 0.;
   pstep = TMath::Max(pstep, dz);
   Double_t theta, phi;
   Int_t idaughter;
   fTimer->Start();
   for (Int_t i = 0; i < 1000000; i++) {
      lpt[0] = ox - dx + 2 * dx * gRandom->Rndm();
      lpt[1] = oy - dy + 2 * dy * gRandom->Rndm();
      lpt[2] = oz - dz + 2 * dz * gRandom->Rndm();
      fGeoManager->GetCurrentMatrix()->LocalToMaster(lpt, point);
      fGeoManager->SetCurrentPoint(point[0], point[1], point[2]);
      phi = 2 * TMath::Pi() * gRandom->Rndm();
      theta = TMath::ACos(1. - 2. * gRandom->Rndm());
      ldir[0] = TMath::Sin(theta) * TMath::Cos(phi);
      ldir[1] = TMath::Sin(theta) * TMath::Sin(phi);
      ldir[2] = TMath::Cos(theta);
      fGeoManager->GetCurrentMatrix()->LocalToMasterVect(ldir, dir);
      fGeoManager->SetCurrentDirection(dir);
      fGeoManager->SetStep(pstep);
      fGeoManager->ResetState();
      // dist = TGeoShape::Big();
      if (!vol->IsAssembly()) {
         Bool_t inside = vol->Contains(lpt);
         if (!inside) {
            // dist = vol->GetShape()->DistFromOutside(lpt,ldir,3,pstep);
            // if (dist>=pstep) continue;
         } else {
            vol->GetShape()->DistFromInside(lpt, ldir, 3, pstep);
         }
 
         if (!vol->GetNdaughters())
            vol->GetShape()->Safety(lpt, inside);
      }
      if (vol->GetNdaughters()) {
         fGeoManager->Safety();
         fGeoManager->FindNextDaughterBoundary(point, dir, idaughter, kFALSE);
      }
   }
   fTimer->Stop();
   Double_t time_per_track = fTimer->CpuTime();
   Int_t uid = vol->GetNumber();
   Int_t ncrossings = (Int_t)fVal1[uid];
   if (!vol->GetNdaughters())
      printf("Time for volume %s (shape=%s): %g [ms] ndaughters=%d ncross=%d\n", vol->GetName(),
             vol->GetShape()->GetName(), time_per_track, vol->GetNdaughters(), ncrossings);
   else
      printf("Time for volume %s (assemb=%d): %g [ms] ndaughters=%d ncross=%d\n", vol->GetName(), vol->IsAssembly(),
             time_per_track, vol->GetNdaughters(), ncrossings);
   return time_per_track;
}
 
////////////////////////////////////////////////////////////////////////////////
/// 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.
 
void TGeoChecker::CheckGeometry(Int_t nrays, Double_t startx, Double_t starty, Double_t startz) const
{
   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);
      array1 = 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);
      array2 = 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;
   delete pma; // markers are drawn on the pad
}
 
////////////////////////////////////////////////////////////////////////////////
/// Clean-up the mesh of pcon/pgon from useless points
 
void TGeoChecker::CleanPoints(Double_t *points, Int_t &numPoints) const
{
   Int_t ipoint = 0;
   Int_t j, k = 0;
   Double_t rsq;
   for (Int_t i = 0; i < numPoints; i++) {
      j = 3 * i;
      rsq = points[j] * points[j] + points[j + 1] * points[j + 1];
      if (rsq < 1.e-10)
         continue;
      points[k] = points[j];
      points[k + 1] = points[j + 1];
      points[k + 2] = points[j + 2];
      ipoint++;
      k = 3 * ipoint;
   }
   numPoints = ipoint;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Compute overlap/extrusion for given candidate using mesh points of the shapes.
 
Bool_t TGeoChecker::ComputeOverlap(const TGeoOverlapCandidate &c, TGeoOverlapResult &out) const
{
   out = TGeoOverlapResult{}; // reset
 
   TGeoVolume *vol1 = c.fVol1;
   TGeoVolume *vol2 = c.fVol2;
   if (!vol1 || !vol2)
      return kFALSE;
 
   // Keep legacy early outs
   if (vol1->IsAssembly() || vol2->IsAssembly())
      return kFALSE;
 
   if (dynamic_cast<TGeoTessellated *>(vol1->GetShape()) || dynamic_cast<TGeoTessellated *>(vol2->GetShape()))
      return kFALSE;
 
   TGeoShape *shape1 = vol1->GetShape();
   TGeoShape *shape2 = vol2->GetShape();
 
   const auto &m1 = GetMeshPoints(shape1);
   const auto &m2 = GetMeshPoints(shape2);
   if (m1.fNPoints <= 0 || m2.fNPoints <= 0)
      return kFALSE;
 
   const Double_t *points1 = m1.fPoints.data();
   const Double_t *points2 = m2.fPoints.data();
   const Int_t numPoints1 = m1.fNPoints;
   const Int_t numPoints2 = m2.fNPoints;
 
   const TGeoHMatrix mat1 = c.fMat1;
   const TGeoHMatrix mat2 = c.fMat2;
 
   const Double_t tolerance = TGeoShape::Tolerance();
   const Double_t ovlp = c.fOvlp;
 
   auto add_point = [&](const Double_t p[3]) {
      if (out.fPoints.size() < 100)
         out.fPoints.push_back({{p[0], p[1], p[2]}});
   };
 
   auto start_result_if_needed = [&]() {
      if (out.fVol1)
         return; // already started
      out.fName = c.fName;
      out.fVol1 = vol1;
      out.fVol2 = vol2;
      out.fMat1 = mat1;
      out.fMat2 = mat2;
      out.fIsOverlap = c.fIsOverlap;
      out.fMaxOverlap = 0.0;
   };
 
   // Scratch
   Double_t local[3], local1[3], point[3];
   Double_t safety = TGeoShape::Big();
 
   // ---------------------------
   // Extrusion case (c.fIsOverlap == false):
   // "Test vol2 extrudes vol1" and also the legacy “mother points inside daughter” pass.
   // ---------------------------
   if (!c.fIsOverlap) {
      Bool_t found = kFALSE;
 
      // loop all points of the daughter (vol2)
      for (Int_t ip = 0; ip < numPoints2; ++ip) {
         memcpy(local, &points2[3 * ip], 3 * sizeof(Double_t));
         if (TMath::Abs(local[0]) < tolerance && TMath::Abs(local[1]) < tolerance)
            continue;
 
         mat2.LocalToMaster(local, point);
         mat1.MasterToLocal(point, local);
 
         Bool_t extrude = !shape1->Contains(local);
         if (extrude) {
            safety = shape1->Safety(local, kFALSE);
            if (safety < ovlp)
               extrude = kFALSE;
         }
 
         if (extrude) {
            start_result_if_needed();
            found = kTRUE;
            if (safety > out.fMaxOverlap)
               out.fMaxOverlap = safety;
            add_point(point);
         }
      }
 
      // loop all points of the mother (vol1)
      for (Int_t ip = 0; ip < numPoints1; ++ip) {
         memcpy(local, &points1[3 * ip], 3 * sizeof(Double_t));
         if (local[0] < 1e-10 && local[1] < 1e-10)
            continue;
 
         mat1.LocalToMaster(local, point);
         mat2.MasterToLocal(point, local1);
 
         Bool_t extrude = shape2->Contains(local1);
         if (extrude) {
            safety = shape1->Safety(local, kTRUE);
            if (safety > 1E-6) {
               extrude = kFALSE;
            } else {
               safety = shape2->Safety(local1, kTRUE);
               if (safety < ovlp)
                  extrude = kFALSE;
            }
         }
 
         if (extrude) {
            start_result_if_needed();
            found = kTRUE;
            if (safety > out.fMaxOverlap)
               out.fMaxOverlap = safety;
            add_point(point);
         }
      }
 
      return found;
   }
 
   // ---------------------------
   // Overlap case
   // ---------------------------
   Bool_t found = kFALSE;
 
   // loop all points of first candidate (vol1) in vol2
   for (Int_t ip = 0; ip < numPoints1; ++ip) {
      memcpy(local, &points1[3 * ip], 3 * sizeof(Double_t));
      if (local[0] < 1e-10 && local[1] < 1e-10)
         continue;
 
      mat1.LocalToMaster(local, point);
      mat2.MasterToLocal(point, local);
 
      Bool_t overlap = shape2->Contains(local);
      if (overlap) {
         safety = shape2->Safety(local, kTRUE);
         if (safety < ovlp)
            overlap = kFALSE;
      }
 
      if (overlap) {
         start_result_if_needed();
         found = kTRUE;
         if (safety > out.fMaxOverlap)
            out.fMaxOverlap = safety;
         add_point(point);
      }
   }
 
   // loop all points of second candidate (vol2) in vol1
   for (Int_t ip = 0; ip < numPoints2; ++ip) {
      memcpy(local, &points2[3 * ip], 3 * sizeof(Double_t));
      if (local[0] < 1e-10 && local[1] < 1e-10)
         continue;
 
      mat2.LocalToMaster(local, point);
      mat1.MasterToLocal(point, local);
 
      Bool_t overlap = shape1->Contains(local);
      if (overlap) {
         safety = shape1->Safety(local, kTRUE);
         if (safety < ovlp)
            overlap = kFALSE;
      }
 
      if (overlap) {
         start_result_if_needed();
         found = kTRUE;
         if (safety > out.fMaxOverlap)
            out.fMaxOverlap = safety;
         add_point(point);
      }
   }
 
   return found;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Check illegal overlaps for volume VOL within a limit OVLP by sampling npoints
/// inside the volume shape.
 
void TGeoChecker::CheckOverlapsBySampling(const TGeoVolume *vol, Double_t ovlp, Int_t npoints) const
{
   Int_t nd = vol->GetNdaughters();
   if (nd < 2)
      return;
   TGeoVoxelFinder *voxels = vol->GetVoxels();
   if (!voxels)
      return;
   TGeoBBox *box = (TGeoBBox *)vol->GetShape();
   TGeoShape *shape;
   TGeoNode *node;
   Double_t dx = box->GetDX();
   Double_t dy = box->GetDY();
   Double_t dz = box->GetDZ();
   Double_t pt[3];
   Double_t local[3];
   Int_t *check_list = nullptr;
   Int_t ncheck = 0;
   const Double_t *orig = box->GetOrigin();
   Int_t ipoint = 0;
   Int_t itry = 0;
   Int_t iovlp = 0;
   Int_t id = 0, id0 = 0, id1 = 0;
   Bool_t in, incrt;
   Double_t safe;
   TString name1 = "";
   TString name2 = "";
   TGeoOverlap **flags = nullptr;
   TGeoNode *node1, *node2;
   Int_t novlps = 0;
   TGeoHMatrix mat1, mat2;
   //   Int_t tid = TGeoManager::ThreadId();
   TGeoNavigator *nav = fGeoManager->GetCurrentNavigator();
   TGeoStateInfo &td = *nav->GetCache()->GetInfo();
   while (ipoint < npoints) {
      // Shoot randomly in the bounding box.
      pt[0] = orig[0] - dx + 2. * dx * gRandom->Rndm();
      pt[1] = orig[1] - dy + 2. * dy * gRandom->Rndm();
      pt[2] = orig[2] - dz + 2. * dz * gRandom->Rndm();
      if (!vol->Contains(pt)) {
         itry++;
         if (itry > 10000 && !ipoint) {
            Error("CheckOverlapsBySampling", "No point inside volume!!! - aborting");
            break;
         }
         continue;
      }
      // Check if the point is inside one or more daughters
      in = kFALSE;
      ipoint++;
      check_list = voxels->GetCheckList(pt, ncheck, td);
      if (!check_list || ncheck < 2)
         continue;
      for (id = 0; id < ncheck; id++) {
         id0 = check_list[id];
         node = vol->GetNode(id0);
         // Ignore MANY's
         if (node->IsOverlapping())
            continue;
         node->GetMatrix()->MasterToLocal(pt, local);
         shape = node->GetVolume()->GetShape();
         incrt = shape->Contains(local);
         if (!incrt)
            continue;
         if (!in) {
            in = kTRUE;
            id1 = id0;
            continue;
         }
         // The point is inside 2 or more daughters, check safety
         safe = shape->Safety(local, kTRUE);
         if (safe < ovlp)
            continue;
         // We really have found an overlap -> store the point in a container
         iovlp++;
         if (!novlps) {
            if (flags)
               delete[] flags; // should never happen
            flags = new TGeoOverlap *[nd * nd];
            memset(flags, 0, nd * nd * sizeof(TGeoOverlap *));
         }
         TGeoOverlap *nodeovlp = flags[nd * id1 + id0];
         if (!nodeovlp) {
            novlps++;
            node1 = vol->GetNode(id1);
            name1 = node1->GetName();
            mat1 = node1->GetMatrix();
            Int_t cindex = node1->GetVolume()->GetCurrentNodeIndex();
            while (cindex >= 0) {
               node1 = node1->GetVolume()->GetNode(cindex);
               name1 += TString::Format("/%s", node1->GetName());
               mat1.Multiply(node1->GetMatrix());
               cindex = node1->GetVolume()->GetCurrentNodeIndex();
            }
            node2 = vol->GetNode(id0);
            name2 = node2->GetName();
            mat2 = node2->GetMatrix();
            cindex = node2->GetVolume()->GetCurrentNodeIndex();
            while (cindex >= 0) {
               node2 = node2->GetVolume()->GetNode(cindex);
               name2 += TString::Format("/%s", node2->GetName());
               mat2.Multiply(node2->GetMatrix());
               cindex = node2->GetVolume()->GetCurrentNodeIndex();
            }
            nodeovlp = new TGeoOverlap(
               TString::Format("Volume %s: node %s overlapping %s", vol->GetName(), name1.Data(), name2.Data()),
               node1->GetVolume(), node2->GetVolume(), &mat1, &mat2, kTRUE, safe);
            flags[nd * id1 + id0] = nodeovlp;
            fGeoManager->AddOverlap(nodeovlp);
         }
         // Max 100 points per marker
         if (nodeovlp->GetPolyMarker()->GetN() < 100)
            nodeovlp->SetNextPoint(pt[0], pt[1], pt[2]);
         if (nodeovlp->GetOverlap() < safe)
            nodeovlp->SetOverlap(safe);
      }
   }
   nav->GetCache()->ReleaseInfo();
   if (flags)
      delete[] flags;
   if (!novlps)
      return;
   Double_t capacity = vol->GetShape()->Capacity();
   capacity *= Double_t(iovlp) / Double_t(npoints);
   Double_t err = 1. / TMath::Sqrt(Double_t(iovlp));
   Info("CheckOverlapsBySampling", "#Found %d overlaps adding-up to %g +/- %g [cm3] for daughters of %s", novlps,
        capacity, err * capacity, vol->GetName());
}
 
/// @brief Materialize a TGeoOverlapResult into a TGeoOverlap
/// @param r The result to materialize
void TGeoChecker::MaterializeOverlap(const TGeoOverlapResult &r)
{
   auto *ov = new TGeoOverlap(r.fName.Data(), r.fVol1, r.fVol2, &r.fMat1, &r.fMat2, r.fIsOverlap, r.fMaxOverlap);
 
   for (const auto &p : r.fPoints)
      ov->SetNextPoint(p[0], p[1], p[2]);
 
   fGeoManager->AddOverlap(ov);
}
 
////////////////////////////////////////////////////////////////////////////////
/// Create a policy stamp for overlap checking
 
TGeoChecker::MeshPolicyStamp TGeoChecker::CurrentMeshPolicyStamp() const
{
   MeshPolicyStamp s;
   s.fNmeshPoints = fNmeshPoints;
   s.fNsegments = gGeoManager ? gGeoManager->GetNsegments() : -1;
   return s;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Compute number of overlaps combinations to check per volume
 
Int_t TGeoChecker::NChecksPerVolume(TGeoVolume *vol)
{
   if (vol->GetFinder())
      return 0;
   UInt_t nd = vol->GetNdaughters();
   if (!nd)
      return 0;
   Bool_t is_assembly = vol->IsAssembly();
   TGeoIterator next1(vol);
   TGeoIterator next2(vol);
   Int_t nchecks = 0;
   TGeoNode *node;
   UInt_t id;
   if (!is_assembly) {
      while ((node = next1())) {
         if (node->IsOverlapping()) {
            next1.Skip();
            continue;
         }
         if (!node->GetVolume()->IsAssembly()) {
            nchecks++;
            next1.Skip();
         }
      }
   }
   // now check if the daughters overlap with each other
   if (nd < 2)
      return nchecks;
   TGeoVoxelFinder *vox = vol->GetVoxels();
   if (!vox)
      return nchecks;
 
   TGeoNode *node1, *node01, *node02;
   Int_t novlp;
   Int_t *ovlps;
   Int_t ko;
   UInt_t io;
   for (id = 0; id < nd; id++) {
      node01 = vol->GetNode(id);
      if (node01->IsOverlapping())
         continue;
      vox->FindOverlaps(id);
      ovlps = node01->GetOverlaps(novlp);
      if (!ovlps)
         continue;
      for (ko = 0; ko < novlp; ko++) { // loop all possible overlapping candidates
         io = ovlps[ko];               // (node1, shaped, matrix1, points, fBuff1)
         if (io <= id)
            continue;
         node02 = vol->GetNode(io);
         if (node02->IsOverlapping())
            continue;
 
         // We have to check node against node1, but they may be assemblies
         if (node01->GetVolume()->IsAssembly()) {
            next1.Reset(node01->GetVolume());
            while ((node = next1())) {
               if (!node->GetVolume()->IsAssembly()) {
                  if (node02->GetVolume()->IsAssembly()) {
                     next2.Reset(node02->GetVolume());
                     while ((node1 = next2())) {
                        if (!node1->GetVolume()->IsAssembly()) {
                           nchecks++;
                           next2.Skip();
                        }
                     }
                  } else {
                     nchecks++;
                  }
                  next1.Skip();
               }
            }
         } else {
            // node not assembly
            if (node02->GetVolume()->IsAssembly()) {
               next2.Reset(node02->GetVolume());
               while ((node1 = next2())) {
                  if (!node1->GetVolume()->IsAssembly()) {
                     nchecks++;
                     next2.Skip();
                  }
               }
            } else {
               // node1 also not assembly
               nchecks++;
            }
         }
      }
      node01->SetOverlaps(nullptr, 0);
   }
   return nchecks;
}
 
const TGeoChecker::MeshPointsEntry &TGeoChecker::GetMeshPoints(const TGeoShape *shape) const
{
   static const TGeoChecker::MeshPointsEntry kEmpty;
   if (!shape)
      return kEmpty;
 
   std::shared_lock lock(fMeshCacheMutex);
   auto it = fMeshPointsCache.find(shape);
   if (it != fMeshPointsCache.end())
      return it->second;
 
   // Strict mode: caller forgot to include this shape in BuildMeshPointsCache
   ::Error("TGeoChecker::GetMeshPoints",
           "Mesh cache miss for shape %s. BuildMeshPointsCache did not cover all candidates.", shape->ClassName());
   return kEmpty;
}
 
void TGeoChecker::BuildMeshPointsCache(const std::vector<TGeoOverlapCandidate> &candidates)
{
   std::unique_lock lock(fMeshCacheMutex);
 
   const auto nowStamp = CurrentMeshPolicyStamp();
   const bool policyChanged = PolicyChanged(fMeshPolicyStamp, nowStamp);
 
   if (policyChanged) {
      fMeshPointsCache.clear();
      fMeshPointsCache.reserve(2 * candidates.size());
      fMeshPolicyStamp = nowStamp;
   }
 
   auto ensureShape = [&](const TGeoShape *shape) {
      if (!shape)
         return;
      if (fMeshPointsCache.find(shape) != fMeshPointsCache.end())
         return;
      MeshPointsEntry entry;
 
      // IMPORTANT: meshN is the *only* correct size for SetPoints() buffer.
      Int_t meshN = 0, meshS = 0, meshP = 0;
      const_cast<TGeoShape *>(shape)->GetMeshNumbers(meshN, meshS, meshP);
 
      // 1) Base mesh points via SetPoints (exact sizing!)
      if (meshN > 0) {
         entry.fPoints.resize(3 * meshN);
         const_cast<TGeoShape *>(shape)->SetPoints(entry.fPoints.data());
         entry.fNPoints = meshN;
      } else {
         entry.fNPoints = 0;
      }
 
      // 2) Extra points on segments: fill into exactly-sized array, then append
      const Int_t nSegWanted = fNmeshPoints;
      if (nSegWanted > 0) {
         std::vector<Double_t> seg(3 * nSegWanted);
         Bool_t usedSeg = const_cast<TGeoShape *>(shape)->GetPointsOnSegments(nSegWanted, seg.data());
         if (usedSeg) {
            const auto oldSize = entry.fPoints.size();
            entry.fPoints.resize(oldSize + seg.size());
            std::memcpy(entry.fPoints.data() + oldSize, seg.data(), seg.size() * sizeof(Double_t));
            entry.fNPoints += nSegWanted;
         }
      }
 
      // 3) Hardening: ensure consistency between count and storage.
      // (No dedup; just guard against inconsistencies.)
      if ((Int_t)entry.fPoints.size() != 3 * entry.fNPoints) {
         // This should never happen with the logic above, but keep a guard.
         const Int_t n3 = entry.fPoints.size() / 3;
         entry.fNPoints = n3;
         entry.fPoints.resize(3 * n3);
      }
 
      fMeshPointsCache.emplace(shape, std::move(entry));
   };
 
   // Ensure all shapes referenced by *this* candidate set exist in the cache.
   for (const auto &c : candidates) {
      if (c.fVol1)
         ensureShape(c.fVol1->GetShape());
      if (c.fVol2)
         ensureShape(c.fVol2->GetShape());
   }
}
 
Int_t TGeoChecker::FillMeshPoints(TBuffer3D &buff, const TGeoShape *shape, const TGeoShape *&lastShape, Int_t &cachedN,
                                  Int_t nMeshPoints, Double_t *&points) const
{
   // Cache hit: reuse last filled content for this shape
   if (lastShape == shape && cachedN > 0) {
      points = buff.fPnts;
      return cachedN;
   }
 
   // Refill: compute mesh numbers only now
   Int_t meshN = buff.NbPnts();
   Int_t meshS = buff.NbSegs();
   Int_t meshP = buff.NbPols();
   shape->GetMeshNumbers(meshN, meshS, meshP);
 
   const Int_t nraw = TMath::Max(meshN, nMeshPoints);
   buff.SetRawSizes(nraw, 3 * nraw, 0, 0, 0, 0);
   points = buff.fPnts;
 
   const Bool_t useSeg = shape->GetPointsOnSegments(nMeshPoints, points);
 
   Int_t filledN = 0;
   if (useSeg) {
      filledN = nMeshPoints;
   } else {
      shape->SetPoints(points);
      filledN = meshN;
   }
 
   lastShape = shape;
   cachedN = filledN;
   return filledN;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Stage1: Enumerate all overlap candidates for volume VOL within a limit OVLP.
 
Int_t TGeoChecker::EnumerateOverlapCandidates(const TGeoVolume *vol, Double_t ovlp, Option_t *option,
                                              std::vector<TGeoOverlapCandidate> &out)
{
   if (vol->GetFinder())
      return 0;
   UInt_t nd = vol->GetNdaughters();
   if (!nd)
      return 0;
 
   const Bool_t is_assembly = vol->IsAssembly();
 
   TString opt(option);
   opt.ToLower();
   fFullCheck = opt.Contains("f");
 
   Int_t ncand = 0;
 
   // ---- EXTRUSIONS (only for daughters of a non-assembly volume)
   if (!is_assembly) {
      TGeoIterator next1((TGeoVolume *)vol);
      TGeoNode *node = nullptr;
      TGeoNode *nodecheck = nullptr;
      TString path;
      Int_t level;
 
      while ((node = next1())) {
         if (node->IsOverlapping()) {
            next1.Skip();
            continue;
         }
         if (!node->GetVolume()->IsAssembly()) {
            if (fSelectedNode) {
               // only overlaps of the selected node OR real daughters if selected is assembly
               if ((fSelectedNode != node) && (!fSelectedNode->GetVolume()->IsAssembly())) {
                  next1.Skip();
                  continue;
               }
               if (node != fSelectedNode) {
                  level = next1.GetLevel();
                  while ((nodecheck = next1.GetNode(level--))) {
                     if (nodecheck == fSelectedNode)
                        break;
                  }
                  if (!nodecheck) {
                     next1.Skip();
                     continue;
                  }
               }
            }
 
            next1.GetPath(path);
            ncand++;
            PushCandidate(out, TString::Format("%s extruded by: %s", vol->GetName(), path.Data()), (TGeoVolume *)vol,
                          node->GetVolume(), gGeoIdentity, next1.GetCurrentMatrix(), kFALSE, ovlp);
 
            next1.Skip();
         }
      }
   }
 
   // ---- OVERLAPS between daughters
   if (nd < 2)
      return ncand;
 
   TGeoVoxelFinder *vox = vol->GetVoxels();
   if (!vox) {
      Warning("EnumerateOverlapCandidates", "Volume %s with %i daughters but not voxelized", vol->GetName(), nd);
      return ncand;
   }
 
   TGeoIterator next1((TGeoVolume *)vol);
   TGeoIterator next2((TGeoVolume *)vol);
 
   TGeoNode *node1 = nullptr;
   TGeoNode *node01 = nullptr;
   TGeoNode *node02 = nullptr;
   TGeoNode *nodecheck = nullptr;
 
   TGeoHMatrix hmat1, hmat2;
   TGeoMatrix const *tr1, *tr2;
   TGeoBBox const *box1, *box2;
   TString path, path1;
   UInt_t id, io;
   Int_t level;
 
   for (id = 0; id < nd; id++) {
      node01 = vol->GetNode(id);
      if (node01->IsOverlapping())
         continue;
 
      next1.SetTopName(node01->GetName());
      path = node01->GetName();
      //=== ANY <-> ANY ===//
      for (io = id + 1; io < nd; io++) {
         // check node01 against node02 daughters of volume
         node02 = vol->GetNode(io);
         if (node02->IsOverlapping())
            continue;
         box1 = (const TGeoBBox *)node01->GetVolume()->GetShape();
         tr1 = node01->GetMatrix();
         box2 = (const TGeoBBox *)node02->GetVolume()->GetShape();
         tr2 = node02->GetMatrix();
         // OBB check
         if (!TGeoBBox::MayIntersect(box1, tr1, box2, tr2))
            continue;
 
         next2.SetTopName(node02->GetName());
         path1 = node02->GetName();
 
         // We have to check node01 against node02, but they may be assemblies
         if (node01->GetVolume()->IsAssembly()) {
            // left node assembly - make a visitor
            next1.Reset(node01->GetVolume());
            TGeoNode *node = nullptr; // will iterate components of node01
            //=== ASSEMBLY/ANY <-> ANY ===//
            while ((node = next1())) {
               hmat1 = node01->GetMatrix();
               hmat1 *= *next1.GetCurrentMatrix();
               box1 = (const TGeoBBox *)node->GetVolume()->GetShape();
               tr1 = &hmat1;
               // OBB check
               if (!TGeoBBox::MayIntersect(box1, tr1, box2, tr2)) {
                  // No intersection, skip the full left branch
                  next1.Skip();
                  continue;
               }
 
               if (!node->GetVolume()->IsAssembly()) {
                  //=== ASSEMBLY/REAL <-> ANY ===//
                  // Current daughter of node01 is real, get its path and transformation
                  next1.GetPath(path);
 
                  if (node02->GetVolume()->IsAssembly()) {
                     next2.Reset(node02->GetVolume());
                     //=== ASSEMBLY/REAL <-> ASSEMBLY/ANY ===//
                     while ((node1 = next2())) {
                        hmat2 = node02->GetMatrix();
                        hmat2 *= *next2.GetCurrentMatrix();
                        box2 = (const TGeoBBox *)node1->GetVolume()->GetShape();
                        // OBB check
                        if (!TGeoBBox::MayIntersect(box1, &hmat1, box2, &hmat2)) {
                           // No intersection, skip the full right branch
                           next2.Skip();
                           continue;
                        }
                        if (!node1->GetVolume()->IsAssembly()) {
                           //=== ASSEMBLY/REAL <-> ASSEMBLY/REAL ===//
                           // Selected node skip logic
                           if (fSelectedNode) {
                              if ((fSelectedNode != node) && (fSelectedNode != node1) &&
                                  (!fSelectedNode->GetVolume()->IsAssembly())) {
                                 next2.Skip();
                                 continue;
                              }
                              if ((fSelectedNode != node1) && (fSelectedNode != node)) {
                                 level = next2.GetLevel();
                                 while ((nodecheck = next2.GetNode(level--))) {
                                    if (nodecheck == fSelectedNode)
                                       break;
                                 }
                                 if (node02 == fSelectedNode)
                                    nodecheck = node02;
                                 if (!nodecheck) {
                                    level = next1.GetLevel();
                                    while ((nodecheck = next1.GetNode(level--))) {
                                       if (nodecheck == fSelectedNode)
                                          break;
                                    }
                                 }
                                 if (node01 == fSelectedNode)
                                    nodecheck = node01;
                                 if (!nodecheck) {
                                    next2.Skip();
                                    continue;
                                 }
                              }
                           }
 
                           next2.GetPath(path1);
 
                           ncand++;
                           PushCandidate(out,
                                         TString::Format("%s/%s overlapping %s/%s", vol->GetName(), path.Data(),
                                                         vol->GetName(), path1.Data()),
                                         node->GetVolume(), node1->GetVolume(), &hmat1, &hmat2, kTRUE, ovlp);
 
                           next2.Skip();
                        }
                     }
 
                  } else {
                     //=== ASSEMBLY/REAL <-> REAL
                     // Selected node skip logic
                     if (fSelectedNode) {
                        if ((fSelectedNode != node) && (fSelectedNode != node02) &&
                            (!fSelectedNode->GetVolume()->IsAssembly())) {
                           next1.Skip();
                           continue;
                        }
                        if ((fSelectedNode != node) && (fSelectedNode != node02)) {
                           level = next1.GetLevel();
                           while ((nodecheck = next1.GetNode(level--))) {
                              if (nodecheck == fSelectedNode)
                                 break;
                           }
                           if (node01 == fSelectedNode)
                              nodecheck = node01;
                           if (!nodecheck) {
                              next1.Skip();
                              continue;
                           }
                        }
                     }
 
                     ncand++;
                     PushCandidate(out,
                                   TString::Format("%s/%s overlapping %s/%s", vol->GetName(), path.Data(),
                                                   vol->GetName(), path1.Data()),
                                   node->GetVolume(), node02->GetVolume(), &hmat1, node02->GetMatrix(), kTRUE, ovlp);
                  }
 
                  next1.Skip();
               }
            }
 
         } else {
            if (node02->GetVolume()->IsAssembly()) {
               next2.Reset(node02->GetVolume());
               //=== REAL <-> ASSEMBLY/ANY ===//
               while ((node1 = next2())) {
                  hmat2 = node02->GetMatrix();
                  hmat2 *= *next2.GetCurrentMatrix();
                  box2 = (const TGeoBBox *)node1->GetVolume()->GetShape();
                  // OBB check
                  if (!TGeoBBox::MayIntersect(box1, node01->GetMatrix(), box2, &hmat2)) {
                     // No intersection, skip the entire right branch
                     next2.Skip();
                     continue;
                  }
 
                  if (!node1->GetVolume()->IsAssembly()) {
                     // Selected node skip logic
                     if (fSelectedNode) {
                        if ((fSelectedNode != node1) && (fSelectedNode != node01) &&
                            (!fSelectedNode->GetVolume()->IsAssembly())) {
                           next2.Skip();
                           continue;
                        }
                        if ((fSelectedNode != node1) && (fSelectedNode != node01)) {
                           level = next2.GetLevel();
                           while ((nodecheck = next2.GetNode(level--))) {
                              if (nodecheck == fSelectedNode)
                                 break;
                           }
                           if (node02 == fSelectedNode)
                              nodecheck = node02;
                           if (!nodecheck) {
                              next2.Skip();
                              continue;
                           }
                        }
                     }
 
                     ncand++;
                     next2.GetPath(path1);
                     PushCandidate(out,
                                   TString::Format("%s/%s overlapping %s/%s", vol->GetName(), path.Data(),
                                                   vol->GetName(), path1.Data()),
                                   node01->GetVolume(), node1->GetVolume(), node01->GetMatrix(), &hmat2, kTRUE, ovlp);
 
                     next2.Skip();
                  }
               }
 
            } else {
               // both not assembly
               if (fSelectedNode && (fSelectedNode != node01) && (fSelectedNode != node02))
                  continue;
 
               ncand++;
               PushCandidate(
                  out,
                  TString::Format("%s/%s overlapping %s/%s", vol->GetName(), path.Data(), vol->GetName(), path1.Data()),
                  node01->GetVolume(), node02->GetVolume(), node01->GetMatrix(), node02->GetMatrix(), kTRUE, ovlp);
            }
         }
      }
   }
   return ncand;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Helper to fill candidates list
 
void TGeoChecker::PushCandidate(std::vector<TGeoOverlapCandidate> &out, const TString &name, TGeoVolume *vol1,
                                TGeoVolume *vol2, const TGeoMatrix *mat1, const TGeoMatrix *mat2, Bool_t isovlp,
                                Double_t ovlp) const
{
   TGeoOverlapCandidate c;
   c.fName = name;
   c.fVol1 = vol1;
   c.fVol2 = vol2;
   c.fMat1 = TGeoHMatrix(*mat1);
   c.fMat2 = TGeoHMatrix(*mat2);
   c.fIsOverlap = isovlp;
   c.fOvlp = ovlp;
   out.emplace_back(std::move(c));
}
 
////////////////////////////////////////////////////////////////////////////////
/// Print the current list of overlaps held by the manager class.
 
void TGeoChecker::PrintOverlaps() const
{
   TIter next(fGeoManager->GetListOfOverlaps());
   TGeoOverlap *ov;
   printf("=== Overlaps for %s ===\n", fGeoManager->GetName());
   while ((ov = (TGeoOverlap *)next()))
      ov->PrintInfo();
}
 
////////////////////////////////////////////////////////////////////////////////
/// Draw point (x,y,z) over the picture of the daughters 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.
 
void TGeoChecker::CheckPoint(Double_t x, Double_t y, Double_t z, Option_t *, Double_t safety)
{
   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 = (safety > 0.) ? safety : 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();
}
 
////////////////////////////////////////////////////////////////////////////////
/// Test for shape navigation methods. Summary for test numbers:
///  - 1: DistFromInside/Outside. Sample points inside the shape. Generate
///    directions randomly in cos(theta). Compute DistFromInside and move the
///    point with bigger distance. Compute DistFromOutside back from new point.
///    Plot d-(d1+d2)
///  - 2: Safety test. Sample points inside the bounding and compute safety. Generate
///    directions randomly in cos(theta) and compute distance to boundary. Check if
///    distance to boundary is bigger than safety
 
void TGeoChecker::CheckShape(TGeoShape *shape, Int_t testNo, Int_t nsamples, Option_t *option)
{
   switch (testNo) {
   case 1: ShapeDistances(shape, nsamples, option); break;
   case 2: ShapeSafety(shape, nsamples, option); break;
   case 3: ShapeNormal(shape, nsamples, option); break;
   default: Error("CheckShape", "Test number %d not existent", testNo);
   }
}
 
////////////////////////////////////////////////////////////////////////////////
///  Test TGeoShape::DistFromInside/Outside. Sample points inside the shape. Generate
/// directions randomly in cos(theta). Compute d1 = DistFromInside and move the
/// point on the boundary. Compute DistFromOutside and propagate with d2 making sure that
/// the shape is not re-entered. Swap direction and call DistFromOutside that
/// should fall back on the same point on the boundary (at d2). Propagate back on boundary
/// then compute DistFromInside that should be bigger than d1.
///    Plot d-(d1+d2)
 
void TGeoChecker::ShapeDistances(TGeoShape *shape, Int_t nsamples, Option_t *)
{
   Double_t dx = ((TGeoBBox *)shape)->GetDX();
   Double_t dy = ((TGeoBBox *)shape)->GetDY();
   Double_t dz = ((TGeoBBox *)shape)->GetDZ();
   // the box might be displaced
   auto origin_p = ((TGeoBBox *)shape)->GetOrigin();
   Double_t ox = origin_p[0];
   Double_t oy = origin_p[1];
   Double_t oz = origin_p[2];
 
   Double_t dmax = 2. * TMath::Sqrt(dx * dx + dy * dy + dz * dz);
   Double_t d1, d2, dmove, dnext;
   Int_t itot = 0;
   // Number of tracks shot for every point inside the shape
   const Int_t kNtracks = 1000;
   Int_t n10 = nsamples / 10;
   Int_t i, j;
   Double_t point[3], pnew[3];
   Double_t dir[3], dnew[3];
   Double_t theta, phi, delta;
   TPolyMarker3D *pmfrominside = nullptr;
   TPolyMarker3D *pmfromoutside = nullptr;
   TH1D *hist = new TH1D("hTest1", "Residual distance from inside/outside", 200, -20, 0);
   hist->GetXaxis()->SetTitle("delta[cm] - first bin=overflow");
   hist->GetYaxis()->SetTitle("count");
   hist->SetMarkerStyle(kFullCircle);
 
   if (!fTimer)
      fTimer = new TStopwatch();
   fTimer->Reset();
   fTimer->Start();
   while (itot < nsamples) {
      Bool_t inside = kFALSE;
      while (!inside) {
         point[0] = ox + gRandom->Uniform(-dx, dx);
         point[1] = oy + gRandom->Uniform(-dy, dy);
         point[2] = oz + gRandom->Uniform(-dz, dz);
         inside = shape->Contains(point);
      }
      itot++;
      if (n10) {
         if ((itot % n10) == 0)
            printf("%i percent\n", Int_t(100 * itot / nsamples));
      }
      for (i = 0; i < kNtracks; i++) {
         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);
         dmove = dmax;
         // We have track direction, compute distance from inside
         d1 = shape->DistFromInside(point, dir, 3);
         if (d1 > dmove || d1 < TGeoShape::Tolerance()) {
            // Bad distance or bbox size, to debug
            new TCanvas("shape01", Form("Shape %s (%s)", shape->GetName(), shape->ClassName()), 1000, 800);
            shape->Draw();
            printf("DistFromInside: (%19.15f, %19.15f, %19.15f, %19.15f, %19.15f, %19.15f) %f/%f(max)\n", point[0],
                   point[1], point[2], dir[0], dir[1], dir[2], d1, dmove);
            pmfrominside = new TPolyMarker3D(2);
            pmfrominside->SetMarkerColor(kRed);
            pmfrominside->SetMarkerStyle(24);
            pmfrominside->SetMarkerSize(0.4);
            pmfrominside->SetNextPoint(point[0], point[1], point[2]);
            for (j = 0; j < 3; j++)
               pnew[j] = point[j] + d1 * dir[j];
            pmfrominside->SetNextPoint(pnew[0], pnew[1], pnew[2]);
            pmfrominside->Draw();
            return;
         }
         // Propagate BEFORE the boundary and make sure that DistFromOutside
         // does not return 0 (!!!)
         // Check if there is a second crossing
         for (j = 0; j < 3; j++)
            pnew[j] = point[j] + (d1 - TGeoShape::Tolerance()) * dir[j];
         dnext = shape->DistFromOutside(pnew, dir, 3);
         if (d1 + dnext < dmax)
            dmove = d1 + 0.5 * dnext;
         // Move point and swap direction
         for (j = 0; j < 3; j++) {
            pnew[j] = point[j] + dmove * dir[j];
            dnew[j] = -dir[j];
         }
         // Compute now distance from outside
         d2 = shape->DistFromOutside(pnew, dnew, 3);
         delta = dmove - d1 - d2;
         if (TMath::Abs(delta) > 1E-6 || dnext < 2. * TGeoShape::Tolerance()) {
            // Error->debug this
            new TCanvas("shape01", Form("Shape %s (%s)", shape->GetName(), shape->ClassName()), 1000, 800);
            shape->Draw();
            printf("Error: (%19.15f, %19.15f, %19.15f, %19.15f, %19.15f, %19.15f) d1=%f d2=%f dmove=%f\n", point[0],
                   point[1], point[2], dir[0], dir[1], dir[2], d1, d2, dmove);
            if (dnext < 2. * TGeoShape::Tolerance()) {
               printf(" (*)DistFromOutside(%19.15f, %19.15f, %19.15f, %19.15f, %19.15f, %19.15f)  dnext = %f\n",
                      point[0] + (d1 - TGeoShape::Tolerance()) * dir[0],
                      point[1] + (d1 - TGeoShape::Tolerance()) * dir[1],
                      point[2] + (d1 - TGeoShape::Tolerance()) * dir[2], dir[0], dir[1], dir[2], dnext);
            } else {
               printf("   DistFromOutside(%19.15f, %19.15f, %19.15f, %19.15f, %19.15f, %19.15f)  dnext = %f\n",
                      point[0] + d1 * dir[0], point[1] + d1 * dir[1], point[2] + d1 * dir[2], dir[0], dir[1], dir[2],
                      dnext);
            }
            printf("   DistFromOutside(%19.15f, %19.15f, %19.15f, %19.15f, %19.15f, %19.15f)  = %f\n", pnew[0], pnew[1],
                   pnew[2], dnew[0], dnew[1], dnew[2], d2);
            pmfrominside = new TPolyMarker3D(2);
            pmfrominside->SetMarkerStyle(24);
            pmfrominside->SetMarkerSize(0.4);
            pmfrominside->SetMarkerColor(kRed);
            pmfrominside->SetNextPoint(point[0], point[1], point[2]);
            for (j = 0; j < 3; j++)
               point[j] += d1 * dir[j];
            pmfrominside->SetNextPoint(point[0], point[1], point[2]);
            pmfrominside->Draw();
            pmfromoutside = new TPolyMarker3D(2);
            pmfromoutside->SetMarkerStyle(20);
            pmfromoutside->SetMarkerStyle(7);
            pmfromoutside->SetMarkerSize(0.3);
            pmfromoutside->SetMarkerColor(kBlue);
            pmfromoutside->SetNextPoint(pnew[0], pnew[1], pnew[2]);
            for (j = 0; j < 3; j++)
               pnew[j] += d2 * dnew[j];
            if (d2 < 1E10)
               pmfromoutside->SetNextPoint(pnew[0], pnew[1], pnew[2]);
            pmfromoutside->Draw();
            return;
         }
         // Compute distance from inside which should be bigger than d1
         for (j = 0; j < 3; j++)
            pnew[j] += (d2 - TGeoShape::Tolerance()) * dnew[j];
         dnext = shape->DistFromInside(pnew, dnew, 3);
         if (dnext < d1 - TGeoShape::Tolerance() || dnext > dmax) {
            new TCanvas("shape01", Form("Shape %s (%s)", shape->GetName(), shape->ClassName()), 1000, 800);
            shape->Draw();
            printf("Error DistFromInside(%19.15f, %19.15f, %19.15f, %19.15f, %19.15f, %19.15f) d1=%f d1p=%f\n", pnew[0],
                   pnew[1], pnew[2], dnew[0], dnew[1], dnew[2], d1, dnext);
            pmfrominside = new TPolyMarker3D(2);
            pmfrominside->SetMarkerStyle(24);
            pmfrominside->SetMarkerSize(0.4);
            pmfrominside->SetMarkerColor(kRed);
            pmfrominside->SetNextPoint(point[0], point[1], point[2]);
            for (j = 0; j < 3; j++)
               point[j] += d1 * dir[j];
            pmfrominside->SetNextPoint(point[0], point[1], point[2]);
            pmfrominside->Draw();
            pmfromoutside = new TPolyMarker3D(2);
            pmfromoutside->SetMarkerStyle(20);
            pmfromoutside->SetMarkerStyle(7);
            pmfromoutside->SetMarkerSize(0.3);
            pmfromoutside->SetMarkerColor(kBlue);
            pmfromoutside->SetNextPoint(pnew[0], pnew[1], pnew[2]);
            for (j = 0; j < 3; j++)
               pnew[j] += dnext * dnew[j];
            if (d2 < 1E10)
               pmfromoutside->SetNextPoint(pnew[0], pnew[1], pnew[2]);
            pmfromoutside->Draw();
            return;
         }
         if (TMath::Abs(delta) < 1E-20)
            delta = 1E-30;
         hist->Fill(TMath::Max(TMath::Log(TMath::Abs(delta)), -20.));
      }
   }
   fTimer->Stop();
   fTimer->Print();
   new TCanvas("Test01", "Residuals DistFromInside/Outside", 800, 600);
   hist->Draw();
}
 
////////////////////////////////////////////////////////////////////////////////
/// Check of validity of safe distance for a given shape.
/// Sample points inside the 2x bounding box and compute safety. Generate
/// directions randomly in cos(theta) and compute distance to boundary. Check if
/// distance to boundary is bigger than safety.
 
void TGeoChecker::ShapeSafety(TGeoShape *shape, Int_t nsamples, Option_t *)
{
   Double_t dx = ((TGeoBBox *)shape)->GetDX();
   Double_t dy = ((TGeoBBox *)shape)->GetDY();
   Double_t dz = ((TGeoBBox *)shape)->GetDZ();
   // the box might be displaced
   auto origin_p = ((TGeoBBox *)shape)->GetOrigin();
   Double_t ox = origin_p[0];
   Double_t oy = origin_p[1];
   Double_t oz = origin_p[2];
 
   // Number of tracks shot for every point inside the shape
   const Int_t kNtracks = 1000;
   Int_t n10 = nsamples / 10;
   Int_t i;
   Double_t dist;
   Double_t point[3];
   Double_t dir[3];
   Double_t theta, phi;
   TPolyMarker3D *pm1 = nullptr;
   TPolyMarker3D *pm2 = nullptr;
   if (!fTimer)
      fTimer = new TStopwatch();
   fTimer->Reset();
   fTimer->Start();
   Int_t itot = 0;
   while (itot < nsamples) {
      Bool_t inside = kFALSE;
      point[0] = ox + gRandom->Uniform(-2 * dx, 2 * dx);
      point[1] = oy + gRandom->Uniform(-2 * dy, 2 * dy);
      point[2] = oz + gRandom->Uniform(-2 * dz, 2 * dz);
      inside = shape->Contains(point);
      Double_t safe = shape->Safety(point, inside);
      itot++;
      if (n10) {
         if ((itot % n10) == 0)
            printf("%i percent\n", Int_t(100 * itot / nsamples));
      }
      for (i = 0; i < kNtracks; i++) {
         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);
         if (inside)
            dist = shape->DistFromInside(point, dir, 3);
         else
            dist = shape->DistFromOutside(point, dir, 3);
         if (dist < safe) {
            printf("Error safety (%19.15f, %19.15f, %19.15f, %19.15f, %19.15f, %19.15f) safe=%f  dist=%f\n", point[0],
                   point[1], point[2], dir[0], dir[1], dir[2], safe, dist);
            new TCanvas("shape02", Form("Shape %s (%s)", shape->GetName(), shape->ClassName()), 1000, 800);
            shape->Draw();
            pm1 = new TPolyMarker3D(2);
            pm1->SetMarkerStyle(24);
            pm1->SetMarkerSize(0.4);
            pm1->SetMarkerColor(kRed);
            pm1->SetNextPoint(point[0], point[1], point[2]);
            pm1->SetNextPoint(point[0] + safe * dir[0], point[1] + safe * dir[1], point[2] + safe * dir[2]);
            pm1->Draw();
            pm2 = new TPolyMarker3D(1);
            pm2->SetMarkerStyle(7);
            pm2->SetMarkerSize(0.3);
            pm2->SetMarkerColor(kBlue);
            pm2->SetNextPoint(point[0] + dist * dir[0], point[1] + dist * dir[1], point[2] + dist * dir[2]);
            pm2->Draw();
            return;
         }
      }
   }
   fTimer->Stop();
   fTimer->Print();
}
 
////////////////////////////////////////////////////////////////////////////////
/// Check of validity of the normal for a given shape.
/// Sample points inside the shape. Generate directions randomly in cos(theta)
/// and propagate to boundary. Compute normal and safety at crossing point, plot
/// the point and generate a random direction so that (dir) dot (norm) <0.
 
void TGeoChecker::ShapeNormal(TGeoShape *shape, Int_t nsamples, Option_t *)
{
   Double_t dx = ((TGeoBBox *)shape)->GetDX();
   Double_t dy = ((TGeoBBox *)shape)->GetDY();
   Double_t dz = ((TGeoBBox *)shape)->GetDZ();
   // the box might be displaced
   auto origin_p = ((TGeoBBox *)shape)->GetOrigin();
   Double_t ox = origin_p[0];
   Double_t oy = origin_p[1];
   Double_t oz = origin_p[2];
 
   Double_t dmax = 2. * TMath::Sqrt(dx * dx + dy * dy + dz * dz);
   // Number of tracks shot for every point inside the shape
   const Int_t kNtracks = 1000;
   Int_t n10 = nsamples / 10;
   Int_t itot = 0, errcnt = 0, errsame = 0;
   Int_t i;
   Double_t dist, olddist, safe, dot;
   Double_t theta, phi, ndotd;
   Double_t *spoint = new Double_t[3 * nsamples];
   Double_t *sdir = new Double_t[3 * nsamples];
   while (itot < nsamples) {
      Bool_t inside = kFALSE;
      while (!inside) {
         spoint[3 * itot] = ox + gRandom->Uniform(-dx, dx);
         spoint[3 * itot + 1] = oy + gRandom->Uniform(-dy, dy);
         spoint[3 * itot + 2] = oz + gRandom->Uniform(-dz, dz);
         inside = shape->Contains(&spoint[3 * itot]);
      }
      phi = 2 * TMath::Pi() * gRandom->Rndm();
      theta = TMath::ACos(1. - 2. * gRandom->Rndm());
      sdir[3 * itot] = TMath::Sin(theta) * TMath::Cos(phi);
      sdir[3 * itot + 1] = TMath::Sin(theta) * TMath::Sin(phi);
      sdir[3 * itot + 2] = TMath::Cos(theta);
      itot++;
   }
   Double_t point[3], newpoint[3], oldpoint[3];
   Double_t dir[3], olddir[3];
   Double_t norm[3], newnorm[3], oldnorm[3];
   TCanvas *errcanvas = nullptr;
   TPolyMarker3D *pm1 = nullptr;
   TPolyMarker3D *pm2 = nullptr;
   pm2 = new TPolyMarker3D();
   //   pm2->SetMarkerStyle(24);
   pm2->SetMarkerSize(0.2);
   pm2->SetMarkerColor(kBlue);
   if (!fTimer)
      fTimer = new TStopwatch();
   fTimer->Reset();
   fTimer->Start();
   for (itot = 0; itot < nsamples; itot++) {
      if (n10) {
         if ((itot % n10) == 0)
            printf("%i percent\n", Int_t(100 * itot / nsamples));
      }
      oldnorm[0] = oldnorm[1] = oldnorm[2] = 0.;
      olddist = 0.;
      for (Int_t j = 0; j < 3; j++) {
         oldpoint[j] = point[j] = spoint[3 * itot + j];
         olddir[j] = dir[j] = sdir[3 * itot + j];
      }
      for (i = 0; i < kNtracks; i++) {
         if (errcnt > 0)
            break;
         dist = shape->DistFromInside(point, dir, 3);
         if (dist > dmax) {
            // distance is bigger than from bounding box
            // should not happen for point inside
            printf("Error DistFromInside Too large(%19.15f, %19.15f, %19.15f, %19.15f, %19.15f, %19.15f) =%g\n",
                   point[0], point[1], point[2], dir[0], dir[1], dir[2], dist);
         }
         for (Int_t j = 0; j < 3; j++) {
            newpoint[j] = point[j] + dist * dir[j];
         }
         shape->ComputeNormal(newpoint, dir, newnorm);
 
         dot = olddir[0] * oldnorm[0] + olddir[1] * oldnorm[1] + olddir[2] * oldnorm[2];
         if (!shape->Contains(point) && shape->Safety(point, kFALSE) > 1.E-3) {
            errcnt++;
            printf("Error point outside (%19.15f, %19.15f, %19.15f, %19.15f, %19.15f, %19.15f) =%g olddist=%g\n",
                   point[0], point[1], point[2], dir[0], dir[1], dir[2], dist, olddist);
            printf("         old point: (%19.15f, %19.15f, %19.15f, %19.15f, %19.15f, %19.15f)\n", oldpoint[0],
                   oldpoint[1], oldpoint[2], olddir[0], olddir[1], olddir[2]);
            if (!errcanvas)
               errcanvas =
                  new TCanvas("shape_err03", Form("Shape %s (%s)", shape->GetName(), shape->ClassName()), 1000, 800);
            if (!pm1) {
               pm1 = new TPolyMarker3D();
               pm1->SetMarkerStyle(24);
               pm1->SetMarkerSize(0.4);
               pm1->SetMarkerColor(kRed);
            }
            pm1->SetNextPoint(point[0], point[1], point[2]);
            pm1->SetNextPoint(oldpoint[0], oldpoint[1], oldpoint[2]);
            break;
         }
         if ((dist < TGeoShape::Tolerance() && olddist * dot > 1.E-3) || dist > dmax) {
            errsame++;
            if (errsame > 1) {
               errcnt++;
               printf("Error DistFromInside(%19.15f, %19.15f, %19.15f, %19.15f, %19.15f, %19.15f) =%g olddist=%g\n",
                      point[0], point[1], point[2], dir[0], dir[1], dir[2], dist, olddist);
               printf("         new norm:  (%g, %g, %g)\n", newnorm[0], newnorm[1], newnorm[2]);
               printf("         old point: (%19.15f, %19.15f, %19.15f, %19.15f, %19.15f, %19.15f)\n", oldpoint[0],
                      oldpoint[1], oldpoint[2], olddir[0], olddir[1], olddir[2]);
               printf("         old norm:  (%g, %g, %g)\n", oldnorm[0], oldnorm[1], oldnorm[2]);
               if (!errcanvas)
                  errcanvas =
                     new TCanvas("shape_err03", Form("Shape %s (%s)", shape->GetName(), shape->ClassName()), 1000, 800);
               if (!pm1) {
                  pm1 = new TPolyMarker3D();
                  pm1->SetMarkerStyle(24);
                  pm1->SetMarkerSize(0.4);
                  pm1->SetMarkerColor(kRed);
               }
               pm1->SetNextPoint(point[0], point[1], point[2]);
               pm1->SetNextPoint(oldpoint[0], oldpoint[1], oldpoint[2]);
               break;
            }
         } else
            errsame = 0;
 
         olddist = dist;
         for (Int_t j = 0; j < 3; j++) {
            oldpoint[j] = point[j];
            point[j] += dist * dir[j];
         }
         safe = shape->Safety(point, kTRUE);
         if (safe > 1.E-3) {
            errcnt++;
            printf("Error safety (%19.15f, %19.15f, %19.15f) safe=%g\n", point[0], point[1], point[2], safe);
            if (!errcanvas)
               errcanvas =
                  new TCanvas("shape_err03", Form("Shape %s (%s)", shape->GetName(), shape->ClassName()), 1000, 800);
            if (!pm1) {
               pm1 = new TPolyMarker3D();
               pm1->SetMarkerStyle(24);
               pm1->SetMarkerSize(0.4);
               pm1->SetMarkerColor(kRed);
            }
            pm1->SetNextPoint(point[0], point[1], point[2]);
            break;
         }
         // Compute normal
         shape->ComputeNormal(point, dir, norm);
         memcpy(oldnorm, norm, 3 * sizeof(Double_t));
         memcpy(olddir, dir, 3 * sizeof(Double_t));
         while (true) {
            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);
            ndotd = dir[0] * norm[0] + dir[1] * norm[1] + dir[2] * norm[2];
            if (ndotd <= 0)
               break; // backwards, still inside shape
         }
         if ((itot % 10) == 0)
            pm2->SetNextPoint(point[0], point[1], point[2]);
      }
   }
   fTimer->Stop();
   fTimer->Print();
   if (errcanvas) {
      shape->Draw();
      pm1->Draw();
   }
 
   new TCanvas("shape03", Form("Shape %s (%s)", shape->GetName(), shape->ClassName()), 1000, 800);
   pm2->Draw();
   delete[] spoint;
   delete[] sdir;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Generate a lego plot fot the top volume, according to option.
 
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)
{
   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 = nullptr;
         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 == nullptr && 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;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Draw random points in the bounding box of a volume.
 
void TGeoChecker::RandomPoints(TGeoVolume *vol, Int_t npoints, Option_t *option)
{
   if (!vol)
      return;
   vol->VisibleDaughters(kTRUE);
   vol->Draw();
   TString opt = option;
   opt.ToLower();
   TObjArray *pm = new TObjArray(128);
   TPolyMarker3D *marker = nullptr;
   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 = nullptr;
   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 = 1;
      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);
   if (igen)
      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;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Randomly shoot nrays from point (startx,starty,startz) and plot intersections
/// with surfaces for current top node.
 
void TGeoChecker::RandomRays(Int_t nrays, Double_t startx, Double_t starty, Double_t startz, const char *target_vol,
                             Bool_t check_norm)
{
   TObjArray *pm = new TObjArray(128);
   TString starget = target_vol;
   TPolyLine3D *line = nullptr;
   TPolyLine3D *normline = nullptr;
   TGeoVolume *vol = fGeoManager->GetTopVolume();
   //   vol->VisibleDaughters(kTRUE);
 
   Bool_t random = kFALSE;
   if (nrays <= 0) {
      nrays = 100000;
      random = kTRUE;
   }
   Double_t start[3];
   Double_t dir[3];
   const 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, inull;
   Int_t itot = 0;
   Int_t n10 = nrays / 10;
   Double_t theta, phi, step, normlen;
   Double_t ox = ((TGeoBBox *)vol->GetShape())->GetOrigin()[0];
   Double_t oy = ((TGeoBBox *)vol->GetShape())->GetOrigin()[1];
   Double_t oz = ((TGeoBBox *)vol->GetShape())->GetOrigin()[2];
   Double_t dx = ((TGeoBBox *)vol->GetShape())->GetDX();
   Double_t dy = ((TGeoBBox *)vol->GetShape())->GetDY();
   Double_t dz = ((TGeoBBox *)vol->GetShape())->GetDZ();
   normlen = TMath::Max(dx, dy);
   normlen = TMath::Max(normlen, dz);
   normlen *= 0.05;
   while (itot < nrays) {
      itot++;
      inull = 0;
      ipoint = 0;
      if (n10) {
         if ((itot % n10) == 0)
            printf("%i percent\n", Int_t(100 * itot / nrays));
      }
      if (random) {
         start[0] = ox - dx + 2 * dx * gRandom->Rndm();
         start[1] = oy - dy + 2 * dy * gRandom->Rndm();
         start[2] = oz - dz + 2 * dz * gRandom->Rndm();
      } else {
         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 = nullptr;
      if (fGeoManager->IsOutside())
         startnode = nullptr;
      vis1 = kFALSE;
      if (target_vol) {
         if (startnode && starget == startnode->GetVolume()->GetName())
            vis1 = kTRUE;
      } else {
         if (startnode && startnode->IsOnScreen())
            vis1 = kTRUE;
      }
      if (vis1) {
         line = new TPolyLine3D(2);
         line->SetLineColor(startnode->GetVolume()->GetLineColor());
         line->SetPoint(ipoint++, start[0], start[1], start[2]);
         i++;
         pm->Add(line);
      }
      while ((endnode = fGeoManager->FindNextBoundaryAndStep())) {
         step = fGeoManager->GetStep();
         if (step < TGeoShape::Tolerance())
            inull++;
         else
            inull = 0;
         if (inull > 5)
            break;
         const Double_t *normal = nullptr;
         if (check_norm) {
            normal = fGeoManager->FindNormalFast();
            if (!normal)
               break;
         }
         vis2 = kFALSE;
         if (target_vol) {
            if (starget == endnode->GetVolume()->GetName())
               vis2 = kTRUE;
         } else if (endnode->IsOnScreen())
            vis2 = kTRUE;
         if (ipoint > 0) {
            // old visible node had an entry point -> finish segment
            line->SetPoint(ipoint, point[0], point[1], point[2]);
            if (!vis2 && check_norm) {
               normline = new TPolyLine3D(2);
               normline->SetLineColor(kBlue);
               normline->SetLineWidth(1);
               normline->SetPoint(0, point[0], point[1], point[2]);
               normline->SetPoint(1, point[0] + normal[0] * normlen, point[1] + normal[1] * normlen,
                                  point[2] + normal[2] * normlen);
               pm->Add(normline);
            }
            ipoint = 0;
            line = nullptr;
         }
         if (vis2) {
            // create new segment
            line = new TPolyLine3D(2);
            line->SetLineColor(endnode->GetVolume()->GetLineColor());
            line->SetPoint(ipoint++, point[0], point[1], point[2]);
            i++;
            if (check_norm) {
               normline = new TPolyLine3D(2);
               normline->SetLineColor(kBlue);
               normline->SetLineWidth(2);
               normline->SetPoint(0, point[0], point[1], point[2]);
               normline->SetPoint(1, point[0] + normal[0] * normlen, point[1] + normal[1] * normlen,
                                  point[2] + normal[2] * normlen);
            }
            pm->Add(line);
            if (!random)
               pm->Add(normline);
         }
      }
   }
   // 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;
}
 
////////////////////////////////////////////////////////////////////////////////
/// shoot npoints randomly in a box of 1E-5 around current point.
/// return minimum distance to points outside
/// make sure that path to current node is updated
/// get the response of tgeo
 
TGeoNode *TGeoChecker::SamplePoints(Int_t npoints, Double_t &dist, Double_t epsil, const char *g3path)
{
   TGeoNode *node = fGeoManager->FindNode();
   TGeoNode *nodegeo = nullptr;
   TGeoNode *nodeg3 = nullptr;
   TGeoNode *solg3 = nullptr;
   if (!node) {
      dist = -1;
      return nullptr;
   }
   Bool_t hasg3 = kFALSE;
   if (strlen(g3path))
      hasg3 = kTRUE;
   TString geopath = fGeoManager->GetPath();
   dist = 1E10;
   TString common = "";
   // cd to common path
   Double_t point[3];
   Double_t closest[3];
   TGeoNode *node1 = nullptr;
   TGeoNode *node_close = nullptr;
   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;
   const 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 nullptr;
         if (!nodeg3)
            return nullptr;
         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;
   }
 
   // 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
      fGeoManager->SetCurrentPoint(point[0] - eps[0] + 2 * eps[0] * gRandom->Rndm(),
                                   point[1] - eps[1] + 2 * eps[1] * gRandom->Rndm(),
                                   point[2] - eps[2] + 2 * eps[2] * gRandom->Rndm());
      node1 = fGeoManager->FindNode();
      // check if new node1 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
   fGeoManager->FindNode(point[0], point[1], point[2]); // really needed ?
   if (!node_close)
      dist = -1;
   return node_close;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Shoot one ray from start point with direction (dirx,diry,dirz). Fills input array
/// with points just after boundary crossings.
 
Double_t *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
{
   //   Int_t array_dimension = 3*dim;
   nelem = 0;
   Int_t istep = 0;
   if (!dim) {
      printf("empty input array\n");
      return array;
   }
   //   fGeoManager->CdTop();
   const 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 array;
   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 array;
         }
      }
      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 == nullptr && step > 1E10) {
            //            printf("exit : NULL endnode. nelem=%i\n", nelem);
            return array;
         }
         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 array;
            }
            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++;
      }
      fGeoManager->FindNextBoundary();
      step = fGeoManager->GetStep();
      //      printf("---next at step=%f\n", step);
      endnode = fGeoManager->Step();
      is_entering = fGeoManager->IsEntering();
   }
   return array;
   //   printf("exit : INFINITE step. nelem=%i\n", nelem);
}
 
////////////////////////////////////////////////////////////////////////////////
/// Check time of finding "Where am I" for n points.
 
void TGeoChecker::Test(Int_t npoints, Option_t *option)
{
   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;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Geometry overlap checker based on sampling.
 
void TGeoChecker::TestOverlaps(const char *path)
{
   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
   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 = nullptr;
   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;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Estimate weight of top level volume with a precision SIGMA(W)/W
/// better than PRECISION. Option can be "v" - verbose (default).
 
Double_t TGeoChecker::Weight(Double_t precision, Option_t *option)
{
   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));
   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;
}
////////////////////////////////////////////////////////////////////////////////
/// count voxel timing
 
Double_t TGeoChecker::CheckVoxels(TGeoVolume *vol, TGeoVoxelFinder *voxels, Double_t *xyz, Int_t npoints)
{
   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;
   TGeoNavigator *nav = fGeoManager->GetCurrentNavigator();
   TGeoStateInfo &td = *nav->GetCache()->GetInfo();
   timer.Start();
   for (Int_t i = 0; i < npoints; i++) {
      point = xyz + 3 * i;
      if (!shape->Contains(point))
         continue;
      checklist = voxels->GetCheckList(point, ncheck, td);
      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;
      }
   }
   nav->GetCache()->ReleaseInfo();
   time = timer.CpuTime();
   return time;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Returns optimal voxelization type for volume vol.
///  - kFALSE - cartesian
///  - kTRUE  - cylindrical
 
Bool_t TGeoChecker::TestVoxels(TGeoVolume * /*vol*/, Int_t /*npoints*/)
{
   return kFALSE;
}