GGSRootFileService.cpp

/*
 * GGSRootFileService.cpp
 *
 *  Created on: 03 Jun 2011
 *      Author: Nicola Mori
 */

#include "montecarlo/services/GGSRootFileService.h"
#include "geometry/pluginmanagers/GGSGeoPluginManager.h"
#include "montecarlo/dataobjs/GGSTGeoParams.h"
#include "montecarlo/dataobjs/GGSTHitDetInfo.h"
#include "montecarlo/dataobjs/GGSTHitVolInfo.h"
#include "montecarlo/dataobjs/GGSTParameters.h"
#include "montecarlo/mccore/GGSRunManager.h"
#include "utils/GGSSmartLog.h"

#include "G4LogicalVolume.hh"
#include "G4Run.hh"
#include "G4SystemOfUnits.hh"

#include "ROOT/TBufferMerger.hxx"
#include "TTree.h"

#include <sstream>

// Definition of the std::hash<GGSRootFileService::VolumeKey> specialization
namespace std {
// The standard library does not provide an equivalent of boost::hash_combine, so
// the implementation from Boost 1.75.0 for 64 bit machines is copy-pasted here.
inline void hash_combine_impl(uint64_t &h, uint64_t k) {
  constexpr uint64_t m = UINT64_C(0xc6a4a7935bd1e995);
  constexpr int r = 47;

  k *= m;
  k ^= k >> r;
  k *= m;

  h ^= k;
  h *= m;

  // Completely arbitrary number, to prevent 0's
  // from hashing to 0.
  h += 0xe6546b64;
}
template <class T> inline void hash_combine(uint64_t &seed, T const &v) {
  std::hash<T> hasher;
  hash_combine_impl(seed, hasher(v));
}

// With the above definition of hash_combine, define the hash struct for GGSRootFileService::VolumeKey
template <> struct hash<GGSRootFileService::VolumeKey> {
  uint64_t operator()(GGSRootFileService::VolumeKey const &key) const noexcept {
    uint64_t seed = 0;
    hash_combine(seed, key.name);
    hash_combine(seed, key.physVol);
    hash_combine(seed, key.position[0]);
    hash_combine(seed, key.position[1]);
    hash_combine(seed, key.position[2]);
    return seed;
  }
};
} // namespace std

// Variables used for performance in StoreVolume
G4ThreadLocal TFile *GGSRootFileService::_currVolStorageTFile = nullptr;
G4ThreadLocal std::string GGSRootFileService::_currDetVolName = "";
G4ThreadLocal GGSRootFileService::FileInfoContainer::iterator GGSRootFileService::_currVolStorageFileInfo;
G4ThreadLocal GGSRootFileService::HitDetMap::iterator GGSRootFileService::_currDetector;
G4ThreadLocal GGSTHitDetInfo *GGSRootFileService::_currPersDetector = nullptr;

GGSRootFileService &GGSRootFileService::GetInstance() {
  static GGSRootFileService instance;
  return instance;
}

GGSRootFileService::GGSRootFileService()
    : _suffix(""), _defaultOutBase("GGSRootOutput"), _singleFileForMT(true), _messenger(this, "/GGS/rootOutput/") {
  _messenger.DeclareMethod("setMultiFileForMT", &GGSRootFileService::_SetMultiFileMode,
                           "Set the output mode for MT mode to one file per thread");
}

