CorsikaShowerFileParticleIterator.cc

Go to the documentation of this file.

#include <io/CorsikaShowerFileParticleIterator.h>

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

#include <evt/ShowerSimData.h>

#include <utl/AugerUnits.h>
#include <utl/CoordinateSystem.h>
#include <utl/ErrorLogger.h>
#include <utl/Particle.h>
#include <utl/PhysicalConstants.h>
#include <utl/Point.h>
#include <utl/Vector.h>
#include <utl/GeometryUtilities.h>
#include <utl/ReferenceEllipsoid.h>
#include <utl/UTMPoint.h>
#include <utl/Math.h>

#include <utl/RandomEngine.h>
#include <fwk/LocalCoordinateSystem.h>
#include <fwk/RandomEngineRegistry.h>

#include <CLHEP/Random/Randomize.h>

#include <det/Detector.h> // to retrieve references-CS

#include <sstream>
#include <cmath>

using namespace evt;
using namespace utl;
using namespace std;

using namespace io;
using ::io::Corsika::Block;
using ::io::Corsika::BlockUnthinned;

using CLHEP::RandFlat;


CorsikaShowerFileParticleIterator::
CorsikaShowerFileParticleIterator(Corsika::RawFile& rawFile,
                                  const io::Corsika::EventHeader& header,
                                  const Corsika::AtmosphereParameters& params,
                                  const double timeOffset,
                                  const PositionType& startPosition,
                                  const unsigned int observationLevel,
                                  const bool isThinned,
                                  const double wmaxEM,
                                  const bool CERFIL,
                                  const double CERwlMin,
                                  const double CERwlMax) :
  fRawFile(rawFile),
  fAtmPars(params),
  fTimeOffset(timeOffset),
  fStartPosition(startPosition),
  fVersion(header.fVersion),
  fObservationLevel(observationLevel),
  fCurrentPosition(startPosition),
  fCurrentParticle(utl::Particle::eUndefined,
                   utl::Particle::eShower,
                   utl::Point(0, 0, 0, CoordinateSystem::GetRootCoordinateSystem()),
                   utl::Vector(0, 0, 0, CoordinateSystem::GetRootCoordinateSystem()),
                   utl::TimeInterval(0),
                   0.),
  fIsThinned(isThinned),
  fWMaxEM(wmaxEM),
  fCERFIL(CERFIL),
  fMCERFI(header.fCerenkovOutputFlag),
  fWlLowerLimit(CERwlMin),
  fWlUpperLimit(CERwlMax),
  fCurvedObsLevel(header.fCurvedObsLevel),
  fCurved(header.fFlagCurved),
  fObsLevelHeight(header.fObservationHeight[observationLevel-1]*cm),
  fParent(utl::Particle(utl::Particle::eUndefined,
                        utl::Particle::eShower,
                        utl::Point(0, 0, 0, CoordinateSystem::GetRootCoordinateSystem()),
                        utl::Vector(0, 0, 0, CoordinateSystem::GetRootCoordinateSystem()),
                        utl::TimeInterval(0),
                        0.)),
  fGrandParent(utl::Particle(utl::Particle::eUndefined,
                             utl::Particle::eShower,
                             utl::Point(0, 0, 0, CoordinateSystem::GetRootCoordinateSystem()),
                             utl::Vector(0, 0, 0, CoordinateSystem::GetRootCoordinateSystem()),
                             utl::TimeInterval(0),
                             0.))
{
  fGrandParent.SetValid(false);
  fParent.SetValid(false);
}


void
CorsikaShowerFileParticleIterator::Rewind()
{
  fCurrentPosition = fStartPosition;
  fParticleInBlock = 0;
  fBlockBufferValid = false;
  fIteratorValid = true;
}


Particle::Source
CorsikaShowerFileParticleIterator::GetEMParticleSource(const float wt, const short int gen)
  const
{
  if (fVersion < 6.71)
    return Particle::eShower;
  else if (gen > 50)
    return Particle::eShowerFromMuonDecay;
  else if (fIsThinned) {
    // A rectangular cut
    const double wlimits[2] = { 0, 2000.*fWMaxEM/10000 };
    const double glimits[2] = { 5, 11 };
    const double w = double(wlimits[1] - wlimits[0]) / (glimits[1] - glimits[0]) * (gen - glimits[0]) + wlimits[0];
    if (wt < w)
      return Particle::eShowerFromLocalHadron;
    return Particle::eShower;
  }
  return Particle::eShower;
}


