AMVF algorithm

[rubihason]

Hello,
I am trying to get information about the tracks associated with each reconstructed vertex.
To do so, I write a CSV file from the AdaptiveMultiVertexFinderAlgorithm.cpp file as follows in the code.
I think something wrong here. When I compare the tracks number from the ambu to the tracks number from AMVF I get a bigger number of tracks inside the AMVF.
It doesn't make sense that the vertex algorithm has a bigger number of tracks than those in ambu.
for example I add 3 CSV files of the vertex information from AMVF, tracks associates AMVF and tracks from ambu.

vertex0.csv
tracks0.csv
event000000000-tracks_ambi.csv

// This file is part of the Acts project.
//
// Copyright (C) 2020-2021 CERN for the benefit of the Acts project
//
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at http://mozilla.org/MPL/2.0/.

#include <iostream>
#include <fstream>



#include "ActsExamples/Vertexin

g/AdaptiveMultiVertexFinderAlgorithm.hpp"

#include "Acts/Definitions/Algebra.hpp"
#include "Acts/Definitions/Units.hpp"
#include "Acts/Geometry/GeometryContext.hpp"
#include "Acts/MagneticField/MagneticFieldContext.hpp"
#include "Acts/Propagator/EigenStepper.hpp"
#include "Acts/Propagator/Propagator.hpp"
#include "Acts/Surfaces/PerigeeSurface.hpp"
#include "Acts/Utilities/Helpers.hpp"
#include "Acts/Utilities/Logger.hpp"
#include "Acts/Vertexing/AdaptiveMultiVertexFinder.hpp"
#include "Acts/Vertexing/AdaptiveMultiVertexFitter.hpp"
#include "Acts/Vertexing/HelicalTrackLinearizer.hpp"
#include "Acts/Vertexing/ImpactPointEstimator.hpp"
#include "Acts/Vertexing/LinearizedTrack.hpp"
#include "Acts/Vertexing/TrackDensityVertexFinder.hpp"
#include "Acts/Vertexing/Vertex.hpp"
#include "Acts/Vertexing/VertexingOptions.hpp"
#include "ActsExamples/EventData/ProtoVertex.hpp"
#include "ActsExamples/EventData/Track.hpp"
#include "ActsExamples/EventData/Trajectories.hpp"
#include "ActsExamples/Framework/ProcessCode.hpp"
#include "ActsExamples/Framework/WhiteBoard.hpp"

#include <chrono>
#include <memory>

#include "VertexingHelpers.hpp"

#include <filesystem>
#include <string>
#include <iostream>
#include <string.h>
#include <vector>
#include <algorithm>
#include <fstream>
using namespace std;
namespace fs = std::filesystem;

ActsExamples::AdaptiveMultiVertexFinderAlgorithm::
    AdaptiveMultiVertexFinderAlgorithm(const Config& config,
                                       Acts::Logging::Level level)
    : ActsExamples::IAlgorithm("AdaptiveMultiVertexFinder", level),
      m_cfg(config) {
  if (m_cfg.inputTrackParameters.empty() == m_cfg.inputTrajectories.empty()) {
    throw std::invalid_argument(
        "You have to either provide track parameters or trajectories");
  }
  if (m_cfg.outputProtoVertices.empty()) {
    throw std::invalid_argument("Missing output proto vertices collection");
  }
  if (m_cfg.outputVertices.empty()) {
    throw std::invalid_argument("Missing output vertices collection");
  }
  if (m_cfg.outputTime.empty()) {
    throw std::invalid_argument("Missing output reconstruction time");
  }

  m_inputTrackParameters.maybeInitialize(m_cfg.inputTrackParameters);
  m_inputTrajectories.maybeInitialize(m_cfg.inputTrajectories);

  m_outputProtoVertices.initialize(m_cfg.outputProtoVertices);
  m_outputVertices.initialize(m_cfg.outputVertices);
  m_outputTime.initialize(m_cfg.outputTime);
}


