CONEXFile.cc

Go to the documentation of this file.

#include <io/CONEXFile.h>
#include <io/CONEXStructures.h>
#include <io/CONEXIOException.h>

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

#include <utl/ErrorLogger.h>
#include <utl/AugerUnits.h>
#include <utl/SaveCurrentTDirectory.h>

#include <cmath>
#include <sstream>
#include <string>
#include <sstream>
#include <iostream>
#include <iomanip>

#include <TFile.h>
#include <TTree.h>

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


namespace io {

  CONEXFile::~CONEXFile()
  {
    delete fFile;
    delete fShower;
    delete fHeader;
  }


  CONEXFile::CONEXFile(const std::string& fileName, const Mode mode, utl::Branch* const b) :
    VEventFile(fileName, mode)
  {
    Open(fileName, mode, b);
  }


  void
  CONEXFile::Close()
  {
    if (fFile) {
      if (fFile->IsOpen())
        fFile->Close();
      delete fFile;
      fFile = nullptr;
      delete fShower;
      fShower = nullptr;
      delete fHeader;
      fHeader = nullptr;
    }
  }


  void
  CONEXFile::Open(const std::string& fileName, const Mode mode, utl::Branch* const /*b*/)
  {
    if (fFile)
      Close();

    if (mode != eRead) {
      ostringstream msg;
      msg << "Cannot open file " << fileName << " - CONEX files are read-only";
      ERROR(msg);
      throw CONEXIOException(msg.str());
    }

    const SaveCurrentTDirectory save;
    fFile = TFile::Open(fileName.c_str());

    if (!fFile) {
      ostringstream msg;
      msg << "Cannot open file " << fileName << " - ERROR";
      ERROR(msg);
      throw CONEXIOException(msg.str());
    }

    if (!fFile->IsOpen()) {
      ostringstream msg;
      msg << "Cannot open file " << fileName << " - ERROR";
      ERROR(msg);
      throw CONEXIOException(msg.str());
    }

    fHeaderTree = (TTree*)fFile->Get("Header");
    fShowerTree = (TTree*)fFile->Get("Shower");

    if (!fHeaderTree || !fShowerTree) {
      ostringstream msg;
      msg << "Cannot open trees in file " << fileName << " - ERROR";
      ERROR(msg);
      throw CONEXIOException(msg.str());
    }

    // check CONEX version
    float version = 0;
    fHeaderTree->SetBranchAddress("Version", &version);
    fHeaderTree->GetEntry(0);
    fHeaderTree->ResetBranchAddresses();

    // check CXROOT version (more important)
    fROOTVersion = 0;
    TBranch* rVbranch = fHeaderTree->GetBranch("OutputVersion");
    if (rVbranch) {
      rVbranch->SetAddress(&fROOTVersion);
      rVbranch->GetEntry(0);
    }
    fHeaderTree->ResetBranchAddresses();

    ostringstream infoVersion;
    infoVersion << " CONEX version " << version << " Output version " << fROOTVersion;
    INFO(infoVersion);

    // read number of slant depth slices
    if (version > 0 && fROOTVersion < 2.0) {
      const int nX = int(fHeaderTree->GetMaximum("nX"));
      if (fHeaderTree->GetMaximum("delX") != fHeaderTree->GetMinimum("delX") ) {
        const string msg = "Cannot read from merged CONEXFile with different grammage binnings !!!!! Check your CONEX input.";
        ERROR(msg);
        throw CONEXIOException(msg);
      }
      float grammage[nX];
      fHeaderTree->SetBranchAddress("X", grammage);
      for (int iHeader = 0; iHeader < fHeaderTree->GetEntries(); ++iHeader) {
        fHeaderTree->GetEntry(iHeader);
        if (!iHeader) {
          fFirstGrammageBin = grammage[0];
        } else {
          if (grammage[0] != fFirstGrammageBin) {
            const string msg =
              "Cannot read from merged CONEXFile with incompatible shower content !!!!! "
              "Check your CONEX merging or don't use merged CONEX files !!!";
            ERROR(msg);
            throw CONEXIOException(msg);
          }
        }
      }

      fHeaderTree->ResetBranchAddresses();
      fShower = new CONEXShower(nX);
      MapShower(fROOTVersion);
    }

    fHeader = new CONEXHeader();
    fHeader->fFirstDepth = fFirstGrammageBin;
    MapHeader(fROOTVersion);
    fHeaderTree->GetEntry(0);
  }


  Status
  CONEXFile::Read(evt::Event& event)
  {
    if (!fFile)
      return eFail;

    if (event.HasSimShower()) {
      ERROR("Event not cleared - has SimShower. Cannot read CONEX.");
      return eFail;
    }
    event.MakeSimShower(DefaultShowerGeometryProducer());

    if (!fFile->IsOpen() || fCurrentPosition >= (unsigned int)GetNEvents())
      return eEOF;

    const int entry = NextEntry();

    if (fROOTVersion >= 2.0) {
      int nX_shower = 0;
      fShowerTree->ResetBranchAddresses();
      fShowerTree->SetBranchAddress("nX", &nX_shower);
      fShowerTree->GetEntry(entry);
      fShowerTree->ResetBranchAddresses();
      if (!fShower)
        fShower = new CONEXShower(0);
      fShower->SetVectorLength(nX_shower);
      MapShower(fROOTVersion);
    }

    fShowerTree->GetEntry(entry);

    ShowerSimData& shower = event.GetSimShower();
    FillShowerSimDataFromConex(*fHeader, *fShower, shower);
    FillShowerProfileDataFromConex(*fHeader, *fShower, fLPvector, fROOTVersion, shower);


    shower.SetShowerNumber(entry);
    ostringstream runIdentifier;
    runIdentifier << "cx_" << fHeader->fVersion
                  << "_" << fHeader->fHighEnergyModel
                  << "_" << fHeader->fSeed;
    shower.SetShowerRunId(runIdentifier.str());

    ostringstream info;
    info << "\n"
            "\t ----------- CONEX -------------\n"
            "\t event:      " << entry << "\n"
            "\t xmax/g*cm2: " << shower.GetGHParameters().GetXMax()/g*cm2 << "\n"
            "\t energy/eV:  " << shower.GetEnergy()/eV << "\n"
            "\t zenith/deg: " << fShower->fZenith << "\n"
            "\t ----------- CONEX -------------";

    INFO(info);


    return eSuccess;
  }