namespace {
G4Mutex mutexForGetFile = G4MUTEX_INITIALIZER;
}
TFile *GGSRootFileService::GetFileForThisRun(const std::filesystem::path &baseName, const G4Run *run) {

  G4AutoLock lock(mutexForGetFile);

  std::filesystem::path absPath;

  // Set default if no base name has been provided
  if (baseName.string() == "") {
    absPath = _defaultOutBase;
  } else {
    absPath = baseName;
  }

  // Retrieve absolute path
  absPath = _GetAbsolutePath(absPath);

  // Strip extension and append suffix and extension
  absPath = _AppendSuffixAndExt(absPath, run);

  // Lambda that creates a new TFile for current thread and for the given file info
  auto CreateNewTFile = [this, &absPath](FileInfoContainer::iterator fileInfoIter) -> TFile * {
    std::shared_ptr<TFile> newTFile = nullptr;
    if (G4Threading::IsMultithreadedApplication() && _singleFileForMT) {
      // Create the TFile with TBufferMerger in MT single-file mode
      if (!fileInfoIter->bufferMerger) {
        fileInfoIter->bufferMerger = std::make_unique<ROOT::TBufferMerger>(absPath.string().c_str());
      }
      newTFile = fileInfoIter->bufferMerger->GetFile();
    } else {
      // Create the TFile normally in ST mode or in MT multi-file mode
      if (fileInfoIter->tFilesInfo.size() != 0) {
        G4Exception("GGSRootFileService::GetFileForThisRun", "TFile clash", G4ExceptionSeverity::FatalException,
                    ("A TFile is already existing for file " + absPath.string() + " in ST mode; can't create it again")
                        .c_str());
      }
      newTFile = std::make_shared<TFile>(absPath.string().c_str(), "RECREATE");
    }
    fileInfoIter->tFilesInfo.emplace_back();
    auto &tFileInfo = fileInfoIter->tFilesInfo.back();
    tFileInfo.tFile = newTFile;
    tFileInfo.defaultEventsTree = nullptr;
    tFileInfo.threadNo = G4Threading::G4GetThreadId();
    tFileInfo.nRequests = 1;

    return newTFile.get();
  };

  // Check if the file info is already present
  int iThread = G4Threading::G4GetThreadId();
  auto fileInfoIter = _filesInfo.begin();
  for (; fileInfoIter != _filesInfo.end(); ++fileInfoIter) {
    if (absPath == fileInfoIter->absPath) {
      // File info found. Check for the TFile
      for (auto &tFileInfo : fileInfoIter->tFilesInfo) {
        if (tFileInfo.threadNo == iThread) {
          ++(tFileInfo.nRequests);
          return tFileInfo.tFile.get();
        }
      }
      // TFile not found. Create and return it
      return CreateNewTFile(fileInfoIter);
    }
  }

  // File info has not been found, so create it
  _filesInfo.emplace_back();
  fileInfoIter = --(_filesInfo.end());
  fileInfoIter->absPath = absPath;

  // Create the TFile
  TFile *newTFile = CreateNewTFile(fileInfoIter);

  // Set storage of detector to nullptr values to force proper initialization in StoreVolume
  _currVolStorageTFile = nullptr;
  _currDetVolName = "";

  return newTFile;
}

namespace {
G4Mutex mutexForCloseFile = G4MUTEX_INITIALIZER;
}
void GGSRootFileService::CloseFileForThisRun(const std::filesystem::path &baseName, const G4Run *run) {

  G4AutoLock lock(mutexForCloseFile);

  std::filesystem::path absPath;

  // Set default if no base name has been provided
  if (baseName.string() == "") {
    absPath = _defaultOutBase;
  } else
    absPath = baseName;
  // Retrieve absolute path
  absPath = _GetAbsolutePath(absPath);

  // Strip extension and append suffix and extension
  absPath = _AppendSuffixAndExt(absPath, run);

  // Helper lambda that writes the general simulation information on the output file
  auto WriteGeneralSimInfo = [this](TFile *outFile) {
    // Write simulation informations
    outFile->WriteObjectAny(&_simInfo, "GGSTSimInfo", _simInfo.GetName());

    // Write geometry informations
    GGSVGeometryConstruction *geoCons = GGSGeoPluginManager::GetInstance().GetGeoConstruction();
    GGSTGeoParams geoParams;
    geoParams.SetIntGeoParams(geoCons->GetIntParameters());
    geoParams.SetBoolGeoParams(geoCons->GetBoolParameters());
    geoParams.SetRealGeoParams(geoCons->GetRealParameters());
    geoParams.SetStringGeoParams(geoCons->GetStringParameters());
    geoParams.SetVectIntGeoParams(geoCons->GetVectIntParameters());
    geoParams.SetVectBoolGeoParams(geoCons->GetVectBoolParameters());
    geoParams.SetVectRealGeoParams(geoCons->GetVectRealParameters());
    geoParams.SetVectStringGeoParams(geoCons->GetVectStringParameters());
    outFile->WriteObject(&geoParams, "GGSGeoParams");

    // Write generator informations
    // These are written only for GGS generators that export a parameter object containing the "generator" name
    // setting.
    auto *genAction =
        dynamic_cast<const GGSGeneratorAction *>(GGSRunManager::GetRunManager()->GetUserPrimaryGeneratorAction());
    if (genAction) {
      GGSTParameters generatorParams(genAction->GetParameters());
      try {
        generatorParams.GetParam<std::string>("generator"); // Throws if "generator" parameter is not present
        outFile->WriteObject(&generatorParams, "GGSGenParams");
      } catch (...) {
      }
    }
  };

  int iThread = G4Threading::G4GetThreadId();
  for (auto fileInfoIter = _filesInfo.begin(); fileInfoIter != _filesInfo.end(); ++fileInfoIter) {
    if (fileInfoIter->absPath == absPath) {
      for (auto tFileInfoIter = fileInfoIter->tFilesInfo.begin(); tFileInfoIter != fileInfoIter->tFilesInfo.end();
           ++tFileInfoIter) {
        if (tFileInfoIter->threadNo == iThread) {
          --(tFileInfoIter->nRequests);
          if (tFileInfoIter->nRequests == 0) {
            if (tFileInfoIter->defaultEventsTree) {
              tFileInfoIter->defaultEventsTree->Write(); // The tree is already associated to the right TFile
            }
            // Write the detectors array and the general simulation information in the file, but if there is a buffer
            // merger then do it later just on the merged file
            if (!(fileInfoIter->bufferMerger)) {
              if (fileInfoIter->detectorsArray) {
                tFileInfoIter->tFile->WriteObject(fileInfoIter->detectorsArray, "GGSHitDetInfo", "SingleKey");
              }
              WriteGeneralSimInfo(tFileInfoIter->tFile.get());
            }
            // Mute due to this:
            // https://root-forum.cern.ch/t/warning-message-when-merging-files-with-tbuffermerger/42349
            // WARNING: this will suppress any output message from other threads occurring at the same time.
            GGSSmartLog::MuteOutput();
            tFileInfoIter->tFile->Write();
            GGSSmartLog::UnmuteOutput();
            tFileInfoIter = --(fileInfoIter->tFilesInfo.erase(tFileInfoIter));
          } else {
            return;
          }
        }
      }
      if (fileInfoIter->tFilesInfo.size() == 0) {
        if (fileInfoIter->bufferMerger) {
          // Delete the buffer merger
          fileInfoIter->bufferMerger.reset();
          // Open the file for writing global objects
          auto outFile = std::unique_ptr<TFile>(TFile::Open(fileInfoIter->absPath.string().c_str(), "UPDATE"));
          if (fileInfoIter->detectorsArray) {
            outFile->WriteObject(fileInfoIter->detectorsArray, "GGSHitDetInfo", "SingleKey");
          }
          WriteGeneralSimInfo(outFile.get());
          outFile->Close();
        }
        fileInfoIter = --(_filesInfo.erase(fileInfoIter));
      }
    }
  }
}

