HDFFileReader.cc

Go to the documentation of this file.

#include "dune-raw-data/HDFUtils/HDFFileReader.hh"
#include "dune-raw-data/Overlays/FragmentType.hh"
#include "dune-raw-data/Overlays/FelixFragment.hh"
#include "dune-raw-data/Overlays/SSPFragment.hh"
#include "artdaq-core/Utilities/TimeUtils.hh"
#include <boost/algorithm/string.hpp>
#include <regex>
#include <set>
#include <sstream>

#include <sys/time.h>
#include <sys/resource.h>

namespace dune {
  namespace HDFFileReader {

    // *****************************************
    // *** General-purpose data and routines ***
    // *****************************************

    HDFFileInfoPtr openFile(const std::string& fileName)
    {
      HDFFileInfoPtr hdfFileInfoPtr(new HDFFileInfo());
      hdfFileInfoPtr->filePtr = H5Fopen(fileName.data(), H5F_ACC_RDONLY, H5P_DEFAULT);
      hdfFileInfoPtr->bytesWritten = 0;
      hdfFileInfoPtr->fileName = fileName;
      hdfFileInfoPtr->runNumber = 0;   // don't know it yet.  Fill it when we get fragments.
      hdfFileInfoPtr->dataFormatVersion = 0;

      hid_t grp = H5Gopen(hdfFileInfoPtr->filePtr,"/", H5P_DEFAULT);
      if (attrExists(grp,"data_format_version"))
        {
          hid_t ga = H5Aopen_name(grp,"data_format_version");
          H5Aread(ga, H5Aget_type(ga), &hdfFileInfoPtr->dataFormatVersion);
          H5Aclose(ga);
        }
      H5Gclose(grp);
      return hdfFileInfoPtr;
    }

    HDFFileInfoPtr reopenFile(HDFFileInfoPtr oldHdfFileInfoPtr)
    {
      H5Fclose(oldHdfFileInfoPtr->filePtr);
      oldHdfFileInfoPtr->filePtr = 0;

      HDFFileInfoPtr newHdfFileInfoPtr(new HDFFileInfo());
      newHdfFileInfoPtr->filePtr = H5Fopen(oldHdfFileInfoPtr->fileName.data(), H5F_ACC_RDONLY, H5P_DEFAULT);
      newHdfFileInfoPtr->bytesWritten = 0;
      newHdfFileInfoPtr->fileName = oldHdfFileInfoPtr->fileName;
      newHdfFileInfoPtr->runNumber = oldHdfFileInfoPtr->runNumber; // gets overwritten but fill it anyway
      newHdfFileInfoPtr->dataFormatVersion = oldHdfFileInfoPtr->dataFormatVersion;
      return newHdfFileInfoPtr;
    }

    void closeFile(HDFFileInfoPtr hdfFileInfoPtr)
    {
      H5Fclose(hdfFileInfoPtr->filePtr);
      hdfFileInfoPtr->filePtr = 0;
    }

    std::list<std::string> getTopLevelGroupNames(HDFFileInfoPtr& hdfFileInfoPtr)
    {
      hid_t grp = H5Gopen(hdfFileInfoPtr->filePtr,"/", H5P_DEFAULT);
      std::list<std::string> theList = getMidLevelGroupNames(grp);
      H5Gclose(grp);
      return theList;
    }

    std::list<std::string> getMidLevelGroupNames(hid_t grp)
    {
      std::list<std::string> theList;
      hsize_t nobj = 0;
      H5Gget_num_objs(grp,&nobj);
      for (hsize_t idx = 0; idx < nobj; ++idx)
        {
          hsize_t len = H5Gget_objname_by_idx(grp, idx, NULL, 0 );
          char *memb_name = new char(len+1);
          H5Gget_objname_by_idx(grp, idx, memb_name, len+1 );
          theList.emplace_back(memb_name);
          delete[] memb_name;
        }
      return theList;
    }

    hid_t getGroupFromPath(HDFFileInfoPtr& hdfFileInfoPtr, const std::string &path)
    {
      hid_t grp = H5Gopen(hdfFileInfoPtr->filePtr,path.data(),H5P_DEFAULT);
      return grp;
    }

    // check to see if an attribute exists -- quiet the error messages when probing
    // using example from the HDF5 documentation

