SimpleEMComponent.cc

Go to the documentation of this file.

#include "SimpleEMComponent.h"

#include <cmath>

#include <utl/AugerUnits.h>

using namespace tls;
using namespace utl;
using namespace std;

SimpleEMComponent::SimpleEMComponent(utl::Branch /*branch*/)
{
}

SimpleEMComponent::~SimpleEMComponent()
{
}

double
SimpleEMComponent::SignalRatio(double x, double y, double, double theta, double phi)
{
  if (theta < GetThetaMin()) theta = GetThetaMin();
  else if (theta > GetThetaMax()) theta = GetThetaMax();

  //Ines parameterization. 7/11/2005
  //Modified access by Victor. 3/3/2006
  double xp, yp, zp;
  double xShower, yShower/*, zShower*/;

  xp = cos(phi)*x + sin(phi)*y;
  yp =-sin(phi)*x + cos(phi)*y;
  zp = 0.;

  xShower = cos(theta)*xp - sin(theta)*zp;
  yShower = yp;
  //zShower = sin(theta)*xp + cos(theta)*zp;

  double r = sqrt(xShower*xShower + yShower*yShower);

  if (r>4000*m) r = 4000*m;

  double thetaDeg = theta/deg;

  double p0 = -28.6294;
  double p1 = 1.13353;
  double p2 = -0.0149130;
  double p3 = 6.53067e-05;

  double B = p0 + p1*thetaDeg + p2*pow(thetaDeg,2) + p3*pow(thetaDeg,3.);
  double C = 0.204301 + 3.17347*B;

  double q0 = 2.40906e+01;
  double q1 = -3.13028e-01;
  double q2 = -4.95571e+02;
  double q3 = 1.15072e+01;
  double q4 = -6.42068e-02;

  double Rp=exp(q0 + thetaDeg*q1) + q2 + thetaDeg*q3+ pow(thetaDeg,2.)*q4;

  double ratio = Rp/r + C - B*log10(r);

  return ratio;
}


// Configure (x)emacs for this file ...
// Local Variables:
// mode:c++
// compile-command: "make -C .. -k"
// End: