#include <string>
#include <iostream>
#include <sstream>

#include "globals.hh"
#include "geomdefs.hh"
#include "G4PhysicalConstants.hh"
#include "G4SystemOfUnits.hh"

#include "G4ThreeVector.hh"
#include "G4RotationMatrix.hh"
#include "G4Transform3D.hh"
#include "G4CutTubs.hh"

void CreateCutTubsChain(const G4String& Name, G4double R,
                        const std::vector<G4ThreeVector>& points,
                        std::vector<G4VSolid *>& tubs,
                        std::vector<G4Transform3D>& transforms)
{
  // Clear output vectors and check parameters
  tubs.clear();
  transforms.clear();
  G4int psize = points.size();
  if (psize < 2) return;
  for (G4int i = 0; i < psize-1; ++i) { if (points[i] == points[i+1]) return; }

  // Main loop along points
  for (G4int i = 0; i < psize-1; ++i)
  {
    // end points of current line sector
    G4ThreeVector pa = points[i];
    G4ThreeVector pb = points[i+1];
    // local basis vectors
    G4ThreeVector e3 = (pb - pa).unit();
    G4ThreeVector e2 = (e3.orthogonal()).unit();
    G4ThreeVector e1 = e2.cross(e3);
    // set rotation and transformation
    G4RotationMatrix rot(e1, e2, e3);
    transforms.push_back(G4Transform3D(rot, 0.5*(pa + pb)));
    // calculate normals
    rot.invert();
    G4ThreeVector na = rot*((i == 0) ? -e3 : ((points[i-1] - pa).unit() - e3).unit());
    G4ThreeVector nb = rot*((i == psize-2) ? e3 : ((points[i+2] - pb).unit() + e3).unit());
    // generate name of the segment
    std::stringstream buf;
    buf << Name << "_" << i + 1;
    // create cut tubs solid
    tubs.push_back(new G4CutTubs(buf.str(), 0., R, 0.5*(pb - pa).mag(), 0., CLHEP::twopi, na, nb));
  }
}
  
int main()
{
  std::vector<G4ThreeVector> points;
  points.push_back(G4ThreeVector(1,1,1));
  points.push_back(G4ThreeVector(1,1,2));
  points.push_back(G4ThreeVector(2,2,2));
  points.push_back(G4ThreeVector(3,3,1));

  std::vector<G4VSolid *> tubs;
  std::vector<G4Transform3D> transforms;
  CreateCutTubsChain("Test", 0.1, points, tubs, transforms);

  for (int i=0; i<tubs.size(); ++i)
  {
    tubs[i]->DumpInfo() ;
  }
  return 0;
}