// define useful units
static const Double_t  m = 1.;
static const Double_t cm = 1e-2*m;
static const Double_t mm = 1e-3*m;
static const Double_t um = 1e-6*m;

void geotest()
{
  TGeoManager* manager = new TGeoManager("manager", "SC");

  // Make the world
  TGeoBBox* worldbox = new TGeoBBox("worldbox", 30*m, 30*m, 30*m);
  TGeoVolume* world = new TGeoVolume("world", worldbox);
  manager->SetTopVolume(world);

  // Top volume 
  TGeoBBox* topbox = new TGeoBBox("topbox", 30*m, 30*m, 30*m);
  TGeoVolume* top = new TGeoVolume("top", topbox);

  // Parameters in Table 1
  const Double_t kk1 = -19.75*mm;
  const Double_t kk2 =  0.503*mm;

  // Dummy transparent material
  TGeoMaterial* mat = new TGeoMaterial("mat", 0, 0, 0);
  mat->SetTransparency(70); // needed in OpenGL view
  TGeoMedium* med = new TGeoMedium("med", 1, mat);

  // Make a single MAPMT photocathode
  const Double_t kPixelSize = 5.8*mm; // width of a MAPMT pixel
  const Double_t kPixelPitch = 6.08*mm; // pitch of MAPMT channels
  TGeoBBox* mapmtCathodeV = new TGeoBBox("mapmtCathodeV", kPixelSize/2., kPixelSize/2., 100*um); // very thin box
  TGeoVolume* mapmtCathode = new TGeoVolume("mapmtCathode", mapmtCathodeV);
  mapmtCathode->SetLineColor(2);

  // Make a single MAPMT
  const Double_t kMAPMTWidth = 52.0*mm;
  const Double_t kMAPMTLength = 32.7*mm; // between input window and anode pins
  const Double_t kInputWindowThickness = 1.5*mm;
  TGeoBBox* mapmtV = new TGeoBBox("mapmtV", kMAPMTWidth/2., kMAPMTWidth/2.,
				  kMAPMTLength/2.);
  TGeoVolume* mapmt = new TGeoVolume("mapmt", mapmtV, med);
  //  mapmt->SetVisibility(0);
  TGeoBBox* mapmtInputWindowV = new TGeoBBox("mapmtInputWindowV",
					     kMAPMTWidth/2., kMAPMTWidth/2.,
					     kInputWindowThickness/2.);

  TGeoVolume* mapmtInputWindow = new TGeoVolume("mapmtInputWindow", mapmtInputWindowV, med);
  mapmtInputWindow->SetLineColor(4);
  mapmt->AddNodeOverlap(mapmtInputWindow, 1, new TGeoTranslation(0, 0, kMAPMTLength/2. - kInputWindowThickness/2.));

  Int_t n = 1;
  for(Int_t i = 0; i < 8; i++){
    Double_t dx = (i - 3.5)*kPixelPitch;
    for(Int_t j = 0; j < 8; j++){
      Double_t dy = (j - 3.5)*kPixelPitch;
      mapmt->AddNodeOverlap(mapmtCathode, n, new TGeoTranslation(dx, dy, kMAPMTLength/2. - kInputWindowThickness - 100*um));
      n++;
    } // j
  } // i

  TGeoBBox* mapmtBackObsV = new TGeoBBox("mapmtBackObsV",
					kMAPMTWidth/2., kMAPMTWidth/2.,
					15*mm);
  TGeoVolume* mapmtBackObs = new TGeoVolume("mapmtBackObs", mapmtBackObsV);
  mapmt->AddNode(mapmtBackObs, 1, new TGeoTranslation(0, 0, -kMAPMTLength/2. + 15*mm));

  // Make the focal plane
  n = 1;
  for(Int_t i = -7; i <= 7; i++){
    Double_t dx = i*kMAPMTWidth;
    for(Int_t j = -7; j <= 7; j++){
      Double_t dy = j*kMAPMTWidth;
      Double_t dz = kk1*((i/7.)*(i/7.) + (j/7.)*(j/7.)) +
	kk2*((i/7.)*(i/7.) + (j/7.)*(j/7.))*((i/7.)*(i/7.) + (j/7.)*(j/7.));
      top->AddNodeOverlap(mapmt, n, new TGeoTranslation(dx, dy, kMAPMTLength/2. + dz));
      //      if(n == 151) goto skip;
      n++;
    } // y
  } // x
 skip:
  world->AddNode(top, 1);

  manager->CloseGeometry();

  TCanvas* canGeometry = new TCanvas("canGeometry", "canGeometry", 800, 800);
  top->Draw();
}