  void
  CONEXFile::Write(const evt::Event& /*event*/)
  {
    const string msg = "Cannot write to CONEXFile - read-only file";
    ERROR(msg);
    throw CONEXIOException(msg);
  }


  Status
  CONEXFile::FindEvent(const unsigned int eventId)
  {
    // CONEX shower don't have event ids
    return GotoPosition(eventId);
  }


  Status
  CONEXFile::GotoPosition(const unsigned int position)
  {
    if (position > (unsigned int)GetNEvents())
      return eFail;
    else {
      fCurrentPosition = position;
      return eSuccess;
    }
  }


  int
  CONEXFile::NextEntry()
  {
    ++fCurrentPosition;
    return fCurrentPosition - 1;
  }


  int
  CONEXFile::GetNEvents()
  {
    if (fShowerTree)
      return fShowerTree->GetEntries();
    else {
      const string msg = "Cannot request number of events from NULL tree";
      ERROR(msg);
      throw CONEXIOException(msg);
    }
  }


  void
  CONEXFile::MapHeader(const float /*version*/)
  {
    // no version checks yet

    fHeaderTree->SetBranchAddress("Seed1", &(fHeader->fSeed));
    fHeaderTree->SetBranchAddress("HEModel", &(fHeader->fHighEnergyModel));

    fHeaderTree->SetBranchAddress("Particle", &(fHeader->fParticleID));
    fHeaderTree->SetBranchAddress("hadCut", &(fHeader->fECutHadrons));
    fHeaderTree->SetBranchAddress("emCut", &(fHeader->fECutEM));
    fHeaderTree->SetBranchAddress("hadThr", &(fHeader->fEThrHadrons));
    fHeaderTree->SetBranchAddress("muThr", &(fHeader->fEThrMuons));
    fHeaderTree->SetBranchAddress("emThr", &(fHeader->fEThrEM));
    if (fROOTVersion < 2.0)
      fHeaderTree->SetBranchAddress("delX", &(fHeader->fDelX));
    else if (fROOTVersion >= 2.5) {
      fHeaderTree->SetBranchAddress("LEModel", &(fHeader->fLowEnergyModel));
      fHeaderTree->SetBranchAddress("HiLowEgy", &(fHeader->fHighLowEnergy));
      fHeaderTree->SetBranchAddress("haCut", &(fHeader->fECutLongProfHadrons));
      fHeaderTree->SetBranchAddress("muCut", &(fHeader->fECutLongProfMuons));
      fHeaderTree->SetBranchAddress("elCut", &(fHeader->fECutLongProfElectrons));
      fHeaderTree->SetBranchAddress("gaCut", &(fHeader->fECutLongProfPhotons));
    }
  }


  void
  CONEXFile::MapShower(const float /*version*/)
  {
    // no version checks yet

    fShowerTree->SetBranchAddress("lgE", &(fShower->fLogE));
    fShowerTree->SetBranchAddress("zenith", &(fShower->fZenith));
    fShowerTree->SetBranchAddress("Xfirst", &(fShower->fXfirst));
    fShowerTree->SetBranchAddress("X0", &(fShower->fX0));
    fShowerTree->SetBranchAddress("Xmax", &(fShower->fXmax));
    fShowerTree->SetBranchAddress("Nmax", &(fShower->fNmax));
    fShowerTree->SetBranchAddress("p1", &(fShower->fP1));
    fShowerTree->SetBranchAddress("p2", &(fShower->fP2));
    fShowerTree->SetBranchAddress("p3", &(fShower->fP3));
    fShowerTree->SetBranchAddress("chi2", &(fShower->fChi2));
    fShowerTree->SetBranchAddress("nX", &(fShower->fNX));
    fShowerTree->SetBranchAddress("N", fShower->fN);
    fShowerTree->SetBranchAddress("dEdX", fShower->fdEdX);
    fShowerTree->SetBranchAddress("Mu", fShower->fNmu);

    if (fROOTVersion >= 1.1) {
      fShowerTree->SetBranchAddress("Gamma", fShower->fNgam);
      fShowerTree->SetBranchAddress("Electrons", fShower->fNele);
      fShowerTree->SetBranchAddress("dMu", fShower->fMuProd);
    }

    fShowerTree->SetBranchAddress("EGround", fShower->fGroundEnergy);

    if (fROOTVersion >= 1.2)
      fShowerTree->SetBranchAddress("azimuth", &(fShower->fAzimuth));

    if (fROOTVersion >= 2.0)
      fShowerTree->SetBranchAddress("X", fShower->fX);
  }

}