BmnRun9DM.cxx

Go to the documentation of this file.

/*
    BM@N alignment routine
    BM@N experiment at NICA complex, JINR, 2025
 
    Department: Math & Soft Group of HEP lab
    Author:     Igor Polev, polev@jinr.int
 
    BmnRun9DM class implementation (out-of-line vtable anchor and geometry)
 
    The enumeration itself is inline in the header - only the destructor and the
    geometry loading are kept here, so that the vtable and RTTI are emitted in
    this translation unit alone instead of once per includer, and the detector
    configuration headers are not pulled into every includer of the model.
*/
 
#include "BmnRun9DM.h"
 
#include "BmnGemStripStationSet.h"
#include "BmnSiliconStationSet.h"
#include "BmnVSPStationSet.h"
#include "TGeoSystemOfUnits.h"
 
#include <algorithm>
#include <fstream>
#include <memory>
#include <sstream>
#include <string>
 
/*
 *  The three XML configurations of Run9 (VSP_Run9, SiliconRun9, GemRun9) are
 *  already parsed by the detector classes, and the same files are used by the
 *  hit makers - see config/bmnconf.json. Only the module envelope (the zone
 *  envelope for GEM) is taken from them: pitch, strip angles and dead zones are
 *  not needed here.
 *
 *  The target is described apart from them, by a FairRoot ASCII geometry file
 *  (geometry/target_CsI.geo). Reading it through FairGeoSet would pull in a
 *  media table and a TGeoManager the alignment builds nowhere else, while the
 *  format is fixed and only a few of its numbers are wanted, so the few lines
 *  below read them directly.
 */
