#include <stdio.h>
//////////////functions////////////////
double trace(double A[2][2]){
 double ans;
 ans=A[0][0]+A[1][1];
 return ans;
}
double det(double A[2][2]){
 double ans;
 ans=A[0][0]*A[1][1]-A[0][1]*A[1][0];
 return ans;
}
int TtoT(double A[2][2],double B[2][2],double C[2][2]){
 for(int i=0;i<2;++i){
    for(int j=0;j<2;++j){
       C[i][j]=0;
       for(int k=0;k<2;++k){
          C[i][j]=C[i][j]+A[i][k]*B[k][j];
       }
       printf("%lf ",C[i][j]);
    }
   printf("\n");
 }
 return 1;
}
int SumTwo(double A[2][2], double B[2][2],double C[2][2]){
 for(int i=0;i<2;++i){
    for(int j=0;j<2;++j){
       C[i][j]=0;
       C[i][j]=A[i][j]+B[i][j];
    }
 }
 return 1;
}
int DiffTwo(double A[2][2], double B[2][2],double C[2][2]){
 for(int i=0;i<2;++i){
    for(int j=0;j<2;++j){
       C[i][j]=0;
       C[i][j]=A[i][j]-B[i][j];
    }
 }
 return 1;
}
int FtoF(double A[4][4], double B[4][4],double C[4][4]){
 printf("FtoF-A*B---\n");

 printf("---A-------------\n");
 printf("%lf	%lf	%lf	%lf\n",A[0][0],A[0][1],A[0][2],A[0][3]);
 printf("%lf	%lf	%lf	%lf\n",A[1][0],A[1][1],A[1][2],A[1][3]);
 printf("%lf	%lf	%lf	%lf\n",A[2][0],A[2][1],A[2][2],A[2][3]);
 printf("%lf	%lf	%lf	%lf\n",A[3][0],A[3][1],A[3][2],A[3][3]);
 printf("---B-------------\n");
 printf("%lf	%lf	%lf	%lf\n",B[0][0],B[0][1],B[0][2],B[0][3]);
 printf("%lf	%lf	%lf	%lf\n",B[1][0],B[1][1],B[1][2],B[1][3]);
 printf("%lf	%lf	%lf	%lf\n",B[2][0],B[2][1],B[2][2],B[2][3]);
 printf("%lf	%lf	%lf	%lf\n",B[3][0],B[3][1],B[3][2],B[3][3]);

 
 printf("---A*B-------------\n");
 for(int i=0;i<4;++i){
    for(int j=0;j<4;++j){
       C[i][j]=0;
       for(int k=0;k<4;++k){
          C[i][j]=C[i][j]+A[i][k]*B[k][j];
       }
       printf("%lf ",C[i][j]);
    }
   printf("\n");
 }
 return 1;
}

