inline Float_t dotprodZ(Int_t icha, Int_t ichb){return Rx[icha].x*Rx[ichb].x + Rx[icha].y*Rx[ichb].y + Rx[icha].z*Rx[ichb].z;}

Float_t dotprod(Float_t* a, Float_t* b, int length){
  Float_t sum=0;
  for (int cnt=0; cnt<length; cnt++)
    sum+=(*(a+cnt))*(*(b+cnt));	
  return sum;
}

const Int_t MSOLN=2048;
void vertex4Z(Float_t refrac,Float_t drxmin){
  //int numhits, double refrac, float rlvtx[3], double *a, double *t, int *j, float *volt, double *vrms, double& t01, double& t02, double& confit){
  //#include <fstream.h>
  
#define BETA	1
  ////////////////////////////////////////////////////
  //Jeffs Test Calculation
  ////////////////////////////////////////////////////
  //  const Int_t MCH=20;
  vector<Float_t>phiChZ4V1[MCH];
  vector<Float_t>phiChZ4V2[MCH];
  vector<Float_t>theChZ4V1[MCH];
  vector<Float_t>theChZ4V2[MCH];

  Float_t t01[MSOLN],rV1[MSOLN],phiV1[MSOLN],theV1[MSOLN],xV1[MSOLN],yV1[MSOLN],zV1[MSOLN];
  Float_t t02[MSOLN],rV2[MSOLN],phiV2[MSOLN],theV2[MSOLN],xV2[MSOLN],yV2[MSOLN],zV2[MSOLN];
  TH1F *hVphi=new TH1F("hVphi","hVphi",180,-180,180);
  TH1F *hVthe=new TH1F("hVthe","hVthe",90,0,180);
  Float_t xbar, ybar, zbar, phiAvgCh, theAvgCh;
  Bool_t LUseChZ4Vtx[MCH][MSOLN]={{false}};
  Int_t isoln=0;
  Int_t nsoln=0;
  Int_t nsolnt0=0;
  Int_t nsolnrV=0;
  Int_t nCombCh[MCH]={0};
  Float_t t01Ch[MCH]={0};
  Float_t t02Ch[MCH]={0};
  Float_t jx12,jx13,jx14,jy12,jy13,jy14,jz12,jz13,jz14,jt12,jt13,jt14,jt1sq,jt2sq,jt3sq,jt4sq,ja1sq,ja2sq,ja3sq,ja4sq,jdt12,jdt13,jdt14;
  Float_t jvtx1[3],jvtx2[3];
  Float_t jvec[3],jdt[3],jtdif[3],jmatrix[3][3],jvecmult[3];
  Float_t jray[3],jtev1vtx[4],jtev2vtx[4];
  Float_t ja,jb,jc;
  Float_t ja12[3],ja13[3],ja14[3];
  Float_t C=TMath::C()/refrac;

  Float_t v1vtx[3];
  Float_t v2vtx[3];
  Float_t vxvtx;
  Float_t vyvtx;
  Float_t vzvtx;
  //  Float_t a1z[3],a2z[3],a3z[3],a1zmag,a2zmag,a3zmag,amatrx[3][3],bvec[3],cx1mod,cx2mod,cx3mod,cx4mod;  
  Float_t antloc[4][3];
  Float_t times[4];
  Float_t volts[4];
  Float_t rms[4];
  
  Float_t VMOD(Float_t vmA[], int vmLength);
  void VSUB(Float_t vA[], Float_t vB[], Float_t* vC, int vLength);
  void QUAD(Float_t a, Float_t b, Float_t c, Float_t& xq1, Float_t& xq2);
  
  void matrxvecmult(Float_t* matrx, Float_t* vect, Float_t* dest);
  Int_t linsystem(Float_t* matrx, Float_t* vect);
  Int_t invert(Float_t* matrx);
  Float_t dotprod(Float_t* a, Float_t* b, int length);
  //  Float_t CTHETC=1/(BETA*refrac);
  
  Vvtx4Z1.t0=-NNN;
  Vvtx4Z1.r=-NNN;
  Vvtx4Z1.phi=-NNN;
  Vvtx4Z1.the=-NNN;
  Vvtx4Z1.dphi=-NNN;
  Vvtx4Z1.dthe=-NNN;

  Vvtx4Z2.t0=-NNN;
  Vvtx4Z2.r=-NNN;
  Vvtx4Z2.phi=-NNN;
  Vvtx4Z2.the=-NNN;
  Vvtx4Z2.dphi=-NNN;
  Vvtx4Z2.dthe=-NNN; 
  Int_t nfaildrmin=0;
  Int_t nfailNDim=0;
  Int_t ichi=0; 
  Int_t ichf=7;
  Int_t iRxCodeV;
  if(iPolVtx==-1 && false){
    ichi=0;
    ichf=15;
  }
  if(iPolVtx==1){
    ichi=8;
    ichf=15;
  }

  Float_t chisqVtxZ1, chisqVtxZ2, AchisqVtxZ1, AchisqVtxZ2, nHitCode;
  
  for(Int_t ich0=ichi; ich0<ichf; ich0++){
    if(!Rx[ich0].LVtx)continue;
    antloc[0][0]=Rx[ich0].x;   antloc[0][1]=Rx[ich0].y;   antloc[0][2]=Rx[ich0].z; times[0]=Rx[ich0].HitTime*1e-9; volts[0]=Rx[ich0].VMax; rms[0]=Rx[ich0].Vrms;
    for(Int_t ich1=ich0+1; ich1<ichf; ich1++){
      if(!Rx[ich1].LVtx)continue;
      antloc[1][0]=Rx[ich1].x;   antloc[1][1]=Rx[ich1].y;   antloc[1][2]=Rx[ich1].z; times[1]=Rx[ich1].HitTime*1e-9; volts[1]=Rx[ich1].VMax; rms[1]=Rx[ich1].Vrms;
      for(Int_t ich2=ich1+1; ich2<ichf; ich2++){
	if(!Rx[ich2].LVtx)continue;
	antloc[2][0]=Rx[ich2].x;   antloc[2][1]=Rx[ich2].y;   antloc[2][2]=Rx[ich2].z; times[2]=Rx[ich2].HitTime*1e-9; volts[2]=Rx[ich2].VMax; rms[2]=Rx[ich2].Vrms;
	for(Int_t ich3=ich2+1; ich3<ichf; ich3++){	
	  if(!Rx[ich3].LVtx)continue;
	  if(!LnRxDimXYZ(ich0,ich1,ich2,ich3)){
	    nfailNDim++;
	    continue;
	  }
	  if(!LdrRx(ich0,ich1,ich2,ich3,drxmin)){
	    nfaildrmin++;
	    continue; //compromised
	  }
	  if(!LMC && (iPolVtx==0 && !LParPairsEvV[ich0][ich1][ich2][ich3]))continue; // && !LParPairsEvH[ich0][ich1][ich2][ich3])continue;
	  if(isoln+1>=MSOLN)break;

	  antloc[3][0]=Rx[ich3].x;   antloc[3][1]=Rx[ich3].y;   antloc[3][2]=Rx[ich3].z; times[3]=Rx[ich3].HitTime*1e-9; volts[3]=Rx[ich3].VMax; rms[3]=Rx[ich3].Vrms;
	  jx12=Rx[ich0].x-Rx[ich1].x; //antloc[0][0]-antloc[1][0];
	  jy12=Rx[ich0].y-Rx[ich1].y; //antloc[0][1]-antloc[1][1];
	  jz12=Rx[ich0].z-Rx[ich1].z; //antloc[0][2]-antloc[1][2];
	  jx13=Rx[ich0].x-Rx[ich2].x; //antloc[0][0]-antloc[2][0];
	  jy13=Rx[ich0].y-Rx[ich2].y; //antloc[0][1]-antloc[2][1];
	  jz13=Rx[ich0].z-Rx[ich2].z; //antloc[0][2]-antloc[2][2];
	  jx14=Rx[ich0].x-Rx[ich3].x; //antloc[0][0]-antloc[3][0];
	  jy14=Rx[ich0].y-Rx[ich3].y; //antloc[0][1]-antloc[3][1];
	  jz14=Rx[ich0].z-Rx[ich3].z; //antloc[0][2]-antloc[3][2];
	  jt12=times[0]-times[1];
	  jt13=times[0]-times[2];
	  jt14=times[0]-times[3];
	  ja1sq=dotprodZ(ich0,ich0); //&antloc[0][0],&antloc[0][0],3);
	  ja2sq=dotprodZ(ich1,ich1); //(&antloc[1][0],&antloc[1][0],3);
	  ja3sq=dotprodZ(ich2,ich2); //(&antloc[2][0],&antloc[2][0],3);
	  ja4sq=dotprodZ(ich3,ich3); //(&antloc[3][0],&antloc[3][0],3);
	  jt1sq=C*C*times[0]*times[0];
	  jt2sq=C*C*times[1]*times[1];
	  jt3sq=C*C*times[2]*times[2];
	  jt4sq=C*C*times[3]*times[3];
	  jdt12=(ja1sq-ja2sq)-(jt1sq-jt2sq);
	  jdt13=(ja1sq-ja3sq)-(jt1sq-jt3sq);
	  jdt14=(ja1sq-ja4sq)-(jt1sq-jt4sq);
	  ja12[0]=jx12;
	  ja12[1]=jy12;
	  ja12[2]=jz12;
	  ja13[0]=jx13;
	  ja13[1]=jy13;
	  ja13[2]=jz13;
	  ja14[0]=jx14;
	  ja14[1]=jy14;
	  ja14[2]=jz14;  
  
	  jmatrix[0][0]=2*jx12;
	  jmatrix[0][1]=2*jy12;
	  jmatrix[0][2]=2*jz12;
	  jmatrix[1][0]=2*jx13;
	  jmatrix[1][1]=2*jy13;
	  jmatrix[1][2]=2*jz13;
	  jmatrix[2][0]=2*jx14;
	  jmatrix[2][1]=2*jy14;
	  jmatrix[2][2]=2*jz14;
  
	  jdt[0]=jdt12;
	  jdt[1]=jdt13;
	  jdt[2]=jdt14;
	  
	  jtdif[0]=jt12;
	  jtdif[1]=jt13;
	  jtdif[2]=jt14;
	  
	  jvec[0]=jdt[0];
	  jvec[1]=jdt[1];
	  jvec[2]=jdt[2];
	  
	  linsystem(&jmatrix[0][0],&jvec[0]);
	  
	  jvecmult[0]=jmatrix[0][0]*jtdif[0]+jmatrix[0][1]*jtdif[1]+jmatrix[0][2]*jtdif[2];
	  jvecmult[1]=jmatrix[1][0]*jtdif[0]+jmatrix[1][1]*jtdif[1]+jmatrix[1][2]*jtdif[2];
	  jvecmult[2]=jmatrix[2][0]*jtdif[0]+jmatrix[2][1]*jtdif[1]+jmatrix[2][2]*jtdif[2];
  
	  jray[0]=2*C*C*jvecmult[0];
	  jray[1]=2*C*C*jvecmult[1];
	  jray[2]=2*C*C*jvecmult[2];
	  
	  ja=dotprod(&jray[0],&jray[0],3)-C*C;
	  jb=2*C*C*times[0]-2*dotprod(&jray[0],&antloc[0][0],3)+2*dotprod(&jray[0],&jvec[0],3);
	  jc=dotprod(&jvec[0],&jvec[0],3)-2*dotprod(&jvec[0],&antloc[0][0],3)+dotprod(&antloc[0][0],&antloc[0][0],3)-C*C*times[0]*times[0];
	  QUAD(ja,jb,jc,jtev1vtx[0],jtev2vtx[0]);
	  
	  ja=dotprod(&jray[0],&jray[0],3)-C*C;
	  jb=2*C*C*times[1]-2*dotprod(&jray[0],&antloc[1][0],3)+2*dotprod(&jray[0],&jvec[0],3);
	  jc=dotprod(&jvec[0],&jvec[0],3)-2*dotprod(&jvec[0],&antloc[1][0],3)+dotprod(&antloc[1][0],&antloc[1][0],3)-C*C*times[1]*times[1];
	  QUAD(ja,jb,jc,jtev1vtx[1],jtev2vtx[1]);
  
	  ja=dotprod(&jray[0],&jray[0],3)-C*C;
	  jb=2*C*C*times[2]-2*dotprod(&jray[0],&antloc[2][0],3)+2*dotprod(&jray[0],&jvec[0],3);
	  jc=dotprod(&jvec[0],&jvec[0],3)-2*dotprod(&jvec[0],&antloc[2][0],3)+dotprod(&antloc[2][0],&antloc[2][0],3)-C*C*times[2]*times[2];
	  QUAD(ja,jb,jc,jtev1vtx[2],jtev2vtx[2]);
  
	  ja=dotprod(&jray[0],&jray[0],3)-C*C;
	  jb=2*C*C*times[3]-2*dotprod(&jray[0],&antloc[3][0],3)+2*dotprod(&jray[0],&jvec[0],3);
	  jc=dotprod(&jvec[0],&jvec[0],3)-2*dotprod(&jvec[0],&antloc[3][0],3)+dotprod(&antloc[3][0],&antloc[3][0],3)-C*C*times[3]*times[3];
	  QUAD(ja,jb,jc,jtev1vtx[3],jtev2vtx[3]);
  
	  jvtx1[0]=jvec[0]+jtev1vtx[0]*jray[0];
	  jvtx1[1]=jvec[1]+jtev1vtx[0]*jray[1];
	  jvtx1[2]=jvec[2]+jtev1vtx[0]*jray[2];
	  jvtx2[0]=jvec[0]+jtev2vtx[0]*jray[0];
	  jvtx2[1]=jvec[1]+jtev2vtx[0]*jray[1];
	  jvtx2[2]=jvec[2]+jtev2vtx[0]*jray[2];

	  Float_t jraymag;
	  jraymag=sqrt(jray[0]*jray[0]+jray[1]*jray[1]+jray[2]*jray[2]);
	  VNorm2Avg(jvtx1[0],jvtx1[1],jvtx1[2],xV1[isoln],yV1[isoln],zV1[isoln]);
	  Cart2SpherDeg(jvtx1[0],jvtx1[1],jvtx1[2],rV1[isoln],phiV1[isoln],theV1[isoln]);
	  v1vtx[0]=jvtx1[0];
	  v1vtx[1]=jvtx1[1];
	  v1vtx[2]=jvtx1[2];

	  VNorm2Avg(jvtx2[0],jvtx2[1],jvtx2[2],xV2[isoln],yV2[isoln],zV2[isoln]);
	  Cart2SpherDeg(jvtx2[0],jvtx2[1],jvtx2[2],rV2[isoln],phiV2[isoln],theV2[isoln]);
	  v2vtx[0]=jvtx2[0];
	  v2vtx[1]=jvtx2[1];
	  v2vtx[2]=jvtx2[2];
	  vxvtx=v1vtx[0];
	  vyvtx=v1vtx[1];
	  vzvtx=v1vtx[2];
	  
	  t01[isoln]=jtev1vtx[0]*1e9;
	  t02[isoln]=jtev2vtx[0]*1e9;
	  t01Ch[ich0]+=t01[isoln];
	  t01Ch[ich1]+=t01[isoln];
	  t01Ch[ich2]+=t01[isoln];
	  t01Ch[ich3]+=t01[isoln];
	  t02Ch[ich0]+=t02[isoln];
	  t02Ch[ich1]+=t02[isoln];
	  t02Ch[ich2]+=t02[isoln];
	  t02Ch[ich3]+=t02[isoln];
	  nCombCh[ich0]++;
	  nCombCh[ich1]++;
	  nCombCh[ich2]++;
	  nCombCh[ich3]++;

	  if(fabs(t01[isoln]-t02[isoln])/(t01[isoln]+t02[isoln])>0.001){
	    phiV2[isoln]-=180; //noted empircally
	    xV2[isoln]=rV2[isoln]*cos(TMath::DegToRad()*phiV2[isoln])*sin(TMath::DegToRad()*theV2[isoln]);
	    yV2[isoln]=rV2[isoln]*sin(TMath::DegToRad()*phiV2[isoln])*sin(TMath::DegToRad()*theV2[isoln]);
	  }
	  phiV1[isoln]=Map2pm180deg(phiV1[isoln]);
	  phiV2[isoln]=Map2pm180deg(phiV2[isoln]);
	  if(LDebug)printf("vtx4Z Soln1: t0=%g/r=%g/phi=%g/the=%g // Soln2: t0=%g/r=%g/phi=%g/the=%g\n",t01[isoln],rV1[isoln],phiV1[isoln],theV1[isoln],t02[isoln],rV2[isoln],phiV2[isoln],theV2[isoln]);
	  nHitCode=ich0+100*ich1+10000*ich2+1000000*ich3;
	  iRxCodeV=pow(2,ich0)+pow(2,ich1)+pow(2,ich2)+pow(2,ich3);
	  
	  phiChZ4V1[ich0].push_back(phiV1[isoln]);	  phiChZ4V1[ich1].push_back(phiV1[isoln]);	  phiChZ4V1[ich2].push_back(phiV1[isoln]);	  phiChZ4V1[ich3].push_back(phiV1[isoln]);
	  theChZ4V1[ich0].push_back(theV1[isoln]);	  theChZ4V1[ich1].push_back(theV1[isoln]);	  theChZ4V1[ich2].push_back(theV1[isoln]);	  theChZ4V1[ich3].push_back(theV1[isoln]);
	  phiChZ4V2[ich0].push_back(phiV2[isoln]);	  phiChZ4V2[ich1].push_back(phiV2[isoln]);	  phiChZ4V2[ich2].push_back(phiV2[isoln]);	  phiChZ4V2[ich3].push_back(phiV2[isoln]);
	  theChZ4V2[ich0].push_back(theV2[isoln]);	  theChZ4V2[ich1].push_back(theV2[isoln]);	  theChZ4V2[ich2].push_back(theV2[isoln]);	  theChZ4V2[ich3].push_back(theV2[isoln]);

	  chisqVtxZ1=chisqVtxPt(ich0,ich1,ich2,ich3,999,phiV1[isoln],theV1[isoln],1.);
	  chisqVtxZ2=chisqVtxPt(ich0,ich1,ich2,ich3,999,phiV2[isoln],theV2[isoln],1.);
	  if(false)AchisqVtxZ1=AchisqVtxPt(ich0,ich1,ich2,ich3,phiV1[isoln],theV1[isoln]);
	  if(false)AchisqVtxZ2=AchisqVtxPt(ich0,ich1,ich2,ich3,phiV2[isoln],theV2[isoln]);
	  NFIT4ZSOLNS->Fill(nHitCode,iRxCodeV,isoln,t01[isoln],rV1[isoln],phiV1[isoln],theV1[isoln],t02[isoln],rV2[isoln],phiV2[isoln],theV2[isoln],chisqVtxZ1,chisqVtxZ2,AchisqVtxZ1,AchisqVtxZ2);
	  nsoln++;
	  if(t01[isoln]==0 || t02[isoln]==0)continue;
	  nsolnt0++;
	  if(rV1[isoln]>2048 || rV2[isoln]>2048)continue;
	  nsolnrV++;

	  if(Rx[ich0].x==Rx[ich1].x||Rx[ich0].x==Rx[ich2].x||Rx[ich0].x==Rx[ich3].x)continue;
	  if(Rx[ich0].y==Rx[ich1].y||Rx[ich0].y==Rx[ich2].y||Rx[ich0].y==Rx[ich3].y)continue;
	  if(Rx[ich0].z==Rx[ich1].z||Rx[ich0].z==Rx[ich2].z||Rx[ich0].z==Rx[ich3].z)continue;

	  if(Rx[ich1].x==Rx[ich2].x||Rx[ich1].x==Rx[ich3].x)continue;
	  if(Rx[ich1].y==Rx[ich2].y||Rx[ich1].y==Rx[ich3].y)continue;
	  if(Rx[ich1].z==Rx[ich2].z||Rx[ich1].z==Rx[ich3].z)continue;

	  if(Rx[ich2].x==Rx[ich3].x)continue;
	  if(Rx[ich2].y==Rx[ich3].y)continue;
	  if(Rx[ich2].z==Rx[ich3].z)continue;

	  hVphi->Fill(phiV1[isoln]);
	  hVphi->Fill(phiV2[isoln]);
	  hVthe->Fill(theV1[isoln]);
	  hVthe->Fill(theV2[isoln]);
	  if(ich3<8)h2ZphiVRxCode->Fill(phiV1[isoln],iRxCodeV);
	  LUseChZ4Vtx[ich0][isoln]=true; LUseChZ4Vtx[ich1][isoln]=true; LUseChZ4Vtx[ich2][isoln]=true; LUseChZ4Vtx[ich3][isoln]=true;
	  isoln++;
	}
      }
    }
  }
  if(isoln==0){
    delete hVphi;
    delete hVthe;
    return;
  }
  idum=hVphi->GetMaximumBin();
  Vvtx4Z1.phiTH1=hVphi->GetBinCenter(idum);
  Vvtx4Z1.dphiTH1=hVphi->GetRMS();
  hVphi->SetBinContent(idum,0);
  printf("vtx4Z max at %g/entries=%d /2nd max=%g/entries=%d nfailNDim=%d ifaildrmin=%d nsoln=%d/nsolnrt0=%d/nsolnV=%d/isoln=%d\n",Vvtx4Z1.phiTH1,idum,GetxMax(hVphi),int(hVphi->GetMaximum()),nfailNDim,nfaildrmin,nsoln,nsolnt0,nsolnrV,isoln);
  delete hVphi;

  idum=hVthe->GetMaximumBin();
  Vvtx4Z1.theTH1=hVthe->GetBinCenter(idum);
  Vvtx4Z1.dtheTH1=hVthe->GetRMS();
  
  delete hVthe;

  Vvtx4Z1.t0=TMath::Mean(isoln,t01);
  Vvtx4Z1.r=TMath::Mean(isoln,rV1);
  Vvtx4Z1.dphi=TMath::RMS(isoln,phiV1);
  Vvtx4Z1.dthe=TMath::RMS(isoln,theV1);
  xbar=TMath::Mean(isoln,xV1);
  ybar=TMath::Mean(isoln,yV1);
  zbar=TMath::Mean(isoln,zV1);
  Cart2SpherDeg(xbar,ybar,zbar,fdum,fdum1,Vvtx4Z1.the);
  Vvtx4Z1.phi=Map2pm180deg(fdum1);

  Vvtx4Z2.t0=TMath::Mean(isoln,t02);
  Vvtx4Z2.r=TMath::Mean(isoln,rV2);
  Vvtx4Z2.dphi=TMath::RMS(isoln,phiV2);
  Vvtx4Z2.dthe=TMath::RMS(isoln,theV2);
  xbar=TMath::Mean(isoln,xV2);
  ybar=TMath::Mean(isoln,yV2);
  zbar=TMath::Mean(isoln,zV2);
  Cart2SpherDeg(xbar,ybar,zbar,fdum,fdum1,Vvtx4Z2.the);
  Vvtx4Z2.phi=Map2pm180deg(fdum1);
  printf("vtx4Z <hVphi>=%g+/-%g / <hVthe>=%g+/-%g \n<Soln1> r=%g/phi=%g+/-%g/the=%g+/-%g/t0=%g\n<Soln2> r=%g/phi=%g+/-%g/the=%g+/-%g/t0=%g\n",Vvtx4Z1.phiTH1,Vvtx4Z1.dphiTH1,Vvtx4Z1.theTH1,Vvtx4Z1.dtheTH1,Vvtx4Z1.r,Vvtx4Z1.phi,Vvtx4Z1.dphi,Vvtx4Z1.the,Vvtx4Z1.dthe,Vvtx4Z1.t0,Vvtx4Z2.r,Vvtx4Z2.phi,Vvtx4Z2.dphi,Vvtx4Z2.the,Vvtx4Z2.dthe,Vvtx4Z2.t0);
  for(Int_t ich=0; ich<MCH; ich++){
    phiAvgCh=TMath::Mean(sizeof(phiChZ4V1[ich]),phiChZ4V1[ich]);
    theAvgCh=TMath::Mean(sizeof(theChZ4V1[ich]),theChZ4V1[ich]);
    if(nCombCh[ich]==0)continue;
    h2residZ4VCh->Fill(t01Ch[ich]/nCombCh[ich]-Vvtx4Z1.t0,ich);
    h2residZ4VCh->Fill(t02Ch[ich]/nCombCh[ich]-Vvtx4Z2.t0,ich);
      abort();
    for(Int_t jsoln=0; jsoln<isoln; jsoln++){
      if(!LUseChZ4Vtx[ich][isoln]&&jsoln<sizeof(phiChZ4V1[ich])){
	h2dphiZ41VCh->Fill(phiAvgCh-phiChZ4V1[ich][jsoln],ich);
	h2dtheZ41VCh->Fill(theAvgCh-theChZ4V1[ich][jsoln],ich);
	h2dphiZ42VCh->Fill(phiAvgCh-phiChZ4V2[ich][jsoln],ich);
	h2dtheZ42VCh->Fill(theAvgCh-theChZ4V2[ich][jsoln],ich);
      }
    }
  }
}
	
