Skip to content

Commit

Permalink
pull changes in
Browse files Browse the repository at this point in the history
  • Loading branch information
@benjaminhuth
benjaminhuth committed Jul 9, 2024
1 parent e4f3d90 commit 22386f4
Show file tree
Hide file tree
Showing 2 changed files with 156 additions and 23 deletions.
30 changes: 25 additions & 5 deletions Examples/Io/Root/include/ActsExamples/Io/Root/RootAthenaDumpReader.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,13 +7,16 @@
// file, You can obtain one at http://mozilla.org/MPL/2.0/.

#pragma once

#include "Acts/Utilities/Logger.hpp"
#include "ActsExamples/EventData/Measurement.hpp"
#include "ActsExamples/EventData/SimSpacePoint.hpp"
#include "ActsExamples/Framework/DataHandle.hpp"
#include "ActsExamples/Framework/IReader.hpp"
#include "ActsExamples/Framework/ProcessCode.hpp"
#include <ActsExamples/EventData/Cluster.hpp>
#include <ActsExamples/EventData/SimParticle.hpp>
#include <ActsExamples/EventData/Track.hpp>

#include <map>
#include <memory>
Expand Down Expand Up @@ -43,15 +46,26 @@ class RootAthenaDumpReader : public IReader {
// Name of inputfile
std::string inputfile;
// name of the output measurements
std::string outputMeasurements = "ath_meas";
std::string outputMeasurements = "athena_measurements";
// name of the output pixel space points
std::string outputPixelSpacePoints = "outputPixelSpacepoints";
std::string outputPixelSpacePoints = "athena_pixel_spacepoints";
// name of the output strip space points
std::string outputStripSpacePoints = "outputStripSpacepoints";
std::string outputStripSpacePoints = "athena_strip_spacepoints";
// name of the output space points
std::string outputSpacePoints = "output_spacepoints";
std::string outputSpacePoints = "athena_spacepoints";
// name of the output clusters
std::string outputClusters = "output_clusters";
std::string outputClusters = "athena_clusters";
// name of the output particles
std::string outputParticles = "athena_particles";
// name of the simhit map
std::string outputMeasurementParticlesMap = "athena_meas_parts_map";
// name of the track parameters (fitted by athena?)
std::string outputTrackParameters = "athena_track_parameters";

/// Only extract particles that passed the tracking requirements, for
/// details see:
/// https://gitlab.cern.ch/atlas/athena/-/blob/main/InnerDetector/InDetGNNTracking/src/DumpObjects.cxx?ref_type=heads#L1363
bool onlyPassedParticles = false;
};

