testREASFile.cc

Go to the documentation of this file.

#include <io/REASFile.h>
#include <io/REASIOException.h>
#include <io/IoCodes.h>

#include <evt/Event.h>
#include <evt/ShowerSimData.h>
#include <evt/VGaisserHillasParameter.h>
#include <evt/RadioSimulation.h>
#include <evt/SimRadioPulse.h>

#include <utl/AugerUnits.h>
#include <utl/Particle.h>

#include <fwk/CentralConfig.h>

#include <det/Detector.h>

#include <tst/Verify.h>
#include <cppunit/extensions/HelperMacros.h>

#include <iostream>

using io::REASFile;
using namespace fwk;
using namespace utl;
using namespace tst;


class testREASFile : public CppUnit::TestFixture {

  CPPUNIT_TEST_SUITE (testREASFile);
  CPPUNIT_TEST (testConstruct);
  //CPPUNIT_TEST_EXCEPTION (testConstructException, io::REASIOException);
  //CPPUNIT_TEST (testOpen);
  //CPPUNIT_TEST_EXCEPTION (testOpenNonExisting, io::REASIOException);
  //CPPUNIT_TEST_EXCEPTION (testOpenForWrite, io::REASIOException);
  //CPPUNIT_TEST (testClose);
  //CPPUNIT_TEST_EXCEPTION (testWrite, io::REASIOException);
  //CPPUNIT_TEST (testFindEvent);
  //CPPUNIT_TEST (testGotoPosition);
  CPPUNIT_TEST (testGetNEvents);
  CPPUNIT_TEST (testRead);
  //CPPUNIT_TEST_EXCEPTION (testGetNEventsClosedException, io::REASIOException);
  //CPPUNIT_TEST_EXCEPTION(testIDThrow, io::REASIOException);
  //CPPUNIT_TEST(testSizeBlock);
  CPPUNIT_TEST_SUITE_END();

public:
  // shared data for tests

  REASFile* f;
  evt::Event event;

  void
  setUp()
  {
    f = new REASFile();
    event = evt::Event();
  }

  void
  tearDown()
  {
    delete f;
  }

  void
  testConstruct()
  {
    REASFile* file;
    CPPUNIT_ASSERT(file = new REASFile(TESTFILE));
    delete file;
  }

  /*
  void
  testConstructException()
  {
    Expected();
    REASFile f(TESTFILE, io::eWrite);
  }

  void
  testOpen()
  {
    f->Open(TESTFILE);
  }

  void testOpenNonExisting()
  {
    Expected();
    f->Open("NoSuchREASShowerFile.root");
  }

  void
  testOpenForWrite()
  {
    Expected();
    f->Open(TESTFILE, io::eWrite);
  }

  void
  testClose()
  {
    f->Open(TESTFILE);
    f->Close();
  }

  void
  testFindEvent()
  {
    f->Open(TESTFILE);
    CPPUNIT_ASSERT(f->FindEvent(2) == io::eFail);
    CPPUNIT_ASSERT(f->FindEvent(0) != io::eFail);
  }

  void
  testGotoPosition()
  {
    f->Open(TESTFILE);
    CPPUNIT_ASSERT(f->GotoPosition(0) != io::eFail);
    CPPUNIT_ASSERT(f->GotoPosition(2) == io::eFail);
  }

  void
  testWrite()
  {
    f->Open(TESTFILE);
    Expected();
    f->Write(event);
  }
*/

  void
  testGetNEvents()
  {
    f->Open(TESTFILE);
    CPPUNIT_ASSERT(f->GetNEvents() == 1);
  }

  /*
  void
  testGetNEventsClosedException()
  {
    Expected();
    f->GetNEvents();
  }
*/

