#include <iostream>
#include <iomanip>
#include <cstring>
#include <fstream>
#include <cmath>
#include <unistd.h>
#include <TFile.h>
#include <TTree.h>
#include <TRandom.h>
#include <sys/stat.h>
using namespace std;



float mass,temp;
const float Kb=TMath::K()*1e+7;
const float Pi=TMath::Pi();



	float Pv (float v)
	{
 	 float c1=4/sqrt(Pi);
	 float c2=mass/(2*Kb*temp);
	 float c3=pow(c2,1.5);
	 float vv=v*v;
	 float Pv=c1*c3*vv*exp(-c2*vv);
	 return Pv;
	}



	struct Rec1
	{
	 int Nmol;
	 char Moltype[20];
	 float mass, temp, vmax;
	};

	struct Recn
	{
	 float x, y, z, vx, vy, vz;
	};

 void root()
 {
  Rec1 rec1;
  Recn recn;


  float vmean=0, vrms=0, tmean=0, tmax=0;
  char Moltype[20];
  int Nmol;
  float x, y, z, vx, vy, vz, vmax;
  string out_name,root_name,in_name;
  struct stat controlmaxd, controlmaxr, controltab;
  ofstream maxb,out;
  Nmol=1000000;
  strcpy(Moltype, "hydrogen");
  mass=3.35;
  temp=300;
  vmax=5;
/*  cout << "type Nmol: ";
  cin >> Nmol;
  cout << "type Molecule kind (max 20 chars & no blanks): ";
  cin >> Moltype;
  cout << "type Molecule mass (units 10ˆ-24g): ";
  cin >> mass;
  cout << "typetemperature (K): ";
  cin>>temp;
  cout<<"type vmax ( units 10ˆ5 cm/s):";
  cin >> vmax;
    if(cin.fail())
    {
     cout << "Error!!! Wrong format of input data, STOP NOW!\n";
     exit(0);
    }*/

  mass*=1e-24;
  vmax*=pow(10,5);
  float c0=Kb*(temp/mass);
  float vpeak=sqrt(2*c0);
  float Pvmax=Pv(vpeak);
  float vmeanth=sqrt(8*c0/Pi);
  float vrmsth=sqrt(3*c0);

  rec1.Nmol=Nmol;
  strncpy(rec1.Moltype, Moltype, 20);
  rec1.mass=mass;
  rec1.temp=temp;
  rec1.vmax=vmax;


     cout << "Your data will write, in binary, on a file named ``Maxwell.dat'', if you are ok push 0, else write your file's name: ";
 //    cin >> in_name;
in_name="0";
     if(in_name=="0")
     in_name="Maxwell.dat";
     if(stat(in_name.c_str(), &controlmaxd)!=-1)
     {
      char a;
      cerr << "Maxwell.dat already exist, you want to overwrite?[y,n] ";
      cin >> a;
	 if(a=='n' || a=='0')
	 {
	  exit(0);
         }
 	 else
      	   goto start;
	 
     }
     else
     goto start;

  start:
  maxb.open(in_name,ios::binary);
  maxb.write((char*)&rec1,sizeof(rec1));

     cout << "Your root model will write on a file named ``Maxwell.root'', if you are ok push 0, else write your file's name: ";
//     cin >> root_name;
root_name="0";
     if(root_name=="0")
     root_name="Maxwell.root";
     if(stat(root_name.c_str(), &controlmaxr)!=-1)
     {
      char a;
      cerr << "Maxwell.root already exist, you want to overwrite?[y,n] ";
      cin >> a;
	 if(a=='n'|| a=='0')
	 {
	  exit(0);
         }
     }
  

  TFile maxR(root_name.c_str(),"RECREATE","ftitle");
  TTree *tree=new TTree("tree","root tree with one branch");
  tree->Branch("b1",&recn.x,"x:y:z:vx:vy:vz");
  int igood=1;
  int draws=0;
  TRandom3 nrand;
  const float oneof3=1/3.;

	  while(igood<=Nmol)
	  {
	   ++draws;
	   float v=gRandom->Uniform(vmax);
	   float Pvv=Pv(v);
	   float ratio=Pvv/Pvmax;
	   float test=gRandom->Uniform();
	     switch(ratio>=test)
	     {
	      case 0:
               break;
	      case 1:
	      {
	       double RR[5];
	       nrand.RndmArray(5,RR);
	       float rho=pow(RR[0],oneof3);
	       float phi=RR[1]*2*Pi;
	       float costh=(RR[2]*2)-1.;
	       float sinth=sqrt(1-costh*costh);
	       x=rho*sinth*cos(phi);
	       y=rho*sinth*sin(phi);
	       z=rho*costh;
	       //coseni direttori indipendenti dalle posizioni
	       phi=RR[3]*2*Pi;
	       vz=RR[4]*2-1;
	       sinth=sqrt(1-vz*vz);
	       vx=sinth*cos(phi);
	       vy=sinth*sin(phi);
	       float dis=x*vx+y*vy+z*vz;
	       dis=-dis+sqrt(dis*dis-x*x-y*y-z*z+1);
	       float time=dis/v;
	       tmean+=time;
	       vmean+=v;
	       vrms+=v*v;
	       tmax=(tmax>time)?tmax:time;
	       vx*=v;
	       vy*=v;
	       vz*=v;	       
	       //recnassignment
	       recn.x=x;
	       recn.y=y;
	       recn.z=z;
	       recn.vx=vx;
	       recn.vy=vy;
	       recn.vz=vz;
	       tree->Fill();
	       maxb.write((char*)&recn,sizeof(recn));
	       ++igood;
	       }
	      }
	   }

  float d=1./Nmol;
  tmean*=d;
  vmean*=d;
  vrms=sqrt(vrms*d);
  maxR.Write();
  maxR.Close();
  maxb.close();

     cout << "Your data will write on a file named ``TABLE.dat'', if  you are ok push 0, else write your file's name: ";
     cin >> out_name;
     if(out_name=="0")
     out_name="TABLE.dat";
     if(stat(out_name.c_str(), &controltab)!=-1)
     {
      char a;
      cerr << "TABLE.dat already exist, want you to overwrite?[y,n] ";
      cin >> a;
      if(a=='y' || a=='1')
      out.open(out_name.c_str());
         else if(a=='n' || a=='0')
              {char b;
	       cout<<"Want you to write in append mode?[y,n] ";
               cin >> b;
               if(b=='n' || b=='0')
                exit(0);
               else 
                out.open(out_name.c_str(),ios::app);
	      }
     }
     else
     out.open(out_name.c_str());

  string TIT="Nmol\t\t\tMoleculeName\t\tMass(g)\t\t\tTemp(K)\t\t\tVmax(cm/s)";
  string TIT2="Vmeanth\t\t\tVmeansqth\t\tVmeangen\t\tVmeansqgen";
  string TIT3="Tmean\t\t\tTmax";

  out.flags(ios::scientific);
  out<<"Input data:"<<endl;
  out<<setw(10)<<setprecision(3)<<TIT<<endl;
  out<<Nmol<<"\t\t\t"<<Moltype<<"\t\t"<<mass<<"\t\t"<<temp<<"\t\t"<<vmax<<endl;
  out<<endl;
  out<<"efficiency Hit-or-Miss:"<<(float(Nmol)/draws)*100<<"%"<<endl;
  out<<endl;
  out<<TIT2<<endl;
  out<<vmeanth<<"\t\t"<<vrmsth<<"\t\t"<<vmean<<"\t\t"<<vrms<<endl;
  out<<endl;
  out<<TIT3<<endl;
  out<<tmean<<"\t\t"<<tmax<<endl;
  out.close();
  
 }