#include "include.h"

int toro(){
//int main(){

	gStyle->SetCanvasPreferGL(true);
	gSystem->Load("libGeom");
	gSystem->Load("libPhysics.so");
	
	/*if (gGeoManager)
		delete gGeoManager;*/
	TGeoManager *geom = new TGeoManager("cylindric", "Cylindric Light Guide");
	
	//-- Define Shape Variables
	Bool_t DEBUG = kFALSE;	// Debug mode to show message
	Bool_t DrawRay = kTRUE;	// Flag to show rays or not
	Bool_t ROUGH = kTRUE;	// Add roughness to the shape's surface
	Int_t Nray = 100;	// Number of rays generated
	Double_t toroR = 5;	// Toroid radius
	Double_t tubR  = 1;	// External raius of the tub formed by toroid
	Double_t toroA = 90;	// Final angle of the toroid
	Double_t apertureA = 20;// Aperture angle for light source; from 0 to 89.99
	Double_t TIR_A = 1/1.49;// Snell Law
	Double_t roughMax = 0.01;// Max deviation by roughness
	Double_t roughApertureA = 10; // Aperture of the roughtness
	Double_t degrad = TMath::Pi()/180;
	
	//--- Initializate the viewer
	TCanvas *myCan = new TCanvas("myCan","Toroidal Light Guide Test", 800, 600);

	//--- define some materials
	TGeoMaterial *matVacuum = new TGeoMaterial("Vacuum", 0,0,0);
	TGeoMaterial *matPoly = new TGeoMaterial("Poly",0,0,0);
	
	//--- define some media
	TGeoMedium *Vacuum = new TGeoMedium("Vacuum",1, matVacuum);
	TGeoMedium *Plexi = new TGeoMedium ("Plexiglass",3,matPoly);

	//--- make the top container volume
	TGeoVolume *top = geom->MakeBox("TOP", Vacuum, 250, 250, 250);
	geom->SetTopVolume(top);
	top->SetVisibility(kTRUE);

	//--- Define some Shapes
	TGeoTorus *tor = new TGeoTorus("tor",toroR,0,tubR,0,toroA);// The toroidal part
	TGeoVolume *cub = geom->MakeBox("cub", Plexi, tubR/TMath::Sqrt(2), tubR/TMath::Sqrt(2), tubR/TMath::Sqrt(2));
	TGeoTrap *trapShape = new TGeoTrap("trap1",tubR,0,0,tubR/TMath::Sqrt(2),tubR/TMath::Sqrt(2),tubR/TMath::Sqrt(2),0,tubR,tubR,tubR,0);
	TGeoTube *cil1Shape = new TGeoTube("cil1",0,tubR,tubR);
	
	//--- Define some Transformation
	TGeoTranslation *t1 = new TGeoTranslation("t1",0,0,tubR/TMath::Sqrt(2));
	TGeoTranslation *t2 = new TGeoTranslation("t2",0,0 , tubR*(1+2/TMath::Sqrt(2)));
	TGeoTranslation *t3 = new TGeoTranslation("t3",-toroR,0 , 2*tubR*(1+3/TMath::Sqrt(2)));
	TGeoCombiTrans *c1 = new TGeoCombiTrans("c1",-toroR,0,2*tubR+tubR*TMath::Sqrt(2) -1E-15, new TGeoRotation("r1",0,90,0));
	t1->RegisterYourself();
	t2->RegisterYourself();
	c1->RegisterYourself();
	
	//--- Defining the whole volume
	TGeoCompositeShape *wholeShape = new TGeoCompositeShape("WS1","(cil1*trap1):t2 + cub:t1 + tor:c1");
	TGeoVolume *lGuide1 = new TGeoVolume("lguide1",wholeShape,Plexi);
	lGuide1->SetLineColor(kGreen);
	top->AddNode(lGuide1,1);
	
	geom->CloseGeometry();
        geom->SetVisLevel(4);
        top->Draw();

	//--- Positioning and Directing the ray
	Bool_t cont;
	Int_t dummy,i,j,trackId,outCounter;
	Double_t x,z,tet,costet,phy,TIR_B;
	TVector3 dir,n,roughness;
	TGeoNode *current;
	TGeoNode *otro;
	Double_t *temp;
	Double_t *pos = new Double_t[3];
	outCounter = 0;
	std::string path = "/TOP_1/lguide1_1";
	
	for(i=0;i<Nray;i++){
		pos[0] = (2*gRandom->Rndm(gRandom->Rndm(dummy)) - 1)*tubR/TMath::Sqrt(2);// X
		pos[1] = (2*gRandom->Rndm(gRandom->Rndm(dummy)) - 1)*tubR/TMath::Sqrt(2);// Y
		pos[2] = 0;// Z
		
		costet = gRandom->Uniform(TMath::Sin(degrad*(90-apertureA)),1);
		tet = TMath::ACos(costet);
		phy = 2*TMath::Pi()*gRandom->Rndm(dummy);
		dir.SetXYZ(TMath::Sin(tet)*TMath::Cos(phy),TMath::Sin(tet)*TMath::Sin(phy),costet);
		dir.SetMag(1);
		
		TGeoTrack *track = new TGeoTrack(i+1,i+1);
		cont = kTRUE;
		j=0;
		while(cont){
			j++;
			trackId = geom->AddTrack(j,1);
			geom->SetCurrentTrack(trackId);
			track->AddPoint(pos[0],pos[1],pos[2],0);
			geom->InitTrack(pos[0],pos[1],pos[2],dir.X(),dir.Y(),dir.Z());
			if(path.compare(gGeoManager->GetPath())!=0){
				if(DEBUG) std::cout << "Track is outside the shape" << std::endl;
				break;
			}
			otro = geom->FindNextBoundary();
			otro = geom->Step(kTRUE,kFALSE);
			temp = geom->FindNormal(kFALSE);
			if(temp == NULL){
				if(DEBUG) std::cout << "Normal Vector is NULL:" << std::endl;
				break;
			}
			n.SetXYZ(temp[0],temp[1],temp[2]);
			pos = (Double_t *) geom->GetCurrentPoint();
			if(pos == NULL){
				if(DEBUG) std::cout << "Current Position is NULL" << std::endl;
				break;
			}
			//--- Check if the normal correspond to the output face
			if(n.X() == -TMath::Sin(degrad*toroA) && n.Y()<1E-15 && n.Z()==TMath::Cos(degrad*toroA)){
				if(DEBUG) std::cout << "End face of the Geometry" << std::endl;
				outCounter++;
				cont = kFALSE;
				if(DrawRay){
					track->SetLineColor(kBlue);
				}
			}else{
				//--- Adding some roughness to the surface
				if(ROUGH){
				//--- The code for add roughness need some evaluation in order to
				//--- add the roughness in the direction of normal.
				//--- I don't know if this is accomplished in the actual way.
						costet = gRandom->Uniform(TMath::Sin(degrad*(90-roughApertureA)),1);
						tet = TMath::ACos(costet);
						phy = 2*TMath::Pi()*gRandom->Rndm(dummy);
						roughness.SetXYZ(TMath::Sin(tet)*TMath::Cos(phy),TMath::Sin(tet)*TMath::Sin(phy),costet);
						roughness.SetMag(roughMax);
						n += roughness;
				}
				TIR_B = TMath::Sin(TMath::ACos(n.Dot(dir)));
				if(TIR_B < TIR_A){
					if(DEBUG) std::cout << "Snell Law not accomplished" << std::endl;
					cont = kFALSE;
				}else{
					dir = (dir-2*n.Dot(dir)*n).Unit();
				}
			}
			if(DrawRay)
				track->Draw();
		}
		if(i%100==0)
			std::cout << "Iteración n°: " << i << std::endl;
	}
	std::cout << "Losses ratio: " << (Double_t) outCounter/Nray << std::endl;
	//myCan->SaveAs("torus.root");
	return 0;
}