  void
  testRead()
  {
    fwk::CentralConfig::GetInstance(BOOTSTRAPFILE);

    f->Open(TESTFILE);
    int nEventsInFile = f->GetNEvents();
    int nEventsSeen = 0;

    while (f->Read(event) != io::eEOF) {

      ++nEventsSeen;
      evt::ShowerSimData& shower = event.GetSimShower();

      const CoordinateSystemPtr refCS = det::Detector::GetInstance().GetReferenceCoordinateSystem();
      shower.MakeGeometry(utl::Point(0,0,0,refCS));

      // test if CORSIKA variables are set up correctly
      CPPUNIT_ASSERT(Verify<Equal>(shower.GetPrimaryParticle(), int(Particle::eProton)));
      CPPUNIT_ASSERT(Verify<CloseTo>(shower.GetEnergy(), 1.36772e+18*eV));

      const utl::CoordinateSystemPtr localCS = shower.GetLocalCoordinateSystem();
      const double zenith = (-shower.GetDirection()).GetTheta(localCS);
      const double azimuth = (-shower.GetDirection()).GetPhi(localCS); 

      // THETAP  58.7811 58.7811
      // PHIP    23.3532 23.3532
      // declination: 2.714
      
      CPPUNIT_ASSERT(Verify<CloseTo>(zenith/deg, 58.7811));
      CPPUNIT_ASSERT(Verify<CloseTo>(azimuth/deg+360, 23.3532 + 180 - 2.714 + 90.0)); // the 2.714 degrees are the magnetic declination
      CPPUNIT_ASSERT(Verify<Equal>(shower.GetShowerNumber(), nEventsSeen));
      CPPUNIT_ASSERT(Verify<Equal>(shower.GetShowerRunId(), std::string("1")));
      CPPUNIT_ASSERT(Verify<Equal>(shower.HasGroundParticles(), false)); // no particle file
      CPPUNIT_ASSERT(Verify<Equal>(shower.HasGHParameters(), true));
      CPPUNIT_ASSERT(Verify<CloseTo>(shower.GetGHParameters().GetXMax()/g*cm2, 6.7450E+02));
      CPPUNIT_ASSERT(Verify<CloseTo>(shower.GetGHParameters().GetNMax(), 8.3636E+08));
      CPPUNIT_ASSERT(Verify<CloseTo>(shower.GetDistanceOfShowerMaximum()/cm, 123456.7));

      // test if REAS contents are set up correctly
      CPPUNIT_ASSERT(Verify<Equal>(shower.HasRadioSimulation(), true));
      evt::RadioSimulation& radsim = shower.GetRadioSimulation();
      CPPUNIT_ASSERT(Verify<Equal>(radsim.GetNumPulses(), (long)14));
      CPPUNIT_ASSERT(Verify<Equal>(radsim.HasCorePosition(), true));
      bool ok;
      const evt::SimRadioPulse& simpulse = radsim.GetNextSimRadioPulse(ok);
      const unsigned int frontpadding = 0;  // samples used for padding of the timeseries front
      const double binning = 1.4*nanosecond;
      CPPUNIT_ASSERT(Verify<Equal>(ok, true));
      CPPUNIT_ASSERT(Verify<Equal>(simpulse.GetEfieldTimeSeries().GetSize(), (TraceV3D::SizeType)2012+frontpadding));
      CPPUNIT_ASSERT(Verify<CloseTo>(simpulse.GetEfieldTimeSeries().GetBinning(), binning));
      CPPUNIT_ASSERT(Verify<CloseTo>(simpulse.GetBinning(), binning));
      CPPUNIT_ASSERT(Verify<CloseTo>(simpulse.GetStartTime(), 5.026e-07*second-frontpadding*binning));

      //const double conversionfactor = 2.99792458e10; // conversion factor from cgs to microVolt/meter
      //const unsigned int samplenumber = 10; // test sample number 10
      //Vector3D sample = simpulse.GetEfieldTimeSeries()[samplenumber+frontpadding-1];

      // cannot easily test amplitudes any more because of electric field rotation due to magnetic declination
      //CPPUNIT_ASSERT (Verify<CloseTo>(sample[0]/(micro*volt)*meter, 5.7855425609404e-12*conversionfactor)); // east component (-column 2 from REAS)
      //CPPUNIT_ASSERT (Verify<CloseTo>(sample[1]/(micro*volt)*meter, -3.375375296975e-11*conversionfactor)); // north component (column 1 from REAS)
      //CPPUNIT_ASSERT (Verify<CloseTo>(sample[2]/(micro*volt)*meter, -5.5331174800789e-11*conversionfactor));  // vertical up component (column 3 from REAS)

    }

    CPPUNIT_ASSERT(Verify<Equal>(nEventsInFile, nEventsSeen));
  }

};


CPPUNIT_TEST_SUITE_REGISTRATION(testREASFile);


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