Energy deposition in a calorimeter

From: Paolo MAESTRO (paolo.maestro@pi.infn.it)
Date: Mon May 10 2004 - 17:40:41 CEST

  • Next message: Giuseppe Battistoni: "Re: FLUKA Energy deposition in a calorimeter"

    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


  • Next message: Giuseppe Battistoni: "Re: FLUKA Energy deposition in a calorimeter"

    This archive was generated by hypermail 2.1.6 : Tue May 11 2004 - 00:06:51 CEST