void QUAD(Float_t a, Float_t b, Float_t c, Float_t& xq1, Float_t& xq2){
  Float_t arg;
  xq1=xq2=0;
  arg=b*b-4*a*c;
  if (arg<0){
    arg=0;
  }
  xq1=(-b+sqrt(arg))/(2*a);
  xq2=(-b-sqrt(arg))/(2*a);
  if (a > 1E20 || b > 1E20 || c > 1E20){
    xq1=xq2=0;
  }
}


void VSUB(Float_t vA[3], Float_t vB[3], Float_t* vC, int vLength){
  for (int vCnt=0; vCnt<vLength; vCnt++)
    *(vC+vCnt)=*(vA+vCnt)-*(vB+vCnt);
}


Float_t VMOD(Float_t vmA[], int vmLength){
  Float_t vmSum=0;
  for (int vmCnt=0; vmCnt<vmLength; vmCnt++)
    vmSum+=vmA[vmCnt]*vmA[vmCnt];
  return sqrt(vmSum);
}


int linsystem(Float_t* matrx, Float_t* vect){
  Float_t destvec[3]={0,0,0};
  int determ;
  void matrxvecmult(Float_t* matrx, Float_t* vect, Float_t* dest);
  int invert(Float_t* matrx);
  
  determ=invert(matrx);
  if (determ==1){
    matrxvecmult(matrx, vect, &destvec[0]);
    *vect=destvec[0];
    *(vect+1)=destvec[1];
    *(vect+2)=destvec[2];
    return 1;
  }
  else
    return 0;
}

