macro matrixcorrect sys=4 part=pim ybn0=y0 ybin=b2 yin=b1 yout=b1

hi/del *
*vec/del *

*************************************************************************
*This kumac calculates the renormalized Transpose of the 
*TPC Response Matrix, based on the number of non-zero data points in
*the measured data.  The reason for this choice is to prevent
*edge effects in the proposed "corrected spectra".
*The renormalized transpose is then combined with a paramterized
*loss vector to fully correct for momentum skew and tracking losses.
*************************************************************************

hi/file 31 hbk/tpc-and-geo-accept_[sys]gev_[part]_[ybn0].hbk

*vec/cre par(4) r 1.0 0.5 .02 .08
hi/fit 6000(1:45) lossfit2.for ! 4 par
vec/cre loss0(40) r
hi/get_vect/func 6000 loss0
close 31
hi/del *

hi/file 30 hbk/tpc-and-geo-accept_[sys]gev_[part]_[ybin].hbk

*vec/cre par(4) r 1.0 0.5 .02 .08
hi/fit 6000(1:45) lossfit2.for ! 4 par
vec/cre loss(40) r
hi/get_vect/func 6000 loss

do i=1,40
  vec/op/vadd loss0([i]) loss([i]) loss([i])
  vec/op/vscale loss([I]) 0.5 loss([i])
enddo

vec/cre matrix(80,80) r
vec/cre rmtm0(80) r
vec/cre skewoutput(40) r
vec/cre lossoutput(40) r
vec/cre output(40) r

hi/get_vect/con 1000 matrix
hi/get_vect/con 55 rmtm0 

vec/cre vectocorrect(2,40) r
vec/read vectocorrect ~cebra/e895/specmt/[sys]gev_pip_[yin]_raw.dat

*Find out where the data stops
hidatabin=40
do i = 1,40
  data = vectocorrect(1,[i])
  if ([data].gt.0) then
   hidatabin = [i]
  endif
enddo

*Compute the renormalization factor using the hidatabin
vec/cre renorm([hidatabin]) r
do j = 1,[hidatabin]
 sum = 0
 do i = 1,[hidatabin]
   mat = matrix([i],[j])
   sum = [sum] + [mat]
 enddo
 vec/in renorm([j]) [sum]
enddo

vec/cre NormMatrix([hidatabin],[hidatabin]) r
*ReNormalize mtm0 matrix to output mtm0
do j=1,[hidatabin]
 do i=1,[hidatabin]
   renorm = renorm([i])
   if ([renorm].gt.0) then
     vec/op/vdivide matrix([i],[j]) renorm([i]) NormMatrix([i],[j])
   else
     vec/in NormMatrix([i],[j]) 0
   endif
 enddo
enddo

*Take Transpose of renormalizied mtm0 matrix:
do i = 1,[hidatabin]
  do j = 1,[hidatabin]
     matij = NormMatrix([i],[j])
     matji = NormMatrix([j],[i])
     vec/in NormMatrix([i],[j]) [matji]
     vec/in NormMatrix([j],[i]) [matij]
  enddo
enddo

do j = 1,[hidatabin]
   ov = 0
   do i = 1,[hidatabin]
      rmt = vectocorrect(1,[i])
      mat = NormMatrix([i],[j])
      ov = [ov] + [rmt]*[mat]
   enddo
   if ([ov].lt.0) then
      vec/in skewouput([j]) 0
   else
      vec/in skewoutput([j]) [ov]
   endif
   vec/op/vdivide vectocorrect(1,[j]) loss([j]) lossoutput([j])
   vec/op/vdivide skewoutput([j]) loss([j]) output([j])
enddo

1d 500 'Uncorrected-output' 40 0 1 
hi/put_vect/con 500 vectocorrect(1)
hi/put_vect/err 500 vectocorrect(2)

1d 501 'Skew-corrected-output' 40 0 1 
hi/put_vect/con 501 skewoutput
hi/put_vect/err 501 vectocorrect(2)

1d 502 'Loss-corrected-output' 40 0 1 
hi/put_vect/con 502 lossoutput
hi/put_vect/err 502 vectocorrect(2)

1d 503 'Skew-AND-Loss-corrected-output' 40 0 1 
hi/put_vect/con 503 output
hi/put_vect/err 503 vectocorrect(2)

1d 504 'Skew-correction-only' 40 0 1 
vec/cre skewonly(40) r
do i = 1,[hidatabin]
  vec/op/vdivide skewoutput([i]) vectocorrect(1,[i]) skewonly([i])
enddo
hi/put_vect/con 504 skewonly

1d 505 'Loss-correction-only' 40 0 1 
vec/cre lossCF(40) r
vec/cre unity(40) r 40*1.0
do i = 1,40
  vec/op/vdivide unity([i]) loss([i]) lossCF([i])
enddo
hi/put_vect/con 505 lossCF

1d 506 'Skew-AND-Loss-correction-combined' 40 0 1 
vec/cre CF(40) r
do i = 1,[hidatabin]
  vec/op/vdivide output([i]) vectocorrect(1,[i]) CF([i])
enddo
hi/put_vect/con 506 CF

2d 1001 'Renormalized-transpose-Correction-Matrix' [hidatabin] 0 $eval([hidatabin]*0.025) [hidatabin] 0 $eval([hidatabin]*0.025)
hi/put_vect/con 1001 NormMatrix

*Write NORMALIZED vectors to data files
*vec/cre spec(40) r
*vec/cre erro(40) r
*hi/get_vec/con 501 spec
*hi/get_vec/err 501 erro
*vec/write spec,erro dat/[sys]gev_[part]_[yout]_cskw.dat _
*'(40(F11.5,3x,F10.5))' 'oc'

vec/cre spec(40) r
vec/cre erro(40) r
hi/get_vec/con 502 spec
hi/get_vec/err 502 erro
vec/write spec,erro dat/[sys]gev_pip_[yout]_clss.dat _
'(40(F11.5,3x,F10.5))' 'oc'

*vec/cre spec(40) r
*vec/cre erro(40) r
*hi/get_vec/con 503 spec
*hi/get_vec/err 503 erro
*vec/write spec,erro dat/[sys]gev_[part]_[yout]_csal.dat _
*'(40(F11.5,3x,F10.5))' 'oc'


hi/file 29 hbk/[sys]gev_[part]_[yout]_corrected.hbk ! n

*mtm0 data:
hi/hrout 500
hi/hrout 501
hi/hrout 502
hi/hrout 503
hi/hrout 504
hi/hrout 505
hi/hrout 506
*mtm0 correction info:
hi/hrout 1000
hi/hrout 1001
hi/hrout 6000
close 29
close 30

return