#include "iostream"
#include "stdio.h"
#include "TRandom2.h"
#include "TMath.h"
#include "TCanvas.h"
#include "TH2.h"
//#include "TNtupleD.h"
//#include "tipp_tab.h"

//using namespace TMath;

// Coordonnées de la source ponctuelle

// Caractéristiques du cylindre détecteur

double dist(double x1,double x2,double y1,double y2,double z1=0,double z2=0);
int calc_limits(double HCS,double PCS,double alpha0CS,double R0,double DL,
	double &alphaminCS,double &alphamaxCS,double &thetaminCS,double &thetamaxCS);

void tirMGauss(
	//char *filename, //nom du fichier de sauvegarde des données sur les rayons
	int N=20000, //Nombre de rayons générés (différent du nombre de rayons qui intersectent)
	double sigma=1, //écart-type de la source gaussienne
	double xd=5, //position x du cylindre (centre du bas du cylindre)
	double yd=0, //position y du cylindre (centre du bas du cylindre)
	double zd=10, //position z du cylindre (centre du bas du cylindre)
	double dl=8, //hauteur du cylindre
	double r0=2) //rayon du cylindre
{
  FILE *f=fopen("rays.txt","w");//filename,"w");
  TCanvas *myCan = new TCanvas("myCan","Mon Canvas");
  double pi = TMath::Pi();
  double deuxpi= 2*pi;
  double piontwo = pi/2.0;
  double omega=0; // Angle solide à renvoyer
  double err=0; // Erreur associée à l'angle solide calculé
  
//  int N=200000; // Nombre de tirages

  // Définition de la source gaussienne
  double Sx=0,Sy=0,Sz=0;
  
  TRandom2 Sgauss;
  Sgauss.SetSeed(45);
  double sigmax=sigma; // cm
  double mux=Sx; // cm
  double sigmay=sigma;//cm
  double muy=Sy; //cm
  
  // Définition du détecteur
  double Oprimex=xd;
  double Oprimey=yd;
  double Oprimez=zd;
  double DL=dl;
  double R0=r0;

  // Histogramme permettant la visualisation des tirages
  
  TH2D *HistoAlphaTheta = new TH2D("Histo2D","HistoAlphaTheta",100,-pi,pi,100,0,pi);

  HistoAlphaTheta->GetXaxis()->SetTitle("Alpha");
  HistoAlphaTheta->GetYaxis()->SetTitle("Theta");

  TH2D *HistoSource = new TH2D("Source","Source",100,-7*sigmax,7*sigmax,100,-7*sigmay,7*sigmay);

  HistoSource->GetXaxis()->SetTitle("X (cm)");
  HistoSource->GetYaxis()->SetTitle("Y (cm)");

  TRandom2 rand;
  int seed = 20;
  rand.SetSeed(seed);
  
  double *tabtheta,*tabalpha,*tabsx,*tabsy,*tabsz,*imptheta,*impalpha,*impsx,*impsy;
  tabtheta = new double[N];
  tabalpha = new double[N];
  tabsx = new double[N];
  tabsy = new double[N];
  tabsz = new double[N];
  
  int M=0; // Nombre de faisceaux ayant touchés le détecteur
  imptheta = new double[N];
  impalpha = new double[N];
  impsx = new double[N];
  impsy = new double[N];
  
  for (int i=0;i<N;i++){
    tabalpha[i]=rand.Uniform(-pi,pi);
    tabtheta[i]=rand.Uniform(0,pi);
    tabsx[i]=Sgauss.Gaus(mux,sigmax);
    tabsy[i]=Sgauss.Gaus(muy,sigmay);
    tabsz[i]=0;

  }
  
  //calculate angle limits
  double alphaminCS=0;
  double alphamaxCS=0;
  double thetaminCS=0;
  double thetamaxCS=0;

  double PCS = dist(Sx,Oprimex,Sy,Oprimey); // cm
  double alpha0CS=TMath::ATan2(yd-Sy,xd-Sx);
  double HCS=Oprimez-Sz;

  calc_limits(HCS,PCS,alpha0CS,R0,DL,alphaminCS,alphamaxCS,thetaminCS,thetamaxCS);
  
  //determine intersecting rays
  for (int i=0;i<N;i++){
    bool save=0;
    double sx,sy,sz;
    // Coordonnées de la source ponctuelle-> deux tirs gaussiens
    sx = tabsx[i];
    sy = tabsy[i];
    sz = tabsz[i];
    //double sx=Sgauss.Gaus(mux,sigmax);
    //double sy=Sgauss.Gaus(muy,sigmay);
    HistoSource->Fill(sx,sy);
    //double sz = Sz;
   
    // Calcul des distances caractéristiques
    double H=Oprimez-sz; // cm
    double P=dist(sx,Oprimex,sy,Oprimey); // cm
   
    double R0carre = R0*R0; // pour optimisation
    double Pcarre = P*P;
    double alphamax;
    double thetacr;
    double thetamin;
    double thetamax;
    double alphamin; 
    double alpha,theta;
    //alpha=rand.Uniform(-pi,pi);
    //theta=rand.Uniform(0,pi);
    
    alpha=tabalpha[i];
    theta=tabtheta[i];
    double alpha0=TMath::ATan2(yd-sy,xd-sx);
    alpha= alpha-alpha0;

    //Cas S1    
    if (H>0 && P>R0){

      alphamax = TMath::ASin(R0/P);
      alphamin = -alphamax;
      double HDL = H + DL;
      double W1 = alphamax/pi;
      
      float sinalpha = TMath::Sin(alpha);
      if (alpha<=alphamax && alpha>=alphamin){
	thetamax = TMath::ATan((P*TMath::Cos(alpha)+TMath::Sqrt(R0carre-Pcarre*sinalpha*sinalpha))/H);
	thetamin = TMath::ATan((P*TMath::Cos(alpha)-TMath::Sqrt(R0carre-Pcarre*sinalpha*sinalpha))/(HDL));
	
	double W2 = (TMath::Cos(thetamin)-TMath::Cos(thetamax))/2;
	theta=rand.Uniform(0,pi);
	double W = W1*W2;
	omega+= W;
	err += W*W;
	if (theta<=thetamax && theta>=thetamin){
      alpha=alpha+alpha0;
	  HistoAlphaTheta->Fill(alpha,theta);
	  save=1;
      //cout<<"S1 : "<< alphamin<<"<"<<alpha-alpha0<<"<"<<alphamax<<endl;
	}
      }
    }// fin cas S1
    
    //Cas S2
    if(P>R0 && H<0){

      alphamax = TMath::ASin(R0/P);
      alphamin = - alphamax;
      double W1 = alphamax/pi;
      float sinalpha = TMath::Sin(alpha);
      if (alpha<=alphamax && alpha>=alphamin){
	thetamax=piontwo + TMath::ATan(-H/(P*TMath::Cos(alpha)-TMath::Sqrt(R0carre-Pcarre*sinalpha*sinalpha)));
	thetamin = TMath::ATan((P*TMath::Cos(alpha)-TMath::Sqrt(R0carre-Pcarre*sinalpha*sinalpha))/(DL+H));	
	double W2 = (TMath::Cos(thetamin)-TMath::Cos(thetamax))/2;
	double W = W1*W2;
	omega+= W;
	err += W*W;
	if (theta<=thetamax && theta>=thetamin){
      alpha=alpha+alpha0;
	  HistoAlphaTheta->Fill(alpha,theta);
	  save=1;
	  //cout<<"S2 : "<<thetamin<<"<"<<theta<<"<"<<thetamax<<endl; 
	}
      }
      
    }// fin cas S2

    // Cas S3
    if(P<R0 && H>0){

      thetamax = TMath::ATan((R0+P)/H);
      thetamin = 0;
      thetacr = TMath::ATan((R0-P)/H);
      double W2 = (1-TMath::Cos(thetamax))/2; 
      
      if (theta>= thetamin && theta<=thetamax){
	double W1;
	if (theta<=thetacr){ 
	  W1 = 1;
	  // Dans ce cas, tous les alphas passent dans le détecteur
      alpha=alpha+alpha0;
	  HistoAlphaTheta->Fill(alpha,theta);
	  save=1;
	  //cout<<"S3"<<endl; 
	}

	else{ 
	  
	  double tantheta = TMath::Tan(theta);
	  alphamax = TMath::ACos((Pcarre+H*H*tantheta*tantheta-R0carre)/(2*H*P*tantheta));
	  alphamin = -alphamax;
	  W1 = alphamax/pi;

	  if(alpha>=alphamin && alpha<=alphamax){
        alpha=alpha+alpha0;                         
	    HistoAlphaTheta->Fill(alpha,theta);	
        save=1;    
	  }
	}
	double W = W1*W2;
	omega+= W;
	err += W*W;
      }
      
    }// fin cas S3
    if (save){
       impsx[M]=sx;
       impsy[M]=sy;
       imptheta[M]=theta;
       impalpha[M]=alpha;
       M++;
       //fprintf(f,"%lf %lf %lf %lf\n",sx,sy,theta,alpha);
    }
  }// fin for

  // Fin du calcul de l'angle solide et de l'erreur :
  omega = omega/N;
  double k = N-1;
  double Nn = N*k;
  
  err = TMath::Sqrt((err-N*omega*omega)/(Nn));
  omega = 4*pi*omega;
  err = 4*pi*err;
  cout <<"Nb de rayons : "<<N<<endl;
  cout <<"Nb de rayons intersectants : "<<M<<endl;
  cout <<"Angle solide : "<<omega<<endl;
  cout <<"Erreur : "<<err<<endl;

  HistoAlphaTheta->SetFillColor(65);
  myCan->Divide(2,2);
  myCan->cd(1);
  HistoAlphaTheta->Draw("COLZ");
  myCan->cd(2);
  HistoSource->Draw("COLZ");

//activate/deactivate rays
//M=0;
  
  // nb de pts tirés, nb de pts intersectant, EMS source, params cylindre (pos,taille),angle solide effectif,erreur sur l'angle.
  fprintf(f,"%d %d %lf %lf %lf %lf %lf %lf %lf %lf\n\n",N,M+4,sigma,xd,yd,zd,dl,r0,omega,err);
  
  fprintf(f,"%lf %lf %lf %lf\n",Sx,Sy,pi/2,alphaminCS);
  fprintf(f,"%lf %lf %lf %lf\n",Sx,Sy,pi/2,alphamaxCS);
  fprintf(f,"%lf %lf %lf %lf\n",Sx,Sy,thetaminCS,alpha0CS);
  fprintf(f,"%lf %lf %lf %lf\n",Sx,Sy,thetamaxCS,alpha0CS);
  
  for(int i=0;i<M;i++){
    fprintf(f,"%lf %lf %lf %lf\n",impsx[i],impsy[i],imptheta[i],impalpha[i]);
  }
  
  fclose(f);
}

