#include "SiCapCal.h"


ClassImp(SiCapCal)

#define NRANSI
#define TINY 1.0e-20;
 void SiCapCal::ludcmp(float **a, int n, int *indx, float *d)
{
	int i,imax,j,k;
	float big,dum,sum,temp;
	float *vv;
	
	vv=vector(1,n);
	*d=1.0;
	for (i=1;i<=n;i++) {
		big=0.0;
		for (j=1;j<=n;j++)
			if ((temp=TMath::Abs(a[i][j])) > big) big=temp;
		if (big == 0.0) nrerror("Singular matrix in routine ludcmp");
		vv[i]=1.0/big;
	}
	for (j=1;j<=n;j++) {
		for (i=1;i<j;i++) {
			sum=a[i][j];
			for (k=1;k<i;k++) sum -= a[i][k]*a[k][j];
			a[i][j]=sum;
		}
		big=0.0;
		for (i=j;i<=n;i++) {
			sum=a[i][j];
			for (k=1;k<j;k++)
				sum -= a[i][k]*a[k][j];
			a[i][j]=sum;
			if ( (dum=vv[i]*TMath::Abs(sum)) >= big) {
				big=dum;
				imax=i;
			}
		}
		if (j != imax) {
			for (k=1;k<=n;k++) {
				dum=a[imax][k];
				a[imax][k]=a[j][k];
				a[j][k]=dum;
			}
			*d = -(*d);
			vv[imax]=vv[j];
		}
		indx[j]=imax;
		if (a[j][j] == 0.0) a[j][j]=TINY;
		if (j != n) {
			dum=1.0/(a[j][j]);
			for (i=j+1;i<=n;i++) a[i][j] *= dum;
		}
	}
	free_vector(vv,1,n);
}
#undef TINY
#undef NRANSI

 void SiCapCal::lubksb(float **a, int n, int *indx, float b[])
{
	int i,ii=0,ip,j;
	float sum;

	for (i=1;i<=n;i++) {
		ip=indx[i];
		sum=b[ip];
		b[ip]=b[i];
		if (ii)
			for (j=ii;j<=i-1;j++) sum -= a[i][j]*b[j];
		else if (sum) ii=i;
		b[i]=sum;
	}
	for (i=n;i>=1;i--) {
		sum=b[i];
		for (j=i+1;j<=n;j++) sum -= a[i][j]*b[j];
		b[i]=sum/a[i][i];
	}
}

 SiCapCal::SiCapCal(Int_t x1, Int_t x2, Double_t x3, Double_t x4,Double_t x5,Double_t x6,Int_t x7,Int_t x8)
{
Nread=x1; 
Niter=x2;
Pitch=x3;
Width=x4;
Thick=x5;
Lenght=x6;
Nwir=x7;
Sigstr=x8;
 ReturnCap=0;
}

 Double_t *SiCapCal::CapCal()
{
Double_t *x=CapCal(Nread,Niter,Pitch,Width,Thick,Lenght,Nwir,Sigstr);
return x;
}


 Double_t *SiCapCal::CapCal(int nread,int ninter,double pitch,double width,double thick,double lenght, int nwir, int sigstr)
{
    //     nread 	//!  no. of readout strips 
    //     ninter 	//!  no. of interpolation strips  
    //     pitch        //!  readout pitch
    //     width        //!  strip width
    //     thick        //!  thickness of detector
    //     lenght       //!  strip lenght
    //     sigstr       //!  strip with signal
    //       nwir       //! number of wires used in the approximation
  int nstrips;  // total number of strips
  int dim;      // dimension of arrays
  int ii, kk, ir, is, iw, nr, ns, nw, ifail; // indexes in loops
  int nwzero, maxwir, line, col, hitstr; 
  double wird;     //distance between wires
  double thick2;   //thiskness squared
  double eps;      //some phy. constants multiplied
 

  double rwir; 
  double *q,dist, arg, capt, capi, capb, rd;
  int *indx;
  Float_t totsum=0;
  Float_t **a,*d,*b,*v;



    eps=1/(3.14159*(1.+11.7)*8.854);  //!  1/(pi*(1+eps[Si])*eps0[pF/m])
    nstrips=(nread-1)*(ninter+1) + 1; //total number of strips
    dim=nstrips*nwir+(nread-1)*ninter; //!  dimension of arrays 
    //    Creation of an array
    if(dim>2000) {printf("Too many interpolation strips - not enough memory.n"); return(NULL);} 

    indx=new int [dim+1];
    q=new double [dim+1];
    v=new float [dim+1];
    d=new float [dim+1];
    b=new float [dim+1];
    a=matrix(0,dim+1,0,dim+1); 
    

    pitch = pitch/(ninter+1); // ! Strip pitch
    if(pitch<width) {printf("Pitch should be greater than width.n"); return(NULL);}

    rd=1.0;
    
    if(sigstr<0 || sigstr>ninter)  {printf("Can not set signal strip!!n"); return(NULL);}


  

    wird=width/nwir; //	    !  dist. betwen wires 	
    nwzero=int(nread/2)*nwir;
    thick2=thick*thick;
    rwir=rd*wird/2;
  

  // Print stage debug.
    printf("number of wires/strip=%dn",nwir);
    printf("precision of calculation=%en",eps);
    printf("number of strips=%dn",nstrips);
    printf("dimension of array=%dn",dim);
    printf("Readout pitch=%fn",pitch);
    printf("Interpolation strip with signal=%dn",sigstr);
    printf("Distance between the wires=%fn",wird);
    printf("Thickness=%fn",thick);
    
  
    for(ii=0;ii<dim;ii++) v[ii+1]=0;
    if(sigstr==0) 
      {
	nw=int(nread/2)*(ninter+1)*nwir;
	for(iw=nw; iw<nw+nwir; iw++) v[iw+1]=1;
      }
    else
	v[nstrips*nwir+int(nread/2.-1)*ninter+sigstr]=1;
       
    //    for(ii=0;ii<dim;ii++) printf("v[%d]=%fn",ii+1,v[ii+1]); // Control print
 
    for(ns=0;ns<nstrips;ns++)
       for(nw=0;nw<nwir;nw++)
	  for(is=0;is<nstrips;is++)
	     for(iw=0;iw<nwir;iw++)
	       if((iw==nw) && (is==ns))
		 a[nwir*ns+nw+1][nwir*is+iw+1] = -eps*TMath::Log(rwir/(2*thick));
	       else
		 {
	 	     dist=TMath::Abs((ns-is)*pitch + (nw-iw)*wird);
        	     arg=dist/TMath::Sqrt(dist*dist + 4*thick2);
		     a[nwir*ns+nw+1][nwir*is+iw+1] = -eps*TMath::Log(arg);
		 }
 
    

 for(nr=0;nr<nread-1;nr++)
       for(ns=0;ns<ninter;ns++)
	 {
	   line = nstrips*nwir+nr*ninter+ns;
	    for(iw=0;iw<dim;iw++)
	      if (((nr*(ninter+1)+ns+1)*nwir-1<iw) && ((nr*(ninter+1)+ns+2)*nwir-1>=iw)) 
		a[line+1][iw+1] = 1.; else a[line+1][iw+1] = 0.;	    
	 }

 for(nr=0;nr<nread-1;nr++)
       for(ns=0;ns<ninter;ns++)
	 {
	   col = nstrips*nwir+nr*ninter+ns;
	    for(iw=0;iw<nstrips*nwir;iw++)
	      if (((nr*(ninter+1)+ns+1)*nwir-1<iw) && ((nr*(ninter+1)+ns+2)*nwir-1>=iw))
		a[iw+1][col+1]= 1.; else a[iw+1][col+1] = 0.;	    
	 }


 //   for(ii=0;ii<dim;ii++) { printf("n"); for(kk=0;kk<dim;kk++); printf("%4.2e,",a[ii+1][kk+1]); } // Cont. Print
   printf("nCalculation startn");

      ludcmp(a,dim,indx,d);  
      lubksb(a,dim,indx,v);

  //   for(ii=0;ii<dim;ii++) { printf("n"); for(kk=0;kk<dim;kk++); printf("%4.2e,",a[ii+1][kk+1]); } // Cont. Print
 


  // for(ii=0;ii<dim;ii++) printf("b[%d]=%fn",ii,v[ii+1]); // Cont. Print

 if (sigstr==0)
   {
     capt = 0.;
     capi = 0.;
     hitstr=int(nread/2)*(ninter+1)+1;
     nw=(hitstr-1)*nwir;
     for(iw=nw;iw<nw+nwir;iw++) capt=capt+v[iw+1]; //  ! total readout strip capacitance
     printf("n");
      for(ns=0;ns<nstrips;ns+=ninter+1)
	{
	q[ns]=0;
	for(nw=0;nw<nwir;nw++) q[ns]=q[ns]+v[ns*nwir+nw+1];
	if (ns+1!=hitstr) capi = capi+q[ns];
	}
      capt = capt*lenght/100.;
      capi = -capi*lenght/100.;
      capb = capt - capi;
      printf("n Interstrip capacitance = %9.4f  pF",capi);
      printf("n Backplane capacitance of one strip = %8.3f  pF",capb);
      printf("n Total capacitance = %9.4f  pFn",capt);
   }
 else
   {
 for(ns=0;ns<nstrips;ns+=ninter+1)
	{
	q[ns]=0;
	for(nw=0;nw<nwir;nw++) q[ns]=q[ns]+v[ns*nwir+nw+1]; 
	printf("n Strip: %3.0d , Signal: %7.3f",ns+1,q[ns]);
	totsum+=q[ns];
	}
  printf("nTotal charge induced = %fn",totsum);
   }


//   TH1F *his=new TH1F("chsh","chsh",nstrips,0,nstrips);
//   his->SetMaximum(1);
//   for(ns=0;ns<nstrips;ns+=ninter+1) his->SetBinContent(ns,TMath::Abs(q[ns]));
//   his->SetFillColor(1);
//   his->DrawCopy(); his->Reset(); his->SetFillColor(2);
//   his->SetBinContent((int)((nread/2.-1)*(ninter+1))+sigstr,1);
//   his->DrawCopy("SAME");
//   // free_matrix(a, dim+1,  dim+1,  dim+1,  dim+1);
 delete [] indx;
 delete [] v;
 delete [] d;
 delete [] b;
 
 if(ReturnCap==1) {q[0]=capt; q[1]=capb; q[2]=capi;}
 return q;
}