#include<iostream>
#include<cstdlib>
#include<ctime>
#include<vector>
#include<math.h>
//#include</usr/local/include/root/TCanvas.h>
//#include</usr/local/include/root/TGraph.h>
#include<TCanvas.h>
#include<TGraph.h>

using namespace std;

int main()
{
  int NSAMPLES=10000, x=0, y=0;
  double X2=0, X4=0, Y2=0, Y4=0;
  vector<double> TIME, R2X, R2Y, SR2X, SR2Y;
  vector<int> NSTEP;
  NSTEP.push_back(6);
  NSTEP.push_back(10);
  NSTEP.push_back(25);
  NSTEP.push_back(50);
  NSTEP.push_back(75);
  NSTEP.push_back(100);
  NSTEP.push_back(1000);
  NSTEP.push_back(1500);
  NSTEP.push_back(1800);
  NSTEP.push_back(2000);
  /*NSTEP.push_back(2500);
  NSTEP.push_back(3000);
  NSTEP.push_back(3500);
  NSTEP.push_back(5000);
  NSTEP.push_back(10000);
  NSTEP.push_back(20000);
  NSTEP.push_back(50000);
  NSTEP.push_back(100000);
  NSTEP.push_back(500000);
  NSTEP.push_back(1000000);
  NSTEP.push_back(5000000);*/


  srand(time(0));
  //for(int i=0;i<NSTEP;i++) cout << rand()%10 << endl;
  
  //cout << "extern" << endl;
  //ensemble averages vs step length
  for(unsigned int k=0;k<NSTEP.size();k++)
    {
      X2=0; X4=0; Y2=0; Y4=0;
      //one ensemble average with fixed step length NSTEP[k]
      for(int i=0;i<NSAMPLES;i++)
	{
	  x=0; y=0; 
	  //single system evolution 
	  for(int j=0;j<NSTEP[k];j++)
	    {
	      int RND=0;
	      RND = rand()%4;
	      //cout << "RND " << RND << endl;
	      if(RND==0)      {x+=1;}// cout << x << " flag0" << endl;} 
	      else if(RND==1) {y+=1;}// cout << y << " flag1" << endl;}
	      else if(RND==2) {x-=1;}// cout << x << " flag2" << endl;}
	      else if(RND==3) {y-=1;}// cout << y << " flag3" << endl;}
	    }
	  //cout << "x single system "<< x << endl;
	  //cout << "y single system "<< y << endl;
	  //cout << "x^2 single system "<< pow(x,2) << endl;
	  //cout << "y^2 single system "<< pow(y,2) << endl;
	  
	  X2 += pow(x,2);
	  //cout << "X^2 of the ensemble "<< X2 << endl;
	  X4 += pow(x,4);
	  //cout << "X^4 of the ensemble "<< X4 << endl;
	  Y2 += pow(y,2);
	  //cout << "Y^2 of the ensemble "<< Y2 << endl;
	  Y4 += pow(y,4);
	  //cout << "Y^4 of the ensemble "<< Y4 << endl;
	}
      //cout << "NSAMPLES " << NSAMPLES << endl;
      //cout << "X2 " << (X2/NSAMPLES) << endl << "pow " << pow((X2/NSAMPLES),2) << endl;
      R2X.push_back(X2/NSAMPLES);
      X4/=NSAMPLES;
      //cout << "X4 " << X4 << endl;
      R2Y.push_back(Y2/NSAMPLES);
      Y4/=NSAMPLES;
      //cout << "X4 " << X4 << endl;
      SR2X.push_back((X4-pow((X2/NSAMPLES),2))/sqrt(NSAMPLES));
      SR2Y.push_back((Y4-pow((Y2/NSAMPLES),2))/sqrt(NSAMPLES));
    }
  
  TCanvas* c1 = new TCanvas("c1","Graph",200,10,700,500);
  Int_t n = NSTEP.size();
  Double_t ord[n], absc[n];
  for(Int_t i=0;i<n;i++)
    {
      ord[i]=R2X[i];
      absc[i]=NSTEP[i];
    }
  TGraph* gr = new TGraph(n,ord,absc);
  gr->Draw("ACP");
  c1->Update();
  c1->Print("Rw.png");
  
  //for(unsigned int i=0;i<R2X.size();i++)
  // cout << NSTEP[i] <<" "<< R2X[i] <<" " << /*SR2X[i] << " " <<*/ R2Y[i] << " " <</* SR2Y[i] <<*/ endl;
  
  return 0;
}