//======================================================================
// File and Version Information:
// 	$Id: $
//
// Description:
//       Class EmcStructure
//
//======================================================================


//-----------------------
// This Class's Header --
//-----------------------
#include "EmcStructure.h"
#include "EmcMapper.h"		
#include "CbmEmcReader.h"

#include "TwoCoordIndex.h"
#include "EmcXtal.h"
#include "TGeoArb8.h"
#include "TVector3.h"		
#include "TRotation.h"
#include "TMath.h"		

//---------------
// C++ Headers --
//---------------
#include "stdlib.h"
#include <iostream>
#include <string>
#include <cassert>		


//-------------------------------
// Collaborating Class Headers --
//-------------------------------
EmcStructure* EmcStructure::_instance = 0;

EmcStructure::~EmcStructure(){}
		
// It is aasumed that first call of EmcStructure provide
// geometry file name, but for subsequent calls filename can be empty
EmcStructure* EmcStructure::Instance (string fileGeo)
{
  if (_instance == 0) {
    if (fileGeo.empty())
      {
	cout<<"Geometry file is empty, programm aborted"<<endl;
	abort(); 
      }
    else
      _instance = new EmcStructure(fileGeo);
  }
  return _instance;
}
		
EmcStructure::EmcStructure(string fileGeo)
{
  TString work = getenv("VMCWORKDIR");
  work = work + "/geometry/" + fileGeo;
  cout << "-I- <EmcStructure> Emc geometry loaded from: " << work << endl;
  Int_t detId = -1;
  CbmEmcReader read(work);
  
  // Instantiate mapper to convert from detId to TCI 
  // EmcMapper should be instatiated before first call here with proper input index
  EmcMapper *fEmcMap=EmcMapper::Instance(0);
  
  for(Int_t module=1; module<=read.GetMaxModules(); module++) 
    for(Int_t row=1; row<=read.GetMaxRows(module); row++)
      for(Int_t crystal=1; crystal<=read.GetMaxCrystals(module,row); crystal++) {
	DataG4 data = read.GetData(module,row,crystal);
	
	 Int_t copyNo;
	 if (module==3 || module==4 || module==5) 
	   copyNo=4; 
	 else
	   copyNo=16;
	
	 
	 for(Int_t copy=1; copy<=copyNo; copy++)
	   { 
	     TVector3 ref(1., 1., 1.); // reflection matrix
	     detId =  module*100000000 + row*1000000 + copy*10000 + crystal; 
	     
	     TRotation rot1;
	     emcTheta[detId] = data.theta; 
	     emcTau[detId] = data.tau;
	     
	     if ((module==1) || (module==2))
	       {
		 emcPhi[detId] = fmod(data.phi+22.5*(copy-1),360);
		 rot1.RotateZ(TMath::Pi()*22.5/180.*(copy-1));
	       }
	     
	     TVector3 pos(0., 0., 0.);
	     if ((module==1) || (module==2)) 
	       {
		 pos.SetXYZ(data.posX,data.posY,data.posZ+37.); // shift of 37 mm for the barrel EMC
	       }
	     else
	       {
		 if (module==5)
		   pos.SetXYZ(data.posX,data.posY,data.posZ+0.5*300*(2.*data.pDz+data.pAlp1+data.pAlp2)); //added half of the pad size
		 else
		   pos.SetXYZ(data.posX,data.posY,data.posZ);
		 if (copy==2) ref.SetXYZ( 1., -1, 1.);
		 if (copy==3) ref.SetXYZ(-1., -1, 1.);
		 if (copy==4) ref.SetXYZ(-1.,  1, 1.);
	       }
	     
	     pos*=rot1;
	     emcX[detId] = pos.X() * ref.X(); 
	     emcY[detId] = pos.Y() * ref.Y();
	     emcZ[detId] = pos.Z() * ref.Z();
	     
	     // Fill map between Two Coordinate Index (TCI) of crystal and EmcXtal object 
	     //which containe information about crystal geometry and 
	     TwoCoordIndex *tci=fEmcMap->GetTCI(detId);
	     
	     Double_t trapParams[11];
	     
	     if (module==5) 
	       {
		 trapParams[0] = 300.*(2.*data.pDz+data.pAlp1+data.pAlp2)/10.;
		 trapParams[1] = 0.;
		 trapParams[2] = 0.;
		 trapParams[3] = TMath::Max(data.pDy1, data.pDy2)/10.;
		 trapParams[4] = TMath::Max(data.pDx1, data.pDx2)/10.;
		 trapParams[5] = TMath::Max(data.pDx1, data.pDx2)/10.;
		 trapParams[6] = 0.;
		 trapParams[7] = TMath::Max(data.pDy1, data.pDy2)/10.;
		 trapParams[8] = TMath::Max(data.pDx1, data.pDx2)/10.;
		 trapParams[9] = TMath::Max(data.pDx1, data.pDx2)/10.;
		 trapParams[10] = 0.;
	       }
	     else
	       {
		 trapParams[0] = data.pDz/10.;
		 trapParams[1] = data.pTheta;
		 trapParams[2] = data.pPhi;
		 trapParams[3] = data.pDy1/10.;
		 trapParams[4] = data.pDx1/10.;
		 trapParams[5] = data.pDx2/10.;
		 trapParams[6] = data.pAlp1;
		 trapParams[7] = data.pDy2/10.;
		 trapParams[8] = data.pDx3/10.;
		 trapParams[9] = data.pDx4/10.;
		 trapParams[10] = data.pAlp2;
	       }
	     
	     TGeoTrap trap;
	     trap.SetDimensions(trapParams);
	     
	     TRotation rot;
	     rot.RotateZ(TMath::Pi()*data.tau/180);
	     rot.RotateY(TMath::Pi()*data.theta/180);
	     rot.RotateZ(TMath::Pi()*data.phi/180);
	     rot.RotateZ(TMath::Pi()*22.5*(copy-1)/180);
	     
	     EmcXtal *xtal=new EmcXtal(tci,trap,pos,rot, ref);
	     
	     fTciXtalMap[tci]=xtal;
	     
	   }
      }
}


TwoCoordIndex* 
EmcStructure::locateIndex( double theta, double phi ) const
{
	// Find crystal with minimal angular difference with given direction
	TVector3 vec(0,0,10);
	vec.SetTheta(theta);
	vec.SetPhi(phi);
	TwoCoordIndex *tci;
	double diff=1000;
     
	std::map <TwoCoordIndex*, EmcXtal*>::const_iterator iter=fTciXtalMap.begin();
	while(iter!=(fTciXtalMap).end())
	{
		TVector3 vec2= ((*iter).second)->frontCentre();
		double tmpDiff=vec2.Angle(vec);

		if(tmpDiff<diff)
		{
			tci=(*iter).first;
			diff=tmpDiff;
		}
		
		iter++;
	}

	return tci;
}

ClassImp(EmcStructure)