*------------------------------------------------------------------------ * MATRIX find the reconstruction matrix * -======- * * finds the coefficients of the reconstruction matrix up to a given order * - reads the data from a PAW/hbook ntuple written by hms.f * - writtes the results (coeffs & errors) to a text file * * based on Chris Cothran's mcrecon * widely modified by Markus Muehlbauer *------------------------------------------------------------------------ IMPLICIT NONE INTEGER nwpawc,istat PARAMETER (nwpawc = 1000000) REAL paw(nwpawc) COMMON /pawc/paw ! commonblocks holding the event data and the reconstruction matrix INTEGER NTERMS,num PARAMETER (NTERMS = 300) INTEGER e(NTERMS,4) REAL c(NTERMS,4),err(NTERMS,4),chisq(4) COMMON /matrix/ c,err,e,num,chisq INTEGER NEVENTS,cnt PARAMETER (NEVENTS = 100000) REAL u(NEVENTS,6),v(NEVENTS,4) COMMON /events/ u,v,cnt REAL event(12) INTEGER i,j,minord,maxord,sym CHARACTER filename1*80,filename2*80,name*80 print *,'Matrix: calculates a reconstruction matrix' print *,'======= ' print * PRINT '(X,A,$)','Enter the name of the event file: ' READ (5,'(A)') filename1 PRINT '(X,A,$)','Enter the name of the matrix: ' READ (5,'(A)') filename2 PRINT '(X,A,$)','Enter the mimimal order of the matrix: ' READ (5,*) minord PRINT '(X,A,$)','Enter the maximal order of the matrix: ' READ (5,*) maxord PRINT '(X,A,$)','Assume midplane symmetry [Yes=CR,No]: ' READ (5,'(A)') name sym = 1 IF ((name(1:1) .EQ. 'n') .OR. (name(1:1) .EQ. 'N')) sym = 0 name = filename1(:index(filename1,' ')-1)//'.hbook' PRINT *,'Reading the event data '//name(:index(name,' ')) CALL hlimit(nwpawc) CALL hropen(20,'Input',name,' ',4096,istat) CALL hcdir ('//Input',' ') CALL hgnpar(20,'MAIN') CALL hnoent(20,cnt) cnt = min(cnt,NEVENTS) DO i=1,cnt CALL HGNF (20,i,event,j) DO j=1,6 u(i,j)=event(j) ENDDO DO j=1,4 v(i,j)=event(j+6) ENDDO ENDDO CALL hrend('Input') print *,'Creating the matrix' CALL initExponents (minord,maxord) print *,'Fitting the coeff''s for the delta reconstruction' CALL findCoeffs (6,4,sym*2) print *,'Fitting the coeff''s for the dx/dz reconstruction' CALL findCoeffs (2,1,sym*2) print *,'Fitting the coeff''s for the y reconstruction' CALL findCoeffs (3,2,sym*1) print *,'Fitting the coeff''s for the dy/dz reconstruction' CALL findCoeffs (4,3,sym*1) name = filename2(:index(filename2,' ')-1)//'.map' print *,'Saving the matrix '//name(:index(name,' ')) OPEN (UNIT=97,FILE=name,STATUS='UNKNOWN') DO i=1,num WRITE (97,'(4(E14.7E2,X),4I1)') > c(i,1),c(i,2),c(i,3),c(i,4), > e(i,1),e(i,2),e(i,3),e(i,4) ENDDO CLOSE (97) name = filename2(:index(filename2,' ')-1)//'.error' print *,'Saving the error matrix '//name(:index(name,' ')) OPEN (UNIT=97,FILE=name,STATUS='UNKNOWN') WRITE (97,'(A,I16)') 'NEVENTS=',cnt WRITE (97,'(A)') 'CHISQ =' WRITE (97,'(X,4(E14.7E2,X))') > chisq(1),chisq(2),chisq(3),chisq(4) WRITE (97,'(A)') 'ERRORS =' DO i=1,num WRITE (97,'(X,4(E14.7E2,X),4I1)') > err(i,1),err(i,2),err(i,3),err(i,4), > e(i,1), e(i,2), e(i,3), e(i,4) ENDDO CLOSE (97) STOP END SUBROUTINE findCoeffs (i1,i2,sym) IMPLICIT NONE INTEGER i1,i2,sym * -- find the coeffs ! commonblocks holding the event data and the reconstruction matrix INTEGER NTERMS,num PARAMETER (NTERMS = 300) INTEGER e(NTERMS,4) REAL c(NTERMS,4),err(NTERMS,4),chisq(4) COMMON /matrix/ c,err,e,num,chisq INTEGER NEVENTS,cnt PARAMETER (NEVENTS = 100000) REAL u(NEVENTS,6),v(NEVENTS,4) COMMON /events/ u,v,cnt ! other variables INTEGER order INTEGER e1(NTERMS), e2(NTERMS), e3(NTERMS), e4(NTERMS) COMMON /selected/ order,e1,e2,e3,e4 REAL x(NEVENTS),y(NEVENTS),sigma(NEVENTS),a(NTERMS) REAL uu(NEVENTS,NTERMS),vv(NTERMS,NTERMS),ww(NTERMS) REAL cvm(NTERMS,NTERMS) INTEGER i,j,n EXTERNAL funct ! create the input vector DO i=1,cnt x(i) = i y(i) = u(i,i1) sigma(i) = 1. ENDDO ! apply a given symmetry n = 0 order = 0 DO i=1,num j = e(i,1)+e(i,2)+e(i,3)+e(i,4) IF ((sym .EQ. 0) .OR. (j .EQ. 0) .OR. > (mod(sym,2) .EQ. mod(e(i,3)+e(i,4),2))) THEN n = n+1 e1(n) = e(i,1) e2(n) = e(i,2) e3(n) = e(i,3) e4(n) = e(i,4) order = max (order,j) ENDIF ENDDO ! calculate the matrix coefficients CALL svdfit(x,y,sigma,cnt,a,n, > uu,vv,ww,NEVENTS,NTERMS,chisq(i2),funct) CALL svdvar(vv,n,NTERMS,ww,cvm,NTERMS) ! and fill them into the matrix array according the symmetry n = 0 DO i=1,num IF ((sym .EQ. 0) .OR. > (e(i,1)+e(i,2)+e(i,3)+e(i,4) .EQ. 0) .OR. > (mod(sym,2) .EQ. mod(e(i,3)+e(i,4),2))) THEN n = n+1 c (i,i2) = a(n) err (i,i2) = SQRT(cvm(n,n)) ELSE c (i,i2) = 0 err (i,i2) = 0 ENDIF ENDDO RETURN END SUBROUTINE funct(ix,a,m) IMPLICIT NONE INTEGER m REAL ix,a(m) * -- function passed to the fit routine ! commonblocks holding the event data and the reconstruction matrix INTEGER NTERMS PARAMETER (NTERMS = 300) INTEGER order INTEGER e1(NTERMS), e2(NTERMS), e3(NTERMS), e4(NTERMS) COMMON /selected/ order,e1,e2,e3,e4 INTEGER NEVENTS,cnt PARAMETER (NEVENTS = 100000) REAL u(NEVENTS,6),v(NEVENTS,4) COMMON /events/ u,v,cnt ! other variables INTEGER i REAL vv1(0:10),vv2(0:10),vv3(0:10),vv4(0:10) i = INT(ix) vv1(0) = 1. vv2(0) = 1. vv3(0) = 1. vv4(0) = 1. vv1(1) = v(i,1) vv2(1) = v(i,2) vv3(1) = v(i,3) vv4(1) = v(i,4) DO i=2,order vv1(i)=vv1(i-1)*vv1(1) vv2(i)=vv2(i-1)*vv2(1) vv3(i)=vv3(i-1)*vv3(1) vv4(i)=vv4(i-1)*vv4(1) ENDDO DO i=1,m a(i) = vv1(e1(i)) * vv2(e2(i)) * vv3(e3(i)) * vv4(e4(i)) ENDDO RETURN END SUBROUTINE initExponents (minord,maxord) IMPLICIT NONE INTEGER minord,maxord * -- creates the exponent table for the reconstruction matrix * * minord I : minimal order of the matrix * maxord O : maximal order of matrix elements ! commonblock holding the reconstruction matrix INTEGER NTERMS,num PARAMETER (NTERMS = 300) INTEGER e(NTERMS,4) REAL c(NTERMS,4),err(NTERMS,4),chisq(4) COMMON /matrix/ c,err,e,num,chisq ! other variables INTEGER ee(4),i,j,k num = 0 DO j=minord,maxord DO i=1,4 ee(i) = j ENDDO DO WHILE (ee(1) .GE. 0) k=0 DO i=1,4 k=k+ee(i) ENDDO IF (k .EQ. j) THEN num=num+1 DO i=1,4 e(num,i)=ee(i) ENDDO ENDIF i = 4 ee(i)=ee(i) - 1 DO WHILE ((i .GT. 1) .AND. (ee(i) .LT. 0)) ee(i) = j i = i-1 ee(i) = ee(i)-1 ENDDO ENDDO ENDDO RETURN END