SimMockEvent.cc

Go to the documentation of this file.

#include "SimMockEvent.h"

#include <vector>
#include <map>

#include <utl/TabulatedFunction.h>

#include <utl/UTMPoint.h>
#include <utl/ReferenceEllipsoid.h>
#include <utl/Vector.h>
#include <utl/AxialVector.h>
#include <utl/TimeStamp.h>

#include <utl/MathConstants.h>
#include <utl/PhysicalConstants.h>
#include <utl/AugerUnits.h>
#include <utl/CoordinateSystem.h>
#include <utl/AugerCoordinateSystem.h>

#include <utl/Trace.h>

#include <evt/Event.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerFRecData.h>
#include <evt/ShowerSimData.h>
#include <evt/GaisserHillas4Parameter.h>
#include <evt/DefaultShowerGeometryProducer.h>

#include <fevt/FEvent.h>
#include <fevt/Eye.h>
#include <fevt/EyeHeader.h>
#include <fevt/Telescope.h>
#include <fevt/EyeRecData.h>
#include <fevt/TelescopeRecData.h>
#include <fevt/TelescopeSimData.h>
#include <fevt/PixelRecData.h>
#include <fevt/Pixel.h>
#include <fevt/FdConstants.h>

#include <det/Detector.h>
#include <atm/Atmosphere.h>

#include <fdet/FDetector.h>
#include <fdet/Eye.h>
#include <fdet/Pixel.h>
#include <fdet/Channel.h>
#include <fdet/Camera.h>
#include <fdet/Telescope.h>

#include <utl/AugerException.h>
#include <utl/Vector.h>

#include <sstream>

// modules to be remote controlled
#include <FdReconstruction/FdProfileReconstructorKG/FdProfileReconstructor.h>

using namespace FdLightCollectionEfficiencyKG;
using namespace fwk;
using namespace utl;
using namespace evt;
using namespace fevt;
using namespace det;
using namespace fdet;
using namespace std;


/*************************************************************************/
SimMockEvent::SimMockEvent(const double xBinning) :
  fVerbosity(0),
  fXBinning(xBinning)
{
  fEvent.MakeFEvent();
  fEvent.MakeSimShower(evt::DefaultShowerGeometryProducer());
}


/*************************************************************************/
void
SimMockEvent::FillSimEvent(const evt::ShowerFRecData& recData)
{
  if (fVerbosity > 0)
    cout << "  --[FdLightCollectionEfficiency]--> Preparing Mock SimEvent..." << endl;

  if (!fEvent.HasSimShower())
    fEvent.MakeSimShower(evt::DefaultShowerGeometryProducer());

  ShowerSimData& simShower = fEvent.GetSimShower();

  const Point& position = recData.GetCorePosition();
  const Vector axis = recData.GetAxis();

  const ReferenceEllipsoid& e = ReferenceEllipsoid::GetWGS84();
  const CoordinateSystemPtr cs = AugerCoordinateSystem::Create(position, e, e.GetECEF());
  
  const double zenith = axis.GetTheta(cs);
  const double azimuth = axis.GetPhi(cs);
  simShower.SetGroundParticleCoordinateSystemZenith(zenith);
  simShower.SetGroundParticleCoordinateSystemAzimuth(azimuth);  
  simShower.MakeGeometry(position);
  
  // Need to do this before SetLongitudinalProfilesFromGH
  simShower.MakeGHParameters(recData.GetGHParameters());
  simShower.SetEnergyCutoff(recData.GetEnergyCutoff());
  SetLongitudinalProfilesFromGH(simShower);
}


/*************************************************************************/
void
SimMockEvent::SetLongitudinalProfilesFromGH(evt::ShowerSimData& simShower)
{
  const VGaisserHillasParameter& ghParam = simShower.GetGHParameters();
  const GaisserHillas4Parameter& gh4Param =
    dynamic_cast<const GaisserHillas4Parameter&>(ghParam);

  simShower.SetXFirst(fmax(0.2*g/cm2, gh4Param.GetShapeParameter(gh::eX0)));
  if (fVerbosity > 0)
    cout << "  X1 is: " << simShower.GetXFirst()/g*cm2 << endl;

  const double enCut = simShower.GetEnergyCutoff();

  TabulatedFunction dEdX;
  TabulatedFunction Ne;

  double depth = 0.2*g/cm2;
  const double xMax = ghParam.GetXMax();

  for (;;) {
    const double age = utl::ShowerAge(depth, xMax);

    double thisdEdX = ghParam.Eval(depth);
    if (age < 1. && thisdEdX < 2.*MeV/g*cm2)
      thisdEdX = 2.*MeV/g*cm2;

    const double thisNe = thisdEdX / utl::EnergyDeposit(age, enCut); // lives in utl::PhysicalFunctions

    if (age > 2. || depth > 3000.*g/cm2)
      break;

    dEdX.PushBack(depth, thisdEdX);
    Ne.PushBack(depth, thisNe);

    depth += fXBinning;
  }

  simShower.MakedEdX(dEdX);
  simShower.MakeLongitudinalProfile(Ne);
}


