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


//-----------------------
// This Class's Header --
//-----------------------
#include "PndEmcStructure.h"
#include "PndEmcMapper.h"		
#include "PndEmcTwoCoordIndex.h"
#include "PndEmcXtal.h"

#include "TGeoMatrix.h"		
#include "TGeoManager.h"			
#include "TGeoArb8.h"
#include "TVector3.h"		
#include "TRotation.h"
#include "TMath.h"		
#include "TPRegexp.h"
#include "TString.h"
#include "TObjString.h"

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

using std::cout;
using std::endl;

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

PndEmcStructure::~PndEmcStructure(){}
		
PndEmcStructure* PndEmcStructure::Instance (TGeoManager *geoMan)
{
	if (_instance == 0) {
		if (geoMan==0){
			cout<<"Empty geometry passed to PndEmcStructure"<<endl;
			abort(); 
      }
		else
			_instance = new PndEmcStructure(geoMan);
  }
  return _instance;
}

// If PndEmcStructure is called without input parameters
// it means that it should have been instantiated with TGeoManager before
PndEmcStructure* PndEmcStructure::Instance ()
{
	if (_instance == 0) {
		cout<<"Wrong call for PndEmcStructure. It should be instantiated first with TGeoManager"<<endl;
		abort();
	}
	else
		return _instance;
}
		
PndEmcStructure::PndEmcStructure(TGeoManager *geoMan)
{
  
	// Instantiate mapper to convert from detId to TCI 
	// PndEmcMapper should be instatiated before first call here with proper input index
	PndEmcMapper *fEmcMap=PndEmcMapper::Instance(0);
  
	int module,copy,row,crystal;
	TString node_path;
	const char* crystal_name;
	const TGeoMatrix *crystal_matrix;
	const double *trans;
	TGeoRotation geoRot;
	TRotation rot;
	Double_t phi, theta, psi;
	int detId;
	
	TGeoIterator next(geoMan->GetTopVolume());
   TGeoNode *node;
   while ((node=next())) {
		next.GetPath(node_path);
		
		bool isEmcModule=crystal_name_analysis(node_path,module,copy,row,crystal);
		if (!isEmcModule) continue;
		
		detId =  module*100000000 + row*1000000 + copy*10000 + crystal; 
		PndEmcTwoCoordIndex *tci=fEmcMap->GetTCI(detId);
		
		// Obtaine TGeoMatrix matrix for the current node, later it is factorised to translation aand rotation
		crystal_matrix=next.GetCurrentMatrix();
		
		trans=crystal_matrix->GetTranslation();
		TVector3 pos(trans[0],trans[1],trans[2]);
		emcX[detId] = pos.X(); 
		emcY[detId] = pos.Y();
		emcZ[detId] = pos.Z();
		
		geoRot.SetMatrix(crystal_matrix->GetRotationMatrix());
		
		PndEmcXtal *xtal;
		if (module==5) {
			TGeoBBox *box = (TGeoBBox *) node->GetVolume()->GetShape();
			
			TGeoTrap crystal(box->GetDZ(), 0, 0,
			box->GetDY(), box->GetDX(), box->GetDX(), 0,
			box->GetDY(), box->GetDX(), box->GetDX(), 0);

			xtal=new PndEmcXtal(tci,crystal,pos,geoRot);
		}
		else{
			
			TGeoTrap *trap1 = (TGeoTrap*) node->GetVolume()->GetShape();
			
			TGeoTrap crystal(trap1->GetDz(), trap1->GetTheta(), trap1->GetPhi(),
			trap1->GetH1(), trap1->GetBl1(), trap1->GetTl1(), trap1->GetAlpha1(),
			trap1->GetH2(), trap1->GetBl2(), trap1->GetTl2(), trap1->GetAlpha2());
			xtal=new PndEmcXtal(tci,crystal,pos,geoRot);
		}
		
		fTciXtalMap[tci]=xtal;
   }
}

