CorsikaShowerFile.cc

Go to the documentation of this file.

#include <io/CorsikaShowerFile.h>
#include <io/CorsikaShowerFileParticleIterator.h>
#include <io/CorsikaShowerFileGeometryProducer.h>
#include <io/CorsikaIOException.h>
#include <io/RawCorsikaFile.h>
#include <io/CorsikaBlock.h>
#include <io/CorsikaUtilities.h>

#include <evt/Event.h>
#include <evt/ShowerSimData.h>
#include <evt/GaisserHillas6Parameter.h>
#include <evt/AtmosphereParameters.h>
#include <evt/DefaultShowerGeometryProducer.h>

#include <utl/TabularStream.h>
#include <utl/ErrorLogger.h>
#include <utl/AugerUnits.h>
#include <utl/TabulatedFunction.h>
#include <utl/MathConstants.h>
#include <utl/Math.h>
#include <utl/PhysicalConstants.h>
#include <utl/PhysicalFunctions.h>
#include <utl/Reader.h>
#include <utl/String.h>
#include <utl/VParticleProperties.h>
#include <utl/ShowerParticleList.h>
#include <utl/ParticlePropertiesFactory.h>
#include <utl/Point.h>
#include <utl/UTMPoint.h>

#include <fwk/CentralConfig.h>
#include <fwk/CoordinateSystemRegistry.h>

#include <boost/format.hpp>
#include <boost/lexical_cast.hpp>
#include <boost/algorithm/string/trim.hpp>
#include <boost/algorithm/string/replace.hpp>
#include <boost/algorithm/string.hpp>
#include <boost/algorithm/string/case_conv.hpp>
#include <boost/filesystem.hpp>
#include <boost/regex.hpp>

#include <TMinuit.h>
#include <TMath.h>

#include <sstream>
#include <string>
#include <cmath>
#include <iostream>

using namespace io;
using namespace evt;
using namespace utl;
using boost::format;
// do not import, namespce pollution! using namespace std;
using std::string;
using std::vector;
using std::ostringstream;
using std::stringstream;
using std::ifstream;
using std::cerr;
using std::endl;


namespace io {

  namespace GHFit {

    // store the necessary variables outside the loop they are defined in
    struct Profile {
      vector<double> fSlantDepth;
      vector<double> fEnergyDeposit;
      double fElectronCutoff;
    };


    // global data for Minuit fitting
    Profile gProfile;


    // The 2 functions below and the structure are used for the correct calculation of GH parameters of up-going events

    inline
    double
    GaisserHillas4ParameterFunction(const double x, const double dEdXmax, const double x0, const double xmax, const double lambda)
    {
      const double xm0 = xmax - x0;
      return dEdXmax * pow((x - x0) / xm0, xm0 / lambda) * exp((xmax - x)/lambda);
    }


    // function for TMinuit for dE/dX GH parameters -> get dEdXmax, xMax, x0 and lambda
    // fit done for dE/dx in PeV/(g/cm2)
    void
    GH4ParamFitFunctiondEdxGH(int& /*npar*/, double* /*gin*/, double& f, double* par, int /*iflag*/)
    {
      double chi2 = 0;

      const double dEdXmax = par[0];
      const double x0 = par[1];
      const double xMax = par[2];
      const double lambda = par[3];

      for (unsigned int i = 0, n = GHFit::gProfile.fSlantDepth.size(); i < n; ++i) {
        const double x = GHFit::gProfile.fSlantDepth[i] / (g/cm2);
        const double dEdX = GHFit::gProfile.fEnergyDeposit[i] / (PeV / (g/cm2));
        //const double dEdXerr = 0.01 * data.EnergyDepositStruct[i] + 0.0003;  // not used as we don't care about chi2 value
        const double dEdXFct = GaisserHillas4ParameterFunction(x, dEdXmax, x0, xMax, lambda);
        const double residual = dEdX - dEdXFct;
        chi2 += Sqr(residual);
      }

      f = chi2;
    }

  }


  inline
  bool
  Contains(const string& str, const string& what)
  {
    return str.find(what) != string::npos;
  }


  vector<string>
  FixedColumnWidthSplit(const string& line,
                        const unsigned int nColumns,
                        const unsigned int firstColumnWidth, const unsigned int columnWidth)
  {
    using boost::algorithm::trim;

    if (!nColumns || line.length() != firstColumnWidth + (nColumns - 1)*columnWidth) {
      const string err = "Line is too short to fix it.";
      ERROR(err);
      throw CorsikaIOException(err);
    }
    vector<string> res(nColumns);
    res[0] = line.substr(0, firstColumnWidth);
    trim(res[0]);
    unsigned int pos = firstColumnWidth;
    for (unsigned int i = 1; i < nColumns; ++i) {
      res[i] = line.substr(pos, columnWidth);
      trim(res[i]);
      pos += columnWidth;
    }
    return res;
  }


  CorsikaShowerFile::~CorsikaShowerFile()
  {
    delete fParticleIterator;
  }


  CorsikaShowerFile::CorsikaShowerFile(const string& fileName,
                                       const Mode mode,
                                       utl::Branch* const branch)
  {
    if (!branch || !*branch) {
      WARNING("No Branch was passed for configuration.");
    } else {
      fObservationLevel = 1;
      if (branch->GetChild("CorsikaObservationLevel"))
        branch->GetChild("CorsikaObservationLevel").GetData(fObservationLevel);
      if (fObservationLevel < 1 || fObservationLevel > 10) {
        ostringstream err;
        err << "CORSIKA observation level must be within [1,10] and " << fObservationLevel << " is not allowed";
        ERROR(err);
        throw CorsikaIOException(err.str());
      }
      ostringstream info;
      info << "Reading particles at observation level: " << fObservationLevel;
      INFO(info);

      info.str("");
      if (branch->GetChild("UseRealisticMagneticDeclination"))
        branch->GetChild("UseRealisticMagneticDeclination").GetData(fUseRealisticMagneticDeclination);
      info << "Using magnetic field declination from time-dependent model "
              "(If False, a hardcoded value (-4.233 deg) is taken or the value"
              " specified in the .reas file (if used)): "
           << fUseRealisticMagneticDeclination;
      INFO(info);

      Branch bRef = branch->GetChild("ReferenceCoordinateSystem");
      if (bRef) {
        fRefCoordinateSystemName = bRef.Get<string>();
        ostringstream inf;
        inf << "Core reference system: " << fRefCoordinateSystemName;
        INFO(inf);
      }
      Branch bMulti = branch->GetChild("MultiCoreSelection");
      if (bMulti) {
        ostringstream info;
        info << "Multi-core selection: ";
        string entry = bMulti.Get<string>();
        boost::trim(entry);
        if (boost::to_lower_copy(entry) == "all") {
          info << "all";
        } else {
          vector<string> strs;
          vector<string> range;
          boost::split(strs, entry, boost::is_any_of(","));
          string comma = "";
          for (const auto& it : strs) {
            boost::split(range, it, boost::is_any_of("-"));
            auto x = range.begin();
            const int low = boost::lexical_cast<int>(range[0]);
            int high = low;
            ++x;
            if (x != range.end())
              high = boost::lexical_cast<int>(range[1]);
            for (int i = low; i <= high; ++i) {
              if (i < 1) {
                const string err = "CORIKSA (multi) core counting start at 1.";
                ERROR(err);
                throw CorsikaIOException(err);
              }
              fMultiCoreSelection.push_back(i);
              info << comma << i;
              comma = ", ";
            }
          }
        }
        INFO(info);
      }
    }
    Open(fileName, mode, branch);
  }