Particle*
CorsikaShowerFileParticleIterator::GetOneParticle(const utl::CoordinateSystemPtr& groundCS)
{
  for (;;) {

    if (fIsThinned) {

      const Corsika::ParticleData* const corsikaParticle = GetOneParticleRecord();

      // end of particle list
      if (!corsikaParticle)
        return 0;

      if (fCERFIL) {

        const Corsika::CherenkovData* const cherenkovPhoton = (const Corsika::CherenkovData*)corsikaParticle;

        const double weight = cherenkovPhoton->fPhotonsInBunch *
          (fMCERFI == 1 ? cherenkovPhoton->fWavelength : 1); // "Wavelength" is still "thinning weight"

        if (weight == 0 || (cherenkovPhoton->fX == 0 && cherenkovPhoton->fY == 0)) // empty block is all zeros
          return 0;

        const double x = cherenkovPhoton->fX*cm;
        const double y = cherenkovPhoton->fY*cm;
        const double z = 0;
        /*
          Note: Cherenkov photon output is NEVER on curved observation level !!!!
                RU, Do 15. Dez 15:49:12 CET 2016
        */
        const double sintheta_cosphi = cherenkovPhoton->fU;
        const double sintheta_sinphi = cherenkovPhoton->fV;
        const double costheta = -sqrt(1 - Sqr(sintheta_sinphi) - Sqr(sintheta_cosphi)); // downwards

        const Point posOnGround(x, y, z, groundCS);
        const Vector direction(sintheta_cosphi, sintheta_sinphi, costheta, groundCS);

        double lambda = (fMCERFI > 1) ? cherenkovPhoton->fWavelength * nanometer : 0; // for MCERFI==1 no wavelength is stored!
        if (!lambda) {
          // zero also means that CORSIKA did not generate wavelength already
          utl::RandomEngine* const randomEngine = &fwk::RandomEngineRegistry::GetInstance().Get(fwk::RandomEngineRegistry::eDetector);
          lambda = 1 / (1/fWlLowerLimit - RandFlat::shoot(&randomEngine->GetEngine(), 0, 1) * (1/fWlLowerLimit - 1/fWlUpperLimit));
        }

        // for MCERFI<=2, fDistance is PRODUCTION HEIGHT, for MCERFI>2 it is DISTANCE
        const double distance = cherenkovPhoton->fDistance * (fMCERFI > 2 ? 1. : -1/costheta) * cm;
        const Point posAtProduction = posOnGround - distance*direction;

        const double energy = utl::kPlanck / lambda;

        fCurrentParticle =
          Particle(utl::Particle::ePhoton,
                   Particle::eShower,
                   posAtProduction, direction,
                   TimeInterval(cherenkovPhoton->fT*ns - fTimeOffset - distance/utl::kSpeedOfLight),
                   weight, energy);

        fProductionPoint.SetValid(false);
        fGrandParent.SetValid(false);
        fParent.SetValid(false);

        return &fCurrentParticle;
      }

      const int info = corsikaParticle->fDescription / 1000;
      const int particleId = Corsika::CorsikaToPDG(abs(info));

      // muon additional info (75,76: regular muons. 85,86: decaying muons at production point)
      if (info == 75 || info == 76 || info == 85 || info == 86) {
        if (!fCurved) {
          fProductionPoint = Point(corsikaParticle->fX*cm,
                                   corsikaParticle->fY*cm,
                                   corsikaParticle->fTorZ*cm - fObsLevelHeight,
                                   groundCS);
          continue;
        }
      }

      // negative particle description denotes mother and grandmother particles
      if (corsikaParticle->fDescription < 0) {

        if (particleId == Particle::eUndefined) {
          ostringstream warn;
          warn << "Unknown corsika particle code " << abs(int(corsikaParticle->fDescription/1000));
          WARNING(warn);
        }

        if (!fParent.IsValid()) {
          if (!fCurved) {
            fParent = Particle(particleId,
                               Particle::eShower,
                               Point(corsikaParticle->fX*cm,
                                     corsikaParticle->fY*cm,
                                     corsikaParticle->fTorZ*cm - fObsLevelHeight,
                                     groundCS),
                               Vector(corsikaParticle->fPx*GeV,
                                      corsikaParticle->fPy*GeV,
                                      -corsikaParticle->fPz*GeV,
                                      groundCS),
                               TimeInterval(),
                               corsikaParticle->fWeight);
          }
        } else if (!fGrandParent.IsValid()) {
          if (!fCurved) {
            const Vector direction(corsikaParticle->fPx*GeV,
                                   corsikaParticle->fPy*GeV,
                                   -corsikaParticle->fPz*GeV,
                                   groundCS);
            const double cosTheta = direction.GetCosTheta(groundCS);
            const double parent_h = corsikaParticle->fTorZ*cm;
            // Subtract traversed depth along the momentum direction
            const double grand_parent_depth = GetDepthAtHeight(parent_h) - corsikaParticle->fY*g/cm2*cosTheta;
            const double d = (GetHeightAtDepth(grand_parent_depth) - parent_h)/cosTheta;
            const Point position = fParent.Get().GetPosition() - d*direction/direction.GetMag();
            fGrandParent = Particle(particleId,
                                    Particle::eShower,
                                    position,
                                    direction,
                                    TimeInterval(),
                                    corsikaParticle->fWeight);
          }
        }
        continue;
      }

      // skip unknown particles...
      if (particleId == Particle::eUndefined) {
        fProductionPoint.SetValid(false);
        fGrandParent.SetValid(false);
        fParent.SetValid(false);
        continue;
      }

      // these are not defined in ehistory particles
      const short int hadGen = fmod(corsikaParticle->fDescription, 1000) / 10;
      short unsigned int obsLevel = fmod(corsikaParticle->fDescription, 10);
      //short unsigned int fateIndex = 0;

      // 95,96: decaying muons at point of decay. observation level is replaced by fate index
      if (info == 95 || info == 96) {
        obsLevel = fObservationLevel;
        //fateIndex = fmod(corsikaParticle->fDescription, 10);
      }

      if (obsLevel != fObservationLevel) {
        fProductionPoint.SetValid(false);
        fGrandParent.SetValid(false);
        fParent.SetValid(false);
        continue;
      }

      const float weight = corsikaParticle->fWeight;

      const Particle::Source source =
        (particleId == Particle::eElectron ||
         particleId == Particle::ePositron ||
         particleId == Particle::ePhoton) ?
         GetEMParticleSource(weight, hadGen) : Particle::eShower;

      double x = corsikaParticle->fX*cm;
      double y = corsikaParticle->fY*cm;
      double z = 0;
      if ((fVersion == 6.98) ||
          (fVersion >= 6.981 && fVersion < 6.990 && fCurvedObsLevel)) {
        const double theta = sqrt(x*x + y*y) / Corsika::kEarthRadius;
        const double phi = atan2(y, x);
        const double r = Corsika::kEarthRadius + fObsLevelHeight;
        const double d = r * sin(theta);
        x = d * cos(phi);
        y = d * sin(phi);
        z = r * (cos(theta) - 1);
      } else if (fVersion >= 6.990 && fCurvedObsLevel) {
        const double r = Corsika::kEarthRadius + fObsLevelHeight;
        const double theta = sqrt(x*x + y*y) / r;
        const double phi = atan2(y, x);
        const double d = r * sin(theta);
        x = d * cos(phi);
        y = d * sin(phi);
        z = r * (cos(theta) - 1);
      }

      fCurrentParticle =
        Particle(particleId,
                 source,
                 Point(x,
                       y,
                       z,
                       groundCS),
                 Vector(corsikaParticle->fPx*GeV,
                        corsikaParticle->fPy*GeV,
                        -corsikaParticle->fPz*GeV,
                        groundCS),
                 TimeInterval(corsikaParticle->fTorZ*ns - fTimeOffset),
                 weight);

      if (fProductionPoint.IsValid())
        fCurrentParticle.SetProductionPoint(fProductionPoint.Get());
      if (fParent.IsValid())
        fCurrentParticle.SetParent(fParent.Get());
      if (fGrandParent.IsValid())
        fCurrentParticle.GetParent().SetParent(fGrandParent.Get());

      fProductionPoint.SetValid(false);
      fGrandParent.SetValid(false);
      fParent.SetValid(false);

      return &fCurrentParticle;

    } else { // NOT THINNED

      const Corsika::ParticleDataUnthinned* const corsikaParticle = GetOneParticleRecordUnthinned();

      // end of particle list
      if (!corsikaParticle)
        return 0;

      if (fCERFIL) {

        const Corsika::CherenkovDataUnthinned* const cherenkovPhoton = (const Corsika::CherenkovDataUnthinned*)corsikaParticle;

        const double weight = cherenkovPhoton->fPhotonsInBunch;
        if (weight == 0 || (cherenkovPhoton->fX == 0 && cherenkovPhoton->fY == 0)) // empty block is all zeros
          return 0;

        const double x = cherenkovPhoton->fX*cm;
        const double y = cherenkovPhoton->fY*cm;
        const double z = 0;
        /*
          Note: Cherenkov photon output is NEVER on curved observation level !!!!
                RU, Do 15. Dez 15:49:12 CET 2016
        */
        const double sintheta_cosphi = cherenkovPhoton->fU; // corsika x-direction
        const double sintheta_sinphi = cherenkovPhoton->fV; // corsika y-direction
        const double costheta = -sqrt(1 - Sqr(sintheta_sinphi) - Sqr(sintheta_cosphi)); // downwards

        const Point posOnGround(x, y, z, groundCS);
        const Vector direction(sintheta_cosphi, sintheta_sinphi, costheta, groundCS);

        // for MCERFI<=2, fDistance is PRODUCTION HEIGHT, for MCERFI>2 it is DISTANCE
        const double distance = cherenkovPhoton->fDistance * (fMCERFI > 2 ? 1. : 1/costheta) * cm;
        const Point posAtProduction = posOnGround - distance*direction;

        utl::RandomEngine* const randomEngine = &fwk::RandomEngineRegistry::GetInstance().Get(fwk::RandomEngineRegistry::eDetector);

        // if (since) wavelength is not stored, sample one:
        const double lambda = 1 / (1/fWlLowerLimit - RandFlat::shoot(&randomEngine->GetEngine(), 0, 1) * (1/fWlLowerLimit - 1/fWlUpperLimit));
        const double energy = utl::kPlanck / (lambda );

        fCurrentParticle =
          Particle(utl::Particle::ePhoton,
                   Particle::eShower,
                   posAtProduction, direction,
                   TimeInterval(cherenkovPhoton->fT*ns - fTimeOffset - distance/utl::kSpeedOfLight),
                   weight, energy);

        fProductionPoint.SetValid(false);
        fGrandParent.SetValid(false);
        fParent.SetValid(false);

        return &fCurrentParticle;
      }

      const int info = corsikaParticle->fDescription / 1000;
      const int particleId = Corsika::CorsikaToPDG(abs(info));

      // muon additional info (75,76: regular muons. 85,86: decaying muons at production point)
      if (info == 75 || info == 76 || info == 85 || info == 86) {
        if (!fCurved) {
          fProductionPoint = Point(corsikaParticle->fX*cm,
                                   corsikaParticle->fY*cm,
                                   corsikaParticle->fTorZ*cm - fObsLevelHeight,
                                   groundCS);
          continue;
        }
      }

      // negative particle description denotes mother and grandmother particles
      if (corsikaParticle->fDescription < 0) {

        if (particleId == Particle::eUndefined) {
          ostringstream what;
          what << "Unknown corsika particle code " << abs(int(corsikaParticle->fDescription/1000));
          WARNING(what);
        }

        if (!fParent.IsValid()) {
          if (!fCurved) {
            fParent = Particle(particleId,
                               Particle::eShower,
                               Point(corsikaParticle->fX*cm,
                                     corsikaParticle->fY*cm,
                                     corsikaParticle->fTorZ*cm - fObsLevelHeight,
                                     groundCS),
                               Vector(corsikaParticle->fPx*GeV,
                                      corsikaParticle->fPy*GeV,
                                      -corsikaParticle->fPz*GeV,
                                      groundCS),
                               TimeInterval(),
                               1);
          }
        } else if (!fGrandParent.IsValid()) {
          if (!fCurved) {
            const Vector direction(corsikaParticle->fPx*GeV,
                                   corsikaParticle->fPy*GeV,
                                   -corsikaParticle->fPz*GeV,
                                   groundCS);
            const double cosTheta = direction.GetCosTheta(groundCS);
            const double parent_h = corsikaParticle->fTorZ*cm;
            // Subtract traversed depth along the momentum direction
            const double grand_parent_depth = GetDepthAtHeight(parent_h) - corsikaParticle->fY*g/cm2 * cosTheta;
            const double d = (GetHeightAtDepth(grand_parent_depth) - parent_h) / cosTheta;
            const Point position = fParent.Get().GetPosition() - d*direction/direction.GetMag();
            fGrandParent = Particle(particleId,
                                    Particle::eShower,
                                    position,
                                    direction,
                                    TimeInterval(),
                                    1);
          }
        }
        continue;
      }

      // skip unknown particles...
      if (particleId == Particle::eUndefined) {
        fProductionPoint.SetValid(false);
        fGrandParent.SetValid(false);
        fParent.SetValid(false);
        continue;
      }

      // these are not defined in ehistory particles
      const short int hadGen = fmod(corsikaParticle->fDescription, 1000) / 10;
      short unsigned int obsLevel = fmod(corsikaParticle->fDescription, 10);
      //short unsigned int fateIndex = 0;

      // 95,96: decaying muons at point of decay. observation level is replaced by fate index
      if (info == 95 || info == 96) {
        obsLevel = fObservationLevel;
        //fateIndex = fmod(corsikaParticle->fDescription, 10);
      }

      if (obsLevel != fObservationLevel) {
        fProductionPoint.SetValid(false);
        fGrandParent.SetValid(false);
        fParent.SetValid(false);
        continue;
      }

      const float weight = 1;

      const Particle::Source source =
        (fVersion >= 6.71 &&
         (particleId == Particle::eElectron ||
          particleId == Particle::ePositron ||
          particleId == Particle::ePhoton) &&
         hadGen > 50) ?
        Particle::eShowerFromMuonDecay : Particle::eShower;

      double x = corsikaParticle->fX*cm;
      double y = corsikaParticle->fY*cm;
      double z = 0;
      if ((fVersion == 6.98) ||
          (fVersion >= 6.981 && fVersion < 6.990 && fCurvedObsLevel)) {
        const double theta = sqrt(x*x + y*y) / Corsika::kEarthRadius;
        const double phi = atan2(y, x);
        const double r = Corsika::kEarthRadius + fObsLevelHeight;
        const double d = r * sin(theta);
        x = d * cos(phi);
        y = d * sin(phi);
        z = r * (cos(theta) - 1);
      } else if (fVersion >= 6.990 && fCurvedObsLevel) {
        const double r = Corsika::kEarthRadius + fObsLevelHeight;
        const double theta = sqrt(x*x + y*y) / r;
        const double phi = atan2(y, x);
        const double d = r * sin(theta);
        x = d * cos(phi);
        y = d * sin(phi);
        z = r * (cos(theta) - 1);
      }

      fCurrentParticle =
        Particle(particleId,
                 source,
                 Point(x,
                       y,
                       z,
                       groundCS),
                 Vector(corsikaParticle->fPx*GeV,
                        corsikaParticle->fPy*GeV,
                        -corsikaParticle->fPz*GeV,
                        groundCS),
                 TimeInterval(corsikaParticle->fTorZ*ns - fTimeOffset),
                 weight);
      if (fProductionPoint.IsValid())
        fCurrentParticle.SetProductionPoint(fProductionPoint.Get());
      if (fParent.IsValid())
        fCurrentParticle.SetParent(fParent.Get());
      if (fGrandParent.IsValid())
        fCurrentParticle.GetParent().SetParent(fGrandParent.Get());

      fProductionPoint.SetValid(false);
      fGrandParent.SetValid(false);
      fParent.SetValid(false);

      return &fCurrentParticle;

    }
  }
}


