CorsikaBlock.h

Go to the documentation of this file.

#ifndef _io_Corsika_Block_h_
#define _io_Corsika_Block_h_

#include <cstddef>


namespace io {

  namespace Corsika {

    typedef float RWORD;  // 4 byte real


    class CherenkovFlags {

      /* C  ---------------------------------------------------------
         C  ELEMENT 77 OF EVENT HEADER HAS THE FOLLOWING CONTENTS IF
         C  CONVERTED TO AN INTEGER WITH SUITABLE ROUNDING APPLIED:
         C  BIT  1: CERENKOV OPTION COMPILED IN
         C       2: IACT OPTION COMPILED IN
         C       3: CEFFIC OPTION COMPILED IN
         C       4: ATMEXT OPTION COMPILED IN
         C       5: ATMEXT OPTION USED WITH REFRACTION ENABLED
         C       6: VOLUMEDET OPTION COMPILED IN
         C       7: CURVED OPTION COMPILED IN (SEE ALSO EVTH(79))
         C       8: USED BY THE IACT INTERFACE
         C       9: SLANT OPTION
         C      10: UNUSED
         C   11-21: TABLE NUMBER FOR EXTERNAL ATMOSPHERE TABLE
         C          (BUT LIMITED TO 1023 IF THE NUMBER IS LARGER)
         C  -------------------------------------------------------- */

      enum Flags {
        eCHERENKOV = 0x00001,
        eIACT      = 0x00002,
        eCEFFIC    = 0x00004,
        eATMEXT    = 0x00008,
        eREFRACT   = 0x00010,
        eVOLUMEDET = 0x00020,
        eCURVED    = 0x00040,
        //           0x00080 used by IACT interface
        eSLANT     = 0x00100,
        // unused    0x00200
        eATMO_ID   = 0xffc00,
        eATMO_IDShift = 10
      };

    public:
      CherenkovFlags(const unsigned int d) : fData(d) { }

      bool IsCHERENKOV() const { return fData & eCHERENKOV; }
      bool IsIACT() const { return fData & eIACT; }
      bool IsCEFFIC() const { return fData & eCEFFIC; }
      bool IsATMEXT() const { return fData & eATMEXT; }
      bool IsREFRACT() const { return fData & eREFRACT; }
      bool IsVOLUMEDET() const { return fData & eVOLUMEDET; }
      bool IsCURVED() const { return fData & eCURVED; }
      bool IsSLANT() const { return fData & eSLANT; }

      int GetATMO_ID() const { return (fData & eATMO_ID) >> eATMO_IDShift; }

    private:
      unsigned int fData = 0;
    };


    static const unsigned int kMaxObservationLevels = 10;
    static const unsigned int kParticlesInBlock  = 39;


    class BlockID {
    public:
      static const size_t kLength = 4;

      void SetID(const char* const theID);

      const char* GetID() const { return fID; }

      bool Is(const char* const theID) const;

    private:
      char fID[kLength];  // hold strings like RUNH etc to
                          // identify the block type
    };


    struct RunHeader {
      BlockID fID;
      RWORD fRunNumber;
      RWORD fDateStart;
      RWORD fVersion;

      RWORD fObservationLevels;
      RWORD fObservationHeight[kMaxObservationLevels]; // in cm

      RWORD fSpectralSlope;
      RWORD fEMin, fEMax;     // in GeV

      RWORD fFlagEGS4;
      RWORD fFlagNKG;

      RWORD fCutoffHadrons;   // in GeV
      RWORD fCutoffMuons;     // in GeV
      RWORD fCutoffElectrons; // in GeV
      RWORD fCutoffPhotons;   // in GeV

      RWORD fConstC[50];
      RWORD fConstUNUSED1a[5]; // inclined observation level parameters. Check Guide.
      RWORD fCoordinateRotation; // c.f. ARRANG option of CORSIKA
      RWORD fConstUNUSED1b[14];
      RWORD fConstCKA[40];
      RWORD fConstCETA[5];
      RWORD fConstCSTRBA[11];
      RWORD fConstUNUSED2[97];
      RWORD fConstXSCATT;
      RWORD fConstYSCATT;
      RWORD fConstHLAY[5];
      RWORD fConstAATM[5];
      RWORD fConstBATM[5];
      RWORD fConstCATM[5];
      RWORD fConstNFLAIN;
      RWORD fConstNFLDIF;
      RWORD fConstNFLPI;
      RWORD fConstNFLCHE;
    };


