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);
      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;

  //gGeoManager->FindNextBoundary (99999.);
  //gGeoManager->FindNextBoundaryAndStep(.5, kTRUE) ;
  gGeoManager->FindNextDaughterBoundary (currpos, currdir, nodeid);

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

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

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

  cout << "SafeDistance = " << fSafeDistance << "\t";
  cout << "fStep = " << fStep << "\t";
  cout << "NodeID = " << nodeid << "\t";
  cout << "Next point = (" << position[0] << ", "<< position[1] << ", " << position[2] << ")" << endl;
}

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

      gGeoManager->AddTrack(i, 11);
      gGeoManager->InitTrack (0.,0.,0., 0.,0.,1.);

      while (!gGeoManager->IsOutside ()) 
	DoTargetStep();  
    }
}