const Corsika::ParticleData*
CorsikaShowerFileParticleIterator::GetOneParticleRecord()
{
  using std::ostringstream;

  if (!fIteratorValid)
    throw CorsikaIOException("CorsikaShowerFileParticleIterator not valid.");

  if (!fBlockBufferValid) {
    fRawFile.SeekTo(fCurrentPosition);
    if (!fRawFile.GetNextBlock(fCurrentBlock)) {
      ostringstream err;
      err << "Error reading block " << fCurrentPosition << " in CORSIKA file.";
      ERROR(err);
      throw CorsikaIOException(err.str());
    }

    if (fCurrentBlock.IsControl()) { // end of particle records
      // cout << " ------------------- CONTROL BLOCK ----------" << endl;
      fIteratorValid = false;
      return 0;
    }
  }

  const Corsika::ParticleData* const currentRecord =
    fCurrentBlock.AsParticleBlock().fParticle + fParticleInBlock;
  ++fParticleInBlock;

  if (fParticleInBlock >= io::Corsika::kParticlesInBlock) {
    ++fCurrentPosition;
    fParticleInBlock = 0;
    fBlockBufferValid = false;
  }

  return currentRecord;
}


const Corsika::ParticleDataUnthinned*
CorsikaShowerFileParticleIterator::GetOneParticleRecordUnthinned()
{
  using std::ostringstream;

  if (!fIteratorValid)
    throw CorsikaIOException("CorsikaShowerFileParticleIterator not valid.");

  if (!fBlockBufferValid) {
    fRawFile.SeekTo(fCurrentPosition);
    if (!fRawFile.GetNextBlockUnthinned(fCurrentBlockUnthinned)) {
      ostringstream err;
      err << "Error reading block " << fCurrentPosition << " in CORSIKA file.";
      ERROR(err);
      throw CorsikaIOException(err.str());
    }

    if (fCurrentBlockUnthinned.IsControl()) { // end of particle records
      fIteratorValid = false;
      return 0;
    }
  }

  const Corsika::ParticleDataUnthinned* const currentRecord =
    fCurrentBlockUnthinned.AsParticleBlock().fParticle + fParticleInBlock;
  ++fParticleInBlock;

  if (fParticleInBlock >= io::Corsika::kParticlesInBlock) {
    ++fCurrentPosition;
    fParticleInBlock = 0;
    fBlockBufferValid = false;
  }

  return currentRecord;
}