  void
  CorsikaShowerFile::Open(const string& fileName,
                          const Mode mode,
                          utl::Branch* const branch)
  {
    if (fRawFile.IsOpen())
      Close();
    for (auto& cherFile : fCherenkovFile)
      if (cherFile.second.IsOpen())
        cherFile.second.Close();

    if (mode != eRead) {
      const string msg = "Cannot open file '" + fileName + "' - Corsika files are read-only";
      ERROR(msg);
      throw CorsikaIOException(msg);
    }

    fOpenedAnything = false;
    bool foundShowerDefinition = false;

    if (!branch) {

      ostringstream warn;
      warn << "No Branch was passed for configuration. Try to open file(s) based on " << fileName << ".";
      WARNING(warn);

      // base name should be as close as possible to "DATxxxxxx" as possible
      string base = fileName.substr(0, fileName.rfind(".")); // check for any extension and remove
      boost::replace_last(base, "RUN", "DAT");
      boost::replace_last(base, "SIM", "DAT");
      // check if fileName is a ground particle file, or w/ extension
      if (ScanGroundFile(base, fRawFile, fPositionToRaw, fPositionOfEventTrailer, !foundShowerDefinition, true) == eSuccess) {
        foundShowerDefinition = true;
        fOpenedAnything = true;
        INFO(string("Found: ") + base);
      } else if (ScanGroundFile(base + ".part", fRawFile, fPositionToRaw, fPositionOfEventTrailer, !foundShowerDefinition, true) == eSuccess) {
        foundShowerDefinition = true;
        fOpenedAnything = true;
        INFO(string("Found: ") + base + ".part");
      } else if (ScanGroundFile(base + ".grdptcls", fRawFile, fPositionToRaw, fPositionOfEventTrailer, !foundShowerDefinition, true) == eSuccess) {
        foundShowerDefinition = true;
        fOpenedAnything = true;
        INFO(string("Found: ") + base + ".grdptcls");
      } else if (ScanGroundFile(fileName, fRawFile, fPositionToRaw, fPositionOfEventTrailer, !foundShowerDefinition, true) == eSuccess) {
        foundShowerDefinition = true;
        fOpenedAnything = true;
        INFO(string("Found: ") + fileName);
      } /*else if (boost::filesystem::exists(fileName)) {
        fRawFile.Open(fileName);
        ScanGroundFile(fRawFile, fPositionToRaw, fPositionOfEventTrailer, !foundShowerDefinition);
        foundShowerDefinition = true;
        fOpenedAnything = true;
      }*/

      // check if there is a long file
      if (ScanLongFile(base + ".long") == eSuccess) {
        fLongFile = base + ".long";
        fReadLongFile = true;
        fOpenedAnything = true;
        INFO(string("Found profile: ") + base + ".long");
      }

      // check if ther is an INP steering file
      if (ScanSteerFile(base + ".inp") == eSuccess) {
        if (!fNumberOfEvents)
          fNumberOfEvents = 1;
        fSteerFile = base + ".inp";
        fReadSteerFile = true;
        fOpenedAnything = true;
        foundShowerDefinition = true; // at least maybe...
        INFO(string("Found steerfile: ") + base + ".inp");
      } else {
        string testBase = base;
        boost::replace_last(testBase, "DAT", "SIM");
        if (ScanSteerFile(testBase + ".inp") == eSuccess) {
          if (!fNumberOfEvents)
            fNumberOfEvents = 1;
          fSteerFile = testBase + ".inp";
          fReadSteerFile = true;
          fOpenedAnything = true;
          foundShowerDefinition = true; // at least maybe...
          INFO(string("Found steerfile: ") + testBase + ".inp");
        } else {
          string testBase = base;
          boost::replace_last(testBase, "DAT", "RUN");
          if (ScanSteerFile(testBase + ".inp") == eSuccess) {
            if (!fNumberOfEvents)
              fNumberOfEvents = 1;
            fSteerFile = testBase + ".inp";
            fReadSteerFile = true;
            fOpenedAnything = true;
            foundShowerDefinition = true; // at least maybe...
            INFO(string("Found steerfile: ") + testBase + ".inp");
          }
        }
      }
      // check if there is an LST output file
      if (boost::filesystem::exists(base + ".lst")) {
        // not implemented yet
        /*fListFile = base + ".lst";
        fReadListFile = true;
        fOpenedAnything = true;
        foundShowerDefinition = true;
        */
      }

    } else {

      utl::Branch mapping = branch->GetChild("InputMapping");
      utl::Branch bRegex = branch->GetChild("PatternRegex");
      utl::Branch bPath = branch->GetChild("PatternPath");

      boost::regex re;
      if (bRegex) {

        try {
          re = boost::regex(bRegex.Get<string>());
        } catch (const boost::regex_error& e) {
          const string err = "Please make sure to specify valid regex expression in PatternRegex!";
          ERROR(err);
          throw CorsikaIOException(err);
        }

      } else {

        if (!bPath) {
          const string err = "InputMapping requires one of pattern_regex or pattern_path attribute!";
          ERROR(err);
          throw CorsikaIOException(err);
        }

        auto pattern = bPath.Get<string>();
        boost::replace_all(pattern, "\\", "\\\\");
        boost::replace_all(pattern, "^", "\\^");
        boost::replace_all(pattern, ".", "\\.");
        boost::replace_all(pattern, "$", "\\$");
        boost::replace_all(pattern, "|", "\\|");
        boost::replace_all(pattern, "(", "\\(");
        boost::replace_all(pattern, ")", "\\)");
        boost::replace_all(pattern, "[", "\\[");
        boost::replace_all(pattern, "]", "\\]");
        boost::replace_all(pattern, "*", "\\*");
        boost::replace_all(pattern, "+", "\\+");
        boost::replace_all(pattern, "?", "\\?");
        boost::replace_all(pattern, "/", "\\/");
        // I hope this is enough .... please extend if needed
        boost::replace_all(pattern, "\\?", "."); // -> group !!
        boost::replace_all(pattern, "\\*", "(.*)"); // -> group
        re = boost::regex(pattern);
      }

      vector<string> patterns;
      try {
        boost::smatch what;
        if (boost::regex_match(fileName, what, re)) {
          for (int i = 1, n = what.size(); i < n; ++i)
            patterns.push_back(what[i]);
        } else {
          const string err = "Please make sure you specify valid regex/path expression! No match was found.";
          ERROR(err);
          throw CorsikaIOException(err);
        }
      } catch (const boost::regex_error& e) {
        const string err = "Please make sure you specify valid regex/path expression!";
        ERROR(err);
        throw CorsikaIOException(err);
      }

      for (mapping = mapping.GetFirstChild(); mapping; mapping = mapping.GetNextSibling()) {

        // determine pattern(s) data
        const string name = mapping.GetName();
        string path = mapping.Get<string>();

        for (unsigned int i = 0; i < patterns.size(); ++i) {
          ostringstream marker;
          marker << '(' << i+1 << ')';
          boost::replace_all(path, marker.str(), patterns[i]);
        }

        {
          ostringstream inf;
          inf << "Read " << name << " from file " << path;
          INFO(inf);
        }

        if (name == "ProfileData") {

          fReadLongFile = true;
          fLongFile = path;
          fOpenedAnything = true;

        } else if (name == "GroundData") {

          ScanGroundFile(path, fRawFile, fPositionToRaw, fPositionOfEventTrailer, !foundShowerDefinition);
          foundShowerDefinition = true;
          fOpenedAnything = true;

        } else if (name == "ShowerINFO") { // .lst files

          foundShowerDefinition = true;
          fOpenedAnything = true;

          ifstream info(path);
          string line;
          while (getline(info, line)) {
            if (line.empty())
              continue;
            if (line.find("energy_prim") == 0)
              fShowerEnergy.push_back(stod(line.substr(13, line.length() - 13)));
            if (line.find("theta_prim") == 0)
              fShowerTheta.push_back(stod(line.substr(12, line.length() - 12)));
            if (line.find("phi_prim") == 0)
              fShowerPhi.push_back(stod(line.substr(11, line.length() - 11)));
          }

          const string err = "this is not yet completely implemented";
          ERROR(err);
          throw CorsikaIOException(err);

        } else if (name == "ShowerDBASE") {

          foundShowerDefinition = true;
          fOpenedAnything = true;

          ifstream dbase(path);
          string line;
          while (getline(dbase, line)) {
            if (line.empty())
              continue;
            if (line.find("#energy_prim") == 0) {
              vector<string> strs;
              boost::split(strs, line, boost::is_any_of("#"));
              fShowerEnergy.push_back(stod(strs[1]));
              fShowerTheta.push_back(stod(strs[3]));
              fShowerPhi.push_back(stod(strs[5]));
            }
          }
          dbase.close();

          const string err = "this is not yet completely implemented";
          ERROR(err);
          throw CorsikaIOException(err);

        } else if (name == "TelescopeData") {

          fOpenedAnything = true;

          const AttributeMap am = mapping.GetAttributes();
          if (fCherenkovFile.count(am)) {
            ostringstream err;
            err << "Already defined TelescopeData entry for attributes: ";
            if (am.empty()) {
              err << "[none]";
            } else {
              string comma = "";
              for (const auto& iatt : am) {
                err << comma << '"' << iatt.first << "\"=\"" << iatt.second << '"';
                comma = ", ";
              }
            }
            ERROR(err);
            throw CorsikaIOException(err.str());
          }
          ScanGroundFile(path, fCherenkovFile[am], fPositionToCherenkovRaw[am], fPositionOfCherenkovEventTrailer[am], !foundShowerDefinition);
          foundShowerDefinition = true;

        } else {
          ERROR("unknown input mapping");
          throw CorsikaIOException("unknown input mapping" );
        }
      }
    }

    if (!fOpenedAnything) {
      const string err = "Did not open anything! Stop here.";
      WARNING(err);
      throw CorsikaIOException(err);
    }

    if (!foundShowerDefinition) {
      const string err = "No suited input file with primary shower data found. Cannot produce full MC evnets.";
      WARNING(err);
    }
    GotoPosition(0);
  }


  // Helper function to set runheader data for both thinned/unthinned
  void
  CorsikaShowerFile::SetRunHeader(const bool first, const io::Corsika::RunHeader& header)
  {
    if (first) {
      fRunNumber = int(header.fRunNumber);
      // Prior to this version, the h parameters were not saved.
      if (header.fVersion >= 6.9) {
        for (int i = 0; i < 5; ++i)
          fAtmPars.fHLAY[i] = header.fConstHLAY[i] * cm;
        if (fAtmPars.fHLAY[0] < 0)
          fAtmPars.fHLAY[0] = 0;
      } else {
        ostringstream msg;
        msg << "The current CORSIKA file was created using version " << header.fVersion << ". "
               "The atmospheric layers are not defined! Using the values from US std. atmosphere.";
        WARNING(msg);
        fAtmPars.fHLAY[0] = 0*km;
        fAtmPars.fHLAY[1] = 4*km;
        fAtmPars.fHLAY[2] = 10*km;
        fAtmPars.fHLAY[3] = 40*km;
        fAtmPars.fHLAY[4] = 100*km;
      }
      for (int i = 0; i < 5; ++i) {
        fAtmPars.fAATM[i] = header.fConstAATM[i] * g/cm2;
        fAtmPars.fBATM[i] = header.fConstBATM[i] * g/cm2;
        fAtmPars.fCATM[i] = header.fConstCATM[i] * cm;
        fAtmPars.fZLAY[i] = GetDepthVsHeight(fAtmPars.fHLAY[i]);
      }
      INFO("Loaded atmosphere parameters:");
      TabularStream tab(" l | l | l | l ");
      tab << "h[km]" << endc << "a[g/cm2]" << endc << "b[g/cm2]" << endc << "c[km]" << endr
          << hline;
      for (int i = 0; i < 5; ++i) {
        if (i == 4)
          tab << hline;
        tab << fAtmPars.fHLAY[i] / km << endc
            << fAtmPars.fAATM[i] / (g/cm2) << endc
            << fAtmPars.fBATM[i] / (g/cm2) << endc
            << fAtmPars.fCATM[i] / km << endr;
      }
      tab << delr;
      cerr << tab << endl;
    } else {
      // if NOT first -> consistency check
      if (fRunNumber != int(header.fRunNumber)) {
        const string msg = "Inconsistent RunNumber found";
        ERROR(msg);
        throw CorsikaIOException(msg);
      }
      // Prior to this version, the h parameters were not saved.
      if (header.fVersion >= 6.9) {
        for (int i = 0; i < 5; ++i) {
          if (fAtmPars.fHLAY[i] != header.fConstHLAY[i] * cm) {
            const string msg = "Inconsistent Atmosphere HLAY found";
            ERROR(msg);
            throw CorsikaIOException(msg);
          }
        }
      }
      for (int i = 0; i < 5; ++i) {
        if (fAtmPars.fAATM[i] != header.fConstAATM[i] * g/cm2) {
          const string msg = "Inconsistent Atmosphere AATM found";
          ERROR(msg);
          throw CorsikaIOException(msg);
        }
        if (fAtmPars.fBATM[i] != header.fConstBATM[i] * g/cm2) {
          const string msg = "Inconsistent Atmosphere BATM found";
          ERROR(msg);
          throw CorsikaIOException(msg);
        }
        if (fAtmPars.fCATM[i] != header.fConstCATM[i] * cm) {
          const string msg = "Inconsistent Atmosphere CATM found";
          ERROR(msg);
          throw CorsikaIOException(msg);
        }
        /*if (fAtmPars.fZLAY[i] != GetDepthVsHeight(fAtmPars.fHLAY[i])) {
          const string msg = "Inconsistent Atmosphere found";
          ERROR(msg);
          throw CorsikaIOException(msg);
        }*/
      }
    }
  }