ActsExamples::ProcessCode
ActsExamples::AdaptiveMultiVertexFinderAlgorithm::execute(
    const ActsExamples::AlgorithmContext& ctx) const {
  // Set up EigenStepper
  Acts::EigenStepper<> stepper(m_cfg.bField);

  // Set up the propagator
  auto propagator = std::make_shared<Propagator>(stepper);

  // Set up ImpactPointEstimator
  IPEstimator::Config ipEstimatorCfg(m_cfg.bField, propagator);
  IPEstimator ipEstimator(ipEstimatorCfg);

  // Set up the helical track linearizer
  Linearizer::Config ltConfig(m_cfg.bField, propagator);
  Linearizer linearizer(ltConfig, logger().cloneWithSuffix("HelLin"));

  // Set up deterministic annealing with user-defined temperatures
  Acts::AnnealingUtility::Config annealingConfig;
  annealingConfig.setOfTemperatures = {1.};
  Acts::AnnealingUtility annealingUtility(annealingConfig);

  // Set up the vertex fitter with user-defined annealing
  Fitter::Config fitterCfg(ipEstimator);
  fitterCfg.annealingTool = annealingUtility;
  fitterCfg.minWeight = 0.001;
  fitterCfg.doSmoothing = true;
  Fitter fitter(fitterCfg, logger().cloneWithSuffix("AMVFitter"));

  // Set up the vertex seed finder
  SeedFinder seedFinder;

  Finder::Config finderConfig(std::move(fitter), seedFinder, ipEstimator,
                              std::move(linearizer), m_cfg.bField);
  // We do not want to use a beamspot constraint here
  finderConfig.useBeamSpotConstraint = false;
  finderConfig.tracksMaxZinterval = 1. * Acts::UnitConstants::mm;

  // Instantiate the finder
  Finder finder(finderConfig, logger().cloneWithSuffix("AMVFinder"));

  // retrieve input tracks and convert into the expected format

  auto [inputTrackParameters, inputTrackPointers] =
      makeParameterContainers(ctx, m_inputTrackParameters, m_inputTrajectories);

  //////////////////////////////////////////////
  /* Full tutorial example code for reference */
  //////////////////////////////////////////////

  // The vertex finder state
  Finder::State state;

  // Default vertexing options, this is where e.g. a constraint could be set
  using VertexingOptions = Acts::VertexingOptions<Acts::BoundTrackParameters>;
  VertexingOptions finderOpts(ctx.geoContext, ctx.magFieldContext);

  VertexCollection vertices;

  auto t1 = std::chrono::high_resolution_clock::now();

  if (inputTrackParameters.empty()) {
    ACTS_DEBUG("Empty track parameter collection found, skipping vertexing");
  } else {
    ACTS_DEBUG("Have " << inputTrackParameters.size()
                       << " input track parameters, running vertexing");
    // find vertices and measure elapsed time
    auto result = finder.find(inputTrackPointers, finderOpts, state);

    if (result.ok()) {
      vertices = std::move(result.value());
    } else {
      ACTS_ERROR("Error in vertex finder: " << result.error().message());
    }
  }

  auto t2 = std::chrono::high_resolution_clock::now();

  // show some debug output
  ACTS_INFO("Found " << vertices.size() << " vertices in event");
  for (const auto& vtx : vertices) {
    ACTS_INFO("Found vertex at " << vtx.fullPosition().transpose() << " with "
                                  << vtx.tracks().size() << " tracks.");
  }

int i = 0;
int j = 0;


  std::ofstream myfile,myfile1;
  std::string vertex = "Vertecis/vertex";
  std::string tracks = "Tracks/tracks";
  std::string file_type = ".csv";

  std::string full_vertex_name = vertex + std::to_string(ctx.eventNumber) + file_type;
  std::string full_tracks_name = tracks + std::to_string(ctx.eventNumber) + file_type;
  
  ACTS_INFO("full_vertex_name:")
  ACTS_INFO(full_vertex_name)
  ACTS_INFO("full_tracks_name:")
  ACTS_INFO(full_tracks_name)

  myfile1.open(full_vertex_name);
  myfile.open(full_tracks_name);
  myfile1 << "Number,position,chi squared,number of degrees of freedom\n";
  myfile << "Number,vertex number,compatibility,weight,number of degrees of freedom,chi2,absoluteMomentum,transverseMomentum,charge,Phi,Theta,loc0,loc1\n";
  for (const auto& vtx : vertices) 
  {
      

    i = i + 1; 
    

     myfile1  << i <<"," <<vtx.position()[2] <<","<<  to_string(vtx.m_chiSquared) << ","<< to_string(vtx.m_numberDoF) <<",\n";
 
    
    for (const auto& trx : vtx.tracks()) 
    {
        j = j + 1;
        myfile << "Track " << j << "," << i << "," << trx.vertexCompatibility <<","<<  trx.trackWeight << ","<< trx.ndf << ","<< trx.chi2Track << ","<< trx.originalParams->absoluteMomentum()<< ","<< trx.originalParams->transverseMomentum() <<","<< trx.originalParams->charge() <<","<< trx.originalParams->phi() << ","<< trx.originalParams->Theta() << ","<< trx.originalParams->loc0() << ","<< trx.originalParams->loc1() << ",\n";
    }
  }

  myfile.close();
  myfile1.close();
  // store proto vertices extracted from the found vertices
  m_outputProtoVertices(ctx, makeProtoVertices(inputTrackParameters, vertices));

  // store found vertices
  m_outputVertices(ctx, std::move(vertices));

  // time in milliseconds
  int timeMS =
      std::chrono::duration_cast<std::chrono::milliseconds>(t2 - t1).count();
  // store reconstruction time
  m_outputTime(ctx,
               std::move(timeMS));  // NOLINT(performance-move-const-arg)

  return ActsExamples::ProcessCode::SUCCESS;
}

[andiwand]

are you sure the AMVF is receiving the tracks from the ambi solver and not from the CKF directly?

[rubihason]

I think so, this is the code:

