#include "TGeoManager.h"
#include "TROOT.h"
#include "TStyle.h"
#include "TMath.h"

using namespace TMath;

void PMTs()
{
  TGeoManager *geom = new TGeoManager("PMTs", "inner PMT array");

  //--- define some materials
  TGeoMaterial *matVacuum = new TGeoMaterial("Vacuum", 0,0,0);
  TGeoMaterial *matAl = new TGeoMaterial("Al", 26.98,13,2.7);

  //--- define some media
  TGeoMedium *Vacuum = new TGeoMedium("Vacuum",1, matVacuum);
  TGeoMedium *Al = new TGeoMedium("Root Material",2, matAl);

  //--- make one hexagonal PMT
  TGeoVolume *HexPMT = geom->MakePgon("Hexagonal PMT", Al, 0., 360., 6, 2);
  // DefineSection(sectionID,z,innerR, RofInscribedCircle)
  ((TGeoPgon*)(HexPMT->GetShape()))->DefineSection(0,-0.01,0,5.84/2.);
  ((TGeoPgon*)(HexPMT->GetShape()))->DefineSection(1, 0.01,0,5.84/2.);
  HexPMT->SetLineColor(kBlue);

  //--- make one round PMT
  TGeoVolume *roundPMT = geom->MakeTube("Round PMT", Al, 0., 4.7/2., 0.1);
  roundPMT->SetLineColor(kRed);

  //--- dimentions from technical drawings. Unit: cm
  Double_t lb = 31.03165; // length of the bottom of the triangular PMT holder
  Double_t lh = 23.15; // length of the height of the triangular PMT holder
  Double_t ls = Sqrt(lh*lh+lb/2.*lb/2.); // length of the side of the holder
  Double_t pz = .1; // FIXME:z position of PMT surface w.r.o. holder's inner surface

  // angle between triangular holder and the combined pentagon
  Double_t thetaTP = ACos(0.5*lb/lh/Tan(Pi()/5.));

  // distance between the pentagon center and the triangle top
  Double_t dc = lh*Sin(thetaTP);
  Double_t dv = 5.94; // distance of the centers of two PMTs vertically
  Double_t dh = 5.35; // distance of the centers of two PMTs horizontally
  Double_t dt = 5.57; // distance of the top PMT to the holder top
  Double_t ds =14.36; // distance of the side PMTs to the holder top

  //--- Make the assembly of PMTs in one triangle hoder
  TGeoVolume *triangle = new TGeoVolumeAssembly("Triangular Holder");
  triangle->AddNode(HexPMT, 1, new TGeoTranslation(0,-dt,       -pz));
  triangle->AddNode(HexPMT, 2, new TGeoTranslation(0,-dt-dv,    -pz));
  triangle->AddNode(HexPMT, 3, new TGeoTranslation(0,-dt-dv-dv, -pz));
  triangle->AddNode(HexPMT, 4, new TGeoTranslation(-dh, -ds,    -pz));
  triangle->AddNode(HexPMT, 5, new TGeoTranslation( dh, -ds,    -pz));
  triangle->AddNode(HexPMT, 6, new TGeoTranslation(-dh, -ds-dv, -pz));
  triangle->AddNode(HexPMT, 7, new TGeoTranslation( dh, -ds-dv, -pz));
  //triangle->AddNode(roundPMT, 2, new TGeoTranslation(0, 0, 0));

  //--- Make the assembly of one pentagon
  TGeoVolume *pentagon = new TGeoVolumeAssembly("Pentagon");
  for (Int_t i=0; i<5; i++)
    pentagon->AddNode(triangle, i+1, new TGeoCombiTrans(0,0,dc, 
          new TGeoRotation(Form("rTri%d",i+1), 360./5.*i, thetaTP/Pi()*180., 0)));

  //--- Make the assembly of half sphere
  TGeoVolume *halfSphere = new TGeoVolumeAssembly("Half Sphere");

  Double_t gr = (1.+Sqrt(5.))/2.; // golden ratio
  // half edge of the dodecahedron's dual icosahedron
  Double_t he = lb/2./Sqrt(5.)*gr*gr; 

  TGeoTranslation tra1(0., he, he*gr);
  TGeoTranslation tra2(0.,-he, he*gr);
  TGeoTranslation tra3(he*gr, 0., he);
  TGeoTranslation tra4(-he*gr,0., he);
  TGeoTranslation tra5(he,-he*gr, 0.);
  TGeoTranslation tra6(-he,-he*gr,0.);

  TGeoRotation rot1("rot1"); 
  rot1.RotateX(-ATan(gr/(1.+gr))*180./Pi()); // z -> y

  TGeoRotation rot2("rot2"); 
  rot2.RotateZ(180.); // x -> y
  rot2.RotateX(ATan(gr/(1.+gr))*180./Pi()); // y -> z

  TGeoRotation rot3("rot3"); 
  rot3.RotateZ(90.); // x -> y
  rot3.RotateY(ATan((1.+gr)/gr)*180./Pi()); // z -> x

  TGeoRotation rot4("rot4"); 
  rot4.RotateZ(-90.); // y -> x
  rot4.RotateY(-ATan((1.+gr)/gr)*180./Pi()); // x -> z

  TGeoRotation rot5("rot5"); 
  rot5.RotateY(90.); // z -> x
  rot5.RotateZ(-ATan((1.+gr)/gr)*180./Pi()); // y -> x

  TGeoRotation rot6("rot6"); 
  rot6.RotateY(-90.); // x -> z
  rot6.RotateZ(ATan((1.+gr)/gr)*180./Pi()); // x -> y

  TGeoCombiTrans *com1 = new TGeoCombiTrans(tra1, rot1);
  TGeoCombiTrans *com2 = new TGeoCombiTrans(tra2, rot2);
  TGeoCombiTrans *com3 = new TGeoCombiTrans(tra3, rot3);
  TGeoCombiTrans *com4 = new TGeoCombiTrans(tra4, rot4);
  TGeoCombiTrans *com5 = new TGeoCombiTrans(tra5, rot5);
  TGeoCombiTrans *com6 = new TGeoCombiTrans(tra6, rot6);

  halfSphere->AddNode(pentagon, 1, com1);
  halfSphere->AddNode(pentagon, 2, com2);
  halfSphere->AddNode(pentagon, 3, com3);
  halfSphere->AddNode(pentagon, 4, com4);
  halfSphere->AddNode(pentagon, 5, com5);
  halfSphere->AddNode(pentagon, 6, com6);

  //--- Make the assembly of the whole Sphere
  TGeoVolume *sphere = new TGeoVolumeAssembly("sphere");
  sphere->AddNode(halfSphere, 1);
  TGeoRotation *rS = new TGeoRotation("rotatePiAlongX"); rS->RotateX(180.);
  sphere->AddNode(halfSphere, 2, rS);

  TGeoVolume *inner = geom->MakeBox("Inner Detector", Vacuum, 1000., 1000., 1000.);
  geom->SetTopVolume(inner);

  inner->AddNode(sphere, 1, new TGeoTranslation(0., 0., 0.));

  //--- close the geometry
  // MT: request node-ids to be calculated.
  geom->CloseGeometry("i");

  //--- draw the ROOT box.
  // by default the picture will appear in the standard ROOT TPad.
  //if you have activated the following line in system.rootrc, 
  //it will appear in the GL viewer
  //#Viewer3D.DefaultDrawOption:   ogl

  // geom->SetVisLevel(18);
  // inner->Draw(); // use ROOT Browser
  // sphere->Draw("ogl"); // use OpenGL
}