Line data    Source code

       1             : #ifndef CRPROPA_EMISSION_MAP_H
       2             : #define CRPROPA_EMISSION_MAP_H
       3             : 
       4             : #include "Referenced.h"
       5             : #include "Candidate.h"
       6             : 
       7             : namespace crpropa {
       8             : 
       9             : /**
      10             :  @class CylindricalProjectionMap
      11             :  @brief 2D histogram of spherical coordinates in equal-area projection
      12             :  */
      13             : class CylindricalProjectionMap : public Referenced {
      14             : private:
      15             :         size_t nPhi, nTheta;
      16             :         double sPhi, sTheta;
      17             :         mutable bool dirty;
      18             :         std::vector<double> pdf;
      19             :         mutable std::vector<double> cdf;
      20             : 
      21             :         /** Calculate the cdf from the pdf */
      22             :         void updateCdf() const;
      23             : 
      24             : public:
      25             :         CylindricalProjectionMap();
      26             :         /** constructur
      27             :          * @param nPhi number of bins for phi (0-2pi)
      28             :          * @param nTheta number of bins for theta (0-pi)
      29             :          */
      30             :         CylindricalProjectionMap(size_t nPhi, size_t nTheta);
      31             : 
      32             :         /** Increment the bin value in direction by weight. */
      33             :         void fillBin(const Vector3d& direction, double weight = 1.);
      34             : 
      35             :         /** Increment the bin value by weight. */
      36             :         void fillBin(size_t bin, double weight = 1.);
      37             : 
      38             :         /** Draw a random vector from the distribution. */
      39             :         Vector3d drawDirection() const;
      40             : 
      41             :         /** Check if the direction has a non zero propabiliy. */
      42             :         bool checkDirection(const Vector3d &direction) const;
      43             : 
      44             :         const std::vector<double>& getPdf() const;
      45             :         std::vector<double>& getPdf();
      46             : 
      47             :         const std::vector<double>& getCdf() const;
      48             : 
      49             :         size_t getNPhi();
      50             :         size_t getNTheta();
      51             : 
      52             :         /** Calculate the bin from a direction */
      53             :         size_t binFromDirection(const Vector3d& direction) const;
      54             : 
      55             :         /** Calculate a random vector inside the bin boundaries */
      56             :         Vector3d directionFromBin(size_t bin) const;
      57             : };
      58             : 
      59             : /**
      60             :  @class EmissionMap
      61             :  @brief Particle Type and energy binned emission maps.
      62             : 
      63             :  Use SourceEmissionMap to suppress directions at the source. Use EmissionMapFiller to create EmissionMap from Observer.
      64             :  */
      65           2 : class EmissionMap : public Referenced {
      66             : public:
      67             :         typedef std::pair<int, size_t> key_t;
      68             :         typedef std::map<key_t, ref_ptr<CylindricalProjectionMap> > map_t;
      69             : 
      70             :         EmissionMap();
      71             :         /**
      72             :          * @param nPhi number of bins for phi (0-2pi)
      73             :          * @param nTheta number of bins for theta (0-pi)
      74             :          * @param nEnergy number of bins for energy (1e-4 - 1e4 EeV)
      75             :          */
      76             :         EmissionMap(size_t nPhi, size_t nTheta, size_t nEnergy);
      77             : 
      78             :         /**
      79             :          * @param nPhi number of bins for phi (0-2pi)
      80             :          * @param nTheta number of bins for theta (0-pi)
      81             :          * @param nEnergy number of bins for energy (1e-4 - 1e4 EeV)
      82             :          * @param minEnergy minimum energy for binning
      83             :          * @param maxEnergy maximum energy for binning
      84             :          */
      85             :         EmissionMap(size_t nPhi, size_t nTheta, size_t nEnergy, double minEnergy, double maxEnergy);
      86             : 
      87             :         /** Calculate energy from bin */
      88             :         double energyFromBin(size_t bin) const;
      89             : 
      90             :         /** Calculate bin from energy */
      91             :         size_t binFromEnergy(double energy) const;
      92             : 
      93             :         map_t &getMaps();
      94             :         const map_t &getMaps() const;
      95             : 
      96             :         /** Increment the value for particle type, energy and direction by weight. */
      97             :         void fillMap(int pid, double energy, const Vector3d& direction, double weight = 1.);
      98             :         /** Increment the value for the particle state by weight. */
      99             :         void fillMap(const ParticleState& state, double weight = 1.);
     100             : 
     101             :         /** Draw a random vector from the distribution. */
     102             :         bool drawDirection(int pid, double energy, Vector3d& direction) const;
     103             :         /** Draw a random vector from the distribution. */
     104             :         bool drawDirection(const ParticleState& state, Vector3d& direction) const;
     105             : 
     106             :         /** Check if the direction has a non zero propabiliy. */
     107             :         bool checkDirection(int pid, double energy, const Vector3d& direction) const;
     108             :         /** Check if the direction has a non zero propabiliy. */
     109             :         bool checkDirection(const ParticleState& state) const;
     110             : 
     111             :         /** Check if a valid map exists */
     112             :         bool hasMap(int pid, double energy);
     113             : 
     114             :         /** Get the map for the specified pid and energy */
     115             :         ref_ptr<CylindricalProjectionMap> getMap(int pid, double energy);
     116             : 
     117             :         /** Save the content of the maps into a text file */
     118             :         void save(const std::string &filename);
     119             :         /** Load the content of the maps from a text file */
     120             :         void load(const std::string &filename);
     121             : 
     122             :         /** Merge other maps, add pdfs */
     123             :         void merge(const EmissionMap *other);
     124             : 
     125             :         /** Merge maps from file */
     126             :         void merge(const std::string &filename);
     127             : 
     128             : protected:
     129             :         double minEnergy, maxEnergy, logStep;
     130             :         size_t nPhi, nTheta, nEnergy;
     131             :         map_t maps;
     132             : };
     133             : 
     134             : } // namespace crpropa
     135             : 
     136             : #endif // CRPROPA_EMISSION_MAP_H