ShowerSimData.h

Go to the documentation of this file.

#ifndef _evt_ShowerSimData_h_
#define _evt_ShowerSimData_h_

#include <utl/CoordinateSystemPtr.h>
#include <utl/ShowerParticleIterator.h>
#include <utl/ShowerParticleList.h>
#include <utl/VParticleProperties.h>
#include <utl/MultiTabulatedFunction.h>
#include <utl/Vector.h>
#include <utl/Point.h>
#include <utl/TimeStamp.h>
#include <utl/AugerUnits.h>
#include <utl/LameShadowPtr.h>
#include <utl/config.h>
#include <utl/ShadowPtr.h>
#include <utl/Branch.h>
#include <utl/IteratorRange.h>

#include <evt/VShowerGeometryProducer.h>
#include <evt/VGaisserHillasParameter.h>
#include <evt/AtmosphereParameters.h>
#include <evt/GenParticle.h>

#include <string>
#include <cmath>


namespace io {
  class ShowerSimData_ROOT;
}

namespace evt {

  class LaserData;
  class RadioSimulation;
  class Event;


  class ShowerSimData  {

  public:
    enum ProfileType {
      eCharged = 0,
      eEnergyDeposit,
      ePhoton,
      eMuon,
      eElectron,
      eHadron,
      eMuonProduction,
      eCherenkov
    };

    int GetShowerNumber() const { return fShowerNumber; }
    void SetShowerNumber(const int sid) { fShowerNumber = sid; }

    std::string GetShowerRunId() const { return fShowerRunId; }
    void SetShowerRunId(const std::string srid) { fShowerRunId = srid; }

    int GetPrimaryParticle() const { return fPrimaryProperties->GetType(); }
    void SetPrimaryParticle(const int type);

    double GetEnergy() const { return fEnergy; }
    void SetEnergy(const double theEnergy) { fEnergy = theEnergy; }

    double GetMuonNumber() const { return fMuonNumber; }
    void SetMuonNumber(const double nmuon) { fMuonNumber = nmuon; }

    double GetMuonWeightScale() const { return fMuonWeightScale; }
    void SetMuonWeightScale(const double scale) { fMuonWeightScale = scale; }

    double GetCalorimetricEnergy() const { return fCalorimetricEnergy; }
    void SetCalorimetricEnergy(const double energy)
    { fCalorimetricEnergy = energy; }

    double GetElectromagneticEnergy() const { return fElectromagneticEnergy; }
    void SetElectromagneticEnergy(const double energy)
    { fElectromagneticEnergy = energy; }

    double GetEnergyCutoff(const ProfileType type = eElectron) const;
    void SetEnergyCutoff(const double energy, const ProfileType type = eElectron);

    double GetMinCherenkovWavelength() const { return fMinCherenkovWavelength; }
    void SetMinCherenkovWavelength(const double wl) { fMinCherenkovWavelength = wl; }

    double GetMaxCherenkovWavelength() const { return fMaxCherenkovWavelength; }
    void SetMaxCherenkovWavelength(const double wl) { fMaxCherenkovWavelength = wl; }

    void SetGroundParticleCoordinateSystemAzimuth(const double azimuth);

    double GetGroundParticleCoordinateSystemAzimuth() const { return fGroundParticleCoordinateSystemAzimuth; }

    void SetGroundParticleCoordinateSystemZenith(const double zenith);

    double GetGroundParticleCoordinateSystemZenith() const { return fGroundParticleCoordinateSystemZenith; }

    bool HasGeometry() const { return bool(fParticleCS); }
    void MakeGeometry(const utl::Point& pointOnShowerAxis);
    void CheckGeometry() const;


    double GetMinRadiusCut() const { return fMinRadiusCut; }
    void SetMinRadiusCut(const double minR) { fMinRadiusCut = minR; }


    double GetMaxRadiusCut() const { return fMaxRadiusCut; }
    // Set the maximum radius cut
    void SetMaxRadiusCut(const double maxR) { fMaxRadiusCut = maxR; }


    double GetMagneticFieldZenith() const { return 90*utl::deg - std::fabs(fMagneticFieldInclination); }

    double GetMagneticFieldAzimuth() const { return fMagneticFieldDeclination + 90*utl::degree; }
    void SetMagneticFieldAzimuth(const double magneticFieldAzimuth) { fMagneticFieldDeclination = magneticFieldAzimuth - 90*utl::degree; }

