Installation

Get the latest version of Corsika from here (right now this is at least 7.64): ​https://www.ikp.kit.edu/corsika/

Extract Corsika

[0] tar xvf corsika-76400.tar.bz2

Call coconut

[1] cd corsika-76400
[2] ./coconut

You can use the defaults (just press return) for

• Compile in 32 or 64bit mode? [64 bit]
• Which high energy hadronic interaction model do you want to use? [QGSJET 01C]
• Which low energy hadronic interaction model do you want to use? [GHEISHA 2002d]
• Which detector geometry do you have? [horizontal flat]

Now you should get a long list of possible options

Which additional CORSIKA program options do you need?

Enter

(multiple selections accepted, leading '-' removes option): 1a 1b 1c 1e 7c 9 

There options correspond to: CERENKOV, IACT, CEFFIC, TRAJECT, VIEWCONE and ATMEXT.

Now you will get asked more question, answer them as follows

• Cherenkov light vertical (longitudinal) distribution option? [No Cherenkov light distribution at all]
• Do you want Cherenkov light emission angle wavelength dependence? [depending on wavelength] {This corresponds to the CERWLEN option}
• IACTEXT external output file option? [not stored] {no eventio format}

Now, coconut can create the input files, configure the make system and make. To proceed, select *** Finish selection ***.

You should see

 Are you sure you want to continue with these current option selection:
    BERNLOHRDIR IACT TRAJECT CEFFIC CERENKOV CERWLEN IACTDIR VIEWCONE ATMEXT
    yes or no ? (default: yes) >

Which you acknowledge. Now, go ahead!

Options description

CERENKOV[1a]: Enables production of Cherenkov photons

IACT[1b]: Redirects photons to eventio (telescope.dat) file (CER files exists, but only for one telescope and it is empty), extends wavelength range to 2000nm.

CERWLEN: Calculates a wavelength dependent refractive index for each step and replaces the production height (8-th) in the photon output by the wavelength.

IACTEXT: Write also particles to eventio file

CEFFIC[1c]: All three look-up-tables (QE, Atmabs, Mirror) are limited to 105 values between 180nm and 700nm -> This turns the extension to 2000nm in the IACT option off

TRAJECT[1e]: To be able to simulate along a trajetory (zd, az, orientation to magnetic field) with "TRAFLG T" (Turn off with "TRAFLG F")

THIN[2a]: Thinning is mainly interesting for >PeV showers and calculates average results for high densities. With CEFFIC and/or CERWLEN, the wavelength is written 8th value to the Cherenkov bunch into the CER-file instead of the weight. To turn the thinning off for the production use "THIN 0 1 0" (The production should behave like with no THIN option). To turn thinning only off for electromagnetic-particles use "THINEM 0 1"

VIEWCONE[7c]: Logically allows to simulate a cone rather than a rectangle in Alt/Az

ATMEXT[9]: To use an atmosphere from an input file rather than a built-in atmosphere

This is from the tracking of the electrons through the Earth's magnetic field

C  LIMITING FACTOR FOR STEP SIZE OF ELECTRON IN MAGNETIC FIELD
#if __CERENKOV__ && __IACT__
C  LIMIT IN DEFLECTION ANGLE IS 2.5 MILLIRADIAN = 0.143 DEG
C  WE USE A LIMIT OF ABOUT 0.05 DEG (APPROX. 1 MILLIRAD)
      BLIMIT   = 0.001D0 / BNORM
#else
C  WE USE A LIMIT OF ABOUT 11.4 DEG (0.2 RAD)
      BLIMIT   = 0.2D0 / BNORM
#endif

This is from the muon tracking (MUTRAC)

#if __CERENKOV__ && __IACT__
C  SCATTERING ANGLES OF MUONS SHOULD BE SMALLER THAN THE PIXEL SIZE.
*     AUX = MIN( 1.D0, 0.015D0*GAMMA )
C  THE SAME SHOULD HOLD FOR DEFLECTION IN THE GEOMAGNETIC FIELD.
C  HERE USING A MAXIMUM RMS SCATTERING / DEFLECTION ANGLE OF 0.05 DEG
C  AND APPROXIMATE ALL BETA*GAMMA TERMS BY GAMMA.
Cxx      Write(*,*) 'mu step old-style step=',MIN( 1.D0,0.015D0*GAMMA )
C  FOR A MEAN SCATTERING ANGLE THETA WE HAVE A STEP LENGTH OF ABOUT
C  (THETA / (13.6 MEV/(BETA*C*P))**2 RADIATION LENGTHS (PDG),
C  NOT TAKING INTO ACCOUNT THE NON-GAUSSIAN PART OF THE DISTIBUTION.
C  FOR THE MOMENT DON''T CARE ABOUT THE DIFFERENCE BETWEEN THE
C  'COULOMB SCATTERING LENGTH' 37.7 G/CM**2 (=C(21)) AND THE
C  RADIATION LENGTH OF 36.66 OR 36.62 G/CM**2 IN AIR.
C  NOTE: PI/180/(13.6 MEV/(BETA*C*P)) APPROX 0.136*GAMMA FOR MUONS.
      AUX = MIN( 1.D0, ((0.05*0.136)*GAMMA)**2 )
      IF ( BNORMC .GT. 0.D0 ) THEN
C  NOTE: PI/180*PAMA(5)*BETA*GAMMA APPROX 0.00185*GAMMA
         AUX = MIN( AUX, (0.05*0.00185)*GAMMA*RHOF(H)/(BNORMC*C(21)) )
      ENDIF
Cxx      Write(*,*) 'mu step new-style step=',AUX
#else
      AUX = MIN( 10.D0, 0.015D0*GAMMA )
#endif

For me this seems only relevant if someone wants to backtrack a single muon, i.e. determin its arrival direction.

An interesting possibility of the IACT option is:

j) Starting with version 1.25, the package has been prepared for importance sampling ofcore position offsets. This would mean that actual core offsets can be generated in anon-uniform distribution and can extend to different distances, depending on primaryenergy, primary type, zenith angle and so on. This package, however, does not provide areal implementation of importance sampling (other than for testing that the later stagesof the processing properly get the weights for each event). If you do nothing about it,you will get uniformly distributed core offsets as before. If you plan to make use ofimportance sampling, you have to replace the file ’sampling.c’ with an implementationof your choice.

