AiresShowerFile.cc

Go to the documentation of this file.

#include <io/AiresShowerFile.h>
#include <io/AiresUtilities.h>
#include <io/AiresWrapper.h>
#include <io/AiresShowerFileParticleIterator.h>
#include <evt/Event.h>
#include <evt/DefaultShowerGeometryProducer.h>
#include <evt/ShowerSimData.h>
#include <evt/GaisserHillas2Parameter.h>
#include <evt/GaisserHillas4Parameter.h>
#include <evt/GaisserHillasTypes.h>
#include <utl/AugerUnits.h>
#include <utl/AugerCoordinateSystem.h>
#include <utl/CoordinateSystem.h>
#include <utl/ErrorLogger.h>
#include <utl/Point.h>
#include <utl/TabulatedFunction.h>
#include <utl/Vector.h>
#include <utl/Particle.h>
#include <boost/tuple/tuple.hpp>
#include <boost/tuple/tuple_io.hpp>
#include <boost/filesystem/operations.hpp>
#include <boost/filesystem/path.hpp>
#include <cmath>
#include <iostream>
#include <sstream>
#include <fstream>
#include <vector>

using namespace std;
using namespace io;
using namespace evt;
using namespace utl;
namespace fs = boost::filesystem;


AiresShowerFile::AiresShowerFile() :
  VEventFile()
{
  Init();
}


AiresShowerFile::AiresShowerFile(const string& theFileName,
                                 const Mode theMode,
                                 const utl::Branch* const branch) :
  VEventFile(theFilename, eRead)
{
  Init();
  Open(theFilename, eRead, branch);
}


AiresShowerFile::~AiresShowerFile()
{
  Close();
  delete fParticleIterator;
}


void
AiresShowerFile::Open(const string& theFileName, const Mode theMode, const utl::Branch* const branch)
{
  if (fIsOpen) {
    Close();
    Init();
  }

  // If trying to specify a mode other than read throw an exception
  if (theMode != eRead) {
    FATAL("AIRES data file can only be opened in read-only mode!");
    throw utl::IOFailureException();
  }

  DefaultOpen(theFilename, theMode);

  // We have to split TheFilename up into its path and actual filename since
  // AiresWrapper::opencrofilec takes both the directory and filename as
  // arguments.
#if BOOST_FILESYSTEM_VERSION < 3
  fs::path fullPath(fFilename.c_str(), fs::native);
  string theDirectory = fullPath.branch_path().native_file_string();
  string theFile = fullPath.leaf();
#else
  fs::path fullPath(fFilename.c_str());
  string theDirectory = fullPath.branch_path().c_str();
  string theFile = fullPath.leaf().c_str();
#endif

  // Initialise AIRES
  AiresWrapper::ciorinit(&fInitLevel, &fCodsys, &fVerbosity, &fIrc);

  // Throw the information to AIRES and see if we can open the file. AIRES
  // will fail if the directory and/or file do not exist
  AiresWrapper::opencrofilec(theDirectory.c_str(), theFile.c_str(),
                             &fSkipHeader, &fLogBase, &fVerbosity,
                             &fChannel, &fReturnCode);

  // Set the indices into the file for the various record information
  SetRecordIndices();

  // Set the data record information
  SetRecordInformation();

  // Read and fill the file header information - this is information
  // directly connected with the file and so should reside within what
  // effectively amounts to the in-memory file descriptor
  int intData[15], showerPrimaryCode;
  double doubleData[37], showerPrimaryWeight;

  fill(intData, intData + 15, 0);
  fill(doubleData, doubleData + 37, 0);

  AiresWrapper::croinputdata0(intData, doubleData,
                              &showerPrimaryCode, &showerPrimaryWeight);

  fNumberDiffPrimaries = intData[0];
  fPrimaryEnergyDist = intData[4];
  fZenithDist = intData[5];
  fAzimuthDist = intData[6];
  fNumObsLvls = intData[7];
  fAtmosphereModel = intData[8];
  fFirstShowerNumber = intData[14];
  fShowerPrimaryCode = showerPrimaryCode;

  fMinPrimaryEnergy = doubleData[0]*GeV;
  fMaxPrimaryEnergy = doubleData[1]*GeV;
  fEnergyExponent = doubleData[2];
  fMinZenith = doubleData[3]*deg;
  fMaxZenith = doubleData[4]*deg;
  fMinAzimuth = doubleData[5]*deg;
  fMaxAzimuth = doubleData[6]*deg;
  fThinningParameter = doubleData[7];
  fInjectionAltitude = doubleData[8]*m;
  fInjectionDepth = doubleData[9]*g/cm2;
  fGroundAltitude = doubleData[10]*m;
  fGroundDepth = doubleData[11]*g/cm2;
  fSiteLatitude = doubleData[19]*deg;
  fSiteLongitude = doubleData[20]*deg;
  fGeomagneticFieldStrength = doubleData[21]*(1.e-9*tesla);
  fGeomagneticFieldInclination = doubleData[22]*deg;
  fGeomagneticFieldDeclination = doubleData[23]*deg;
  fMinRadiusCut = doubleData[35]*meter;
  fMaxRadiusCut = doubleData[36]*meter;

  // Set the declination to ensure we get the transformation to the
  // Auger coordinate system correct
  io::Aires::kMagneticFieldDeclination = fGeomagneticFieldDeclination;

  fShowerPrimaryWeight = showerPrimaryWeight;

  fCurrentShowerLocation = fDataLocations.begin();

  fIsOpen = true;
}