double
CorsikaShowerFileParticleIterator::GetHeightAtDepth(const double depth)
  const
{
  const double d = (depth < 0) ? 0 : depth;

  // depth values in ZLAY[] are in descending order
  for (int i = 0; i < 4; ++i)
    if (fAtmPars.fZLAY[i+1] < d && d <= fAtmPars.fZLAY[i])
      return -fAtmPars.fCATM[i] * log((d - fAtmPars.fAATM[i]) / fAtmPars.fBATM[i]);

  if (d <= fAtmPars.fZLAY[4])
    return fAtmPars.fCATM[4] * (fAtmPars.fAATM[4] - d) / fAtmPars.fBATM[4];

  return 0;
}


double
CorsikaShowerFileParticleIterator::GetDepthAtHeight(const double height)
  const
{
  // height values in HLAY[] are in ascending order

//#warning why use fabs(height)?
  const double h = fabs(height);

  // figure out in what layer of the atmosphere model the first interaction happened
  int section = 0;
  if (h > fAtmPars.fHLAY[4])
    section = 4;
  else if (h > fAtmPars.fHLAY[3])
    section = 3;
  else if (h > fAtmPars.fHLAY[2])
    section = 2;
  else if (h > 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(-h / fAtmPars.fCATM[section]);
  return fAtmPars.fAATM[section] -
    fAtmPars.fBATM[section] * (h / fAtmPars.fCATM[section]);
}