Line data    Source code

       1             : #include "crpropa/massDistribution/Massdistribution.h"
       2             : #include <sstream>
       3             : namespace crpropa {
       4             : 
       5           2 : void DensityList::addDensity(ref_ptr<Density> dens) {
       6           2 :         DensityList.push_back(dens);
       7           2 : }
       8             : 
       9           1 : double DensityList::getDensity(const Vector3d &position) const {
      10             :         double n = 0.;
      11           3 :         for (int i = 0; i < DensityList.size(); i++)
      12           2 :                 n += DensityList[i]->getDensity(position);
      13           1 :         return n;
      14             : }
      15             : 
      16           1 : double DensityList::getHIDensity(const Vector3d &position) const {
      17             :         double n = 0.;
      18           3 :         for (int i = 0; i < DensityList.size(); i++)
      19           2 :                 n += DensityList[i]->getHIDensity(position);
      20           1 :         return n;
      21             : }
      22             : 
      23           1 : double DensityList::getHIIDensity(const Vector3d &position) const {
      24             :         double n = 0.;
      25           3 :         for (int i = 0; i < DensityList.size(); i++)
      26           2 :                 n += DensityList[i]->getHIIDensity(position);
      27           1 :         return n;
      28             : }
      29             : 
      30           1 : double DensityList::getH2Density(const Vector3d &position) const {
      31             :         double n = 0.;
      32           3 :         for (int i = 0; i < DensityList.size(); i++)
      33           2 :                 n += DensityList[i]->getH2Density(position);
      34           1 :         return n;
      35             : }
      36             : 
      37           1 : double DensityList::getNucleonDensity(const Vector3d &position) const {
      38             :         double n = 0.;
      39           3 :         for (int i = 0; i < DensityList.size(); i++)
      40           2 :                 n += DensityList[i]->getNucleonDensity(position);
      41           1 :         return n;
      42             : }
      43             : 
      44           0 : std::string DensityList::getDescription() {
      45           0 :         std::stringstream ss; 
      46           0 :         ss << "DensityList with " << DensityList.size() << " modules: \n";
      47           0 :         for (int i = 0; i < DensityList.size(); i++) {
      48           0 :                 ss << "density " << i + 1 << ": " << DensityList[i] -> getDescription();
      49             :         }
      50             :         
      51           0 :         return ss.str();
      52           0 : }
      53             : 
      54             : // ----------- DensityGrid -----------------------------------------------------------------
      55             : 
      56           2 : DensityGrid::DensityGrid(ref_ptr<Grid1f> grid, bool isForHI, bool isForHII, bool isForH2) : 
      57           2 :         grid(grid), isForHI(isForHI), isForHII(isForHII), isForH2(isForH2) {
      58           2 :                 checkAndWarn();
      59           2 :         }
      60             : 
      61           5 : void DensityGrid::checkAndWarn() {
      62           5 :         bool allDeactivated = (isForHI == false) && (isForHII == false) && (isForH2 == false);
      63             :         if (allDeactivated) {
      64           0 :                 KISS_LOG_WARNING << "DensityGrid has all types deactivated."
      65           0 :                         << "In this case all output will be n = 0. \n"
      66           0 :                         << "Please activate the intended particle type. \n";
      67             :         }
      68           5 : }
      69             : 
      70           3 : double DensityGrid::getHIDensity(const Vector3d &position) const {
      71           3 :         if (isForHI)
      72           3 :                 return grid -> interpolate(position);
      73             :         else 
      74             :                 return 0.;
      75             : }
      76             : 
      77           3 : double DensityGrid::getHIIDensity(const Vector3d &position) const {
      78           3 :         if (isForHII) 
      79           0 :                 return grid -> interpolate(position);
      80             :         else
      81             :                 return 0.;
      82             : }
      83             : 
      84           3 : double DensityGrid::getH2Density(const Vector3d &position) const {
      85           3 :         if (isForH2)
      86           0 :                 return grid -> interpolate(position);
      87             :         else
      88             :                 return 0.;
      89             : }
      90             : 
      91           1 : double DensityGrid::getDensity(const Vector3d &position) const {
      92             :         double n = 0;
      93           1 :         n += getHIDensity(position);
      94           1 :         n += getHIIDensity(position);
      95           1 :         n += getH2Density(position);
      96             : 
      97           1 :         return n;
      98             : }
      99             : 
     100           1 : double DensityGrid::getNucleonDensity(const Vector3d &position) const {
     101             :         double n = 0;
     102           1 :         n += getHIDensity(position);
     103           1 :         n += getHIIDensity(position);
     104           1 :         n += 2 * getH2Density(position);
     105             : 
     106           1 :         return n;
     107             : }
     108             : 
     109           2 : bool DensityGrid::getIsForHI() {
     110           2 :         return isForHI;
     111             : }
     112             : 
     113           2 : bool DensityGrid::getIsForHII() {
     114           2 :         return isForHII;
     115             : }
     116             : 
     117           2 : bool DensityGrid::getIsForH2() {
     118           2 :         return isForH2;
     119             : }
     120             : 
     121           1 : void DensityGrid::setIsForHI(bool b) {
     122           1 :         isForHI = b;
     123           1 :         checkAndWarn();
     124           1 : }
     125             : 
     126           1 : void DensityGrid::setIsForHII(bool b) {
     127           1 :         isForHII = b;
     128           1 :         checkAndWarn();
     129           1 : }
     130             : 
     131           1 : void DensityGrid::setIsForH2(bool b) {
     132           1 :         isForH2 = b;
     133           1 :         checkAndWarn();
     134           1 : }
     135             : 
     136           0 : void DensityGrid::setGrid(ref_ptr<Grid1f> grid) {
     137           0 :         this->grid = grid;
     138           0 : }
     139             : 
     140           0 : std::string DensityGrid::getDescription() {
     141           0 :         std::stringstream ss;
     142           0 :         ss << "Density in a given grid \n"; 
     143           0 :         return ss.str();
     144           0 : }
     145             : 
     146             : } //namespace crpropa