void geomBuild (Bool_t v=kFALSE, Double_t fRadii = 2.5, Double_t fLength = 10.0)
{
  // world specs.
  const Double_t Wr = fRadii*2;
  const Double_t Wl = fLength*2;
	
  fWorldMat = new TGeoMaterial ("Vacuum");
  fWorldMed = new TGeoMedium ("Vacuum", 1, fWorldMat);
  fWorldVolume = gGeoManager->MakeTube ("world", fWorldMed, 
					0., Wr, Wl/2);
  gGeoManager->SetTopVolume (fWorldVolume);

  // Phantom volume
  fTargetMat = new TGeoMaterial ("Water", 0, 0, 1);
  fTargetMed = new TGeoMedium ("Water", 1, fTargetMat);	
  fTargetVolume	= gGeoManager->MakeTube ("phantom_phys", fTargetMed, 
					 0., fRadii, fLength/2);
  fTargetVolume->SetLineColor (kBlue);
  fWorldVolume->AddNode (fTargetVolume, 1, new TGeoTranslation (0., 0., 0.));

  // Voxelizing
  if (v)
    {
      sliceZ = fTargetVolume->Divide ("sliceZ", 3,50, 0, 0);
      // A divided volume never has "voxels"
      // fTargetVolume->PrintVoxels();
    }
	
  gGeoManager->CloseGeometry();
  gGeoManager->SetVisLevel(4);
  gGeoManager->SetTopVisible();
  fWorldVolume->Draw (); 
}

void DoTargetStep ()
{
  Double_t *currpos = gGeoManager->GetCurrentPoint ();
  Double_t *currdir = gGeoManager->GetCurrentDirection ();
//  Int_t nodeid = 0;
  TGeoNode *node;
  TVirtualGeoTrack *track; // Use the virtual interface

  // This computes distance to next boundary + safety within a range WITHOUT propagating the 
  // current point. You should propagate the point yourself using Step().
  // I advice you not using this when trying to propagate and cross boundaries.
//  node = gGeoManager->FindNextBoundary (99999.);
  
  // This finds the distance to next boundary within the proposed range AND
  // propagates the point. This is the new improved version of FindNextBoundary()
  // dealing with boundary crossings.
  node = gGeoManager->FindNextBoundaryAndStep(.5, kTRUE) ;
  
  // This is just used internally by the upper 2 methods (it should become protected)
  // Do not use it unless you know exactly what you are doing
  //gGeoManager->FindNextDaughterBoundary (currpos, currdir, nodeid);

  Double_t fSafeDistance = gGeoManager->GetSafeDistance ();
  Double_t fStep = gGeoManager->GetStep();

/*     // Not needed anymore
  if (fStep == 0.)
    {
      gGeoManager->SetStep (0.3);
      gGeoManager->Step (kFALSE);
      cout << "NULL STEP" << endl;
    }
  else
    gGeoManager->Step ();
*/

  Double_t *position = gGeoManager->GetCurrentPoint ();
  track = gGeoManager->GetLastTrack ();
  track->AddPoint (position[0], position[1], position[2], 0.);
//  fGeoTrack->Draw ("/*");

  if (gGeoManager->IsOnBoundary()) cout << "On boundary" << "\t";
  else cout << "SafeDistance = " << fSafeDistance << "\t";
  cout << "fStep = " << fStep << "\t";
  if (node) cout << "NodeID = " << node->GetName() << "\t";
  cout << "Next point = (" << position[0] << ", "<< position[1] << ", " << position[2] << ")" << endl;
}

void geomtracktest (UInt_t histories = 2, Bool_t v=kTRUE)
{
  geomBuild (v);
  TGeoNode *node;
  for (UInt_t i = 0; i< histories; i++)
    {
      cout << " ===========> Event: " << i  << endl;

      gGeoManager->AddTrack(i, 11);
      node = gGeoManager->InitTrack (0.,0.,0., 0.,0.,1.);
      if (node) printf("Start: %s\n", node->GetName());

      while (!gGeoManager->IsOutside ()) DoTargetStep();  
      gGeoManager->GetLastTrack ()->Draw ("/*"); // Draw the track when finished
    }
}