    struct EventHeader {
      BlockID fID;
      RWORD fEventNumber;
      RWORD fParticleId;
      RWORD fEnergy;
      RWORD fStartingAltitude;        // g cm^-2
      RWORD fFirstTarget;     //
      RWORD fZFirst;          // cm
      RWORD fPx, fPy, fPz;    // GeV
      RWORD fTheta;           // zenith in radians
      RWORD fPhi;             // azimuth in radians

      struct Seed {
        RWORD fSeed;
        RWORD fInitialCallsMod;
        RWORD fInitialCallsDiv;
      };
      RWORD fRandomSequences;
      Seed fSeeds[10];

      RWORD fRunNumber;
      RWORD fDateStart;
      RWORD fVersion;

      RWORD fObservationLevels;
      RWORD fObservationHeight[kMaxObservationLevels]; // in cm

      RWORD fSpectralSlope;   //
      RWORD fEMin, fEMax;     // in GeV

      RWORD fCutoffHadrons;   // in GeV
      RWORD fCutoffMuons;     // in GeV
      RWORD fCutoffElectrons; // in GeV
      RWORD fCutoffPhotons;   // in GeV

      RWORD fNFLAIN;
      RWORD fNFLDIF;
      RWORD fNFLPI0;
      RWORD fNFLPIF;
      RWORD fNFLCHE;
      RWORD fNFRAGM;

      RWORD fBx, fBz;         // magnetic field in mu T
      RWORD fFlagEGS4;
      RWORD fFlagNKG;

      RWORD fFlagLowEnergyModel;
      RWORD fFlagHighEnergyModel;
      RWORD fFlagCerenkov;
      RWORD fFlagNeutrino;
      RWORD fFlagCurved;
      RWORD fFlagComputer;    // 1: IBM, 2: Transputer, 3: DEC/UNIX, 4: Mac,
                              // 5: VAX/VMS, 6: GNU/Linux

      RWORD fThetaMin, fThetaMax; // degrees
      RWORD fPhiMin, fPhiMax; // degrees

      RWORD fCerenkovBunch;
      RWORD fCerenkovNumberX, fCerenkovNumberY;
      RWORD fCerenkovGridX, fCerenkovGridY; // cm
      RWORD fCerenkovDetectorX, fCerenkovDetectorY; // cm
      RWORD fCerenkovOutputFlag;

      RWORD fArrayRotation;
      RWORD fFlagMuonAdditionalInformation;

      RWORD fMultipleScatteringStep;
      RWORD fCerenkovBandwidthMin, fCerenkovBandwidthMax; // nm
      RWORD fUsesOfCerenkovEvent;
      RWORD fCoreX[20], fCoreY[20]; // cm

      RWORD fFlagSIBYLL, fFlagSIBYLLCross;
      RWORD fFlagQGSJET, fFlagQGSJETCross;
      RWORD fFlagDPMJET, fFlagDPMJETCross;
      RWORD fFlagVENUSCross;
      RWORD fFlagMuonMultiple; // 0: Gauss, 1: Moilere
      RWORD fNKGRadialRange;  // cm
      RWORD fEFractionThinningH; // Energy fraction of thinning
                                   // level
      // These are in the CORSIKA manual but not in Lukas's original code
      RWORD fEFractionThinningEM; // Energy fraction of thinning level EM
      RWORD fWMaxHadronic, fWMaxEM;
      RWORD fRMaxThinning;
      RWORD fInnerAngle, fOuterAngle;

      RWORD fHighLowEnergyTransition;
      RWORD fSkimmingIncidence;
      RWORD fSkimmingAxisAltitude;
      RWORD fStartingHeight;
      RWORD fFlagCharm;
      RWORD fFlagEMAdditionalInformation;
      RWORD fUnused[7];
      RWORD fCurvedObsLevel;

      // Padding adjusted according to additions described above
      RWORD fPad[105];
    };


    struct RunTrailer {
      BlockID fID;
      RWORD fRunNumber;
      RWORD fEventsProcessed;

      RWORD fPad[270];
    };


    struct EventTrailer {
      BlockID fID;
      RWORD fEventNumber;

      RWORD fPhotons;
      RWORD fElectrons;
      RWORD fHadrons;
      RWORD fMuons;
      RWORD fParticles;

      // NKG output
      RWORD fLateral1X[21], fLateral1Y[21]; // cm^-2
      RWORD fLateral1XY[21], fLateral1YX[21]; // cm^-2

      RWORD fLateral2X[21], fLateral2Y[21]; // cm^-2
      RWORD fLateral2XY[21], fLateral2YX[21]; // cm^-2

