//----------------------------------------------------------------------------
//
//    File     : CellCorners.cpp
//    Purpose  : Calculate the corners of a cell presumably for graphics
//               purposes.  A coarse and fine description are avaialable.
//    Created  : 24-FEB-1999  by Stephen Kahn
//    
//----------------------------------------------------------------------------

//  include files

#include <list>
#include "calorimeter_geometry/utils/CellCorners.hpp"
#include "geometry_system/components/GeometryElement.hpp"
#include "geometry_system/components/CartesianCoordinate.hpp"
#include "caladdress/CALparams.h"
#include "PhysicsVectors/Rotation.h"

//  namespaces

using  namespace std;
using  namespace dgs;

//  constructors

CellCorners::CellCorners( const CellAddress& ca) : _ca(ca)  
{  fill_map(); 
   fill_edges();
}

//  accessors

//  methods

bool  CellCorners::fill_map()
{   CALcell cc( _ca.ieta(), _ca.iphi(), _ca.layer());
    if( !cc.valid_cell()) return false;

    Ctype ctype = cc.cell_type();

    return cell_corners(); 
}

bool  CellCorners::cell_corners()
{   CellShape  cs( _ca, -1);
    double  xl, yl, zl;
    float  drmid, rmid, h1, dz, r, rmin, rmax, zz;
    double dphi;
    short  isgn[2] = { -1, 1};
    InOut   inout[2] = { INNER, OUTER };
    InOut   frtbck[4] = { INNER_FRONT, OUTER_FRONT, INNER_BACK, OUTER_BACK};
    int  i, j, k, ikey, nzpts;
    GeometryXform xform( cs.position(), cs.orientation());

    switch (cs.shape_type())
      {  case SLTUBS :
	 { dphi = cs.half_ang_width( INNER);
	   drmid = 0.5*(cs.radius(INNER)+cs.radius(OUTER))*cos(dphi);
	   rmid = 0.5*(cs.radius(INNER)+cs.radius(OUTER));
	   for( i=0; i<2; i++)  // inner-outer loop
	    {  InOut io = inout[i];
	       for( j=0; j<2; j++)  // lower-higher azimuth
		 {  yl = isgn[j]*cs.radius(io)*sin(dphi);
		    xl = cs.radius( io) - drmid;
		    for( k=0; k<2; k++)
		      {  zl = cs.cell_local_z(xl) + isgn[k]*cs.half_length();
		         CartesianCoordinate lpos( xl, yl, zl);
			 CartesianCoordinate gpos = xform.apply(lpos);
			 SpaceVector pos( gpos.x(), gpos.y(), gpos.z());
			 ikey = corner_key( i, j, k);
			 add_corner( ikey, pos);
		      }
		 }
	    }
	   break;
	 }

         case TUBS :
	 { dphi = cs.half_ang_width( INNER);
	   drmid = 0.5*(cs.radius(INNER)+cs.radius(OUTER))*cos(dphi);
	   rmid = 0.5*(cs.radius(INNER)+cs.radius(OUTER));
	   for( i=0; i<2; i++)  // inner-outer loop
	    {  InOut io = inout[i];
	       for( j=0; j<2; j++)  // lower-higher azimuth
		 {  yl = isgn[j]*cs.radius(io)*sin(dphi);
		    xl = cs.radius( io) - drmid;
		    for( k=0; k<2; k++)
		      {  zl = isgn[k]*cs.half_length();
		         CartesianCoordinate lpos( xl, yl, zl);
			 CartesianCoordinate gpos = xform.apply(lpos);
			 SpaceVector pos( gpos.x(), gpos.y(), gpos.z());
			 ikey = corner_key( i, j, k);
			 add_corner( ikey, pos);
		      }
		 }
	    }
	   break;
	 }

         case CONS :
	 { dphi = cs.half_ang_width( INNER);
	   drmid = 0.25*(cs.radius(INNER_FRONT) + cs.radius(OUTER_FRONT)
		+  cs.radius(INNER_BACK) + cs.radius(OUTER_BACK))*cos(dphi);
	   rmid =  0.25*(cs.radius(INNER_FRONT) + cs.radius(OUTER_FRONT)
		+  cs.radius(INNER_BACK) + cs.radius(OUTER_BACK));
	   for( k=0; k<2; k++)  // front-back loop
	    {  zl = isgn[k]*cs.half_length();
	       for( i=0; i<2; i++)  //  inner-outer loop
	         {  InOut io = frtbck[2*k+i];
        	    for( j=0; j<2; j++)  // lower-higher azimuth
  			 {  yl = isgn[j]*cs.radius(io)*sin(dphi);
			    xl = cs.radius( io) - drmid;
		            CartesianCoordinate lpos( xl, yl, zl);
			    CartesianCoordinate gpos = xform.apply(lpos);
			    SpaceVector pos( gpos.x(), gpos.y(), gpos.z());
			    ikey = corner_key( i, j, k);
			    add_corner( ikey, pos);
		      }
		 }
	    }
	   break;
	 }

         case TRAP1 :
	 { h1 = cs.half_height();
	   dz = cs.half_length();
	   for( i=0; i<2; i++)
	     {  xl = isgn[i]*h1;
	        for( j=0; j<2; j++)
		  {  yl = cs.cell_local_x( xl) - cs.cell_local_x(0.);
		     yl += isgn[j]*cs.cell_dx( xl);
		     for( k=0; k<2; k++)
		       {  zl = isgn[k]*cs.half_length();
		          CartesianCoordinate lpos( xl, yl, zl);
			  CartesianCoordinate gpos = xform.apply(lpos);
			  SpaceVector pos( gpos.x(), gpos.y(), gpos.z());
			  int ikey = corner_key( i, j, k);
			  add_corner( ikey, pos);
		       }
		  }
	     }
	   break;
	 }

         case TRAP2 :
	 { h1 = cs.half_height();
	   for( i=0; i<2; i++)
	     {  xl = isgn[i]*h1;
	        for( j=0; j<2; j++)
		  {  yl = cs.cell_local_x( xl) - cs.cell_local_x(0.);
		     yl += isgn[j]*cs.cell_dx( xl);
		     for( k=0; k<2; k++)
		       {  zl = cs.cell_local_z( xl) + isgn[k]*cs.cell_dz(xl);
		          CartesianCoordinate lpos( xl, yl, zl);
			  CartesianCoordinate gpos = xform.apply(lpos);
			  SpaceVector pos( gpos.x(), gpos.y(), gpos.z());
			  ikey = corner_key( i, j, k);
   			 add_corner( ikey, pos);
		       }
		  }
	     }
	   break;
	 }

	case  PCON :
	{   nzpts = static_cast<int> (cs.param(2));
	    const float tol = 4.0;
	    int ieta = _ca.ieta();
	    double cell_azimuth = atan2( cs.position().y(), cs.position().x());
	    for ( j=0; j<2; j++)    // azimuthal loop
		{  float phi = cs.center_angle() + isgn[j]*cs.half_ang_width();
		   phi += cell_azimuth - cs.center_angle(); //shortcut to 
		                                            // avoid xform 
		   for( k=0; k<nzpts; k++)   //  z positions loop
		    {  zz = cs.param( 3*k + 3);
		       rmin = cs.param( 3*k + 4);
		       rmax = cs.param( 3*k + 5);
		       int kk = k;
		       if( ieta < 0) kk = nzpts - k -1;
			if( nzpts == 3 )  // triangular (prizm) shapes
			  {  r = rmin;
     			     if( k == 1 && ((cs.param(4) == cs.param(7) &&
				 ieta>0) || (cs.param(10) == cs.param(7) &&
				 ieta<0)))   r = rmax;
			     xl = r*cos(phi);
			     yl = r*sin(phi);
			     SpaceVector lpos( xl, yl, zz); // global position
			     ikey = corner_key( 0, j, kk, 1);
			     add_corner( ikey, lpos);
			  }
			else if ( nzpts == 4) // trapezoidal
			  {  r = rmin;
			     if( (ieta >  0 && k == 2)||(ieta < 0 && k == 1))
				 r = rmax;
			     xl = r*cos(phi);
			     yl = r*sin(phi);
			     SpaceVector lpos( xl, yl, zz); // global position
			     ikey = corner_key( 0, j, kk, 2);
			     add_corner( ikey, lpos);
			  }

			else if ( nzpts == 5 && (cs.param(3) == cs.param(6) 
			      && ieta>0) || (cs.param(15) == cs.param(12) &&
			      ieta<0))  // ICH or MCH pentagon 
			  {  r = rmax;
			     if((ieta>0 && (k == 0 || k==3)) || 
				(ieta<0 && (k == 4 || k==1))) r = rmin;
			     xl = r*cos(phi);
			     yl = r*sin(phi);
			     SpaceVector lpos( xl, yl, zz);
			     ikey = corner_key( 0, j, kk, 1);
			     add_corner( ikey, lpos);
			  }

			else if (nzpts == 5 && ((ieta>0 && fabs(cs.param(11) -
			     cs.param(14)) > tol) || ( ieta<0 && 
			     fabs(cs.param(11) - cs.param(8)) > tol)))
			  {  r = rmin;   // ECOH lower pentagon
			     if( (k == 3 && ieta>0) || (k == 1 && ieta<0 ))
			       r = rmax;
			     xl = r*cos(phi);
			     yl = r*sin(phi);
			     SpaceVector lpos( xl, yl, zz);
			     ikey = corner_key( 0, j, kk, 1);
			     add_corner( ikey, lpos);
			  }

			else if (nzpts == 5) // ECOH upper pentagon
			  {  r = rmax;   
			     if( (ieta>0 && k == 1) || (ieta<0 && k == 3) )
			       r = rmin;
			     xl = r*cos(phi);
			     yl = r*sin(phi);
			     SpaceVector lpos( xl, yl, zz);
			     ikey = corner_key( 0, j, kk, 3);
			     add_corner( ikey, lpos);
			  }
		    }
		}
	}	
      } 
    return true;
}

