#include <stdlib.h>
#include <stdio.h>
#include <math.h>

/*  

Event generator for kinematics of inelastic scattering in two-body kinematics
 1.) Find CM system, assuming target is at rest
 2.) set direction of beam = direction of CM and remember angular deviation of ion 
 3.) do Lorentz transformation to CM
 4.) Calculate E,p of recoils considering different mass
 5.) fill angles theta and phi with random values 
 6.) project momenta to px, py, pz in CM, 
 7.) switch for which ion, mass03 or mass04, definition of forward changes
 8.) Lorentz transformation back to LAB system
 9.) Rotate coord. system back to system with Z = optical axis
10.) Convert to ion-optical coordinates
 variable tof() is filled with scattering angle theta_CM in deg.
 It can be used later as weighting function, 
 e.g. uwfunc=1/sin((tof(2))/180*pi) for homogenous distribution in theta,
   or uwfunc=1/sin(tof(2)*pi/180/2)^4 for Rutherford d_sigma/d_theta
       H.Weick 19.Aug 2009   

------------ how to compile --------------------     
make shared library object
>gcc -fPIC -c inelastic.c  
>gcc -shared -Wl,-soname,inelastic.so -o inelastic.so inelastic.o


---------- how to use in MOCADI -----------------     
* param[0]= mass target nucleus [u],
* param[1]= mass projectile recoil [u],
* param[2]= mass target recoil [u],
* param[3]= switch which recoil to look at (1=target, others=beam)
* param[4]= charge projectile recoil [e],
* more parameters are optional to limit angular range in CM 
* param[5]=theta_min [deg], param[6]=theta_max [deg] 
CALL /u/weick/gicosy/esr/mocadi/inelastic.so ext_beam
1.0073 55.94213 1.0073 0 28 0 180
parameters

*/


