Line data    Source code

       1             : #ifndef CRPROPA_PARTICLEMAPSCONTAINER_HH
       2             : #define CRPROPA_PARTICLEMAPSCONTAINER_HH
       3             : 
       4             : #include <map>
       5             : #include <vector>
       6             : #include "crpropa/magneticLens/Pixelization.h"
       7             : #include "crpropa/magneticLens/MagneticLens.h"
       8             : 
       9             : #include "crpropa/Vector3.h"
      10             : 
      11             : namespace crpropa {
      12             : /**
      13             :  * \addtogroup MagneticLenses 
      14             :  * @{
      15             :  */
      16             : 
      17             : /** 
      18             :  @class ParticleMapsContainer
      19             :  @brief A container for particle maps
      20             : 
      21             :  The maps are stored with discrete energies on a logarithmic scale. The
      22             :  default energy width is 0.02 with an energy bin from 10**17.99 - 10**18.01 eV.
      23             :  */
      24             : class ParticleMapsContainer {
      25             : private:
      26             :         std::map<int, std::map <int, double*> > _data;
      27             :         Pixelization _pixelization;
      28             :         double _deltaLogE;
      29             :         double _bin0lowerEdge;
      30             : 
      31             :         // get the bin number of the energy
      32             :         int energy2Idx(double energy) const;
      33             :         double idx2Energy(int idx) const;
      34             : 
      35             :         // weights of the particles
      36             :         double _sumOfWeights;
      37             :         std::map<int, double > _weightsPID;
      38             :         std::map<int, map<int, double> > _weights_pidEnergy;
      39             : 
      40             :         // lazy update of weights
      41             :         bool _weightsUpToDate;
      42             :         void _updateWeights();
      43             : 
      44             : public:
      45             :         /** Constructor.
      46             :          @param deltaLogE               width of logarithmic energy bin [in eV]
      47             :          @param bin0lowerEdge   logarithm of energy of the lower edge of first bin [in log(eV)]
      48             :          */
      49           2 :         ParticleMapsContainer(double deltaLogE = 0.02, double bin0lowerEdge = 17.99) : _deltaLogE(deltaLogE), _bin0lowerEdge(bin0lowerEdge), _pixelization(6), _weightsUpToDate(false), _sumOfWeights(0) {
      50           2 :         }
      51             :         /** Destructor.
      52             :          */
      53             :         ~ParticleMapsContainer();
      54             : 
      55             :         size_t getNumberOfPixels() {
      56           0 :                 return _pixelization.getNumberOfPixels();
      57             :         }
      58             : 
      59             :         /** Get the map for the particleId with the given energy.
      60             :          @param particleId              id of the particle following the PDG numbering scheme
      61             :          @param energy                  the energy of the particle [in Joules]
      62             :          @returns The map for a given particleId with a given energy
      63             :          */
      64             :         double *getMap(const int particleId, double energy);
      65             : 
      66             :         /** Adds a particle to the map container.
      67             :          @param particleId                      id of the particle following the PDG numbering scheme
      68             :          @param energy                          the energy of the particle [in Joules]
      69             :          @param galacticLongitude       galactic longitude [radians]
      70             :          @param galacticLatitude        galactic latitude [radians]
      71             :          @param weight                          relative weight for the specific particle
      72             :         */
      73             :         void addParticle(const int particleId, double energy, double galacticLongitude, double galacticLatitude, double weight = 1);
      74             :         /** Adds a particle to the map container.
      75             :          @param particleId                      id of the particle following the PDG numbering scheme
      76             :          @param energy                          the energy of the particle [in Joules]
      77             :          @param v                                       vector containing the arrival directions of a particle
      78             :          @param weight                          relative weight for the specific particle
      79             :         */
      80             :         void addParticle(const int particleId, double energy, const Vector3d &v, double weight = 1);
      81             : 
      82             :         /** Get all particle ids in the map.
      83             :          @returns Vector of all ids.
      84             :          */
      85             :         std::vector<int> getParticleIds();
      86             : 
      87             :         /** Get energies in map.
      88             :          @param pid     id of the particle following the PDG numbering scheme
      89             :          @returns Energies are returned in units of eV (unlike in other CRPropa modules) for performance reasons
      90             :          */
      91             :         std::vector<double> getEnergies(int pid);
      92             : 
      93             :         void applyLens(MagneticLens &lens);;
      94             : 
      95             :         /** Get random particles from map.
      96             :          The arguments are the vectors where the information will be stored.
      97             :          @param N                                       number of particles to be selected
      98             :          @param particleId                      id of the particle following the PDG numbering scheme
      99             :          @param energy                          energy of interest [in eV]
     100             :          @param galacticLongitudes      longitude in the interval [-pi, pi] [in radians]
     101             :          @param galacticLatitudes       latitude in the interval [-pi/2, pi/2] [in radians]
     102             :          */
     103             :         void getRandomParticles(size_t N, vector<int> &particleId,
     104             :                 vector<double> &energy, vector<double> &galacticLongitudes,
     105             :                 vector<double> &galacticLatitudes);
     106             : 
     107             :         /** Places a particle with given id and energy according to the  probability maps. 
     108             :          @param pid                                     id of the particle following the PDG numbering scheme
     109             :          @param energy                          energy of interest [in eV]
     110             :          @param galacticLongitude       longitude in the interval [-pi, pi] [in radians]
     111             :          @param galacticLatitude        latitude in the interval [-pi/2, pi/2] [in radians]
     112             :          @returns Returns false if operation not possible; true otherwise.
     113             :          */
     114             :         bool placeOnMap(int pid, double energy, double &galacticLongitude, double &galacticLatitude);
     115             : 
     116             :         /** Force weight update prior to getting random particles. 
     117             :          Only necessary when reusing pointer to maps after calculating weights.
     118             :         */
     119             :         void forceWeightUpdate();
     120             : 
     121             :         /** Get sum of weights of the maps.
     122             :          @returns Sum of all weights.
     123             :          */
     124             :         double getSumOfWeights() {
     125           0 :                 if (!_weightsUpToDate)
     126           0 :                         _updateWeights();
     127           0 :                 return _sumOfWeights;
     128             :         }
     129             : 
     130             :         /** Get weight for a given particle and energy
     131             :          @param pid                                     id of the particle following the PDG numbering scheme
     132             :          @param energy                          energy of interest [in eV]
     133             :          @returns Weight for the chosen particle and energy.
     134             :          */
     135           0 :         double getWeight(int pid, double energy) {
     136           0 :                 if (!_weightsUpToDate)
     137           0 :                         _updateWeights();
     138           0 :                 return _weights_pidEnergy[pid][energy2Idx(energy)];
     139             :         }
     140             : };
     141             : /** @}*/
     142             : 
     143             : } // namespace crpropa
     144             : 
     145             : #endif // CRPROPA_PARTICLEMAPSCONTAINER_HH