// Your own function to generate and print the histogram
// Inputs :
//		1. Data array (pointer to float)
//		2. Size of data array (cont int)
//		3. Number of bins (int)
//		4. Lower limit of the histogram (float)
//		5. Upper limit of the histogram (float)
//		
void myhist(float *data, const int N, const int H, float low, float high){

  float range = high - low;			// Range of the histogram
  float width = range/float(H);		// Width of each bins
  float  w[H+2]={0}, X[H+2]={0};		// Arrays containing counts for bins and
									// bin centers

  for(int ii=0; ii<N; ii++){		// Looping over all data points
    int jj = (data[ii]-low)/width;	// Finding the bin for data point
    if (data[ii]<low) jj = -1;		// Underflow check
    if(data[ii]>high)  jj=H;    	// Over flow check
    w[jj+1] = w[jj+1] + 1.;			// Increase the bin count
  }
  for(int jj=0; jj<H+2; jj++){		// Printing the bin entries
    cout<<"myhist bin "<<jj<<" entries = "<<w[jj]<<endl;
  }
//return w;		// Returning the pointer to bin count
}
//----------------------------------------------------------------------------------

//function to compute moments of a distribution
void moments(){
  const int N = 10000;
  const int H = 100;

  float r1[N], w[H+2], X[H+2];
  const float pi = 3.1415;
  TH1F *hr = new TH1F("hr","uniform random", H,0.0,1.0);

  gRandom->SetSeed(43540);

// Generation of the data and filling the ROOT histogram 
  float mean = 0.0;
  for(int ii=0; ii<N; ii++){
    r1[ii] = gRandom->Uniform();
    mean +=r1[ii];     
    hr->Fill(r1[ii]);
  }
  mean = mean/double(N); // Mean calculated from raw data 
  cout<<"Mean = "<<mean<<endl;
hr->Draw();

// Generating the historam using myhist
  cout<<"dumping myhist"<<endl;
  myhist(r1, N, H, 0.0, 1.0);
  cout<<"end myhist"<<endl;


// Calculating mean for deviation from mean
  float x1 = 0.0;
  float x2=0.0;
  float x3=0.0;
  float x4=0.0;
  for(int ii=0; ii<N; ii++){
    x1 += (r1[ii] - mean);
	x2 += pow(r1[ii]-mean,2);
	x3 += pow(r1[ii]-mean,3);
	x4 += pow(r1[ii]-mean,4);
  }
  x1 = x1/double(N); 
  x2 = x2/double(N); 
  x3 = x3/double(N); 
  x4 = x4/double(N); 
  cout<<"Moments from raw data"<<endl;
  cout<<mean<<"   "<<x1<<"   "<<x2<<"   "<<x3<<"   "<<x4<<endl;


  float hmean=0.0,d1 = 0.0,d2 = 0.0, d3 = 0.0, d4 = 0.0;
  for(int j=0; j<H; j++){
    w[j] = hr->GetBinContent(j);  //ROOT command 
    X[j] = hr->GetBinCenter(j);   //ROOT command
    hmean += w[j]*X[j];
  }
  hmean = hmean/double(N);

  for(int j=0; j<H; j++){
    d1 += w[j]*(X[j] - hmean);
    d2 += w[j]*pow((X[j] - hmean),2);
    d3 += w[j]*pow((X[j] - hmean),3);
    d4 += w[j]*pow((X[j] - hmean),4);
  }
  d1 = d1/double(N); 
  d2 = d2/double(N); 
  d3 = d3/double(N); 
  d4 = d4/double(N); 
  cout<<"Moments from histogram data"<<endl;
  cout<<hmean<<"   "<<d1<<"   "<<d2<<"   "<<d3<<"   "<<d4<<endl;
  //return; 
}