CbmOnlineParWrite.cxx

Go to the documentation of this file.

/* Copyright (C) 2024-2025 FIAS Frankfurt Institute for Advanced Studies, Frankfurt / Main
   SPDX-License-Identifier: GPL-3.0-only
   Authors: Felix Weiglhofer [committer], Sergei Zharko */
 
#include "CbmOnlineParWrite.h"
 
#include "CbmCaParametersHandler.h"
#include "CbmMcbmUtils.h"
#include "CbmMuchGeoScheme.h"
#include "CbmRecoSetupManager.h"
#include "CbmRunDatabaseContainer.h"
#include "CbmSetup.h"
#include "CbmSinkDummy.h"
#include "CbmSourceDummy.h"
#include "CbmStsAddress.h"
#include "CbmStsParAsic.h"
#include "CbmStsParModule.h"
#include "CbmStsParSensor.h"
#include "CbmStsParSensorCond.h"
#include "CbmStsParSetSensor.h"
#include "CbmStsSetup.h"
#include "CbmTaskStsHitFinderParWrite.h"
#include "CbmTaskTofClusterizerParWrite.h"
#include "CbmTaskTrdHitFinderParWrite.h"
#include "CbmTaskTrdUnpackParWrite.h"
#include "CbmTrdParSetAsic.h"
#include "CbmTrdParSetDigi.h"
 
#include <FairParAsciiFileIo.h>
#include <FairParRootFileIo.h>
#include <FairRootFileSink.h>
#include <FairRunAna.h>
#include <FairRuntimeDb.h>
 
#include <TGeoManager.h>
#include <TObjString.h>
#include <TStopwatch.h>
 
#include <boost/filesystem.hpp>
 
#include <array>
#include <fstream>
#include <iostream>
#include <sstream>
 
using namespace cbm::algo;
 
namespace fs = boost::filesystem;
 
// ---------------------------------------------------------------------------------------------------------------------
//
void CbmOnlineParWrite::AddDetectors()
{
  // Add detectors here
  AddTrd();
  AddTof();  // NOTE: BMON as well
  AddSts();
}
 
// ---------------------------------------------------------------------------------------------------------------------
//
void CbmOnlineParWrite::AddTrd()
{
  // Copied from macro/beamtime/mcbm2022/trd_hitfinder_run.C
  fSetup->SetActive(ECbmModuleId::kTrd, kTRUE);
  fSetup->SetActive(ECbmModuleId::kTrd2d, kTRUE);
 
  // ----- TRD digitisation parameters -------------------------------------
  TString geoTagTrd;
  if (!fSetup->IsActive(ECbmModuleId::kTrd)) {
    throw std::runtime_error{"TRD not enabled in current geometry!"};
  }
 
  if (!fSetup->GetGeoTag(ECbmModuleId::kTrd, geoTagTrd)) {
    throw std::runtime_error{"Failed to get TRD geo tag from CbmSetup!"};
  }
 
  TString paramFilesTrd(Form("%s/parameters/trd/trd_%s", fSrcDir.Data(), geoTagTrd.Data()));
  std::vector<TString> paramFilesVecTrd = {"asic", "digi", "gas", "gain"};
  for (auto parIt : paramFilesVecTrd) {
    fParList->Add(new TObjString(Form("%s.%s.par", paramFilesTrd.Data(), parIt.Data())));
  }
 
  for (auto parFileVecIt : *fParList) {
    std::cout << Form("TrdParams - %s - added to parameter file list", parFileVecIt->GetName()) << std::endl;
  }
 
  // ----- TRD task ---------------------------------------------------------
  auto* trdHitfinderPar = new CbmTaskTrdHitFinderParWrite{};
  trdHitfinderPar->SetRecoParOutputDir(fOutDirReco.string());
  fRun->AddTask(trdHitfinderPar);
 
  fParFiles.timesliceReco->trd = std::make_optional(RecoParFiles::Trd{
    .hitfinder = "TrdHitfinderPar.yaml",
  });
 
  fParFiles.timesliceReco->trd2d = std::make_optional(RecoParFiles::Trd2d{
    .hitfinder = "TrdHitfinder2DPar.yaml",
  });
 
  // Initialize input files
  // FairParAsciiFileIo asciiInput;
  // std::string digiparfile = Form("%s/parameters/trd/trd_%s.digi.par", fSrcDir.Data(), geoTagTrd.Data());
  // std::string asicparfile = Form("%s/parameters/trd/trd_%s.asic.par", fSrcDir.Data(), geoTagTrd.Data());
 
  // // Read the .digi file and store result
  // auto* digiparset = new CbmTrdParSetDigi{};
  // if (asciiInput.open(digiparfile.data())) {
  //   digiparset->init(&asciiInput);
  // }
  // asciiInput.close();
 
  // // Read the .asic file and store result
  // auto* asicparset = new CbmTrdParSetAsic{};
  // if (asciiInput.open(asicparfile.data())) {
  //   asicparset->init(&asciiInput);
  // }
  // asciiInput.close();
 
  // // Unpack
  // CbmTaskTrdUnpackParWrite::Pars parFilesUnpack{.asic = asicparset, .digi = digiparset, .runId = fConfig.runId};
  // auto* trdUnpackPar = new CbmTaskTrdUnpackParWrite{parFilesUnpack};
  // fRun->AddTask(trdUnpackPar);
}
 