#!/usr/bin/env python3
from pathlib import Path
import os, argparse, pathlib, contextlib, acts, acts.examples
from acts.examples import Sequencer, GenericDetector, RootParticleReader
import os, argparse, pathlib, contextlib, acts, acts.examples
from acts.examples.simulation import (
    addPythia8,
    addFatras,
    addDigitization,
)
from acts.examples.reconstruction import (
    addSeeding,
    TruthSeedRanges,
    ParticleSmearingSigmas,
    addCKFTracks,
    SeedFinderConfigArg,
    CKFPerformanceConfig,
    addSeeding,
    TruthSeedRanges,
    TrackSelectorRanges,
    ParticleSmearingSigmas,
    SeedFinderConfigArg,
    CKFPerformanceConfig,
    addCKFTracks,
    addAmbiguityResolutionML,
    AmbiguityResolutionMLConfig,
    addAmbiguityResolution,
    AmbiguityResolutionConfig,
    addVertexFitting,
    VertexFinder,
)

from pathlib import Path
from typing import Optional, Union
from collections import namedtuple

from acts.examples import Sequencer, GenericDetector, RootParticleReader

import acts

from acts import UnitConstants as u
srcdir = Path(__file__).resolve().parent.parent.parent.parent
# PreliminariesonMode23:onMode=off
rnd = acts.examples.RandomNumbers()
outputDir = Path("/root/Thesis/acts/Examples/Scripts/Python/Data")
detector, trackingGeometry, decorators = GenericDetector.create()
field = acts.ConstantBField(acts.Vector3(0, 0, 2 * u.T))
digiConfigFile=srcdir / "Examples/Algorithms/Digitization/share/default-smearing-config-generic.json"
geometrySelection=srcdir/ "Examples/Algorithms/TrackFinding/share/geoSelection-genericDetector.json"
    # Sequencer
s =  acts.examples.Sequencer(
        events=1, numThreads=-1, logLevel=acts.logging.INFO
)
addPythia8(
        s,
        rnd=rnd,
        nhard=0,
        npileup=1,
        cmsEnergy=13.6 * acts.UnitConstants.TeV,
        hardProcess=["WeakSingleBoson:ffbar2ffbar(s:gmZ)=on" ,"WeakZ0:gmZmode=2" ,"23:onMode=off" ,"23:onIfAny=13","Beams:allowVertexSpread =on","Beams:sigmaVertexZ=1"],
        pileupProcess=["SoftQCD:all = on","Beams:allowVertexSpread =on","Beams:sigmaVertexZ=1"],
        outputDirRoot = "/root/Thesis/acts/Examples/Scripts/Python/Data",
        outputDirCsv=Path("/root/Thesis/acts/Examples/Scripts/Python/Pythia8") ,
        #outputDirRoot=outputDir if outputRoot else None,
    )

addFatras(
        s,
        trackingGeometry,
        field,
        rnd=rnd,
        # outputDirCsv=outputDir,
    )

addDigitization(
        s,
        trackingGeometry,
        field,
        digiConfigFile=digiConfigFile,
        rnd=rnd,
        # outputDirCsv=outputDir,
    )

addSeeding(
        s,
        trackingGeometry,
        field,
        TruthSeedRanges(pt=(500.0 * u.MeV, None), nHits=(9, None)),
        ParticleSmearingSigmas(pRel=0.01),  # only used by SeedingAlgorithm.TruthSmeared
        SeedFinderConfigArg(
            r=(None, 200 * u.mm),  # rMin=default, 33mm
            deltaR=(1 * u.mm, 60 * u.mm),
            collisionRegion=(-250 * u.mm, 250 * u.mm),
            z=(-2000 * u.mm, 2000 * u.mm),
            maxSeedsPerSpM=1,
            sigmaScattering=5,
            radLengthPerSeed=0.1,
            minPt=500 * u.MeV,
            impactMax=3 * u.mm,
        ),

        geoSelectionConfigFile=geometrySelection,
        outputDirRoot=outputDir,
        rnd=rnd,  # only used by SeedingAlgorithm.TruthSmeared
    )

addCKFTracks(
    s,
        trackingGeometry,
        field,
        CKFPerformanceConfig(ptMin=400.0 * u.MeV, nMeasurementsMin=6),
        TrackSelectorRanges(pt=(0.5* u.GeV, None),absEta=(None, 2.5)),
        outputDirRoot=outputDir,
        outputDirCsv=Path("/root/Thesis/acts/Examples/Scripts/Python/CKFTracks"),

    )

    



addAmbiguityResolution(
            s,
            AmbiguityResolutionConfig(
                maximumSharedHits=3, maximumIterations=10000, nMeasurementsMin=7
            ),
            CKFPerformanceConfig(0.0,nMeasurementsMin=7
            ),
            outputDirRoot=outputDir,
            outputDirCsv= Path("/root/Thesis/acts/Examples/Scripts/Python/AmbiguityResolution"),
        )

addVertexFitting(
        s,
        field,
        vertexFinder=VertexFinder.AMVF,
        outputDirRoot=outputDir,
    )

s.run()

[rubihason]

I just lookup on the file tracksummary_ckf.root and I see a different number of tracks from the csv file.
maybe there is some problem with the csv output?

[andiwand]

which CSV writer are you using? also be aware that there might be some filtering applied to one of the output but not to the other

[rubihason]

I using the default writer for addCKFTracks.
you say that in the csv file may have some filters inside the writer that are different from the root writer?