int Inverse(double A[4][4],double B[4][4]){
int i,j,k;
double tmp_a[4][4];
double inv_a[4][4];
double buf;
int n=4;

for(i=0;i<n;i++){
  for(j=0;j<n;j++){
    tmp_a[i][j]=A[i][j];
  }
}
for(i=0;i<n;i++){
  for(j=0;j<n;j++){
    inv_a[i][j]=(i==j)?1.0:0.0;
  }
}
for(i=0;i<n;i++){
  buf=1/tmp_a[i][i];
  for(j=0;j<n;j++){
    tmp_a[i][j]*=buf;
    inv_a[i][j]*=buf;
  }
  for(j=0;j<n;j++){
    if(i!=j){
      buf=tmp_a[j][i];
      for(k=0;k<n;k++){
        tmp_a[j][k]-=tmp_a[i][k]*buf;
        inv_a[j][k]-=inv_a[i][k]*buf;
      }
    }
  }
}
for(i=0;i<n;i++){
  for(j=0;j<n;j++){
    B[i][j]=inv_a[i][j];
  }
}
 return 1;
}
int FtoV(double A[4][4], double V[4],double VV[4]){
 int i,j;
 for(i=0;i<4;++i){
   VV[i]=0;  
   for(j=0;j<4;++j){
     VV[i]=VV[i]+A[i][j]*V[j];
   }
 }
 return 1;
}
int TransT4(double A[4][4], double B[4][4]){
 int i,j;
 for(i=0;i<4;++i){
  for(j=0;j<4;++j){
    B[j][i]=A[i][j];
  }
 }
 return 1;
}
double ellipseTheta(double A, double B, double C){
 double ans;
  if(C!=A){
    ans=0.5*atan(2*B/(double)(C-A));
  }else{
    printf("ellipse ERROE\n");//getchar();
  }
 return ans;
}
double slope(double A, double B, double C){
 double ans;
  if(C!=A){
    ans=2*B/(double)(C-A);
  }else{
    printf("ellipse slope ERROR\n");getchar();
  }
 return ans;
}
double ellipseABC(double inpA,double inpB, double inpC, double inpD,double *A, double *B, double *C){
  double ans;
  *A=pow(inpD,2)+pow(inpC,2); 
  *B=-1*(inpD*inpB+inpA*inpC); 
  *C=pow(inpA,2)+pow(inpB,2); 
 return ans;

}
void calcTwiss2021_miniSol(){
 char iname[256];
 int i,j,k,idmy;
 double pi=TMath::Pi();
 int muMinus=0;
///////////////////////////
  FILE *inputFile;
  char parName[256];

 double ax,ay;
 double a1,a2;

//HIATCHI2016_2020Jan typ11
 //ax=-0.9176;//rad/m
 //ay=-0.685;
 //a1=0.59*1;
 //a2=0.354535*1;
 //
//minisol_2021Jan 
ax=-2.160473;//rad/m
ay=-0.716764;
a1=0.295926;
a2=0.098059;

FILE *infile;
  infile=fopen("parm.txt","r");
   fscanf(infile,"%s",&parName);
   fscanf(infile,"%lf,%lf,%lf,%lf",&ax,&ay,&a1,&a2);
   fclose(infile);

double alphaX0=-2.7;
double alphaY0=-0.159;
double betaX0=0.955;
double betaY0=1.92;
alphaX0=0;
alphaY0=0;
betaX0=5;
betaY0=5;
double JJ[2][2],JJ4[4][4];
  JJ[0][0]=0;JJ[0][1]=1;
  JJ[1][0]=-1;JJ[1][1]=0;
  JJ4[0][0]=0; JJ4[0][1]=1;JJ4[0][2]=0; JJ4[0][3]=0;
  JJ4[1][0]=-1;JJ4[1][1]=0;JJ4[1][2]=0; JJ4[1][3]=0; 
  JJ4[2][0]=0; JJ4[2][1]=0;JJ4[2][2]=0; JJ4[2][3]=1;
  JJ4[3][0]=0; JJ4[3][1]=0;JJ4[3][2]=-1;JJ4[3][3]=0; 
double MAT[4][4],MATR[4][4],MATRcheck[4][4],tMAT[4][4],sMAT[4][4],dmy4[4][4];
double PP[2][2],TT[2][2],SS[2][2],QQ[2][2];
double AA[2][2],BB[2][2],DD[4][4];
double RR[4][4],RRinv[4][4],tRR[4][4],sRR[4][4];
double RRini[4][4],RRini_inv[4][4];


////////////////////////////////////
//calc UMat

 double opera_x=ax;
 double opera_y=ay;
 double goalX=-2*opera_x/(double)(pow(opera_x,2)-1);
 double goalY=-2*opera_y/(double)(pow(opera_y,2)-1);
 double opera_xy=a1;
 double opera_xpyp=a2;
 double goalXY=-2*opera_xy/(double)(pow(opera_xy,2)-1);
 double goalXPYP=-2*opera_xpyp/(double)(pow(opera_xpyp,2)-1);

 double fDD[4][4];
//test
fDD[0][0]=1;fDD[0][1]=0;fDD[0][2]=0;fDD[0][3]=0;
//fDD[1][0]=ax;fDD[1][1]=0;fDD[1][2]=0;fDD[1][3]=0;//typA
fDD[1][0]=ax;fDD[1][1]=1;fDD[1][2]=0;fDD[1][3]=0;//typB
fDD[2][0]=0;fDD[2][1]=0;fDD[2][2]=1;fDD[2][3]=0;
//fDD[3][0]=0;fDD[3][1]=0;fDD[3][2]=ay;fDD[3][3]=0;//typA
fDD[3][0]=0;fDD[3][1]=0;fDD[3][2]=ay;fDD[3][3]=1;//typB

double alphaXp=-1*ax*betaX0;
double alphaYp=-1*ay*betaY0;
double betaXp=betaX0*1;
double betaYp=betaY0*1;

//from Text
printf("ax=%lf ay=%lf\n",ax,ay);
printf("alphaX0=%lf alphaY0=%lf betaX0=%lf betaY0=%lf\n",alphaX0,alphaY0,betaX0,betaY0);
printf("alphaXp=%lf alphaYp=%lf betaXp=%lf betaYp=%lf\n",alphaXp,alphaYp,betaXp,betaYp);
getchar();
 int i_1,i_2,i_3,i_4;
 double t_r1,t_r2,t_r3,t_r4,t_mu;
 double inpA,inpB,inpC,inpD,outA,outB,outC;
 double ans,judgexxp,judgeyyp,judgexy,judgexpyp,judge;
 double xxp_r1[100],xxp_r2[100],xxp_r3[100],xxp_r4[100];
 double yyp_r1[100],yyp_r2[100],yyp_r3[100],yyp_r4[100];
 double xy_r1[100],xy_r2[100],xy_r3[100],xy_r4[100];
 double xpyp_r1[100],xpyp_r2[100],xpyp_r3[100];
 int num1,num2,num3,num4;
 num1=50;num2=50;num3=50;num4=50;
 double range[4]={1 ,1 ,1 ,1};//13,5
// double range[4]={0.500 ,0.500 ,0.500 ,0.500};//13,5
double sR[4]= {
//-8.73,-10.99,8.14,10.21//2018Feb3 redo 
//-0.385,-1.37277,-2.954,-2.5518//2018Feb3 redo 
//0.075000,0.207230,0.866000,1.128200
//1.85000,6.65230,3.200000,3.966200 
//-0.385289*-1,-1.37277*-1,-2.891678*-1,-2.489342*-1
//2.085289,3.772770,2.491678,1.689342
//2.865289,5.272770,3.451678,2.289342
//
//-3.8,-5.1,-3.1,1.9
//-6.1,-2.6,-0.6,-0.4
//-5.6,-2.3,-1.12,-0.2//<--keep 20210106
//-5.83,-2.74,-0.83,-0.53//<--keep 20210106
//5.83,2.74,0.83,0.53//<--keep 20210106
7.55,3.27,2.83,1.11//<--keep 20210106
};

 double ResultsR[4]= {0, 0, 0, 0};
 double ResultMu;
  t_r1=sR[0]-range[0]*0.5;
  t_r2=sR[1]-range[1]*0.5;
  t_r3=sR[2]-range[2]*0.5;
  t_r4=sR[3]-range[3]*0.5;
 double dtr1=range[0]/(double)num1;
 double dtr2=range[1]/(double)num2;
 double dtr3=range[2]/(double)num3;
 double dtr4=range[3]/(double)num4;
 double match=0.50E01;
 double minimum=10000;
 int good_frag=-1;
  printf("initial %lf %lf %lf %lf\n",t_r1,t_r2,t_r3,t_r4);
 for(i_1=0;i_1<num1;++i_1){
  t_r1=t_r1+dtr1;
       printf("i_1= %d\n ",i_1);
  t_r2=sR[1]-range[0]*0.5;
//  printf("i_1=%d:: %lf %lf %lf %lf\n",i_1,t_r1,t_r2,t_r3,t_r4);
 for(i_2=0;i_2<num2;++i_2){
  t_r2=t_r2+dtr2;
  //     printf("R2: %lf ",R2);
  //     printf(":%d : %lf\n",i_2,t_r2);getchar();
  t_r3=sR[2]-range[1]*0.5;
  t_r4=sR[3]-range[2]*0.5;
 for(i_3=0;i_3<num3;++i_3){
  t_r3=t_r3+dtr3;
  t_r4=sR[3]-range[3]*0.5;
 for(i_4=0;i_4<num4;++i_4){
  t_r4=t_r4+dtr4;

//  printf("i_1,2,3,4=%d %d %d %d\n",i_1,i_2,i_3,i_4);



   good_frag=-1;
//x-x'
   double symp=t_r1*t_r4-t_r2*t_r3;
   if(symp<1){
   //if(symp<1 && symp>=0){
     t_mu=1-symp;
     t_mu=sqrt(t_mu);
     if(muMinus==1)t_mu=-1*t_mu;
//     printf("t_mu=%lf\n",t_mu);
//x-x'
     inpA=t_mu*fDD[0][0];
     inpB=t_r4*fDD[2][2]-t_r2*fDD[3][2];
     inpC=t_mu*fDD[1][0];
     inpD=t_r1*fDD[3][2]-t_r3*fDD[2][2];
//     printf("inpA,B,C,D=%lf %lf %lf %lf\n",inpA,inpB,inpC,inpD);
     ellipseABC(inpA,inpB,inpC,inpD,&outA,&outB,&outC);
//     printf("outA,B,C=%lf %lf %lf\n",outA,outB,outC);
     //ans= slope(outA,outB,outC);
     ans= ellipseTheta(outA,outB,outC);
     ans=ans*-1;
     ans=tan(ans);
     judgexxp=(ax-ans)/(double)ax;
     if(sqrt(pow(judgexxp,2))<match){
//       printf("R1,R2,R3,R4: %lf %lf %lf %lf\n",R1,R2,R3,R4);
//       printf("good sample X-X':%d %d %d %d: %lf %lf %lf %lf mu=%lf\n",i_1,i_2,i_3,i_4,t_r1,t_r2,t_r3,t_r4,t_mu);
         good_frag=0;
     }else{
  //     printf("Not good sample X-X' %lf %lf %lf %lf mu=%lf\n",t_r1,t_r2,t_r3,t_r4,t_mu);
     }
//y-y'
     inpA=-1*(t_r1*fDD[0][0]+t_r2*fDD[1][0]);
     inpB=t_mu*fDD[2][2];
     inpC=-1*(t_r3*fDD[0][0]+t_r4*fDD[1][0]);
     inpD=t_mu*fDD[3][2];
     ellipseABC(inpA,inpB,inpC,inpD,&outA,&outB,&outC);
     //ans= slope(outA,outB,outC);
     ans= ellipseTheta(outA,outB,outC);
     ans=ans*-1;
     ans=tan(ans);
     //judgeyyp=(goalY-ans)/(double)goalY;
     judgeyyp=(ay-ans)/(double)ay;
     if(sqrt(pow(judgeyyp,2))<match){
       //printf("good sample Y-Y':%d %d %d %d:  %lf %lf %lf %lf mu=%lf\n",i_1,i_2,i_3,i_4,t_r1,t_r2,t_r3,t_r4,t_mu);
     }else{
         good_frag=-1;
   //    printf("Not good sample Y-Y' %lf %lf %lf %lf mu=%lf\n",t_r1,t_r2,t_r3,t_r4,t_mu);
     }
//x-y
     inpA=t_mu*fDD[0][0];
     inpB=t_r4*fDD[2][2]-t_r2*fDD[3][2];
     inpC=-1*(t_r1*fDD[0][0]+t_r2*fDD[1][0]);
     inpD=t_mu*fDD[2][2];
     ellipseABC(inpA,inpB,inpC,inpD,&outA,&outB,&outC);
     //ans= slope(outA,outB,outC);
     ans= ellipseTheta(outA,outB,outC);
     ans=ans*-1;
     ans=tan(ans);
     judgexy=(a1-ans)/(double)a1;
     if(sqrt(pow(judgexy,2))<match){
       //printf("good sample X-Y :%d %d %d %d: %lf %lf %lf %lf mu=%lf\n",i_1,i_2,i_3,i_4,t_r1,t_r2,t_r3,t_r4,t_mu);
     }else{
         good_frag=-1;
    //   printf("Not good sample X-Y %lf %lf %lf %lf mu=%lf\n",t_r1,t_r2,t_r3,t_r4,t_mu);
     }
//x'-y'
     inpA=t_mu*fDD[1][0];
     inpB=t_r1*fDD[3][2]-t_r3*fDD[2][2];
     inpC=-1*(t_r3*fDD[0][0]+t_r4*fDD[1][0]);
     inpD=t_mu*fDD[3][2];
     ellipseABC(inpA,inpB,inpC,inpD,&outA,&outB,&outC);
     //ans= slope(outA,outB,outC);
     ans= ellipseTheta(outA,outB,outC);
     ans=ans*-1;
     ans=tan(ans);
     //judgexpyp=(goalXPYP-ans)/(double)goalXPYP;
     judgexpyp=(a2-ans)/(double)a2;
     if(sqrt(pow(judgexpyp,2))<match){
       //printf("good sample X'-Y' :%d %d %d %d: %lf %lf %lf %lf mu=%lf\n",i_1,i_2,i_3,i_4,t_r1,t_r2,t_r3,t_r4,t_mu);
     }else{
         good_frag=-1;
     //  printf("Not good sample X'-Y' %lf %lf %lf %lf mu=%lf\n",t_r1,t_r2,t_r3,t_r4,t_mu);
     }
     if(good_frag>-1){
       //judge=sqrt(pow(judgeyyp,2)+pow(judgexy,2)); //2018Feb3-keep
       //judge=sqrt(pow(judgexy,2)); //2018Feb3-keep
      judge=sqrt(pow(judgexxp,2)+pow(judgeyyp,2)+pow(judgexy,2)); //2018Feb3-keep
      // judge=sqrt(pow(judgexxp,2)+pow(judgeyyp,2)+pow(judgexy,2)+pow(judgexpyp,2));
       //judge=sqrt(pow(judgexy,2)+pow(judgexpyp,2));
       //judge=sqrt(pow(judgeyyp,2)+pow(judgexy,2)+pow(judgexpyp,2));
       //judge=sqrt(pow(judgexxp,2)+pow(judgeyyp,2));
       //judge=sqrt(pow(judgexy,2));
       if(minimum>judge){
         printf("good sample  :%d %d %d %d: %lf %lf %lf %lf mu=%lf r1r4-r2r3=%lf\n",i_1,i_2,i_3,i_4,t_r1,t_r2,t_r3,t_r4,t_mu,symp);
         printf("judge  :: %lf %lf %lf %lf \n",judgexxp,judgeyyp,judgexy,judgexpyp);
         minimum=judge;
         ResultsR[0]=t_r1;
         ResultsR[1]=t_r2;
         ResultsR[2]=t_r3;
         ResultsR[3]=t_r4;
         ResultMu=t_mu;
       }
     }
     if(i_1==(num1/2-1) && i_2==(num2/2-1) && i_3==(num3/2-1) && i_4==(num4/2-1)){
       printf("R1,R2,R3,R4: %lf %lf %lf %lf\n",sR[0],sR[1],sR[2],sR[3]);
       printf(":%d %d %d %d: %lf %lf %lf %lf mu=%lf\n",i_1,i_2,i_3,i_4,t_r1,t_r2,t_r3,t_r4,t_mu);
       printf("judge=%lf minimum=%lf\n",judge,minimum);
       getchar();
     }
   }//end if(symp<1)
 }
 }
 }
 }
 
//////////RESULT///////////////
printf("---GET RESULT R1-R4\n");
printf("R=%lf %lf %lf %lf mu=%lf\n",ResultsR[0],ResultsR[1],ResultsR[2],ResultsR[3],ResultMu);
printf("R2/R1=%lf R3/R1=%lf R4/R1=%lf\n",ResultsR[1]/(double)ResultsR[0],ResultsR[2]/(double)ResultsR[0],ResultsR[3]/(double)ResultsR[0]);

//MAIN END
//input r1-r4 and mu
double R1,R2,R3,R4;
R1=ResultsR[0];R2= ResultsR[1];R3=ResultsR[2];R4= ResultsR[3];//2018Feb3 keep
//R1=-0.19;R2= -1.0;R3=-2.26;R4= -1.66;//2018Feb3 keep
double  Rmu=1-(R1*R4-R2*R3);
  Rmu=sqrt(Rmu);
printf("R1-4 =%lf %lf %lf %lf\n",R1,R2,R3,R4);
////////////////////////////
 RR[0][0]=Rmu;RR[0][1]=0;  RR[0][2]=-1*R4;RR[0][3]=R2;
 RR[1][0]=0;  RR[1][1]=Rmu;RR[1][2]=R3;   RR[1][3]=-1*R1;
 RR[2][0]=R1; RR[2][1]=R2; RR[2][2]=Rmu;  RR[2][3]=0;
 RR[3][0]=R3; RR[3][1]=R4; RR[3][2]=0;    RR[3][3]=Rmu;

 RRini[0][0]=1; RRini[0][1]=0;RRini[0][2]=0; RRini[0][3]=0;
 RRini[1][0]=0; RRini[1][1]=1;RRini[1][2]=0; RRini[1][3]=0;
 RRini[2][0]=0; RRini[2][1]=0;RRini[2][2]=1; RRini[2][3]=0;
 RRini[3][0]=0; RRini[3][1]=0;RRini[3][2]=0; RRini[3][3]=1;
///symplectic check///
 TransT4(RR,tRR);
 FtoF(tRR,JJ4,dmy4);
 FtoF(dmy4,RR,sRR);
printf("symplectic? RR---\n");
 for(int ii=0;ii<4;++ii){
   for(int jj=0;jj<4;++jj){
     printf("%lf ",sRR[ii][jj]);
   }
   printf("\n");
 } 
if(pow(sRR[0][0],2)>1)getchar();

//////////////////
 double RR_inv[4][4];
 Inverse(RR,RR_inv);
printf(" RR---\n");
 for(int ii=0;ii<4;++ii){
   for(int jj=0;jj<4;++jj){
     printf("%lf ",RR[ii][jj]);
   }
   printf("\n");
 } 
printf(" RR_inv---\n");
 for(int ii=0;ii<4;++ii){
   for(int jj=0;jj<4;++jj){
     printf("%lf ",RR_inv[ii][jj]);
   }
   printf("\n");
 } 
getchar();
 Inverse(RRini,RRini_inv);
/////////////////////////////
 FtoF(RR_inv,fDD,MAT);
///symplectic check///
 TransT4(MAT,tMAT);
 FtoF(tMAT,JJ4,dmy4);
 FtoF(dmy4,MAT,sMAT);
printf("symplectic? MAT---\n");
 for(int ii=0;ii<4;++ii){
   for(int jj=0;jj<4;++jj){
     printf("%lf ",sMAT[ii][jj]);
   }
   printf("\n");
 } 
if(pow(sMAT[0][0],2)>1)getchar();


////check plot//////////
double VecIn[4],VecOut[4],theta;
double cos_alphaX,cos_alphaY;
double sin_alphaX,sin_alphaY;
cos_alphaX=ax;
sin_alphaX=sqrt(1-cos_alphaX*cos_alphaX);
cos_alphaY=ay;
sin_alphaY=sqrt(1-cos_alphaY*cos_alphaY);
double fx[128],fy[128],fz[128];
double fxx[128],fyy[128];
 double fVx[128],fVy[128],fVz[128];
 double fVxx[128],fVyy[128];
 double fy2[128],fVy2[128];
 double rang=0.02;//3cm
 for(int ii=0;ii<128;++ii){
   fx[ii]=(ii-64)/(double)64*rang;
   fz[ii]=0;

   fy[ii]=fx[ii]*a1;
   fVy[ii]=fy[ii]*ay;
   fVx[ii]=fx[ii]*ax;

   fVyy[ii]=fVxx[ii]*a2;
   fyy[ii]=fxx[ii]*a1;
 }

 TGraph *gr_testx;
 TGraph *gr_testy;
 TGraph *gr_testxy,*gr_testxpyp;

 gr_testx = new TGraph(128,fx,fVx);
 gr_testy = new TGraph(128,fy,fVy);
 gr_testxy = new TGraph(128,fx,fy);
 gr_testxpyp = new TGraph(128,fVx,fVy);
////////////CIRCLE/////////////
 int evtNum=100;
 double xc[100],yc[100],pxc[100],pyc[100];
 double xc1[100],yc1[100],pxc1[100],pyc1[100];
 double xc2[100],yc2[100],pxc2[100],pyc2[100];
 double xc0[100],yc0[100],pxc0[100],pyc0[100];
 double xc20[100],yc20[100],pxc20[100],pyc20[100];
 double ranx=2E-3;
 double rany=2E-3;
 double ranpx=7E-5;
 double ranpy=7E-5;

 srand((unsigned) time(NULL));
 double dd1,dd2,dd3,dd4; 
 double theta1,theta2;
 double theta3,theta4;
//cirmod0 x-x' theta1, y-y' theta1
//cirmod01,02 x-x' theta1, y-y' -theta1
//cirmod5,03 x-x' theta1, y-y' theta2
//cirmod6,04 x-y theta1, x'-y' theta2
//cirmod7,05 x theta1, x' theta2 y theta3 y' theta4
 for(int i=0;i<evtNum;++i){
  theta=i*0.5*pi/(double)evtNum*4;
  theta1=2*pi*rand()/(double)RAND_MAX;
  dd1=rand()/(double)RAND_MAX;
  xc[i]=ranx*dd1*cos(theta1);
  xc0[i]=ranx*cos(theta);
  dd2=1;//rand()/(double)RAND_MAX;
  theta3=2*pi*rand()/(double)RAND_MAX;
  yc[i]=rany*dd2*sin(theta2);
  yc0[i]=rany*sin(theta);
  dd3=1;//rand()/(double)RAND_MAX;
  theta2=2*pi*rand()/(double)RAND_MAX;
  pxc[i]=ranpx*dd3*sin(theta1);
  pxc0[i]=ranpx*sin(theta);
  dd4=1;//rand()/(double)RAND_MAX;
  theta4=2*pi*rand()/(double)RAND_MAX;
  pyc[i]=ranpy*dd4*cos(theta2);
  pyc0[i]=ranpy*cos(theta);
// printf("i=%d dd1=%e %e %e %e \n",i,dd1,dd2,dd3,dd4);getchar();
 }
 for(int i=0;i<evtNum;++i){

    VecIn[0]=xc[i];
    VecIn[1]=pxc[i];
    VecIn[2]=yc[i];
    VecIn[3]=pyc[i];
    FtoV(MAT,VecIn,VecOut); 
    xc2[i]=VecOut[0];
    pxc2[i]=VecOut[1];
    yc2[i]=VecOut[2];
    pyc2[i]=VecOut[3];

 }
 for(int i=0;i<evtNum;++i){

    VecIn[0]=xc0[i];
    VecIn[1]=pxc0[i];
    VecIn[2]=yc0[i];
    VecIn[3]=pyc0[i];
    FtoV(MAT,VecIn,VecOut); 
   // FtoV(fDD,VecIn,VecOut); 
    xc20[i]=VecOut[0];
    pxc20[i]=VecOut[1];
    yc20[i]=VecOut[2];
    pyc20[i]=VecOut[3];

 }
 for(int i=0;i<evtNum;++i){

    VecIn[0]=xc[i];
    VecIn[1]=pxc[i];
    VecIn[2]=yc[i];
    VecIn[3]=pyc[i];
    FtoV(fDD,VecIn,VecOut); 
    xc1[i]=VecOut[0];
    pxc1[i]=VecOut[1];
    yc1[i]=VecOut[2];
    pyc1[i]=VecOut[3];

 }

 TGraph *gr_xc,*gr_yc,*gr_xyc,*gr_xpypc;
 TGraph *gr_xc0,*gr_yc0,*gr_xyc0,*gr_xpypc0;
 TGraph *gr_xc1,*gr_yc1,*gr_xyc1,*gr_xpypc1;
 TGraph *gr_xc2,*gr_yc2,*gr_xyc2,*gr_xpypc2;
 TGraph *gr_xc20,*gr_yc20,*gr_xyc20,*gr_xpypc20,*gr_yxpc20,*gr_xypc20;
 gr_xc2 = new TGraph(100,xc2,pxc2);
 gr_yc2 = new TGraph(100,yc2,pyc2);
 gr_xyc2 = new TGraph(100,xc2,yc2);
 gr_xpypc2 = new TGraph(100,pxc2,pyc2);
 gr_xc = new TGraph(100,xc,pxc);
 gr_yc = new TGraph(100,yc,pyc);
 gr_xyc = new TGraph(100,xc,yc);
 gr_xpypc = new TGraph(100,pxc,pyc);
 gr_xc1 = new TGraph(100,xc1,pxc1);
 gr_yc1 = new TGraph(100,yc1,pyc1);
 gr_xyc1 = new TGraph(100,xc1,yc1);
 gr_xpypc1 = new TGraph(100,pxc1,pyc1);

 gr_xc20 = new TGraph(100,xc20,pxc20);
 gr_yc20 = new TGraph(100,yc20,pyc20);
 gr_xyc20 = new TGraph(100,xc20,yc20);
 gr_xpypc20 = new TGraph(100,pxc20,pyc20);
 gr_xypc20 = new TGraph(100,xc20,pyc20);
 gr_yxpc20 = new TGraph(100,yc20,pxc20);
 gr_xc0 = new TGraph(100,xc0,pxc0);
 gr_yc0 = new TGraph(100,yc0,pyc0);
 gr_xyc0 = new TGraph(100,xc0,yc0);
 gr_xpypc0 = new TGraph(100,pxc0,pyc0);
////check plot END//////////

//////////////
gROOT->SetStyle("Plain");
gROOT->ForceStyle();
gStyle->SetOptStat(0);
gStyle->SetCanvasBorderMode(0);
gStyle->SetPadColor(0);
gStyle->SetTitleFillColor(0);
gStyle->SetStatColor(0);

int color=0;
 TCanvas *c2 = new TCanvas("c2","Acceptance c2",200,100,600,600);
 c2->Divide(2,2);
///
 c2->cd(1);
 c2->cd(1)->SetLeftMargin(0.18);
// TH2F* frame2_1 = new TH2F("frame2_1"," ",10, -6, 6, 10, -6, 6);
 TH2F* frame2_1 = new TH2F("frame2_1"," ",10, -2E-2, 2E-2, 10, -2*1E-2, 2*1E-2);
 frame2_1->SetTitle("");
 frame2_1->GetXaxis()->SetLabelSize(0.04);
 frame2_1->GetXaxis()->SetTitleSize(0.05);
 frame2_1->GetYaxis()->SetLabelSize(0.04);
 frame2_1->GetYaxis()->SetTitleSize(0.05);
 frame2_1->GetXaxis()->SetTitle("X (m)");
 frame2_1->GetYaxis()->SetTitle("X' (rad)");
 frame2_1->GetYaxis()->SetTitleOffset(1.3);
 frame2_1->SetStats(0);
 frame2_1->Draw();

   gr_xc0->SetMarkerStyle(7);
   gr_xc0->SetMarkerColor(2);
   gr_xc0->SetLineColor(2);
   gr_xc0->Draw("pl");
   gr_xc1->SetMarkerStyle(7);
   gr_xc1->SetMarkerColor(3);
   gr_xc1->SetLineColor(3);
//   gr_xc1->Draw("pl");
   gr_xc20->SetMarkerStyle(7);
   gr_xc20->SetMarkerColor(4);
   gr_xc20->SetLineColor(4);
   gr_xc20->Draw("p");

   gr_testx->SetMarkerStyle(7);
   gr_testx->SetMarkerColor(3);
   gr_testx->SetLineColor(6);
   gr_testx->Draw("l");


c2->cd(1)->SetGridx();
c2->cd(1)->SetGridy();

 c2->cd(2);
 c2->cd(2)->SetLeftMargin(0.18);
 TH2F* frame2_2 = new TH2F("frame2_2"," ",10, -2E-2, 2E-2, 10, -2E-2, 2E-2);
 frame2_2->SetTitle("");
 frame2_2->GetXaxis()->SetLabelSize(0.04);
 frame2_2->GetXaxis()->SetTitleSize(0.05);
 frame2_2->GetYaxis()->SetLabelSize(0.04);
 frame2_2->GetYaxis()->SetTitleSize(0.05);
 frame2_2->GetXaxis()->SetTitle("Y (m)");
 frame2_2->GetYaxis()->SetTitle("Y' (rad)");
 frame2_2->GetYaxis()->SetTitleOffset(1.3);
 frame2_2->SetStats(0);
 frame2_2->Draw();

   gr_yc0->SetMarkerStyle(7);
   gr_yc0->SetMarkerColor(2);
   gr_yc0->SetLineColor(2);
   gr_yc0->Draw("pl");
   gr_yc1->SetMarkerStyle(7);
   gr_yc1->SetMarkerColor(3);
   gr_yc1->SetLineColor(3);
//   gr_yc1->Draw("pl");
   gr_yc20->SetMarkerStyle(7);
   gr_yc20->SetMarkerColor(4);
   gr_yc20->SetLineColor(4);
   gr_yc20->Draw("pl");
   //

   gr_testy->SetMarkerStyle(7);
   gr_testy->SetMarkerColor(3);
   gr_testy->SetLineColor(6);
   gr_testy->Draw("l");

c2->cd(2)->SetGridx();
c2->cd(2)->SetGridy();

///
///
 c2->cd(3);
 c2->cd(3)->SetLeftMargin(0.18);
 TH2F* frame2_3 = new TH2F("frame2_3"," ",10, -2E-2, 2E-2, 10, -2E-2, 2E-2);
 frame2_3->SetTitle("");
 frame2_3->GetXaxis()->SetLabelSize(0.04);
 frame2_3->GetXaxis()->SetTitleSize(0.05);
 frame2_3->GetYaxis()->SetLabelSize(0.04);
 frame2_3->GetYaxis()->SetTitleSize(0.05);
 frame2_3->GetXaxis()->SetTitle("X (m)");
 frame2_3->GetYaxis()->SetTitle("Y (m)");
 frame2_3->GetYaxis()->SetTitleOffset(1.3);
 frame2_3->SetStats(0);
 frame2_3->Draw();

   gr_xyc0->SetMarkerStyle(7);
   gr_xyc0->SetMarkerColor(2);
   gr_xyc0->SetLineColor(2);
   gr_xyc0->Draw("lp");
   gr_xyc1->SetMarkerStyle(9);
   gr_xyc1->SetMarkerColor(3);
   gr_xyc1->SetLineColor(3);
//   gr_xyc1->Draw("p");
   gr_xyc20->SetMarkerStyle(7);
   gr_xyc20->SetMarkerColor(4);
   gr_xyc20->SetLineColor(4);
   gr_xyc20->Draw("lp");
   gr_xyc2->SetMarkerStyle(7);
   gr_xyc2->SetMarkerColor(7);
   gr_xyc2->SetLineColor(4);
//   gr_xyc2->Draw("p");
//

   gr_testxy->SetMarkerStyle(7);
   gr_testxy->SetMarkerColor(3);
   gr_testxy->SetLineColor(6);
   gr_testxy->Draw("l");

   color=1;


c2->cd(3)->SetGridx();
c2->cd(3)->SetGridy();

 c2->cd(4);
 c2->cd(4)->SetLeftMargin(0.18);
 TH2F* frame2_4 = new TH2F("frame2_4"," ",10, -2E-2, 2E-2, 10, -2E-2, 2E-2);
 frame2_4->SetTitle("");
 frame2_4->GetXaxis()->SetLabelSize(0.04);
 frame2_4->GetXaxis()->SetTitleSize(0.05);
 frame2_4->GetYaxis()->SetLabelSize(0.04);
 frame2_4->GetYaxis()->SetTitleSize(0.05); 
 frame2_4->GetXaxis()->SetTitle("X' (rad)"); 
 frame2_4->GetYaxis()->SetTitle("Y' (rad)"); 
 frame2_4->GetYaxis()->SetTitleOffset(1.3);
 frame2_4->SetStats(0);
 frame2_4->Draw();

   gr_xpypc0->SetMarkerStyle(7);
   gr_xpypc0->SetMarkerColor(2);
   gr_xpypc0->SetLineColor(2);
   gr_xpypc0->Draw("pl");
   gr_xpypc1->SetMarkerStyle(9);
   gr_xpypc1->SetMarkerColor(3);
   gr_xpypc1->SetLineColor(3);
//   gr_xpypc1->Draw("p");
//   gr_xpypc->SetMarkerStyle(9);
//   gr_xpypc->SetMarkerColor(kGray);
//   gr_xpypc->SetLineColor(kGray);
//   gr_xpypc->Draw("pl");
   gr_xpypc20->SetMarkerStyle(7);
   gr_xpypc20->SetMarkerColor(4);
   gr_xpypc20->SetLineColor(4);
   gr_xpypc20->Draw("pl");

   gr_testxpyp->SetMarkerStyle(7);
   gr_testxpyp->SetMarkerColor(3);
   gr_testxpyp->SetLineColor(6);
   gr_testxpyp->Draw("l");
   color=1;

c2->cd(4)->SetGridx();
c2->cd(4)->SetGridy();

////check plot END//////////
printf("R=%lf %lf %lf %lf mu=%lf\n",ResultsR[0],ResultsR[1],ResultsR[2],ResultsR[3],ResultMu);

}