CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C H M O M R E S C C C C COLLECTION OF HELP-SUBROUTINES FOR MAIN PROGRAM M O M R E S C C C C DATE: AUG-1985 C C C C CHANGES: UNIX VERSION 04/27/2000 C C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C SUBROUTINE I G L FOR MAIN PROGRAM I M P A U F C C C C PURPOSE: SOLVES AN INHOMOGENEOUS SYSTEM OF 3 LINEAR EQUATIONS C C OF THE TYPE A * X = B C C C C WRITTEN BY: B.A.M. C C C C CHANGES: C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C 1 2 3 4 5 6 7 8 C2345678901234567890123456789012345678901234567890123456789012345678901234567890 C SUBROUTINE IGL(A,X,B,D) SAVE DIMENSION A(3,3),X(3),B(3),AV(3,3) D=DET(A) IF(D.EQ.0.) RETURN DO 1 K=1,3 DO 2 I=1,3 DO 3 J=1,3 AV(I,J)=A(I,J) 3 IF(J.EQ.K) AV(I,J)=B(I) 2 CONTINUE X(K)=DET(AV)/D 1 CONTINUE RETURN END C C FUNCTION DET(A) C CALCULATES THE DETERMINANT OF A 3*3 MATRIX ACCORDING TO SARRUS' RULE DIMENSION A(3,3) DET=A(1,1)*A(2,2)*A(3,3) + 1 A(1,2)*A(2,3)*A(3,1) + 2 A(1,3)*A(2,1)*A(3,2) - 3 A(1,3)*A(2,2)*A(3,1) - 4 A(1,1)*A(2,3)*A(3,2) - 5 A(1,2)*A(2,1)*A(3,3) RETURN END C C SUBROUTINE DGLABW(PIMP,NSTEP,SS,F123,IMIN,IMAX,ICONT) SAVE C C CALCULATES THE SOLUTION OF THE DIFFERENTIAL EQUATION FOR THE DEVIATIONS C FROM THE REFERENCE TRAJECTORY IN AN INHOMOGENEOUS MAGNETIC FIELD C C THIS IS THE MOST PRIMITIVE WAY TO TREAT A SECOND ORDER DIFFERENTIAL C EQUATION (THE RUNGE-KUTTA METHOD WOULD BE FAR MORE ACCURATE) C C FOR THE DIFFERENTIAL EQUATION SEE C MECKING'S LOGBOOK M2/? C OR NOELDEKE: LECTURE ON PARTICLE ACCELERATORS C OR STEFFEN : HIGH ENERGY BEAM OPTICS C OR CALCULATE IT YOURSELF C DIMENSION NSTEP(50),SS(0:4,0:50),F123(3,0:50),XAKT(3) DIMENSION XSAKT(3),A(3),ABL2(3),XA(3) C C THE MATRIX SS(I,J) CONTAINS THE PARTICLE POSITION FOR STEP # J C I=0 TRACK LENGTH C I=1,3 SPATIAL COORDINATES C I=4 INVERSE BENDING RADIUS C GOTO(100,200,300),ICONT C 100 CONTINUE C SET THE INITIAL CONDITIONS XAKT(1)=1. XAKT(2)=0. XAKT(3)=0. XSAKT(1)=0. XSAKT(2)=1. XSAKT(3)=0. A(1)=0. A(2)=0. A(3)=1. C C STEPWISE INTEGRATION OF THE DIFFERENTIAL EQUATION SA=0. C DO 102 I=IMIN,IMAX FNS=FLOAT(NSTEP(I)) DS=(SS(0,I)-SS(0,I-1))/FNS DS22=DS**2/2. DO 103 K=1,NSTEP(I) SA=SA+DS C INTERPOLATION BETWEEN SUBSEQUENT INTERFACES DO 111 I1=1,3 IF(FNS.GE.1) 1 XA(I1)=SS(I1,I-1)+(SS(I1,I)-SS(I1,I-1))*FLOAT(K)/FNS IF(FNS.LT.1) XA(I1)=SS(I1,I) 111 CONTINUE C IF(FNS.GE.1) 1 RADI=SS(4,I-1)+(SS(4,I)-SS(4,I-1))*FLOAT(K)/FNS IF(FNS.LT.1) RADI=SS(4,I) FACB0=1. C DO 101 J=1,3 ABL2(J)=(A(J)-XAKT(J)*RADI)*RADI XAKT(J)=XAKT(J)+XSAKT(J)*DS+ABL2(J)*DS22 XSAKT(J)=XSAKT(J)+ABL2(J)*DS 101 CONTINUE 103 CONTINUE C DO 104 J=1,3 104 F123(J,I)=XAKT(J) C 102 CONTINUE C RETURN C 200 CONTINUE RETURN C 300 CONTINUE c WRITE(3,1002) c DO 105 I=1,IMAX c105 WRITE(3,1001) I,(F123(J,I),J=1,3) RETURN 1000 FORMAT(1H ,10X,3I4,8F10.4) 1001 FORMAT(1H ,5X,I2,5X,8F12.6) 1002 FORMAT(1H-,' MATRIX-COEFFICIENTS FOR X-DEVIATION ',//, 1 ' NUMBER MATRIX-COEFFICIENTS X(X0,THETA0,DP/P) ',/) END C C SUBROUTINE RESOLX(AMAT,F123,SIGX,XAUF,IMAX) SAVE C C SUBPROGRAM FOR SR CHIMP OR IMPAUF C CALCULATES THE RESOLUTION IN X0,THETA0 AND DELTA(P)/P FOR A C GIVEN POSITION RESOLUTION SIGX C DOUBLE PRECISION AMAT,BMAT DIMENSION AMAT(3,3),BMAT(3,3),F123(3,0:50),XAUF(3) DIMENSION SIGX(50) C C MATRIX-INVERSION BMAT = AMAT**-1 DO 10 I=1,3 DO 10 J=1,3 10 BMAT(I,J)=AMAT(I,J) C CALL MATINV(BMAT,3,DET) C DO 1 K=1,3 XK2=0. DO 2 I=1,IMAX SU=0. DO 3 J=1,3 IF(SIGX(I).LE.0.) GOTO 3 SU=SU+BMAT(K,J)*F123(J,I)/SIGX(I) 3 CONTINUE XK2=XK2+SU**2 2 CONTINUE XAUF(K)=SQRT(XK2) 1 CONTINUE RETURN END C C SUBROUTINE RESOLV(AM12,AMAT,F123,SIGX2,SS,NSTEP,X0,P, 1 BETA,XAUF,IMAX) SAVE C C SUBPROGRAM FOR SR CHIMP OR IMPAUF C CALCULATES THE RESOLUTION IN X0,THETA0 AND DELTA(P)/P C CAUSED BY MULTIPLE SCATTERING C C P MOMENTUM IN GEV/C C C DOUBLE PRECISION AMAT DIMENSION AMAT(3,3),B(3),X(3),F123(3,0:50),SIGX2(50) DIMENSION XAUF(3),NSTEP(50),SS(0:4,0:50),X0(50) DIMENSION X2(3),AM12(50,50) ckm test for pion: is it correct betta and momentum? real m_pion,imp,skor m_pion=.139 imp=P scor=imp/sqrt(imp*imp+m_pion*m_pion) if(abs(BETA-scor).gt.1.e-7)then print*,'RESOLV:P,BETTA,scor,M',P,BETTA,scor,m_pion endif C SA=0. DO 110 I=1,3 110 X2(I)=0. C C DETERMINATION OF THE CORRECTION FACTOR FOR THE SCATTERING ANGLE XTOT=0. C DO 112 I=1,IMAX IF(X0(I).GT.0.) XTOT=XTOT+(SS(0,I)-SS(0,I-1))/X0(I) 112 CONTINUE C FCOR=1. IF(XTOT.GT.0.) FCOR=1.+ALOG10(XTOT/2.)/9. C DO 102 I=1,IMAX DS=(SS(0,I)-SS(0,I-1))/FLOAT(NSTEP(I)) C C CALCULATION OF DELTA(THETA**2) DTH=0. IF((DS.GT.0.).AND.(X0(I).GT.0.)) 1 DTH=(.014/(P*BETA))*SQRT(DS/X0(I))*FCOR !!!! SIGMA !!! C DO 103 I1=1,NSTEP(I) SA=SA+DS C C SET VECTOR B DO 200 K=1,3 S1=0. DO 201 IM=1,IMAX IF((SS(0,IM)-SS(0,0)).LT.SA) GOTO 201 IF(SIGX2(IM).LE.0.) GOTO 201 C C APPRROXIMATE DETERMINATION OF THE COEFFICIENT M12 BETWEEN THE POSITION C SA AND THE DATA POINT #IM BY INTERPOLATION F2=AM12(I,IM) IF(DS.GT.0.) F2=AM12(I,IM)+(SA+SS(0,0)-SS(0,I-1))* 1 (AM12(I+1,IM)-AM12(I,IM))/(SS(0,I)-SS(0,I-1)) C S1=S1+F2*F123(K,IM)/SIGX2(IM) 201 CONTINUE! END IM B(K)=S1*DTH 200 CONTINUE C C SOLVE THE INHOMOGENEOUS SYSTEM OF LINEAR EQUATIONS AMAT*X=B C CALL IGL(AMAT,X,B,DET) C DO 202 K=1,3 X2(K)=X2(K)+X(K)**2 202 CONTINUE C 103 CONTINUE 102 CONTINUE DO 101 K=1,3 XAUF(K)=SQRT(X2(K)) 101 CONTINUE RETURN END C C SUBROUTINE MATINV(ARRAY,NORDER,DET) SAVE C C DOUBLE PRECISION ARRAY,AMAX,SAVE DIMENSION ARRAY(3,3),IK(10),JK(10) C DET=1. DO 100 K=1,NORDER C C FIND LARGEST ELEMENT ARRAY(I,J) IN REST OF MATRIX C AMAX=0. 21 DO 30 I=K,NORDER DO 30 J=K,NORDER IF(ABS(AMAX)-ABS(ARRAY(I,J))) 24,24,30 C IF(DABS(AMAX)-DABS(ARRAY(I,J))) 24,24,30 24 AMAX=ARRAY(I,J) IK(K)=I JK(K)=J 30 CONTINUE C C INTERCHANGE ROWS AND COLUMNS TO PUT AMAX IN ARRAY(K,K) C IF(AMAX) 41,32,41 32 DET=0. GOTO 140 41 I=IK(K) IF(I-K) 21,51,43 43 DO 50 J=1,NORDER SAVE=ARRAY(K,J) ARRAY(K,J)=ARRAY(I,J) 50 ARRAY(I,J)=-SAVE 51 J=JK(K) IF(J-K) 21,61,53 53 DO 60 I=1,NORDER SAVE=ARRAY(I,K) ARRAY(I,K)=ARRAY(I,J) 60 ARRAY(I,J)=-SAVE C C ACCUMULATE ELEMENTS OF INVERSE MATRIX C 61 DO 70 I=1,NORDER IF(I-K) 63,70,63 63 ARRAY(I,K)=-ARRAY(I,K)/AMAX 70 CONTINUE DO 80 I=1,NORDER DO 80 J=1,NORDER IF(I-K) 74,80,74 74 IF(J-K) 75,80,75 75 ARRAY(I,J)=ARRAY(I,J)+ARRAY(I,K)*ARRAY(K,J) 80 CONTINUE DO 90 J=1,NORDER IF(J-K) 83,90,83 83 ARRAY(K,J)=ARRAY(K,J)/AMAX 90 CONTINUE ARRAY(K,K)=1./AMAX 100 DET=DET*AMAX C C RESTORE ORDERING OF MATRIX C DO 130 L=1,NORDER K=NORDER-L+1 J=IK(K) IF(J-K) 111,111,105 105 DO 110 I=1,NORDER SAVE=ARRAY(I,K) ARRAY(I,K)=-ARRAY(I,J) 110 ARRAY(I,J)=SAVE 111 I=JK(K) IF(I-K) 130,130,113 113 DO 120 J=1,NORDER SAVE=ARRAY(K,J) ARRAY(K,J)=-ARRAY(I,J) 120 ARRAY(I,J)=SAVE 130 CONTINUE C 140 RETURN END