// Simulation of Nuclear decays as a function of time
// Probability of decay for each nucleus per unit time is
//		P(t)*dt = 1/tau  * dt
// Between time t and t+dt, go to each nucleus and check if it decays
//
void NuclearDecay(){
	const int Nsteps = 1000;
	double tau=1.0,dt=0.001, t[Nsteps], N[Nsteps], dp;
	int N0,N1;
	TRandom2 r;

	N0 = N1 = 1000;
	dp = dt/tau;

	for(int i=0;i<Nsteps;i++){	// Loop over time
		for(int j=0;j<N0;j++){	// for this time step loop over each nucleus
			if(r.Uniform() < dp) N1--;	// Check if it decays, reduce count
		}
		N0 = N1;
		t[i] = i*dt;
		N[i] = N1;
		cout<< i*dt << "   "<<N1<<endl;
	}
	TGraph *gr  = new TGraph(Nsteps,t,N);
	gr->SetTitle("N(t) vs t");
	gr->GetXaxis()->SetTitle("time");
	gr->GetYaxis()->SetTitle("N(t)");
	gr->Draw("ACP");
}