TPhoton_o.h

Go to the documentation of this file.

#ifndef RAD_TPhoton_o
#define RAD_TPhoton_o
 
#include "TLorentzVector.h"
#include "TRandom.h"
namespace rb {
  class TPhoton {
  protected:
    TLorentzVector fk;
    double         fknorm;
    double         fNorm;
    double         fE;
    double         fP0;       // auxilary parameter 0
    double         fP1;       // auxilary parameter 1
    double         fX;        // emitted energy
    double         fCosTheta; // cosine of photon on a big angle
    double         fLnD;      // ln((1-beta*c1)/(1+beta*c1))
    double         fBetaI;    // velocity of final particles
    double         fiBetaI;   // inverse velocity of final particles
    double         fLnD_n;    // ln((1-beta*c1)/(1+beta*c1)*(1+beta*c2)/(1-beta*c2))
    double         fBt;       // (1+beta*c1)/(1-beta*c1)
  public:
    TPhoton() {}
    TPhoton( double E );
    TPhoton( double E, double Emin, double Emax );
    ~TPhoton() {}
 
    void SetEnergyRange( double Emin, double Emax );
    void SetThetaRange( double thmin, double thmax );
    void SetCosThetaRange( double c1, double c2 );
 
    const TLorentzVector& GetPhoton() { return fk; }
 
    const TLorentzVector& GetNewPhoton();
    inline double         GetPhotNorm() { return fknorm; }
    inline double         GetEnergy() {
      fX = fP0 * exp( gRandom->Rndm() * fP1 );
      return fX;
    }
    inline double GetENorm() { return fX * fNorm; }
    inline double GetCosThetaF() {
      double d  = fBt * exp( fLnD_n * gRandom->Rndm() );
      fCosTheta = ( d - 1 ) / ( fBetaI * ( d + 1 ) );
      return fCosTheta;
    }
    inline double GetThNormF() {
      double k = fBetaI * fCosTheta;
      return 0.5 * fiBetaI * fLnD_n * ( 1 - k * k );
    }
    inline double GetPhi() { return 2 * M_PI * gRandom->Rndm(); }
    inline double GetPhiNorm() { return 2 * M_PI; }
  };
 
  class TPhotonD : public TPhoton {
  public:
    TPhotonD() : TPhoton(){};
    void Init( const double& delta, const double& beta ) {
      // Initializing photons in final state
      // by distribution 1/x
      fP0   = delta;
      fP1   = 2 / beta;
      fNorm = pow( delta, beta / 2 );
    }
    inline double GetEnergy() {
      double t1 = log( gRandom->Rndm() ) * fP1;
      fX        = fP0 * exp( t1 );
      return fX;
    }
    inline double GetENorm() { return fNorm; }
  };
} // namespace rb
#endif // RAD_TPhoton