A problem with the IACT option might be this:

m) Starting with version 1.38, the dynamic range in a telescope is basically unlimited dueto automatic thinning of bunches. When a detector sphere is hit by more than the givenmaximum number of bunches, the actual number of bunches is reduced by increasingpowers of two, by discarding every second bunch and increasing the bunch size of theremaining bunches by factors of two. Very large eventio buffers are now possible on64-bit machines. They were formerly limited to less than 1 GiB per telescope array andevent and now to 2 GiB on 32-bit machines but can be increased to 4 Terabytes on 64-bitmachines.

The secret of wavelength 999 and 9999 in eventio:

IACTEXT The interface to theTELOUTfunction is extended by parameters describing theemitting particle. This extended information is stored as an additional photon bunch(after the normal one) with mass, charge, energy, and emission time replacing thecx,cy,photons, andzemfields, respectively, and are identified by a wavelength of 9999.The compact output format is disabled for making that possible. In addition, all particlesarriving at the CORSIKA observation level are included in the eventio format output file,in a photon-bunch like block identified by array and detector numbers 999.

The 700nm problem

This is were Corsika checks the limits

#if __IACT__ && !__CEFFIC__
      IF ( WAVLGL .LT. 100.D0  .OR.  WAVLGU .GT. 2000.D0
     *                         .OR.  WAVLGL .GE. WAVLGU ) THEN
#else
      IF ( WAVLGL .LT. 100.D0  .OR.  WAVLGU .GT. 700.D0
     *                         .OR.  WAVLGL .GE. WAVLGU ) THEN
#endif

This is were the data is read from the files

C=======================================================================

      SUBROUTINE TPDINI( OBLECE )

C-----------------------------------------------------------------------
[...]
C  READ THE TABLE OF ATMOSPHERIC EXTINCTION AND STORE IT IN THE
C  ATMABS MATRIX.
[...]
        DO  I = 1, 105
          READ(MCERABS,21) WLT
 21       FORMAT(I4)
          READ(MCERABS,*,ERR=995,END=995) (ATMABS(I,J), J=0,50)
        ENDDO
[...]
C  CALCULATE THE ATM. EXTINCTION FOR OBS. LEVEL INTERPOLATING THE TABLE
        DO  IWL = 1, 105
          FX0 = ATMABS(IWL,X0)
          FX1 = ATMABS(IWL,X1)
          ATEOBS(IWL) = LINEAR(X,X0,X1,FX0,FX1)
        ENDDO
[...]
C  READ THE QUANTUM EFFECIENCY OF THE PMT
        OPEN(UNIT=MCERQEF,FILE=DATDIR(1:INDEX(DATDIR,' ')-1)//
     *          'quanteff.dat',STATUS='OLD',FORM='FORMATTED',ERR=996)
        READ(MCERQEF,20) TEXTQEF
        READ(MCERQEF,23,ERR=997,END=997) (QUAEFF(I),I = 1,105)
[...]
C  READ THE MIRROR REFLECTIVITY
        OPEN(UNIT=MCERMIR,FILE=DATDIR(1:INDEX(DATDIR,' ')-1)//
     *          'mirreff.dat',STATUS='OLD',FORM='FORMATTED',ERR=998)
        READ(MCERMIR,20) TEXTREF
        READ(MCERMIR,23,END=999,ERR=999) (MIRREF(I), I=1,105)

This is where the index or the y-value and the corresponding x-values are calculated for mirror and pde efficiency.

C=====================================================================

      SUBROUTINE TELEFF( ABSORB )

C-----------------------------------------------------------------------
[...]
C  CALCULATE THE REFERENCE WL AND INDEX OF WL FOR THE INTERPOLATIONS
      RIWL = 1 + INT( (WL-180.D0)/5.D0 )
      WLI0 = RIWL*5 + 175
      WLI1 = RIWL*5 + 175 + 5

This is the atmospheric absorption

*-- Author : V. de Souza Filho, Uni. Campinas   22/06/1999
C=======================================================================

      SUBROUTINE ATABSO( ABSORB )

C-----------------------------------------------------------------------
[...]
C  CALCULATE THE REFERENCE WL AND INDEX OF WL FOR THE INTERPOLATIONS
      RIWL = 1 + INT( (WL-180.D0)/5.D0 )
      WLI0 = RIWL*5 + 175
      WLI1 = RIWL*5 + 175 + 5

Conclusion

If we want to do any absorption in Corsika (CEFFIC), we have to patch Corsika to allow wavelength larger than 700nm or accept that we will miss photons above 700nm.

To be able to define telescope positions ("TELESCOPE"), the IACT option is not required. Plain Corsika can nowadays do that.

Enabling IACT might make sense, due to the code snippets shown above, the automatic thinning (faster) and the work on improving runtime by faster interpolation algorithms.

The disadvantage of the IACT option is that teh data of all telescopes end up in the same file and that *eventio* format has to be used. Although, the iACT options writes a temporary file for each individual telescope, they are sorted together into one file later.

Using IACT means that we have to understand the automatic thinning and implement it and to make sure eventio works.

I have the impression that eventio files are larger than the corsika output (which is a bit strange)... this has to be checked with long runs and high statistics!