Program Listing for File PhotoPionProduction.h

↰ Return to documentation for file (include/crpropa/module/PhotoPionProduction.h)

#ifndef CRPROPA_PHOTOPIONPRODUCTION_H
#define CRPROPA_PHOTOPIONPRODUCTION_H

#include "crpropa/Module.h"
#include "crpropa/PhotonBackground.h"

#include <vector>

namespace crpropa {
struct SophiaEventOutput {
        int nParticles;
        std::vector<double> energy;
        std::vector<int> id;
};

class PhotoPionProduction: public Module {

protected:
        ref_ptr<PhotonField> photonField;
        std::vector<double> tabLorentz;
        std::vector<double> tabRedshifts;
        std::vector<double> tabProtonRate;
        std::vector<double> tabNeutronRate;
        double limit;
        bool havePhotons;
        bool haveNeutrinos;
        bool haveElectrons;
        bool haveAntiNucleons;
        bool haveRedshiftDependence;
        std::string interactionTag = "PPP";

        // called by: sampleEps
        // - input: s [GeV^2]
        // - output: (s-p^2) * sigma_(nucleon/gamma) [GeV^2 * mubarn]
        double functs(double s, bool onProton) const;

        // called by: sampleEps, gaussInt
        // - input: photon energy eps [eV], Ein [GeV]
        // - output: probability to encounter photon of energy eps
        double probEps(double eps, bool onProton, double Ein, double z) const;

        double epsMinInteraction(bool onProton, double Ein) const;

        double momentum(bool onProton, double Ein) const;

        // called by: functs
        // - input: photon energy [eV]
        // - output: crossection of nucleon-photon-interaction [mubarn]
        double crossection(double eps, bool onProton) const;

        // called by: crossection
        // - input: photon energy [eV], threshold [eV], max [eV], unknown [no unit]
        // - output: unknown [no unit]
        double Pl(double eps, double xth, double xMax, double alpha) const;

        // called by: crossection
        // - input: photon energy [eV], threshold [eV], unknown [eV]
        // - output: unknown [no unit]
        double Ef(double eps, double epsTh, double w) const;

        // called by: crossection
        // - input: cross section [µbarn], width [GeV], mass [GeV/c^2], rest frame photon energy [GeV]
        // - output: Breit-Wigner crossection of a resonance of width Gamma
        double breitwigner(double sigma0, double gamma, double DMM, double epsPrime, bool onProton) const;

        // called by: probEps, crossection, breitwigner, functs
        // - input: is proton [bool]
        // - output: mass [Gev/c^2]
        double mass(bool onProton) const;

        // - output: [GeV^2] head-on collision
        double sMin() const;

        bool sampleLog = true;
        double correctionFactor = 1.6; // increeses the maximum of the propability function


public:
        PhotoPionProduction(
                ref_ptr<PhotonField> photonField,
                bool photons = false,
                bool neutrinos = false,
                bool electrons = false,
                bool antiNucleons = false,
                double limit = 0.1,
                bool haveRedshiftDependence = false);

        // set the target photon field
        void setPhotonField(ref_ptr<PhotonField> photonField);

        // decide if secondary photons are added to the simulation
        void setHavePhotons(bool b);

        // decide if secondary neutrinos are added to the simulation
        void setHaveNeutrinos(bool b);

        // decide if secondary electrons are added to the simulation
        void setHaveElectrons(bool b);

        // decide if secondary anti nucleons are added to the simulation
        void setHaveAntiNucleons(bool b);

        // decide if redshift dependent tabulated loss rates are used
        void setHaveRedshiftDependence(bool b);

        void setLimit(double limit);

        void setInteractionTag(std::string tag);

        void initRate(std::string filename);

        double nucleonMFP(double gamma, double z, bool onProton) const;

        double nucleiModification(int A, int X) const;
        void process(Candidate *candidate) const;
        void performInteraction(Candidate *candidate, bool onProton) const;

        double lossLength(int id, double gamma, double z = 0);

        SophiaEventOutput sophiaEvent(bool onProton, double Ein, double eps) const;

        double sampleEps(bool onProton, double E, double z) const;

        double probEpsMax(bool onProton, double Ein, double z, double epsMin, double epsMax) const;

        // using log or lin spacing of photons in the range between epsMin and
        // epsMax for computing the maximum probability of photons in field
        void setSampleLog(bool log);

        // given the discrete steps to compute the maximum interaction probability pEpsMax
        // of photons in field, the real pEpsMax may lie between the descrete tested photon energies.
        // A correction factor can be set to increase pEpsMax by that factor
        void setCorrectionFactor(double factor);

        ref_ptr<PhotonField> getPhotonField() const;
        bool getHavePhotons() const;
        bool getHaveNeutrinos() const;
        bool getHaveElectrons() const;
        bool getHaveAntiNucleons() const;
        bool getHaveRedshiftDependence() const;
        double getLimit() const;
        bool getSampleLog() const;
        double getCorrectionFactor() const;
        std::string getInteractionTag() const;
};
} // namespace crpropa

#endif // CRPROPA_PHOTOPIONPRODUCTION_H