Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Improve MC truth information in TICL, Sim-Reco Associators, TICLDumper #42518

Merged
merged 12 commits into from
Oct 27, 2023
9 changes: 9 additions & 0 deletions DataFormats/HGCalReco/interface/TICLCandidate.h
Original file line number Diff line number Diff line change
Expand Up @@ -53,7 +53,16 @@ class TICLCandidate : public reco::LeafCandidate {
return idProbabilities_[(int)type];
}

inline const std::array<float, 8>& idProbabilities() const { return idProbabilities_; }

void zeroProbabilities() {
for (auto& p : idProbabilities_) {
p = 0.f;
}
}

void setIdProbabilities(const std::array<float, 8>& idProbs) { idProbabilities_ = idProbs; }
inline void setIdProbability(ParticleType type, float value) { idProbabilities_[int(type)] = value; }

private:
float time_;
Expand Down
76 changes: 36 additions & 40 deletions DataFormats/HGCalReco/interface/Trackster.h
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@
namespace ticl {
class Trackster {
public:
typedef math::XYZVector Vector;
typedef math::XYZVectorF Vector;

enum IterationIndex { TRKEM = 0, EM, TRKHAD, HAD, MIP, SIM, SIM_CP };

Expand All @@ -37,21 +37,20 @@ namespace ticl {
enum class PCAOrdering { ascending = 0, descending };

Trackster()
: iterationIndex_(0),
seedIndex_(-1),
time_(0.f),
timeError_(-1.f),
: barycenter_({0.f, 0.f, 0.f}),
regressed_energy_(0.f),
raw_energy_(0.f),
time_(0.f),
timeError_(-1.f),
raw_em_energy_(0.f),
id_probabilities_{},
raw_pt_(0.f),
raw_em_pt_(0.f),
barycenter_({0., 0., 0.}),
eigenvalues_{{0.f, 0.f, 0.f}},
sigmas_{{0.f, 0.f, 0.f}},
sigmasPCA_{{0.f, 0.f, 0.f}} {
zeroProbabilities();
}
seedIndex_(-1),
eigenvalues_{},
sigmas_{},
sigmasPCA_{},
iterationIndex_(0) {}

inline void setIteration(const Trackster::IterationIndex i) { iterationIndex_ = i; }
std::vector<unsigned int> &vertices() { return vertices_; }
Expand All @@ -73,6 +72,9 @@ namespace ticl {
inline void setRawPt(float value) { raw_pt_ = value; }
inline void setRawEmPt(float value) { raw_em_pt_ = value; }
inline void setBarycenter(Vector value) { barycenter_ = value; }
inline void setTrackIdx(int index) { track_idx_ = index; }
int trackIdx() const { return track_idx_; }

inline void fillPCAVariables(Eigen::Vector3d &eigenvalues,
Eigen::Matrix3d &eigenvectors,
Eigen::Vector3d &sigmas,
Expand Down Expand Up @@ -147,20 +149,23 @@ namespace ticl {
}

private:
// TICL iteration producing the trackster
uint8_t iterationIndex_;
Vector barycenter_;
float regressed_energy_;
float raw_energy_;
// -99, -1 if not available. ns units otherwise
float time_;
float timeError_;
float raw_em_energy_;

// trackster ID probabilities
std::array<float, 8> id_probabilities_;

// The vertices of the DAG are the indices of the
// 2d objects in the global collection
std::vector<unsigned int> vertices_;
std::vector<float> vertex_multiplicity_;

// The edges connect two vertices together in a directed doublet
// ATTENTION: order matters!
// A doublet generator should create edges in which:
// the first element is on the inner layer and
// the outer element is on the outer layer.
std::vector<std::array<unsigned int, 2> > edges_;
float raw_pt_;
float raw_em_pt_;

// Product ID of the seeding collection used to create the Trackster.
// For GlobalSeeding the ProductID is set to 0. For track-based seeding
Expand All @@ -173,31 +178,22 @@ namespace ticl {
// created the trackster. For track-based seeding the pointer to the track
// can be cooked using the previous ProductID and this index.
int seedIndex_;
int track_idx_ = -1;

// We also need the pointer to the original seeding region ??
// something like:
// int seedingRegionIdx;

// -99, -1 if not available. ns units otherwise
float time_;
float timeError_;

// regressed energy
float regressed_energy_;
float raw_energy_;
float raw_em_energy_;
float raw_pt_;
float raw_em_pt_;

// PCA Variables
Vector barycenter_;
std::array<float, 3> eigenvalues_;
std::array<Vector, 3> eigenvectors_;
std::array<float, 3> eigenvalues_;
std::array<float, 3> sigmas_;
std::array<float, 3> sigmasPCA_;

// trackster ID probabilities
std::array<float, 8> id_probabilities_;
// The edges connect two vertices together in a directed doublet
// ATTENTION: order matters!
// A doublet generator should create edges in which:
// the first element is on the inner layer and
// the outer element is on the outer layer.
std::vector<std::array<unsigned int, 2> > edges_;

// TICL iteration producing the trackster
uint8_t iterationIndex_;
};

typedef std::vector<Trackster> TracksterCollection;
Expand Down
3 changes: 2 additions & 1 deletion DataFormats/HGCalReco/src/classes_def.xml
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@

<lcgdict>
<class name="ticl::Trackster" ClassVersion="9">
<class name="ticl::Trackster" ClassVersion="10">
<version ClassVersion="10" checksum="556627704"/>
<version ClassVersion="9" checksum="1001808235"/>
<version ClassVersion="8" checksum="749412802"/>
<version ClassVersion="7" checksum="1072014319"/>
Expand Down
15 changes: 14 additions & 1 deletion RecoHGCal/Configuration/python/RecoHGCal_EventContent_cff.py
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,12 @@
'keep *_ticlTrackstersHFNoseMIP_*_*',
'keep *_ticlTrackstersHFNoseHAD_*_*',
'keep *_ticlTrackstersHFNoseMerge_*_*',] +
['keep *_pfTICL_*_*']
['keep *_pfTICL_*_*'] +
['keep CaloParticles_mix_*_*', 'keep SimClusters_mix_*_*'] +
['keep *_layerClusterSimClusterAssociationProducer_*_*','keep *_layerClusterCaloParticleAssociationProducer_*_*', 'keep *_layerClusterSimTracksterAssociationProducer_*_*'] +
['keep *_tracksterSimTracksterAssociationLinking_*_*' ,'keep *_tracksterSimTracksterAssociationPR_*_*'] +
['keep *_tracksterSimTracksterAssociationLinkingPU_*_*' ,'keep *_tracksterSimTracksterAssociationPRPU_*_*'] +
['keep *_tracksterSimTracksterAssociationLinkingbyCLUE3D_*_*', 'keep *_tracksterSimTracksterAssociationPRbyCLUE3D_*_*']
)
)
TICL_RECO.outputCommands.extend(TICL_AOD.outputCommands)
Expand All @@ -28,6 +33,7 @@
TICL_FEVT = cms.PSet(
outputCommands = cms.untracked.vstring(
'keep *_ticlSimTracksters_*_*',
'keep *_ticlSimTICLCandidates_*_*',
'keep *_ticlSimTrackstersFromCP_*_*',
)
)
Expand All @@ -48,6 +54,13 @@ def cleanOutputAndSet(outputModule, ticl_outputCommads):
'keep *_layerClusterSimClusterAssociationProducer_*_*',
'keep *_layerClusterCaloParticleAssociationProducer_*_*',
'keep *_randomEngineStateProducer_*_*',
'keep *_layerClusterSimTracksterAssociationProducer_*_*',
'keep *_tracksterSimTracksterAssociationLinking_*_*',
'keep *_tracksterSimTracksterAssociationPR_*_*',
'keep *_tracksterSimTracksterAssociationLinkingPU_*_*',
'keep *_tracksterSimTracksterAssociationPRPU_*_*',
'keep *_tracksterSimTracksterAssociationLinkingbyCLUE3D_*_*',
'keep *_tracksterSimTracksterAssociationPRbyCLUE3D_*_*',
])