    bool attrExists(hid_t object, const std::string &attrname)
    {
      // Save old error handler 
      H5E_auto_t old_func;
      void *old_client_data;
      H5Eget_auto(H5E_DEFAULT,&old_func, &old_client_data);

      // Turn off error handling */
      H5Eset_auto(H5E_DEFAULT,NULL, NULL);

      // Probe. On failure, retval is supposed to be negative

      hid_t retval = H5Aopen_name(object, attrname.data());

      // Restore previous error handler 
      H5Eset_auto(H5E_DEFAULT,old_func, old_client_data);

      bool result = (retval >= 0);
      return result;
    }

    // ************************************************
    // *** Trigger-Record-related data and routines ***
    // ************************************************

    // We don't know the run number until we read some fragments in, so this method modifies the run number
    // in hdfFileInfoPtr

    FragmentListsByType getFragmentsForEvent(HDFFileInfoPtr& hdfFileInfoPtr, const std::string& topLevelGroupName)
    {
      FragmentListsByType fragmentMap;
      artdaq::Fragment::fragment_id_t sampleFragmentID = 1001;

      hid_t requestedGroup = getGroupFromPath(hdfFileInfoPtr, topLevelGroupName);
      {
        std::list<std::string> detectorTypeNames = getMidLevelGroupNames(requestedGroup);
        for (auto& detectorTypeName : detectorTypeNames)
          {
            if (detectorTypeName != "TPC") continue;  // skip non-TPC fragments for now
            std::string dtnameloc = detectorTypeName;
            if (dtnameloc == "TPC") 
              {
                dtnameloc = "FELIX";  // label the fragments as FELIX fragments for later consumption
              }
            fragmentMap[dtnameloc] = std::make_unique<artdaq::Fragments>();

            std::string subdetGroupPath = topLevelGroupName + "/" + detectorTypeName;
            hid_t subdetGroup = getGroupFromPath(hdfFileInfoPtr, subdetGroupPath);

            // this loop doesn't make sense for TriggerRecordHeader datasets

            std::list<std::string> subdetGeoNames = getMidLevelGroupNames(subdetGroup);
            for (auto& subdetGeoName : subdetGeoNames)
              {
                std::string geoGroupPath = subdetGroupPath + "/" + subdetGeoName;
                //std::cout << "     Geo path: " << geoGroupPath << std::endl;
                hid_t geoGroup = getGroupFromPath(hdfFileInfoPtr,geoGroupPath);
                std::list<std::string> dataSetNames = getMidLevelGroupNames(geoGroup);
                for (auto& dataSetName : dataSetNames)  
                  {
                    std::string dataSetPath = geoGroupPath + "/" + dataSetName;
                    //std::cout << "Dataset path: " << dataSetPath << std::endl;

                    // read a dataset and put it in a buffer so we can extract metadata
                    // Future optimization -- just read in the first part of the dataset to extract
                    // the header, and read in the rest of it directly into the fragment.
                    // need to know header size to do that.  80 bytes for v2 TPC data.

                    hid_t datasetid = H5Dopen(geoGroup,dataSetName.data(),H5P_DEFAULT);
                    hsize_t ds_size = H5Dget_storage_size(datasetid);
                    std::vector<char> ds_data(ds_size);   
                    H5Dread(datasetid, H5T_STD_I8LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, ds_data.data());
                    H5Dclose(datasetid);

                    artdaq::FragmentPtr frag;
                    artdaq::Fragment::type_t type = artdaq::Fragment::InvalidFragmentType;
                    uint32_t payloadOffset = 0;
                    if (detectorTypeName == TPC_GROUP_NAME)
                      {
                        if (ds_size < 81) continue;  // skip if there isn't even enough for a header
                        frag = artdaq::Fragment::FragmentBytes(ds_size - 80);
                        type = dune::detail::FELIX;
                        frag->setUserType(type);
                        size_t numframes = (ds_size - 80)/sizeof(dune::FelixFrame);

                        //int magic_word = 0;
                        //memcpy(&magic_word,&ds_data[0],4);
                        //std::cout << "   Magic word: 0x" << std::hex << magic_word << std::dec << std::endl;

                        //int version = 0;
                        //memcpy(&version, &ds_data[4],4);
                        //std::cout << "   Version: " << std::dec << version << std::dec << std::endl;

                        //uint64_t fragsize=0;
                        //memcpy(&fragsize, &ds_data[8],8);
                        //std::cout << "   Frag Size: " << std::dec << fragsize << std::dec << std::endl;

                        uint64_t trignum=0;
                        memcpy(&trignum, &ds_data[16],8);
                        //std::cout << "   Trig Num: " << std::dec << trignum << std::dec << std::endl;
                        frag->setSequenceID(trignum);

                        uint64_t triggerTime = 0;
                        memcpy(&triggerTime, &ds_data[24],8);
                        //std::cout << "   Trig Timestamp: " << std::dec << trig_timestamp << std::dec << std::endl;

                        uint64_t windowStart=0;
                        memcpy(&windowStart, &ds_data[32],8);
                        //std::cout << "   Window Begin:   " << std::dec << windowbeg << std::dec << std::endl;

                        uint64_t windowEnd=0;
                        memcpy(&windowEnd, &ds_data[40],8);
                        //std::cout << "   Window End:     " << std::dec << windowend << std::dec << std::endl;

                        int runno=0;
                        memcpy(&runno, &ds_data[48], 4);
                        //std::cout << "   Run Number: " << std::dec << runno << std::endl;
                        hdfFileInfoPtr->runNumber = runno;  // keep overwriting this with each dataset

                        //int errbits=0;
                        //memcpy(&errbits, &ds_data[52], 4);
                        //std::cout << "   Error bits: " << std::dec << errbits << std::endl;

                        //int fragtype=0;
                        //memcpy(&fragtype, &ds_data[56], 4);
                        //std::cout << "   Fragment type: " << std::dec << fragtype << std::endl;

                        //int fragpadding=0;
                        //memcpy(&fragtype, &ds_data[60], 4);
                        //std::cout << "   Fragment padding: " << std::dec << fragpadding << std::endl;

                        //int geoidversion=0;
                        //memcpy(&geoidversion, &ds_data[64], 4);
                        //std::cout << "   GeoID version: " << std::dec << geoidversion << std::endl;

                        //unsigned short geoidtype;
                        //memcpy(&geoidtype, &ds_data[70], 1);
                        //std::cout << "   GeoID type: " << geoidtype << std::endl;

                        //unsigned short geoidregion=0;
                        //memcpy(&geoidregion, &ds_data[71], 1);
                        //std::cout << "   GeoID region: " << std::dec << geoidregion << std::endl;

                        //int geoidelement=0;
                        //memcpy(&geoidelement, &ds_data[72], 4);
                        //std::cout << "   GeoID element: " << std::dec << geoidelement << std::endl;

                        //int geoidpadding=0;
                        //memcpy(&geoidpadding, &ds_data[76], 4);
                        //std::cout << "   GeoID padding: " << std::dec << geoidpadding << std::endl;

                        dune::FelixFragmentBase::Metadata felixMetadata;
                        felixMetadata.control_word = 0xabc;
                        felixMetadata.version = 1;
                        felixMetadata.reordered = 0;
                        felixMetadata.compressed = 0;
                        felixMetadata.offset_frames = (uint32_t) ((triggerTime - windowStart) / 25);
                        felixMetadata.window_frames = (uint32_t) ((windowEnd - windowStart) / 25);
                        felixMetadata.num_frames = numframes;

                        frag->setMetadata(felixMetadata);
                        frag->setTimestamp(windowStart);

                        memcpy(frag->dataBeginBytes() + payloadOffset, &ds_data[80], ds_size-80);
                      }

                    uint64_t pdsPayloadHeaderWord = 0;
                    if (detectorTypeName == PDS_GROUP_NAME)  
                      {
                        payloadOffset = sizeof(uint64_t);
                        frag = artdaq::Fragment::FragmentBytes(ds_size + sizeof(uint64_t));
                        type = dune::detail::PHOTON;

                        // fill in the SSPFragment metadata, as much as we can
                        SSPDAQ::MillisliceHeader pdsMetadata;
                        pdsMetadata.triggerTime = 0;
                        pdsMetadata.triggerType = 0;

                        if (attrExists(requestedGroup,"trigger_timestamp"))
                          {
                            uint64_t timestamp;
                            hid_t a1 = H5Aopen_name(requestedGroup,"trigger_timestamp");
                            H5Aread(a1, H5Aget_type(a1), &timestamp);
                            H5Aclose(a1);
                            pdsMetadata.triggerTime = timestamp * 3;
                          }

                        if (attrExists(requestedGroup,"timeslice_start"))
                          {
                            uint64_t timestamp;
                            hid_t a1 = H5Aopen_name(requestedGroup,"timeslice_start");
                            H5Aread(a1, H5Aget_type(a1), &timestamp);
                            H5Aclose(a1);
                            pdsMetadata.startTime = timestamp * 3;
                            if (pdsMetadata.triggerTime == 0) {pdsMetadata.triggerTime = timestamp * 3;}
                          }
                        else
                          {
                            if (attrExists(datasetid,"timestamp_of_first_packet"))
                              {
                                uint64_t timestamp;
                                hid_t a1 = H5Aopen_name(datasetid,"timestamp_of_first_packet");
                                H5Aread(a1, H5Aget_type(a1), &timestamp);
                                H5Aclose(a1);
                                pdsMetadata.startTime = timestamp * 3;
                              }
                          }


                        if (attrExists(requestedGroup,"timeslice_end"))
                          {
                            uint64_t timestamp;
                            hid_t a1 = H5Aopen_name(requestedGroup,"timeslice_end");
                            H5Aread(a1, H5Aget_type(a1), &timestamp);
                            H5Aclose(a1);
                            pdsMetadata.endTime = timestamp * 3;
                          }
                        else
                          {
                            if (attrExists(datasetid,"timestamp_of_last_packet"))
                              {
                                uint64_t timestamp;
                                hid_t a1 = H5Aopen_name(datasetid,"timestamp_of_last_packet");
                                H5Aread(a1, H5Aget_type(a1), &timestamp);
                                H5Aclose(a1);
                                pdsMetadata.endTime = timestamp * 3;
                              }
                          }


                        if (attrExists(datasetid,"number_of_packets"))
                          {
                            uint32_t packetCount;
                            hid_t a1 = H5Aopen_name(datasetid,"number_of_packets");
                            H5Aread(a1, H5Aget_type(a1), &packetCount);
                            H5Aclose(a1);
                            pdsMetadata.nTriggers = packetCount;
                          }

                        if (attrExists(datasetid,"size_in_bytes"))
                          {
                            uint32_t sizeBytes;
                            hid_t a1 = H5Aopen_name(datasetid,"size_in_bytes");
                            H5Aread(a1, H5Aget_type(a1), &sizeBytes);
                            H5Aclose(a1);
                            pdsMetadata.length = 9 + (sizeBytes / sizeof(uint32_t));
                            pdsPayloadHeaderWord |= ((2 + (sizeBytes / sizeof(uint32_t))) & 0xffffffff);
                          }
                        else
                          {
                            pdsMetadata.length = 9 + (ds_size / sizeof(uint32_t));
                            pdsPayloadHeaderWord |= ((2 + (ds_size / sizeof(uint32_t))) & 0xffffffff);
                          }

                        frag->setMetadata(pdsMetadata);
                        pdsPayloadHeaderWord |= 0x000003e700000000;  // run number always seems to be set to 999
                      }

                    artdaq::Fragment::sequence_id_t seqID;
                    artdaq::Fragment::fragment_id_t fragID;

                    std::stringstream seqStr(topLevelGroupName);
                    seqStr >> seqID;
                    // sequence ID in artdaq::Fragment header is only 48 bits
                    seqID &= 0xffffffffffff;

                    fragID = sampleFragmentID++;

                    frag->setFragmentID(fragID);
                    //frag->setUserType(type);  // set this in the relevant block
                    //frag->setSequenceID(seqID);  // to be read out of each fragment's header
                    //frag->setTimestamp(headerTimestamp);  // to be read out of each fragment's header

                    if (detectorTypeName == PDS_GROUP_NAME)
                      {
                        memcpy(frag->dataBeginBytes(), &pdsPayloadHeaderWord, sizeof(uint64_t));
                      }

                    fragmentMap[dtnameloc]->emplace_back(std::move(*(frag.release())));

                  }
              }
            H5Gclose(subdetGroup);
          }
      }

      H5Gclose(requestedGroup);

      //hdfFileInfoPtr->filePtr->flush(H5F_SCOPE_GLOBAL);
      return fragmentMap;
    }

  } // namespace HDFFileReader
}   // namespace dune