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            1 : 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