void cube()
{
// Cube example
   TGeoManager *geom = new TGeoManager("simple", "A cube");
   TGeoMaterial *matVacuum = new TGeoMaterial("Vacuum", 0,0,0);
   TGeoMedium *Vacuum = new TGeoMedium("Vacuum",1, matVacuum);
   // The world
   TGeoVolume *top = geom->MakeBox("WORLD", Vacuum, 1000., 1000., 1000.);
   geom->SetTopVolume(top);
   // And cube volume in the center
   TGeoVolume *cube = geom->MakeBox("CUBE", Vacuum,10.,10.,10.);
   cube->SetLineColor(kBlue);
   top->AddNode(cube,0);
   geom->CloseGeometry();
   // Draw the world
   top->Draw();
   new TRandom3();
   shoot();
}

void shoot()
{   
   // Shoot a track towards the cube (from outside the cube to the center)
   Double_t point[3] = {0.,0.,0.};
   Double_t dir[3];
   Int_t i;
   TPolyLine3D *pl = new TPolyLine3D(2);
   pl->SetLineColor(kRed);
   TGeoVolume *cube = gGeoManager->GetVolume("CUBE");
   while (cube->Contains(point)) {
      for (i=0; i<3; i++) point[i] = -20 + 40.*gRandom->Rndm();
   }
   Double_t r = TMath::Sqrt(point[0]*point[0]+point[1]*point[1]+point[2]*point[2]);
   for (i=0; i<3; i++) dir[i] = -point[i]/r;
   printf("Shooting track from : (%g, %g, %g), direction: (%g, %g, %g)\n",
          point[0],point[1],point[2],dir[0],dir[1],dir[2]);
   gGeoManager->InitTrack(point,dir);
   gGeoManager->DrawCurrentPoint();
   // Compute distance to entering
   Double_t length = 0.;
   gGeoManager->FindNextBoundaryAndStep();
   Double_t step = gGeoManager->GetStep();
   length += step;
   const Double_t *crt = gGeoManager->GetCurrentPoint();
   pl->SetNextPoint(crt[0],crt[1],crt[2]);
   printf("Entring cube at: (%g, %g, %g) length=%g\n", crt[0],crt[1],crt[2],length);
   // Now exit
   gGeoManager->FindNextBoundaryAndStep();  
   step = gGeoManager->GetStep();
   length += step;
   pl->SetNextPoint(crt[0],crt[1],crt[2]);
   printf("Exiting cube at: (%g, %g, %g) length=%g\n", crt[0],crt[1],crt[2],length);
   pl->Draw("same");
   printf("Type shoot() to shoot again...\n");
}