    double GetMagneticFieldInclination() const { return fMagneticFieldInclination; }
    void SetMagneticFieldInclination(const double magneticFieldInclination) { fMagneticFieldInclination = magneticFieldInclination; }

    double GetMagneticFieldDeclination() const { return fMagneticFieldDeclination; }
    void SetMagneticFieldDeclination(const double magneticFieldDeclination) { fMagneticFieldDeclination = magneticFieldDeclination; }

    double GetMagneticFieldStrength() const { return fMagneticFieldStrength; }
    void SetMagneticFieldStrength(const double magneticFieldStrength) { fMagneticFieldStrength = magneticFieldStrength; }

    void AddSimCore(const utl::Point& pos);
    int GetNSimCores() const { return fSimCores.size(); }
    const utl::Point& GetSimCore(const int i) const { return fSimCores[i]; }

    const utl::Vector& GetDirection() const;
    bool HasDirection() const { return HasGeometry(); }

    const utl::Point& GetPosition() const;
    bool HasPosition() const { return HasGeometry(); }

    utl::CoordinateSystemPtr GetLocalCoordinateSystem() const { CheckGeometry(); return fLocalCS; }

    utl::CoordinateSystemPtr GetShowerCoordinateSystem() const { CheckGeometry(); return fShowerCS; }

    utl::CoordinateSystemPtr GetGroundParticleCoordinateSystem() const { CheckGeometry(); return fParticleCS; }

    const utl::TabulatedFunction& GetLongitudinalProfile(const ProfileType type = eCharged) const
    { return fLongitudinalProfile.GetTabulatedFunction(type); }
    utl::TabulatedFunction& GetLongitudinalProfile(const ProfileType type = eCharged)
    { return fLongitudinalProfile.GetTabulatedFunction(type); }

    bool HasLongitudinalProfile(const ProfileType type = eCharged) const
    { return fLongitudinalProfile.HasLabel(type); }

    void MakeLongitudinalProfile(const utl::TabulatedFunction& lp, const ProfileType type = eCharged)
    { fLongitudinalProfile.AddTabulatedFunction(lp, type); }
    void MakeLongitudinalProfile(const ProfileType type = eCharged)
    { fLongitudinalProfile.AddTabulatedFunction(type); }

    const utl::TabulatedFunction& GetdEdX() const { return GetLongitudinalProfile(eEnergyDeposit); }
    bool HasdEdX() const { return HasLongitudinalProfile(eEnergyDeposit); }
    void MakedEdX(const utl::TabulatedFunction& dEdX) { MakeLongitudinalProfile(dEdX, eEnergyDeposit); }

    const utl::TabulatedFunction& GetFluorescencePhotons(const int wavelength) const
    { return fFluorescencePhotons.GetTabulatedFunction(wavelength); }
    utl::TabulatedFunction& GetFluorescencePhotons(const int wavelength)
    { return fFluorescencePhotons.GetTabulatedFunction(wavelength); }
    bool HasFluorescencePhotons(const int wavelength) const
    { return fFluorescencePhotons.HasLabel(wavelength); }
    void AddFluorescencePhotons(const utl::TabulatedFunction& fp, const int wavelength)
    { fFluorescencePhotons.AddTabulatedFunction(fp, wavelength); }

    const utl::TabulatedFunction& GetCherenkovPhotons(const int wavelength) const
    { return fCherenkovPhotons.GetTabulatedFunction(wavelength); }
    utl::TabulatedFunction& GetCherenkovPhotons(const int wavelength)
    { return fCherenkovPhotons.GetTabulatedFunction(wavelength); }
    bool HasCherenkovPhotons(const int wavelength) const
    { return fCherenkovPhotons.HasLabel(wavelength); }
    void AddCherenkovPhotons(const utl::TabulatedFunction& cp, const int wavelength)
    { fCherenkovPhotons.AddTabulatedFunction(cp, wavelength); }