void matrxvecmult(Float_t* matrx, Float_t* vect, Float_t* dest){
  for (int cnt=0; cnt<3; cnt++)
    for (int cnt2=0; cnt2<3; cnt2++) 
      *(dest+cnt)+=*(matrx+3*cnt+cnt2)*(*(vect+cnt2));
}

int invert(Float_t* matrx){
  Float_t adj[3][3], adjt[3][3], det;
  adj[0][0]=(*(matrx+4)*(*(matrx+8)))-(*(matrx+5)*(*(matrx+7)));
  adj[0][1]=-((*(matrx+3)*(*(matrx+8)))-(*(matrx+5)*(*(matrx+6))));
  adj[0][2]=(*(matrx+3)*(*(matrx+7)))-(*(matrx+4)*(*(matrx+6)));
  adj[1][0]=-((*(matrx+1)*(*(matrx+8)))-(*(matrx+2)*(*(matrx+7))));
  adj[1][1]=(*(matrx)*(*(matrx+8)))-(*(matrx+2)*(*(matrx+6)));
  adj[1][2]=-((*(matrx)*(*(matrx+7)))-(*(matrx+1)*(*(matrx+6))));
  adj[2][0]=(*(matrx+1)*(*(matrx+5)))-(*(matrx+2)*(*(matrx+4)));
  adj[2][1]=-((*(matrx)*(*(matrx+5)))-(*(matrx+2)*(*(matrx+3))));
  adj[2][2]=(*(matrx)*(*(matrx+4)))-(*(matrx+1)*(*(matrx+3)));
  
  adjt[0][0]=adj[0][0];
  adjt[0][1]=adj[1][0];
  adjt[0][2]=adj[2][0];
  adjt[1][0]=adj[0][1];
  adjt[1][1]=adj[1][1];
  adjt[1][2]=adj[2][1];
  adjt[2][0]=adj[0][2];
  adjt[2][1]=adj[1][2];
  adjt[2][2]=adj[2][2];
  
  det=*(matrx)*adj[0][0]+*(matrx+1)*adj[0][1]+*(matrx+2)*adj[0][2];
  
  if (det!=0){
    for (int cnt=0; cnt<3; cnt++){
      for (int cnt2=0; cnt2<3; cnt2++){
	adjt[cnt][cnt2]*=1/det;
	*(matrx+3*cnt+cnt2)=adjt[cnt][cnt2];
      }
    }
    return 1;
  }
  else
    return 0;
}