AiresShowerFileParticleIterator.cc

Go to the documentation of this file.

#include <io/AiresShowerFileParticleIterator.h>

#include <io/AiresShowerFile.h>
#include <io/AiresUtilities.h>
#include <io/AiresWrapper.h>
using namespace io;

#include <evt/ShowerSimData.h>
using namespace evt;

#include <utl/AugerUnits.h>
#include <utl/AugerException.h>
#include <utl/CoordinateSystem.h>
#include <utl/ErrorLogger.h>
#include <utl/MathConstants.h>
#include <utl/Point.h>
#include <utl/Vector.h>
using namespace utl;

#include <boost/tuple/tuple.hpp>
#include <boost/tuple/tuple_io.hpp>

#include <cmath>
#include <iostream>
#include <sstream>

using namespace std;

AiresShowerFileParticleIterator::
AiresShowerFileParticleIterator(AiresShowerFile& theFile) :
  fFile(theFile),
  fCurrentParticle(utl::Particle::eUndefined,
                   utl::Particle::eShower,
                   utl::Point(0.0, 0.0, 0.0,
                              CoordinateSystem::GetRootCoordinateSystem()),
                   utl::Vector(0.0, 0.0, 0.0,
                               CoordinateSystem::GetRootCoordinateSystem()),
                   utl::TimeInterval(0.0),
                   0.0),
  fChannel(theFile.fChannel),
  fVerbosity(theFile.fVerbosity),
  fIrc(theFile.fIrc),
  fHeadOfList(-1), fEndOfList(-1),
  fAltType(theFile.fAltType),
  fGroundParticleCodeIndex(theFile.fGroundParticleCodeIndex),
  fLogRIndex(theFile.fLogRIndex),
  fLogEIndex(theFile.fLogEIndex),
  fThetaIndex(theFile.fThetaIndex),
  fUxIndex(theFile.fUxIndex),
  fUyIndex(theFile.fUyIndex),
  fTimeIndex(theFile.fTimeIndex),
  fWeightIndex(theFile.fWeightIndex),
  fCurrentRecordNumber(-1),
  fIteratorValid(false) {

  fHeadOfList = fFile.fCurrentShowerLocation->second.first + 1;
  fEndOfList = fFile.fCurrentShowerLocation->second.second;

  for (unsigned int i = 0; i < 99; ++i) {

    fIntData[i] = 0;
    fDoubleData[i] = 0.0;

  }

  io::Aires::kMagneticFieldDeclination = theFile.fGeomagneticFieldDeclination;

}

AiresShowerFileParticleIterator::~AiresShowerFileParticleIterator() {

}

Particle *AiresShowerFileParticleIterator::GetOneParticle(const utl::CoordinateSystemPtr& cs) {

  if (!fIteratorValid)
    throw IOFailureException("AiresShowerFileParticleIterator not valid.");

  if (fCurrentRecordNumber != fEndOfList) {

    int recordNumber = fCurrentRecordNumber - 1;
    AiresWrapper::crogotorec(&fChannel, &recordNumber, &fVerbosity, &fIrc);

    if (!AiresWrapper::getcrorecord(&fChannel, fIntData, fDoubleData,
                                    &fAltType, &fVerbosity, &fIrc)) {

      ERROR("Error attempting to retrieve record from AIRES file.");
      fIteratorValid = false;
      return NULL;

    }

    if (fIrc != 0 && fIrc != 1 && fIrc != 2) {

      ++fCurrentRecordNumber;
      return GetOneParticle(cs);

    }

    const double r = pow(10.0, fDoubleData[fLogRIndex])*m;
    const double theta = fDoubleData[fThetaIndex]*rad;
    const double z = 0.0;

    const double pX = fDoubleData[fUxIndex];
    const double pY = fDoubleData[fUyIndex];
    const double pZ2 = 1.0 - pX*pX - pY*pY;
    const double pZ = pZ2 > 0 ? -sqrt(pZ2) : 0.0;

    fCurrentParticle =
      Particle(Aires::AiresToPDG(int(fIntData[fGroundParticleCodeIndex])),
               utl::Particle::eShower,
               Point(r*cos(theta), r*sin(theta), z, cs),
               Vector(pX, pY, pZ, cs),
               utl::TimeInterval(fDoubleData[fTimeIndex]*ns),
               fDoubleData[fWeightIndex],
               pow(10.0, fDoubleData[fLogEIndex])*GeV);

    ++fCurrentRecordNumber;

    return &fCurrentParticle;

  }

  fIteratorValid = false;
  return NULL;

}

void AiresShowerFileParticleIterator::Rewind() {

  AiresWrapper::crogotorec(&fChannel, &fHeadOfList, &fVerbosity, &fIrc);

  fCurrentRecordNumber =
    AiresWrapper::crorecnumber(&fChannel, &fVerbosity, &fIrc);

  fIteratorValid = true;

}

CoordinateSystemPtr
AiresShowerFileParticleIterator::ComputeExternalShowerCoordinateSystem(
  const CoordinateSystemPtr& theCoordinateSystem)
{
  const double angle = io::Aires::AiresAzimuthToAuger(0.0*deg);

  return CoordinateSystem::RotationZ(angle, theCoordinateSystem);
}


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