Belle II Software development
KinkVariables.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9// Own header.
10#include <analysis/variables/KinkVariables.h>
11
12// dataobjects from the analysis
13#include <analysis/dataobjects/Particle.h>
14
15// dataobjects from the framework
16#include <framework/dataobjects/Helix.h>
17
18// dataobjects from the MDST
19#include <mdst/dataobjects/Track.h>
20#include <mdst/dataobjects/Kink.h>
21#include <mdst/dataobjects/MCParticle.h>
22#include <mdst/dataobjects/TrackFitResult.h>
23#include <mdst/dataobjects/HitPatternCDC.h>
24#include <mdst/dataobjects/HitPatternVXD.h>
25
26// framework aux
27#include <framework/logging/Logger.h>
28
29// extra
30#include <algorithm>
31#include <cmath>
32#include <Math/Boost.h>
33
34using namespace std;
35
36namespace Belle2 {
41 namespace Variable {
42
43 //-----------------------------------------------
44 // Helper Functions
45
46 Helix kinkDaughterHelixAtKinkVertex(const Kink* kink)
47 {
48 Helix daughterHelixAtKinkVertex = kink->getDaughterTrackFitResult()->getHelix();
49 daughterHelixAtKinkVertex.passiveMoveBy(kink->getFittedVertexX(), kink->getFittedVertexY(), 0);
50 return daughterHelixAtKinkVertex;
51 }
52
53 ROOT::Math::XYZVector kinkDaughterMomentumAtKinkVertex(const Kink* kink)
54 {
55 const double BzAtKinkVertex = BFieldManager::getFieldInTesla(
56 {kink->getFittedVertexX(), kink->getFittedVertexY(), kink->getFittedVertexZ()}).Z();
57 Helix daughterHelixAtKinkVertex = kink->getDaughterTrackFitResult()->getHelix();
58 daughterHelixAtKinkVertex.passiveMoveBy(kink->getFittedVertexX(), kink->getFittedVertexY(), 0);
59 return daughterHelixAtKinkVertex.getMomentum(BzAtKinkVertex);
60 }
61
62 Helix kinkMotherHelixAtKinkVertex(const Kink* kink)
63 {
64 Helix motherHelixAtKinkVertex = kink->getMotherTrackFitResultEnd()->getHelix();
65 motherHelixAtKinkVertex.passiveMoveBy(kink->getFittedVertexX(), kink->getFittedVertexY(), 0);
66 return motherHelixAtKinkVertex;
67 }
68
69 ROOT::Math::XYZVector kinkMotherMomentumAtKinkVertex(const Kink* kink)
70 {
71 const double BzAtKinkVertex = BFieldManager::getFieldInTesla(
72 {kink->getFittedVertexX(), kink->getFittedVertexY(), kink->getFittedVertexZ()}).Z();
73 Helix motherHelixAtKinkVertex = kink->getMotherTrackFitResultEnd()->getHelix();
74 motherHelixAtKinkVertex.passiveMoveBy(kink->getFittedVertexX(), kink->getFittedVertexY(), 0);
75 return motherHelixAtKinkVertex.getMomentum(BzAtKinkVertex);
76 }
77
78 double kinkDaughterMomentumAndCosThetaInMotherRF(const Particle* part, Const::ChargedStable motherType,
79 Const::ChargedStable daughterType, bool returnCosTheta)
80 {
81 const Kink* kink = part->getKink();
82 if (!kink) return Const::doubleNaN;
83
84 ROOT::Math::XYZVector motherMomentumAtKinkVertex = kinkMotherMomentumAtKinkVertex(kink);
85 ROOT::Math::XYZVector daughterMomentumAtKinkVertex = kinkDaughterMomentumAtKinkVertex(kink);
86
87 double motherEnergy = sqrt(motherMomentumAtKinkVertex.Mag2() + motherType.getMass() * motherType.getMass());
88 double daughterEnergy = sqrt(daughterMomentumAtKinkVertex.Mag2() + daughterType.getMass() * daughterType.getMass());
89
90 ROOT::Math::PxPyPzEVector mother4MomentumAtKinkVertex(motherMomentumAtKinkVertex.X(),
91 motherMomentumAtKinkVertex.Y(),
92 motherMomentumAtKinkVertex.Z(),
93 motherEnergy);
94 int signDiff = kink->getMotherTrackFitResultEnd()->getChargeSign() * kink->getDaughterTrackFitResult()->getChargeSign();
95 ROOT::Math::PxPyPzEVector daughter4MomentumAtKinkVertex(signDiff * daughterMomentumAtKinkVertex.X(),
96 signDiff * daughterMomentumAtKinkVertex.Y(),
97 signDiff * daughterMomentumAtKinkVertex.Z(),
98 daughterEnergy);
99
100 ROOT::Math::XYZVector motherBoostAtKinkVertex = mother4MomentumAtKinkVertex.BoostToCM();
101 daughter4MomentumAtKinkVertex = ROOT::Math::Boost(motherBoostAtKinkVertex) * daughter4MomentumAtKinkVertex;
102
103 if (returnCosTheta)
104 return daughter4MomentumAtKinkVertex.Vect().Unit().Dot(mother4MomentumAtKinkVertex.Vect().Unit());
105 else
106 return sqrt(daughter4MomentumAtKinkVertex.Vect().Mag2());
107 }
108
109 ROOT::Math::PxPyPzEVector kinkMotherMCP4AtDecayVertex(const MCParticle* p)
110 {
111 ROOT::Math::PxPyPzEVector P4(0, 0, 0, 0);
112 ROOT::Math::XYZVector mcDecayVertex = p->getDecayVertex();
113 std::vector<MCParticle*> mcDaughters = p->getDaughters();
114 for (std::vector<MCParticle*>::iterator daughterIter = mcDaughters.begin();
115 daughterIter != mcDaughters.end(); ++daughterIter) {
116 ROOT::Math::XYZVector mcDaughterVertex = (*daughterIter)->getVertex();
117 if ((mcDaughterVertex - mcDecayVertex).Mag2() < 1)
118 P4 += (*daughterIter)->get4Vector();
119 }
120 return P4;
121 }
122
123 double kinkDaughterMCMomentumAndCosThetaInMotherRF(const Particle* part, bool returnCosTheta)
124 {
125 const Kink* kink = part->getKink();
126 if (!kink) return Const::doubleNaN;
127
128 const Track* motherTrack = kink->getMotherTrack();
129 const Track* daughterTrack = kink->getDaughterTrack();
130 if (motherTrack == daughterTrack) return Const::doubleNaN;
131
132 const MCParticle* motherMCParticle = motherTrack->getRelated<MCParticle>();
133 const MCParticle* daughterMCParticle = daughterTrack->getRelated<MCParticle>();
134
135 if (!motherMCParticle || !daughterMCParticle) return Const::doubleNaN;
136
137 ROOT::Math::PxPyPzEVector mother4MomentumAtDecayVertex = kinkMotherMCP4AtDecayVertex(motherMCParticle);
138 ROOT::Math::PxPyPzEVector daughter4MomentumAtProductionVertex = daughterMCParticle->get4Vector();
139
140 ROOT::Math::XYZVector motherBoostAtDecayVertex = mother4MomentumAtDecayVertex.BoostToCM();
141 ROOT::Math::PxPyPzEVector daughter4MomentumAtDecayVertex = ROOT::Math::Boost(motherBoostAtDecayVertex) *
142 daughter4MomentumAtProductionVertex;
143
144 if (returnCosTheta)
145 return daughter4MomentumAtDecayVertex.Vect().Unit().Dot(mother4MomentumAtDecayVertex.Vect().Unit());
146 else
147 return sqrt(daughter4MomentumAtDecayVertex.Vect().Mag2());
148 }
149
150 //-----------------------------------------------
151 // MEASURED VARIABLES
152 // Kink General Variables
153
154 double kinkVertexX(const Particle* part)
155 {
156 const Kink* kink = part->getKink();
157 if (!kink) return Const::doubleNaN;
158 return kink->getFittedVertexX();
159 }
160
161 double kinkVertexY(const Particle* part)
162 {
163 const Kink* kink = part->getKink();
164 if (!kink) return Const::doubleNaN;
165 return kink->getFittedVertexY();
166 }
167
168 double kinkVertexZ(const Particle* part)
169 {
170 const Kink* kink = part->getKink();
171 if (!kink) return Const::doubleNaN;
172 return kink->getFittedVertexZ();
173 }
174
175 double kinkFilterID(const Particle* part)
176 {
177 const Kink* kink = part->getKink();
178 if (!kink) return Const::doubleNaN;
179 return kink->getPrefilterFlag();
180 }
181
182 double kinkCombinedFitResultFlag(const Particle* part)
183 {
184 const Kink* kink = part->getKink();
185 if (!kink) return Const::doubleNaN;
186 if (kinkFilterID(part) < 3)
187 return kink->getCombinedFitResultFlag();
188 else
189 return Const::doubleNaN;
190 }
191
192 double kinkCombinedFitResultFlagBit1(const Particle* part)
193 {
194 double flag = kinkCombinedFitResultFlag(part);
195 if (flag < 16)
196 return static_cast<int>(flag) & 0b0001;
197 else
198 return Const::doubleNaN;
199 }
200
201 double kinkCombinedFitResultFlagBit2(const Particle* part)
202 {
203 double flag = kinkCombinedFitResultFlag(part);
204 if (flag < 16)
205 return static_cast<bool>(static_cast<int>(flag) & 0b0010);
206 else
207 return Const::doubleNaN;
208 }
209
210 double kinkCombinedFitResultFlagBit3(const Particle* part)
211 {
212 double flag = kinkCombinedFitResultFlag(part);
213 if (flag < 16)
214 return static_cast<bool>(static_cast<int>(flag) & 0b0100);
215 else
216 return Const::doubleNaN;
217 }
218
219 double kinkCombinedFitResultFlagBit4(const Particle* part)
220 {
221 double flag = kinkCombinedFitResultFlag(part);
222 if (flag < 16)
223 return static_cast<bool>(static_cast<int>(flag) & 0b1000);
224 else
225 return Const::doubleNaN;
226 }
227
228 double kinkSplitTrackDistanceAtVertexFlag(const Particle* part)
229 {
230 const Kink* kink = part->getKink();
231 if (!kink) return Const::doubleNaN;
232 if (kinkFilterID(part) < 3)
233 return Const::doubleNaN;
234 else
235 return kink->getSplitTrackDistanceAtVertexFlag();
236 }
237
238 double kinkNumberOfReassignedHits(const Particle* part)
239 {
240 const Kink* kink = part->getKink();
241 if (!kink) return Const::doubleNaN;
242 return kink->getNumberOfReassignedHits();
243 }
244
245 double kinkIsSameCharge(const Particle* part)
246 {
247 const Kink* kink = part->getKink();
248 if (!kink) return Const::doubleNaN;
249
250 return kink->getMotherTrackFitResultEnd()->getChargeSign() *
251 kink->getDaughterTrackFitResult()->getChargeSign() > 0 ? true : false;
252 }
253
254 // Kink Decay Kinematics
255
256 double kinkDaughterMomentumInMotherRF(const Particle* part)
257 {
258 if (!(part->hasExtraInfo("kinkDaughterPDGCode")))
259 return Const::doubleNaN;
260 Const::ChargedStable daughterType(abs(part->getExtraInfo("kinkDaughterPDGCode")));
261 Const::ChargedStable motherType(abs(part->getPDGCode()));
262 return kinkDaughterMomentumAndCosThetaInMotherRF(part, motherType, daughterType, false);
263 }
264
265 double kinkDaughterCosThetaInMotherRF(const Particle* part)
266 {
267 if (!(part->hasExtraInfo("kinkDaughterPDGCode")))
268 return Const::doubleNaN;
269 Const::ChargedStable daughterType(abs(part->getExtraInfo("kinkDaughterPDGCode")));
270 Const::ChargedStable motherType(abs(part->getPDGCode()));
271 return kinkDaughterMomentumAndCosThetaInMotherRF(part, motherType, daughterType, true);
272 }
273
274 double kinkDaughterMomentumInMotherRFKPi(const Particle* part)
275 {
276 return kinkDaughterMomentumAndCosThetaInMotherRF(part, Const::kaon, Const::pion, false);
277 }
278
279 double kinkDaughterCosThetaInMotherRFKPi(const Particle* part)
280 {
281 return kinkDaughterMomentumAndCosThetaInMotherRF(part, Const::kaon, Const::pion, true);
282 }
283
284 double kinkDaughterMomentumInMotherRFKMu(const Particle* part)
285 {
286 return kinkDaughterMomentumAndCosThetaInMotherRF(part, Const::kaon, Const::muon, false);
287 }
288
289 double kinkDaughterCosThetaInMotherRFKMu(const Particle* part)
290 {
291 return kinkDaughterMomentumAndCosThetaInMotherRF(part, Const::kaon, Const::muon, true);
292 }
293
294 double kinkDaughterMomentumInMotherRFPiMu(const Particle* part)
295 {
296 return kinkDaughterMomentumAndCosThetaInMotherRF(part, Const::pion, Const::muon, false);
297 }
298
299 double kinkDaughterCosThetaInMotherRFPiMu(const Particle* part)
300 {
301 return kinkDaughterMomentumAndCosThetaInMotherRF(part, Const::pion, Const::muon, true);
302 }
303
304 double kinkDaughterMomentumInMotherRFMuE(const Particle* part)
305 {
306 return kinkDaughterMomentumAndCosThetaInMotherRF(part, Const::muon, Const::electron, false);
307 }
308
309 double kinkDaughterCosThetaInMotherRFMuE(const Particle* part)
310 {
311 return kinkDaughterMomentumAndCosThetaInMotherRF(part, Const::muon, Const::electron, true);
312 }
313
314 // Kink Daughter Measured Track Parameters
315
316 Manager::FunctionPtr kinkDaughterTrack(const std::vector<std::string>& arguments)
317 {
318 if (arguments.size() == 1) {
319 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
320 auto func = [var](const Particle * particle) -> double {
321 if (particle->getParticleSource() != Particle::EParticleSourceObject::c_Kink) return Const::doubleNaN;
322 const Kink* kink = particle->getKink();
323 if (!kink) return Const::doubleNaN;
324 if (not particle->hasExtraInfo("kinkDaughterPDGCode")) return Const::doubleNaN;
325 Particle tmpParticle(kink, Const::ChargedStable(abs(particle->getExtraInfo("kinkDaughterPDGCode"))), kink->getTrackFitResultIndexDaughter());
326 auto var_result = var->function(&tmpParticle);
327 if (std::holds_alternative<double>(var_result))
328 {
329 return std::get<double>(var_result);
330 } else if (std::holds_alternative<int>(var_result))
331 {
332 return std::get<int>(var_result);
333 } else if (std::holds_alternative<bool>(var_result))
334 {
335 return std::get<bool>(var_result);
336 } else
337 {
338 return Const::doubleNaN;
339 }
340 };
341 return func;
342 } else {
343 B2FATAL("Wrong number of arguments for meta function kinkDaughterTrack");
344 }
345 }
346
347 Manager::FunctionPtr kinkDaughterInitTrack(const std::vector<std::string>& arguments)
348 {
349 if (arguments.size() == 1) {
350 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
351 auto func = [var](const Particle * particle) -> double {
352 if (particle->getParticleSource() != Particle::EParticleSourceObject::c_Kink) return Const::doubleNaN;
353 const Kink* kink = particle->getKink();
354 if (!kink) return Const::doubleNaN;
355 if (not particle->hasExtraInfo("kinkDaughterPDGCode")) return Const::doubleNaN;
356 Particle tmpParticle(kink->getDaughterTrack()->getArrayIndex(), kink->getDaughterTrack()->getTrackFitResultWithClosestMass(Const::pion), Const::ChargedStable(abs(particle->getExtraInfo("kinkDaughterPDGCode"))));
357 auto var_result = var->function(&tmpParticle);
358 if (std::holds_alternative<double>(var_result))
359 {
360 return std::get<double>(var_result);
361 } else if (std::holds_alternative<int>(var_result))
362 {
363 return std::get<int>(var_result);
364 } else if (std::holds_alternative<bool>(var_result))
365 {
366 return std::get<bool>(var_result);
367 } else
368 {
369 return Const::doubleNaN;
370 }
371 };
372 return func;
373 } else {
374 B2FATAL("Wrong number of arguments for meta function kinkDaughterInitTrack");
375 }
376 }
377
378 double kinkDaughterTrackD0AtKinkVertex(const Particle* part)
379 {
380 const Kink* kink = part->getKink();
381 if (!kink) return Const::doubleNaN;
382
383 Helix daughterHelixAtKinkVertex = kinkDaughterHelixAtKinkVertex(kink);
384 return daughterHelixAtKinkVertex.getD0();
385 }
386
387 double kinkDaughterTrackZ0AtKinkVertex(const Particle* part)
388 {
389 const Kink* kink = part->getKink();
390 if (!kink) return Const::doubleNaN;
391
392 Helix daughterHelixAtKinkVertex = kinkDaughterHelixAtKinkVertex(kink);
393 return daughterHelixAtKinkVertex.getZ0();
394 }
395
396 double kinkDaughterPtAtKinkVertex(const Particle* part)
397 {
398 const Kink* kink = part->getKink();
399 if (!kink) return Const::doubleNaN;
400
401 ROOT::Math::XYZVector daughterMomentumAtKinkVertex = kinkDaughterMomentumAtKinkVertex(kink);
402 return sqrt(daughterMomentumAtKinkVertex.Perp2());
403 }
404
405 double kinkDaughterPzAtKinkVertex(const Particle* part)
406 {
407 const Kink* kink = part->getKink();
408 if (!kink) return Const::doubleNaN;
409
410 ROOT::Math::XYZVector daughterMomentumAtKinkVertex = kinkDaughterMomentumAtKinkVertex(kink);
411 return daughterMomentumAtKinkVertex.Z();
412 }
413
414 double kinkDaughterPAtKinkVertex(const Particle* part)
415 {
416 const Kink* kink = part->getKink();
417 if (!kink) return Const::doubleNaN;
418
419 ROOT::Math::XYZVector daughterMomentumAtKinkVertex = kinkDaughterMomentumAtKinkVertex(kink);
420 return sqrt(daughterMomentumAtKinkVertex.Mag2());
421 }
422
423 // Kink Mother Measured Track Parameters
424
425 Manager::FunctionPtr kinkMotherInitTrack(const std::vector<std::string>& arguments)
426 {
427 if (arguments.size() == 1) {
428 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
429 auto func = [var](const Particle * particle) -> double {
430 if (particle->getParticleSource() != Particle::EParticleSourceObject::c_Kink) return Const::doubleNaN;
431 const Kink* kink = particle->getKink();
432 if (!kink) return Const::doubleNaN;
433 Particle tmpParticle(kink->getMotherTrack()->getArrayIndex(), kink->getMotherTrack()->getTrackFitResultWithClosestMass(Const::pion), Const::pion);
434 auto var_result = var->function(&tmpParticle);
435 if (std::holds_alternative<double>(var_result))
436 {
437 return std::get<double>(var_result);
438 } else if (std::holds_alternative<int>(var_result))
439 {
440 return std::get<int>(var_result);
441 } else if (std::holds_alternative<bool>(var_result))
442 {
443 return std::get<bool>(var_result);
444 } else
445 {
446 return Const::doubleNaN;
447 }
448 };
449 return func;
450 } else {
451 B2FATAL("Wrong number of arguments for meta function kinkMotherInitTrack");
452 }
453 }
454
455 double kinkMotherTrackD0AtKinkVertex(const Particle* part)
456 {
457 const Kink* kink = part->getKink();
458 if (!kink) return Const::doubleNaN;
459
460 Helix motherHelixAtKinkVertex = kinkMotherHelixAtKinkVertex(kink);
461 return motherHelixAtKinkVertex.getD0();
462 }
463
464 double kinkMotherTrackZ0AtKinkVertex(const Particle* part)
465 {
466 const Kink* kink = part->getKink();
467 if (!kink) return Const::doubleNaN;
468
469 Helix motherHelixAtKinkVertex = kinkMotherHelixAtKinkVertex(kink);
470 return motherHelixAtKinkVertex.getZ0();
471 }
472
473 double kinkMotherPtAtKinkVertex(const Particle* part)
474 {
475 const Kink* kink = part->getKink();
476 if (!kink) return Const::doubleNaN;
477
478 ROOT::Math::XYZVector motherMomentumAtKinkVertex = kinkMotherMomentumAtKinkVertex(kink);
479 return sqrt(motherMomentumAtKinkVertex.Perp2());
480 }
481
482 double kinkMotherPzAtKinkVertex(const Particle* part)
483 {
484 const Kink* kink = part->getKink();
485 if (!kink) return Const::doubleNaN;
486
487 ROOT::Math::XYZVector motherMomentumAtKinkVertex = kinkMotherMomentumAtKinkVertex(kink);
488 return motherMomentumAtKinkVertex.Z();
489 }
490
491 double kinkMotherPAtKinkVertex(const Particle* part)
492 {
493 const Kink* kink = part->getKink();
494 if (!kink) return Const::doubleNaN;
495
496 ROOT::Math::XYZVector motherMomentumAtKinkVertex = kinkMotherMomentumAtKinkVertex(kink);
497 return sqrt(motherMomentumAtKinkVertex.Mag2());
498 }
499
500 //-----------------------------------------------
501 // Check if the particle is a part of any Kink
502
503 double particleIsInKink(const Particle* part)
504 {
505 const Track* track = part->getTrack();
506 if (!track) return Const::doubleNaN;
507
508 const short trackIndex = track->getArrayIndex();
509 StoreArray<Kink> kinks;
510
511 for (int i = 0; i < kinks.getEntries(); i++) {
512 const Kink* kink = kinks[i];
513 const short motherTrackIndex = kink->getMotherTrackIndex();
514 const short daughterTrackIndex = kink->getDaughterTrackIndex();
515 if (trackIndex == motherTrackIndex || trackIndex == daughterTrackIndex)
516 return true;
517 }
518 return false;
519 }
520
521
522 double particleIsMotherInKink(const Particle* part)
523 {
524 const Track* track = part->getTrack();
525 if (!track) return Const::doubleNaN;
526
527 const short trackIndex = track->getArrayIndex();
528 StoreArray<Kink> kinks;
529
530 for (int i = 0; i < kinks.getEntries(); i++) {
531 const Kink* kink = kinks[i];
532 const short motherTrackIndex = kink->getMotherTrackIndex();
533 if (trackIndex == motherTrackIndex)
534 return true;
535 }
536 return false;
537 }
538
539
540 double particleIsDaughterInKink(const Particle* part)
541 {
542 const Track* track = part->getTrack();
543 if (!track) return Const::doubleNaN;
544
545 const short trackIndex = track->getArrayIndex();
546 StoreArray<Kink> kinks;
547
548 for (int i = 0; i < kinks.getEntries(); i++) {
549 const Kink* kink = kinks[i];
550 const short daughterTrackIndex = kink->getDaughterTrackIndex();
551 if (trackIndex == daughterTrackIndex)
552 return true;
553 }
554 return false;
555 }
556
557 double particleIsSplitKink(const Particle* part)
558 {
559 const Track* track = part->getTrack();
560 if (!track) return Const::doubleNaN;
561
562 const short trackIndex = track->getArrayIndex();
563 StoreArray<Kink> kinks;
564
565 for (int i = 0; i < kinks.getEntries(); i++) {
566 const Kink* kink = kinks[i];
567 const short motherTrackIndex = kink->getMotherTrackIndex();
568 const short daughterTrackIndex = kink->getDaughterTrackIndex();
569 if (trackIndex == motherTrackIndex && motherTrackIndex == daughterTrackIndex)
570 return true;
571 }
572 return false;
573 }
574
575 //-----------------------------------------------
576 // MC VARIABLES
577 // Kink from track pair MC variables
578
579 double kinkPairIsMCRelated(const Particle* part)
580 {
581 const Kink* kink = part->getKink();
582 if (!kink) return Const::doubleNaN;
583
584 const Track* motherTrack = kink->getMotherTrack();
585 const Track* daughterTrack = kink->getDaughterTrack();
586
587 if (motherTrack == daughterTrack) return Const::doubleNaN;
588
589 const MCParticle* motherMCParticle = motherTrack->getRelated<MCParticle>();
590 const MCParticle* daughterMCParticle = daughterTrack->getRelated<MCParticle>();
591
592 if (!motherMCParticle || !daughterMCParticle) return false;
593 if (motherMCParticle == daughterMCParticle) return true;
594 if (daughterMCParticle->getMother() && daughterMCParticle->getMother() == motherMCParticle) return true;
595 return false;
596 }
597
598 double kinkPairIsClone(const Particle* part)
599 {
600 const Kink* kink = part->getKink();
601 if (!kink) return Const::doubleNaN;
602
603 const Track* motherTrack = kink->getMotherTrack();
604 const Track* daughterTrack = kink->getDaughterTrack();
605
606 if (motherTrack == daughterTrack) return Const::doubleNaN;
607
608 const MCParticle* motherMCParticle = motherTrack->getRelated<MCParticle>();
609 const MCParticle* daughterMCParticle = daughterTrack->getRelated<MCParticle>();
610
611 if (!motherMCParticle || !daughterMCParticle) return false;
612 if (motherMCParticle == daughterMCParticle) return true;
613 return false;
614 }
615
616 double kinkPairIsReal(const Particle* part)
617 {
618 const Kink* kink = part->getKink();
619 if (!kink) return Const::doubleNaN;
620
621 const Track* motherTrack = kink->getMotherTrack();
622 const Track* daughterTrack = kink->getDaughterTrack();
623
624 if (motherTrack == daughterTrack) return Const::doubleNaN;
625
626 const MCParticle* motherMCParticle = motherTrack->getRelated<MCParticle>();
627 const MCParticle* daughterMCParticle = daughterTrack->getRelated<MCParticle>();
628
629 if (!motherMCParticle || !daughterMCParticle) return false;
630 if (daughterMCParticle->getMother() && daughterMCParticle->getMother() == motherMCParticle) return true;
631 return false;
632 }
633
634 double kinkPairIsDecayInFlight(const Particle* part)
635 {
636 const Kink* kink = part->getKink();
637 if (!kink) return Const::doubleNaN;
638
639 const Track* motherTrack = kink->getMotherTrack();
640 const Track* daughterTrack = kink->getDaughterTrack();
641
642 if (motherTrack == daughterTrack) return Const::doubleNaN;
643
644 const MCParticle* motherMCParticle = motherTrack->getRelated<MCParticle>();
645 const MCParticle* daughterMCParticle = daughterTrack->getRelated<MCParticle>();
646
647 if (!motherMCParticle || !daughterMCParticle) return false;
648 if (daughterMCParticle->getMother() &&
649 daughterMCParticle->getMother() == motherMCParticle &&
650 daughterMCParticle->getSecondaryPhysicsProcess() == 201) return true;
651 return false;
652 }
653
654 double kinkPairIsHadronScattering(const Particle* part)
655 {
656 const Kink* kink = part->getKink();
657 if (!kink) return Const::doubleNaN;
658
659 const Track* motherTrack = kink->getMotherTrack();
660 const Track* daughterTrack = kink->getDaughterTrack();
661
662 if (motherTrack == daughterTrack) return Const::doubleNaN;
663
664 const MCParticle* motherMCParticle = motherTrack->getRelated<MCParticle>();
665 const MCParticle* daughterMCParticle = daughterTrack->getRelated<MCParticle>();
666
667 if (!motherMCParticle || !daughterMCParticle) return false;
668 if (daughterMCParticle->getMother() &&
669 daughterMCParticle->getMother() == motherMCParticle &&
670 daughterMCParticle->getSecondaryPhysicsProcess() == 121) return true;
671 return false;
672 }
673
674 // Kink from track pair MC kinematics
675
676 double kinkDaughterMomentumInMotherRFMC(const Particle* part)
677 {
678 return kinkDaughterMCMomentumAndCosThetaInMotherRF(part, false);
679 }
680
681 double kinkDaughterCosThetaInMotherRFMC(const Particle* part)
682 {
683 return kinkDaughterMCMomentumAndCosThetaInMotherRF(part, true);
684 }
685
686 // Kink from track pair daughter MC variables
687
688 Manager::FunctionPtr kinkPairDaughterMC(const std::vector<std::string>& arguments)
689 {
690 if (arguments.size() == 1) {
691 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
692 auto func = [var](const Particle * particle) -> double {
693 if (particle->getParticleSource() != Particle::EParticleSourceObject::c_Kink) return Const::doubleNaN;
694 const Kink* kink = particle->getKink();
695 if (!kink) return Const::doubleNaN;
696 const Track* daughterTrack = kink->getDaughterTrack();
697 if (daughterTrack == kink->getMotherTrack()) return Const::doubleNaN;
698 Particle tmpParticle(daughterTrack->getRelated<MCParticle>());
699 auto var_result = var->function(&tmpParticle);
700 if (std::holds_alternative<double>(var_result))
701 {
702 return std::get<double>(var_result);
703 } else if (std::holds_alternative<int>(var_result))
704 {
705 return std::get<int>(var_result);
706 } else if (std::holds_alternative<bool>(var_result))
707 {
708 return std::get<bool>(var_result);
709 } else
710 {
711 return Const::doubleNaN;
712 }
713 };
714 return func;
715 } else {
716 B2FATAL("Wrong number of arguments for meta function kinkPairDaughterMC");
717 }
718 }
719
720 // Kink from track pair mother MC variables
721
722 Manager::FunctionPtr kinkPairMotherMC(const std::vector<std::string>& arguments)
723 {
724 if (arguments.size() == 1) {
725 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
726 auto func = [var](const Particle * particle) -> double {
727 if (particle->getParticleSource() != Particle::EParticleSourceObject::c_Kink) return Const::doubleNaN;
728 const Kink* kink = particle->getKink();
729 if (!kink) return Const::doubleNaN;
730 const Track* motherTrack = kink->getMotherTrack();
731 if (motherTrack == kink->getDaughterTrack()) return Const::doubleNaN;
732 Particle tmpParticle(motherTrack->getRelated<MCParticle>());
733 auto var_result = var->function(&tmpParticle);
734 if (std::holds_alternative<double>(var_result))
735 {
736 return std::get<double>(var_result);
737 } else if (std::holds_alternative<int>(var_result))
738 {
739 return std::get<int>(var_result);
740 } else if (std::holds_alternative<bool>(var_result))
741 {
742 return std::get<bool>(var_result);
743 } else
744 {
745 return Const::doubleNaN;
746 }
747 };
748 return func;
749 } else {
750 B2FATAL("Wrong number of arguments for meta function kinkPairMotherMC");
751 }
752 }
753
754 const MCParticle* kinkPairMotherMCParticle(const Particle* part)
755 {
756 const Kink* kink = part->getKink();
757 if (!kink) return nullptr;
758
759 const Track* motherTrack = kink->getMotherTrack();
760 const Track* daughterTrack = kink->getDaughterTrack();
761
762 if (motherTrack == daughterTrack) return nullptr;
763
764 return motherTrack->getRelated<MCParticle>();
765 }
766
767 double kinkPairMotherMCPXAtDecayVertex(const Particle* part)
768 {
769 const MCParticle* motherMCParticle = kinkPairMotherMCParticle(part);
770 if (!motherMCParticle) return Const::doubleNaN;
771 return kinkMotherMCP4AtDecayVertex(motherMCParticle).X();
772 }
773
774 double kinkPairMotherMCPYAtDecayVertex(const Particle* part)
775 {
776 const MCParticle* motherMCParticle = kinkPairMotherMCParticle(part);
777 if (!motherMCParticle) return Const::doubleNaN;
778 return kinkMotherMCP4AtDecayVertex(motherMCParticle).Y();
779 }
780
781 double kinkPairMotherMCPZAtDecayVertex(const Particle* part)
782 {
783 const MCParticle* motherMCParticle = kinkPairMotherMCParticle(part);
784 if (!motherMCParticle) return Const::doubleNaN;
785 return kinkMotherMCP4AtDecayVertex(motherMCParticle).Z();
786 }
787
788 double kinkPairMotherMCPTAtDecayVertex(const Particle* part)
789 {
790 const MCParticle* motherMCParticle = kinkPairMotherMCParticle(part);
791 if (!motherMCParticle) return Const::doubleNaN;
792 return sqrt(kinkMotherMCP4AtDecayVertex(motherMCParticle).Vect().Perp2());
793 }
794
795 double kinkPairMotherMCPAtDecayVertex(const Particle* part)
796 {
797 const MCParticle* motherMCParticle = kinkPairMotherMCParticle(part);
798 if (!motherMCParticle) return Const::doubleNaN;
799 return sqrt(kinkMotherMCP4AtDecayVertex(motherMCParticle).Vect().Mag2());
800 }
801
802 double kinkPairMotherMCEAtDecayVertex(const Particle* part)
803 {
804 const MCParticle* motherMCParticle = kinkPairMotherMCParticle(part);
805 if (!motherMCParticle) return Const::doubleNaN;
806 return kinkMotherMCP4AtDecayVertex(motherMCParticle).E();
807 }
808
809 //-----------------------------------------------
810 // REGISTER VARIABLES
811
812 VARIABLE_GROUP("Kink");
813
814 //-----------------------------------------------
815 // MEASURED VARIABLES
816 // Kink General Variables
817 REGISTER_VARIABLE("kinkVertexX", kinkVertexX, "x coordinate of kink vertex");
818 REGISTER_VARIABLE("kinkVertexY", kinkVertexY, "y coordinate of kink vertex");
819 REGISTER_VARIABLE("kinkVertexZ", kinkVertexZ, "z coordinate of kink vertex");
820 REGISTER_VARIABLE("kinkFilterID", kinkFilterID, "Filter ID with which kink was preselected");
821 REGISTER_VARIABLE("kinkCombinedFitResultFlag", kinkCombinedFitResultFlag, "Flag of the combined kink fit result");
822 REGISTER_VARIABLE("kinkCombinedFitResultFlagB1", kinkCombinedFitResultFlagBit1,
823 "The first bit of the flag of the combined kink fit result");
824 REGISTER_VARIABLE("kinkCombinedFitResultFlagB2", kinkCombinedFitResultFlagBit2,
825 "The second bit of the flag of the combined kink fit result");
826 REGISTER_VARIABLE("kinkCombinedFitResultFlagB3", kinkCombinedFitResultFlagBit3,
827 "The third bit of the flag of the combined kink fit result");
828 REGISTER_VARIABLE("kinkCombinedFitResultFlagB4", kinkCombinedFitResultFlagBit4,
829 "The fourth bit of the flag of the combined kink fit result");
830 REGISTER_VARIABLE("kinkSplitTrackDistanceAtVertexFlag", kinkSplitTrackDistanceAtVertexFlag,
831 "Flag showing if the split kink failed the distance criteria at the kink vertex");
832 REGISTER_VARIABLE("kinkNumberOfReassignedHits", kinkNumberOfReassignedHits,
833 "Number of reassigned hits between kink mother and daughter tracks");
834 REGISTER_VARIABLE("kinkIsSameCharge", kinkIsSameCharge,
835 "Check if charges of mother and daughter tracks are the same");
836
837 // Kink Decay Kinematics
838 REGISTER_VARIABLE("kinkDaughterMomentumInMotherRestFrame", kinkDaughterMomentumInMotherRF,
839 "Kink daughter momentum in mother rest frame with default pair of mass hypotheses,"
840 " set by a user in the decay string");
841 REGISTER_VARIABLE("kinkDaughterCosAngleInMotherRestFrame", kinkDaughterCosThetaInMotherRF,
842 "Kink daughter direction in mother rest frame with respect to mother momentum direction in the lab frame"
843 " with default pair of mass hypotheses, set by a user in the decay string");
844 REGISTER_VARIABLE("kinkDaughterMomentumInMotherRestFrameKPiHypothesis", kinkDaughterMomentumInMotherRFKPi,
845 "Kink daughter momentum in mother rest frame with pion and kaon mass hypotheses");
846 REGISTER_VARIABLE("kinkDaughterCosAngleInMotherRestFrameKPiHypothesis", kinkDaughterCosThetaInMotherRFKPi,
847 "Kink daughter direction in mother rest frame with respect to mother momentum direction in the lab frame"
848 " with pion and kaon mass hypotheses");
849 REGISTER_VARIABLE("kinkDaughterMomentumInMotherRestFrameKMuHypothesis", kinkDaughterMomentumInMotherRFKMu,
850 "Kink daughter momentum in mother rest frame with muon and kaon mass hypotheses");
851 REGISTER_VARIABLE("kinkDaughterCosAngleInMotherRestFrameKMuHypothesis", kinkDaughterCosThetaInMotherRFKMu,
852 "Kink daughter direction in mother rest frame with respect to mother momentum direction in the lab frame"
853 " with muon and kaon mass hypotheses");
854 REGISTER_VARIABLE("kinkDaughterMomentumInMotherRestFramePiMuHypothesis", kinkDaughterMomentumInMotherRFPiMu,
855 "Kink daughter momentum in mother rest frame with muon and pion mass hypotheses");
856 REGISTER_VARIABLE("kinkDaughterCosAngleInMotherRestFramePiMuHypothesis", kinkDaughterCosThetaInMotherRFPiMu,
857 "Kink daughter direction in mother rest frame with respect to mother momentum direction in the lab frame"
858 " with muon and pion mass hypotheses");
859 REGISTER_VARIABLE("kinkDaughterMomentumInMotherRestFrameMuEHypothesis", kinkDaughterMomentumInMotherRFMuE,
860 "Kink daughter momentum in mother rest frame with electron and muon mass hypotheses");
861 REGISTER_VARIABLE("kinkDaughterCosAngleInMotherRestFrameMuEHypothesis", kinkDaughterCosThetaInMotherRFMuE,
862 "Kink daughter direction in mother rest frame with respect to mother momentum direction in the lab frame"
863 " with electron and muon mass hypotheses");
864
865 // Kink Daughter Measured Track Parameters
866 REGISTER_METAVARIABLE("kinkDaughterTrack(variable)", kinkDaughterTrack,
867 "Returns variable for the kink daughter track", Manager::VariableDataType::c_double);
868 REGISTER_METAVARIABLE("kinkDaughterInitTrack(variable)", kinkDaughterInitTrack,
869 "Returns variable for the initial kink daughter track", Manager::VariableDataType::c_double);
870 REGISTER_VARIABLE("kinkDaughterTrackD0AtKinkVertex", kinkDaughterTrackD0AtKinkVertex,
871 "D0 impact parameter of kink daughter track at kink vertex");
872 REGISTER_VARIABLE("kinkDaughterTrackZ0AtKinkVertex", kinkDaughterTrackZ0AtKinkVertex,
873 "Z0 impact parameter of kink daughter track at kink vertex");
874 REGISTER_VARIABLE("kinkDaughterPtAtKinkVertex", kinkDaughterPtAtKinkVertex,
875 "Pt of kink daughter track at kink vertex");
876 REGISTER_VARIABLE("kinkDaughterPzAtKinkVertex", kinkDaughterPzAtKinkVertex,
877 "Pz of kink daughter track at kink vertex");
878 REGISTER_VARIABLE("kinkDaughterPAtKinkVertex", kinkDaughterPAtKinkVertex,
879 "P of kink daughter track at kink vertex");
880
881 // Kink Mother Measured Track Parameters
882 REGISTER_METAVARIABLE("kinkMotherInitTrack(variable)", kinkMotherInitTrack,
883 "Returns variable for the initial kink mother track fit", Manager::VariableDataType::c_double);
884 REGISTER_VARIABLE("kinkMotherTrackD0AtKinkVertex", kinkMotherTrackD0AtKinkVertex,
885 "D0 impact parameter of kink mother track at kink vertex");
886 REGISTER_VARIABLE("kinkMotherTrackZ0AtKinkVertex", kinkMotherTrackZ0AtKinkVertex,
887 "Z0 impact parameter of kink mother track at kink vertex");
888 REGISTER_VARIABLE("kinkMotherPtAtKinkVertex", kinkMotherPtAtKinkVertex,
889 "Pt of kink mother track at kink vertex");
890 REGISTER_VARIABLE("kinkMotherPzAtKinkVertex", kinkMotherPzAtKinkVertex,
891 "Pz of kink mother track at kink vertex");
892 REGISTER_VARIABLE("kinkMotherPAtKinkVertex", kinkMotherPAtKinkVertex,
893 "P of kink mother track at kink vertex");
894
895 //-----------------------------------------------
896 // Check if the particle is a part of Kink
897 REGISTER_VARIABLE("particleIsInKink", particleIsInKink,
898 "Particle is used in a Kink object");
899 REGISTER_VARIABLE("particleIsMotherInKink", particleIsMotherInKink,
900 "Particle is a mother in a Kink object");
901 REGISTER_VARIABLE("particleIsDaughterInKink", particleIsDaughterInKink,
902 "Particle is a daughter in a Kink object");
903 REGISTER_VARIABLE("particleIsSplitKink", particleIsSplitKink,
904 "Particle is a split track in a Kink object");
905
906 //-----------------------------------------------
907 // MC VARIABLES
908
909 // Kink from track pair MC variables
910 REGISTER_VARIABLE("kinkIsMCRelated", kinkPairIsMCRelated,
911 "Mother and daughter tracks have MC relations and two of them are related");
912 REGISTER_VARIABLE("kinkIsClone", kinkPairIsClone,
913 "Mother and daughter tracks have the same MCParticle relation");
914 REGISTER_VARIABLE("kinkIsReal", kinkPairIsReal,
915 "Mother and daughter tracks are mother/daughter related");
916 REGISTER_VARIABLE("kinkIsDecayInFlight", kinkPairIsDecayInFlight,
917 "Kink is a decay-in-flight");
918 REGISTER_VARIABLE("kinkIsHadronScattering", kinkPairIsHadronScattering,
919 "Kink is a hadron scattering");
920
921 // Kink from track pair MC kinematics
922 REGISTER_VARIABLE("kinkDaughterMCMomentumInMotherRestFrame", kinkDaughterMomentumInMotherRFMC,
923 "Kink daughter momentum in mother rest frame MC (works only for kink created from a track pair)\n"
924 "Makes sense only for real decays-in-flight; however, it is not checked here");
925 REGISTER_VARIABLE("kinkDaughterMCCosAngleInMotherRestFrame", kinkDaughterCosThetaInMotherRFMC,
926 "Kink daughter direction in mother rest frame with respect to mother momentum direction in the lab frame"
927 " MC (works only for kink created from a track pair)\n"
928 "Makes sense only for real decays-in-flight; however, it is not checked here");
929
930 // Kink from track pair daughter MC variables
931 REGISTER_METAVARIABLE("kinkPairDaughterMC(variable)", kinkPairDaughterMC,
932 "Returns MC variable for the kink daughter for kinks created from two separate tracks", Manager::VariableDataType::c_double);
933
934 // Kink from track pair mother MC variables
935 REGISTER_METAVARIABLE("kinkPairMotherMC(variable)", kinkPairMotherMC,
936 "Returns MC variable for the kink mother for kinks created from two separate tracks", Manager::VariableDataType::c_double);
937 REGISTER_VARIABLE("kinkMotherMCPXAtDV", kinkPairMotherMCPXAtDecayVertex,
938 "Generated PX of the kink mother at the decay vertex for kinks created from two separate tracks");
939 REGISTER_VARIABLE("kinkMotherMCPYAtDV", kinkPairMotherMCPYAtDecayVertex,
940 "Generated PY of the kink mother at the decay vertex for kinks created from two separate tracks");
941 REGISTER_VARIABLE("kinkMotherMCPZAtDV", kinkPairMotherMCPZAtDecayVertex,
942 "Generated PZ of the kink mother at the decay vertex for kinks created from two separate tracks");
943 REGISTER_VARIABLE("kinkMotherMCPTAtDV", kinkPairMotherMCPTAtDecayVertex,
944 "Generated PT of the kink mother at the decay vertex for kinks created from two separate tracks");
945 REGISTER_VARIABLE("kinkMotherMCPAtDV", kinkPairMotherMCPAtDecayVertex,
946 "Generated P of the kink mother at the decay vertex for kinks created from two separate tracks");
947 REGISTER_VARIABLE("kinkMotherMCEAtDV", kinkPairMotherMCEAtDecayVertex,
948 "Generated E of the kink mother at the decay vertex for kinks created from two separate tracks");
949 }
951}
Belle2::BFieldManager::getFieldInTesla
static ROOT::Math::XYZVector getFieldInTesla(const ROOT::Math::XYZVector &pos)
return the magnetic field at a given position in Tesla.
Definition: BFieldManager.h:60
Belle2::Const::muon
static const ChargedStable muon
muon particle
Definition: Const.h:660
Belle2::Const::pion
static const ChargedStable pion
charged pion particle
Definition: Const.h:661
Belle2::Const::doubleNaN
static const double doubleNaN
quiet_NaN
Definition: Const.h:703
Belle2::Const::kaon
static const ChargedStable kaon
charged kaon particle
Definition: Const.h:662
Belle2::Const::electron
static const ChargedStable electron
electron particle
Definition: Const.h:659
Belle2::RelationsInterface::getRelated
T * getRelated(const std::string &name="", const std::string &namedRelation="") const
Get the object to or from which this object has a relation.
Definition: RelationsObject.h:278
Belle2::Variable::Manager::FunctionPtr
std::function< VarVariant(const Particle *)> FunctionPtr
functions stored take a const Particle* and return VarVariant.
Definition: Manager.h:112
Belle2::Variable::Manager::getVariable
const Var * getVariable(std::string name)
Get the variable belonging to the given key.
Definition: Manager.cc:58
Belle2::Variable::Manager::Instance
static Manager & Instance()
get singleton instance.
Definition: Manager.cc:26
Belle2::sqrt
double sqrt(double a)
sqrt for double
Definition: beamHelpers.h:28
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:17
std
STL namespace.