#include "PndTrkCTGeometryCalculations.h"
#include "PndTrkMergeSort.h"
#include <iostream>
#include <cmath>

// Root includes
#include "TROOT.h"

#define PI		3.141592654

using namespace std;

//----------begin of function PndTrkCTGeometryCalculations::CalculateArcLength

Double_t PndTrkCTGeometryCalculations::CalculateArcLength(
	Double_t Oxx,
	Double_t Oyy,
	Double_t Rr,
	Short_t charge,
	Double_t *Xcross, // entrance-exit point
	Double_t *Ycross // entrance-exit point
	)
{
	Double_t	dis,
			theta1,
			theta2;

	theta1 = atan2(Ycross[0]-Oyy,  Xcross[0]- Oxx);
	theta2 = atan2(Ycross[1]-Oyy,  Xcross[1]- Oxx);


	if(charge>0){  // the rotation was clockwise.
		dis = theta1-theta2;
	} else {  // the rotation was counterclockwise.
		dis = theta2-theta1;
	}

	if(dis<0.) dis += 2.*PI;
	if(dis<0.) dis =0.;

	dis *= Rr;

	return dis;
}
//----------end of function PndTrkCTGeometryCalculations::CalculateArcLength




//----------begin of function PndTrkCTGeometryCalculations::CalculateCircleThru3Points

bool PndTrkCTGeometryCalculations::CalculateCircleThru3Points(
	Double_t x1,
	Double_t y1,
	Double_t x2,
	Double_t y2,
	Double_t x3,
	Double_t y3,
	Double_t *o_x,
	Double_t *o_y,
	Double_t *r_r
	)
{
	// inspiration from  http://mathworld.wolfram.com/Circle.html

	Double_t a,
		 d,
		 e,
		 q1,
		 q2,
		 q3;

	a = x1*y2 + y1*x3 + x2*y3 - x3*y2 - x2*y1 - x1*y3;

	if( fabs(a)<1.e-10 ) return false;

	q1 = x1*x1 + y1*y1;
	q2 = x2*x2 + y2*y2;
	q3 = x3*x3 + y3*y3;

	d = -(q1*y2 + y1*q3 + q2*y3 - q3*y2 - q2*y1 - q1*y3);

	e = -(x1*q2 + q1*x3 + x2*q3 - x3*q2 - x2*q1 - x1*q3);

	*o_x = -0.5*d/a;
	*o_y = -0.5*e/a;
	*r_r = sqrt( (x1-*o_x)*(x1-*o_x) + (y1-*o_y)*(y1-*o_y) );


	return true;
}
//----------end of function PndTrkCTGeometryCalculations::CalculateCircleThru3Points




//-------------------  begin of function   PndTrkCTGeometryCalculations::calculateintersections

void PndTrkCTGeometryCalculations::calculateintersections(
	Double_t Oxx,
	Double_t Oyy,
	Double_t Rr,
	Double_t C0x,
	Double_t C0y,
	Double_t C0z,
	Double_t r,
	Double_t vx,
	Double_t vy,
	Double_t vz,
	Int_t *STATUS,
	Double_t* POINTS)
{


//-------------------------------------------


  Double_t P1x, P1y, P1z, P2x, P2y, P2z;

  Double_t AAA, DELTA, ax, ay, aaa;




/*
 INPUTS :

  Oxx, Oyy        = abscissa and ordinate of the center of the circular trajectory of
                  the particle;
  Rr             = radius of such trajectory;
  C0x, C0y, Coz = x, y, z coordinates of a point belonging to the axis of the
                  skewed straw;
  r  = radius of equidrift of such skewed straw;
  vx, vy, vz    =  versor of the direction along which the skewed straw lies.


 OUTPUTS :

  P1x, P1y, P1z  =  x, y, z coordinates of the point intersection between the
                    particle trajectory circle and the equidrift cylinder of
                    the skewed straw calculated as a function of theta (first
                    solution);
  P2x, P2y, P2z  =  x, y, z coordinates of the point intersection between the
                    particle trajectory circle and the equidrift cylinder of
                    the skewed straw calculated as a function of theta (second
                    solution);

 P1x, P1y, .... , P2z  are stored in the array POINTS[0-5];  the status is stored
 in *STATUS.

*/




//--------------------------------------------------------------------------


// from the resolving formula on page 23 of Gianluigi's notes, setting  r=0.


  ax = C0x - Oxx;
  ay = C0y - Oyy;
  DELTA = Rr*Rr*(vx*vx+vy*vy) - (vx*ay - vy*ax)*(vx*ay - vy*ax);
  AAA = vx*vx+vy*vy;


  if ( DELTA < 0. ) {
      *STATUS=-2;
  } else if (AAA == 0.){
      *STATUS=-3;
  } else if (DELTA == 0.){
      *STATUS=1;
      POINTS[0] = C0x - vx*(vx*ax + vy*ay)/AAA;
      POINTS[1] = C0y - vy*(vx*ax + vy*ay)/AAA;
      POINTS[2] = C0z - vz*(vx*ax + vy*ay)/AAA;
  }else {
      *STATUS=0;
      DELTA = sqrt(DELTA);
      POINTS[0] = C0x - vx*(vx*ax + vy*ay - DELTA)/AAA;
      POINTS[1] = C0y - vy*(vx*ax + vy*ay - DELTA)/AAA;
      POINTS[2] = C0z - vz*(vx*ax + vy*ay - DELTA)/AAA;
      POINTS[3] = C0x - vx*(vx*ax + vy*ay + DELTA)/AAA;
      POINTS[4] = C0y - vy*(vx*ax + vy*ay + DELTA)/AAA;
      POINTS[5] = C0z - vz*(vx*ax + vy*ay + DELTA)/AAA;
  }

   return;


}



//----------end of function PndTrkCTGeometryCalculations::calculateintersections



//------------------------- begin of function  PndTrkCTGeometryCalculations::CalculateSandZ

void PndTrkCTGeometryCalculations::CalculateSandZ(
	Double_t Oxx,
	Double_t Oyy,
	Double_t Rr,
	Double_t STRAWRESOLUTION,
	Short_t skewnum,
	Double_t info[][7],
	Double_t *WDX,
	Double_t *WDY,
	Double_t *WDZ,
	Double_t S[2],
	Double_t Z[2], //  Zcoordinate of the central wire.
	Double_t Zdrift[2],// drift radius projected onto the Helix.
	Double_t Zerror[2] //  150 micron projected onto the Helix.
	)
{


	Int_t STATUS;
	Short_t i,
		j,
		ii;

	Double_t aaa,
		bbb,
		distance,
		Aellipsis1,
		Bellipsis1,
		LL,
		SkewInclWithRespectToS,
		Rx,
		Ry,
		vx1,
		vy1,
		vz1,
		C0x1,
		C0y1,
		C0z1,
		POINTS1[6];


	 Z[0]= 999999.;
	 Z[1]= 999999.;
         i = skewnum ;

         aaa = sqrt(WDX[i]*WDX[i]+WDY[i]*WDY[i]+ WDZ[i]*WDZ[i]);
         vx1 = WDX[i]/aaa;
         vy1 = WDY[i]/aaa;
         vz1 = WDZ[i]/aaa;
         C0x1 = info[i][0];
         C0y1 = info[i][1];
         C0z1 = info[i][2];


       calculateintersections(Oxx,Oyy,Rr,C0x1,C0y1,C0z1,info[i][3],
                              vx1,vy1,vz1,
                              &STATUS,POINTS1);

       if(STATUS < 0 ) return;



       for( ii=0; ii<2; ii++){
        j=3*ii;
        distance = sqrt(
                  (POINTS1[j]-C0x1)*(POINTS1[j]-C0x1) + 
                  (POINTS1[1+j]-C0y1)*(POINTS1[1+j]-C0y1) + 
                  (POINTS1[2+j]-C0z1)*(POINTS1[2+j]-C0z1) 
                            );


        Rx = POINTS1[j]-Oxx ;   //  x component Radial vector of cylinder of trajectory
        Ry = POINTS1[1+j]-Oyy ;   //  y direction Radial vector of cylinder of trajectory

        aaa = sqrt(Rx*Rx+Ry*Ry);
        SkewInclWithRespectToS = (-Ry*vx1 + Rx*vy1)/aaa ;
        SkewInclWithRespectToS /= Rr;
        bbb = sqrt( SkewInclWithRespectToS*SkewInclWithRespectToS + vz1*vz1);



        //  the tilt direction of this ellipse is (1,0)  when major axis along Z direction



        LL = fabs(vx1*Rx + vy1*Ry);
        if( LL < 1.e-10) continue;
        Aellipsis1 = info[i][3]*aaa/LL;

        Bellipsis1 = info[i][3]/Rr;
        if( distance >= info[i][4]  + Aellipsis1) continue;

//--------------------------
        S[ii] = atan2(POINTS1[j+1]-Oyy, POINTS1[j]-Oxx) ;  // atan2 returns radians in (-pi and +pi]
        if( S[ii] < 0.) S[ii] += 2.*PI;


	Z[ii] = POINTS1[j+2];
	Zdrift[ii] =  Aellipsis1;
	Zerror[ii] = STRAWRESOLUTION*aaa/LL;



//	Double_t rotation1 = 180.*atan2(Tiltdirection1[1],Tiltdirection1[0])/PI;


   }    //  end of    for( ii=0; ii<2; ii++)

	return;
}

