L1InitManager.cxx

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

/************************************************************************************************************
 * @file L1InitManager.cxx
 * @bried Input data management class for L1Algo
 * @since 19.01.2022
 ***********************************************************************************************************/

#include "L1InitManager.h"

#include <algorithm>
#include <sstream>

#include "L1Assert.h"

//-----------------------------------------------------------------------------------------------------------------------
//
void L1InitManager::AddStation(const L1BaseStationInfo& inStation)
{
  // Check if other fields were defined already
  // Active detector IDs
  L1MASSERT(0, fInitController.GetFlag(EInitKey::kActiveDetectorIDs),
            "Attempt to add a station info before the active detetors set had been initialized");

  // Number of stations check
  L1MASSERT(0, fInitController.GetFlag(EInitKey::kStationsNumberCrosscheck),
            "Attempt to add a station info before the numbers of stations for each detector had been initialized");

  // Field function
  L1MASSERT(0, fInitController.GetFlag(EInitKey::kFieldFunction),
            "Attempt to add a station info before the magnetic field function had been intialized");

  // Check activeness of this station type
  bool isStationActive =
    inStation.GetTrackingStatus() && fActiveDetectorIDs.find(inStation.GetDetectorID()) != fActiveDetectorIDs.end();
  if (isStationActive) {
    // initialize magnetic field slice
    L1BaseStationInfo inStationCopy = L1BaseStationInfo(inStation);  // make a copy of station so it can be initialized
    inStationCopy.SetFieldFunction(fFieldFunction);

    // check, if material map is used
    if (!inStationCopy.GetInitController().GetFlag(L1BaseStationInfo::EInitKey::kThicknessMap)) {
      LOG(warn) << "Station material map was not set for detectorID = "
                << static_cast<int>(inStationCopy.GetDetectorID()) << ", stationID = " << inStationCopy.GetStationID()
                << ". Homogenious material budget will be used: " << inStationCopy.GetRadThick()[0];
      L1Material material;
      material.SetBins(1, 100);
      material.SetRadThick(0, 0, inStationCopy.GetRadThick()[0]);
      inStationCopy.SetMaterialMap(std::move(material));
    }

    // Check station init
    LOG(debug) << "L1InitManager::AddStation:(original) L1BaseStationInfo " << inStation.GetInitController().ToString();
    LOG(debug) << "L1InitManager::AddStation:(copy)     L1BaseStationInfo " << inStation.GetInitController().ToString();
    if (!inStationCopy.GetInitController().IsFinalized()) {
      LOG(fatal) << "Attempt to add an incompletely initialized station info object (detectorID = "
                 << static_cast<int>(inStationCopy.GetDetectorID()) << ", stationID = " << inStationCopy.GetStationID()
                 << ")";
    }
    // insert the station in a set
    auto insertionResult = fStationsInfo.insert(std::move(inStationCopy));
    if (!insertionResult.second) {
      LOG(fatal) << "Attempt to add a dublicating station info object (detectorID = "
                 << static_cast<int>(inStationCopy.GetDetectorID()) << ", stationID = " << inStationCopy.GetStationID()
                 << ")";
    }

    int index = fStationsInfo.size() - 1
                + (fNstationsGeometry[fNstationsGeometry.size() - 1] - fNstationsActive[fNstationsActive.size() - 1]);
    fParameters.fActiveStationGlobalIDs[index] = fStationsInfo.size() - 1;
  }
  else {
    int index = fStationsInfo.size()
                + (fNstationsGeometry[fNstationsGeometry.size() - 1] - fNstationsActive[fNstationsActive.size() - 1]);
    fParameters.fActiveStationGlobalIDs[index] = -1;
    fNstationsActive[static_cast<L1DetectorID_t>(inStation.GetDetectorID())]--;
    fNstationsActive[fNstationsActive.size() - 1]--;
  }
  LOG(debug) << "L1InitManager: adding a station with stationID = " << inStation.GetStationID()
             << " and detectorID = " << static_cast<int>(inStation.GetDetectorID())
             << ". Is active: " << isStationActive;
}

//-----------------------------------------------------------------------------------------------------------------------
//
void L1InitManager::CheckInit()
{
  this->CheckCAIterationsInit();
  this->CheckStationsInfoInit();
}