// ---------------------------------------------------------------------------------------------------------------------
//
void CbmOnlineParWrite::AddTof()
{
  // Copied from macro/tools/tof_hitfinder_run.C
  fSetup->SetActive(ECbmModuleId::kTof, kTRUE);
 
  TString geoTag;
  if (fSetup->IsActive(ECbmModuleId::kTof)) {
    fSetup->GetGeoTag(ECbmModuleId::kTof, geoTag);
    TObjString* tofBdfFile = new TObjString(fSrcDir + "/parameters/tof/tof_" + geoTag + ".digibdf.par");
    fParList->Add(tofBdfFile);
    std::cout << "-I- TOF: Using parameter file " << tofBdfFile->GetString() << std::endl;
  }
 
  // -----   TOF defaults ------------------------
  Int_t calMode      = 93;
  Int_t calSel       = 1;
  Double_t dDeadtime = 50.;
 
  TString TofFileFolder = fSrcDir + "/parameters/mcbm/";
  bool doCalibration    = true;
  TString cCalId        = "490.100.5.0";
  Int_t iCalSet         = 30040500;  // calibration settings
  TString cFname        = "";
 
  // Select calibration setting (NOTE: legacy support)
  if (fConfig.runId >= 2350 && fConfig.runId < 2724) {  // mCBM2022
    iCalSet = 22002500;
  }
  else if (fConfig.runId >= 2724 && fConfig.runId < 2918) {  // mCBM2024_03
    iCalSet = 012032500;
  }
 
  std::string tofCalTag = cbm::RunDatabaseContainer::Instance().GetTofCalibrationTag();
  cFname                = Form("%s/%s1_TofCal.hst.root", TofFileFolder.Data(), tofCalTag.c_str());
 
  if (cFname.IsNull() && doCalibration) {
    cFname =
      Form("%s/%s_set%09d_%02d_%01dtofClust.hst.root", TofFileFolder.Data(), cCalId.Data(), iCalSet, calMode, calSel);
  }
 
  auto* tofCluster = new CbmTaskTofClusterizerParWrite("Task TOF Clusterizer", 0, 1);
  tofCluster->SetCalParFileName(cFname);
  tofCluster->SetCalMode(calMode);
  tofCluster->SetTotMax(20.);                 // Tot upper limit for walk corection
  tofCluster->SetTotMin(0.);                  //(12000.);  // Tot lower limit for walk correction
  tofCluster->SetTotMean(5.);                 // Tot calibration target value in ns
  tofCluster->SetMaxTimeDist(1.0);            // default cluster range in ns
  tofCluster->SetChannelDeadtime(dDeadtime);  // artificial deadtime in ns
  tofCluster->PosYMaxScal(0.75);              //in % of length
  tofCluster->SetRunId(fConfig.runId);
  tofCluster->SetRecoParOutputDir(fOutDirReco.string());
  fRun->AddTask(tofCluster);
 
  fParFiles.timesliceReco->tof = std::make_optional(RecoParFiles::Tof{
    .calibrate = "TofCalibratePar.yaml",
    .hitfinder = "TofHitfinderPar.yaml",
  });
 
  fParFiles.timesliceReco->bmon = std::make_optional(RecoParFiles::Bmon{
    .calibrate = "BmonCalibratePar.yaml",
    .hitfinder = "BmonHitfinderPar.yaml",
  });
}
 