    const utl::TabulatedFunction& GetCherenkovBeamPhotons(const int wavelength) const
    { return fCherenkovBeamPhotons.GetTabulatedFunction(wavelength); }
    utl::TabulatedFunction& GetCherenkovBeamPhotons(const int wavelength)
    { return fCherenkovBeamPhotons.GetTabulatedFunction(wavelength); }
    bool HasCherenkovBeamPhotons(const int wavelength) const
    { return fCherenkovBeamPhotons.HasLabel(wavelength); }
    void AddCherenkovBeamPhotons(const utl::TabulatedFunction& cp, const int wavelength)
    { fCherenkovBeamPhotons.AddTabulatedFunction(cp, wavelength); }

    const utl::TabulatedFunction& GetCherenkovBeamProductionPhotons(const int wavelength) const
    { return fCherenkovBeamProductionPhotons.GetTabulatedFunction(wavelength); }
    utl::TabulatedFunction& GetCherenkovBeamProductionPhotons(const int wavelength)
    { return fCherenkovBeamProductionPhotons.GetTabulatedFunction(wavelength); }
    bool HasCherenkovBeamProductionPhotons(const int wavelength) const
    { return fCherenkovBeamProductionPhotons.HasLabel(wavelength); }
    void AddCherenkovBeamProductionPhotons(const utl::TabulatedFunction& fp, const int wavelength)
    { fCherenkovBeamProductionPhotons.AddTabulatedFunction(fp, wavelength); }

#ifdef PRIM_CHER
    const utl::TabulatedFunction& GetPrimaryCherenkovPhotons(const int wavelength) const
    { return fPrimaryCherenkovPhotons.GetTabulatedFunction(wavelength); }
    utl::TabulatedFunction& GetPrimaryCherenkovPhotons(const int wavelength)
    { return fPrimaryCherenkovPhotons.GetTabulatedFunction(wavelength); }
    bool HasPrimaryCherenkovPhotons(const int wavelength) const
    { return fPrimaryCherenkovPhotons.HasLabel(wavelength); }
    void AddPrimaryCherenkovPhotons(const utl::TabulatedFunction& cp, const int wavelength)
    { fPrimaryCherenkovPhotons.AddTabulatedFunction(cp, wavelength); }

    const utl::TabulatedFunction& GetPrimaryCherenkovBeamPhotons(const int wavelength) const
    { return fPrimaryCherenkovBeamPhotons.GetTabulatedFunction(wavelength); }
    utl::TabulatedFunction& GetPrimaryCherenkovBeamPhotons(const int wavelength)
    { return fPrimaryCherenkovBeamPhotons.GetTabulatedFunction(wavelength); }
    bool HasPrimaryCherenkovBeamPhotons(const int wavelength) const
    { return fPrimaryCherenkovBeamPhotons.HasLabel(wavelength); }
    void AddPrimaryCherenkovBeamPhotons(const utl::TabulatedFunction& cp, const int wavelength)
    { fPrimaryCherenkovBeamPhotons.AddTabulatedFunction(cp, wavelength); }
#endif

    typedef utl::MultiTabulatedFunction::Iterator PhotonIterator;
    typedef utl::MultiTabulatedFunction::ConstIterator ConstPhotonIterator;

    PhotonIterator FirstFluorescenceWavelength()
    { return fFluorescencePhotons.TabulatedFunctionsBegin(); }
    ConstPhotonIterator FirstFluorescenceWavelength() const
    { return fFluorescencePhotons.TabulatedFunctionsBegin(); }
    PhotonIterator LastFluorescenceWavelength()
    { return fFluorescencePhotons.TabulatedFunctionsEnd(); }
    ConstPhotonIterator LastFluorescenceWavelength() const
    { return fFluorescencePhotons.TabulatedFunctionsEnd(); }

    PhotonIterator FirstCherenkovWavelength()
    { return fCherenkovPhotons.TabulatedFunctionsBegin(); }
    ConstPhotonIterator FirstCherenkovWavelength() const
    { return fCherenkovPhotons.TabulatedFunctionsBegin(); }
    PhotonIterator LastCherenkovWavelength()
    { return fCherenkovPhotons.TabulatedFunctionsEnd(); }
    ConstPhotonIterator LastCherenkovWavelength() const
    { return fCherenkovPhotons.TabulatedFunctionsEnd(); }

