96{
97
98
99
100
101 const EventContext& context = Gaudi::Hive::currentContext();
102 SG::ReadHandle<xAOD::TruthParticleContainer>
103 xTruthParticleContainer(
105 if (!xTruthParticleContainer.isValid()) {
106 ATH_MSG_ERROR(
"Could not retrieve xAOD::TruthParticleGenContainer with key:"
108
109 return StatusCode::FAILURE;
110 }
111
113 double primx = 0.;
114 double primy = 0.;
115 int NumMuons = 0;
116 int NumDstars = 0;
117
118 std::vector<const xAOD::TruthParticle *> Muons;
119
120 unsigned int nPart = xTruthParticleContainer->size();
121 ATH_MSG_DEBUG(
"xAODMuDstarFilter:number of particles " << nPart);
122
123
125
130
131 ATH_MSG_DEBUG(
"xAODMuDstarFilter: PV x, y = " << primx <<
" , " << primy);
132 }
133
135 continue;
136
137
138 if (pitr->isMuon())
139 {
143 {
144 NumMuons++;
145 Muons.push_back(pitr);
146 }
147 }
148 }
149
150 if (NumMuons == 0){
151 setFilterPassed(false);
152 return StatusCode::SUCCESS;
153 }
154
155 ATH_MSG_DEBUG(
"xAODMuDstarFilter: NumMuons = " << NumMuons );
156
159 continue;
160
161
162 if (std::abs(pitr->pdgId()) ==
MC::DSTAR)
163 {
167 {
168
169
170 if (!(pitr->decayVtx()))
171 continue;
172 double Rxy = std::hypot(pitr->decayVtx()->x() - primx, pitr->decayVtx()->y() - primy);
173
175
177 continue;
178
179 auto firstChild = pitr->decayVtx()->outgoingParticle(0);
180 if (firstChild->pdgId() == pitr->pdgId())
181 continue;
182
183 TLorentzVector p4_pis;
184
185 int pis_pdg = 0;
186 int K_pdg = 0;
187
188 int NumD0 = 0;
189 int NumPis = 0;
190
191 for (size_t thisChild_id = 0; thisChild_id < pitr->decayVtx()->nOutgoingParticles(); thisChild_id++)
192 {
193
194 auto thisChild = pitr->decayVtx()->outgoingParticle(thisChild_id);
196 continue;
197
198 if (std::abs(thisChild->pdgId()) ==
MC::PIPLUS)
199 {
203 {
204 NumPis++;
205
206 p4_pis.SetPtEtaPhiM(thisChild->pt(), thisChild->eta(), thisChild->phi(),
m_PionMass);
207 pis_pdg = thisChild->pdgId();
208 }
209 }
210 }
211
212 if (NumPis == 0)
213 continue;
214
216
218
219 int NumChildD0 = 0;
220 int NumChildD0Charged = 0;
221 int NumChildD0neutrinos = 0;
222 int NumChildD0gammas = 0;
223 int ChargeD0Child1 = 0;
224 int ChargeD0Child2 = 0;
225 int NumChildD0K = 0;
226 int NumChildD0pi = 0;
227 int NumChildD0mu = 0;
228
229 for (size_t thisChild_id = 0; thisChild_id < pitr->decayVtx()->nOutgoingParticles(); thisChild_id++)
230 {
231 auto thisChild = pitr->decayVtx()->outgoingParticle(thisChild_id);
233 continue;
234
235 if (std::abs(thisChild->pdgId()) ==
MC::D0)
236 {
237 if (! thisChild->decayVtx()) continue;
238
239 for (size_t thisChild1_id = 0; thisChild1_id < thisChild->decayVtx()->nOutgoingParticles(); thisChild1_id++)
240 {
241 auto thisChild1 = thisChild->decayVtx()->outgoingParticle(thisChild1_id);
243 continue;
244
245 if (thisChild1->isElectron() || thisChild1->isMuon() ||
247 {
248
249 NumChildD0++;
253 {
254 NumChildD0Charged++;
255
256 if (NumChildD0Charged == 1)
257 {
258 D0Child1 = thisChild1;
259 ChargeD0Child1 = D0Child1->
charge();
260 }
261 if (NumChildD0Charged == 2)
262 {
263 D0Child2 = thisChild1;
264 ChargeD0Child2 = D0Child2->charge();
265 }
266 if (thisChild1->isMuon())
267 {
268 NumChildD0mu++;
269 D0ChildMu = thisChild1;
270 }
271 if (std::abs(thisChild1->pdgId()) ==
MC::PIPLUS)
272 NumChildD0pi++;
273 if (std::abs(thisChild1->pdgId()) ==
MC::KPLUS)
274 {
275 NumChildD0K++;
276 K_pdg = thisChild1->pdgId();
277 }
278 }
279 }
280 else if (std::abs(thisChild1->pdgId()) ==
MC::PI0)
281 {
282 NumChildD0++;
283 }
284 else if (thisChild1->isNeutrino())
285 {
286 NumChildD0neutrinos++;
287 }
288 else if (thisChild1->isPhoton())
289 {
290 NumChildD0gammas++;
291 }
292 else if (std::abs(thisChild1->pdgId()) ==
MC::K0 || std::abs(thisChild1->pdgId()) ==
MC::K0L ||
293 std::abs(thisChild1->pdgId()) ==
MC::K0S)
294 {
295 NumChildD0++;
296 NumChildD0++;
297 }
298 else if (thisChild1->decayVtx())
299 {
300 for (size_t thisChild2_id = 0; thisChild2_id < thisChild1->decayVtx()->nOutgoingParticles(); thisChild2_id++)
301 {
302
303 auto thisChild2 = thisChild1->decayVtx()->outgoingParticle(thisChild2_id);
305 continue;
306
307 if (thisChild2->isElectron() || thisChild2->isMuon() ||
309 {
310 NumChildD0++;
311
315 {
316 NumChildD0Charged++;
317
318 if (NumChildD0Charged == 1)
319 {
320 D0Child1 = thisChild2;
321 ChargeD0Child2 = D0Child1->
charge();
322 }
323 if (NumChildD0Charged == 2)
324 {
325 D0Child2 = thisChild2;
326 ChargeD0Child2 = D0Child2->charge();
327 }
328 if (thisChild2->isMuon())
329 {
330 NumChildD0mu++;
331 D0ChildMu = thisChild2;
332 }
333 }
334 }
335 else if (std::abs(thisChild2->pdgId()) ==
MC::PI0)
336 {
337 NumChildD0++;
338 }
339 else if (thisChild2->isNeutrino())
340 {
341 NumChildD0neutrinos++;
342 }
343 else if (thisChild2->isPhoton())
344 {
345 NumChildD0gammas++;
346 }
347 else if (std::abs(thisChild2->pdgId()) ==
MC::K0 || std::abs(thisChild2->pdgId()) ==
MC::K0L ||
348 std::abs(thisChild2->pdgId()) ==
MC::K0S)
349 {
350 NumChildD0++;
351 NumChildD0++;
352 }
353 else if (thisChild2->decayVtx())
354 {
355 NumChildD0++;
356 NumChildD0++;
357 }
358 else
359 {
360 NumChildD0++;
361 NumChildD0++;
362 ATH_MSG_DEBUG(
"xAODMuDstarFilter: unexpected D0 granddaughter = " << thisChild2->pdgId());
363 }
364 }
365 }
366 else
367 {
368 NumChildD0++;
369 NumChildD0++;
370 ATH_MSG_DEBUG(
"xAODMuDstarFilter: unexpected D0 daughter = " << thisChild1->pdgId());
371 }
372 }
373
374 ATH_MSG_DEBUG(
"xAODMuDstarFilter: NumChildD0, NumChildD0Charged = " << NumChildD0 <<
" , " << NumChildD0Charged);
375
376 if (NumChildD0 <= 3 && NumChildD0Charged == 2 && ChargeD0Child1 * ChargeD0Child2 < 0)
377 {
379 {
380 if (NumChildD0 == 2 && NumChildD0K == 1 && NumChildD0pi == 1){
381 NumD0++;
382 }
383 }
384 else
385 {
386 NumD0++;
387 }
388 }
389 }
390
392
393 if (NumD0 == 1)
394 {
395
396 if (pis_pdg * ChargeD0Child1 > 0)
std::swap(D0Child1, D0Child2);
397 TLorentzVector p4_K;
399 TLorentzVector p4_pi;
400 p4_pi.SetPtEtaPhiM(D0Child2->pt(), D0Child2->eta(), D0Child2->phi(),
m_PionMass);
401
402 TLorentzVector p4_D0 = p4_K + p4_pi;
403 double mKpi = p4_D0.M();
404
406
408 {
409
410 TLorentzVector p4_Dstar = p4_D0 + p4_pis;
411
412 double delta_m = p4_Dstar.M() - mKpi;
413
415
417 {
418 NumDstars++;
419
420 ATH_MSG_DEBUG(
"xAODMuDstarFilter: NumDstars = " << NumDstars);
421
422 for (
size_t i = 0;
i < Muons.size(); ++
i)
423 {
424
425 if (NumChildD0mu == 1)
426 {
427 if (std::fabs(Muons[i]->
pt() - D0ChildMu->pt()) < std::numeric_limits<double>::epsilon())
428 continue;
429 ATH_MSG_DEBUG(
"xAODMuDstarFilter: Mu(pT), D0Mu(pT) = " << Muons[i]->
pt() <<
" , " << D0ChildMu->pt());
430 }
431
432 TLorentzVector p4_Mu;
434
435 TLorentzVector p4_DstarMu = p4_Dstar + p4_Mu;
436
437 ATH_MSG_DEBUG(
"xAODMuDstarFilter: p4_DstarMu.M() = " << p4_DstarMu.M());
438
440 {
441
443 ATH_MSG_DEBUG(
"xAODMuDstarFilter: p4_DstarMu.M() = " << p4_DstarMu.M());
444 ATH_MSG_DEBUG(
"xAODMuDstarFilter: NumChildD0, NumChildD0Charged = " << NumChildD0 <<
" , " << NumChildD0Charged);
445 ATH_MSG_DEBUG(
"xAODMuDstarFilter: NumChildD0K, NumChildD0pi, NumChildD0mu = " << NumChildD0K <<
" , " << NumChildD0pi <<
" , " << NumChildD0mu);
446
447 if (NumChildD0mu == 1)
448 {
449
450 ATH_MSG_DEBUG(
"xAODMuDstarFilter: Mu(pT), D0Mu(pT) = " << Muons[i]->
pt() <<
" , " << D0ChildMu->pt());
451 }
452
453 ATH_MSG_DEBUG(
"xAODMuDstarFilter: NumChildD0neutrinos, NumChildD0gammas = " << NumChildD0neutrinos <<
" , " << NumChildD0gammas);
454 ATH_MSG_DEBUG(
"xAODMuDstarFilter: pis_pdg, K_pdg, ChargeD0Child1, ChargeD0Child2 = " << pis_pdg <<
" , " << K_pdg <<
" , " << ChargeD0Child1 <<
" , " << ChargeD0Child2);
455
456 setFilterPassed(true);
457 return StatusCode::SUCCESS;
458 }
459 }
460
461 }
462 }
463 }
464 }
465 }
466
467 }
468
469 }
470
471
472
473
474
475
476 setFilterPassed(false);
477 return StatusCode::SUCCESS;
478}
Scalar eta() const
pseudorapidity method
Scalar phi() const
phi method
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.