// ---------------------------------------------------------------------------------------------------------------------
//
void CbmOnlineParWrite::AddSts()
{
  // Copied from macro/beamtime/mcbm2022/mcbm_reco.C
 
  auto* recoSts = new CbmTaskStsHitFinderParWrite{};
  // recoSts->SetMode(ECbmRecoMode::EventByEvent);
 
  // recoSts->SetTimeCutDigisAbs(20.0);     // cluster finder: time cut in ns
  // recoSts->SetTimeCutClustersAbs(20.0);  // hit finder: time cut in ns
 
  TString geoTag;
  fSetup->GetGeoTag(ECbmModuleId::kTof, geoTag);
  LOG(info) << "STS geo tag: " << geoTag;
 
  // Sensor params
  CbmStsParSensor sensor6cm(CbmStsSensorClass::kDssdStereo);
  sensor6cm.SetPar(0, 6.2092);  // Extension in x
  sensor6cm.SetPar(1, 6.2);     // Extension in y
  sensor6cm.SetPar(2, 0.03);    // Extension in z
  sensor6cm.SetPar(3, 5.9692);  // Active size in y
  sensor6cm.SetPar(4, 1024.);   // Number of strips front side
  sensor6cm.SetPar(5, 1024.);   // Number of strips back side
  sensor6cm.SetPar(6, 0.0058);  // Strip pitch front side
  sensor6cm.SetPar(7, 0.0058);  // Strip pitch back side
  sensor6cm.SetPar(8, 0.0);     // Stereo angle front side
  sensor6cm.SetPar(9, 7.5);     // Stereo angle back side
 
  CbmStsParSensor sensor12cm(sensor6cm);  // copy all parameters, change then only the y size
  sensor12cm.SetPar(1, 12.4);             // Extension in y
  sensor12cm.SetPar(3, 12.1692);          // Active size in y
 
  // --- Addresses for sensors
  // --- They are defined in each station as sensor 1, module 1, halfladderD (2), ladder 1
  //  Int_t GetAddress(UInt_t unit = 0, UInt_t ladder = 0, UInt_t halfladder = 0, UInt_t module = 0, UInt_t sensor = 0,
  //                   UInt_t side = 0, UInt_t version = kCurrentVersion);
 
  // --- Now we can define the sensor parameter set and tell recoSts to use it
  auto sensorParSet = new CbmStsParSetSensor("CbmStsParSetSensor", "STS sensor parameters"
                                                                   "mcbm2021");
 
  // TODO: is it possible to read these values from a parameter file?
  if (fConfig.runId < 2918) {
 
    Int_t stsAddress01 = CbmStsAddress::GetAddress(0, 0, 1, 0, 0, 0);  // U0 L0 M0  6 cm
    Int_t stsAddress02 = CbmStsAddress::GetAddress(0, 0, 1, 1, 0, 0);  // U0 L0 M1  6 cm
    Int_t stsAddress03 = CbmStsAddress::GetAddress(0, 1, 1, 0, 0, 0);  // U0 L1 M0  6 cm
    Int_t stsAddress04 = CbmStsAddress::GetAddress(0, 1, 1, 1, 0, 0);  // U0 L1 M1  6 cm
    Int_t stsAddress05 = CbmStsAddress::GetAddress(1, 0, 1, 0, 0, 0);  // U1 L0 M0  6 cm
    Int_t stsAddress06 = CbmStsAddress::GetAddress(1, 0, 1, 1, 0, 0);  // U1 L0 M1 12 cm
    Int_t stsAddress07 = CbmStsAddress::GetAddress(1, 1, 1, 0, 0, 0);  // U1 L1 M0  6 cm
    Int_t stsAddress08 = CbmStsAddress::GetAddress(1, 1, 1, 1, 0, 0);  // U1 L1 M1 12 cm
    Int_t stsAddress09 = CbmStsAddress::GetAddress(1, 2, 1, 0, 0, 0);  // U1 L2 M0  6 cm
    Int_t stsAddress10 = CbmStsAddress::GetAddress(1, 2, 1, 1, 0, 0);  // U1 L2 M1  6 cm
    Int_t stsAddress11 = CbmStsAddress::GetAddress(1, 2, 1, 2, 0, 0);  // U1 L2 M2  6 cm
 
    LOG(info) << "STS address01 " << CbmStsAddress::ToString(stsAddress01);
    LOG(info) << "STS address02 " << CbmStsAddress::ToString(stsAddress02);
    LOG(info) << "STS address03 " << CbmStsAddress::ToString(stsAddress03);
    LOG(info) << "STS address04 " << CbmStsAddress::ToString(stsAddress04);
    LOG(info) << "STS address05 " << CbmStsAddress::ToString(stsAddress05);
    LOG(info) << "STS address06 " << CbmStsAddress::ToString(stsAddress06);
    LOG(info) << "STS address07 " << CbmStsAddress::ToString(stsAddress07);
    LOG(info) << "STS address08 " << CbmStsAddress::ToString(stsAddress08);
    LOG(info) << "STS address09 " << CbmStsAddress::ToString(stsAddress09);
    LOG(info) << "STS address10 " << CbmStsAddress::ToString(stsAddress10);
    LOG(info) << "STS address11 " << CbmStsAddress::ToString(stsAddress11);
 
    sensorParSet->SetParSensor(stsAddress01, sensor6cm);
    sensorParSet->SetParSensor(stsAddress02, sensor6cm);
    sensorParSet->SetParSensor(stsAddress03, sensor6cm);
    sensorParSet->SetParSensor(stsAddress04, sensor6cm);
    sensorParSet->SetParSensor(stsAddress05, sensor6cm);
    sensorParSet->SetParSensor(stsAddress06, sensor12cm);
    sensorParSet->SetParSensor(stsAddress07, sensor6cm);
    sensorParSet->SetParSensor(stsAddress08, sensor12cm);
    sensorParSet->SetParSensor(stsAddress09, sensor6cm);
    sensorParSet->SetParSensor(stsAddress10, sensor6cm);
    sensorParSet->SetParSensor(stsAddress11, sensor6cm);
  }
  else {
    uint32_t addr01 = 0x10008012;
    uint32_t addr02 = 0x10018012;
    uint32_t addr03 = 0x10008412;
    uint32_t addr04 = 0x10018412;
    uint32_t addr05 = 0x10008422;
    uint32_t addr06 = 0x10018422;
    uint32_t addr07 = 0x10008822;
    uint32_t addr08 = 0x10018822;
    uint32_t addr09 = 0x10028822;
    uint32_t addr10 = 0x10008022;
    uint32_t addr11 = 0x10018022;
    uint32_t addr00 = 0x10000002;  // New station 0 in mCBM2024
 
    LOG(info) << "STS address01 " << CbmStsAddress::ToString(addr01);
    LOG(info) << "STS address02 " << CbmStsAddress::ToString(addr02);
    LOG(info) << "STS address03 " << CbmStsAddress::ToString(addr03);
    LOG(info) << "STS address04 " << CbmStsAddress::ToString(addr04);
    LOG(info) << "STS address05 " << CbmStsAddress::ToString(addr05);
    LOG(info) << "STS address06 " << CbmStsAddress::ToString(addr06);
    LOG(info) << "STS address07 " << CbmStsAddress::ToString(addr07);
    LOG(info) << "STS address08 " << CbmStsAddress::ToString(addr08);
    LOG(info) << "STS address09 " << CbmStsAddress::ToString(addr09);
    LOG(info) << "STS address10 " << CbmStsAddress::ToString(addr10);
    LOG(info) << "STS address11 " << CbmStsAddress::ToString(addr11);
    LOG(info) << "STS address00 " << CbmStsAddress::ToString(addr00);
 
    sensorParSet->SetParSensor(addr01, sensor6cm);
    sensorParSet->SetParSensor(addr02, sensor6cm);
    sensorParSet->SetParSensor(addr03, sensor6cm);
    sensorParSet->SetParSensor(addr04, sensor6cm);
    sensorParSet->SetParSensor(addr05, sensor6cm);
    sensorParSet->SetParSensor(addr06, sensor12cm);
    sensorParSet->SetParSensor(addr07, sensor6cm);
    sensorParSet->SetParSensor(addr08, sensor6cm);
    sensorParSet->SetParSensor(addr09, sensor6cm);
    sensorParSet->SetParSensor(addr10, sensor6cm);
    sensorParSet->SetParSensor(addr11, sensor12cm);
    sensorParSet->SetParSensor(addr00, sensor6cm);
  }
  //else {
  //  throw std::runtime_error("STS: Unknown setup type");
  //}
 
  recoSts->UseSensorParSet(sensorParSet);
 
  // ASIC params: #ADC channels, dyn. range, threshold, time resol., dead time,
  // noise RMS, zero-threshold crossing rate
  auto parAsic = new CbmStsParAsic(128, 31, 75000., 3000., 5., 800., 1000., 3.9789e-3);
 
  // Module params: number of channels, number of channels per ASIC
  auto parMod = new CbmStsParModule(2048, 128);
  parMod->SetAllAsics(*parAsic);
  recoSts->UseModulePar(parMod);
 
  // Sensor conditions: full depletion voltage, bias voltage, temperature,
  // coupling capacitance, inter-strip capacitance
  auto sensorCond = new CbmStsParSensorCond(70., 140., 268., 17.5, 1.);
  recoSts->UseSensorCond(sensorCond);
  recoSts->SetRecoParOutputDir(fOutDirReco.string());
 
  fRun->AddTask(recoSts);
 
  fParFiles.timesliceReco->sts = std::make_optional(RecoParFiles::Sts{
    .hitfinder = "StsHitfinder.yaml",
  });
}
 