int calc_limits(double HCS,double PCS,double alpha0CS,double R0,double DL,
	double &alphaminCS,double &alphamaxCS,double &thetaminCS,double &thetamaxCS)
{
	double pi = TMath::Pi();
	double deuxpi= 2*pi;
	double piontwo = pi/2.0;

	//limits based on source center
	//case S1
	if (HCS>0 && PCS>R0)
	{
		alphamaxCS = alpha0CS+TMath::ASin(R0/PCS);
		alphaminCS = alpha0CS-TMath::ASin(R0/PCS);
		thetamaxCS = TMath::ATan((PCS+R0)/HCS);
		thetaminCS = TMath::ATan((PCS-R0)/(HCS+DL));
	}
	//case S2
	if(PCS>R0 && HCS<0)
	{
		alphamaxCS = alpha0CS+TMath::ASin(R0/PCS);
		alphaminCS = alpha0CS-TMath::ASin(R0/PCS);
		thetamaxCS = piontwo + TMath::ATan(-HCS/(PCS-R0));
		thetaminCS = TMath::ATan((PCS-R0)/(DL+HCS));
	}
	//case S3
	if(PCS<R0 && HCS>0)
	{
		alphamaxCS = deuxpi;
		alphaminCS = 0;
		thetamaxCS = TMath::ATan((R0+PCS)/(double)HCS);
		thetaminCS = -TMath::ATan((R0-PCS)/(double)HCS);
	}
	return(0);
}