//-----------------------------------------------------------------------------------------------------------------------
// NOTE: this function should be called once in the TransferParametersContainer
void L1InitManager::FormParametersContainer()
{
  // Check initialization
  this->CheckInit();

  if (!fInitController.IsFinalized()) {
    LOG(fatal) << "Attempt to form parameters container before all necessary fields were initialized"
               << fInitController.ToString();
  }

  // Copy station numbers
  std::copy(fNstationsActive.begin(), fNstationsActive.end(), fParameters.fNstationsActive.begin());
  std::copy(fNstationsGeometry.begin(), fNstationsGeometry.end(), fParameters.fNstationsGeometry.begin());

  // Form array of stations
  auto destinationArrayIterator = fParameters.fStations.begin();
  for (const auto& item : fStationsInfo) {
    *destinationArrayIterator = item.GetL1Station();
    ++destinationArrayIterator;
  }

  // Form array of material map
  auto thickMapIt = fParameters.fThickMap.begin();
  for (auto it = fStationsInfo.begin(); it != fStationsInfo.end(); ++it) {
    auto node   = fStationsInfo.extract(it);
    *thickMapIt = std::move(node.value().TakeMaterialMap());
    fStationsInfo.insert(std::move(node));
    ++thickMapIt;
  }

  // Check the consistency of the parameters object. If object inconsistent, it throws std::logic_error
  fParameters.CheckConsistency();
}

//-----------------------------------------------------------------------------------------------------------------------
//
void L1InitManager::InitTargetField(double zStep)
{
  if (fInitController.GetFlag(EInitKey::kPrimaryVertexField)) {
    LOG(warn) << "L1InitManager::InitTargetField: attempt to reinitialize the field value and field region "
              << "near target. Ignore";
    return;
  }

  // Check for field function
  L1MASSERT(0, fInitController.GetFlag(EInitKey::kFieldFunction),
            "Attempt to initialze the field value and field region near target before initializing field function");

  // Check for target defined
  L1MASSERT(
    0, fInitController.GetFlag(EInitKey::kTargetPos),
    "Attempt to initialize the field value and field region near target before the target position initialization");

  constexpr int nDimensions {3};
  constexpr int nPointsNodal {3};

  std::array<double, nPointsNodal> inputNodalZ {fTargetZ, fTargetZ + zStep, fTargetZ + 2. * zStep};
  std::array<L1FieldValue, nPointsNodal> B {};
  std::array<fvec, nPointsNodal> z {};
  // loop over nodal points
  for (int idx = 0; idx < nPointsNodal; ++idx) {
    double point[nDimensions] {0., 0., inputNodalZ[idx]};
    double field[nDimensions] {};
    fFieldFunction(point, field);
    z[idx]   = inputNodalZ[idx];
    B[idx].x = field[0];
    B[idx].y = field[1];
    B[idx].z = field[2];
  }  // loop over nodal points: end
  fParameters.fVertexFieldRegion.Set(B[0], z[0], B[1], z[1], B[2], z[2]);
  fParameters.fVertexFieldValue = B[0];

  fInitController.SetFlag(EInitKey::kPrimaryVertexField);
}

//-----------------------------------------------------------------------------------------------------------------------
//
void L1InitManager::PrintCAIterations(int verbosityLevel) const
{
  for (const auto& iteration : fParameters.fCAIterations) {
    iteration.Print(verbosityLevel);
  }
}

//-----------------------------------------------------------------------------------------------------------------------
//
void L1InitManager::PrintStations(int verbosityLevel) const
{
  for (const auto& station : fStationsInfo) {
    station.Print(verbosityLevel);
  }
}

//-----------------------------------------------------------------------------------------------------------------------
//
void L1InitManager::PushBackCAIteration(const L1CAIteration& iteration)
{
  // TODO: probably some checks must be inserted here (S.Zharko)
  bool control = fInitController.GetFlag(EInitKey::kCAIterationsNumberCrosscheck);
  //std::cout << "L1InitManager::PushBackCAIteration " << control << '\n';
  L1MASSERT(0, control,  //fInitController.GetFlag(EInitKey::kCAIterationsNumberCrosscheck),
            "Attempt to push back a CA track finder iteration before the number of iterations was defined");

  fParameters.fCAIterations.push_back(iteration);
}

