95 {
96 const EventContext &ctx = Gaudi::Hive::currentContext();
97 SG::ReadHandle<xAOD::VertexContainer> vertices(
m_vertexInKey, ctx);
99
100 SG::ReadHandle<xAOD::PhotonContainer> photonsIn(
m_photonsInKey, ctx);
102
105 SG::ReadHandle<xAOD::MuonContainer> muonsIn(
m_muonsInKey, ctx);
107 SG::ReadHandle<xAOD::JetContainer> jetsIn(
m_jetsInKey, ctx);
109 SG::ReadHandle<xAOD::EventInfo> eventInfo(
m_eventInKey,ctx);
111
112 SG::WriteDecorHandle<xAOD::VertexContainer, std::vector<ElementLink<xAOD::PhotonContainer>>> dec_photonLinks(
m_photonLinksKey, ctx);
113 SG::WriteDecorHandle<xAOD::VertexContainer, std::vector<ElementLink<xAOD::JetContainer>>> dec_jetLinks(
m_jetLinksKey, ctx);
114 SG::WriteDecorHandle<xAOD::VertexContainer, std::vector<ElementLink<xAOD::ElectronContainer>>> dec_electronLinks(
m_electronLinksKey, ctx);
115 SG::WriteDecorHandle<xAOD::VertexContainer, std::vector<ElementLink<xAOD::MuonContainer>>> dec_muonLinks(
m_muonLinksKey, ctx);
116 SG::WriteDecorHandle<xAOD::VertexContainer, std::vector<ElementLink<xAOD::CompositeParticleContainer>>> dec_multiPhotonLinks(
m_multiPhotonLinksKey, ctx);
117
118 SG::ReadDecorHandle<xAOD::VertexContainer, float> acc_deltaZ(
m_deltaZKey, ctx);
119 SG::ReadDecorHandle<xAOD::VertexContainer, float> acc_deltaPhi(
m_deltaPhiKey, ctx);
120 SG::ReadDecorHandle<xAOD::PhotonContainer, float> acc_caloPointingZ(
m_caloPointingZKey, ctx);
121 SG::ReadDecorHandle<xAOD::PhotonContainer, float> acc_zCommon(
m_zCommonKey, ctx);
122 SG::ReadDecorHandle<xAOD::PhotonContainer, float> acc_zCommonError(
m_zCommonErrorKey, ctx);
123
124
125 SG::WriteDecorHandle<xAOD::VertexContainer, int> dec_ntrk(
m_mDecor_ntrk, ctx);
126 SG::WriteDecorHandle<xAOD::VertexContainer, float> dec_sumPt(
m_mDecor_sumPt,ctx);
128 SG::WriteDecorHandle<xAOD::VertexContainer, float> dec_z_asym(
m_mDecor_z_asym,ctx);
131 SG::WriteDecorHandle<xAOD::VertexContainer, float> dec_z_skew(
m_mDecor_z_skew,ctx);
136
137 auto mpCont = std::make_unique<xAOD::CompositeParticleContainer>();
138 auto mpAux = std::make_unique<xAOD::CompositeParticleAuxContainer>();
139 mpCont->setStore(mpAux.get());
140
143 bool haveMP = false;
144
145
146 const bool haveZCommonDecor = acc_zCommon.isAvailable() && acc_zCommonError.isAvailable();
147 const bool haveCaloPointingDecor = acc_caloPointingZ.isAvailable();
148
149 if (photonsIn->size() > 1 && (!haveZCommonDecor || !haveCaloPointingDecor))
150 {
151 ATH_MSG_ERROR(
"photonsIn has size > 1 but required decorations are missing: "
154 return StatusCode::FAILURE;
155 }
156
157 if (photonsIn->size() > 1)
158 {
159
160 double sumW = 0.0;
161 double sumWZ = 0.0;
162
163 TLorentzVector p4sum;
164 int nUsed = 0;
166 {
167 const float zCommon = acc_zCommon(*ph);
168 const float zCaloPointing = acc_caloPointingZ(*ph);
169
170 const bool useZCommonValue = std::isfinite(zCommon);
171 const float z = useZCommonValue ? zCommon : zCaloPointing;
172 if (!std::isfinite(
z)) {
173 ATH_MSG_ERROR(
"Non-finite photon pointing values found in both "
175 return StatusCode::FAILURE;
176 }
177
178
179
181 if (useZCommonValue) {
182 const float s = acc_zCommonError(*ph);
183 if (std::isfinite(s) && s > 0.0f)
w = 1.0 / (
double(s) *
double(s));
184 }
185
188 p4sum += ph->p4();
189 ++nUsed;
190 }
191
192 if (nUsed > 1 && sumW > 0.0)
193 {
194 const float zPoint = sumWZ / sumW;
195
196
197 float bestAbsDZ = 1e30;
199 {
201
202 const float dz = zPoint - vtx->z();
203 const float adz = std::abs(dz);
204 if (adz < bestAbsDZ)
205 {
206 bestAbsDZ = adz;
207 bestVtxForMP = vtx;
208 }
209 }
210
211 if (bestVtxForMP)
212 {
213
215 mpCont->push_back(mp);
216
217 mp->setP4(p4sum);
218
219 SG::AuxElement::Decorator<float> dec_mp_deltaZ("deltaZ");
220 SG::AuxElement::Decorator<float> dec_mp_deltaPhi("deltaPhi");
221 SG::AuxElement::Decorator<int> dec_mp_nPhotons("nPhotons");
222 SG::AuxElement::Decorator<float> dec_mp_zPointing("zPointing");
223
224 const float deltaZ = zPoint - bestVtxForMP->
z();
225 dec_mp_deltaZ(*mp) =
deltaZ;
226 dec_mp_zPointing(*mp) = zPoint;
227 dec_mp_nPhotons(*mp) = static_cast<int>(photonsIn->size());
228
229 float dphi = 0.0f;
230 if (acc_deltaPhi.isAvailable())
231 {
232 dphi = acc_deltaPhi(*bestVtxForMP);
233 if (!std::isfinite(dphi)) dphi = 0.0f;
234 }
235 dec_mp_deltaPhi(*mp) = dphi;
238
239 haveMP = true;
240 }
241 }
242 }
243
244
245 ATH_CHECK(mpHandle.record(std::move(mpCont), std::move(mpAux)));
246
247 std::map< const xAOD::Vertex*, std::vector<ElementLink<xAOD::JetContainer>> > jetsInVertex;
248 std::map< const xAOD::Jet*, std::map< const xAOD::Vertex*, int> > jetVertexPt;
249
250
252 jetsInVertex[
vertex] = {};
254 jetVertexPt[jet][
vertex] = 0;
255 }
256 }
257
258
260
261
264 if( !jtrk ) continue;
266 if(jetTrackVertex) jetVertexPt[jet][*jetTrackVertex] += jtrk->pt();
267 }
268
269
270 float maxPtFrac = -1;
274 if(jetVertexPt[jet][vertex] > maxPtFrac){
275 maxPtFrac = jetVertexPt[jet][
vertex];
276 uniqueVertexAddress =
vertex;
277 }
278 }
279
280
281 ElementLink<xAOD::JetContainer> jetLink;
284 jetsInVertex[uniqueVertexAddress].push_back(jetLink);
285 }
286
288 {
290 continue;
291
292 if (
vertex->nTrackParticles() < 2)
293 continue;
294
295 dec_actualInterPerXing(*vertex) = eventInfo->actualInteractionsPerCrossing();
296
297
298 float z_asym = 0;
299 float sumDZ = 0;
301 float modsumDZ = 0;
302 float weighted_sumDZ = 0;
303 float weighted_deltaZ = 0;
304 float weighted_modsumDZ = 0;
305 float weighted_z_asym = 0;
306
307
308 std::vector<float> track_deltaZ;
309
310 for (
size_t i = 0;
i <
vertex->nTrackParticles();
i++) {
311
313
314 if(!trackTmp) continue;
315
317 track_deltaZ.push_back(deltaZ);
318
319 float trk_weight =
vertex->trackWeight(i);
320 weighted_deltaZ =
deltaZ * trk_weight;
321
323 modsumDZ += std::abs(deltaZ);
324 weighted_sumDZ += weighted_deltaZ;
325 weighted_modsumDZ += std::abs(weighted_deltaZ);
326 }
327
328 if (modsumDZ > 0) {
329 z_asym = sumDZ / modsumDZ;
330 }
331 if (weighted_modsumDZ > 0) {
332 weighted_z_asym = weighted_sumDZ / weighted_modsumDZ;
333 }
334
335 const float number_tracks = track_deltaZ.size();
336 const float mean_Dz =
337 number_tracks > 0 ? sumDZ / number_tracks : 0.F;
338
339 float z_skew = 0;
340 float z_kurt = 0;
341 float z_var = 0;
342
343 for (auto i : track_deltaZ)
344 {
345 float z_zbar = (
i - mean_Dz);
346 z_var += std::pow(z_zbar, 2);
347 z_skew += std::pow(z_zbar, 3);
348 z_kurt += std::pow(z_zbar, 4);
349 }
350 if (number_tracks > 1 && z_var > 0) {
351 z_var /= (number_tracks - 1);
352 float z_sd = std::sqrt(z_var);
353 const float skew_denom = (number_tracks - 1) * std::pow(z_sd, 3);
354 const float kurt_denom = (number_tracks - 1) * std::pow(z_sd, 4);
355 if (std::isfinite(skew_denom) && skew_denom != 0.F) {
356 z_skew /= skew_denom;
357 } else {
358 z_skew = 0.F;
359 }
360 if (std::isfinite(kurt_denom) && kurt_denom != 0.F) {
361 z_kurt /= kurt_denom;
362 } else {
363 z_kurt = 0.F;
364 }
365 }
366 else
367 {
369 z_skew = 0.;
370 z_kurt = 0.;
371 }
372
373 dec_ntrk(*vertex) = number_tracks;
374
375 if(!dec_sumPt.isAvailable()){
377 }
378 const float numberDoF =
vertex->numberDoF();
379 dec_chi2Over_ndf(*vertex) =
380 numberDoF > 0.F ?
vertex->chiSquared() / numberDoF : 0.F;
381 dec_z_asym(*vertex) = z_asym;
382 dec_weighted_z_asym(*vertex) = weighted_z_asym;
383 dec_z_kurt(*vertex) = z_kurt;
384 dec_z_skew(*vertex) = z_skew;
385
386 if (acc_deltaZ.isAvailable() && photonsIn->size() > 0 ) {
387
388 if (!std::isfinite(acc_deltaZ(*vertex))) {
390 dec_photon_deltaz(*vertex) = -1;
391 }
392 else{
393 dec_photon_deltaz(*vertex) = acc_deltaZ(*vertex);
394 }
395 }
396 else{
397 dec_photon_deltaz(*vertex) = -1;
398 }
399 if (acc_deltaPhi.isAvailable() && photonsIn->size() > 0) {
400
401 if (!std::isfinite(acc_deltaPhi(*vertex))) {
403 dec_photon_deltaPhi(*vertex) = -1;
404 }
405 else{
406 dec_photon_deltaPhi(*vertex) = acc_deltaPhi(*vertex);
407 }
408 }
409 else{
410 dec_photon_deltaPhi(*vertex) = -1;
411 }
412
413
414 std::vector<ElementLink<xAOD::ElectronContainer>> electronLinks;
417 if(!id_trk) continue;
419 if(!eleVertex) continue;
420 ElementLink<xAOD::ElectronContainer> elLink;
421 if(*eleVertex == vertex){
424 electronLinks.push_back(elLink);
425 }
426 }
427 dec_electronLinks(*vertex) = electronLinks;
428
429 std::vector<ElementLink<xAOD::PhotonContainer>> photonLinks;
431 ElementLink<xAOD::PhotonContainer> phLink;
434 photonLinks.push_back(phLink);
435 }
436 dec_photonLinks(*vertex) = photonLinks;
437
438
439 std::vector<ElementLink<xAOD::CompositeParticleContainer>> mpLinks;
440 if (haveMP && vertex == bestVtxForMP)
441 {
442 mpLinks.push_back(mpLink);
443 }
444 dec_multiPhotonLinks(*vertex) = mpLinks;
445
446
447 dec_jetLinks(*vertex) = jetsInVertex[
vertex];
448
449 std::vector<ElementLink<xAOD::MuonContainer>> muonLinks;
451 const auto *tp =
muon->trackParticle(xAOD::Muon::InnerDetectorTrackParticle);
452 if(!tp) continue;
453 try{
455 if(!muonVertex) continue;
456 ElementLink<xAOD::MuonContainer> muonLink;
457 if(*muonVertex == vertex){
460 muonLinks.push_back(muonLink);
461 }
462 }catch(...) {
463 ATH_MSG_DEBUG(
"Skipping muon as the track is not associated to any PV ");
465 }
466
467 }
468 dec_muonLinks(*vertex) = muonLinks;
469
470
472 }
473
474 return StatusCode::SUCCESS;
475 }
#define ATH_CHECK
Evaluate an expression and check for errors.
ElementLink()
Default constructor.
#define ATH_MSG_WARNING(x)
bool setElement(ElementType element)
Set link to point to an Element (slowest).
bool setStorableObject(BaseConstReference data, bool replace=false)
Set link storable to data object pointed by data (slower).
SG::WriteDecorHandleKey< xAOD::VertexContainer > m_mDecor_photon_deltaPhi
ToolHandle< GNNTool > m_gnnTool
SG::WriteDecorHandleKey< xAOD::VertexContainer > m_mDecor_z_skew
SG::WriteDecorHandleKey< xAOD::VertexContainer > m_muonLinksKey
SG::WriteDecorHandleKey< xAOD::VertexContainer > m_mDecor_sumPt
SG::WriteHandleKey< xAOD::CompositeParticleContainer > m_multiPhotonsOutKey
SG::WriteDecorHandleKey< xAOD::VertexContainer > m_mDecor_chi2Over_ndf
SG::ReadDecorHandleKey< xAOD::VertexContainer > m_deltaZKey
SG::WriteDecorHandleKey< xAOD::VertexContainer > m_photonLinksKey
SG::ReadHandleKey< xAOD::MuonContainer > m_muonsInKey
SG::WriteDecorHandleKey< xAOD::VertexContainer > m_mDecor_weighted_z_asym
SG::ReadDecorHandleKey< xAOD::PhotonContainer > m_zCommonKey
SG::WriteDecorHandleKey< xAOD::VertexContainer > m_electronLinksKey
SG::ReadHandleKey< xAOD::EventInfo > m_eventInKey
SG::WriteDecorHandleKey< xAOD::VertexContainer > m_mDecor_ntrk
SG::WriteDecorHandleKey< xAOD::VertexContainer > m_mDecor_actualInterPerXing
ToolHandle< CP::TrackVertexAssociationTool > m_trkVtxAssociationTool
SG::ReadDecorHandleKey< xAOD::VertexContainer > m_deltaPhiKey
SG::WriteDecorHandleKey< xAOD::VertexContainer > m_mDecor_weighted_z_kurt
SG::ReadDecorHandleKey< xAOD::PhotonContainer > m_zCommonErrorKey
SG::WriteDecorHandleKey< xAOD::VertexContainer > m_mDecor_photon_deltaz
SG::ReadDecorHandleKey< xAOD::PhotonContainer > m_caloPointingZKey
SG::WriteDecorHandleKey< xAOD::VertexContainer > m_multiPhotonLinksKey
SG::WriteDecorHandleKey< xAOD::VertexContainer > m_mDecor_z_asym
SG::ReadHandleKey< xAOD::VertexContainer > m_vertexInKey
SG::ReadHandleKey< xAOD::PhotonContainer > m_photonsInKey
SG::WriteDecorHandleKey< xAOD::VertexContainer > m_jetLinksKey
SG::ReadHandleKey< xAOD::ElectronContainer > m_electronsInKey
SG::ReadHandleKey< xAOD::JetContainer > m_jetsInKey
float z0() const
Returns the parameter.
float vz() const
The z origin for the parameters.
float z() const
Returns the z position.
const xAOD::TrackParticle * getOriginalTrackParticle(const xAOD::Electron *el)
Helper function for getting the "Original" Track Particle (i.e before GSF) via the electron.
float getVertexSumPt(const xAOD::Vertex *vertex, int power=1, bool useAux=true)
Loop over track particles associated with vertex and return scalar sum of pT^power in GeV (from auxda...
@ NoVtx
Dummy vertex. TrackParticle was not used in vertex fit.
Jet_v1 Jet
Definition of the current "jet version".
TrackParticle_v1 TrackParticle
Reference the current persistent version:
Vertex_v1 Vertex
Define the latest version of the vertex class.
CompositeParticle_v1 CompositeParticle
Define the latest version of the composite particle class.
Muon_v1 Muon
Reference the current persistent version:
Photon_v1 Photon
Definition of the current "egamma version".
Electron_v1 Electron
Definition of the current "egamma version".