double dist(double x1,double x2,double y1,double y2,double z1,double z2)
{
  double x,y,z;
  x = x2-x1;
  y = y2-y1;
  z = z2-z1;
  return (TMath::Sqrt(x*x+y*y+z*z));
}

/*
double calc_alphamin()
{
      alphamax = TMath::ASin(R0/P);
      alphamin = -alphamax;
      double HDL = H + DL;
      double W1 = alphamax/pi;
      
      float sinalpha = TMath::Sin(alpha);
      if (alpha<=alphamax && alpha>=alphamin){
	thetamax = TMath::ATan((P*TMath::Cos(alpha)+TMath::Sqrt(R0carre-Pcarre*sinalpha*sinalpha))/H);
	thetamin = TMath::ATan((P*TMath::Cos(alpha)+TMath::Sqrt(R0carre-Pcarre*sinalpha*sinalpha))/(HDL));
///////////////////
      alphamax = TMath::ASin(R0/P);
      alphamin = - alphamax;
      double W1 = alphamax/pi;
      float sinalpha = TMath::Sin(alpha);
      if (alpha<=alphamax && alpha>=alphamin){
	thetamax=piontwo + TMath::ATan(-H/(P*TMath::Cos(alpha)-TMath::Sqrt(R0carre-Pcarre*sinalpha*sinalpha)));
	thetamin = TMath::ATan((P*TMath::Cos(alpha)-TMath::Sqrt(R0carre-Pcarre*sinalpha*sinalpha))/(DL+H));	
//////////////////////
    if(P<R0 && H>0){
      thetamax = TMath::ATan((R0+P)/H);
      thetamin = 0;
      thetacr = TMath::ATan((R0-P)/H);
      double W2 = (1-TMath::Cos(thetamax))/2; 
      
      if (theta>= thetamin && theta<=thetamax){
	double W1;
	if (theta<=thetacr){ 
	  W1 = 1;
	  // Dans ce cas, tous les alphas passent dans le détecteur
      alpha=alpha+alpha0;
	  HistoAlphaTheta->Fill(alpha,theta);
	  save=1;
	  //cout<<"S3"<<endl; 
	}

	else{ 
	  
	  double tantheta = TMath::Tan(theta);
	  alphamax = TMath::ACos((Pcarre+H*H*tantheta*tantheta-R0carre)/(2*H*P*tantheta));
	  alphamin = -alphamax;
	  W1 = alphamax/pi;

	  if(alpha>=alphamin && alpha<=alphamax){
        alpha=alpha+alpha0;                         
	    HistoAlphaTheta->Fill(alpha,theta);	

}
double calc_alphamax()
{}
double calc_thetamin()
{}
double calc_thetamax()
{}
*/