  // Helper function to set runheader data for both thinned/unthinned
  void
  CorsikaShowerFile::SetEventHeader(const bool first, const unsigned int eventsSoFar, const io::Corsika::EventHeader& eh)
  {
    if (first) { // first read and record
      fCoordinateRotation = eh.fArrayRotation * radian;
      fIdToPosition[int(eh.fEventNumber)] = eventsSoFar;
      const Corsika::CherenkovFlags chf(eh.fFlagCerenkov);
      fIsSlant = chf.IsSLANT();
      if (chf.IsCHERENKOV()) {
        if (chf.IsCEFFIC()) { // || chf.IsCURVED() ||
                              //blockUnth.AsEventHeader().fCerenkovOutputFlag only support option 1 !!
          ostringstream msg;
          msg << "NOTE: Offline does currently not support to read Cherenkov photons from ground particle files with optinos: "
                 "CEFFIC = " << chf.IsCEFFIC() << ", "
                 "CURVED = " << chf.IsCURVED() << ", "
                 "OUTPUT = " << eh.fCerenkovOutputFlag;
          ERROR(msg);
          throw CorsikaIOException(msg.str());
        } else {
          fCherenkovBandwidthMin = eh.fCerenkovBandwidthMin * nanometer;
          fCherenkovBandwidthMax = eh.fCerenkovBandwidthMax * nanometer;
          ostringstream inf;
          inf << "Cherenkov option is available: "
                 "CHERENKOV = " << chf.IsCHERENKOV() << ", "
                 "CEFFIC = " << chf.IsCEFFIC() << ", "
                 "CURVED = " << chf.IsCURVED() << ", "
                 "Ch-Output = " << eh.fCerenkovOutputFlag << ", "
                 "lambda = [" << fCherenkovBandwidthMin / nanometer << ", "
                            << fCherenkovBandwidthMax / nanometer << "]";
          INFO(inf);
        }
      }
    } else {
      // if NOT first -> check consistency
      if (fCoordinateRotation != eh.fArrayRotation*radian) {
        const string msg = "Inconsistent AARANG found";
        ERROR(msg);
        throw CorsikaIOException(msg);
      }
      if (fIdToPosition[int(eh.fEventNumber)] != eventsSoFar) {
        const string msg = "Inconsistent EventID found";
        ERROR(msg);
        throw CorsikaIOException(msg);
      }
      const Corsika::CherenkovFlags chf(eh.fFlagCerenkov);
      if (fIsSlant != chf.IsSLANT()) {
        const string msg = "Inconsistent SLANT";
        ERROR(msg);
        throw CorsikaIOException(msg);
      }
      if (chf.IsCHERENKOV()) {
        if (chf.IsCEFFIC()) { // || chf.IsCURVED() ||
          ;
        } else {
          if (fCherenkovBandwidthMin != eh.fCerenkovBandwidthMin * nanometer) {
            const string msg = "Inconsistent CherWLmin found";
            ERROR(msg);
            throw CorsikaIOException(msg);
          }
          if (fCherenkovBandwidthMax != eh.fCerenkovBandwidthMax * nanometer) {
            const string msg = "Inconsistent CherWLmax found";
            ERROR(msg);
            throw CorsikaIOException(msg);
          }
        }
      }
    }
  };


  Status
  CorsikaShowerFile::ScanGroundFile(const std::string& name,
                                    io::Corsika::RawFile& raw,
                                    PositionVector& positionToRaw,
                                    PositionVector& positionOfEventTrailer,
                                    const bool first,
                                    const bool noException)
  {
    if (!boost::filesystem::exists(name))
      return eFail;

    if (!raw.Open(name, noException))
      return eFail;

    raw.SeekTo(0, true); // reset in order to start fresh for sure
    unsigned int eventsSoFar = 0;

    unsigned int blockIndex = 0;
    bool foundEventHeader = false;
    bool foundRunHeader = false;
    fIsSlant = true;
    fIsThinned = true;

    Corsika::BlockUnthinned blockUnth;
    while (raw.GetNextBlockUnthinned(blockUnth) && !blockUnth.IsRunTrailer()) {
      ++blockIndex;
      if (blockIndex == 1 && !blockUnth.IsRunHeader()) // check if this is not un-thinned
        break;
      if (blockIndex == 2 && !blockUnth.IsEventHeader()) // check if this is not un-thinned
        break;
      const Corsika::RawFile::PositionType rawPosition = raw.GetNextPosition();
      if (blockUnth.IsEventHeader()) {
        fIsThinned = false;
        foundEventHeader = true;
        positionToRaw.push_back(rawPosition - 1);
        const Corsika::EventHeader& eh = blockUnth.AsEventHeader();
        SetEventHeader(first, eventsSoFar, eh);
        ++eventsSoFar;
      } else {
        if (blockUnth.IsEventTrailer())
          positionOfEventTrailer.push_back(rawPosition - 1);
        else if (blockUnth.IsRunHeader()) {
          const Corsika::RunHeader& RUNH = blockUnth.AsRunHeader();
          SetRunHeader(first, RUNH);
          foundRunHeader = true;
        } // end runheader
      }
      if (blockIndex > 400 && !foundRunHeader) {
        const string msg = "Error scanning Corsika ground file: could not find run header";
        ERROR(msg);
        throw CorsikaIOException(msg);
      }
      if (blockIndex > 400 && !foundEventHeader)
        break;
    }

    if (!foundEventHeader) { // now try thinned reading

      eventsSoFar = 0;
      blockIndex = 0;
      raw.SeekTo(0, true); // also "reset" in order to cleanly go from unthinned to thinned mode

      Corsika::Block blockTh;
      while (raw.GetNextBlock(blockTh) && !blockTh.IsRunTrailer()) {
        ++blockIndex;
        if (blockIndex == 1 && !blockTh.IsRunHeader()) // this is not un-thinned
          break;
        if (blockIndex == 2 && !blockTh.IsEventHeader()) // this is not un-thinned
          break;
        const Corsika::RawFile::PositionType rawPosition = raw.GetNextPosition();
        if (blockTh.IsEventHeader()) {
          foundEventHeader = true;
          fIsThinned = true;
          positionToRaw.push_back(rawPosition - 1);
          const Corsika::EventHeader& eh = blockTh.AsEventHeader();
          SetEventHeader(first, eventsSoFar, eh);
          ++eventsSoFar;
        } else if (blockTh.IsEventTrailer())
          positionOfEventTrailer.push_back(rawPosition - 1);
        else if (blockTh.IsRunHeader()) {
          const Corsika::RunHeader& rh = blockTh.AsRunHeader();
          if (!foundRunHeader)
            SetRunHeader(first, rh);
          foundRunHeader = true;
        }
        if (blockIndex > 400) {
          if (!foundRunHeader) {
            const string err = "Error scanning thinned Corsika ground file: "
                               "could not find run header";
            ERROR(err);
            throw CorsikaIOException(err);
          }
          if (!foundEventHeader) {
            const string err = "Error scanning Corsika ground file: "
                               "could not find Event header";
            ERROR(err);
            throw CorsikaIOException(err);
          }
        }
      }

      if (!blockTh.IsRunTrailer()) {
        const string err = "Error scanning Corsika ground file: could not find run end";
        ERROR(err);
        throw CorsikaIOException(err);
      }
    }

    if (positionToRaw.size() != positionOfEventTrailer.size()) {
      const string err = "Found different number of event-headers and -trailers";
      ERROR(err);
      throw CorsikaIOException(err);
    }

    if (!first) {
      // consistency checking mode: .cher files vs. DAT file
      if (eventsSoFar != fPositionToRaw.size()) {
        ostringstream err;
        err << "Found different number of events in Cherenkov file (" << eventsSoFar << ") "
          "compared to ground particle file (" << fPositionToRaw.size() << ")";
        ERROR(err);
        throw CorsikaIOException(err.str());
      }
    }

    fNumberOfEvents = positionToRaw.size();

    return eSuccess;
  }


  void
  CorsikaShowerFile::Close()
  {
    if (fRawFile.IsOpen())
      fRawFile.Close();

    fPositionToRaw.clear();
    fIdToPosition.clear();

    delete fParticleIterator;
    fParticleIterator = nullptr;
  }


