//2004.9.20 updated
#include <iostream.h>
#include <stdlib.h>
#include <math.h>
using namespace std;

const double pi = 3.141592;
const double twopi = 2.0 * pi;
const double eta_min = -5.0;
const double eta_max = 5.0;
const double phi_min = 0.0;
const double phi_max = twopi;

int main()
{
// define a box object
// one can change these parameters, make sure your box is inside the
// scanned eta-phi region.
     double xmin = -10.;
     double ymin = 20.;
     double zmin = -20.;
     double xmax = 15.0;
     double ymax = 30.0;
     double zmax = 50.0;
     double volume = (xmax-xmin)*(ymax-ymin)*(zmax-zmin);
     double density = 2.70;
     double weight_expected = density * volume;
     cout << "Monte Carlo program to calculate weight"<< endl;
     cout << "weight expected = " << weight_expected << " gram" << endl;

// set seed number for randum number generator
     long zzz=12459888;
     srand48(zzz);
     cout << "set seed number to : " << zzz << endl;

// strat start Monte Carlo
     double total_weight = 0.0;
     int itry = 0;
     for (int n=0; n<30; n++){
        int nmax = 10000;
        double step = 0.1;
//  repeat nmax times
        for(int i=0; i<nmax; i++){
           itry++;
           double eta = (eta_max-eta_min)*drand48() + eta_min ;
           double phi = drand48() * twopi ;
           double theta = 2.0*atan(exp(-eta));
	   double fx = sin(theta) * cos(phi);
           double fy = sin(theta) * sin(phi);
           double fz = cos(theta);
           if(itry<=5) cout << i << ": eta, phi, theta=" << eta << 
                " , " << phi << " , " << theta << endl;
           double length = 0;

// step along the line from 0.1 cm to 100 cm
           for(int j=0; j<1000 ; j++){
              length = length + step;
              double x = length * fx;
              double y = length * fy;
              double rxy2 = x*x + y*y;
              double z = length * fz;
              double totall = sqrt(x*x + y*y + z*z);
              if(itry <= 1 && j < 5)cout << "x,y,z,l=" << x <<" , "<< y << 
                  " , "<<z << " , "<<totall<<endl;

// check if the step point is now inside the box
              if( ((xmin < x ) && (x <= xmax)) &&
                  ((ymin < y ) && (y <= ymax)) &&
                  ((zmin < z ) && (z <= zmax)) )
              total_weight = total_weight + density*rxy2*step;
           }
        }
        double weight_obtained = total_weight * 
                     (eta_max-eta_min)  * (phi_max-phi_min) / itry ;
        double ratio = weight_obtained / weight_expected;
        cout << " itry=" << itry << " total_weight=" << weight_obtained << 
             " ratio=" << ratio << endl;
     }
     return 0;
}