cmpEvent.cc

Go to the documentation of this file.

#include "cmpEvent.h"

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

#include <TFile.h>
#include <TDirectory.h>

#include <adst/RecEventFile.h>
#include <adst/DetectorGeometry.h>
#include <adst/FileInfo.h>
#include <adst/RecEvent.h>
#include <adst/FDEvent.h>
#include <adst/GenShower.h>

#include "testlib.h"
#include "adstReading.h"
#include "cmpFDEvent.h"
#include "cmpSDEvent.h"

using namespace std;

void
cmpRecEventFile(RecEventFile& refFile, RecEventFile& testFile)
{
  cmpDetectorGeometry(refFile, testFile);
  cmpFileInfo(refFile, testFile);

  const unsigned int nRef  = refFile.GetNEvents();
  const unsigned int nTest = testFile.GetNEvents();
  isDiag(nTest, nRef, "Same no. of events");

  RecEvent* refRE  = new RecEvent();
  RecEvent* testRE = new RecEvent();
  ok(refFile.SetBuffers(&refRE) == RecEventFile::eSuccess, "reference file: SetBuffers == eSuccess");
  ok(testFile.SetBuffers(&testRE) == RecEventFile::eSuccess, "test file: SetBuffers == eSuccess");

  for (unsigned int iEvent = 0; iEvent < nRef && iEvent < nTest; ++iEvent)
  {
    ostringstream eventNoStr;
    eventNoStr << iEvent;
    string eventNo = string(" (event ") +  eventNoStr.str() + string(")");

    diag(string("Testing event ") + eventNoStr.str());

    ok(refFile.ReadEvent(iEvent) == RecEventFile::eSuccess, string("Read event from reference file") + eventNo);
    ok(testFile.ReadEvent(iEvent) == RecEventFile::eSuccess, string("Read event from test file") + eventNo);

    diag(string("This event has auger id ") + refRE->GetAugerId());

    cmpEvent(*refRE, *testRE);
  }

  delete refRE;
  delete testRE;
}

void
cmpFileInfo(const RecEventFile& refFile, const RecEventFile& testFile)
{
  FileInfo* refInfo  = readFileInfo(refFile, "reference");
  FileInfo* testInfo = readFileInfo(testFile, "test");
  if (!refInfo || !testInfo)
    return;


  isDiag(testInfo->HasMC(), refInfo->HasMC(), "FileInfo->HasMC");
  isDiag(testInfo->HasVEMTraces(), refInfo->HasVEMTraces(), "FileInfo->HasVEMTraces");
  isDiag(testInfo->HasSdAllTraces(), refInfo->HasSdAllTraces(), "FileInfo->HasSdAllTraces");
  isDiag(testInfo->HasFdTraces(), refInfo->HasFdTraces(), "FileInfo->HasFdTraces");
  isDiag(testInfo->HasFdSpotRecTraces(), refInfo->HasFdSpotRecTraces(), "FileInfo->HasFdSpotRecTraces");
  ok(testInfo->GetHost() != string(""), "FileInfo->GetHost");
  ok(testInfo->GetUser() != string(""), "FileInfo->GetUser");
  ok(testInfo->GetOfflineConfiguration() != string(""), "FileInfo->GetOfflineConfiguration");
  ok(testInfo->GetRecEventVersion() != string(""), "FileInfo->GetRecEventVersion");
  isDiag(testInfo->GetLDFType(), refInfo->GetLDFType(), "FileInfo->GetLDFType");
  isDiag(testInfo->GetLDFType(), refInfo->GetLDFType(), "FileInfo->GetLDFType");

  delete refInfo;
  delete testInfo;
}