//-------------------------  end of function  PndTrkCTGeometryCalculations::CalculateSandZ



//----------star of function PndTrkCTGeometryCalculations::ChooseEntranceExitbis
void PndTrkCTGeometryCalculations::ChooseEntranceExitbis(
	Double_t Oxx,
	Double_t Oyy,
	Short_t  Charge,
	Double_t FiStart,
	Short_t nIntersections,
	Double_t *XintersectionList, //  second index =1 -->inner polygon;
	Double_t *YintersectionList, //  second index =2 -->outer polygon.
	Double_t Xcross[2],	// output
	Double_t Ycross[2]	// output
					)
{

 Short_t
	i,
	j;

 PndTrkMergeSort MergeSort;

// this method works under the hypothesis that there are at least 2 intersections.
	if (nIntersections<2) return;

	  if(Charge > 0) {
		Int_t auxIndex[100];
		Double_t fi[100];
		for( i=0;i<nIntersections;i++){
		  fi[i] = atan2(YintersectionList[i]-Oyy,
				XintersectionList[i]-Oxx);
		  if( fi[i] < 0.) fi[i]  += 2.*PI;
		  if( fi[i] > FiStart) fi[i]  -= 2.*PI;
		  if( fi[i] > FiStart) fi[i] = FiStart;
		  auxIndex[i]=i;
		} // end of for( i=0;i<nIntersections[j];i++)
		MergeSort.Merge_Sort( nIntersections, fi, auxIndex);
		Xcross[0] = XintersectionList[ auxIndex[nIntersections-1] ];
		Ycross[0] = YintersectionList[ auxIndex[nIntersections-1] ];
		Xcross[1] = XintersectionList[ auxIndex[nIntersections-2] ];
		Ycross[1] = YintersectionList[ auxIndex[nIntersections-2] ];

	  } else { // case in which Charge is negative.
		Int_t auxIndex[nIntersections];
		Double_t fi[nIntersections];
		for( i=0;i<nIntersections;i++){
		  fi[i] = atan2(YintersectionList[i]-Oyy,
				XintersectionList[i]-Oxx);
		  if( fi[i] < 0.) fi[i]  += 2.*PI;
		  if( fi[i] < FiStart) fi[i]  += 2.*PI;
		  if( fi[i] < FiStart) fi[i] += FiStart;
		  auxIndex[i]=i;
		} // end of for( i=0;i<nIntersections;i++)
		MergeSort.Merge_Sort( nIntersections, fi, auxIndex);
		Xcross[0] = XintersectionList[ auxIndex[0] ];
		Ycross[0] = YintersectionList[ auxIndex[0] ];
		Xcross[1] = XintersectionList[ auxIndex[1] ];
		Ycross[1] = YintersectionList[ auxIndex[1] ];
	  }




}
//----------end of function PndTrkCTGeometryCalculations::ChooseEntranceExitbis



//-------------------------  begin of function  PndTrkCTGeometryCalculations::Dist_SZ

Double_t PndTrkCTGeometryCalculations::Dist_SZ(
	Double_t Rr,
	Double_t KAPPA,
	Double_t FI0,
	Double_t ZED,
	Double_t S,
	Int_t *nrounds
	)
{

//	Defining :	ZZ = (S-FI0)/KAPPA
//	this method returns the distance (WITH ITS SIGN ) :  ZZ - ZED.  Therefore this number
//	can be negative.
//	Care is taken to calculate this distance properly taking into
//	account that we are dealing with the function  FI = mod(KAPPA*Z + FI0, 2*3.14). 

// the limits of an Int_t are : -2147483648 <= n <= 2147483647

	Double_t aaa,
		ABSdis1,
		dis1,
		dis2,
		dis_segments,
		gap;

	if(fabs(KAPPA) < 1.e-10){
		return -999999999.;
	} else if (fabs(KAPPA)>1.e10) {
		return -ZED;
	}

//	gap = fabs(2.*PI/KAPPA);
	
	aaa = (KAPPA*ZED)/(2.*PI);
	if( aaa<-2147483648.) {
		*nrounds = -2147483647;
	} else if (aaa > 2147483647.) {
		*nrounds = 2147483647;
	} else {
		*nrounds = (Int_t) aaa;
	}

	dis_segments = 2.*PI*Rr/sqrt(1.+KAPPA*KAPPA*Rr*Rr); // distance between
						// two consecutive segments
						//  of trajectory.

	dis1 = fabs( Rr*S -Rr*KAPPA*ZED - Rr*FI0)/sqrt( 1.+KAPPA*KAPPA*Rr*Rr);
	dis1 = fmod(dis1,dis_segments);
	dis2 = dis_segments-dis1; if(dis2<0.)dis2=0.;


	if( dis1 < dis2 )
	{
		return dis1;
	} else {
		return dis2;
	}

}

//-------------------------  end of function  PndTrkCTGeometryCalculations::Dist_SZ



//------------------------- begin of function  PndTrkCTGeometryCalculations::FindDistance

Double_t PndTrkCTGeometryCalculations::FindDistance(
	Double_t Oxx,	//  center from wich distance is calculated
	Double_t Oyy,	//  center from wich distance is calculated
	Double_t Rr,
	Double_t tanlow,
	Double_t tanmid,
	Double_t tanup,
	Double_t alfa,	//  intersection circumference parameter
	Double_t beta,	//  intersection circumference parameter
	Double_t gamma	//  intersection circumference parameter
			)
{


	Short_t i,
		 n;

	Double_t Delta,
		 m[3],
		 q,
		 dist,
		 dist1,
		 dist2,
		 distlow,
		 distmid,
		 distup,
		 totaldist,
		 x1,
		 x2,
		 y1,
		 y2;

int nevento=1;





	m[0] = tanlow;
	m[1] = tanmid;
	m[2] = tanup;

	n=0;
	totaldist=0.;


    for(i=0;i<3;i++){

	if( tanlow<999998.) {
		q = Oyy-m[i]*Oxx;
		Delta = alfa*alfa + beta*beta*m[i]*m[i] - 4.*gamma*m[i]*m[i] - 4.*q*q + 4.*m[i]*q*alfa +
			+ 2.*alfa*beta*m[i] - 4.*beta*q - 4.*gamma;
		if( Delta < 0.){
			dist = -1.;
		} else if (Delta==0.){
			x1 = 0.5*(-alfa - 2.*m[i]*q - beta*m[i] )/(1.+m[i]*m[i]);
			y1 = m[i]*x1+q;
			dist = fabs( sqrt( (Oxx-x1)*(Oxx-x1) + (Oyy-y1)*(Oyy-y1) ) - Rr);
		} else {
			Delta = sqrt(Delta);
			x1 = 0.5*(-alfa - 2.*m[i]*q - beta*m[i] - Delta)/(1.+m[i]*m[i]);
			x2 = 0.5*(-alfa - 2.*m[i]*q - beta*m[i] + Delta)/(1.+m[i]*m[i]);
			y1 = m[i]*x1+q;
			y2 = m[i]*x2+q;
			dist1 = fabs( sqrt( (Oxx-x1)*(Oxx-x1) + (Oyy-y1)*(Oyy-y1) ) - Rr);
			dist2 = fabs( sqrt( (Oxx-x2)*(Oxx-x2) + (Oyy-y2)*(Oyy-y2) ) - Rr);
			if(dist1<dist2){
				dist = dist1;
			} else{
				dist = dist2;
			}
		}
	} else {
		Delta =  beta*beta - 4.*Oxx*Oxx - 4.*Oxx*alfa - 4.*gamma;

		if( Delta < 0.){
			dist = -1.;
		} else if (Delta==0.){
			dist = fabs( fabs(Oyy+0.5*beta) - Rr);
		} else {
			Delta = sqrt(Delta);
			y1 = -0.5*beta + Delta/2.;
			y2 = -0.5*beta - Delta/2.;
			dist1 = fabs( fabs(Oyy+0.5*beta + Delta/2.) - Rr);
			dist2 = fabs( fabs(Oyy+0.5*beta - Delta/2.) - Rr);
			if(dist1<dist2) dist = dist1;
			else  dist = dist2;
		}
	}

	if( dist>-0.5) {
		totaldist+=dist;
		n++;
	}


    }	//   end of  for(i=0;i<3;i++)


	if(n!=3)  totaldist = -1.;
	else	totaldist = totaldist / n;

	return totaldist;

}