int ext_beam(double *in, double *out, double *param, char *option)
{

/*
in[0]=X [cm]   in[4]=energy[MeV/u] in[8] =electrons    in[12]=deltaE[MeV]
in[1]=X'[mrad] in[5]=time [us]     in[9] =nf/nsf       in[13]=reserved
in[2]=Y [cm]   in[6]=mass [amu]    in[10]=range[mg/c2] in[14]=reserved
in[3]=Y'[mrad] in[7]=Z             in[11]=tof [us]
the variable range is varied after the "stop" keyword
the variable deltaE is varied behind energy loss materials (matter, wedge etc.)
*/

double m01, m02, m03, m04 ;     /* rest masses in amu */
double E1, E2, E3, E4 ;         /* energies in MeV    */
double pc1,  pc2,  pc3,  pc4 ;  /* absolute momentum in MeV/c */
double pc1x, pc2x, pc3x, pc4x ; /* momentum in x in MeV/c */
double pc1y, pc2y, pc3y, pc4y ; /* momentum in y in MeV/c */
double pc1z, pc2z, pc3z, pc4z ; /* momentum in z in MeV/c */
double Z1,   Z2,   Z3,   Z4 ;   /* charge in e */
double gamma1, beta_CM, gamma_CM ; /* gamma of m01 in LAB, velocity of CM vs. LAB */
double theta, phi,pcr ;         /* scattering angle in CM, radial momentum in CM */
double help,min,max ;
double alpha, beta ;            /* angle of ion direction vs. optical axis */
double amu = 931.494028 ;       /* amu in MeV */
int i;

/* sort input, masses in MeV */
 m01 = in[6]   *amu ;  
 m02 = param[0]*amu ;
 m03 = param[1]*amu ;
 m04 = param[2]*amu ;
 Z1 = in[7] ;
 Z3 = param[4] ;
 //printf("parameter=%f %f %f %f %f %f\n",param[0],param[1],param[2],param[3],param[4],param[5]);

/* find beta_CM, projectile moving in z, target at rest */
 gamma1 = 1. + in[4]/amu ;
 E1 = m01 *gamma1 ;
 pc1= m01 *sqrt(gamma1*gamma1-1.) ;
 E2 = m02 ;
 pc2 = 0. ;
 beta_CM = (pc1+pc2) / (E1+E2) ; 
 gamma_CM = 1./sqrt(1.-beta_CM*beta_CM) ;
 //printf("m01=%f m02=%f E1=%f beta_CM=%f gamma_CM=%f \n",m01,m02,E1,beta_CM,gamma_CM);

/* remember angle and set momentum in rotated coord. system */
 alpha = atan(in[1]/1000.) ;
 beta  = atan(in[3]/1000.) ;
 pc1x = 0 ;
 pc1y = 0 ;
 pc1z = pc1 ;
 //printf("pc1x=%f pc1y=%f pc1z=%f E1=%f \n",pc1x,pc1y,pc1z,E1);

/* Lorentz transformation to CM */
 help = pc1z ;
 pc1x = pc1x ;
 pc1y = pc1y ;
 pc1z = gamma_CM * (pc1z - beta_CM * E1) ; 
 E1   = gamma_CM * (E1   - beta_CM * help) ;
 E2   = gamma_CM * (E2   - beta_CM * 0) ;
 //printf("pc1x=%f pc1y=%f pc1z=%f E1=%f E2=%f\n",pc1x,pc1y,pc1z,E1,E2);

/* throw dice for scattering angle distributed on surface of sphere */
 help = rand()/((double)RAND_MAX+1) ;
 theta = acos(2.*help -1) ;
 if(param[5]>=0&&param[5]<param[6]&&param[6]<=180) {    /* limited theta range */
   min = (cos(param[5]/180*3.14159265)+1)/2;
   max = (cos(param[6]/180*3.14159265)+1)/2;
   help = help *(max-min) + min ;
   theta = acos(2.*help -1) ;
 }
 else printf("*** warning inconsistent angular range for theta, use 0-180deg.\n");
 help = rand()/((double)RAND_MAX+1) ;
 phi = help *2.*3.14159265 ;
 //printf("theta=%f phi=%f \n",theta,phi);

/* recoil momentum vector in radial coordinates */
 help= pow(E1+E2,4) + pow(m03*m03-m04*m04,2) - 2*pow(E1+E2,2)*(m03*m03+m04*m04) ;
 if(help<0) printf("*** error in kinematics, argument of sqrt<0 \n");
 pcr = 0.5 * sqrt(help) / (E1+E2); /* could be plus or minus */
 //printf("E1=%f E2=%f m03=%f m04=%f \n",E1,E2,m03,m04);
 //printf("pr=%f %f \n",pcr,pc1z);

/* recoil momentum vector in cartesian coordinates */
 pc3x = pcr * cos(phi) * sin(theta) ;
 pc3y = pcr * sin(phi) * sin(theta) ;
 pc3z = pcr * cos(theta) ;
 E3   = sqrt(pcr * pcr + m03 * m03) ;
 pc4x = -pc3x ;
 pc4y = -pc3y ;
 pc4z = -pc3z ;
 E4   = sqrt(pcr * pcr + m04 * m04) ;
 //printf("CM:  pc3x=%f pc3z=%f E3=%f theta=%f pcr=%f \n",pc3x,pc3z,E3,theta*180/3.14159,pcr);

 /* look at scattered target ion or at beam ion */
 if ((0.999<param[3])&&(param[3]<1.001)) {
   pc3x = -pc4x ;
   pc3y = -pc4y ;
   pc3z =  pc4z ;
   E3   =  E4 ; 
   m03  =  m04 ;
 }

/* Lorentz transformation back to LAB */
 help = pc3z ;
 pc3x = pc3x ;
 pc3y = pc3y ;
 pc3z = gamma_CM * (pc3z + beta_CM * E3) ; 
 E3   = gamma_CM * (E3   + beta_CM * help) ;
 help = pc4z ;
 pc4x = pc4x ;
 pc4y = pc4y ;
 pc4z = gamma_CM * (pc4z + beta_CM * E4) ; 
 E4   = gamma_CM * (E4   + beta_CM * help) ;
 //printf("pc3x=%f pc3y=%f pc3z=%f E3=%f E4=%f\n",pc3x,pc3y,pc3z,E3,E4);
 //printf("LAB: pc3x=%f pc3z=%f E3=%f alpha=%f \n",pc3x,pc3z,E3,atan(pc3x/pc3z)*180/3.14159265);

/* rotate beam to original direction, first in x than in y */
 help = pc3x ;
 pc3x = cos(alpha)*help +sin(alpha)*pc3z ;
 pc3z =-sin(alpha)*help +cos(alpha)*pc3z ;
 help = pc3y ;
 pc3y = cos(beta) *help +sin(beta) *pc3z ;
 pc3z =-sin(beta) *help +cos(beta) *pc3z ;
 if(pc3z<0) {
   printf("*** warning backward scattered ions, MOCADI can't handle opposite direction.\n");
 }

/* transform to ion-optical coordinates, angle is tangens*1000 */
 for(i=0;i<14;i++) out[i]=in[i];
 out[1] = pc3x / pc3z *1000 ;
 out[3] = pc3y / pc3z *1000 ;
 out[4] = (E3-m03)/(m03/amu) ;
 out[6] = m03/amu ;
 out[7] = Z3 ;
 out[11] = theta*180/3.14159265 ;
 //printf(" E-kin3=%f A=%f, B=%f \n",out[4],out[1],out[3]);

return(0);
}
 



/* ************ function to rotate coordinate system ******** */
/*  Rotation in x,z and y,z plane                             */
/*       H.Weick 15.June 2009                                 */
/* ************************************************************/

