%L %E !INDENT M 4; "INDENT MORTRAN LISTING 4 PER NESTING LEVEL" !INDENT F 2; "INDENT FORTRAN OUTPUT 2 PER NESTING LEVEL" "This line is 80 characters long, use it to set up the screen width" "23456789|123456789|123456789|123456789|123456789|123456789|123456789|123456789" "*********************************************************************" " " " ******************** " " * * " " * COMPTON_PET.MOR * " " * * " " ******************** " " " " " " An EGS4 user code which records all interactions in the animal " " PET instrument consisting Si and BGO detectors when a 511 keV " " annihilation photons is incident isotropically " " " " " "*********************************************************************" "---------------------------------------------------------------------" "STEP 1: USER-OVERRIDE-OF-EGS4-MACROS " "---------------------------------------------------------------------" REPLACE {$MXMED} WITH {3} "only 3 medium in the problem(default 10)" REPLACE {$MXREG} WITH {9} "only 9 geometric regions (default 2000)" REPLACE {$MXSTACK} WITH {15}"less than 15 particles on stack at once" REPLACE {$MXCYLS} WITH {4} REPLACE {;OUTPUT61#;#;} WITH {;{SETR A =@LG} WRITE(6,{COPY A}){P1};{COPY A}FORMAT{P2};} REPLACE {;OUTPUT69#;#;} WITH {;{SETR A =@LG} WRITE(6,{COPY A}){P1};WRITE(9,{COPY A}){P1};{COPY A}FORMAT{P2};} REPLACE {;COMIN/RANDOM/;} WITH {;COMMON/RANDOM/IXX;} REPLACE {$RANDOMSET#;} WITH {IXX=IXX*663608941;{P1}=0.5 + IXX*0.23283064E-09;} REPLACE {;COMIN/SCORE/;} WITH {;COMMON/SCORE/IWATCH;} REPLACE {;COMIN/DSOURCE/;} WITH {;COMMON/DSOURCE/DISK_R(6),DISK_X(6),DISK_Y(6),DISK_Z(6), A_DISK(6),F_DISK(6),SUM_F(7);} REPLACE {;COMIN/STACK/;} WITH {;COMMON/STACK/$LGN(E,X,Y,Z,U,V,W,DNEAR,WT, IQ,IR,LATCH($MXSTACK)),LATCHI,NP; $ENERGY PRECISION E;} REPLACE {$TRANSFER PROPERTIES TO # FROM#;} WITH {X{P1}=X{P2};Y{P1}=Y{P2};Z{P1}=Z{P2};IR{P1}=IR{P2};WT{P1}=WT{P2}; DNEAR{P1}=DNEAR{P2};LATCH{P1}=LATCH{P2};} ;COMIN/RANDOM,BOUNDS,EPCONT,MEDIA,MISC,SCORE,STACK,PLADTA,THRESH, CYLDTA,ELECIN,USER,USEFUL,DSOURCE/; DATA PI /3.1415926536/; %E "---------------------------------------------------------------------" "STEP 2 PRE-HATCH-CALL-INITIALIZATION " "---------------------------------------------------------------------" $TYPE MEDARR(24,3) /$S'SI', 22*' ', $S'AIR',21*' ', $S'BGO',21*' '/; NMED=3; NREG=9; DO J=1,NMED [DO I=1,24 [MEDIA(I,J)=MEDARR(I,J);]] DO K=1,NREG [MED(K)=2;] MED(2)=1; MED(5)=3; DO I=1,NREG [ECUT(I)=0.01;PCUT(I)=0.01;] DO I=1,NREG [IRAYLR(I)=1;] " TURN ON RAYLEIGH SCATTERING IN THE SLAB " "---------------------------------------------------------------------" "STEP 3 HATCH-CALL " "---------------------------------------------------------------------" CALL HATCH; "---------------------------------------------------------------------" "STEP 4 INITIALIZATION-FOR-HOWFAR " "---------------------------------------------------------------------" "1st PLANE" PCOORD(1,1)= 0.0; PCOORD(2,1)= 0.0; PCOORD(3,1) =-2.0; PNORM(1,1) = 0.0; PNORM(2,1) = 0.0; PNORM(3,1) = 1.0; PCOORD(1,2)= 0.0; PCOORD(2,2)= 0.0; PCOORD(3,2) = 2.0; PNORM(1,2) = 0.0; PNORM(2,2) = 0.0; PNORM(3,2) = 1.0; PCOORD(1,3)= 0.0; PCOORD(2,3)= 0.0; PCOORD(3,3) =-8.0; PNORM(1,3) = 0.0; PNORM(2,3) = 0.0; PNORM(3,3) = 1.0; PCOORD(1,4)= 0.0; PCOORD(2,4)= 0.0; PCOORD(3,4) = 8.0; PNORM(1,4) = 0.0; PNORM(2,4) = 0.0; PNORM(3,4) = 1.0; CYRAD2(1)= 4.00; CYRAD2(2)= 12.96; CYRAD2(3)= 77.44; CYRAD2(4)= 116.64; "---------------------------------------------------------------------" "STEP 5 INITIALIZATION-FOR-AUSGAB " "---------------------------------------------------------------------" /IAUSFL(18),IAUSFL(24)/=1; IWATCH = 1; CALL WATCH(-99,IWATCH); "---------------------------------------------------------------------" "STEP 6 DETERMINATION-OF-INICIDENT-PARTICLE-PARAMETERS " "---------------------------------------------------------------------" IQIN=0;" INCIDENT CHARGE - ELECTRONS" EIN=0.511;" 511 KEV KINETIC ENERGY" /XIN,YIN,ZIN/=0.0;" INCIDENT AT ORIGIN" /UIN,VIN,WIN/=0.0;" MOVING ALONG Z AXIS" IRIN=1;" STARTS IN REGION 1" WTIN=1.0;" WEIGHT = 1 SINCE NO VARIANCE REDUCTION USED" LATCHI=0;" LATCH SET TO ZERO AT START OF EACH HISTORY" %E "---------------------------------------------------------------------" "STEP 7 SHOWER-CALL " "---------------------------------------------------------------------" NSHOWER=0; NESCAPE=0; NCASE = 100000; " <--- TOTAL NUMBER OF 511 keV PHOTONS " NPAIR = NCASE/2; OUTPUT69 NPAIR; (I13); DO I=1,NCASE [ "**************************************************************" DISK_S = 0; "FOR DISK SOURCE" IF(DISK_S.EQ.1) [ CALL DISK_SOURCE(XIN,YIN,ZIN) ; ] "**************************************************************" IWATCH=1; IF(MOD(I,2)=1) [ $RANDOMSET RTHE; ANGLE1=RTHE*PI; $RANDOMSET RFAI; ANGLE2=RFAI*PI*2.0; UIN=SIN(ANGLE1)*COS(ANGLE2); VIN=SIN(ANGLE1)*SIN(ANGLE2); IF (ANGLE1<(PI/2.0)) [ WIN=SQRT(1-UIN**2.0-VIN**2.0);] ELSE [ WIN=-SQRT(1-UIN**2.0-VIN**2.0);] ] ELSEIF(MOD(I,2)=0) [ UIN=UIN*(-1); VIN=VIN*(-1); WIN=WIN*(-1);] "=== CHECK BIASED SOURCE ==== " " ANGLE_START=PI/2.0-ATAN(2.24/2.0); " " ANGLE_END =PI/2.0+ATAN(2.24/2.0); " " " " BIASOK = 0; " " IF((ANGLE_START <= ANGLE1) & (ANGLE1 <= ANGLE_END))[ " " BIASOK = 1 ;] " " " " IF(BIASOK=0)[ " " NESCAPE=NESCAPE+1;] " BIASOK = 1; IF(BIASOK=1)[ CALL WATCH(-1,IWATCH); CALL SHOWER(IQIN,EIN,XIN,YIN,ZIN,UIN,VIN,WIN,IRIN,WTIN); CALL WATCH(99,IWATCH); NSHOWER = NSHOWER+1 ;] ] ;OPEN(9,FILE='/tmp/comp_pet.dat',STATUS='OLD'); OUTPUT69 ;(23X,' 0 0 0 0 0 0 0 0'); "---------------------------------------------------------------------" "STEP 8 OUTPUT-OF-RESULTS " "---------------------------------------------------------------------" STOP;END;"END OF TUTOR5 MAIN ROUTINE" %E "**********************************************************************" " " SUBROUTINE AUSGAB(IARG); " " "**********************************************************************" ;COMIN/SCORE,STACK/; CALL WATCH(IARG,IWATCH); RETURN;END;"END OF AUSGAB" %E "*********************************************************************" " " SUBROUTINE HOWFAR; " " "*********************************************************************" COMIN/EPCONT,PLADTA,CYLDTA,STACK/; IRL=IR(NP); "SET LOCAL VARIABLE" IF(IRL.EQ.7) [IDISC=1; "TERMINATE THIS HISTORY"] IF(IRL.EQ.8) [IDISC=1;] IF(IRL.EQ.9) [IDISC=1;] ELSEIF(IRL.EQ.5) ["WE ARE IN THE BGO DETECTOR" $PLAN2P(3,7,-1,4,8,1); $CYLNDR(3,0,IHIT,TCYL); IF(IHIT.EQ.1) [ $FINVAL(TCYL,XF,YF,ZF); IF(ZF.LE.-2) [IRNXT=4;] ELSEIF(ZF.GE.2) [IRNXT=6;] ELSE [IRNXT=3;] $CHGTR(TCYL,IRNXT); ] ELSEIF(IHIT.EQ.0) [ $CYLNDR(4,1,IHIT,TCYL); IF(IHIT.EQ.1) [$CHGTR(TCYL,9);] ] ] ELSEIF(IRL.EQ.3) ["WE ARE IN THE AIR BTW SI AND BGO" $PLAN2P(1,4,-1,2,6,1); $CYL2(2,2,3,5); ] ELSEIF(IRL.EQ.2) ["WE ARE IN THE SI DETECTOR" $PLAN2P(1,4,-1,2,6,1); $CYL2(1,1,2,3); ] ELSEIF(IRL.EQ.1) ["WEHN WE IN THE AIR BTW TWO SI DETECTORS" $CYLNDR(1,1,IHIT,TCYL); IF(IHIT.EQ.1) [$CHGTR(TCYL,2);] $PLAN2P(1,4,-1,2,6,1); ] ELSEIF(IRL.EQ.4) ["WE ARE IN THE AIR LEFT OF SI" $CYLNDR(3,1,IHIT,TCYL); IF(IHIT.EQ.1) [$CHGTR(TCYL,5);] $PLANE1(1,1,IHIT,TPLN); IF(IHIT.EQ.1) [ $FINVAL(TPLN,XF,YF,ZF); IF((XF*XF + YF*YF).GE.CYRAD2(2)) [IRNXT=3;] ELSEIF((XF*XF + YF*YF).LE.CYRAD2(1)) [IRNXT=1;] ELSE [IRNXT=2;] $CHGTR(TPLN,IRNXT); ] ELSEIF(IHIT.EQ.0) [ $PLANE1(3,-1,IHIT,TPLN); IF(IHIT.EQ.1) [$CHGTR(TPLN,7);] ] ] ELSEIF(IRL.EQ.6) ["WE ARE IN THE AIR RIGHT OF SI" $CYLNDR(3,1,IHIT,TCYL); IF(IHIT.EQ.1) [$CHGTR(TCYL,5);] $PLANE1(2,-1,IHIT,TPLN); IF(IHIT.EQ.1) [ $FINVAL(TPLN,XF,YF,ZF); IF((XF*XF + YF*YF).GE.CYRAD2(2)) [IRNXT=3; ] ELSEIF((XF*XF + YF*YF).LE.CYRAD2(1)) [IRNXT=1;] ELSE [IRNXT=2;] $CHGTR(TPLN,IRNXT); ] ELSEIF(IHIT.EQ.0) [ $PLANE1(4,1,IHIT,TPLN); IF(IHIT.EQ.1) [$CHGTR(TPLN,8);] ] ] RETURN; END;" END OF SUBROUTINE HOWFAR" "**********************************************************************" " " " WATCH " " " ;SUBROUTINE WATCH(IARG,IWATCH); " " " A GENERAL PURPOSE AUXILIARY ROUTINE FOR USE WITH THE EGS4 SYSTEM " " IT PRINTS OUT INFORMATION ABOUT THE PARTICLE TRANSPORT " " " " FOR IWATCH = 1 IT PRINTS INFORMATION ABOUT EACH DISCRETE INTERACTION " " FOR IWATCH = 2 OR 3 IT PRINTS INFORMATION ABOUT EACH STEP AS WELL " " FOR IWATCH = 4 IT PRINTS GRAPHING DATA " " " " " " ROUTINE IS USED VIA TWO MANDATORY AND 1 OPTIONAL CALL FROM THE USER'S " " CODE " " " " 1)THE ROUTINE MUST BE INITIALIZED BY A CALL WITH IARG=-99 BEFORE THE FIRST" " CALL TO SHOWER " " 2)THE ROUTINE MUST BE CALLED NEAR THE BEGINNING OF THE AUSGAB SUBROUTINE " " IF (IWATCH > 0 ) CALL WATCH(IARG,IWATCH); " " 3)THE ROUTINE MAY BE CALLED AT THE END OF EACH HISTORY WITH IARG = - 1 SO " " A MESSAGE WILL GET PRINTED " " " " " " THE ROUTINE USES UP TO 132 COLUMNS FOR OUTPUT AND WILL INITIALIZE " " AN PRINTER LINEPRITER VIA A CALL TO SUBROUTINNE PRINTER. IF NOT AVAILABLE" " REPLACE THE CALL TO PRNTER WTH: OUTPUT;('1'); " " " " JAN 1984 GENERALIZED VERSION WITH INITIALIZATION " " DAVE ROGERS NRCC " " JUN 1987 PUT IN IWATCH = 4 OPTION AFB " " JUL 1988 COMPATIBLE WITH X-RAY FLUORESCENCE DWOR " " SEP 1990 ADDED ENERGY OUTPUT TO IWATCH = 4 OPTION AFB " " OCT 1990 UNIX compatible carriage control DWOR " " "*****************************************************************************" "DEFINE A LOCAL MACRO" REAL A; DIMENSION A(100,8); REPLACE {;OUTPUT69#;#;} WITH {;{SETR A =@LG} WRITE(9,{COPY A}){P1};{COPY A}FORMAT{P2};} REPLACE {;OUTPUT61#;#;} WITH {;{SETR A =@LG} WRITE(6,{COPY A}){P1};{COPY A}FORMAT{P2};} REPLACE {$CNTOUT(#);(#);} WITH { ICOUNT=ICOUNT+1; IF( ((0 <= X({P1})**2) & (0 <= Y({P1})**2 )) &( 0 <= Z({P1})**2 ) ) [ IF(IQ({P1}).EQ.0) [ KK=KK+1; A(KK,1) = NGAMMA; A(KK,2) = IR({P1}); A(KK,3) = {P1}; A(KK,4) = {P2}; A(KK,5) = EKE; A(KK,6) = X({P1}); A(KK,7) = Y({P1}); A(KK,8) = Z({P1}); "OUTPUT69 " "A(KK,1),A(KK,2),A(KK,3),A(KK,4),A(KK,5),A(KK,6),A(KK,7),A(KK,8);" "(F3.1,2X,F3.1,2X,F3.1,2X,F5.1,F7.3,1X,3(1X,1PE12.5));" ]] } ;COMIN/STACK,EPCONT/; ;OPEN(9,FILE='/tmp/comp_pet.dat'); DATA ICOUNT/0/,JHSTRY/1/; IF(MOD(JHSTRY,2)=1) [ NGAMMA = 2;] ELSE [ NGAMMA = 1;] IF(IARG = -99) [ "INITIALIZE FLAGS SO WE WILL GET CALLS THRU AUSGAB" DO J=1,25[IAUSFL(J)=1;]; /IAUSFL(6),IAUSFL(22),IAUSFL(23),IAUSFL(24)/=0; IF(IWATCH = 4) IAUSFL(6)= 0; ] IF(IARG = -1) ["MAIN IS ASSUMED TO CALL AUSGAB WITH IARG=-1 AT END OF HISTORY" IF(IWATCH = 4) [ WRITE(13,:GRAPHICS_FORMAT:) 0,0,0,0.0,0.0,0.0,0.0,JHSTRY; JHSTRY=JHSTRY+1; ] ELSE[ IF(IWATCH=10)[ JHSTRY=JHSTRY+1;ICOUNT=ICOUNT+2;RETURN;] ELSE[ IF(NGAMMA.EQ.2) [ N_PAIRS = (JHSTRY+1)/2 ; OUTPUT69 N_PAIRS ;(23X,' 0 0 0 0 0 0 0 ', I12); KK = 0; ] JHSTRY=JHSTRY+1;ICOUNT=ICOUNT+2;RETURN;] ] ] IF((IARG.EQ.99).AND.(NGAMMA.EQ.2)) [ DO II=1, KK [OUTPUT69 A(II,1),A(II,2),A(II,3),A(II,4),A(II,5),A(II,6),A(II,7),A(II,8); (F3.1,2X,F3.1,2X,F3.1,2X,F5.1,F7.3,1X,3(1X,1PE12.5)); ] ] IF( (IWATCH.NE.4).AND.((ICOUNT.EQ.0).OR.(IARG.EQ.-99)) )[ "PRINT HEADER" ICOUNT=1;] IF((IWATCH = 4) & (IARG >= 0)) [ "GRAPHICS OUTPUT" WRITE(13,:GRAPHICS_FORMAT:) NP,IQ(NP),IR(NP),X(NP),Y(NP),Z(NP),E(NP); :GRAPHICS_FORMAT:FORMAT(3I4,4G15.8,I12);] IF(IARG = 5 | IARG < 0) RETURN; IF(IWATCH = 4) RETURN; "NONE OF THE REST NEEDED FOR GRAPHICS OUTPUT" EKE=E(NP);IF(IQ(NP).NE.0)[EKE=E(NP)-0.511;] IF(IARG = 0 & IWATCH = 2)[$CNTOUT(NP);(0);] ELSEIF(IARG = 0)[RETURN;] IF(IARG = 1)[$CNTOUT(NP);(1);] ELSEIF(IARG = 2)[$CNTOUT(NP);(2);] ELSEIF(IARG = 3)[ IF((NP=1 & EKE = 0.511) & (IQ(1)=0 & IR(1)=1)) [IWATCH = 10;] ELSE[$CNTOUT(NP);(3);] ] ELSEIF(IARG = 4)[$CNTOUT(NP);(4);] " ELSEIF(IARG = 6)[$CNTOUT(NP);(6);]" " ELSEIF(IARG = 7)[$CNTOUT(NP);(7);]" ELSEIF(IARG = 8)[$CNTOUT(NP);(8);] ELSEIF(IARG = 9)[$CNTOUT(NP);(9);] ELSEIF(IARG = 10)[$CNTOUT(NP);(10);] ELSEIF(IARG = 11)[$CNTOUT(NP);(11);] ELSEIF(IARG = 12)[$CNTOUT(NP);(12);] ELSEIF(IARG = 13)[$CNTOUT(NP);(13);] ELSEIF(IARG = 14)[$CNTOUT(NP);(14);] ELSEIF(IARG = 15)[$CNTOUT(NP);(15);] ELSEIF(IARG = 16)[$CNTOUT(NP);(16);] " ELSEIF(IARG = 17)[$CNTOUT(NP);(17);] " ELSEIF(IARG = 18)[$CNTOUT(NP);(18);] " ELSEIF(IARG = 19)[$CNTOUT(NP);(19);]" " ELSEIF(IARG = 20)[ " " if photon on top of stack, it must be the fluorescent photon" " it may have 0 energy if incident photon was below binding energy" " IF(IQ(NP).EQ.0)[$CNTOUT(NP);(201);] " " ELSE [$CNTOUT(NP);(202);] " " ] " ELSEIF(IARG = 24)[$CNTOUT(NP);(24);] IF(IARG = 0 & IWATCH = 2)[OUTPUT61 TUSTEP,VSTEP,TVSTEP,EDEP; (T11,'TUSTEP,VSTEP,TVSTEP,EDEP',T36,':',4(1PE15.7));ICOUNT=ICOUNT+1; ] IF (NP=1 | IARG=0) [RETURN;] IF(IARG = 7 | IARG = 9 | IARG = 11 | IARG = 13 | IARG = 14 | IARG = 16 | IARG = 18 | (IARG = 20 & IQ(NP)=0 & E(NP) ~= 0) | IARG <= 3)[ "IARG condition was without brackets before and IQ(NP)=0" "it was wrong since IQ(NP)=0 occured when there was no real" "event (below binding) as well as when there was a fluorescent" "photon. This condition means a fluorescent" "photon was created and above it is the ejected electron" N=NP-1;EKE= E(N); IF(IQ(N).NE.0) EKE=E(N)-0.511; IF(IARG <= 3) [ ;] " [$CNTOUT(N);(100);]" ELSEIF(IARG = 20)[$CNTOUT(N);(203);] ELSE [$CNTOUT(N);(204);] ] RETURN;END; %I4 %C80 %Q1 %E "**************************************************************************" " REAL FUNCTION TAN(X); " " " GIVEN X IN RADIANS THIS CALCULATES THE TANGENT OF X " "**************************************************************************" A=COS(X); IF(ABS(A).GT.1.E-10) [TAN = SIN(X)/A;] ELSE [TAN = 1.E10;] RETURN;END; "**************************************************************************" " " " " ;SUBROUTINE DISK_SOURCE(XIN,YIN,ZIN); " " " " " " "**************************************************************************" COMIN/DSOURCE,RANDOM/; DATA PI /3.1415926536/; NDISK = 6; DO L=1,NDISK [DISK_Z(L)=0;] DISK_R(1)= 0.5; DISK_X(1)= 1.0; DISK_Y(1)= 0.0; DISK_R(2)= 0.05; DISK_X(2)= 0.5; DISK_Y(2)= SIN(PI/3); DISK_R(3)= 0.4; DISK_X(3)= -0.5; DISK_Y(3)= SIN(PI/3); DISK_R(4)= 0.2; DISK_X(4)= -1.0; DISK_Y(4)= 0.0; DISK_R(5)= 0.3; DISK_X(5)= -0.5; DISK_Y(5)= SIN(PI/3)*(-1); DISK_R(6)= 0.1; DISK_X(6)= 0.5; DISK_Y(6)= SIN(PI/3)*(-1); DO M=1,NDISK [ A_DISK(M) = PI*DISK_R(M)*DISK_R(M);] TOT_DISK = 0.0; DO I=1,NDISK [ TOT_DISK = TOT_DISK + A_DISK(I);] DO J=1,NDISK [ F_DISK(J) = A_DISK(J)/TOT_DISK;] SUM_F(1) = 0.0; DO K=1,NDISK [ SUM_F(K+1) = SUM_F(K) + F_DISK(K);] $RANDOMSET POS_DISK ; DO L=1,NDISK [ IF( SUM_F(L).LE.POS_DISK .AND. POS_DISK.LE.SUM_F(L+1) ) [ DISK_NUM = L ; ] ] $RANDOMSET RAND_R ; POS_R = RAND_R * DISK_R(DISK_NUM) ; $RANDOMSET RAND_THETA ; POS_THETA = RAND_THETA * 2.0 * PI; XIN = DISK_X(DISK_NUM) + POS_R * COS(POS_THETA) ; YIN = DISK_Y(DISK_NUM) + POS_R * SIN(POS_THETA) ; ZIN = 0.0; RETURN;END;