// ---------------------------------------------------------------------------------------------------------------------
//
void CbmOnlineParWrite::AddCa()
{
  using cbm::algo::ca::ParametersIO;
  using cbm::algo::kf::EFieldMode;
  using cbm::algo::kf::FloatTag;
 
  // ----- MuCh digi parameters initialization --------------------------------------
  if (CbmSetup::Instance()->IsActive(ECbmModuleId::kMuch)) {
    // Parameter file name
    TString geoTag;
    CbmSetup::Instance()->GetGeoTag(ECbmModuleId::kMuch, geoTag);
    Int_t muchFlag  = (geoTag.Contains("mcbm") ? 1 : 0);
    TString parFile = gSystem->Getenv("VMCWORKDIR");
    parFile += "/parameters/much/much_" + geoTag(0, 4) + "_digi_sector.root";
 
    // Initialization of the geometry scheme
    auto muchGeoScheme = CbmMuchGeoScheme::Instance();
    if (!muchGeoScheme->IsInitialized()) {
      muchGeoScheme->Init(parFile, muchFlag);
    }
  }
  // --------------------------------------------------------------------------------
 
  std::string caParFilename = Form("%s.ca.par", fGeoSetupTag.Data());
  fs::path mainCfgPath      = fs::absolute(fs::weakly_canonical(fOutDirReco)) / "MainConfig.yaml";
  fs::path caParPath        = fs::absolute(fs::weakly_canonical(fOutDirReco)) / caParFilename;
  if (!fs::exists(mainCfgPath)) {
    throw std::runtime_error("Could not locate MainConfig.yaml");
  }
 
  // Tracking parameters handler
  auto& pCaParHandler = cbm::ca::ParametersHandler::Instance();
  pCaParHandler.SetIgnoreHitPresence(true);
  pCaParHandler.SetMainConfig(mainCfgPath.string());
  pCaParHandler.StoreParameters(caParPath.string());
  pCaParHandler.SetFieldMode(cbm::algo::kf::EFieldMode::Interpolated);
  fRun->AddTask(&pCaParHandler);
 
  // Update the config
  fParFiles.timesliceReco->tracking = std::make_optional(RecoParFiles::Tracking{.capar = caParFilename});
}
 
