636 {
638 ATH_MSG_WARNING(
"Requested soft track element links, but no track selection tool supplied.");
639 }
641
642 if(
helper.map().empty()) {
643 ATH_MSG_WARNING(
"Incomplete association map received. Cannot rebuild MET.");
644 ATH_MSG_WARNING(
"Note: ColumnarMETMaker should only be run on events containing at least one PV");
645 return StatusCode::SUCCESS;
646 }
647
648 if(doJetJVT &&
m_JvtWP ==
"None"){
650 doJetJVT = false;
651 }
652
654 if(!metSoftClus && !metSoftTrk) {
655 ATH_MSG_WARNING(
"Neither soft cluster nor soft track term has been supplied!");
656 return StatusCode::SUCCESS;
657 }
659 std::optional<columnar::MetHelpers::ObjectWeightHandle<columnar::ContainerId::mutableMet,columnar::ContainerId::jet>> metSoftClusLinks;
660 if(metSoftClus) {
662 if(!coreSoftClus) {
663 ATH_MSG_ERROR(
"Soft cluster term provided without a core term!");
664 return StatusCode::FAILURE;
665 }
671
672
673
674 if(softConst && acc_softConst.isAvailable(*coreSoftClus.getXAODObject())) {
675 for(const auto& constit : acc_softConst(*coreSoftClus.getXAODObject())) {
676 softConst->push_back(*constit);
677 }
678 ATH_MSG_DEBUG(softConst->size() <<
" soft constituents from core term");
679 }
680 }
681 std::optional<columnar::MetHelpers::ObjectWeightHandle<columnar::ContainerId::mutableMet,columnar::ContainerId::jet>> metSoftTrkLinks;
682 if(metSoftTrk) {
684 if(!coreSoftTrk) {
685 ATH_MSG_ERROR(
"Soft track term provided without a core term!");
686 return StatusCode::FAILURE;
687 }
693 if(softConst && acc_softConst.isAvailable(*coreSoftTrk.getXAODObject()) && !
m_doPFlow && !
m_doSoftTruth) {
694 for(const auto& constit : acc_softConst(*coreSoftTrk.getXAODObject())) {
695 softConst->push_back(*constit);
696 }
697 ATH_MSG_DEBUG(softConst->size() <<
" soft constituents from trk core term");
698 }
699 }
700
701 columnar::MetHelpers::ObjectWeightHandle<columnar::ContainerId::mutableMet,columnar::ContainerId::jet> metJetWeights(*
this,
m_jetOutputMetWeightDecRegular,metJet,jets);
702
703
704
705
706
707
708
709
711 for(auto jet : jets) {
712 auto originalJet = jetsOriginals.getOriginal(jet);
713 auto assoc =
helper.getJetAssociation(originalJet);
716 continue;
717 }
718
720
724 }
725 else
726 ATH_MSG_ERROR(
"No nominal calibrated jet available for jet " << jet <<
". Cannot simplify overlap removal!");
727 }
729
731 bool JVT_reject(false);
732 bool isMuFSRJet(false);
733
734
736 JVT_reject = true;
737
738 if(doJetJVT) {
739
741 ATH_MSG_VERBOSE(
"Jet " << (JVT_reject ?
"fails" :
"passes") <<
" JVT selection");
742 }
743
744
746 bool hardJet(false);
748 bool caloverlap = false;
749 caloverlap = calvec.
ce()>0;
750 ATH_MSG_DEBUG(
"Jet " << jet <<
" is " << ( caloverlap ?
"" :
"non-") <<
"overlapping");
751
753 for(
const auto object :
m_assocAcc.objects(*assoc)) {
754
756 }
757 }
758
760 double constSF(1);
762 constjet =
m_assocAcc.getAlternateConstVec(*assoc);
763 } else {
768 calvec *= constSF;
769 }
770 double jpx = constjet.Px();
771 double jpy = constjet.Py();
772 double jpt = constjet.Pt();
773 double opx = jpx - calvec.
cpx();
774 double opy = jpy - calvec.
cpy();
775
778
780 if(!met_muonEloss) {
781 ATH_MSG_WARNING(
"Attempted to apply muon Eloss correction, but corresponding MET term does not exist!");
782 return StatusCode::FAILURE;
783 }
784 }
785
786 float total_eloss(0);
788 std::vector<columnar::MuonId> muons_in_jet;
789 std::vector<columnar::ElectronId> electrons_in_jet;
790 bool passJetForEl=false;
794 return StatusCode::FAILURE;
795 }
798 ATH_MSG_ERROR(
"Invalid element link to ghost muon! Quitting.");
799 return StatusCode::FAILURE;
800 }
801 muons_in_jet.push_back(*
static_cast<const xAOD::Muon*
>(*el));
802 }
803 }
804 for(
const auto obj :
m_assocAcc.objects(*assoc)) {
805 if(!obj) continue;
807 auto mu_test =
obj.tryGetVariant<columnar::ContainerId::muon>().
value();
812 }
813 if(
helper.objSelected(mu_test)) {
814 muons_in_jet.push_back(mu_test);
816 }
817 }
819 auto el_test =
obj.tryGetVariant<columnar::ContainerId::electron>().
value();
823 }
824 if(
helper.objSelected(*assoc,el_test)){
826 electrons_in_jet.push_back(el_test);
828 }
829 }
830 }
831 }
834 float jet_ORtrk_sumpt =
helper.overlapTrkVec(*assoc).sumpt();
835 float jet_all_trk_pt = initialTrkMom.
sumpt();
836 float jet_unique_trk_pt = jet_all_trk_pt - jet_ORtrk_sumpt;
839 for(const auto& elec : electrons_in_jet) {
841 el_trkvec +=
m_assocAcc.trkVec(*assoc,&elec.getXAODObject());
842 }
843 float el_cal_pt = el_calvec.
cpt();
844 float el_trk_pt = el_trkvec.
cpt();
846 << " jetalltrk: " << jet_all_trk_pt
847 << " jetORtrk: " << jet_ORtrk_sumpt
848 << " electrk-jetORtrk: " << (el_trk_pt-jet_ORtrk_sumpt)
849 << " elec cal: " << el_cal_pt
850 << " jetalltrk-electrk: " << (jet_all_trk_pt-el_trk_pt)
851 << " jetalltrk-jetORtrk: " << (jet_all_trk_pt-jet_ORtrk_sumpt) );
852
853
854
855 if(el_trk_pt>1
e-9 && jet_unique_trk_pt>10.0e3) passJetForEl=
true;
856 }
857
858 for(auto mu_in_jet : muons_in_jet) {
859 float mu_Eloss =
acc_Eloss(mu_in_jet.getXAODObject());
860
861 if(!JVT_reject) {
863
864 float mu_id_pt = mu_in_jet.getXAODObject().trackParticle(xAOD::Muon::InnerDetectorTrackParticle) ? mu_in_jet.getXAODObject().trackParticle(xAOD::Muon::InnerDetectorTrackParticle)->pt() : 0.;
866
867
868 if(0.9999*mu_id_pt>jet_trk_sumpt)
869 jet_trk_sumpt+=mu_id_pt;
874 if(jet_from_muon) {
876 JVT_reject = true;
877 }
878 }
879
881
882 float mu_id_pt = mu_in_jet.getXAODObject().trackParticle(xAOD::Muon::InnerDetectorTrackParticle) ? mu_in_jet.getXAODObject().trackParticle(xAOD::Muon::InnerDetectorTrackParticle)->pt() : 0.;
884
885 if(0.9999*mu_id_pt>jet_trk_sumpt)
886 jet_trk_sumpt+=mu_id_pt;
888
889 float jet_psE = 0.;
894 } else {
895 ATH_MSG_ERROR(
"Jet PS fraction or sampling energy must be available to calculate MET with doSetMuonJetEMScale");
896 return StatusCode::FAILURE;
897 }
898
901 ATH_MSG_VERBOSE(
"Jet has trk sumpt " << jet_trk_sumpt <<
", trk N " << jet_trk_N <<
", PS E " << jet_psE <<
", width " <<
m_acc_width(jet) <<
", emfrac " <<
m_acc_emf(jet));
902
903 if(jet_from_muon) {
904 ATH_MSG_VERBOSE(
"Jet is from muon -- set to EM scale and subtract Eloss.");
905
906
907 ATH_MSG_VERBOSE(
"Jet e: " << constjet.E() <<
", mu Eloss: " << mu_Eloss);
908 float elosscorr = mu_Eloss >= constjet.e() ? 0. : 1.-mu_Eloss/constjet.e();
909
910
911 opx *= elosscorr;
912 opy *= elosscorr;
913 ATH_MSG_VERBOSE(
" Jet eloss factor " << elosscorr <<
", final pt: " << sqrt(opx*opx+opy*opy));
914
915 isMuFSRJet = true;
916 }
917 }
918 }
919
920 switch(mu_in_jet.getXAODObject().energyLossType()) {
921 case xAOD::Muon::Parametrized:
922 case xAOD::Muon::MOP:
923 case xAOD::Muon::Tail:
924 case xAOD::Muon::FSRcandidate:
925 case xAOD::Muon::NotIsolated:
926
927
928
929 total_eloss += mu_Eloss;
930 muons_selflags |= (1<<
m_assocAcc.findIndex(*assoc,mu_in_jet));
931 }
932 }
935
937
938 for(
size_t iKey = 0; iKey <
m_assocAcc.sizeCal(*assoc); iKey++) {
940 if(selector) mu_calovec +=
m_assocAcc.calVec(*assoc,iKey);
942 }
944 if(
m_muEloss) mu_calovec *= std::max<float>(0.,1-(total_eloss/mu_calovec.
ce()));
946
947
949 ATH_MSG_VERBOSE(
"Jet " << jet <<
" const pT after OR " << sqrt(opx*opx+opy*opy));
950 opx += mu_calovec.
cpx();
951 opy += mu_calovec.
cpy();
952 double opt = sqrt( opx*opx+opy*opy );
953 ATH_MSG_VERBOSE(
"Jet " << jet <<
" const pT diff after OR readding muon clusters " << opt-jpt);
954 double uniquefrac = 1. - (calvec.
ce() - mu_calovec.
ce()) / constjet.E();
955 ATH_MSG_VERBOSE(
"Jet constscale px, py, pt, E = " << jpx <<
", " << jpy <<
", " << jpt <<
", " << constjet.E() );
957 ATH_MSG_VERBOSE(
"Jet OR px, py, pt, E = " << opx <<
", " << opy <<
", " << opt <<
", " << constjet.E() - calvec.
ce() );
958
959 if(isMuFSRJet) {
960 if(!met_muonEloss){
961 ATH_MSG_ERROR(
"Attempted to apply muon Eloss correction, but corresponding MET term does not exist!");
962 return StatusCode::FAILURE;
963 }
965 continue;
966 }
967
968 if(selected && !JVT_reject) {
969 if(!caloverlap) {
970
971 hardJet = true;
972 if (!tracksForHardJets) {
975 else
977 }
978 }
980
981 hardJet = true;
982 if(!tracksForHardJets) {
986 else {
989 }
990 } else {
992 m_outputMetMomAcc.addParticle(metJet,uniquefrac*jpx,uniquefrac*jpy,uniquefrac*jpt);
993 else{
996 else
998 }
999 }
1000 }
1001 }
1002 }
1003
1004 if(hardJet){
1006 metJetWeights.emplace_back(jet, uniquefrac);
1007 } else {
1008 if(metSoftClus && !JVT_reject) {
1009
1012 metSoftClusLinks->emplace_back (jet, uniquefrac);
1014 }
1015
1016
1017
1018
1019
1020
1021 if(softConst) {
1022 for(size_t iConst=0; iConst<jet.getXAODObject().numConstituents(); ++iConst) {
1023 const IParticle* constit = jet.getXAODObject().rawConstituent(iConst);
1024 softConst->push_back(constit);
1025 }
1026 }
1027 }
1028 }
1029
1030 if(!metSoftTrk || (hardJet && !tracksForHardJets)) continue;
1031
1032
1033
1036 if(jettrkvec.
ce()>1
e-9) {
1037 jpx = jettrkvec.
cpx();
1038 jpy = jettrkvec.
cpy();
1039 jpt = jettrkvec.
sumpt();
1040 jettrkvec -= trkvec;
1041 opx = jettrkvec.
cpx();
1042 opy = jettrkvec.
cpy();
1044 ATH_MSG_VERBOSE(
"Jet track px, py, sumpt = " << jpx <<
", " << jpy <<
", " << jpt );
1045 ATH_MSG_VERBOSE(
"Jet OR px, py, sumpt = " << opx <<
", " << opy <<
", " << opt );
1046 } else {
1047 opx = opy =
opt = 0;
1049 }
1053
1054 if(!metSoftClus) {
1055 if (metSoftTrkLinks) metSoftTrkLinks->emplace_back (jet, uniquefrac);
1056 }
1057
1058
1059
1060
1061
1062
1064 std::vector<const IParticle*> jettracks;
1066 for(size_t iConst=0; iConst<jettracks.size(); ++iConst) {
1068 if (
acceptTrack(pTrk,pv)) softConst->push_back(pTrk);
1069 }
1070 }
1071 }
1072 }
1073
1074 ATH_MSG_DEBUG(
"Number of selected jets: " << metJetWeights.size());
1075
1076 if(metSoftTrk) {
1077 ATH_MSG_DEBUG(
"Number of softtrk jets: " << metSoftTrkLinks->size());
1078 }
1079
1080 if(metSoftClus) {
1081 ATH_MSG_DEBUG(
"Number of softclus jets: " << metSoftClusLinks->size());
1082 }
1083
1084 if(softConst)
ATH_MSG_DEBUG(softConst->size() <<
" soft constituents from core term + jets");
1085
1086 auto assoc =
helper.getMiscAssociation();
1087 if(!assoc) return StatusCode::SUCCESS;
1088
1089 if(metSoftTrk) {
1090
1091
1093 double opx = trkvec.
cpx();
1094 double opy = trkvec.
cpy();
1095 double osumpt = trkvec.
sumpt();
1096 ATH_MSG_VERBOSE(
"Misc track px, py, sumpt = " << opx <<
", " << opy <<
", " << osumpt );
1101 }
1102
1103 if(metSoftClus) {
1104
1105
1106 float total_eloss(0.);
1109 double opx = calvec.
cpx();
1110 double opy = calvec.
cpy();
1111 double osumpt = calvec.
sumpt();
1112 for(
const auto objId :
m_assocAcc.objects(*assoc)) {
1113 auto *
obj = objId.getXAODObject();
1117 if(
helper.objSelected(mu_test)) {
1119 switch(mu_test->energyLossType()) {
1120 case xAOD::Muon::Parametrized:
1121 case xAOD::Muon::MOP:
1122 case xAOD::Muon::Tail:
1123 case xAOD::Muon::FSRcandidate:
1124 case xAOD::Muon::NotIsolated:
1125
1126
1127
1128 total_eloss += mu_Eloss;
1129 muons_selflags |= (1<<
m_assocAcc.findIndex(*assoc,mu_test));
1130 }
1132 }
1133 }
1136
1138
1139 for(
size_t iKey = 0; iKey <
m_assocAcc.sizeCal(*assoc); iKey++) {
1142 <<
" this calvec E: " <<
m_assocAcc.calVec(*assoc,iKey).ce());
1143 if(selector) mu_calovec +=
m_assocAcc.calVec(*assoc,iKey);
1144 }
1146 mu_calovec *= std::max<float>(0.,1-(total_eloss/mu_calovec.
ce()));
1147 opx += mu_calovec.
cpx();
1148 opy += mu_calovec.
cpy();
1149 osumpt += mu_calovec.
sumpt();
1150 }
1152
1153 ATH_MSG_VERBOSE(
"Misc cluster px, py, sumpt = " << opx <<
", " << opy <<
", " << osumpt );
1158 }
1159
1160 return StatusCode::SUCCESS;
1161 }
#define ATH_MSG_VERBOSE(x)
SG::ConstAccessor< T, ALLOC > ConstAccessor
size_t index() const
Return the index of this element within its container.
bool acceptTrack(const xAOD::TrackParticle *trk, const xAOD::Vertex *vx) const
columnar::MetHelpers::ObjectWeightDecorator< columnar::ContainerId::mutableMet, columnar::ContainerId::jet > m_jetOutputMetWeightDecSoft
columnar::JetAccessor< float > m_acc_width
columnar::JetAccessor< std::vector< int > > m_acc_trkN
columnar::MetHelpers::MetMomentumAccessors< columnar::ContainerId::mutableMet > m_outputMetMomAcc
columnar::MutableMetAccessor< std::string > m_outputMetNameAcc
columnar::JetAccessor< float > m_acc_emf
columnar::JetAccessor< std::vector< float > > m_acc_trksumpt
columnar::Met1Accessor< MissingETBase::Types::bitmask_t > m_inputMetSourceAcc
columnar::MetHelpers::MetMomentumAccessors< columnar::ContainerId::met1 > m_inputMetMomAcc
columnar::JetAccessor< std::vector< float > > m_acc_sampleE
columnar::JetAccessor< float > m_acc_psf
const xAOD::Vertex * getPV() const
columnar::ElectronAccessor< columnar::RetypeColumn< double, float > > m_electronPtAcc
columnar::MetHelpers::InputMomentumAccessors< columnar::ContainerId::jet > m_jetMomAcc
columnar::MetHelpers::ObjectWeightDecorator< columnar::ContainerId::mutableMet, columnar::ContainerId::jet > m_jetOutputMetWeightDecRegular
float ce() const
Returns .
float sumpt() const
Returns sum of component pt.
float cpt() const
Returns .
float cpx() const
Returns .
float cpy() const
Returns .
xAOD::MissingETAssociation_v1::ConstVec constvec_t
Type for constituent vector.
uint64_t bitmask_t
Type for status word bit mask.
OriginalObjectHandle(const asg::AsgTool &, ObjectRange< CI, CM >) -> OriginalObjectHandle< CI, CM >
OptObjectId< ContainerId::mutableMet > OptMutableMetId
static const SG::AuxElement::ConstAccessor< iplink_t > acc_nominalObject("nominalObjectLink")
static const SG::AuxElement::ConstAccessor< iplink_t > acc_originalObject("originalObjectLink")
static const SG::AuxElement::ConstAccessor< float > acc_Eloss("EnergyLoss")
static const SG::AuxElement::ConstAccessor< std::vector< iplink_t > > acc_ghostMuons("GhostMuon")
@ Photon
The object is a photon.
@ Muon
The object is a muon.
Jet_v1 Jet
Definition of the current "jet version".
TrackParticle_v1 TrackParticle
Reference the current persistent version:
Muon_v1 Muon
Reference the current persistent version:
setBGCode setTAP setLVL2ErrorBits bool
Electron_v1 Electron
Definition of the current "egamma version".
ROOT::Math::LorentzVector< ROOT::Math::PtEtaPhiM4D< double > > JetFourMom_t
Base 4 Momentum type for Jet.
@ Central
Indicator for MET contribution from the central region.
static constexpr bool isXAOD
Whether this is the xAOD mode.