/********************************* constants *********************************/

#define pi 3.14159265358979

/****************************** global variables *****************************/

int n_nodes = 0;

TGeoVolume *LArg;

/****************************** auxiliar functions ***************************/

Double_t eta2rad(Double_t eta)
{
	return pi/2.0 - acos(tanh(eta));
}

void BuilLArg()
{
	// variables

	Double_t v[16], dx1, dx2, dx;
	Double_t phi;

	// materials

	LArgMaterial = new TGeoMaterial("LArg",0,0,0);
	CellMaterial = new TGeoMaterial("Cell",0,0,0);

	LArgMedium = new TGeoMedium  ("LArg",1,LArgMaterial);
	CellMedium = new TGeoMedium  ("Cell",1,CellMaterial);

	LArgMaterial->SetTransparency(80);

	// top volume

	LArg = gGeoManager->MakePcon("LArg", LArgMedium, 0.0, 360, 2);
	TGeoPcon   *pcon = (TGeoPcon*)(LArg->GetShape());
	pcon->DefineSection(0,-4000,1500,1970);
	pcon->DefineSection(1,+4000,1500,1970);
	gGeoManager->SetNsegments(64);

	//TGeoVolume *LArg = gGeoManager->MakeTube("LArg", LArgMedium, 1500.0, 1970.0, 3400.0);
   // rotations
   TObjArray *rot = new TObjArray(256);

	// external camade

   TGeoVolumeAssembly *extcam = new TGeoVolumeAssembly("ExtCam");
	Double_t h1   = 1927.84753363;
	Double_t h2   = 1970;
	Double_t deta = 0.05;
	Double_t dphi = 2.0*pi/128.0;

	Double_t dy1  = h1*tan(dphi/2);
	Double_t dy2  = h2*tan(dphi/2);

	Double_t eta = -1.350; // maximum angle in eta
   phi = dphi/2 + pi/2;
	for (int j = 0; j < 128; j++) {
		rot->AddAt(new TGeoRotation("rot", -phi*360.0/(2*pi), 90.0, 90.0),j);
		phi += dphi;
   }      
	for (int i = 0; i < 54; i++)
	{

		dx1 = h1/tan(pi/2 - eta2rad(eta+deta)) - h1/tan(pi/2 - eta2rad(eta));
		dx2 = h2/tan(pi/2 - eta2rad(eta+deta)) - h2/tan(pi/2 - eta2rad(eta));
		dx  = (h2-h1)/tan(pi/2 - eta2rad(eta));

		v[ 0] =  0.0; v[ 1] = -dy1;
		v[ 2] = -dx1; v[ 3] = -dy1;
		v[ 4] = -dx1; v[ 5] = +dy1;
		v[ 6] =  0.0; v[ 7] = +dy1;

		v[ 8] = -dx;     v[ 9] = -dy2;
		v[10] = -dx-dx2; v[11] = -dy2;
		v[12] = -dx-dx2; v[13] = +dy2;
		v[14] = -dx;     v[15] = +dy2;
      TGeoVolumeAssembly *etaslice = new TGeoVolumeAssembly(Form("extEta%d",i));
   	TGeoVolume *cell = gGeoManager->MakeArb8("cell", CellMedium, (h2-h1)/2, v);
		cell->SetLineColor(kRed);


		phi = dphi/2 + pi/2;
		for (int j = 0; j < 128; j++)
		{

			TGeoCombiTrans *comb = new TGeoCombiTrans((-h1-(h2-h1)/2)*sin(phi), (-h1-(h2-h1)/2)*cos(phi), -h1/tan(pi/2 - eta2rad(eta)), (TGeoRotation*)rot->At(j));

			etaslice->AddNode(cell, j, comb);

			phi += dphi;
		}
		extcam->AddNode(etaslice, i);
		eta += deta;
	}
   LArg->AddNode(extcam, 1);

	// medium camade
   TGeoVolumeAssembly *medcam = new TGeoVolumeAssembly("MedCam");

	Double_t h1   = 1590.62780269;
	Double_t h2   = 1927.84753363;
	Double_t deta = 0.025;
	Double_t dphi = 2.0*pi/256.0;

	Double_t dy1  = h1*tan(dphi/2);
	Double_t dy2  = h2*tan(dphi/2);

	Double_t eta = -1.450; // maximum angle in eta
   phi = dphi/2 + pi/2;
	for (int j = 0; j < 256; j++) {
		rot->AddAt(new TGeoRotation("rot",  - phi*360.0/(2*pi), 90.0, 90.0),j);
		phi += dphi;
   }      
	for (int i = 0; i < 116; i++)
	{

		dx1 = h1/tan(pi/2 - eta2rad(eta+deta)) - h1/tan(pi/2 - eta2rad(eta));
		dx2 = h2/tan(pi/2 - eta2rad(eta+deta)) - h2/tan(pi/2 - eta2rad(eta));
		dx  = (h2-h1)/tan(pi/2 - eta2rad(eta));

		v[ 0] =  0.0; v[ 1] = -dy1;
		v[ 2] = -dx1; v[ 3] = -dy1;
		v[ 4] = -dx1; v[ 5] = +dy1;
		v[ 6] =  0.0; v[ 7] = +dy1;

		v[ 8] = -dx;     v[ 9] = -dy2;
		v[10] = -dx-dx2; v[11] = -dy2;
		v[12] = -dx-dx2; v[13] = +dy2;
		v[14] = -dx;     v[15] = +dy2;
      TGeoVolumeAssembly *etaslice = new TGeoVolumeAssembly(Form("medEta%d",i));
   	TGeoVolume *cell = gGeoManager->MakeArb8("cell", CellMedium, (h2-h1)/2, v);
		cell->SetLineColor(kBlue);

		phi = dphi/2 + pi/2;
		for (int j = 0; j < 256; j++)
		{

			TGeoCombiTrans *comb = new TGeoCombiTrans((-h1-(h2-h1)/2)*sin(phi), (-h1-(h2-h1)/2)*cos(phi), -h1/tan(pi/2 - eta2rad(eta)),(TGeoRotation*) rot->At(j));

			etaslice->AddNode(cell, j, comb);

			phi += dphi;
		}
		medcam->AddNode(etaslice, i);
		eta += deta;
	}
   LArg->AddNode(medcam, 1);

	// inner camade
   TGeoVolumeAssembly *incam = new TGeoVolumeAssembly("InCam");

	Double_t h1   = 1500;
	Double_t h2   = 1590;
	Double_t deta = 0.025/8;
	Double_t dphi = 2.0*pi/64.0;

	Double_t dy1  = h1*tan(dphi/2);
	Double_t dy2  = h2*tan(dphi/2);

	//Double_t eta = -1.4640625; // maximum angle in eta
	Double_t eta = -1.45; // maximum angle in eta
   phi = dphi/2 + pi/2;
	for (int j = 0; j < 64; j++) {
		rot->AddAt(new TGeoRotation("rot",  - phi*360.0/(2*pi), 90.0, 90.0),j);
		phi += dphi;
   }      
	for (int i = 0; i < 928; i++)
	{

		dx1 = h1/tan(pi/2 - eta2rad(eta+deta)) - h1/tan(pi/2 - eta2rad(eta));
		dx2 = h2/tan(pi/2 - eta2rad(eta+deta)) - h2/tan(pi/2 - eta2rad(eta));
		dx  = (h2-h1)/tan(pi/2 - eta2rad(eta));

		v[ 0] =  0.0; v[ 1] = -dy1;
		v[ 2] = -dx1; v[ 3] = -dy1;
		v[ 4] = -dx1; v[ 5] = +dy1;
		v[ 6] =  0.0; v[ 7] = +dy1;

		v[ 8] = -dx;     v[ 9] = -dy2;
		v[10] = -dx-dx2; v[11] = -dy2;
		v[12] = -dx-dx2; v[13] = +dy2;
		v[14] = -dx;     v[15] = +dy2;
      TGeoVolumeAssembly *etaslice = new TGeoVolumeAssembly(Form("intEta%d",i));
		TGeoVolume *cell = gGeoManager->MakeArb8("cell", CellMedium, (h2-h1)/2, v);
		cell->SetLineColor(kYellow);

		phi = dphi/2 + pi/2;
		for (int j = 0; j < 64; j++)
		{

			TGeoCombiTrans *comb = new TGeoCombiTrans((-h1-(h2-h1)/2)*sin(phi), (-h1-(h2-h1)/2)*cos(phi), -h1/tan(pi/2 - eta2rad(eta)),(TGeoRotation*) rot->At(j));

			etaslice->AddNode(cell, j, comb);

			phi += dphi;
		}
		incam->AddNode(etaslice, i);

		cout << i << " " << eta << endl;

		eta += deta;
	}
   delete rot;
   LArg->AddNode(incam, 1);
}

void CaloBuild()
{
	BuilLArg();

	gGeoManager->SetTopVolume(LArg);
	gGeoManager->CloseGeometry();
	gGeoManager->SetVisLevel(3);

	gGeoManager->Export("CaloGeometry.root");
}