Line data    Source code

       1              : #ifndef CRPROPA_PHOTODISINTEGRATION_H
       2              : #define CRPROPA_PHOTODISINTEGRATION_H
       3              : 
       4              : #include "crpropa/Module.h"
       5              : #include "crpropa/PhotonBackground.h"
       6              : 
       7              : #include <vector>
       8              : #include <map>
       9              : 
      10              : namespace crpropa {
      11              : /**
      12              :  * \addtogroup EnergyLosses
      13              :  * @{
      14              :  */
      15              :  
      16              : /**
      17              :  @class PhotoDisintegration
      18              :  @brief Photodisintegration of nuclei by background photons.
      19              :  */
      20              : class PhotoDisintegration: public Module {
      21              : private:
      22              :         ref_ptr<PhotonField> photonField;
      23              :         double limit; // fraction of mean free path for limiting the next step
      24              :         bool havePhotons;
      25              :         std::string interactionTag = "PD";
      26              : 
      27       454129 :         struct Branch {
      28              :                 int channel; // number of emitted (n, p, H2, H3, He3, He4)
      29              :                 std::vector<double> branchingRatio; // branching ratio as function of nucleus Lorentz factor
      30              :         };
      31              : 
      32      1003680 :         struct PhotonEmission {
      33              :                 double energy; // energy of emitted photon [J]
      34              :                 std::vector<double> emissionProbability; // emission probability as function of nucleus Lorentz factor
      35              :         };
      36              : 
      37              :         std::vector<std::vector<double> > pdRate; // pdRate[Z * NUCLEAR_NSTRIDE + N] = total interaction rate
      38              :         std::vector<std::vector<Branch> > pdBranch; // pdBranch[Z * NUCLEAR_NSTRIDE + N] = branching ratios
      39              :         mutable std::map<int, std::vector<PhotonEmission> > pdPhoton; // map of emitted photon energies and photon emission probabilities
      40              : 
      41              :         static const double lgmin; // minimum log10(Lorentz-factor)
      42              :         static const double lgmax; // maximum log10(Lorentz-factor)
      43              :         static const size_t nlg; // number of Lorentz-factor steps
      44              : 
      45              : public:
      46              :         /** Constructor.
      47              :          @param photonField  target photon field
      48              :          @param havePhotons  if true, add secondary photons as candidates
      49              :          @param limit        step size limit as fraction of mean free path
      50              :          @param superheavy   if true, also load rate/branching _superheavy files extending
      51              :                              coverage from A≥56 up to Pb; files are loaded if present and
      52              :                              silently skipped with a warning if not yet generated
      53              :          */
      54              :         PhotoDisintegration(ref_ptr<PhotonField> photonField, bool havePhotons = false, double limit = 0.1, bool superheavy = false);
      55              : 
      56              :         /** Set the target photon field and reload all interaction tables.
      57              :          @param photonField  target photon field
      58              :          @param superheavy   if true, also load rate_<field>_superheavy.txt and
      59              :                              branching_<field>_superheavy.txt (logged as warning if absent)
      60              :          */
      61              :         void setPhotonField(ref_ptr<PhotonField> photonField, bool superheavy = false);
      62              : 
      63              :         // decide if secondary photons are added to the simulation
      64              :         void setHavePhotons(bool havePhotons);
      65              : 
      66              :         /** Limit the propagation step to a fraction of the mean free path
      67              :          * @param limit fraction of the mean free path
      68              :          */
      69              :         void setLimit(double limit);
      70              : 
      71              :         /** set a custom interaction tag to trace back this interaction
      72              :          * @param tag string that will be added to the candidate and output
      73              :          */
      74              :         void setInteractionTag(std::string tag);
      75              :         std::string getInteractionTag() const;
      76              : 
      77              :         void initRate(std::string filename);
      78              :         void initBranching(std::string filename);
      79              :         void initPhotonEmission(std::string filename);
      80              : 
      81              :         void process(Candidate *candidate) const;
      82              :         void performInteraction(Candidate *candidate, int channel) const;
      83              : 
      84              :         /**
      85              :          Calculates the loss length E dx/dE in [m] physical distance.
      86              :          This is not used in the simulation.
      87              :          @param id              PDG particle id
      88              :          @param gamma   Lorentz factor of particle
      89              :          @param z               redshift
      90              :          @returns E dx/dE [in meters]
      91              :          */
      92              :         double lossLength(int id, double gamma, double z = 0);
      93              : };
      94              : 
      95              : /** @}*/
      96              : } // namespace crpropa
      97              : 
      98              : #endif // CRPROPA_PHOTODISINTEGRATION_H