Belle II Software development
MetaVariables.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/MetaVariables.h>
11#include <analysis/variables/MCTruthVariables.h>
12
13#include <analysis/VariableManager/Utility.h>
14#include <analysis/dataobjects/Particle.h>
15#include <analysis/dataobjects/ParticleList.h>
16#include <analysis/dataobjects/EventKinematics.h>
17#include <analysis/utility/PCmsLabTransform.h>
18#include <analysis/utility/ReferenceFrame.h>
19#include <analysis/utility/EvtPDLUtil.h>
20#include <analysis/utility/ParticleCopy.h>
21#include <analysis/utility/ValueIndexPairSorting.h>
22#include <analysis/ClusterUtility/ClusterUtils.h>
23#include <analysis/variables/VariableFormulaConstructor.h>
24
25#include <framework/logging/Logger.h>
26#include <framework/datastore/StoreArray.h>
27#include <framework/datastore/StoreObjPtr.h>
28#include <framework/dataobjects/EventExtraInfo.h>
29#include <framework/utilities/Conversion.h>
30#include <framework/utilities/MakeROOTCompatible.h>
31#include <framework/gearbox/Const.h>
32
33#include <mdst/dataobjects/Track.h>
34#include <mdst/dataobjects/MCParticle.h>
35#include <mdst/dataobjects/ECLCluster.h>
36#include <mdst/dataobjects/TrackFitResult.h>
37
38#include <boost/algorithm/string.hpp>
39#include <limits>
40
41#include <cmath>
42#include <stdexcept>
43#include <memory>
44#include <regex>
45
46#include <TDatabasePDG.h>
47#include <Math/Vector4D.h>
48#include <Math/VectorUtil.h>
49
50namespace Belle2 {
55 namespace Variable {
56
57 Manager::FunctionPtr useRestFrame(const std::vector<std::string>& arguments)
58 {
59 if (arguments.size() == 1) {
60 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
61 auto func = [var](const Particle * particle) -> double {
62 UseReferenceFrame<RestFrame> frame(particle);
63 double result = std::get<double>(var->function(particle));
64 return result;
65 };
66 return func;
67 } else {
68 B2FATAL("Wrong number of arguments for meta function useRestFrame");
69 }
70 }
71
72 Manager::FunctionPtr useCMSFrame(const std::vector<std::string>& arguments)
73 {
74 if (arguments.size() == 1) {
75 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
76 auto func = [var](const Particle * particle) -> double {
77 UseReferenceFrame<CMSFrame> frame;
78 double result = std::get<double>(var->function(particle));
79 return result;
80 };
81 return func;
82 } else {
83 B2FATAL("Wrong number of arguments for meta function useCMSFrame");
84 }
85 }
86
87 Manager::FunctionPtr useLabFrame(const std::vector<std::string>& arguments)
88 {
89 if (arguments.size() == 1) {
90 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
91 auto func = [var](const Particle * particle) -> double {
92 UseReferenceFrame<LabFrame> frame;
93 double result = std::get<double>(var->function(particle));
94 return result;
95 };
96 return func;
97 } else {
98 B2FATAL("Wrong number of arguments for meta function useLabFrame");
99 }
100 }
101
102 Manager::FunctionPtr useTagSideRecoilRestFrame(const std::vector<std::string>& arguments)
103 {
104 if (arguments.size() == 2) {
105 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
106
107 int daughterIndexTagB = 0;
108 try {
109 daughterIndexTagB = convertString<int>(arguments[1]);
110 } catch (std::invalid_argument&) {
111 B2FATAL("Second argument of useTagSideRecoilRestFrame meta function must be integer!");
112 }
113
114 auto func = [var, daughterIndexTagB](const Particle * particle) -> double {
115 if (particle->getPDGCode() != 300553)
116 {
117 B2ERROR("Variable should only be used on a Upsilon(4S) Particle List!");
118 return Const::doubleNaN;
119 }
120
121 PCmsLabTransform T;
122 ROOT::Math::PxPyPzEVector pSigB = T.getBeamFourMomentum() - particle->getDaughter(daughterIndexTagB)->get4Vector();
123 Particle tmp(pSigB, -particle->getDaughter(daughterIndexTagB)->getPDGCode());
124
125 UseReferenceFrame<RestFrame> frame(&tmp);
126 double result = std::get<double>(var->function(particle));
127 return result;
128 };
129
130 return func;
131 } else {
132 B2FATAL("Wrong number of arguments for meta function useTagSideRecoilRestFrame");
133 }
134 }
135
136 Manager::FunctionPtr useParticleRestFrame(const std::vector<std::string>& arguments)
137 {
138 if (arguments.size() == 2) {
139 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
140 std::string listName = arguments[1];
141 auto func = [var, listName](const Particle * particle) -> double {
142 StoreObjPtr<ParticleList> list(listName);
143 unsigned listSize = list->getListSize();
144 if (listSize == 0)
145 return Const::doubleNaN;
146 if (listSize > 1)
147 B2WARNING("The selected ParticleList contains more than 1 Particles in this event. The variable useParticleRestFrame will use only the first candidate, and the result may not be the expected one."
148 << LogVar("ParticleList", listName)
149 << LogVar("Number of candidates in the list", listSize));
150 const Particle* p = list->getParticle(0);
151 UseReferenceFrame<RestFrame> frame(p);
152 double result = std::get<double>(var->function(particle));
153 return result;
154 };
155 return func;
156 } else {
157 B2FATAL("Wrong number of arguments for meta function useParticleRestFrame.");
158 }
159 }
160
161 Manager::FunctionPtr useRecoilParticleRestFrame(const std::vector<std::string>& arguments)
162 {
163 if (arguments.size() == 2) {
164 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
165 std::string listName = arguments[1];
166 auto func = [var, listName](const Particle * particle) -> double {
167 StoreObjPtr<ParticleList> list(listName);
168 unsigned listSize = list->getListSize();
169 if (listSize == 0)
170 return Const::doubleNaN;
171 if (listSize > 1)
172 B2WARNING("The selected ParticleList contains more than 1 Particles in this event. The variable useParticleRestFrame will use only the first candidate, and the result may not be the expected one."
173 << LogVar("ParticleList", listName)
174 << LogVar("Number of candidates in the list", listSize));
175 const Particle* p = list->getParticle(0);
176 PCmsLabTransform T;
177 ROOT::Math::PxPyPzEVector recoil = T.getBeamFourMomentum() - p->get4Vector();
178 /* Let's use 0 as PDG code to avoid wrong assumptions. */
179 Particle pRecoil(recoil, 0);
180 pRecoil.setVertex(particle->getVertex());
181 UseReferenceFrame<RestFrame> frame(&pRecoil);
182 double result = std::get<double>(var->function(particle));
183 return result;
184 };
185 return func;
186 } else {
187 B2FATAL("Wrong number of arguments for meta function useParticleRestFrame.");
188 }
189 }
190
191 Manager::FunctionPtr useDaughterRestFrame(const std::vector<std::string>& arguments)
192 {
193 if (arguments.size() >= 2) {
194 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
195 auto func = [var, arguments](const Particle * particle) -> double {
196
197 // Sum of the 4-momenta of all the selected daughters
198 ROOT::Math::PxPyPzEVector pSum(0, 0, 0, 0);
199
200 for (unsigned int i = 1; i < arguments.size(); i++)
201 {
202 auto generalizedIndex = arguments[i];
203 const Particle* dauPart = particle->getParticleFromGeneralizedIndexString(generalizedIndex);
204 if (dauPart)
205 pSum += dauPart->get4Vector();
206 else
207 return Const::doubleNaN;
208 }
209 Particle tmp(pSum, 0);
210 UseReferenceFrame<RestFrame> frame(&tmp);
211 double result = std::get<double>(var->function(particle));
212 return result;
213 };
214 return func;
215 } else {
216 B2FATAL("Wrong number of arguments for meta function useDaughterRestFrame.");
217 }
218 }
219
220 Manager::FunctionPtr useDaughterRecoilRestFrame(const std::vector<std::string>& arguments)
221 {
222 if (arguments.size() >= 2) {
223 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
224 auto func = [var, arguments](const Particle * particle) -> double {
225
226 // Sum of the 4-momenta of all the selected daughters
227 ROOT::Math::PxPyPzEVector pSum(0, 0, 0, 0);
228
229 for (unsigned int i = 1; i < arguments.size(); i++)
230 {
231 auto generalizedIndex = arguments[i];
232 const Particle* dauPart = particle->getParticleFromGeneralizedIndexString(generalizedIndex);
233 if (dauPart)
234 pSum += dauPart->get4Vector();
235 else
236 return Const::doubleNaN;
237 }
238 PCmsLabTransform T;
239 ROOT::Math::PxPyPzEVector recoil = T.getBeamFourMomentum() - pSum;
240 /* Let's use 0 as PDG code to avoid wrong assumptions. */
241 Particle pRecoil(recoil, 0);
242 UseReferenceFrame<RestFrame> frame(&pRecoil);
243 double result = std::get<double>(var->function(particle));
244 return result;
245 };
246 return func;
247 } else {
248 B2FATAL("Wrong number of arguments for meta function useDaughterRecoilRestFrame.");
249 }
250 }
251
252 Manager::FunctionPtr useMCancestorBRestFrame(const std::vector<std::string>& arguments)
253 {
254 if (arguments.size() == 1) {
255 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
256 auto func = [var](const Particle * particle) -> double {
257 int index = ancestorBIndex(particle);
258 if (index < 0) return Const::doubleNaN;
259 StoreArray<MCParticle> mcparticles;
260 Particle temp(mcparticles[index]);
261 UseReferenceFrame<RestFrame> frame(&temp);
262 double result = std::get<double>(var->function(particle));
263 return result;
264 };
265 return func;
266 } else {
267 B2FATAL("Wrong number of arguments for meta function useMCancestorBRestFrame.");
268 }
269 }
270
271 Manager::FunctionPtr extraInfo(const std::vector<std::string>& arguments)
272 {
273 if (arguments.size() == 1) {
274 auto extraInfoName = arguments[0];
275 auto func = [extraInfoName](const Particle * particle) -> double {
276 if (particle == nullptr)
277 {
278 B2WARNING("Returns NaN because the particle is nullptr! If you want EventExtraInfo variables, please use eventExtraInfo() instead");
279 return Const::doubleNaN;
280 }
281 if (particle->hasExtraInfo(extraInfoName))
282 {
283 return particle->getExtraInfo(extraInfoName);
284 } else
285 {
286 return Const::doubleNaN;
287 }
288 };
289 return func;
290 } else {
291 B2FATAL("Wrong number of arguments for meta function extraInfo");
292 }
293 }
294
295 Manager::FunctionPtr eventExtraInfo(const std::vector<std::string>& arguments)
296 {
297 if (arguments.size() == 1) {
298 auto extraInfoName = arguments[0];
299 auto func = [extraInfoName](const Particle*) -> double {
300 StoreObjPtr<EventExtraInfo> eventExtraInfo;
301 if (not eventExtraInfo.isValid())
302 return Const::doubleNaN;
303 if (eventExtraInfo->hasExtraInfo(extraInfoName))
304 {
305 return eventExtraInfo->getExtraInfo(extraInfoName);
306 } else
307 {
308 return Const::doubleNaN;
309 }
310 };
311 return func;
312 } else {
313 B2FATAL("Wrong number of arguments for meta function extraInfo");
314 }
315 }
316
317 Manager::FunctionPtr eventCached(const std::vector<std::string>& arguments)
318 {
319 if (arguments.size() == 1) {
320 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
321 std::string key = std::string("__") + MakeROOTCompatible::makeROOTCompatible(var->name);
322 auto func = [var, key](const Particle*) -> double {
323
324 StoreObjPtr<EventExtraInfo> eventExtraInfo;
325 if (not eventExtraInfo.isValid())
326 eventExtraInfo.create();
327 if (eventExtraInfo->hasExtraInfo(key))
328 {
329 return eventExtraInfo->getExtraInfo(key);
330 } else
331 {
332 double value = Const::doubleNaN;
333 auto var_result = var->function(nullptr);
334 if (std::holds_alternative<double>(var_result)) {
335 value = std::get<double>(var_result);
336 } else if (std::holds_alternative<int>(var_result)) {
337 return std::get<int>(var_result);
338 } else if (std::holds_alternative<bool>(var_result)) {
339 return std::get<bool>(var_result);
340 }
341 eventExtraInfo->addExtraInfo(key, value);
342 return value;
343 }
344 };
345 return func;
346 } else {
347 B2FATAL("Wrong number of arguments for meta function eventCached");
348 }
349 }
350
351 Manager::FunctionPtr particleCached(const std::vector<std::string>& arguments)
352 {
353 if (arguments.size() == 1) {
354 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
355 std::string key = std::string("__") + MakeROOTCompatible::makeROOTCompatible(var->name);
356 auto func = [var, key](const Particle * particle) -> double {
357
358 if (particle->hasExtraInfo(key))
359 {
360 return particle->getExtraInfo(key);
361 } else
362 {
363 double value = std::get<double>(var->function(particle));
364 // Remove constness from Particle pointer.
365 // The extra-info is used as a cache in our case,
366 // indicated by the double-underscore in front of the key.
367 // One could implement the cache as a separate property of the particle object
368 // and mark it as mutable, however, this would only lead to code duplication
369 // and an increased size of the particle object.
370 // Thus, we decided to use the extra-info field and cast away the const in this case.
371 const_cast<Particle*>(particle)->addExtraInfo(key, value);
372 return value;
373 }
374 };
375 return func;
376 } else {
377 B2FATAL("Wrong number of arguments for meta function particleCached");
378 }
379 }
380
381 // Formula of other variables, going to require a space between all operators and operations.
382 // Later can add some check for : (colon) trailing + or - to distinguish between particle lists
383 // and operations, but for now cbf.
384 Manager::FunctionPtr formula(const std::vector<std::string>& arguments)
385 {
386 if (arguments.size() != 1) B2FATAL("Wrong number of arguments for meta function formula");
387 FormulaParser<VariableFormulaConstructor> parser;
388 try {
389 return parser.parse(arguments[0]);
390 } catch (std::runtime_error& e) {
391 B2FATAL(e.what());
392 }
393 }
394
395 Manager::FunctionPtr nCleanedTracks(const std::vector<std::string>& arguments)
396 {
397 if (arguments.size() <= 1) {
398
399 std::string cutString;
400 if (arguments.size() == 1)
401 cutString = arguments[0];
402 std::shared_ptr<Variable::Cut> cut = std::shared_ptr<Variable::Cut>(Variable::Cut::compile(cutString));
403 auto func = [cut](const Particle*) -> int {
404
405 int number_of_tracks = 0;
406 StoreArray<Track> tracks;
407 for (const auto& track : tracks)
408 {
409 const TrackFitResult* trackFit = track.getTrackFitResultWithClosestMass(Const::pion);
410 if (!trackFit) continue;
411 if (trackFit->getChargeSign() == 0) {
412 // Ignore track
413 } else {
414 Particle particle(&track, Const::pion);
415 if (cut->check(&particle))
416 number_of_tracks++;
417 }
418 }
419
420 return number_of_tracks;
421
422 };
423 return func;
424 } else {
425 B2FATAL("Wrong number of arguments for meta function nCleanedTracks");
426 }
427 }
428
429 Manager::FunctionPtr nCleanedECLClusters(const std::vector<std::string>& arguments)
430 {
431 if (arguments.size() <= 1) {
432
433 std::string cutString;
434 if (arguments.size() == 1)
435 cutString = arguments[0];
436 std::shared_ptr<Variable::Cut> cut = std::shared_ptr<Variable::Cut>(Variable::Cut::compile(cutString));
437 auto func = [cut](const Particle*) -> int {
438
439 int number_of_clusters = 0;
440 StoreArray<ECLCluster> clusters;
441 for (const auto& cluster : clusters)
442 {
443 // look only at momentum of N1 (n photons) ECLClusters
444 if (!cluster.hasHypothesis(ECLCluster::EHypothesisBit::c_nPhotons))
445 continue;
446
447 Particle particle(&cluster);
448 if (cut->check(&particle))
449 number_of_clusters++;
450 }
451
452 return number_of_clusters;
453
454 };
455 return func;
456 } else {
457 B2FATAL("Wrong number of arguments for meta function nCleanedECLClusters");
458 }
459 }
460
461 Manager::FunctionPtr passesCut(const std::vector<std::string>& arguments)
462 {
463 if (arguments.size() == 1) {
464 std::string cutString = arguments[0];
465 std::shared_ptr<Variable::Cut> cut = std::shared_ptr<Variable::Cut>(Variable::Cut::compile(cutString));
466 auto func = [cut](const Particle * particle) -> bool {
467 if (cut->check(particle))
468 return 1;
469 else
470 return 0;
471 };
472 return func;
473 } else {
474 B2FATAL("Wrong number of arguments for meta function passesCut");
475 }
476 }
477
478 Manager::FunctionPtr passesEventCut(const std::vector<std::string>& arguments)
479 {
480 if (arguments.size() == 1) {
481 std::string cutString = arguments[0];
482 std::shared_ptr<Variable::Cut> cut = std::shared_ptr<Variable::Cut>(Variable::Cut::compile(cutString));
483 auto func = [cut](const Particle*) -> bool {
484 if (cut->check(nullptr))
485 return 1;
486 else
487 return 0;
488 };
489 return func;
490 } else {
491 B2FATAL("Wrong number of arguments for meta function passesEventCut");
492 }
493 }
494
495 Manager::FunctionPtr varFor(const std::vector<std::string>& arguments)
496 {
497 if (arguments.size() == 2) {
498 int pdgCode = 0;
499 try {
500 pdgCode = convertString<int>(arguments[0]);
501 } catch (std::invalid_argument&) {
502 B2FATAL("The first argument of varFor meta function must be a positive integer!");
503 }
504 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[1]);
505 auto func = [pdgCode, var](const Particle * particle) -> double {
506 if (std::abs(particle->getPDGCode()) == std::abs(pdgCode))
507 {
508 auto var_result = var->function(particle);
509 if (std::holds_alternative<double>(var_result)) {
510 return std::get<double>(var_result);
511 } else if (std::holds_alternative<int>(var_result)) {
512 return std::get<int>(var_result);
513 } else if (std::holds_alternative<bool>(var_result)) {
514 return std::get<bool>(var_result);
515 } else return Const::doubleNaN;
516 } else return Const::doubleNaN;
517 };
518 return func;
519 } else {
520 B2FATAL("Wrong number of arguments for meta function varFor");
521 }
522 }
523
524 Manager::FunctionPtr varForMCGen(const std::vector<std::string>& arguments)
525 {
526 if (arguments.size() == 1) {
527 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
528 auto func = [var](const Particle * particle) -> double {
529 if (particle->getMCParticle())
530 {
531 if (particle->getMCParticle()->getStatus(MCParticle::c_PrimaryParticle)
532 && (! particle->getMCParticle()->getStatus(MCParticle::c_IsVirtual))
533 && (! particle->getMCParticle()->getStatus(MCParticle::c_Initial))) {
534 auto var_result = var->function(particle);
535 if (std::holds_alternative<double>(var_result)) {
536 return std::get<double>(var_result);
537 } else if (std::holds_alternative<int>(var_result)) {
538 return std::get<int>(var_result);
539 } else if (std::holds_alternative<bool>(var_result)) {
540 return std::get<bool>(var_result);
541 } else return Const::doubleNaN;
542 } else return Const::doubleNaN;
543 } else return Const::doubleNaN;
544 };
545 return func;
546 } else {
547 B2FATAL("Wrong number of arguments for meta function varForMCGen");
548 }
549 }
550
551 Manager::FunctionPtr nParticlesInList(const std::vector<std::string>& arguments)
552 {
553 if (arguments.size() == 1) {
554 std::string listName = arguments[0];
555 auto func = [listName](const Particle * particle) -> int {
556
557 (void) particle;
558 StoreObjPtr<ParticleList> listOfParticles(listName);
559
560 if (!(listOfParticles.isValid())) B2FATAL("Invalid Listname " << listName << " given to nParticlesInList");
561
562 return listOfParticles->getListSize();
563
564 };
565 return func;
566 } else {
567 B2FATAL("Wrong number of arguments for meta function nParticlesInList");
568 }
569 }
570
571 Manager::FunctionPtr isInList(const std::vector<std::string>& arguments)
572 {
573 // unpack arguments, there should be only one: the name of the list we're checking
574 if (arguments.size() != 1) {
575 B2FATAL("Wrong number of arguments for isInList");
576 }
577 auto listName = arguments[0];
578
579 auto func = [listName](const Particle * particle) -> bool {
580
581 // check the list exists
582 StoreObjPtr<ParticleList> list(listName);
583 if (!(list.isValid()))
584 {
585 B2FATAL("Invalid Listname " << listName << " given to isInList");
586 }
587
588 // is the particle in the list?
589 return list->contains(particle);
590
591 };
592 return func;
593 }
594
595 Manager::FunctionPtr sourceObjectIsInList(const std::vector<std::string>& arguments)
596 {
597 // unpack arguments, there should be only one: the name of the list we're checking
598 if (arguments.size() != 1) {
599 B2FATAL("Wrong number of arguments for sourceObjectIsInList");
600 }
601 auto listName = arguments[0];
602
603 auto func = [listName](const Particle * particle) -> int {
604
605 // check the list exists
606 StoreObjPtr<ParticleList> list(listName);
607 if (!(list.isValid()))
608 {
609 B2FATAL("Invalid Listname " << listName << " given to sourceObjectIsInList");
610 }
611
612 // this only makes sense for particles that are *not* composite and come
613 // from some mdst object (tracks, clusters..)
614 Particle::EParticleSourceObject particlesource = particle->getParticleSource();
615 if (particlesource == Particle::EParticleSourceObject::c_Composite
616 or particlesource == Particle::EParticleSourceObject::c_Undefined)
617 return -1;
618
619 // it *is* possible to have a particle list from different sources (like
620 // hadrons from the ECL and KLM) so we have to check each particle in
621 // the list individually
622 for (unsigned i = 0; i < list->getListSize(); ++i)
623 {
624 Particle* iparticle = list->getParticle(i);
625 if (particle->getMdstSource() == iparticle->getMdstSource())
626 return 1;
627 }
628 return 0;
629
630 };
631 return func;
632 }
633
634 Manager::FunctionPtr mcParticleIsInMCList(const std::vector<std::string>& arguments)
635 {
636 // unpack arguments, there should be only one: the name of the list we're checking
637 if (arguments.size() != 1) {
638 B2FATAL("Wrong number of arguments for mcParticleIsInMCList");
639 }
640 auto listName = arguments[0];
641
642 auto func = [listName](const Particle * particle) -> bool {
643
644 // check the list exists
645 StoreObjPtr<ParticleList> list(listName);
646 if (!(list.isValid()))
647 B2FATAL("Invalid Listname " << listName << " given to mcParticleIsInMCList");
648
649 // this can only be true for mc-matched particles or particles are created from MCParticles
650 const MCParticle* mcp = particle->getMCParticle();
651 if (mcp == nullptr) return false;
652
653 // check every particle in the input list is not matched to (or created from) the same MCParticle
654 for (unsigned i = 0; i < list->getListSize(); ++i)
655 {
656 const MCParticle* imcp = list->getParticle(i)->getMCParticle();
657 if ((imcp != nullptr) and (mcp->getArrayIndex() == imcp->getArrayIndex()))
658 return true;
659 }
660 return false;
661 };
662 return func;
663 }
664
665 Manager::FunctionPtr isDaughterOfList(const std::vector<std::string>& arguments)
666 {
667 B2WARNING("isDaughterOfList is outdated and replaced by isDescendantOfList.");
668 std::vector<std::string> new_arguments = arguments;
669 new_arguments.push_back(std::string("1"));
670 return isDescendantOfList(new_arguments);
671 }
672
673 Manager::FunctionPtr isGrandDaughterOfList(const std::vector<std::string>& arguments)
674 {
675 B2WARNING("isGrandDaughterOfList is outdated and replaced by isDescendantOfList.");
676 std::vector<std::string> new_arguments = arguments;
677 new_arguments.push_back(std::string("2"));
678 return isDescendantOfList(new_arguments);
679 }
680
681 Manager::FunctionPtr isDescendantOfList(const std::vector<std::string>& arguments)
682 {
683 if (arguments.size() > 0) {
684 auto listNames = arguments;
685 auto func = [listNames](const Particle * particle) -> bool {
686 bool output = false;
687 int generation_flag = -1;
688 try
689 {
690 generation_flag = convertString<int>(listNames.back());
691 } catch (std::exception& e) {}
692
693 for (auto& iListName : listNames)
694 {
695 try {
696 convertString<int>(iListName);
697 continue;
698 } catch (std::exception& e) {}
699
700 // Creating recursive lambda
701 auto list_comparison = [](auto&& self, const Particle * m, const Particle * p, int flag)-> bool {
702 bool result = false;
703 for (unsigned i = 0; i < m->getNDaughters(); ++i)
704 {
705 const Particle* daughter = m->getDaughter(i);
706 if ((flag == 1.) or (flag < 0)) {
707 if (p->isCopyOf(daughter)) {
708 return true;
709 }
710 }
711
712 if (flag != 1.) {
713 if (daughter->getNDaughters() > 0) {
714 result = self(self, daughter, p, flag - 1);
715 if (result == 1) {
716 return true;
717 }
718 }
719 }
720 }
721 return result;
722 };
723
724 StoreObjPtr<ParticleList> listOfParticles(iListName);
725
726 if (!(listOfParticles.isValid())) B2FATAL("Invalid Listname " << iListName << " given to isDescendantOfList");
727
728 for (unsigned i = 0; i < listOfParticles->getListSize(); ++i) {
729 Particle* iParticle = listOfParticles->getParticle(i);
730 output = list_comparison(list_comparison, iParticle, particle, generation_flag);
731 if (output) {
732 return output;
733 }
734 }
735 }
736 return output;
737 };
738 return func;
739 } else {
740 B2FATAL("Wrong number of arguments for meta function isDescendantOfList");
741 }
742 }
743
744 Manager::FunctionPtr isMCDescendantOfList(const std::vector<std::string>& arguments)
745 {
746 if (arguments.size() > 0) {
747 auto listNames = arguments;
748 auto func = [listNames](const Particle * particle) -> bool {
749 bool output = false;
750 int generation_flag = -1;
751 try
752 {
753 generation_flag = convertString<int>(listNames.back());
754 } catch (std::exception& e) {}
755
756 if (particle->getMCParticle() == nullptr)
757 {
758 return false;
759 }
760
761 for (auto& iListName : listNames)
762 {
763 try {
764 std::stod(iListName);
765 continue;
766 } catch (std::exception& e) {}
767 // Creating recursive lambda
768 auto list_comparison = [](auto&& self, const Particle * m, const Particle * p, int flag)-> bool {
769 bool result = false;
770 for (unsigned i = 0; i < m->getNDaughters(); ++i)
771 {
772 const Particle* daughter = m->getDaughter(i);
773 if ((flag == 1.) or (flag < 0)) {
774 if (daughter->getMCParticle() != nullptr) {
775 if (p->getMCParticle()->getArrayIndex() == daughter->getMCParticle()->getArrayIndex()) {
776 return true;
777 }
778 }
779 }
780 if (flag != 1.) {
781 if (daughter->getNDaughters() > 0) {
782 result = self(self, daughter, p, flag - 1);
783 if (result) {
784 return true;
785 }
786 }
787 }
788 }
789 return result;
790 };
791
792 StoreObjPtr<ParticleList> listOfParticles(iListName);
793
794 if (!(listOfParticles.isValid())) B2FATAL("Invalid Listname " << iListName << " given to isMCDescendantOfList");
795
796 for (unsigned i = 0; i < listOfParticles->getListSize(); ++i) {
797 Particle* iParticle = listOfParticles->getParticle(i);
798 output = list_comparison(list_comparison, iParticle, particle, generation_flag);
799 if (output) {
800 return output;
801 }
802 }
803 }
804 return output;
805 };
806 return func;
807 } else {
808 B2FATAL("Wrong number of arguments for meta function isMCDescendantOfList");
809 }
810 }
811
812 Manager::FunctionPtr daughterProductOf(const std::vector<std::string>& arguments)
813 {
814 if (arguments.size() == 1) {
815 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
816 auto func = [var](const Particle * particle) -> double {
817 double product = 1.0;
818 if (particle->getNDaughters() == 0)
819 {
820 return Const::doubleNaN;
821 }
822 if (std::holds_alternative<double>(var->function(particle->getDaughter(0))))
823 {
824 for (unsigned j = 0; j < particle->getNDaughters(); ++j) {
825 product *= std::get<double>(var->function(particle->getDaughter(j)));
826 }
827 } else if (std::holds_alternative<int>(var->function(particle->getDaughter(0))))
828 {
829 for (unsigned j = 0; j < particle->getNDaughters(); ++j) {
830 product *= std::get<int>(var->function(particle->getDaughter(j)));
831 }
832 } else return Const::doubleNaN;
833 return product;
834 };
835 return func;
836 } else {
837 B2FATAL("Wrong number of arguments for meta function daughterProductOf");
838 }
839 }
840
841 Manager::FunctionPtr daughterSumOf(const std::vector<std::string>& arguments)
842 {
843 if (arguments.size() == 1) {
844 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
845 auto func = [var](const Particle * particle) -> double {
846 double sum = 0.0;
847 if (particle->getNDaughters() == 0)
848 {
849 return Const::doubleNaN;
850 }
851 if (std::holds_alternative<double>(var->function(particle->getDaughter(0))))
852 {
853 for (unsigned j = 0; j < particle->getNDaughters(); ++j) {
854 sum += std::get<double>(var->function(particle->getDaughter(j)));
855 }
856 } else if (std::holds_alternative<int>(var->function(particle->getDaughter(0))))
857 {
858 for (unsigned j = 0; j < particle->getNDaughters(); ++j) {
859 sum += std::get<int>(var->function(particle->getDaughter(j)));
860 }
861 } else return Const::doubleNaN;
862 return sum;
863 };
864 return func;
865 } else {
866 B2FATAL("Wrong number of arguments for meta function daughterSumOf");
867 }
868 }
869
870 Manager::FunctionPtr daughterLowest(const std::vector<std::string>& arguments)
871 {
872 if (arguments.size() == 1) {
873 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
874 auto func = [var](const Particle * particle) -> double {
875 double min = Const::doubleNaN;
876 if (particle->getNDaughters() == 0)
877 {
878 return Const::doubleNaN;
879 }
880 if (std::holds_alternative<double>(var->function(particle->getDaughter(0))))
881 {
882 for (unsigned j = 0; j < particle->getNDaughters(); ++j) {
883 double iValue = std::get<double>(var->function(particle->getDaughter(j)));
884 if (std::isnan(iValue)) continue;
885 if (std::isnan(min)) min = iValue;
886 if (iValue < min) min = iValue;
887 }
888 } else if (std::holds_alternative<int>(var->function(particle->getDaughter(0))))
889 {
890 for (unsigned j = 0; j < particle->getNDaughters(); ++j) {
891 int iValue = std::get<int>(var->function(particle->getDaughter(j)));
892 if (std::isnan(min)) min = iValue;
893 if (iValue < min) min = iValue;
894 }
895 }
896 return min;
897 };
898 return func;
899 } else {
900 B2FATAL("Wrong number of arguments for meta function daughterLowest");
901 }
902 }
903
904 Manager::FunctionPtr daughterHighest(const std::vector<std::string>& arguments)
905 {
906 if (arguments.size() == 1) {
907 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
908 auto func = [var](const Particle * particle) -> double {
909 double max = Const::doubleNaN;
910 if (particle->getNDaughters() == 0)
911 {
912 return Const::doubleNaN;
913 }
914 if (std::holds_alternative<double>(var->function(particle->getDaughter(0))))
915 {
916 for (unsigned j = 0; j < particle->getNDaughters(); ++j) {
917 double iValue = std::get<double>(var->function(particle->getDaughter(j)));
918 if (std::isnan(iValue)) continue;
919 if (std::isnan(max)) max = iValue;
920 if (iValue > max) max = iValue;
921 }
922 } else if (std::holds_alternative<int>(var->function(particle->getDaughter(0))))
923 {
924 for (unsigned j = 0; j < particle->getNDaughters(); ++j) {
925 int iValue = std::get<int>(var->function(particle->getDaughter(j)));
926 if (std::isnan(max)) max = iValue;
927 if (iValue > max) max = iValue;
928 }
929 }
930 return max;
931 };
932 return func;
933 } else {
934 B2FATAL("Wrong number of arguments for meta function daughterHighest");
935 }
936 }
937
938 Manager::FunctionPtr daughterDiffOf(const std::vector<std::string>& arguments)
939 {
940 if (arguments.size() == 3) {
941 auto func = [arguments](const Particle * particle) -> double {
942 if (particle == nullptr)
943 return Const::doubleNaN;
944 const Particle* dau_i = particle->getParticleFromGeneralizedIndexString(arguments[0]);
945 const Particle* dau_j = particle->getParticleFromGeneralizedIndexString(arguments[1]);
946 auto variablename = arguments[2];
947 if (dau_i == nullptr || dau_j == nullptr)
948 {
949 B2ERROR("One of the first two arguments doesn't specify a valid (grand-)daughter!");
950 return Const::doubleNaN;
951 }
952 const Variable::Manager::Var* var = Manager::Instance().getVariable(variablename);
953 auto result_j = var->function(dau_j);
954 auto result_i = var->function(dau_i);
955 double diff = Const::doubleNaN;
956 if (std::holds_alternative<double>(result_j) && std::holds_alternative<double>(result_i))
957 {
958 diff = std::get<double>(result_j) - std::get<double>(result_i);
959 } else if (std::holds_alternative<int>(result_j) && std::holds_alternative<int>(result_i))
960 {
961 diff = std::get<int>(result_j) - std::get<int>(result_i);
962 } else
963 {
964 throw std::runtime_error("Bad variant access");
965 }
966 if (variablename == "phi" or variablename == "clusterPhi" or std::regex_match(variablename, std::regex("use.*Frame\\(phi\\)"))
967 or std::regex_match(variablename, std::regex("use.*Frame\\(clusterPhi\\)")))
968 {
969 if (fabs(diff) > M_PI) {
970 if (diff > M_PI) {
971 diff = diff - 2 * M_PI;
972 } else {
973 diff = 2 * M_PI + diff;
974 }
975 }
976 }
977 return diff;
978 };
979 return func;
980 } else {
981 B2FATAL("Wrong number of arguments for meta function daughterDiffOf");
982 }
983 }
984
985 Manager::FunctionPtr mcDaughterDiffOf(const std::vector<std::string>& arguments)
986 {
987 if (arguments.size() == 3) {
988 auto func = [arguments](const Particle * particle) -> double {
989 if (particle == nullptr)
990 return Const::doubleNaN;
991 const Particle* dau_i = particle->getParticleFromGeneralizedIndexString(arguments[0]);
992 const Particle* dau_j = particle->getParticleFromGeneralizedIndexString(arguments[1]);
993 auto variablename = arguments[2];
994 if (dau_i == nullptr || dau_j == nullptr)
995 {
996 B2ERROR("One of the first two arguments doesn't specify a valid (grand-)daughter!");
997 return Const::doubleNaN;
998 }
999 const MCParticle* iMcDaughter = dau_i->getMCParticle();
1000 const MCParticle* jMcDaughter = dau_j->getMCParticle();
1001 if (iMcDaughter == nullptr || jMcDaughter == nullptr)
1002 return Const::doubleNaN;
1003 Particle iTmpPart(iMcDaughter);
1004 Particle jTmpPart(jMcDaughter);
1005 const Variable::Manager::Var* var = Manager::Instance().getVariable(variablename);
1006 auto result_j = var->function(&jTmpPart);
1007 auto result_i = var->function(&iTmpPart);
1008 double diff = Const::doubleNaN;
1009 if (std::holds_alternative<double>(result_j) && std::holds_alternative<double>(result_i))
1010 {
1011 diff = std::get<double>(result_j) - std::get<double>(result_i);
1012 } else if (std::holds_alternative<int>(result_j) && std::holds_alternative<int>(result_i))
1013 {
1014 diff = std::get<int>(result_j) - std::get<int>(result_i);
1015 } else
1016 {
1017 throw std::runtime_error("Bad variant access");
1018 }
1019 if (variablename == "phi" or std::regex_match(variablename, std::regex("use.*Frame\\(phi\\)")))
1020 {
1021 if (fabs(diff) > M_PI) {
1022 if (diff > M_PI) {
1023 diff = diff - 2 * M_PI;
1024 } else {
1025 diff = 2 * M_PI + diff;
1026 }
1027 }
1028 }
1029 return diff;
1030 };
1031 return func;
1032 } else {
1033 B2FATAL("Wrong number of arguments for meta function mcDaughterDiffOf");
1034 }
1035 }
1036
1037 Manager::FunctionPtr grandDaughterDiffOf(const std::vector<std::string>& arguments)
1038 {
1039 if (arguments.size() == 5) {
1040 try {
1041 convertString<int>(arguments[0]);
1042 convertString<int>(arguments[1]);
1043 convertString<int>(arguments[2]);
1044 convertString<int>(arguments[3]);
1045 } catch (std::invalid_argument&) {
1046 B2FATAL("First four arguments of grandDaughterDiffOf meta function must be integers!");
1047 }
1048 std::vector<std::string> new_arguments;
1049 new_arguments.push_back(std::string(arguments[0] + ":" + arguments[2]));
1050 new_arguments.push_back(std::string(arguments[1] + ":" + arguments[3]));
1051 new_arguments.push_back(arguments[4]);
1052 return daughterDiffOf(new_arguments);
1053 } else {
1054 B2FATAL("Wrong number of arguments for meta function grandDaughterDiffOf");
1055 }
1056 }
1057
1058 Manager::FunctionPtr daughterDiffOfPhi(const std::vector<std::string>& arguments)
1059 {
1060 std::vector<std::string> new_arguments = arguments;
1061 new_arguments.push_back(std::string("phi"));
1062 return daughterDiffOf(new_arguments);
1063 }
1064
1065 Manager::FunctionPtr mcDaughterDiffOfPhi(const std::vector<std::string>& arguments)
1066 {
1067 std::vector<std::string> new_arguments = arguments;
1068 new_arguments.push_back(std::string("phi"));
1069 return mcDaughterDiffOf(new_arguments);
1070 }
1071
1072 Manager::FunctionPtr grandDaughterDiffOfPhi(const std::vector<std::string>& arguments)
1073 {
1074 std::vector<std::string> new_arguments = arguments;
1075 new_arguments.push_back(std::string("phi"));
1076 return grandDaughterDiffOf(new_arguments);
1077 }
1078
1079 Manager::FunctionPtr daughterDiffOfClusterPhi(const std::vector<std::string>& arguments)
1080 {
1081 std::vector<std::string> new_arguments = arguments;
1082 new_arguments.push_back(std::string("clusterPhi"));
1083 return daughterDiffOf(new_arguments);
1084 }
1085
1086 Manager::FunctionPtr grandDaughterDiffOfClusterPhi(const std::vector<std::string>& arguments)
1087 {
1088 std::vector<std::string> new_arguments = arguments;
1089 new_arguments.push_back(std::string("clusterPhi"));
1090 return grandDaughterDiffOf(new_arguments);
1091 }
1092
1093 Manager::FunctionPtr daughterDiffOfPhiCMS(const std::vector<std::string>& arguments)
1094 {
1095 std::vector<std::string> new_arguments = arguments;
1096 new_arguments.push_back(std::string("useCMSFrame(phi)"));
1097 return daughterDiffOf(new_arguments);
1098 }
1099
1100 Manager::FunctionPtr mcDaughterDiffOfPhiCMS(const std::vector<std::string>& arguments)
1101 {
1102 std::vector<std::string> new_arguments = arguments;
1103 new_arguments.push_back(std::string("useCMSFrame(phi)"));
1104 return mcDaughterDiffOf(new_arguments);
1105 }
1106
1107 Manager::FunctionPtr daughterDiffOfClusterPhiCMS(const std::vector<std::string>& arguments)
1108 {
1109 std::vector<std::string> new_arguments = arguments;
1110 new_arguments.push_back(std::string("useCMSFrame(clusterPhi)"));
1111 return daughterDiffOf(new_arguments);
1112 }
1113
1114 Manager::FunctionPtr daughterNormDiffOf(const std::vector<std::string>& arguments)
1115 {
1116 if (arguments.size() == 3) {
1117 auto func = [arguments](const Particle * particle) -> double {
1118 if (particle == nullptr)
1119 return Const::doubleNaN;
1120 const Particle* dau_i = particle->getParticleFromGeneralizedIndexString(arguments[0]);
1121 const Particle* dau_j = particle->getParticleFromGeneralizedIndexString(arguments[1]);
1122 if (!(dau_i && dau_j))
1123 {
1124 B2ERROR("One of the first two arguments doesn't specify a valid (grand-)daughter!");
1125 return Const::doubleNaN;
1126 }
1127 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[2]);
1128 double iValue, jValue;
1129 if (std::holds_alternative<double>(var->function(dau_j)))
1130 {
1131 iValue = std::get<double>(var->function(dau_i));
1132 jValue = std::get<double>(var->function(dau_j));
1133 } else if (std::holds_alternative<int>(var->function(dau_j)))
1134 {
1135 iValue = std::get<int>(var->function(dau_i));
1136 jValue = std::get<int>(var->function(dau_j));
1137 } else return Const::doubleNaN;
1138 return (jValue - iValue) / (jValue + iValue);
1139 };
1140 return func;
1141 } else {
1142 B2FATAL("Wrong number of arguments for meta function daughterNormDiffOf");
1143 }
1144 }
1145
1146 Manager::FunctionPtr daughterMotherDiffOf(const std::vector<std::string>& arguments)
1147 {
1148 if (arguments.size() == 2) {
1149 int daughterNumber = 0;
1150 try {
1151 daughterNumber = convertString<int>(arguments[0]);
1152 } catch (std::invalid_argument&) {
1153 B2FATAL("First argument of daughterMotherDiffOf meta function must be integer!");
1154 }
1155 auto variablename = arguments[1];
1156 auto func = [variablename, daughterNumber](const Particle * particle) -> double {
1157 if (particle == nullptr)
1158 return Const::doubleNaN;
1159 if (daughterNumber >= int(particle->getNDaughters()))
1160 return Const::doubleNaN;
1161 else
1162 {
1163 const Variable::Manager::Var* var = Manager::Instance().getVariable(variablename);
1164 auto result_mother = var->function(particle);
1165 auto result_daughter = var->function(particle->getDaughter(daughterNumber));
1166 double diff = Const::doubleNaN;
1167 if (std::holds_alternative<double>(result_mother) && std::holds_alternative<double>(result_daughter)) {
1168 diff = std::get<double>(result_mother) - std::get<double>(result_daughter);
1169 } else if (std::holds_alternative<int>(result_mother) && std::holds_alternative<int>(result_daughter)) {
1170 diff = std::get<int>(result_mother) - std::get<int>(result_daughter);
1171 } else {
1172 throw std::runtime_error("Bad variant access");
1173 }
1174
1175 if (variablename == "phi" or variablename == "useCMSFrame(phi)") {
1176 if (fabs(diff) > M_PI) {
1177 if (diff > M_PI) {
1178 diff = diff - 2 * M_PI;
1179 } else {
1180 diff = 2 * M_PI + diff;
1181 }
1182 }
1183 }
1184 return diff;
1185 }
1186 };
1187 return func;
1188 } else {
1189 B2FATAL("Wrong number of arguments for meta function daughterMotherDiffOf");
1190 }
1191 }
1192
1193 Manager::FunctionPtr daughterMotherNormDiffOf(const std::vector<std::string>& arguments)
1194 {
1195 if (arguments.size() == 2) {
1196 int daughterNumber = 0;
1197 try {
1198 daughterNumber = convertString<int>(arguments[0]);
1199 } catch (std::invalid_argument&) {
1200 B2FATAL("First argument of daughterMotherDiffOf meta function must be integer!");
1201 }
1202 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[1]);
1203 auto func = [var, daughterNumber](const Particle * particle) -> double {
1204 if (particle == nullptr)
1205 return Const::doubleNaN;
1206 if (daughterNumber >= int(particle->getNDaughters()))
1207 return Const::doubleNaN;
1208 else
1209 {
1210 double daughterValue = 0.0, motherValue = 0.0;
1211 if (std::holds_alternative<double>(var->function(particle))) {
1212 daughterValue = std::get<double>(var->function(particle->getDaughter(daughterNumber)));
1213 motherValue = std::get<double>(var->function(particle));
1214 } else if (std::holds_alternative<int>(var->function(particle))) {
1215 daughterValue = std::get<int>(var->function(particle->getDaughter(daughterNumber)));
1216 motherValue = std::get<int>(var->function(particle));
1217 }
1218 return (motherValue - daughterValue) / (motherValue + daughterValue);
1219 }
1220 };
1221 return func;
1222 } else {
1223 B2FATAL("Wrong number of arguments for meta function daughterMotherNormDiffOf");
1224 }
1225 }
1226
1227 Manager::FunctionPtr angleBetweenDaughterAndRecoil(const std::vector<std::string>& arguments)
1228 {
1229 if (arguments.size() >= 1) {
1230
1231 auto func = [arguments](const Particle * particle) -> double {
1232 if (particle == nullptr)
1233 return Const::doubleNaN;
1234
1235 const auto& frame = ReferenceFrame::GetCurrent();
1236
1237 ROOT::Math::PxPyPzEVector pSum(0, 0, 0, 0);
1238 for (auto& generalizedIndex : arguments)
1239 {
1240 const Particle* dauPart = particle->getParticleFromGeneralizedIndexString(generalizedIndex);
1241 if (dauPart) pSum += frame.getMomentum(dauPart);
1242 else {
1243 B2WARNING("Trying to access a daughter that does not exist. Index = " << generalizedIndex);
1244 return Const::doubleNaN;
1245 }
1246 }
1247
1248 PCmsLabTransform T;
1249 ROOT::Math::PxPyPzEVector pIN = T.getBeamFourMomentum(); // Initial state (e+e- momentum in LAB)
1250 ROOT::Math::PxPyPzEVector pRecoil = frame.getMomentum(pIN - particle->get4Vector());
1251
1252 return ROOT::Math::VectorUtil::Angle(pRecoil, pSum);
1253 };
1254 return func;
1255 } else {
1256 B2FATAL("Wrong number of arguments for meta function angleBetweenDaughterAndRecoil");
1257 }
1258 }
1259
1260 Manager::FunctionPtr angleBetweenDaughterAndMissingMomentum(const std::vector<std::string>& arguments)
1261 {
1262 if (arguments.size() >= 1) {
1263 auto func = [arguments](const Particle * particle) -> double {
1264 if (particle == nullptr)
1265 return Const::doubleNaN;
1266
1267 StoreObjPtr<EventKinematics> evtShape;
1268 if (!evtShape)
1269 {
1270 B2WARNING("Cannot find missing momentum information, did you forget to run EventKinematicsModule?");
1271 return Const::doubleNaN;
1272 }
1273 ROOT::Math::XYZVector missingMomentumCMS = evtShape->getMissingMomentumCMS();
1274 ROOT::Math::PxPyPzEVector missingTotalMomentumCMS(missingMomentumCMS.X(),
1275 missingMomentumCMS.Y(),
1276 missingMomentumCMS.Z(),
1277 evtShape->getMissingEnergyCMS());
1278 PCmsLabTransform T;
1279 ROOT::Math::PxPyPzEVector missingTotalMomentumLab = T.rotateCmsToLab() * missingTotalMomentumCMS;
1280
1281 const auto& frame = ReferenceFrame::GetCurrent();
1282 ROOT::Math::PxPyPzEVector pMiss = frame.getMomentum(missingTotalMomentumLab); // transform from lab to reference frame
1283
1284 ROOT::Math::PxPyPzEVector pSum(0, 0, 0, 0);
1285 for (auto& generalizedIndex : arguments)
1286 {
1287 const Particle* dauPart = particle->getParticleFromGeneralizedIndexString(generalizedIndex);
1288 if (dauPart) pSum += frame.getMomentum(dauPart);
1289 else {
1290 B2WARNING("Trying to access a daughter that does not exist. Index = " << generalizedIndex);
1291 return Const::doubleNaN;
1292 }
1293 }
1294
1295 return ROOT::Math::VectorUtil::Angle(pMiss, pSum);
1296 };
1297 return func;
1298 } else {
1299 B2FATAL("Wrong number of arguments for meta function angleBetweenDaughterAndMissingMomentum");
1300 }
1301 }
1302
1303 Manager::FunctionPtr daughterAngle(const std::vector<std::string>& arguments)
1304 {
1305 if (arguments.size() == 2 || arguments.size() == 3) {
1306
1307 auto func = [arguments](const Particle * particle) -> double {
1308 if (particle == nullptr)
1309 return Const::doubleNaN;
1310
1311 std::vector<ROOT::Math::PxPyPzEVector> pDaus;
1312 const auto& frame = ReferenceFrame::GetCurrent();
1313
1314 // Parses the generalized indexes and fetches the 4-momenta of the particles of interest
1315 for (auto& generalizedIndex : arguments)
1316 {
1317 const Particle* dauPart = particle->getParticleFromGeneralizedIndexString(generalizedIndex);
1318 if (dauPart)
1319 pDaus.push_back(frame.getMomentum(dauPart));
1320 else {
1321 B2WARNING("Trying to access a daughter that does not exist. Index = " << generalizedIndex);
1322 return Const::doubleNaN;
1323 }
1324 }
1325
1326 // Calculates the angle between the selected particles
1327 if (pDaus.size() == 2)
1328 return ROOT::Math::VectorUtil::Angle(pDaus[0], pDaus[1]);
1329 else
1330 return ROOT::Math::VectorUtil::Angle(pDaus[2], pDaus[0] + pDaus[1]);
1331 };
1332 return func;
1333 } else {
1334 B2FATAL("Wrong number of arguments for meta function daughterAngle");
1335 }
1336 }
1337
1338 double grandDaughterDecayAngle(const Particle* particle, const std::vector<double>& arguments)
1339 {
1340 if (arguments.size() == 2) {
1341
1342 if (!particle)
1343 return Const::doubleNaN;
1344
1345 int daughterIndex = std::lround(arguments[0]);
1346 if (daughterIndex >= int(particle->getNDaughters()))
1347 return Const::doubleNaN;
1348 const Particle* dau = particle->getDaughter(daughterIndex);
1349
1350 int grandDaughterIndex = std::lround(arguments[1]);
1351 if (grandDaughterIndex >= int(dau->getNDaughters()))
1352 return Const::doubleNaN;
1353
1354 ROOT::Math::XYZVector boost = dau->get4Vector().BoostToCM();
1355
1356 ROOT::Math::PxPyPzEVector motherMomentum = - particle->get4Vector();
1357 motherMomentum = ROOT::Math::Boost(boost) * motherMomentum;
1358
1359 ROOT::Math::PxPyPzEVector grandDaughterMomentum = dau->getDaughter(grandDaughterIndex)->get4Vector();
1360 grandDaughterMomentum = ROOT::Math::Boost(boost) * grandDaughterMomentum;
1361
1362 return ROOT::Math::VectorUtil::Angle(motherMomentum, grandDaughterMomentum);
1363
1364 } else {
1365 B2FATAL("The variable grandDaughterDecayAngle needs exactly two integers as arguments!");
1366 }
1367 }
1368
1369 Manager::FunctionPtr mcDaughterAngle(const std::vector<std::string>& arguments)
1370 {
1371 if (arguments.size() == 2 || arguments.size() == 3) {
1372
1373 auto func = [arguments](const Particle * particle) -> double {
1374 if (particle == nullptr)
1375 return Const::doubleNaN;
1376
1377 std::vector<ROOT::Math::PxPyPzEVector> pDaus;
1378 const auto& frame = ReferenceFrame::GetCurrent();
1379
1380 // Parses the generalized indexes and fetches the 4-momenta of the particles of interest
1381 if (particle->getParticleSource() == Particle::EParticleSourceObject::c_MCParticle) // Check if MCParticle
1382 {
1383 for (auto& generalizedIndex : arguments) {
1384 const MCParticle* mcPart = particle->getMCParticle();
1385 if (mcPart == nullptr)
1386 return Const::doubleNaN;
1387 const MCParticle* dauMcPart = mcPart->getParticleFromGeneralizedIndexString(generalizedIndex);
1388 if (dauMcPart == nullptr)
1389 return Const::doubleNaN;
1390
1391 pDaus.push_back(frame.getMomentum(dauMcPart->get4Vector()));
1392 }
1393 } else
1394 {
1395 for (auto& generalizedIndex : arguments) {
1396 const Particle* dauPart = particle->getParticleFromGeneralizedIndexString(generalizedIndex);
1397 if (dauPart == nullptr)
1398 return Const::doubleNaN;
1399
1400 const MCParticle* dauMcPart = dauPart->getMCParticle();
1401 if (dauMcPart == nullptr)
1402 return Const::doubleNaN;
1403
1404 pDaus.push_back(frame.getMomentum(dauMcPart->get4Vector()));
1405 }
1406 }
1407
1408 // Calculates the angle between the selected particles
1409 if (pDaus.size() == 2)
1410 return ROOT::Math::VectorUtil::Angle(pDaus[0], pDaus[1]);
1411 else
1412 return ROOT::Math::VectorUtil::Angle(pDaus[2], pDaus[0] + pDaus[1]);
1413 };
1414 return func;
1415 } else {
1416 B2FATAL("Wrong number of arguments for meta function mcDaughterAngle");
1417 }
1418 }
1419
1420 double daughterClusterAngleInBetween(const Particle* particle, const std::vector<double>& daughterIndices)
1421 {
1422 if (daughterIndices.size() == 2) {
1423 int daughterIndexi = std::lround(daughterIndices[0]);
1424 int daughterIndexj = std::lround(daughterIndices[1]);
1425 if (std::max(daughterIndexi, daughterIndexj) >= int(particle->getNDaughters())) {
1426 return Const::doubleNaN;
1427 } else {
1428 const ECLCluster* clusteri = particle->getDaughter(daughterIndexi)->getECLCluster();
1429 const ECLCluster* clusterj = particle->getDaughter(daughterIndexj)->getECLCluster();
1430 if (clusteri and clusterj) {
1431 const auto& frame = ReferenceFrame::GetCurrent();
1432 const ECLCluster::EHypothesisBit clusteriBit = (particle->getDaughter(daughterIndexi))->getECLClusterEHypothesisBit();
1433 const ECLCluster::EHypothesisBit clusterjBit = (particle->getDaughter(daughterIndexj))->getECLClusterEHypothesisBit();
1434 ClusterUtils clusutils;
1435 ROOT::Math::PxPyPzEVector pi = frame.getMomentum(clusutils.Get4MomentumFromCluster(clusteri, clusteriBit));
1436 ROOT::Math::PxPyPzEVector pj = frame.getMomentum(clusutils.Get4MomentumFromCluster(clusterj, clusterjBit));
1437 return ROOT::Math::VectorUtil::Angle(pi, pj);
1438 }
1439 return Const::doubleNaN;
1440 }
1441 } else if (daughterIndices.size() == 3) {
1442 int daughterIndexi = std::lround(daughterIndices[0]);
1443 int daughterIndexj = std::lround(daughterIndices[1]);
1444 int daughterIndexk = std::lround(daughterIndices[2]);
1445 if (std::max(std::max(daughterIndexi, daughterIndexj), daughterIndexk) >= int(particle->getNDaughters())) {
1446 return Const::doubleNaN;
1447 } else {
1448 const ECLCluster* clusteri = (particle->getDaughter(daughterIndices[0]))->getECLCluster();
1449 const ECLCluster* clusterj = (particle->getDaughter(daughterIndices[1]))->getECLCluster();
1450 const ECLCluster* clusterk = (particle->getDaughter(daughterIndices[2]))->getECLCluster();
1451 if (clusteri and clusterj and clusterk) {
1452 const auto& frame = ReferenceFrame::GetCurrent();
1453 const ECLCluster::EHypothesisBit clusteriBit = (particle->getDaughter(daughterIndices[0]))->getECLClusterEHypothesisBit();
1454 const ECLCluster::EHypothesisBit clusterjBit = (particle->getDaughter(daughterIndices[1]))->getECLClusterEHypothesisBit();
1455 const ECLCluster::EHypothesisBit clusterkBit = (particle->getDaughter(daughterIndices[2]))->getECLClusterEHypothesisBit();
1456 ClusterUtils clusutils;
1457 ROOT::Math::PxPyPzEVector pi = frame.getMomentum(clusutils.Get4MomentumFromCluster(clusteri, clusteriBit));
1458 ROOT::Math::PxPyPzEVector pj = frame.getMomentum(clusutils.Get4MomentumFromCluster(clusterj, clusterjBit));
1459 ROOT::Math::PxPyPzEVector pk = frame.getMomentum(clusutils.Get4MomentumFromCluster(clusterk, clusterkBit));
1460 return ROOT::Math::VectorUtil::Angle(pk, pi + pj);
1461 }
1462 return Const::doubleNaN;
1463 }
1464 } else {
1465 B2FATAL("Wrong number of arguments for daughterClusterAngleInBetween!");
1466 }
1467 }
1468
1469 Manager::FunctionPtr daughterInvM(const std::vector<std::string>& arguments)
1470 {
1471 if (arguments.size() > 1) {
1472 auto func = [arguments](const Particle * particle) -> double {
1473 const auto& frame = ReferenceFrame::GetCurrent();
1474 ROOT::Math::PxPyPzEVector pSum;
1475
1476 for (auto& generalizedIndex : arguments)
1477 {
1478 const Particle* dauPart = particle->getParticleFromGeneralizedIndexString(generalizedIndex);
1479 if (dauPart)
1480 pSum += frame.getMomentum(dauPart);
1481 else {
1482 return Const::doubleNaN;
1483 }
1484 }
1485 return pSum.M();
1486 };
1487 return func;
1488 } else {
1489 B2FATAL("Wrong number of arguments for meta function daughterInvM. At least two integers are needed.");
1490 }
1491 }
1492
1493 Manager::FunctionPtr modulo(const std::vector<std::string>& arguments)
1494 {
1495 if (arguments.size() == 2) {
1496 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
1497 int divideBy = 1;
1498 try {
1499 divideBy = convertString<int>(arguments[1]);
1500 } catch (std::invalid_argument&) {
1501 B2FATAL("Second argument of modulo meta function must be integer!");
1502 }
1503 auto func = [var, divideBy](const Particle * particle) -> int {
1504 auto var_result = var->function(particle);
1505 if (std::holds_alternative<double>(var_result))
1506 {
1507 return int(std::get<double>(var_result)) % divideBy;
1508 } else if (std::holds_alternative<int>(var_result))
1509 {
1510 return std::get<int>(var_result) % divideBy;
1511 } else if (std::holds_alternative<bool>(var_result))
1512 {
1513 return int(std::get<bool>(var_result)) % divideBy;
1514 } else return 0;
1515 };
1516 return func;
1517 } else {
1518 B2FATAL("Wrong number of arguments for meta function modulo");
1519 }
1520 }
1521
1522 Manager::FunctionPtr isNAN(const std::vector<std::string>& arguments)
1523 {
1524 if (arguments.size() == 1) {
1525 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
1526
1527 auto func = [var](const Particle * particle) -> bool { return std::isnan(std::get<double>(var->function(particle))); };
1528 return func;
1529 } else {
1530 B2FATAL("Wrong number of arguments for meta function isNAN");
1531 }
1532 }
1533
1534 Manager::FunctionPtr ifNANgiveX(const std::vector<std::string>& arguments)
1535 {
1536 if (arguments.size() == 2) {
1537 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
1538 double defaultOutput;
1539 try {
1540 defaultOutput = convertString<double>(arguments[1]);
1541 } catch (std::invalid_argument&) {
1542 B2FATAL("The second argument of ifNANgiveX meta function must be a number!");
1543 }
1544 auto func = [var, defaultOutput](const Particle * particle) -> double {
1545 double output = std::get<double>(var->function(particle));
1546 if (std::isnan(output)) return defaultOutput;
1547 else return output;
1548 };
1549 return func;
1550 } else {
1551 B2FATAL("Wrong number of arguments for meta function ifNANgiveX");
1552 }
1553 }
1554
1555 Manager::FunctionPtr isInfinity(const std::vector<std::string>& arguments)
1556 {
1557 if (arguments.size() == 1) {
1558 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
1559
1560 auto func = [var](const Particle * particle) -> bool { return std::isinf(std::get<double>(var->function(particle))); };
1561 return func;
1562 } else {
1563 B2FATAL("Wrong number of arguments for meta function isInfinity");
1564 }
1565 }
1566
1567 Manager::FunctionPtr unmask(const std::vector<std::string>& arguments)
1568 {
1569 if (arguments.size() >= 2) {
1570 // get the function pointer of variable to be unmasked
1571 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
1572
1573 // get the final mask which summarize all the input masks
1574 int finalMask = 0;
1575 for (size_t i = 1; i < arguments.size(); ++i) {
1576 try {
1577 finalMask |= convertString<int>(arguments[i]);
1578 } catch (std::invalid_argument&) {
1579 B2FATAL("The input flags to meta function unmask() should be integer!");
1580 return nullptr;
1581 }
1582 }
1583
1584 // unmask the variable
1585 auto func = [var, finalMask](const Particle * particle) -> double {
1586 int value = 0;
1587 auto var_result = var->function(particle);
1588 if (std::holds_alternative<double>(var_result))
1589 {
1590 // judge if the value is nan before unmasking
1591 if (std::isnan(std::get<double>(var_result))) {
1592 return Const::doubleNaN;
1593 }
1594 value = int(std::get<double>(var_result));
1595 } else if (std::holds_alternative<int>(var_result))
1596 {
1597 value = std::get<int>(var_result);
1598 }
1599
1600 // apply the final mask
1601 value &= (~finalMask);
1602
1603 return value;
1604 };
1605 return func;
1606
1607 } else {
1608 B2FATAL("Meta function unmask needs at least two arguments!");
1609 }
1610 }
1611
1612 Manager::FunctionPtr conditionalVariableSelector(const std::vector<std::string>& arguments)
1613 {
1614 if (arguments.size() == 3) {
1615
1616 std::string cutString = arguments[0];
1617 std::shared_ptr<Variable::Cut> cut = std::shared_ptr<Variable::Cut>(Variable::Cut::compile(cutString));
1618
1619 const Variable::Manager::Var* variableIfTrue = Manager::Instance().getVariable(arguments[1]);
1620 const Variable::Manager::Var* variableIfFalse = Manager::Instance().getVariable(arguments[2]);
1621
1622 auto func = [cut, variableIfTrue, variableIfFalse](const Particle * particle) -> double {
1623 if (particle == nullptr)
1624 return Const::doubleNaN;
1625 if (cut->check(particle))
1626 {
1627 auto var_result = variableIfTrue->function(particle);
1628 if (std::holds_alternative<double>(var_result)) {
1629 return std::get<double>(var_result);
1630 } else if (std::holds_alternative<int>(var_result)) {
1631 return std::get<int>(var_result);
1632 } else if (std::holds_alternative<bool>(var_result)) {
1633 return std::get<bool>(var_result);
1634 } else return Const::doubleNaN;
1635 } else
1636 {
1637 auto var_result = variableIfFalse->function(particle);
1638 if (std::holds_alternative<double>(var_result)) {
1639 return std::get<double>(var_result);
1640 } else if (std::holds_alternative<int>(var_result)) {
1641 return std::get<int>(var_result);
1642 } else if (std::holds_alternative<bool>(var_result)) {
1643 return std::get<bool>(var_result);
1644 } else return Const::doubleNaN;
1645 }
1646 };
1647 return func;
1648
1649 } else {
1650 B2FATAL("Wrong number of arguments for meta function conditionalVariableSelector");
1651 }
1652 }
1653
1654 Manager::FunctionPtr pValueCombination(const std::vector<std::string>& arguments)
1655 {
1656 if (arguments.size() > 0) {
1657 std::vector<const Variable::Manager::Var*> variables;
1658 for (auto& argument : arguments)
1659 variables.push_back(Manager::Instance().getVariable(argument));
1660
1661 auto func = [variables, arguments](const Particle * particle) -> double {
1662 double pValueProduct = 1.;
1663 for (auto variable : variables)
1664 {
1665 double pValue = std::get<double>(variable->function(particle));
1666 if (pValue < 0)
1667 return -1;
1668 else
1669 pValueProduct *= pValue;
1670 }
1671 double pValueSum = 1.;
1672 double factorial = 1.;
1673 for (unsigned int i = 1; i < arguments.size(); ++i)
1674 {
1675 factorial *= i;
1676 pValueSum += pow(-std::log(pValueProduct), i) / factorial;
1677 }
1678 return pValueProduct * pValueSum;
1679 };
1680 return func;
1681 } else {
1682 B2FATAL("Wrong number of arguments for meta function pValueCombination");
1683 }
1684 }
1685
1686 Manager::FunctionPtr pValueCombinationOfDaughters(const std::vector<std::string>& arguments)
1687 {
1688 if (arguments.size() == 1) {
1689 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
1690 auto func = [var](const Particle * particle) -> double {
1691 double pValueProduct = 1.;
1692 if (particle->getNDaughters() == 0)
1693 {
1694 return Const::doubleNaN;
1695 }
1696
1697 for (unsigned j = 0; j < particle->getNDaughters(); ++j)
1698 {
1699 double pValue = std::get<double>(var->function(particle->getDaughter(j)));
1700 if (pValue < 0) return -1;
1701 else pValueProduct *= pValue;
1702 }
1703
1704 double pValueSum = 1.;
1705 double factorial = 1.;
1706 for (unsigned int i = 1; i < particle->getNDaughters(); ++i)
1707 {
1708 factorial *= i;
1709 pValueSum += pow(-std::log(pValueProduct), i) / factorial;
1710 }
1711 return pValueProduct * pValueSum;
1712 };
1713 return func;
1714 } else {
1715 B2FATAL("Wrong number of arguments for meta function pValueCombinationOfDaughters");
1716 }
1717 }
1718
1719 Manager::FunctionPtr abs(const std::vector<std::string>& arguments)
1720 {
1721 if (arguments.size() == 1) {
1722 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
1723 auto func = [var](const Particle * particle) -> double {
1724 auto var_result = var->function(particle);
1725 if (std::holds_alternative<double>(var_result))
1726 {
1727 return std::abs(std::get<double>(var_result));
1728 } else if (std::holds_alternative<int>(var_result))
1729 {
1730 return std::abs(std::get<int>(var_result));
1731 } else return Const::doubleNaN;
1732 };
1733 return func;
1734 } else {
1735 B2FATAL("Wrong number of arguments for meta function abs");
1736 }
1737 }
1738
1739 Manager::FunctionPtr max(const std::vector<std::string>& arguments)
1740 {
1741 if (arguments.size() == 2) {
1742 const Variable::Manager::Var* var1 = Manager::Instance().getVariable(arguments[0]);
1743 const Variable::Manager::Var* var2 = Manager::Instance().getVariable(arguments[1]);
1744
1745 if (!var1 or !var2)
1746 B2FATAL("One or both of the used variables doesn't exist!");
1747
1748 auto func = [var1, var2](const Particle * particle) -> double {
1749 double val1, val2;
1750 auto var_result1 = var1->function(particle);
1751 auto var_result2 = var2->function(particle);
1752 if (std::holds_alternative<double>(var_result1))
1753 {
1754 val1 = std::get<double>(var_result1);
1755 } else if (std::holds_alternative<int>(var_result1))
1756 {
1757 val1 = std::get<int>(var_result1);
1758 }
1759 if (std::holds_alternative<double>(var_result2))
1760 {
1761 val2 = std::get<double>(var_result2);
1762 } else if (std::holds_alternative<int>(var_result2))
1763 {
1764 val2 = std::get<int>(var_result2);
1765 }
1766 return std::max(val1, val2);
1767 };
1768 return func;
1769 } else {
1770 B2FATAL("Wrong number of arguments for meta function max");
1771 }
1772 }
1773
1774 Manager::FunctionPtr min(const std::vector<std::string>& arguments)
1775 {
1776 if (arguments.size() == 2) {
1777 const Variable::Manager::Var* var1 = Manager::Instance().getVariable(arguments[0]);
1778 const Variable::Manager::Var* var2 = Manager::Instance().getVariable(arguments[1]);
1779
1780 if (!var1 or !var2)
1781 B2FATAL("One or both of the used variables doesn't exist!");
1782
1783 auto func = [var1, var2](const Particle * particle) -> double {
1784 double val1, val2;
1785 auto var_result1 = var1->function(particle);
1786 auto var_result2 = var2->function(particle);
1787 if (std::holds_alternative<double>(var_result1))
1788 {
1789 val1 = std::get<double>(var_result1);
1790 } else if (std::holds_alternative<int>(var_result1))
1791 {
1792 val1 = std::get<int>(var_result1);
1793 }
1794 if (std::holds_alternative<double>(var_result2))
1795 {
1796 val2 = std::get<double>(var_result2);
1797 } else if (std::holds_alternative<int>(var_result2))
1798 {
1799 val2 = std::get<int>(var_result2);
1800 }
1801 return std::min(val1, val2);
1802 };
1803 return func;
1804 } else {
1805 B2FATAL("Wrong number of arguments for meta function min");
1806 }
1807 }
1808
1809 Manager::FunctionPtr sin(const std::vector<std::string>& arguments)
1810 {
1811 if (arguments.size() == 1) {
1812 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
1813 auto func = [var](const Particle * particle) -> double {
1814 auto var_result = var->function(particle);
1815 if (std::holds_alternative<double>(var_result))
1816 return std::sin(std::get<double>(var_result));
1817 else if (std::holds_alternative<int>(var_result))
1818 return std::sin(std::get<int>(var_result));
1819 else return Const::doubleNaN;
1820 };
1821 return func;
1822 } else {
1823 B2FATAL("Wrong number of arguments for meta function sin");
1824 }
1825 }
1826
1827 Manager::FunctionPtr asin(const std::vector<std::string>& arguments)
1828 {
1829 if (arguments.size() == 1) {
1830 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
1831 auto func = [var](const Particle * particle) -> double {
1832 auto var_result = var->function(particle);
1833 if (std::holds_alternative<double>(var_result))
1834 return std::asin(std::get<double>(var_result));
1835 else if (std::holds_alternative<int>(var_result))
1836 return std::asin(std::get<int>(var_result));
1837 else return Const::doubleNaN;
1838 };
1839 return func;
1840 } else {
1841 B2FATAL("Wrong number of arguments for meta function asin");
1842 }
1843 }
1844
1845 Manager::FunctionPtr cos(const std::vector<std::string>& arguments)
1846 {
1847 if (arguments.size() == 1) {
1848 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
1849 auto func = [var](const Particle * particle) -> double {
1850 auto var_result = var->function(particle);
1851 if (std::holds_alternative<double>(var_result))
1852 return std::cos(std::get<double>(var_result));
1853 else if (std::holds_alternative<int>(var_result))
1854 return std::cos(std::get<int>(var_result));
1855 else return Const::doubleNaN;
1856 };
1857 return func;
1858 } else {
1859 B2FATAL("Wrong number of arguments for meta function cos");
1860 }
1861 }
1862
1863 Manager::FunctionPtr acos(const std::vector<std::string>& arguments)
1864 {
1865 if (arguments.size() == 1) {
1866 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
1867 auto func = [var](const Particle * particle) -> double {
1868 auto var_result = var->function(particle);
1869 if (std::holds_alternative<double>(var_result))
1870 return std::acos(std::get<double>(var_result));
1871 else if (std::holds_alternative<int>(var_result))
1872 return std::acos(std::get<int>(var_result));
1873 else return Const::doubleNaN;
1874 };
1875 return func;
1876 } else {
1877 B2FATAL("Wrong number of arguments for meta function acos");
1878 }
1879 }
1880
1881 Manager::FunctionPtr tan(const std::vector<std::string>& arguments)
1882 {
1883 if (arguments.size() == 1) {
1884 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
1885 auto func = [var](const Particle * particle) -> double { return std::tan(std::get<double>(var->function(particle))); };
1886 return func;
1887 } else {
1888 B2FATAL("Wrong number of arguments for meta function tan");
1889 }
1890 }
1891
1892 Manager::FunctionPtr atan(const std::vector<std::string>& arguments)
1893 {
1894 if (arguments.size() == 1) {
1895 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
1896 auto func = [var](const Particle * particle) -> double { return std::atan(std::get<double>(var->function(particle))); };
1897 return func;
1898 } else {
1899 B2FATAL("Wrong number of arguments for meta function atan");
1900 }
1901 }
1902
1903 Manager::FunctionPtr exp(const std::vector<std::string>& arguments)
1904 {
1905 if (arguments.size() == 1) {
1906 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
1907 auto func = [var](const Particle * particle) -> double {
1908 auto var_result = var->function(particle);
1909 if (std::holds_alternative<double>(var_result))
1910 return std::exp(std::get<double>(var_result));
1911 else if (std::holds_alternative<int>(var_result))
1912 return std::exp(std::get<int>(var_result));
1913 else return Const::doubleNaN;
1914 };
1915 return func;
1916 } else {
1917 B2FATAL("Wrong number of arguments for meta function exp");
1918 }
1919 }
1920
1921 Manager::FunctionPtr log(const std::vector<std::string>& arguments)
1922 {
1923 if (arguments.size() == 1) {
1924 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
1925 auto func = [var](const Particle * particle) -> double {
1926 auto var_result = var->function(particle);
1927 if (std::holds_alternative<double>(var_result))
1928 return std::log(std::get<double>(var_result));
1929 else if (std::holds_alternative<int>(var_result))
1930 return std::log(std::get<int>(var_result));
1931 else return Const::doubleNaN;
1932 };
1933 return func;
1934 } else {
1935 B2FATAL("Wrong number of arguments for meta function log");
1936 }
1937 }
1938
1939 Manager::FunctionPtr log10(const std::vector<std::string>& arguments)
1940 {
1941 if (arguments.size() == 1) {
1942 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
1943 auto func = [var](const Particle * particle) -> double {
1944 auto var_result = var->function(particle);
1945 if (std::holds_alternative<double>(var_result))
1946 return std::log10(std::get<double>(var_result));
1947 else if (std::holds_alternative<int>(var_result))
1948 return std::log10(std::get<int>(var_result));
1949 else return Const::doubleNaN;
1950 };
1951 return func;
1952 } else {
1953 B2FATAL("Wrong number of arguments for meta function log10");
1954 }
1955 }
1956
1957 Manager::FunctionPtr originalParticle(const std::vector<std::string>& arguments)
1958 {
1959 if (arguments.size() == 1) {
1960 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
1961 auto func = [var](const Particle * particle) -> double {
1962 if (particle == nullptr)
1963 return Const::doubleNaN;
1964
1965 StoreArray<Particle> particles;
1966 if (!particle->hasExtraInfo("original_index"))
1967 return Const::doubleNaN;
1968
1969 auto originalParticle = particles[particle->getExtraInfo("original_index")];
1970 if (!originalParticle)
1971 return Const::doubleNaN;
1972 auto var_result = var->function(originalParticle);
1973 if (std::holds_alternative<double>(var_result))
1974 {
1975 return std::get<double>(var_result);
1976 } else if (std::holds_alternative<int>(var_result))
1977 {
1978 return std::get<int>(var_result);
1979 } else if (std::holds_alternative<bool>(var_result))
1980 {
1981 return std::get<bool>(var_result);
1982 } else return Const::doubleNaN;
1983 };
1984 return func;
1985 } else {
1986 B2FATAL("Wrong number of arguments for meta function originalParticle");
1987 }
1988 }
1989
1990 Manager::FunctionPtr daughter(const std::vector<std::string>& arguments)
1991 {
1992 if (arguments.size() == 2) {
1993 int daughterNumber = 0;
1994 try {
1995 daughterNumber = convertString<int>(arguments[0]);
1996 } catch (std::invalid_argument&) {
1997 B2FATAL("First argument of daughter meta function must be integer!");
1998 }
1999 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[1]);
2000 auto func = [var, daughterNumber](const Particle * particle) -> double {
2001 if (particle == nullptr)
2002 return Const::doubleNaN;
2003 if (daughterNumber >= int(particle->getNDaughters()))
2004 return Const::doubleNaN;
2005 else
2006 {
2007 auto var_result = var->function(particle->getDaughter(daughterNumber));
2008 if (std::holds_alternative<double>(var_result)) {
2009 return std::get<double>(var_result);
2010 } else if (std::holds_alternative<int>(var_result)) {
2011 return std::get<int>(var_result);
2012 } else if (std::holds_alternative<bool>(var_result)) {
2013 return std::get<bool>(var_result);
2014 } else return Const::doubleNaN;
2015 }
2016 };
2017 return func;
2018 } else {
2019 B2FATAL("Wrong number of arguments for meta function daughter");
2020 }
2021 }
2022
2023 Manager::FunctionPtr originalDaughter(const std::vector<std::string>& arguments)
2024 {
2025 if (arguments.size() == 2) {
2026 int daughterNumber = 0;
2027 try {
2028 daughterNumber = convertString<int>(arguments[0]);
2029 } catch (std::invalid_argument&) {
2030 B2FATAL("First argument of daughter meta function must be integer!");
2031 }
2032 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[1]);
2033 auto func = [var, daughterNumber](const Particle * particle) -> double {
2034 if (particle == nullptr)
2035 return Const::doubleNaN;
2036 if (daughterNumber >= int(particle->getNDaughters()))
2037 return Const::doubleNaN;
2038 else
2039 {
2040 StoreArray<Particle> particles;
2041 if (!particle->getDaughter(daughterNumber)->hasExtraInfo("original_index"))
2042 return Const::doubleNaN;
2043 auto originalDaughter = particles[particle->getDaughter(daughterNumber)->getExtraInfo("original_index")];
2044 if (!originalDaughter)
2045 return Const::doubleNaN;
2046
2047 auto var_result = var->function(originalDaughter);
2048 if (std::holds_alternative<double>(var_result)) {
2049 return std::get<double>(var_result);
2050 } else if (std::holds_alternative<int>(var_result)) {
2051 return std::get<int>(var_result);
2052 } else if (std::holds_alternative<bool>(var_result)) {
2053 return std::get<bool>(var_result);
2054 } else return Const::doubleNaN;
2055 }
2056 };
2057 return func;
2058 } else {
2059 B2FATAL("Wrong number of arguments for meta function daughter");
2060 }
2061 }
2062
2063 Manager::FunctionPtr mcDaughter(const std::vector<std::string>& arguments)
2064 {
2065 if (arguments.size() == 2) {
2066 int daughterNumber = 0;
2067 try {
2068 daughterNumber = convertString<int>(arguments[0]);
2069 } catch (std::invalid_argument&) {
2070 B2FATAL("First argument of mcDaughter meta function must be integer!");
2071 }
2072 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[1]);
2073 auto func = [var, daughterNumber](const Particle * particle) -> double {
2074 if (particle == nullptr)
2075 return Const::doubleNaN;
2076 if (particle->getMCParticle()) // has MC match or is MCParticle
2077 {
2078 if (daughterNumber >= int(particle->getMCParticle()->getNDaughters())) {
2079 return Const::doubleNaN;
2080 }
2081 Particle tempParticle = Particle(particle->getMCParticle()->getDaughters().at(daughterNumber));
2082 auto var_result = var->function(&tempParticle);
2083 if (std::holds_alternative<double>(var_result)) {
2084 return std::get<double>(var_result);
2085 } else if (std::holds_alternative<int>(var_result)) {
2086 return std::get<int>(var_result);
2087 } else if (std::holds_alternative<bool>(var_result)) {
2088 return std::get<bool>(var_result);
2089 } else return Const::doubleNaN;
2090 } else
2091 {
2092 return Const::doubleNaN;
2093 }
2094 };
2095 return func;
2096 } else {
2097 B2FATAL("Wrong number of arguments for meta function mcDaughter");
2098 }
2099 }
2100
2101 Manager::FunctionPtr mcMother(const std::vector<std::string>& arguments)
2102 {
2103 if (arguments.size() == 1) {
2104 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
2105 auto func = [var](const Particle * particle) -> double {
2106 if (particle == nullptr)
2107 return Const::doubleNaN;
2108 if (particle->getMCParticle()) // has MC match or is MCParticle
2109 {
2110 if (particle->getMCParticle()->getMother() == nullptr) {
2111 return Const::doubleNaN;
2112 }
2113 Particle tempParticle = Particle(particle->getMCParticle()->getMother());
2114 auto var_result = var->function(&tempParticle);
2115 if (std::holds_alternative<double>(var_result)) {
2116 return std::get<double>(var_result);
2117 } else if (std::holds_alternative<int>(var_result)) {
2118 return std::get<int>(var_result);
2119 } else if (std::holds_alternative<bool>(var_result)) {
2120 return std::get<bool>(var_result);
2121 } else return Const::doubleNaN;
2122 } else
2123 {
2124 return Const::doubleNaN;
2125 }
2126 };
2127 return func;
2128 } else {
2129 B2FATAL("Wrong number of arguments for meta function mcMother");
2130 }
2131 }
2132
2133 Manager::FunctionPtr genParticle(const std::vector<std::string>& arguments)
2134 {
2135 if (arguments.size() == 2) {
2136 int particleNumber = 0;
2137 try {
2138 particleNumber = convertString<int>(arguments[0]);
2139 } catch (std::invalid_argument&) {
2140 B2FATAL("First argument of genParticle meta function must be integer!");
2141 }
2142 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[1]);
2143
2144 auto func = [var, particleNumber](const Particle*) -> double {
2145 StoreArray<MCParticle> mcParticles("MCParticles");
2146 if (particleNumber >= mcParticles.getEntries())
2147 {
2148 return Const::doubleNaN;
2149 }
2150
2151 MCParticle* mcParticle = mcParticles[particleNumber];
2152 Particle part = Particle(mcParticle);
2153 auto var_result = var->function(&part);
2154 if (std::holds_alternative<double>(var_result))
2155 {
2156 return std::get<double>(var_result);
2157 } else if (std::holds_alternative<int>(var_result))
2158 {
2159 return std::get<int>(var_result);
2160 } else if (std::holds_alternative<bool>(var_result))
2161 {
2162 return std::get<bool>(var_result);
2163 } else return Const::doubleNaN;
2164 };
2165 return func;
2166 } else {
2167 B2FATAL("Wrong number of arguments for meta function genParticle");
2168 }
2169 }
2170
2171 Manager::FunctionPtr genUpsilon4S(const std::vector<std::string>& arguments)
2172 {
2173 if (arguments.size() == 1) {
2174 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
2175
2176 auto func = [var](const Particle*) -> double {
2177 StoreArray<MCParticle> mcParticles("MCParticles");
2178 if (mcParticles.getEntries() == 0)
2179 {
2180 return Const::doubleNaN;
2181 }
2182
2183 MCParticle* mcUpsilon4S = mcParticles[0];
2184 if (mcUpsilon4S->isInitial()) mcUpsilon4S = mcParticles[2];
2185 if (mcUpsilon4S->getPDG() != 300553)
2186 {
2187 return Const::doubleNaN;
2188 }
2189
2190 Particle upsilon4S = Particle(mcUpsilon4S);
2191 auto var_result = var->function(&upsilon4S);
2192 if (std::holds_alternative<double>(var_result))
2193 {
2194 return std::get<double>(var_result);
2195 } else if (std::holds_alternative<int>(var_result))
2196 {
2197 return std::get<int>(var_result);
2198 } else if (std::holds_alternative<bool>(var_result))
2199 {
2200 return std::get<bool>(var_result);
2201 } else return Const::doubleNaN;
2202 };
2203 return func;
2204 } else {
2205 B2FATAL("Wrong number of arguments for meta function genUpsilon4S");
2206 }
2207 }
2208
2209 Manager::FunctionPtr getVariableByRank(const std::vector<std::string>& arguments)
2210 {
2211 if (arguments.size() == 4) {
2212 std::string listName = arguments[0];
2213 std::string rankedVariableName = arguments[1];
2214 std::string returnVariableName = arguments[2];
2215 std::string extraInfoName = rankedVariableName + "_rank";
2216 int rank = 1;
2217 try {
2218 rank = convertString<int>(arguments[3]);
2219 } catch (std::invalid_argument&) {
2220 B2ERROR("3rd argument of getVariableByRank meta function (Rank) must be an integer!");
2221 return nullptr;
2222 }
2223
2224 const Variable::Manager::Var* var = Manager::Instance().getVariable(returnVariableName);
2225 auto func = [var, rank, extraInfoName, listName](const Particle*)-> double {
2226 StoreObjPtr<ParticleList> list(listName);
2227
2228 const unsigned int numParticles = list->getListSize();
2229 for (unsigned int i = 0; i < numParticles; i++)
2230 {
2231 const Particle* p = list->getParticle(i);
2232 if (p->getExtraInfo(extraInfoName) == rank) {
2233 auto var_result = var->function(p);
2234 if (std::holds_alternative<double>(var_result)) {
2235 return std::get<double>(var_result);
2236 } else if (std::holds_alternative<int>(var_result)) {
2237 return std::get<int>(var_result);
2238 } else if (std::holds_alternative<bool>(var_result)) {
2239 return std::get<bool>(var_result);
2240 } else return Const::doubleNaN;
2241 }
2242 }
2243 // return 0;
2244 return std::numeric_limits<double>::signaling_NaN();
2245 };
2246 return func;
2247 } else {
2248 B2FATAL("Wrong number of arguments for meta function getVariableByRank");
2249 }
2250 }
2251
2252 Manager::FunctionPtr countInList(const std::vector<std::string>& arguments)
2253 {
2254 if (arguments.size() == 1 or arguments.size() == 2) {
2255
2256 std::string listName = arguments[0];
2257 std::string cutString = "";
2258
2259 if (arguments.size() == 2) {
2260 cutString = arguments[1];
2261 }
2262
2263 std::shared_ptr<Variable::Cut> cut = std::shared_ptr<Variable::Cut>(Variable::Cut::compile(cutString));
2264
2265 auto func = [listName, cut](const Particle*) -> int {
2266
2267 StoreObjPtr<ParticleList> list(listName);
2268 int sum = 0;
2269 for (unsigned int i = 0; i < list->getListSize(); i++)
2270 {
2271 const Particle* particle = list->getParticle(i);
2272 if (cut->check(particle)) {
2273 sum++;
2274 }
2275 }
2276 return sum;
2277 };
2278 return func;
2279 } else {
2280 B2FATAL("Wrong number of arguments for meta function countInList");
2281 }
2282 }
2283
2284 Manager::FunctionPtr veto(const std::vector<std::string>& arguments)
2285 {
2286 if (arguments.size() == 2 or arguments.size() == 3) {
2287
2288 std::string roeListName = arguments[0];
2289 std::string cutString = arguments[1];
2290 int pdgCode = Const::electron.getPDGCode();
2291 if (arguments.size() == 2) {
2292 B2INFO("Use pdgCode of electron as default in meta variable veto, other arguments: " << roeListName << ", " << cutString);
2293 } else {
2294 try {
2295 pdgCode = convertString<int>(arguments[2]);;
2296 } catch (std::invalid_argument&) {
2297 B2FATAL("Third argument of veto meta function must be integer!");
2298 }
2299 }
2300
2301 auto flavourType = (EvtPDLUtil::hasAntiParticle(pdgCode)) ? Particle::c_Flavored : Particle::c_Unflavored;
2302 std::shared_ptr<Variable::Cut> cut = std::shared_ptr<Variable::Cut>(Variable::Cut::compile(cutString));
2303
2304 auto func = [roeListName, cut, pdgCode, flavourType](const Particle * particle) -> bool {
2305 StoreObjPtr<ParticleList> roeList(roeListName);
2306 ROOT::Math::PxPyPzEVector vec = particle->get4Vector();
2307 for (unsigned int i = 0; i < roeList->getListSize(); i++)
2308 {
2309 const Particle* roeParticle = roeList->getParticle(i);
2310 if (not particle->overlapsWith(roeParticle)) {
2311 ROOT::Math::PxPyPzEVector tempCombination = roeParticle->get4Vector() + vec;
2312 std::vector<int> indices = { particle->getArrayIndex(), roeParticle->getArrayIndex() };
2313 Particle tempParticle = Particle(tempCombination, pdgCode, flavourType, indices, particle->getArrayPointer());
2314 if (cut->check(&tempParticle)) {
2315 return 1;
2316 }
2317 }
2318 }
2319 return 0;
2320 };
2321 return func;
2322 } else {
2323 B2FATAL("Wrong number of arguments for meta function veto");
2324 }
2325 }
2326
2327 Manager::FunctionPtr countDaughters(const std::vector<std::string>& arguments)
2328 {
2329 if (arguments.size() == 1) {
2330 std::string cutString = arguments[0];
2331 std::shared_ptr<Variable::Cut> cut = std::shared_ptr<Variable::Cut>(Variable::Cut::compile(cutString));
2332 auto func = [cut](const Particle * particle) -> int {
2333 int n = 0;
2334 for (auto& daughter : particle->getDaughters())
2335 {
2336 if (cut->check(daughter))
2337 ++n;
2338 }
2339 return n;
2340 };
2341 return func;
2342 } else {
2343 B2FATAL("Wrong number of arguments for meta function countDaughters");
2344 }
2345 }
2346
2347 Manager::FunctionPtr countFSPDaughters(const std::vector<std::string>& arguments)
2348 {
2349 if (arguments.size() == 1) {
2350 std::string cutString = arguments[0];
2351 std::shared_ptr<Variable::Cut> cut = std::shared_ptr<Variable::Cut>(Variable::Cut::compile(cutString));
2352 auto func = [cut](const Particle * particle) -> int {
2353
2354 std::vector<const Particle*> fspDaughters;
2355 particle->fillFSPDaughters(fspDaughters);
2356
2357 int n = 0;
2358 for (auto& daughter : fspDaughters)
2359 {
2360 if (cut->check(daughter))
2361 ++n;
2362 }
2363 return n;
2364 };
2365 return func;
2366 } else {
2367 B2FATAL("Wrong number of arguments for meta function countFSPDaughters");
2368 }
2369 }
2370
2371 Manager::FunctionPtr countDescendants(const std::vector<std::string>& arguments)
2372 {
2373 if (arguments.size() == 1) {
2374 std::string cutString = arguments[0];
2375 std::shared_ptr<Variable::Cut> cut = std::shared_ptr<Variable::Cut>(Variable::Cut::compile(cutString));
2376 auto func = [cut](const Particle * particle) -> int {
2377
2378 std::vector<const Particle*> allDaughters;
2379 particle->fillAllDaughters(allDaughters);
2380
2381 int n = 0;
2382 for (auto& daughter : allDaughters)
2383 {
2384 if (cut->check(daughter))
2385 ++n;
2386 }
2387 return n;
2388 };
2389 return func;
2390 } else {
2391 B2FATAL("Wrong number of arguments for meta function countDescendants");
2392 }
2393 }
2394
2395 Manager::FunctionPtr numberOfNonOverlappingParticles(const std::vector<std::string>& arguments)
2396 {
2397
2398 auto func = [arguments](const Particle * particle) -> int {
2399
2400 int _numberOfNonOverlappingParticles = 0;
2401 for (const auto& listName : arguments)
2402 {
2403 StoreObjPtr<ParticleList> list(listName);
2404 if (not list.isValid()) {
2405 B2FATAL("Invalid list named " << listName << " encountered in numberOfNonOverlappingParticles.");
2406 }
2407 for (unsigned int i = 0; i < list->getListSize(); i++) {
2408 const Particle* p = list->getParticle(i);
2409 if (not particle->overlapsWith(p)) {
2410 _numberOfNonOverlappingParticles++;
2411 }
2412 }
2413 }
2414 return _numberOfNonOverlappingParticles;
2415 };
2416
2417 return func;
2418
2419 }
2420
2421 void appendDaughtersRecursive(Particle* mother)
2422 {
2423
2424 auto* mcmother = mother->getRelated<MCParticle>();
2425
2426 if (!mcmother)
2427 return;
2428
2429 std::vector<MCParticle*> mcdaughters = mcmother->getDaughters();
2430 StoreArray<Particle> particles;
2431
2432 for (auto& mcdaughter : mcdaughters) {
2433 if (!mcdaughter->hasStatus(MCParticle::c_PrimaryParticle)) continue;
2434 Particle tmp_daughter(mcdaughter);
2435 Particle* new_daughter = particles.appendNew(tmp_daughter);
2436 new_daughter->addRelationTo(mcdaughter);
2437 mother->appendDaughter(new_daughter, false);
2438
2439 if (mcdaughter->getNDaughters() > 0)
2440 appendDaughtersRecursive(new_daughter);
2441 }
2442
2443 }
2444
2445 Manager::FunctionPtr matchedMC(const std::vector<std::string>& arguments)
2446 {
2447 if (arguments.size() == 1) {
2448 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
2449 auto func = [var](const Particle * particle) -> double {
2450 const MCParticle* mcp = particle->getMCParticle();
2451 if (!mcp) // Has no MC match and is no MCParticle
2452 {
2453 return Const::doubleNaN;
2454 }
2455 Particle tmpPart(mcp);
2456 StoreArray<Particle> particles;
2457 Particle* newPart = particles.appendNew(tmpPart);
2458 newPart->addRelationTo(mcp);
2459
2460 appendDaughtersRecursive(newPart);
2461
2462 auto var_result = var->function(newPart);
2463 if (std::holds_alternative<double>(var_result))
2464 {
2465 return std::get<double>(var_result);
2466 } else if (std::holds_alternative<int>(var_result))
2467 {
2468 return std::get<int>(var_result);
2469 } else if (std::holds_alternative<bool>(var_result))
2470 {
2471 return std::get<bool>(var_result);
2472 } else return Const::doubleNaN;
2473 };
2474 return func;
2475 } else {
2476 B2FATAL("Wrong number of arguments for meta function matchedMC");
2477 }
2478 }
2479
2480 Manager::FunctionPtr clusterBestMatchedMCParticle(const std::vector<std::string>& arguments)
2481 {
2482 if (arguments.size() == 1) {
2483 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
2484
2485 auto func = [var](const Particle * particle) -> double {
2486
2487 const ECLCluster* cluster = particle->getECLCluster();
2488 if (!cluster) return Const::doubleNaN;
2489
2490 auto mcps = cluster->getRelationsTo<MCParticle>();
2491 if (mcps.size() == 0) return Const::doubleNaN;
2492
2493 std::vector<std::pair<double, int>> weightsAndIndices;
2494 for (unsigned int i = 0; i < mcps.size(); ++i)
2495 weightsAndIndices.emplace_back(mcps.weight(i), i);
2496
2497 // sort descending by weight
2498 std::sort(weightsAndIndices.begin(), weightsAndIndices.end(),
2499 ValueIndexPairSorting::higherPair<decltype(weightsAndIndices)::value_type>);
2500
2501 // cppcheck-suppress containerOutOfBounds
2502 const MCParticle* mcp = mcps.object(weightsAndIndices[0].second);
2503
2504 Particle tmpPart(mcp);
2505 StoreArray<Particle> particles;
2506 Particle* newPart = particles.appendNew(tmpPart);
2507 newPart->addRelationTo(mcp);
2508
2509 appendDaughtersRecursive(newPart);
2510
2511 auto var_result = var->function(newPart);
2512 if (std::holds_alternative<double>(var_result))
2513 {
2514 return std::get<double>(var_result);
2515 } else if (std::holds_alternative<int>(var_result))
2516 {
2517 return std::get<int>(var_result);
2518 } else if (std::holds_alternative<bool>(var_result))
2519 {
2520 return std::get<bool>(var_result);
2521 } else
2522 {
2523 return Const::doubleNaN;
2524 }
2525 };
2526
2527 return func;
2528 } else {
2529 B2FATAL("Wrong number of arguments for meta function clusterBestMatchedMCParticle");
2530 }
2531 }
2532
2533 Manager::FunctionPtr clusterBestMatchedMCKlong(const std::vector<std::string>& arguments)
2534 {
2535 if (arguments.size() == 1) {
2536 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
2537
2538 auto func = [var](const Particle * particle) -> double {
2539
2540 const ECLCluster* cluster = particle->getECLCluster();
2541 if (!cluster) return Const::doubleNaN;
2542
2543 auto mcps = cluster->getRelationsTo<MCParticle>();
2544 if (mcps.size() == 0) return Const::doubleNaN;
2545
2546 std::map<int, double> mapMCParticleIndxAndWeight;
2547 getKlongWeightMap(particle, mapMCParticleIndxAndWeight);
2548
2549 // Klong is not found
2550 if (mapMCParticleIndxAndWeight.size() == 0)
2551 return Const::doubleNaN;
2552
2553 // find max totalWeight
2554 auto maxMap = std::max_element(mapMCParticleIndxAndWeight.begin(), mapMCParticleIndxAndWeight.end(),
2555 [](const auto & x, const auto & y) { return x.second < y.second; }
2556 );
2557
2558 StoreArray<MCParticle> mcparticles;
2559 const MCParticle* mcKlong = mcparticles[maxMap->first];
2560
2561 Particle tmpPart(mcKlong);
2562 auto var_result = var->function(&tmpPart);
2563 if (std::holds_alternative<double>(var_result))
2564 {
2565 return std::get<double>(var_result);
2566 } else if (std::holds_alternative<int>(var_result))
2567 {
2568 return std::get<int>(var_result);
2569 } else if (std::holds_alternative<bool>(var_result))
2570 {
2571 return std::get<bool>(var_result);
2572 } else
2573 {
2574 return Const::doubleNaN;
2575 }
2576 };
2577
2578 return func;
2579 } else {
2580 B2FATAL("Wrong number of arguments for meta function clusterBestMatchedMCKlong");
2581 }
2582 }
2583
2584 double matchedMCHasPDG(const Particle* particle, const std::vector<double>& pdgCode)
2585 {
2586 if (pdgCode.size() != 1) {
2587 B2FATAL("Too many arguments provided to matchedMCHasPDG!");
2588 }
2589 int inputPDG = std::lround(pdgCode[0]);
2590
2591 const MCParticle* mcp = particle->getMCParticle();
2592 if (!mcp)
2593 return Const::doubleNaN;
2594
2595 return std::abs(mcp->getPDG()) == inputPDG;
2596 }
2597
2598 Manager::FunctionPtr totalEnergyOfParticlesInList(const std::vector<std::string>& arguments)
2599 {
2600 if (arguments.size() == 1) {
2601 std::string listName = arguments[0];
2602 auto func = [listName](const Particle * particle) -> double {
2603
2604 (void) particle;
2605 StoreObjPtr<ParticleList> listOfParticles(listName);
2606
2607 if (!(listOfParticles.isValid())) B2FATAL("Invalid Listname " << listName << " given to totalEnergyOfParticlesInList");
2608 double totalEnergy = 0;
2609 int nParticles = listOfParticles->getListSize();
2610 for (int i = 0; i < nParticles; i++)
2611 {
2612 const Particle* part = listOfParticles->getParticle(i);
2613 const auto& frame = ReferenceFrame::GetCurrent();
2614 totalEnergy += frame.getMomentum(part).E();
2615 }
2616 return totalEnergy;
2617
2618 };
2619 return func;
2620 } else {
2621 B2FATAL("Wrong number of arguments for meta function totalEnergyOfParticlesInList");
2622 }
2623 }
2624
2625 Manager::FunctionPtr totalPxOfParticlesInList(const std::vector<std::string>& arguments)
2626 {
2627 if (arguments.size() == 1) {
2628 std::string listName = arguments[0];
2629 auto func = [listName](const Particle*) -> double {
2630 StoreObjPtr<ParticleList> listOfParticles(listName);
2631
2632 if (!(listOfParticles.isValid())) B2FATAL("Invalid Listname " << listName << " given to totalPxOfParticlesInList");
2633 double totalPx = 0;
2634 int nParticles = listOfParticles->getListSize();
2635 const auto& frame = ReferenceFrame::GetCurrent();
2636 for (int i = 0; i < nParticles; i++)
2637 {
2638 const Particle* part = listOfParticles->getParticle(i);
2639 totalPx += frame.getMomentum(part).Px();
2640 }
2641 return totalPx;
2642 };
2643 return func;
2644 } else {
2645 B2FATAL("Wrong number of arguments for meta function totalPxOfParticlesInList");
2646 }
2647 }
2648
2649 Manager::FunctionPtr totalPyOfParticlesInList(const std::vector<std::string>& arguments)
2650 {
2651 if (arguments.size() == 1) {
2652 std::string listName = arguments[0];
2653 auto func = [listName](const Particle*) -> double {
2654 StoreObjPtr<ParticleList> listOfParticles(listName);
2655
2656 if (!(listOfParticles.isValid())) B2FATAL("Invalid Listname " << listName << " given to totalPyOfParticlesInList");
2657 double totalPy = 0;
2658 int nParticles = listOfParticles->getListSize();
2659 const auto& frame = ReferenceFrame::GetCurrent();
2660 for (int i = 0; i < nParticles; i++)
2661 {
2662 const Particle* part = listOfParticles->getParticle(i);
2663 totalPy += frame.getMomentum(part).Py();
2664 }
2665 return totalPy;
2666 };
2667 return func;
2668 } else {
2669 B2FATAL("Wrong number of arguments for meta function totalPyOfParticlesInList");
2670 }
2671 }
2672
2673 Manager::FunctionPtr totalPzOfParticlesInList(const std::vector<std::string>& arguments)
2674 {
2675 if (arguments.size() == 1) {
2676 std::string listName = arguments[0];
2677 auto func = [listName](const Particle*) -> double {
2678 StoreObjPtr<ParticleList> listOfParticles(listName);
2679
2680 if (!(listOfParticles.isValid())) B2FATAL("Invalid Listname " << listName << " given to totalPzOfParticlesInList");
2681 double totalPz = 0;
2682 int nParticles = listOfParticles->getListSize();
2683 const auto& frame = ReferenceFrame::GetCurrent();
2684 for (int i = 0; i < nParticles; i++)
2685 {
2686 const Particle* part = listOfParticles->getParticle(i);
2687 totalPz += frame.getMomentum(part).Pz();
2688 }
2689 return totalPz;
2690 };
2691 return func;
2692 } else {
2693 B2FATAL("Wrong number of arguments for meta function totalPzOfParticlesInList");
2694 }
2695 }
2696
2697 Manager::FunctionPtr invMassInLists(const std::vector<std::string>& arguments)
2698 {
2699 if (arguments.size() > 0) {
2700
2701 auto func = [arguments](const Particle * particle) -> double {
2702
2703 ROOT::Math::PxPyPzEVector total4Vector;
2704 // To make sure particles in particlesList don't overlap.
2705 std::vector<Particle*> particlePool;
2706
2707 (void) particle;
2708 for (const auto& argument : arguments)
2709 {
2710 StoreObjPtr <ParticleList> listOfParticles(argument);
2711
2712 if (!(listOfParticles.isValid())) B2FATAL("Invalid Listname " << argument << " given to invMassInLists");
2713 int nParticles = listOfParticles->getListSize();
2714 for (int i = 0; i < nParticles; i++) {
2715 bool overlaps = false;
2716 Particle* part = listOfParticles->getParticle(i);
2717 for (auto poolPart : particlePool) {
2718 if (part->overlapsWith(poolPart)) {
2719 overlaps = true;
2720 break;
2721 }
2722 }
2723 if (!overlaps) {
2724 total4Vector += part->get4Vector();
2725 particlePool.push_back(part);
2726 }
2727 }
2728 }
2729 double invariantMass = total4Vector.M();
2730 return invariantMass;
2731
2732 };
2733 return func;
2734 } else {
2735 B2FATAL("Wrong number of arguments for meta function invMassInLists");
2736 }
2737 }
2738
2739 Manager::FunctionPtr totalECLEnergyOfParticlesInList(const std::vector<std::string>& arguments)
2740 {
2741 if (arguments.size() == 1) {
2742 std::string listName = arguments[0];
2743 auto func = [listName](const Particle * particle) -> double {
2744
2745 (void) particle;
2746 StoreObjPtr<ParticleList> listOfParticles(listName);
2747
2748 if (!(listOfParticles.isValid())) B2FATAL("Invalid Listname " << listName << " given to totalEnergyOfParticlesInList");
2749 double totalEnergy = 0;
2750 int nParticles = listOfParticles->getListSize();
2751 for (int i = 0; i < nParticles; i++)
2752 {
2753 const Particle* part = listOfParticles->getParticle(i);
2754 const ECLCluster* cluster = part->getECLCluster();
2755 const ECLCluster::EHypothesisBit clusterHypothesis = part->getECLClusterEHypothesisBit();
2756 if (cluster != nullptr) {
2757 totalEnergy += cluster->getEnergy(clusterHypothesis);
2758 }
2759 }
2760 return totalEnergy;
2761
2762 };
2763 return func;
2764 } else {
2765 B2FATAL("Wrong number of arguments for meta function totalECLEnergyOfParticlesInList");
2766 }
2767 }
2768
2769 Manager::FunctionPtr maxPtInList(const std::vector<std::string>& arguments)
2770 {
2771 if (arguments.size() == 1) {
2772 std::string listName = arguments[0];
2773 auto func = [listName](const Particle*) -> double {
2774 StoreObjPtr<ParticleList> listOfParticles(listName);
2775
2776 if (!(listOfParticles.isValid())) B2FATAL("Invalid Listname " << listName << " given to maxPtInList");
2777 int nParticles = listOfParticles->getListSize();
2778 const auto& frame = ReferenceFrame::GetCurrent();
2779 double maxPt = 0;
2780 for (int i = 0; i < nParticles; i++)
2781 {
2782 const Particle* part = listOfParticles->getParticle(i);
2783 const double Pt = frame.getMomentum(part).Pt();
2784 if (Pt > maxPt) maxPt = Pt;
2785 }
2786 return maxPt;
2787 };
2788 return func;
2789 } else {
2790 B2FATAL("Wrong number of arguments for meta function maxPtInList");
2791 }
2792 }
2793
2794 Manager::FunctionPtr eclClusterTrackMatchedWithCondition(const std::vector<std::string>& arguments)
2795 {
2796 if (arguments.size() <= 1) {
2797
2798 std::string cutString;
2799 if (arguments.size() == 1)
2800 cutString = arguments[0];
2801 std::shared_ptr<Variable::Cut> cut = std::shared_ptr<Variable::Cut>(Variable::Cut::compile(cutString));
2802 auto func = [cut](const Particle * particle) -> double {
2803
2804 if (particle == nullptr)
2805 return Const::doubleNaN;
2806
2807 const ECLCluster* cluster = particle->getECLCluster();
2808
2809 if (cluster)
2810 {
2811 auto tracks = cluster->getRelationsFrom<Track>();
2812
2813 for (const auto& track : tracks) {
2814 Particle trackParticle(&track, Const::pion);
2815
2816 if (cut->check(&trackParticle))
2817 return 1;
2818 }
2819 return 0;
2820 }
2821 return Const::doubleNaN;
2822 };
2823 return func;
2824 } else {
2825 B2FATAL("Wrong number of arguments for meta function eclClusterSpecialTrackMatched");
2826 }
2827 }
2828
2829 Manager::FunctionPtr averageValueInList(const std::vector<std::string>& arguments)
2830 {
2831 if (arguments.size() == 2) {
2832 std::string listName = arguments[0];
2833 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[1]);
2834
2835 auto func = [listName, var](const Particle*) -> double {
2836 StoreObjPtr<ParticleList> listOfParticles(listName);
2837
2838 if (!(listOfParticles.isValid())) B2FATAL("Invalid list name " << listName << " given to averageValueInList");
2839 int nParticles = listOfParticles->getListSize();
2840 if (nParticles == 0)
2841 {
2842 return Const::doubleNaN;
2843 }
2844 double average = 0;
2845 if (std::holds_alternative<double>(var->function(listOfParticles->getParticle(0))))
2846 {
2847 for (int i = 0; i < nParticles; i++) {
2848 average += std::get<double>(var->function(listOfParticles->getParticle(i))) / nParticles;
2849 }
2850 } else if (std::holds_alternative<int>(var->function(listOfParticles->getParticle(0))))
2851 {
2852 for (int i = 0; i < nParticles; i++) {
2853 average += std::get<int>(var->function(listOfParticles->getParticle(i))) / nParticles;
2854 }
2855 } else return Const::doubleNaN;
2856 return average;
2857 };
2858 return func;
2859 } else {
2860 B2FATAL("Wrong number of arguments for meta function averageValueInList");
2861 }
2862 }
2863
2864 Manager::FunctionPtr medianValueInList(const std::vector<std::string>& arguments)
2865 {
2866 if (arguments.size() == 2) {
2867 std::string listName = arguments[0];
2868 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[1]);
2869
2870 auto func = [listName, var](const Particle*) -> double {
2871 StoreObjPtr<ParticleList> listOfParticles(listName);
2872
2873 if (!(listOfParticles.isValid())) B2FATAL("Invalid list name " << listName << " given to medianValueInList");
2874 int nParticles = listOfParticles->getListSize();
2875 if (nParticles == 0)
2876 {
2877 return Const::doubleNaN;
2878 }
2879 std::vector<double> valuesInList;
2880 if (std::holds_alternative<double>(var->function(listOfParticles->getParticle(0))))
2881 {
2882 for (int i = 0; i < nParticles; i++) {
2883 valuesInList.push_back(std::get<double>(var->function(listOfParticles->getParticle(i))));
2884 }
2885 } else if (std::holds_alternative<int>(var->function(listOfParticles->getParticle(0))))
2886 {
2887 for (int i = 0; i < nParticles; i++) {
2888 valuesInList.push_back(std::get<int>(var->function(listOfParticles->getParticle(i))));
2889 }
2890 } else return Const::doubleNaN;
2891 std::sort(valuesInList.begin(), valuesInList.end());
2892 if (nParticles % 2 != 0)
2893 {
2894 return valuesInList[nParticles / 2];
2895 } else
2896 {
2897 return 0.5 * (valuesInList[nParticles / 2] + valuesInList[nParticles / 2 - 1]);
2898 }
2899 };
2900 return func;
2901 } else {
2902 B2FATAL("Wrong number of arguments for meta function medianValueInList");
2903 }
2904 }
2905
2906 Manager::FunctionPtr sumValueInList(const std::vector<std::string>& arguments)
2907 {
2908 if (arguments.size() == 2) {
2909 std::string listName = arguments[0];
2910 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[1]);
2911
2912 auto func = [listName, var](const Particle*) -> double {
2913 StoreObjPtr<ParticleList> listOfParticles(listName);
2914
2915 if (!(listOfParticles.isValid())) B2FATAL("Invalid list name " << listName << " given to sumValueInList");
2916 int nParticles = listOfParticles->getListSize();
2917 if (nParticles == 0)
2918 {
2919 return Const::doubleNaN;
2920 }
2921 double sum = 0;
2922 if (std::holds_alternative<double>(var->function(listOfParticles->getParticle(0))))
2923 {
2924 for (int i = 0; i < nParticles; i++) {
2925 sum += std::get<double>(var->function(listOfParticles->getParticle(i)));
2926 }
2927 } else if (std::holds_alternative<int>(var->function(listOfParticles->getParticle(0))))
2928 {
2929 for (int i = 0; i < nParticles; i++) {
2930 sum += std::get<int>(var->function(listOfParticles->getParticle(i)));
2931 }
2932 } else return Const::doubleNaN;
2933 return sum;
2934 };
2935 return func;
2936 } else {
2937 B2FATAL("Wrong number of arguments for meta function sumValueInList");
2938 }
2939 }
2940
2941 Manager::FunctionPtr productValueInList(const std::vector<std::string>& arguments)
2942 {
2943 if (arguments.size() == 2) {
2944 std::string listName = arguments[0];
2945 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[1]);
2946
2947 auto func = [listName, var](const Particle*) -> double {
2948 StoreObjPtr<ParticleList> listOfParticles(listName);
2949
2950 if (!(listOfParticles.isValid())) B2FATAL("Invalid list name " << listName << " given to productValueInList");
2951 int nParticles = listOfParticles->getListSize();
2952 if (nParticles == 0)
2953 {
2954 return Const::doubleNaN;
2955 }
2956 double product = 1;
2957 if (std::holds_alternative<double>(var->function(listOfParticles->getParticle(0))))
2958 {
2959 for (int i = 0; i < nParticles; i++) {
2960 product *= std::get<double>(var->function(listOfParticles->getParticle(i)));
2961 }
2962 } else if (std::holds_alternative<int>(var->function(listOfParticles->getParticle(0))))
2963 {
2964 for (int i = 0; i < nParticles; i++) {
2965 product *= std::get<int>(var->function(listOfParticles->getParticle(i)));
2966 }
2967 } else return Const::doubleNaN;
2968 return product;
2969 };
2970 return func;
2971 } else {
2972 B2FATAL("Wrong number of arguments for meta function productValueInList");
2973 }
2974 }
2975
2976 Manager::FunctionPtr angleToClosestInList(const std::vector<std::string>& arguments)
2977 {
2978 // expecting the list name
2979 if (arguments.size() != 1)
2980 B2FATAL("Wrong number of arguments for meta function angleToClosestInList");
2981
2982 std::string listname = arguments[0];
2983
2984 auto func = [listname](const Particle * particle) -> double {
2985 // get the list and check it's valid
2986 StoreObjPtr<ParticleList> list(listname);
2987 if (not list.isValid())
2988 B2FATAL("Invalid particle list name " << listname << " given to angleToClosestInList");
2989
2990 // check the list isn't empty
2991 if (list->getListSize() == 0)
2992 return Const::doubleNaN;
2993
2994 // respect the current frame and get the momentum of our input
2995 const auto& frame = ReferenceFrame::GetCurrent();
2996 const auto p_this = frame.getMomentum(particle);
2997
2998 // find the particle index with the smallest opening angle
2999 double minAngle = 2 * M_PI;
3000 for (unsigned int i = 0; i < list->getListSize(); ++i)
3001 {
3002 const Particle* compareme = list->getParticle(i);
3003 const auto p_compare = frame.getMomentum(compareme);
3004 double angle = ROOT::Math::VectorUtil::Angle(p_compare, p_this);
3005 if (minAngle > angle) minAngle = angle;
3006 }
3007 return minAngle;
3008 };
3009 return func;
3010 }
3011
3012 Manager::FunctionPtr closestInList(const std::vector<std::string>& arguments)
3013 {
3014 // expecting the list name and a variable name
3015 if (arguments.size() != 2)
3016 B2FATAL("Wrong number of arguments for meta function closestInList");
3017
3018 std::string listname = arguments[0];
3019
3020 // the requested variable and check it exists
3021 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[1]);
3022
3023 auto func = [listname, var](const Particle * particle) -> double {
3024 // get the list and check it's valid
3025 StoreObjPtr<ParticleList> list(listname);
3026 if (not list.isValid())
3027 B2FATAL("Invalid particle list name " << listname << " given to closestInList");
3028
3029 // respect the current frame and get the momentum of our input
3030 const auto& frame = ReferenceFrame::GetCurrent();
3031 const auto p_this = frame.getMomentum(particle);
3032
3033 // find the particle index with the smallest opening angle
3034 double minAngle = 2 * M_PI;
3035 int iClosest = -1;
3036 for (unsigned int i = 0; i < list->getListSize(); ++i)
3037 {
3038 const Particle* compareme = list->getParticle(i);
3039 const auto p_compare = frame.getMomentum(compareme);
3040 double angle = ROOT::Math::VectorUtil::Angle(p_compare, p_this);
3041 if (minAngle > angle) {
3042 minAngle = angle;
3043 iClosest = i;
3044 }
3045 }
3046
3047 // final check that the list wasn't empty (or some other problem)
3048 if (iClosest == -1) return Const::doubleNaN;
3049 auto var_result = var->function(list->getParticle(iClosest));
3050 if (std::holds_alternative<double>(var_result))
3051 {
3052 return std::get<double>(var_result);
3053 } else if (std::holds_alternative<int>(var_result))
3054 {
3055 return std::get<int>(var_result);
3056 } else if (std::holds_alternative<bool>(var_result))
3057 {
3058 return std::get<bool>(var_result);
3059 } else return Const::doubleNaN;
3060 };
3061 return func;
3062 }
3063
3064 Manager::FunctionPtr angleToMostB2BInList(const std::vector<std::string>& arguments)
3065 {
3066 // expecting the list name
3067 if (arguments.size() != 1)
3068 B2FATAL("Wrong number of arguments for meta function angleToMostB2BInList");
3069
3070 std::string listname = arguments[0];
3071
3072 auto func = [listname](const Particle * particle) -> double {
3073 // get the list and check it's valid
3074 StoreObjPtr<ParticleList> list(listname);
3075 if (not list.isValid())
3076 B2FATAL("Invalid particle list name " << listname << " given to angleToMostB2BInList");
3077
3078 // check the list isn't empty
3079 if (list->getListSize() == 0)
3080 return Const::doubleNaN;
3081
3082 // respect the current frame and get the momentum of our input
3083 const auto& frame = ReferenceFrame::GetCurrent();
3084 const auto p_this = frame.getMomentum(particle);
3085
3086 // find the most back-to-back (the largest opening angle before they
3087 // start getting smaller again!)
3088 double maxAngle = 0;
3089 for (unsigned int i = 0; i < list->getListSize(); ++i)
3090 {
3091 const Particle* compareme = list->getParticle(i);
3092 const auto p_compare = frame.getMomentum(compareme);
3093 double angle = ROOT::Math::VectorUtil::Angle(p_compare, p_this);
3094 if (maxAngle < angle) maxAngle = angle;
3095 }
3096 return maxAngle;
3097 };
3098 return func;
3099 }
3100
3101 Manager::FunctionPtr deltaPhiToMostB2BPhiInList(const std::vector<std::string>& arguments)
3102 {
3103 // expecting the list name
3104 if (arguments.size() != 1)
3105 B2FATAL("Wrong number of arguments for meta function deltaPhiToMostB2BPhiInList");
3106
3107 std::string listname = arguments[0];
3108
3109 auto func = [listname](const Particle * particle) -> double {
3110 // get the list and check it's valid
3111 StoreObjPtr<ParticleList> list(listname);
3112 if (not list.isValid())
3113 B2FATAL("Invalid particle list name " << listname << " given to deltaPhiToMostB2BPhiInList");
3114
3115 // check the list isn't empty
3116 if (list->getListSize() == 0)
3117 return Const::doubleNaN;
3118
3119 // respect the current frame and get the momentum of our input
3120 const auto& frame = ReferenceFrame::GetCurrent();
3121 const auto phi_this = frame.getMomentum(particle).Phi();
3122
3123 // find the most back-to-back in phi (largest absolute value of delta phi)
3124 double maxAngle = 0;
3125 for (unsigned int i = 0; i < list->getListSize(); ++i)
3126 {
3127 const Particle* compareme = list->getParticle(i);
3128 const auto phi_compare = frame.getMomentum(compareme).Phi();
3129 double angle = std::abs(phi_compare - phi_this);
3130 if (angle > M_PI) {angle = 2 * M_PI - angle;}
3131 if (maxAngle < angle) maxAngle = angle;
3132 }
3133 return maxAngle;
3134 };
3135 return func;
3136 }
3137
3138 Manager::FunctionPtr mostB2BInList(const std::vector<std::string>& arguments)
3139 {
3140 // expecting the list name and a variable name
3141 if (arguments.size() != 2)
3142 B2FATAL("Wrong number of arguments for meta function mostB2BInList");
3143
3144 std::string listname = arguments[0];
3145
3146 // the requested variable and check it exists
3147 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[1]);
3148
3149 auto func = [listname, var](const Particle * particle) -> double {
3150 // get the list and check it's valid
3151 StoreObjPtr<ParticleList> list(listname);
3152 if (not list.isValid())
3153 B2FATAL("Invalid particle list name " << listname << " given to mostB2BInList");
3154
3155 // respect the current frame and get the momentum of our input
3156 const auto& frame = ReferenceFrame::GetCurrent();
3157 const auto p_this = frame.getMomentum(particle);
3158
3159 // find the most back-to-back (the largest opening angle before they
3160 // start getting smaller again!)
3161 double maxAngle = -1.0;
3162 int iMostB2B = -1;
3163 for (unsigned int i = 0; i < list->getListSize(); ++i)
3164 {
3165 const Particle* compareme = list->getParticle(i);
3166 const auto p_compare = frame.getMomentum(compareme);
3167 double angle = ROOT::Math::VectorUtil::Angle(p_compare, p_this);
3168 if (maxAngle < angle) {
3169 maxAngle = angle;
3170 iMostB2B = i;
3171 }
3172 }
3173
3174 // final check that the list wasn't empty (or some other problem)
3175 if (iMostB2B == -1) return Const::doubleNaN;
3176 auto var_result = var->function(list->getParticle(iMostB2B));
3177 if (std::holds_alternative<double>(var_result))
3178 {
3179 return std::get<double>(var_result);
3180 } else if (std::holds_alternative<int>(var_result))
3181 {
3182 return std::get<int>(var_result);
3183 } else if (std::holds_alternative<bool>(var_result))
3184 {
3185 return std::get<bool>(var_result);
3186 } else return Const::doubleNaN;
3187 };
3188 return func;
3189 }
3190
3191 Manager::FunctionPtr maxOpeningAngleInList(const std::vector<std::string>& arguments)
3192 {
3193 if (arguments.size() == 1) {
3194 std::string listName = arguments[0];
3195 auto func = [listName](const Particle*) -> double {
3196 StoreObjPtr<ParticleList> listOfParticles(listName);
3197
3198 if (!(listOfParticles.isValid())) B2FATAL("Invalid Listname " << listName << " given to maxOpeningAngleInList");
3199 int nParticles = listOfParticles->getListSize();
3200 // return NaN if number of particles is less than 2
3201 if (nParticles < 2) return Const::doubleNaN;
3202
3203 const auto& frame = ReferenceFrame::GetCurrent();
3204 double maxOpeningAngle = -1;
3205 for (int i = 0; i < nParticles; i++)
3206 {
3207 ROOT::Math::PxPyPzEVector v1 = frame.getMomentum(listOfParticles->getParticle(i));
3208 for (int j = i + 1; j < nParticles; j++) {
3209 ROOT::Math::PxPyPzEVector v2 = frame.getMomentum(listOfParticles->getParticle(j));
3210 const double angle = ROOT::Math::VectorUtil::Angle(v1, v2);
3211 if (angle > maxOpeningAngle) maxOpeningAngle = angle;
3212 }
3213 }
3214 return maxOpeningAngle;
3215 };
3216 return func;
3217 } else {
3218 B2FATAL("Wrong number of arguments for meta function maxOpeningAngleInList");
3219 }
3220 }
3221
3222 Manager::FunctionPtr daughterCombination(const std::vector<std::string>& arguments)
3223 {
3224 // Expect 2 or more arguments.
3225 if (arguments.size() >= 2) {
3226 // First argument is the variable name
3227 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
3228
3229 // Core function: calculates a variable combining an arbitrary number of particles
3230 auto func = [var, arguments](const Particle * particle) -> double {
3231 if (particle == nullptr)
3232 {
3233 B2WARNING("Trying to access a daughter that does not exist. Skipping");
3234 return Const::doubleNaN;
3235 }
3236 const auto& frame = ReferenceFrame::GetCurrent();
3237
3238 // Sum of the 4-momenta of all the selected daughters
3239 ROOT::Math::PxPyPzEVector pSum(0, 0, 0, 0);
3240
3241 // Loop over the arguments. Each one of them is a generalizedIndex,
3242 // pointing to a particle in the decay tree.
3243 for (unsigned int iCoord = 1; iCoord < arguments.size(); iCoord++)
3244 {
3245 auto generalizedIndex = arguments[iCoord];
3246 const Particle* dauPart = particle->getParticleFromGeneralizedIndexString(generalizedIndex);
3247 if (dauPart)
3248 pSum += frame.getMomentum(dauPart);
3249 else {
3250 B2WARNING("Trying to access a daughter that does not exist. Index = " << generalizedIndex);
3251 return Const::doubleNaN;
3252 }
3253 }
3254
3255 // Make a dummy particle out of the sum of the 4-momenta of the selected daughters
3256 Particle sumOfDaughters(pSum, 100); // 100 is one of the special numbers
3257
3258 auto var_result = var->function(&sumOfDaughters);
3259 // Calculate the variable on the dummy particle
3260 if (std::holds_alternative<double>(var_result))
3261 {
3262 return std::get<double>(var_result);
3263 } else if (std::holds_alternative<int>(var_result))
3264 {
3265 return std::get<int>(var_result);
3266 } else if (std::holds_alternative<bool>(var_result))
3267 {
3268 return std::get<bool>(var_result);
3269 } else return Const::doubleNaN;
3270 };
3271 return func;
3272 } else
3273 B2FATAL("Wrong number of arguments for meta function daughterCombination");
3274 }
3275
3276 Manager::FunctionPtr useAlternativeDaughterHypothesis(const std::vector<std::string>& arguments)
3277 {
3278 /*
3279 `arguments` contains the variable to calculate and a list of colon-separated index-particle pairs.
3280 Overall, it looks like {"M", "0:K+", "1:p+", "3:e-"}.
3281 The code is thus divided in two parts:
3282 1) Parsing. A loop over the elements of `arguments` that first separates the variable from the rest, and then splits all the index:particle
3283 pairs, filling a std::vector with the indexes and another one with the new mass values.
3284 2) Replacing: A loop over the particle's daughters. We take the 4-momentum of each of them, recalculating it with a new mass if needed, and then we calculate
3285 the variable value using the sum of all the 4-momenta, both updated and non-updated ones.
3286 */
3287
3288 // Expect 2 or more arguments.
3289 if (arguments.size() >= 2) {
3290
3291 //----
3292 // 1) parsing
3293 //----
3294
3295 // First argument is the variable name
3296 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
3297
3298 // Parses the other arguments, which are in the form of index:particleName pairs,
3299 // and stores indexes and pdgs in std::unordered_map
3300 std::unordered_map<unsigned int, int> mapOfReplacedDaughters;
3301
3302 // Loop over the arguments to parse them
3303 for (unsigned int iCoord = 1; iCoord < arguments.size(); iCoord++) {
3304 auto replacedDauString = arguments[iCoord];
3305 // Split the string in index and new mass
3306 std::vector<std::string> indexAndMass;
3307 boost::split(indexAndMass, replacedDauString, boost::is_any_of(":"));
3308
3309 // Checks that the index:particleName pair is properly formatted.
3310 if (indexAndMass.size() > 2) {
3311 B2WARNING("The string indicating which daughter's mass should be replaced contains more than two elements separated by a colon. Perhaps you tried to pass a generalized index, which is not supported yet for this variable. The offending string is "
3312 << replacedDauString << ", while a correct syntax looks like 0:K+.");
3313 return nullptr;
3314 }
3315
3316 if (indexAndMass.size() < 2) {
3317 B2WARNING("The string indicating which daughter's mass should be replaced contains only one colon-separated element instead of two. The offending string is "
3318 << replacedDauString << ", while a correct syntax looks like 0:K+.");
3319 return nullptr;
3320 }
3321
3322 // indexAndMass[0] is the daughter index as string. Try to convert it
3323 int dauIndex = 0;
3324 try {
3325 dauIndex = convertString<int>(indexAndMass[0]);
3326 } catch (std::invalid_argument&) {
3327 B2FATAL("Found the string " << indexAndMass[0] << "instead of a daughter index.");
3328 }
3329
3330 // Determine PDG code corresponding to indexAndMass[1] using the particle names defined in evt.pdl
3331 TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(indexAndMass[1].c_str());
3332 if (!particlePDG) {
3333 B2WARNING("Particle not in evt.pdl file! " << indexAndMass[1]);
3334 return nullptr;
3335 }
3336
3337 // Stores the indexes and the pdgs in the map that will be passed to the lambda function
3338 int pdgCode = particlePDG->PdgCode();
3339 mapOfReplacedDaughters[dauIndex] = pdgCode;
3340 } // End of parsing
3341
3342 // Check the size of mapOfReplacedDaughters
3343 if (mapOfReplacedDaughters.size() != arguments.size() - 1)
3344 B2FATAL("Overlapped daughter's index is detected in the meta-variable useAlternativeDaughterHypothesis");
3345
3346 //----
3347 // 2) replacing
3348 //----
3349
3350 // Core function: creates a new particle from the original one changing
3351 // some of the daughters' masses
3352 auto func = [var, mapOfReplacedDaughters](const Particle * particle) -> double {
3353 if (particle == nullptr)
3354 {
3355 B2WARNING("Trying to access a particle that does not exist. Skipping");
3356 return Const::doubleNaN;
3357 }
3358
3359 const auto& frame = ReferenceFrame::GetCurrent();
3360
3361 // Create a dummy particle from the given particle to overwrite its kinematics
3362 Particle* dummy = ParticleCopy::copyParticle(particle);
3363
3364 // Sum of the 4-momenta of all the daughters with the new mass assumptions
3365 ROOT::Math::PxPyPzMVector pSum(0, 0, 0, 0);
3366
3367 for (unsigned int iDau = 0; iDau < particle->getNDaughters(); iDau++)
3368 {
3369 const Particle* dauPart = particle->getDaughter(iDau);
3370 if (not dauPart) {
3371 B2WARNING("Trying to access a daughter that does not exist. Index = " << iDau);
3372 return Const::doubleNaN;
3373 }
3374
3375 ROOT::Math::PxPyPzMVector dauMom = ROOT::Math::PxPyPzMVector(frame.getMomentum(dauPart));
3376
3377 int pdgCode;
3378 try {
3379 pdgCode = mapOfReplacedDaughters.at(iDau);
3380 } catch (std::out_of_range&) {
3381 // iDau is not in mapOfReplacedDaughters
3382 pSum += dauMom;
3383 continue;
3384 }
3385
3386 // overwrite the daughter's kinematics
3387 double p_x = dauMom.Px();
3388 double p_y = dauMom.Py();
3389 double p_z = dauMom.Pz();
3390 dauMom.SetCoordinates(p_x, p_y, p_z, TDatabasePDG::Instance()->GetParticle(pdgCode)->Mass());
3391 const_cast<Particle*>(dummy->getDaughter(iDau))->set4VectorDividingByMomentumScaling(ROOT::Math::PxPyPzEVector(dauMom));
3392
3393 // overwrite the daughter's pdg
3394 const int charge = dummy->getDaughter(iDau)->getCharge();
3395 if (TDatabasePDG::Instance()->GetParticle(pdgCode)->Charge() / 3.0 == charge)
3396 const_cast<Particle*>(dummy->getDaughter(iDau))->setPDGCode(pdgCode);
3397 else
3398 const_cast<Particle*>(dummy->getDaughter(iDau))->setPDGCode(-1 * pdgCode);
3399
3400 pSum += dauMom;
3401 } // End of loop over number of daughter
3402
3403 // overwrite the particle's kinematics
3404 dummy->set4Vector(ROOT::Math::PxPyPzEVector(pSum));
3405
3406 auto var_result = var->function(dummy);
3407
3408 // Calculate the variable on the dummy particle
3409 if (std::holds_alternative<double>(var_result))
3410 {
3411 return std::get<double>(var_result);
3412 } else if (std::holds_alternative<int>(var_result))
3413 {
3414 return std::get<int>(var_result);
3415 } else if (std::holds_alternative<bool>(var_result))
3416 {
3417 return std::get<bool>(var_result);
3418 } else return Const::doubleNaN;
3419 }; // end of lambda function
3420 return func;
3421 }// end of check on number of arguments
3422 else
3423 B2FATAL("Wrong number of arguments for meta function useAlternativeDaughterHypothesis");
3424 }
3425
3426 Manager::FunctionPtr varForFirstMCAncestorOfType(const std::vector<std::string>& arguments)
3427 {
3428 if (arguments.size() == 2) {
3429 int pdg_code = -1;
3430 std::string arg = arguments[0];
3431 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[1]);
3432 TParticlePDG* part = TDatabasePDG::Instance()->GetParticle(arg.c_str());
3433
3434 if (part != nullptr) {
3435 pdg_code = std::abs(part->PdgCode());
3436 } else {
3437 try {
3438 pdg_code = convertString<int>(arg);
3439 } catch (std::exception& e) {}
3440 }
3441
3442 if (pdg_code == -1) {
3443 B2FATAL("Ancestor " + arg + " is not recognised. Please provide valid PDG code or particle name.");
3444 }
3445
3446 auto func = [pdg_code, var](const Particle * particle) -> double {
3447 const Particle* p = particle;
3448
3449 int ancestor_level = std::get<double>(Manager::Instance().getVariable("hasAncestor(" + std::to_string(pdg_code) + ", 0)")->function(p));
3450 if ((ancestor_level <= 0) or (std::isnan(ancestor_level)))
3451 {
3452 return Const::doubleNaN;
3453 }
3454
3455 const MCParticle* i_p = p->getMCParticle();
3456
3457 for (int a = 0; a < ancestor_level ; a = a + 1)
3458 {
3459 i_p = i_p->getMother();
3460 }
3461
3462 Particle m_p(i_p);
3463 StoreArray<Particle> particles;
3464 Particle* newPart = particles.appendNew(m_p);
3465 newPart->addRelationTo(i_p);
3466
3467 appendDaughtersRecursive(newPart);
3468
3469 auto var_result = var->function(newPart);
3470 if (std::holds_alternative<double>(var_result))
3471 {
3472 return std::get<double>(var_result);
3473 } else if (std::holds_alternative<int>(var_result))
3474 {
3475 return std::get<int>(var_result);
3476 } else if (std::holds_alternative<bool>(var_result))
3477 {
3478 return std::get<bool>(var_result);
3479 } else return Const::doubleNaN;
3480 };
3481 return func;
3482 } else {
3483 B2FATAL("Wrong number of arguments for meta function varForFirstMCAncestorOfType (expected 2: type and variable of interest)");
3484 }
3485 }
3486
3487 Manager::FunctionPtr nTrackFitResults(const std::vector<std::string>& arguments)
3488 {
3489 if (arguments.size() != 1) {
3490 B2FATAL("Number of arguments for nTrackFitResults must be 1, particleType or PDGcode");
3491 }
3492
3493 std::string arg = arguments[0];
3494 TParticlePDG* part = TDatabasePDG::Instance()->GetParticle(arg.c_str());
3495 int absPdg;
3496 if (part != nullptr) {
3497 absPdg = std::abs(part->PdgCode());
3498 } else {
3499 try {
3500 absPdg = convertString<int>(arg);
3501 } catch (std::exception& e) {}
3502 }
3503
3504 auto func = [absPdg](const Particle*) -> int {
3505
3506 Const::ChargedStable type(absPdg);
3507 StoreArray<Track> tracks;
3508
3509 int nTrackFitResults = 0;
3510
3511 for (const auto& track : tracks)
3512 {
3513 const TrackFitResult* trackFit = track.getTrackFitResultWithClosestMass(type);
3514
3515 if (!trackFit) continue;
3516 if (trackFit->getChargeSign() == 0) continue;
3517
3518 nTrackFitResults++;
3519 }
3520
3521 return nTrackFitResults;
3522
3523 };
3524 return func;
3525 }
3526
3527
3528 Manager::FunctionPtr convertToInt(const std::vector<std::string>& arguments)
3529 {
3530 if (arguments.size() == 2) {
3531 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
3532 int default_val = convertString<int>(arguments[1]);
3533 auto func = [var, default_val](const Particle * particle) -> int {
3534 auto var_result = var->function(particle);
3535 if (std::holds_alternative<double>(var_result))
3536 {
3537 double value = std::get<double>(var_result);
3538 if (value > std::numeric_limits<int>::max())
3539 value = std::numeric_limits<int>::max();
3540 if (value < std::numeric_limits<int>::min())
3541 value = std::numeric_limits<int>::min();
3542 if (std::isnan(value))
3543 value = default_val;
3544 return static_cast<int>(value);
3545 } else if (std::holds_alternative<int>(var_result))
3546 return std::get<int>(var_result);
3547 else if (std::holds_alternative<bool>(var_result))
3548 return static_cast<int>(std::get<bool>(var_result));
3549 else return default_val;
3550 };
3551 return func;
3552 } else {
3553 B2FATAL("Wrong number of arguments for meta function int, please provide variable name and replacement value for NaN!");
3554 }
3555 }
3556
3557 VARIABLE_GROUP("MetaFunctions");
3558 REGISTER_METAVARIABLE("nCleanedECLClusters(cut)", nCleanedECLClusters,
3559 "[Eventbased] Returns the number of clean Clusters in the event\n"
3560 "Clean clusters are defined by the clusters which pass the given cut assuming a photon hypothesis.",
3561 Manager::VariableDataType::c_int);
3562 REGISTER_METAVARIABLE("nCleanedTracks(cut)", nCleanedTracks,
3563 "[Eventbased] Returns the number of clean Tracks in the event\n"
3564 "Clean tracks are defined by the tracks which pass the given cut assuming a pion hypothesis.", Manager::VariableDataType::c_int);
3565 REGISTER_METAVARIABLE("formula(v1 + v2 * [v3 - v4] / v5^v6)", formula, R"DOCSTRING(
3566Returns the result of the given formula, where v1 to vN are variables or floating
3567point numbers. Currently the only supported operations are addition (``+``),
3568subtraction (``-``), multiplication (``*``), division (``/``) and power (``^``
3569or ``**``). Parenthesis can be in the form of square brackets ``[v1 * v2]``
3570or normal brackets ``(v1 * v2)``. It will work also with variables taking
3571arguments. Operator precedence is taken into account. For example ::
3572
3573 (daughter(0, E) + daughter(1, E))**2 - p**2 + 0.138
3574
3575.. versionchanged:: release-03-00-00
3576 now both, ``[]`` and ``()`` can be used for grouping operations, ``**`` can
3577 be used for exponent and float literals are possible directly in the
3578 formula.
3579)DOCSTRING", Manager::VariableDataType::c_double);
3580 REGISTER_METAVARIABLE("useRestFrame(variable)", useRestFrame,
3581 "Returns the value of the variable using the rest frame of the given particle as current reference frame.\n"
3582 "E.g. ``useRestFrame(daughter(0, p))`` returns the total momentum of the first daughter in its mother's rest-frame", Manager::VariableDataType::c_double);
3583 REGISTER_METAVARIABLE("useCMSFrame(variable)", useCMSFrame,
3584 "Returns the value of the variable using the CMS frame as current reference frame.\n"
3585 "E.g. ``useCMSFrame(E)`` returns the energy of a particle in the CMS frame.", Manager::VariableDataType::c_double);
3586 REGISTER_METAVARIABLE("useLabFrame(variable)", useLabFrame, R"DOC(
3587Returns the value of ``variable`` in the *lab* frame.
3588
3589.. tip::
3590 The lab frame is the default reference frame, usually you don't need to use this meta-variable.
3591 E.g. ``useLabFrame(E)`` returns the energy of a particle in the Lab frame, same as just ``E``.
3592
3593Specifying the lab frame is useful in some corner-cases. For example:
3594``useRestFrame(daughter(0, formula(E - useLabFrame(E))))`` which is the difference of the first daughter's energy in the rest frame of the mother (current particle) with the same daughter's lab-frame energy.
3595)DOC", Manager::VariableDataType::c_double);
3596 REGISTER_METAVARIABLE("useTagSideRecoilRestFrame(variable, daughterIndexTagB)", useTagSideRecoilRestFrame,
3597 "Returns the value of the variable in the rest frame of the recoiling particle to the tag side B meson.\n"
3598 "The variable should only be applied to an Upsilon(4S) list.\n"
3599 "E.g. ``useTagSideRecoilRestFrame(daughter(1, daughter(1, p)), 0)`` applied on a Upsilon(4S) list (``Upsilon(4S)->B+:tag B-:sig``) returns the momentum of the second daughter of the signal B meson in the signal B meson rest frame.", Manager::VariableDataType::c_double);
3600 REGISTER_METAVARIABLE("useParticleRestFrame(variable, particleList)", useParticleRestFrame,
3601 "Returns the value of the variable in the rest frame of the first Particle contained in the given ParticleList.\n"
3602 "It is strongly recommended to pass a ParticleList that contains at most only one Particle in each event. "
3603 "When more than one Particle is present in the ParticleList, only the first Particle in the list is used for "
3604 "computing the rest frame and a warning is thrown. If the given ParticleList is empty in an event, it returns NaN.", Manager::VariableDataType::c_double);
3605 REGISTER_METAVARIABLE("useRecoilParticleRestFrame(variable, particleList)", useRecoilParticleRestFrame,
3606 "Returns the value of the variable in the rest frame of recoil system against the first Particle contained in the given ParticleList.\n"
3607 "It is strongly recommended to pass a ParticleList that contains at most only one Particle in each event. "
3608 "When more than one Particle is present in the ParticleList, only the first Particle in the list is used for "
3609 "computing the rest frame and a warning is thrown. If the given ParticleList is empty in an event, it returns NaN.", Manager::VariableDataType::c_double);
3610 REGISTER_METAVARIABLE("useDaughterRestFrame(variable, daughterIndex_1, [daughterIndex_2, ... daughterIndex_3])", useDaughterRestFrame,
3611 "Returns the value of the variable in the rest frame of the selected daughter particle.\n"
3612 "The daughter is identified via generalized daughter index, e.g. ``0:1`` identifies the second daughter (1) "
3613 "of the first daughter (0). If the daughter index is invalid, it returns NaN.\n"
3614 "If two or more indices are given, the rest frame of the sum of the daughters is used.",
3615 Manager::VariableDataType::c_double);
3616 REGISTER_METAVARIABLE("useDaughterRecoilRestFrame(variable, daughterIndex_1, [daughterIndex_2, ... daughterIndex_3])", useDaughterRecoilRestFrame,
3617 "Returns the value of the variable in the rest frame of the recoil of the selected daughter particle.\n"
3618 "The daughter is identified via generalized daughter index, e.g. ``0:1`` identifies the second daughter (1) "
3619 "of the first daughter (0). If the daughter index is invalid, it returns NaN.\n"
3620 "If two or more indices are given, the rest frame of the sum of the daughters is used.",
3621 Manager::VariableDataType::c_double);
3622 REGISTER_METAVARIABLE("useMCancestorBRestFrame(variable)", useMCancestorBRestFrame,
3623 "Returns the value of the variable in the rest frame of the ancestor B MC particle.\n"
3624 "If no B or no MC-matching is found, it returns NaN.", Manager::VariableDataType::c_double);
3625 REGISTER_METAVARIABLE("passesCut(cut)", passesCut,
3626 "Returns 1 if particle passes the cut otherwise 0.\n"
3627 "Useful if you want to write out if a particle would have passed a cut or not.", Manager::VariableDataType::c_bool);
3628 REGISTER_METAVARIABLE("passesEventCut(cut)", passesEventCut,
3629 "[Eventbased] Returns 1 if event passes the cut otherwise 0.\n"
3630 "Useful if you want to select events passing a cut without looping into particles, such as for skimming.\n", Manager::VariableDataType::c_bool);
3631 REGISTER_METAVARIABLE("countDaughters(cut)", countDaughters,
3632 "Returns number of direct daughters which satisfy the cut.\n"
3633 "Used by the skimming package (for what exactly?)", Manager::VariableDataType::c_int);
3634 REGISTER_METAVARIABLE("countFSPDaughters(cut)", countDescendants,
3635 "Returns number of final-state daughters which satisfy the cut.",
3636 Manager::VariableDataType::c_int);
3637 REGISTER_METAVARIABLE("countDescendants(cut)", countDescendants,
3638 "Returns number of descendants for all generations which satisfy the cut.",
3639 Manager::VariableDataType::c_int);
3640 REGISTER_METAVARIABLE("varFor(pdgCode, variable)", varFor,
3641 "Returns the value of the variable for the given particle if its abs(pdgCode) agrees with the given one.\n"
3642 "E.g. ``varFor(11, p)`` returns the momentum if the particle is an electron or a positron.", Manager::VariableDataType::c_double);
3643 REGISTER_METAVARIABLE("varForMCGen(variable)", varForMCGen,
3644 "Returns the value of the variable for the given particle if the MC particle related to it is primary, not virtual, and not initial.\n"
3645 "If no MC particle is related to the given particle, or the MC particle is not primary, virtual, or initial, NaN will be returned.\n"
3646 "E.g. ``varForMCGen(PDG)`` returns the PDG code of the MC particle related to the given particle if it is primary, not virtual, and not initial.", Manager::VariableDataType::c_double);
3647 REGISTER_METAVARIABLE("nParticlesInList(particleListName)", nParticlesInList,
3648 "[Eventbased] Returns number of particles in the given particle List.", Manager::VariableDataType::c_int);
3649 REGISTER_METAVARIABLE("isInList(particleListName)", isInList,
3650 "Returns 1 if the particle is in the list provided, 0 if not. Note that this only checks the particle given. For daughters of composite particles, please see :b2:var:`isDaughterOfList`.", Manager::VariableDataType::c_bool);
3651 REGISTER_METAVARIABLE("isDaughterOfList(particleListNames)", isDaughterOfList,
3652 "Returns 1 if the given particle is a daughter of at least one of the particles in the given particle Lists.", Manager::VariableDataType::c_bool);
3653 REGISTER_METAVARIABLE("isDescendantOfList(particleListName[, anotherParticleListName][, generationFlag = -1])", isDescendantOfList, R"DOC(
3654 Returns 1 if the given particle appears in the decay chain of the particles in the given ParticleLists.
3655
3656 Passing an integer as the last argument, allows to check if the particle belongs to the specific generation:
3657
3658 * ``isDescendantOfList(<particle_list>,1)`` returns 1 if particle is a daughter of the list,
3659 * ``isDescendantOfList(<particle_list>,2)`` returns 1 if particle is a granddaughter of the list,
3660 * ``isDescendantOfList(<particle_list>,3)`` returns 1 if particle is a great-granddaughter of the list, etc.
3661 * Default value is ``-1`` that is inclusive for all generations.
3662 )DOC", Manager::VariableDataType::c_bool);
3663 REGISTER_METAVARIABLE("isMCDescendantOfList(particleListName[, anotherParticleListName][, generationFlag = -1])", isMCDescendantOfList, R"DOC(
3664 Returns 1 if the given particle is linked to the same MC particle as any reconstructed daughter of the decay lists.
3665
3666 Passing an integer as the last argument, allows to check if the particle belongs to the specific generation:
3667
3668 * ``isMCDescendantOfList(<particle_list>,1)`` returns 1 if particle is matched to the same particle as any daughter of the list,
3669 * ``isMCDescendantOfList(<particle_list>,2)`` returns 1 if particle is matched to the same particle as any granddaughter of the list,
3670 * ``isMCDescendantOfList(<particle_list>,3)`` returns 1 if particle is matched to the same particle as any great-granddaughter of the list, etc.
3671 * Default value is ``-1`` that is inclusive for all generations.
3672
3673 It makes only sense for lists created with `fillParticleListFromMC` function with ``addDaughters=True`` argument.
3674 )DOC", Manager::VariableDataType::c_bool);
3675
3676 REGISTER_METAVARIABLE("sourceObjectIsInList(particleListName)", sourceObjectIsInList, R"DOC(
3677Returns 1 if the underlying mdst object (e.g. track, or cluster) was used to create a particle in ``particleListName``, 0 if not.
3678
3679.. note::
3680 This only makes sense for particles that are not composite. Returns -1 for composite particles.
3681)DOC", Manager::VariableDataType::c_int);
3682
3683 REGISTER_METAVARIABLE("mcParticleIsInMCList(particleListName)", mcParticleIsInMCList, R"DOC(
3684Returns 1 if the particle's matched MC particle is also matched to a particle in ``particleListName``
3685(or if either of the lists were filled from generator level `modularAnalysis.fillParticleListFromMC`.)
3686
3687.. seealso:: :b2:var:`isMCDescendantOfList` to check daughters.
3688)DOC", Manager::VariableDataType::c_bool);
3689
3690 REGISTER_METAVARIABLE("isGrandDaughterOfList(particleListNames)", isGrandDaughterOfList,
3691 "Returns 1 if the given particle is a grand daughter of at least one of the particles in the given particle Lists.", Manager::VariableDataType::c_bool);
3692 REGISTER_METAVARIABLE("originalParticle(variable)", originalParticle, R"DOC(
3693 Returns value of variable for the original particle from which the given particle is copied.
3694
3695 The copy of particle is created, for example, when the vertex fit updates the daughters and `modularAnalysis.copyParticles` is called.
3696 Returns NaN if the given particle is not copied and so there is no original particle.
3697 )DOC", Manager::VariableDataType::c_double);
3698 REGISTER_METAVARIABLE("daughter(i, variable)", daughter, R"DOC(
3699 Returns value of variable for the i-th daughter. E.g.
3700
3701 * ``daughter(0, p)`` returns the total momentum of the first daughter.
3702 * ``daughter(0, daughter(1, p)`` returns the total momentum of the second daughter of the first daughter.
3703
3704 Returns NaN if particle is nullptr or if the given daughter-index is out of bound (>= amount of daughters).
3705 )DOC", Manager::VariableDataType::c_double);
3706 REGISTER_METAVARIABLE("originalDaughter(i, variable)", originalDaughter, R"DOC(
3707 Returns value of variable for the original particle from which the i-th daughter is copied.
3708
3709 The copy of particle is created, for example, when the vertex fit updates the daughters and `modularAnalysis.copyParticles` is called.
3710 Returns NaN if the daughter is not copied and so there is no original daughter.
3711
3712 Returns NaN if particle is nullptr or if the given daughter-index is out of bound (>= amount of daughters).
3713 )DOC", Manager::VariableDataType::c_double);
3714 REGISTER_METAVARIABLE("mcDaughter(i, variable)", mcDaughter, R"DOC(
3715 Returns the value of the requested variable for the i-th Monte Carlo daughter of the particle.
3716
3717 Returns NaN if the particle is nullptr, if the particle is not matched to an MC particle,
3718 or if the i-th MC daughter does not exist.
3719
3720 E.g. ``mcDaughter(0, PDG)`` will return the PDG code of the first MC daughter of the matched MC
3721 particle of the reconstructed particle the function is applied to.
3722
3723 The meta variable can also be nested: ``mcDaughter(0, mcDaughter(1, PDG))``.
3724 )DOC", Manager::VariableDataType::c_double);
3725 REGISTER_METAVARIABLE("mcMother(variable)", mcMother, R"DOC(
3726 Returns the value of the requested variable for the Monte Carlo mother of the particle.
3727
3728 Returns NaN if the particle is nullptr, if the particle is not matched to an MC particle,
3729 or if the MC mother does not exist.
3730
3731 E.g. ``mcMother(PDG)`` will return the PDG code of the MC mother of the matched MC
3732 particle of the reconstructed particle the function is applied to.
3733
3734 The meta variable can also be nested: ``mcMother(mcMother(PDG))``.
3735 )DOC", Manager::VariableDataType::c_double);
3736 REGISTER_METAVARIABLE("genParticle(index, variable)", genParticle, R"DOC(
3737[Eventbased] Returns the ``variable`` for the ith generator particle.
3738The arguments of the function must be the ``index`` of the particle in the MCParticle Array,
3739and ``variable``, the name of the function or variable for that generator particle.
3740If ``index`` goes beyond the length of the MCParticles array, NaN will be returned.
3741
3742E.g. ``genParticle(0, p)`` returns the total momentum of the first MCParticle, which in a generic decay up to MC15 is
3743the Upsilon(4S) and for MC16 and beyond the initial electron.
3744)DOC", Manager::VariableDataType::c_double);
3745 REGISTER_METAVARIABLE("genUpsilon4S(variable)", genUpsilon4S, R"DOC(
3746[Eventbased] Returns the ``variable`` evaluated for the generator-level :math:`\Upsilon(4S)`.
3747If no generator level :math:`\Upsilon(4S)` exists for the event, NaN will be returned.
3748
3749E.g. ``genUpsilon4S(p)`` returns the total momentum of the :math:`\Upsilon(4S)` in a generic decay.
3750``genUpsilon4S(mcDaughter(1, p))`` returns the total momentum of the second daughter of the
3751generator-level :math:`\Upsilon(4S)` (i.e. the momentum of the second B meson in a generic decay).
3752)DOC", Manager::VariableDataType::c_double);
3753 REGISTER_METAVARIABLE("daughterProductOf(variable)", daughterProductOf,
3754 "Returns product of a variable over all daughters.\n"
3755 "E.g. ``daughterProductOf(extraInfo(SignalProbability))`` returns the product of the SignalProbabilitys of all daughters.", Manager::VariableDataType::c_double);
3756 REGISTER_METAVARIABLE("daughterSumOf(variable)", daughterSumOf,
3757 "Returns sum of a variable over all daughters.\n"
3758 "E.g. ``daughterSumOf(nDaughters)`` returns the number of grand-daughters.", Manager::VariableDataType::c_double);
3759 REGISTER_METAVARIABLE("daughterLowest(variable)", daughterLowest,
3760 "Returns the lowest value of the given variable among all daughters.\n"
3761 "E.g. ``useCMSFrame(daughterLowest(p))`` returns the lowest momentum in CMS frame.", Manager::VariableDataType::c_double);
3762 REGISTER_METAVARIABLE("daughterHighest(variable)", daughterHighest,
3763 "Returns the highest value of the given variable among all daughters.\n"
3764 "E.g. ``useCMSFrame(daughterHighest(p))`` returns the highest momentum in CMS frame.", Manager::VariableDataType::c_double);
3765 REGISTER_METAVARIABLE("daughterDiffOf(daughterIndex_i, daughterIndex_j, variable)", daughterDiffOf, R"DOC(
3766 Returns the difference of a variable between the two given daughters.
3767 E.g. ``useRestFrame(daughterDiffOf(0, 1, p))`` returns the momentum difference between first and second daughter in the rest frame of the given particle.
3768 (That means that it returns :math:`p_j - p_i`)
3769
3770 The daughters can be provided as generalized daughter indexes, which are simply colon-separated
3771 lists of daughter indexes, ordered starting from the root particle. For example, ``0:1``
3772 identifies the second daughter (1) of the first daughter (0) of the mother particle.
3773
3774 )DOC", Manager::VariableDataType::c_double);
3775 REGISTER_METAVARIABLE("mcDaughterDiffOf(i, j, variable)", mcDaughterDiffOf,
3776 "MC matched version of the `daughterDiffOf` function.", Manager::VariableDataType::c_double);
3777 REGISTER_METAVARIABLE("grandDaughterDiffOf(i, j, variable)", grandDaughterDiffOf,
3778 "Returns the difference of a variable between the first daughters of the two given daughters.\n"
3779 "E.g. ``useRestFrame(grandDaughterDiffOf(0, 1, p))`` returns the momentum difference between the first daughters of the first and second daughter in the rest frame of the given particle.\n"
3780 "(That means that it returns :math:`p_j - p_i`)", Manager::VariableDataType::c_double);
3781 MAKE_DEPRECATED("grandDaughterDiffOf", false, "light-2402-ocicat", R"DOC(
3782 The difference between any combination of (grand-)daughters can be calculated with the more general variable :b2:var:`daughterDiffOf`
3783 by using generalized daughter indexes.)DOC");
3784 REGISTER_METAVARIABLE("daughterDiffOfPhi(i, j)", daughterDiffOfPhi,
3785 "Returns the difference in :math:`\\phi` between the two given daughters. The unit of the angle is ``rad``.\n"
3786 "The difference is signed and takes account of the ordering of the given daughters.\n"
3787 "The function returns :math:`\\phi_j - \\phi_i`.", Manager::VariableDataType::c_double);
3788 MAKE_DEPRECATED("daughterDiffOfPhi", true, "release-06-00-00", R"DOC(
3789 The difference of the azimuthal angle :math:`\\phi` of two daughters can be calculated with the generic variable :b2:var:`daughterDiffOf`.)DOC");
3790 REGISTER_METAVARIABLE("mcDaughterDiffOfPhi(i, j)", mcDaughterDiffOfPhi,
3791 "MC matched version of the `daughterDiffOfPhi` function. The unit of the angle is ``rad``", Manager::VariableDataType::c_double);
3792 MAKE_DEPRECATED("mcDaughterDiffOfPhi", true, "release-06-00-00", R"DOC(
3793 The difference of the azimuthal angle :math:`\\phi` of the MC partners of two daughters can be calculated with the generic variable :b2:var:`mcDaughterDiffOf`.)DOC");
3794 REGISTER_METAVARIABLE("grandDaughterDiffOfPhi(i, j)", grandDaughterDiffOfPhi,
3795 "Returns the difference in :math:`\\phi` between the first daughters of the two given daughters. The unit of the angle is ``rad``.\n"
3796 "The difference is signed and takes account of the ordering of the given daughters.\n"
3797 "The function returns :math:`\\phi_j - \\phi_i`.\n", Manager::VariableDataType::c_double);
3798 MAKE_DEPRECATED("grandDaughterDiffOfPhi", true, "release-06-00-00", R"DOC(
3799 The difference of the azimuthal angle :math:`\\phi` of two granddaughters can be calculated with the generic variable :b2:var:`grandDaughterDiffOf`.)DOC");
3800 REGISTER_METAVARIABLE("daughterDiffOfClusterPhi(i, j)", daughterDiffOfClusterPhi,
3801 "Returns the difference in :math:`\\phi` between the ECLClusters of two given daughters. The unit of the angle is ``rad``.\n"
3802 "The difference is signed and takes account of the ordering of the given daughters.\n"
3803 "The function returns :math:`\\phi_j - \\phi_i`.\n"
3804 "The function returns NaN if at least one of the daughters is not matched to or not based on an ECLCluster.", Manager::VariableDataType::c_double);
3805 MAKE_DEPRECATED("daughterDiffOfClusterPhi", true, "release-06-00-00", R"DOC(
3806 The difference of the azimuthal angle :math:`\\phi` of the related ECL clusters of two daughters can be calculated with the generic variable :b2:var:`daughterDiffOf`.)DOC");
3807 REGISTER_METAVARIABLE("grandDaughterDiffOfClusterPhi(i, j)", grandDaughterDiffOfClusterPhi,
3808 "Returns the difference in :math:`\\phi` between the ECLClusters of the daughters of the two given daughters. The unit of the angle is ``rad``.\n"
3809 "The difference is signed and takes account of the ordering of the given daughters.\n"
3810 "The function returns :math:`\\phi_j - \\phi_i`.\n"
3811 "The function returns NaN if at least one of the daughters is not matched to or not based on an ECLCluster.\n", Manager::VariableDataType::c_double);
3812 MAKE_DEPRECATED("grandDaughterDiffOfClusterPhi", true, "release-06-00-00", R"DOC(
3813 The difference of the azimuthal angle :math:`\\phi` of the related ECL clusters of two granddaughters can be calculated with the generic variable :b2:var:`grandDaughterDiffOf`.)DOC");
3814 REGISTER_METAVARIABLE("daughterDiffOfPhiCMS(i, j)", daughterDiffOfPhiCMS,
3815 "Returns the difference in :math:`\\phi` between the two given daughters in the CMS frame. The unit of the angle is ``rad``.\n"
3816 "The difference is signed and takes account of the ordering of the given daughters.\n"
3817 "The function returns :math:`\\phi_j - \\phi_i`.", Manager::VariableDataType::c_double);
3818 MAKE_DEPRECATED("daughterDiffOfPhiCMS", true, "release-06-00-00", R"DOC(
3819 The difference of the azimuthal angle :math:`\\phi` of two daughters in the CMS frame can be calculated with the generic variable :b2:var:`daughterDiffOf`.)DOC");
3820 REGISTER_METAVARIABLE("mcDaughterDiffOfPhiCMS(i, j)", daughterDiffOfPhiCMS,
3821 "MC matched version of the `daughterDiffOfPhiCMS` function. The unit of the angle is ``rad``", Manager::VariableDataType::c_double);
3822 MAKE_DEPRECATED("mcDaughterDiffOfPhiCMS", true, "release-06-00-00", R"DOC(
3823 The difference of the azimuthal angle :math:`\\phi` of the MC partners of two daughters in the CMS frame can be calculated with the generic variable :b2:var:`mcDaughterDiffOf`.)DOC");
3824 REGISTER_METAVARIABLE("daughterDiffOfClusterPhiCMS(i, j)", daughterDiffOfClusterPhiCMS,
3825 "Returns the difference in :math:`\\phi` between the ECLClusters of two given daughters in the CMS frame. The unit of the angle is ``rad``.\n"
3826 "The difference is signed and takes account of the ordering of the given daughters.\n"
3827 "The function returns :math:`\\phi_j - \\phi_i``.\n"
3828 "The function returns NaN if at least one of the daughters is not matched to or not based on an ECLCluster.", Manager::VariableDataType::c_double);
3829 MAKE_DEPRECATED("daughterDiffOfClusterPhiCMS", true, "release-06-00-00", R"DOC(
3830 The difference of the azimuthal angle :math:`\\phi` of the related ECL clusters of two daughters in the CMS frame can be calculated with the generic variable :b2:var:`daughterDiffOf`.)DOC");
3831 REGISTER_METAVARIABLE("daughterNormDiffOf(i, j, variable)", daughterNormDiffOf,
3832 "Returns the normalized difference of a variable between the two given daughters.\n"
3833 "E.g. ``daughterNormDiffOf(0, 1, p)`` returns the normalized momentum difference between first and second daughter in the lab frame.", Manager::VariableDataType::c_double);
3834 REGISTER_METAVARIABLE("daughterMotherDiffOf(i, variable)", daughterMotherDiffOf,
3835 "Returns the difference of a variable between the given daughter and the mother particle itself.\n"
3836 "E.g. ``useRestFrame(daughterMotherDiffOf(0, p))`` returns the momentum difference between the given particle and its first daughter in the rest frame of the mother.", Manager::VariableDataType::c_double);
3837 REGISTER_METAVARIABLE("daughterMotherNormDiffOf(i, variable)", daughterMotherNormDiffOf,
3838 "Returns the normalized difference of a variable between the given daughter and the mother particle itself.\n"
3839 "E.g. ``daughterMotherNormDiffOf(1, p)`` returns the normalized momentum difference between the given particle and its second daughter in the lab frame.", Manager::VariableDataType::c_double);
3840 REGISTER_METAVARIABLE("angleBetweenDaughterAndRecoil(daughterIndex_1, daughterIndex_2, ... )", angleBetweenDaughterAndRecoil, R"DOC(
3841 Returns the angle between the momentum recoiling against the particle and the sum of the momenta of the given daughters.
3842 The unit of the angle is ``rad``.
3843
3844 The particles are identified via generalized daughter indexes, which are simply colon-separated lists of
3845 daughter indexes, ordered starting from the root particle. For example, ``0:1:3`` identifies the fourth
3846 daughter (3) of the second daughter (1) of the first daughter (0) of the mother particle. ``1`` simply
3847 identifies the second daughter of the root particle.
3848
3849 At least one generalized index has to be given to ``angleBetweenDaughterAndRecoil``.
3850
3851 .. tip::
3852 ``angleBetweenDaughterAndRecoil(0)`` will return the angle between pRecoil and the momentum of the first daughter.
3853
3854 ``angleBetweenDaughterAndRecoil(0, 1)`` will return the angle between pRecoil and the sum of the momenta of the first and second daughter.
3855
3856 ``angleBetweenDaughterAndRecoil(0:0, 3:0)`` will return the angle between pRecoil and the sum of the momenta of the: first daughter of the first daughter, and
3857 the first daughter of the fourth daughter.)DOC", Manager::VariableDataType::c_double);
3858 REGISTER_METAVARIABLE("angleBetweenDaughterAndMissingMomentum(daughterIndex_1, daughterIndex_2, ... )", angleBetweenDaughterAndMissingMomentum, R"DOC(
3859 Returns the angle between the missing momentum in the event and the sum of the momenta of the given daughters.
3860 The unit of the angle is ``rad``. EventKinematics module has to be called to use this.
3861
3862 The particles are identified via generalized daughter indexes, which are simply colon-separated lists of
3863 daughter indexes, ordered starting from the root particle. For example, ``0:1:3`` identifies the fourth
3864 daughter (3) of the second daughter (1) of the first daughter (0) of the mother particle. ``1`` simply
3865 identifies the second daughter of the root particle.
3866
3867 At least one generalized index has to be given to ``angleBetweenDaughterAndMissingMomentum``.
3868
3869 .. tip::
3870 ``angleBetweenDaughterAndMissingMomentum(0)`` will return the angle between missMom and the momentum of the first daughter.
3871
3872 ``angleBetweenDaughterAndMissingMomentum(0, 1)`` will return the angle between missMom and the sum of the momenta of the first and second daughter.
3873
3874 ``angleBetweenDaughterAndMissingMomentum(0:0, 3:0)`` will return the angle between missMom and the sum of the momenta of the: first daughter of the first daughter, and
3875 the first daughter of the fourth daughter.)DOC", Manager::VariableDataType::c_double);
3876 REGISTER_METAVARIABLE("daughterAngle(daughterIndex_1, daughterIndex_2[, daughterIndex_3])", daughterAngle, R"DOC(
3877 Returns the angle in between any pair of particles belonging to the same decay tree.
3878 The unit of the angle is ``rad``.
3879
3880 The particles are identified via generalized daughter indexes, which are simply colon-separated lists of
3881 daughter indexes, ordered starting from the root particle. For example, ``0:1:3`` identifies the fourth
3882 daughter (3) of the second daughter (1) of the first daughter (0) of the mother particle. ``1`` simply
3883 identifies the second daughter of the root particle.
3884
3885 Both two and three generalized indexes can be given to ``daughterAngle``. If two indices are given, the
3886 variable returns the angle between the momenta of the two given particles. If three indices are given, the
3887 variable returns the angle between the momentum of the third particle and a vector which is the sum of the
3888 first two daughter momenta.
3889
3890 .. tip::
3891 ``daughterAngle(0, 3)`` will return the angle between the first and fourth daughter.
3892 ``daughterAngle(0, 1, 3)`` will return the angle between the fourth daughter and the sum of the first and
3893 second daughter.
3894 ``daughterAngle(0:0, 3:0)`` will return the angle between the first daughter of the first daughter, and
3895 the first daughter of the fourth daughter.
3896
3897 )DOC", Manager::VariableDataType::c_double);
3898 REGISTER_METAVARIABLE("mcDaughterAngle(daughterIndex_1, daughterIndex_2, [daughterIndex_3])", mcDaughterAngle,
3899 "MC matched version of the `daughterAngle` function. Also works if applied directly to MC particles. The unit of the angle is ``rad``", Manager::VariableDataType::c_double);
3900 REGISTER_VARIABLE("grandDaughterDecayAngle(i, j)", grandDaughterDecayAngle,
3901 "Returns the decay angle of the granddaughter in the daughter particle's rest frame.\n"
3902 "It is calculated with respect to the reverted momentum vector of the particle.\n"
3903 "Two arguments representing the daughter and granddaughter indices have to be provided as arguments.\n\n", "rad");
3904 REGISTER_VARIABLE("daughterClusterAngleInBetween(i, j)", daughterClusterAngleInBetween,
3905 "Returns the angle between clusters associated to the two daughters."
3906 "If two indices given: returns the angle between the momenta of the clusters associated to the two given daughters."
3907 "If three indices given: returns the angle between the momentum of the third particle's cluster and a vector "
3908 "which is the sum of the first two daughter's cluster momenta."
3909 "Returns nan if any of the daughters specified don't have an associated cluster."
3910 "The arguments in the argument vector must be integers corresponding to the ith and jth (and kth) daughters.\n\n", "rad");
3911 REGISTER_METAVARIABLE("daughterInvM(i[, j, ...])", daughterInvM, R"DOC(
3912 Returns the invariant mass adding the Lorentz vectors of the given daughters. The unit of the invariant mass is GeV/:math:`\text{c}^2`
3913 E.g. ``daughterInvM(0, 1, 2)`` returns the invariant Mass :math:`m = \sqrt{(p_0 + p_1 + p_2)^2}` of the first, second and third daughter.
3914
3915 Daughters from different generations of the decay tree can be combined using generalized daughter indexes,
3916 which are simply colon-separated daughter indexes for each generation, starting from the root particle. For
3917 example, ``0:1:3`` identifies the fourth daughter (3) of the second daughter (1) of the first daughter(0) of
3918 the mother particle.
3919
3920 Returns NaN if the given daughter-index is out of bound (>= number of daughters))DOC", Manager::VariableDataType::c_double);
3921 REGISTER_METAVARIABLE("extraInfo(name)", extraInfo,
3922 "Returns extra info stored under the given name.\n"
3923 "The extraInfo has to be set by a module first.\n"
3924 "E.g. ``extraInfo(SignalProbability)`` returns the SignalProbability calculated by the ``MVAExpert`` module.\n"
3925 "If nothing is set under the given name or if the particle is a nullptr, NaN is returned.\n"
3926 "In the latter case please use `eventExtraInfo` if you want to access an EventExtraInfo variable.", Manager::VariableDataType::c_double);
3927 REGISTER_METAVARIABLE("eventExtraInfo(name)", eventExtraInfo,
3928 "[Eventbased] Returns extra info stored under the given name in the event extra info.\n"
3929 "The extraInfo has to be set first by another module like MVAExpert in event mode.\n"
3930 "If nothing is set under this name, NaN is returned.", Manager::VariableDataType::c_double);
3931 REGISTER_METAVARIABLE("eventCached(variable)", eventCached,
3932 "[Eventbased] Returns value of event-based variable and caches this value in the EventExtraInfo.\n"
3933 "The result of second call to this variable in the same event will be provided from the cache.\n"
3934 "It is recommended to use this variable in order to declare custom aliases as event-based. This is "
3935 "necessary if using the eventwise mode of variablesToNtuple).", Manager::VariableDataType::c_double);
3936 REGISTER_METAVARIABLE("particleCached(variable)", particleCached,
3937 "Returns value of given variable and caches this value in the ParticleExtraInfo of the provided particle.\n"
3938 "The result of second call to this variable on the same particle will be provided from the cache.", Manager::VariableDataType::c_double);
3939 REGISTER_METAVARIABLE("modulo(variable, n)", modulo,
3940 "Returns rest of division of variable by n.", Manager::VariableDataType::c_int);
3941 REGISTER_METAVARIABLE("abs(variable)", abs,
3942 "Returns absolute value of the given variable.\n"
3943 "E.g. abs(mcPDG) returns the absolute value of the mcPDG, which is often useful for cuts.", Manager::VariableDataType::c_double);
3944 REGISTER_METAVARIABLE("max(var1,var2)", max, "Returns max value of two variables.\n", Manager::VariableDataType::c_double);
3945 REGISTER_METAVARIABLE("min(var1,var2)", min, "Returns min value of two variables.\n", Manager::VariableDataType::c_double);
3946 REGISTER_METAVARIABLE("sin(variable)", sin, "Returns sine value of the given variable.", Manager::VariableDataType::c_double);
3947 REGISTER_METAVARIABLE("asin(variable)", asin, "Returns arcsine of the given variable. The unit of the asin() is ``rad``", Manager::VariableDataType::c_double);
3948 REGISTER_METAVARIABLE("cos(variable)", cos, "Returns cosine value of the given variable.", Manager::VariableDataType::c_double);
3949 REGISTER_METAVARIABLE("acos(variable)", acos, "Returns arccosine value of the given variable. The unit of the acos() is ``rad``", Manager::VariableDataType::c_double);
3950 REGISTER_METAVARIABLE("tan(variable)", tan, "Returns tangent value of the given variable.", Manager::VariableDataType::c_double);
3951 REGISTER_METAVARIABLE("atan(variable)", atan, "Returns arctangent value of the given variable. The unit of the atan() is ``rad``", Manager::VariableDataType::c_double);
3952 REGISTER_METAVARIABLE("exp(variable)", exp, "Returns exponential evaluated for the given variable.", Manager::VariableDataType::c_double);
3953 REGISTER_METAVARIABLE("log(variable)", log, "Returns natural logarithm evaluated for the given variable.", Manager::VariableDataType::c_double);
3954 REGISTER_METAVARIABLE("log10(variable)", log10, "Returns base-10 logarithm evaluated for the given variable.", Manager::VariableDataType::c_double);
3955 REGISTER_METAVARIABLE("int(variable, nan_replacement)", convertToInt, R"DOC(
3956 Casts the output of the variable to an integer value.
3957
3958 .. note::
3959 Overflow and underflow are clipped at maximum and minimum values, respectively. NaN values are replaced with the value of the 2nd argument.
3960
3961 )DOC", Manager::VariableDataType::c_int);
3962 REGISTER_METAVARIABLE("isNAN(variable)", isNAN,
3963 "Returns true if variable value evaluates to nan (determined via std::isnan(double)).\n"
3964 "Useful for debugging.", Manager::VariableDataType::c_bool);
3965 REGISTER_METAVARIABLE("ifNANgiveX(variable, x)", ifNANgiveX,
3966 "Returns x (has to be a number) if variable value is nan (determined via std::isnan(double)).\n"
3967 "Useful for technical purposes while training MVAs.", Manager::VariableDataType::c_double);
3968 REGISTER_METAVARIABLE("isInfinity(variable)", isInfinity,
3969 "Returns true if variable value evaluates to infinity (determined via std::isinf(double)).\n"
3970 "Useful for debugging.", Manager::VariableDataType::c_bool);
3971 REGISTER_METAVARIABLE("unmask(variable, flag1, flag2, ...)", unmask,
3972 "unmask(variable, flag1, flag2, ...) or unmask(variable, mask) sets certain bits in the variable to zero.\n"
3973 "For example, if you want to set the second, fourth and fifth bits to zero, you could call \n"
3974 "``unmask(variable, 2, 8, 16)`` or ``unmask(variable, 26)``.\n"
3975 "", Manager::VariableDataType::c_double);
3976 REGISTER_METAVARIABLE("conditionalVariableSelector(cut, variableIfTrue, variableIfFalse)", conditionalVariableSelector,
3977 "Returns one of the two supplied variables, depending on whether the particle passes the supplied cut.\n"
3978 "The first variable is returned if the particle passes the cut, and the second variable is returned otherwise.", Manager::VariableDataType::c_double);
3979 REGISTER_METAVARIABLE("pValueCombination(p1, p2, ...)", pValueCombination,
3980 "Returns the combined p-value of the provided p-values according to the formula given in `Nucl. Instr. and Meth. A 411 (1998) 449 <https://doi.org/10.1016/S0168-9002(98)00293-9>`_ .\n"
3981 "If any of the p-values is invalid, i.e. smaller than zero, -1 is returned.", Manager::VariableDataType::c_double);
3982 REGISTER_METAVARIABLE("pValueCombinationOfDaughters(variable)", pValueCombinationOfDaughters,
3983 "Returns the combined p-value of the daughter p-values according to the formula given in `Nucl. Instr. and Meth. A 411 (1998) 449 <https://doi.org/10.1016/S0168-9002(98)00293-9>`_ .\n"
3984 "If any of the p-values is invalid, i.e. smaller than zero, -1 is returned.", Manager::VariableDataType::c_double);
3985 REGISTER_METAVARIABLE("veto(particleList, cut, pdgCode = 11)", veto,
3986 "Combines current particle with particles from the given particle list and returns 1 if the combination passes the provided cut. \n"
3987 "For instance one can apply this function on a signal Photon and provide a list of all photons in the rest of event and a cut \n"
3988 "around the neutral Pion mass (e.g. ``0.130 < M < 0.140``). \n"
3989 "If a combination of the signal Photon with a ROE photon fits this criteria, hence looks like a neutral pion, the veto-Metavariable will return 1", Manager::VariableDataType::c_bool);
3990 REGISTER_METAVARIABLE("matchedMC(variable)", matchedMC,
3991 "Returns variable output for the matched MCParticle by constructing a temporary Particle from it.\n"
3992 "This may not work too well if your variable requires accessing daughters of the particle.\n"
3993 "E.g. ``matchedMC(p)`` returns the total momentum of the related MCParticle.\n"
3994 "Returns NaN if no matched MCParticle exists.", Manager::VariableDataType::c_double);
3995 REGISTER_METAVARIABLE("clusterBestMatchedMCParticle(variable)", clusterBestMatchedMCParticle,
3996 "Returns variable output for the MCParticle that is best-matched with the ECLCluster of the given Particle.\n"
3997 "E.g. To get the energy of the MCParticle that matches best with an ECLCluster, one could use ``clusterBestMatchedMCParticle(E)``\n"
3998 "When the variable is called for ``gamma`` and if the ``gamma`` is matched with MCParticle, it works same as `matchedMC`.\n"
3999 "If the variable is called for ``gamma`` that fails to match with an MCParticle, it provides the mdst-level MCMatching information abouth the ECLCluster.\n"
4000 "Returns NaN if the particle is not matched to an ECLCluster, or if the ECLCluster has no matching MCParticles", Manager::VariableDataType::c_double);
4001 REGISTER_METAVARIABLE("varForBestMatchedMCKlong(variable)", clusterBestMatchedMCKlong,
4002 "Returns variable output for the Klong MCParticle which has the best match with the ECLCluster of the given Particle.\n"
4003 "Returns NaN if the particle is not matched to an ECLCluster, or if the ECLCluster has no matching Klong MCParticle", Manager::VariableDataType::c_double);
4004
4005 REGISTER_METAVARIABLE("countInList(particleList, cut='')", countInList, "[Eventbased] "
4006 "Returns number of particle which pass given in cut in the specified particle list.\n"
4007 "Useful for creating statistics about the number of particles in a list.\n"
4008 "E.g. ``countInList(e+, isSignal == 1)`` returns the number of correctly reconstructed electrons in the event.\n"
4009 "The variable is event-based and does not need a valid particle pointer as input.", Manager::VariableDataType::c_int);
4010 REGISTER_METAVARIABLE("getVariableByRank(particleList, rankedVariableName, variableName, rank)", getVariableByRank, R"DOC(
4011 [Eventbased] Returns the value of ``variableName`` for the candidate in the ``particleList`` with the requested ``rank``.
4012
4013 .. note::
4014 The `BestCandidateSelection` module available via `rankByHighest` / `rankByLowest` has to be used before.
4015
4016 .. warning::
4017 The first candidate matching the given rank is used.
4018 Thus, it is not recommended to use this variable in conjunction with ``allowMultiRank`` in the `BestCandidateSelection` module.
4019
4020 The suffix ``_rank`` is automatically added to the argument ``rankedVariableName``,
4021 which either has to be the name of the variable used to order the candidates or the selected outputVariable name without the ending ``_rank``.
4022 This means that your selected name for the rank variable has to end with ``_rank``.
4023
4024 An example of this variable's usage is given in the tutorial `B2A602-BestCandidateSelection <https://gitlab.desy.de/belle2/software/basf2/-/tree/main/analysis/examples/tutorials/B2A602-BestCandidateSelection.py>`_
4025 )DOC", Manager::VariableDataType::c_double);
4026 REGISTER_VARIABLE("matchedMCHasPDG(PDGCode)", matchedMCHasPDG,
4027 "Returns if the absolute value of the PDGCode of the MCParticle related to the Particle matches a given PDGCode."
4028 "Returns 0/NAN/1 if PDGCode does not match/is not available/ matches");
4029 REGISTER_METAVARIABLE("numberOfNonOverlappingParticles(pList1, pList2, ...)", numberOfNonOverlappingParticles,
4030 "Returns the number of non-overlapping particles in the given particle lists"
4031 "Useful to check if there is additional physics going on in the detector if one reconstructed the Y4S", Manager::VariableDataType::c_int);
4032 REGISTER_METAVARIABLE("totalEnergyOfParticlesInList(particleListName)", totalEnergyOfParticlesInList,
4033 "[Eventbased] Returns the total energy of particles in the given particle List. The unit of the energy is ``GeV``", Manager::VariableDataType::c_double);
4034 REGISTER_METAVARIABLE("totalPxOfParticlesInList(particleListName)", totalPxOfParticlesInList,
4035 "[Eventbased] Returns the total momentum Px of particles in the given particle List. The unit of the momentum is ``GeV/c``", Manager::VariableDataType::c_double);
4036 REGISTER_METAVARIABLE("totalPyOfParticlesInList(particleListName)", totalPyOfParticlesInList,
4037 "[Eventbased] Returns the total momentum Py of particles in the given particle List. The unit of the momentum is ``GeV/c``", Manager::VariableDataType::c_double);
4038 REGISTER_METAVARIABLE("totalPzOfParticlesInList(particleListName)", totalPzOfParticlesInList,
4039 "[Eventbased] Returns the total momentum Pz of particles in the given particle List. The unit of the momentum is ``GeV/c``", Manager::VariableDataType::c_double);
4040 REGISTER_METAVARIABLE("invMassInLists(pList1, pList2, ...)", invMassInLists,
4041 "[Eventbased] Returns the invariant mass of the combination of particles in the given particle lists. The unit of the invariant mass is GeV/:math:`\\text{c}^2` ", Manager::VariableDataType::c_double);
4042 REGISTER_METAVARIABLE("totalECLEnergyOfParticlesInList(particleListName)", totalECLEnergyOfParticlesInList,
4043 "[Eventbased] Returns the total ECL energy of particles in the given particle List. The unit of the energy is ``GeV``", Manager::VariableDataType::c_double);
4044 REGISTER_METAVARIABLE("maxPtInList(particleListName)", maxPtInList,
4045 "[Eventbased] Returns maximum transverse momentum Pt in the given particle List. The unit of the transverse momentum is ``GeV/c``", Manager::VariableDataType::c_double);
4046 REGISTER_METAVARIABLE("eclClusterSpecialTrackMatched(cut)", eclClusterTrackMatchedWithCondition,
4047 "Returns if at least one Track that satisfies the given condition is related to the ECLCluster of the Particle.", Manager::VariableDataType::c_double);
4048 REGISTER_METAVARIABLE("averageValueInList(particleListName, variable)", averageValueInList,
4049 "[Eventbased] Returns the arithmetic mean of the given variable of the particles in the given particle list.", Manager::VariableDataType::c_double);
4050 REGISTER_METAVARIABLE("medianValueInList(particleListName, variable)", medianValueInList,
4051 "[Eventbased] Returns the median value of the given variable of the particles in the given particle list.", Manager::VariableDataType::c_double);
4052 REGISTER_METAVARIABLE("sumValueInList(particleListName, variable)", sumValueInList,
4053 "[Eventbased] Returns the sum of the given variable of the particles in the given particle list.", Manager::VariableDataType::c_double);
4054 REGISTER_METAVARIABLE("productValueInList(particleListName, variable)", productValueInList,
4055 "[Eventbased] Returns the product of the given variable of the particles in the given particle list.", Manager::VariableDataType::c_double);
4056 REGISTER_METAVARIABLE("angleToClosestInList(particleListName)", angleToClosestInList,
4057 "Returns the angle between this particle and the closest particle (smallest opening angle) in the list provided. The unit of the angle is ``rad`` ", Manager::VariableDataType::c_double);
4058 REGISTER_METAVARIABLE("closestInList(particleListName, variable)", closestInList,
4059 "Returns `variable` for the closest particle (smallest opening angle) in the list provided.", Manager::VariableDataType::c_double);
4060 REGISTER_METAVARIABLE("angleToMostB2BInList(particleListName)", angleToMostB2BInList,
4061 "Returns the angle between this particle and the most back-to-back particle (closest opening angle to 180) in the list provided. The unit of the angle is ``rad`` ", Manager::VariableDataType::c_double);
4062 REGISTER_METAVARIABLE("deltaPhiToMostB2BPhiInList(particleListName)", deltaPhiToMostB2BPhiInList,
4063 "Returns the abs(delta phi) between this particle and the most back-to-back particle in phi (closest opening angle to 180) in the list provided. The unit of the angle is ``rad`` ", Manager::VariableDataType::c_double);
4064 REGISTER_METAVARIABLE("mostB2BInList(particleListName, variable)", mostB2BInList,
4065 "Returns `variable` for the most back-to-back particle (closest opening angle to 180) in the list provided.", Manager::VariableDataType::c_double);
4066 REGISTER_METAVARIABLE("maxOpeningAngleInList(particleListName)", maxOpeningAngleInList,
4067 "[Eventbased] Returns maximum opening angle in the given particle List. The unit of the angle is ``rad`` ", Manager::VariableDataType::c_double);
4068 REGISTER_METAVARIABLE("daughterCombination(variable, daughterIndex_1, daughterIndex_2 ... daughterIndex_n)", daughterCombination,R"DOC(
4069Returns a ``variable`` function only of the 4-momentum calculated on an arbitrary set of (grand)daughters.
4070
4071.. warning::
4072 ``variable`` can only be a function of the daughters' 4-momenta.
4073
4074Daughters from different generations of the decay tree can be combined using generalized daughter indexes, which are simply colon-separated
4075the list of daughter indexes, starting from the root particle: for example, ``0:1:3`` identifies the fourth
4076daughter (3) of the second daughter (1) of the first daughter (0) of the mother particle.
4077
4078.. tip::
4079 ``daughterCombination(M, 0, 3, 4)`` will return the invariant mass of the system made of the first, fourth and fifth daughter of particle.
4080 ``daughterCombination(M, 0:0, 3:0)`` will return the invariant mass of the system made of the first daughter of the first daughter and the first daughter of the fourth daughter.
4081
4082)DOC", Manager::VariableDataType::c_double);
4083 REGISTER_METAVARIABLE("useAlternativeDaughterHypothesis(variable, daughterIndex_1:newMassHyp_1, ..., daughterIndex_n:newMassHyp_n)", useAlternativeDaughterHypothesis,R"DOC(
4084Returns a ``variable`` calculated using new mass hypotheses for (some of) the particle's daughters.
4085
4086.. warning::
4087 ``variable`` can only be a function of the particle 4-momentum, which is re-calculated as the sum of the daughters' 4-momenta, and the daughters' 4-momentum.
4088 This means that if you made a kinematic fit without updating the daughters' momenta, the result of this variable will not reflect the effect of the kinematic fit.
4089 Also, the track fit is not performed again: the variable only re-calculates the 4-vectors using different mass assumptions.
4090 In the variable, a copy of the given particle is created with daughters' alternative mass assumption (i.e. the original particle and daughters are not changed).
4091
4092.. warning::
4093 Generalized daughter indexes are not supported (yet!): this variable can be used only on first-generation daughters.
4094
4095.. tip::
4096 ``useAlternativeDaughterHypothesis(M, 0:K+, 2:pi-)`` will return the invariant mass of the particle assuming that the first daughter is a kaon and the third is a pion, instead of whatever was used in reconstructing the decay.
4097 ``useAlternativeDaughterHypothesis(mRecoil, 1:p+)`` will return the recoil mass of the particle assuming that the second daughter is a proton instead of whatever was used in reconstructing the decay.
4098
4099)DOC", Manager::VariableDataType::c_double);
4100 REGISTER_METAVARIABLE("varForFirstMCAncestorOfType(type, variable)",varForFirstMCAncestorOfType,R"DOC(Returns requested variable of the first ancestor of the given type.
4101Ancestor type can be set up by PDG code or by particle name (check evt.pdl for valid particle names))DOC", Manager::VariableDataType::c_double);
4102
4103 REGISTER_METAVARIABLE("nTrackFitResults(particleType)", nTrackFitResults,
4104 "[Eventbased] Returns the total number of TrackFitResults for a given particleType. The argument can be the name of particle (e.g. pi+) or PDG code (e.g. 211).",
4105 Manager::VariableDataType::c_int);
4106
4107 }
4109}
int getPDGCode() const
PDG code.
Definition Const.h:473
static const ChargedStable pion
charged pion particle
Definition Const.h:661
static const double doubleNaN
quiet_NaN
Definition Const.h:703
static const ChargedStable electron
electron particle
Definition Const.h:659
EHypothesisBit
The hypothesis bits for this ECLCluster (Connected region (CR) is split using this hypothesis.
Definition ECLCluster.h:31
@ c_nPhotons
CR is split into n photons (N1)
Definition ECLCluster.h:41
static std::unique_ptr< GeneralCut > compile(const std::string &cut)
Definition GeneralCut.h:84
@ c_Initial
bit 5: Particle is initial such as e+ or e- and not going to Geant4
Definition MCParticle.h:57
@ c_PrimaryParticle
bit 0: Particle is primary particle.
Definition MCParticle.h:47
@ c_IsVirtual
bit 4: Particle is virtual and not going to Geant4.
Definition MCParticle.h:55
static std::string makeROOTCompatible(std::string str)
Remove special characters that ROOT dislikes in branch names, e.g.
EParticleSourceObject
particle source enumerators
Definition Particle.h:83
@ c_Flavored
Is either particle or antiparticle.
Definition Particle.h:98
static const ReferenceFrame & GetCurrent()
Get current rest frame.
std::function< VarVariant(const Particle *)> FunctionPtr
functions stored take a const Particle* and return VarVariant.
Definition Manager.h:112
const Var * getVariable(std::string name)
Get the variable belonging to the given key.
Definition Manager.cc:58
static Manager & Instance()
get singleton instance.
Definition Manager.cc:26
#define MAKE_DEPRECATED(name, make_fatal, version, description)
Registers a variable as deprecated.
Definition Manager.h:448
T convertString(const std::string &str)
Converts a string to type T (one of float, double, long double, int, long int, unsigned long int).
double charge(int pdgCode)
Returns electric charge of a particle with given pdg code.
Definition EvtPDLUtil.cc:44
bool hasAntiParticle(int pdgCode)
Checks if the particle with given pdg code has an anti-particle or not.
Definition EvtPDLUtil.cc:12
Particle * copyParticle(const Particle *original)
Function takes argument Particle and creates a copy of it and copies of all its (grand-)^n-daughters.
Abstract base class for different kinds of events.
const std::vector< double > v2
MATLAB generated random vector.
const std::vector< double > v1
MATLAB generated random vector.