  void
  CorsikaShowerFile::SetHeaderTrailer(evt::Event& event, double& timeShift,
                                      GaisserHillas6Parameter& gh, bool& haveGH, double& zenith,
                                      const Corsika::EventHeader& header, const Corsika::EventTrailer& trailer)
  {
    // if fUseRealisticMagneticDeclination is true the realistic declination is optained by a model
    // inside CorsikaShowerFileGeometryProducer
    const double magneticNorthDeclination = fUseRealisticMagneticDeclination ? 0 : fMagneticDeclination;
    const auto refCoordinateSystem = fwk::CoordinateSystemRegistry::Get(fRefCoordinateSystemName);
    const bool multiCore = header.fUsesOfCerenkovEvent > 0;
    event.MakeSimShower(
      io::CorsikaShowerFileGeometryProducer(
        header.fTheta*radian, header.fPhi*radian + fCoordinateRotation,
        header.fObservationHeight[fObservationLevel-1]*cm,
        multiCore,
        refCoordinateSystem,
        fCoordinateRotation + magneticNorthDeclination,
        fUseRealisticMagneticDeclination
      )
    );
    auto& shower = event.GetSimShower();

    zenith = header.fTheta * radian;
    const double cosZenith = cos(zenith);
    const double xmax = trailer.fLongitudinalPar[2] * g/cm2;
    const double x0 = trailer.fLongitudinalPar[1] * g/cm2;
    const double nmax = trailer.fLongitudinalPar[0];
    const double rMin = header.fRMaxThinning * cm;  // max radius for radial thinning

    const bool hasValidGHfit = (trailer.fChi2 > 0 && trailer.fChi2 < 1e13) && (nmax || xmax || x0);

    shower.SetPrimaryParticle(Corsika::CorsikaToPDG(int(header.fParticleId)));
    shower.SetEnergy(header.fEnergy * GeV);
    shower.SetMuonNumber(trailer.fMuons);
    shower.SetGroundParticleCoordinateSystemZenith(zenith);
    // phi always relative to magnetic north in CORSIKA. This is the definition we also use (JUST) here!
    shower.SetGroundParticleCoordinateSystemAzimuth(
      header.fPhi*radian
      + 180*degree  // incoming vs. outgoing
      - fCoordinateRotation  // eventual array rotation
    );
    shower.SetMinRadiusCut(rMin);

    shower.SetShowerNumber(int(header.fEventNumber));
    shower.SetShowerRunId(boost::lexical_cast<string>(header.fRunNumber));

    shower.SetEnergyCutoff(header.fCutoffElectrons * GeV, ShowerSimData::eElectron);
    shower.SetEnergyCutoff(header.fCutoffMuons * GeV, ShowerSimData::eMuon);

    const double bx = header.fBx * micro*tesla;
    const double bz = header.fBz * micro*tesla;
    shower.SetMagneticFieldInclination(atan2(bz, bx));
    shower.SetMagneticFieldStrength(sqrt(bx * bx + bz * bz));

    shower.SetMinCherenkovWavelength(header.fCerenkovBandwidthMin*nanometer);
    shower.SetMaxCherenkovWavelength(header.fCerenkovBandwidthMax*nanometer);

    // Corsika starts at the top of the atmosphere, not
    const double heightObsLevel = header.fObservationHeight[fObservationLevel - 1] * cm;
    const double heightFirstInt = fabs(header.fZFirst) * cm;

    const double hAtmBoundary = 112.8292*km;
    // for the SLANT and CURVED options, clock starts at the margin of
    // the atmosphere. This is indicated by fZFirst < 0
    double hReference = (header.fZFirst < 0) ? hAtmBoundary : heightFirstInt;

    // calculating slant depth of first interaction point
    if (cosZenith > 0) {
      shower.SetXFirst(header.fFlagCurved ? GetSlantDepthCurved(heightObsLevel, zenith, heightFirstInt) :
                                            GetDepthVsHeight(heightFirstInt) / fabs(cosZenith));
    } else if (cosZenith < 0) {
      // for up-going showers atmosphere starts at the bottom and ends at the top.amount of
      // atmosphere at 1.4 km = 875.539 g/cm2 -> set to 875.6
      shower.SetXFirst(header.fFlagCurved ? GetSlantDepthCurved(heightObsLevel, zenith, heightFirstInt) :
                                     (875.6*(g/cm2) - GetDepthVsHeight(heightFirstInt)) / fabs(cosZenith));
    }
    {
      ostringstream info;
      const utl::ParticlePropertiesPtr& pid = utl::ParticlePropertiesFactory::Create(Corsika::CorsikaToPDG(int(header.fParticleId)));
      info << "Primary = " << pid.get()->GetName() << ", "
              "Energy = " << header.fEnergy* GeV / eV << " eV, "
              "Height of first interaction = " << heightFirstInt / km << " km, "
              "X_1 = " << shower.GetXFirst() / (g/cm2) << " g/cm2, "
              "Zenith = " << zenith / deg << " deg, "
              "Azimuth = " << (header.fPhi + 180*deg)/ deg << " deg, "
              "cutoff_e = " << header.fCutoffElectrons << " GeV, "
              "cutoff_mu = " << header.fCutoffMuons << " GeV, "
              "r_Min^thin = " << rMin / cm << " cm";
      INFO(info);
      GHFit::gProfile.fElectronCutoff = header.fCutoffElectrons*GeV;
    }

    if (fCoordinateRotation) {
      ostringstream info;
      info << "CORSIKA output is rotated by " << fCoordinateRotation / deg << " deg with respect to CORSIKA USER GUIDE";
      INFO(info);
    }

    if (header.fCurvedObsLevel) {
      ostringstream info;
      info << "CORSIKA output is on curved surface" << endl;
    }

    for (int iCore = 0; iCore < header.fUsesOfCerenkovEvent; ++iCore) {
      ostringstream info;
      info << "Found simulated multi-core #" << iCore+1 << " "
              "x=" << header.fCoreX[iCore] / 100. << " m, "
              "y=" << header.fCoreY[iCore] / 100. << " m, "
              "height=" << heightObsLevel / m << " m ";
      if (fMultiCoreSelection.empty() || // in this case all cores are selected
          find(fMultiCoreSelection.begin(), fMultiCoreSelection.end(), iCore+1) != fMultiCoreSelection.end()) {
        info << " [selected]";
        const UTMPoint utm(Point(0,0,0, refCoordinateSystem), ReferenceEllipsoid::eWGS84);
        const double origin_height = utm.GetHeight();
        const Point mcCore(header.fCoreX[iCore]*cm, header.fCoreY[iCore]*cm, heightObsLevel-origin_height, refCoordinateSystem);
        shower.AddSimCore(mcCore);
      } else {
        info << " [ignore]";
      }
      INFO(info);
    }

    ostringstream info;
    if (header.fFlagCurved) {

      info << "CURVED version: ";

      // value taken from CORSIKA
      const double kREarth = 6.371315e6*m; // radius of the earth

      timeShift = (Sqr((kREarth + heightObsLevel) * cosZenith) +
                   Sqr(hReference - heightObsLevel) +
                   2 * (kREarth + heightObsLevel) * (hReference - heightObsLevel));
      timeShift = sqrt(timeShift);
      timeShift -= (kREarth + heightObsLevel) * cosZenith;
      timeShift /= kSpeedOfLight;

    } else {
      timeShift = (hReference - heightObsLevel) / (cosZenith * kSpeedOfLight);
    }
    info << "TimeShift to core: " << timeShift / ns << " ns";
    INFO(info);

    if (fObservationLevel > header.fObservationLevels) {
      ostringstream info;
      info << "The requested observation level: " << fObservationLevel
           << " does not exist (max obs. level: "
           << header.fObservationLevels << "), "
              "switching to level 1.";
      fObservationLevel = 1;
      INFO(info);
    }

    if (fIsThinned && header.fWMaxEM / header.fWMaxHadronic < 100) {
      ostringstream msg;
      msg << "The ratio of EM to hadron maximum weight is " << header.fWMaxEM / header.fWMaxHadronic << ". "
             "This is smaller than required to separate the EM contribution from local hadrons and "
             "the rest of the cascade.";
      WARNING(msg);
    }

    if (header.fVersion >= 6.9 && !shower.HasAtmosphereParameters()) {
      shower.MakeAtmosphereParameters(
        AtmosphereParameters(fAtmPars.fHLAY, fAtmPars.fAATM, fAtmPars.fBATM, fAtmPars.fCATM)
      );
    }

    if (!haveGH) {
      haveGH = true;
      const double a = trailer.fLongitudinalPar[3] * g/cm2;
      const double b = trailer.fLongitudinalPar[4];
      const double c = trailer.fLongitudinalPar[5] / (g/cm2);
      if (hasValidGHfit) {
        if (fIsSlant) {
          gh.SetXMax(xmax, 0);
          gh.SetNMax(nmax, 0);
          gh.SetXZero(x0, 0);
        } else {
          gh.SetXMax(xmax / fabs(cosZenith), 0); //cosZenith<0 for Zenith>90
          gh.SetNMax(nmax, 0);
          gh.SetXZero(x0 / cosZenith, 0);
        }
        gh.SetA(a, 0);
        gh.SetB(b, 0);
        gh.SetC(c, 0);
        gh.SetChiSquare(trailer.fChi2, 1);
      } else {
        // FS: Changed ERROR to WARNING, the GH Fit is read out from the .long file in case invalid values are detected
        // For parallel corsika simulations the GH Fit will be always invalid in the DAT* file
        ostringstream warn;
        warn << "CORISKA shower with invalid GH fit: Xmax=" << xmax / cosZenith / (g/cm2) << " g/cm2, "
                "Nmax=" << nmax << ", X0=" << x0 / cosZenith / (g/cm2) << " g/cm2, "
                "chi2/ndf=" << trailer.fChi2;
        WARNING(warn);
      }
    }
  }


