KfSetup.h

Go to the documentation of this file.

/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
   SPDX-License-Identifier: GPL-3.0-only
   Authors: Sergei Zharko [committer] */

 
#pragma once  // include this header only once per compilation unit
 
#include "KfActiveLayer.h"
#include "KfDefs.h"
#include "KfField.h"
#include "KfMaterialMap.h"
#include "KfModuleIndexMap.h"
#include "KfTarget.h"
#include "KfVector.h"
 
#include <boost/serialization/access.hpp>
 
#include <fstream>
#include <numeric>
#include <optional>
#include <sstream>
#include <string>
#include <unordered_map>
#include <unordered_set>
#include <vector>
 
namespace cbm::algo::kf
{
  template<typename F>
  class alignas(VcMemAlign) Setup {
    template<typename>
    friend class Setup;
 
   public:
    // FIXME: hide to private
    Setup() = default;

    Setup(const ModuleIndexMap& moduleIndexMap,            //
          const std::vector<MaterialMap>& materialLayers,  //
          std::vector<ActiveLayer<F>> activeLayers,        //
          const Field<F>& field,                           //
          const Target<F>& target)
      : fModuleIndexMap(moduleIndexMap)
      , fvMaterialLayers(materialLayers)
      , fvActiveLayers(std::move(activeLayers))
      , fField(field)
      , fTarget(target)
    {
      if (!Validate()) {
        throw std::logic_error("Attempt to create a kf::Setup from inconsistent set of parameters. Please, use the "
                               "kf::SetupBuilder to create a valid setup.");
      }
    }

    template<typename FSrc>
    Setup(const Setup<FSrc>& other)
      : fModuleIndexMap(other.fModuleIndexMap)
      , fvMaterialLayers(other.fvMaterialLayers)
      , fField(std::move(Field<F>::CopyConvert(other.fField)))
      , fTarget(other.fTarget)
    {
      fvActiveLayers.reserve(other.fvActiveLayers.size());
      for (const auto& layer : other.fvActiveLayers) {
        fvActiveLayers.emplace_back(layer);
      }
    }

    Setup(const Setup&) = default;

    Setup(Setup&&) = default;

    ~Setup() = default;

    Setup& operator=(const Setup& other) = default;

    Setup& operator=(Setup&&) = default;

    template<class EDetID>
    Setup CreateSetupWithDisabledLayers(const std::vector<std::pair<EDetID, int>>& toDisable) const;

    template<class EDetID>
    void DisableLayer(EDetID iDet, int iLoc);

    const ActiveLayer<F>& GetActiveLayer(int iLayer) const { return fvActiveLayers[iLayer]; }

    const auto& GetActiveLayers() const { return fvActiveLayers; }

    template<class EDetID>
    const ActiveLayer<F>& GetActiveLayer(EDetID iDet, int iLoc) const
    {
      return fvActiveLayers[fModuleIndexMap.LocalToGlobal(iDet, iLoc)];
    }

    const Field<F>& GetField() const { return fField; }

    const MaterialMap& GetMaterial(int iLayer) const { return fvMaterialLayers[iLayer]; }

    template<class EDetID>
    const MaterialMap& GetMaterial(EDetID iDet, int iLoc) const
    {
      return fvMaterialLayers[fModuleIndexMap.LocalToGlobal(iDet, iLoc)];
    }

    const ModuleIndexMap& GetIndexMap() const { return fModuleIndexMap; }

    int GetNofLayers() const { return fModuleIndexMap.GetNofLayers(); }

    template<class EDetID>
    int GetNofLayers(EDetID iDet) const
    {
      return fModuleIndexMap.GetNofLayers(iDet);
    }

    const Target<F>& GetTarget() const { return fTarget; }

    std::string ToString(int verbosity = 1, int indentLevel = 0) const;
 
   private:
    ModuleIndexMap fModuleIndexMap{};              
    std::vector<MaterialMap> fvMaterialLayers{};   
    std::vector<ActiveLayer<F>> fvActiveLayers{};  
    Field<F> fField{};                             
    Target<F> fTarget{};                           
 

    bool Validate() const;

