SenecaShowerFile.cc

Go to the documentation of this file.

#include "SenecaShowerFile.h"

#include <evt/Event.h>
#include <evt/ShowerSimData.h>
#include <evt/DefaultShowerGeometryProducer.h>

#include <utl/AugerException.h>
#include <utl/ErrorLogger.h>
#include <utl/Math.h>

#include <string>
#include <fstream>
#include <iostream>

using namespace evt;
using namespace io;
using namespace utl;
using namespace std;


namespace io {

  SenecaShowerFile::SenecaShowerFile()
  {
    fSenecaParticleIterator.SetSenecaFileParser(&fFileParser);
  }


  SenecaShowerFile::SenecaShowerFile(const std::string& fileName, const io::Mode mode, utl::Branch* const b)
  {
    fSenecaParticleIterator.SetSenecaFileParser(&fFileParser);
    Open(fileName, mode, b);
  }


  void
  SenecaShowerFile::Open(const std::string& fileName, const io::Mode mode, utl::Branch* const /*b*/)
  {
    const io::Status check = fFileParser.OpenSenecaFile(fileName, mode);
    if (check != eSuccess) {
      const string msg = "Seneca IO Exception: Could not open " + fileName;
      throw utl::IOFailureException(msg);
    }
  }


  void
  SenecaShowerFile::Close()
  {
    fFileParser.CloseSenecaFile();
  }


  Status
  SenecaShowerFile::Read(evt::Event& event)
  {
    Status check = fFileParser.FileStatus();
    if (check != eSuccess)
      return check;

    if (!event.HasSimShower()) {
      const double rotationZ = utl::kPi;  // Seneca phi-angle is "incoming" not "outgoing"
      event.MakeSimShower(DefaultShowerGeometryProducer(rotationZ));
    } else {
      ERROR("ShowerSimData already exists. Cannot read seneca.");
      return eFail;
    }

    evt::ShowerSimData& shower = event.GetSimShower();
    check = fFileParser.GetNextEvent();
    if (check != eSuccess)
      return check;

    SenecaEventHeaderRecord* const eventHeader = fFileParser.GetCurrentEventHeader();
    if (!eventHeader)
      return eFail;

    shower.SetPrimaryParticle(SenecaParticleIDToAuger(eventHeader->ID0));
    shower.SetEnergy(eventHeader->E0 * GeV);
    shower.SetGroundParticleCoordinateSystemZenith(eventHeader->Theta);
    shower.SetGroundParticleCoordinateSystemAzimuth(eventHeader->Phi);

    ostringstream info;
    info << "SenecaShowerFile" << eventHeader->Theta/degree << ' ' << eventHeader->Phi/degree;
    INFO(info);
    shower.SetShowerNumber(fFileParser.GetCurrentEvent());
    shower.SetShowerRunId("1");

    if (!shower.HasGroundParticles())
      shower.MakeGroundParticles();

    if (CreateSenecaFileParticleIterator() != eSuccess)
      return eFail;

    shower.GetGroundParticles().SetFileInterface(&fSenecaParticleIterator);

    return eSuccess;
  }


  void
  SenecaShowerFile::Write(const evt::Event& /*event*/)
  {
    const string msg = "Seneca IO Exception: Cannot write to Seneca files.\n";
    throw utl::IOFailureException(msg);
  }


  Status
  SenecaShowerFile::FindEvent(const unsigned int eventId)
  {
    return fFileParser.GotoEvent(eventId);
  }


  Status
  SenecaShowerFile::GotoPosition(const unsigned int position)
  {
    return fFileParser.GotoPosition(position);
  }


  int
  SenecaShowerFile::GetNEvents()
  {
    return fFileParser.GetNumberEvents();
  }


  Status
  SenecaShowerFile::CreateSenecaFileParticleIterator()
  {
    return fSenecaParticleIterator.SetSenecaFileParser(&fFileParser);
  }


  int
  SenecaShowerFile::SenecaParticleIDToAuger(const int senecaID)
  {
    int particleID = utl::Particle::eUndefined;
    switch (senecaID) {
    case 10:
      particleID = utl::Particle::ePhoton;
      break;
    case 12:
      particleID = utl::Particle::eElectron;
      break;
    case -12:
      particleID = utl::Particle::ePositron;
      break;
    case 14:
      particleID = utl::Particle::eMuon;
      break;
    case -14:
      particleID = utl::Particle::eAntiMuon;
      break;
    case 110:
      particleID = utl::Particle::ePiZero;
      break;
    case 120:
      particleID = utl::Particle::ePiPlus;
      break;
    case -120:
      particleID = utl::Particle::ePiMinus;
      break;
    case 220:
      particleID = utl::Particle::eEta;
      break;
    case 130:
      particleID = utl::Particle::eKaonPlus;
      break;
    case -130:
      particleID = utl::Particle::eKaonMinus;
      break;
    case 230:
      particleID = utl::Particle::eUndefined;
      break;
    case -230:
      particleID = utl::Particle::eUndefined;
      break;
    case 330:
      particleID = utl::Particle::eUndefined;
      break;
    case 1120:
      particleID = utl::Particle::eProton;
      break;
    case 1220:
      particleID = utl::Particle::eNeutron;
      break;
    case 1130:
      particleID = utl::Particle::eUndefined;
      break;
    case 1230:
      particleID = utl::Particle::eUndefined;
      break;
    case 2130:
      particleID = utl::Particle::eLambda;
      break;
    }
    return particleID;
  }

}