Re: FLUKA: Access to models ?
Hi Laurent,
Please find attached the mgdraw which I used to compute the damage to the
NLC positron production target. It dumps recoil momenta of residual
nuclei, heavy fragments and low-energy neutrons. I carefully compute the
energy balance to make sure that PXRES,PYRES,PZRES,ERES,IBRES,ICRES,AMMRES
are consistent with the present event.
Please let me know if you find mistakes in this routine as compared to the
one Alfredo is going to send you.
Stefan
___________________________________
Stefan Roesler
SLAC
P.O. Box 4349, MS 48
Stanford, CA 94309
Tel: 650-926-2048
Fax: 650-926-3569
E-mail: sroesler@slac.stanford.edu
Stefan.Roesler@cern.ch
On Fri, 6 Apr 2001, Laurent APHECETCHE wrote:
> Alfredo Ferrari wrote:
> >
> > Hi all
> >
> > the excitation energy CANNOT be accessed at mgdraw level. TVCMS if
> > everything is going correctly should come out =0 (no residual excitation
> > energy). The excitation energy after PEANUT cannot be accessed and it
> > would not be physically meaningful, since it depends a lot on the
> > termination conditions and thresholds of the preequilibrium part (where
> > the code assumes equilibrium is reached). The excitation energy after the
> > cascade part which perhaps was the one Laurent was referring to, is also
> > of course strongly dependent on the boundary between cascade and
> > prequilibrium and does not bear any real physics meaning.
> >
> > Also Ekres,Px,y,zres are not necessarily updated correctly: their status
> > depends on the physics input cards set by the user and if
>fission/fragmentation
> > occurred etc. DON'T rely on any of these variables.
> >
> > Ciao
> > Alfredo
>
> Hi Alfredo,
>
> Thanks for the reply, even if it seems discouraging...
> What I would like to have is the distribution of nuclei recoil momenta,
> in order to be able to compute damage energies. Is there a way to do
> that, e.g. some combination of input cards insuring that Ekres,Px,y,zres
> _are_ meaningful, or is it hopeless ?
>
> Best regards,
>
> --
> Dr. Laurent APHECETCHE (mailto:aphecetc@in2p3.fr) (IN2P3-CNRS)
> SUBATECH-EMN-4 rue Alfred Kastler-BP 20722-44307 NANTES cedex 03
> TEL (+33/0) 2 51 85 84 17 - FAX (+33/0) 2 51 85 84 24 (France)
> Collaborations PHENIX http://www.phenix.bnl.gov/~aphecetc et MEGAPIE.
>
*$ CREATE MGDRAW.FOR
*COPY MGDRAW
* *
*=== mgdraw ===========================================================*
* *
SUBROUTINE MGDRAW ( ICODE, MREG )
INCLUDE '(DBLPRC)'
INCLUDE '(DIMPAR)'
INCLUDE '(IOUNIT)'
*
*----------------------------------------------------------------------*
* *
* MaGnetic field trajectory DRAWing: actually this entry manages *
* all trajectory dumping for *
* drawing *
* *
* Created by Alfredo Ferrari *
* INFN - Milan *
* last change 19-oct-98 by Alfredo Ferrari *
* INFN - Milan *
* *
*----------------------------------------------------------------------*
*
INCLUDE '(CASLIM)'
INCLUDE '(COMPUT)'
INCLUDE '(EPISOR)'
INCLUDE '(FINUC)'
INCLUDE '(MGDDCM)'
INCLUDE '(PAPROP)'
INCLUDE '(STACK)'
INCLUDE '(STARS)'
INCLUDE '(TRACKR)'
*
INCLUDE '(RESNUC)'
INCLUDE '(FHEAVY)'
*
CHARACTER*20 FILNAM
LOGICAL LFCOPE,LFIRST
SAVE LFCOPE,LFIRST
DATA LFCOPE,LFIRST / .FALSE.,.TRUE. /
*
*----------------------------------------------------------------------*
* *
* 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 =
& 'UNFORMATTED' )
END IF
C WRITE (IODRAW) NTRACK, MTRACK, JTRACK, SNGL (ETRACK),
C & SNGL (WTRACK)
C WRITE (IODRAW) ( SNGL (XTRACK (I)), SNGL (YTRACK (I)),
C & SNGL (ZTRACK (I)), I = 0, NTRACK ),
C & ( SNGL (DTRACK (I)), I = 1,MTRACK ),
C & SNGL (CTRACK)
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 *
* Icode = 2x: call from Emfsco *
* 20: local energy deposition (i.e. photoelectric) *
* 21: below threshold, iarg=1 *
* 22: below threshold, iarg=2 *
* 23: escape *
* Icode = 3x: call from Kasneu *
* 30: target recoil *
* 31: below threshold *
* 32: escape *
* Icode = 4x: call from Kashea *
* 40: escape *
* Icode = 5x: call from Kasoph *
* 50: optical photon absorption *
* 51: escape *
* *
*----------------------------------------------------------------------*
* *
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 =
& 'UNFORMATTED' )
END IF
C WRITE (IODRAW) 0, ICODE, JTRACK, SNGL (ETRACK), SNGL (WTRACK)
C WRITE (IODRAW) SNGL (XSCO), SNGL (YSCO), SNGL (ZSCO), SNGL (RULL)
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 =
& 'UNFORMATTED' )
END IF
C WRITE (IODRAW) -NCASE, LSTACK, LSTMAX, SNGL (TKESUM),
C & SNGL (WEIPRI)
C WRITE (IODRAW) ( ILO (I), SNGL (TKE(I)+AM(ILO(I))), SNGL (WT(I)),
C & SNGL (XA(I)), SNGL (YA(I)), SNGL (ZA(I)),
C & SNGL (TX(I)), SNGL (TY(I)), SNGL (TZ(I)),
C & I = 1, LSTACK )
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 *
* 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 *
* 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 FINUC common (kp=1,np) *
* but for KASHEA delta ray generation where only the secondary elec- *
* tron is present and stacked on STACK common for kp=lstack *
* *
*======================================================================*
*
ENTRY USDRAW ( ICODE, MREG, XSCO, YSCO, ZSCO )
*
IF (LFIRST) THEN
OPEN(50,FILE='recoil.dbg',STATUS='UNKNOWN')
OPEN(51,FILE='recoil.out',STATUS='UNKNOWN')
OPEN(52,FILE='neutron.out',STATUS='UNKNOWN')
OPEN(53,FILE='heavy.out',STATUS='UNKNOWN')
LFIRST = .FALSE.
ENDIF
*
* zero the momentum-balance accumulators
PXBAL = ZERZER
PYBAL = ZERZER
PZBAL = ZERZER
*
* select elastic and inelastic interactions
IF ((ICODE.LE.101).OR.((ICODE.EQ.300).AND.(ISPUSR(1).EQ.999)))
& THEN
IF ((ICODE.EQ.300).AND.(ISPUSR(1).EQ.999)) ISPUSR(1) = 0
c WRITE(50,1000) INT(WEIPRI)
1000 FORMAT(/,1X,'Beam particle #',I7)
c WRITE(50,1001) ICODE,XSCO,YSCO,ZSCO
1001 FORMAT(1X,'Interaction code ',I5,' Coordinates ',3F8.4)
*
* incoming particle
PXBAL = PXBAL+CXTRCK*PTRACK
PYBAL = PYBAL+CYTRCK*PTRACK
PZBAL = PZBAL+CZTRCK*PTRACK
c WRITE(50,1002) JTRACK,WTRACK,ETRACK-AM(JTRACK),
c & CXTRCK*PTRACK,CYTRCK*PTRACK,CZTRCK*PTRACK
1002 FORMAT(1X,'Incoming particle ',I4,1P,2E12.4,2X,3F8.4)
*
* just to make sure that no incoming neutron in the low-energy module
* is above threshold
IF ((ICODE.EQ.300).AND.((ETRACK-AM(JTRACK)).GE.1.9600E-02))
& WRITE(50,*) ' WARNING ! neutron above threshold !'
*
* outgoing particles
DO 10 IP=1,NP
PXBAL = PXBAL-CXR(IP)*PLR(IP)
PYBAL = PYBAL-CYR(IP)*PLR(IP)
PZBAL = PZBAL-CZR(IP)*PLR(IP)
IF (IP.EQ.1) THEN
c WRITE(50,1003) KPART(IP),WEI(IP),TKI(IP),
c & CXR(IP)*PLR(IP),CYR(IP)*PLR(IP),CZR(IP)*PLR(IP)
1003 FORMAT(1X,'Outgoing particles',I4,1P,2E12.4,2X,3F8.4)
ELSE
c WRITE(50,1004) KPART(IP),WEI(IP),TKI(IP),
c & CXR(IP)*PLR(IP),CYR(IP)*PLR(IP),CZR(IP)*PLR(IP)
1004 FORMAT(19X,I4,1P,2E12.4,2X,3F8.4)
ENDIF
*
* flag neutrons produced by the high-energy module
IF ((ICODE.LE.101).AND.(KPART(IP).EQ.8)) THEN
ISPUSR(1) = 999
* dump neutrons below 19.6MeV
IF (TKI(IP).LT.1.9600E-02) THEN
WRITE(52,1011)
& INT(WEIPRI),XSCO,YSCO,ZSCO,WEI(IP),TKI(IP),
& CXR(IP)*PLR(IP),CYR(IP)*PLR(IP),CZR(IP)*PLR(IP)
1011 FORMAT(I6,3F9.4,1P,5E12.4)
ENDIF
ENDIF
10 CONTINUE
*
* outgoing heavy particles
DO 11 IP=1,NPHEAV
PXBAL = PXBAL-CXHEAV(IP)*PHEAVY(IP)
PYBAL = PYBAL-CYHEAV(IP)*PHEAVY(IP)
PZBAL = PZBAL-CZHEAV(IP)*PHEAVY(IP)
IF (IP.EQ.1) THEN
c WRITE(50,1005) KHEAVY(IP),WHEAVY(IP),TKHEAV(IP),
c & CXHEAV(IP)*PHEAVY(IP),CYHEAV(IP)*PHEAVY(IP),
c & CZHEAV(IP)*PHEAVY(IP)
1005 FORMAT(1X,'Outgoing heavies ',I4,2E12.4,2X,3F8.4)
ELSE
c WRITE(50,1004) KHEAVY(IP),WHEAVY(IP),TKHEAV(IP),
c & CXHEAV(IP)*PHEAVY(IP),CYHEAV(IP)*PHEAVY(IP),
c & CZHEAV(IP)*PHEAVY(IP)
ENDIF
* dump heavies
IF (ICODE.LE.101) THEN
WRITE(53,1012)
& INT(WEIPRI),XSCO,YSCO,ZSCO,WHEAVY(IP),TKHEAV(IP),
& CXHEAV(IP)*PHEAVY(IP),CYHEAV(IP)*PHEAVY(IP),
& CZHEAV(IP)*PHEAVY(IP),ICHEAV(KHEAVY(IP)),
& IBHEAV(KHEAVY(IP))
1012 FORMAT(I6,3F9.4,1P,5E12.4,2I4)
ENDIF
11 CONTINUE
*
* momentum balance
c WRITE(50,1006) PXBAL,PYBAL,PZBAL
1006 FORMAT(1X,'Momentum balance',28X,3F8.4)
*
* patch for elastic interactions:
* eres is not updated for some reasons in fluka, so we have to
* calculate it here
IF (ICODE.EQ.100) THEN
BRES = DBLE(IBRES)
ZRES = DBLE(ICRES)
XERES = IBRES*AMUGEV+EMVGEV*ENERGY(BRES,ZRES)+EKRES
ELSE
XERES = ERES
ENDIF
*
* dump residual nucleus
IF ((ICODE.LE.101).AND.(IBRES.GT.0)) THEN
WRITE(51,1012)
& INT(WEIPRI),XSCO,YSCO,ZSCO,WTRACK,EKRES,
& PXRES,PYRES,PZRES,ICRES,IBRES
ENDIF
*
* residual nucleus from common resnuc
c WRITE(50,1007) ICRES,IBRES,AMMRES,
c & XERES,EKRES,PXRES,PYRES,PZRES
1007 FORMAT(1X,'Residual nucleus ',2I4,F10.4,
& /,25X,F8.4,E12.4,3F8.4)
*
* residual nucleus: mass check
AMCHEK = SQRT(ABS(XERES**2-PXRES**2-PYRES**2-PZRES**2))
c IF (ABS(AMMRES-AMCHEK).GT.1.D-2)
c & WRITE(50,1008) AMCHEK,AMMRES
1008 FORMAT(1X,'---> Residual nucleus: inconsistent mass ',2E12.4)
*
* residual nucleus: kinetic energy check
C IF (ABS(EKRES-TVRECL).GT.1.D-6) WRITE(50,1009) EKRES,TVRECL
C1009 FORMAT(1X,'---> Residual nucleus: inconsistent energy ',2E12.4)
*
* momentum balance check
PXBAL = PXBAL-PXRES
PYBAL = PYBAL-PYRES
PZBAL = PZBAL-PZRES
c IF ((ABS(PXBAL)+ABS(PYBAL)+ABS(PZBAL)).GT.5.D-2)
c & WRITE(50,1010) ICODE,PXBAL,PYBAL,PZBAL
1010 FORMAT(1X,'---> Inconsistent momentum balance',3X,I4,3X,3F8.4)
*
ENDIF
*
RETURN
*=== End of subrutine Mgdraw ==========================================*
END
Partial thread listing:
- Re: FLUKA: Access to models ?, (continued)
FLUKA: new (only bug fixes) release,
Alfredo Ferrari
FLUKA: Re: FLUKA
Stefan Roesler