int ext_rotate(double *in, double *out, double *param, char *option)
{

/*
in[0]=X [cm]   in[4]=energy[MeV/u] in[8] =electrons    in[12]=deltaE[MeV]
in[1]=X'[mrad] in[5]=time [us]     in[9] =nf/nsf       in[13]=reserved
in[2]=Y [cm]   in[6]=mass [amu]    in[10]=range[mg/c2] in[14]=reserved
in[3]=Y'[mrad] in[7]=Z             in[11]=tof [us]
the variable range is varied after the "stop" keyword
the variable deltaE is varied behind energy loss materials (matter, wedge etc.)
*/
 double alpha, beta ;             /* angle of rotation vs. optical axis in degree */
 double pc1x, pc1y, pc1z, pc1 ;   /* momentum vector in MeV */
 double pc2x, pc2y, pc2z;         /* momentum vector after rotation1 in MeV */
 double gamma, m0, E ;            /* Lorentz factor, rest mass in MeV, total energy in MeV */
 double amu = 931.494028 ;        /* amu in MeV */
 double help ;
 int i ;

 /* find three momentum vector */
 m0 = in[6]*amu;
 gamma = 1. + in[4]/amu;
 E = gamma * m0;
 pc1= m0 *sqrt(gamma*gamma-1.) ;
 pc1z = pc1 / sqrt(1. +pow(in[1]/1000,2) +pow(in[3]/1000,2)) ;
 pc1x = pc1z * in[1]/1000;
 pc1y = pc1z * in[3]/1000;
 //printf("pc1x=%f pc1y=%f pc1z=%f pc1=%f gamma=%f\n",pc1x,pc1y,pc1z,pc1,gamma);

 /* rotate vector by angle alpha, beta */
 alpha= param[0]/180*3.14159265;
 beta = param[1]/180*3.14159265;
 pc2x = cos(alpha)*pc1x +sin(alpha)*pc1z ;
 pc2y = pc1y ;
 pc2z =-sin(alpha)*pc1x +cos(alpha)*pc1z ;
 help = pc2y ;
 pc2y = cos(beta) *help +sin(beta) *pc2z ;
 pc2z =-sin(beta) *help +cos(beta) *pc2z ;
 //printf("pc2x=%f pc2y=%f pc2z=%f\n",pc2x,pc2y,pc2z);

 /* convert back to ion-optical coordinates for output*/
 for(i=0;i<14;i++) out[i]=in[i];
 out[1] = pc2x / pc2z *1000;
 out[3] = pc2y / pc2z *1000;
 // printf("A=%f B=%f alpha=%f beta=%f A'=%f, B'=%f \n",in[1],in[3],alpha,beta,out[1],out[3]);

return(0);
}
 




/**************************************************************/
/* user function to put all individual ions on z-axis         */
/* then they will pass following matter always perpendicular  */
/*       H.Weick 16.June 2009                                 */
/* ************************************************************/

int ext_straight(double *in, double *out, double *param, char *option)
{
/*
in[0]=X [cm]   in[4]=energy[MeV/u] in[8] =electrons    in[12]=deltaE[MeV]
in[1]=X'[mrad] in[5]=time [us]     in[9] =nf/nsf       in[13]=reserved
in[2]=Y [cm]   in[6]=mass [amu]    in[10]=range[mg/c2] in[14]=reserved
in[3]=Y'[mrad] in[7]=Z             in[11]=tof [us]
the variable range is varied after the "stop" keyword
the variable deltaE is varied behind energy loss materials (matter, wedge etc.)
*/
  int i ;

 /* convert back to ion-optical coordinates for output*/
 for(i=0;i<14;i++) out[i]=in[i];
 out[1] = 0.0 ;
 out[3] = 0.0 ;

return(0);
}
 




/**************************************************************/
/* user function for a special collimator                     */
/* A rectangular catcher which blocks a part of the beam      */
/* The parameters are: x0, y0, x-width, y-width               */
/*       H.Weick 08.July 2009                                 */
/* ************************************************************/

int ext_boxcatcher(double *in, double *out, double *param, char *option)
{
/*
in[0]=X [cm]   in[4]=energy[MeV/u] in[8] =electrons    in[12]=deltaE[MeV]
in[1]=X'[mrad] in[5]=time [us]     in[9] =nf/nsf       in[13]=reserved
in[2]=Y [cm]   in[6]=mass [amu]    in[10]=range[mg/c2] in[14]=reserved
in[3]=Y'[mrad] in[7]=Z             in[11]=tof [us]
the variable range is varied after the "stop" keyword
the variable deltaE is varied behind energy loss materials (matter, wedge etc.)
A negative return value means the particle is lost.
*/
  int i ;
  int flag ;
  double x,y,x0,y0,xw,yw;  

  x = in[0];
  y = in[2];
  x0 = param[0];
  y0 = param[1];
  xw = param[2]/2.;
  yw = param[3]/2.;

  flag = 0 ; 
  if((x>x0-xw&&x<x0+xw)&&(y>y0-yw&&y<y0+yw)) flag=-1;
  for(i=0;i<14;i++) out[i]=in[i];

return(flag);
}