SUBROUTINE CALC_PATHLEN(ICELL,ILAY,ISLAY,ITRACK,LPATH,PATHLEN) IMPLICIT NONE $$INCLUDE 'KOFIA$LIB:UTC_GEOM.CMN' $$INCLUDE 'KOFIA$LIB:LUNS.CMN' $$INCLUDE 'KOFIA$LIB:DC_TRACK.PAR' $$INCLUDE 'KOFIA$LIB:DC_TRACK.CMN' $$INCLUDE 'KOFIA$LIB:INFO.CMN' INTEGER nintersect INTEGER ICELL,ILAY,ISLAY INTEGER ITRACK REAL*4 PATHLEN, pathlenxy REAL*4 rinner, router, phiedge1, phiedge2, phioffset REAL*4 XCEN,YCEN,PHICEN,LCEN,RADIUS REAL*4 PHI,arg,I REAL*4 rtest,phitest REAL*4 rintersect(4),phiintersect(4) REAL*4 angleint(4),angle(4),anglespan REAL*4 PI LOGICAL*4 LPATH PI = 2.0*ASIN(1.0) nintersect = 0 c WRITE(LOUT,*) ' ' c WRITE (LLOG,*) ' ' c WRITE (LOUT,*) 'SLAY,LAY,CELL',ISLAY,ILAY,ICELL c WRITE (LLOG,*) 'SLAY,LAY,CELL',ISLAY,ILAY,ICELL c- To calculate the path length, we find the intersection of the track c- with the boundaries of the cell. c- The variables describing the positions of the boundaries c- as well as any corrections to this geometry are found c- in utc_geom.cmn, and the data is stored in utc_init.inc. c- The track data is found in the DC_TRACK.CMN common block. c ------ c Fetch the relevant geometry c ------ IF (ILAY .eq. 1) THEN rinner = RINL(ISLAY) + DELTAR_W(ISLAY) ELSE rinner= RW((ILAY-1)*2,ISLAY) + DELTAR_W(ISLAY) ENDIF IF (ILAY .eq. NLAY(ISLAY)) THEN router=ROUTL(ISLAY) + DELTAR_W(ISLAY) ELSE router= RW(ILAY*2,ISLAY) + DELTAR_W(ISLAY) ENDIF c the edge of the cell is at phiedge degrees from the x axis c We ensure than phiedge is between -pi and pi c Depending on the layer, the first cell will have its edge on the c x-axis, or its centre on the x-axis. C in the even layers, the cells are centred on the x-axis IF (MOD(ILAY,2).EQ.0) THEN phioffset = .5 ELSE phioffset = 0 ENDIF phiedge1 = (((2*PI)/NCELL(ILAY,ISLAY))*(ICELL-1-phioffset)) + & ROT_W(ISLAY) IF (phiedge1 .gt. PI) phiedge1 = phiedge1-(2*PI) phiedge2 = (((2*PI)/NCELL(ILAY,ISLAY))*(ICELL-phioffset)) + & ROT_W(ISLAY) IF (phiedge2 .gt. PI) phiedge2 = phiedge2-(2*PI) c write(lout,*) 'rinner',rinner c write(lout,*) 'router',router c write(lout,*) 'phiedge1',phiedge1*(180/PI) c write(lout,*) 'phiedge2',phiedge2*(180/PI) c write(llog,*) 'rinner',rinner c write(llog,*) 'router',router c write(llog,*) 'phiedge1',phiedge1*(180/PI) c write(llog,*) 'phiedge2',phiedge2*(180/PI) c -------- c Compute some of the track parameters c -------- c In Cartesian coordinates w.r.t. the center of the superlayer, C (X-XCEN)**2+(Y-YCEN)**2=RADIUS**2 XCEN=XCEN_D(ITRACK)-TRAN_W_X(ISLAY) YCEN=YCEN_D(ITRACK)-TRAN_W_Y(ISLAY) c The centre of curvature of the track in polar coord. is LCEN = SQRT(XCEN**2+YCEN**2) PHICEN = ATAN2(YCEN,XCEN) RADIUS = RADIUS_D(ITRACK) c write(lout,*) 'lcen',lcen c write(lout,*) 'phicen',phicen*(180/pi) c write(lout,*) 'radius',radius c write(llog,*) 'lcen',lcen c write(llog,*) 'phicen',phicen*(180/pi) c write(llog,*) 'radius',radius c ------- C Check for intersections at the side edges of the cell c ------- c Using the cosine law, we find radial position of the pts. of intersection c w.r.t. the center of the superlayer C phiedge1 PHI=ABS(PHICEN-phiedge1) C Ensure that the angle between the track center and c phiedge1 is acute IF (PHI.GT.PI) PHI=2*PI-PHI IF (RADIUS**2-((LCEN*SIN(PHI))**2).GE.0) THEN rtest = LCEN*COS(PHI) + SQRT(RADIUS**2-((LCEN*SIN(PHI))**2)) c WRITE(LOUT,*) 'rtest',rtest c WRITE(LLOG,*) 'rtest',rtest IF (rtest .le. router .and. rtest .ge. rinner) THEN nintersect=nintersect+1 rintersect(nintersect)=rtest phiintersect(nintersect)=phiedge1 ENDIF rtest = LCEN*COS(PHI) - SQRT(RADIUS**2-((LCEN*SIN(PHI))**2)) c WRITE(LOUT,*) 'rtest',rtest c WRITE(LLOG,*) 'rtest',rtest IF (rtest .le. router .and. rtest .ge. rinner) THEN nintersect=nintersect+1 rintersect(nintersect)=rtest phiintersect(nintersect)=phiedge1 ENDIF ENDIF C phiedge2 PHI=ABS(PHICEN-phiedge2) C Ensure that the angle between the track center and c phiedge2 is acute IF (PHI.GT.PI) PHI=2*PI-PHI C Ensure that the track crosses phiedge2 IF (RADIUS**2-((LCEN*SIN(PHI))**2).GE.0) THEN rtest = LCEN*COS(PHI) + SQRT(RADIUS**2-((LCEN*SIN(PHI))**2)) c WRITE(LOUT,*) 'rtest',rtest c WRITE(LLOG,*) 'rtest',rtest IF (rtest .le. router .and. rtest .ge. rinner) THEN nintersect=nintersect+1 rintersect(nintersect)=rtest phiintersect(nintersect)=phiedge2 ENDIF rtest = LCEN*COS(PHI) - SQRT(RADIUS**2-((LCEN*SIN(PHI))**2)) c WRITE(LOUT,*) 'rtest',rtest c WRITE(LLOG,*) 'rtest',rtest IF (rtest .le. router .and. rtest .ge. rinner) THEN nintersect=nintersect+1 rintersect(nintersect)=rtest phiintersect(nintersect)=phiedge2 ENDIF ENDIF c -------- C Check for intersections at the inner and outer edges of the cell c -------- c To find the phi of the track at a given radius, we again use the c polar coordinates of the track's center of curvature and the c law of cosines c rinner arg=((rinner**2)+(LCEN**2)-(RADIUS**2))/(2*rinner*LCEN) IF (ABS(arg).LE.1) THEN PHI=ACOS(arg) phitest = PHICEN + PHI IF (phitest.GT.PI) phitest=phitest-2*PI IF (phitest.LT.-PI) phitest=phitest+2*PI c WRITE(LOUT,*) 'phitest',phitest*(180/PI) c WRITE(LLOG,*) 'phitest',phitest*(180/PI) IF (((phitest .ge. phiedge1) .and. & (phitest .le. phiedge2)) & .or. & ((phiedge1 .gt. 0) .and. (phiedge2 .lt. 0) .and. & ((phitest .ge. phiedge1) .or. (phitest .le. phiedge2))))THEN nintersect = nintersect+1 IF (nintersect .le. 4) THEN rintersect(nintersect)=rinner phiintersect(nintersect)=phitest ENDIF ENDIF phitest = PHICEN - PHI IF (phitest.GT.PI) phitest=phitest-2*PI IF (phitest.LT.-PI) phitest=phitest+2*PI c WRITE(LOUT,*) 'phitest',phitest*(180/PI) c WRITE(LLOG,*) 'phitest',phitest*(180/PI) IF (((phitest .ge. phiedge1) .and. & (phitest .le. phiedge2)) & .or. & ((phiedge1 .gt. 0) .and. (phiedge2 .lt. 0) .and. & ((phitest .ge. phiedge1) .or. (phitest .le. phiedge2))))THEN nintersect = nintersect+1 IF (nintersect .le. 4) THEN rintersect(nintersect)=rinner phiintersect(nintersect)=phitest ENDIF ENDIF ENDIF c router arg=((router**2)+(LCEN**2)-(RADIUS**2))/(2*router*LCEN) IF (ABS(arg).LE.1)THEN PHI=ACOS(arg) phitest = PHICEN + PHI IF (phitest.GT.PI) phitest=phitest-2*PI IF (phitest.LT.-PI) phitest=phitest+2*PI c WRITE(LOUT,*) 'phitest',phitest*(180/PI) c WRITE(LLOG,*) 'phitest',phitest*(180/PI) IF (phitest.GT.PI) phitest=2*PI-phitest IF (phitest.lt.-PI) phitest=phitest+2*PI IF (((phitest .ge. phiedge1) .and. & (phitest .le. phiedge2)) & .or. & ((phiedge1 .gt. 0) .and. (phiedge2 .lt. 0) .and. & ((phitest .ge. phiedge1) .or. (phitest .le. phiedge2))))THEN nintersect = nintersect+1 IF (nintersect .le. 4) THEN rintersect(nintersect)=router phiintersect(nintersect)=phitest ENDIF ENDIF phitest = PHICEN - PHI IF (phitest.GT.PI) phitest=phitest-2*PI IF (phitest.LT.-PI) phitest=phitest+2*PI c WRITE(LOUT,*) 'phitest',phitest*(180/PI) c WRITE(LLOG,*) 'phitest',phitest*(180/PI) IF (phitest.GT.PI) phitest=2*PI-phitest IF (phitest.lt.-PI) phitest=phitest+2*PI IF (((phitest .ge. phiedge1) .and. & (phitest .le. phiedge2)) & .or. & ((phiedge1 .gt. 0) .and. (phiedge2 .lt. 0) .and. & ((phitest .ge. phiedge1) .or. (phitest .le. phiedge2))))THEN nintersect = nintersect+1 IF (nintersect .le. 4) THEN rintersect(nintersect)=router phiintersect(nintersect)=phitest ENDIF ENDIF ENDIF c ------- c If there is a track through the cell, Calculate the pathlength c ------- C the number of intersections in a cell is 2 or 4 c (It can be 4 when the pt. at which the track is tangent to the c cell boundary is just outside that boundary) c We find the angle of the intersection pts in the coord system c centered on the track's center of curvature (angleint(1-4)). c The pathlength in xy for 2 intersection pts. is then given by c pathlenxy = RADIUS * anglespan c where anglespan is the angular diff between intersection pts. C for 4 intersection pts. we add the pathlength in xy between the c first two angles and the last two angles. c The total path length is therefore given by c pathlen = pathlenxy/sin(phiz) c where phiz is the angle the track makes to the z axis DO I=1,nintersect angleint(I)=ACOS(((rintersect(I)**2)-(RADIUS**2)-(LCEN**2))/ & (2*RADIUS*LCEN)) IF (phiintersect(I).LT.PHICEN) angleint(I)=-angleint(I) c WRITE(LOUT,*) 'rintersect',rintersect(I) c WRITE(LOUT,*) 'phiintersect',phiintersect(I)*(180/PI) c WRITE(LOUT,*) 'angleint',angleint(I)*(180/PI) c WRITE(LLOG,*) 'rintersect',rintersect(I) c WRITE(LLOG,*) 'phiintersect',phiintersect(I)*(180/PI) c WRITE(LLOG,*) 'angleint',angleint(I)*(180/PI) ENDDO IF (nintersect.EQ.2) THEN LPATH = .TRUE. anglespan = ABS(angleint(1)-angleint(2)) pathlenxy = RADIUS * anglespan PATHLEN = pathlenxy/SQRT(1-COS3_D(ITRACK)**2) ELSEIF (nintersect.EQ.4) THEN LPATH=.TRUE. angle(1)=-PI angle(2)=-PI angle(3)=-PI angle(4)=-PI DO I=1,4 IF (angleint(I).GT.angle(1)) angle(1)=angleint(I) ENDDO DO I=1,4 IF (angleint(I).NE.angle(1)) THEN IF (angleint(I).GT.angle(2)) angle(2)=angleint(I) ENDIF ENDDO DO I=1,4 IF ((angleint(I).NE.angle(1)).AND. & (angleint(I).NE.angle(2))) THEN IF (angleint(I).GT.angle(3)) angle(3)=angleint(I) ENDIF ENDDO DO I=1,4 IF ((angleint(I).NE.angle(1)).AND. & (angleint(I).NE.angle(2)).AND. & (angleint(I).NE.angle(3))) THEN angle(4)=angleint(I) ENDIF ENDDO anglespan = ABS(angle(1)-angle(2))+ABS(angle(3)-angle(4)) pathlenxy = RADIUS * anglespan PATHLEN = pathlenxy/SQRT(1-COS3_D(ITRACK)**2) ELSEIF (nintersect.NE.0) THEN WRITE(LOUT,*) 'Funky number of intersection points'// & ' between track and cell boudaries. ' WRITE(LOUT,*) 'nintersect=',nintersect WRITE(LOUT,*) 'event',NEVT WRITE(LOUT,*) 'Superlayer, layer, cell=',ISLAY,ILAY,ICELL WRITE(LLOG,*) 'Funky number of intersection points'// & ' between track and cell boudaries. ' WRITE(LLOG,*) 'nintersect=',nintersect WRITE(LLOG,*) 'event',NEVT WRITE(LLOG,*) 'Superlayer, layer, cell=',ISLAY,ILAY,ICELL write(lout,*) 'rinner',rinner write(lout,*) 'router',router write(lout,*) 'phiedge1',phiedge1*(180/PI) write(lout,*) 'phiedge2',phiedge2*(180/PI) write(llog,*) 'rinner',rinner write(llog,*) 'router',router write(llog,*) 'phiedge1',phiedge1*(180/PI) write(llog,*) 'phiedge2',phiedge2*(180/PI) write(lout,*) 'lcen',lcen write(lout,*) 'phicen',phicen*(180/pi) write(lout,*) 'radius',radius write(llog,*) 'lcen',lcen write(llog,*) 'phicen',phicen*(180/pi) write(llog,*) 'radius',radius WRITE (LOUT,*) 'rintersect #1=',rintersect(1) WRITE (LOUT,*) 'phiintersect #1=',phiintersect(1)*(180/PI) WRITE (LOUT,*) 'rintersect #2=',rintersect(2) WRITE (LOUT,*) 'phiintersect #2=',phiintersect(2)*(180/PI) WRITE (LOUT,*) 'rintersect #3=',rintersect(3) WRITE (LOUT,*) 'phiintersect #3=',phiintersect(3)*(180/PI) WRITE (LLOG,*) 'rintersect #1=',rintersect(1) WRITE (LLOG,*) 'phiintersect #1=',phiintersect(1)*(180/PI) WRITE (LLOG,*) 'rintersect #2=',rintersect(2) WRITE (LLOG,*) 'phiintersect #2=',phiintersect(2)*(180/PI) WRITE (LLOG,*) 'rintersect #3=',rintersect(3) WRITE (LLOG,*) 'phiintersect #3=',phiintersect(3)*(180/PI) GOTO 400 ENDIF c write(lout,*) 'nintersect', nintersect c write(llog,*) 'nintersect', nintersect c WRITE (LOUT,*) 'rintersect #1=',rintersect(1) c WRITE (LOUT,*) 'phiintersect #1=',phiintersect(1)*(180/PI) c WRITE (LOUT,*) 'rintersect #2=',rintersect(2) c WRITE (LOUT,*) 'phiintersect #2=',phiintersect(2)*(180/PI) c WRITE (LLOG,*) 'rintersect #1=',rintersect(1) c WRITE (LLOG,*) 'phiintersect #1=',phiintersect(1)*(180/PI) c WRITE (LLOG,*) 'rintersect #2=',rintersect(2) c WRITE (LLOG,*) 'phiintersect #2=',phiintersect(2)*(180/PI) c WRITE (LOUT,*) 'PATHLENXY', PATHLENXY c WRITE (LOUT,*) 'PATHLEN', PATHLEN c WRITE (LLOG,*) 'PATHLENXY', PATHLENXY c WRITE (LLOG,*) 'PATHLEN', PATHLEN IF (COS3_D(ITRACK) .EQ. 0.) THEN WRITE (LOUT,*) 'COS3_D',COS3_D(ITRACK), & 'event',NEVT, & 'Superlayer, layer, cell=',ISLAY,ILAY,ICELL WRITE (LLOG,*) 'COS3_D',COS3_D(ITRACK), & 'event',NEVT, & 'Superlayer, layer, cell=',ISLAY,ILAY,ICELL ENDIF 400 CONTINUE RETURN END