      RWORD fElectronNumber[10]; // in steps of 100 g cm^-2
      RWORD fAge[10];         // in steps of 100 g cm^-2
      RWORD fDistances[10];   // cm
      RWORD fLocalAge1[10];

      RWORD fLevelHeightMass[10]; // g cm^-2
      RWORD fLevelHeightDistance[10]; // cm
      RWORD fDistanceBinsAge[10]; // cm
      RWORD fLocalAge2[10];

      // Longitudinal distribution
      RWORD fLongitudinalPar[6];
      RWORD fChi2;
      // Added according to the CORSIKA manual
      RWORD fWeightedPhotons;
      RWORD fWeightedElectrons;
      RWORD fWeightedHadrons;
      RWORD fWeightedMuons;
      RWORD fNumberEMFromPreshower;

      RWORD fPad[6];
    };


    //


    template<class T, unsigned int n>
    class Padded : public T {
      RWORD fPadding[n];
    };


    struct ParticleDataUnthinned {
      bool IsParticle() const
      { return 0 < fDescription && fDescription < 100000; }
      bool IsNucleus() const
      { return 100000 <= fDescription && fDescription < 9900000; }
      bool IsCherenkov() const
      { return 9900000 <= fDescription; }

      RWORD fDescription;     // Particle ID
      RWORD fPx, fPy, fPz;    // GeV
      RWORD fX, fY;           // cm
      RWORD fTorZ;            // ns or cm
    };


    struct ParticleData : public ParticleDataUnthinned {
      RWORD fWeight;  // absent if CORSIKA compiled w/o thinning
    };


    struct CherenkovDataUnthinned {
      RWORD fPhotonsInBunch;  // includes weight
      RWORD fX, fY;           // cm
      RWORD fU, fV;           // cos to X and Y axis
      RWORD fT;               // ns
      //RWORD fProductionHeight;  // cm
      RWORD fDistance;        // cm
    };


    struct CherenkovData : public CherenkovDataUnthinned {
      RWORD fWavelength;      // absent if CORSIKA compiled w/o thinning
    };


    template<class ParticleType, class CherenkovType,
             unsigned int nParticles, unsigned int padding>
    class BasicBlock {
    public:
      static const unsigned int kPadding = padding;

      struct ParticleBlock {
        ParticleType fParticle[nParticles];
      };

      struct CherenkovBlock {
        CherenkovType fParticle[nParticles];
      };

      typedef Padded<RunHeader, padding> RunHeaderPadded;
      typedef Padded<RunTrailer, padding> RunTrailerPadded;
      typedef Padded<EventHeader, padding> EventHeaderPadded;
      typedef Padded<EventTrailer, padding> EventTrailerPadded;

      union SubBlock {
        RunHeaderPadded fRunHeader;
        RunTrailerPadded fRunTrailer;
        EventHeaderPadded fEventHeader;
        EventTrailerPadded fEventTrailer;
        ParticleBlock fParticleBlock;
        CherenkovBlock fCherenkovBlock;
      };

      bool IsRunHeader() const { return fSubBlock.fRunHeader.fID.Is("RUNH"); }
      bool IsRunTrailer() const { return fSubBlock.fRunTrailer.fID.Is("RUNE"); }
      bool IsEventHeader() const { return fSubBlock.fEventHeader.fID.Is("EVTH"); }
      bool IsEventTrailer() const { return fSubBlock.fEventTrailer.fID.Is("EVTE"); }
      bool IsControl() const
      { return IsRunHeader() || IsRunTrailer() || IsEventHeader() || IsEventTrailer(); }

      const RunHeader& AsRunHeader() const { return fSubBlock.fRunHeader; }
      const RunTrailer& AsRunTrailer() const { return fSubBlock.fRunTrailer; }
      const EventHeader& AsEventHeader() const { return fSubBlock.fEventHeader; }
      const EventTrailer& AsEventTrailer() const { return fSubBlock.fEventTrailer; }
      const ParticleBlock& AsParticleBlock() const { return fSubBlock.fParticleBlock; }
      const CherenkovBlock& AsCherenkovBlock() const { return fSubBlock.fCherenkovBlock; }

    private:
      SubBlock fSubBlock;

    };


    // Padding for thinning (in RWORDS): 39 (thinning) or 0 (no thinning)

    typedef BasicBlock<ParticleData, CherenkovData, kParticlesInBlock, 39> Block;

    typedef BasicBlock<ParticleDataUnthinned, CherenkovDataUnthinned, kParticlesInBlock, 0> BlockUnthinned;

  }

}


#endif