CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C MAIN : M O M R E S C C C C PURPOSE: CALCULATES AND PLOTS MOMENTUM RESOLUTION C C USING SR MORE C C C C WRITTEN BY: B.A.M. C C C C DATE: 23-AUG-85 C C C C CHANGES: WHAT WHEN WHO C C C C UNIX VERSION 4/27/0 BAM C c linux version 2/01/05 kmikhail@jlab.org c add hbook 2/01/05 kmikhail@jlab.org C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C 1 2 3 4 5 6 7 8 C2345678901234567890123456789012345678901234567890123456789012345678901234567890 C PROGRAM TO CALCULATE POSITION, ANGULAR AND MOMENTUM RESOLUTION C IN INHOMOGENEOUS MAGNETIC FIELDS C program momres ckm CHARACTER*16 FNAME ckm CHARACTER*40 FNAME CHARACTER*60 FNAME CHARACTER*9 DATUM CHARACTER*5 TEKO(20),SNAME(50,4),DUMMY(30) DIMENSION NSTEP(50),SS(0:4,0:50) DIMENSION SIGX(50),X0(50),BVEK(3) DIMENSION RESO(3,3),XA(3),RI(3) real cent_forw !angle theta real Bxdl equivalence (icent_forw,cent_forw) ckm hbook: integer iw common/pawc/iw(1000000) call hlimit(1000000) print*,'+-----------------------------------------------------+' print*,'| momres.for |' print*,'+-----------------------------------------------------+' C C INPUT DATA C WRITE(6,1033) ckm READ(5,2002) FNAME c fname='DAT_ANGLE/5thin.dat' * fname='DAT_ANGLE/10thin.dat' c fname='DAT_ANGLE/15thin.dat' c fname='DAT_ANGLE/20thin.dat' c fname='DAT_ANGLE/25thin.dat' c fname='DAT_ANGLE/30thin.dat' c fname='DAT_ANGLE/35thin.dat' c fname='DAT_ANGLE/40thin.dat' c fname='dc90.dat' c fname='sv90.dat' c fname='cdc_pion.dat' fname='momres.dat' c * * fname='DAT_ANGLE/20thin1p5.dat' * fname='DAT_ANGLE/test.dat' *--- Plaing with different configurations: * move all dat in DAT directory * fname='20gail.dat' * fname='20gail-stv-dc.dat' * fname='20gail-1stv-dc-cc.dat' * fname='20gail-fix-target-dc-cc.dat' c----------------------- print*,'filename:',fname ckm read vertex ON(1) or OFF(0) ckm WRITE(6,1040) ckm READ(5,1000) NVERTEX NVERTEX=0 c NVERTEX=1 if(NVERTEX.EQ.1)print*,'WE USE VERTEX, NVERTEX=',NVERTEX if(NVERTEX.EQ.0)print*,'WE DO NOT USE VERTEX, NVERTEX=',NVERTEX if(NVERTEX.LT.0.OR.NVERTEX.GT.1)then print*,'ERR: Wrong vertex set!!!,NVERTEX=',NVERTEX stop endif ckm Number of bins for momenta region ckm ISTEP=10 ISTEP=100 C OPEN(UNIT=1,ACCESS='SEQUENTIAL',STATUS='OLD',FILE=FNAME) C READ THE COMMENT READ(1,1001) (TEKO(IX),IX=1,20) WRITE(6,1014) (TEKO(IX),IX=1,20) READ(1,1001) (DUMMY(IX),IX=1,30) C READ(1,1003) AM,PMIN,PMAX,icent_forw changed to real ES 10/03/08, see equiv above READ(1,1003) AM,PMIN,PMAX,cent_forw call myhbook(pmin,pmax,istep) C c CALL DATE(DATUM) datum='March2005' C C READ DATA FOR PARTICLE TRACKING (MAX. 50) READ(1,1001) (DUMMY(IX),IX=1,10) READ(1,1008) NST1,SS(0,0) IF(NVERTEX.NE.1) READ(1,1001) (DUMMY(IX),IX=1,10) !skip vertex line SUMX0=0. C DO 101 I=1,50 !50 -> max items of material in dat file READ(1,1008) NSTEP(I),SS(0,I),SIGX(I),X0(I),(SNAME(I,J),J=1,4) ckm test print: print*,i,NSTEP(I),SS(0,I),SIGX(I),X0(I),(SNAME(I,J),J=1,4) C WRITE(3,2008) NSTEP(I),SS(0,I),SIGX(I),X0(I),(SNAME(I,J),J=1,4) C SS(2,I)=SS(0,I) C C NSTEP NUMBER OF STEPS (TO SOLVE THE DIFFERENTIAL EQUATION) C SS(0,I) POSITION OF DATA POINT ON THE TRACK C THE ORIGIN IS AT SS(0,0) C C SIGX POSITION RESOLUTION IN METER C (NULL OR NEGATIV TO MARK THE INTERFACE BETWEEN TWO DIFFERENT C MATERIALS; NO POSITION INFORMATION AT THIS POINT C C X0 RADIATION LENGTH (IN METER) IF(NSTEP(I).LE.0) Then !break the loop print*,'Number of materials in dat file:',i,IMAX GOTO 102 EndIf IF(X0(I).GT.0.) then call hf1(1,SS(0,I),(SS(0,I)-SS(0,I-1))/X0(I)) SUMX0=SUMX0+(SS(0,I)-SS(0,I-1))/X0(I) print*,'i=',i,' SS(0,I),(SS(0,I)-SS(0,I-1))/X0(I)', . SS(0,I),(SS(0,I)-SS(0,I-1))/X0(I) else call hf1(1,SS(0,I),4.) endif IMAX=I 101 CONTINUE C 102 CONTINUE CLOSE(UNIT=1) C DP=(PMAX-PMIN)/FLOAT(ISTEP-1) C 300 CONTINUE C WRITE(3,1015) WRITE(3,1024) DATUM WRITE(3,1014) (TEKO(IX),IX=1,20) WRITE(3,1016) WRITE(3,1010) C print*,'+-----------------------------------------------------+' print*,'| end of reading input file |' print*,'+-----------------------------------------------------+' c if(cent_forw.eq.20.) c & CALL FIELD(XA,BVEK,RI,1.,1.,1.,0,1) !initialization of B-field c CALL FIELD(XA,BVEK,RI,1.,1.,1.,0,1) !initialization of B-field c if(cent_forw.eq.90.) c & CALL FIELD_CONST(XA,BVEK,RI,1.,1.,1.,0,1) !initialization of B-field c c use constant field option for halld 10/01/06 c after renaning the files, temporarily fieldr has the constant field option. c CALL FIELD (XA,BVEK,RI,1.,1.,1.,0,1) !initialization of B-field c CALL FIELD_CONST(XA,BVEK,RI,1.,1.,1.,0,2) !print B-field ckm XA(3) 3-vector coordinates of point on the track ckm our S is radius-vector, and theta=20degree and phi=0degree ckm fill B-field histograms XA(1)=0. XA(2)=0. do i=1,1000 XA(3)=(float(i-1)/1000.*5.)+0.001 CALL FIELD(XA,BVEK,RI,1.,1.,1.,0,2) call hf1(2,XA(3),BVEK(1)) !X call hf1(3,XA(3),BVEK(2)) !Y call hf1(4,XA(3),BVEK(3)) !Z enddo ckm put B-field in array BVEK XA(1)=0. XA(2)=0. DO 301 I=0,IMAX XA(3)=SS(0,I) CALL FIELD(XA,BVEK,RI,1.,1.,1.,0,2) ckm test print:XA,BVEK,RI c print 11038,XA,' |',BVEK,' |',RI C IF(ABS(RI(3)).GT.1.0E-05) THEN RAD=1./ABS(RI(3)) ELSE RAD=0. ENDIF C IF(I.EQ.0) THEN WRITE(3,1111) I,SS(0,I),BVEK(3),RAD WRITE(6,1111) I,SS(0,I),BVEK(3),RAD GOTO 301 ENDIF C IF(SIGX(I).GT.0.) WRITE(3,1011) I,NSTEP(I),SS(0,I),SIGX(I), 1 BVEK(3),RAD,X0(I),(SNAME(I,J),J=1,4) C IF(SIGX(I).LE.0.) WRITE(3,1012) I,NSTEP(I),SS(0,I),BVEK(3), 1 RAD,X0(I),(SNAME(I,J),J=1,4) IF(SIGX(I).GT.0.) WRITE(6,1011) I,NSTEP(I),SS(0,I),SIGX(I), 1 BVEK(3),RAD,X0(I),(SNAME(I,J),J=1,4) C IF(SIGX(I).LE.0.) WRITE(6,1012) I,NSTEP(I),SS(0,I),BVEK(3), 1 RAD,X0(I),(SNAME(I,J),J=1,4) 301 CONTINUE C CALL B_INT(NSTEP,SS,1,IMAX,BI) WRITE(3,1041) SUMX0,BI WRITE(6,1041) SUMX0,BI Bxdl=BI C DO 201 IP=1,ISTEP P=PMIN+FLOAT(IP-1)*DP C C INITIALIZE FIELD SUBROUTINE AND TRANSFER THE MOMENTUM P CALL MORE(P,AM,NSTEP,SS,SIGX,X0,IMAX,RESO,0,2) C BETA=P/SQRT(P**2+AM**2) C DO 221 IO1=0,1 IF((IO1.EQ.1).AND.(IOUT.EQ.5)) GOTO 221 IF(IO1.EQ.0) IOU=3 IF(IO1.EQ.1) IOU=6 WRITE(IOU,1016) WRITE(IOU,1013) P,AM,BETA WRITE(IOU,1028) WRITE(IOU,1029) (RESO(1,K),K=1,3) WRITE(IOU,1030) (RESO(2,K),K=1,3) WRITE(IOU,1031) (RESO(3,K),K=1,3) 221 CONTINUE BETAINV=1./BETA WRITE(10,*)P,BETAINV,(RESO(K,1),K=1,3) WRITE(11,*)P,BETAINV,(RESO(K,2),K=1,3) WRITE(12,*)P,BETAINV,(RESO(K,3),K=1,3) ckm dP/P call hf1(11,p,reso(1,3)) !pos (resolution) call hf1(12,p,reso(2,3)) !mult scat call hf1(10,p,reso(3,3)) !sum + pos+scat c print*,'-- -- p,reso(3,3)', p,reso(3,3)*100. ckm dTheta (Radians) call hf1(21,p,reso(1,2)) !pos (resolution) call hf1(22,p,reso(2,2)) !mult scat call hf1(20,p,reso(3,2)) !sum + pos+scat ckm dX (microns) call hf1(31,p,reso(1,1)*1000000.) !pos (resolution) call hf1(32,p,reso(2,1)*1000000.) !mult scat call hf1(30,p,reso(3,1)*1000000.) !sum + pos+scat CALL MORE(P,AM,NSTEP,SS,SIGX,X0,IMAX,RESO,0,3) C 201 CONTINUE ckm hbook finish: CALL hrput(0,'momres.hbook','T') print *, 'Histogramms are saved into the file : ' & ,'momres.hbook' C C CALL DGLABW(PIMP,NSTEP,SS,F123,IMIN,IMAX,3) 400 CONTINUE print*,'filename:',fname,'Bxdl=',Bxdl 1000 FORMAT(I10) 1001 FORMAT(10A5) 2001 FORMAT(' ',10A5) 1002 FORMAT(3I1) ckm 1003 FORMAT(8F10.4) 1003 FORMAT(8F10.6) 2003 FORMAT(' ',8F10.4) 1004 FORMAT(I2) 1007 FORMAT(1I10) ckm it doesn't matter which Fy.x format we use ckm 1008 FORMAT(I3,2X,3F10.4,4A5) 1008 FORMAT(I3,2X,3F10.8,4A5) 1009 FORMAT(F10.5) C' DATA POINTS FOR PARTICLE TRACKING ',//, 1010 FORMAT(1H , 1 ' I NSTEP POSITION SIGMA(X) FIELD ', 2 'RADIUS X0 ',/) 1011 FORMAT(1H ,I2,2X,I4,F10.4,F10.6,F9.4,2F10.2,2X,4A5) 1012 FORMAT(1H ,I2,2X,I4,F10.4,10X,F9.4,2F10.2,2X,4A5) 1013 FORMAT(//,' P = ',F7.3,' GeV/c Mass = ',F7.4,' GeV ', 1 ' Beta = ',F6.3,/) 1014 FORMAT(' ',10A5) 1015 FORMAT(' O U T P U T OF M A I N M O M R E S ',/) 1016 FORMAT(1H-) 1017 FORMAT(' MOMENTUM (IN GEV/C) ',$) 1018 FORMAT(1H ,/,' MOMENTUM =',F8.3, 1 ' GEV/C , MASS = ',F8.3,' GEV '/) 1019 FORMAT(' MASS IN GEV ',$) 1020 FORMAT(' PRINTOUT MATRIX COEFFICIENTS ? YES=1 ',$) 1021 FORMAT(10A10) 1023 FORMAT(1H ,8(E12.6,2X)) 1024 FORMAT(' DATE: ',A9) 1028 FORMAT(1H ,' RESOLUTION (sigma) IN X0 THETA0 ', 1 ' DP/P'/) 1029 FORMAT(1H ,' POSITION ',3(E12.6,2X),/) 1030 FORMAT(1H ,' SCATTERING ',3(E12.6,2X),/) 1031 FORMAT(1H ,' SUM ',3(E12.6,2X),/) 1032 FORMAT(' ',5(E12.6,2X)) 1033 FORMAT(' INPUT FILE FOR MAIN MOMRES, FILENAME: ',$) 1034 FORMAT(/,' OUTPUT DEVICE # (return -> NO OUTPUT) ',$) 1035 FORMAT(/,' OUTPUT FOR ALL MOMENTA = 1 ',/, 1 ' ONLY FOR 1. AND LAST MOMENTUM = 0 ',$) 1036 FORMAT(/,' MAIN M O M R E S READING DATA FILE "MOMRES" ',/) 1037 FORMAT(/,' PRINTOUT MATRIX COEFFICIENTS? YES=1 ',$) 1040 FORMAT(/,' USE VERTEX (FIRST LINE) YES=1 ',$) 1041 FORMAT(/,' Total Thickness ',F10.4,' X0 ', 1 ' Int B x dl = ',F10.4,' Tm ') 1111 FORMAT(1H ,I2,6X,F10.4,10X,F9.4,2F10.2) 2002 FORMAT(A16) 2008 FORMAT(I3,2X,3F10.4,2X,4A5) END C C SUBROUTINE B_INT(NSTEP,SS,IMIN,IMAX,BI) SAVE C C CALCULATES THE INTEGRAL B x dl along the trajectory C DIMENSION NSTEP(50),SS(0:4,0:50),XA(3),BVEK(3),RI(3) C C THE MATRIX SS(I,J) CONTAINS THE PARTICLE POSITION FOR STEP # J C I=0 TRACK LENGTH C CALL FIELD(XA,BVEK,RI,1.,1.,1.,0,1) XA(1)=0. XA(2)=0. BI=0. C SA=SS(0,0) DO 102 I=IMIN,IMAX FNS=FLOAT(NSTEP(I)) DS=(SS(0,I)-SS(0,I-1))/FNS DO 103 K=1,NSTEP(I) SA=SA+DS XA(3)=SA CALL FIELD(XA,BVEK,RI,1.,1.,1.,0,2) BI=BI+BVEK(3)*DS 103 CONTINUE 102 CONTINUE RETURN END