//------------------------- end of function  PndTrkCTGeometryCalculations::FindDistance



//----------start  function PndTrkCTGeometryCalculations::FindingParallelTrackAngularRange

void  PndTrkCTGeometryCalculations::FindingParallelTrackAngularRange(
	Double_t oX,
	Double_t oY,
	Double_t Rr,
	Short_t  Charge,
	Double_t *Fi_low_limit,
	Double_t *Fi_up_limit,
	Short_t * status,
	Double_t Rmi,	// Rmin of cylindrical volume intersected by track;
	Double_t Rma	// Rmax of cylindrical volume intersected by track;
	)
{
// -------------- calculate the maximum fi and minimum fi spanned by this track,

// see logbook pag.270; by using the Rmin and Rmax of the straw detector.
// this function works also when circular trajectory in XY doesn't pass through (0,0).

	bool	intersection_inner,
		intersection_outer;
	Double_t	teta1,
			teta2,
			tetavertex,
			a,
			cosT,
			cost,
			cosFi,
			cosfi,
			Fi,
			fi,
			FI0,
			Px,
			Py,
			tmp;


	Rma += 1. ; // add a safety margin.
	Rmi -= 1. ; // add a safety margin.




	a = sqrt(oX*oX+oY*oY);

	//  preliminary condition
	if(a + Rr <= Rmi )	// in this case there might be hits at radius < Rmi.
		 { *status = -1 ;return;}
	if( a >= Rr + Rma || Rr >= a + Rma)  // in this case there can be no hits with radius < Rma.
		 { *status = -2;return;}

	if( a - Rr >= Rmi ) intersection_inner = false; else intersection_inner = true;

	if( a + Rr <= Rma || a - Rr >= Rma  )
		 intersection_outer = false; else intersection_outer = true;

	if( (! intersection_inner) && (! intersection_outer) ){
		*Fi_low_limit = 0.;
		*Fi_up_limit = 2.*PI;
		*status = 1;
		return;
	}

//	now the calculation

	FI0 = atan2(-oY,-oX);
	if( intersection_outer ){
		cosFi = (a*a + Rr*Rr - Rma*Rma)/(2.*Rr*a);
		if(cosFi<-1.) cosFi=-1.; else if(cosFi>1.) cosFi=1.;
		Fi = acos(cosFi);
	}

	if( intersection_inner ){
		cosfi = (a*a + Rr*Rr - Rmi*Rmi)/(2.*Rr*a);
		if(cosfi<-1.) cosfi=-1.; else if(cosfi>1.) cosfi=1.;
		fi = acos(cosfi);
	}


	if( Charge < 0.){ // this particle rotates counterclockwise when looking into the beam
		if( intersection_outer && intersection_inner){
			*Fi_low_limit=FI0 + fi;
			*Fi_up_limit= FI0 +Fi;

		} else if (intersection_inner) {
			*Fi_low_limit=FI0 + fi;
			*Fi_up_limit= FI0 - fi;
		} else {
			*Fi_low_limit=FI0 - Fi;
			*Fi_up_limit= FI0 + Fi;
		}	// end of    if( intersection_outer && intersection_inner



	} else {	// continuation of   if( Charge < 0.)

		if( intersection_outer && intersection_inner){
			*Fi_low_limit=FI0 - Fi;
			*Fi_up_limit= FI0 - fi;

		} else if (intersection_inner) {
			*Fi_low_limit=FI0 + fi;	// must invert because low limit must be < up limit
			*Fi_up_limit= FI0 - fi;
		} else {
			*Fi_low_limit=FI0 - Fi;
			*Fi_up_limit= FI0 + Fi;
		}	// end of    if( intersection_outer && intersection_inner


	}	// end of  if( Charge < 0.)




	if(*Fi_low_limit<0.) {
		*Fi_low_limit=fmod(*Fi_low_limit,2.*PI);
		*Fi_low_limit += 2.*PI;
	} else if (*Fi_low_limit>=2.*PI){
		*Fi_low_limit=fmod(*Fi_low_limit,2.*PI);
	}
	if(*Fi_up_limit<0.) {
		*Fi_up_limit=fmod(*Fi_up_limit,2.*PI);
		*Fi_up_limit += 2.*PI;
	} else if (*Fi_up_limit>=2.*PI){
		*Fi_up_limit=fmod(*Fi_up_limit,2.*PI);
	}

	//	Modify *Fi_up_limit by adding
	//	2PI if it is the case, in order to make *Fi_up_limit > *Fi_low_limit.
	if( *Fi_up_limit < *Fi_low_limit ) *Fi_up_limit += 2.*PI;
	if( *Fi_up_limit < *Fi_low_limit ) *Fi_up_limit = *Fi_low_limit;

	*status = 0;



      return;
}


//---------- end of  function PndTrkCTGeometryCalculations::FindingParallelTrackAngularRange



//----------begin of function PndTrkCTGeometryCalculations::FindIntersectionsOuterCircle

Short_t PndTrkCTGeometryCalculations::FindIntersectionsOuterCircle(
	Double_t oX,
	Double_t oY,
	Double_t Rr,
	Double_t Rma,
	Double_t Xcross[2],
	Double_t Ycross[2]
	)
{

	// return -1 --> non intersection;
	// return 0  --> 2 intersections.

	Double_t	a,
			cosFi,
			Fi,
			FI0;
	a = sqrt(oX*oX+oY*oY);

	// case with no intersections.
	if( a >= Rr + Rma || Rr >= a + Rma ||  a + Rr <= Rma) return -1;


	FI0 = atan2(-oY,-oX);
	cosFi = (a*a + Rr*Rr - Rma*Rma)/(2.*Rr*a);
	if(cosFi<-1.) cosFi=-1.; else if(cosFi>1.) cosFi=1.;
	Fi = acos(cosFi);

	Xcross[0] = oX + Rr*cos(FI0+Fi);
	Ycross[0] = oY + Rr*sin(FI0+Fi);
	Xcross[1] = oX + Rr*cos(FI0-Fi);
	Ycross[1] = oY + Rr*sin(FI0-Fi);

	return 0;
}
//----------end of function PndTrkCTGeometryCalculations::FindIntersectionsOuterCircle


//----------begin of function PndTrkCTGeometryCalculations::FindTrackEntranceExitbiHexagonLeft

Short_t PndTrkCTGeometryCalculations::FindTrackEntranceExitbiHexagonLeft(
	Double_t vgap,
	Double_t Oxx,
	Double_t Oyy,
	Double_t Rr,
	Short_t  Charge,
	Double_t Start[3],
	Double_t ApotemaMin, // Apotema=distance Hexagon side from (0,0).
	Double_t ApotemaMax,
	Double_t Xcross[2],
	Double_t Ycross[2]
	)
{
	Short_t flag;
	Short_t	nIntersections;
	Double_t	FiStart,
			XintersectionList[16], // all the possible intersections
			YintersectionList[16]; // (up to 16 intersections).

// The following is the form of the Left (looking from downstream into the beam) biHexagon
// geometrical shape considered in this method :
//
/*

	      /|
	     / |
	    /  |
	   /  /
	  /  / 
	 /  /  
	/  /   
	|  |   
	|  |   
	|  |   
	|  |   
	\  \   
	 \  \  
	  \  \ 
	   \  \
	    \  |
	     \ |
	      \|

*/
// finding all possible intersections with inner parallel straw region.
// The inner parallel straw region is delimited by two hexagons.

	// flag meaning :
	// -1 -->  track outside outer perimeter;
	// 0 -->  at least 1 intersection with polygon, therefore a possible entry and an exit;
	// 1 -->  track contained completely between the two polygons;


	flag=IntersectionsWithClosedbiHexagonLeft(
		vgap,
		Oxx,
		Oyy,
		Rr,
		ApotemaMin,	// Apotema of the inner Hexagon,
		ApotemaMax,// Apotema of the outer Hexagon.
		&nIntersections,
		XintersectionList, // XintersectionList[..][0] --> inner polygon,
				   // XintersectionList[..][1] --> outer polygon.
		YintersectionList
					);

	// IMPORTANT :
	// this is true because here it is assumed that the track comes from (0,0,0)
	// otherwise the code must be changed!

	if (!(flag == 0)) return flag;
	if( nIntersections<2) return -1;

//-------  among all  possible intersection find the entrance point (Xcross[0], Ycross[0])
//	   and the exit point (Xcross[1], Ycross[1])  of this track.

//-------- the starting point of the track.
	FiStart = atan2( Start[1]-Oyy,Start[0]-Oxx);
	if(FiStart<0.) FiStart+= 2.*PI;
	if(FiStart<0.) FiStart =0.;

	// this method selects the entrance and exit points of the trajectory among all
	// geometrical intersections of the circular trajectory with the straw particular
	// volume.

	// so at this point, the intersections are at least 2.

	ChooseEntranceExitbis(
			Oxx,
			Oyy,
			Charge,
			FiStart,
			nIntersections,
			XintersectionList,
			YintersectionList,
			Xcross,	// output
			Ycross	// output
					);


//----------------------

	return flag;

}

