10#include "GaudiKernel/SystemOfUnits.h"
49 float result(std::numeric_limits<float>::quiet_NaN());
52 if (truthMatchProbabilityAcc.isAvailable(trackParticle)) {
53 result = truthMatchProbabilityAcc(trackParticle);
60 safelyGetEta(
const T& pTrk,
const float safePtThreshold = 0.1) {
61 return (pTrk->pt() > safePtThreshold) ? (pTrk->eta()) : std::nan(
"");
67 inRange(
const T& value,
const T& minVal,
const T& maxVal) {
68 return not ((value < minVal)or(value > maxVal));
73 inRange(
const T& value,
const T& absoluteMax) {
74 return not (std::abs(value) > absoluteMax);
82 binIndex(
const T& val,
const std::vector<T>& partitions) {
89 return nf ?
i :
i - 1;
94 const float x(vtx->
x()),
y(vtx->
y()),
z(vtx->
z());
95 const float vrtR2 = (
x *
x +
y *
y);
97 return inRange(
z, 500.f) and not (vrtR2 > 1);
101 2.7, 3.5, std::numeric_limits<float>::infinity()
108 const IInterface* parent) :
126 {not m_vertexContainerName.key().empty() &&
127 not m_hardScatterSelectionTool.empty()}));
137 (
nullptr,
static_cast<std::string
>(
m_dirName) +
static_cast<std::string
>(
m_folder),
148 std::vector<std::string> trk_decorations;
164 std::vector<std::string> required_float_track_decorations {
"d0",
"hitResiduals_residualLocX",
"d0err"};
165 std::vector<std::string> required_int_track_decorations {};
166 std::vector<std::string> required_float_truth_decorations {
"d0"};
167 std::vector<std::string> required_int_truth_decorations {};
168 std::vector<std::string> required_int_jet_decorations {
"HadronConeExclTruthLabelID"};
173 m_datfile <<
"EvtNumber,PdgID,Px,Py,Pz,E,Pt,Eta,Phi,Mass,numPU,numvtx,MatchProb"<<std::endl;
175 std::string empty_prefix;
189 return StatusCode::SUCCESS;
316 if (not trackHandle.
isValid()) {
318 return StatusCode::FAILURE;
329 ATH_MSG_WARNING(
"Shouldn't happen - EventInfo is buggy, setting mu to 0");
337 return StatusCode::SUCCESS;
341 std::vector<const xAOD::TruthParticle*> truthParticlesVec =
getTruthParticles(ctx);
350 unsigned int truthMu = 0;
353 truthMu =
static_cast<int>( truthPileupEventContainer->size() );
355 if(pie.
isValid()) actualMu = pie->actualInteractionsPerCrossing();
358 float puEvents = truthMu>0 ? truthMu : actualMu;
361 unsigned int nVertices = 0;
362 float beamSpotWeight = 1;
365 ATH_MSG_DEBUG(
"Getting number of pu interactings per event");
371 nVertices = not vertices->empty() ? vertices->size() : 0;
372 beamSpotWeight = pie->beamSpotWeight();
373 ATH_MSG_DEBUG(
"beamSpotWeight is equal to " << beamSpotWeight);
377 if(readDecorHandle.
isAvailable()) prwWeight = readDecorHandle(*pie);
378 ATH_MSG_DEBUG(
"Applying pileup weight equal to " << prwWeight);
379 beamSpotWeight *= prwWeight;
382 if (vertices.
isValid() and not vertices->empty()) {
383 ATH_MSG_DEBUG(
"Number of vertices retrieved for this event " << vertices->size());
386 for (
const auto& vertex : *vertices){
388 primaryvertex = vertex;
399 ATH_MSG_DEBUG(
"Failed to find a hard scatter vertex in this event.");
405 std::pair<std::vector<const xAOD::TruthVertex*>, std::vector<const xAOD::TruthVertex*>> truthVertices =
getTruthVertices(ctx);
406 std::vector<const xAOD::TruthVertex*> truthHSVertices = truthVertices.first;
407 std::vector<const xAOD::TruthVertex*> truthPUVertices = truthVertices.second;
414 m_monPlots->fill(*vertices, primaryvertex, truthHSVertices, truthPUVertices, actualMu, beamSpotWeight);
418 m_monPlots->fill(*vertices, truthMu, actualMu, beamSpotWeight);
435 const auto& stdVertexContainer = truthVrt->stdcont();
437 auto findVtx = std::find_if(stdVertexContainer.rbegin(), stdVertexContainer.rend(), acceptTruthVertex);
438 truthVertex = (findVtx == stdVertexContainer.rend()) ? nullptr : *findVtx;
442 if (not truthVertex)
ATH_MSG_INFO (
"Truth vertex did not pass cuts");
447 unsigned int nSelectedTruthTracks(0), nSelectedRecoTracks(0), nSelectedMatchedTracks(0), nAssociatedTruth(0), nMissingAssociatedTruth(0), nTruths(0);
456 std::map<const xAOD::TruthParticle*, std::vector<const xAOD::TrackParticle*>> cacheTruthMatching {};
458 std::vector<const xAOD::TrackParticle*> selectedTracks {};
459 selectedTracks.reserve(tracks->
size());
460 unsigned int nTrackTOT = 0;
461 unsigned int nTrackCentral = 0;
462 unsigned int nTrackPt1GeV = 0;
463 for (
const auto *
const thisTrack: *tracks) {
469 selectedTracks.push_back(thisTrack);
471 nSelectedRecoTracks++;
475 if (thisTrack->pt() >= (1 * Gaudi::Units::GeV))
477 if (std::abs(thisTrack->eta()) < 2.5)
480 m_monPlots->fill(*thisTrack, puEvents, nVertices, beamSpotWeight);
482 float prob = getMatchingProbability(*thisTrack);
485 if (associatedTruth) {
496 nSelectedMatchedTracks++;
497 bool truthIsFromB =
false;
501 m_monPlots->fill(*thisTrack, *associatedTruth, truthIsFromB, puEvents, beamSpotWeight);
505 const bool isAssociatedTruth = associatedTruth !=
nullptr;
506 const bool isFake = not std::isnan(prob) ? (prob <
m_lowProb) :
true;
508 if(!isAssociatedTruth) nMissingAssociatedTruth++;
509 m_monPlots->fillFakeRate(*thisTrack, isFake, puEvents, beamSpotWeight);
515 if (isAssociatedTruth) {
520 auto cachedAssoc = cacheTruthMatching.find(associatedTruth);
522 if (cachedAssoc == cacheTruthMatching.end()) {
524 cacheTruthMatching[associatedTruth] = {thisTrack};
528 cachedAssoc->second.push_back(thisTrack);
533 m_monPlots->fillNtuple(*thisTrack, primaryvertex);
541 for (
auto& cachedAssoc: cacheTruthMatching) {
548 std::sort(cachedAssoc.second.begin(), cachedAssoc.second.end(),
554 for (
int itrack = 0; itrack < (int) cachedAssoc.second.size(); itrack++) {
558 m_monPlots->fillNtuple(*thisTrack, *thisTruth, primaryvertex, itrack);
564 m_monPlots->fill(nTrackTOT, nTrackCentral, nTrackPt1GeV, truthMu, actualMu, nVertices, beamSpotWeight);
570 m_truthCutFlow.merge(std::move(tmp_truth_cutflow));
576 for (
int itruth = 0; itruth < (int) truthParticlesVec.size(); itruth++) {
583 ++nSelectedTruthTracks;
584 bool isEfficient(
false);
585 float matchingProbability{};
589 auto cachedAssoc = cacheTruthMatching.find(thisTruth);
591 if (cachedAssoc == cacheTruthMatching.end()) {
593 cacheTruthMatching[thisTruth] = {};
595 m_monPlots->fillDuplicate(*thisTruth, cacheTruthMatching[thisTruth], beamSpotWeight);
602 for (
const auto& thisTrack: selectedTracks) {
604 if (associatedTruth && associatedTruth == thisTruth) {
605 float prob = getMatchingProbability(*thisTrack);
606 if (not std::isnan(prob) && prob >
m_lowProb) {
607 matchingProbability = prob;
609 matchedTrack = thisTrack;
616 <<thisTruth->
py()/ Gaudi::Units::GeV<<
","<<thisTruth->
pz()/ Gaudi::Units::GeV<<
","
617 <<thisTruth->
e()/ Gaudi::Units::GeV<<
","<<thisTruth->
pt()/ Gaudi::Units::GeV<<
","
618 <<thisTruth->
eta()<<
","<<thisTruth->
phi()<<
","<<thisTruth->
m()/ Gaudi::Units::GeV<<
","
619 <<puEvents<<
","<<nVertices<<
","<<matchingProbability<<std::endl;
622 ATH_MSG_ERROR(
"An error occurred: Truth particle for tracking efficiency calculation is a nullptr");
624 else if (isEfficient && !matchedTrack){
625 ATH_MSG_ERROR(
"Something went wrong - we log a truth particle as reconstructed, but the reco track is a nullptr! Bailing out... ");
628 ATH_MSG_DEBUG(
"Filling efficiency plots info monitoring plots");
629 m_monPlots->fillEfficiency(*thisTruth, matchedTrack, isEfficient, truthMu, actualMu, beamSpotWeight);
631 ATH_MSG_DEBUG(
"Filling technical efficiency plots info monitoring plots");
635 m_monPlots->fillTechnicalEfficiency(*thisTruth, isEfficient,
636 truthMu, actualMu, beamSpotWeight);
639 ATH_MSG_DEBUG(
"Cannot fill technical efficiency. Missing si hit information for truth particle.");
657 if (nSelectedRecoTracks == nMissingAssociatedTruth) {
662 ATH_MSG_DEBUG(nAssociatedTruth <<
" tracks out of " << tracks->
size() <<
" had associated truth.");
681 return StatusCode::SUCCESS;
707 ATH_MSG_DEBUG(
"Selected by pileup switch decoration requested from a truth particle but not available");
714 for (
const auto& thisTruth: truthParticles) {
723 std::vector<HistData> hists =
m_monPlots->retrieveBookedHistograms();
724 for (
const auto& hist : hists) {
728 std::vector<EfficiencyData> effs =
m_monPlots->retrieveBookedEfficiencies();
729 for (
auto& eff : effs) {
736 std::vector<TreeData> trees =
m_monPlots->retrieveBookedTrees();
737 for (
auto& t : trees) {
742 return StatusCode::SUCCESS;
767 return StatusCode::SUCCESS;
770const std::vector<const xAOD::TruthParticle*>
773 std::vector<const xAOD::TruthParticle*> tempVec {};
780 if (not truthParticleContainer.
isValid()) {
783 tempVec.insert(tempVec.begin(), truthParticleContainer->begin(), truthParticleContainer->end());
793 const auto& links =
event->truthParticleLinks();
794 tempVec.reserve(event->nTruthParticles());
795 for (
const auto& link : links) {
797 tempVec.push_back(*link);
806 if (truthPileupEventContainer.
isValid()) {
807 const unsigned int nPileup = truthPileupEventContainer->size();
808 tempVec.reserve(nPileup * 200);
809 for (
unsigned int i(0); i != nPileup; ++i) {
810 const auto *eventPileup = truthPileupEventContainer->at(i);
812 int ntruth = eventPileup->nTruthParticles();
813 ATH_MSG_VERBOSE(
"Adding " << ntruth <<
" truth particles from TruthPileupEvents container");
814 const auto& links = eventPileup->truthParticleLinks();
815 for (
const auto& link : links) {
817 tempVec.push_back(*link);
832std::pair<std::vector<const xAOD::TruthVertex*>, std::vector<const xAOD::TruthVertex*>>
835 std::vector<const xAOD::TruthVertex*> truthHSVertices = {};
836 truthHSVertices.reserve(5);
837 std::vector<const xAOD::TruthVertex*> truthPUVertices = {};
838 truthPUVertices.reserve(100);
861 if (truthEventContainer.
isValid()) {
862 for (
const auto *
const evt : *truthEventContainer) {
863 truthVtx = evt->signalProcessVertex();
865 truthHSVertices.push_back(truthVtx);
879 if (truthPileupEventContainer.
isValid()) {
880 for (
const auto *
const evt : *truthPileupEventContainer) {
885 size_t i_vtx = 0;
size_t n_vtx = evt->nTruthVertices();
886 while(!truthVtx && i_vtx<n_vtx){
887 truthVtx = evt->truthVertex(i_vtx);
892 truthPUVertices.push_back(truthVtx);
902 return {std::move(truthHSVertices), std::move(truthPUVertices)};
912 std::vector<int>& cutFlow) {
914 if (cutFlow.empty()) {
915 names.emplace_back(
"preCut");
916 cutFlow.push_back(0);
917 for (
unsigned int i = 0; i != accept.getNCuts(); ++i) {
918 cutFlow.push_back(0);
919 names.push_back((std::string) accept.getCutName(i));
924 bool cutPositive =
true;
925 for (
unsigned int i = 0; i != (accept.getNCuts() + 1); ++i) {
929 if (accept.getCutResult(i)) {
938 double absEta = std::abs(truth.
eta());
941 ATH_MSG_INFO(
"Requesting cut value outside of configured eta range: clamping eta = "
942 << std::abs(truth.
eta()) <<
" to eta= " << absEta);
945 const auto pVal = std::lower_bound(
m_etaBins.value().begin(),
m_etaBins.value().end(), absEta);
946 const int bin = std::distance(
m_etaBins.value().begin(), pVal) - 1;
947 ATH_MSG_DEBUG(
"Checking (abs(eta)/bin) = (" << absEta <<
"," <<
bin <<
")");
953 const std::vector<const xAOD::TruthParticle*>& truthParticles,
956 float beamSpotWeight) {
961 if (truthParticles.empty()) {
962 ATH_MSG_WARNING(
"No entries in TruthParticles truth particle container. Skipping jet plots.");
963 return StatusCode::SUCCESS;
969 return StatusCode::SUCCESS;
974 for (
const xAOD::Jet *
const thisJet: *jets) {
982 isBjet = (btagLabel(*thisJet) == 5);
990 if (not el.isValid())
continue;
996 if (std::find(truthParticles.begin(), truthParticles.end(), truth) == truthParticles.end()) {
1007 if(!accept)
continue;
1009 bool isEfficient(
false);
1011 for (
const auto *thisTrack: tracks) {
1017 if (associatedTruth and associatedTruth == truth) {
1018 float prob = getMatchingProbability(*thisTrack);
1019 if (not std::isnan(prob) && prob >
m_lowProb) {
1026 bool truthIsFromB =
false;
1028 truthIsFromB =
true;
1030 m_monPlots->fillEfficiency(*truth, *thisJet, isEfficient, isBjet, truthIsFromB, beamSpotWeight);
1040 if (thisJet->p4().DeltaR(thisTrack->p4()) >
m_maxTrkJetDR) {
1044 float prob = getMatchingProbability(*thisTrack);
1045 if(std::isnan(prob)) prob = 0.0;
1048 const bool isFake = (associatedTruth && prob <
m_lowProb);
1049 bool truthIsFromB =
false;
1051 truthIsFromB =
true;
1053 m_monPlots->fill(*thisTrack, *thisJet, isBjet, isFake, truthIsFromB, beamSpotWeight);
1054 if (associatedTruth){
1055 m_monPlots->fillFakeRate(*thisTrack, *thisJet, isFake, isBjet, truthIsFromB, beamSpotWeight);
1061 return StatusCode::SUCCESS;
1066 const float jetPt =
jet.pt();
1067 const float jetEta = std::abs(
jet.eta());
#define ATH_CHECK
Evaluate an expression and check for errors.
#define ATH_MSG_DEBUG(x,...)
#define ATH_MSG_ERROR(x,...)
#define ATH_MSG_WARNING(x,...)
#define ATH_MSG_VERBOSE(x,...)
#define ATH_MSG_INFO(x,...)
Helper class to provide constant type-safe access to aux data.
virtual void lock()=0
Interface to allow an object to lock itself when made const in SG.
header file for class of same name
bool inRange(const double *boundaries, const double value, const double tolerance=0.02)
ServiceHandle< StoreGateSvc > & evtStore()
size_type size() const noexcept
Returns the number of elements in the collection.
ElementLink implementation for ROOT usage.
static void neededTrackParticleDecorations(std::vector< std::string > &decorations)
const xAOD::TruthParticle * getTruth(const xAOD::TrackParticle *const trackParticle)
Helper class to provide constant type-safe access to aux data.
bool isAvailable(const ELT &e) const
Test to see if this variable exists in the store.
Handle class for reading a decoration on an object.
bool isAvailable()
Test to see if this variable exists in the store, for the referenced object.
virtual bool isValid() override final
Can the handle be successfully dereferenced?
const_pointer_type cptr()
Dereference the pointer.
const_pointer_type get() const
Dereference the pointer, but don't cache anything.
@ IS_SIMULATION
true: simulation, false: data
uint64_t eventNumber() const
The current event's event number.
virtual double m() const override final
The mass of the particle.
int pdgId() const
PDG ID code.
float px() const
The x component of the particle's momentum.
virtual double e() const override final
The total energy of the particle.
virtual double pt() const override final
The transverse momentum ( ) of the particle.
float py() const
The y component of the particle's momentum.
virtual double eta() const override final
The pseudorapidity ( ) of the particle.
virtual double phi() const override final
The azimuthal angle ( ) of the particle.
virtual FourMom_t p4() const override final
The full 4-momentum of the particle.
float pz() const
The z component of the particle's momentum.
float z() const
Vertex longitudinal distance along the beam line form the origin.
float y() const
Vertex y displacement.
float x() const
Vertex x displacement.
void addReadDecoratorHandleKeys(T_Parent &parent, const SG::ReadHandleKey< T_Cont > &container_key, const std::string &prefix, const std::vector< std::string > &decor_names, std::vector< SG::ReadDecorHandleKey< T_Cont > > &decor_out)
float safelyGetEta(const T &pTrk, const float safePtThreshold=0.1)
Safely get eta.
unsigned int binIndex(const T &val, const std::vector< T > &partitions)
general utility function to return bin index given a value and the upper endpoints of each bin
HardScatterType classifyHardScatter(const xAOD::VertexContainer &vxContainer)
void sort(typename DataModel_detail::iterator< DVL > beg, typename DataModel_detail::iterator< DVL > end)
Specialization of sort for DataVector/List.
Jet_v1 Jet
Definition of the current "jet version".
EventInfo_v1 EventInfo
Definition of the latest event info version.
TruthVertex_v1 TruthVertex
Typedef to implementation.
TrackParticle_v1 TrackParticle
Reference the current persistent version:
Vertex_v1 Vertex
Define the latest version of the vertex class.
TruthEvent_v1 TruthEvent
Typedef to implementation.
TruthParticle_v1 TruthParticle
Typedef to implementation.
TrackParticleContainer_v1 TrackParticleContainer
Definition of the current "TrackParticle container version".
JetContainer_v1 JetContainer
Definition of the current "jet container version".
implementation file for function of same name
helper struct - steer the configuration from the parent tool's side
bool doTrkInJetPlots_matched_bjets
bool doHitsFakeTracksPlots
int detailLevel
detail level (kept for compatibility)
bool doNtupleTruthToReco
Ntuple functionality.
bool doEfficienciesPerAuthor
per author plots
bool doHitsRecoTracksPlots
bool doHitsRecoTracksPlotsPerAuthor
bool doEffPlots
Efficiency and duplicate plots - require truth, optionally matching reco.
bool doTrkInJetPlots_bjets
bool doTrkInJetPlots_matched
bool doFakePlots
Fake plots.
bool doResolutionPlotPrim
Resolution and "matched track" plots - filled if both reco and truth exist.
bool doTrkInJetPlots_fake_bjets
bool doVertexTruthMatchingPlots
Vertexing plots - truth requirement.
bool doVertexPlots
Vertexing plots - no truth requirement.
bool doHardScatterVertexTruthMatchingPlots
bool doTrkInJetPlots_truthFromB
bool doTrkInJetPlots_fake
bool doTrackParametersPerAuthor
bool doTrackParameters
Plots for (selected) tracks, not necessarily truth matched.
bool doTrkInJetPlots
Plots for tracks in jets.
bool doResolutionPlotPrim_truthFromB
bool doHardScatterVertexPlots
bool doResolutionPlotSecd
bool doResolutionsPerAuthor
bool doHitsMatchedTracksPlots