TTree *GGSRootFileService::GetDefaultTree(TFile *file) {
  // Don't allow other threads to ask for a file in the meantime since this might invalidate the references in the
  // loops
  G4AutoLock lock(mutexForGetFile);

  for (auto &fileInfo : _filesInfo) {
    for (auto &tFileInfo : fileInfo.tFilesInfo) {
      if (file == tFileInfo.tFile.get()) {
        if (tFileInfo.defaultEventsTree == nullptr) {
          tFileInfo.defaultEventsTree = new TTree("GGSEventsTree", "GGS events tree. Info: ", 99, file);
        }
        return tFileInfo.defaultEventsTree;
      }
    }
  }
  return nullptr;
}

std::filesystem::path GGSRootFileService::_GetAbsolutePath(const std::filesystem::path &baseName) {

  std::filesystem::path absPath(baseName);

  absPath = absolute(baseName);

  return absPath;
}

std::filesystem::path GGSRootFileService::_AppendSuffixAndExt(const std::filesystem::path &baseName, const G4Run *run) {

  std::filesystem::path absPath(baseName);
  absPath.replace_extension(); // Remove the extension

  // Append thread no. in MT multi-file mode
  if (G4Threading::IsMultithreadedApplication() && !_singleFileForMT) {
    absPath += "_Thread" + std::to_string(G4Threading::G4GetThreadId());
  }

  // if runID > 0 append _RunID
  int runID = run->GetRunID();
  if (runID > 0) {
    absPath += "_Run" + std::to_string(runID);
  }

  // Append suffix
  if (_suffix != "") {
    absPath += "_" + _suffix;
  }

  // Append extension
  absPath += ".root";

  return absPath;
}

void GGSRootFileService::_SetMultiFileMode() {
  const std::string routineName("GGSRootFileService::_SetMultiFileMode");
  if (G4Threading::IsMultithreadedApplication()) {
    _singleFileForMT = false;
  } else {
    GGSCOUT(WARNING) << "ST application: ignoring multi-file output setting" << GGSENDL;
  }
}

