126 {
127
128
129
131
133
134 for (itr = events_const(ctx) -> begin(); itr != events_const(ctx) -> end(); ++itr) {
135
136 double primx = 0.;
137 double primy = 0.;
138
140
141 primx = vprim -> position().x();
142 primy = vprim -> position().y();
143
144 ATH_MSG_DEBUG(
"MuDstarFilter: PV x, y = " << primx <<
" , " << primy);
145
146 int NumMuons = 0;
147 int NumDstars = 0;
148
149 std::vector < HepMC::ConstGenParticlePtr >
Muons;
150
151
153 for (const auto & pitr: * genEvt) {
154
156
157
159 if ((pitr -> momentum().
perp() >= m_PtMinMuon) &&
160 (pitr -> momentum().
perp() < m_PtMaxMuon) &&
161 (std::abs(pitr -> momentum().
eta()) < m_EtaRangeMuon)) {
162 ++NumMuons;
163 Muons.push_back(pitr);
164 }
165 }
166 }
167
169
170 if (NumMuons == 0) break;
171 for (const auto & pitr: * genEvt) {
172
174
175
176 if (std::abs(pitr -> pdg_id()) ==
MC::DSTAR) {
179 (std::abs(pitr ->
momentum().
eta()) < m_EtaRangeDstar)) {
180
181
182 if (!(pitr -> end_vertex())) continue;
183
184 double Rxy = std::sqrt(
pow(pitr -> end_vertex() ->
position().
x() - primx, 2) +
pow(pitr -> end_vertex() ->
position().
y() - primy, 2));
185
187
188 if (Rxy < m_RxyMinDstar) continue;
189
190
191 auto firstChild = pitr -> end_vertex() -> particles_out().begin();
192 auto endChild = pitr -> end_vertex() -> particles_out().end();
193 auto thisChild = firstChild;
194
195 if (( * firstChild) -> pdg_id() == pitr -> pdg_id()) continue;
196
197 TLorentzVector p4_K;
198 TLorentzVector p4_pi;
199 TLorentzVector p4_pis;
200 TLorentzVector p4_D0;
201 TLorentzVector p4_Dstar;
202 TLorentzVector p4_Mu;
203 TLorentzVector p4_DstarMu;
204
205 int pis_pdg = 0;
206 int K_pdg = 0;
207
208 int NumD0 = 0;
209 int NumPis = 0;
210
211 for (; thisChild != endChild; ++thisChild) {
212
214
215 if (std::abs(( * thisChild) -> pdg_id()) ==
MC::PIPLUS) {
216 if ((( * thisChild) ->
momentum().
perp() >= m_PtMinPis) &&
218 (std::abs(( * thisChild) ->
momentum().
eta()) < m_EtaRangePis)) {
219 ++NumPis;
220
222 pis_pdg = ( * thisChild) -> pdg_id();
223 }
224 }
225 }
226
228
229 if (NumPis != 1) continue;
230
234
235 int NumChildD0 = 0;
236 int NumChildD0Charged = 0;
237 int NumChildD0neutrinos = 0;
238 int NumChildD0gammas = 0;
239 int ChargeD0Child1 = 0;
240 int ChargeD0Child2 = 0;
241 int NumChildD0K = 0;
242 int NumChildD0pi = 0;
243 int NumChildD0mu = 0;
244
245 for (thisChild = firstChild; thisChild != endChild; ++thisChild) {
246
248
249 if (std::abs(( * thisChild) -> pdg_id()) ==
MC::D0) {
250 if ((( * thisChild) -> end_vertex())) {
251 auto firstChild1 = ( * thisChild) -> end_vertex() -> particles_out().begin();
252 auto endChild1 = ( * thisChild) -> end_vertex() -> particles_out().end();
253 auto thisChild1 = firstChild1;
254
255 for (; thisChild1 != endChild1; ++thisChild1) {
256
258
259 if (std::abs(( * thisChild1) -> pdg_id()) ==
MC::ELECTRON || std::abs(( * thisChild1) -> pdg_id()) ==
MC::MUON ||
260 std::abs(( * thisChild1) -> pdg_id()) ==
MC::PIPLUS || std::abs(( * thisChild1) -> pdg_id()) ==
MC::KPLUS) {
261
262 ++NumChildD0;
263
264 if ((( * thisChild1) ->
momentum().
perp() >= m_PtMinKpi) &&
266 (std::abs(( * thisChild1) ->
momentum().
eta()) < m_EtaRangeKpi)) {
267 ++NumChildD0Charged;
268
269 if (NumChildD0Charged == 1) {
270 D0Child1 = ( * thisChild1);
273 ChargeD0Child1 = +1;
274 } else {
275 ChargeD0Child1 = -1;
276 }
277 }
278 if (NumChildD0Charged == 2) {
279 D0Child2 = ( * thisChild1);
282 ChargeD0Child2 = +1;
283 } else {
284 ChargeD0Child2 = -1;
285 }
286 }
287 if (std::abs(( * thisChild1) -> pdg_id()) ==
MC::MUON) {
288 ++NumChildD0mu;
289 D0ChildMu = ( * thisChild1);
290 }
291 if (std::abs(( * thisChild1) -> pdg_id()) ==
MC::PIPLUS) ++NumChildD0pi;
292 if (std::abs(( * thisChild1) -> pdg_id()) ==
MC::KPLUS) {
293 ++NumChildD0K;
294 K_pdg = ( * thisChild1) -> pdg_id();
295 }
296 }
297
298 }
else if (std::abs(( * thisChild1) -> pdg_id()) ==
MC::PI0) {
299 ++NumChildD0;
300 }
else if (std::abs(( * thisChild1) -> pdg_id()) ==
MC::NU_E || std::abs(( * thisChild1) -> pdg_id()) ==
MC::NU_MU) {
301 ++NumChildD0neutrinos;
302 }
else if (std::abs(( * thisChild1) -> pdg_id()) ==
MC::PHOTON) {
303 ++NumChildD0gammas;
304 }
else if (std::abs(( * thisChild1) -> pdg_id()) ==
MC::K0 || std::abs(( * thisChild1) -> pdg_id()) ==
MC::K0L ||
305 std::abs(( * thisChild1) -> pdg_id()) ==
MC::K0S) {
306 ++NumChildD0;
307 ++NumChildD0;
308 } else if ((( * thisChild1) -> end_vertex())) {
309
310
311 auto firstChild2 = ( * thisChild1) -> end_vertex() -> particles_out().begin();
312 auto endChild2 = ( * thisChild1) -> end_vertex() -> particles_out().end();
313 auto thisChild2 = firstChild2;
314 for (; thisChild2 != endChild2; ++thisChild2) {
315
317
318 if (std::abs(( * thisChild2) -> pdg_id()) ==
MC::ELECTRON || std::abs(( * thisChild2) -> pdg_id()) ==
MC::MUON ||
319 std::abs(( * thisChild2) -> pdg_id()) ==
MC::PIPLUS || std::abs(( * thisChild2) -> pdg_id()) ==
MC::KPLUS) {
320
321 ++NumChildD0;
322
323 if ((( * thisChild2) ->
momentum().
perp() >= m_PtMinKpi) &&
325 (std::abs(( * thisChild2) ->
momentum().
eta()) < m_EtaRangeKpi)) {
326 ++NumChildD0Charged;
327
328 if (NumChildD0Charged == 1) {
329 D0Child1 = ( * thisChild2);
332 ChargeD0Child1 = +1;
333 } else {
334 ChargeD0Child1 = -1;
335 }
336 }
337 if (NumChildD0Charged == 2) {
338 D0Child2 = ( * thisChild2);
341 ChargeD0Child2 = +1;
342 } else {
343 ChargeD0Child2 = -1;
344 }
345 }
346 if (std::abs(( * thisChild2) -> pdg_id()) ==
MC::MUON) {
347 ++NumChildD0mu;
348 D0ChildMu = ( * thisChild2);
349 }
350 }
351
352 }
else if (std::abs(( * thisChild2) -> pdg_id()) ==
MC::PI0) {
353 ++NumChildD0;
354 }
else if (std::abs(( * thisChild2) -> pdg_id()) ==
MC::NU_E || std::abs(( * thisChild2) -> pdg_id()) ==
MC::NU_MU) {
355 ++NumChildD0neutrinos;
356 }
else if (std::abs(( * thisChild2) -> pdg_id()) ==
MC::PHOTON) {
357 ++NumChildD0gammas;
358 }
else if (std::abs(( * thisChild2) -> pdg_id()) ==
MC::K0 || std::abs(( * thisChild2) -> pdg_id()) ==
MC::K0L ||
359 std::abs(( * thisChild2) -> pdg_id()) ==
MC::K0S) {
360 ++NumChildD0;
361 ++NumChildD0;
362 } else if ((( * thisChild2) -> end_vertex())) {
363 ++NumChildD0;
364 ++NumChildD0;
365 } else {
366 ++NumChildD0;
367 ++NumChildD0;
368 ATH_MSG_DEBUG(
"MuDstarFilter: unexpected D0 granddaughter = " << (*thisChild2)->pdg_id() );
369 }
370 }
371
372 } else {
373 ++NumChildD0;
374 ++NumChildD0;
375 ATH_MSG_DEBUG(
"MuDstarFilter: unexpected D0 daughter = " << (*thisChild1)->pdg_id() );
376 }
377 }
378
379 ATH_MSG_DEBUG(
"MuDstarFilter: NumChildD0, NumChildD0Charged = " << NumChildD0 <<
" , " << NumChildD0Charged);
380
381 if (NumChildD0 <= 3 && NumChildD0Charged == 2 && ChargeD0Child1 * ChargeD0Child2 < 0) {
382 if (m_D0Kpi_only) {
383 if (NumChildD0 == 2 && NumChildD0K == 1 && NumChildD0pi == 1) ++NumD0;
384 } else {
385 ++NumD0;
386 }
387 }
388 }
389 }
390
392
393 if (NumD0 == 1) {
394
395 if (pis_pdg * ChargeD0Child1 < 0) {
396
399
400 } else {
401
404 }
405
406 p4_D0 = p4_K + p4_pi;
407 double mKpi = p4_D0.M();
408
410
411 if (mKpi >= m_mKpiMin && mKpi <= m_mKpiMax) {
412
413 p4_Dstar = p4_D0 + p4_pis;
414
415 double delta_m = p4_Dstar.M() - mKpi;
416
418
419 if (delta_m <= m_delta_m_Max) {
420 ++NumDstars;
421
423
424 for (
size_t i = 0;
i < Muons.size(); ++
i) {
425
426 if (NumChildD0mu == 1) {
427 ATH_MSG_DEBUG(
"MuDstarFilter: Mu(pT), D0Mu(pT) = " << Muons[i]->
momentum().
perp() <<
" , " << D0ChildMu->momentum().perp());
428 if (std::fabs(Muons[i] ->
momentum().
perp() - D0ChildMu ->
momentum().
perp())< std::numeric_limits<double>::epsilon())
continue;
429 ATH_MSG_DEBUG(
"MuDstarFilter: Mu(pT), D0Mu(pT) = " << Muons[i]->
momentum().
perp() <<
" , " << D0ChildMu->momentum().perp() ) ;
430 }
431
433
434 p4_DstarMu = p4_Dstar + p4_Mu;
435
436 ATH_MSG_DEBUG(
"MuDstarFilter: p4_DstarMu.M() = " << p4_DstarMu.M());
437
438 if (p4_DstarMu.M() <= m_DstarMu_m_Max) {
439
440 ATH_MSG_INFO(
"MuDstarFilter: MuDstar candidate found" );
441 ATH_MSG_INFO(
"MuDstarFilter: p4_DstarMu.M() = " << p4_DstarMu.M() );
442 ATH_MSG_INFO(
"MuDstarFilter: NumChildD0, NumChildD0Charged = " << NumChildD0 <<
" , " << NumChildD0Charged );
443 ATH_MSG_INFO(
"MuDstarFilter: NumChildD0K, NumChildD0pi, NumChildD0mu = " << NumChildD0K <<
" , " << NumChildD0pi <<
" , " << NumChildD0mu );
444
445 if (NumChildD0mu == 1) {
446
447 ATH_MSG_DEBUG(
"MuDstarFilter: Mu(pT), D0Mu(pT) = " << Muons[i]->
momentum().
perp() <<
" , " << D0ChildMu->momentum().perp() );
448 }
449
450 ATH_MSG_INFO(
"MuDstarFilter: NumChildD0neutrinos, NumChildD0gammas = " << NumChildD0neutrinos <<
" , " << NumChildD0gammas );
451 ATH_MSG_INFO(
"MuDstarFilter: pis_pdg, K_pdg, ChargeD0Child1, ChargeD0Child2 = " << pis_pdg <<
" , " << K_pdg <<
" , " << ChargeD0Child1 <<
" , " << ChargeD0Child2 );
452
454 return StatusCode::SUCCESS;
455 }
456 }
457
458 }
459 }
460 }
461 }
462 }
463 }
464 }
465 }
466
467
468
469
471 return StatusCode::SUCCESS;
472}
Scalar eta() const
pseudorapidity method
Scalar perp() const
perp method - perpendicular length
Scalar phi() const
phi method
virtual void setFilterPassed(bool state, const EventContext &ctx) const
DataModel_detail::const_iterator< DataVector > const_iterator
HepMC3::ConstGenParticlePtr ConstGenParticlePtr
HepMC3::ConstGenVertexPtr ConstGenVertexPtr
HepMC3::GenEvent GenEvent
static const int ELECTRON
bool isPhysical(const T &p)
Identify if the particle is physical, i.e. is stable or decayed.
const Amg::Vector3D & position() const
Method to retrieve the position of the Intersection.