From: Paolo MAESTRO (paolo.maestro@pi.infn.it)
Date: Mon May 10 2004 - 17:40:41 CEST
Hi ,
I'm a beginner with FLUKA , and I'm trying to understand some basic
things with a simple example.
I defined a block of tungsten 20 radiation lengths thick,
and I want to calculate the total energy ETOT_CAL deposited by a 20 GeV
electron
entering the calorimeter and record it into a histogram.
I modified the mgdraw.f routine so to calculate ETOT_CAL as the sum
of two different contributions:
1. a loop over DTRACK(J) J=1,...Mtrack i.e. the energy deposition along the tracks
2. the RULL variable (local energy depositions) in the entry ENDRAW
I send you the code in attachment.
Doing this way , it seems that the output ETOT_CAL has the value I
expect, i.e. about 19.4 GeV since the shower is almost completely
contained.
Then I tested the method changing both the beam
energy and the particle type, and it seems it works well.
Anyhow, since I came to my procedure by several empirical attempts ,
I wonder if that is correct and if there are alternative ways to do the
same thing at run time, I mean without analyzing the output files
produced by fluka.
I also put in attachment the input file I used.
Thanks
Paolo Maestro
TITLE
Test of simple calorimeter
*
* Sets the standards defaults of FLUKA (see user manual)
*
DEFAULTS CALORIME
*
*
*
*23456789 123456789 123456789 123456789 123456789 123456789 123456789 123456789
*BEAM 20. RAY
*BEAM 20. ELECTRON
*BEAM 20. PROTON
***
*** Heavy ion case: BEAM WHAT(1) is kin energy (momentum) per nucleon NOT total
***
BEAM 20. HEAVYION
HI-PROPE 5.0 10.0 0.0
EVENTYPE 0. 0. 2. 0. 0. 0. DPMJET
*
* Particles must injected in vacuum or material
* not in black hole
* direction +z
BEAMPOS 0.0 0.0 -50.0
SOURCE
********************
*
* Switching off Electromagnetic Showering
*EMF EMF-OFF
*
* Do not transport electrons, positrons, photons, all charged and
* neutral hadrons, neutrinos
*
*DISCARD 3. 4. 7. 218. 223.
*DISCARD 5. 6. 27. 28.
*********************
*
* Geometry description
*
GEOBEGIN 0.01 COMBINAT
water layer
*123 12345 123456789 123456789 123456789 123456789 123456789 123456789
SPH 1 0.0 0.0 0.0 +5000.0
RPP 2 -500.0 +500.0 -500.0 +500.0 -500.0 +500.0
RPP 3 -25.0 +25.0 -25.0 +25.0 -7.0 0.0
RPP 4 -0.75 +0.75 -0.75 +0.75 9.965 10.0
* XYP 3 -3.5
* XYP 4 +3.5
END
* black hole
1 5 1 -2
* vacuum
2 5 2 -3 -4
* cal block
3 5 3
* si pixel
4 5 4
END
GEOEND
*
* Here defines materials and Compounds of sea water (Antares 1)
*
* 1) Hydrogen
MATERIAL 1.0 1.0079 .0000899 3.0 1.0 HYDROGEN
* 2) Oxygen
MATERIAL 8.0 15.999 0.001429 8.0 OXYGEN
* 3) Magnesium
MATERIAL 12.0 24.305 1.738 9.0 MAGNESIU
* 4) Potassium
MATERIAL 19.0 39.102 0.031165 11.0 POTASSIU
* 5) Calcium
MATERIAL 20.0 40.08 1.54 12.0 CALCIUM
* 6) Sodium
MATERIAL 11.0 22.990 0.038349 13.0 SODIUM
* 7) Chlorine
MATERIAL 17.0 35.4529 0.0029947 14.0 CHLORINE
* 8) Sulphur
MATERIAL 16.0 32.066 2.070 15.0 SULFUR
* 8) Tungsten
MATERIAL 74.0 183.84 19.3 23.0 TUNGSTEN
* 8) Carbon
MATERIAL 6.0 12.011 2.265 6.0 CARBON
* 9) Silicon
MATERIAL 14.0 28.085 2.33 26.0 SILICON
* 9) Sea water compound (Antares 1) using atom relative contents
MATERIAL 1.0341 18.0 SEAWATER
COMPOUND 2.00000 3.0 1.008800 8.0 0.00943 13.0 SEAWATER
COMPOUND 0.000209 11.0 0.001087 9.0 0.000209 12.0 SEAWATER
COMPOUND 0.01106 14.0 0.005820 15.0 SEAWATER
**
*23456789 123456789 123456789 123456789 123456789 123456789 123456789 123456789
EMFCUT 0.0001 0.0001 0. 3. 3.
EMFCUT 0.00001 0.00001 0. 4. 4.
*
* Here assigns materials to regions
* tungsten in region 3
ASSIGNMAT 23.0 3.0
* External Black Hole in region 1
ASSIGNMAT 1.0 1.0
* Vacuum at beginning (region2)
ASSIGNMAT 2.0 2.0
* Silicon at end (region 4)
ASSIGNMAT 26.0 4.0
*
* Setting SCORE and Ouptut levels.
*
SCORE 208.0 211.0
OUTLEVEL 1.0 7.0 1.
*EVENTDAT 11.0 PIPPO
*
* Activates Boundary Crossing Scoring. In this example this
* is useful since it enables the Entry GEOBDX in GEODEN routine.
* There ntupla of survived muons is filled
*
*23456789 123456789 123456789 123456789 123456789 123456789 123456789 123456789
USERDUMP 101.0 0. 1.0 SEAMU
*
* Activates Random Number initialization
*
RANDOMIZE 1.
*
* Activates user dependent initialization
*
USRICALL
*
* Starts Job asking for a given number of events
*
START 5. 99999999. 0.0 0.0
*
* Activates user dependent output
*
USROCALL
STOP
*$ 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 10-nov-02 by Alfredo Ferrari *
* INFN - Milan *
* *
*----------------------------------------------------------------------*
*
INCLUDE '(CASLIM)'
INCLUDE '(COMPUT)'
INCLUDE '(EPISOR)'
INCLUDE '(FHEAVY)'
INCLUDE '(FINUC)'
INCLUDE '(MGDDCM)'
INCLUDE '(PAPROP)'
INCLUDE '(STACK)'
INCLUDE '(STARS)'
INCLUDE '(TRACKR)'
INCLUDE 'MYVAR.inc'
*
CHARACTER*20 FILNAM
LOGICAL LFCOPE
SAVE LFCOPE
DATA LFCOPE / .FALSE. /
PARAMETER (NNTUPLE=8)
REAL*4 XTUP(nntuple)
*
*----------------------------------------------------------------------*
* *
* 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
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)
*******************
IF(MREG.EQ.3) THEN
WRITE (*,*) 'TRACK ', JTRACK,' ',ICODE,' ',MREG, ' '
DO I = 1,MTRACK
ETOT_CAL = ETOT_CAL + DTRACK(I)
WRITE (*,*) I,' ',DTRACK(I)
ENDDO
ENDIF
IF(MREG.EQ.4) THEN
WRITE (*,*) 'TRACK ', JTRACK,' ',ICODE,' ',MREG, ' '
DO I = 1,MTRACK
ETOT_SCD = ETOT_SCD + DTRACK(I)
WRITE (*,*) I,' ',DTRACK(I)
ENDDO
ENDIF
********************
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 )
C-STA
WRITE(*,*) 'BXDRAW ',icode,mreg,newreg,jtrack
$ ,etrack,am(jtrack)
IF(MREG.EQ.3.AND.NEWREG.EQ.4) THEN ! Select the desired boundary
*IF( JTRACK.EQ.10.OR.JTRACK.EQ.11 ) THEN ! Select muons
IF(ETRACK.GT.AM(JTRACK)) THEN ! Muon has survived
XTUP(1) = JTRACK
XTUP(2) = ETRACK-AM(JTRACK)
XTUP(3) = XSCO
XTUP(4) = YSCO
XTUP(5) = ZSCO
XTUP(6) = CXTRCK
XTUP(7) = CYTRCK
XTUP(8) = CZTRCK
write(*,*) "Fortran says : ",nusrbx
* CALL treefill(JTRACK,
* & (ETRACK-AM(JTRACK)),XSCO,YSCO,ZSCO,
* & CXTRCK,CYTRCK,CZTRCK) ! <---- Function Call Added
* CALL HFN(100,XTUP)
ENDIF
* ENDIF
ENDIF
C-END
RETURN
*
*======================================================================*
* *
* Event End DRAWing: *
* *
*======================================================================*
* *
ENTRY EEDRAW ( ICODE )
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 =
& 'UNFORMATTED' )
END IF
WRITE (IODRAW) 0, ICODE, JTRACK, SNGL (ETRACK), SNGL (WTRACK)
WRITE (IODRAW) SNGL (XSCO), SNGL (YSCO), SNGL (ZSCO), SNGL (RULL)
**
IF( MREG.EQ.3 ) THEN
WRITE(*,*) 'HIT region=', MREG,' edep=',RULL,
&' coo=',XSCO,',',YSCO,',',ZSCO,' code=',ICODE,
&' type=',JTRACK,' ',ETRACK
ETOT_CAL = ETOT_CAL + RULL
ENDIF
IF( MREG.EQ.4 ) THEN
WRITE(*,*) 'scd HIT region=', MREG,' edep=',RULL,
&' coo=',XSCO,',',YSCO,',',ZSCO,' code=',ICODE,
&' type=',JTRACK,' ',ETRACK
ETOT_SCD = ETOT_SCD + RULL
END IF
**
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
WRITE (IODRAW) -NCASE, LSTACK, LSTMAX, SNGL (TKESUM),
& SNGL (WEIPRI)
* +-------------------------------------------------------------------*
* | Patch for heavy ions: it works only for 1 source particle on
* | the stack for the time being
IF ( ABS (ILO(LSTACK)) .GE. 10000 ) THEN
IONID = ILO (LSTACK)
CALL DCDION ( IONID )
WRITE (IODRAW) ( IONID,SNGL(TKE(I)+AMNHEA(-IONID)),SNGL(WT(I)),
& SNGL (XA(I)), SNGL (YA(I)), SNGL (ZA(I)),
& SNGL (TX(I)), SNGL (TY(I)), SNGL (TZ(I)),
& I = 1, LSTACK )
* |
* +-------------------------------------------------------------------*
* |
ELSE
WRITE (IODRAW) ( ILO(I), SNGL (TKE(I)+AM(ILO(I))), SNGL(WT(I)),
& SNGL (XA(I)), SNGL (YA(I)), SNGL (ZA(I)),
& SNGL (TX(I)), SNGL (TY(I)), SNGL (TZ(I)),
& I = 1, LSTACK )
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 *
* 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 ( .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
* No output by default:
RETURN
*=== End of subrutine Mgdraw ==========================================*
END
This archive was generated by hypermail 2.1.6 : Tue May 11 2004 - 00:06:51 CEST