//Calculate OPERA Euler for 3D

#include "TFile.h"
#include "TTree.h"
#include "TCanvas.h"
#include "TFrame.h"
#include "TH1F.h"
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

//void initialTryLsol(){
//	int tt=0;
void initialTryLsol(int tt){

 double gamma_mu,beta_mu;
 double gamma_e,beta_e;
 double e_mass=0.510998928;//MeV
 //double e_mass=0.511;//MeV
 double mu_mass=105.66;//MeV
 double pi=TMath::Pi();
 double light_c=TMath::C();
 double q=1.60217733*1E-19;//[c]
 double GeVkG=1E-26/(double)5.61;
 char comment[256];


 double vx,vy,vz,vv1,vv2;
 double x,y,z,x0,y0,z0;
 double mom;

 FILE *ofile;
 char oname[256];
 //sprintf(oname,"Lsol_test_2021Jan_XY_line.txt");
 //sprintf(oname,"Lsol_test_2021Jan_XY4_8.txt");
 //sprintf(oname,"Lsol_test_2021Jan_XY.txt");
 //sprintf(oname,"Lsol_test_2021Jan_nonXY4_8.txt");
// sprintf(oname,"Lsol_test_2021Jan_base.txt");
 //sprintf(oname,"Lsol_test_2021Jan_base_inv.txt");
 //sprintf(oname,"Lsol_test_2021Jan_base_inv_down.txt");
 //sprintf(oname,"Lsol_test_2021Jan_base_inv_XYZminus2.txt");
 //sprintf(oname,"Lsol_test_2021Jan_base_nonon8.txt");
// sprintf(oname,"Lsol_test_2021Mar_15_old5.txt");
// sprintf(oname,"Lsol_test_2021Mar_15_site5.txt");
 //sprintf(oname,"Lsol_test_2021Mar_22_site.txt");
 //sprintf(oname,"Lsol_test_2021Sep_27_XY.txt");
 sprintf(oname,"Lsol_test.txt");

 if(tt==0)ofile=fopen(oname,"w");
 if(tt!=0)ofile=fopen(oname,"a");

TH1D *h1= new TH1D("h1","check",201,-0.5,1.5);

//31.275588,60.010399,16.426613,1.125152e+10,1.000928e+10,2.601769e+08

x=31.275588;
y=60.010399;
z=16.426613;
vx=1.125152e+10*-1;//injection
vy=1.000928e+10*-1;
vz=2.601769e+08*-1;

double dMom=0;
double dP=0;

//////////////////////
int kmax=10;
//int kmax=1;
//int kmax=12;
 if(tt==0)fprintf(ofile,"%d\n",kmax); //2013May17
FILE *ifileT;
ifileT=fopen("/home/hiromi/OPERA_model_2021_Jan/miniSol/InputTest_Lsol_XY_line4.txt","r");
char sdmy[256];
fscanf(ifileT,"%s",sdmy);
int ii=0;
int idmy;
double ldmy;
double vx0,vy0,vz0;
double cos_theta,sin_theta;
double cos_phi,sin_phi;
double theta,phi;
for(int k=0;k<kmax;++k){

  fscanf(ifileT,"%d,%d,%d,%lf,%lf,%lf,",&ii,&idmy,&idmy,&x,&y,&z);
  fscanf(ifileT,"%lf,%lf,%lf,%lf,%lf,",&vx,&vy,&vz,&ldmy,&ldmy);
  fscanf(ifileT,"%d,%lf,%lf,%lf,",&idmy,&x0,&y0,&z0);
  fscanf(ifileT,"%lf,%lf,%lf\n",&vx0,&vy0,&vz0);
  printf("k=%d x=%lf y=%lf z=%lf\n",k,x,y,z);
  //x=z*-1;
  //z=x*-1;
//  vx=vx*-1;
//  y=y*-1;
//  vy=vy*-1;
//  vz=vz*-1;
  if(pow(k-tt,2)<0.1)break;
  
}//end for

   vv1=sqrt(vx*vx+vy*vy+vz*vz);
   vv2=sqrt(vx*vx+vz*vz);
   printf("---->input vv1 vv2 vy=%9.6e %9.6e %9.6e \n",vv1,vv2,vy);
   printf("---->input vx vy vz=%lf %lf %lf \n",vx/(double)vv1,vy/(double)vv1,vz/(double)vv1);
   printf("---->input vx vy vz=%e %e %e \n",vx,vy,vz);


 cos_theta=vz/(double)vv1;
 sin_theta=sqrt(1-cos_theta*cos_theta);

 printf("vv1=%lf sin_theta=%lf cos_theta=%lf\n",vv1,sin_theta,cos_theta);
 printf("theta=%lf \n",acos(cos_theta)*180/(double)pi);
 printf("theta=%lf \n",asin(sin_theta)*180/(double)pi);
 
 if(sin_theta!=0){
   cos_phi=vx/(double)(vv1*sin_theta);
   sin_phi=-vy/(double)(vv1*sin_theta);
   theta=acos(cos_theta);
   phi=acos(cos_phi);
   if(sin_phi<0)phi=phi*-1;//3rd

   phi=phi*-1;// needed
  fprintf(ofile,"%lf %lf %lf %9.8e %9.8e %9.8e %e\n",
   x,y,z,(phi)*180/(double)pi,(theta)*180/(double)pi,0.0,dMom);//from y=0 (positive mu) and  from y>0(negative mu) 2013Apr22
 }else{
   theta=0;
   phi=0;
  fprintf(ofile,"%lf %lf %lf %9.8e %9.8e %9.8e %e\n",
   x,y,z,(phi-pi)*180/(double)pi,(pi-theta)*180/(double)pi,0.,dMom);
//   printf("sin_theta=0, theta=0, phi is arbitorary\n");
 }



//}//end for
fclose(ofile);
}