testRDetector.cc

Go to the documentation of this file.

// Test giving out infos on Stations

#include <iostream>
#include <sstream>
#include <string>
#include <vector>
#include <det/Detector.h>
#include <rdet/RDetector.h>
#include <rdet/Station.h>
#include <rdet/Channel.h>
#include <sdet/SDetector.h>
#include <sdet/Station.h>
#include <fwk/CentralConfig.h>
#include <utl/TimeStamp.h>
#include <utl/UTCDateTime.h>
#include <cppunit/extensions/HelperMacros.h>
#include <tst/Verify.h>
#include <utl/Point.h>
#include <utl/ReferenceEllipsoid.h>
#include <fwk/CoordinateSystemRegistry.h>
#include <fwk/LocalCoordinateSystem.h>
#include <TROOT.h>
#include <TH2F.h>
#include <TFile.h>
#include <TCanvas.h>
#include <TPaveText.h>
#include <TAttText.h>

using namespace det;
using namespace utl;
using namespace tst;
using namespace rdet;
using std::cout;
using std::endl;

#define ASSERT_CLOSE(x, y, eps) CPPUNIT_ASSERT(Verify<CloseTo>(x, y, eps))
#define ASSERT_EQUAL(x, y) CPPUNIT_ASSERT(Verify<Equal>(x, y))

class RDetectorTest : public CppUnit::TestFixture  {

  CPPUNIT_TEST_SUITE(RDetectorTest);
  CPPUNIT_TEST(testStationListManager);
  CPPUNIT_TEST(testRModel);
  CPPUNIT_TEST(testAntennaType);
  CPPUNIT_TEST_SUITE_END();

public:
  void
  setUp()
  {
    ErrorLogger::GetInstance().SetVerbosity(Verbosity::eVerbose);

    // CentralConfig must be created first.
    //
    /*CentralConfig* fCentralConfig =*/
    fwk::CentralConfig::GetInstance(BOOTSTRAPFILE);

    // Create a detector
    //
    /*Detector& theDetector =*/
    Detector::GetInstance();
  }

  void tearDown() { }

  void
  testAntennaType()
  {
    cout << "Test of the AntennaTypes of RDetector" << endl;
    
    Detector& detector = Detector::GetInstance();
    const TimeStamp detTime = UTCDateTime(2008, 11, 10, 16, 45, 12).GetTimeStamp();
    detector.Update(detTime);

    const rdet::RDetector& rDetector = det::Detector::GetInstance().GetRDetector();
    const Station& rStation = rDetector.GetStation(102);
    std::pair<std::complex<double>, std::complex<double>> cEffectiveAntennaHeight;

    cout << "test channel 1" << endl;
    const Channel& rChannel1  = rStation.GetChannel(1);
    for(unsigned int i = 0; i < 100; i+=5)
    {
        cEffectiveAntennaHeight =
                  rChannel1.GetElectricFieldResponse(0.2, 0.3, double(i) * megahertz, "LinearInterpolation");
        cout << double(i) * megahertz / (giga * hertz) << " GHz \t-> e_theta_amp : "
             << abs(cEffectiveAntennaHeight.first) / meter << " m, e_theta_phase : "
             << arg(cEffectiveAntennaHeight.first) / degree << " degrees" << endl;
    }
    cEffectiveAntennaHeight =
                  rChannel1.GetElectricFieldResponse(0.2, 0.3, 50 * megahertz, "LinearInterpolation");
    CPPUNIT_ASSERT(Verify<Equal>(abs(cEffectiveAntennaHeight.first) > 0, true));

    cout << "test channel 2" << endl;
    const Channel& rChannel2 = rStation.GetChannel(2);
    for(unsigned int i = 0; i < 100; i+=5)
    {
        cEffectiveAntennaHeight =
                  rChannel2.GetElectricFieldResponse(0.2, 0.3, double(i) * megahertz, "LinearInterpolation");
        cout << double(i) * megahertz / (giga * hertz) << " GHz \t-> e_theta_amp : "
             << abs(cEffectiveAntennaHeight.first) / meter << " m, e_theta_phase : "
             << arg(cEffectiveAntennaHeight.first) / degree << " degrees" << endl;
    }
    cEffectiveAntennaHeight =
                  rChannel2.GetElectricFieldResponse(0.2, 0.3, 50 * megahertz, "LinearInterpolation");
    CPPUNIT_ASSERT(Verify<Equal>(abs(cEffectiveAntennaHeight.first) > 0, true));

  };