  Status
  CorsikaShowerFile::Read(evt::Event& event)
  {
    if (event.HasSimShower()) {
      ERROR("Event not cleared - has SimShower. Cannot read CORSIKA.");
      return eFail;
    }

    bool haveGH = false;
    GaisserHillas6Parameter gh;
    bool setSimShowerAlready = false;
    double timeShift = 0;
    double zenith = 0;

    if (fRawFile.IsOpen()) { // if DAT file open

      if (fCurrentPosition >= fPositionToRaw.size()) {
        INFO("End of file");
        return eEOF;
      }
      fRawFile.SeekTo(fPositionToRaw[fCurrentPosition]);

      Corsika::EventHeader header;

      if (fIsThinned) {
        Corsika::Block headerBlock;
        if (!fRawFile.GetNextBlock(headerBlock)) {
          ostringstream err;
          err << "Cannot read CORSIKA shower header for position "
              << fCurrentPosition;
          FATAL(err);
          return eFail;
        }
        if (!headerBlock.IsEventHeader()) {
          ostringstream err;
          err << "First block at position " << fCurrentPosition
              << " is not event header";
          FATAL(err);
          return eFail;
        }
        header = headerBlock.AsEventHeader();

        fRawFile.SeekTo(fPositionOfEventTrailer[fCurrentPosition]);
        Corsika::Block trailerBlock;
        if (!fRawFile.GetNextBlock(trailerBlock)) {
          ostringstream err;
          err << "Cannot read CORSIKA shower trailer for position "
              << fCurrentPosition;
          FATAL(err);
          return eFail;
        }
        if (!trailerBlock.IsEventTrailer()) {
          ostringstream err;
          err << "Block at position " << fCurrentPosition
              << " is not event trailer";
          FATAL(err);
          return eFail;
        }
        const Corsika::EventTrailer& trailer = trailerBlock.AsEventTrailer();

        if (!setSimShowerAlready) {
          SetHeaderTrailer(event, timeShift, gh, haveGH, zenith, header, trailer);
          setSimShowerAlready = true;
        }

      } else { // switch to UN-THINNED version of DAT files

        Corsika::BlockUnthinned headerBlock;
        if (!fRawFile.GetNextBlockUnthinned(headerBlock)) {
          ostringstream msg;
          msg << "Cannot read CORSIKA shower header for position "
              << fCurrentPosition;
          FATAL(msg);
          return eFail;
        }
        if (!headerBlock.IsEventHeader()) {
          ostringstream msg;
          msg << "First block at position " << fCurrentPosition
              << " is not event header";
          FATAL(msg);
          return eFail;
        }
        header = headerBlock.AsEventHeader();

        fRawFile.SeekTo(fPositionOfEventTrailer[fCurrentPosition]);
        Corsika::BlockUnthinned trailerBlock;
        if (!fRawFile.GetNextBlockUnthinned(trailerBlock)) {
          ostringstream msg;
          msg << "Cannot read CORSIKA shower trailer for position "
              << fCurrentPosition;
          FATAL(msg);
          return eFail;
        }
        if (!trailerBlock.IsEventTrailer()) {
          ostringstream msg;
          msg << "BlockUnthinned at position " << fCurrentPosition
              << " is not event trailer";
          FATAL(msg);
          return eFail;
        }
        const Corsika::EventTrailer& trailer = trailerBlock.AsEventTrailer();

        if (!setSimShowerAlready) {
          SetHeaderTrailer(event, timeShift, gh, haveGH, zenith, header, trailer);
          setSimShowerAlready = true;
        }

      } // end thinned/unthinned block

      delete fParticleIterator;
      fParticleIterator =
        new CorsikaShowerFileParticleIterator(fRawFile,
                                              header,
                                              fAtmPars,
                                              timeShift,
                                              fPositionToRaw[fCurrentPosition] + 1,
                                              fObservationLevel,
                                              fIsThinned,
                                              header.fWMaxEM);

      ShowerSimData& shower = event.GetSimShower();
      if (!shower.HasGroundParticles())
        shower.MakeGroundParticles();
      shower.GetGroundParticles().SetFileInterface(fParticleIterator);

    } // DAT file reading is finished

    for (auto& it : fCherenkovFile) {

      const AttributeMap& am = it.first;
      Corsika::RawFile& file = it.second;

      if (!file.IsOpen()) {
        FATAL("Try to read data from file that is not open.");
        return eFail;
      }

      if (fCurrentPosition >= fPositionToCherenkovRaw[am].size()) {
        INFO("End of file (Cherenkov)");
        return eEOF;
      }
      file.SeekTo(fPositionToCherenkovRaw[am][fCurrentPosition]);

      Corsika::EventHeader header;

      if (fIsThinned) { // thinning version
        Corsika::Block headerBlock;
        if (!file.GetNextBlock(headerBlock)) {
          ostringstream err;
          err << "Cannot read Cherenkov CORSIKA shower header for position "
              << fCurrentPosition;
          FATAL(err);
          return eFail;
        }
        if (!headerBlock.IsEventHeader()) {
          ostringstream err;
          err << "First block at position " << fCurrentPosition
              << " is not cherenkov event header";
          FATAL(err);
          return eFail;
        }
        header = headerBlock.AsEventHeader();

        file.SeekTo(fPositionOfCherenkovEventTrailer[am][fCurrentPosition]);
        Corsika::Block trailerBlock;
        if (!file.GetNextBlock(trailerBlock)) {
          ostringstream err;
          err << "Cannot read CORSIKA shower cherenkov trailer for position "
              << fCurrentPosition;
          FATAL(err);
          return eFail;
        }
        if (!trailerBlock.IsEventTrailer()) {
          ostringstream err;
          err << "Block at position " << fCurrentPosition
              << " is not cherenkov event trailer";
          FATAL(err);
          return eFail;
        }
        const Corsika::EventTrailer& trailer = trailerBlock.AsEventTrailer();

        if (!setSimShowerAlready) {
          SetHeaderTrailer(event, timeShift, gh, haveGH, zenith, header, trailer);
          setSimShowerAlready = true;
        }

      } else { // switch to UN-THINNED version of Cherenkov files

        Corsika::BlockUnthinned headerBlock;
        if (!file.GetNextBlockUnthinned(headerBlock)) {
          ostringstream err;
          err << "Cannot read Cherenkov CORSIKA shower header for position "
              << fCurrentPosition;
          FATAL(err);
          return eFail;
        }
        if (!headerBlock.IsEventHeader()) {
          ostringstream err;
          err << "First block at position " << fCurrentPosition
              << " is not chernkov event header";
          FATAL(err);
          return eFail;
        }
        header = headerBlock.AsEventHeader();

        file.SeekTo(fPositionOfEventTrailer[fCurrentPosition]);

        Corsika::BlockUnthinned trailerBlock;
        if (!file.GetNextBlockUnthinned(trailerBlock)) {
          ostringstream err;
          err << "Cannot read CORSIKA shower cherenkov trailer for position "
              << fCurrentPosition;
          FATAL(err);
          return eFail;
        }
        if (!trailerBlock.IsEventTrailer()) {
          ostringstream err;
          err << "Block at position " << fCurrentPosition
              << " is not cherenkov event trailer";
          FATAL(err);
          return eFail;
        }
        const Corsika::EventTrailer& trailer = trailerBlock.AsEventTrailer();

        if (!setSimShowerAlready) {
          SetHeaderTrailer(event, timeShift, gh, haveGH, zenith, header, trailer);
          setSimShowerAlready = true;
        }

      }

      if (fCherenkovIterator.count(am))
        delete fCherenkovIterator[am];

      fCherenkovIterator[am] =
        new CorsikaShowerFileParticleIterator(file,
                                              header,
                                              fAtmPars,
                                              timeShift,
                                              fPositionToCherenkovRaw[am][fCurrentPosition] + 1,
                                              1,
                                              fIsThinned,
                                              header.fWMaxEM,
                                              true, fCherenkovBandwidthMin, fCherenkovBandwidthMax);
      ShowerSimData& shower = event.GetSimShower();
      if (!shower.HasGroundCherenkov(am))
        shower.MakeGroundCherenkov(am);
      shower.GetGroundCherenkov(am).SetFileInterface(fCherenkovIterator[am]);

    } // end loop Cherenkov files

    // check if ShowerSimData is already defined (from header/trailer from binary files)
    if (!setSimShowerAlready) {
      if (fReadSteerFile) {
        setSimShowerAlready = true;
        const Status status = ReadSteerFile(event);
        if (status != eSuccess) {
          INFO("Problem reading steering input");
          return status;
        }
      }
    }

    // if still no shower-sim-data: --> Default geometry producer
    // BUT this is not 100% filled. If the user does not use special module-sequence this will produce crab...
    if (!setSimShowerAlready) {
      WARNING("No geometry information found while reading CORSIKA. ShowerSimData is filled only partially!");
      //const double field_declination = 0; // don't know better here
      const double fieldDeclination = fUseRealisticMagneticDeclination ? 0 : fMagneticDeclination;
      const double corsikaRotation = Corsika::CorsikaAzimuthToAuger(fCoordinateRotation, fieldDeclination);
      event.MakeSimShower(DefaultShowerGeometryProducer(corsikaRotation, fUseRealisticMagneticDeclination));
    }

    ShowerSimData& shower = event.GetSimShower();

    if (fReadLongFile) {
      const Status status = ReadProfile(shower, zenith); // only needed for non-SLANT: avoid
      if (status != eSuccess) {
        INFO("End of file (profile)");
        return status;
      }
    }

    if (!shower.HasGHParameters()) {
      if (!haveGH) {
        ERROR("Found no (valid) GH fit in DAT file! No profile data available!");
      } else {

        const double xmax = gh.GetXMax();
        const double nmax = gh.GetNMax();
        const double x0 = gh.GetXZero();

        ostringstream msg;
        /*msg << "Take GH parameters from CORSIKA event trailer!\n"
          "CORSIKA SLANT option cannot be recognized! Use .long file if you use SLANT showers!";
          WARNING(msg);
          msg.str("");*/
        msg << "Nmax = " << nmax << ", "
               "Xmax = " << xmax / (g/cm2) << " g/cm2, "
               "X0 = " << x0 / (g/cm2) << " g/cm2";
        INFO(msg);
        shower.MakeGHParameters(gh);
      }
    }

    ++fCurrentPosition;
    return eSuccess;
  }