void
cmpDetector(const RecEvent& ref, const RecEvent& test)
{
//   const TBits& GetMirrorsInDAQ(int iEye)
//   isDiag(test.GetDetector().\&(), ref.GetDetector().\&(), "detector.\&");

   isDiag(test.GetDetector().GetTemperature(),
          ref.GetDetector().GetTemperature(),
          "detector.GetTemperature");
   isDiag(test.GetDetector().GetPressure(),
          ref.GetDetector().GetPressure(),
          "detector.GetPressure");
  isDiag(test.GetDetector().GetTemperaturePressureGPSSeconds(),
         ref.GetDetector().GetTemperaturePressureGPSSeconds(),
         "detector.GetTemperaturePressureGPSSeconds");
  isDiag(test.GetDetector().HasWeatherStationData(),
 ref.GetDetector().HasWeatherStationData(),
         "detector.HasWeatherStationData");
  isDiag(test.GetDetector().HasMieDatabase(),
         ref.GetDetector().HasMieDatabase(),
         "detector.HasMieDatabase");
  isDiag(test.GetDetector().HasVAODs(),
         ref.GetDetector().HasVAODs(),
         "detector.HasVAODs");
  isDiag(test.GetDetector().HasOverallQualityDatabase(),
         ref.GetDetector().HasOverallQualityDatabase(),
         "detector.HasOverallQualityDatabase");
  isDiag(test.GetDetector().GetQDBHorizonalUniformity(),
         ref.GetDetector().GetQDBHorizonalUniformity(),
         "detector.GetQDBHorizonalUniformity");
  isDiag(test.GetDetector().GetQDBCloudCoverage(),
         ref.GetDetector().GetQDBCloudCoverage(),
         "detector.GetQDBCloudCoverage");
  isDiag(test.GetDetector().GetQDBMinCloudBaseHeight(),
         ref.GetDetector().GetQDBMinCloudBaseHeight(),
         "detector.GetQDBMinCloudBaseHeight");
  isDiag(test.GetDetector().GetQDBMinCloudBaseDepth(),
         ref.GetDetector().GetQDBMinCloudBaseDepth(),
         "detector.GetQDBMinCloudBaseDepth");
  isDiag(test.GetDetector().HasAerosolData(),
         ref.GetDetector().HasAerosolData(),
         "detector.HasAerosolData");
  isDiag(test.GetDetector().HasLidarData(),
         ref.GetDetector().HasLidarData(),
         "detector.HasLidarData");
  isEpsDiag(test.GetDetector().GetVAODsAtReferenceHeight().size(),
            ref.GetDetector().GetVAODsAtReferenceHeight().size(),
            1.e-9,
            "detector.GetVAODsAtReferenceHeight().size");

  const TBits refEyes = ref.GetDetector().GetActiveEyes();
  const TBits testEyes = test.GetDetector().GetActiveEyes();
  isDiag(testEyes.GetNbits(), refEyes.GetNbits(), "# of telescopes");
  for (unsigned int eye = 0; eye < refEyes.GetNbits(); ++eye) {
    ostringstream eyeNoS;
    eyeNoS << eye;
    diag("Testing eye " + eyeNoS.str());
    if (refEyes.TestBitNumber(eye)) {
      isDiag(test.GetDetector().GetSDFDTimeOffset(eye),
             ref.GetDetector().GetSDFDTimeOffset(eye),
             "GetSDFDTimeOffset");
      isDiag(test.GetDetector().HasAerosolData(eye),
             ref.GetDetector().HasAerosolData(eye),
             "detector.HasAerosolData(eye)");
      isEpsDiag(test.GetDetector().GetVAODAtReferenceHeight(eye),
                ref.GetDetector().GetVAODAtReferenceHeight(eye),
                1.e-9,
                "detector.GetVAODAtReferenceHeight(eye)");
      isDiag(test.GetDetector().HasVAOD(eye),
             ref.GetDetector().HasVAOD(eye),
             "detector.HasVAOD(eye)");
      isDiag(test.GetDetector().GetPointingIdsForEye(eye).size(),
             ref.GetDetector().GetPointingIdsForEye(eye).size(),
             "detector.GetPointingIdsForEye(eye)");
      isDiag(test.GetDetector().HasLidarData(eye),
             ref.GetDetector().HasLidarData(eye),
             "detector.HasLidarData(eye)");

//   bool IsInDAQ (eye,  tel)
//   TString GetPointingId(eye, tel)
//   bool HasCorrectorRing(eye, tel)

//   bool IsMismatchedPixel(eyeId, telId, pixelId)
//   bool IsBadPixel(eyeId, telId, pixelId)

      if (test.GetDetector().HasLidarData(eye) &&
          ref.GetDetector().HasLidarData(eye)) {
        isDiag(test.GetDetector().GetAerosolData(eye).GetMieAttenuationLength(),
               ref.GetDetector().GetAerosolData(eye).GetMieAttenuationLength(),
               "aerosols.att_length");
        isDiag(test.GetDetector().GetAerosolData(eye).GetMieAttenuationHeight(),
               ref.GetDetector().GetAerosolData(eye).GetMieAttenuationHeight(),
               "aerosols.att_height");
      }

      if (test.GetDetector().HasLidarData(eye) &&
          ref.GetDetector().HasLidarData(eye)) {
         isDiag(test.GetDetector().GetLidarData(eye).GetCloudCoverage(),
                ref.GetDetector().GetLidarData(eye).GetCloudCoverage(),
                "lidarData.GetCloudCoverage");
         isDiag(test.GetDetector().GetLidarData(eye).GetCloudHeight(),
                ref.GetDetector().GetLidarData(eye).GetCloudHeight(),
                "lidarData.GetCloudHeight");
         isDiag(test.GetDetector().GetLidarData(eye).GetCloudThickness(),
                ref.GetDetector().GetLidarData(eye).GetCloudThickness(),
                "lidarData.GetCloudThickness");
         isDiag(test.GetDetector().GetLidarData(eye).GetCloudDepth(),
                ref.GetDetector().GetLidarData(eye).GetCloudDepth(),
                "lidarData.GetCloudDepth");
         isDiag(test.GetDetector().GetLidarData(eye).GetCloudSlantThickness(),
                ref.GetDetector().GetLidarData(eye).GetCloudSlantThickness(),
                "lidarData.GetCloudSlantThickness");
         isDiag(test.GetDetector().GetLidarData(eye).GetCloudVAOD(),
                ref.GetDetector().GetLidarData(eye).GetCloudVAOD(),
                "lidarData.GetCloudVAOD");
         isDiag(test.GetDetector().GetLidarData(eye).GetLidarRange(),
                ref.GetDetector().GetLidarData(eye).GetLidarRange(),
                "lidarData.GetLidarRange");
      }
    }
  }
  isDiag(test.GetDetector().GetBadPixelSet().size(),
         ref.GetDetector().GetBadPixelSet().size(),
         "pixel set");


  const TBits refStations = ref.GetDetector().GetActiveStations();
  const TBits testStations = test.GetDetector().GetActiveStations();

  const std::vector<SdRecStation>& stations = ref.GetSDEvent().GetStationVector();
  for (unsigned int i = 0; i != stations.size(); ++i) {
    if (stations[i].IsDense()) {
      TVector3 ref_pos = ref.GetDetector().GetStationPosition(stations[i].GetId());
      TVector3 test_pos = test.GetDetector().GetStationPosition(stations[i].GetId());
      ostringstream name;
      name << "station " << stations[i].GetId() << " pos ";
      isDiag(ref_pos.x(),test_pos .x(), name.str());
      isDiag(ref_pos.y(),test_pos .y(), name.str());
    }
  }

  for (unsigned int id = 0; id < refStations.GetNbits(); ++id) {
    ostringstream name;
    isDiag(testStations.TestBitNumber(id), refStations.TestBitNumber(id), name.str());
  }
}