void
AiresShowerFile::Close()
{
  // Get rid of the invalid particle iterator and close the file
  delete fParticleIterator;
  fParticleIterator = NULL;

  AiresWrapper::cioclose1(&fChannel);
  AiresWrapper::ciorshutdown();

  fIsOpen = false;
}


void
AiresShowerFile::Write(const evt::Event& /*event*/)
{
  // This should never get called
  ERROR("Attempt to write AIRES data file!");
  throw utl::IOFailureException();
}


Status
AiresShowerFile::Read(evt::Event& theEvent)
{
  if (fCurrentShowerLocation == fDataLocations.end())
    return eEOF;

  // crogotorec positions at the end of the record we request to go to.
  // So subtract one from the record we want.

  // We want to make sure we've got both the header and trailer record
  // set up before we even begin to make the shower and extract the
  // information. If we don't get both of these correct, return with failure.

  // Read and try to extract from header first.

  int headerRecordNumber = fCurrentShowerLocation->second.first - 1;

  if (!AiresWrapper::crogotorec(&fChannel, &headerRecordNumber,
                                &fVerbosity, &fIrc)) {

    ERROR("Error in positioning with AIRES data file.");
    return eFail;

  }

  int headerIntData[99];
  double headerDoubleData[99];

  fill(headerIntData, headerIntData + 99, 0);
  fill(headerDoubleData, headerDoubleData + 99, 0);

  if (!AiresWrapper::getcrorecord(&fChannel, headerIntData, headerDoubleData,
                                  &fAltType, &fVerbosity,
                                  &fReturnCode)) {

    ERROR("Error attempting to retrieve header for "
          "shower from AIRES data file.");

    return eFail;

  }

  // Trailer extraction comes next

  int trailerRecordNumber = fCurrentShowerLocation->second.second - 1;

  if (!AiresWrapper::crogotorec(&fChannel, &trailerRecordNumber,
                                &fVerbosity, &fIrc)) {
    ERROR("Error in positioning with AIRES data file.");
    return eFail;
  }

  int trailerIntData[99];
  double trailerDoubleData[99];

  fill(trailerIntData, trailerIntData + 99, 0);
  fill(trailerDoubleData, trailerDoubleData + 99, 0);

  if (!AiresWrapper::getcrorecord(&fChannel, trailerIntData, trailerDoubleData,
                                  &fAltType, &fVerbosity,
                                  &fReturnCode)) {

    ERROR("Error attempting to retrieve header for "
          "shower from AIRES data file.");

    return eFail;

  }

  // Okay since we got past those stages we can go ahead and start the
  // creation of all the in-memory junk we need to use the aires shower

  if (!theEvent.HasSimShower())
    theEvent.MakeSimShower(evt::DefaultShowerGeometryProducer());

  ShowerSimData& theShower = theEvent.GetSimShower();

  // Get the Aires task id.  This id is used to set the
  // shower run id portion of the event identifier
  //
  char taskNameChar[999];
  int taskNameCharLength;
  int taskNameVersion;
  char taskDate[20];  // size set by Aires convention
  AiresWrapper::crotaskidc(taskNameChar, &taskNameCharLength,
                           &taskNameVersion, taskDate);
  string taskName(taskNameChar, taskNameCharLength);
  theShower.SetShowerRunId(taskName);

  // Set the Min and Max radii cuts used to generate the shower
  //


  // Reading the longitudinal profile of the shower if available
  //
#if BOOST_FILESYSTEM_VERSION < 3
  fs::path fullPath(fFilename.c_str(), fs::native);
  const string theDirectory = fullPath.branch_path().native_file_string();
  const string theFile =
    fullPath.leaf().substr(0, fullPath.leaf().find(string(".")));
#else
  fs::path fullPath(fFilename.c_str());
  const string theDirectory = fullPath.branch_path().c_str();
  const string tmpPath = fullPath.leaf().c_str();
  const string theFile =
    tmpPath.substr(0, string(fullPath.leaf().c_str()).find(string(".")));
#endif

  stringstream filename;
  filename << theDirectory << '/' << theFile;

  // If ExportPerShower option is invoked during shower generation
  // , a _s appendix is included in the filename
  //
  int showerNumber = headerIntData[fShowerNumberIndex];
  stringstream showerNumberAppendixMultipleShowers;
  showerNumberAppendixMultipleShowers << "_s";
  if (showerNumber < 10000)
    showerNumberAppendixMultipleShowers.width(4);
  else
    showerNumberAppendixMultipleShowers.width(6);
  showerNumberAppendixMultipleShowers.fill('0');
  showerNumberAppendixMultipleShowers << showerNumber;

  // Change for Bug 352 - TAP 21/09/2005
  const string totalFilenameMultiple = filename.str() + showerNumberAppendixMultipleShowers.str();

  double depth, data, dummy;
  int index;
  TabulatedFunction profile;

  // Longitudinal charge profile

  profile.Clear();

  const double cosTheta = cos(headerDoubleData[fPrimaryZenithIndex]*deg);
  const double auxStep = fGroundDepth/(g/cm2)/fNumObsLvls/cosTheta;

  ifstream longitudinalDataFile;

  // Keep track of whether we are dealing with single or multiple
  // shower(s)/file.  The filenames (and formats) of the export
  // files are different in the two cases.
  //
  bool singleShower = true;

  // first attempt to open : case of multiple showers / file

  ostringstream infoOpeningLongProf;
  infoOpeningLongProf << "Reading profile from: ";

  if (fs::exists(totalFilenameMultiple + fChargeTableAppendix)) {
    longitudinalDataFile.open((totalFilenameMultiple + fChargeTableAppendix).c_str());
    singleShower = false;
    infoOpeningLongProf << totalFilenameMultiple;
  } else if (fs::exists(filename.str() + fChargeTableAppendix)) {
    longitudinalDataFile.open((filename.str() + fChargeTableAppendix).c_str());
    singleShower = true;
    infoOpeningLongProf << filename.str();
  }

  infoOpeningLongProf << fChargeTableAppendix;
  INFO(infoOpeningLongProf);

  if (longitudinalDataFile.is_open()) {

    char onechar;
    longitudinalDataFile.get(onechar);

    while (onechar) {

      if (onechar != '#') {

        longitudinalDataFile.unget();

        for (int i = 0; i < fNumObsLvls; ++i) {

          if (singleShower) // conform to different formats for single shower vs. multiple shower profiles
            longitudinalDataFile >> index >> depth >> data >> dummy >> dummy >> dummy >> dummy;
          else
            longitudinalDataFile >> index >> depth >> data;

          //BRD 17/2/05 Store profile in terms of slant depth. Here we
          //assume Aires tables are written with the default vertical
          //depth.  We apply a cosTheta correction here, really only
          //valid for a flat Earth.  If precision is required for
          //very inclined showers we suggest writing tables from Aires
          //with the slant depth option, as Aires does have a model
          //for a curved Earth.
          profile.PushBack(depth/cosTheta*g/cm2, data);
        }
        break;

      } else {
        char line[128];
        longitudinalDataFile.getline(line, 128);
        longitudinalDataFile.get(onechar);
      }
    }

    if (!theShower.HasLongitudinalProfile(ShowerSimData::eCharged))
      theShower.MakeLongitudinalProfile(profile, ShowerSimData::eCharged);

    // TAP - 09/02/2006 as per the resolution of bug 242 we do not attempt
    // to fit using aires: we obtain this information directly from the
    // aires trailer itself.

    // Fit to extract Gaisser-Hillas parameters from Aires profile

    double depth1[2000], nch[2000], weight[2000];

    if (index > 2000)
      index = 2000;

    for (int i = 0; i < index; ++i) {
      //BRD use real grams rather than Auger units for the Aires fitter
      depth1[i] = profile[i].X()/(g/cm2);
      nch[i] = profile[i].Y();
      weight[i] = 1.0;

    }

    int bodata = 1, eodata = index - 1, ws = 2;
    double minmax = 0.0, nminratio = 100.0;

    //output
    int bodataeff, eodataeff;
    double  nmax, xmax, x0, lambda, sq;
    int irc;

    AiresWrapper::fitghf(&bodata, &eodata, depth1, nch, weight,
                         &ws, &minmax, &nminratio, &bodataeff, &eodataeff,
                         &nmax, &xmax, &x0, &lambda, &sq, &irc);

    // End of fit

    // Adding GH information to simulated shower

    if (irc) {

      ostringstream message;
      message << "Aires fit for the Gaisser-Hillas function failed! \n"
              << "GHParameters object from simulated shower will \n"
              << "NOT be made.\n";
      WARNING(message);

    } else {

      GaisserHillas4Parameter gh(gh::eClassic);

      unsigned ndof = eodataeff - bodataeff - 4;  // 4-parameter fit I believe

      //BRD convert to Auger units before setting

      gh.SetXMax(xmax*(g/cm/cm), 0);
      gh.SetNMax(nmax, 0);
      gh.SetShapeParameter(gh::eX0, x0*(g/cm/cm), 0);
      gh.SetShapeParameter(gh::eLambda, lambda*(g/cm/cm), 0);
      gh.SetChiSquare(ndof*sq, ndof);

      if (!theShower.HasGHParameters())
          theShower.MakeGHParameters(gh);

    }

    longitudinalDataFile.close();

  } else {

    ostringstream message;
    message << "Unable to find Aires charge profile file ( neither "
            << totalFilenameMultiple << fChargeTableAppendix << " nor "
            << filename.str() << fChargeTableAppendix
            << ")";
    INFO(message);

  }

  // Longitudinal energy deposit profile

  profile.Clear();

  ostringstream infoOpeningEnergyDep;
  infoOpeningEnergyDep << "Reading energy deposit from: ";

  ifstream energyDepositDataFile;

  if (!singleShower) {
    energyDepositDataFile.open((totalFilenameMultiple + fdEdXTableAppendix).c_str());
    infoOpeningEnergyDep << totalFilenameMultiple;
  }
  else {
    energyDepositDataFile.open((filename.str() + fdEdXTableAppendix).c_str());
    infoOpeningEnergyDep << filename.str();
  }

  infoOpeningEnergyDep << fdEdXTableAppendix;
  INFO(infoOpeningEnergyDep);

  if (energyDepositDataFile.is_open()) {

    char onechar;
    energyDepositDataFile.get(onechar);

    while (onechar) {

      if (onechar != '#') {

        energyDepositDataFile.unget();

        for (int i = 0; i < fNumObsLvls; ++i) {

          if (singleShower)
            energyDepositDataFile >> index >> depth >> data >> dummy >> dummy >> dummy >> dummy;
          else
            energyDepositDataFile >> index >> depth >> data;

          //BRD 17/2/05 Store profile in terms of slant depth. Here we
          //assume Aires tables are written with the default vertical
          //depth.  We apply a cosTheta correction here, really only
          //valid for a flat Earth.  If precision is required for
          //very inclined showers we suggest writing tables from Aires
          //with the slant depth option, as Aires does have a model
          //for a curved Earth.
          profile.PushBack(depth/cosTheta*g/cm2,
                           data/auxStep*(GeV/(g/cm2)));

        }
        break;

      } else {

        char line[128];
        energyDepositDataFile.getline(line, 128);
        energyDepositDataFile.get(onechar);

      }
    }

    if (!theShower.HasdEdX())
      theShower.MakedEdX(profile);

    energyDepositDataFile.close();

  } else {
    ostringstream message;
    message << "Unable to find Aires energy deposit profile file (neither "
            << totalFilenameMultiple << fdEdXTableAppendix << " nor "
            << filename.str() << fdEdXTableAppendix
            << ")";
    INFO(message);
  }


  // --------------------------------------------------------------------
  // Longitudinal photon profile

  profile.Clear();

  ostringstream infoOpeningGammaProf;
  infoOpeningGammaProf << "Reading gamma profile from: ";

  ifstream gammaProfileDataFile;

  if (!singleShower) {
    gammaProfileDataFile.open((totalFilenameMultiple + fGammaTableAppendix).c_str());
    infoOpeningGammaProf << totalFilenameMultiple;
  }
  else {
    gammaProfileDataFile.open((filename.str() + fGammaTableAppendix).c_str());
    infoOpeningGammaProf << filename.str();
  }

  infoOpeningGammaProf << fGammaTableAppendix;
  INFO(infoOpeningGammaProf);

  if (gammaProfileDataFile.is_open()) {

    char onechar;
    gammaProfileDataFile.get(onechar);

    while (onechar) {

      if (onechar != '#') {

        gammaProfileDataFile.unget();

        for (int i = 0; i < fNumObsLvls; ++i) {

          if (singleShower)
            gammaProfileDataFile >> index >> depth >> data >> dummy >> dummy >> dummy >> dummy;
          else
            gammaProfileDataFile >> index >> depth >> data;

          profile.PushBack(depth/cosTheta*g/cm2,
                           data/auxStep*(GeV/(g/cm2)));

        }
        break;

      } else {

        char line[128];
        gammaProfileDataFile.getline(line, 128);
        gammaProfileDataFile.get(onechar);

      }
    }

    if (!theShower.HasLongitudinalProfile(ShowerSimData::ePhoton))
      theShower.MakeLongitudinalProfile(profile, ShowerSimData::ePhoton);

    gammaProfileDataFile.close();

  } else {
    // Unable to find Aires gamma profile file
  }


  // --------------------------------------------------------------------
  // Longitudinal electron profile

  profile.Clear();

  ostringstream infoOpeningElectronProf;
  infoOpeningElectronProf << "Reading electron profile from: ";

  ifstream electronProfileDataFile;

  if (!singleShower) {
    electronProfileDataFile.open((totalFilenameMultiple + fElectronTableAppendix).c_str());
    infoOpeningElectronProf << totalFilenameMultiple;
  }
  else {
    electronProfileDataFile.open((filename.str() + fElectronTableAppendix).c_str());
    infoOpeningElectronProf << filename.str();
  }

  infoOpeningElectronProf << fElectronTableAppendix;
  INFO(infoOpeningElectronProf);

  if (electronProfileDataFile.is_open()) {

    char onechar;
    electronProfileDataFile.get(onechar);

    while (onechar) {

      if (onechar != '#') {

        electronProfileDataFile.unget();

        for (int i = 0; i < fNumObsLvls; ++i) {

          if (singleShower)
            electronProfileDataFile >> index >> depth >> data >> dummy >> dummy >> dummy >> dummy;
          else
            electronProfileDataFile >> index >> depth >> data;

          profile.PushBack(depth/cosTheta*g/cm2,
                           data/auxStep*(GeV/(g/cm2)));

        }
        break;

      } else {

        char line[128];
        electronProfileDataFile.getline(line, 128);
        electronProfileDataFile.get(onechar);

      }
    }

    if (!theShower.HasLongitudinalProfile(ShowerSimData::eElectron))
      theShower.MakeLongitudinalProfile(profile, ShowerSimData::eElectron);

    electronProfileDataFile.close();

  } else {
    // Unable to find Aires electron profile file
  }


  // --------------------------------------------------------------------
  // Longitudinal muon profile

  profile.Clear();

  ostringstream infoOpeningMuonProf;
  infoOpeningMuonProf << "Reading muon profile from: ";

  ifstream muonProfileDataFile;

  if (!singleShower) {
    muonProfileDataFile.open((totalFilenameMultiple + fMuonTableAppendix).c_str());
    infoOpeningMuonProf << totalFilenameMultiple;
  }
  else {
    muonProfileDataFile.open((filename.str() + fMuonTableAppendix).c_str());
    infoOpeningMuonProf << filename.str();
  }

  infoOpeningMuonProf << fMuonTableAppendix;
  INFO(infoOpeningMuonProf);

  if (muonProfileDataFile.is_open()) {

    char onechar;
    muonProfileDataFile.get(onechar);

    while (onechar) {

      if (onechar != '#') {

        muonProfileDataFile.unget();

        for (int i = 0; i < fNumObsLvls; ++i) {

          if (singleShower)
            muonProfileDataFile >> index >> depth >> data >> dummy >> dummy >> dummy >> dummy;
          else
            muonProfileDataFile >> index >> depth >> data;

          profile.PushBack(depth/cosTheta*g/cm2,
                           data/auxStep*(GeV/(g/cm2)));

        }
        break;

      } else {

        char line[128];
        muonProfileDataFile.getline(line, 128);
        muonProfileDataFile.get(onechar);

      }
    }

    if (!theShower.HasLongitudinalProfile(ShowerSimData::eMuon))
      theShower.MakeLongitudinalProfile(profile, ShowerSimData::eMuon);

    muonProfileDataFile.close();

  } else {
    // Unable to find Aires muon profile file
  }

  const int code = headerIntData[fPrimaryCodeIndex];

  theShower.SetPrimaryParticle(io::Aires::AiresToPDG(code));
  theShower.SetEnergy(std::pow(10.0, headerDoubleData[fPrimaryLogEIndex])*GeV);
  theShower.SetXFirst(headerDoubleData[fFirstInteractionDepthIndex]*g/cm2);

  const double zenith = headerDoubleData[fPrimaryZenithIndex]*deg;
  const double azimuth = headerDoubleData[fPrimaryAzimuthIndex]*deg;

  theShower.SetGroundParticleCoordinateSystemAzimuth(azimuth);
  theShower.SetGroundParticleCoordinateSystemZenith(zenith);
  theShower.SetMinRadiusCut(fMinRadiusCut);
  theShower.SetMaxRadiusCut(fMaxRadiusCut);
  theShower.SetShowerNumber(headerIntData[fShowerNumberIndex]);

  // This stuff comes from the trailer.

  const int longFitReturnCode = trailerIntData[fLongFitReturnCodeIndex];
  const double xMax = trailerDoubleData[fShowerMaxIndex];
  const double nMax = trailerDoubleData[fTotalChargedParticlesIndex];

  if (!theShower.HasGHParameters() && !longFitReturnCode) {

    // longitudinal fit converged
    evt::GaisserHillas2Parameter gh;
    gh.SetXMax(xMax * (g/cm2), 0.);
    gh.SetNMax(nMax, 0.);
    theShower.MakeGHParameters(gh);

  }

  // Set up the information related to the ParticleList for this event
  delete fParticleIterator;

  fParticleIterator = new AiresShowerFileParticleIterator(*this);

  if (fMinRadiusCut <= 1.001*meter) {
    ostringstream info;
    info << "Inner radius cut r = " << fMinRadiusCut/meter
         << " m is small [max: 1 meter]:"
            " computing total number of muons at ground level";
    INFO(info);

    // make throw-away ShowerFileParticleIterator
    AiresShowerFileParticleIterator iter(*this);
    iter.Rewind();

    double muonWeight = 0;
    while (true) {
      Particle* const p = iter.GetOneParticle(CoordinateSystem::GetRootCoordinateSystem());
      if (!p)
        break;
      if (abs(p->GetType()) == Particle::eMuon)
        muonWeight += p->GetWeight();
    }
    theShower.SetMuonNumber(muonWeight);
  }

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

  theShower.GetGroundParticles().SetFileInterface(fParticleIterator);

  ++fCurrentShowerLocation;

  return eSuccess;
}


