96{
97
98
99
100
101 SG::ReadHandle<xAOD::TruthParticleContainer>
102 xTruthParticleContainer(
104 if (!xTruthParticleContainer.isValid()) {
105 ATH_MSG_ERROR(
"Could not retrieve xAOD::TruthParticleGenContainer with key:"
107
108 return StatusCode::FAILURE;
109 }
110
112 double primx = 0.;
113 double primy = 0.;
114 int NumMuons = 0;
115 int NumDstars = 0;
116
117 std::vector<const xAOD::TruthParticle *> Muons;
118
119 unsigned int nPart = xTruthParticleContainer->size();
120 ATH_MSG_DEBUG(
"xAODMuDstarFilter:number of particles " << nPart);
121
122
124
129
130 ATH_MSG_DEBUG(
"xAODMuDstarFilter: PV x, y = " << primx <<
" , " << primy);
131 }
132
134 continue;
135
136
137 if (pitr->isMuon())
138 {
142 {
143 NumMuons++;
144 Muons.push_back(pitr);
145 }
146 }
147 }
148
149 if (NumMuons == 0){
151 return StatusCode::SUCCESS;
152 }
153
154 ATH_MSG_DEBUG(
"xAODMuDstarFilter: NumMuons = " << NumMuons );
155
158 continue;
159
160
161 if (std::abs(pitr->pdgId()) ==
MC::DSTAR)
162 {
166 {
167
168
169 if (!(pitr->decayVtx()))
170 continue;
171 double Rxy = std::hypot(pitr->decayVtx()->x() - primx, pitr->decayVtx()->y() - primy);
172
174
176 continue;
177
178 auto firstChild = pitr->decayVtx()->outgoingParticle(0);
179 if (firstChild->pdgId() == pitr->pdgId())
180 continue;
181
182 TLorentzVector p4_pis;
183
184 int pis_pdg = 0;
185 int K_pdg = 0;
186
187 int NumD0 = 0;
188 int NumPis = 0;
189
190 for (size_t thisChild_id = 0; thisChild_id < pitr->decayVtx()->nOutgoingParticles(); thisChild_id++)
191 {
192
193 auto thisChild = pitr->decayVtx()->outgoingParticle(thisChild_id);
195 continue;
196
197 if (std::abs(thisChild->pdgId()) ==
MC::PIPLUS)
198 {
202 {
203 NumPis++;
204
205 p4_pis.SetPtEtaPhiM(thisChild->pt(), thisChild->eta(), thisChild->phi(),
m_PionMass);
206 pis_pdg = thisChild->pdgId();
207 }
208 }
209 }
210
211 if (NumPis == 0)
212 continue;
213
215
217
218 int NumChildD0 = 0;
219 int NumChildD0Charged = 0;
220 int NumChildD0neutrinos = 0;
221 int NumChildD0gammas = 0;
222 int ChargeD0Child1 = 0;
223 int ChargeD0Child2 = 0;
224 int NumChildD0K = 0;
225 int NumChildD0pi = 0;
226 int NumChildD0mu = 0;
227
228 for (size_t thisChild_id = 0; thisChild_id < pitr->decayVtx()->nOutgoingParticles(); thisChild_id++)
229 {
230 auto thisChild = pitr->decayVtx()->outgoingParticle(thisChild_id);
232 continue;
233
234 if (std::abs(thisChild->pdgId()) ==
MC::D0)
235 {
236 if (! thisChild->decayVtx()) continue;
237
238 for (size_t thisChild1_id = 0; thisChild1_id < thisChild->decayVtx()->nOutgoingParticles(); thisChild1_id++)
239 {
240 auto thisChild1 = thisChild->decayVtx()->outgoingParticle(thisChild1_id);
242 continue;
243
244 if (thisChild1->isElectron() || thisChild1->isMuon() ||
246 {
247
248 NumChildD0++;
252 {
253 NumChildD0Charged++;
254
255 if (NumChildD0Charged == 1)
256 {
257 D0Child1 = thisChild1;
258 ChargeD0Child1 = D0Child1->
charge();
259 }
260 if (NumChildD0Charged == 2)
261 {
262 D0Child2 = thisChild1;
263 ChargeD0Child2 = D0Child2->charge();
264 }
265 if (thisChild1->isMuon())
266 {
267 NumChildD0mu++;
268 D0ChildMu = thisChild1;
269 }
270 if (std::abs(thisChild1->pdgId()) ==
MC::PIPLUS)
271 NumChildD0pi++;
272 if (std::abs(thisChild1->pdgId()) ==
MC::KPLUS)
273 {
274 NumChildD0K++;
275 K_pdg = thisChild1->pdgId();
276 }
277 }
278 }
279 else if (std::abs(thisChild1->pdgId()) ==
MC::PI0)
280 {
281 NumChildD0++;
282 }
283 else if (thisChild1->isNeutrino())
284 {
285 NumChildD0neutrinos++;
286 }
287 else if (thisChild1->isPhoton())
288 {
289 NumChildD0gammas++;
290 }
291 else if (std::abs(thisChild1->pdgId()) ==
MC::K0 || std::abs(thisChild1->pdgId()) ==
MC::K0L ||
292 std::abs(thisChild1->pdgId()) ==
MC::K0S)
293 {
294 NumChildD0++;
295 NumChildD0++;
296 }
297 else if (thisChild1->decayVtx())
298 {
299 for (size_t thisChild2_id = 0; thisChild2_id < thisChild1->decayVtx()->nOutgoingParticles(); thisChild2_id++)
300 {
301
302 auto thisChild2 = thisChild1->decayVtx()->outgoingParticle(thisChild2_id);
304 continue;
305
306 if (thisChild2->isElectron() || thisChild2->isMuon() ||
308 {
309 NumChildD0++;
310
314 {
315 NumChildD0Charged++;
316
317 if (NumChildD0Charged == 1)
318 {
319 D0Child1 = thisChild2;
320 ChargeD0Child2 = D0Child1->
charge();
321 }
322 if (NumChildD0Charged == 2)
323 {
324 D0Child2 = thisChild2;
325 ChargeD0Child2 = D0Child2->charge();
326 }
327 if (thisChild2->isMuon())
328 {
329 NumChildD0mu++;
330 D0ChildMu = thisChild2;
331 }
332 }
333 }
334 else if (std::abs(thisChild2->pdgId()) ==
MC::PI0)
335 {
336 NumChildD0++;
337 }
338 else if (thisChild2->isNeutrino())
339 {
340 NumChildD0neutrinos++;
341 }
342 else if (thisChild2->isPhoton())
343 {
344 NumChildD0gammas++;
345 }
346 else if (std::abs(thisChild2->pdgId()) ==
MC::K0 || std::abs(thisChild2->pdgId()) ==
MC::K0L ||
347 std::abs(thisChild2->pdgId()) ==
MC::K0S)
348 {
349 NumChildD0++;
350 NumChildD0++;
351 }
352 else if (thisChild2->decayVtx())
353 {
354 NumChildD0++;
355 NumChildD0++;
356 }
357 else
358 {
359 NumChildD0++;
360 NumChildD0++;
361 ATH_MSG_DEBUG(
"xAODMuDstarFilter: unexpected D0 granddaughter = " << thisChild2->pdgId());
362 }
363 }
364 }
365 else
366 {
367 NumChildD0++;
368 NumChildD0++;
369 ATH_MSG_DEBUG(
"xAODMuDstarFilter: unexpected D0 daughter = " << thisChild1->pdgId());
370 }
371 }
372
373 ATH_MSG_DEBUG(
"xAODMuDstarFilter: NumChildD0, NumChildD0Charged = " << NumChildD0 <<
" , " << NumChildD0Charged);
374
375 if (NumChildD0 <= 3 && NumChildD0Charged == 2 && ChargeD0Child1 * ChargeD0Child2 < 0)
376 {
378 {
379 if (NumChildD0 == 2 && NumChildD0K == 1 && NumChildD0pi == 1){
380 NumD0++;
381 }
382 }
383 else
384 {
385 NumD0++;
386 }
387 }
388 }
389
391
392 if (NumD0 == 1)
393 {
394
395 if (pis_pdg * ChargeD0Child1 > 0)
std::swap(D0Child1, D0Child2);
396 TLorentzVector p4_K;
398 TLorentzVector p4_pi;
399 p4_pi.SetPtEtaPhiM(D0Child2->pt(), D0Child2->eta(), D0Child2->phi(),
m_PionMass);
400
401 TLorentzVector p4_D0 = p4_K + p4_pi;
402 double mKpi = p4_D0.M();
403
405
407 {
408
409 TLorentzVector p4_Dstar = p4_D0 + p4_pis;
410
411 double delta_m = p4_Dstar.M() - mKpi;
412
414
416 {
417 NumDstars++;
418
419 ATH_MSG_DEBUG(
"xAODMuDstarFilter: NumDstars = " << NumDstars);
420
421 for (
size_t i = 0;
i < Muons.size(); ++
i)
422 {
423
424 if (NumChildD0mu == 1)
425 {
426 if (std::fabs(Muons[i]->
pt() - D0ChildMu->pt()) < std::numeric_limits<double>::epsilon())
427 continue;
428 ATH_MSG_DEBUG(
"xAODMuDstarFilter: Mu(pT), D0Mu(pT) = " << Muons[i]->
pt() <<
" , " << D0ChildMu->pt());
429 }
430
431 TLorentzVector p4_Mu;
433
434 TLorentzVector p4_DstarMu = p4_Dstar + p4_Mu;
435
436 ATH_MSG_DEBUG(
"xAODMuDstarFilter: p4_DstarMu.M() = " << p4_DstarMu.M());
437
439 {
440
442 ATH_MSG_DEBUG(
"xAODMuDstarFilter: p4_DstarMu.M() = " << p4_DstarMu.M());
443 ATH_MSG_DEBUG(
"xAODMuDstarFilter: NumChildD0, NumChildD0Charged = " << NumChildD0 <<
" , " << NumChildD0Charged);
444 ATH_MSG_DEBUG(
"xAODMuDstarFilter: NumChildD0K, NumChildD0pi, NumChildD0mu = " << NumChildD0K <<
" , " << NumChildD0pi <<
" , " << NumChildD0mu);
445
446 if (NumChildD0mu == 1)
447 {
448
449 ATH_MSG_DEBUG(
"xAODMuDstarFilter: Mu(pT), D0Mu(pT) = " << Muons[i]->
pt() <<
" , " << D0ChildMu->pt());
450 }
451
452 ATH_MSG_DEBUG(
"xAODMuDstarFilter: NumChildD0neutrinos, NumChildD0gammas = " << NumChildD0neutrinos <<
" , " << NumChildD0gammas);
453 ATH_MSG_DEBUG(
"xAODMuDstarFilter: pis_pdg, K_pdg, ChargeD0Child1, ChargeD0Child2 = " << pis_pdg <<
" , " << K_pdg <<
" , " << ChargeD0Child1 <<
" , " << ChargeD0Child2);
454
456 return StatusCode::SUCCESS;
457 }
458 }
459
460 }
461 }
462 }
463 }
464 }
465
466 }
467
468 }
469
470
471
472
473
474
476 return StatusCode::SUCCESS;
477}
Scalar eta() const
pseudorapidity method
Scalar phi() const
phi method
virtual void setFilterPassed(bool state, const EventContext &ctx) const
Gaudi::Property< double > m_EtaRangeKpi
Gaudi::Property< double > m_EtaRangeDstar
Gaudi::Property< double > m_delta_m_Max
Gaudi::Property< double > m_EtaRangePis
Gaudi::Property< double > m_PtMinKpi
Gaudi::Property< double > m_DstarMu_m_Max
Gaudi::Property< double > m_PtMinMuon
Gaudi::Property< double > m_mKpiMin
Gaudi::Property< double > m_PtMaxMuon
Gaudi::Property< bool > m_D0Kpi_only
SG::ReadHandleKey< xAOD::TruthParticleContainer > m_xaodTruthParticleContainerNameGenKey
Gaudi::Property< double > m_RxyMinDstar
Gaudi::Property< double > m_PtMaxPis
Gaudi::Property< double > m_EtaRangeMuon
Gaudi::Property< double > m_mKpiMax
Gaudi::Property< double > m_PtMinDstar
Gaudi::Property< double > m_PtMaxKpi
Gaudi::Property< double > m_PtMaxDstar
Gaudi::Property< double > m_PtMinPis
virtual double pt() const override final
The transverse momentum ( ) of the particle.
double charge() const
Physical charge.
virtual double eta() const override final
The pseudorapidity ( ) of the particle.
virtual double phi() const override final
The azimuthal angle ( ) of the particle.
float y() const
Vertex y displacement.
float x() const
Vertex x displacement.
bool isBeam(const T &p)
Identify if the particle is beam particle.
bool isPhysical(const T &p)
Identify if the particle is physical, i.e. is stable or decayed.
void swap(ElementLinkVector< DOBJ > &lhs, ElementLinkVector< DOBJ > &rhs)
TruthVertex_v1 TruthVertex
Typedef to implementation.
TruthParticle_v1 TruthParticle
Typedef to implementation.