JPowerLawFlux.hh

Go to the documentation of this file.

#ifndef __JAANET__JPOWERLAWFLUX__
#define __JAANET__JPOWERLAWFLUX__
 
#include "JLang/JEquals.hh"
#include "JLang/JClonable.hh"
 
#include "Jeep/JProperties.hh"
 
#include "JAAnet/JFlux.hh"
#include "JAAnet/JAAnetToolkit.hh"
 
 
/**
 * \author bjung
 */
 
namespace JAANET {}
namespace JPP { using namespace JAANET; }
 
namespace JAANET {
 
  using JLANG::JEquals;
  using JLANG::JClonable;
 
  using JEEP::JProperties;
  
 
  /**
   * Example function object for computing power-law flux.
   */
  struct JPowerLawFlux final :
    public JEquals<JPowerLawFlux>,
    public JClonable<JFlux, JPowerLawFlux>
  {
    /**
     * Default constructor.
     */
    JPowerLawFlux() :
      normalisation(1.0),
      spectralIndex(0.0)
    {}
    
    
    /**
     * Constructor.
     *
     * \param  normalisation         normalisation   [GeV * m^-2 * sr^-1 * s^-1]
     * \param  spectralIndex         spectral index
     */
    JPowerLawFlux(const double normalisation,
                  const double spectralIndex) :
      normalisation(normalisation),
      spectralIndex(spectralIndex)
    {}
 
 
    /**
     * Check whether this power-law flux is valid.
     *
     * \return                 true if valid; else false
     */
    bool is_valid() const override final
    {
      return normalisation > 0.0;
    }
    
    
    /**
     * Get flux of given event.
     *
     * \param  evt             event
     * \return                 flux                  [GeV * m^-2 * sr^-1 * s^-1]
     */
    double getFactor(const Evt& evt) const override final
    {
      using namespace std;
      using namespace JPP;
 
      const Trk& primary = get_primary(evt);
        
      return normalisation * pow(primary.E, -spectralIndex);
    }
 
 
    /**
     * Check if this flux is identical to the given flux.
     *
     * \param  object                power-law flux object
     * \return                       true if this flux is identical to given flux; else flase
     */
    bool equals(const JPowerLawFlux& object) const
    {
      return (object.normalisation == this->normalisation &&
              object.spectralIndex == this->spectralIndex);
    }
 
 
    /**
     * Get properties of this class.
     *
     * \param  eqpars          equation parameters
     */
    JProperties getProperties(const JEquationParameters& eqpars = JEvtWeightFactor::getEquationParameters()) override final
    {
      return JPowerLawFluxHelper(*this, eqpars);
    }
 
 
    /**
     * Get properties of this class.
     *
     * \param  eqpars          equation parameters
     */
    JProperties getProperties(const JEquationParameters& eqpars = JEvtWeightFactor::getEquationParameters()) const override final
    {
      return JPowerLawFluxHelper(*this, eqpars);      
    }
 
 
    /**
     * Stream input.
     *
     * \param  in              input stream
     * \return                 input stream
     */
    std::istream& read(std::istream& in) override final
    {
      using namespace std;
 
      streampos pos = in.tellg();
      
      if (!(in >> normalisation >> spectralIndex)) {
 
        in.clear();
        in.seekg(pos);
 
        JProperties properties = getProperties();
 
        in >> properties;
      }
 
      check_validity();
 
      return in;
    }
 
 
  private:
 
    /**
     * Auxiliary class for I/O of power-law flux.
     */
    struct JPowerLawFluxHelper :
      public JProperties
    {
      /**
       * Constructor.
       *
       * \param  flux          power-law flux
       * \param  eqpars        equation parameters 
       */
      template<class JPowerLawFlux_t>
      JPowerLawFluxHelper(JPowerLawFlux_t&           flux,
                          const JEquationParameters& eqpars) :
        JProperties(eqpars, 1)
      {
        (*this)[JEvtWeightFactor::getTypeKey()] = "power-law flux";
        
        this->insert(gmake_property(flux.normalisation));
        this->insert(gmake_property(flux.spectralIndex));
      }
    };
    
    
    double normalisation;   //!< normalisation      [GeV * m^-2 * sr^-1 * s^-1]
    double spectralIndex;   //!< spectral index     >= 0
  };
}
 
#endif