#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 cylinder 
     double rmin = 50.0;
     double rmax = 55.0;
     double zmin = -500.0;
     double zmax = 500.0;
     double volume = pi*(rmax*rmax-rmin*rmin)*(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
     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;   // step length = 0.1 cm

//  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));
        if(itry<=5) cout << i << ": eta, phi, theta=" << eta << 
                " , " << phi << " , " << theta << endl;
        double length = 0;

// step along the line from 0.1 cm to 600 cm
        for(int j=0; j<6000 ; j++){
           length = length + step;
           double x = length * sin(theta) * cos(phi);
           double y = length * sin(theta) * sin(phi);
           double rxy = sqrt(x*x + y*y);
           double z = length * cos(theta);
           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 inside the cylinder material
           if( ((rmin < rxy ) && (rxy <= rmax)) &&
               ((zmin < z ) && (z <= zmax)) )
           total_weight = total_weight + density*rxy*rxy*step;
        }
     }
 
//  results
    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;
}