VectorStack.hpp

/*
 * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
 *
 * This software is distributed under the terms of the GNU General Public
 * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
 * the license.
 */

#pragma once

#include <corsika/framework/core/ParticleProperties.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/stack/Stack.hpp>

#include <corsika/framework/geometry/Point.hpp>
#include <corsika/framework/geometry/Vector.hpp>
#include <corsika/framework/geometry/PhysicalGeometry.hpp>

#include <string>
#include <tuple>
#include <vector>

namespace corsika {

  template <typename TStackIterator>
  class ParticleInterface : public ParticleBase<TStackIterator> {

  private:
    typedef ParticleBase<TStackIterator> super_type;

  public:
    typedef std::tuple<Code, HEPEnergyType, DirectionVector, Point, TimeType>
        particle_data_type;

    typedef std::tuple<Code, MomentumVector, Point, TimeType> particle_data_momentum_type;

    std::string asString() const;

    void setParticleData(particle_data_type const& v);

    void setParticleData(ParticleInterface<TStackIterator> const& p,
                         particle_data_type const& v);

    void setParticleData(particle_data_momentum_type const& v);

    void setParticleData(ParticleInterface<TStackIterator> const& p,
                         particle_data_momentum_type const& v);

    void setPID(Code const id) {
      super_type::getStackData().setPID(super_type::getIndex(), id);
    }

    void setEnergy(HEPEnergyType const& e) {
      super_type::getStackData().setKineticEnergy(super_type::getIndex(),
                                                  e - this->getMass());
    }

    void setKineticEnergy(HEPEnergyType const& ekin) {
      super_type::getStackData().setKineticEnergy(super_type::getIndex(), ekin);
    }

    void setMomentum(MomentumVector const& v) {
      HEPMomentumType const P = v.getNorm();
      if (P == 0_eV) {
        super_type::getStackData().setKineticEnergy(super_type::getIndex(), 0_eV);
        super_type::getStackData().setDirection(
            super_type::getIndex(), DirectionVector(v.getCoordinateSystem(), {0, 0, 0}));
      } else {
        super_type::getStackData().setKineticEnergy(
            super_type::getIndex(),
            sqrt(square(getMass()) + square(P)) - this->getMass());
        super_type::getStackData().setDirection(super_type::getIndex(), v / P);
      }
    }
    void setDirection(DirectionVector const& v) {
      super_type::getStackData().setDirection(super_type::getIndex(), v);
    }
    void setPosition(Point const& v) {
      super_type::getStackData().setPosition(super_type::getIndex(), v);
    }
    void setTime(TimeType const& v) {
      super_type::getStackData().setTime(super_type::getIndex(), v);
    }

    Code getPID() const {
      return super_type::getStackData().getPID(super_type::getIndex());
    }
    PDGCode getPDG() const { return get_PDG(getPID()); }
    HEPEnergyType getKineticEnergy() const {
      return super_type::getStackData().getKineticEnergy(super_type::getIndex());
    }
    DirectionVector getDirection() const {
      return super_type::getStackData().getDirection(super_type::getIndex());
    }
    Point getPosition() const {
      return super_type::getStackData().getPosition(super_type::getIndex());
    }
    TimeType getTime() const {
      return super_type::getStackData().getTime(super_type::getIndex());
    }
    VelocityVector getVelocity() const {
      return this->getMomentum() / this->getEnergy() * constants::c;
    }
    MomentumVector getMomentum() const {
      auto const P = sqrt(square(getEnergy()) - square(this->getMass()));
      return super_type::getStackData().getDirection(super_type::getIndex()) * P;
    }
    HEPMassType getMass() const { return get_mass(this->getPID()); }

    ElectricChargeType getCharge() const { return get_charge(this->getPID()); }

    HEPEnergyType getEnergy() const { return this->getKineticEnergy() + this->getMass(); }

    int16_t getChargeNumber() const { return get_charge_number(this->getPID()); }
  };

  class VectorStackImpl {

  public:
    typedef std::vector<Code> code_vector_type;
    typedef std::vector<HEPEnergyType> kinetic_energy_vector_type;
    typedef std::vector<Point> point_vector_type;
    typedef std::vector<TimeType> time_vector_type;
    typedef std::vector<DirectionVector> direction_vector_type;

    VectorStackImpl() = default;

    VectorStackImpl(VectorStackImpl const& other) = default;

    VectorStackImpl(VectorStackImpl&& other) = default;

    VectorStackImpl& operator=(VectorStackImpl const& other) = default;

    VectorStackImpl& operator=(VectorStackImpl&& other) = default;

    void dump() const {}

    void clear();

    unsigned int getSize() const { return dataPID_.size(); }
    unsigned int getCapacity() const { return dataPID_.size(); }

    void setPID(size_t i, Code const id) { dataPID_[i] = id; }
    void setKineticEnergy(size_t i, HEPEnergyType const& e) { dataEkin_[i] = e; }
    void setDirection(size_t i, DirectionVector const& v) { direction_[i] = v; }
    void setPosition(size_t i, Point const& v) { position_[i] = v; }
    void setTime(size_t i, TimeType const& v) { time_[i] = v; }

    Code getPID(size_t i) const { return dataPID_[i]; }
    HEPEnergyType getKineticEnergy(size_t i) const { return dataEkin_[i]; }
    DirectionVector getDirection(size_t i) const { return direction_[i]; }
    Point getPosition(size_t i) const { return position_[i]; }
    TimeType getTime(size_t i) const { return time_[i]; }

    void copy(size_t i1, size_t i2);

    void swap(size_t i1, size_t i2);

    void incrementSize();
    void decrementSize();

  private:
    code_vector_type dataPID_;
    kinetic_energy_vector_type dataEkin_;
    direction_vector_type direction_;
    point_vector_type position_;
    time_vector_type time_;

  }; // end class VectorStackImpl

  typedef Stack<VectorStackImpl, ParticleInterface> VectorStack;

} // namespace corsika

#include <corsika/detail/stack/VectorStack.inl>