//-----------------------------------------------------------------------------------------------------------------------
//
void L1InitManager::SetActiveDetectorIDs(const L1DetectorIDSet_t& detectorIDs)
{
  // TODO: To think about redifinition possibilities: should it be allowed or not? (S.Zh.)
  fActiveDetectorIDs = detectorIDs;
  fInitController.SetFlag(EInitKey::kActiveDetectorIDs);
}

//-----------------------------------------------------------------------------------------------------------------------
//
void L1InitManager::SetCAIterationsNumberCrosscheck(int nIterations)
{
  fCAIterationsNumberCrosscheck                = nIterations;
  L1Vector<L1CAIteration>& iterationsContainer = fParameters.fCAIterations;

  // NOTE: should be called to prevent multiple copyings of objects between the memory realocations
  iterationsContainer.reserve(nIterations);
  fInitController.SetFlag(EInitKey::kCAIterationsNumberCrosscheck);
}

//-----------------------------------------------------------------------------------------------------------------------
//
void L1InitManager::SetFieldFunction(const L1FieldFunction_t& fieldFunction)
{
  if (!fInitController.GetFlag(EInitKey::kFieldFunction)) {
    fFieldFunction = fieldFunction;
    fInitController.SetFlag(EInitKey::kFieldFunction);
  }
  else {
    LOG(warn) << "L1InitManager::SetFieldFunction: attempt to reinitialize the field function. Ignored";
  }
}

//-----------------------------------------------------------------------------------------------------------------------
//
void L1InitManager::SetGhostSuppression(int ghostSuppression)
{
  if (fInitController.GetFlag(EInitKey::kGhostSuppression)) {
    LOG(warn) << "L1InitManager::SetGhostSuppression: attempt of reinitializating the ghost suppresion flag. Ignore";
    return;
  }
  fGhostSuppression = ghostSuppression;
  fInitController.SetFlag(EInitKey::kGhostSuppression);
}

//-----------------------------------------------------------------------------------------------------------------------
//
void L1InitManager::SetMomentumCutOff(float momentumCutOff)
{
  if (fInitController.GetFlag(EInitKey::kMomentumCutOff)) {
    LOG(warn) << "L1InitManager::SetMomentumCutOff: attempt of reinitializating the momentum cutoff value. Ignore";
    return;
  }
  fMomentumCutOff = momentumCutOff;
  fInitController.SetFlag(EInitKey::kMomentumCutOff);
}

//-----------------------------------------------------------------------------------------------------------------------
//
void L1InitManager::SetNstations(L1DetectorID detectorID, int nStations)
{
  L1MASSERT(0, fInitController.GetFlag(EInitKey::kActiveDetectorIDs),
            "Attempt to set crosscheck number of stations before the active detetors set had been initialized");

  // NOTE: We add and check only those detectors which will be active (?)
  // For INACTIVE detectors the initialization code for it inside CbmL1/BmnL1 can (and must) be still in,
  // but it will be ignored inside L1InitManager.
  if (fActiveDetectorIDs.find(detectorID) != fActiveDetectorIDs.end()) {
    if (nStations) {
      fNstationsGeometry[static_cast<L1DetectorID_t>(detectorID)] = nStations;
      fNstationsActive[static_cast<L1DetectorID_t>(detectorID)]   = nStations;
    }
    else {
      // TODO: Probably it is better to replace fatal with warn and remove the detectorID from active detectors (S.Zharko)
      LOG(fatal) << "L1InitManager::SetNstations: attempt to initialize zero stations for active detector: "
                 << static_cast<L1DetectorID_t>(detectorID);
    }
  }

  // Check if all the station numbers for active detectors are initialized now:
  LOG(debug) << "L1InitManager::SetNstations called for detectorID = " << static_cast<int>(detectorID);
  if (!fInitController.GetFlag(EInitKey::kStationsNumberCrosscheck)) {
    bool ifInitialized = true;
    for (auto item : fActiveDetectorIDs) {
      if (fNstationsGeometry[static_cast<L1DetectorID_t>(item)] == 0) {
        ifInitialized = false;
        break;
      }
    }
    fInitController.SetFlag(EInitKey::kStationsNumberCrosscheck, ifInitialized);
  }
  if (fInitController.GetFlag(EInitKey::kStationsNumberCrosscheck)) {
    fNstationsGeometry[L1Constants::size::kMaxNdetectors] =
      std::accumulate(fNstationsGeometry.begin(), fNstationsGeometry.end() - 1, 0);
    fNstationsActive[L1Constants::size::kMaxNdetectors] = fNstationsGeometry[L1Constants::size::kMaxNdetectors];
  }
}