Status
AiresShowerFile::FindEvent(const unsigned int eventId)
{
  for (DataLocationsType::iterator dataIt = fDataLocations.begin();
       dataIt != fDataLocations.end(); ++dataIt)
    if (dataIt->first == int(eventId)) {
      fCurrentShowerLocation = dataIt;
      return eSuccess;
    }

  return eFail;
}


Status
AiresShowerFile::GotoPosition(const unsigned int position)
{
  if (position < fDataLocations.size()) {

    unsigned int i = 0;

    for (DataLocationsType::iterator dataIt = fDataLocations.begin();
         dataIt != fDataLocations.end(); ++dataIt, ++i)
      if (i == position) {
        fCurrentShowerLocation = dataIt;
        return eSuccess;
      }
  }

  return eFail;
}


int
AiresShowerFile::GetNEvents()
{
  return fDataLocations.size(); //fNumRecords[1];
}


void
AiresShowerFile::Init()
{
  fIsOpen = false;
  fParticleIterator = NULL;
  fChannel = 0;
  fVerbosity = 0;
  fReturnCode = -1;
  fInitLevel = 0;
  fCodsys = 0;
  fIrc = 0;
  fSkipHeader = 0;
  fLogBase = 10;
  fAltType = false;
  fGroundParticleCodeIndex = 0;
  fLogRIndex = 0;
  fLogEIndex = 0;
  fThetaIndex = 0;
  fUxIndex = 0;
  fUyIndex = 0;
  fTimeIndex = 0;
  fWeightIndex = 0;
  fPrimaryCodeIndex = 0;
  fShowerNumberIndex = 0;
  fPrimaryLogEIndex = 0;
  fPrimaryZenithIndex = 0;
  fPrimaryAzimuthIndex = 0;
  fPrimaryThinningIndex = 0;
  fFirstInteractionDepthIndex = 0;
  fInjectionAltitudeIndex = 0;
  fGlobalTimeShiftIndex = 0;
  fTotalParticlesInShowerIndex = 0;
  fNumberGroundParticlesIndex = 0;
  fLongFitReturnCodeIndex = 0;
  fShowerMaxIndex = 0;
  fTotalChargedParticlesIndex = 0;
  fNumberDiffPrimaries = 0;
  fPrimaryEnergyDist = 0;
  fZenithDist = 0;
  fAzimuthDist = 0;
  fNumObsLvls = 0;
  fAtmosphereModel = 0;
  fFirstShowerNumber = 0;
  fShowerPrimaryCode = 0;
  fMinPrimaryEnergy = 0;
  fMaxPrimaryEnergy = 0;
  fEnergyExponent = 0;
  fMinZenith = 0;
  fMaxZenith = 0;
  fMinAzimuth = 0;
  fMaxAzimuth = 0;
  fThinningParameter = 0;
  fInjectionAltitude = 0;
  fInjectionDepth = 0;
  fGroundAltitude = 0;
  fGroundDepth = 0;
  fSiteLatitude = 0;
  fSiteLongitude = 0;
  fShowerPrimaryWeight = 0;
  fGeomagneticFieldStrength = 0;
  fGeomagneticFieldInclination = 0;
  fGeomagneticFieldDeclination = 0;
  fMinRadiusCut = 0.;
  fMaxRadiusCut = 0.;
  fCurrentShowerNumber = -1;
  fChargeTableAppendix = ".t1291";
  fdEdXTableAppendix = ".t7991";
  fGammaTableAppendix = ".t1001";
  fElectronTableAppendix = ".t1205";
  fMuonTableAppendix = ".t1207";

  for (int i = 0; i < 5; ++i)
    fNumRecords[i] = 0;
}