// ---------------------------------------------------------------------------------------------------------------------
//
void CbmOnlineParWrite::CopyParameters()
{
 
  std::string srcDirName = TString(gSystem->Getenv("VMCWORKDIR")).Data();
  if (srcDirName.empty()) {
    LOG(fatal) << "Cannot copy parameters: the VMCWORKDIR environmental variable is not specified";
  }
 
  const auto& cbmRunDb = cbm::RunDatabaseContainer::Instance();
  fs::path srcDir      = srcDirName + "/" + cbmRunDb.GetRecoParDir();
  if (!fs::exists(srcDir)) {
    LOG(fatal) << "Source directory " << srcDir.string()
               << " for the parameter files does not exist (runId=" << fConfig.runId << "). "
               << "Maybe, the run ID was not defined in the CbmRunDatabase.yaml";
  }
 
  fs::path dstDir = fs::absolute(fs::weakly_canonical(fOutDirReco));
  if (!fs::exists(dstDir)) {
    fs::create_directories(dstDir);
  }
  LOG(info) << "Copying parameters from " << srcDir.string() << " to " << dstDir.string();
  for (fs::directory_iterator it(srcDir); it != fs::directory_iterator{}; ++it) {
    fs::path srcFile = it->path();
    fs::path dstFile = fs::absolute(fs::weakly_canonical(dstDir / srcFile.filename()));
    if (fs::exists(dstFile)) {
      fs::remove(dstFile);
    }
    fs::copy_file(srcFile, dstFile);
    LOG(info) << " - copying " << srcFile.string() << " -> " << dstFile.string();
  }
}
 