void
cmpDetectorGeometry(const RecEventFile& refFile, const RecEventFile& testFile)
{
  DetectorGeometry* refGeo  = readDetectorGeometry(refFile, "reference");
  DetectorGeometry* testGeo = readDetectorGeometry(testFile, "test");
  if (!refGeo || !testGeo)
    return;

  {
    //fail("Need to implement more DetectorGeometry comparison");
  }

  isDiag(testGeo->GetNEyes(), refGeo->GetNEyes(), "Same no of eyes in DetectorGeometry");

  for (DetectorGeometry::ConstEyeIterator ieye = refGeo->EyesBegin(),
       end = refGeo->EyesEnd(); ieye != end; ++ieye)
  {
    ok(testGeo->HasEye( ieye->first ), "Has correct eye in detgeo");
    if (!testGeo->HasEye(ieye->first))
        continue;

    isDiag(testGeo->GetEye(ieye->first).GetNumberOfMirrors(),
           refGeo->GetEye(ieye->first).GetNumberOfMirrors(),
           "eyeGeom.GetNumberOfMirrors");
    isEpsDiag(testGeo->GetEye(ieye->first).GetEyePos(),
              refGeo->GetEye(ieye->first).GetEyePos(), 1.e-6,
              "eyeGeom.GetEyePos");
    isDiag(testGeo->GetEye(ieye->first).GetEyeName(),
           refGeo->GetEye(ieye->first).GetEyeName(),
           "eyeGeom.GetEyeName");
    isDiag(testGeo->GetEye(ieye->first).GetEyeNameAbbr(),
           refGeo->GetEye(ieye->first).GetEyeNameAbbr(),
           "eyeGeom.GetEyeNameAbbr");
    isDiag(testGeo->GetEye(ieye->first).GetBackWallAngle(),
           refGeo->GetEye(ieye->first).GetBackWallAngle(),
           "eyeGeom.GetBackWallAngle");
    isDiag(testGeo->GetEye(ieye->first).GetEyePhiZ(),
           refGeo->GetEye(ieye->first).GetEyePhiZ(),
           "eyeGeom.GetEyePhiZ");
    isDiag(testGeo->GetEye(ieye->first).GetEyeThetaZ(),
           refGeo->GetEye(ieye->first).GetEyeThetaZ(),
           "eyeGeom.GetEyeThetaZ");
    isDiag(testGeo->GetEye(ieye->first).GetTelescopeIDs(),
           refGeo->GetEye(ieye->first).GetTelescopeIDs(),
           "eyeGeom.GetTelescopeIDs");

    for (EyeGeometry::TelescopeIterator telescope = refGeo->GetEye(ieye->first).TelescopesBegin(),
            end = refGeo->GetEye(ieye->first).TelescopesEnd();
         telescope != end; ++telescope) {
      ok(testGeo->GetEye(ieye->first).HasTelescope(telescope->first),
         "Has correct telescope in detgeo");
      if (!testGeo->GetEye(ieye->first).HasTelescope(telescope->first))
        continue;

      isDiag(testGeo->GetEye(ieye->first).GetTelescope(telescope->first).GetNumberOfPixels(),
             refGeo->GetEye(ieye->first).GetTelescope(telescope->first).GetNumberOfPixels(),
             "telescopeGeom.GetNumberOfPixels");
      isDiag(testGeo->GetEye(ieye->first).GetTelescope(telescope->first).GetFADCSize(),
             refGeo->GetEye(ieye->first).GetTelescope(telescope->first).GetFADCSize(),
             "telescopeGeom.GetFADCSize");
      isDiag(testGeo->GetEye(ieye->first).GetTelescope(telescope->first).GetFADCBinning(),
             refGeo->GetEye(ieye->first).GetTelescope(telescope->first).GetFADCBinning(),
             "telescopeGeom.GetFADCBinning");
      isDiag(testGeo->GetEye(ieye->first).GetTelescope(telescope->first).GetPointingIds(),
             refGeo->GetEye(ieye->first).GetTelescope(telescope->first).GetPointingIds(),
             "telescopeGeom.GetPointingIds");

      vector<TString> pointing =
         refGeo->GetEye(ieye->first).GetTelescope(telescope->first).GetPointingIds();

      for (vector<TString>::iterator it = pointing.begin(); it != pointing.end(); ++it) {
        isEpsDiag(testGeo->GetEye(ieye->first).GetTelescope(telescope->first).GetElevation(*it),
                  refGeo->GetEye(ieye->first).GetTelescope(telescope->first).GetElevation(*it),
                  1.e-9,
                  "telescopeGeom.GetElevation");
        isDiag(testGeo->GetEye(ieye->first).GetTelescope(telescope->first).GetPixelMinPhi(*it),
               refGeo->GetEye(ieye->first).GetTelescope(telescope->first).GetPixelMinPhi(*it),
               "telescopeGeom.GetPixelMinPhi");
        isDiag(testGeo->GetEye(ieye->first).GetTelescope(telescope->first).GetPixelMaxPhi(*it),
               refGeo->GetEye(ieye->first).GetTelescope(telescope->first).GetPixelMaxPhi(*it),
               "telescopeGeom.GetPixelMaxPhi");
        isDiag(testGeo->GetEye(ieye->first).GetTelescope(telescope->first).GetPixelMinOmega(*it),
               refGeo->GetEye(ieye->first).GetTelescope(telescope->first).GetPixelMinOmega(*it),
               "telescopeGeom.GetPixelMinOmega");
        isDiag(testGeo->GetEye(ieye->first).GetTelescope(telescope->first).GetPixelMaxOmega(*it),
               refGeo->GetEye(ieye->first).GetTelescope(telescope->first).GetPixelMaxOmega(*it),
               "telescopeGeom.GetPixelMaxOmega");
        isDiag(testGeo->GetEye(ieye->first).GetTelescope(telescope->first).GetPixelsPhi(*it),
               refGeo->GetEye(ieye->first).GetTelescope(telescope->first).GetPixelsPhi(*it),
               "telescopeGeom.GetPixelsPhi");
        isDiag(testGeo->GetEye(ieye->first).GetTelescope(telescope->first).GetPixelsOmega(*it),
               refGeo->GetEye(ieye->first).GetTelescope(telescope->first).GetPixelsOmega(*it),
               "telescopeGeom.GetPixelsOmega");

        const TelescopeGeometry::PixelList_t& pixels_phi = telescope->second.GetPixelsPhi(*it);
        for (unsigned int i = 0; i != pixels_phi.size(); ++i) {
          isDiag(testGeo->GetEye(ieye->first).GetTelescope(telescope->first).HasPixel(pixels_phi[i]),
                 refGeo->GetEye(ieye->first).GetTelescope(telescope->first).HasPixel(pixels_phi[i]),
                 "telescopeGeom.HasPixel");
          if (!testGeo->GetEye(ieye->first).GetTelescope(telescope->first).HasPixel(pixels_phi[i]))
            continue;

          isDiag(testGeo->GetEye(ieye->first).GetTelescope(telescope->first).GetPixelPhi(pixels_phi[i], *it),
                 refGeo->GetEye(ieye->first).GetTelescope(telescope->first).GetPixelPhi(pixels_phi[i], *it),
                 "telescopeGeom.GetPixelPhi");
          isDiag(testGeo->GetEye(ieye->first).GetTelescope(telescope->first).GetPixelOmega(pixels_phi[i], *it),
                 refGeo->GetEye(ieye->first).GetTelescope(telescope->first).GetPixelOmega(pixels_phi[i], *it),
                 "telescopeGeom.GetPixelOmega");
        }

        const TelescopeGeometry::PixelList_t& pixels_omega = telescope->second.GetPixelsOmega(*it);
        for (unsigned int i = 0; i != pixels_omega.size(); ++i) {
          isDiag(testGeo->GetEye(ieye->first).GetTelescope(telescope->first).HasPixel(pixels_omega[i]),
                 refGeo->GetEye(ieye->first).GetTelescope(telescope->first).HasPixel(pixels_omega[i]),
                 "telescopeGeom.HasPixel");
          if (!testGeo->GetEye(ieye->first).GetTelescope(telescope->first).HasPixel(pixels_omega[i]))
            continue;

          isDiag(testGeo->GetEye(ieye->first).GetTelescope(telescope->first).GetPixelPhi(pixels_omega[i], *it),
                 refGeo->GetEye(ieye->first).GetTelescope(telescope->first).GetPixelPhi(pixels_omega[i], *it),
                 "telescopeGeom.GetPixelPhi");
          isDiag(testGeo->GetEye(ieye->first).GetTelescope(telescope->first).GetPixelOmega(pixels_omega[i], *it),
                 refGeo->GetEye(ieye->first).GetTelescope(telescope->first).GetPixelOmega(pixels_omega[i], *it),
                 "telescopeGeom.GetPixelOmega");
        }
      }
    }
  }

  isDiag(testGeo->GetNStations(), refGeo->GetNStations(), "Same no of stations in DetectorGeometry");

  DetectorGeometry::ConstStationPosMap pos = testGeo->GetStations();
  for (DetectorGeometry::StationPosMapConstIterator is = refGeo->GetStationsBegin(),
       end = refGeo->GetStationsEnd(); is != end; ++is)
  {
    const unsigned int id = is->first;
    if (!ok(pos.find(id) != pos.end(), "station exists in test geo")) {
      ostringstream msg;
      msg << "Station " << id << " is missing from test DetectorGeometry";
      diag(msg);
      continue;
    }
    isDiag(testGeo->GetStationName(id), refGeo->GetStationName(id), "Station name identical");
    isEpsDiag(testGeo->GetStationPosition(id), refGeo->GetStationPosition(id), 0.001, "Station position identical");
    isDiag(testGeo->GetHexagon(id), refGeo->GetHexagon(id), "Station hexagon");
    isDiag(testGeo->GetFdTriggerReadoutSection(id), refGeo->GetFdTriggerReadoutSection(id), "Station's FD readout section");
  }


  delete refGeo;
  delete testGeo;
}


