CRSource.H

Go to the documentation of this file.

 
#ifndef _CRSOURCE_H_
#define _CRSOURCE_H_
 
#include <cmath>
#include <gsl/gsl_sf.h>
#include "CRPacket.H"
#include "../Utils/Constants.H"
#include "../Utils/PartTypes.H"
#include "../Utils/RngThread.H"
#include "../Utils/SR.H"
#include "../Utils/Types.H"
 
namespace criptic {
 
  class CRSource {
 
  public:
 
    // Data
    partTypes::pType type;   
    IdxType uniqueID;        
    Real p0;                 
    Real p1;                 
    Real q;                  
    Real k;                  
#ifdef TRACK_PITCH_ANGLE
    Real mu0;                
    Real mu1;                
#endif
    RealVec v;               
    RealVec a;               
    Real T0() const {
      return sr::T_from_p(type, p0);
    }
 
    Real T1() const {
      return sr::T_from_p(type, p1);
    }
 
    void set_T0(const Real T0) {
      p0 = sr::p_from_T(type, T0);
    }
    
    void set_T1(const Real T1) {
      p1 = sr::p_from_T(type, T1);
    }
 
    Real E0() const {
      return sr::E_from_p(type, p0);
    }
 
    Real E1() const {
      return sr::E_from_p(type, p1);
    }
 
    Real R0() const {
      return p0 / partTypes::charge[type];
    }
    
    Real R1() const {
      return p1 / partTypes::charge[type];
    }
    
    Real nRate() const {
      if (p1 <= p0) return k;
      else if (q == -1.0) return k * std::log(p1/p0);
      else return k * (std::pow(p1, q+1) - std::pow(p0, q+1)) / (q+1);
    }
 
    Real pRate() const {
      if (p1 <= p0) return k * p0;
      else if (q == -2.0) return k * std::log(p1/p0);
      else return k * (std::pow(p1, q+2) - std::pow(p0, q+2)) / (q+2);
    }
 
    Real TRate() const {
      if (p1 <= p0) return k * sr::T_from_p(type, p0);
      else {
    Real m = partTypes::mass[type];
    Real m2 = m*m;
    Real p02 = p0*p0;
    Real p12 = p1*p1;
    if (q == -2.0) {
      return k * m * (1/p1 - 1/p0) +
        std::sqrt(1 + m2/p02) - std::sqrt(1 + m2/p12) +
        std::asinh(p1/m) - std::asinh(p0/m);
    } else if (q == 0.0) {
      // Use series expansions to avoid underflows
      Real p1fac = (m2/p12 > 1e-6) ?
        std::log( p1 * (std::sqrt(1 + m2/p12) - 1) ) :
        std::log( m2 / (2 * p1) - m2*m2 / (8 * p1 * p12));
      Real p0fac = (m2/p12 > 1e-6) ?
        std::log( p0 * (std::sqrt(1 + m2/p02) - 1) ) :
        std::log( m2 / (2 * p0) - m2*m2 / (8 * p0 * p02));
      return k * 0.5 * (p1 * std::sqrt(p12 + m2) -
                p0 * std::sqrt(p02 + m2) -
                2 * m * (p1 - p0) -
                m2 * (p1fac - p0fac) );
    } else {
      Real a = -1-q/2.0;
      Real b = (1+q)/2.0;
      Real term1 = m/(q+1) *
        (std::pow(p0, q+1) - std::pow(p1, q+1));
      Real term2 = 0.5 * std::pow(m, q+2) * gsl_sf_beta(a, b) *
        ( gsl_sf_beta_inc( a, b, 1/(1+p02/m2) ) -
          gsl_sf_beta_inc( a, b, 1/(1+p12/m2) ) );
      return k * (term1 + term2);
    }
      }
    }
 
    Real pBar() const {
      return pRate() / nRate();
    }
 
    Real TBar() const {
      return TRate() / nRate();
    }
 
    Real wgt(const Real qSamp) const {
      if (p1 <= p0)
    return k / std::pow(p0, qSamp);
      else if (qSamp == -1.0)
    return k * std::log(p1/p0) / std::pow(p0, qSamp);
      else
    return k * (std::pow(p1, qSamp+1) - std::pow(p0, qSamp+1)) /
      ((qSamp+1) * std::pow(p0, qSamp));
    }
 
    void setLum(const Real lum, const bool codeUnits = false) {
      k = 1.0;
      if (codeUnits) k = lum / TRate();
      else k = (lum / constants::mp_c2) / TRate();
    }
 
    std::vector<CRPacket> draw(const IdxType n,
                   const Real t0,
                   const Real t1,
                   const Real qSamp,
                   const RngThread& rng) const;
  };
 
}
 
#endif
// _CRSOURCE_H_