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

#define pi 3.1415926535897932384626433832795

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

TGeoVolume *Calo;
TGeoVolume *LArg;
TGeoVolume *Tile;

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

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

void BuildLArg()
{
	// variables

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

   new TGeoManager("calorimeter", "ATLAS calorimeter");
	TObjArray *rot = new TObjArray(256);

   printf("Building LArg...\n");
	// 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,-3400,1500,1970);
	pcon->DefineSection(1,+3400,1500,1970);
	gGeoManager->SetNsegments(64);

	// external camade
   printf("   ExtCam...\n");

	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/256.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 < 256; j++) 
	{
		rot->AddAt(new TGeoRotation("rot", 180+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;

   	TGeoArb8 *arb8 = new TGeoArb8("arb8", (h2-h1)/2.0, v);
		TGeoBBox *box  = new TGeoBBox("box", 3400, 3400, 3400);
		TGeoVolumeAssembly *etaslice = new TGeoVolumeAssembly(Form("extEta%d",i));

		phi = dphi/2 + pi/2;
		for (int j = 0; j < 256; j++)
		{
			TGeoVolume *cell = gGeoManager->MakeArb8("cell", CellMedium, (h2-h1)/2.0, v);
			cell->SetLineColor(kRed);

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

			if (i == 0 || i == 53)
			{
				comb->RegisterYourself();

            TGeoIntersection *boolnode = new TGeoIntersection(arb8, box, comb, NULL);
            TGeoCompositeShape *cs   = new TGeoCompositeShape("cs", boolnode);
            // The constructor below looks by name for all components. This is slow in
            // a loop and it even fails in your case since it takes always the first "comb"
//				TGeoCompositeShape *cs   = new TGeoCompositeShape("cs","arb8:comb * box");	
				TGeoVolume         *cell2 = new TGeoVolume("cell2",cs);

				etaslice->AddNode(cell2, j, new TGeoTranslation(0,0,0));
			}
			else
			{
				etaslice->AddNode(cell, j, comb);
			}

			phi += dphi;
		}

		extcam->AddNode(etaslice, i);
		eta += deta;
	}

	LArg->AddNode(extcam, 1);

	// medium camade
   printf("   MedCam...\n");

	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
	Double_t eta  = -1.50 +0.025;

	phi = dphi/2 + pi/2;
	for (int j = 0; j < 256; j++) 
	{
		rot->AddAt(new TGeoRotation("rot", 180+phi*360.0/(2*pi), 90.0, -90.0),j);
		phi += dphi;
	}
	
	for (int i = 0; i < 116; i++)
	{
		if (i == 0 || i == 115) 
			deta = 3*0.025; 
		else if (i == 1 || i == 114) deta = 0;
		else deta = 0.025;

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

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

			TGeoVolume *cell = gGeoManager->MakeArb8("cell", CellMedium, (h2-h1)/2, v);
			cell->SetLineColor(kBlue);

			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
   printf("   InCam...\n");

	TGeoVolumeAssembly *incam = new TGeoVolumeAssembly("InCam");

	Double_t h1   = 1500;
	Double_t h2   = 1589;
	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.50 + 0.025;
	//Double_t eta  = -1.46562492847442630;
	//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", 180+phi*360.0/(2*pi), 90.0, -90.0),j);
		phi += dphi;
	}
	
	for (int i = 0; i < 938; i++)
	{

		if (i > 920 || i < 17)
		{
			if (i == 921)
			{
				eta -= 0.025/2.0;
			}

			if (i == 921 || i == 929 || i == 937 || i == 0 || i == 8 || i == 16)
			{
				deta = 0.025;
			}
			else {deta = 0.0;}
		}
		else {deta = 0.025/8.0;}

		if (i == 17) eta -= 0.025/2.0;

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

		phi = dphi/2 + pi/2;
		for (int j = 0; j < 64; j++)
		{
			TGeoVolume *cell = gGeoManager->MakeArb8("cell", CellMedium, (h2-h1)/2, v);
			cell->SetLineColor(kYellow);

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

		eta += deta;
	}

	delete rot;

	LArg->AddNode(incam, 1);
}

void BuildTile()
{
	// constants

	const float cell_d_w[07] = {432.245, 411.442, 383.705, 365.212, 383.705, 411.442, 432.245};
	const float cell_d_t[07] = {+2387.750, +1544.063, +748.916, 0.0, -748.916, -1544.063, -2387.750};

	const float cell_c_w[16] = {185.120, 197.462, 185.120, 182.035, 169.694, 166.608, 163.523, 160.438,
		                        160.438, 163.523, 166.608, 169.694, 182.035, 185.120, 197.462, 185.120};
	const float cell_c_t[16] = {+2634.880, +2252.298, +1869.716, +1502.561, +1150.832, +814.530, +484.399, +160.438,
	                            -160.438, -484.399, -814.530, -1150.832, -1502.561, -1869.716, -2252.298, -2634.880};

	const float cell_b_w[18] = {166.245, 181.638, 169.323, 163.167, 157.009, 147.772, 141.616, 144.694, 138.537,
	                            138.537, 144.694, 141.616, 147.772, 157.009, 163.167, 169.323, 181.638, 166.245};
	const float cell_b_t[18] = {+2653.757, +2305.874, +1954.913, +1622.423, +1302.247, +997.466, +708.078, +421.768, +138.537,
	                            -138.537, -421.768, -708.078, -997.466, -1302.247, -1622.423, -1954.913, -2305.874, -2653.757};

	const float cell_a_w[20] = {148.096, 169.694, 157.352, 154.267, 141.925, 135.755, 129.584, 129.584, 123.413, 120.328,
	                            120.328, 123.413, 129.584, 129.584, 135.755, 141.925, 154.267, 157.352, 169.694, 148.096};
	const float cell_a_t[20] = {+2671.900, +2354.110, +2027.064, +1715.445, +1419.253, +1141.573, +876.234, +617.066, +364.069, +120.328,
	                            -120.328, -364.069, -617.066, -876.234, -1141.573, -1419.253, -1715.445, -2027.064, -2354.110, -2671.900}

	// variables

	Double_t phi;

   printf("Building Tile...\n");
	TObjArray *rot = new TObjArray(256);

	TGeoTrd1 *cellb[18][64];
	TGeoTrd1 *cellc[16][64];

	TGeoCombiTrans *combb[18][64];
	TGeoCombiTrans *combc[16][64];

	// materials

	TileMaterial = new TGeoMaterial("Tile",0,0,0);
	CellMaterial = new TGeoMaterial("Cell",0,0,0);

	TileMedium = new TGeoMedium("Tile",1,TileMaterial);
	CellMedium = new TGeoMedium("Cell",1,CellMaterial);

	TileMaterial->SetTransparency(80);

	// top volume

	Tile = gGeoManager->MakePcon("Tile", TileMedium, 0.0, 360, 2);
	TGeoPcon   *pcon = (TGeoPcon*)(Tile->GetShape());
	pcon->DefineSection(0,-6110.0,2280,3886.352);
	pcon->DefineSection(1,+6110.0,2280,3886.352);
	gGeoManager->SetNsegments(64);

	// external camade

   printf("   ExtCam...\n");
	TGeoVolumeAssembly *extcam = new TGeoVolumeAssembly("ExtCam");

	Double_t h1   = 3396.702;
	Double_t h2   = 3886.352;
	Double_t deta = 0.2;
	Double_t dphi = 2.0*pi/64.0;

	Double_t eta  = -0.7; // maximum angle in eta

	Double_t phi = dphi/2 - pi/2;
	for (int j = 0; j < 64; j++) 
	{
		rot->AddAt(new TGeoRotation("rot", phi*360.0/(2*pi), 90.0, 0.0),j);
		phi += dphi;
	}
	
	for (int i = 0; i < 7; i++)
	{
		TGeoVolumeAssembly *etaslice = new TGeoVolumeAssembly(Form("extEta%d",i));

		phi = dphi/2 - pi/2;
		for (int j = 0; j < 64; j++)
		{
			TGeoVolume *cell = gGeoManager->MakeTrd1("cell", CellMedium, 165.491, 188.965, cell_d_w[i], 244.825);
			cell->SetLineColor(kRed);

			TGeoCombiTrans *comb = new TGeoCombiTrans("comb", -(-h1-(h2-h1)/2)*sin(phi), (-h1-(h2-h1)/2)*cos(phi), -cell_d_t[i], (TGeoRotation*)rot->At(j));
			etaslice->AddNode(cell, j, comb);

			phi += dphi;
		}

		extcam->AddNode(etaslice, i);
		eta += deta;
	}

	Tile->AddNode(extcam, 1);

	// medium camade out

   printf("   MedCamOut...\n");
	TGeoVolumeAssembly *medcam = new TGeoVolumeAssembly("MedCamOut");

	Double_t h1 = 3396.702 - 2*216.797;
	Double_t h2 = 3396.702;
	Double_t deta = 0.1;
	Double_t dphi = 2.0*pi/64.0;

	Double_t eta  = -0.8; // 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, 0.0),j);
		phi += dphi;
	}
	
	for (int i = 0; i < 16; i++)
	{
		TGeoVolumeAssembly *etaslice = new TGeoVolumeAssembly(Form("medEta%d",i));

		phi = dphi/2 - pi/2;
		for (int j = 0; j < 64; j++)
		{
			cellc[i][j] = new TGeoTrd1(Form("medc_%d_%d",i,j), 144.704, 165.491, cell_c_w[i], 216.797);
			TGeoVolume *cell = gGeoManager->MakeTrd1("cell", CellMedium, 144.704, 165.491, cell_c_w[i], 216.797);
			cell->SetLineColor(kRed);

			TGeoCombiTrans *comb = new TGeoCombiTrans(Form("combc_%d_%d",i,j), -(-h1-(h2-h1)/2)*sin(phi), (-h1-(h2-h1)/2)*cos(phi), -cell_c_t[i], (TGeoRotation*)rot->At(j));
			etaslice->AddNode(cell, j, comb);
			combc[i][j] = comb;

			phi += dphi;
		}

		medcam->AddNode(etaslice, i);
		eta += deta;
	}

	Tile->AddNode(medcam, 1);

	// medium camade in

   printf("   MedCamIn...\n");
	TGeoVolumeAssembly *medcam = new TGeoVolumeAssembly("MedCamIn");

	Double_t h1 = 3396.702 - 2*216.797 - 2*188.431;
	//Double_t h2 = 3396.702 - 2*216.797;
	Double_t h2 = 3396.702 - 2*216.8;
	Double_t deta = 0.1;
	Double_t dphi = 2.0*pi/64.0;

	Double_t eta  = -0.9; // 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, 0.0),j);
		phi += dphi;
	}
	
	for (int i = 0; i < 18; i++)
	{
		TGeoVolumeAssembly *etaslice = new TGeoVolumeAssembly(Form("medEta%d",i));

		phi = dphi/2 - pi/2;
		for (int j = 0; j < 64; j++)
		{
			cellb[i][j] = new TGeoTrd1(Form("medb_%d_%d",i,j),126.638, 144.704, cell_b_w[i], 188.431);
			TGeoVolume *cell = gGeoManager->MakeTrd1("cell", CellMedium, 126.638, 144.704, cell_b_w[i], 188.431);
			cell->SetLineColor(kRed);

			TGeoCombiTrans *comb = new TGeoCombiTrans(Form("combb_%d_%d",i,j), -(-h1-(h2-h1)/2)*sin(phi), (-h1-(h2-h1)/2)*cos(phi), -cell_b_t[i], (TGeoRotation*)rot->At(j));
			etaslice->AddNode(cell, j, comb);
			combb[i][j] = comb;

			phi += dphi;
		}

		medcam->AddNode(etaslice, i);
		eta += deta;
	}

	Tile->AddNode(medcam, 2);

	/*// joing medium camade

	Double_t h1 = 3396.702 - 2*216.797 - 2*188.431;
	Double_t h2 = 3396.702 - 2*216.797;
	Double_t h3 = 3396.702;
	Double_t dphi = 2.0*pi/64.0;

	TGeoVolumeAssembly *medcam = new TGeoVolumeAssembly("MedCam");

	for (int i = 0; i < 16; i++)
	{
		TGeoVolumeAssembly *etaslice = new TGeoVolumeAssembly(Form("medjEta%d",i));

		phi = dphi/2 - pi/2;
		for (int j = 0; j < 64; j++)
		{
			TGeoCombiTrans *tr1 = new TGeoCombiTrans(Form("tr1_%d_%d",i+1,j), 0, (-h1-(h2-h1)/2.0), -cell_b_t[i+1], new TGeoRotation("rot", 0, 90.0, 0.0));
			tr1->RegisterYourself();
			TGeoCombiTrans *tr2 = new TGeoCombiTrans(Form("tr2_%d_%d",i,j), 0, (-h2-(h3-h2)/2.0)-0.1, -cell_c_t[i], new TGeoRotation("rot", 0, 90.0, 0.0));
			tr2->RegisterYourself();

			TGeoCompositeShape *cs = new TGeoCompositeShape("cs",Form("medb_%d_%d:tr1_%d_%d + medc_%d_%d:tr2_%d_%d", i+1,j,i+1,j,i,j,i,j));
			TGeoVolume *cell = new TGeoVolume("cell", cs);
			etaslice->AddNode(cell, j, new TGeoRotation("rot", 0.0, 0.0, phi*360.0/(2*pi)));

			phi += dphi;
		}

		medcam->AddNode(etaslice, i);
	}
	
	Tile->AddNode(medcam, 3);*/

	// joing medium camade

  // Not clear what you are doing here... You make unions with objects already positioned

   printf("   MedCamJnt...\n");
	TGeoVolumeAssembly *medcam = new TGeoVolumeAssembly("MedCamJnt");

	for (int i = 0; i < 16; i++)
	{
		TGeoVolumeAssembly *etaslice = new TGeoVolumeAssembly(Form("medjEta%d",i));

		for (int j = 0; j < 64; j++)
		{
         TGeoUnion *boolnode = new TGeoUnion(cellb[i+1][j], cellb[i][j], combb[i+1][j], combb[i][j]);
//			TGeoCompositeShape *cs = new TGeoCompositeShape("cs",Form("medb_%d_%d:combb_%d_%d + medc_%d_%d:combc_%d_%d", i+1,j,i+1,j,i,j,i,j));
         TGeoCompositeShape *cs = new TGeoCompositeShape("cs",boolnode);
			TGeoVolume *cell = new TGeoVolume("cell", cs);
			etaslice->AddNode(cell, j);
		}

		medcam->AddNode(etaslice, i);
	}
	
	Tile->AddNode(medcam, 3);

	// internal camade

   printf("   IntCam...\n");
	TGeoVolumeAssembly *intcam = new TGeoVolumeAssembly("IntCam");

	Double_t h1   = 3396.702 - 2*216.797 - 2*188.431 - 2*151.623;
	Double_t h2   = 3396.702 - 2*216.797 - 2*188.431;
	Double_t deta = 0.1;
	Double_t dphi = 2.0*pi/64.0;

	Double_t eta  = -1.0; // 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, 0.0),j);
		phi += dphi;
	}
	
	for (int i = 0; i < 20; i++)
	{
		TGeoVolumeAssembly *etaslice = new TGeoVolumeAssembly(Form("intEta%d",i));

		phi = dphi/2 - pi/2;
		for (int j = 0; j < 64; j++)
		{
			TGeoVolume *cell = gGeoManager->MakeTrd1("cell", CellMedium, 112.100, 126.638, cell_a_w[i], 151.623);
			cell->SetLineColor(kRed);

			TGeoCombiTrans *comb = new TGeoCombiTrans("comb", -(-h1-(h2-h1)/2)*sin(phi), (-h1-(h2-h1)/2)*cos(phi), -cell_a_t[i], (TGeoRotation*)rot->At(j));
			etaslice->AddNode(cell, j, comb);

			phi += dphi;
		}

		intcam->AddNode(etaslice, i);
		eta += deta;
	}

	Tile->AddNode(intcam, 4);
}

void CaloBuild()
{
	BuildLArg();
	BuildTile();

	// materials

	CaloMaterial = new TGeoMaterial("Calo",0,0,0);
	CaloMedium   = new TGeoMedium  ("Calo",1,CaloMaterial);

	CaloMaterial->SetTransparency(80);

	// top volume

	Calo = gGeoManager->MakePcon("Calo", CaloMedium, 0.0, 360, 6);
	TGeoPcon   *pcon = (TGeoPcon*)(Calo->GetShape());
	pcon->DefineSection(0,-6110.0,2280,3886.352);
	pcon->DefineSection(1,-3400  ,2280,3886.352);
	pcon->DefineSection(2,-3400  ,1500,3886.352);
	pcon->DefineSection(3,+3400  ,1500,3886.352);
	pcon->DefineSection(4,+3400  ,2280,3886.352);
	pcon->DefineSection(5,+6110.0,2280,3886.352);
	gGeoManager->SetNsegments(64);

	Calo->AddNode(LArg, 1);
	Calo->AddNode(Tile, 1);

	gGeoManager->SetTopVolume(Calo);
	gGeoManager->CloseGeometry();
	gGeoManager->SetTopVisible();
   gGeoManager->SetVisLevel(5);

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