RootAthenaDumpReader(const RootAthenaDumpReader &) = delete;
Expand Down Expand Up @@ -91,6 +105,12 @@ class RootAthenaDumpReader : public IReader {
WriteDataHandle<SimSpacePointContainer> m_outputSpacePoints{
this, "output_spacepoints"};
WriteDataHandle<ClusterContainer> m_outputClusters{this, "output_clusters"};
WriteDataHandle<SimParticleContainer> m_outputParticles{this,
"output_particles"};
WriteDataHandle<MeasurementContainer> m_outputMeasurements{
this, "output_measurements"};
WriteDataHandle<IndexMultimap<ActsFatras::Barcode>> m_outputMeasParticleMap{
this, "output_meas_part_map"};

std::unique_ptr<const Acts::Logger> m_logger;
std::mutex m_read_mutex;
Expand Down
149 changes: 131 additions & 18 deletions Examples/Io/Root/src/RootAthenaDumpReader.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,10 +14,60 @@
#include "ActsExamples/EventData/Cluster.hpp"
#include "ActsExamples/EventData/GeometryContainers.hpp"
#include "ActsExamples/EventData/IndexSourceLink.hpp"
#include <ActsExamples/Digitization/MeasurementCreation.hpp>

#include <cmath>

#include <TChain.h>
#include <boost/container/static_vector.hpp>

class BarcodeConstructor {
/// Particles with barcodes larger then this value are considered to be
/// secondary particles
constexpr static int s_maxBarcodeForPrimary = 200000;

std::uint16_t m_primaryCount = 0;
std::uint16_t m_secondaryCount = 0;
std::unordered_map<std::uint64_t, ActsFatras::Barcode> m_barcodeMap;

static std::uint64_t concatInts(int a, int b) {
auto va = static_cast<std::uint32_t>(a);
auto vb = static_cast<std::uint32_t>(b);
std::uint64_t value = (static_cast<std::uint64_t>(va) << 32) | vb;
return value;
}

public:
ActsFatras::Barcode getBarcode(int barcode, int evtnumber) {
auto v = concatInts(barcode, evtnumber);
auto found = m_barcodeMap.find(v);
if (found != m_barcodeMap.end()) {
return found->second;
}

auto primary = (barcode < s_maxBarcodeForPrimary);

ActsFatras::Barcode fBarcode;

// vertex primary shouldn't be zero for a valid particle
fBarcode.setVertexPrimary(1);
if (primary) {
fBarcode.setVertexSecondary(0);
fBarcode.setParticle(m_primaryCount);
assert(m_primaryCount < std::numeric_limits<std::uint16_t>::max());
m_primaryCount++;
} else {
fBarcode.setVertexSecondary(1);
fBarcode.setParticle(m_secondaryCount);
assert(m_primaryCount < std::numeric_limits<std::uint16_t>::max());
m_secondaryCount++;
}

m_barcodeMap[v] = fBarcode;
return fBarcode;
}
};

enum SpacePointType { ePixel = 1, eStrip = 2 };

ActsExamples::RootAthenaDumpReader::RootAthenaDumpReader(
Expand All @@ -39,6 +89,9 @@ ActsExamples::RootAthenaDumpReader::RootAthenaDumpReader(
m_outputStripSpacePoints.initialize(m_cfg.outputStripSpacePoints);
m_outputSpacePoints.initialize(m_cfg.outputSpacePoints);
m_outputClusters.initialize(m_cfg.outputClusters);
m_outputParticles.initialize(m_cfg.outputParticles);
m_outputMeasParticleMap.initialize(m_cfg.outputMeasurementParticlesMap);
m_outputMeasurements.initialize(m_cfg.outputMeasurements);

// Set the branches

Expand Down Expand Up @@ -214,23 +267,45 @@ ActsExamples::ProcessCode ActsExamples::RootAthenaDumpReader::read(

m_inputchain->GetEntry(entry);

// Loop on clusters (measurements)
ACTS_DEBUG("Found " << nSP << " space points");
ACTS_DEBUG("Found " << nCL << " clusters / measurements");
// Concat the two 32bit integers from athena to a Fatras barcode

SimParticleContainer particles;
BarcodeConstructor barcodeConstructor;

for (auto ip = 0; ip < nPartEVT; ++ip) {
if (m_cfg.onlyPassedParticles && !static_cast<bool>(Part_passed[ip])) {
continue;
}

auto barcode =
barcodeConstructor.getBarcode(Part_barcode[ip], Part_event_number[ip]);
SimParticle particle(barcode,
static_cast<Acts::PdgParticle>(Part_pdg_id[ip]));

ClusterContainer clusters;
clusters.resize(nCL);
auto p = Acts::Vector3{Part_px[ip], Part_py[ip], Part_pz[ip]};
particle.setAbsoluteMomentum(p.norm());
particle.setDirection(p.normalized());

auto x = Acts::Vector4{Part_vx[ip], Part_vy[ip], Part_vz[ip], 0.0};
particle.setPosition4(x);

particles.insert(particle);
}

ClusterContainer clusters(nCL);
MeasurementContainer measurements;
measurements.reserve(nCL);

IndexMultimap<ActsFatras::Barcode> measPartMap;

for (int im = 0; im < nCL; im++) {
int bec = CLbarrel_endcap[im];
int lydisk = CLlayer_disk[im];
int etamod = CLeta_module[im];
int phimod = CLphi_module[im];
int side = CLside[im];
// ULong64_t moduleID = CLmoduleID [im];
if (!(CLhardware->at(im) == "PIXEL" || CLhardware->at(im) == "STRIP")) {
ACTS_ERROR("hardware is neither 'PIXEL' or 'STRIP'");
return ActsExamples::ProcessCode::ABORT;
}
ACTS_VERBOSE("Cluster " << im << ": " << CLhardware->at(im));

ACTS_VERBOSE(bec << " " << lydisk << " " << etamod << " " << phimod << " "
<< side << " ");
auto type = (CLhardware->at(im) == "PIXEL") ? ePixel : eStrip;

// Make cluster
// TODO refactor ActsExamples::Cluster class so it is not so tedious
Expand Down Expand Up @@ -265,11 +340,46 @@ ActsExamples::ProcessCode ActsExamples::RootAthenaDumpReader::read(
activation);
}

ACTS_VERBOSE("Cluster " << im << ": " << cluster.channels.size()
<< "cells, dimensions: " << cluster.sizeLoc0 << ", "
<< cluster.sizeLoc1);
cluster.globalPosition = {CLx[im], CLy[im], CLz[im]};

ACTS_VERBOSE("CL shape: " << cluster.channels.size()
<< "cells, dimensions: " << cluster.sizeLoc0
<< ", " << cluster.sizeLoc1);

clusters[im] = cluster;

// Measurement creation
ACTS_VERBOSE("CL loc dims:" << CLloc_direction1[im] << ", "
<< CLloc_direction2[im] << ", "
<< CLloc_direction3[im]);
const auto& locCov = CLlocal_cov->at(im);

DigitizedParameters digiPars;
if (type == ePixel) {
digiPars.indices = {Acts::eBoundLoc0, Acts::eBoundLoc1};
digiPars.values = {CLloc_direction1[im], CLloc_direction2[im]};
assert(locCov.size() == 4);
digiPars.variances = {locCov[0], locCov[3]};
} else {
// TODO is this correct ???
digiPars.values = {CLloc_direction2[im]};
digiPars.indices = {Acts::eBoundLoc1};

// why is this 4
// assert(locCov.size() == 1);
digiPars.variances = {locCov[0]};
}

IndexSourceLink sl(Acts::GeometryIdentifier{CLmoduleID[im]}, im);

measurements.push_back(createMeasurement(digiPars, sl));

// Create measurement particles map and particles container
for (const auto& [subevt, bc] : Acts::zip(CLparticleLink_eventIndex->at(im),
CLparticleLink_barcode->at(im))) {
auto barcode = barcodeConstructor.getBarcode(bc, subevt);
measPartMap.insert({im, barcode});
}
}

// Prepare pixel space points
Expand Down Expand Up @@ -323,14 +433,17 @@ ActsExamples::ProcessCode ActsExamples::RootAthenaDumpReader::read(
spacePoints.push_back(sp);
}

ACTS_DEBUG("Created " << particles.size() << " particles");
ACTS_DEBUG("Created " << spacePoints.size() << " overall space points");
ACTS_DEBUG("Created " << pixelSpacePoints.size() << " "
<< " pixel space points");

ACTS_DEBUG("Created " << spacePoints.size() << " overall space points");

m_outputPixelSpacePoints(ctx, std::move(pixelSpacePoints));
m_outputSpacePoints(ctx, std::move(spacePoints));
m_outputClusters(ctx, std::move(clusters));
m_outputParticles(ctx, std::move(particles));
m_outputMeasParticleMap(ctx, std::move(measPartMap));
m_outputMeasurements(ctx, std::move(measurements));

return ProcessCode::SUCCESS;
}

0 comments on commit 22386f4

Please sign in to comment.