Commit a42c5156 authored by Plyushchenko's avatar Plyushchenko
Browse files

undo merge master because it's broken

parent d8368824
......@@ -23,3 +23,4 @@ setup.sh
*.root
.vscode
.idea
output/*
......@@ -24,7 +24,7 @@ set(CMAKE_CXX_FLAGS_DEBUG "-O0 -g -DDEBUG")
set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -lpthread")
# Configuration of CUDA compute architecture
set(CUDA_ARCH "MIN" CACHE STRING "Cuda architecture")
set(CUDA_ARCH "MAX" CACHE STRING "Cuda architecture")
if (CUDA_ARCH STREQUAL "MIN" OR CUDA_ARCH STREQUAL "MAX" OR CUDA_ARCH STREQUAL "COMP")
set(OUTPUTFILE ${PROJECT_BINARY_DIR}/cuda_arch) # No suffix required
......@@ -71,10 +71,7 @@ else(USE_LZ4)
set(LZ4_FOUND OFF)
endif(USE_LZ4)
# Set Kalman single precision.
option(USE_KALMAN_SINGLE_PRECISION
"Use single precision in ParKalman"
ON)
option(USE_KALMAN_SINGLE_PRECISION ON)
# find_package(CUDA REQUIRED)
#set(CUDA_HOST_COMPILER "g++")
......@@ -100,9 +97,9 @@ endif()
# Cuda: Deal with build type
if(${CMAKE_BUILD_TYPE} STREQUAL RelWithDebInfo)
set(CMAKE_CUDA_FLAGS_RELWITHDEBINFO "-O3 -g -DNDEBUG --generate-line-info ")
set(CMAKE_CUDA_FLAGS_RELWITHDEBINFO "-O3 -g -DNDEBUG --generate-line-info")
elseif(${CMAKE_BUILD_TYPE} STREQUAL Release)
set(CMAKE_CUDA_FLAGS_RELEASE "-O3 -DNDEBUG --generate-line-info ")
set(CMAKE_CUDA_FLAGS_RELEASE "-O3 -DNDEBUG --generate-line-info")
elseif(${CMAKE_BUILD_TYPE} STREQUAL Debug)
set(CMAKE_CUDA_FLAGS_DEBUG "-O0 -G -g -DDEBUG ")
endif(${CMAKE_BUILD_TYPE} STREQUAL RelWithDebInfo)
......@@ -146,6 +143,25 @@ include_directories(cuda/velo/mask_clustering/include)
include_directories(cuda/velo/search_by_triplet/include)
include_directories(cuda/velo/simplified_kalman_filter/include)
include_directories(cuda/SciFi/common/include)
include_directories(cuda/SciFi/looking_forward/common/include)
include_directories(cuda/SciFi/looking_forward/calculate_first_layer_window/include)
include_directories(cuda/SciFi/looking_forward/calculate_second_layer_window/include)
include_directories(cuda/SciFi/looking_forward/form_seeds_from_candidates/include)
include_directories(cuda/SciFi/looking_forward/calculate_candidate_extrapolation_window/include)
include_directories(cuda/SciFi/looking_forward/promote_candidates/include)
include_directories(cuda/SciFi/looking_forward/calculate_track_extrapolation_window/include)
include_directories(cuda/SciFi/looking_forward/extend_tracks/include)
include_directories(cuda/SciFi/looking_forward_sbt/search_initial_windows/include)
include_directories(cuda/SciFi/looking_forward_sbt/collect_candidates/include)
include_directories(cuda/SciFi/looking_forward_sbt/triplet_seeding/include)
include_directories(cuda/SciFi/looking_forward_sbt/extend_tracks_x/include)
include_directories(cuda/SciFi/looking_forward_sbt/composite_algorithms/include)
include_directories(cuda/SciFi/looking_forward_sbt/extend_tracks_uv/include)
include_directories(cuda/SciFi/looking_forward_sbt/quality_filter/include)
include_directories(cuda/SciFi/looking_forward_sbt/quality_filter_x/include)
include_directories(cuda/SciFi/looking_forward_sbt/search_uv_windows/include)
include_directories(cuda/SciFi/PrForward/include)
include_directories(cuda/SciFi/consolidate/include)
include_directories(cuda/muon/common/include)
include_directories(cuda/utils/prefix_sum/include)
include_directories(cuda/event_model/velo/include)
......@@ -155,9 +171,9 @@ include_directories(cuda/event_model/common/include)
include_directories(checker/tracking/include)
include_directories(checker/pv/include)
include_directories(stream/sequence/include)
include_directories(x86/SciFi/include)
include_directories(cuda/SciFi/PrForward/include)
include_directories(cuda/SciFi/consolidate/include)
include_directories(x86/SciFi/PrForward/include)
include_directories(x86/SciFi/LookingForward/include)
include_directories(x86/SciFi/MomentumForward/include)
include_directories(cuda/UT/UTDecoding/include)
include_directories(cuda/kalman/ParKalman/include)
include_directories(mdf/include)
......
......@@ -43,22 +43,55 @@ struct CheckerInvoker {
std::tuple<bool, MCEvents> read_mc_folder() const;
template<typename T>
void check(const uint start_event_offset, const std::vector<trackChecker::Tracks>& tracks) const
std::vector<std::vector<std::vector<uint32_t>>> check(
const uint start_event_offset,
const std::vector<Checker::Tracks>& tracks,
std::vector<std::vector<float>>& p_events) const
{
std::vector<std::vector<std::vector<uint32_t>>> scifi_ids_events;
if (is_mc_folder_populated) {
T trackChecker {};
T track_checker {};
#ifdef WITH_ROOT
trackChecker.histos.initHistos(trackChecker.histo_categories());
track_checker.histos.initHistos(track_checker.histo_categories());
#endif
for (int evnum = 0; evnum < selected_mc_events.size(); ++evnum) {
const auto& mc_event = selected_mc_events[evnum];
const auto& event_tracks = tracks[evnum];
const auto& mcps = mc_event.mc_particles<T>();
MCAssociator mcassoc {mcps};
std::vector<uint32_t> matched_mcp_keys =
track_checker(event_tracks, mc_event, get_num_hits_subdetector<typename T::subdetector_t>);
std::vector<std::vector<uint32_t>> scifi_ids_tracks;
std::vector<float> p_tracks;
for (const auto key : matched_mcp_keys) {
std::vector<uint32_t> scifi_ids;
float p = 1e9;
if (!(key == 0xFFFFFF)) { // track was matched to an MCP
// Find this MCP
for (const auto mcp : mc_event.m_mcps) {
if (mcp.key == key) {
// Save momentum and charge of this MCP
p = mcp.p * mcp.charge;
// debug_cout << "Adding particle with PID = " << mcp.pid << " and charge " << charge << std::endl;
// Find SciFi IDs of this MCP
if (mcp.isLong) { // found matched long MCP
for (const auto id : mcp.hits) {
const uint32_t detector_id = (id >> 20) & 0xFFF;
if (detector_id == 0xa00) { // hit in the SciFi
scifi_ids.push_back(id);
}
}
}
}
}
}
scifi_ids_tracks.push_back(scifi_ids);
p_tracks.push_back(p);
}
trackChecker(event_tracks, mcassoc, mcps);
scifi_ids_events.push_back(scifi_ids_tracks);
p_events.push_back(p_tracks);
// Check all tracks for duplicate LHCb IDs
for (int i_track = 0; i_track < event_tracks.size(); ++i_track) {
......@@ -78,5 +111,6 @@ struct CheckerInvoker {
}
}
}
return scifi_ids_events;
}
};
......@@ -15,12 +15,21 @@
#include "LHCbID.h"
namespace trackChecker {
class Track {
namespace Checker {
namespace Subdetector {
struct Velo;
struct UT;
struct SciFi;
} // namespace Subdetector
struct TruthCounter {
uint n_velo {0};
uint n_ut {0};
uint n_scifi {0};
};
public:
SomeLHCbIDs allids;
struct Track {
LHCbIDs allids;
// Kalman information.
float z, x, y, tx, ty, qop;
float first_qop, best_qop;
......@@ -32,14 +41,16 @@ namespace trackChecker {
float velo_docaz;
float long_ip, long_ip_chi2, long_ipx, long_ipy;
std::size_t n_matched_total = 0;
float p;
float p, pt, eta;
float muon_catboost_output;
bool is_muon;
void addId(LHCbID id) { allids.push_back(id); }
SomeLHCbIDs ids() const { return allids; }
LHCbIDs ids() const { return allids; }
int nIDs() const { return allids.size(); }
};
using Tracks = std::vector<Track>;
} // namespace trackChecker
} // namespace Checker
......@@ -7,7 +7,7 @@
#include "Logger.h"
#include "MCAssociator.h"
#include "MCEvent.h"
#include "Tracks.h"
#include "CheckerTypes.h"
#ifdef WITH_ROOT
#include "TTree.h"
......@@ -16,5 +16,5 @@
void checkKalmanTracks(
const uint start_event_offset,
const std::vector<trackChecker::Tracks>& tracks,
const std::vector<Checker::Tracks>& tracks,
const MCEvents selected_mc_events);
......@@ -13,6 +13,7 @@
#include <array>
#include <cstdint>
#include <cmath>
#include <vector>
/// encapsulate an LHCbID
class LHCbID {
......@@ -22,6 +23,8 @@ private:
public:
constexpr LHCbID(uint32_t id) : m_id(id) {}
enum channelIDtype { Velo = 0x8, FT = 0xa, UT = 0xb };
LHCbID() = default;
LHCbID(const LHCbID& other) = default;
LHCbID(LHCbID&& other) = default;
......@@ -61,8 +64,11 @@ public:
/// ordering of LHCbIDs
constexpr bool operator>=(const LHCbID& other) const noexcept { return m_id >= other.m_id; }
// FIXME: ultimately, more methods are needed to e.g. get Velo sensor
// numbers for hits etc.
// get subdetector
inline uint32_t detectorType() const { return (unsigned int) ((m_id & 0xF0000000L) >> 28); }
bool isVelo() const { return (Velo == detectorType()); };
bool isUT() const { return (UT == detectorType()); };
bool isSciFi() const { return (FT == detectorType()); };
};
typedef std::vector<LHCbID> SomeLHCbIDs;
typedef std::vector<LHCbID> LHCbIDs;
......@@ -18,10 +18,10 @@
#include "LHCbID.h"
#include "MCParticle.h"
#include "Logger.h"
#include "CheckerTypes.h"
/// simple MC associator
class MCAssociator {
private:
struct MCAssociator {
using LHCbIDWithIndex = std::pair<LHCbID, uint>;
using AssocMap = std::vector<LHCbIDWithIndex>;
......@@ -29,19 +29,34 @@ private:
struct MCParticleWithWeight {
std::size_t m_idx;
float m_w;
MCParticleWithWeight(std::size_t idx, float w) : m_idx(idx), m_w(w) {}
uint m_counter_sum;
MCParticleWithWeight(std::size_t idx, float w, uint counter_sum) : m_idx(idx), m_w(w), m_counter_sum(counter_sum) {}
MCParticleWithWeight(const MCParticleWithWeight&) = default;
MCParticleWithWeight(MCParticleWithWeight&&) = default;
MCParticleWithWeight& operator=(const MCParticleWithWeight&) = default;
MCParticleWithWeight& operator=(MCParticleWithWeight&&) = default;
};
// internal structure with index into tracks matched to an MCP
struct TrackWithWeight {
int m_idx;
float m_w;
int m_counter_subdetector;
TrackWithWeight(int idx, float w, int counter_subdetector) :
m_idx(idx), m_w(w), m_counter_subdetector(counter_subdetector)
{}
TrackWithWeight(const TrackWithWeight&) = default;
TrackWithWeight(TrackWithWeight&&) = default;
TrackWithWeight& operator=(const TrackWithWeight&) = default;
TrackWithWeight& operator=(TrackWithWeight&&) = default;
};
const MCParticles& m_mcps; // keep a reference to MCParticles
AssocMap m_map; // association LHCbID -> MCParticle index
// little helper which does the hard work
AssocMap::const_iterator find(LHCbID id) const noexcept;
AssocMap::const_iterator find_id(const LHCbID&) const noexcept;
AssocMap::const_iterator find_id(const LHCbID&, const AssocMap::const_iterator& begin) const noexcept;
std::vector<AssocMap::const_iterator> find_ids(const LHCbID&) const noexcept;
public:
MCAssociator(const MCParticles& mcps);
MCAssociator(const MCAssociator&) = default;
MCAssociator(MCAssociator&&) = default;
......@@ -52,6 +67,7 @@ public:
class MCAssocResult {
private:
using AssocVector = std::vector<MCParticleWithWeight>;
using AssocTable = std::map<MCParticle, std::vector<TrackWithWeight>>;
AssocVector m_assoc;
const MCParticles& m_mcps;
......@@ -136,30 +152,30 @@ public:
const auto& b = reinterpret_cast<const AssocVector::const_iterator&>(other);
return a >= b;
}
std::pair<MCParticles::const_reference, float> operator*() const noexcept
std::tuple<MCParticles::const_reference, float, int> operator*() const noexcept
{
AssocVector::const_reference ref = reinterpret_cast<const AssocVector::const_iterator&>(*this).operator*();
return {m_mcps[ref.m_idx], ref.m_w};
return {m_mcps[ref.m_idx], ref.m_w, ref.m_counter_sum};
}
};
iterator begin() const noexcept { return iterator(m_assoc.begin(), m_mcps); }
iterator end() const noexcept { return iterator(m_assoc.end(), m_mcps); }
std::pair<MCParticles::const_reference, float> front() const noexcept { return *begin(); }
std::pair<MCParticles::const_reference, float> back() const noexcept { return *--end(); }
std::tuple<MCParticles::const_reference, float, int> front() const noexcept { return *begin(); }
std::tuple<MCParticles::const_reference, float, int> back() const noexcept { return *--end(); }
};
private:
// private:
using AssocPreResult = std::map<std::size_t, std::size_t>;
/// little helper for the final step of multi-MCP association
MCAssocResult buildResult(const AssocPreResult& assocmap, std::size_t total) const noexcept;
public:
// public:
/// associate a single LHCbID
MCAssocResult operator()(LHCbID id) const noexcept
{
auto it = find(id);
auto it = find_id(id);
if (m_map.end() == it) return MCAssocResult({}, m_mcps);
return MCAssocResult({{it->second, 1.f}}, m_mcps);
return MCAssocResult({{it->second, 1.f, 1}}, m_mcps);
}
/// associate a range of LHCbIDs
template<typename IT>
......@@ -171,9 +187,9 @@ public:
// and how many hits of the reconstructed track are matched
// to the MCP
for (; last != first; ++first) {
auto it = find(*first);
if (m_map.end() == it) continue;
// std::cout << "Matched LHCbID to MCP: " << std::hex << *first << std::endl;
const auto it = find_id(*first);
if (it == m_map.end()) continue;
// std::cout << "Matched LHCbID to MCP: " << *first << " to " << it->second << std::endl;
++n_matched_total;
++assoc[it->second];
}
......
......@@ -18,38 +18,29 @@
#include <string>
#include <vector>
#include <algorithm>
#include "MCParticle.h"
#include "Common.h"
#include "Logger.h"
#include "MCParticle.h"
#include "TrackChecker.h"
#include "LHCbID.h"
#include "CheckerTypes.h"
struct MCEvent {
MCParticles velo_mcps;
MCParticles ut_mcps;
MCParticles scifi_mcps;
MCParticles m_mcps;
uint32_t size;
// Constructor
MCEvent() {};
MCEvent(const std::vector<char>& _event, const bool checkFile = true);
/**
* @brief Print a specific set of MC particles.
*/
void print(const MCParticles& mcps) const;
// Checks if a LHCb ID is in a particular subdetector
bool is_subdetector_impl(const LHCbIDs& vector, const LHCbID& id) const;
// Subdetector-specialized check
template<typename T>
const MCParticles& mc_particles() const;
bool is_subdetector(const LHCbID& id) const;
/**
* @brief Print the set of MC particles of type T.
*/
template<typename T>
void print() const
{
print(mc_particles<T>());
}
// Checks an MCP does not contain invalid values
void check_mcp(const MCParticle& mcp);
};
using MCEvents = std::vector<MCEvent>;
......@@ -12,6 +12,8 @@
#pragma once
#include "CheckerTypes.h"
// Monte Carlo information
struct MCParticle {
uint32_t key;
......@@ -28,6 +30,11 @@ struct MCParticle {
bool fromBeautyDecay;
bool fromCharmDecay;
bool fromStrangeDecay;
uint32_t motherKey;
float charge;
uint32_t velo_num_hits;
uint32_t ut_num_hits;
uint32_t scifi_num_hits;
uint32_t numHits;
uint32_t nPV; // # of reconstructible primary vertices in event
std::vector<uint32_t> hits;
......@@ -36,4 +43,10 @@ struct MCParticle {
bool inEta2_5() const { return (eta < 5. && eta > 2.); };
};
template<typename T>
uint32_t get_num_hits(const MCParticle& mc_particle);
template<typename T>
uint32_t get_num_hits_subdetector(const MCParticle& mc_particle);
using MCParticles = std::vector<MCParticle>;
#pragma once
#include <vector>
#include "Tracks.h"
#include "CheckerTypes.h"
#include "Logger.h"
#include "InputTools.h"
#include "UTDefinitions.cuh"
......@@ -25,16 +25,16 @@ float ipChi2Kalman(const ParKalmanFilter::FittedTrack& track, const PV::Vertex&
float kalmanDOCAz(const ParKalmanFilter::FittedTrack& track, const PV::Vertex& vertex);
// Velo tracks.
float ipVelo(const Velo::Consolidated::States& velo_kalman_states, const uint state_index, const PV::Vertex& vertex);
float ipxVelo(const Velo::Consolidated::States& velo_kalman_states, const uint state_index, const PV::Vertex& vertex);
float ipyVelo(const Velo::Consolidated::States& velo_kalman_states, const uint state_index, const PV::Vertex& vertex);
float ipVelo(const Velo::Consolidated::KalmanStates& velo_kalman_states, const uint state_index, const PV::Vertex& vertex);
float ipxVelo(const Velo::Consolidated::KalmanStates& velo_kalman_states, const uint state_index, const PV::Vertex& vertex);
float ipyVelo(const Velo::Consolidated::KalmanStates& velo_kalman_states, const uint state_index, const PV::Vertex& vertex);
float ipChi2Velo(
const Velo::Consolidated::States& velo_kalman_states,
const Velo::Consolidated::KalmanStates& velo_kalman_states,
const uint state_index,
const PV::Vertex& vertex);
float veloDOCAz(const Velo::Consolidated::States& velo_kalman_states, const uint state_index, const PV::Vertex& vertex);
float veloDOCAz(const Velo::Consolidated::KalmanStates& velo_kalman_states, const uint state_index, const PV::Vertex& vertex);
std::vector<trackChecker::Tracks> prepareKalmanTracks(
std::vector<Checker::Tracks> prepareKalmanTracks(
const uint* velo_track_atomics,
const uint* velo_track_hit_number,
const char* velo_track_hits,
......
......@@ -5,13 +5,16 @@
#include "Logger.h"
#include "UTDefinitions.cuh"
#include "SciFiDefinitions.cuh"
#include "SciFiEventModel.cuh"
#include "UTEventModel.cuh"
#include "MiniState.cuh"
#include "States.cuh"
float eta_from_rho(const float rho);
/**
* @brief Prepares tracks for Velo consolidated datatypes.
*/
std::vector<trackChecker::Tracks> prepareVeloTracks(
std::vector<Checker::Tracks> prepareVeloTracks(
const uint* track_atomics,
const uint* track_hit_number,
const char* track_hits,
......@@ -20,10 +23,11 @@ std::vector<trackChecker::Tracks> prepareVeloTracks(
/**
* @brief Prepares tracks for Velo, UT consolidated datatypes.
*/
std::vector<trackChecker::Tracks> prepareUTTracks(
std::vector<Checker::Tracks> prepareUTTracks(
const uint* velo_track_atomics,
const uint* velo_track_hit_number,
const char* velo_track_hits,
const char* kalman_velo_states,
const int* ut_track_atomics,
const uint* ut_track_hit_number,
const char* ut_track_hits,
......@@ -34,10 +38,11 @@ std::vector<trackChecker::Tracks> prepareUTTracks(
/**
* @brief Prepares tracks for Velo, UT, SciFi consolidated datatypes.
*/
std::vector<trackChecker::Tracks> prepareSciFiTracks(
std::vector<Checker::Tracks> prepareSciFiTracks(
const uint* velo_track_atomics,
const uint* velo_track_hit_number,
const char* velo_track_hits,
const char* kalman_velo_states,
const int* ut_track_atomics,
const uint* ut_track_hit_number,
const char* ut_track_hits,
......@@ -51,4 +56,21 @@ std::vector<trackChecker::Tracks> prepareSciFiTracks(
const MiniState* scifi_states,
const char* host_scifi_geometry,
const std::array<float, 9>& host_inv_clus_res,
const float* muon_catboost_output,
const bool* is_muon,
const uint number_of_events);
std::vector<Checker::Tracks> prepareSciFiTracks(
const uint* velo_track_atomics,
const uint* velo_track_hit_number,
const char* velo_track_hits,
const char* kalman_velo_states,
const int* ut_track_atomics,
const uint* ut_track_hit_number,
const char* ut_track_hits,
const uint* ut_track_velo_indices,
const float* ut_qop,
const std::vector<std::vector<SciFi::TrackHits>>& scifi_tracks,
const SciFi::Hits& scifi_hits,
const uint* host_scifi_hit_count,
const uint number_of_events);
......@@ -20,7 +20,8 @@
#include <vector>
#include "Logger.h"
#include "MCAssociator.h"
#include "Tracks.h"
#include "CheckerTypes.h"
#include "MCEvent.h"
#ifdef WITH_ROOT
#include "TDirectory.h"
......@@ -31,6 +32,8 @@
class TrackChecker {
protected:
bool m_print = false;
using AcceptFn = std::function<bool(MCParticles::const_reference&)>;
struct TrackEffReport {
std::string m_name;
......@@ -44,7 +47,11 @@ protected:
float m_effperevt = 0.f;
float m_hitpur = 0.f;
float m_hiteff = 0.f;
std::set<uint32_t> m_keysseen;
std::size_t m_naccept_per_event = 0;
std::size_t m_nfound_per_event = 0;
std::size_t m_nclones_per_event = 0;
float m_eff_per_event = 0.f;
float m_number_of_events = 0.f;
/// no default construction
TrackEffReport() = delete;
......@@ -67,10 +74,14 @@ protected:
/// register MC particles
void operator()(const MCParticles& mcps);
/// register track and its MC association
void
operator()(trackChecker::Tracks::const_reference& track, MCParticles::const_reference& mcp, const float weight);
/// notify of end of event
void evtEnds();
void operator()(
const std::vector<MCAssociator::TrackWithWeight> tracks,
MCParticles::const_reference& mcp,
const std::function<uint32_t(const MCParticle&)>& get_num_hits_subdetector);
void event_start();
void event_done();
/// free resources, and print result
~TrackEffReport();
};
......@@ -78,7 +89,6 @@ protected:
struct HistoCategory {
std::string m_name;
AcceptFn m_accept;
std::set<uint32_t> m_keysseen;
/// construction from name and accept criterion for eff. denom.
template<typename F>
......@@ -88,8 +98,6 @@ protected:
template<typename F>
HistoCategory(std::string&& name, F&& accept) : m_name(std::move(name)), m_accept(std::move(accept))
{}
/// notify of end of event
void evtEnds();
};
std::vector<TrackEffReport> m_categories;
......@@ -116,6 +124,12 @@ protected:
TH2D* h_qop_resolution;
TH2D* h_dqop_versus_qop;
TH1D* h_momentum_matched;
TH1D* h_muon_catboost_output_matched_muon;
TH1D* h_muon_catboost_output_matched_notMuon;
TH1D* h_muon_catboost_output_matched_muon_ismuon_true;
TH1D* h_muon_catboost_output_matched_notMuon_ismuon_true;
TH1D* h_is_muon_matched_muon;
TH1D* h_is_muon_matched_notMuon;