    friend class boost::serialization::access;
    template<class Archive>
    void serialize(Archive& ar, const unsigned int /*version*/)
    {
      ar& fModuleIndexMap;
      ar& fvMaterialLayers;
      ar& fvActiveLayers;
      ar& fTarget;
      ar& fField;
    }
  };
 
 
  // -------------------------------------------------------------------------------------------------------------------
  //
  template<typename F>
  template<class EDetID>
  Setup<F> Setup<F>::CreateSetupWithDisabledLayers(const std::vector<std::pair<EDetID, int>>& toDisableInp) const
  {
    Setup<F> res(*this);
    // We also take into account a situation, when some of the elements in the input vector make no sense in terms of
    // the geometry setup, so we create a new container for keeping the inactive layers of only existing stations and
    // sort it in the direction of geoId decrease:
 
    using Inp_t = typename std::remove_cv_t<std::remove_reference_t<decltype(toDisableInp)>>;
 
    // Build a temporary array, which contains only valid {iDet, iLoc} to be disabled
    Inp_t sorted;
    {
      // Collect the {iDet: locIds} map: actual iLoc for each iDet to be disabled
      // This approach addresses the following:
      //   1) filters {iDet, iLoc(> -1)} pairs, which are passed with toDisable vector, which do not exist in the setup
      //   2) translates iLoc == -1 into a sequence of iLoc [0, NofLayers(detId)), which are to be disabled
      std::unordered_map<EDetID, std::unordered_set<int>> filtered;
      // Lookup the input toDisable, fill the toDisableTmp map:
      for (auto [iDet, iLoc] : toDisableInp) {
        // Disable a single layer (with locId > -1), if it exists in this setup
        if (iLoc > -1 && fModuleIndexMap.LocalToGlobal(iDet, iLoc) > -1) {
          filtered[iDet].insert(iLoc);
        }
 
        // Disable all layers of the detector subsystem, it at has at least one layer
        if (iLoc == -1 && fModuleIndexMap.GetNofLayers(iDet) > 0) {
          auto& locIds = filtered[iDet];
          for (int iL = 0; iL < fModuleIndexMap.GetNofLayers(iDet); ++iL) {
            locIds.insert(iL);
          }
        }
      }
 
      // Fill the temporary array
      sorted.reserve(std::accumulate(filtered.begin(), filtered.end(), 0,
                                     [](int sum, const auto& item) { return sum + item.second.size(); }));
      for (const auto& [iDet, locIds] : filtered) {
        // NOTE: [iDet = iDet] is a solution to pass structure bindings in lambdas for some compilers
        std::transform(locIds.begin(), locIds.end(), std::back_inserter(sorted), [iDet = iDet](int iLoc) {
          return std::pair{iDet, iLoc};
        });
      }
    }
 
    // Sort the prepared array in the decreasing order of global Ids
    std::sort(sorted.begin(), sorted.end(), [&](auto l, auto r) {
      return fModuleIndexMap.LocalToGlobal(l.first, l.second) > fModuleIndexMap.LocalToGlobal(r.first, r.second);
    });
 
    // Disable the selected layers
    for (auto [iDet, iLoc] : sorted)
      res.DisableLayer(iDet, iLoc);
    return res;
  }
 
  // -------------------------------------------------------------------------------------------------------------------
  //
  template<typename F>
  template<class EDetID>
  void Setup<F>::DisableLayer(EDetID iDet, int iLoc)
  {
    int iLayer{fModuleIndexMap.LocalToGlobal(iDet, iLoc)};
    if (iLayer < 0) {
      return;
    }
 
    // Remove material layer and add it to the next one
    if (iLayer < static_cast<int>(fvMaterialLayers.size() - 1)) {
      fvMaterialLayers[iLayer + 1].Add(fvMaterialLayers[iLayer], utils::simd::Cast<F, float>(fTarget.GetZ()));
    }
    fvMaterialLayers.erase(fvMaterialLayers.begin() + iLayer);
 
    // Remove active layer
    fvActiveLayers.erase(fvActiveLayers.begin() + iLayer);
 
    // Remove field slice
    fField.RemoveSlice(iLayer);
 
    // Disable layer in the index map
    fModuleIndexMap.Disable(iDet, iLoc);
  }
 
 
}  // namespace cbm::algo::kf