SUBROUTINE CALC_PATHLEN_B0(ICELL,ILAY,ISLAY,ITRACK,LPATH,PATHLEN) c This routine is called for cells which are crossed by a straigth track. c If the pathlength through the cell is successfully calculated, LPATH c is returned as .TRUE., and PATHLEN returns the pathlength. 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 PATHLEN, pathlenpolar REAL rinner, router, phiedge1, phiedge2, phioffset REAL PHIPCA, PHITRACK, LOWERLIM, UPPERLIM, DIFFLIM REAL rintersect(2),phiintersect(2) REAL xintersect(2),yintersect(2), xintdiff, yintdiff REAL*4 a,b REAL*4 rtest,x,y,phitest 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 parameters are found in UTCBLOCK. 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 Fetch the track data c ------ C DIST0_D(ITRACK) is the distance of closest approach to the origin PHITRACK=ATAN2(COS2_D(ITRACK),COS1_D(ITRACK)) c PHITRACK is the angle of the track w.r.t the x-axis c in the x-y plane C Note that ATAN2 will yield PHITRACK to be between -PI and PI PHIPCA=ATAN2(Y0_D(ITRACK),X0_D(ITRACK)) C PHIPCA is the angle of the pt. of closest approach w.r.t. the origin c in polar coordinates the equation of the track w.r.t the origin is c r = b / (sin(phi)-a*cos(phi)) C Make sure the slope isn't infinite IF (PHITRACK .EQ. PI/2) GOTO 400 c Find the equation paramaters a=tan(PHITRACK) b=DIST0_D(ITRACK)*(sin(PHIPCA)-a*cos(PHIPCA)) IF (PHIPCA .GT. 0) b = ABS(b) IF (PHIPCA .LT. 0) b = -ABS(b) c The sign of the offset is the same as the sign of PHIPCA C If phipca is 0 then phitrack was pi/2, and was not considered c Alternatively SIGN(b,PHIPCA) could have been used c Express the equation of the track w.r.t the centre of the SL b = b + a*TRAN_W_X(ISLAY) - TRAN_W_Y(ISLAY) c ------- C Check for intersections at the side edges of the cell c ------- c Avoid division by zero or out of range error by requiring radius c to be within the drift chamber IF ((sin(phiedge1)-a*cos(phiedge1))/b & .gt. 1/ROUTL(NSLAY_MAX)) THEN rtest = b / (sin(phiedge1)-a*cos(phiedge1)) 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 IF ((sin(phiedge2)-a*cos(phiedge2))/b & .gt. 1/ROUTL(NSLAY_MAX)) THEN rtest = b / (sin(phiedge2)-a*cos(phiedge2)) 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 solve for the c x and y coordinates of the line through the track at that c radius, then use ATAN2 to obtain phi (-pi to pi). Since the c line through the track crosses the radius twice, there are two c points. We choose the one with a phi that falls between the c angles spanned by the cell. x= (-a*b + SQRT(-b**2 + (1+a**2)*rinner**2))/(1+a**2) y= a*x+b phitest = ATAN2(y,x) 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. 2) THEN rintersect(nintersect)=rinner phiintersect(nintersect)=phitest ENDIF ENDIF x= (-a*b - SQRT(-b**2 + (1+a**2)*rinner**2))/(1+a**2) y= a*x+b phitest = ATAN2(y,x) 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. 2) THEN rintersect(nintersect)=rinner phiintersect(nintersect)=phitest ENDIF ENDIF x= (-a*b + SQRT(-b**2 + (1+a**2)*router**2))/(1+a**2) y= a*x+b phitest = ATAN2(y,x) 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. 2) THEN rintersect(nintersect)=router phiintersect(nintersect)=phitest ENDIF ENDIF x= (-a*b - SQRT(-b**2 + (1+a**2)*router**2))/(1+a**2) y= a*x+b phitest = ATAN2(y,x) 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. 2) THEN rintersect(nintersect)=router phiintersect(nintersect)=phitest ENDIF ENDIF c ------- c If there is a track through the cell, Calculate the pathlength c ------- C the number of intersections in a cell should be 0 or 2 C The path length is given in polar coord by c pathlenpolar = sqrt(r2^2 + r1^2 - r1*r2*cos(phi1)*cos(phi2)) c The total path length is therefore given by c pathlen = pathlenpolar/sin(phiz) c where phiz is the angle the track makes to the z axis c write(lout,*) 'nintersect', nintersect c write(llog,*) 'nintersect', nintersect IF (nintersect .eq. 2) THEN LPATH = .TRUE. c pathlenpolar = SQRT(rintersect(1)**2+rintersect(2)**2 c & -2*rintersect(1) *rintersect(2) c & *COS(phiintersect(1))*COS(phiintersect(2))) xintersect(1)=rintersect(1)*Cos(phiintersect(1)) yintersect(1)=rintersect(1)*Sin(phiintersect(1)) xintersect(2)=rintersect(2)*Cos(phiintersect(2)) yintersect(2)=rintersect(2)*Sin(phiintersect(2)) xintdiff = xintersect(2)-xintersect(1) yintdiff = yintersect(2)-yintersect(1) pathlenpolar = SQRT(xintdiff**2 + yintdiff**2) PATHLEN = pathlenpolar/SQRT(1-COS3_D(ITRACK)**2) 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,*) 'PATHLENPOLAR', PATHLENPOLAR c WRITE (LOUT,*) 'PATHLEN', PATHLEN c WRITE (LLOG,*) 'PATHLENPOLAR', PATHLENPOLAR 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 ELSEIF (nintersect.ne.0) THEN WRITE(LOUT,*) 'Funky number of intersection points'// & ' between track and cell boudaries. '// & 'nintersect=',nintersect,'', & 'event',NEVT,'', & 'Superlayer, layer, cell=',ISLAY,ILAY,ICELL WRITE(LLOG,*) 'Funky number of intersection points'// & ' between track and cell boudaries. '// & 'nintersect=',nintersect,'', & 'event',NEVT,'', & 'Superlayer, layer, cell=',ISLAY,ILAY,ICELL c CALL KERROR(-3,0,'DPLOT',MSG) GOTO 400 ENDIF 400 CONTINUE RETURN END