CorsikaShowerFile.h

Go to the documentation of this file.

#ifndef _io_CorsikaShowerFile_h_
#define _io_CorsikaShowerFile_h_

#include <io/RawCorsikaFile.h>
#include <io/CorsikaShowerFileParticleIterator.h>
#include <io/CorsikaIOException.h>
#include <io/VEventFile.h>

#include <utl/CoordinateSystemPtr.h>
#include <utl/Branch.h>

#include <string>
#include <fstream>
#include <map>
#include <vector>


namespace evt {
  class Event;
  class ShowerSimData;
  class GaisserHillas6Parameter;
}


namespace io {

  class CorsikaShowerFile : public VEventFile {

  public:
    CorsikaShowerFile() { }
    CorsikaShowerFile(const std::string& fileName, const Mode mode = eRead, utl::Branch* const b = nullptr);
    virtual ~CorsikaShowerFile();

    void Open(const std::string& fileName, const Mode mode = eRead, utl::Branch* const b = nullptr) override;
    void Close() override;
    Status Read(evt::Event& event) override;
    void Write(const evt::Event& event) override;
    Status FindEvent(const unsigned int eventId) override;
    Status GotoPosition(const unsigned int position) override;
    int GetNEvents() override;

    void SetMagneticFieldDeclination(const double d) { fMagneticDeclination = d; }
    void UseRealisticMagneticDeclination(const bool v) { fUseRealisticMagneticDeclination = v; }
    void SetReadLongFile(const bool f) { fReadLongFile = f; }

    CorsikaShowerFileParticleIterator* GetCorsikaShowerFileParticleIterator() const;

  private:
    typedef std::vector<Corsika::RawFile::PositionType> PositionVector;
    typedef std::map<unsigned int, unsigned int> IdToPositionMap;

    void SetRunHeader(const bool first, const Corsika::RunHeader& header);
    void SetEventHeader(const bool first, const unsigned int eventsSoFar, const Corsika::EventHeader& eh);
    void SetHeaderTrailer(evt::Event& event, double& timeShift,
                          evt::GaisserHillas6Parameter& gh, bool& haveGH, double& zenith,
                          const Corsika::EventHeader& header, const Corsika::EventTrailer& trailer);

    Status ScanGroundFile(const std::string& name, Corsika::RawFile& raw, PositionVector& posEVTH, PositionVector& posEVTE,
                          const bool first, const bool noException = false);

    Status ScanLongFile(const std::string& name);
    Status ReadProfile(evt::ShowerSimData& shower, const double zenith);

    Status ScanSteerFile(const std::string& name);

    Status ReadSteerFile(evt::Event& event);

    double GetDensityVsHeight(const double height) const;
    double GetDepthVsHeight(const double h) const;
    double GetSlantDepthCurved(const double hObs, const double theta, const double h) const;

    // data members
    bool fOpenedAnything = false;
    unsigned int fNumberOfEvents = 0;
    io::Corsika::RawFile fRawFile;
    std::map<utl::AttributeMap, io::Corsika::RawFile> fCherenkovFile;

    std::string fRefCoordinateSystemName = "PampaAmarilla";
    //utl::CoordinateSystemPtr fRefCoordinateSystem;
    std::vector<int> fMultiCoreSelection;

    std::string fLongFile;
    bool fReadLongFile = false;

    std::string fSteerFile;
    bool fReadSteerFile = false;

    std::string fListFile;
    bool fReadListFile = false;

    PositionVector fPositionToRaw;
    PositionVector fPositionOfEventTrailer;
    std::map<utl::AttributeMap, PositionVector> fPositionToCherenkovRaw;
    std::map<utl::AttributeMap, PositionVector> fPositionOfCherenkovEventTrailer;
    IdToPositionMap fIdToPosition;
    unsigned int fCurrentPosition = 0;
    int fRunNumber = 0;
    double fCoordinateRotation = 0;
    unsigned int fObservationLevel = 1;
    bool fIsThinned = true;
    bool fIsSlant = true;

    bool fHasCherenkov = false;
    double fCherenkovBandwidthMin = 0;
    double fCherenkovBandwidthMax = 0;

    double fMagneticDeclination = -4.233*utl::degree; // some history halve-wrong value [WEST=-EAST]
    bool fUseRealisticMagneticDeclination = false;

    Corsika::AtmosphereParameters fAtmPars;

    CorsikaShowerFileParticleIterator* fParticleIterator = nullptr;
    std::map<utl::AttributeMap, CorsikaShowerFileParticleIterator*> fCherenkovIterator;

    std::vector<double> fShowerEnergy;
    std::vector<double> fShowerTheta;
    std::vector<double> fShowerPhi;

    friend class REASFile;  // necessary so that REASFile can call ReadProfile()
    friend class REASH5File;  // necessary so that REASFile can call ReadProfile()

  };

}


#endif