// ---------------------------------------------------------------------------------------------------------------------
//
void CbmOnlineParWrite::Run(const Config& config)
{
  // Copied and adjusted from macro/beamtime/mcbm2022/trd_hitfinder_run.C
 
  static bool callOnce = true;
  if (!callOnce) {
    throw std::runtime_error("CbmOnlineParWrite::Run() can only be called once at the moment!");
  }
  callOnce = false;
 
  fConfig = config;
 
  // -----   Environment   --------------------------------------------------
  fSrcDir = gSystem->Getenv("VMCWORKDIR");  // top source directory
  // ------------------------------------------------------------------------
 
  // --- CBM run database -----------
  auto& cbmRunDb = cbm::RunDatabaseContainer::Instance();
  cbmRunDb.LoadRun(fConfig.runId);
 
  // Terminate, if alignment was required, but no alignment was found for the run ID
  if (config.doAlignment && cbmRunDb.GetAlignmentTag().empty()) {
    std::stringstream msg;
    msg << "LOGIC ERROR: alignment application was required, but no alignment tag was defined in the "
        << "data-base. Please, re-run the application with \"--no-alignment\" option";
    throw std::logic_error(msg.str());
  }
 
  // --- Output path ----------------
  fOutDir = fs::absolute(fs::weakly_canonical(fConfig.outputDir));
  if (!fConfig.experimental && !cbmRunDb.GetRecoParTag().empty()) {
    fOutDirReco = fOutDir / cbmRunDb.GetRecoParTag();
    fs::create_directories(fOutDirReco);
  }
  else if (fConfig.experimental) {
    fOutDirReco = fOutDir;
  }
  else {
    LOG(warn) << "A tag for reconstruction parameters was not defined in the CBM run database. The parameters "
              << "will be stored in the output directory, but please provide a valid tag for this run ID";
    fOutDirReco = fOutDir;
  }
 
  TString geoSetupTag = cbm::RunDatabaseContainer::Instance().GetGeoSetupTag();
 
  LOG(info) << "Using geometry setup: " << geoSetupTag;
  fGeoSetupTag = geoSetupTag;
 
  auto geoFilePath = fs::absolute(fs::weakly_canonical(fConfig.geoFileDir)) / (geoSetupTag + ".geo.root").Data();
  if (!fs::exists(geoFilePath)) {
    throw std::runtime_error{"Geometry file " + geoFilePath.string() + " does not exist"};
  }
 
  TString geoFile = geoFilePath.string();
  LOG(info) << "Using geometry file: " << geoFile;
  fSetup = CbmSetup::Instance();
  fSetup->LoadSetup(geoSetupTag);
 
  //-----  Load Parameters --------------------------------------------------
  fParList = new TList();
 
  // -----   FairRunAna   ---------------------------------------------------
  fRun = new FairRunAna();
 
  // Dummy source required for the FairRunAna as it will crash without a source
  fRun->SetSource(new CbmSourceDummy{});
  fRun->SetSink(new CbmSinkDummy{});
 
 
  // =========================================================================
  // ===                   Alignment Correction                            ===
  // =========================================================================
  if (fConfig.doAlignment) {
    TString alignmentTag            = cbm::RunDatabaseContainer::Instance().GetAlignmentTag();
    TString alignmentMatrixFileName = fSrcDir + "/parameters/mcbm/AlignmentMatrices_" + alignmentTag + ".root";
    if (alignmentMatrixFileName.Length() == 0) {
      throw std::runtime_error{"Alignment matrix file name is empty"};
    }
 
    LOG(info) << "Applying alignment for file '" << alignmentMatrixFileName << "'";
 
    // Define the basic structure which needs to be filled with information
    // This structure is stored in the output file and later passed to the
    // FairRoot framework to do the (miss)alignment
    std::map<std::string, TGeoHMatrix>* matrices{nullptr};
 
    // read matrices from disk
    LOG(info) << "Filename: " << alignmentMatrixFileName;
    TFile* misalignmentMatrixRootfile = new TFile(alignmentMatrixFileName, "READ");
    if (misalignmentMatrixRootfile->IsOpen()) {
      gDirectory->GetObject("MisalignMatrices", matrices);
      misalignmentMatrixRootfile->Close();
    }
    else {
      throw std::runtime_error{"Could not open alignment matrix file: " + alignmentMatrixFileName};
    }
 
    if (matrices) {
      fRun->AddAlignmentMatrices(*matrices);
    }
    else {
      throw std::runtime_error{"Could not read alignment matrices from file: " + alignmentMatrixFileName};
    }
  }
 
  // -----   Create ParFiles instance -----------------------------------
  fParFiles.timesliceReco = std::make_optional<cbm::algo::RecoParFiles>();
  fRecoParTopDir = Form("recoPar/%s/", cbm::RunDatabaseContainer::Instance().GetRecoParTag().c_str());
 
 
  // -----   Add detectors   ------------------------------------------------
  if (fConfig.experimental) {
    CopyParameters();
  }
 
  AddDetectors();
 
  fParFiles.eventReco     = fParFiles.timesliceReco;
  fParFiles.recoSetupPath = std::make_optional<std::string>(Form("%s.crs", fGeoSetupTag.Data()));
 
  auto* pRecoSetupManager = cbm::RecoSetupManager::Instance();
  pRecoSetupManager->StoreSetup((fOutDirReco / *(fParFiles.recoSetupPath)).string());
  fRun->AddTask(pRecoSetupManager);
 
  AddCa();
 
  if (fConfig.experimental) {
    UpdateParFiles();
  }
 
  // -----  Parameter database   --------------------------------------------
  FairRuntimeDb* rtdb        = fRun->GetRuntimeDb();
  FairParAsciiFileIo* parIo2 = new FairParAsciiFileIo();
  parIo2->setAutoWritable(false);
  parIo2->open(fParList, "in");
  rtdb->setSecondInput(parIo2);
  // ------------------------------------------------------------------------
 
  // -----   Run initialisation   -------------------------------------------
  fRun->SetGeomFile(geoFile);
  fRun->Init();
 
 
 
  // ----- Clean up ---------------------------------------------------------
  gGeoManager->GetListOfVolumes()->Delete();
  gGeoManager->GetListOfShapes()->Delete();
  delete gGeoManager;
 
  // Delete files created by FairRun that I don't know how disable otherwise
  gSystem->Exec("rm all_*.par *.mat.kf.bin");
}
 
 
// ---------------------------------------------------------------------------------------------------------------------
//
void CbmOnlineParWrite::UpdateParFiles()
{
  namespace fs = boost::filesystem;
 
  fs::path mainCfgPath = fs::absolute(fs::weakly_canonical(fOutDirReco)) / "MainConfig.yaml";
  if (!fs::exists(mainCfgPath)) {
    throw std::runtime_error("Cannot find the main config " + mainCfgPath.string() + " in the output directory ");
  }
 
  // The following approach is needed to preserve the comments and formatting in the MainConfig.yaml for the
  // rest of the nodes. It reads the file in three streams:
  //  iStream == 0 (head): everything before the parFiles node
  //  iStream == 1 (body): the parFiles node itself, which will be parsed and rewritten
  //  iStream == 2 (tail): everything after the parFiles node
  std::array<std::stringstream, 3> content;
  {
    int iStream{0};
    int nFound{0};  // number of lines, which contain "parFiles:" entry in order to catch a potential error
    std::ifstream in(mainCfgPath.string());
    for (std::string line; std::getline(in, line, '\n');) {
      bool bFound = line.find("parFiles:") != std::string::npos;
      if (bFound) ++nFound;
      if (nFound > 1) {
        throw std::runtime_error("CbmOnlineParWrite::UpdateParFiles(): found more than one lines with \"parFiles:\"");
      }
      if (iStream == 1 && line[0] != ' ' && line[0] != '\t') iStream = 2;
      if (iStream == 0 && bFound) iStream = 1;
      content[iStream] << line << '\n';
    }
  }
 
  // Read the ParFiles node and update it with timesliceReco, eventReco and recoSetupPath
  auto parFilesNode      = YAML::Load(content[1])["parFiles"];
  auto parFiles          = yaml::Read<ParFiles>(parFilesNode);
  parFiles.timesliceReco = fParFiles.timesliceReco;
  parFiles.eventReco     = fParFiles.eventReco;
  parFiles.recoSetupPath = fParFiles.recoSetupPath;
 
  // Store the updated config
  std::ofstream out(mainCfgPath.string());
  out << content[0].str() << '\n';
  {
    cbm::util::yaml::Dump dump;
    YAML::Emitter yamlEmitter;
    yamlEmitter << YAML::BeginMap;
    yamlEmitter << YAML::Key << std::string{"parFiles"};
    yamlEmitter << YAML::Value;
    dump.DoDump(parFiles, yamlEmitter);
    yamlEmitter << YAML::EndMap;
    out << yamlEmitter.c_str() << '\n';
  }
  out << content[2].str();
  out.close();
}