//----------end of function PndTrkCTGeometryCalculations::FindTrackEntranceExitbiHexagonLeft




//----------begin of function PndTrkCTGeometryCalculations::FindTrackEntranceExitbiHexagonRight

Short_t PndTrkCTGeometryCalculations::FindTrackEntranceExitbiHexagonRight(
	Double_t vgap,
	Double_t Oxx,
	Double_t Oyy,
	Double_t Rr,
	Short_t  Charge,
	Double_t Start[3],
	Double_t ApotemaMin, // Apotema=distance Hexagon side from (0,0).
	Double_t ApotemaMax,
	Double_t Xcross[2],
	Double_t Ycross[2]
	)
{
	Short_t flag;
	Short_t	nIntersections;
	Double_t	FiStart,
			XintersectionList[16], // all the possible intersections
			YintersectionList[16]; // (up to 16 intersections).

// The following is the form of the Left (looking from downstream into the beam) biHexagon
// geometrical shape considered in this method :
//
/*
	|\
	| \
	|  \
	 \  \
	  \  \
	   |  |
	   |  |
	   /  /
	  /  /
	 /  /
	|  /
	| /
	|/
*/

// finding all possible intersections with inner parallel straw region.
// The inner parallel straw region is delimited by two hexagons.

	// flag meaning :
	// -1 -->  track outside outer perimeter;
	// 0 -->  at least 1 intersection with polygon, therefore a possible entry and an exit;
	// 1 -->  track contained completely between the two polygons;


	flag=IntersectionsWithClosedbiHexagonRight(
		vgap,
		Oxx,
		Oyy,
		Rr,
		ApotemaMin,	// Apotema of the inner Hexagon,
		ApotemaMax,// Apotema of the outer Hexagon.
		&nIntersections,
		XintersectionList,
		YintersectionList
					);




	// IMPORTANT :
	// this is true because here it is assumed that the track comes from (0,0,0)
	// otherwise the code must be changed!

	if (!(flag == 0)) return flag;
	if( nIntersections<2) return -1;

//-------  among all  possible intersection find the entrance point (Xcross[0], Ycross[0])
//	   and the exit point (Xcross[1], Ycross[1])  of this track.

//-------- the starting point of the track.
	FiStart = atan2( Start[1]-Oyy,Start[0]-Oxx);
	if(FiStart<0.) FiStart+= 2.*PI;
	if(FiStart<0.) FiStart =0.;

	// this method selects the entrance and exit points of the trajectory among all
	// geometrical intersections of the circular trajectory with the straw particular
	// volume.

	// so at this point, the intersections are at least 2.

	ChooseEntranceExitbis(
			Oxx,
			Oyy,
			Charge,
			FiStart,
			nIntersections,
			XintersectionList,
			YintersectionList,
			Xcross,	// output
			Ycross	// output
					);


//----------------------

	return flag;

}


//----------end of function PndTrkCTGeometryCalculations::FindTrackEntranceExitbiHexagonRight



//----------begin function PndTrkCTGeometryCalculations::FindTrackEntranceExitHexagonCircleLeft

Short_t PndTrkCTGeometryCalculations::FindTrackEntranceExitHexagonCircleLeft(
	Double_t Oxx,
	Double_t Oyy,
	Double_t Rr,
	Short_t  Charge,
	Double_t Start[3],
	Double_t ApotemaMin, // Apotema=distance Hexagon side from (0,0).
	Double_t Rma, // outer radius of the Circle.
	Double_t GAP,
	Double_t Xcross[2],
	Double_t Ycross[2]
	)
{

//  This methods finds the intersections between a trajectory coming from (0,0) and
//  parameters  Oxx, Oyy, Rr   with the closed geometrical figure (in XY) formed by the STT
//  external Left semicircle and the Outer STT parallel Left straw semi-Hexagon + Gap for
//  the pellet target target.
//  It returns -1 if there are 0 or 1 intersections, 0 if there are at least 2 intersections.


	Double_t	cosFi,
			theta1,
			theta2,
			Theta1,
			Theta2,
			aaa,
			Fi,
			FI0,
			x1,
			x2,
			y1,
			y2;
//------------------

	Short_t	nIntersectionsCircle,
			nIntersections;
	Double_t	FiStart,
			XintersectionList[12],
			YintersectionList[12]; // all the possible intersections (up to 12 intersections).

// finding all possible intersections with inner parallel straw region.


	Double_t Side_x[] = {	-GAP/2., -GAP/2. , -ApotemaMin, -ApotemaMin, -GAP/2., -GAP/2. },
		 Side_y[] = {	sqrt(Rma*Rma-GAP*GAP/4.),  (2.*ApotemaMin-0.5*GAP)/sqrt(3.),
				ApotemaMin/sqrt(3.),
				-ApotemaMin/sqrt(3.),
				-(2.*ApotemaMin-0.5*GAP)/sqrt(3.), -sqrt(Rma*Rma-GAP*GAP/4.)},
		 a[] =	{1.,		-1./sqrt(3.),	1.,	1./sqrt(3.),	1.},
		 b[] =	{0.,		1.,		0.,	1.,		0.},
		 c[] =	{GAP/2., -2.*ApotemaMin/sqrt(3.),ApotemaMin, 2.*ApotemaMin/sqrt(3.), GAP/2.};

	nIntersections=IntersectionsWithOpenPolygon(
		Oxx,
		Oyy,
		Rr,
		5, //  n. Sides of open Polygon.
		a, //  coefficient of formula :  aX + bY + c = 0 defining the Polygon sides.
		b,
		c,
		Side_x,  // X,Y coordinate of the Sides vertices (in sequence, following
		Side_y,  // the Polygon along.
		XintersectionList, // XintersectionList
		YintersectionList // YintersectionList.
					);




//-------------------------------------------------------------------------
// finding intersections of trajectory [assumed to originate from (0,0) ]
// with outer semicircle, the Left part.

	nIntersectionsCircle=IntersectionsWithGapSemicircle(
		Oxx, // input from trajectory
		Oyy, // input from trajectory
		Rr, // input from trajectory
		GAP, // input, vertical gap in XY plane of STT detector.
		true, // true --> Left semicircle, false --> Right semicircle.
		Rma, // radius of external Stt containment.
		&XintersectionList[nIntersections], // output, X list of intersections (maximal 2).
		&YintersectionList[nIntersections]
						);
	nIntersections += nIntersectionsCircle;

//-------- the starting point of the track.

	if(nIntersections<2) return -1;

	FiStart = atan2( Start[1]-Oyy,Start[0]-Oxx);
	if(FiStart<0.) FiStart+= 2.*PI;
	if(FiStart<0.) FiStart =0.;

	// this method selects the entrance and exit points of the trajectory among all
	// geometrical intersections of the circular trajectory with the straw particular
	// volume.


	ChooseEntranceExitbis(
			Oxx,
			Oyy,
			Charge,
			FiStart,
			nIntersections,
			XintersectionList,
			YintersectionList,
			Xcross,	// output
			Ycross	// output
					);


	return 0;

}

//----------end of function PndTrkCTGeometryCalculations::FindTrackEntranceExitHexagonCircleLeft



//----------begin of function PndTrkCTGeometryCalculations::FindTrackEntranceExitHexagonCircleRight