/*************************************************************************/
bool
SimMockEvent::PrepareEvent(const unsigned int eyeId,
                           const evt::ShowerFRecData& recData,
                           const evt::Event& event,
                           const fevt::Eye& eye)
{
  ostringstream output;
  string oindent = "                                 ";
  if (fVerbosity > -1)
    output << "  --[FdLightCollectionEfficiency]--> Preparing Mock FEvent...\n";

  const fdet::FDetector& detFD = det::Detector::GetInstance().GetFDetector();
  const fdet::Eye& theDetEye = detFD.GetEye(eyeId);

  // Create the event structures
  if (!event.GetFEvent().HasEye(eyeId, ComponentSelector::eHasData))
    return false;

  const int eyeDaqStatus = theDetEye.GetDAQStatus();
  if (eyeDaqStatus == 0) {
    if (fVerbosity > 0)
      output << oindent << "Skipping eye " << eyeId << " (not in DAQ)\n";
    return false;
  }

  if (fVerbosity > 0)
    output << oindent << "Creating eye " << eyeId << " (Telescopes in DAQ:";

  fEvent.GetFEvent().MakeEye(eyeId, fevt::ComponentSelector::eInDAQ);
  fevt::Eye& eyeEvent = fEvent.GetFEvent().GetEye(eyeId, fevt::ComponentSelector::eInDAQ);

  bool hasTelescope = false;
  // Only what is written into the FEvent HERE will be considered by any Modules
  // further donwstream in the efficiency-calculation-sequence!
  for (fdet::Eye::TelescopeIterator iTel = theDetEye.TelescopesBegin(), end = theDetEye.TelescopesEnd();
       iTel != end; ++iTel) {

    const unsigned int telId = iTel->GetId();

    const int telDaqStatus = iTel->GetDAQStatus();
    if (telDaqStatus == 0)
      continue;

    hasTelescope = true;

    if (fVerbosity > 0)
      output << " " << telId;

    eyeEvent.MakeTelescope(telId, fevt::ComponentSelector::eInDAQ);
    fevt::Telescope& tel = eyeEvent.GetTelescope(telId, fevt::ComponentSelector::eInDAQ);
    if (event.GetFEvent().HasEye(eyeId, fevt::ComponentSelector::eHasData)) {
      const fevt::Eye& refEye = event.GetFEvent().GetEye(eyeId, fevt::ComponentSelector::eHasData);
      if (refEye.HasTelescope(telId, fevt::ComponentSelector::eHasData)) {
        // FIXME: It is entirely unclear whether this makes any sense at all.
        //tel.SetTracesStartTime( refEye.GetTelescope(telId, fevt::ComponentSelector::eHasData).GetTracesStartTime() );
        tel.SetTimeOffset( refEye.GetTelescope(telId, fevt::ComponentSelector::eHasData).GetTimeOffset() );
      }
    }

  } // end for telescopes
  output << ")\n";

  cout << output.str() << endl;

  // Direction, core position, etc.
  FillSimEvent(recData);

  // The following code computes the time-at-core which is needed
  // for the simShower. (See also: OfflineKG's SetAxis does the opposite)
  ShowerSimData& simShower = fEvent.GetSimShower();

  const Point& corePosition = simShower.GetPosition();
  const Vector axis = simShower.GetDirection();

  const double T0 = eye.GetRecData().GetTZero();
  utl::TimeStamp eyeTriggerTime =
    eye.GetHeader().GetTimeStamp() - utl::TimeInterval(100.*microsecond);

  const utl::Point& eyePos = detFD.GetEye(eye).GetPosition();
  const Vector eyeToCore = corePosition-eyePos;
  const double distCoreEye = eyeToCore.GetMag();
  const double beta = Angle(-axis, eyeToCore);

  const utl::TimeStamp timeAtCore =
    eyeTriggerTime + TimeInterval(T0) - TimeInterval(distCoreEye/kSpeedOfLight * cos(beta));

  /*cout << "Eye Trigger time: " << eyeTriggerTime << endl;
  cout << "time at core.     " << timeAtCore << endl;
  cout << "T0                " << T0 << endl;
  cout << "beta              " << beta << endl;
  cout << "distcoreeye       " << distCoreEye << endl;
  */

  simShower.MakeTimeStamp(timeAtCore);

  return hasTelescope;
}