if hasattr(process, 'FEVTDEBUGEventContent'):
Expand Down
54 changes: 54 additions & 0 deletions RecoHGCal/TICL/interface/commons.h
Original file line number Diff line number Diff line change
Expand Up @@ -4,12 +4,28 @@
#include "DataFormats/Common/interface/Ref.h"
#include "DataFormats/CaloRecHit/interface/CaloClusterFwd.h"
#include "DataFormats/HGCalReco/interface/Trackster.h"
#include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h"

namespace ticl {

//constants
constexpr double mpion = 0.13957;
constexpr float mpion2 = mpion * mpion;
typedef math::XYZVectorF Vector;
enum LayerType {

CE_E_120 = 0,
CE_E_200 = 1,
CE_E_300 = 2,
CE_H_120_F = 3,
CE_H_200_F = 4,
CE_H_300_F = 5,
CE_H_200_C = 6,
CE_H_300_C = 7,
CE_H_SCINT_C = 8,
EnumSize = 9

};

inline Trackster::ParticleType tracksterParticleTypeFromPdgId(int pdgId, int charge) {
if (pdgId == 111) {
Expand Down Expand Up @@ -40,6 +56,44 @@ namespace ticl {
// verbosity levels for ticl algorithms
enum VerbosityLevel { None = 0, Basic, Advanced, Expert, Guru };

inline int returnClusterType(DetId& lc_seed, const hgcal::RecHitTools& rhtools_) {
auto layer_number = rhtools_.getLayerWithOffset(lc_seed);
auto thickness = rhtools_.getSiThickIndex(lc_seed);
auto isEELayer = (layer_number <= rhtools_.lastLayerEE(false));
auto isScintillator = rhtools_.isScintillator(lc_seed);
auto isFine = (layer_number <= rhtools_.lastLayerEE(false) + 7);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

7?
No hard coded numbers: these will cause troubles in the long term, especially if related to the geometry.
On more general grounds, maybe all this should be moved to RecHittools?
In the end, the type of cluster is derived from its seed which is a RecHit.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also, why returning -1? If, for whatever reason, we reach that statement, all downstream code has to check if the type returned is valid.
A nice assert is safer, IMHO. At least you die where you were supposed to.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@waredjeb Do you want to comment here?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@srimanob Yes, we discussed it offline, I moved this piece to the RecHitTools and added an assert to avoid returning -1

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry, it's still confusing why we have a 7 there? Can this be clarified?


if (isScintillator) {
return CE_H_SCINT_C;
}
if (isEELayer) {
if (thickness == 0) {
return CE_E_120;
} else if (thickness == 1) {
return CE_E_200;
} else if (thickness == 2) {
return CE_E_300;
}
} else {
if (isFine) {
if (thickness == 0) {
return CE_H_120_F;
} else if (thickness == 1) {
return CE_H_200_F;
} else if (thickness == 2) {
return CE_H_300_F;
}
} else {
if (thickness == 1) {
return CE_H_200_C;
} else if (thickness == 2) {
return CE_H_300_C;
}
}
}
return -1;
};

} // namespace ticl

#endif
2 changes: 2 additions & 0 deletions RecoHGCal/TICL/plugins/BuildFile.xml
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,8 @@
<use name="TrackingTools/GeomPropagators"/>
<use name="TrackingTools/Records"/>
<use name="TrackingTools/TrajectoryState"/>
<use name="PhysicsTools/UtilAlgos"/>
<use name="FWCore/ServiceRegistry"/>
<use name="fastjet"/>
<library file="*.cc" name="RecoHGCalTICLPlugins">
<flags EDM_PLUGIN="1"/>
Expand Down
9 changes: 5 additions & 4 deletions RecoHGCal/TICL/plugins/LinkingAlgoByDirectionGeometric.cc
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"

#include "RecoParticleFlow/PFProducer/interface/PFMuonAlgo.h"
#include "DataFormats/HGCalReco/interface/Trackster.h"

using namespace ticl;

Expand All @@ -36,10 +37,10 @@ void LinkingAlgoByDirectionGeometric::initialize(const HGCalDDDConstants *hgcons
propagator_ = propH;
}

math::XYZVector LinkingAlgoByDirectionGeometric::propagateTrackster(const Trackster &t,
const unsigned idx,
float zVal,
std::array<TICLLayerTile, 2> &tracksterTiles) {
Vector LinkingAlgoByDirectionGeometric::propagateTrackster(const Trackster &t,
const unsigned idx,
float zVal,
std::array<TICLLayerTile, 2> &tracksterTiles) {
// needs only the positive Z co-ordinate of the surface to propagate to
// the correct sign is calculated inside according to the barycenter of trackster
Vector const &baryc = t.barycenter();
Expand Down
11 changes: 6 additions & 5 deletions RecoHGCal/TICL/plugins/LinkingAlgoByDirectionGeometric.h
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@
#include "CommonTools/Utils/interface/StringCutObjectSelector.h"

#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
#include "DataFormats/HGCalReco/interface/Trackster.h"

namespace ticl {
class LinkingAlgoByDirectionGeometric final : public LinkingAlgoBase {
Expand All @@ -44,14 +45,14 @@ namespace ticl {
static void fillPSetDescription(edm::ParameterSetDescription &desc);

private:
typedef math::XYZVector Vector;
using Vector = ticl::Trackster::Vector;

void buildLayers();

math::XYZVector propagateTrackster(const Trackster &t,
const unsigned idx,
float zVal,
std::array<TICLLayerTile, 2> &tracksterTiles);
Vector propagateTrackster(const Trackster &t,
const unsigned idx,
float zVal,
std::array<TICLLayerTile, 2> &tracksterTiles);

void findTrackstersInWindow(const std::vector<std::pair<Vector, unsigned>> &seedingCollection,
const std::array<TICLLayerTile, 2> &tracksterTiles,
Expand Down
Loading