void
AiresShowerFile::SetRecordIndices()
{
  int datatype;

  int rectype = 0;

  fGroundParticleCodeIndex =
    AiresWrapper::crofieldindex(&fChannel, &rectype, "Particle code",
                                &fVerbosity, &datatype, &fReturnCode) - 1;
  fLogRIndex =
    AiresWrapper::crofieldindex(&fChannel, &rectype, "Distance from the core",
                                &fVerbosity, &datatype, &fReturnCode) - 1;
  fLogEIndex =
    AiresWrapper::crofieldindex(&fChannel, &rectype, "Energy",
                                &fVerbosity, &datatype, &fReturnCode) - 1;
  fThetaIndex =
    AiresWrapper::crofieldindex(&fChannel, &rectype, "Ground plane polar angle",
                                &fVerbosity, &datatype, &fReturnCode) - 1;
  fUxIndex =
    AiresWrapper::crofieldindex(&fChannel, &rectype, "Direction of motion (x component)",
                                &fVerbosity, &datatype, &fReturnCode) - 1;
  fUyIndex =
    AiresWrapper::crofieldindex(&fChannel, &rectype, "Direction of motion (y component)",
                                &fVerbosity, &datatype, &fReturnCode) - 1;
  fTimeIndex =
    AiresWrapper::crofieldindex(&fChannel, &rectype, "Arrival time delay",
                                &fVerbosity, &datatype, &fReturnCode) - 1;
  fWeightIndex =
    AiresWrapper::crofieldindex(&fChannel, &rectype, "Particle weight",
                                &fVerbosity, &datatype, &fReturnCode) - 1;

  rectype = 1;

  fPrimaryCodeIndex =
    AiresWrapper::crofieldindex(&fChannel, &rectype, "Primary particle code",
                                &fVerbosity, &datatype, &fReturnCode) - 1;
  fShowerNumberIndex =
    AiresWrapper::crofieldindex(&fChannel, &rectype, "Shower number",
                                &fVerbosity, &datatype, &fReturnCode) - 1;
  fPrimaryLogEIndex =
    AiresWrapper::crofieldindex(&fChannel, &rectype, "Primary energy",
                                &fVerbosity, &datatype, &fReturnCode) - 1;
  fPrimaryZenithIndex =
    AiresWrapper::crofieldindex(&fChannel, &rectype, "Primary zenith angle",
                                &fVerbosity, &datatype, &fReturnCode) - 1;
  fPrimaryAzimuthIndex =
    AiresWrapper::crofieldindex(&fChannel, &rectype, "Primary azimuth angle",
                                &fVerbosity, &datatype, &fReturnCode) - 1;
  fPrimaryThinningIndex =
    AiresWrapper::crofieldindex(&fChannel, &rectype, "Thinning energy",
                                &fVerbosity, &datatype, &fReturnCode) - 1;
  fFirstInteractionDepthIndex =
    AiresWrapper::crofieldindex(&fChannel, &rectype, "First interaction depth",
                                &fVerbosity, &datatype, &fReturnCode) - 1;
  fInjectionAltitudeIndex =
    AiresWrapper::crofieldindex(&fChannel, &rectype, "Central injection altitude",
                                &fVerbosity, &datatype, &fReturnCode) - 1;
  fGlobalTimeShiftIndex =
    AiresWrapper::crofieldindex(&fChannel, &rectype, "Global time shift",
                                &fVerbosity, &datatype, &fReturnCode) - 1;

  rectype = 2;

  fTotalParticlesInShowerIndex =
    AiresWrapper::crofieldindex(&fChannel, &rectype, "Total number of shower particles",
                                &fVerbosity, &datatype, &fReturnCode) - 1;
  fNumberGroundParticlesIndex =
    AiresWrapper::crofieldindex(&fChannel, &rectype, "Number of particles reaching ground",
                                &fVerbosity, &datatype, &fReturnCode) - 1;
  fLongFitReturnCodeIndex =
    AiresWrapper::crofieldindex(&fChannel, &rectype, "Xmax fit return code",
                                &fVerbosity, &datatype, &fReturnCode) - 1;
  fShowerMaxIndex =
    AiresWrapper::crofieldindex(&fChannel, &rectype, "Shower maximum depth",
                                &fVerbosity, &datatype, &fReturnCode) - 1;
  fTotalChargedParticlesIndex =
    AiresWrapper::crofieldindex(&fChannel, &rectype, "Total charged particles at shower maximum",
                                &fVerbosity, &datatype, &fReturnCode) - 1;
}