Short_t PndTrkCTGeometryCalculations::FindTrackEntranceExitHexagonCircleRight(
	Double_t Oxx,
	Double_t Oyy,
	Double_t Rr,
	Short_t  Charge,
	Double_t Start[3],
	Double_t ApotemaMin, // Apotema=distance Hexagon side from (0,0).
	Double_t Rma, // outer radius of the Circle.
	Double_t GAP,
	Double_t Xcross[2],
	Double_t Ycross[2]
	)
{

//  This methods finds the intersections between a trajectory coming from (0,0) and
//  parameters  Oxx, Oyy, Rr   with the closed geometrical figure (in XY) formed by the STT
//  external Left semicircle and the Outer STT parallel Left straw semi-Hexagon + Gap for
//  the pellet target target.
//  It returns -1 if there are 0 or 1 intersections, 0 if there are at least 2 intersections.


	Double_t	cosFi,
			theta1,
			theta2,
			Theta1,
			Theta2,
			aaa,
			Fi,
			FI0,
			x1,
			x2,
			y1,
			y2;
//------------------

	Short_t	nIntersectionsCircle,
			nIntersections;
	Double_t	FiStart,
			XintersectionList[12],
			YintersectionList[12]; // all the possible intersections (up to 10 intersections).

// finding all possible intersections with inner parallel straw region.


	Double_t Side_x[] = {	GAP/2., GAP/2. , ApotemaMin, ApotemaMin, GAP/2., GAP/2. },
		 Side_y[] = {	sqrt(Rma*Rma-GAP*GAP/4.),  (2.*ApotemaMin-GAP)/sqrt(3.),
				ApotemaMin/sqrt(3.),
				-ApotemaMin/sqrt(3.),
				-(2.*ApotemaMin-GAP)/sqrt(3.), -sqrt(Rma*Rma-GAP*GAP/4.)},
		 a[] =	{1.,		1./sqrt(3.),	1.,	-1./sqrt(3.),	1.},
		 b[] =	{0.,		1.,		0.,	1.,		0.},
		 c[] =	{-GAP/2.,-2.*ApotemaMin/sqrt(3.),-ApotemaMin, 2.*ApotemaMin/sqrt(3.), -GAP/2.};

	nIntersections=IntersectionsWithOpenPolygon(
		Oxx,
		Oyy,
		Rr,
		5, //  n. Sides of open Polygon.
		a, //  coefficient of formula :  aX + bY + c = 0 defining the Polygon sides.
		b,
		c,
		Side_x,  // X,Y coordinate of the Sides vertices (in sequence, following
		Side_y,  // the Polygon along.
		XintersectionList, // XintersectionList
		YintersectionList // YintersectionList.
					);



//-------------------------------------------------------------------------
// finding intersections of trajectory [assumed to originate from (0,0) ]
// with outer semicircle, the Left part.

	nIntersectionsCircle=IntersectionsWithGapSemicircle(
		Oxx, // input from trajectory
		Oyy, // input from trajectory
		Rr, // input from trajectory
		GAP, // input, vertical gap in XY plane of STT detector.
		false, // true --> Left semicircle, false --> Right semicircle.
		Rma, // radius of external Stt containment.
		&XintersectionList[nIntersections], // output, X list of intersections (maximal 2).
		&YintersectionList[nIntersections]
						);


	nIntersections += nIntersectionsCircle;

//-------- the starting point of the track.

	if(nIntersections<2) return -1;

	FiStart = atan2( Start[1]-Oyy,Start[0]-Oxx);
	if(FiStart<0.) FiStart+= 2.*PI;
	if(FiStart<0.) FiStart =0.;

	// this method selects the entrance and exit points of the trajectory among all
	// geometrical intersections of the circular trajectory with the straw particular
	// volume.

	ChooseEntranceExitbis(
			Oxx,
			Oyy,
			Charge,
			FiStart,
			nIntersections,
			XintersectionList,
			YintersectionList,
			Xcross,	// output
			Ycross	// output
					);


	return 0;

}

//----------end of function PndTrkCTGeometryCalculations::FindTrackEntranceExitHexagonCircleRight



//----------begin of function PndTrkCTGeometryCalculations::IntersectionCircle_Segment

bool PndTrkCTGeometryCalculations::IntersectionCircle_Segment(
	Double_t a, // coefficients implicit equation.
	Double_t b, // of segment : a*x + b*y + c =0.
	Double_t c,
	Double_t P1x, // point delimiting the segment.
	Double_t P2x, // point delimiting the segment.
	Double_t P1y, // point delimiting the segment.
	Double_t P2y, // point delimiting the segment.
	Double_t Oxx, // center of circle.
	Double_t Oyy,
	Double_t Rr, // Radius of circle.
	Short_t * Nintersections,
	Double_t XintersectionList[2],
	Double_t YintersectionList[2],
	Double_t *distance
	)
{

// this method finds the intersection of a circle with a segment. If the circle
// passes through an endpoint of the segment, that is considered an intersection also.

	bool status;

	Short_t ipossibility;

	Double_t aperp,
		 bperp,
		 cperp,
		 det,
		 distq,
		 dist1,
		 dist2,
		 length,
		 length_segmentq,
		 Rq,
		 x,
		 y,
		 Xintersection,
		 Yintersection;

	*Nintersections=0;
	det = a*a+b*b ;
	if(det < 1.e-20){
		cout<<"from  PndTrkCTGeometryCalculations::IntersectionCircle_Segment :"
			<<" this is not the equation of a segment; a = "<<a<<", b = "
			<<b<<", return false!\n";
		return false;
	}
	//  find if intersection circle - segment is possible.

	// check if there is an intersection (with the squares, it's the same!).
	distq = ( a*Oxx+ b*Oyy + c )*( a*Oxx+ b*Oyy + c )/det;
	*distance = sqrt(distq);
	Rq=Rr*Rr;
	length = Rq - distq;
	if(length <= 0. ) return false; // no intersection between trajectory and this
						// segment.
	length = sqrt(length);
	// coefficients of line perpendicular to input segment and passing for (Oxx, Oyy).
	aperp = -b;
	bperp = a;
	cperp = - aperp*Oxx - bperp*Oyy;

	// find intersection of segment with perpendicular : no need to check if
	// det is different from 0.
	Xintersection = (-bperp*c + b*cperp)/det;
	Yintersection = (-a*cperp + aperp*c)/det;
	det = sqrt(det);
	length_segmentq = (P1x-P2x)*(P1x-P2x) + (P1y-P2y)*(P1y-P2y);

	status = false;
	for(ipossibility=-1;ipossibility<2;ipossibility +=2){
		x = Xintersection + ipossibility*length*b/det ;
		y = Yintersection - ipossibility*length*a/det ;
		if ( (x-P1x)*(x-P1x)+(y-P1y)*(y-P1y) >  length_segmentq
						||
		     (x-P2x)*(x-P2x)+(y-P2y)*(y-P2y) >  length_segmentq
				    )  continue;
		status = true;
		XintersectionList[*Nintersections] = x;
		YintersectionList[*Nintersections] = y;
		(*Nintersections)++;
	}  //  end of  for(ipossibility=-1; ...


	return status;
}

//----------end of function PndTrkCTGeometryCalculations::IntersectionCircle_Segment


//----------begin of function PndTrkCTGeometryCalculations::IntersectionSciTil_Circle
bool PndTrkCTGeometryCalculations::IntersectionSciTil_Circle(
	Double_t DIMENSIONSCITIL,
	Double_t posizSciTilx,
	Double_t posizSciTily,
	Double_t Oxx, // center of circle.
	Double_t Oyy,
	Double_t Rr, // Radius of circle.
	Short_t * Nintersections,
	Double_t XintersectionList[2],
	Double_t YintersectionList[2]
	)
{

 bool intersect;

 Double_t
	distance,
	QQ,
	sqrtRR,
	SIGN;


 QQ = posizSciTilx*posizSciTilx+posizSciTily*posizSciTily;
 sqrtRR=sqrt(QQ);

 if( posizSciTily<0.)  SIGN=-1.;
 else  SIGN=1.;


 intersect = IntersectionCircle_Segment(
	posizSciTilx,
	posizSciTily,
	-QQ,
	posizSciTilx-SIGN*0.5*DIMENSIONSCITIL*posizSciTily/sqrtRR,
	posizSciTilx+SIGN*0.5*DIMENSIONSCITIL*posizSciTily/sqrtRR,
	posizSciTily-SIGN*posizSciTilx*0.5*DIMENSIONSCITIL/sqrtRR,
	posizSciTily+SIGN*posizSciTilx*0.5*DIMENSIONSCITIL/sqrtRR,
	Oxx,
	Oyy,
	Rr,
	Nintersections,  // output
	XintersectionList,  // output
	YintersectionList,  // output
	&distance  // output
						);

	return intersect;
}
//----------end of function PndTrkCTGeometryCalculations::IntersectionSciTil_Circle




//----------start  function PndTrkCTGeometryCalculations::IntersectionsWithClosedbiHexagonLeft

