void Coverage()
{
  double rnd1,rnd2;
  int Nthrows=50000;
  int Ntotal = 300;
  float p=0.15;
  int nblue;
  int ncov=0;
  int ncov2=0;
  for (int it=0;it!=Nthrows;it++) {
    nblue=0;
    for (int i=0;i!=Ntotal;i++)
      {      
	rnd1=rand()/2147483647.;
	//	cout << rnd1 << endl;
	if (rnd1<p) nblue++;
      }	
    double Pmeas = nblue/300.;
    //    cout << nblue << " p = " << Pmeas << " Error " << sqrt(Pmeas*(1-Pmeas)/300.) << endl;
    if (fabs(Pmeas-p)<=sqrt(Pmeas*(1-Pmeas)/300.))
      {
	//	cout << "Value is in!" << endl;
	ncov++;
      }
    if (fabs(Pmeas-p)<=sqrt(p*(1-p)/300.))
      {
	//	cout << "Value is in!" << endl;
	ncov2++;
      }
  }
  
  cout << "Fraction of throws not covered " << (float) ncov/Nthrows << "  " << sqrt(p*(1-p)/300./sqrt(2.*Nthrows))<< endl;
  cout << "Fraction of throws not covered by stat error " << (float) ncov2/Nthrows << endl;


}