  void
  testStationListManager()
  {
    Detector& detector = Detector::GetInstance();
    const TimeStamp detTime = UTCDateTime(2008, 11, 10, 16, 45, 12).GetTimeStamp();
    detector.Update(detTime);

    // try to get the name of the RDetector.Station id=1
    const RDetector& RDetector = detector.GetRDetector();
    const Station& RStation = RDetector.GetStation(101);
    cout << "Name of rdet::Station with Id 1: "
         << RStation.GetName()
         << endl;

    // get an SD detector (Olaia) as a reference
    /*
    const sdet::SDetector& SDetector = detector.GetSDetector();
    const sdet::Station& OlaiaStation = SDetector.GetStation(1325); // get tank Olaia (infill tank close to BLS)
    const utl::Point OlaiaLocation = OlaiaStation.GetPosition();
    */

    /*Oliva: get antennas positions
      (originset to Station with id=1)
      then prepare the histo to be filled
    */
    utl::CoordinateSystemPtr localcoor = RDetector.GetStation(101).GetLocalCoordinateSystem();
    // get reference coordinate system (usually PampaAmarilla)
    utl::CoordinateSystemPtr referenceCS= detector.GetReferenceCoordinateSystem();
    TH2F *histo = new TH2F("histo","Antennas Position", 1000, -27100., -25600., 1000, -6000., -4600.);
    TFile antennasposition("antennasposition.root","recreate");
    TCanvas* RSPlot = new TCanvas("RSPlot","RSPlot");
    RSPlot->cd();
    histo->Draw();

    /*
    // print out coordinates of Olaia relative to radio station 1
    std::cout << " Olaia coordinates relative to first radio station: " << OlaiaLocation.GetX(localcoor) << " " << OlaiaLocation.GetY(localcoor) << " " << OlaiaLocation.GetZ(localcoor) << "\n" << endl;
    */

    //Oliva @brief: iterator over radio stations to get infos on them
    for (rdet::RDetector::StationIterator rdIt = RDetector.StationsBegin(); rdIt != RDetector.StationsEnd(); ++rdIt) {
      std::string StationName = rdIt->GetName();
      utl::Point Position = rdIt->GetPosition();
      histo->Fill(Position.GetX(localcoor), Position.GetY(localcoor));

      cout << " Station Name:                                 " << StationName << endl;
      cout << " Position (x,y,z) meters:                      " 
           << Position.GetX(localcoor) << ", " <<  Position.GetY(localcoor) << ", " << Position.GetZ(localcoor) << endl;
      cout << " Station Spher. coo. (theta, phi):             " 
           << Position.GetTheta(localcoor) << ", " <<  Position.GetPhi(localcoor) <<endl;
      cout << " Radius r (cylind. coo.) from origin (m):      " << Position.GetRho(localcoor) << endl;
      cout << " Radius R (spheri. coo.) from origin (m):      " << Position.GetR(localcoor) << endl;
      /*Oliva: the cylindric. coor. radius r
        must b equal to "dist" defined right down here.
        For a check UNCOMMENT next lines!
      */
      //  double dist =
      //  sqrt(Position.GetX(localcoor)*Position.GetX(localcoor)+Position.GetY(localcoor)*Position.GetY(localcoor));
      //  cout <<" Dist. from origin sqrt(x^2+y^2) -CHECK-  (m): " << dist  <<endl;

      /*Oliva: Creates the RSPlot.png file
        which displays the Radio Station positions
        and their funny names.
      */
      std::stringstream StationText("");
      StationText << "   " << rdIt->GetId() << "  " << StationName;
      TText* RStext = new TText(Position.GetX(referenceCS), Position.GetY(referenceCS), StationText.str().c_str());
      RStext->SetTextAlign(20);
      RStext->SetTextAngle(35);
      RStext->SetTextSize(.04);
      RStext->Draw();
    }

    histo->SetStats(false);
    histo->SetMarkerStyle(2);
    histo->SetFillColor(kWhite);
    histo->SetMarkerSize(1.);
    histo->Draw("same");
    histo->Write();
    RSPlot->SaveAs("RSPlot.png");
  }


