Matrix hexagon(double xc, double yc, double r=1., double theta=0.) 
{
	Matrix hex(2,7);
	hex(0,0) = xc;  hex(1,0) = yc;
	for(int i=1; i<=6; i++) {
	  double tet = (i-1)*M_PI/3.+theta;
	  hex(0,i) = r*cos(tet)+xc;
	  hex(1,i) = r*sin(tet)+yc;
	}
	return hex;
}

Matrix couche(int num, double r=0.)
{
	if (num <= 0) {
      Matrix pos(3,1);
	  pos(0,0) = pos(1,0) = pos(2,0) = 0.;
	  return pos;
	}
	else {
	  if (r <= 0.0001)  r = 5.*num;  	
	  int n = num*6;
	  double tet0 = 0.;
	  tet0 = M_PI/n;
	  double dtet = M_PI*2./(double)n;
	  Matrix pos(3, n);
  	  for(int i=0; i<n; i++) {
	    double tet = i*dtet+tet0;
	    pos(0,i) = r*cos(tet);
	    pos(1,i) = r*sin(tet);
	    pos(2,i) = tet;
	  }
	  return pos;
	}
}

int FillVis(Histo2D& h2, double x, double y, double d=1) 
{
  h2.Add(x, y, 1.);

  double wa = 0.6;  
  double wb = 0.3;  
  for(int i=0; i<6; i++) {
  	double dx1a = 0.35*cos(i*M_PI/3.)*d;
  	double dy1a = 0.35*sin(i*M_PI/3.)*d;
  	double dx1b = 2*dx1a; 
  	double dy1b = 2*dy1a;
  	for(int j=0; j<6; j++) {
  	  double dx2a = 0.35*cos(j*M_PI/3.)*d;
  	  double dy2a = 0.35*sin(j*M_PI/3.)*d;
  	  double dx2b = 2*dx2a; 
  	  double dy2b = 2*dy2a;
  	  h2.Add((x+dx1a-dx2a), (y+dy1a-dy2a), wa); 
  	  h2.Add((x+dx1b-dx2b), (y+dy1b-dy2b), wb);  	
  	}
  }
  return 0;
}