  Status
  CorsikaShowerFile::ReadSteerFile(evt::Event& event)
  {
    ifstream corsikaSteer(fSteerFile.c_str());
    if (!corsikaSteer.is_open()) {
      ostringstream msg;
      msg << "Cannot open file " << fSteerFile << " - ERROR";
      ERROR(msg);
      throw io::CorsikaIOException(msg.str());
    }

    string corsikaDataFileBase;
    string runNo;
    double magneticFieldInclination = 0;
    double magneticFieldStrength = 0;
    double showerAzimuth = 0;
    double showerZenith = 0;
    int primaryParticleType = 0;
    double primaryEnergy = 0;
    double observationHeight = 0;

    string currLine;
    while (getline(corsikaSteer, currLine)) {

      if (currLine.empty())
        continue;

      const char firstChar = currLine.at(0);
      // should check for more variants of commenting out options
      if (firstChar == 'c' || firstChar == 'C')
        continue;

      // not a comment line, parse it
      stringstream ss(currLine);

      if (Contains(currLine, "RUNNR")) {

        string skip;
        int rn = 0;
        if (!(ss >> skip >> rn))
          throw io::CorsikaIOException("RUNNR parse failed");
        corsikaDataFileBase = (format("DAT%06d") % rn).str();
        runNo = boost::lexical_cast<string>(rn);

      } else if (Contains(currLine, "MAGNET")) {

        string skip;
        double bx = 0;
        double bz = 0;
        if (!(ss >> skip >> bx >> bz))
          throw io::CorsikaIOException("MAGNET parse failed");
        bx *= micro * tesla;
        bz *= micro * tesla;
        magneticFieldInclination = atan2(bz, bx);
        magneticFieldStrength = sqrt(bx*bx + bz*bz);

      } else if (Contains(currLine, "PHIP")) {

        // get azimuth angle range, can only be used if singular!
        std::string skip;
        double az2 = 0;
        if (!(ss >> skip >> showerAzimuth >> az2))
          throw io::CorsikaIOException("PHIP parse failed");
        if (showerAzimuth != az2) {
          const string msg =
            "Cannot use CORSIKA azimuth because azimuth angle range in CORSIKA parameter file is not singular";
          ERROR(msg);
          throw io::CorsikaIOException(msg);
        }
        showerAzimuth *= utl::degree;

      } else if (Contains(currLine, "THETAP")) {

        // get zenith angle range, can only be used if singular!
        std::string skip;
        double zen2 = 0;
        if (!(ss >> skip >> showerZenith >> zen2))
          throw io::CorsikaIOException("THETAP parse failed");
        if (showerZenith != zen2) {
          const string msg =
            "Cannot use CORSIKA zenith because zenith angle range in CORSIKA parameter file is not singular";
          ERROR(msg);
          throw io::CorsikaIOException(msg);
        }
        showerZenith *= utl::degree;

      } else if (Contains(currLine, "ERANGE")) {

        // get energy range, can only be used if singular!
        std::string skip;
        double en2 = 0;
        if (!(ss >> skip >> primaryEnergy >> en2))
          throw io::CorsikaIOException("ERANGE parse failed");
        if (primaryEnergy != en2) {
          const string msg =
            "Cannot use CORSIKA energy because energy angle range in CORSIKA parameter file is not singular";
          ERROR(msg);
          throw io::CorsikaIOException(msg);
        }
        primaryEnergy *= utl::GeV;

      } else if (Contains(currLine, "PRMPAR")) {

        std::string skip;
        if (!(ss >> skip >> primaryParticleType))
          throw io::CorsikaIOException("PRMPAR parse failed");

      } else if (Contains(currLine, "ARRANG")) {

        // get corsika array rotation
        std::string skip;
        double skip2 = 0;
        if (!(ss >> skip >> fCoordinateRotation >> skip2))
          throw io::CorsikaIOException("ARRANG parse failed");
        fCoordinateRotation *= utl::degree;

      } else if (Contains(currLine, "OBSLEV")) {

        // read out the observation level; will use last observation level specified in the steering file
        std::string skip;
        if (!(ss >> skip >> observationHeight))
          throw io::CorsikaIOException("OBSLEV parse failed");
        observationHeight *= utl::cm;

      }

    }

    // geometry producer from .inp steering file -> no ground particles, no multi-core, etc.
    const bool multiCore = false;
    const CoordinateSystemPtr& refCoordinateSystem = fwk::CoordinateSystemRegistry::Get(fRefCoordinateSystemName);
    event.MakeSimShower(
      io::CorsikaShowerFileGeometryProducer(
        showerZenith, showerAzimuth + fCoordinateRotation,
        observationHeight,
        multiCore,
        refCoordinateSystem,
        fMagneticDeclination,
        fUseRealisticMagneticDeclination
      )
    );
    /*
#warning RU: check geometry here
    // calculate shower aziumth relative to geographic north (correct for magnetic field declination)
    // and normalize it to [0,2pi)
    showerAzimuth -= magneticfielddeclination;
    if (showerAzimuth < 0)
      showerAzimuth += 2 * kPi;
    // the actual setting of the shower azimuth has to be done after CORSIKA::Read() because
    // the CORSIKA reader has a hard-coded value for the magnetic field declination!
    // CORSIKA read-in is done, now set the shower azimuth
    */

    auto& shower = event.GetSimShower();
    shower.SetPrimaryParticle(Corsika::CorsikaToPDG(primaryParticleType));
    shower.SetEnergy(primaryEnergy);
    shower.SetShowerRunId(runNo);
    shower.SetMagneticFieldInclination(magneticFieldInclination);
    shower.SetMagneticFieldStrength(magneticFieldStrength);
    shower.SetGroundParticleCoordinateSystemZenith(showerZenith);
    shower.SetGroundParticleCoordinateSystemAzimuth(
      showerAzimuth
      + 180*degree  // incoming vs. outgoing
      - fCoordinateRotation  // eventual array rotation
    );

    return eSuccess;
  }


  // check existance
  Status
  CorsikaShowerFile::ScanSteerFile(const string& name)
  {
    if (!boost::filesystem::exists(name))
      return eFail;
    return eSuccess;
  }


  // check integrity of file, get number of events
  Status
  CorsikaShowerFile::ScanLongFile(const string& name)
  {
    ifstream longFile(name);
    if (!longFile.is_open())
      return eFail;

    fNumberOfEvents = 0;

    string line;
    while (getline(longFile, line))
      if (Contains(line, "LONGITUDINAL DISTRIBUTION IN"))
        ++fNumberOfEvents;

    return eSuccess;
  }


