126 {
127
128
129
131
133
134 for (itr = events_const() -> begin(); itr != events_const() -> end(); ++itr) {
135
136 double primx = 0.;
137 double primy = 0.;
138
139 #ifdef HEPMC3
141 #else
143 #endif
144
145 primx = vprim -> position().x();
146 primy = vprim -> position().y();
147
148 ATH_MSG_DEBUG(
"MuDstarFilter: PV x, y = " << primx <<
" , " << primy);
149
150 int NumMuons = 0;
151 int NumDstars = 0;
152
153 std::vector < HepMC::ConstGenParticlePtr >
Muons;
154
155
156 const HepMC::GenEvent * genEvt = ( * itr);
157 for (const auto & pitr: * genEvt) {
158
160
161
163 if ((pitr -> momentum().
perp() >= m_PtMinMuon) &&
164 (pitr -> momentum().
perp() < m_PtMaxMuon) &&
165 (std::abs(pitr -> momentum().
eta()) < m_EtaRangeMuon)) {
166 ++NumMuons;
167 Muons.push_back(pitr);
168 }
169 }
170 }
171
173
174 if (NumMuons == 0) break;
175 for (const auto & pitr: * genEvt) {
176
178
179
180 if (std::abs(pitr -> pdg_id()) ==
MC::DSTAR) {
183 (std::abs(pitr ->
momentum().
eta()) < m_EtaRangeDstar)) {
184
185
186 if (!(pitr -> end_vertex())) continue;
187
188 double Rxy = std::sqrt(
pow(pitr -> end_vertex() -> position().
x() - primx, 2) +
pow(pitr -> end_vertex() -> position().
y() - primy, 2));
189
191
192 if (Rxy < m_RxyMinDstar) continue;
193
194
195 #ifdef HEPMC3
196 auto firstChild = pitr -> end_vertex() -> particles_out().begin();
197 auto endChild = pitr -> end_vertex() -> particles_out().end();
198 auto thisChild = firstChild;
199 #else
200 HepMC::GenVertex::particle_iterator firstChild =
201 pitr -> end_vertex() -> particles_begin(HepMC::children);
202 HepMC::GenVertex::particle_iterator endChild =
203 pitr -> end_vertex() -> particles_end(HepMC::children);
204 HepMC::GenVertex::particle_iterator thisChild = firstChild;
205 #endif
206
207 if (( * firstChild) -> pdg_id() == pitr -> pdg_id()) continue;
208
209 TLorentzVector p4_K;
210 TLorentzVector p4_pi;
211 TLorentzVector p4_pis;
212 TLorentzVector p4_D0;
213 TLorentzVector p4_Dstar;
214 TLorentzVector p4_Mu;
215 TLorentzVector p4_DstarMu;
216
217 int pis_pdg = 0;
218 int K_pdg = 0;
219
220 int NumD0 = 0;
221 int NumPis = 0;
222
223 for (; thisChild != endChild; ++thisChild) {
224
226
227 if (std::abs(( * thisChild) -> pdg_id()) ==
MC::PIPLUS) {
228 if ((( * thisChild) ->
momentum().
perp() >= m_PtMinPis) &&
230 (std::abs(( * thisChild) ->
momentum().
eta()) < m_EtaRangePis)) {
231 ++NumPis;
232
234 pis_pdg = ( * thisChild) -> pdg_id();
235 }
236 }
237 }
238
240
241 if (NumPis != 1) continue;
242
246
247 int NumChildD0 = 0;
248 int NumChildD0Charged = 0;
249 int NumChildD0neutrinos = 0;
250 int NumChildD0gammas = 0;
251 int ChargeD0Child1 = 0;
252 int ChargeD0Child2 = 0;
253 int NumChildD0K = 0;
254 int NumChildD0pi = 0;
255 int NumChildD0mu = 0;
256
257 for (thisChild = firstChild; thisChild != endChild; ++thisChild) {
258
260
261 if (std::abs(( * thisChild) -> pdg_id()) ==
MC::D0) {
262 if ((( * thisChild) -> end_vertex())) {
263 #ifdef HEPMC3
264 auto firstChild1 = ( * thisChild) -> end_vertex() -> particles_out().begin();
265 auto endChild1 = ( * thisChild) -> end_vertex() -> particles_out().end();
266 auto thisChild1 = firstChild1;
267 #else
268 HepMC::GenVertex::particle_iterator firstChild1 =
269 ( * thisChild) -> end_vertex() -> particles_begin(HepMC::children);
270 HepMC::GenVertex::particle_iterator endChild1 =
271 ( * thisChild) -> end_vertex() -> particles_end(HepMC::children);
272 HepMC::GenVertex::particle_iterator thisChild1 = firstChild1;
273
274 #endif
275
276 for (; thisChild1 != endChild1; ++thisChild1) {
277
279
280 if (std::abs(( * thisChild1) -> pdg_id()) ==
MC::ELECTRON || std::abs(( * thisChild1) -> pdg_id()) ==
MC::MUON ||
281 std::abs(( * thisChild1) -> pdg_id()) ==
MC::PIPLUS || std::abs(( * thisChild1) -> pdg_id()) ==
MC::KPLUS) {
282
283 ++NumChildD0;
284
285 if ((( * thisChild1) ->
momentum().
perp() >= m_PtMinKpi) &&
287 (std::abs(( * thisChild1) ->
momentum().
eta()) < m_EtaRangeKpi)) {
288 ++NumChildD0Charged;
289
290 if (NumChildD0Charged == 1) {
291 D0Child1 = ( * thisChild1);
294 ChargeD0Child1 = +1;
295 } else {
296 ChargeD0Child1 = -1;
297 }
298 }
299 if (NumChildD0Charged == 2) {
300 D0Child2 = ( * thisChild1);
303 ChargeD0Child2 = +1;
304 } else {
305 ChargeD0Child2 = -1;
306 }
307 }
308 if (std::abs(( * thisChild1) -> pdg_id()) ==
MC::MUON) {
309 ++NumChildD0mu;
310 D0ChildMu = ( * thisChild1);
311 }
312 if (std::abs(( * thisChild1) -> pdg_id()) ==
MC::PIPLUS) ++NumChildD0pi;
313 if (std::abs(( * thisChild1) -> pdg_id()) ==
MC::KPLUS) {
314 ++NumChildD0K;
315 K_pdg = ( * thisChild1) -> pdg_id();
316 }
317 }
318
319 }
else if (std::abs(( * thisChild1) -> pdg_id()) ==
MC::PI0) {
320 ++NumChildD0;
321 }
else if (std::abs(( * thisChild1) -> pdg_id()) ==
MC::NU_E || std::abs(( * thisChild1) -> pdg_id()) ==
MC::NU_MU) {
322 ++NumChildD0neutrinos;
323 }
else if (std::abs(( * thisChild1) -> pdg_id()) ==
MC::PHOTON) {
324 ++NumChildD0gammas;
325 }
else if (std::abs(( * thisChild1) -> pdg_id()) ==
MC::K0 || std::abs(( * thisChild1) -> pdg_id()) ==
MC::K0L ||
326 std::abs(( * thisChild1) -> pdg_id()) ==
MC::K0S) {
327 ++NumChildD0;
328 ++NumChildD0;
329 } else if ((( * thisChild1) -> end_vertex())) {
330
331
332 #ifdef HEPMC3
333 auto firstChild2 = ( * thisChild1) -> end_vertex() -> particles_out().begin();
334 auto endChild2 = ( * thisChild1) -> end_vertex() -> particles_out().end();
335 auto thisChild2 = firstChild2;
336 #else
337 HepMC::GenVertex::particle_iterator firstChild2 =
338 ( * thisChild1) -> end_vertex() -> particles_begin(HepMC::children);
339 HepMC::GenVertex::particle_iterator endChild2 =
340 ( * thisChild1) -> end_vertex() -> particles_end(HepMC::children);
341 HepMC::GenVertex::particle_iterator thisChild2 = firstChild2;
342
343 #endif
344 for (; thisChild2 != endChild2; ++thisChild2) {
345
347
348 if (std::abs(( * thisChild2) -> pdg_id()) ==
MC::ELECTRON || std::abs(( * thisChild2) -> pdg_id()) ==
MC::MUON ||
349 std::abs(( * thisChild2) -> pdg_id()) ==
MC::PIPLUS || std::abs(( * thisChild2) -> pdg_id()) ==
MC::KPLUS) {
350
351 ++NumChildD0;
352
353 if ((( * thisChild2) ->
momentum().
perp() >= m_PtMinKpi) &&
355 (std::abs(( * thisChild2) ->
momentum().
eta()) < m_EtaRangeKpi)) {
356 ++NumChildD0Charged;
357
358 if (NumChildD0Charged == 1) {
359 D0Child1 = ( * thisChild2);
362 ChargeD0Child1 = +1;
363 } else {
364 ChargeD0Child1 = -1;
365 }
366 }
367 if (NumChildD0Charged == 2) {
368 D0Child2 = ( * thisChild2);
371 ChargeD0Child2 = +1;
372 } else {
373 ChargeD0Child2 = -1;
374 }
375 }
376 if (std::abs(( * thisChild2) -> pdg_id()) ==
MC::MUON) {
377 ++NumChildD0mu;
378 D0ChildMu = ( * thisChild2);
379 }
380 }
381
382 }
else if (std::abs(( * thisChild2) -> pdg_id()) ==
MC::PI0) {
383 ++NumChildD0;
384 }
else if (std::abs(( * thisChild2) -> pdg_id()) ==
MC::NU_E || std::abs(( * thisChild2) -> pdg_id()) ==
MC::NU_MU) {
385 ++NumChildD0neutrinos;
386 }
else if (std::abs(( * thisChild2) -> pdg_id()) ==
MC::PHOTON) {
387 ++NumChildD0gammas;
388 }
else if (std::abs(( * thisChild2) -> pdg_id()) ==
MC::K0 || std::abs(( * thisChild2) -> pdg_id()) ==
MC::K0L ||
389 std::abs(( * thisChild2) -> pdg_id()) ==
MC::K0S) {
390 ++NumChildD0;
391 ++NumChildD0;
392 } else if ((( * thisChild2) -> end_vertex())) {
393 ++NumChildD0;
394 ++NumChildD0;
395 } else {
396 ++NumChildD0;
397 ++NumChildD0;
398 ATH_MSG_DEBUG(
"MuDstarFilter: unexpected D0 granddaughter = " << (*thisChild2)->pdg_id() );
399 }
400 }
401
402 } else {
403 ++NumChildD0;
404 ++NumChildD0;
405 ATH_MSG_DEBUG(
"MuDstarFilter: unexpected D0 daughter = " << (*thisChild1)->pdg_id() );
406 }
407 }
408
409 ATH_MSG_DEBUG(
"MuDstarFilter: NumChildD0, NumChildD0Charged = " << NumChildD0 <<
" , " << NumChildD0Charged);
410
411 if (NumChildD0 <= 3 && NumChildD0Charged == 2 && ChargeD0Child1 * ChargeD0Child2 < 0) {
412 if (m_D0Kpi_only) {
413 if (NumChildD0 == 2 && NumChildD0K == 1 && NumChildD0pi == 1) ++NumD0;
414 } else {
415 ++NumD0;
416 }
417 }
418 }
419 }
420
422
423 if (NumD0 == 1) {
424
425 if (pis_pdg * ChargeD0Child1 < 0) {
426
429
430 } else {
431
434 }
435
436 p4_D0 = p4_K + p4_pi;
437 double mKpi = p4_D0.M();
438
440
441 if (mKpi >= m_mKpiMin && mKpi <= m_mKpiMax) {
442
443 p4_Dstar = p4_D0 + p4_pis;
444
445 double delta_m = p4_Dstar.M() - mKpi;
446
448
449 if (delta_m <= m_delta_m_Max) {
450 ++NumDstars;
451
453
454 for (
size_t i = 0;
i < Muons.size(); ++
i) {
455
456 if (NumChildD0mu == 1) {
457 ATH_MSG_DEBUG(
"MuDstarFilter: Mu(pT), D0Mu(pT) = " << Muons[i]->
momentum().
perp() <<
" , " << D0ChildMu->momentum().perp());
458 if (std::fabs(Muons[i] ->
momentum().
perp() - D0ChildMu ->
momentum().
perp())< std::numeric_limits<double>::epsilon())
continue;
459 ATH_MSG_DEBUG(
"MuDstarFilter: Mu(pT), D0Mu(pT) = " << Muons[i]->
momentum().
perp() <<
" , " << D0ChildMu->momentum().perp() ) ;
460 }
461
463
464 p4_DstarMu = p4_Dstar + p4_Mu;
465
466 ATH_MSG_DEBUG(
"MuDstarFilter: p4_DstarMu.M() = " << p4_DstarMu.M());
467
468 if (p4_DstarMu.M() <= m_DstarMu_m_Max) {
469
470 ATH_MSG_INFO(
"MuDstarFilter: MuDstar candidate found" );
471 ATH_MSG_INFO(
"MuDstarFilter: p4_DstarMu.M() = " << p4_DstarMu.M() );
472 ATH_MSG_INFO(
"MuDstarFilter: NumChildD0, NumChildD0Charged = " << NumChildD0 <<
" , " << NumChildD0Charged );
473 ATH_MSG_INFO(
"MuDstarFilter: NumChildD0K, NumChildD0pi, NumChildD0mu = " << NumChildD0K <<
" , " << NumChildD0pi <<
" , " << NumChildD0mu );
474
475 if (NumChildD0mu == 1) {
476
477 ATH_MSG_DEBUG(
"MuDstarFilter: Mu(pT), D0Mu(pT) = " << Muons[i]->
momentum().
perp() <<
" , " << D0ChildMu->momentum().perp() );
478 }
479
480 ATH_MSG_INFO(
"MuDstarFilter: NumChildD0neutrinos, NumChildD0gammas = " << NumChildD0neutrinos <<
" , " << NumChildD0gammas );
481 ATH_MSG_INFO(
"MuDstarFilter: pis_pdg, K_pdg, ChargeD0Child1, ChargeD0Child2 = " << pis_pdg <<
" , " << K_pdg <<
" , " << ChargeD0Child1 <<
" , " << ChargeD0Child2 );
482
483 setFilterPassed(true);
484 return StatusCode::SUCCESS;
485 }
486 }
487
488 }
489 }
490 }
491 }
492 }
493 }
494 }
495 }
496
497
498
499
500 setFilterPassed(false);
501 return StatusCode::SUCCESS;
502}
Scalar eta() const
pseudorapidity method
Scalar perp() const
perp method - perpendicular length
Scalar phi() const
phi method
constexpr int pow(int base, int exp) noexcept
DataModel_detail::const_iterator< DataVector > const_iterator
HepMC::GenVertex * GenVertexPtr
const GenParticle * ConstGenParticlePtr
const HepMC::GenVertex * ConstGenVertexPtr
static const int ELECTRON
bool isPhysical(const T &p)
Identify if the particle is physical, i.e. is stable or decayed.