ostream&  operator <<(ostream& os, CellCorners& me)
{       char spc[4] = "   ";
	os << "Corners for calorimeter cell " << me._ca << endl;
	
	os << (Corners&) me;
/*  skip
	for (int i = 0; i < 2; i++)
	for (int j = 0; j < 2; j++)
	for (int k = 0; k < 2; k++)
	 {  
	    int ikey  = me.corner_key( i, j, k);
	    SpaceVector pos = me.corner(ikey);
	    os << "Corner " << ikey << spc << pos.x() << spc << pos.y()
		<< spc << pos.z() << endl;
	 }
end skip */  

	return os;
}

int  CellCorners::corner_key( int i, int j, int k, int it)
  {   if( it == 0 ) return Corners::corner_key( i, j, k);
  
      if( it == 1)  // triangle (prizm) forms or ICH, MCH, ECOH lower
	            // pentagons from PCON
	{  int ikey = j + 2*k;
	   return ikey;
	}

      else if (it == 2) // trapezoidal forms from PCON
	{  return Corners::corner_key( k%2, j, k/2);
	}

      else if (it == 3) // ECOH upper pentagons
	{  int kk = k;
	   if ( k == 1 ) kk = 3;
	   else if ( k == 2 ) kk = 1;
	   else if ( k == 3 ) kk = 2;
	   return 2*kk + j;
	}

      return 0;
  }

bool  CellCorners::fill_edges()
{   int ncorn = get_corners().size();
    int i;

    if( ncorn == 8)    // quadralaterial
      {  int ipair[12][2] = { {0,1}, {0,2}, {0,4}, {1,3}, {1,5}, {2,3},
			      {2,6}, {3,7}, {4,5}, {4,6}, {5,7}, {6,7}};
         for ( i=0; i<12; i++)
	   add_edge( ipair[i][0], ipair[i][1]);
      }
	      
    else if ( ncorn == 6)  //  triangular
      {  int ipair[9][2] = { {0,1}, {0,2}, {0,4}, {1,3}, {1,5}, {2,3},
			     {2,4}, {3,5}, {4,5}};
         for ( i=0; i<9; i++)
	   add_edge( ipair[i][0], ipair[i][1]);
      }

    else if ( ncorn == 10) //  pentagon
      {  int ipair[15][2] = { {0,1}, {0,2}, {0,6}, {1,3}, {1,7}, {2,4},
			      {2,3}, {3,5}, {4,8}, {4,5}, {6,8}, {6,7},
			      {8,9}, {5,9}, {7,9}};
         for ( i=0; i<15; i++)
	   add_edge( ipair[i][0], ipair[i][1]);
      }

    return true;
}