RecActsTracking.cpp

#include <iostream>
#include <fstream>
#include <cstdlib>
#include <sstream>

// dependence
#include "RecActsTracking.h"

#include "GearSvc/IGearSvc.h"

// csv parser
// #include "csv2/writer.hpp"

// MC
#include "CLHEP/Units/SystemOfUnits.h"


using namespace Acts::UnitLiterals;

DECLARE_COMPONENT(RecActsTracking)

RecActsTracking::RecActsTracking(const std::string& name, ISvcLocator* svcLoc)
    : GaudiAlgorithm(name, svcLoc)
{
}

StatusCode RecActsTracking::initialize()
{
    chronoStatSvc = service<IChronoStatSvc>("ChronoStatSvc");
    _nEvt = 0;

    if (!std::filesystem::exists(TGeo_path.value())) {
        error() << "CEPC TGeo file: " << TGeo_path.value() << " does not exist!" << endmsg;
        return StatusCode::FAILURE;
    }

    if (!std::filesystem::exists(TGeo_config_path.value())) {
        error() << "CEPC TGeo config file: " << TGeo_config_path.value() << " does not exist!" << endmsg;
        return StatusCode::FAILURE;
    }

    if (!std::filesystem::exists(MaterialMap_path.value())) {
        error() << "CEPC Material map file: " << MaterialMap_path.value() << " does not exist!" << endmsg;
        return StatusCode::FAILURE;
    }

    if(!m_particle.value().empty()){
        info() << "Assume Particle: " << m_particle.value() << endmsg;
        if(m_particle.value() == "muon"){
            particleHypothesis = Acts::ParticleHypothesis::muon();
        }
        else if(m_particle.value() == "pion"){
            particleHypothesis = Acts::ParticleHypothesis::pion();
        }
        else if(m_particle.value() == "electron"){
            particleHypothesis = Acts::ParticleHypothesis::electron();
        }
        else if(m_particle.value() == "kaon"){
            particleHypothesis = Acts::ParticleHypothesis::kaon();
        }
        else if(m_particle.value() == "proton"){
            particleHypothesis = Acts::ParticleHypothesis::proton();
        }
        else if(m_particle.value() == "photon"){
            particleHypothesis = Acts::ParticleHypothesis::photon();
        }
        else if(m_particle.value() == "geantino"){
            particleHypothesis = Acts::ParticleHypothesis::geantino();
        }
        else if(m_particle.value() == "chargedgeantino"){
            particleHypothesis = Acts::ParticleHypothesis::chargedGeantino();
        }
        else{
            info() << "Supported Assumed Particle: " << particleNames[0] << ", " << particleNames[1] << ", " << particleNames[2] << ", "
                                             << particleNames[3] << ", " << particleNames[4] << ", " << particleNames[5] << ", "
                                             << particleNames[6] << ", " << particleNames[7] << endmsg;
            error() << "Unsupported particle name " << m_particle.value() << endmsg;
            return StatusCode::FAILURE;
        }
    }

    TGeo_ROOTFilePath = TGeo_path.value();
    TGeoConfig_jFilePath = TGeo_config_path.value();
    MaterialMap_jFilePath = MaterialMap_path.value();

    chronoStatSvc->chronoStart("read geometry");
    m_geosvc = service<IGeomSvc>("GeomSvc");
    if (!m_geosvc) {
        error() << "Failed to find GeomSvc." << endmsg;
        return StatusCode::FAILURE;
    }

    if(m_geosvc){
        const dd4hep::Direction& field = m_geosvc->lcdd()->field().magneticField(dd4hep::Position(0,0,0));
        m_field = field.z()/dd4hep::tesla;
    }

    m_vtx_surfaces = m_geosvc->getSurfaceMap("VXD");
    debug() << "Surface map size: " << m_vtx_surfaces->size() << endmsg;

    m_sit_surfaces = m_geosvc->getSurfaceMap("SIT");
    debug() << "Surface map size: " << m_sit_surfaces->size() << endmsg;

    m_ftd_surfaces = m_geosvc->getSurfaceMap("FTD");
    debug() << "Surface map size: " << m_ftd_surfaces->size() << endmsg;

    vxd_decoder = m_geosvc->getDecoder("VXDCollection");
    if(!vxd_decoder){
        return StatusCode::FAILURE;
    }

    sit_decoder = m_geosvc->getDecoder("SITCollection");
    if(!sit_decoder){
        return StatusCode::FAILURE;
    }

    ftd_decoder = m_geosvc->getDecoder("FTDCollection");
    if(!ftd_decoder){
        return StatusCode::FAILURE;
    }

    info() << "ACTS Tracking Geometry initialized successfully!" << endmsg;
    // initialize tgeo detector
    auto logger = Acts::getDefaultLogger("TGeoDetector", Acts::Logging::INFO);
    trackingGeometry = buildTGeoDetector(geoContext, detElementStore, TGeo_ROOTFilePath, TGeoConfig_jFilePath, MaterialMap_jFilePath, *logger);

    info() << "Seeding tools initialized successfully!" << endmsg;

    // configure the acts tools
    seed_cfg.seedFinderOptions.bFieldInZ = m_field.value();
    seed_cfg.seedFinderConfig.deltaRMin = SeedDeltaRMin.value();
    seed_cfg.seedFinderConfig.deltaRMax = SeedDeltaRMax.value();
    seed_cfg.seedFinderConfig.rMax = SeedRMax.value();
    seed_cfg.seedFinderConfig.rMin = SeedRMin.value();
    seed_cfg.seedFinderConfig.impactMax = SeedImpactMax.value();
    seed_cfg.seedFinderConfig.useVariableMiddleSPRange = false;
    seed_cfg.seedFinderConfig.rMinMiddle = SeedRMinMiddle.value();
    seed_cfg.seedFinderConfig.rMaxMiddle = SeedRMaxMiddle.value();

    // initialize the acts tools
    seed_cfg.seedFilterConfig = seed_cfg.seedFilterConfig.toInternalUnits();
    seed_cfg.seedFinderConfig.seedFilter =
        std::make_unique<Acts::SeedFilter<SimSpacePoint>>(seed_cfg.seedFilterConfig);
    seed_cfg.seedFinderConfig =
        seed_cfg.seedFinderConfig.toInternalUnits().calculateDerivedQuantities();
    seed_cfg.seedFinderOptions =
        seed_cfg.seedFinderOptions.toInternalUnits().calculateDerivedQuantities(seed_cfg.seedFinderConfig);
    seed_cfg.gridConfig = seed_cfg.gridConfig.toInternalUnits();
    seed_cfg.gridOptions = seed_cfg.gridOptions.toInternalUnits();

    if (std::isnan(seed_cfg.seedFinderConfig.deltaRMaxTopSP)) {
        seed_cfg.seedFinderConfig.deltaRMaxTopSP = seed_cfg.seedFinderConfig.deltaRMax;}
    if (std::isnan(seed_cfg.seedFinderConfig.deltaRMinTopSP)) {
        seed_cfg.seedFinderConfig.deltaRMinTopSP = seed_cfg.seedFinderConfig.deltaRMin;}
    if (std::isnan(seed_cfg.seedFinderConfig.deltaRMaxBottomSP)) {
        seed_cfg.seedFinderConfig.deltaRMaxBottomSP = seed_cfg.seedFinderConfig.deltaRMax;}
    if (std::isnan(seed_cfg.seedFinderConfig.deltaRMinBottomSP)) {
        seed_cfg.seedFinderConfig.deltaRMinBottomSP = seed_cfg.seedFinderConfig.deltaRMin;}
        
    m_bottomBinFinder = std::make_unique<const Acts::GridBinFinder<2ul>>(
        seed_cfg.numPhiNeighbors, seed_cfg.zBinNeighborsBottom);
    m_topBinFinder = std::make_unique<const Acts::GridBinFinder<2ul>>(
        seed_cfg.numPhiNeighbors, seed_cfg.zBinNeighborsTop);

    seed_cfg.seedFinderConfig.seedFilter =
        std::make_unique<Acts::SeedFilter<SimSpacePoint>>(seed_cfg.seedFilterConfig);
    m_seedFinder =
        Acts::SeedFinder<SimSpacePoint, Acts::CylindricalSpacePointGrid<SimSpacePoint>>(seed_cfg.seedFinderConfig);

    // initialize the ckf
    findTracks = makeTrackFinderFunction(trackingGeometry, magneticField);

    info() << "CKF Track Finder initialized successfully!" << endmsg;
    chronoStatSvc->chronoStop("read geometry");

    return GaudiAlgorithm::initialize();
}

StatusCode RecActsTracking::execute()
{
    auto trkCol = _outColHdl.createAndPut();

    SpacePointPtrs.clear();
    sourceLinks.clear();
    measurements.clear();
    initialParameters.clear();
    Selected_Seeds.clear();

    chronoStatSvc->chronoStart("read input hits");

    int successVTX = InitialiseVTX();
    if (successVTX == 0)
    {
        _nEvt++;
        return StatusCode::SUCCESS;
    }

    int successSIT = InitialiseSIT();
    if (successSIT == 0)
    {