void GGSRootFileService::_FillDefaultEventsTrees() {
  // No need to acquire a lock since there's always one TFile (either a plain TFile or a TBufferMergerFile) per
  // thread.
  for (auto &fileInfo : _filesInfo) {
    for (auto &tFileInfo : fileInfo.tFilesInfo) {
      if (tFileInfo.tFile && tFileInfo.defaultEventsTree && tFileInfo.threadNo == G4Threading::G4GetThreadId()) {
        tFileInfo.tFile->cd();
        tFileInfo.defaultEventsTree->Fill();
      }
    }
  }
}

namespace {
G4Mutex mutexForStoreVolume = G4MUTEX_INITIALIZER;
}
int GGSRootFileService::StoreVolume(TFile *tFilePtr, const std::string &detector, const G4VPhysicalVolume *volume,
                                    const G4ThreeVector &position, int id) {
  static const std::string routineName("GGSRootFileService::StoreVolume");

  G4AutoLock lock(mutexForStoreVolume);
  std::string detVolName = detector.substr(0, detector.find_first_of('.'));

  // 1. Retrieve file info structure
  bool fileChanged = false, found = false;
  if (tFilePtr != _currVolStorageTFile) {
    fileChanged = true;
    _currVolStorageTFile = tFilePtr;

    for (_currVolStorageFileInfo = _filesInfo.begin(); _currVolStorageFileInfo != _filesInfo.end();
         ++_currVolStorageFileInfo) {
      for (auto &tFileInfo : _currVolStorageFileInfo->tFilesInfo) {
        if (tFileInfo.tFile.get() == tFilePtr) {
          found = true;
          break;
        }
      }
      if (found) {
        break;
      }
    }
    if (!found) {
      G4Exception("GGSRootFileService::StoreVolume", "TFile not found", G4ExceptionSeverity::FatalException,
                  (std::string("TFile for ") + tFilePtr->GetName() + " is not managed by GGSRootFileService").c_str());
    }
  }

  // 2. Find the detector in the index
  if (fileChanged || _currDetVolName != detVolName) {
    _currDetVolName = detVolName;
    int currPersDetectorIndex = 0;
    for (_currDetector = _currVolStorageFileInfo->detectorsMap.begin();
         _currDetector != _currVolStorageFileInfo->detectorsMap.end(); ++_currDetector, ++currPersDetectorIndex) {
      if (_currDetector->first == _currDetVolName) {
        break;
      }
    }
    if (_currDetector == _currVolStorageFileInfo->detectorsMap.end()) {
      // 2.1 Detector not found. Create it both in index and persistence array
      _currVolStorageFileInfo->detectorsMap.push_back(
          std::make_pair<std::string, HitVolMap>(std::string(_currDetVolName), HitVolMap()));
      _currDetector = _currVolStorageFileInfo->detectorsMap.end(); // Past-the-end iterator
      --_currDetector;                                             // Set the iterator to the last element
      if (_currVolStorageFileInfo->detectorsArray == nullptr) {
        _currVolStorageFileInfo->detectorsArray = new TClonesArray("GGSTHitDetInfo");
      }
      currPersDetectorIndex = _currVolStorageFileInfo->detectorsArray->GetEntries(); // Redundant
      _currPersDetector = new ((*(_currVolStorageFileInfo->detectorsArray))[currPersDetectorIndex]) GGSTHitDetInfo;
      _currPersDetector->detectorName = _currDetVolName.data();
    } else {
      // 2.2 Detector found. Set the current persistent detector
      _currPersDetector = (GGSTHitDetInfo *)(_currVolStorageFileInfo->detectorsArray->At(currPersDetectorIndex));
    }
  }

  // 3. Insert the volume in the index and in the persistency structure

  // Append volume ID to volume name to avoid placing all the replicated volumes at the same place in the index
  std::string volAndId = volume->GetName() + std::to_string(id);

  auto insertResult =
      _currDetector->second.emplace(VolumeKey(volAndId, volume, position), _currPersDetector->volumes.GetEntries());
  if (insertResult.second) {
    // 3.1 No hash collision: this is a new volume. Add it also to the volumes array in detector object...
    GGSTHitVolInfo *volInfo =
        new ((_currPersDetector->volumes)[_currPersDetector->volumes.GetEntries()]) GGSTHitVolInfo;
    volInfo->volumeName = volume->GetName().data();
    for (int i = 0; i < 3; i++) {
      volInfo->volumePos[i] = position[i] / cm;
    }
    volInfo->id = id;
    // ... then return its position in the array
    return _currPersDetector->volumes.GetEntries() - 1;
  } else {
    // 3.2 A hash collision has happened: the volume is already present. So take its position from the index
    return insertResult.first->second;
  }
}