Short_t  PndTrkCTGeometryCalculations::IntersectionsWithClosedbiHexagonLeft(
	//-------- inputs
	Double_t vgap,
	Double_t Oxx,
	Double_t Oyy,
	Double_t Rr,
	Double_t Ami,	// min Apotema of hexagonal  volume intersected by track;
	Double_t Ama,	// max Apotema of hexagonal volume intersected by track;
	//-------- outputs
	Short_t *nIntersections,
	Double_t *XintersectionList,
	Double_t *YintersectionList
	)
{


	// return integer convention :
	// -1 -->  track outside outer perimeter;
	// 0 -->  at least 1 intersection with polygon;
  	// 1 -->  track contained completely between the two polygons;


	//  inner Hexagon --> 0
	//  outer Hexagon --> 1

	bool	internal,
		AtLeast1;

	Short_t i,
		 is,
		 j,
		 Nintersections;

	Double_t aaa,
		 maxdistq,
		 distance,
	//-------------------
	// a,b,c == coefficients of the implicit equations of the 3 sides of the half inner Hexagon
	// plus 3 sides of the half outer Hexagon plus 2 vertical sides corresponding to the Gap :
	//    a*x + b*y +c =0; the numbering of these 8 sides follows the convention of Gianluigi's
	// logbook on page 286.
a[] = {-1./sqrt(3.) ,	1.,	1./sqrt(3.),	1.,	1./sqrt(3.),	1.,	-1./sqrt(3.),	1.},
b[] = {1.,		0.,	1.,		0.,	1.,		0.,	1.,		0.},
c[] = {-2.*Ama/sqrt(3.),Ama,	2.*Ama/sqrt(3.),vgap/2.,2.*Ami/sqrt(3.),Ami,	-2.*Ami/sqrt(3.),vgap/2.},
	//----------------------

		tempX[2],
		tempY[2];

	// side_x and side_y are ordered as side1, side2, ..... in Gianluigi's logbook on page 286.
	Double_t
		side_x[] = { -vgap/2.,	-Ama ,	-Ama,	-vgap/2.,	-vgap/2.,	-Ami,
					-Ami,	-vgap/2.,	-vgap/2.},
		side_y[] = {(-0.5*vgap+2.*Ama)/sqrt(3.),	Ama/sqrt(3.),	-Ama/sqrt(3.),
			    -(-0.5*vgap+2.*Ama)/sqrt(3.),	-(-0.5*vgap+2.*Ami)/sqrt(3.),
			    -Ami/sqrt(3.),	Ami/sqrt(3.),	(-0.5*vgap+2.*Ami)/sqrt(3.),
			    (-0.5*vgap+2.*Ama)/sqrt(3.)	};



//-----------------------

//   find intersections (maximum 16) with the 8 sides.



		AtLeast1 = false;
		*nIntersections =0;
		internal = true;
		maxdistq=-9999.;
		for(is=0; is<8; is++){
			aaa = (side_x[is]-Oxx)*(side_x[is]-Oxx)+(side_y[is]-Oyy)*(side_y[is]-Oyy);
			if(aaa>maxdistq) maxdistq=aaa;
			aaa = (side_x[is+1]-Oxx)*(side_x[is+1]-Oxx)+(side_y[is+1]-Oyy)*(side_y[is+1]-Oyy);
			if(aaa>maxdistq) maxdistq=aaa;
			if ( IntersectionCircle_Segment(
						a[is],
						b[is],
						c[is],
						side_x[is],
						side_x[is+1],
						side_y[is],
						side_y[is+1],
						Oxx,
						Oyy,
						Rr,
						&Nintersections,
						tempX,
						tempY,
						&distance // distance of (Oxx,Oyy) from line
							  // defined by  a*x+b*y+c=0.
							)
			   ){
			   AtLeast1=true;
			   for(j=0;j<Nintersections;j++){
				XintersectionList[ *nIntersections ] =tempX[j];
				YintersectionList[ *nIntersections ] =tempY[j];
				(*nIntersections)++;
			   }
			}	// end of if ( IntersectionCircle_Segment( .....


			// the definition of 'internal' here is when the given Point
			// stays at the same side of the origin (0,0) with respect to
			// the given line of equation   a*x+b*y+c=0.
			if(is<3) {
				internal = internal && IsInternal(Oxx,
								Oyy,
								a[is],
								b[is],
								c[is]
								);
			} else {
				internal = internal && (!IsInternal(Oxx,
								Oyy,
								a[is],
								b[is],
								c[is]
								) );
			}

		} // end of  for(is=0; is<8; is++)


	if( !AtLeast1){
	  if (!internal)  return -1;// trajectory outside polygon.
	  if( maxdistq < Rr*Rr ) return -1;// trajectory outside polygon.
	  return 1;	// trajectory completely contained inside this Polygon.
	}
	return 0;

}

//---------- end of  function PndTrkCTGeometryCalculations::IntersectionsWithClosedbiHexagonLeft







//----------start  function PndTrkCTGeometryCalculations::IntersectionsWithClosedbiHexagonRight

Short_t  PndTrkCTGeometryCalculations::IntersectionsWithClosedbiHexagonRight(
	//-------- inputs
	Double_t vgap,
	Double_t Oxx,
	Double_t Oyy,
	Double_t Rr,
	Double_t Ami,	// min Apotema of hexagonal  volume intersected by track;
	Double_t Ama,	// max Apotema of hexagonal volume intersected by track;
	//-------- outputs
	Short_t *nIntersections,
	Double_t *XintersectionList,
	Double_t *YintersectionList
	)
{


	// return integer convention :
	// -1 -->  track outside outer perimeter;
	// 0 -->  at least 1 intersection with polygon;
  	// 1 -->  track contained completely between the two polygons;


	//  inner Hexagon --> 0
	//  outer Hexagon --> 1

	bool	internal,
		AtLeast1;

	Short_t i,
		 is,
		 j,
		 Nintersections;

	Double_t aaa,
		 maxdistq,
		 distance,
	//-------------------
	// a,b,c == coefficients of the implicit equations of the 3 sides of the half inner Hexagon
	// plus 3 sides of the half outer Hexagon plus 2 vertical sides corresponding to the Gap :
	//    a*x + b*y +c =0; the numbering of these 8 sides foolows the convention of Gianluigi's
	// logbook on page 286.
a[] = {1./sqrt(3.) ,	1.,	-1./sqrt(3.),	1.,	-1./sqrt(3.),	1.,	1./sqrt(3.),	1.},
b[] = {1.,		0.,	1.,		0.,	1.,		0.,	1.,		0.},
c[] = {-2.*Ama/sqrt(3.),-Ama,	2.*Ama/sqrt(3.),-vgap/2.,2.*Ami/sqrt(3.),-Ami,	-2.*Ami/sqrt(3.),-vgap/2.},
	//----------------------

		tempX[2],
		tempY[2];

	// side_x and side_y are ordered as side1, side2, ..... in Gianluigi's logbook on page 286.
	Double_t
		side_x[] = { vgap/2.,	Ama ,	Ama,	vgap/2.,	vgap/2.,	Ami,
					Ami,	vgap/2.,	vgap/2.},
		side_y[] = {(-0.5*vgap+2.*Ama)/sqrt(3.),	Ama/sqrt(3.),	-Ama/sqrt(3.),
			    -(-0.5*vgap+2.*Ama)/sqrt(3.),	-(-0.5*vgap+2.*Ami)/sqrt(3.),
			    -Ami/sqrt(3.),	Ami/sqrt(3.),	(-0.5*vgap+2.*Ami)/sqrt(3.),
			    (-0.5*vgap+2.*Ama)/sqrt(3.)	};



//-----------------------

//   find intersections (maximum 16) with the 8 sides.

		AtLeast1 = false;
		*nIntersections =0;
		internal = true;
		maxdistq=-9999.;
		for(is=0; is<8; is++){
			aaa = (side_x[is]-Oxx)*(side_x[is]-Oxx)+(side_y[is]-Oyy)*(side_y[is]-Oyy);
			if(aaa>maxdistq) maxdistq=aaa;
			aaa = (side_x[is+1]-Oxx)*(side_x[is+1]-Oxx)+(side_y[is+1]-Oyy)*(side_y[is+1]-Oyy);
			if(aaa>maxdistq) maxdistq=aaa;
			if ( IntersectionCircle_Segment(
						a[is],
						b[is],
						c[is],
						side_x[is],
						side_x[is+1],
						side_y[is],
						side_y[is+1],
						Oxx,
						Oyy,
						Rr,
						&Nintersections,
						tempX,
						tempY,
						&distance // distance of (Oxx,Oyy) from line
							  // defined by  a*x+b*y+c=0.
							)
			   )
			{
			   AtLeast1=true;
			   for(j=0;j<Nintersections;j++){
				XintersectionList[ *nIntersections ] =tempX[j];
				YintersectionList[ *nIntersections ] =tempY[j];
				(*nIntersections)++;
			   }
			}	// end of if ( IntersectionCircle_Segment( .....


			// the definition of 'internal' here is when the given Point
			// stays at the same side of the origin (0,0) with respect to
			// the given line of equation   a*x+b*y+c=0.
			if(is<3) {
				internal = internal && IsInternal(Oxx,
								Oyy,
								a[is],
								b[is],
								c[is]
								);
			} else {
				internal = internal && (!IsInternal(Oxx,
								Oyy,
								a[is],
								b[is],
								c[is]
								) );
			}

		} // end of  for(is=0; is<8; is++)


	if( !AtLeast1){
	  if (!internal)  return -1;// trajectory outside polygon.
	  if( maxdistq < Rr*Rr ) return -1;// trajectory outside polygon.
	  return 1;	// trajectory completely contained inside this Polygon.
	}
	return 0;

}