Bool_t BmnRun9DM::LoadGeometry(const char* targetGeoFile)
{
    constexpr Int_t PERIOD{9}, SETUP{0};
 
    // The file names the volume, its mother, the shape and the medium, one per
    // line, then gives a TUBE as the centre of its lower plane, its rmin and
    // rmax, and the centre of its upper plane, and finally the translation that
    // places the volume in its mother. The z span goes to the target layer,
    // unbounded in xy (see BmnDetectorModel::Layer); the tube's xy extent at
    // the middle of that span is the plane of the target element. Like every
    // FairGeo ASCII file it gives lengths in millimetres, while the alignment
    // works in centimetres. Without a file both stay zero.
    fTarget = {0.0, 0.0};
    Plane targetPlane{0.0, 0.0, 0.0, 0.0, 0.0};
    if (nullptr != targetGeoFile) {
        std::ifstream geo{targetGeoFile};
        if (!geo) {
            LOG(error) << " <!> BmnRun9DM : cannot open target geometry " << targetGeoFile;
            return kFALSE;
        }
        std::vector<std::string> fields;
        for (std::string line; std::getline(geo, line);) {
            const size_t first{line.find_first_not_of(" \t\r")};
            if (std::string::npos == first || '/' == line[first])
                continue;   // blank line or comment
            fields.push_back(line.substr(first, line.find_last_not_of(" \t\r") - first + 1));
        }
        if (fields.size() < 8 || "TUBE" != fields[2]) {
            LOG(error) << " <!> BmnRun9DM : " << targetGeoFile
                       << " is not a TUBE target in FairRoot ASCII geometry format";
            return kFALSE;
        }
        Double_t x{0.0}, y{0.0}, xUpper{0.0}, yUpper{0.0}, rMin{0.0}, rMax{0.0}, xPlacement{0.0}, yPlacement{0.0},
            zPlacement{0.0};
        std::istringstream lower{fields[4]}, radii{fields[5]}, upper{fields[6]}, placement{fields[7]};
        if (!(lower >> x >> y >> fTarget.zStart) || !(radii >> rMin >> rMax)
            || !(upper >> xUpper >> yUpper >> fTarget.zEnd) || !(placement >> xPlacement >> yPlacement >> zPlacement))
        {
            LOG(error) << " <!> BmnRun9DM : cannot read the target tube from " << targetGeoFile;
            return kFALSE;
        }
        fTarget.zStart += zPlacement;
        fTarget.zEnd += zPlacement;
        x += xPlacement;
        y += yPlacement;
        // conversion into cm
        fTarget.zStart *= TGeoUnit::millimeter;
        fTarget.zEnd *= TGeoUnit::millimeter;
        x *= TGeoUnit::millimeter;
        y *= TGeoUnit::millimeter;
        rMax *= TGeoUnit::millimeter;
        targetPlane = {x - rMax, x + rMax, y - rMax, y + rMax, 0.5 * (fTarget.zStart + fTarget.zEnd)};
    }
 
    const std::unique_ptr<BmnVSPStationSet> vsp{BmnVSPStationSet ::Create(PERIOD, SETUP)};
    const std::unique_ptr<BmnSiliconStationSet> sil{BmnSiliconStationSet ::Create(PERIOD, SETUP)};
    const std::unique_ptr<BmnGemStripStationSet> gem{BmnGemStripStationSet::Create(PERIOD, SETUP)};
 
    for (Int_t station = 0; station < ALL_STATION_CNT; station++)
        LOG_IF(warning, 0 == StationModules(station)) << "  ! BmnRun9DM : st " << station + STA_NUM_START
                                                      << " is left out of the model, its hits, if any, are ignored";
 
    fPlanes.clear();
    fPlanes.reserve(Alignables());
 
    // The three station sets share no base class, but name their geometry
    // getters alike. Indices are checked here because the sets throw or return
    // a null pointer on an out of range request. A negative zone takes the
    // whole module.
    auto append = [this](auto* set, Int_t station, Int_t module, Int_t zone) {
        if (station < 0 || station >= set->GetNStations())
            return kFALSE;
        auto* sta{set->GetStation(station)};
        if (nullptr == sta || module < 0 || module >= sta->GetNModules())
            return kFALSE;
        auto* mod{sta->GetModule(module)};
        if (nullptr == mod)
            return kFALSE;
        Double_t xMin{mod->GetXMinModule()}, xMax{mod->GetXMaxModule()}, yMin{mod->GetYMinModule()},
            yMax{mod->GetYMaxModule()};
        if (zone >= 0) {
            // zone n is the strip layer pair 2n, 2n+1: the rule by which
            // BmnToCbmHitConverter gives a hit its sector
            if (2 * zone + 1 >= mod->GetNStripLayers())
                return kFALSE;
            auto& lower{mod->GetStripLayer(2 * zone)};
            auto& upper{mod->GetStripLayer(2 * zone + 1)};
            xMin = std::min(lower.GetXMinLayer(), upper.GetXMinLayer());
            xMax = std::max(lower.GetXMaxLayer(), upper.GetXMaxLayer());
            yMin = std::min(lower.GetYMinLayer(), upper.GetYMinLayer());
            yMax = std::max(lower.GetYMaxLayer(), upper.GetYMaxLayer());
        }
        if (xMax <= xMin || yMax <= yMin)
            return kFALSE;
        // the configurations store x mirrored with respect to the hit
        // coordinates, therefore the bounds are negated and swapped
        fPlanes.push_back({-xMax, -xMin, yMin, yMax, mod->GetZStartModulePosition() + 0.5 * mod->GetModuleThickness()});
        return kTRUE;
    };
 
    for (Int_t idx = 0; idx < Alignables(); idx++) {
        const Int_t ID{EncodedID(idx)}, system{ID >> 24}, station{(ID >> 16) & 0xFF}, module{(ID >> 4) & 0xFFF};
 
        if (TargetID() == ID) {
            fPlanes.push_back(targetPlane);
            continue;
        }
 
        // A GEM element is one zone of a configuration module: the
        // BmnToCbmHitConverter sector is 2 * module + 1 + zone.
        const Bool_t known{
            VSP_SYSTEM_ID == system   ? append(vsp.get(), station - VSP_STATION_IDX[kFirst], module - MOD_NUM_START, -1)
            : SIL_SYSTEM_ID == system ? append(sil.get(), station - SIL_STATION_IDX[kFirst], module - MOD_NUM_START, -1)
            : GEM_SYSTEM_ID == system ? append(gem.get(), station - GEM_STATION_IDX[kFirst],
                                               (module - MOD_NUM_START) / 2, (module - MOD_NUM_START) % 2)
                                      : kFALSE};
 
        if (!known) {
            LOG(error) << " <!> BmnRun9DM : no geometry for element " << idx << " (system " << system << ", station "
                       << station << ", module " << module << ")";
            return kFALSE;
        }
    }
    return kTRUE;
}