  Status
  CorsikaShowerFile::ReadProfile(evt::ShowerSimData& shower,
                                 const double zenith) // only needed for non-SLANT: avoid
  {
    //const utl::CoordinateSystemPtr localCS = shower.GetLocalCoordinateSystem(); this is not YET DEFINED here !
    //const double cosZenith = (-shower.GetDirection()).GetCosTheta(localCS); // this is only needed if not SLANT -> always use SLANT !
    const double cosZenith = cos(zenith);

    vector<string> lines;
    {
      // Read CORSIKA profile if available
      ifstream longFile(fLongFile.c_str());

      if (!longFile.is_open()) {
        const string info = "Cannot read longitudinal-profile file '" + fLongFile + "'.";
        ERROR(info);
        throw CorsikaIOException(info);
      }

      unsigned int showerN = 0;
      string line;
      while (getline(longFile, line)) {
        if (line.empty())
          continue;
        if (Contains(line, "LONGITUDINAL DISTRIBUTION IN")) {
          ++showerN;
          if (showerN > fCurrentPosition+1)
            break;
        }
        if (showerN == fCurrentPosition+1)
          lines.push_back(line);
      }
    }
    const unsigned int nLines = lines.size();

    if (nLines < 3) {
      // this is not worth an extra output
      //ostringstream err;
      //err << "Cannot find position " << fCurrentPosition << " in the longitudinal file '" << fLongFile << "'. EOF.";
      //INFO(err);
      return eEOF;
    } else {
      ostringstream info;
      info << "Read " << nLines << " lines from '" << fLongFile << "' at position " << fCurrentPosition;
      INFO(info);
    }

    typedef unsigned int (*ToUIntType)(const string&);
    const ToUIntType ToUInt = &boost::lexical_cast<unsigned int>;
    typedef double (*ToDType)(const string&);
    const ToDType ToD = &boost::lexical_cast<double>;;

    bool isProfileSlant = false;
    unsigned int profileLength = 0;
    unsigned int i = 0;
    // particle distributions
    {
      const vector<string> s1 = String::Split(lines.at(i));
      const vector<string> s2 = String::Split(lines.at(i + 1));
      if (!(s1.size() == 12 && s1[0] == "LONGITUDINAL" && s1[1] == "DISTRIBUTION" &&
            s2.size() == 10 && s2[0] == "DEPTH" && s2[1] == "GAMMAS")) {
        const string err = "Mismatch in particle distribution profile headers.";
        ERROR(err);
        throw CorsikaIOException(err);
      } else {
        const unsigned int n = ToUInt(s1[3]);
        profileLength = n;
        if (n < 1) {
          const string msg = "Particle profiles are too short.";
          ERROR(msg);
          throw CorsikaIOException(msg);
        }
        const string& type = s1[4];
        if (type != "SLANT" && type != "VERTICAL") {
          const string err = "Unknown profile type '" + type + "', (only 'VERTICAL' and 'SLANT' supported)";
          ERROR(err);
          throw CorsikaIOException(err);
        }
        isProfileSlant = (type == "SLANT");
        const double slantFac = isProfileSlant ? 1 : 1 / fabs(cosZenith);  // cosZenith < 0 for Zenith > 90
        const double dX = slantFac * ToD(s1[7]) * (g/cm2);
        i += 2;
        vector<double> particleDistributionDepth;
        vector<double> particleDistributionCharge;
        vector<double> particleDistributionMuons;
        vector<double> particleDistributionGammas;
        vector<double> particleDistributionElectrons;
        vector<double> particleProductionCherenkov;
        for (unsigned int j = 0; j < n; ++j) {
          // the Fortran formats for writing the longitudina energy deposits depend on the options
          // With the SLANT option the format is FORMAT(' ', F7.1, 1P, 9(E12.5), 0P)
          // For the standard output without SLANT option the format is FORMAT(' ', F5.0, 1P, 9(E12.5), 0P)
          vector<string> s = String::Split(lines.at(i));
          if (s.size() != 10) {
            WARNING("Corrupted longitudinal particle distribution table detected; will try to fix it.");
            s = FixedColumnWidthSplit(lines[i], 10, isProfileSlant ? 8 : 6, 12);
            ostringstream warn;
            warn << "Fixed line: ";
            string comma = "";
            for (const auto& sss : s) {
              warn << comma << sss;
              comma = ", ";
            }
            WARNING(warn);
          }
          particleDistributionDepth.push_back(slantFac * ToD(s.at(0)) * (g / cm2));
          const double ch = ToD(s.at(7));
          particleDistributionCharge.push_back(ch);
          const double mu = ToD(s.at(4)) + ToD(s.at(5));
          particleDistributionMuons.push_back(mu);
          particleDistributionGammas.push_back(ToD(s.at(1)));
          const double e = ToD(s.at(2)) + ToD(s.at(3));
          particleDistributionElectrons.push_back(e);
          {
            // this is "production per bin", thus need to normalize per bin-width
            particleProductionCherenkov.push_back(double(ToD(s.at(9)) / dX));
          }
          ++i;
        }
        const double normDepth = particleDistributionDepth.back();
        const double normN = shower.HasGHParameters() ? shower.GetGHParameters().Eval(normDepth) : 1;
        const double normCharge = particleDistributionCharge.back();
        const double normMuons = particleDistributionMuons.back();
        const double normGammas = particleDistributionGammas.back();
        const double normElectrons = particleDistributionElectrons.back();
        const double normCherenkov = particleProductionCherenkov.back();
        const double lastDepth = particleDistributionDepth.back();

        // prolongate profiles by some bins
        if (type == "SLANT" || zenith < 90*deg) {
          // do not do for VERTICAL depth in the case of up-going events.
          for (int j = 0; j < 2; ++j) {
            const double addDepth = lastDepth + dX * (j + 1);
            const double nch = shower.HasGHParameters() ? shower.GetGHParameters().Eval(addDepth) : 1;
            const double fac = nch / normN;
            const double addCharge = normCharge * fac;
            const double addMuons = normMuons * fac;
            const double addGammas = normGammas * fac;
            const double addElectrons = normElectrons * fac;
            const double addCherenkov = normCherenkov * fac;
            particleDistributionDepth.push_back(addDepth);
            particleDistributionCharge.push_back(addCharge);
            particleDistributionMuons.push_back(addMuons);
            particleDistributionGammas.push_back(addGammas);
            particleDistributionElectrons.push_back(addElectrons);
            particleProductionCherenkov.push_back(addCherenkov);
          }
        }
        if (type == "VERTICAL" && zenith > 90*deg) {
          // the vertical depth in the .long file goes from top of the
          // atmosphere(min value) to the bottom of the atmosphere(max value)
          // for UPWARD events, the min value of X_slant corresponds to the
          // bottom of the atmosphere, and max value to the top of the
          // atmosphere. So the vertical profile must be reversed for upgoing
          // events
          reverse(particleDistributionCharge.begin(), particleDistributionCharge.end());
          reverse(particleDistributionMuons.begin(), particleDistributionMuons.end());
          reverse(particleDistributionGammas.begin(), particleDistributionGammas.end());
          reverse(particleDistributionElectrons.begin(), particleDistributionElectrons.end());
          reverse(particleProductionCherenkov.begin(), particleProductionCherenkov.end());
        }
        TabulatedFunction muonProfile;
        TabulatedFunction gammaProfile;
        TabulatedFunction electronProfile;
        TabulatedFunction chargeProfile;  // charged
        TabulatedFunction cherenkovProfile;
        for (unsigned int j = 0; j < n; ++j) {
          const double depth = particleDistributionDepth[j];
          chargeProfile.PushBack(depth, particleDistributionCharge[j]);
          muonProfile.PushBack(depth, particleDistributionMuons[j]);
          gammaProfile.PushBack(depth, particleDistributionGammas[j]);
          electronProfile.PushBack(depth, particleDistributionElectrons[j]);
          cherenkovProfile.PushBack(depth, particleProductionCherenkov[j]);
        }
        if (!shower.HasLongitudinalProfile())
          shower.MakeLongitudinalProfile(chargeProfile);
        // RU Thu May 12 10:03:23 CEST 2005 (added moun profile)
        if (!shower.HasLongitudinalProfile(ShowerSimData::eMuon))
          shower.MakeLongitudinalProfile(muonProfile, ShowerSimData::eMuon);
        if (!shower.HasLongitudinalProfile(ShowerSimData::ePhoton))
          shower.MakeLongitudinalProfile(gammaProfile, ShowerSimData::ePhoton);
        if (!shower.HasLongitudinalProfile(ShowerSimData::eElectron))
          shower.MakeLongitudinalProfile(electronProfile, ShowerSimData::eElectron);
        // RU Mo 26. Sep 23:06:39 CEST 2016 (added differential Cherenkov production profile)
        if (!shower.HasLongitudinalProfile(ShowerSimData::eCherenkov))
          shower.MakeLongitudinalProfile(cherenkovProfile, ShowerSimData::eCherenkov);

      }
    }
    // energy deposit
    {
      const vector<string> s1 = String::Split(lines.at(i));
      const vector<string> s2 = String::Split(lines.at(i + 1));
      if (!(s1.size() == 13 && s1[0] == "LONGITUDINAL" && s1[1] == "ENERGY" &&
            s2.size() == 16 && s2[0] == "DEPTH" && s2[1] == "GAMMA")) {
        const string err = "Mismatch in energy deposit profile headers.";
        ERROR(err);
        throw CorsikaIOException(err);
      } else {
        const unsigned int n = ToUInt(s1[4]);
        if (n < 3) {
          const string msg = "Energy deposit profile is too short.";
          ERROR(msg);
          throw CorsikaIOException(msg);
        }
        const string& type = s1[5];
        if (type != "SLANT" && type != "VERTICAL") {
          const string err = "Unknown profile type '" + type + "', (only 'VERTICAL' and 'SLANT' supported)";
          ERROR(err);
          throw CorsikaIOException(err);
        }
        if (isProfileSlant != (type == "SLANT")) {
          const string err = "Energy deposit and particle profiles do not have the same inclination.";
          ERROR(err);
          throw CorsikaIOException(err);
        }
        const double slantFac = isProfileSlant ? 1 : 1 / fabs(cosZenith); // cosZenith < 0 for Zenith > 90
        const double dX = slantFac * ToD(s1[8]) * g/cm2;
        i += 2;
        // The energy deposit of particles hitting ground is accounted for in"
        //        the last row if not SLANT
        //        the two last rows if SLANT (in good approximation*)
        // *A few (early) particles reaching ground are accounted for in previous rows (Effect seems to be very small,
        // Tanguy thinks that could be "corrected" by counting the last row twice for groundEnergy. This is not done yet here!)
        const unsigned int groundRowsAtTheEnd = (isProfileSlant) ? 2 : 1;
        double energyDepositSum = 0;
        double electromagneticEnergyDeposit = 0;
        vector<double> energyDepositDepth;
        vector<double> energyDeposit;
        for (unsigned int j = 0; j < n; ++j) {
          vector<string> s = String::Split(lines.at(i));
          if (s.size() != 10) {
            WARNING("Corrupted longitudinal energy deposit table detected; will try to fix it.");
            s = FixedColumnWidthSplit(lines[i], 10, isProfileSlant ? 8 : 7, 12);
            ostringstream warn;
            warn << "Fixed line: ";
            string comma = "";
            for (const auto& sss : s) {
              warn << comma << sss;
              comma = ", ";
            }
            WARNING(warn);
          }
          energyDepositDepth.push_back(slantFac * ToD(s.at(0)) * (g / cm2));
          const double muCut = ToD(s.at(5)) * GeV;
          const double hadCut = ToD(s.at(7)) * GeV;
          const double nu = ToD(s.at(8)) * GeV;
          const double sum = ToD(s.at(9)) * GeV;
          // Subtracting neutrino energy and take fraction for muons and hadrons
          // Barbosa et al., Astropart. Phys. 22 (2004) 159.
          const double muonFraction = 0.575;
          const double hadronFraction = 0.;  //0.261;  // <-- Tanguy thinks this should be zero
          const double de = sum - nu - muonFraction * muCut - hadronFraction * hadCut;
          energyDeposit.push_back(de / dX);
          if (j < n - groundRowsAtTheEnd)
            energyDepositSum += de;
          else {
            const double gamma = ToD(s.at(1)) * GeV;
            const double emIon = ToD(s.at(2)) * GeV;
            const double emCut = ToD(s.at(3)) * GeV;
            const double muIon = ToD(s.at(4)) * GeV;
            const double hadIon = ToD(s.at(6)) * GeV;
            const double hadronGroundFraction = 0.390;
            const double groundEnergy =
              (1 - hadronGroundFraction) * hadCut + hadIon + muIon + emIon + emCut + gamma;
            energyDepositSum += groundEnergy;
          }

          // Adding GAMMA, EM IONIZ and EM CUT over whole atmosphere including energy deposit in ground plane
          electromagneticEnergyDeposit += (ToD(s.at(1)) + ToD(s.at(2)) + ToD(s.at(3))) * GeV;
          ++i;
        }
        const double normDepth = energyDepositDepth[n - 3];
        const double normN = shower.HasGHParameters() ? shower.GetGHParameters().Eval(normDepth) : 1;
        // go three bins back, since CORSIKA sometimes has a funny second-last bin
        const double normdEdX = energyDeposit[n - 3];
        // recalculate last dedx bins
        for (int j = 0; j < 2; ++j) {
          const int bin = n - 2 + j;
          const double depth = energyDepositDepth[bin];
          energyDeposit[bin] = normdEdX / normN *
                               (shower.HasGHParameters() ? shower.GetGHParameters().Eval(depth) : 1);
        }
        const double lastDepth = energyDepositDepth.back();
        // prolongate profiles by some bins
        if (type == "SLANT" || zenith < 90*deg) {
          // do not do for VERTICAL depth in the case of up-going events
          for (int j = 0; j < 2; ++j) {
            const double addDepth = lastDepth + dX * (j + 1);
            const double nch = shower.HasGHParameters() ? shower.GetGHParameters().Eval(addDepth) : 1;
            const double adddEdX = normdEdX / normN * nch;
            energyDepositDepth.push_back(addDepth);
            energyDeposit.push_back(adddEdX);
          }
        }
        // the vertical depth in the .long file goes from top of the
        // atmosphere(min value) to the bottom of the atmosphere(max value) for
        // UPWARD events, the min value of X_slant corresponds to the bottom of
        // the atmosphere, and max value to the top of the atmosphere. So the
        // vertical profile must be reversed for upgoing events
        if (type == "VERTICAL" && zenith > 90*deg)
          reverse(energyDeposit.begin(), energyDeposit.end());
        TabulatedFunction dEdX;
        for (unsigned int j = 0; j < n; ++j)
          dEdX.PushBack(energyDepositDepth[j], energyDeposit[j]);

        for (unsigned int i = 0; i < energyDepositDepth.size(); ++i) {
          GHFit::gProfile.fSlantDepth.push_back(energyDepositDepth[i]);
          GHFit::gProfile.fEnergyDeposit.push_back(energyDeposit[i]);
        }

        shower.SetCalorimetricEnergy(energyDepositSum);
        shower.SetElectromagneticEnergy(electromagneticEnergyDeposit);
        // BRD we need a trap here.  Only MakedEdX if such profile is not null
        if (!energyDeposit.empty() && !shower.HasdEdX())
          shower.MakedEdX(dEdX);
      }
    }

    {
      const vector<string> s1 = String::Split(i < nLines ? lines[i] : "");
      const vector<string> s2 = String::Split(i + 1 < nLines ? lines[i + 1] : "");
      if (!(s1.size() > 5 && s1[0] == "FIT" && s1[1] == "OF" && s1[2] == "THE" && s1[3] == "HILLAS" &&
            s2.size() == 7 && s2[0] == "TO" && s2[1] == "LONGITUDINAL" && s2[2] == "DISTRIBUTION")) {
        WARNING("Native GH fit not found.");
      } else {
        const vector<string> s3 = String::Split(lines.at(i + 2));
        const vector<string> s4 = String::Split(lines.at(i + 3));
        //const vector<string> s5 = String::Split(lines.at(i + 4));
        if (s3.size() != 8 || s3[0] != "PARAMETERS" || s3[1] != "=" ||
            s4.size() != 3 || s4[0] != "CHI**2/DOF" || s4[1] != "=" /*||
            s5.size() < 3 || s5[0] != "AV."*/) {
          const string err = "Malformed fit of the Hillas curve.";
          ERROR(err);
          throw CorsikaIOException(err);
        }
        // do not do for VERTICAL depth in the case of up-going events
        // parameters in CORSIKA calculated for atmosphere starting to the top and ending to the bottom
        // we need the parameters correpsonding to the reversed profiles (atmo starts at the bottom and ends at top
        if (isProfileSlant || zenith < 90*deg) {
          const double slantFac = isProfileSlant ? 1 : 1 / fabs(cosZenith); // cosZenith < 0 for Zenith > 90
          const double nMax = ToD(s3.at(2));
          const double x0 = slantFac * ToD(s3.at(3)) * g / cm2;
          const double xMax = slantFac * ToD(s3.at(4)) * g / cm2;
          const double apar = ToD(s3.at(5)) * g / cm2;
          const double bpar = ToD(s3.at(6));
          const double cpar = ToD(s3.at(7)) / (g / cm2);
          const double chi2 = ToD(s4.at(2));
          if (chi2 <= 0 || chi2 * profileLength > 1e15) {
            ostringstream err;
            err << "CORISKA shower with invalid GH fit: Xmax = " << xMax / (g/cm2) << " g/cm2, "
                   "Nmax = " << nMax << ", X0 = " << x0 / (g/cm2) << " g/cm2, chi2/ndf = " << chi2;
            ERROR(err);
          } else {
            const double a = slantFac * apar;
            const double b = slantFac * bpar;
            const double c = slantFac * cpar;
            evt::GaisserHillas6Parameter gh;
            gh.SetXMax(xMax, 0);
            gh.SetNMax(nMax, 0);
            gh.SetXZero(x0, 0);
            gh.SetChiSquare(chi2 * profileLength, profileLength);
            gh.SetA(a, 0);
            gh.SetB(b, 0);
            gh.SetC(c, 0);
            ostringstream info;
            info << "Nmax = " << nMax << ", "
                    "Xmax = " << xMax / (g/cm2) << " g/cm2, "
                    "X0 = " << x0 / (g/cm2) << " g/cm2, "
                    "zenith = " << acos(cosZenith) / deg << " deg"
                 << (isProfileSlant ? " (SLANT DEPTH)" : " (VERTICAL DEPTH)");
            INFO(info);
            if (!shower.HasGHParameters())
              shower.MakeGHParameters(gh);
          }
        } else if (!isProfileSlant && zenith > 90*deg) {
          // recalculating GHParameters for the case of up-going showers and VERTICAL atmosphere
          // so that atmosphere starts at the bottom and ends at the top. Using the parameters
          // resulted from the GH fit of N(x), as explained above

          TMinuit minuitdEdXGHFit(4);  // GH with 4 par for fit
          minuitdEdXGHFit.SetFCN(GHFit::GH4ParamFitFunctiondEdxGH);

          double arglist[2] = { 0 };
          int ierflag;
          // get the index for the corresponding depth to dEdXmax from .long data;
          const double index = std::distance(
            GHFit::gProfile.fEnergyDeposit.begin(),
            max_element(GHFit::gProfile.fEnergyDeposit.begin(), GHFit::gProfile.fEnergyDeposit.end())
          );
          arglist[0] = 1;
          minuitdEdXGHFit.mnexcm("SET ERR", arglist, 1, ierflag);
          // take expectded dEdXmax from the current .long data
          double dEdXmaxFit = *max_element(GHFit::gProfile.fEnergyDeposit.begin(), GHFit::gProfile.fEnergyDeposit.end()) / (PeV / (g/cm2));
          // first guess of xmax
          double xMaxFit = GHFit::gProfile.fSlantDepth[index] / (g/cm2);
          // initial guess for X0 and lambda
          double x0Fit = 10;
          double lambdaFit = 70;

          minuitdEdXGHFit.mnparm(0, "dEdXmax", dEdXmaxFit, 0.1 * fabs(dEdXmaxFit), 0, 0, ierflag);
          minuitdEdXGHFit.mnparm(1, "X0", x0Fit, 1, -1000, GHFit::gProfile.fSlantDepth[0] / (g/cm2), ierflag);
          minuitdEdXGHFit.mnparm(2, "Xmax", xMaxFit, 20, 0, 0, ierflag);
          minuitdEdXGHFit.mnparm(3, "lambda", lambdaFit, 5, 10, 150, ierflag);

          arglist[0] = 1000;
          arglist[1] = 1;

          minuitdEdXGHFit.mnexcm("MIGRAD", arglist, 2, ierflag);

          double amin = 0;
          double edm = 0;
          double errdef = 0;
          int nvpar = 0;
          int nparx = 0;
          int icstat = 0;
          minuitdEdXGHFit.mnstat(amin, edm, errdef, nvpar, nparx, icstat);

          double p[4] = { 0 };
          minuitdEdXGHFit.GetParameter(0, dEdXmaxFit, p[0]);
          minuitdEdXGHFit.GetParameter(1, x0Fit, p[1]);
          minuitdEdXGHFit.GetParameter(2, xMaxFit, p[2]);
          minuitdEdXGHFit.GetParameter(3, lambdaFit, p[3]);

          // get Nmax using xmax and dEdXmax provided from dEdX GHfit
          // utl::EnergyDeposit (const double depthX, const double xMax, const double enCut) calculated in PhysicalFunctions.cc
          // it is really the mean ionization loss rate. Parameters are for energies in MeV / (g/cm2)

          const double nMaxFit =
              dEdXmaxFit * PeV / MeV /
              (utl::EnergyDeposit(xMaxFit * (g / cm2), xMaxFit * (g/cm2), GHFit::gProfile.fElectronCutoff / MeV) / (MeV / (g/cm2)));

          const double x0 = x0Fit * g/cm2;
          const double xMax = xMaxFit * g/cm2;
          //getting lambda from a GH4 fit. Let it with GH6Fit for consistency. However by making a =lambda,b=0, c=0 nothing changes
          const double apar = lambdaFit * g/cm2;
          // b and c param only used in GaisserHillas6Parameters.cc, Eval(depth) - l. 28:
          // lambda = fA + depth*fB + Sqr(depth)*fC;
          const double bpar = 0;
          const double cpar = 0;
          const double chi2 = minuitdEdXGHFit.fAmin;
          if (chi2 <= 0 || chi2 * profileLength > 1e15) {
            ostringstream err;
            err << "CORISKA shower with invalid GH fit: Xmax = " << xMax / (g/cm2) << " g/cm2, "
                   "Nmax = " << nMaxFit << ", X0 = " << x0 / (g/cm2) << " g/cm2, chi2/ndf = " << chi2;
            ERROR(err);
          } else {
            const double a = apar;
            const double b = bpar;
            const double c = cpar;
            evt::GaisserHillas6Parameter gh;
            gh.SetXMax(xMax, 0);
            gh.SetNMax(nMaxFit, 0);
            gh.SetXZero(x0, 0);
            gh.SetChiSquare(chi2 * profileLength, profileLength);
            gh.SetA(a, 0);
            gh.SetB(b, 0);
            gh.SetC(c, 0);
            ostringstream info;
            info << "Nmax = " << nMaxFit << ", "
                    "Xmax = " << xMax / (g/cm2) << " g/cm2, "
                    "X0 = " << x0 / (g/cm2) << " g/cm2, "
                    "zenith = " << acos(cosZenith) / deg << " deg"
                 << (isProfileSlant ? " (SLANT DEPTH)" : " (VERTICAL DEPTH)");
            INFO(info);
            if (!shower.HasGHParameters())
              shower.MakeGHParameters(gh);
          }
        }
      }
    }

    // CORSIKA writes into the last bin of the dEdX profile not only
    // the energy deposit in the atmosphere, but also the energy hitting
    // the ground (observation level). For vertical event this can be
    // remarkable, and should not be used for FD simulations!
    //
    // - Fix this by replacing the entries in the last bin by GH fitted
    //   extrapolations.
    // - Add one more bin in depth (extrapolation) to make sure
    //   that also for a different atmospheric profile model the
    //   shower will hit the ground, and does not just disappear
    //   somewhere in the middle of the atmosphere!

    return eSuccess;
  }


