// class Box,  2004.9.16  T. Kondo
#include <iostream.h>
#include <stdlib.h>
#include <math.h>
using namespace std;
const double pi = 3.141592;
const double twopi = 2.0 * pi;

class Box
{
     private:
	double xcenter[3] ;
        double xHalfLength[3] ;
	double density ;
     public:
	Box(double xh, double yh, double zh, double density0);
        void locateBox(double *x0);
        void showBox();
	double xLowEdge(int i);
	double xHighEdge(int i);
	double getDensity(){return density;}
};

class Particle
{
     private:
	double ptotal;
	double mass;
	double x[3] ;
	double p[3] ;
     public:
	Particle(double *x0, double *p0, double mass);
        void showParticle();
        void stepParticle(double delta);
	int checkBox(Box thisbox);
	double getRxy2(){return x[0]*x[0]+x[1]*x[1];}
};

int main()
{
// define a box object
     double xmin = -10.;
     double ymin = 20.;
     double zmin = -20.;
     double xmax = 15.0;
     double ymax = 30.0;
     double zmax = 50.0;
	double x0[3] ;
	x0[0] = xmax - xmin;
	x0[1] = ymax - ymin;
	x0[2] = zmax - zmin;
     	double density = 2.70;
	Box mybox(0.5*x0[0], 0.5*x0[1], 0.5*x0[2], density);
	double xc[3] ;
	xc[0] = xmin + 0.5*x0[0];
	xc[1] = ymin + 0.5*x0[1];
	xc[2] = zmin + 0.5*x0[2];
        mybox.locateBox(xc);
        mybox.showBox();
     	double volume = x0[0] * x0[1] * x0[2];
        double weight_expected = density * volume;
        cout << "Monte Carlo program to calculate weight"<< endl;
        cout << "weight expected = " << weight_expected << " gram" << endl;
       
//----- start Monte Carlo----------
// one can change these parameters, make sure your box is inside the
// scanned eta-phi region.
     const double eta_min = -5.0;
     const double eta_max = 5.0;
     const double phi_min = 0.0;
     const double phi_max = twopi;
     cout << "eta_min, eta_max = " << eta_min << " , "<< eta_max<<endl;
     cout << "phi_min, phi_max = " << phi_min << " , "<< phi_max<<endl;
     x0[0] = 0.0;  x0[1]=0.0 ; x0[2]=0.;
     double pParticle = 1000.0;
     double p0[3];
     double mass = 0.135;
// 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));
           p0[0] = pParticle * sin(theta) * cos(phi);
           p0[1] = pParticle * sin(theta) * sin(phi);
           p0[2] = pParticle * cos(theta);
           Particle pA(x0, p0, mass);
// step along the line from 0.1 cm to 100 cm
           for(int j=0; j<1000 ; j++)
	   {
              pA.stepParticle(step);
//            pA.showParticle();
// check if the step point is now inside the box
              if(( pA.checkBox(mybox) !=0 ) && (pA.checkBox(mybox) != 0)
                  && (pA.checkBox(mybox) != 0))
              total_weight = total_weight + mybox.getDensity()*
                             pA.getRxy2()*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;
}

Box::Box(double xh, double yh, double zh, double density0)
{
	xHalfLength[0] = xh;
	xHalfLength[1] = yh;
	xHalfLength[2] = zh;
	density = density0;
}  

void Box::locateBox(double *xc) 
{
	for (int i=0; i<3; i++)
	{   xcenter[i] = xc[i];
        }
}  

double Box::xLowEdge(int i) 
{
   return  xcenter[i] - xHalfLength[i];
}  

double Box::xHighEdge(int i) 
{
  return  xcenter[i] + xHalfLength[i];
}  

void Box::showBox()
{
	cout << "Box: (xh,yh,zh) =(" << xHalfLength[0] <<" , " << 
	        xHalfLength[1] <<" , "<< xHalfLength[2] <<")"<<endl;
	cout << "Box: (xc,yc,zc) =(" << xcenter[0] <<" , " << 
	        xcenter[1] <<" , "<< xcenter[2] <<")"<<endl;
}

Particle::Particle(double *x0, double *p0, double m)
{
	for (int i=0; i<3; i++)
	{   x[i] = x0[i];
	    p[i] = p0[i];
        }
	mass = m;
        ptotal = sqrt(p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
}

void Particle::showParticle()
{
    cout << "x=" << x[0] << " , "<< x[1] << " , " << x[2] << endl;
    cout << "p=" << p[0] << " , "<< p[1] << " , " << p[2] << 
         " ptotal, mass=" << ptotal << "," << mass << endl;
}

void Particle::stepParticle(double delta)
{
	for (int i=0; i<3; i++)
	 { x[i] += delta*p[i]/ptotal;}
}

int Particle::checkBox(Box thisbox)
{
	for (int i=0; i<3; i++)
	{
	   if ( (x[i] < thisbox.xLowEdge(i) ) || (thisbox.xHighEdge(i) < x[i]))
 		return 0;
	   else;
	}
	return 1;
}