//---------- end of  function PndTrkCTGeometryCalculations::IntersectionsWithClosedbiHexagonRight



//----------start  function PndTrkCTGeometryCalculations::IntersectionsWithClosedPolygon

Short_t  PndTrkCTGeometryCalculations::IntersectionsWithClosedPolygon(
	//-------- inputs
	Double_t Oxx,
	Double_t Oyy,
	Double_t Rr,
	Double_t Rmi,	// min Apotema of hexagonal  volume intersected by track;
	Double_t Rma,	// max Apotema of hexagonal volume intersected by track;
	//-------- outputs
	Short_t nIntersections[2],
	Double_t XintersectionList[][2],
	Double_t YintersectionList[][2]
	)
{


	// return integer convention :
	// -2 -->  track outside outer polygon;
	// -1 -->  track contained completely within inner polygon;
	// 0 -->  at least 1 intersection with inner polygon, at least 1 with outer polygon;
	// 1 -->  track contained completely between the two polygons;
	// 2 -->  track contained completely by larger polygons, with intersections in the smaller;
	// 3 -->  track completely outsiede the small polygon with intersections in the bigger.


	//  inner Hexagon --> 0
	//  outer Hexagon --> 1

	bool	internal[2],
		AtLeast1[2];

	Short_t i,
		 is,
		 j,
		 Nintersections;

	Double_t mindist[2],
		 distance,
	//-------------------
	// a,b,c == coefficients of the implicit equations of the six sides of the Hexagon
	// centered at 0 :   a*x + b*y +c =0; see Gianluigi's logbook on page 277;
	// the coefficient  c  has to be multiplied by Erre.
	// The first side is 
		a[] = { 1./sqrt(3.) , 1., -1./sqrt(3.), 1./sqrt(3.) , 1. , -1./sqrt(3.) } ,
		b[] = { 1., 0. , 1., 1., 0. , 1. },
		c[] = { -2./sqrt(3.) , -1., 2./sqrt(3.), 2./sqrt(3.), 1., -2./sqrt(3.)},
	//----------------------

		Erre[] = {Rmi, Rma},  //  this is the distance from (0,0) 
					// of the SIDES of the Hexagon delimiting the Skew area
		tempX[2],
		tempY[2];

	// both hexagon_side_xlow and hexagon_side_xup must be multiplied by appropriate Erre;
	// sides of Hexagon ordered as a, b, c ..... in Gianluigi's logbook on page 280.
	Double_t
		hexagon_side_x[] = { 0., 1. , 1., 0., -1., -1., 0. },
		hexagon_side_y[] = { 2./sqrt(3.),1./sqrt(3.),-1./sqrt(3.),
					-2./sqrt(3.),-1./sqrt(3.), 1./sqrt(3.), 2./sqrt(3.)};



//-----------------------

//   find intersection with the 6 sides of the small exhagon delimiting the skew straws zone
//   see on page 277 of Gianluigi's logbook.

	// status  =0, at least 1 intersection with inner Hexagon, at least 1 intersection
	// with the outer Hexagon; =1, track contained completely between the
	// two Hexagons; = -1 track contained within inner Hexagon;
	// =-2 track outside outer Hexagon.


	for(i=0;i<2;i++){	// i=0 --> inner Hexagon, i= 1 --> outer Hexagon.
		AtLeast1[i] = false;
		nIntersections[i]=0;
		internal[i] = true;
		mindist[i]=999999.;
		for(is=0; is<6; is++){
			if ( IntersectionCircle_Segment(a[is],
						b[is],
						c[is]*Erre[i],
						hexagon_side_x[is]*Erre[i],
						hexagon_side_x[is+1]*Erre[i],
						hexagon_side_y[is]*Erre[i],
						hexagon_side_y[is+1]*Erre[i],
						Oxx,
						Oyy,
						Rr,
						&Nintersections,
						tempX,
						tempY,
						&distance // distance of (Oxx,Oyy) from line
							  // defined by  a*x+b*y+c=0.
							)
			   )
			{
			   AtLeast1[i]=true;
			   for(j=0;j<Nintersections;j++){
				XintersectionList[ nIntersections[i] ][i] =tempX[j];
				YintersectionList[ nIntersections[i] ][i] =tempY[j];
				nIntersections[i]++;
			   }
			}	// end of if ( IntersectionCircle_Segment( .....

			if(mindist[i]>distance) mindist[i]=distance;

			// the definition of 'internal' here is when the given Point
			// stays at the same side of the origin (0,0) with respect to
			// the given line of equation   a*x+b*y+c=0.
			 internal[i] = internal[i] && IsInternal(Oxx,
								Oyy,
								a[is],
								b[is],
								c[is]*Erre[i]
								);

		} // end of  for(is=0; is<6; is++)
	}  // end of  for(i=0;i<2;i++)


	if( (!AtLeast1[0]) && (!AtLeast1[1]) ){
	  if (!internal[1])  return -2;	// trajectory outside outer Hexagon.
	  if( Rr > mindist[1]) return -2;
	  if( !internal[0]) return 1; // trajectory contained between inner and outer Hexagon.
	  if( mindist[0] >= Rr)  return -1; //  trajectory contained in inner Hexagon.
	  return 1;	// trajectory contained between inner and outer Hexagon.
	} else if (AtLeast1[0] && AtLeast1[1] ){ // continuation of  if( (!AtLeast1[0]) && ...
	  return  0;
	} else if (AtLeast1[0]){
		return 2;
	}

	return 3;
}

//---------- end of  function PndTrkCTGeometryCalculations::IntersectionsWithClosedPolygon



//----------begin of function PndTrkCTGeometryCalculations::IntersectionsWithGapSemicircle

Short_t PndTrkCTGeometryCalculations::IntersectionsWithGapSemicircle(
	Double_t Oxx,
	Double_t Oyy,
	Double_t Rr,
	Double_t GAP,
	bool left,  // if true --> Left semicircle; false --> Right semicircle.
	Double_t Rma,
	Double_t *XintersectionList,
	Double_t *YintersectionList
	)
{

	Short_t	nIntersectionsCircle;

	Double_t	cosFi,
			theta1,
			theta2,
			Theta1,
			Theta2,
			aaa,
			Fi,
			FI0,
			x1,
			x2,
			y1,
			y2;



	nIntersectionsCircle=0;
	aaa = sqrt(Oxx*Oxx+Oyy*Oyy);

	//  preliminary condition for having intersections between trajectory and  Circle.

	if( !( aaa >= Rr + Rma || Rr >= aaa + Rma) &&
		!( aaa + Rr <= Rma || aaa - Rr >= Rma  ) )	{

//	now the calculation

		FI0 = atan2(-Oyy,-Oxx);
		cosFi = (aaa*aaa + Rr*Rr - Rma*Rma)/(2.*Rr*aaa);
		if(cosFi<-1.) cosFi=-1.; else if(cosFi>1.) cosFi=1.;
		Fi = acos(cosFi);


		//  (x1, y1) and (x2,y2) are the intersections between the trajectory
		//  and the circle, in the laboratory reference frame.
		x1 = Oxx+Rr*cos(FI0 - Fi);
		y1 = Oyy+Rr*sin(FI0 - Fi);
		theta1 = atan2(y1,x1); // in this way theta1 is between -PI and PI.
		x2 = Oxx+Rr*cos(FI0 + Fi);
		y2 = Oyy+Rr*sin(FI0 + Fi);
		theta2 = atan2(y2,x2); // in this way theta2 is between -PI and PI.
		//  Theta1, Theta2 = angle of the edges of the outer circle + Gap in the laboratory frame.
		//  Theta1, Theta2 are angles between -PI/2 and +PI/2.
		if(!left){	//  Right (looking into the beam) Semicircle.
			Theta2 = atan2( sqrt(Rma*Rma-GAP*GAP/4.),GAP/2.);
			Theta1 = atan2( -sqrt(Rma*Rma-GAP*GAP/4.),GAP/2.);
			if( Theta1<= theta1 && theta1 <= Theta2 ){
				XintersectionList[nIntersectionsCircle]=x1;
				YintersectionList[nIntersectionsCircle]=y1;
				nIntersectionsCircle++;
			}
			if( Theta1<= theta2 && theta2 <= Theta2 ){
				XintersectionList[nIntersectionsCircle]=x2;
				YintersectionList[nIntersectionsCircle]=y2;
				nIntersectionsCircle++;
			}
		} else {	//  Left (looking into the beam) Semicircle.
			Theta2 = atan2( -sqrt(Rma*Rma-GAP*GAP/4.),-GAP/2.);
			Theta1 = atan2( sqrt(Rma*Rma-GAP*GAP/4.),-GAP/2.);
			if( Theta1<= theta1 || theta1 <= Theta2 ){
				XintersectionList[nIntersectionsCircle]=x1;
				YintersectionList[nIntersectionsCircle]=y1;
				nIntersectionsCircle++;
			}
			if( Theta1<= theta2 || theta2 <= Theta2 ){
				XintersectionList[nIntersectionsCircle]=x2;
				YintersectionList[nIntersectionsCircle]=y2;
				nIntersectionsCircle++;
			}
		}

	}  // end of   if (!( a >= Rr + Rma || .....

//---------------------------- end of calculation intersection with outer circle.

	return nIntersectionsCircle;


}
//----------end of function PndTrkCTGeometryCalculations::IntersectionsWithGapSemicircle