void
cmpEvent(const RecEvent& ref, const RecEvent& test)
{
  cmpEventGlobalInfo(ref, test);

  isDiag(test.GetNEyes(), ref.GetNEyes(), "Same no. of eyes");

  // Compare FDEvents
  for (RecEvent::ConstEyeIterator iEye = ref.EyesBegin();
       iEye != ref.EyesEnd(); ++iEye)
  {
    if (ok(test.HasEye(iEye->GetEyeId()), "Eye exists")) {
      cmpFDEvent(*iEye, test.GetEye(iEye->GetEyeId()));
    }
    else {
      ostringstream m;
      m << "Missing eye from test event: " << iEye->GetEyeId();
      diag(m);
    }
  }

  // Compare SDEvent
  cmpSDEvent(ref.GetSDEvent(), test.GetSDEvent());

  // Compare gen shower
  cmpGenShower(ref.GetGenShower(), test.GetGenShower());

  cmpDetector(ref, test);
}


void
cmpEventGlobalInfo(const RecEvent& ref, const RecEvent& test)
{
  isDiag(test.GetEventId(), ref.GetEventId(), "RecEvent.GetEventId");
  isDiag(test.GetYYMMDD(), ref.GetYYMMDD(), "RecEvent.GetYYMMDD");
  isDiag(test.GetHHMMSS(), ref.GetHHMMSS(), "RecEvent.GetHHMMSS");
  isDiag(test.GetMJD(), ref.GetMJD(), "RecEvent.GetMJD");

  isDiag(test.GetDetector().HasMieDatabase(), ref.GetDetector().HasMieDatabase(), "Detector.HasMieDatabase");

  // FIXME what else?
}