    PhotonIterator FirstCherenkovBeamWavelength()
    { return fCherenkovBeamPhotons.TabulatedFunctionsBegin(); }
    ConstPhotonIterator FirstCherenkovBeamWavelength() const
    { return fCherenkovBeamPhotons.TabulatedFunctionsBegin(); }
    PhotonIterator LastCherenkovBeamWavelength()
    { return fCherenkovBeamPhotons.TabulatedFunctionsEnd(); }
    ConstPhotonIterator LastCherenkovBeamWavelength() const
    { return fCherenkovBeamPhotons.TabulatedFunctionsEnd(); }

    PhotonIterator FirstCherenkovBeamProductionWavelength()
    { return fCherenkovBeamProductionPhotons.TabulatedFunctionsBegin(); }
    ConstPhotonIterator FirstCherenkovBeamProductionWavelength() const
    { return fCherenkovBeamProductionPhotons.TabulatedFunctionsBegin(); }
    PhotonIterator LastCherenkovBeamProductionWavelength()
    { return fCherenkovBeamProductionPhotons.TabulatedFunctionsEnd(); }
    ConstPhotonIterator LastCherenkovBeamProductionWavelength() const
    { return fCherenkovBeamProductionPhotons.TabulatedFunctionsEnd(); }

    const evt::VGaisserHillasParameter& GetGHParameters() const { return *fGHParameters; }
    bool HasGHParameters() const { return bool(fGHParameters); }
    void MakeGHParameters(const evt::VGaisserHillasParameter& ghPar);

    const evt::AtmosphereParameters& GetAtmosphereParameters() const { return *fAtmosphereParameters; }
    bool HasAtmosphereParameters() const { return bool(fAtmosphereParameters); }
    void MakeAtmosphereParameters(const evt::AtmosphereParameters& par);
    double GetXFirst() const { return fXFirst; }
    void SetXFirst(const double xFirst) { fXFirst = xFirst; }
    double GetXInject() const { return fXInject; }
    void SetXInject(const double xInject) { fXInject = xInject; }

    double GetNmu() const { return fNmu; }
    void SetNmu(const double Nmu) { fNmu = Nmu; }

    double GetXmaxMu() const { return fXmaxMu; }
    void SetXmaxMu(const double XmaxMu) { fXmaxMu = XmaxMu; }

    const utl::TimeStamp& GetTimeStamp() const { return *fTimeStamp; }
    bool HasTimeStamp() const { return bool(fTimeStamp); }
    void MakeTimeStamp(const utl::TimeStamp& ts);

    utl::ShowerParticleList& GetGroundParticles()
    { return *fGroundParticles; }
    const utl::ShowerParticleList& GetGroundParticles() const
    { return *fGroundParticles; }
    bool HasGroundParticles() const { return bool(fGroundParticles); }
    void MakeGroundParticles();
    utl::ShowerParticleIterator GroundParticlesBegin() const
    { return fGroundParticles->ParticlesBegin(GetGroundParticleCoordinateSystem()); }
    utl::ShowerParticleIterator GroundParticlesEnd() const
    { return fGroundParticles->ParticlesEnd(GetGroundParticleCoordinateSystem()); }
    OFFLINE_MAKE_CONST_ITERATOR_RANGE(utl::ShowerParticleIterator, GroundParticles)

    
    std::map<utl::AttributeMap, utl::ShadowPtr<utl::ShowerParticleList>>::const_iterator GetGroundCherenkovDataBegin() const { return fGroundCherenkov.begin(); }
    std::map<utl::AttributeMap, utl::ShadowPtr<utl::ShowerParticleList>>::const_iterator GetGroundCherenkovDataEnd() const { return fGroundCherenkov.end(); }
    utl::ShowerParticleList& GetGroundCherenkov(const utl::AttributeMap& am) { return *fGroundCherenkov[am]; }
    const utl::ShowerParticleList& GetGroundCherenkov(const utl::AttributeMap& am) const { return *fGroundCherenkov.find(am)->second; }
    bool HasGroundCherenkov(const utl::AttributeMap& am) const { return bool(fGroundCherenkov.count(am)) && bool(fGroundCherenkov.find(am)->second); }
    void MakeGroundCherenkov(const utl::AttributeMap& am);
    utl::ShowerParticleIterator GroundCherenkovBegin(const utl::AttributeMap& am) const
    { return fGroundCherenkov.find(am)->second->ParticlesBegin(GetGroundParticleCoordinateSystem()); }
    utl::ShowerParticleIterator GroundCherenkovEnd(const utl::AttributeMap& am) const
    { return fGroundCherenkov.find(am)->second->ParticlesEnd(GetGroundParticleCoordinateSystem()); }