bool PndEmcStructure::crystal_name_analysis(TString node_path,int &module,int &copy,int &row,int &crystal)
{
	// The following code extract information about module, copy, row and crystal number from the path name of the TGeoNode which corresponds to emc crystal
	// Name convention is according to PndEmc.cxx
	if (!(node_path.Contains("emc"))) return false;
	
	// Case of Fsc
	if (node_path.Contains("Fsc")){
		//at the moment all the layers of module in Fsc are not taken into account, only the whole block is selected
		if (node_path.Contains("FscLayer"))
			return false;
		TObjArray *subStrL =
		TPRegexp("^cave/Fsc_(\\d+)/emc(\\d+)r(\\d+)c(\\d+)_0$").MatchS(node_path);
		
		copy  = (((TObjString *)subStrL->At(1))->GetString()).Atoi();
		module = (((TObjString *)subStrL->At(2))->GetString()).Atoi();
		row  = (((TObjString *)subStrL->At(3))->GetString()).Atoi();
		crystal  = (((TObjString *)subStrL->At(4))->GetString()).Atoi();
		return true;
	}
	
	// Case of barrel
	if (node_path.Contains("EmcLayer")){
		TObjArray *subStrL = TPRegexp("^cave/Emc\\d+_(\\d+)/EmcLayer\\d_0/emc(\\d+)r(\\d+)c(\\d+)_0$").MatchS(node_path);
		if(subStrL->GetLast()<4){
			cout<<"crystal name in barrel Emc: "<<node_path<<" missmatch pattern"<<endl;
			return false;
		}
		copy  = (((TObjString *)subStrL->At(1))->GetString()).Atoi();
		module = (((TObjString *)subStrL->At(2))->GetString()).Atoi();
		row  = (((TObjString *)subStrL->At(3))->GetString()).Atoi();
		crystal  = (((TObjString *)subStrL->At(4))->GetString()).Atoi();
	}
	//case of endcups of forward calorimeter
	else {
		TObjArray *subStrL = TPRegexp("^cave/Emc\\d_(\\d+)/emc(\\d+)r(\\d+)c(\\d+)_0$").MatchS(node_path);
		if(subStrL->GetLast()<4){
			cout<<"crystal name in endcup Emc: "<<node_path<<" missmatch pattern"<<endl;
			return false;
		}
		copy  = (((TObjString *)subStrL->At(1))->GetString()).Atoi();
		module = (((TObjString *)subStrL->At(2))->GetString()).Atoi();
		row  = (((TObjString *)subStrL->At(3))->GetString()).Atoi();
		crystal  = (((TObjString *)subStrL->At(4))->GetString()).Atoi();
	}
	
	return true;
}


PndEmcTwoCoordIndex* 
PndEmcStructure::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);
	PndEmcTwoCoordIndex *tci;
	double diff=1000;
	     
	std::map <PndEmcTwoCoordIndex*, PndEmcXtal*>::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;
}


// Print information on angular position for all the crystall into the file
void PndEmcStructure::Print(string filename) const
{
	ofstream f(filename.c_str());
	f<<"detID"<<"\t"<<"ThetaInd"<<"\t"<<"PhiInd"<<"\t"<<"theta_centre"<<"\t"<<"theta_frontface"<<"\t"<<"dTheta"<<"\t"<<"phi_centre"<<"\t"<<"phi_frontface"<<"\t"<<"dPhi"<<endl;
	TVector3 centre, front_centre;
	double theta_c, theta_f, phi_c, phi_f;
	double dTheta, dPhi;
	std::map <PndEmcTwoCoordIndex*, PndEmcXtal*>::const_iterator iter=fTciXtalMap.begin();
	PndEmcTwoCoordIndex *tci;
	while(iter!=(fTciXtalMap).end())
	{
		tci = (*iter).first;
		centre= ((*iter).second)->centre();
		front_centre=((*iter).second)->frontCentre();
		theta_c=centre.Theta()*TMath::RadToDeg();
		phi_c=centre.Phi()*TMath::RadToDeg();
		theta_f=front_centre.Theta()*TMath::RadToDeg();
		phi_f=front_centre.Phi()*TMath::RadToDeg();

		dTheta=theta_c-theta_f;
		dPhi=phi_c-phi_f;

		f<<tci->itsIndex()<<"\t"<<tci->itsXCoord()<<"\t"<<tci->itsYCoord()<<"\t"<<theta_c<<"\t"<<theta_f<<"\t"<<dTheta<<"\t"<<phi_c<<"\t"<<phi_f<<"\t"<<dPhi<<endl;
		iter++;
	}
       
	f.close();
	cout<<"write emc structure to log file"<<endl;

}

ClassImp(PndEmcStructure)