//////////////////////////////////////////////////////////////////////////
//
// Component: Geometry.cc
// Purpose: Implementation of Calorimeter Tower Geometry class.
//
//////////////////////////////////////////////////////////////////////////
#include <iostream>
#include <cstdio>
#include <math.h>
#include "Geometry.hh"

// Constructor
Geometry::Geometry(){  

  // Implement Tower Geometry
  for(int i = 0; i < (TOWER_NETA/2) ; ++i){
    
    _thetaMin[TOWER_NETA-1-i] = TOWER_THETA[i]*DEGRAD;
    _thetaMin[i] = (180-TOWER_THETA[i+1])*DEGRAD;
    
    _thetaMax[TOWER_NETA-1-i] = TOWER_THETA[i+1]*DEGRAD;
    _thetaMax[i] = (180-TOWER_THETA[i])*DEGRAD;
    
    _theta[i] = (_thetaMax[i]+_thetaMin[i])/2;
    _theta[TOWER_NETA-1-i] = ((_thetaMax[TOWER_NETA-1-i]+
			       _thetaMin[TOWER_NETA-1-i])/2);
    
    _eta[i] = -log(tan(_theta[i]/2.));
    _eta[TOWER_NETA-1-i] = -log(tan(_theta[TOWER_NETA-1-i]/2.));
    
    _etaMin[i] = -log(tan(_thetaMax[i]/2.));
    _etaMin[TOWER_NETA-1-i] = -log(tan(_thetaMax[TOWER_NETA-1-i]/2.));
    
    _etaMax[i] = -log(tan(_thetaMin[i]/2.));
    _etaMax[TOWER_NETA-1-i] = -log(tan(_thetaMin[TOWER_NETA-1-i]/2.));
    
    _cotTheta[i] = cos(_theta[i])/sin(_theta[i]);
    _cotTheta[TOWER_NETA-1-i] = cos(_theta[TOWER_NETA-1-i])/
      sin(_theta[TOWER_NETA-1-i]);
    
    _cotThetaMin[i] = cos(_thetaMax[i])/sin(_thetaMax[i]);
    _cotThetaMin[TOWER_NETA-1-i] = cos(_thetaMax[TOWER_NETA-1-i])/
      sin(_thetaMax[TOWER_NETA-1-i]);
    
    _cotThetaMax[i] = cos(_thetaMin[i])/sin(_thetaMin[i]);
    _cotThetaMax[TOWER_NETA-1-i] = cos(_thetaMin[TOWER_NETA-1-i])/
      sin(_thetaMin[TOWER_NETA-1-i]);
    
    _etaGranularity[i] = abs(_etaMax[i]-_etaMin[i]);
    _etaGranularity[TOWER_NETA-1-i] = abs(_etaMax[TOWER_NETA-1-i]
					  -_etaMin[TOWER_NETA-1-i]);
    
    _phiGranularity[i] = 2.0*PI*(1.0/(TOWER_PHI_SEG[TOWER_TYPE[i]]));
    
    _phiGranularity[TOWER_NETA-1-i] = 2.0*PI*(1.0/(TOWER_PHI_SEG[TOWER_TYPE[i]])); 
  }
} // Geometry Constructor

// Geometry Destructor
Geometry::~Geometry(){
}

// Tower Phi from ieta and iphi indices
float Geometry::phi(int IETA, int IPHI) const { 
       return 2*PI*((IPHI+0.5)/TOWER_PHI_SEG[TOWER_TYPE[IETA]]);
}

// Tower Phimin from ieta and iphi indices
float Geometry::phiMin(int IETA, int IPHI) const { 
       return 2*PI*((1.0*IPHI)/(1.0*TOWER_PHI_SEG[TOWER_TYPE[IETA]]));
}
								
// Tower Phimax from ieta and iphi indices
float Geometry::phiMax(int IETA, int IPHI) const { 
       return 2*PI*((1.0*(IPHI+1))/(1.0*TOWER_PHI_SEG[TOWER_TYPE[IETA]]));
} 
		
// tower print spew
void Geometry::spew() { 
     for(int i = 0; i < TOWER_NETA ; ++i){
     std::cout << std::endl;
     std::cout << "Tower eta index = " << i << std::endl;
     std::cout << "Tower theta = " << _theta[i] << std::endl;
     std::cout << "Tower thetaMin = " << _thetaMin[i] << std::endl;
     std::cout << "Tower thetaMax = " << _thetaMax[i] << std::endl;
     std::cout << "Tower eta = " << _eta[i] << std::endl;
     std::cout << "Tower etaMin = " << _etaMin[i] << std::endl;
     std::cout << "Tower etaMax = " << _etaMax[i] << std::endl;
     std::cout << "Tower cotTheta = " << _cotTheta[i] << std::endl;
     std::cout << "Tower cotThetaMin = " << _cotThetaMin[i] << std::endl;
     std::cout << "Tower cotThetaMax = " << _cotThetaMax[i] << std::endl;
     std::cout << "Tower eta Granularity = " << _etaGranularity[i] << std::endl;
     std::cout << "Tower phi Granularity = " << _phiGranularity[i] << std::endl;
	  }
} 

// get ieta index from x,y,z position in tower
size_t Geometry::iEta(float x, float y, float z) const{

  float R     = sqrt(x*x + y*y);
  float cotTheta;
  if (R > 0.0 ){
	cotTheta = z/R;
  } else {
	return 0;
  }
  
  if(cotTheta <=0){
    for(size_t i =0; i< (TOWER_NETA/2); ++i){ 
      if(cotTheta< _cotThetaMax[i]) return i;
    }
  }
  else{
    cotTheta = cotTheta*-1.0;
    for(size_t i =0; i< (TOWER_NETA/2); ++i){
      if(cotTheta<_cotThetaMax[i]) return (TOWER_NETA - 1- i);
    }
  }

  std::cerr << "Geometry: cotTheta out of bounds:  " << cotTheta << std::endl;
  std::cout << "Geometry: cotTheta out of bounds:  " << cotTheta << std::endl;
  std::cerr << "Geometry: x, y, z =  " << x << ",  " << y << ",  " << z << std::endl;
  std::cout << "Geometry: x, y, z =  " << x << ",  " << y << ",  " << z << std::endl;
  return 0;
}

// get ieta index from the real eta
size_t Geometry::iEta(float eta) const{

  if(eta <=0){
    for(size_t i =0; i< (TOWER_NETA/2); ++i){ 
      if(eta < _etaMax[i]) return i;
    }
  }
  else{
    eta = eta*-1.0;
    for(size_t i =0; i< (TOWER_NETA/2); ++i){
      if(eta<_etaMax[i]) return (TOWER_NETA - 1- i);
    }
  }

  std::cerr << "Geomtery: eta out of bounds:  " << eta << std::endl;
  return 0;
}

// get phi index from eta index and real phi
size_t Geometry::iPhi(size_t iEta, float phi) const{

  if (phi>2*PI) { phi -= 2.*PI*int ( phi/(2.*PI)); }
  if (phi<0)    { phi += 2.*PI*int((-phi/(2.*PI))+1); }

  size_t iPhi = int(phi/(2*PI) * TEPSEG[TETTYP[iEta]]);  

 return iPhi;
}