    LaserData& GetLaserData() { return *fLaserData; }
    const LaserData& GetLaserData() const { return *fLaserData; }

    bool HasLaserData() const { return bool(fLaserData); }
    void MakeLaserData();

    RadioSimulation& GetRadioSimulation() { return *fRadioSimulation; }
    const RadioSimulation& GetRadioSimulation() const { return *fRadioSimulation; }

    bool HasRadioSimulation() const { return bool(fRadioSimulation); }

    void MakeRadioSimulation();

    double GetDistanceOfShowerMaximum() const { return fDistanceOfShowerMaximum; }

    void SetDistanceOfShowerMaximum(const double parDistance) { fDistanceOfShowerMaximum = parDistance; }

    const evt::GenParticle& GetParticleTree() const { return fParticleTree; }

    void SetParticleTree(const evt::GenParticle& partTree)
    { fParticleTree = partTree; }

    bool HasParticleTree() const { return bool(fParticleTree.GetParentId()); }

  private:
    ShowerSimData(const evt::VShowerGeometryProducer& p)
      : fShowerGeometryProducer(p.Clone()) { }

    ~ShowerSimData() = default;

    utl::ShadowPtr<evt::VShowerGeometryProducer> fShowerGeometryProducer;  // no root streaming
    utl::ShadowPtr<utl::ShowerParticleList> fGroundParticles;  // no root streaming
    std::map<utl::AttributeMap, utl::ShadowPtr<utl::ShowerParticleList>> fGroundCherenkov;  // no root streaming

    utl::ParticlePropertiesPtr fPrimaryProperties;
    double fEnergy = 0;
    double fMuonNumber = 0;
    double fMuonWeightScale = 0;
    double fCalorimetricEnergy = 0;
    double fElectromagneticEnergy = 0;
    double fGroundParticleCoordinateSystemAzimuth = 0;
    double fGroundParticleCoordinateSystemZenith = 0;
    double fEnergyCutoff = 0;
    double fEnergyCutoffMuon = 0;
    double fMinCherenkovWavelength = 0;
    double fMaxCherenkovWavelength = 0;
    double fMinRadiusCut = 0;
    double fMaxRadiusCut = 0;
    double fMagneticFieldInclination = 0;
    double fMagneticFieldDeclination = 0;
    double fMagneticFieldStrength = 0;
    int fShowerNumber = 0;
    std::string fShowerRunId;
    double fXFirst = 0;
    double fXInject = 0;
    double fDistanceOfShowerMaximum = std::numeric_limits<double>::quiet_NaN();

    std::vector<utl::Point> fSimCores;

    double fNmu = 0;
    double fXmaxMu = 0;

    utl::CoordinateSystemPtr fLocalCS;
    utl::CoordinateSystemPtr fParticleCS;
    utl::CoordinateSystemPtr fShowerCS;

    utl::ShadowPtr<utl::TimeStamp> fTimeStamp;
    utl::MultiTabulatedFunction fLongitudinalProfile;
    utl::MultiTabulatedFunction fFluorescencePhotons;
    utl::MultiTabulatedFunction fCherenkovPhotons;
    utl::MultiTabulatedFunction fCherenkovBeamPhotons;
    utl::MultiTabulatedFunction fCherenkovBeamProductionPhotons;
#ifdef PRIM_CHER
    utl::MultiTabulatedFunction fPrimaryCherenkovPhotons;
    utl::MultiTabulatedFunction fPrimaryCherenkovBeamPhotons;
#endif
    utl::ShadowPtr<evt::VGaisserHillasParameter> fGHParameters;
    utl::ShadowPtr<evt::AtmosphereParameters> fAtmosphereParameters;

    utl::LameShadowPtr<LaserData> fLaserData;

    utl::LameShadowPtr<RadioSimulation> fRadioSimulation;

    evt::GenParticle fParticleTree;

    mutable utl::ShadowPtr<utl::Point> fCachedPosition;
    mutable utl::ShadowPtr<utl::Vector> fCachedDirection;

    friend class evt::Event;
    friend class utl::LameShadowPtr<ShowerSimData>;
    friend class io::ShowerSimData_ROOT;

  };

}


#endif