Line data    Source code

       1             : #include "HepPID/ParticleIDMethods.hh"
       2             : #include "crpropa/Random.h"
       3             : #include "crpropa/magneticLens/ParticleMapsContainer.h"
       4             : #include "crpropa/Units.h"
       5             : 
       6             : #include <iostream>
       7             : #include <fstream>
       8             : 
       9             : namespace crpropa  {
      10             : 
      11           2 : ParticleMapsContainer::~ParticleMapsContainer() {
      12           2 :         for(std::map<int, std::map<int, double*> >::iterator pid_iter = _data.begin(); 
      13           4 :                         pid_iter != _data.end(); ++pid_iter) {
      14           2 :                 for(std::map<int, double*>::iterator energy_iter = pid_iter->second.begin();
      15           4 :                         energy_iter != pid_iter->second.end(); ++energy_iter) {
      16           2 :                         delete[] (energy_iter->second);
      17             :                 }
      18             :         }
      19           2 : }
      20             : 
      21         422 : int ParticleMapsContainer::energy2Idx(double energy) const {
      22         422 :         double lE = log10(energy / eV);
      23         422 :         return int((lE - _bin0lowerEdge) / _deltaLogE);
      24             : }
      25             : 
      26         421 : double ParticleMapsContainer::idx2Energy(int idx) const {
      27         421 :         return pow(10, idx * _deltaLogE + _bin0lowerEdge + _deltaLogE / 2) * eV;
      28             : }
      29             : 
      30             :                 
      31           0 : double* ParticleMapsContainer::getMap(const int particleId, double energy) {
      32           0 :         _weightsUpToDate = false;
      33           0 :         if (_data.find(particleId) == _data.end()) {
      34           0 :                 std::cerr << "No map for ParticleID " << particleId << std::endl;
      35           0 :                 return NULL;
      36             :         }
      37           0 :         int energyIdx   = energy2Idx(energy);
      38           0 :         if (_data[particleId].find(energyIdx) == _data[particleId].end()) {
      39           0 :                 std::cerr << "No map for ParticleID and energy" << energy / eV << " eV" << std::endl;
      40           0 :                 return NULL;
      41             :         }
      42           0 :         return _data[particleId][energy2Idx(energy)];
      43             : }
      44             :                         
      45             :                         
      46           2 : void ParticleMapsContainer::addParticle(const int particleId, double energy, double galacticLongitude, double galacticLatitude, double weight) {
      47           2 :         _weightsUpToDate = false;
      48           2 :         if (_data.find(particleId) == _data.end()) {
      49             :                 map<int, double*> M;
      50           2 :                 _data[particleId] = M;
      51             :         }
      52             : 
      53           2 :         int energyIdx   = energy2Idx(energy);
      54           4 :         if (_data[particleId].find(energyIdx) == _data[particleId].end()) {
      55           2 :                 _data[particleId][energyIdx] = new double[_pixelization.getNumberOfPixels()];
      56           2 :                 std::fill(_data[particleId][energyIdx], _data[particleId][energyIdx] + _pixelization.getNumberOfPixels(), 0);
      57             :         }
      58             : 
      59           2 :         uint32_t pixel = _pixelization.direction2Pix(galacticLongitude, galacticLatitude);
      60           2 :         _data[particleId][energyIdx][pixel] += weight;
      61           2 : }
      62             : 
      63             : 
      64           0 : void ParticleMapsContainer::addParticle(const int particleId, double energy, const Vector3d &p, double weight) {
      65           0 :         double galacticLongitude = atan2(-p.y, -p.x);
      66           0 :         double galacticLatitude =       M_PI / 2 - acos(-p.z / p.getR());
      67           0 :         addParticle(particleId, energy, galacticLongitude, galacticLatitude, weight);
      68           0 : }
      69             : 
      70             : 
      71           1 : std::vector<int> ParticleMapsContainer::getParticleIds() {
      72             :         std::vector<int> ids;
      73           1 :         for(std::map<int, std::map<int, double*> >::iterator pid_iter = _data.begin(); 
      74           2 :                         pid_iter != _data.end(); ++pid_iter) {
      75           1 :                 ids.push_back(pid_iter->first);
      76             :         }
      77           1 :         return ids;
      78             : }
      79             : 
      80             : 
      81           1 : std::vector<double> ParticleMapsContainer::getEnergies(int pid) {
      82             :         std::vector<double> energies;
      83           1 :         if (_data.find(pid) != _data.end()) {
      84           1 :                 for(std::map<int, double*>::iterator iter = _data[pid].begin(); 
      85           2 :                         iter != _data[pid].end(); ++iter) {
      86           1 :                         energies.push_back( idx2Energy(iter->first) / eV );
      87             :                 }
      88             :         }
      89           1 :         return energies;
      90             : }
      91             : 
      92             : 
      93           0 : void ParticleMapsContainer::applyLens(MagneticLens &lens) {
      94             :         // if lens is normalized, this should not be necessary.
      95           0 :         _weightsUpToDate = false;
      96             : 
      97           0 :         for(std::map<int, std::map<int, double*> >::iterator pid_iter = _data.begin(); 
      98           0 :                         pid_iter != _data.end(); ++pid_iter) {
      99           0 :                 for(std::map<int, double*>::iterator energy_iter = pid_iter->second.begin();
     100           0 :                         energy_iter != pid_iter->second.end(); ++energy_iter) {
     101             :                         // transform only nuclei
     102           0 :                         double energy = idx2Energy(energy_iter->first);
     103           0 :                         int chargeNumber = HepPID::Z(pid_iter->first);
     104           0 :                         if (chargeNumber != 0 && lens.rigidityCovered(energy / chargeNumber)) {
     105           0 :                                 lens.transformModelVector(energy_iter->second, energy / chargeNumber);
     106             :                         } else { // still normalize the vectors 
     107           0 :                                 for(size_t j=0; j< _pixelization.getNumberOfPixels() ; j++) {
     108           0 :                                         energy_iter->second[j] /= lens.getNorm();
     109             :                                 }
     110             :                         }
     111             :                 }
     112             :         }
     113           0 : }
     114             : 
     115             : 
     116         421 : void ParticleMapsContainer::_updateWeights() {
     117         421 :         if (_weightsUpToDate)
     118             :                 return;
     119             : 
     120           1 :         for(std::map<int, std::map<int, double*> >::iterator pid_iter = _data.begin(); 
     121           2 :                         pid_iter != _data.end(); ++pid_iter) {
     122           1 :                 _weightsPID[pid_iter->first] = 0;
     123             : 
     124           1 :                 for(std::map<int, double*>::iterator energy_iter = pid_iter->second.begin();
     125           2 :                         energy_iter != pid_iter->second.end(); ++energy_iter)  {
     126             : 
     127           1 :                         _weights_pidEnergy[pid_iter->first][energy_iter->first] = 0;
     128       49153 :                         for(size_t j = 0; j< _pixelization.getNumberOfPixels() ; j++) {
     129       49152 :                                 _weights_pidEnergy[pid_iter->first][energy_iter->first] +=energy_iter->second[j];
     130       49152 :                                 _weightsPID[pid_iter->first]+=energy_iter->second[j];
     131             :                         }
     132           1 :                         _sumOfWeights+=_weights_pidEnergy[pid_iter->first][energy_iter->first];
     133             :                 }
     134             :         }
     135           1 :         _weightsUpToDate = true;
     136             : }
     137             : 
     138             : 
     139           1 : void ParticleMapsContainer::getRandomParticles(size_t N, vector<int> &particleId, 
     140             :         vector<double> &energy, vector<double> &galacticLongitudes,
     141             :         vector<double> &galacticLatitudes) {
     142           1 :         _updateWeights();
     143             : 
     144           1 :         particleId.resize(N);
     145           1 :         energy.resize(N);
     146           1 :         galacticLongitudes.resize(N);
     147           1 :         galacticLatitudes.resize(N);
     148             : 
     149         421 :         for(size_t i=0; i< N; i++) {
     150             :                 //get particle
     151         420 :                 double r = Random::instance().rand() * _sumOfWeights;
     152             :                 std::map<int, double>::iterator iter = _weightsPID.begin();
     153         420 :                 while ((r-= iter->second) > 0) {
     154             :                         ++iter; 
     155             :                 }
     156         420 :                 particleId[i] = iter->first;
     157             :         
     158             :                 //get energy
     159         420 :                 r = Random::instance().rand() * iter->second;
     160         420 :                 iter = _weights_pidEnergy[particleId[i]].begin();
     161         420 :                 while ((r-= iter->second) > 0) {
     162             :                  ++iter; 
     163             :                 }
     164         420 :                 energy[i] = idx2Energy(iter->first) / eV;
     165             : 
     166         420 :                 placeOnMap(particleId[i], energy[i] * eV, galacticLongitudes[i], galacticLatitudes[i]);
     167             :         }
     168           1 : }
     169             : 
     170             : 
     171         420 : bool ParticleMapsContainer::placeOnMap(int pid, double energy, double &galacticLongitude, double &galacticLatitude) {
     172         420 :         _updateWeights();
     173             : 
     174         420 :         if (_data.find(pid) == _data.end()) {
     175             :                 return false;
     176             :         }
     177         420 :         int energyIdx   = energy2Idx(energy);
     178         840 :         if (_data[pid].find(energyIdx) == _data[pid].end()) {
     179             :                 return false;
     180             :         }
     181             : 
     182         420 :         double r = Random::instance().rand() * _weights_pidEnergy[pid][energyIdx];
     183             : 
     184    10268580 :         for(size_t j = 0; j< _pixelization.getNumberOfPixels(); j++) {
     185    10268580 :                 r -= _data[pid][energyIdx][j];
     186    10268580 :                 if (r <= 0) {
     187         420 :                         _pixelization.getRandomDirectionInPixel(j, galacticLongitude, galacticLatitude);
     188         420 :                         return true;
     189             :                 }
     190             :         }
     191             :         return false;
     192             : }
     193             : 
     194             : 
     195           0 : void ParticleMapsContainer::forceWeightUpdate() {
     196           0 :         _weightsUpToDate = false;
     197           0 : }
     198             : 
     199             : } // namespace crpropa