//-----------------------------------------------------------------------------------------------------------------------
//
void L1InitManager::SetTargetPosition(double x, double y, double z)
{
  if (fInitController.GetFlag(EInitKey::kTargetPos)) {
    LOG(warn) << "L1InitManager::SetTargetPosition: attempt to reinitialize the target position. Ignore";
    return;
  }

  /// Fill fvec target fields
  fParameters.fTargetPos[0] = x;
  fParameters.fTargetPos[1] = y;
  fParameters.fTargetPos[2] = z;
  /// Set additional field for z component in double precision
  fTargetZ = z;
  fInitController.SetFlag(EInitKey::kTargetPos);
}

//-----------------------------------------------------------------------------------------------------------------------
//
void L1InitManager::SetTrackingLevel(int trackingLevel)
{
  if (fInitController.GetFlag(EInitKey::kTrackingLevel)) {
    LOG(warn) << "L1InitManager::SetTrackingLevel: attempt of reinitialization the tracking level. Ignore";
    return;
  }
  fTrackingLevel = trackingLevel;
  fInitController.SetFlag(EInitKey::kTrackingLevel);
}

//-----------------------------------------------------------------------------------------------------------------------
//
void L1InitManager::TransferParametersContainer(L1Parameters& destination)
{
  this->FormParametersContainer();
  destination = std::move(fParameters);
  LOG(info) << "Parameters object transfered to L1Algo core";
}

//
// INIT CHECKERS
//

//-----------------------------------------------------------------------------------------------------------------------
//
void L1InitManager::CheckCAIterationsInit()
{
  //
  // 1) Check number of iterations
  //
  bool ifInitPassed = true;
  if (!fInitController.GetFlag(EInitKey::kCAIterations)) {
    int nIterationsActual   = fParameters.fCAIterations.size();
    int nIterationsExpected = fCAIterationsNumberCrosscheck;
    if (nIterationsActual != nIterationsExpected) {
      LOG(warn) << "L1InitManager::CheckCAIterations: incorrect number of iterations registered: " << nIterationsActual
                << " of " << nIterationsExpected << " expected";
      ifInitPassed = false;
    }
  }
  fInitController.SetFlag(EInitKey::kCAIterations, ifInitPassed);
}

//-----------------------------------------------------------------------------------------------------------------------
// TODO: REWRITE! and add const qualifier (S.Zharko)
void L1InitManager::CheckStationsInfoInit()
{
  bool ifInitPassed = true;
  if (!fInitController.GetFlag(EInitKey::kStationsInfo)) {
    //
    // 1) Check numbers of stations passed
    //
    // loop over active detectors
    for (auto itemDetector : fActiveDetectorIDs) {
      auto selectDetector = [&itemDetector](const L1BaseStationInfo& station) {
        return station.GetDetectorID() == itemDetector;
      };
      int nStationsExpected = GetNstationsActive(itemDetector);
      int nStations         = std::count_if(fStationsInfo.begin(), fStationsInfo.end(), selectDetector);
      if (nStations != nStationsExpected) {
        LOG(error) << "L1InitManager::CheckStationsInfoInit: Incorrect number of L1BaseStationInfo objects passed"
                   << " to the L1Manager for L1DetectorID = " << static_cast<int>(itemDetector) << ": " << nStations
                   << " of " << nStationsExpected << " expected";
        ifInitPassed = false;
      }
    }  // loop over active detectors: end
    L1MASSERT(0, ifInitPassed, "Station info initialization failed");

    //
    // 2) Check for maximum allowed number of stations
    //
    int nStationsTotal = fNstationsGeometry[fNstationsGeometry.size() - 1];
    if (nStationsTotal > L1Constants::size::kMaxNstations) {
      LOG(fatal) << "Actual total number of registered stations in geometry (" << nStationsTotal
                 << ") is larger then possible (" << L1Constants::size::kMaxNstations
                 << "). Please, select another set of active tracking detectors or recompile the code with enlarged"
                 << " L1Constants::size::kMaxNstations value";
    }
  }

  fInitController.SetFlag(EInitKey::kStationsInfo, ifInitPassed);
}