70{
71
73 if (!
jets.isValid()) {
75 return StatusCode::FAILURE;
76 }
77 unsigned int nJets(
jets->size());
78 std::vector<const xAOD::Jet*> jetToCheck; jetToCheck.clear();
79
80
82 std::vector<int>
entries = m_parser->evaluateAsVector();
84
85 if (nJets != nEntries ) {
86 ATH_MSG_ERROR(
"Sizes incompatible! Are you sure your selection string used jets??");
87 return StatusCode::FAILURE;
88 } else {
89
90 for (
unsigned int i=0;
i<nJets; ++
i)
if (
entries[i]==1) jetToCheck.push_back((*jets)[i]);
91 }
92 } else {
93 for (
unsigned int i=0;
i<nJets; ++
i) jetToCheck.push_back((*jets)[i]);
94 }
95
96 using WDH = SG::WriteDecorHandle<xAOD::JetContainer, float>;
119
120 std::vector<const xAOD::IParticle*> mergedGhostConstits;
121 std::vector<const xAOD::IParticle*> ghosts;
122 std::vector<fastjet::PseudoJet> constituents;
123
124 for (
const xAOD::Jet* jet : jetToCheck) {
125
126
127 mergedGhostConstits.clear();
129 ghosts.clear();
130 ATH_CHECK(jet->getAssociatedObjects(ghostName, ghosts));
131 mergedGhostConstits.insert(mergedGhostConstits.end(), ghosts.begin(), ghosts.end());
132 }
133
134
135 constituents.clear();
138 for (auto constit : mergedGhostConstits){
139 if (constit==nullptr) {
141 continue;
142 }
143
144 if (constit->e()<0) {
ATH_MSG_INFO(
"################## Negative energy cluster");
continue;}
145 constituents.push_back( fastjet::PseudoJet(constit->p4()) );
146
147
150 double eConstit = caloConstit->
e()* caloConstit->e();
151 time += caloConstit->time()* eConstit;
153 }
154 }
155
156
157 if (eTot==0) {
158 dec_timing(*jet) = 0;
159 } else {
161 }
162
163
164 int nConstits = constituents.size();
165
166
167 fastjet::PseudoJet groomed_jet;
168 if (nConstits!=0) {
169 auto jet_def = fastjet::JetDefinition(fastjet::antikt_algorithm, 1.5);
170 fastjet::ClusterSequence cs(constituents, jet_def);
171 fastjet::PseudoJet recluster_jet = cs.inclusive_jets(0.0).front();
172
173
175 groomed_jet =
m_trimmer->result(recluster_jet);
178 } else {
179 ATH_MSG_DEBUG(
" No grooming requested or wrong one, will not apply grooming");
180 groomed_jet = recluster_jet;
181 }
182
183
184 nConstits = groomed_jet.constituents().size();
185 }
186
187
188 if (nConstits==0) {
189 dec_Qw(*jet) = -999;
190
191 dec_Tau1(*jet) = -999;
192 dec_Tau2(*jet) = -999;
193 dec_Tau3(*jet) = -999;
194 dec_Tau4(*jet) = -999;
195
196 dec_Tau21(*jet) = -999;
197 dec_Tau32(*jet) = -999;
198
199 dec_Split12(*jet) = -999;
200 dec_Split23(*jet) = -999;
201 dec_Split34(*jet) = -999;
202
203 dec_ECF1(*jet) = -999;
204 dec_ECF2(*jet) = -999;
205 dec_ECF3(*jet) = -999;
206 dec_ECF4(*jet) = -999;
207
208 dec_C2(*jet) = -999;
209 dec_D2(*jet) = -999;
210
211 dec_pT(*jet) = 0;
212 dec_m(*jet) = 0;
213 dec_NConstits(*jet) = 0;
214 dec_eta(*jet) = -999;
215 dec_phi(*jet) = -999;
216
217 return StatusCode::SUCCESS;
218 }
219
220
221 dec_pT(*jet) = groomed_jet.pt();
222 dec_m(*jet) = groomed_jet.m();
223 dec_NConstits(*jet) = nConstits;
224 dec_eta(*jet) = groomed_jet.eta();
225 dec_phi(*jet) = groomed_jet.phi() -
M_PI;
226
227
228 JetSubStructureUtils::Qw
qw;
229 dec_Qw(*jet) =
qw.result(groomed_jet);
230
231
232 fastjet::contrib::WTA_KT_Axes wta_kt_axes;
233 fastjet::contrib::NormalizedCutoffMeasure normalized_measure(1.0, 1.0, 1000000);
234 JetSubStructureUtils::Nsubjettiness Tau1_wta(1, wta_kt_axes, normalized_measure);
235 JetSubStructureUtils::Nsubjettiness Tau2_wta(2, wta_kt_axes, normalized_measure);
236 JetSubStructureUtils::Nsubjettiness Tau3_wta(3, wta_kt_axes, normalized_measure);
237 JetSubStructureUtils::Nsubjettiness Tau4_wta(4, wta_kt_axes, normalized_measure);
238 float tau1_wta = Tau1_wta.result(groomed_jet);
239 float tau2_wta = Tau2_wta.result(groomed_jet);
240 float tau3_wta = Tau3_wta.result(groomed_jet);
241 float tau4_wta = Tau4_wta.result(groomed_jet);
242 float tau21_wta = -999;
243 float tau32_wta = -999;
244
245 if( tau1_wta > 1e-8 ) {
246 tau21_wta = tau2_wta / tau1_wta;
247 }
248 if( tau2_wta > 1e-8 ) {
249 tau32_wta = tau3_wta / tau2_wta;
250 }
251
252 dec_Tau1(*jet) = tau1_wta;
253 dec_Tau2(*jet) = tau2_wta;
254 dec_Tau3(*jet) = tau3_wta;
255 dec_Tau4(*jet) = tau4_wta;
256 dec_Tau21(*jet) = tau21_wta;
257 dec_Tau32(*jet) = tau32_wta;
258
259
260 JetSubStructureUtils::KtSplittingScale Split12(1);
261 JetSubStructureUtils::KtSplittingScale Split23(2);
262 JetSubStructureUtils::KtSplittingScale Split34(3);
263 float split12 = Split12.result(groomed_jet);
264 float split23 = Split23.result(groomed_jet);
265 float split34 = Split34.result(groomed_jet);
266
267 dec_Split12(*jet) = split12;
268 dec_Split23(*jet) = split23;
269 dec_Split34(*jet) = split34;
270
271
272 JetSubStructureUtils::EnergyCorrelator ECF1(1, 1.0, JetSubStructureUtils::EnergyCorrelator::pt_R);
273 JetSubStructureUtils::EnergyCorrelator ECF2(2, 1.0, JetSubStructureUtils::EnergyCorrelator::pt_R);
274 JetSubStructureUtils::EnergyCorrelator ECF3(3, 1.0, JetSubStructureUtils::EnergyCorrelator::pt_R);
275 JetSubStructureUtils::EnergyCorrelator ECF4(4, 1.0, JetSubStructureUtils::EnergyCorrelator::pt_R);
276 float ecf1 = ECF1.result(groomed_jet);
277 float ecf2 = ECF2.result(groomed_jet);
278 float ecf3 = ECF3.result(groomed_jet);
279 float ecf4 = ECF4.result(groomed_jet);
282
283 if( ecf2 > 1e-8 ) {
284 c2 = ecf3 * ecf1 /
pow( ecf2, 2.0 );
285 d2 = ecf3 *
pow( ecf1, 3.0 ) /
pow( ecf2, 3.0 );
286 }
287
288 dec_ECF1(*jet) = ecf1;
289 dec_ECF2(*jet) = ecf2;
290 dec_ECF3(*jet) = ecf3;
291 dec_ECF4(*jet) = ecf4;
294
295 }
296
297 return StatusCode::SUCCESS;
298}
#define ATH_CHECK
Evaluate an expression and check for errors.
constexpr int pow(int base, int exp) noexcept
WDHK m_dec_ECF1
Energy correlation factors decorators.
SG::ReadHandleKey< xAOD::JetContainer > m_jetKey
StringProperty m_selectionString
StringProperty m_grooming
WDHK m_dec_pT
Reclustered jets information decorators.
WDHK m_dec_Split12
KtSplittingScale decorators.
WDHK m_dec_Tau1
Nsubjetiness decorators.
std::optional< fastjet::contrib::SoftDrop > m_softdropper
WDHK m_dec_Qw
Qw decorator.
std::optional< fastjet::Filter > m_trimmer
WDHK m_dec_timing
Timing information.
Gaudi::Property< std::vector< std::string > > m_ghostNames
virtual double e() const
The total energy of the particle.
time(flags, cells_name, *args, **kw)
float eTot(const U &p)
Accessor utility function for getting the value of Energy.
@ CaloCluster
The object is a calorimeter cluster.
Jet_v1 Jet
Definition of the current "jet version".
CaloCluster_v1 CaloCluster
Define the latest version of the calorimeter cluster class.