  void
  testRModel()
  {
    Detector& detector = Detector::GetInstance();
    const TimeStamp detTime = UTCDateTime(2008,11,10,16,45,12).GetTimeStamp();
    detector.Update(detTime);

    // try to get RDetector.GetStation
    const RDetector& RDetector = detector.GetRDetector();
    const Station&  RStation  = RDetector.GetStation(103);
    const Channel&  RChannel  = RStation.GetChannel(1);
    cout << "Lower Design Freq. of rdet::Channel of Station id=" << RStation.GetId()
        << " and Channel=" << RChannel.GetId()
        << ": " << RChannel.GetDesignLowerFreq() / utl::megahertz << " MHz"
         << " GetOrientationAzimuth: " << RChannel.GetOrientationAzimuth() / utl::degree
        << endl;

    const Channel& RChannel2 = RStation.GetChannel(2);
    cout << "Lower Design Freq. of rdet::Channel of Station id=" <<  RStation.GetId()
        << " and Channel=" << RChannel2.GetId()
        << ": " << RChannel2.GetDesignLowerFreq() / utl::megahertz << " MHz"
         << " GetOrientationAzimuth/degree: "
         << RChannel2.GetOrientationAzimuth() / utl::degree
         << " GetOrientationZeta/degree: "
         << RChannel2.GetOrientationZeta() / utl::degree
         << " GetPosition/m: "
        << ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84).PointToLatitudeLongitudeHeight(RChannel2.GetPosition()).get<2>()/utl::meter
        << endl;

    ResponseMap rMap = RChannel2.GetResponseMap();

    std::string RHRME("");
    cout << "Got ResponseMap for Station  id=" <<  RStation.GetId()
         << " and Channel=" << RChannel2.GetId()
         << endl;
    for (ResponseMap::const_iterator rMapIter=rMap.Begin(); rMapIter!=rMap.End(); ++rMapIter)
      {
        cout << "\n ResponseID:     " << rMapIter->first
             << "\n ResponseWeight: " << rMapIter->second
             << "\n Nr of Points in Profile: " << RDetector.GetHardwareResponseProfile(rMapIter->first).GetNPoints()
             << endl;
        RHRME = rMapIter->first;
      };

    cout << "Testiing HardwareProfile: " << RHRME << endl;
    TabulatedFunctionComplexLgAmpPhase tComplexLgAmpPhase = RDetector.GetHardwareResponseProfile(RHRME);
    cout << "\nNumber of Points: "
         << tComplexLgAmpPhase.GetNPoints()
         << "\n3rd x-value in MHz: "
         << tComplexLgAmpPhase.GetX(2)
         << " Log of Amplitude: "
         << tComplexLgAmpPhase.GetY(2).GetLgAmplitude()
         << ", Phase: "
         << tComplexLgAmpPhase.GetY(2).GetPhase()
         << "\nLog of Amplitude at 25 MHz: "
         << tComplexLgAmpPhase.Y(25*utl::megahertz).GetLgAmplitude()
         << ", Phase/degree: "
         << tComplexLgAmpPhase.Y(25*utl::megahertz).GetPhase() / utl::degree
         << "\nLog of Amplitude at 30 MHz: "
         << tComplexLgAmpPhase.Y(30*utl::megahertz).GetLgAmplitude()
         << ", Phase/degree: "
         << tComplexLgAmpPhase.Y(30*utl::megahertz).GetPhase() / utl::degree
         << endl;

  }
};

CPPUNIT_TEST_SUITE_REGISTRATION(RDetectorTest);