#define NP 10
#define ND 20

#define DIAMETRE 5 

// Timer tm;

Vector lesx(ND);

/* int spos[NP] = {1,2,3,5,6,7,10,12,14,16};
for(int i=0; i<NP; i++)  lesx(i) = -spos[NP-i-1];  
for(int i=NP; i<2*NP; i++)  lesx(i) = spos[i-NP];  
*/
int spos[ND] = {-15,-13,-12,-10,-7,-6,-5,-3,-2,-1,1,2,3,5,6,7,10,12,14,15};
for(int i=0; i<ND; i++)  lesx(i) = spos[i];  


Matrix pos(2,ND*ND); 

int k = 0;
for(int i=0; i<ND; i++)  {
  for(int j=0; j<ND; j++) {
  	pos(0,k) = lesx(i); pos(1,k) = lesx(j);  k++;
  }
} 


Histo h(-0.5, 31.5, 32);

for(int i=0; i<ND-1; i++) 
  for(int j=i+1; j<ND; j++) 
    h.Add(fabs(lesx(i)-lesx(j)),1.);

KeepObj(h);
KeepObj(lesx);


Histo2D h2(-0.5, 31.5, 32, -0.5, 31.5, 32);
for(int i=0; i<ND*ND-1; i++) 
  for(int j=i+1; j<ND*ND; j++) 
    h2.Add(fabs(pos(0,i)-pos(0,j)),fabs(pos(1,i)-pos(1,j)), 1.);

KeepObj(h2);


Histo2D h2f(-0.5, 31.5, 96, -0.5, 31.5, 96);
for(int i=0; i<ND*ND-1; i++) 
  for(int j=i+1; j<ND*ND; j++) {
    FillVis(h2f, fabs(pos(0,i)-pos(0,j)),fabs(pos(1,i)-pos(1,j)));
  }
KeepObj(h2f);

cout << "h2f OK" << endl;

//---------------------------------------------------------------------
//----- Disposition compacte d'antennes circulaires -------------------
//---------------------------------------------------------------------
{
double DIAM_FA = 5.;
double L_FA = 100.;
double DISTA_FA = DIAM_FA;

int ND1_FA = L_FA/DISTA_FA;
int ND_FA = ND1_FA*ND1_FA;

Matrix pos_fa(2,ND_FA);
int kk=0;
for(int i=0; i<ND1_FA; i++) {
  for(int j=0; j<ND1_FA; j++) {
  	pos_fa(0, kk) = i*DISTA_FA;  pos_fa(1, kk) = j*DISTA_FA;  kk++;
  }
}
KeepObj(pos_fa);
cout << "pos_fa OK" << endl;

Histo2D h2_fa(0., 2*ND1_FA*DIAM_FA, 3*ND1_FA, 0., 2*ND1_FA*DIAM_FA, 3*ND1_FA);
Histo2D h2ns_fa(0., 2*ND1_FA*DIAM_FA, 4*ND1_FA, 0., 2*ND1_FA*DIAM_FA, 4*ND1_FA);

for(int i=0; i<pos_fa.NCols()-1; i++) 
  for(int j=i+1; j<pos_fa.NCols(); j++) {
  	h2ns_fa.Add(fabs(pos_fa(0,i)-pos_fa(0,j)),fabs(pos_fa(1,i)-pos_fa(1,j)),1.);
    FillVis(h2_fa, fabs(pos_fa(0,i)-pos_fa(0,j)),fabs(pos_fa(1,i)-pos_fa(1,j)), DIAM_FA);
  }
KeepObj(h2_fa);
KeepObj(h2ns_fa);
cout << "h2_fa OK" << endl;

}
//---------------------------------------------------------------------
//------------------ Disposition hexagonale ---------------------------
//---------------------------------------------------------------------
{
int nlayer = 4;
int nhex = 1;
for(int i=1; i<=nlayer; i++) {
  nhex += i*6;	
}
int ndet = nhex*7;
cout << " NLayer= " << nlayer << " NHex= " << nhex << " NDet= " << ndet << endl;
 
Matrix posh(2,ndet);

char * names[3] = {"xr", "yr", "sz"};
NTuple recp(3, names);
double xnt[5];
xnt[0] = xnt[1] = xnt[2] = xnt[3] = xnt[4] = 0.;

int kk = 0;
for(int i=0; i<= nlayer; i++) {
  double rlayer = i*4.5*DIAMETRE;
//  if (i == 1) rlayer = i*4.*DIAMETRE;
  Matrix ph = couche(i,rlayer) ; 
  cout << " DEBUG - I=" << i << " / NLayer= " << nlayer << " ph.NCols= " << ph.NCols() << endl;
  for(int j=0; j<ph.NCols(); j++) {
  	Matrix hex = hexagon(ph(0,j), ph(1,j), DIAMETRE, ph(2,i));
  	  for(int k=0; k<hex.NCols(); k++) {
  	  	xnt[0] = posh(0, kk) = hex(0, k); 
  		xnt[1] = posh(1, kk) = hex(1, k); 
  		recp.Fill(xnt);
	    kk++;
  	  }  
  }	
}

KeepObj(posh);
KeepObj(recp);

cout << "posh/recp OK" << endl;

Histo2D h2hex(0., 40.*DIAMETRE, 120, 0., 40.*DIAMETRE, 120);
Histo2D h2nshex(0., 20.*DIAMETRE, 120, 0., 20.*DIAMETRE, 120);

for(int i=0; i<posh.NCols()-1; i++) 
  for(int j=i+1; j<posh.NCols(); j++) {
  	h2nshex.Add(fabs(posh(0,i)-posh(0,j)),fabs(posh(1,i)-posh(1,j)), 1.);
    FillVis(h2hex, fabs(posh(0,i)-posh(0,j)),fabs(posh(1,i)-posh(1,j)), DIAMETRE);
    FillVis(h2hex, fabs(posh(0,i)-posh(0,j))*0.75,fabs(posh(1,i)-posh(1,j)), DIAMETRE);
    FillVis(h2hex, fabs(posh(0,i)-posh(0,j)),fabs(posh(1,i)-posh(1,j))*0.75, DIAMETRE);
  }
KeepObj(h2hex);
KeepObj(h2nshex);
cout << "h2hex OK" << endl;
}