testCorsikaShowerFile.cc

Go to the documentation of this file.

#include <io/CorsikaShowerFile.h>
#include <io/CorsikaBlock.h>
#include <io/CorsikaIOException.h>
#include <io/IoCodes.h>
#include <evt/Event.h>
#include <evt/ShowerSimData.h>
#include <utl/AugerUnits.h>
#include <utl/Particle.h>
#include <utl/ShowerParticleIterator.h>
#include <utl/UTCDateTime.h>

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

#include <det/Detector.h>

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

#include <iostream>
#include <iomanip>
#include <string>

using io::Corsika::Block;
using io::CorsikaShowerFile;
using tst::Expected;
using namespace fwk;
using namespace utl;
using namespace tst;

using std::endl;
using std::cerr;

class testCorsikaShowerFile : public CppUnit::TestFixture {

  CPPUNIT_TEST_SUITE(testCorsikaShowerFile);
  CPPUNIT_TEST(testConstruct);
  CPPUNIT_TEST_EXCEPTION(testConstructException, io::CorsikaIOException);
  CPPUNIT_TEST(testOpen);
  CPPUNIT_TEST_EXCEPTION(testOpenOpen, io::CorsikaIOException);
  CPPUNIT_TEST_EXCEPTION(testOpenNonExisting, io::CorsikaIOException);
  CPPUNIT_TEST_EXCEPTION(testOpenForWrite, io::CorsikaIOException);
  CPPUNIT_TEST(testClose);
  CPPUNIT_TEST(testCloseClose);
  CPPUNIT_TEST_EXCEPTION(testWrite, io::CorsikaIOException);
  CPPUNIT_TEST(testFindEvent);
  CPPUNIT_TEST(testGotoPosition);
  CPPUNIT_TEST(testGetNEvents);
  CPPUNIT_TEST(testRead);
  CPPUNIT_TEST(testParticleRead);
  CPPUNIT_TEST(testCherenkovShower);
  CPPUNIT_TEST_EXCEPTION(testGetNEventsClosedException,
                         io::CorsikaIOException);
  CPPUNIT_TEST_SUITE_END();

public:
  // shared data for tests

  CorsikaShowerFile* f;
  evt::Event event;

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

  void tearDown()
  {
    delete f;
  }

  void testConstruct()
  {
    CorsikaShowerFile* file;
    CentralConfig::GetInstance(BOOTSTRAPFILE);
    det::Detector& theDet = det::Detector::GetInstance();
    theDet.Update(utl::UTCDateTime(2000,4,1,1).GetTimeStamp()); // any reasonable time will do for the magnetic field model
    CPPUNIT_ASSERT(file = new CorsikaShowerFile(TESTFILE));
    delete file;
  }

  void testConstructException()
  {
    Expected();
    CorsikaShowerFile f(TESTFILE, io::eWrite);
  }

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

  void testOpenNonExisting()
  {
    Expected();
    f->Open("NoSuchCorsikaShowerFile.part");
  }

  void testOpenOpen()
  {
    f->Open(TESTFILE);
    Expected();
    f->Open(TESTFILE);
  }

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

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

  void testCloseClose()
  {
    f->Close();
    Expected();
    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() == 2);
  }

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

  void testRead()
  {
    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));

      CPPUNIT_ASSERT(Verify<Equal>(shower.GetPrimaryParticle(),
                                 int(Particle::eProton)));
      CPPUNIT_ASSERT(Verify<CloseTo>(shower.GetEnergy()/eV, 1.e20));
      
      const utl::CoordinateSystemPtr localCS = shower.GetLocalCoordinateSystem();
      const double zenith = (-shower.GetDirection()).GetTheta(localCS);
      const double azimuth = (-shower.GetDirection()).GetPhi(localCS); 

      // Test shower has: Zenith=32.42 deg, Azimuth=325.59 deg,
      // Note: azimuth is already corrected by 180deg (incoming vs. outgoing)
      // AARANG: 0

      CPPUNIT_ASSERT(Verify<CloseTo>(zenith/deg, 32.42));
      CPPUNIT_ASSERT(Verify<CloseTo>(azimuth/deg, 325.59+90.-4.233-360));
      CPPUNIT_ASSERT(Verify<CloseTo>(shower.GetMinRadiusCut(), // set in THIN parameter of CorsikaTest.inp
                                   101.*cm));
      CPPUNIT_ASSERT(Verify<Equal>(shower.GetShowerNumber(),
                                 nEventsSeen));
      CPPUNIT_ASSERT(Verify<Equal>(shower.GetShowerRunId(), std::string("1")));

      event = evt::Event(); // clear
    }

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

  void testParticleRead()
  {

    f->Open(TESTFILE);
    for (int i = 0; i < f->GetNEvents(); ++i) {
      event=evt::Event();       // clear event until readers do it.
      f->Read(event);
      evt::ShowerSimData& shower = event.GetSimShower();

      shower.MakeGeometry(Point(0, 0, 0,
                                det::Detector::GetInstance().GetSiteCoordinateSystem()));

      long double nParticles = 0;
      for (ShowerParticleIterator p = shower.GroundParticlesBegin();
           p != shower.GroundParticlesEnd();
           ++p) {
        nParticles += p->GetWeight();
      }
      CPPUNIT_ASSERT(nParticles > 1e11);
      CPPUNIT_ASSERT(nParticles < 1e12);
    }
  }

  void testCherenkovShower()
  {
    f->Open(TESTFILE2);
    int nEventsInFile = f->GetNEvents();
    int nEventsSeen = 0;

    f->SetReadLongFile(false);
    
    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));

      CPPUNIT_ASSERT(Verify<Equal>(shower.GetPrimaryParticle(),
                                 int(Particle::eProton)));
      CPPUNIT_ASSERT(Verify<CloseTo>(shower.GetEnergy()/GeV, 6952381.7620710116));
      
      const utl::CoordinateSystemPtr localCS = shower.GetLocalCoordinateSystem();
      const double zenith = (-shower.GetDirection()).GetTheta(localCS);
      const double azimuth = (-shower.GetDirection()).GetPhi(localCS); 

      // Zenith=25.2917 deg, Azimuth=239.904 deg
      // Note: azimuth is already corrected by 180deg (incoming vs. outgoing)
      // AARANG: -85.767 deg

      CPPUNIT_ASSERT(Verify<CloseTo>(zenith/deg, 25.291737580210896));
      // azimuth is corrected for: geomagnetic north direction, CORSIKA 180deg rotation of arrival direction
      CPPUNIT_ASSERT(Verify<CloseTo>(azimuth/deg, 239.904454153400025+90-4.233-360));
      CPPUNIT_ASSERT(Verify<Equal>(shower.GetShowerNumber(), nEventsSeen));

      event = evt::Event(); // clear
    }

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

};


CPPUNIT_TEST_SUITE_REGISTRATION(testCorsikaShowerFile);


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