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;
317 if (not trackHandle.
isValid()) {
319 return StatusCode::FAILURE;
330 ATH_MSG_WARNING(
"Shouldn't happen - EventInfo is buggy, setting mu to 0");
338 return StatusCode::SUCCESS;
342 std::vector<const xAOD::TruthParticle*> truthParticlesVec =
getTruthParticles(ctx);
351 unsigned int truthMu = 0;
354 truthMu =
static_cast<int>( truthPileupEventContainer->size() );
356 if(pie.
isValid()) actualMu = pie->actualInteractionsPerCrossing();
359 float puEvents = truthMu>0 ? truthMu : actualMu;
362 unsigned int nVertices = 0;
363 float beamSpotWeight = 1;
366 ATH_MSG_DEBUG(
"Getting number of pu interactings per event");
372 nVertices = not vertices->empty() ? vertices->size() : 0;
373 beamSpotWeight = pie->beamSpotWeight();
374 ATH_MSG_DEBUG(
"beamSpotWeight is equal to " << beamSpotWeight);
378 if(readDecorHandle.
isAvailable()) prwWeight = readDecorHandle(*pie);
379 ATH_MSG_DEBUG(
"Applying pileup weight equal to " << prwWeight);
380 beamSpotWeight *= prwWeight;
383 if (vertices.
isValid() and not vertices->empty()) {
384 ATH_MSG_DEBUG(
"Number of vertices retrieved for this event " << vertices->size());
387 for (
const auto& vertex : *vertices){
389 primaryvertex = vertex;
400 ATH_MSG_DEBUG(
"Failed to find a hard scatter vertex in this event.");
406 std::pair<std::vector<const xAOD::TruthVertex*>, std::vector<const xAOD::TruthVertex*>> truthVertices =
getTruthVertices(ctx);
407 std::vector<const xAOD::TruthVertex*> truthHSVertices = truthVertices.first;
408 std::vector<const xAOD::TruthVertex*> truthPUVertices = truthVertices.second;
415 m_monPlots->fill(*vertices, primaryvertex, truthHSVertices, truthPUVertices, actualMu, beamSpotWeight);
419 m_monPlots->fill(*vertices, truthMu, actualMu, beamSpotWeight);
436 const auto& stdVertexContainer = truthVrt->stdcont();
438 auto findVtx = std::find_if(stdVertexContainer.rbegin(), stdVertexContainer.rend(), acceptTruthVertex);
439 truthVertex = (findVtx == stdVertexContainer.rend()) ? nullptr : *findVtx;
443 if (not truthVertex)
ATH_MSG_INFO (
"Truth vertex did not pass cuts");
448 unsigned int nSelectedTruthTracks(0), nSelectedRecoTracks(0), nSelectedMatchedTracks(0), nAssociatedTruth(0), nMissingAssociatedTruth(0), nTruths(0);
457 std::map<const xAOD::TruthParticle*, std::vector<const xAOD::TrackParticle*>> cacheTruthMatching {};
459 std::vector<const xAOD::TrackParticle*> selectedTracks {};
460 selectedTracks.reserve(tracks->
size());
461 unsigned int nTrackTOT = 0;
462 unsigned int nTrackCentral = 0;
463 unsigned int nTrackPt1GeV = 0;
464 for (
const auto *
const thisTrack: *tracks) {
470 selectedTracks.push_back(thisTrack);
472 nSelectedRecoTracks++;
476 if (thisTrack->pt() >= (1 * Gaudi::Units::GeV))
478 if (std::abs(thisTrack->eta()) < 2.5)
481 m_monPlots->fill(*thisTrack, puEvents, nVertices, beamSpotWeight);
483 float prob = getMatchingProbability(*thisTrack);
486 if (associatedTruth) {
497 nSelectedMatchedTracks++;
498 bool truthIsFromB =
false;
502 m_monPlots->fill(*thisTrack, *associatedTruth, truthIsFromB, puEvents, beamSpotWeight);
506 const bool isAssociatedTruth = associatedTruth !=
nullptr;
507 const bool isFake = not std::isnan(prob) ? (prob <
m_lowProb) :
true;
509 if(!isAssociatedTruth) nMissingAssociatedTruth++;
510 m_monPlots->fillFakeRate(*thisTrack, isFake, puEvents, beamSpotWeight);
516 if (isAssociatedTruth) {
521 auto cachedAssoc = cacheTruthMatching.find(associatedTruth);
523 if (cachedAssoc == cacheTruthMatching.end()) {
525 cacheTruthMatching[associatedTruth] = {thisTrack};
529 cachedAssoc->second.push_back(thisTrack);
534 m_monPlots->fillNtuple(*thisTrack, primaryvertex);
542 for (
auto& cachedAssoc: cacheTruthMatching) {
549 std::sort(cachedAssoc.second.begin(), cachedAssoc.second.end(),
555 for (
int itrack = 0; itrack < (int) cachedAssoc.second.size(); itrack++) {
559 m_monPlots->fillNtuple(*thisTrack, *thisTruth, primaryvertex, itrack);
565 m_monPlots->fill(nTrackTOT, nTrackCentral, nTrackPt1GeV, truthMu, actualMu, nVertices, beamSpotWeight);
571 m_truthCutFlow.merge(std::move(tmp_truth_cutflow));
577 for (
int itruth = 0; itruth < (int) truthParticlesVec.size(); itruth++) {
584 ++nSelectedTruthTracks;
585 bool isEfficient(
false);
586 float matchingProbability{};
590 auto cachedAssoc = cacheTruthMatching.find(thisTruth);
592 if (cachedAssoc == cacheTruthMatching.end()) {
594 cacheTruthMatching[thisTruth] = {};
596 m_monPlots->fillDuplicate(*thisTruth, cacheTruthMatching[thisTruth], beamSpotWeight);
603 for (
const auto& thisTrack: selectedTracks) {
605 if (associatedTruth && associatedTruth == thisTruth) {
606 float prob = getMatchingProbability(*thisTrack);
607 if (not std::isnan(prob) && prob >
m_lowProb) {
608 matchingProbability = prob;
610 matchedTrack = thisTrack;
617 <<thisTruth->
py()/ Gaudi::Units::GeV<<
","<<thisTruth->
pz()/ Gaudi::Units::GeV<<
","
618 <<thisTruth->
e()/ Gaudi::Units::GeV<<
","<<thisTruth->
pt()/ Gaudi::Units::GeV<<
","
619 <<thisTruth->
eta()<<
","<<thisTruth->
phi()<<
","<<thisTruth->
m()/ Gaudi::Units::GeV<<
","
620 <<puEvents<<
","<<nVertices<<
","<<matchingProbability<<std::endl;
623 ATH_MSG_ERROR(
"An error occurred: Truth particle for tracking efficiency calculation is a nullptr");
625 else if (isEfficient && !matchedTrack){
626 ATH_MSG_ERROR(
"Something went wrong - we log a truth particle as reconstructed, but the reco track is a nullptr! Bailing out... ");
629 ATH_MSG_DEBUG(
"Filling efficiency plots info monitoring plots");
630 m_monPlots->fillEfficiency(*thisTruth, matchedTrack, isEfficient, truthMu, actualMu, beamSpotWeight);
632 ATH_MSG_DEBUG(
"Filling technical efficiency plots info monitoring plots");
636 m_monPlots->fillTechnicalEfficiency(*thisTruth, isEfficient,
637 truthMu, actualMu, beamSpotWeight);
640 ATH_MSG_DEBUG(
"Cannot fill technical efficiency. Missing si hit information for truth particle.");
658 if (nSelectedRecoTracks == nMissingAssociatedTruth) {
663 ATH_MSG_DEBUG(nAssociatedTruth <<
" tracks out of " << tracks->
size() <<
" had associated truth.");
682 return StatusCode::SUCCESS;
708 ATH_MSG_DEBUG(
"Selected by pileup switch decoration requested from a truth particle but not available");
715 for (
const auto& thisTruth: truthParticles) {
724 std::vector<HistData> hists =
m_monPlots->retrieveBookedHistograms();
725 for (
const auto& hist : hists) {
729 std::vector<EfficiencyData> effs =
m_monPlots->retrieveBookedEfficiencies();
730 for (
auto& eff : effs) {
737 std::vector<TreeData> trees =
m_monPlots->retrieveBookedTrees();
738 for (
auto& t : trees) {
743 return StatusCode::SUCCESS;
768 return StatusCode::SUCCESS;
771const std::vector<const xAOD::TruthParticle*>
774 std::vector<const xAOD::TruthParticle*> tempVec {};
781 if (not truthParticleContainer.
isValid()) {
784 tempVec.insert(tempVec.begin(), truthParticleContainer->begin(), truthParticleContainer->end());
794 const auto& links =
event->truthParticleLinks();
795 tempVec.reserve(event->nTruthParticles());
796 for (
const auto& link : links) {
798 tempVec.push_back(*link);
807 if (truthPileupEventContainer.
isValid()) {
808 const unsigned int nPileup = truthPileupEventContainer->size();
809 tempVec.reserve(nPileup * 200);
810 for (
unsigned int i(0); i != nPileup; ++i) {
811 const auto *eventPileup = truthPileupEventContainer->at(i);
813 int ntruth = eventPileup->nTruthParticles();
814 ATH_MSG_VERBOSE(
"Adding " << ntruth <<
" truth particles from TruthPileupEvents container");
815 const auto& links = eventPileup->truthParticleLinks();
816 for (
const auto& link : links) {
818 tempVec.push_back(*link);
833std::pair<std::vector<const xAOD::TruthVertex*>, std::vector<const xAOD::TruthVertex*>>
836 std::vector<const xAOD::TruthVertex*> truthHSVertices = {};
837 truthHSVertices.reserve(5);
838 std::vector<const xAOD::TruthVertex*> truthPUVertices = {};
839 truthPUVertices.reserve(100);
862 if (truthEventContainer.
isValid()) {
863 for (
const auto *
const evt : *truthEventContainer) {
864 truthVtx = evt->signalProcessVertex();
866 truthHSVertices.push_back(truthVtx);
880 if (truthPileupEventContainer.
isValid()) {
881 for (
const auto *
const evt : *truthPileupEventContainer) {
886 size_t i_vtx = 0;
size_t n_vtx = evt->nTruthVertices();
887 while(!truthVtx && i_vtx<n_vtx){
888 truthVtx = evt->truthVertex(i_vtx);
893 truthPUVertices.push_back(truthVtx);
903 return {std::move(truthHSVertices), std::move(truthPUVertices)};
913 std::vector<int>& cutFlow) {
915 if (cutFlow.empty()) {
916 names.emplace_back(
"preCut");
917 cutFlow.push_back(0);
918 for (
unsigned int i = 0; i != accept.getNCuts(); ++i) {
919 cutFlow.push_back(0);
920 names.push_back((std::string) accept.getCutName(i));
925 bool cutPositive =
true;
926 for (
unsigned int i = 0; i != (accept.getNCuts() + 1); ++i) {
930 if (accept.getCutResult(i)) {
939 double absEta = std::abs(truth.
eta());
942 ATH_MSG_INFO(
"Requesting cut value outside of configured eta range: clamping eta = "
943 << std::abs(truth.
eta()) <<
" to eta= " << absEta);
946 const auto pVal = std::lower_bound(
m_etaBins.value().begin(),
m_etaBins.value().end(), absEta);
947 const int bin = std::distance(
m_etaBins.value().begin(), pVal) - 1;
948 ATH_MSG_DEBUG(
"Checking (abs(eta)/bin) = (" << absEta <<
"," <<
bin <<
")");
954 const std::vector<const xAOD::TruthParticle*>& truthParticles,
957 float beamSpotWeight) {
962 if (truthParticles.empty()) {
963 ATH_MSG_WARNING(
"No entries in TruthParticles truth particle container. Skipping jet plots.");
964 return StatusCode::SUCCESS;
970 return StatusCode::SUCCESS;
975 for (
const xAOD::Jet *
const thisJet: *jets) {
983 isBjet = (btagLabel(*thisJet) == 5);
991 if (not el.isValid())
continue;
1000 if(!accept)
continue;
1002 bool isEfficient(
false);
1004 for (
const auto *thisTrack: tracks) {
1010 if (associatedTruth and associatedTruth == truth) {
1011 float prob = getMatchingProbability(*thisTrack);
1012 if (not std::isnan(prob) && prob >
m_lowProb) {
1019 bool truthIsFromB =
false;
1021 truthIsFromB =
true;
1023 m_monPlots->fillEfficiency(*truth, *thisJet, isEfficient, isBjet, truthIsFromB, beamSpotWeight);
1033 if (thisJet->p4().DeltaR(thisTrack->p4()) >
m_maxTrkJetDR) {
1037 float prob = getMatchingProbability(*thisTrack);
1038 if(std::isnan(prob)) prob = 0.0;
1041 const bool isFake = (associatedTruth && prob <
m_lowProb);
1042 bool truthIsFromB =
false;
1044 truthIsFromB =
true;
1046 m_monPlots->fill(*thisTrack, *thisJet, isBjet, isFake, truthIsFromB, beamSpotWeight);
1047 if (associatedTruth){
1048 m_monPlots->fillFakeRate(*thisTrack, *thisJet, isFake, isBjet, truthIsFromB, beamSpotWeight);
1054 return StatusCode::SUCCESS;
1059 const float jetPt =
jet.pt();
1060 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