  void
  CorsikaShowerFile::Write(const evt::Event& /*event*/)
  {
    const string msg = "Cannot write to CorsikaShowerFile - read-only file";
    ERROR(msg);
    throw CorsikaIOException(msg);
  }


  Status
  CorsikaShowerFile::FindEvent(const unsigned int eventId)
  {
    const auto it = fIdToPosition.find(eventId);
    if (it == fIdToPosition.end())
      return eFail;
    return GotoPosition(it->second);
  }


  Status
  CorsikaShowerFile::GotoPosition(const unsigned int position)
  {
    if (position >= fNumberOfEvents)
      return eFail;
    fCurrentPosition = position;
    return eSuccess;
  }


  int
  CorsikaShowerFile::GetNEvents()
  {
    if (fOpenedAnything)
      return fNumberOfEvents;

    const string msg = "Cannot request number of events from closed file";
    ERROR(msg);
    throw CorsikaIOException(msg);
  }


  CorsikaShowerFileParticleIterator*
  CorsikaShowerFile::GetCorsikaShowerFileParticleIterator()
    const
  {
    return fParticleIterator;
  }


  double
  CorsikaShowerFile::GetDensityVsHeight(const double height)
  const
  {
    for (int i = 0; i < 4; ++i)
      if (fAtmPars.fHLAY[i] <= height && height < fAtmPars.fHLAY[i + 1])
        return fAtmPars.fBATM[i] * exp(-height / fAtmPars.fCATM[i]) / fAtmPars.fCATM[i];
    if (height >= fAtmPars.fHLAY[4])
      return fAtmPars.fBATM[4] / fAtmPars.fCATM[4];
    return 0;
  }


  double
  CorsikaShowerFile::GetDepthVsHeight(const double height)
  const
  {
    //figure out in what layer of the atmosphere model the first interaction happened
    int section = 0;
    //#warning why use fabs(height)?
    if (fabs(height) > fAtmPars.fHLAY[4])
      section = 4;
    else if (fabs(height) > fAtmPars.fHLAY[3])
      section = 3;
    else if (fabs(height) > fAtmPars.fHLAY[2])
      section = 2;
    else if (fabs(height) > fAtmPars.fHLAY[1])
      section = 1;

    //Calculate X1 (slant depth in g/cm^2) from height (in cm) of first interaction
    if (section < 4)
      return fAtmPars.fAATM[section] + fAtmPars.fBATM[section] * exp(-height / fAtmPars.fCATM[section]);
    return fAtmPars.fAATM[section] - fAtmPars.fBATM[section] * (height / fAtmPars.fCATM[section]);
  }


  double
  CorsikaShowerFile::GetSlantDepthCurved(const double hObs, const double theta, const double h)
  const
  {
    using namespace Corsika;
    const double sinZenith = sin(theta);
    const double sinZenith2 = sinZenith * sinZenith;

    // the resolution
    const double dX = 0.01 * (g/cm2);
    double height = h;
    double depth = 0;

    while (height < hAtmBoundary && height > hObs) {
      height += dX * sqrt(1 - Sqr((Corsika::kEarthRadius + hObs) / (Corsika::kEarthRadius + height)) * sinZenith2) /
                GetDensityVsHeight(height);
      depth += dX;
    }
    return depth;
  }

}