RdStationSimAntennaAlignmentError.cc

Go to the documentation of this file.

#include "RdStationSimAntennaAlignmentError.h"

#include <fwk/CentralConfig.h>
#include <fwk/CoordinateSystemRegistry.h>

#include <utl/ErrorLogger.h>
#include <utl/Reader.h>
#include <utl/config.h>
#include <utl/Trace.h>
#include <utl/TraceAlgorithm.h>
#include <utl/TimeStamp.h>
#include <utl/TimeInterval.h>
#include <utl/AugerUnits.h>
#include <utl/AugerException.h>
#include <utl/FFTDataContainerAlgorithm.h>
#include <utl/Math.h>
#include <utl/Vector.h>
#include <utl/AxialVector.h>
#include <utl/AugerCoordinateSystem.h>
#include <utl/CoordinateSystemPtr.h>
#include <utl/UTMPoint.h>

#include <evt/Event.h>
#include <revt/REvent.h>
#include <revt/Header.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerSimData.h>
#include <evt/ShowerRRecData.h>
#include <evt/ShowerFRecData.h>
#include <evt/ShowerSRecData.h>
#include <revt/Channel.h>

#include <rdet/RDetector.h>
#include <rdet/Station.h>

#include <fwk/LocalCoordinateSystem.h>


#define OUT(x) if ((x) <= fInfoLevel) cerr << "  "

using namespace revt;
using namespace fwk;
using namespace utl;
using namespace std;

namespace RdStationSimAntennaAlignmentError {
    
  RdStationSimAntennaAlignmentError::RdStationSimAntennaAlignmentError()
  : fInfoLevel(0),
    fOnlyForVariedEvents(0), 
    fZenithalAlignmentError(0),
    fAzimuthalAlignmentError(0),
    fRandomSeed(0),
    randgen(),
    normalRandom(0.0, 1.0),
    normalRandomGenerator(randgen, normalRandom),
    uniformRandom(0.0, kPi),
    uniformRandomGenerator(randgen, uniformRandom)
  {}

  RdStationSimAntennaAlignmentError::~RdStationSimAntennaAlignmentError() {}
  
  VModule::ResultFlag RdStationSimAntennaAlignmentError::Init()
  {
    INFO("RdStationSimAntennaAlignmentError::Init()");

    // Read in the configurations of the xml file
    Branch topBranch =
      CentralConfig::GetInstance()->GetTopBranch("RdStationSimAntennaAlignmentError");
    topBranch.GetChild("InfoLevel").GetData(fInfoLevel);
    topBranch.GetChild("OnlyForVariedEvents").GetData(fOnlyForVariedEvents);
    topBranch.GetChild("ZenithalAlignmentError").GetData(fZenithalAlignmentError);
    topBranch.GetChild("AzimuthalAlignmentError").GetData(fAzimuthalAlignmentError);
    topBranch.GetChild("randomSeed").GetData(fRandomSeed);
    
    if (fRandomSeed) {
      OUT(0) << "Seeding the random generator with " << fRandomSeed << endl;
      randgen.seed(fRandomSeed);
    }

    return eSuccess;
  }

  Vector3D RotateZ(Vector3D c, double ang) {
    Vector3D d;
    d[0] =  cos(ang)*c[0] - sin(ang)*c[1]                 ;
    d[1] =  sin(ang)*c[0] + cos(ang)*c[1]                 ;
    d[2] =                                            c[2];
    return d;
  }

  Vector3D RotateX(Vector3D c, double ang) {
    Vector3D d;
    d[0] =            c[0]                                ;
    d[1] =                   cos(ang)*c[1] - sin(ang)*c[2];
    d[2] =                   sin(ang)*c[1] + cos(ang)*c[2];
    return d;
  }

  Vector3D ApplyAlignmentError(Vector3D c, double angZ1, double angZ2, double angX)
  {
    return RotateZ(RotateX(RotateZ(RotateZ(c, angZ1), angZ2), angX), -angZ2);
  }
  

  VModule::ResultFlag RdStationSimAntennaAlignmentError::Run(evt::Event& event)
  {
    INFO("RdStationSimAntennaAlignmentError::Run()");
    
    // Check wether there are events at all
    if(!event.HasREvent()) {
      WARNING("No radio event found!");
      return eContinueLoop;
    }

    REvent& rEvent = event.GetREvent();
    const rdet::RDetector& rDetector =
      det::Detector::GetInstance().GetRDetector();
      
    int the_id = rEvent.GetHeader().GetId();
    // This is a workaround for the SELFAS reader that doesn't implement
    // the reading of the event ID into the rHeader. This should be fixed
    // in the SELFAS reader and then this if clause should be removed.
    if (the_id == 1) {
      the_id = atoi(event.GetHeader().GetId().c_str());
      rEvent.GetHeader().SetId(the_id);
    }
    
    if (fOnlyForVariedEvents and the_id < 99999999) {
      OUT(2) << "Not adding a random rotation: the SD event ID is below 99999999 and the antenna orientation will not be randomized." << endl;
      return eSuccess;
    }
    
    
    // loop through stations
    for (REvent::StationIterator sIt = rEvent.StationsBegin();
         sIt != rEvent.StationsEnd(); ++sIt) {

      double angZ1 = normalRandomGenerator()*fAzimuthalAlignmentError;
      double angZ2 = uniformRandomGenerator();
      double angX  = normalRandomGenerator()*fZenithalAlignmentError;
      
      OUT(3) << "Rotating the trace slightly to account for the error in the alignment of the antenna with the earth magnetic field..." << endl;
      OUT(3) << "Rotation1: about the Z-axis with " << angZ1/degree << " deg drawn from a normal distribution with with sigma=" << fAzimuthalAlignmentError/degree << " deg." << endl;
      OUT(3) << "Rotation2: about the Z-axis with " << angZ2/degree << " deg drawn from a uniform distribution between 0 and Pi." << endl;
      OUT(3) << "Rotation3: about the X-axis with " << angX/degree << " deg drawn from a normal distribution with sigma=" << fZenithalAlignmentError/degree << " deg." << endl;
      OUT(3) << "Rotation4: back around the Z-axis with " << -angZ2/degree << " deg (inverse of rotation2)." << endl;
      
      StationTimeSeries &TimeSeries = sIt->GetFFTDataContainer().GetTimeSeries();
      for (unsigned int i = 0; i < TimeSeries.GetSize(); ++i) {
        TimeSeries[i] = ApplyAlignmentError(TimeSeries[i], angZ1, angZ2, angX);
      }
      
    }
    return eSuccess;
  }

  VModule::ResultFlag
  RdStationSimAntennaAlignmentError::Finish()
  {

    INFO("RdStationSimAntennaAlignmentError::Finish()");

    return eSuccess;
  }

}