void
cmpGenShower(const GenShower& ref, const GenShower& test)
{
// [a-z0-9A-Z]+
//  isDiag(test.GetSdRecShower().\&(), ref.GetSdRecShower().\&(), "sdRecShower.\&");
  isDiag(test.GetPrimary(), ref.GetPrimary(), "genShower.GetPrimary");
  isDiag(test.GetPrimaryName(), ref.GetPrimaryName(), "genShower.GetPrimaryName");
  isDiag(test.GetShortPrimaryName(), ref.GetShortPrimaryName(), "genShower.GetShortPrimaryName");
  isDiag(test.GetElecEnergy(), ref.GetElecEnergy(), "genShower.GetElecEnergy");
  isDiag(test.GetX1(), ref.GetX1(), "genShower.GetX1");
  isDiag(test.GetXmax(), ref.GetXmax(), "genShower.GetXmax");
  isDiag(test.GetdEdXmax(), ref.GetdEdXmax(), "genShower.GetdEdXmax");
  isDiag(test.GetElectrons(), ref.GetElectrons(), "genShower.GetElectrons");
  isDiag(test.GetEnergyDeposit(), ref.GetEnergyDeposit(), "genShower.GetEnergyDeposit");
  isDiag(test.GetDepth(), ref.GetDepth(), "genShower.GetDepth");
  isDiag(test.GetShowerAge(), ref.GetShowerAge(), "genShower.GetShowerAge");
//#warning this gives me a segmentation fault (?)
  //isDiag(test.GetCoreTimeNanoSecond(), ref.GetCoreTimeNanoSecond(), "genShower.GetCoreTimeNanoSecond");
}