/*! a toy Gaussian point source compoment */

struct GaussianPointSource : public Component
{

  Parameter ra = { "right assencion of the source (radians)",
                   0,
                   -pi, pi, 0.1
                 } ;

  Parameter dec = { "declinations of the source (radians)",
                    0,
                    -pi/2, pi/2, 0.1
                  } ;

  Parameter sigma = { "width of the point-spread-function (radians)",
                      1 * pi / 180, // 1 degree
                      0.01 / 57.2, 10 / 57.2, 0.1 / 57.2
                    } ;


  GaussianPointSource(string name = "gauspoint" ) : Component(name) {}

  /*! return dP/dx */
  double Rayleigh( double x, double sigma )
  {
    const double sigma2 = sigma * sigma;
    return x / sigma2 * exp( -x * x / (2 * sigma2) );
  }

  double dP_dlogE( double logE )
  {
    // assume some toy thing
    return TMath::Gaus(logE, 3, 1.5, true);
  }

  virtual double dN_dOmega_dlogE(SearchEvent &sev)
  {
    EquatorialCoords S( ra, dec ) ;
    double angle = sev.coordinates.distance( S );
    double dPdOmega = Rayleigh( angle, sigma ) / max( angle, 1e-10);

    double r = norm.value * dPdOmega * dP_dlogE ( log10(sev.Eproxy) );

    return r;


  }

  virtual SearchEvent generate_event()
  {

    double u = gRandom->Rndm();
    double angle = sigma * sqrt( -2 * log(u) );
    double logE = gRandom->Gaus( 3, 1.5 );

    auto eq = EquatorialCoords(ra, dec);

    Vec v = eq.asVec();
    Mat m(v); // matrix that rotates (0,0,1) to v

    double theta = angle;
    double phi   = 2 * pi * gRandom->Rndm();

    Vec r = m * Vec().set_angles( theta, phi );

    double ra_  = r.phi();
    double dec_ = pi / 2 - r.theta();

    return SearchEvent( ra_, dec_, pow(10, logE) );
   // return SearchEvent();
  }

  // -----boilerplate-----
  virtual GaussianPointSource *clone() const
  {
    return new GaussianPointSource(*this);
  }

  virtual ~GaussianPointSource() {}
  ClassDef( GaussianPointSource, 1)
};