void
AiresShowerFile::SetRecordInformation()
{
  int nrtype;

  AiresWrapper::croreccount(&fChannel, &fVerbosity, &nrtype,
                            fNumRecords, &fReturnCode);

  AiresWrapper::crorewind(&fChannel, &fVerbosity, &fReturnCode);

  fDataLocations.clear();

  // fNumRecords contains the number of different record types within
  // the file. fNumRecords[1] constains the number of shower header
  // records (type 1 records)

  int intype1 = 1;
  int intype2 = 2;
  int infield1;
  int currentRecordNumber;

  int intData[99];
  double doubleData[99];

  fill(intData, intData + 99, 0);
  fill(doubleData, doubleData + 99, 0);

  while (AiresWrapper::crorecfind(&fChannel, &intype1, &fVerbosity, &infield1,
                                  &fIrc)) {

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

     if (!AiresWrapper::regetcrorecord(&fChannel, intData, doubleData,
                                    &fAltType, &fVerbosity,
                                    &fReturnCode)) {

      ERROR("Invalid read from AIRES data file.");
      fDataLocations.clear();
      break;
    }

    fDataLocations[intData[fShowerNumberIndex]].first = currentRecordNumber;

    // Now search for the trailer record
    AiresWrapper::crorecfind(&fChannel, &intype2, &fVerbosity,
                             &infield1, &fIrc);

    fDataLocations[intData[fShowerNumberIndex]].second =
      AiresWrapper::crorecnumber(&fChannel, &fVerbosity, &fIrc);
  }

  AiresWrapper::crorewind(&fChannel, &fVerbosity, &fReturnCode);
}