//----------start  function PndTrkCTGeometryCalculations::IntersectionsWithOpenPolygon

Short_t  PndTrkCTGeometryCalculations::IntersectionsWithOpenPolygon(
	//-------- inputs
	Double_t Oxx, // Track parameter
	Double_t Oyy, // Track parameter
	Double_t Rr, // Track parameter
	Short_t nSides, // input, n. of Sides of open Polygon.
	Double_t *a, //  coefficient of formula :  aX + bY + c = 0 defining
	Double_t *b, //  the Polygon sides.
	Double_t *c,
	Double_t *Side_x,  // X,Y coordinate of the Sides vertices (in sequence, following
	Double_t *Side_y,  // the Polygon along.
	//-------- outputs
	Double_t *XintersectionList, // XintersectionList
	Double_t *YintersectionList // YintersectionList.
	)
{

	// this methods returns the n. of intersections.



	Short_t i,
		 is,
		 j,
		 nIntersections,
		 Nintersections;

	Double_t mindist,
		 distance,
	//-------------------
	// a,b,c == coefficients of the implicit equations of the six sides of the Hexagon
	// centered at 0 :   a*x + b*y +c =0; see Gianluigi's logbook on page 277;
	// the coefficient  c  has to be multiplied by Erre.

		tempX[2],
		tempY[2];




//-----------------------

		nIntersections=0;
		for(is=0; is<nSides; is++){
			if ( IntersectionCircle_Segment(a[is],
						b[is],
						c[is],
						Side_x[is],
						Side_x[is+1],
						Side_y[is],
						Side_y[is+1],
						Oxx,
						Oyy,
						Rr,
						&Nintersections,
						tempX,
						tempY,
						&distance // distance of (Oxx,Oyy) from line
							  // defined by  a*x+b*y+c=0.
							)
			   ){
			   for(j=0;j<Nintersections;j++){
				XintersectionList[ nIntersections ] =tempX[j];
				YintersectionList[ nIntersections ] =tempY[j];
				nIntersections += Nintersections;
			   }
			}	// end of if ( IntersectionCircle_Segment( .....


		} // end of  for(is=0; is<nSides; is++)
//	}  // end of  for(i=0;i<2;i++)



	return nIntersections;
}

//---------- end of  function PndTrkCTGeometryCalculations::IntersectionsWithOpenPolygon



//----------begin of function PndTrkCTGeometryCalculations::IsInsideCircle

bool PndTrkCTGeometryCalculations::IsInsideArc(
	Double_t Oxx,
	Double_t Oyy,
	Short_t Charge,
	Double_t Xcross[2],
	Double_t Ycross[2],
	Double_t f  // f should be between 0 and 2PI.
	)
{

	Double_t f1,
		 f2;


	// Xcross[0],Ycross[0] is the point of entrance.


	f1 = atan2(Ycross[0]-Oyy, Xcross[0]-Oxx);
	if(f1<0.) f1+= 2.*PI;
	if(f1<0.) f1= 0.;
	f2 = atan2(Ycross[1]-Oyy, Xcross[1]-Oxx);
	if(f2<0.) f2+= 2.*PI;
	if(f2<0.) f2= 0.;



	if(Charge<0){
		if(f1 > f2 ) f2 +=2.*PI;
		if(f1 > f2 ) f2 = f1;
		if( f>f1){
			if(f<f2) return true; else return false;
		} else {
			f +=2.*PI;
			if(f<f1) f= f1;
			if(f<f2) return true; else return false;
		}
	} else {  // Charge > 0.
		if(f1 < f2 ) f1 +=2.*PI;
		if(f1 < f2 ) f1 = f2;
		if( f>f2){
			if(f<f1) return true; else return false;
		} else {
			f +=2.*PI;
			if(f<f2) f= f2;
			if(f<f1) return true; else return false;
		}
	}	// end of if(Charge<0)
}

//----------end of function PndTrkCTGeometryCalculations::IsInsideCircle




//----------end of function PndTrkCTGeometryCalculations::IsInTargetPipe

bool PndTrkCTGeometryCalculations::IsInTargetPipe(
	Double_t Oxx,
	Double_t Oyy,
	Double_t Rr,
	Double_t fi0,
	Double_t kappa,
	Short_t Charge,
	Double_t gap
	)
{
	Short_t nintersections;

	Double_t	delta,
			XintersectionList[8],
			YintersectionList[8],
			Xcross[2],
			Ycross[2];

	//  find if this track lays in the target pipe volume
	//  for Radius < 15 (supposidly the maximu radius of the Mvd system)
	//  and consequently it cannot have Mvd hits.

	nintersections=0;

	// intersections of trajectory with plane X = gap;
	delta = Rr*Rr - (gap-Oxx)*(gap-Oxx);
	if(delta>0.) {
		XintersectionList[0] = gap;
		YintersectionList[0] = sqrt(delta) + Oyy;
		XintersectionList[1] = gap;
		YintersectionList[1] = -sqrt(delta) + Oyy;
		nintersections+=2;
	}

	// intersections of trajectory with plane X = -gap;
	delta = Rr*Rr - (-gap-Oxx)*(-gap-Oxx);
	if(delta>0.) {
		XintersectionList[nintersections] = -gap;
		YintersectionList[nintersections] = sqrt(delta) + Oyy;
		XintersectionList[nintersections+1] = -gap;
		YintersectionList[nintersections+1] = -sqrt(delta) + Oyy;
		nintersections+=2;
	}


	// intersections of trajectory with plane Z = gap;
		XintersectionList[nintersections] = Oxx + Rr*cos(fi0+kappa*gap);
		YintersectionList[nintersections] = Oyy + Rr*sin(fi0+kappa*gap);
		nintersections++;

	// intersections of trajectory with plane Z = -gap;
		XintersectionList[nintersections] = Oxx + Rr*cos(fi0-kappa*gap);
		YintersectionList[nintersections] = Oyy + Rr*sin(fi0-kappa*gap);
		nintersections++;

	ChooseEntranceExitbis(
				Oxx,
				Oyy,
				Charge,
				fi0,
				nintersections,// n. intersection in input.
				XintersectionList,
				YintersectionList,
				Xcross,	// output
				Ycross	// output
				);

	// here I suppose that the Mvd system CONSERVATIVELY ends at Rmax=5 cm.
	if( fabs(Ycross[0]) > 4. ) return true;
	else  return false;
}
//----------end of function PndTrkCTGeometryCalculations::IsInTargetPipe


//----------star of function PndTrkCTGeometryCalculations::IsInternal

bool PndTrkCTGeometryCalculations::IsInternal(
	Double_t Px,	// point
	Double_t Py,
	Double_t Xtraslation,
	Double_t Ytraslation,
	Double_t Theta
	)
{


	// for explanations see Gianluigi's logbook on page 278-280.

	if( (Xtraslation-Px)*sin(Theta) + (Py-Ytraslation)*cos(Theta) >= 0. ) return true;
	else return false;


}
//----------end of function PndTrkCTGeometryCalculations::IsInternal



ClassImp(PndTrkCTGeometryCalculations);