*$ CREATE MGDRAW.FOR *COPY MGDRAW * * *=== mgdraw ===========================================================* * * SUBROUTINE MGDRAW ( ICODE, MREG ) INCLUDE '(DBLPRC)' INCLUDE '(DIMPAR)' INCLUDE '(IOUNIT)' * *----------------------------------------------------------------------* * * * Copyright (C) 1990-2013 by Alfredo Ferrari * * All Rights Reserved. * * * * * * MaGnetic field trajectory DRAWing: actually this entry manages * * all trajectory dumping for * * drawing * * * * Created on 01 March 1990 by Alfredo Ferrari * * INFN - Milan * * Last change 12-Nov-13 by Alfredo Ferrari * * INFN - Milan * * * *----------------------------------------------------------------------* * modified Nov 2019 by fluka course teachers INCLUDE '(CASLIM)' INCLUDE '(COMPUT)' INCLUDE '(SOURCM)' INCLUDE '(FHEAVY)' INCLUDE '(FLKSTK)' INCLUDE '(GENSTK)' INCLUDE '(MGDDCM)' INCLUDE '(PAPROP)' INCLUDE '(QUEMGD)' INCLUDE '(SUMCOU)' INCLUDE '(TRACKR)' * DIMENSION DTQUEN ( MXTRCK, MAXQMG ) * CHARACTER*20 FILNAM LOGICAL LFCOPE SAVE LFCOPE DATA LFCOPE / .FALSE. / **** Additional Variables for this exercise: * dimension of arrays == max number of produced photons in had/neutron int * above thresh PARAMETER ( MXPHOT=100 ) * How many? INTEGER NPHOT * How many entering? INTEGER NENTER * deposited energy of entering photons DIMENSION EDPHOT (MXPHOT) * parent of entering photons DIMENSION IPARENT(MXPHOT) COMMON /FOLLOW/ THRESH, NPHOT,ICINT ** local, no need to declare, here for clarity: total edep in detector DOUBLE PRECISION EDTOT SAVE /FOLLOW/ SAVE EDTOT SAVE EDPHOT * region number of detector SAVE IRGDET SAVE NENTER SAVE IPARENT *****!!!!!! remember to initialize at each new event: call sodraw * **** End Additional Variables * *----------------------------------------------------------------------* * * * Icode = 1: call from Kaskad * * Icode = 2: call from Emfsco * * Icode = 3: call from Kasneu * * Icode = 4: call from Kashea * * Icode = 5: call from Kasoph * * * *----------------------------------------------------------------------* * * IF ( .NOT. LFCOPE ) THEN LFCOPE = .TRUE. IF ( KOMPUT .EQ. 2 ) THEN FILNAM = '/'//CFDRAW(1:8)//' DUMP A' ELSE FILNAM = CFDRAW END IF OPEN ( UNIT = IODRAW, FILE = FILNAM, STATUS = 'NEW', FORM = & 'FORMATTED' ) END IF * standard output: we do not need it here. Leave for reference * WRITE (IODRAW) NTRACK, MTRACK, JTRACK, SNGL (ETRACK), * & SNGL (WTRACK) * WRITE (IODRAW) ( SNGL (XTRACK (I)), SNGL (YTRACK (I)), * & SNGL (ZTRACK (I)), I = 0, NTRACK ), * & ( SNGL (DTRACK (I)), I = 1, MTRACK ), * & SNGL (CTRACK) * +-------------------------------------------------------------------* * | Quenching is activated IF ( LQEMGD ) THEN IF ( MTRACK .GT. 0 ) THEN RULLL = ZERZER CALL QUENMG ( ICODE, MREG, RULLL, DTQUEN ) * standard output: we do not need it here. Leave for reference * WRITE (IODRAW) ( ( SNGL (DTQUEN (I,JBK)), I = 1, MTRACK ), * & JBK = 1, NQEMGD ) END IF END IF * | End of quenching * +-------------------------------------------------------------------* * +-------------------------------------------------------------------* * | Store en dep in detec IF ( MREG .EQ. IRGDET .AND. MTRACK .GT. 0) THEN DE = ZERZER DO I =1, NTRACK DE = DE + DTRACK (I) END DO IPHOTON = LLOUSE IF (IPHOTON .GT. 0 ) THEN EDPHOT (IPHOTON) = EDPHOT (IPHOTON) + DE END IF EDTOT = EDTOT + DE END IF * | * +-------------------------------------------------------------------* RETURN * *======================================================================* * * * Boundary-(X)crossing DRAWing: * * * * Icode = 1x: call from Kaskad * * 19: boundary crossing * * Icode = 2x: call from Emfsco * * 29: boundary crossing * * Icode = 3x: call from Kasneu * * 39: boundary crossing * * Icode = 4x: call from Kashea * * 49: boundary crossing * * Icode = 5x: call from Kasoph * * 59: boundary crossing * * * *======================================================================* * * ENTRY BXDRAW ( ICODE, MREG, NEWREG, XSCO, YSCO, ZSCO ) * * +-------------------------------------------------------------------* * | do something if new region is detector, a photon is entering * | and it is above thr. IF ( NEWREG .EQ. IRGDET .AND. JTRACK .EQ. 7 .AND. & ETRACK .GT. THRESH) THEN NENTER = NENTER +1 LLOUSE = NENTER IF ( NENTER .GT. MXPHOT ) THEN WRITE (LUNOUT,*) ' TOO MANY ENTERING PHOTONS' CALL FLABRT('BXDRAW','NENTER') END IF IPARENT (NENTER ) = ISPUSR(5) * | ** energy, position,direction,generation number, time of arrival WRITE (IODRAW,'(A,1x,I3,1p,6(1x,E12.3))') & ' Photon entering ', LTRACK, ETRACK, XSCO, YSCO, ZSCO, & ATRACK*1.0D+09, WTRACK * | Target Z ispusr(1) * | Projectile Z ispusr(2) * | Projectile A ispusr(3) * | Generation level of original photon ispusr(4) * | Ordinal number of original photon ispusr(5) * | Track number of original photon ispusr(6) * | Region of first conversion (if any)! ispusr(7) * | Code of interaction generating orig. phot. ispusr(8) * | Id of projectile generating orig. phot. ispusr(9) * | Region of interaction generating orig. phot. ispusr(10) * | Ordinal number of entered photon llouse * | * | Position of first had interaction (3) spausr(1:3) * | Energy of first had projectile spausr(4) * | Position of gamma production point (3) spausr(5:7) * | Energy of last had projectile spausr(8) * | Energy of original photon (might scatter) spausr(9) * +-------------------------------------------------------------------* * | Who's the parent? WRITE (IODRAW, '( A, I3, 1p, E12.3, A)') & ' From original phot ', ISPUSR(5), SPAUSR(9), ' GeV' WRITE (IODRAW,'(A,I4, A, 3(1x,I3),A,1x,I3,1p,A,E12.3,A)') & ' Generated in ', ISPUSR(8),' by ',ISPUSR(9), ISPUSR(2), & ISPUSR(3), ' on ', ISPUSR(1), ' at ', SPAUSR(8),' GeV' WRITE (IODRAW,'(A,I4,A,1p,3(1x,E12.3))')' in region', & ISPUSR(10),' position ',SPAUSR(5), SPAUSR(6), SPAUSR(7) * | +----------------------------------------------------------------* * | | it was a secondary hadronic interaction IF ( SPAUSR(4) .GT. ANGLGB ) THEN WRITE (IODRAW,*) 'Generation level' ,ISPUSR(4) WRITE (IODRAW,'(A,1p,4(1x,E12.3))') & ' Energy and position at first interaction', & SPAUSR(4), SPAUSR(1), SPAUSR(2), SPAUSR(3) END IF * | | * | +----------------------------------------------------------------* * | +----------------------------------------------------------------* * | | the photon was rescattered IF ( ISPUSR (7) .NE. 0 ) THEN WRITE (IODRAW,*) ' First conversion', ISPUSR (7) END IF * | | * | +----------------------------------------------------------------* * | END IF * | * +-------------------------------------------------------------------* RETURN * *======================================================================* * * * Event End DRAWing: * * * *======================================================================* * * ENTRY EEDRAW ( ICODE ) WRITE (IODRAW , *)'total photons in this event: ',NPHOT WRITE (IODRAW , *)'total entering in this event: ',NENTER * +-------------------------------------------------------------------* * | write energy only if energy inside ( Anglgb = machine "zero") IF ( EDTOT .GT. ANGLGB ) THEN WRITE (IODRAW, *)'Detected energy ',EDTOT IF ( NENTER .GT. 0 ) THEN WRITE (IODRAW,*) ' from each entering photon ' WRITE (IODRAW,*) (EDPHOT(I),I=1,NENTER) EPHTOT = ZERZER DO J=1,NENTER EPHTOT=EPHTOT + EDPHOT(J) ENDDO WRITE (IODRAW,*)'Total energy from above-threshold photons', & EPHTOT END IF END IF * | * +-------------------------------------------------------------------* RETURN * *======================================================================* * * * ENergy deposition DRAWing: * * * * Icode = 1x: call from Kaskad * * 10: elastic interaction recoil * * 11: inelastic interaction recoil * * 12: stopping particle * * 13: pseudo-neutron deposition * * 14: escape * * 15: time kill * * Icode = 2x: call from Emfsco * * 20: local energy deposition (i.e. photoelectric) * * 21: below threshold, iarg=1 * * 22: below threshold, iarg=2 * * 23: escape * * 24: time kill * * Icode = 3x: call from Kasneu * * 30: target recoil * * 31: below threshold * * 32: escape * * 33: time kill * * Icode = 4x: call from Kashea * * 40: escape * * 41: time kill * * 42: delta ray stack overflow * * Icode = 5x: call from Kasoph * * 50: optical photon absorption * * 51: escape * * 52: time kill * * * *======================================================================* * * ENTRY ENDRAW ( ICODE, MREG, RULL, XSCO, YSCO, ZSCO ) IF ( .NOT. LFCOPE ) THEN LFCOPE = .TRUE. IF ( KOMPUT .EQ. 2 ) THEN FILNAM = '/'//CFDRAW(1:8)//' DUMP A' ELSE FILNAM = CFDRAW END IF OPEN ( UNIT = IODRAW, FILE = FILNAM, STATUS = 'NEW', FORM = & 'FORMATTED' ) END IF * +-------------------------------------------------------------------* * | Store en dep in detec IF ( MREG .EQ. IRGDET ) THEN IPHOTON = LLOUSE IF (IPHOTON .GT. 0 ) THEN EDPHOT (IPHOTON) = EDPHOT (IPHOTON) + RULL END IF EDTOT = EDTOT + RULL END IF * | * +-------------------------------------------------------------------* * +-------------------------------------------------------------------* * | Quenching is activated : calculate quenching factor * | and store quenched energy in DTQUEN(1, jbk) IF ( LQEMGD ) THEN RULLL = RULL CALL QUENMG ( ICODE, MREG, RULLL, DTQUEN ) END IF * | end quenching * +-------------------------------------------------------------------* RETURN * *======================================================================* * * * SOurce particle DRAWing: * * * *======================================================================* * ENTRY SODRAW IF ( .NOT. LFCOPE ) THEN LFCOPE = .TRUE. IF ( KOMPUT .EQ. 2 ) THEN FILNAM = '/'//CFDRAW(1:8)//' DUMP A' ELSE FILNAM = CFDRAW END IF OPEN ( UNIT = IODRAW, FILE = FILNAM, STATUS = 'NEW', FORM = & 'FORMATTED' ) END IF * Threshold for entering/produced photons THRESH = 1.0D-04 * write on our dump file that there is a new event. * write energy and type of particle WRITE (IODRAW, *)'#####' * +-------------------------------------------------------------------* * | (Radioactive) isotope: it works only for 1 source particle on * | the stack for the time being IF ( ILOFLK (NPFLKA) .GE. 100000 .AND. LRADDC (NPFLKA) ) THEN IARES = MOD ( ILOFLK (NPFLKA), 100000 ) / 100 IZRES = MOD ( ILOFLK (NPFLKA), 10000000 ) / 100000 IISRES = ILOFLK (NPFLKA) / 10000000 IONID = ILOFLK (NPFLKA) * in this case tkeflk is half life with a minus sign WRITE (IODRAW, *) 'New Event ', & NCASE, IONID, IARES, IZRES, -TKEFLK(NPFLKA) * | * +-------------------------------------------------------------------* * | Patch for heavy ions: it works only for 1 source particle on * | the stack for the time being ELSE IF ( ABS (ILOFLK (NPFLKA)) .GE. 10000 ) THEN IONID = ILOFLK (NPFLKA) CALL DCDION ( IONID ) WRITE (IODRAW, *) 'New Event ', & NCASE, IONID, IBHEAV(-IONID), ICHEAV(-IONID), & TKEFLK(NPFLKA) * | * +-------------------------------------------------------------------* * | Patch for heavy ions: ??? ELSE IF ( ILOFLK (NPFLKA) .LT. -6 ) THEN IONID = ILOFLK(NPFLKA) WRITE (IODRAW, *) 'New Event ', & NCASE, IONID, IBHEAV(-IONID), ICHEAV(-IONID), & TKEFLK(NPFLKA) * | * +-------------------------------------------------------------------* * | ELSE IONID = ILOFLK(NPFLKA) WRITE (IODRAW, *) 'New Event ', & NCASE, IONID, IBARCH(IONID), ICHRGE(IONID), & TKEFLK(NPFLKA) END IF * | * +-------------------------------------------------------------------* ******initialise our variables at each event NPHOT = 0 EDTOT = ZERZER NENTER = 0 DO I =1, MXPHOT EDPHOT(I) = ZERZER END DO ***** initialization of detector region at first event IF ( NCASE .EQ. 1 ) THEN CALL GEON2R( "DETEC ", IRGDET, IERR ) IF ( IERR .NE. 0 ) THEN WRITE (LUNOUT, *) 'Detector region not found' CALL FLABRT ('mgdraw','no detec') ELSE WRITE (LUNOUT, *) 'Detector region ', IRGDET END IF END IF RETURN * *======================================================================* * * * USer dependent DRAWing: * * * * Icode = 10x: call from Kaskad * * 100: elastic interaction secondaries * * 101: inelastic interaction secondaries * * 102: particle decay secondaries * * 103: delta ray generation secondaries * * 104: pair production secondaries * * 105: bremsstrahlung secondaries * * 110: radioactive decay products * * Icode = 20x: call from Emfsco * * 208: bremsstrahlung secondaries * * 210: Moller secondaries * * 212: Bhabha secondaries * * 214: in-flight annihilation secondaries * * 215: annihilation at rest secondaries * * 217: pair production secondaries * * 219: Compton scattering secondaries * * 221: photoelectric secondaries * * 225: Rayleigh scattering secondaries * * 237: mu pair production secondaries * * Icode = 30x: call from Kasneu * * 300: interaction secondaries * * Icode = 40x: call from Kashea * * 400: delta ray generation secondaries * * For all interactions secondaries are put on GENSTK common (kp=1,np) * * but for KASHEA delta ray generation where only the secondary elec- * * tron is present and stacked on FLKSTK common for kp=npflka * * * *======================================================================* * ENTRY USDRAW ( ICODE, MREG, XSCO, YSCO, ZSCO ) IF ( .NOT. LFCOPE ) THEN LFCOPE = .TRUE. IF ( KOMPUT .EQ. 2 ) THEN FILNAM = '/'//CFDRAW(1:8)//' DUMP A' ELSE FILNAM = CFDRAW END IF OPEN ( UNIT = IODRAW, FILE = FILNAM, STATUS = 'NEW', FORM = & 'FORMATTED' ) END IF * Pass the code to stuprf ICINT = ICODE * +-------------------------------------------------------------------* * | First inelastic interaction: select by generation level * | and store relevant variables IF ( LTRACK .EQ. 1 .AND. ICODE .EQ. 101 ) THEN * | Position SPAUSR (1) = XSCO SPAUSR (2) = YSCO SPAUSR (3) = ZSCO * | +----------------------------------------------------------------* * | | Kinetic energy: Beware to ions IF ( JTRACK .GE. -6 ) THEN EKGENP = ETRACK - AM (JTRACK) IDGENP = JTRACK ELSE EKGENP = ETRACK - AMNHEA (-JTRACK) IDGENP = -2 END IF SPAUSR (4) = EKGENP * | | * | +----------------------------------------------------------------* END IF * | * +-------------------------------------------------------------------* * * When a photon is produced in inelastic, * we'll have to set and update other variables, * however some of them will have to be associated * to each single photon, and to photons only. Not here ==> stuprf * * Record photon interaction: must do in stupre, bacause it is called * before usdraw RETURN *=== End of subrutine Mgdraw ==========================================* END