Belle II Software development
FlavorTaggingVariables.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/FlavorTaggingVariables.h>
11#include <analysis/variables/ROEVariables.h>
12
13#include <analysis/variables/MCTruthVariables.h>
14#include <analysis/variables/KLMClusterVariables.h>
15
16#include <analysis/ClusterUtility/ClusterUtils.h>
17
18#include <analysis/utility/MCMatching.h>
19#include <analysis/utility/PCmsLabTransform.h>
20
21// framework - DataStore
22#include <framework/datastore/StoreObjPtr.h>
23
24// dataobjects
25#include <analysis/dataobjects/Particle.h>
26#include <analysis/dataobjects/RestOfEvent.h>
27#include <analysis/dataobjects/ParticleList.h>
28#include <analysis/dataobjects/FlavorTaggerInfo.h>
29#include <analysis/dataobjects/FlavorTaggerInfoMap.h>
30#include <analysis/ContinuumSuppression/Thrust.h>
31
32#include <mdst/dataobjects/MCParticle.h>
33#include <mdst/dataobjects/Track.h>
34#include <mdst/dataobjects/ECLCluster.h>
35#include <mdst/dataobjects/KLMCluster.h>
36#include <mdst/dataobjects/PIDLikelihood.h>
37
38// framework aux
39#include <framework/gearbox/Const.h>
40#include <framework/logging/Logger.h>
41
42#include <Math/Vector3D.h>
43#include <Math/Vector4D.h>
44#include <Math/VectorUtil.h>
45
46#include <algorithm>
47#include <cmath>
48
49using namespace std;
50
51namespace Belle2 {
56 namespace Variable {
57
58 // ############################################## FlavorTagger Variables ###############################################
59
60 // Track Level Variables ---------------------------------------------------------------------------------------------------
61
62 double momentumMissingTagSide(const Particle*)
63 {
64 StoreObjPtr<RestOfEvent> roe("RestOfEvent");
65 if (!roe.isValid()) return 0;
66
67 ROOT::Math::PxPyPzEVector roeCMSVec;
68
69 const auto& roeChargedParticles = roe->getChargedParticles();
70 for (auto roeChargedParticle : roeChargedParticles) {
71 roeCMSVec += PCmsLabTransform::labToCms(roeChargedParticle->get4Vector());
72 }
73
74 double missMom = -roeCMSVec.P();
75 return missMom ;
76 }
77
78 Manager::FunctionPtr momentumMissingTagSideWithMask(const std::vector<std::string>& arguments)
79 {
80 std::string maskName;
81 if (arguments.size() == 0)
83 else if (arguments.size() == 1)
84 maskName = arguments[0];
85 else
86 B2FATAL("At most 1 argument (name of mask) accepted.");
87
88 auto func = [maskName](const Particle*) -> double {
89 StoreObjPtr<RestOfEvent> roe("RestOfEvent");
90 if (!roe.isValid()) return 0;
91
92 ROOT::Math::PxPyPzEVector roeCMSVec;
93
94 const auto& roeChargedParticles = roe->getChargedParticles(maskName);
95 for (auto roeChargedParticle : roeChargedParticles)
96 {
97 roeCMSVec += PCmsLabTransform::labToCms(roeChargedParticle->get4Vector());
98 }
99 double missMom = -roeCMSVec.P();
100 return missMom ;
101 };
102 return func;
103 }
104
105 double cosTPTO(const Particle* part)
106 {
107 StoreObjPtr<RestOfEvent> roeobject("RestOfEvent");
108 const RestOfEvent* roe;
109 if (roeobject.isValid()) { // if in for_each loop over ROE
110 roe = &(*roeobject);
111 } else {
112 roe = getRelatedROEObject(part);
113 if (!roe)
114 return 0;
115 }
116
117 std::vector<ROOT::Math::PxPyPzEVector> p_cms_roe;
118 static const double P_MAX(3.2);
119
120 // Charged tracks
121 const auto& roeTracks = roe->getChargedParticles();
122 for (auto& roeChargedParticle : roeTracks) {
123 // TODO: Add helix and KVF with IpProfile once available. Port from L163-199 of:
124 // /belle/b20090127_0910/src/anal/ekpcontsuppress/src/ksfwmoments.cc
125 ROOT::Math::PxPyPzEVector p_cms = PCmsLabTransform::labToCms(roeChargedParticle->get4Vector());
126 if (p_cms != p_cms) continue;
127 if (p_cms.P() > P_MAX) continue;
128 p_cms_roe.push_back(p_cms);
129 }
130
131 // ECLCluster->Gamma
132 const auto& roePhotons = roe->getPhotons();
133 for (auto& roePhoton : roePhotons) {
134 if (roePhoton->getECLClusterEHypothesisBit() == ECLCluster::EHypothesisBit::c_nPhotons) {
135 ROOT::Math::PxPyPzEVector p_lab = roePhoton->get4Vector();
136 if (p_lab != p_lab) continue;
137 if (p_lab.P() < 0.05) continue;
138 ROOT::Math::PxPyPzEVector p_cms = PCmsLabTransform::labToCms(p_lab);
139 if (p_cms != p_cms) continue;
140 if (p_cms.P() > P_MAX) continue;
141 p_cms_roe.push_back(p_cms);
142 }
143 }
144
145 const auto& roeKlongs = roe->getHadrons();
146 for (auto& roeKlong : roeKlongs) {
147 if (nKLMClusterTrackMatches(roeKlong) == 0 && !(roeKlong->getKLMCluster()->getAssociatedEclClusterFlag())) {
148 ROOT::Math::PxPyPzEVector p_lab = roeKlong->get4Vector();
149 if (p_lab != p_lab) continue;
150 if (p_lab.P() < 0.05) continue;
151 ROOT::Math::PxPyPzEVector p_cms = PCmsLabTransform::labToCms(p_lab);
152 if (p_cms != p_cms) continue;
153 if (p_cms.P() > P_MAX) continue;
154 p_cms_roe.push_back(p_cms);
155 }
156 }
157
158 const ROOT::Math::XYZVector thrustO = Thrust::calculateThrust(p_cms_roe);
159 const ROOT::Math::PxPyPzEVector pAxis = PCmsLabTransform::labToCms(part->get4Vector());
160
161 double result = 0 ;
162 if (pAxis == pAxis) result = abs(ROOT::Math::VectorUtil::CosTheta(pAxis, thrustO));
163
164 return result;
165 }
166
167 Manager::FunctionPtr cosTPTOWithMask(const std::vector<std::string>& arguments)
168 {
169 std::string maskName;
170 if (arguments.size() == 0)
172 else if (arguments.size() == 1)
173 maskName = arguments[0];
174 else
175 B2FATAL("At most 1 argument (name of mask) accepted.");
176
177 auto func = [maskName](const Particle * particle) -> double {
178 StoreObjPtr<RestOfEvent> roeobject("RestOfEvent");
179 const RestOfEvent* roe;
180 if (roeobject.isValid()) // if in for_each loop over ROE
181 {
182 roe = &(*roeobject);
183 } else
184 {
185 roe = getRelatedROEObject(particle);
186 if (!roe)
187 return 0;
188 }
189
190 std::vector<ROOT::Math::PxPyPzEVector> p_cms_roe;
191 // static const double P_MAX(3.2);
192
193 // Charged tracks
194 const auto& roeTracks = roe->getChargedParticles(maskName);
195 for (auto& roeChargedParticle : roeTracks)
196 {
197 // TODO: Add helix and KVF with IpProfile once available. Port from L163-199 of:
198 // /belle/b20090127_0910/src/anal/ekpcontsuppress/src/ksfwmoments.cc
199 ROOT::Math::PxPyPzEVector p_cms = PCmsLabTransform::labToCms(roeChargedParticle->get4Vector());
200 if (p_cms != p_cms) continue;
201 // if (p_cms.P() > P_MAX) continue; // Should not be added without any description.
202 p_cms_roe.push_back(p_cms);
203 }
204
205 // ECLCluster->Gamma
206 const auto& roePhotons = roe->getPhotons(maskName);
207 for (auto& roePhoton : roePhotons)
208 {
209 if (roePhoton->getECLClusterEHypothesisBit() == ECLCluster::EHypothesisBit::c_nPhotons) {
210 ROOT::Math::PxPyPzEVector p_lab = roePhoton->get4Vector();
211 if (p_lab != p_lab) continue;
212 // if (p_lab.P() < 0.05) continue; // Should not be added without any description.
213 ROOT::Math::PxPyPzEVector p_cms = PCmsLabTransform::labToCms(p_lab);
214 if (p_cms != p_cms) continue;
215 // if (p_cms.P() > P_MAX) continue; // Should not be added without any description.
216 p_cms_roe.push_back(p_cms);
217 }
218 }
219
220 // KLMCluster
221 const auto& roeKlongs = roe->getHadrons(maskName);
222 for (auto& roeKlong : roeKlongs)
223 {
224 if (nKLMClusterTrackMatches(roeKlong) == 0 && !(roeKlong->getKLMCluster()->getAssociatedEclClusterFlag())) {
225 ROOT::Math::PxPyPzEVector p_lab = roeKlong->get4Vector();
226 if (p_lab != p_lab) continue;
227 // if (p_lab.P() < 0.05) continue; // Should not be added without any description.
228 ROOT::Math::PxPyPzEVector p_cms = PCmsLabTransform::labToCms(p_lab);
229 if (p_cms != p_cms) continue;
230 // if (p_cms.P() > P_MAX) continue; // Should not be added without any description.
231 p_cms_roe.push_back(p_cms);
232 }
233 }
234
235 const ROOT::Math::XYZVector thrustO = Thrust::calculateThrust(p_cms_roe);
236 const ROOT::Math::PxPyPzEVector pAxis = PCmsLabTransform::labToCms(particle->get4Vector());
237
238 double result = 0 ;
239 if (pAxis == pAxis)
240 result = abs(ROOT::Math::VectorUtil::CosTheta(pAxis, thrustO)); // abs??
241
242 return result;
243
244 };
245 return func;
246 }
247
248 int lambdaFlavor(const Particle* particle)
249 {
250 if (particle->getPDGCode() == Const::Lambda.getPDGCode()) return 1; //Lambda0
251 else if (particle->getPDGCode() == Const::antiLambda.getPDGCode()) return -1; //Anti-Lambda0
252 else return 0;
253 }
254
255 bool isLambda(const Particle* particle)
256 {
257 const MCParticle* mcparticle = particle->getMCParticle();
258 if (!mcparticle) return false;
259 return (abs(mcparticle->getPDG()) == Const::Lambda.getPDGCode());
260 }
261
262 double lambdaZError(const Particle* particle)
263 {
264 //This is a simplistic hack. But I see no other way to get that information.
265 //Should be removed if worthless
266 TMatrixFSym ErrorPositionMatrix = particle->getVertexErrorMatrix();
267 return ErrorPositionMatrix[2][2];
268 }
269
270 double momentumOfSecondDaughter(const Particle* part)
271 {
272 if (!part->getDaughter(1)) return 0.0;
273 return part->getDaughter(1)->getP();
274 }
275
276 double momentumOfSecondDaughterCMS(const Particle* part)
277 {
278 if (!part->getDaughter(1)) return 0.0;
279 ROOT::Math::PxPyPzEVector vec = PCmsLabTransform::labToCms(part->getDaughter(1)->get4Vector());
280 return vec.P();
281 }
282
283 double chargeTimesKaonLiklihood(const Particle*)
284 {
285 StoreObjPtr<ParticleList> KaonList("K+:inRoe");
286 if (!KaonList.isValid()) return 0;
287
288 double maximumKaonid = 0;
289 double maximum_charge = 0;
290 for (unsigned int i = 0; i < KaonList->getListSize(); ++i) {
291 const Particle* p = KaonList->getParticle(i);
292 double Kid = p->getRelatedTo<PIDLikelihood>()->getProbability(Const::kaon, Const::pion);
293 if (Kid > maximumKaonid) {
294 maximumKaonid = Kid;
295 maximum_charge = p->getCharge();
296 }
297 }
298 return maximumKaonid * maximum_charge;
299 }
300
301 double transverseMomentumOfChargeTracksInRoe(const Particle* part)
302 {
303 StoreObjPtr<RestOfEvent> roe("RestOfEvent");
304 if (!roe.isValid()) return 0;
305
306 double sum = 0.0;
307
308 for (const auto& track : roe->getChargedParticles()) {
309 if (part->isCopyOf(track, true)) continue;
310 sum += track->getMomentum().Perp2();
311 }
312
313 return sum;
314
315 }
316
317 Manager::FunctionPtr transverseMomentumOfChargeTracksInRoeWithMask(const std::vector<std::string>& arguments)
318 {
319 std::string maskName;
320 if (arguments.size() == 0)
322 else if (arguments.size() == 1)
323 maskName = arguments[0];
324 else
325 B2FATAL("At most 1 argument (name of mask) accepted.");
326
327 auto func = [maskName](const Particle * particle) -> double {
328 StoreObjPtr<RestOfEvent> roe("RestOfEvent");
329 if (!roe.isValid()) return 0;
330
331 double sum = 0.0;
332
333 for (const auto& track : roe->getChargedParticles(maskName))
334 {
335 if (particle->isCopyOf(track, true)) continue;
336 sum += track->getMomentum().Rho();
337 }
338
339 return sum;
340 };
341 return func;
342 }
343
344 Manager::FunctionPtr transverseMomentumSquaredOfChargeTracksInRoeWithMask(const std::vector<std::string>& arguments)
345 {
346 std::string maskName;
347 if (arguments.size() == 0)
349 else if (arguments.size() == 1)
350 maskName = arguments[0];
351 else
352 B2FATAL("At most 1 argument (name of mask) accepted.");
353
354 auto func = [maskName](const Particle * particle) -> double {
355 StoreObjPtr<RestOfEvent> roe("RestOfEvent");
356 if (!roe.isValid()) return 0;
357
358 double sum = 0.0;
359
360 for (const auto& track : roe->getChargedParticles(maskName))
361 {
362 if (particle->isCopyOf(track, true)) continue;
363 sum += track->getMomentum().Perp2();
364 }
365
366 return sum;
367 };
368 return func;
369 }
370
371 int NumberOfKShortsInRoe(const Particle* particle)
372 {
373 StoreObjPtr<ParticleList> KShortList("K_S0:inRoe");
374 if (!KShortList.isValid())
375 B2FATAL("NumberOfKShortsInRoe cannot be calculated because the required particleList K_S0:inRoe could not be found or is not valid");
376
377 int flag = 0;
378 for (unsigned int i = 0; i < KShortList->getListSize(); ++i) {
379 if (!particle->overlapsWith(KShortList->getParticle(i)))
380 ++flag;
381 }
382 return flag;
383 }
384
385// Event Level Variables --------------------------------------------------------------------------------------------
386
387 bool isInElectronOrMuonCat(const Particle* particle)
388 {
389 // check muons
390 StoreObjPtr<ParticleList> MuonList("mu+:inRoe");
391 const Track* trackTargetMuon = nullptr;
392 if (MuonList.isValid()) {
393 double maximumProbMuon = 0;
394 for (unsigned int i = 0; i < MuonList->getListSize(); ++i) {
395 Particle* pMuon = MuonList->getParticle(i);
396 double probMuon = pMuon->getExtraInfo("isRightTrack(Muon)");
397 if (probMuon > maximumProbMuon) {
398 maximumProbMuon = probMuon;
399 trackTargetMuon = pMuon->getTrack();
400 }
401 }
402 }
403 if (particle->getTrack() == trackTargetMuon)
404 return true;
405
406
407 // check electrons
408 StoreObjPtr<ParticleList> ElectronList("e+:inRoe");
409 const Track* trackTargetElectron = nullptr;
410 if (ElectronList.isValid()) {
411 double maximumProbElectron = 0;
412 for (unsigned int i = 0; i < ElectronList->getListSize(); ++i) {
413 Particle* pElectron = ElectronList->getParticle(i);
414 double probElectron = pElectron->getExtraInfo("isRightTrack(Electron)");
415 if (probElectron > maximumProbElectron) {
416 maximumProbElectron = probElectron;
417 trackTargetElectron = pElectron->getTrack();
418 }
419 }
420 }
421 if (particle->getTrack() == trackTargetElectron)
422 return true;
423
424 return false;
425 }
426
427 // helper function to get flavour of MC B0
428 static int getB0flavourMC(const MCParticle* mcParticle)
429 {
430 while (mcParticle) {
431 if (mcParticle->getPDG() == 511) {
432 return 1;
433 } else if (mcParticle->getPDG() == -511) {
434 return -1;
435 }
436 mcParticle = mcParticle->getMother();
437 }
438 return 0; //no B found
439 }
440
441// Target Variables --------------------------------------------------------------------------------------------------
442
443 bool isMajorityInRestOfEventFromB0(const Particle*)
444 {
445 StoreObjPtr<RestOfEvent> roe("RestOfEvent");
446 if (!roe.isValid()) return 0;
447
448 int vote = 0;
449 for (auto& track : roe->getChargedParticles()) {
450 const MCParticle* mcParticle = track->getMCParticle();
451 vote += getB0flavourMC(mcParticle);
452 }
453
454 return vote > 0;
455 }
456
457 Manager::FunctionPtr isMajorityInRestOfEventFromB0WithMask(const std::vector<std::string>& arguments)
458 {
459 std::string maskName;
460 if (arguments.size() == 0)
462 else if (arguments.size() == 1)
463 maskName = arguments[0];
464 else
465 B2FATAL("At most 1 argument (name of mask) accepted.");
466
467 auto func = [maskName](const Particle*) -> bool {
468 StoreObjPtr<RestOfEvent> roe("RestOfEvent");
469 if (!roe.isValid()) return 0;
470
471 int vote = 0;
472 for (auto& track : roe->getChargedParticles(maskName))
473 {
474 const MCParticle* mcParticle = track->getMCParticle();
475 vote += getB0flavourMC(mcParticle);
476 }
477
478 return vote > 0;
479
480 };
481 return func;
482 }
483
484 bool isMajorityInRestOfEventFromB0bar(const Particle*)
485 {
486 StoreObjPtr<RestOfEvent> roe("RestOfEvent");
487 if (!roe.isValid()) return 0;
488
489 int vote = 0;
490 for (auto& track : roe->getChargedParticles()) {
491 const MCParticle* mcParticle = track->getMCParticle();
492 vote += getB0flavourMC(mcParticle);
493 }
494
495 return vote < 0;
496 }
497
498 Manager::FunctionPtr isMajorityInRestOfEventFromB0barWithMask(const std::vector<std::string>& arguments)
499 {
500 std::string maskName;
501 if (arguments.size() == 0)
503 else if (arguments.size() == 1)
504 maskName = arguments[0];
505 else
506 B2FATAL("At most 1 argument (name of mask) accepted.");
507
508 auto func = [maskName](const Particle*) -> bool {
509 StoreObjPtr<RestOfEvent> roe("RestOfEvent");
510 if (!roe.isValid()) return 0;
511
512 int vote = 0;
513 for (auto& track : roe->getChargedParticles(maskName))
514 {
515 const MCParticle* mcParticle = track->getMCParticle();
516 vote += getB0flavourMC(mcParticle);
517 }
518
519 return vote < 0;
520
521 };
522 return func;
523 }
524
525 bool hasRestOfEventTracks(const Particle* part)
526 {
527 const RestOfEvent* roe = part->getRelatedTo<RestOfEvent>();
528 return (roe && roe-> getNTracks() > 0);
529 }
530
531 Manager::FunctionPtr hasRestOfEventTracksWithMask(const std::vector<std::string>& arguments)
532 {
533 std::string maskName;
534 if (arguments.size() == 0)
536 else if (arguments.size() == 1)
537 maskName = arguments[0];
538 else
539 B2FATAL("At most 1 argument (name of mask) accepted.");
540
541 auto func = [maskName](const Particle*) -> bool {
542 StoreObjPtr<RestOfEvent> roe("RestOfEvent");
543 if (!roe.isValid()) return 0;
544
545 return roe->getNTracks(maskName) > 0;
546 };
547 return func;
548 }
549
550 int isRelatedRestOfEventB0Flavor(const Particle* particle)
551 {
552 const RestOfEvent* roe = particle->getRelatedTo<RestOfEvent>();
553 if (!roe) return 0;
554
555 const MCParticle* BcpMC = particle->getMCParticle();
556 if (!BcpMC) return 0;
557 if (Variable::isSignal(particle) <= 0) return 0;
558
559 const MCParticle* Y4S = BcpMC->getMother();
560 if (!Y4S) return 0;
561
562 int BtagFlavor = 0;
563 int BcpFlavor = 0;
564
565 for (auto& roeChargedParticle : roe->getChargedParticles()) {
566 const MCParticle* mcParticle = roeChargedParticle->getMCParticle();
567 while (mcParticle) {
568 if (mcParticle->getMother() == Y4S) {
569 if (mcParticle == BcpMC) {
570 if (mcParticle->getPDG() > 0) BcpFlavor = 2;
571 else BcpFlavor = -2;
572 } else if (BtagFlavor == 0) {
573 if (abs(mcParticle->getPDG()) == 511 || abs(mcParticle->getPDG()) == 521) {
574 if (mcParticle->getPDG() > 0) BtagFlavor = 1;
575 else BtagFlavor = -1;
576 } else BtagFlavor = 5;
577 }
578 break;
579 }
580 mcParticle = mcParticle->getMother();
581 }
582 if (BcpFlavor != 0 || BtagFlavor == 5) break;
583 }
584
585
586 return (BcpFlavor != 0) ? BcpFlavor : BtagFlavor;
587 }
588
589 Manager::FunctionPtr isRelatedRestOfEventB0FlavorWithMask(const std::vector<std::string>& arguments)
590 {
591 std::string maskName;
592 if (arguments.size() == 0)
594 else if (arguments.size() == 1)
595 maskName = arguments[0];
596 else
597 B2FATAL("At most 1 argument (name of mask) accepted.");
598
599 auto func = [maskName](const Particle * particle) -> int {
600 StoreObjPtr<RestOfEvent> roe("RestOfEvent");
601 if (!roe.isValid()) return 0;
602
603 const MCParticle* BcpMC = particle->getMCParticle();
604 if (!BcpMC) return 0;
605 if (Variable::isSignal(particle) <= 0) return 0;
606
607 const MCParticle* Y4S = BcpMC->getMother();
608 if (!Y4S) return 0;
609
610 int BtagFlavor = 0;
611 int BcpFlavor = 0;
612
613 for (auto& roeChargedParticle : roe->getChargedParticles(maskName))
614 {
615 const MCParticle* mcParticle = roeChargedParticle->getMCParticle();
616 while (mcParticle) {
617 if (mcParticle->getMother() != Y4S) {
618 mcParticle = mcParticle->getMother();
619 continue;
620 }
621
622 if (mcParticle == BcpMC) { // if mcParticle is associated with CP-side unfortunately
623 if (mcParticle->getPDG() > 0)
624 BcpFlavor = 2;
625 else
626 BcpFlavor = -2;
627 } else if (BtagFlavor == 0) { // only first mcParticle is checked.
628 if (abs(mcParticle->getPDG()) == 511 || abs(mcParticle->getPDG()) == 521) {
629 if (mcParticle->getPDG() > 0)
630 BtagFlavor = 1;
631 else
632 BtagFlavor = -1;
633 } else {
634 BtagFlavor = 5;
635 }
636 }
637 break;
638 }
639 if (BcpFlavor != 0 || BtagFlavor == 5) break;
640 }
641
642 return (BcpFlavor != 0) ? BcpFlavor : BtagFlavor;
643 };
644 return func;
645
646 }
647
648 int isRestOfEventB0Flavor(const Particle*)
649 {
650 StoreObjPtr<RestOfEvent> roe("RestOfEvent");
651 if (!roe.isValid()) return 0;
652
653 const Particle* Bcp = roe->getRelatedFrom<Particle>();
654 return Variable::isRelatedRestOfEventB0Flavor(Bcp);
655 }
656
657 int ancestorHasWhichFlavor(const Particle* particle)
658 {
659 StoreObjPtr<RestOfEvent> roe("RestOfEvent");
660 if (!roe.isValid()) return 0;
661
662 const MCParticle* BcpMC = roe->getRelatedFrom<Particle>()->getMCParticle();
663 const MCParticle* Y4S = BcpMC->getMother();
664 const MCParticle* mcParticle = particle->getMCParticle();
665
666 int outputB0tagQ = 0;
667 while (mcParticle) {
668 if (mcParticle->getMother() == Y4S) {
669 if (mcParticle != BcpMC && abs(mcParticle->getPDG()) == 511) {
670 if (mcParticle->getPDG() == 511) outputB0tagQ = 1;
671 else outputB0tagQ = -1;
672 } else if (mcParticle == BcpMC) {
673 if (mcParticle->getPDG() == 511) outputB0tagQ = 2;
674 else outputB0tagQ = -2;
675 } else outputB0tagQ = 5;
676 break;
677 }
678 mcParticle = mcParticle->getMother();
679 }
680
681 return outputB0tagQ;
682 }
683
684 int B0mcErrors(const Particle*)
685 {
686 StoreObjPtr<RestOfEvent> roe("RestOfEvent");
687 if (!roe.isValid()) return -1;
688
689 const Particle* Bcp = roe->getRelatedFrom<Particle>();
690 const MCParticle* BcpMC = roe->getRelatedFrom<Particle>()->getMCParticle();
691 return MCMatching::getMCErrors(Bcp, BcpMC);
692 }
693
694 int isRelatedRestOfEventMajorityB0Flavor(const Particle* part)
695 {
696 const RestOfEvent* roe = part->getRelatedTo<RestOfEvent>();
697 if (!roe) return -2;
698
699 int q_MC = 0; //Flavor of B
700
701 if (roe->getNTracks() > 0) {
702 for (auto& track : roe->getChargedParticles()) {
703 const MCParticle* mcParticle = track->getMCParticle();
704 q_MC += getB0flavourMC(mcParticle);
705 }
706 } else if (roe->getNECLClusters() > 0) {
707 for (auto& cluster : roe->getPhotons()) {
708 if (cluster->getECLClusterEHypothesisBit() != ECLCluster::EHypothesisBit::c_nPhotons) continue;
709 const MCParticle* mcParticle = cluster->getMCParticle();
710 q_MC += getB0flavourMC(mcParticle);
711 }
712 } else if (roe->getNKLMClusters() > 0) {
713 for (auto& klmcluster : roe->getHadrons()) {
714 const MCParticle* mcParticle = klmcluster->getMCParticle();
715 q_MC += getB0flavourMC(mcParticle);
716 }
717 }
718
719 if (q_MC == 0)
720 return -2;
721 else
722 return (q_MC > 0);
723 }
724
725 Manager::FunctionPtr isRelatedRestOfEventMajorityB0FlavorWithMask(const std::vector<std::string>& arguments)
726 {
727 std::string maskName;
728 if (arguments.size() == 0)
730 else if (arguments.size() == 1)
731 maskName = arguments[0];
732 else
733 B2FATAL("At most 1 argument (name of mask) accepted.");
734
735 auto func = [maskName](const Particle*) -> int {
736 StoreObjPtr<RestOfEvent> roe("RestOfEvent");
737 if (!roe.isValid()) return -2;
738
739 int q_MC = 0; //Flavor of B
740
741 if (roe->getNTracks(maskName) > 0)
742 {
743 for (auto& track : roe->getChargedParticles(maskName)) {
744 const MCParticle* mcParticle = track->getMCParticle();
745 q_MC += getB0flavourMC(mcParticle);
746 }
747 } else if (roe->getNECLClusters(maskName) > 0) // only if there are no tracks
748 {
749 for (auto& cluster : roe->getPhotons(maskName)) {
750 if (cluster->getECLClusterEHypothesisBit() != ECLCluster::EHypothesisBit::c_nPhotons) continue;
751 const MCParticle* mcParticle = cluster->getMCParticle();
752 q_MC += getB0flavourMC(mcParticle);
753 }
754 } else if (roe->getNKLMClusters(maskName) > 0) // only if there are no tracks nor ecl-clusters
755 {
756 for (auto& klmcluster : roe->getHadrons(maskName)) {
757 const MCParticle* mcParticle = klmcluster->getMCParticle();
758 q_MC += getB0flavourMC(mcParticle);
759 }
760 }
761
762 if (q_MC == 0)
763 return -2;
764 else
765 return int(q_MC > 0);
766
767 };
768 return func;
769 }
770
771 int isRestOfEventMajorityB0Flavor(const Particle*)
772 {
773 StoreObjPtr<RestOfEvent> roe("RestOfEvent");
774 if (!roe.isValid()) return -2;//gRandom->Uniform(0, 1);
775
776 const Particle* Bcp = roe->getRelatedFrom<Particle>();
777 return Variable::isRelatedRestOfEventMajorityB0Flavor(Bcp);
778 }
779
780 double mcFlavorOfOtherB(const Particle* particle)
781 {
782
783 if (std::abs(particle->getPDGCode()) != 511 && std::abs(particle->getPDGCode()) != 521) {
784 B2ERROR("MCFlavorOfOtherB: this variable works only for B mesons.\n"
785 "The given particle with PDG code " << particle->getPDGCode() <<
786 " is not a B-meson candidate (PDG code 511 or 521). ");
787 return Const::doubleNaN;
788 }
789
790 const MCParticle* mcParticle = particle->getMCParticle();
791 if (!mcParticle) return Const::doubleNaN;
792
793 const MCParticle* mcMother = mcParticle->getMother();
794 if (!mcMother) return Const::doubleNaN;
795
796 if (Variable::isSignal(particle) < 1.0) return 0;
797
798 for (auto& upsilon4SDaughter : mcMother->getDaughters()) {
799 if (upsilon4SDaughter == mcParticle) continue;
800 return (upsilon4SDaughter->getPDG() > 0) ? 1 : -1;
801 }
802
803 return 0;
804
805 };
806
807// ######################################### Meta Variables ##############################################
808
809// Track and Event Level variables ------------------------------------------------------------------------
810
811 Manager::FunctionPtr BtagToWBosonVariables(const std::vector<std::string>& arguments)
812 {
813
814 std::string requestedVariable;
815 std::string maskName;
816 if (arguments.size() == 1) {
817 requestedVariable = arguments[0];
819 } else if (arguments.size() == 2) {
820 requestedVariable = arguments[0];
821 maskName = arguments[1];
822 } else {
823 B2FATAL("Number of arguments must be 1 (requestedVariable) or 2 (requestedVariable, maskName).");
824 }
825
826 const std::vector<string> availableVariables = {"recoilMass",
827 "recoilMassSqrd",
828 "pMissCMS",
829 "cosThetaMissCMS",
830 "EW90"
831 };
832
833 if (std::find(availableVariables.begin(), availableVariables.end(), requestedVariable) == availableVariables.end()) {
834 B2FATAL("Wrong variable " << requestedVariable <<
835 " requested. The possibilities are recoilMass, recoilMassSqrd, pMissCMS, cosThetaMissCMS or EW90");
836 }
837
838 auto func = [requestedVariable, maskName](const Particle * particle) -> double {
839 StoreObjPtr<RestOfEvent> roe("RestOfEvent");
840 if (!roe.isValid())
841 return 0;
842
843 ROOT::Math::PxPyPzEVector momXChargedTracks; //Momentum of charged X tracks in lab-System
844 const auto& roeChargedParticles = roe->getChargedParticles(maskName);
845 for (auto& roeChargedParticle : roeChargedParticles)
846 {
847 if (roeChargedParticle->isCopyOf(particle, true)) continue;
848 momXChargedTracks += roeChargedParticle->get4Vector();
849 }
850
851 ROOT::Math::PxPyPzEVector momXNeutralClusters = roe->get4VectorNeutralECLClusters(maskName); //Momentum of neutral X clusters in lab-System
852 const auto& klongs = roe->getHadrons(maskName);
853 for (auto& klong : klongs)
854 {
855 if (nKLMClusterTrackMatches(klong) == 0 && !(klong->getKLMCluster()->getAssociatedEclClusterFlag())) {
856 momXNeutralClusters += klong->get4Vector();
857 }
858 }
859
860 ROOT::Math::PxPyPzEVector momX = PCmsLabTransform::labToCms(momXChargedTracks + momXNeutralClusters); //Total Momentum of the recoiling X in CMS-System
861 ROOT::Math::PxPyPzEVector momTarget = PCmsLabTransform::labToCms(particle->get4Vector()); //Momentum of Mu in CMS-System
862 ROOT::Math::PxPyPzEVector momMiss = -(momX + momTarget); //Momentum of Anti-v in CMS-System
863
864 double output = 0.0;
865 if (requestedVariable == "recoilMass") output = momX.M();
866 if (requestedVariable == "recoilMassSqrd") output = momX.M2();
867 if (requestedVariable == "pMissCMS") output = momMiss.P();
868 if (requestedVariable == "cosThetaMissCMS") output = momTarget.Vect().Unit().Dot(momMiss.Vect().Unit());
869 if (requestedVariable == "EW90")
870 {
871
872 ROOT::Math::PxPyPzEVector momW = momTarget + momMiss; //Momentum of the W-Boson in CMS
873 float E_W_90 = 0 ; // Energy of all charged and neutral clusters in the hemisphere of the W-Boson
874
875 const auto& photons = roe->getPhotons(maskName);
876 for (auto& photon : photons) {
877 if (PCmsLabTransform::labToCms(photon->get4Vector()).Vect().Dot(momW.Vect()) > 0) {
878 E_W_90 += photon->getECLClusterEnergy();
879 }
880 }
881 for (auto& roeChargedParticle : roeChargedParticles) {
882 if (roeChargedParticle->isCopyOf(particle, true))
883 continue;
884
885 for (const ECLCluster& chargedCluster : roeChargedParticle->getTrack()->getRelationsWith<ECLCluster>()) {
886 // ignore everything except the nPhotons hypothesis
887 if (!chargedCluster.hasHypothesis(ECLCluster::EHypothesisBit::c_nPhotons))
888 continue;
889 float iEnergy = chargedCluster.getEnergy(ECLCluster::EHypothesisBit::c_nPhotons);
890 if (iEnergy == iEnergy) {
891 if (PCmsLabTransform::labToCms(ClusterUtils().Get4MomentumFromCluster(&chargedCluster,
892 ECLCluster::EHypothesisBit::c_nPhotons)).Vect().Dot(momW.Vect()) > 0)
893 E_W_90 += iEnergy;
894 }
895 }
896 }
897
898 output = E_W_90;
899 }
900
901 return output;
902 };
903 return func;
904 }
905
906 Manager::FunctionPtr KaonPionVariables(const std::vector<std::string>& arguments)
907 {
908 if (arguments.size() != 1)
909 B2FATAL("Wrong number of arguments (1 required) for meta function KaonPionVariables");
910
911
912 auto requestedVariable = arguments[0];
913 auto func = [requestedVariable](const Particle * particle) -> double {
914 // StoreObjPtr<ParticleList> KaonList("K+:ROE");
915 StoreObjPtr<ParticleList> SlowPionList("pi+:inRoe");
916
917
918 if ((requestedVariable != "HaveOpositeCharges") && (requestedVariable != "cosKaonPion"))
919 B2FATAL("Wrong variable " << requestedVariable << " requested. The possibilities are cosKaonPion or HaveOpositeCharges");
920
921
922 ROOT::Math::PxPyPzEVector momTargetSlowPion;
923 double chargeTargetSlowPion = 0;
924 if (SlowPionList.isValid())
925 {
926 double maximumProbSlowPion = 0;
927 for (unsigned int i = 0; i < SlowPionList->getListSize(); ++i) {
928 Particle* pSlowPion = SlowPionList->getParticle(i);
929 if (!pSlowPion) continue;
930 if (!pSlowPion->hasExtraInfo("isRightCategory(SlowPion)")) continue;
931
932 double probSlowPion = pSlowPion->getExtraInfo("isRightCategory(SlowPion)");
933 if (probSlowPion > maximumProbSlowPion) {
934 maximumProbSlowPion = probSlowPion;
935 chargeTargetSlowPion = pSlowPion->getCharge();
936 momTargetSlowPion = PCmsLabTransform::labToCms(pSlowPion->get4Vector());
937 }
938 }
939 }
940
941 double output = 0.0;
942
943 double chargeTargetKaon = particle->getCharge();
944 if (requestedVariable == "HaveOpositeCharges")
945 {
946 if (chargeTargetKaon * chargeTargetSlowPion == -1)
947 output = 1;
948 }
949 //TODO: when momTargetSlowPion == momTargetSlowPion fail?
950 else if (requestedVariable == "cosKaonPion")
951 {
952 ROOT::Math::PxPyPzEVector momTargetKaon = PCmsLabTransform::labToCms(particle->get4Vector());
953 if (momTargetKaon == momTargetKaon && momTargetSlowPion == momTargetSlowPion)
954 output = momTargetKaon.Vect().Unit().Dot(momTargetSlowPion.Vect().Unit());
955 }
956
957 return output;
958 };
959 return func;
960 }
961
962 Manager::FunctionPtr FSCVariables(const std::vector<std::string>& arguments)
963 {
964 if (arguments.size() != 1)
965 B2FATAL("Wrong number of arguments (1 required) for meta function FSCVariables");
966
967
968 auto requestedVariable = arguments[0];
969 auto func = [requestedVariable](const Particle * particle) -> double {
970 StoreObjPtr<ParticleList> FastParticleList("pi+:inRoe");
971 if (!FastParticleList.isValid()) return 0;
972
973
974 if ((requestedVariable != "pFastCMS") && (requestedVariable != "cosSlowFast") && (requestedVariable != "cosTPTOFast") && (requestedVariable != "SlowFastHaveOpositeCharges"))
975 B2FATAL("Wrong variable " << requestedVariable << " requested. The possibilities are pFastCMS, cosSlowFast, cosTPTOFast or SlowFastHaveOpositeCharges");
976
977
978 double maximumProbFastest = 0;
979 ROOT::Math::PxPyPzEVector momFastParticle; //Momentum of Fast Pion in CMS-System
980 Particle* TargetFastParticle = nullptr;
981 for (unsigned int i = 0; i < FastParticleList->getListSize(); ++i)
982 {
983 Particle* particlei = FastParticleList->getParticle(i);
984 if (!particlei) continue;
985
986 ROOT::Math::PxPyPzEVector momParticlei = PCmsLabTransform::labToCms(particlei->get4Vector());
987 if (momParticlei != momParticlei) continue;
988
989 double probFastest = momParticlei.P();
990 if (probFastest > maximumProbFastest) {
991 maximumProbFastest = probFastest;
992 TargetFastParticle = particlei;
993 momFastParticle = momParticlei;
994 }
995 }
996
997 // if nothing found
998 if (!TargetFastParticle) return 0;
999
1000
1001 double output = 0.0;
1002
1003 if (requestedVariable == "cosTPTOFast")
1004 output = cosTPTO(TargetFastParticle);
1005
1006 ROOT::Math::PxPyPzEVector momSlowPion = PCmsLabTransform::labToCms(particle->get4Vector()); //Momentum of Slow Pion in CMS-System
1007 if (momSlowPion == momSlowPion) // FIXME
1008 {
1009 if (requestedVariable == "cosSlowFast") {
1010 output = momSlowPion.Vect().Unit().Dot(momFastParticle.Vect().Unit());
1011 } else if (requestedVariable == "SlowFastHaveOpositeCharges") {
1012 if (particle->getCharge()*TargetFastParticle->getCharge() == -1) {
1013 output = 1;
1014 }
1015 } else {
1016 output = momFastParticle.P();
1017 }
1018 }
1019
1020 return output;
1021 };
1022 return func;
1023 }
1024
1025
1026// Target Variables ----------------------------------------------------------------------------------------------
1027
1028 // Lists used in target variables
1029 constexpr std::array<int, 18> charmMesons = {
1030 411 /*D+*/, 413 /*D*+*/, 415/*D_2*+*/, 421 /*D0*/, 423 /*D*0*/, 425/*D_2*0*/, 431/*D_s+*/, 433/*D_s*+*/, 435/*D_s2*+*/,
1031 10411 /*D_0*+*/, 10413 /*D_1+*/, 10421 /*D_0*0*/, 10423 /*D_10*/, 10431 /*D_s0*+*/, 10433/*D'_s1+*/,
1032 20413 /*D'_1+*/, 20423 /*D'_10*/, 20433/*D_s1+*/,
1033 };
1034
1035 constexpr std::array<int, 22> charmBaryons = {
1036 4112 /*Sigma_c0*/, 4114 /*Sigma_c*0*/, 4122 /*Lambda_c+*/, 4132 /*Xi_c0*/,
1037 4212 /*Sigma_c+*/, 4214 /*Sigma_c*+*/, 4222 /*Sigma_c++*/, 4224 /*sigma_c*++*/, 4232 /*Xi_c+*/,
1038 4312 /*Xi'_c0*/, 4314 /*Xi_c*0*/, 4322 /*Xi'_c+*/, 4324 /*Xi_c*+*/, 4332 /*Omega_c0*/, 4334 /*Omega_c*0*/,
1039 4412 /*Xi_cc+*/, 4414 /*Xi_cc*+*/, 4422 /*Xi_cc++*/, 4424 /*Xi_cc*0*/, 4432 /*Omega_cc+*/, 4434 /*Omega_cc*+*/,
1040
1041 4444 /*Omega_ccc++ not in evt.pdl*/
1042 };
1043
1044 constexpr std::array<int, 80> qqbarMesons = {
1045 // light qqbar
1046 111 /*pi0*/, 113 /*rho_0*/, 115 /*a_20*/, 117 /*rho(3)(1690)0*/, 119 /*a_4(1970)0*/,
1047 10111 /*pi(2S)0*/, 10113 /*b_10*/, 10115 /*pi(2)(1670)0*/,
1048 20113 /*a_10*/, 30113 /*rho(3S)0*/,
1049 100111 /*pi0_1300 not in evt.pdl*/, 100113 /*rho(2S)0*/,
1050 9000111 /*a_00*/, 9000113 /*pi(1)(1400)0*/, 9010113 /*pi(1)(1600)0*/, 9040113 /*rho(2150)0*/,
1051
1052 // not in evt.pdl
1053 200111 /*?*/, 9020113 /*a1_1640*/, 9030113 /*rho(1900)0*/, 100115 /*?*/, 9000115 /*?*/,
1054 9000117 /*rho(3)(1990)0*/, 9010117 /*rho(3)(2250)0*/,
1055
1056 // ssbar Mesons
1057 221 /*eta*/, 223 /*omega*/, 225 /*f_2*/, 227 /*omega(3)(1670)*/, 229 /*f(4)(2050)*/,
1058 331 /*eta'*/, 333 /*phi*/, 335 /*f'_2*/, 337 /*phi(3)(1850)*/,
1059 10223 /*h_1*/, 10225 /*eta(2)(1645)*/, 10331 /*f(0)(1710)*/, 10333 /*h'_1*/,
1060 20223 /*f_1*/, 20333 /*f'_1*/, 30223 /*omega(1650)*/,
1061 100223 /*omega(2S)*/, 100333 /*phi(1680)*/,
1062 9000221 /*sigma_0*/, 9010221 /*f_0*/, 9020221 /*eta(1405)*/, 9030221 /*f_0(1500)*/,
1063 9050225 /*f(2)(1950)*/, 9060225 /*f(2)(2010)*/,
1064
1065 // not in evt.pdl
1066 10221 /*f0(980)*/, 10335 /*eta(2)(1870)*/,
1067 100331 /*eta(1470)*/,
1068 9000223 /*f_1(1510)*/, 9000225 /*f_2(1430)*/, 9000229 /*f_J(2220)*/,
1069 9010223 /*h_1(1595)*/, 9010225 /*f_2(1565)*/, 9020225 /*f_2(1640)*/, 9030225 /*f_2(1810)*/,
1070 9040221 /*eta(1760)*/, 9040225 /*f_2(1910)*/, 9050221 /*f_0(2020)*/,
1071 9060221 /*f_0(2100)*/, 9070221 /*f_0(2200)*/, 9070225 /*f_2(2150)*/,
1072
1073 100225 /*?*/, 100335 /*?*/,
1074 200221 /*?*/,
1075 9000339 /*?*/,
1076
1077 // ccbar Mesons
1078 441 /*eta_c*/, 443 /*J/psi*/, 445 /*chi_c2*/,
1079 10441 /*chi_c0*/, 10443 /*h_c*/, 20443 /*chi_c1*/, 30443 /*psi(3770)*/,
1080 100441 /*eta_c(2S)*/, 100443 /*phi(2S)*/,
1081 9000443 /*phi(4040)*/, 9010443 /*psi(4160)*/, 9020443 /*psi(4415)*/,
1082 9000445 /*? not in evt.pdl*/,
1083 };
1084
1085 constexpr std::array<int, 63> flavorConservingMesons = {
1086 // Excited light mesons that can decay into hadrons conserving flavor
1087 213 /*rho+*/, 215 /*a_2+*/, 217 /*rho(3)(1690)+*/, 219 /*a_4(1970)+*/,
1088 10211 /*a(0)(1450)+*/, 10213 /*b_1+*/, 10215 /*pi(2)(1670)+*/,
1089 20213 /*a_1+*/, 30213 /*rho(3S)+*/,
1090 100211 /*pi(2S)+*/, 100213 /*rho(2S)+*/,
1091 9000211 /*a_0+*/, 9000213 /*pi(1)(1400)+*/, 9010213 /*pi(1)(1600)+*/,
1092
1093 // not in evt.pdl
1094 9000215 /*a_2(1700)+*/, 9000217 /*rho_3(1990)+*/,
1095 9010217 /*rho_3(2250)+*/, 9020213 /*a_1(1640)+*/, 9030213 /*rho(1900)+*/, 9040213 /*rho(2150)+*/,
1096
1097 100215 /*?*/, 200211 /*?*/,
1098
1099 // Excited K Mesons that hadronize conserving flavor
1100 313 /*K*0*/, 315 /*K_2*0*/, 317 /*K_3*0*/, 319 /*K_4*0*/,
1101 323 /*K*+*/, 325 /*K_2*+*/, 327 /*K_3*+*/, 329 /*K_4*+*/,
1102 10311 /*K_0*0*/, 10313 /*K_10*/, 10315 /*K(2)(1770)0*/,
1103 10321 /*K_0*+*/, 10323 /*K_1+*/, 10325 /*K(2)(1770)+*/,
1104 20313 /*K'_10*/, 20315 /*K(2)(1820)0*/,
1105 20323 /*K'_1+*/, 20325 /*K(2)(1820)+*/,
1106 30313 /*K''*0*/, 30323 /*K''*+*/, 30343 /*Xsd*/,
1107 100313 /*K'*0*/, 100323 /*K'*+*/,
1108
1109 // not in evt.pdl
1110 100311 /*K(1460)0*/, 100321 /*K(1460)+*/,
1111 9000311 /*K*(700)0*/, 9000313 /*K_1(1650)0*/, 9000315 /*K_2(1580)0*/,
1112 9000321 /*K*(700)+*/, 9000323 /*K_1(1650)+*/, 9000325 /*K_2(1580)+*/,
1113 9010315 /*K_2*(1980)0*/, 9010317 /*K_3(2320)0*/,
1114 9010325 /*K_2*(1980)+*/, 9010327 /*K_3(2320)+*/,
1115
1116 100315 /*?*/, 100325 /*?*/,
1117 200311 /*?*/, 200321 /*?*/,
1118 9000319 /*?*/, 9000329 /*?*/,
1119 };
1120
1121 const std::vector<std::string> availableForIsRightTrack = { "Electron", // 0
1122 "IntermediateElectron", // 1
1123 "Muon", // 2
1124 "IntermediateMuon", // 3
1125 "KinLepton", // 4
1126 "IntermediateKinLepton",// 5
1127 "Kaon", // 6
1128 "SlowPion", // 7
1129 "FastHadron", // 8
1130 "Lambda", // 9
1131 "mcAssociated" // 10
1132 };
1133
1134 const std::vector<std::string> availableForIsRightCategory = { "Electron", // 0
1135 "IntermediateElectron", // 1
1136 "Muon", // 2
1137 "IntermediateMuon", // 3
1138 "KinLepton", // 4
1139 "IntermediateKinLepton",// 5
1140 "Kaon", // 6
1141 "SlowPion", // 7
1142 "FastHadron", // 8
1143 "KaonPion", // 9
1144 "MaximumPstar", // 10
1145 "FSC", // 11
1146 "Lambda", // 12
1147 "mcAssociated" // 13
1148 };
1149
1150
1151
1152 Manager::FunctionPtr isRightTrack(const std::vector<std::string>& arguments)
1153 {
1154 if (arguments.size() != 1) {
1155 B2FATAL("Wrong number of arguments (1 required) for meta function isRightTrack");
1156 }
1157
1158 auto particleName = arguments[0];
1159
1160 unsigned index = std::find(availableForIsRightTrack.begin(), availableForIsRightTrack.end(), particleName)
1161 - availableForIsRightTrack.begin();
1162 if (index == availableForIsRightTrack.size()) {
1163 B2FATAL("isRightTrack: Not available category " << particleName <<
1164 ". The possibilities are Electron, IntermediateElectron, Muon, IntermediateMuon, KinLepton, IntermediateKinLepton, Kaon, SlowPion, FastHadron and Lambda");
1165 }
1166
1167 auto func = [index](const Particle * particle) -> int {
1168
1169 const MCParticle* mcParticle = particle->getMCParticle();
1170 if (!mcParticle) return -2;
1171
1172 int mcPDG = abs(mcParticle->getPDG());
1173
1174 // ---------------------------- Mothers and Grandmothers ----------------------------------
1175 std::vector<int> mothersPDG;
1176 std::vector<const MCParticle*> mothersPointers;
1177
1178 const MCParticle* mcMother = mcParticle->getMother();
1179 while (mcMother)
1180 {
1181 mothersPDG.push_back(abs(mcMother->getPDG()));
1182 mothersPointers.push_back(mcMother);
1183 if (abs(mcMother->getPDG()) == 511) break;
1184 mcMother = mcMother->getMother();
1185 }
1186
1187 if (mothersPDG.size() == 0) return -2;
1188
1189 //has associated mothers up to a B meson
1190 if (index == 10) return 1;
1191
1192 // ---------------- Is D Meson in the decay chain --------------------------------------
1193
1194 bool isCharmedMesonInChain = false;
1195 for (auto& iMCMotherPDG : mothersPDG)
1196 {
1197 if (std::find(charmMesons.begin(), charmMesons.end(), iMCMotherPDG) != charmMesons.end()) {
1198 isCharmedMesonInChain = true;
1199 break;
1200 }
1201 }
1202
1203 // ---------------- Is Charmed Baryon in the decay chain --------------------------------
1204
1205 bool isCharmedBaryonInChain = false;
1206 for (auto& iMCMotherPDG : mothersPDG)
1207 {
1208 if (std::find(charmBaryons.begin(), charmBaryons.end(), iMCMotherPDG) != charmBaryons.end()) {
1209 isCharmedBaryonInChain = true;
1210 break;
1211 }
1212 }
1213
1214 // ---------------- Is neutral qqbar Meson in the decay chain --------------------------------
1215
1216 bool isQQbarMesonInChain = false;
1217 for (auto& iMCMotherPDG : mothersPDG)
1218 {
1219 if (std::find(qqbarMesons.begin(), qqbarMesons.end(), iMCMotherPDG) != qqbarMesons.end()) {
1220 isQQbarMesonInChain = true;
1221 break;
1222 }
1223 }
1224
1225 // -------------- Is the Hadron a descendant of a Meson that conserves flavor --------------------------
1226
1227 bool isB0DaughterConservingFlavor = false;
1228 if (std::find(flavorConservingMesons.begin(), flavorConservingMesons.end(),
1229 mothersPDG.rbegin()[1]) != flavorConservingMesons.end())
1230 {
1231 isB0DaughterConservingFlavor = true;
1232 }
1233
1234 // ----------------------------- Is the Hadron a single daughter of a tau ----- --------------------------
1235
1236 bool isHadronSingleTauDaughter = false;
1237 if (mothersPDG.size() > 1 && mothersPDG.rbegin()[1] == 15)
1238 {
1239 int numberOfChargedDaughters = 0;
1240 for (auto& tauDaughter : mothersPointers.rbegin()[1]->getDaughters()) {
1241 if (tauDaughter->getCharge() != 0)
1242 numberOfChargedDaughters += 1;
1243 }
1244 if (numberOfChargedDaughters == 1)
1245 isHadronSingleTauDaughter = true;
1246 }
1247
1248 if (index == 0 // Electron
1249 && mcPDG == Const::electron.getPDGCode() && mothersPDG[0] == 511)
1250 {
1251 return 1;
1252 } else if (index == 1 // IntermediateElectron
1253 && mcPDG == Const::electron.getPDGCode() && mothersPDG.size() > 1 && isQQbarMesonInChain == false)
1254 {
1255 return 1;
1256 } else if (index == 2 // Muon
1257 && mcPDG == Const::muon.getPDGCode() && mothersPDG[0] == 511)
1258 {
1259 return 1;
1260 } else if (index == 3 // IntermediateMuon
1261 && mcPDG == Const::muon.getPDGCode() && mothersPDG.size() > 1 && isQQbarMesonInChain == false)
1262 {
1263 return 1;
1264 } else if (index == 4 // KinLepton
1265 && (mcPDG == Const::muon.getPDGCode() || mcPDG == Const::electron.getPDGCode()) && mothersPDG[0] == 511)
1266 {
1267 return 1;
1268 } else if (index == 5 //IntermediateKinLepton
1269 && (mcPDG == Const::muon.getPDGCode() || mcPDG == Const::electron.getPDGCode()) && mothersPDG.size() > 1
1270 && isQQbarMesonInChain == false)
1271 {
1272 return 1;
1273 } else if (index == 6 // Kaon
1274 && mcPDG == Const::kaon.getPDGCode() && isQQbarMesonInChain == false
1275 && (isCharmedMesonInChain == true || isCharmedBaryonInChain == true))
1276 {
1277 return 1;
1278 } else if (index == 7 // SlowPion
1279 && mcPDG == Const::pion.getPDGCode() && mothersPDG.size() > 1 && mothersPDG[0] == 413 && mothersPDG[1] == 511)
1280 {
1281 return 1;
1282 } else if (index == 8 // FastHadron
1283 && (mcPDG == Const::pion.getPDGCode() || mcPDG == Const::kaon.getPDGCode())
1284 && isQQbarMesonInChain == false
1285 && (mothersPDG[0] == 511
1286 || (mothersPDG.rbegin()[0] == 511 && (isB0DaughterConservingFlavor == true || isHadronSingleTauDaughter == true))))
1287 {
1288 return 1;
1289 } else if (index == 9 // Lambda
1290 && mcPDG == Const::Lambda.getPDGCode() && isCharmedBaryonInChain == true)
1291 {
1292 return 1;
1293 } else
1294 {
1295 return 0;
1296 }
1297 };
1298 return func;
1299 }
1300
1301 Manager::FunctionPtr isRightCategory(const std::vector<std::string>& arguments)
1302 {
1303 if (arguments.size() != 1) {
1304 B2FATAL("Wrong number of arguments (1 required) for meta function isRightCategory");
1305 }
1306
1307 auto particleName = arguments[0];
1308
1309 unsigned index = find(availableForIsRightCategory.begin(), availableForIsRightCategory.end(), particleName)
1310 - availableForIsRightCategory.begin();
1311 if (index == availableForIsRightCategory.size()) {
1312 B2FATAL("isRightCategory: Not available category " << particleName <<
1313 ". The possibilities are Electron, IntermediateElectron, Muon, IntermediateMuon, KinLepton, IntermediateKinLepton, Kaon, SlowPion, FastHadron, KaonPion, MaximumPstar, FSC and Lambda");
1314 }
1315
1316 auto func = [index](const Particle * particle) -> int {
1317
1318 Particle* nullParticle = nullptr;
1319 double qTarget = particle->getCharge();
1320 double qMC = Variable::isRestOfEventB0Flavor(nullParticle);
1321
1322 const MCParticle* mcParticle = particle->getMCParticle();
1323 if (!mcParticle) return -2;
1324
1325 int mcPDG = abs(mcParticle->getPDG());
1326
1327 // ---------------------------- Mothers and Grandmothers ---------------------------------
1328 std::vector<int> mothersPDG;
1329 std::vector<const MCParticle*> mothersPointers;
1330
1331 const MCParticle* mcMother = mcParticle->getMother();
1332 while (mcMother)
1333 {
1334 mothersPDG.push_back(abs(mcMother->getPDG()));
1335 mothersPointers.push_back(mcMother);
1336 if (abs(mcMother->getPDG()) == 511) break;
1337 mcMother = mcMother->getMother();
1338 }
1339
1340 if (mothersPDG.size() == 0) return -2;
1341 //has associated mothers up to a B meson
1342 if (index == 13) return 1;
1343
1344 // ---------------- Is D Meson in the decay chain --------------------------------------
1345
1346 bool isCharmedMesonInChain = false;
1347 for (auto& iMCMotherPDG : mothersPDG)
1348 {
1349 if (std::find(charmMesons.begin(), charmMesons.end(), iMCMotherPDG) != charmMesons.end()) {
1350 isCharmedMesonInChain = true;
1351 break;
1352 }
1353 }
1354
1355 // ---------------- Is Charmed Baryon in the decay chain --------------------------------
1356
1357 bool isCharmedBaryonInChain = false;
1358 for (auto& iMCMotherPDG : mothersPDG)
1359 {
1360 if (std::find(charmBaryons.begin(), charmBaryons.end(), iMCMotherPDG) != charmBaryons.end()) {
1361 isCharmedBaryonInChain = true;
1362 break;
1363 }
1364 }
1365
1366 // ---------------- Is neutral qqbar Meson in the decay chain --------------------------------
1367
1368 bool isQQbarMesonInChain = false;
1369 for (auto& iMCMotherPDG : mothersPDG)
1370 {
1371 if (std::find(qqbarMesons.begin(), qqbarMesons.end(), iMCMotherPDG) != qqbarMesons.end()) {
1372 isQQbarMesonInChain = true;
1373 break;
1374 }
1375 }
1376
1377 // -------------- Is the Hadron a descendant of a Meson that conserves flavor --------------------------
1378
1379 bool isB0DaughterConservingFlavor = false;
1380 if (mothersPDG.size() > 1)
1381 {
1382 if (std::find(flavorConservingMesons.begin(), flavorConservingMesons.end(),
1383 mothersPDG.rbegin()[1]) != flavorConservingMesons.end()) {
1384 isB0DaughterConservingFlavor = true;
1385 }
1386 }
1387
1388 // ----------------------------- Is the Hadron a single daughter of a tau ----- --------------------------
1389
1390 bool isHadronSingleTauDaughter = false;
1391 if (mothersPDG.size() > 1 && mothersPDG.rbegin()[1] == 15)
1392 {
1393 int numberOfChargedDaughters = 0;
1394 for (auto& tauDaughter : mothersPointers.rbegin()[1]->getDaughters()) {
1395 if (tauDaughter->getCharge() != 0)
1396 numberOfChargedDaughters += 1;
1397 }
1398 if (numberOfChargedDaughters == 1)
1399 isHadronSingleTauDaughter = true;
1400 }
1401
1402 // ---------------------------- For KaonPion Category ------------------------------------
1403 bool haveKaonPionSameMother = false;
1404 if (index == 9) // KaonPion
1405 {
1406 const MCParticle* mcSlowPionMother = nullptr;
1407 StoreObjPtr<ParticleList> SlowPionList("pi+:inRoe");
1408 Particle* targetSlowPion = nullptr;
1409 if (SlowPionList.isValid()) {
1410 double mcProbSlowPion = 0;
1411 for (unsigned int i = 0; i < SlowPionList->getListSize(); ++i) {
1412 Particle* pSlowPion = SlowPionList->getParticle(i);
1413 if (!pSlowPion) continue;
1414 if (pSlowPion->hasExtraInfo("isRightCategory(SlowPion)")) {
1415 double probSlowPion = pSlowPion->getExtraInfo("isRightCategory(SlowPion)");
1416 if (probSlowPion > mcProbSlowPion) {
1417 mcProbSlowPion = probSlowPion;
1418 targetSlowPion = pSlowPion;
1419 }
1420 }
1421 }
1422 if (targetSlowPion != nullptr) {
1423 const MCParticle* mcSlowPion = targetSlowPion ->getMCParticle();
1424 // SlowPion_q = targetSlowPion->getCharge();
1425 if (mcSlowPion != nullptr && mcSlowPion->getMother() != nullptr
1426 && abs(mcSlowPion->getPDG()) == Const::pion.getPDGCode() && abs(mcSlowPion->getMother()->getPDG()) == 413) {
1427 mcSlowPionMother = mcSlowPion->getMother();
1428 }
1429 }
1430 }
1431
1432 if (std::find(mothersPointers.begin(), mothersPointers.end(), mcSlowPionMother) != mothersPointers.end())
1433 haveKaonPionSameMother = true;
1434
1435 }
1436
1437 // ---------------------------- For FastSlowCorrelated Category ----------------------------
1438 int FastParticlePDGMother = 0;
1439 double qFSC = 0;
1440 if (index == 11) // FSC
1441 {
1442 StoreObjPtr<ParticleList> FastParticleList("pi+:inRoe");
1443 Particle* targetFastParticle = nullptr;
1444 if (FastParticleList.isValid()) {
1445 double mcProbFastest = 0;
1446 for (unsigned int i = 0; i < FastParticleList->getListSize(); ++i) {
1447 Particle* particlei = FastParticleList->getParticle(i);
1448 if (!particlei) continue;
1449
1450 ROOT::Math::PxPyPzEVector momParticlei = PCmsLabTransform::labToCms(particlei->get4Vector());
1451 if (momParticlei == momParticlei) {
1452 double probFastest = momParticlei.P();
1453 if (probFastest > mcProbFastest) {
1454 mcProbFastest = probFastest;
1455 targetFastParticle = particlei;
1456 }
1457 }
1458 }
1459 if (targetFastParticle != nullptr) {
1460 const MCParticle* mcFastParticle = targetFastParticle ->getMCParticle();
1461 // FastParticle_q = targetFastParticle->getCharge();
1462 if (mcFastParticle != nullptr && mcFastParticle->getMother() != nullptr) {
1463 FastParticlePDGMother = abs(mcFastParticle->getMother()->getPDG());
1464 qFSC = mcFastParticle->getCharge();
1465 }
1466 }
1467 }
1468 }
1469
1470 // ------------------------------ Outputs -----------------------------------
1471 if (index == 0 // Electron
1472 && qTarget == qMC && mcPDG == Const::electron.getPDGCode() && mothersPDG[0] == 511)
1473 {
1474 return 1;
1475 } else if (index == 1 // IntermediateElectron
1476 && qTarget != qMC && mcPDG == Const::electron.getPDGCode() && mothersPDG.size() > 1
1477 && isQQbarMesonInChain == false)
1478 {
1479 return 1;
1480 } else if (index == 2 // Muon
1481 && qTarget == qMC && mcPDG == Const::muon.getPDGCode() && mothersPDG[0] == 511)
1482 {
1483 return 1;
1484 } else if (index == 3 // IntermediateMuon
1485 && qTarget != qMC && mcPDG == Const::muon.getPDGCode() && mothersPDG.size() > 1
1486 && isQQbarMesonInChain == false)
1487 {
1488 return 1;
1489 } else if (index == 4 // KinLepton
1490 && qTarget == qMC
1491 && (mcPDG == Const::electron.getPDGCode() || mcPDG == Const::muon.getPDGCode()) && mothersPDG[0] == 511)
1492 {
1493 return 1;
1494 } else if (index == 5 // IntermediateKinLepton
1495 && qTarget != qMC && (mcPDG == Const::electron.getPDGCode() || mcPDG == Const::muon.getPDGCode())
1496 && mothersPDG.size() > 1 && isQQbarMesonInChain == false)
1497 {
1498 return 1;
1499 } else if (index == 6// Kaon
1500 && qTarget == qMC && mcPDG == Const::kaon.getPDGCode() && isQQbarMesonInChain == false
1501 && (isCharmedMesonInChain == true || isCharmedBaryonInChain == true))
1502 {
1503 return 1;
1504 } else if (index == 7 // SlowPion
1505 && qTarget != qMC && mcPDG == Const::pion.getPDGCode()
1506 && mothersPDG.size() > 1 && mothersPDG[0] == 413 && mothersPDG[1] == 511)
1507 {
1508 return 1;
1509 } else if (index == 8 // FastHadron
1510 && qTarget == qMC && (mcPDG == Const::pion.getPDGCode() || mcPDG == Const::kaon.getPDGCode())
1511 && isQQbarMesonInChain == false
1512 && (mothersPDG[0] == 511 || (mothersPDG.rbegin()[0] == 511
1513 && (isB0DaughterConservingFlavor == true || isHadronSingleTauDaughter == true))))
1514 {
1515 return 1;
1516 } else if (index == 9 // KaonPion
1517 && qTarget == qMC && mcPDG == Const::kaon.getPDGCode() && haveKaonPionSameMother == true)
1518 {
1519 return 1;
1520 } else if (index == 10 && qTarget == qMC) // MaximumPstar
1521 {
1522 return 1;
1523 } else if (index == 11 // FSC
1524 && qTarget != qMC && mothersPDG.size() > 1 && qFSC == qMC
1525 && mcPDG == Const::pion.getPDGCode() && FastParticlePDGMother == 511 && isQQbarMesonInChain == false)
1526 {
1527 return 1;
1528 } else if (index == 12 // Lambda
1529 && (particle->getPDGCode() / abs(particle->getPDGCode())) != qMC
1530 && mcPDG == Const::Lambda.getPDGCode() && isCharmedBaryonInChain == true)
1531 {
1532 return 1;
1533 } else
1534 {
1535 return 0;
1536 }
1537 };
1538
1539 return func;
1540 }
1541
1542 Manager::FunctionPtr hasHighestProbInCat(const std::vector<std::string>& arguments)
1543 {
1544 if (arguments.size() != 2) {
1545 B2FATAL("Wrong number of arguments (2 required) for meta function hasHighestProbInCat");
1546 }
1547
1548 auto particleListName = arguments[0];
1549 auto extraInfoName = arguments[1];
1550
1551 bool isAvailable = false;
1552 for (const auto& name : availableForIsRightTrack) {
1553 if (extraInfoName == "isRightTrack(" + name + ")") {
1554 isAvailable = true;
1555 break;
1556 }
1557 }
1558 for (const auto& name : availableForIsRightCategory) {
1559 if (extraInfoName == "isRightCategory(" + name + ")") {
1560 isAvailable = true;
1561 break;
1562 }
1563 }
1564 if (extraInfoName == "isRightTrack(MaximumPstar)")
1565 isAvailable = true;
1566
1567
1568 if (!isAvailable) {
1569 string strAvailableForIsRightTrack;
1570 for (const auto& name : availableForIsRightTrack)
1571 strAvailableForIsRightTrack += name + " ";
1572 string strAvailableForIsRightCategory;
1573 for (const auto& name : availableForIsRightCategory)
1574 strAvailableForIsRightCategory += name + " ";
1575
1576 B2FATAL("hasHighestProbInCat: Not available category" << extraInfoName <<
1577 ". The possibilities for isRightTrack() are " << endl << strAvailableForIsRightTrack << " MaximumPstar" << endl <<
1578 "The possibilities for isRightCategory() are " << endl << strAvailableForIsRightCategory);
1579 }
1580
1581 auto func = [particleListName, extraInfoName](const Particle * particle) -> bool {
1582 StoreObjPtr<ParticleList> ListOfParticles(particleListName);
1583 if (!ListOfParticles.isValid()) return 0;
1584
1585 double maximumProb = 0;
1586 for (unsigned int i = 0; i < ListOfParticles->getListSize(); ++i)
1587 {
1588 const Particle* particlei = ListOfParticles->getParticle(i);
1589 if (!particlei) continue;
1590
1591 double prob = 0;
1592 if (extraInfoName == "isRightTrack(MaximumPstar)") {
1593 ROOT::Math::PxPyPzEVector momParticlei = PCmsLabTransform::labToCms(particlei->get4Vector());
1594 if (momParticlei == momParticlei) {
1595 prob = momParticlei.P();
1596 }
1597 } else {
1598 if (particlei->hasExtraInfo(extraInfoName)) {
1599 prob = particlei->getExtraInfo(extraInfoName);
1600 }
1601 }
1602 if (prob > maximumProb) {
1603 maximumProb = prob;
1604 }
1605
1606 }
1607
1608 bool output = false;
1609 if ((extraInfoName == "isRightTrack(MaximumPstar)") && (PCmsLabTransform::labToCms(particle->get4Vector()).P() == maximumProb))
1610 {
1611 output = true;
1612 } else if (extraInfoName != "isRightTrack(MaximumPstar)" && particle->hasExtraInfo(extraInfoName))
1613 {
1614 if (particle->getExtraInfo(extraInfoName) == maximumProb) output = true;
1615 }
1616
1617 return output;
1618 };
1619 return func;
1620 }
1621
1622 Manager::FunctionPtr HighestProbInCat(const std::vector<std::string>& arguments)
1623 {
1624 if (arguments.size() != 2) {
1625 B2FATAL("Wrong number of arguments (2 required) for meta function HighestProbInCat");
1626 }
1627
1628 auto particleListName = arguments[0];
1629 auto extraInfoName = arguments[1];
1630
1631 bool isAvailable = false;
1632 for (const auto& name : availableForIsRightTrack) {
1633 if (extraInfoName == "isRightTrack(" + name + ")") {
1634 isAvailable = true;
1635 break;
1636 }
1637 }
1638 for (const auto& name : availableForIsRightCategory) {
1639 if (extraInfoName == "isRightCategory(" + name + ")") {
1640 isAvailable = true;
1641 break;
1642 }
1643 }
1644 if (extraInfoName == "isRightTrack(MaximumPstar)")
1645 isAvailable = true;
1646
1647
1648 if (!isAvailable) {
1649 string strAvailableForIsRightTrack;
1650 for (const auto& name : availableForIsRightTrack)
1651 strAvailableForIsRightTrack += name + " ";
1652 string strAvailableForIsRightCategory;
1653 for (const auto& name : availableForIsRightCategory)
1654 strAvailableForIsRightCategory += name + " ";
1655
1656 B2FATAL("HighestProbInCat: Not available category" << extraInfoName <<
1657 ". The possibilities for isRightTrack() are " << endl << strAvailableForIsRightTrack << " MaximumPstar" << endl <<
1658 "The possibilities for isRightCategory() are " << endl << strAvailableForIsRightCategory);
1659 }
1660
1661 auto func = [particleListName, extraInfoName](const Particle*) -> double {
1662 StoreObjPtr<ParticleList> ListOfParticles(particleListName);
1663 if (!ListOfParticles.isValid()) return 0;
1664
1665 double maximumProb = 0;
1666 for (unsigned int i = 0; i < ListOfParticles->getListSize(); ++i)
1667 {
1668 const Particle* particlei = ListOfParticles->getParticle(i);
1669 if (!particlei) continue;
1670
1671 double prob = 0;
1672 if (extraInfoName == "isRightTrack(MaximumPstar)") {
1673 ROOT::Math::PxPyPzEVector momParticlei = PCmsLabTransform::labToCms(particlei->get4Vector());
1674 if (momParticlei == momParticlei) {
1675 prob = momParticlei.P();
1676 }
1677 } else {
1678 if (particlei->hasExtraInfo(extraInfoName)) {
1679 prob = particlei->getExtraInfo(extraInfoName);
1680 }
1681 }
1682 maximumProb = max(maximumProb, prob);
1683 }
1684
1685 return maximumProb;
1686 };
1687 return func;
1688 }
1689
1690
1691 // List of available extrainfos used in QpOf, weightedQpOf and variableOfTarget.
1692 const std::vector<std::string> availableExtraInfos = {"isRightTrack(Electron)", // 0
1693 "isRightTrack(IntermediateElectron)", // 1
1694 "isRightTrack(Muon)", // 2
1695 "isRightTrack(IntermediateMuon)", // 3
1696 "isRightTrack(KinLepton)", // 4
1697 "isRightTrack(IntermediateKinLepton)",// 5
1698 "isRightTrack(Kaon)", // 6
1699 "isRightTrack(SlowPion)", // 7
1700 "isRightTrack(FastHadron)", // 8
1701 "isRightTrack(MaximumPstar)", // 9
1702 "isRightTrack(Lambda)", // 10
1703 "isRightCategory(Electron)", // 11
1704 "isRightCategory(IntermediateElectron)", // 12
1705 "isRightCategory(Muon)", // 13
1706 "isRightCategory(IntermediateMuon)", // 14
1707 "isRightCategory(KinLepton)", // 15
1708 "isRightCategory(IntermediateKinLepton)",// 16
1709 "isRightCategory(Kaon)", // 17
1710 "isRightCategory(SlowPion)", // 18
1711 "isRightCategory(FastHadron)", // 19
1712 "isRightCategory(MaximumPstar)", // 20
1713 "isRightCategory(Lambda)", // 21
1714 "isRightCategory(KaonPion)", // 22
1715 "isRightCategory(FSC)", // 23
1716 };
1717
1718 Manager::FunctionPtr QpOf(const std::vector<std::string>& arguments)
1719 {
1720 if (arguments.size() != 3) {
1721 B2FATAL("Wrong number of arguments (3 required) for meta function QpOf");
1722 }
1723
1724 auto particleListName = arguments[0];
1725 auto outputExtraInfo = arguments[1];
1726 auto rankingExtraInfo = arguments[2];
1727
1728 unsigned indexRanking = find(availableExtraInfos.begin(), availableExtraInfos.end(),
1729 rankingExtraInfo) - availableExtraInfos.begin();
1730 unsigned indexOutput = find(availableExtraInfos.begin(), availableExtraInfos.end(),
1731 outputExtraInfo) - availableExtraInfos.begin();
1732
1733 if (indexRanking == availableExtraInfos.size() or indexOutput == availableExtraInfos.size()) {
1734 string strAvailableForIsRightTrack;
1735 for (const auto& name : availableForIsRightTrack)
1736 strAvailableForIsRightTrack += name + " ";
1737 string strAvailableForIsRightCategory;
1738 for (const auto& name : availableForIsRightCategory)
1739 strAvailableForIsRightCategory += name + " ";
1740
1741 B2FATAL("QpOf: Not available category " << rankingExtraInfo <<
1742 ". The possibilities for isRightTrack() are " << endl << strAvailableForIsRightTrack << " MaximumPstar" << endl <<
1743 "The possibilities for isRightCategory() are " << endl << strAvailableForIsRightCategory);
1744 }
1745
1746 auto func = [particleListName, indexOutput, indexRanking](const Particle*) -> double {
1747 StoreObjPtr<ParticleList> ListOfParticles(particleListName);
1748 if (!ListOfParticles.isValid()) return 0;
1749
1750 Particle* target = nullptr; //Particle selected as target
1751 double maximumTargetProb = 0; //Probability of being the target track from the track level
1752 for (unsigned int i = 0; i < ListOfParticles->getListSize(); ++i)
1753 {
1754 Particle* particlei = ListOfParticles->getParticle(i);
1755 if (!particlei) continue;
1756
1757 double target_prob = 0;
1758 if (indexRanking == 9 || indexRanking == 20) { // MaximumPstar
1759 ROOT::Math::PxPyPzEVector momParticlei = PCmsLabTransform::labToCms(particlei->get4Vector());
1760 if (momParticlei == momParticlei) {
1761 target_prob = momParticlei.P();
1762 }
1763 } else {
1764 if (particlei->hasExtraInfo(availableExtraInfos[indexRanking])) {
1765 target_prob = particlei->getExtraInfo(availableExtraInfos[indexRanking]);
1766 }
1767 }
1768
1769 if (target_prob > maximumTargetProb) {
1770 maximumTargetProb = target_prob;
1771 target = particlei;
1772 }
1773 }
1774
1775 // nothing found
1776 if (!target) return 0;
1777
1778 double qTarget = 0; //Flavor of the track selected as target
1779 // Get the flavor of the track selected as target
1780 if (indexRanking == 10 || indexRanking == 21) // Lambda
1781 {
1782 qTarget = (-1) * target->getPDGCode() / abs(target->getPDGCode());
1783 // IntermediateElectron IntermediateMuon IntermediateKinLepton SlowPion
1784 } else if (indexRanking == 1 || indexRanking == 3 || indexRanking == 5 || indexRanking == 7 ||
1785 indexRanking == 12 || indexRanking == 14 || indexRanking == 16 || indexRanking == 18)
1786 {
1787 qTarget = (-1) * target->getCharge();
1788 } else
1789 {
1790 qTarget = target->getCharge();
1791 }
1792
1793 //Get the probability of being right classified flavor from event level
1794 double prob = target->getExtraInfo(availableExtraInfos[indexOutput]);
1795
1796 //float r = abs(2 * prob - 1); //Definition of the dilution factor */
1797 //return 0.5 * (qTarget * r + 1);
1798 return qTarget * prob;
1799 };
1800 return func;
1801 }
1802
1803 Manager::FunctionPtr weightedQpOf(const std::vector<std::string>& arguments)
1804 {
1805 if (arguments.size() != 3) {
1806 B2FATAL("Wrong number of arguments (3 required) for meta function weightedQpOf");
1807 }
1808
1809 //used by simple_flavor_tagger
1810
1811 auto particleListName = arguments[0];
1812 auto outputExtraInfo = arguments[1];
1813 auto rankingExtraInfo = arguments[2];
1814
1815
1816 unsigned indexRanking = find(availableExtraInfos.begin(), availableExtraInfos.end(),
1817 rankingExtraInfo) - availableExtraInfos.begin();
1818 unsigned indexOutput = find(availableExtraInfos.begin(), availableExtraInfos.end(),
1819 outputExtraInfo) - availableExtraInfos.begin();
1820
1821
1822 if (indexRanking == availableExtraInfos.size() or indexOutput == availableExtraInfos.size()) {
1823 string strAvailableForIsRightTrack;
1824 for (const auto& name : availableForIsRightTrack)
1825 strAvailableForIsRightTrack += name + " ";
1826 string strAvailableForIsRightCategory;
1827 for (const auto& name : availableForIsRightCategory)
1828 strAvailableForIsRightCategory += name + " ";
1829
1830 B2FATAL("weightedQpOf: Not available category " << rankingExtraInfo <<
1831 ". The possibilities for isRightTrack() are " << endl << strAvailableForIsRightTrack << " MaximumPstar" << endl <<
1832 "The possibilities for isRightCategory() are " << endl << strAvailableForIsRightCategory);
1833 }
1834
1835
1836 auto func = [particleListName, indexOutput, indexRanking, rankingExtraInfo](const Particle*) -> double {
1837
1838 StoreObjPtr<ParticleList> ListOfParticles(particleListName);
1839 if (!ListOfParticles) return 0;
1840 if (ListOfParticles->getListSize() == 0) return 0;
1841
1842
1843 auto compare = [rankingExtraInfo](const Particle * part1, const Particle * part2)-> bool {
1844 double info1 = 0;
1845 double info2 = 0;
1846 if (part1->hasExtraInfo(rankingExtraInfo)) info1 = part1->getExtraInfo(rankingExtraInfo);
1847 if (part2->hasExtraInfo(rankingExtraInfo)) info2 = part2->getExtraInfo(rankingExtraInfo);
1848 return (info1 > info2);
1849 };
1850
1851 auto compareMomentum = [rankingExtraInfo](const Particle * part1, const Particle * part2)-> bool {
1852 double info1 = PCmsLabTransform::labToCms(part1->get4Vector()).P();
1853 double info2 = PCmsLabTransform::labToCms(part2->get4Vector()).P();
1854 return (info1 > info2);
1855 };
1856
1857 std::vector<const Particle*> ParticleVector(ListOfParticles->getListSize());
1858 for (unsigned int i = 0; i < ListOfParticles->getListSize(); ++i)
1859 {
1860 ParticleVector[i] = ListOfParticles->getParticle(i);
1861 }
1862
1863 if (indexRanking == 9 || indexRanking == 20)
1864 std::sort(ParticleVector.begin(), ParticleVector.end(), compareMomentum); // MaximumPstar
1865 else
1866 std::sort(ParticleVector.begin(), ParticleVector.end(), compare);
1867
1868
1869 double final_value = 0.0;
1870 if (ParticleVector.size() != 0) final_value = 1.0;
1871
1872 //Loop over K+ vector until 3 or empty
1873 int Limit = min(3, int(ParticleVector.size()));
1874 double val1 = 1.0;
1875 double val2 = 1.0;
1876 for (int i = 0; i < Limit; ++i)
1877 {
1878 if (ParticleVector[i]->hasExtraInfo(availableExtraInfos[indexOutput])) {
1879 double flavor = 0.0;
1880 if (indexRanking == 10 || indexRanking == 21) { // Lambda
1881 flavor = - copysign(1, ParticleVector[i]->getPDGCode());
1882 // IntermediateElectron IntermediateMuon IntermediateKinLepton SlowPion
1883 } else if (indexRanking == 1 || indexRanking == 3 || indexRanking == 5 || indexRanking == 7 ||
1884 indexRanking == 12 || indexRanking == 14 || indexRanking == 16 || indexRanking == 18) {
1885 flavor = - ParticleVector[i]->getCharge();
1886 } else {
1887 flavor = + ParticleVector[i]->getCharge();
1888 }
1889
1890 double p = ParticleVector[i]->getExtraInfo(availableExtraInfos[indexOutput]);
1891 // B2INFO("Right Track:" << ParticleVector[i]->getExtraInfo(availableExtraInfos[indexRanking]));
1892 // B2INFO("Right Cat:" << ParticleVector[i]->getExtraInfo(availableExtraInfos[indexOutput]));
1893 double qp = (flavor * p);
1894 val1 *= 1 + qp;
1895 val2 *= 1 - qp;
1896 }
1897 }
1898 final_value = (val1 - val2) / (val1 + val2);
1899
1900 return final_value;
1901 };
1902 return func;
1903 }
1904
1905
1906 Manager::FunctionPtr QpTrack(const std::vector<std::string>& arguments)
1907 {
1908 if (arguments.size() != 3) {
1909 B2FATAL("Wrong number of arguments (3 required) for meta function QpTrack");
1910 }
1911
1912 auto particleListName = arguments[0];
1913 auto outputExtraInfo = arguments[1];
1914 auto rankingExtraInfo = arguments[2];
1915
1916 unsigned indexRanking = find(availableExtraInfos.begin(), availableExtraInfos.end(),
1917 rankingExtraInfo) - availableExtraInfos.begin();
1918 unsigned indexOutput = find(availableExtraInfos.begin(), availableExtraInfos.end(),
1919 outputExtraInfo) - availableExtraInfos.begin();
1920
1921 if (indexRanking == availableExtraInfos.size() or indexOutput == availableExtraInfos.size()) {
1922 string strAvailableForIsRightTrack;
1923 for (const auto& name : availableForIsRightTrack)
1924 strAvailableForIsRightTrack += name + " ";
1925 string strAvailableForIsRightCategory;
1926 for (const auto& name : availableForIsRightCategory)
1927 strAvailableForIsRightCategory += name + " ";
1928
1929 B2FATAL("QpTrack: Not available category " << rankingExtraInfo <<
1930 ". The possibilities for isRightTrack() are " << endl << strAvailableForIsRightTrack << " MaximumPstar" << endl <<
1931 "The possibilities for isRightCategory() are " << endl << strAvailableForIsRightCategory);
1932 }
1933
1934 auto func = [particleListName, indexOutput, indexRanking](const Particle * particle) -> double {
1935 StoreObjPtr<ParticleList> ListOfParticles(particleListName);
1936 if (!ListOfParticles.isValid()) return 0;
1937
1938 const auto mdstSource = particle->getMdstSource();
1939
1940 Particle* target = nullptr; //Particle selected as target
1941 for (unsigned int i = 0; i < ListOfParticles->getListSize(); ++i)
1942 {
1943 Particle* particlei = ListOfParticles->getParticle(i);
1944 if (!particlei)
1945 continue;
1946
1947 if (particlei->getMdstSource() == mdstSource) {
1948 target = particlei;
1949 break;
1950 }
1951 }
1952
1953 // nothing found
1954 if (!target) return 0;
1955
1956 double qTarget = 0; //Flavor of the track selected as target
1957 // Get the flavor of the track selected as target
1958 if (indexRanking == 10 || indexRanking == 21) // Lambda
1959 {
1960 qTarget = (-1) * target->getPDGCode() / abs(target->getPDGCode());
1961 // IntermediateElectron IntermediateMuon IntermediateKinLepton SlowPion
1962 } else if (indexRanking == 1 || indexRanking == 3 || indexRanking == 5 || indexRanking == 7 ||
1963 indexRanking == 12 || indexRanking == 14 || indexRanking == 16 || indexRanking == 18)
1964 {
1965 qTarget = (-1) * target->getCharge();
1966 } else
1967 {
1968 qTarget = target->getCharge();
1969 }
1970
1971 //Get the probability of being right classified flavor from event level
1972 double prob = target->getExtraInfo(availableExtraInfos[indexOutput]);
1973
1974 return qTarget * prob;
1975 };
1976 return func;
1977 }
1978
1979 Manager::FunctionPtr variableOfTarget(const std::vector<std::string>& arguments)
1980 {
1981
1982 if (arguments.size() != 3)
1983 B2FATAL("Wrong number of arguments (3 required) for meta function variableOfTarget");
1984
1985 std::string particleListName = arguments[0];
1986 std::string inputVariable = arguments[1];
1987 std::string rankingExtraInfo = arguments[2];
1988
1989 int indexRanking = -1;
1990
1991 for (unsigned i = 0; i < availableExtraInfos.size(); ++i) {
1992 if (rankingExtraInfo == availableExtraInfos[i]) {indexRanking = i; break;}
1993 }
1994
1995 if (indexRanking == -1) {
1996 string strAvailableForIsRightTrack;
1997 for (const auto& name : availableForIsRightTrack)
1998 strAvailableForIsRightTrack += name + " ";
1999 string strAvailableForIsRightCategory;
2000 for (const auto& name : availableForIsRightCategory)
2001 strAvailableForIsRightCategory += name + " ";
2002
2003 B2FATAL("variableOfTarget: Not available category " << rankingExtraInfo <<
2004 ". The possibilities for isRightTrack() are " << endl << strAvailableForIsRightTrack << " MaximumPstar" << endl <<
2005 "The possibilities for isRightCategory() are " << endl << strAvailableForIsRightCategory);
2006 }
2007
2008 const Variable::Manager::Var* var = Manager::Instance().getVariable(inputVariable);
2009 auto func = [particleListName, var, indexRanking](const Particle*) -> double {
2010 StoreObjPtr<ParticleList> ListOfParticles(particleListName);
2011 if (!ListOfParticles.isValid()) return Const::doubleNaN;
2012
2013 Particle* target = nullptr; //Particle selected as target
2014
2015 double maximumTargetProb = 0; //Probability of being the target track from the track level
2016 for (unsigned int i = 0; i < ListOfParticles->getListSize(); ++i)
2017 {
2018 Particle* particlei = ListOfParticles->getParticle(i);
2019 if (!particlei) continue;
2020
2021 double target_prob = 0;
2022 if (indexRanking == 9 || indexRanking == 20) { // MaximumPstar
2023 ROOT::Math::PxPyPzEVector momParticlei = PCmsLabTransform::labToCms(particlei->get4Vector());
2024 if (momParticlei == momParticlei) {
2025 target_prob = momParticlei.P();
2026 }
2027 } else {
2028 if (particlei->hasExtraInfo(availableExtraInfos[indexRanking])) {
2029 target_prob = particlei->getExtraInfo(availableExtraInfos[indexRanking]);
2030 }
2031 }
2032 if (target_prob > maximumTargetProb) {
2033 maximumTargetProb = target_prob;
2034 target = particlei;
2035 }
2036 }
2037
2038 // no target found
2039 if (!target) return Const::doubleNaN;
2040
2041 if (std::holds_alternative<double>(var->function(target)))
2042 {
2043 return std::get<double>(var->function(target));
2044 } else if (std::holds_alternative<int>(var->function(target)))
2045 {
2046 return std::get<int>(var->function(target));
2047 } else if (std::holds_alternative<bool>(var->function(target)))
2048 {
2049 return std::get<bool>(var->function(target));
2050 } else return Const::doubleNaN;
2051 };
2052 return func;
2053 }
2054
2055 Manager::FunctionPtr hasTrueTarget(const std::vector<std::string>& arguments)
2056 {
2057 if (arguments.size() != 1) {
2058 B2FATAL("Wrong number of arguments (1 required) for meta function hasTrueTarget");
2059 }
2060
2061 auto categoryName = arguments[0];
2062
2063 bool isAvailable = false;
2064 for (const auto& name : availableForIsRightCategory) {
2065 if (categoryName == name) {
2066 isAvailable = true;
2067 break;
2068 }
2069 }
2070 if (categoryName == "mcAssociated")
2071 isAvailable = false;
2072
2073 if (!isAvailable) {
2074 string strAvailableForIsRightCategory;
2075 for (const auto& name : availableForIsRightCategory) {
2076 if (name == "mcAssociated") continue;
2077 strAvailableForIsRightCategory += name + " ";
2078 }
2079 B2FATAL("hasTrueTarget: Not available category" << categoryName <<
2080 ". The possibilities for the category name are " << endl << strAvailableForIsRightCategory);
2081 }
2082
2083 auto func = [categoryName](const Particle*) -> double {
2084 std::string particleListName;
2085 std::string trackTargetName = categoryName;
2086
2087 if (categoryName == "Electron" || categoryName == "IntermediateElectron") particleListName = "e+:inRoe";
2088 else if (categoryName == "Muon" || categoryName == "IntermediateMuon" || categoryName == "KinLepton" || categoryName == "IntermediateKinLepton") particleListName = "mu+:inRoe";
2089 else if (categoryName == "Kaon" || categoryName == "KaonPion") {particleListName = "K+:inRoe"; trackTargetName = "Kaon";}
2090 else if (categoryName == "Lambda") particleListName = "Lambda0:inRoe";
2091 else particleListName = "pi+:inRoe";
2092
2093 if (categoryName == "FSC") trackTargetName = "SlowPion";
2094
2095 StoreObjPtr<ParticleList> ListOfParticles(particleListName);
2096 if (!ListOfParticles.isValid()) return Const::doubleNaN;
2097
2098 Variable::Manager& manager = Variable::Manager::Instance();
2099
2100 bool particlesHaveMCAssociated = false;
2101 int nTargets = 0;
2102 for (unsigned int i = 0; i < ListOfParticles->getListSize(); ++i)
2103 {
2104 Particle* iParticle = ListOfParticles->getParticle(i);
2105 if (!iParticle) continue;
2106
2107 if (categoryName == "MaximumPstar") {
2108 static Manager::FunctionPtr selectionFunction =
2109 hasHighestProbInCat({"pi+:inRoe", "isRightTrack(MaximumPstar)"});
2110 bool targetFlag = std::get<bool>(selectionFunction(iParticle));
2111 if (targetFlag) {
2112 particlesHaveMCAssociated = true;
2113 ++nTargets;
2114 }
2115 } else {
2116 int targetFlag = std::get<int>(manager.getVariable("isRightTrack(" + trackTargetName + ")")->function(iParticle));
2117 if (targetFlag != -2) particlesHaveMCAssociated = true;
2118 if (targetFlag == 1) ++nTargets;
2119 }
2120 }
2121
2122 if (!particlesHaveMCAssociated) return Const::doubleNaN;
2123 return (nTargets > 0);
2124 };
2125 return func;
2126 }
2127
2128 Manager::FunctionPtr isTrueCategory(const std::vector<std::string>& arguments)
2129 {
2130 if (arguments.size() != 1) {
2131 B2FATAL("Wrong number of arguments (1 required) for meta function isTrueCategory");
2132 }
2133 auto categoryName = arguments[0];
2134
2135 bool isAvailable = false;
2136 for (const auto& name : availableForIsRightCategory) {
2137 if (categoryName == name) {
2138 isAvailable = true;
2139 break;
2140 }
2141 }
2142 if (categoryName == "mcAssociated")
2143 isAvailable = false;
2144
2145 if (!isAvailable) {
2146 string strAvailableForIsRightCategory;
2147 for (const auto& name : availableForIsRightCategory) {
2148 if (name == "mcAssociated") continue;
2149 strAvailableForIsRightCategory += name + " ";
2150 }
2151 B2FATAL("isTrueCategory: Not available category" << categoryName <<
2152 ". The possibilities for the category name are " << endl << strAvailableForIsRightCategory);
2153 }
2154
2155 auto func = [categoryName](const Particle*) -> double {
2156 std::string particleListName;
2157 std::string trackTargetName = categoryName;
2158
2159 if (categoryName == "Electron" || categoryName == "IntermediateElectron") particleListName = "e+:inRoe";
2160 else if (categoryName == "Muon" || categoryName == "IntermediateMuon" || categoryName == "KinLepton" || categoryName == "IntermediateKinLepton") particleListName = "mu+:inRoe";
2161 else if (categoryName == "Kaon" || categoryName == "KaonPion") {particleListName = "K+:inRoe"; trackTargetName = "Kaon";}
2162 else if (categoryName == "Lambda") particleListName = "Lambda0:inRoe";
2163 else particleListName = "pi+:inRoe";
2164
2165 if (categoryName == "FSC") trackTargetName = "SlowPion";
2166
2167 StoreObjPtr<ParticleList> ListOfParticles(particleListName);
2168 if (!ListOfParticles.isValid()) return Const::doubleNaN;
2169
2170 std::vector<Particle*> targetParticles;
2171 std::vector<Particle*> targetParticlesCategory;
2172 Variable::Manager& manager = Variable::Manager::Instance();
2173
2174
2175 double output = 0;
2176 for (unsigned int i = 0; i < ListOfParticles->getListSize(); ++i)
2177 {
2178 Particle* iParticle = ListOfParticles->getParticle(i);
2179 if (!iParticle) continue;
2180
2181 if (categoryName == "MaximumPstar") {
2182 static Manager::FunctionPtr selectionFunction =
2183 hasHighestProbInCat({"pi+:inRoe", "isRightTrack(MaximumPstar)"});
2184 if (std::get<bool>(selectionFunction(iParticle)))
2185 targetParticles.push_back(iParticle);
2186 } else if (std::get<int>(manager.getVariable("isRightTrack(" + trackTargetName + ")")->function(iParticle))) {
2187 targetParticles.push_back(iParticle);
2188 }
2189 }
2190
2191 for (const auto& targetParticle : targetParticles)
2192 {
2193 int isTargetOfRightCategory = std::get<int>(manager.getVariable("isRightCategory(" + categoryName + ")")->function(
2194 targetParticle));
2195 if (isTargetOfRightCategory == 1) {
2196 output = 1;
2197 targetParticlesCategory.push_back(targetParticle);
2198 } else if (isTargetOfRightCategory == -2 && output != 1)
2199 output = Const::doubleNaN;
2200 }
2201
2202 return output;
2203 };
2204 return func;
2205 }
2206
2207 Manager::FunctionPtr qrOutput(const std::vector<std::string>& arguments)
2208 {
2209 if (arguments.size() != 1)
2210 B2FATAL("Wrong number of arguments for meta function qrOutput");
2211
2212 std::string combinerMethod = arguments[0];
2213 auto func = [combinerMethod](const Particle * particle) -> double {
2214
2215 double output = Const::doubleNaN;
2216 auto* flavorTaggerInfo = particle->getRelatedTo<FlavorTaggerInfo>();
2217
2218 if (flavorTaggerInfo && flavorTaggerInfo->getUseModeFlavorTagger() == "Expert")
2219 output = flavorTaggerInfo->getMethodMap(combinerMethod)->getQrCombined();
2220
2221 return output;
2222 };
2223 return func;
2224 }
2225
2226 Manager::FunctionPtr qOutput(const std::vector<std::string>& arguments)
2227 {
2228 if (arguments.size() != 1)
2229 B2FATAL("Wrong number of arguments for meta function qOutput");
2230
2231 std::string combinerMethod = arguments[0];
2232 auto func = [combinerMethod](const Particle * particle) -> double {
2233
2234 double output = Const::doubleNaN;
2235 auto* flavorTaggerInfo = particle->getRelatedTo<FlavorTaggerInfo>();
2236
2237 if (flavorTaggerInfo && flavorTaggerInfo->getUseModeFlavorTagger() == "Expert")
2238 output = TMath::Sign(1, flavorTaggerInfo->getMethodMap(combinerMethod)->getQrCombined());
2239
2240 return output;
2241 };
2242 return func;
2243 }
2244
2245 Manager::FunctionPtr rBinBelle(const std::vector<std::string>& arguments)
2246 {
2247 if (arguments.size() != 1)
2248 B2FATAL("Wrong number of arguments for meta function rBinBelle");
2249
2250
2251 std::string combinerMethod = arguments[0];
2252 auto func = [combinerMethod](const Particle * particle) -> int {
2253
2254 int output = 0;
2255 auto* flavorTaggerInfo = particle->getRelatedTo<FlavorTaggerInfo>();
2256
2257 if (flavorTaggerInfo && flavorTaggerInfo->getUseModeFlavorTagger() == "Expert")
2258 {
2259 double r = std::abs(flavorTaggerInfo->getMethodMap(combinerMethod)->getQrCombined());
2260 if (r < 0.1) output = 0;
2261 if (r > 0.1 && r < 0.25) output = 1;
2262 if (r > 0.25 && r < 0.5) output = 2;
2263 if (r > 0.5 && r < 0.625) output = 3;
2264 if (r > 0.625 && r < 0.75) output = 4;
2265 if (r > 0.75 && r < 0.875) output = 5;
2266 if (r > 0.875 && r < 1.10) output = 6;
2267 }
2268
2269 return output;
2270 };
2271 return func;
2272 }
2273
2274 Manager::FunctionPtr qpCategory(const std::vector<std::string>& arguments)
2275 {
2276 if (arguments.size() != 1)
2277 B2FATAL("Wrong number of arguments for meta function qpCategory");
2278
2279 std::string categoryName = arguments[0];
2280 auto func = [categoryName](const Particle * particle) -> double {
2281
2282 double output = Const::doubleNaN;
2283 auto* flavorTaggerInfo = particle->getRelatedTo<FlavorTaggerInfo>();
2284
2285 if (flavorTaggerInfo && flavorTaggerInfo->getUseModeFlavorTagger() == "Expert")
2286 {
2287 std::map<std::string, float> iQpCategories = flavorTaggerInfo->getMethodMap("FBDT")->getQpCategory();
2288 if (iQpCategories.find(categoryName) != iQpCategories.end()) output = iQpCategories.at(categoryName);
2289 else if (iQpCategories.size() != 0) B2FATAL("qpCategory: Category with name " << categoryName
2290 << " not found. Check the official category names or if this category is included in the flavor tagger categories list.");
2291 }
2292 return output;
2293 };
2294 return func;
2295 }
2296
2297 Manager::FunctionPtr isTrueFTCategory(const std::vector<std::string>& arguments)
2298 {
2299 if (arguments.size() != 1)
2300 B2FATAL("Wrong number of arguments for meta function isTrueFTCategory");
2301
2302 std::string categoryName = arguments[0];
2303 auto func = [categoryName](const Particle * particle) -> double {
2304
2305 double output = Const::doubleNaN;
2306 auto* flavorTaggerInfo = particle->getRelatedTo<FlavorTaggerInfo>();
2307
2308 if (flavorTaggerInfo && flavorTaggerInfo->getUseModeFlavorTagger() == "Expert")
2309 {
2310 std::map<std::string, float> iIsTrueCategories = flavorTaggerInfo->getMethodMap("FBDT")->getIsTrueCategory();
2311 if (iIsTrueCategories.find(categoryName) != iIsTrueCategories.end()) output = iIsTrueCategories.at(categoryName);
2312 else if (iIsTrueCategories.size() != 0) B2FATAL("isTrueFTCategory: Category with name " << categoryName
2313 << " not found. Check the official category names or if this category is included in the flavor tagger categories list.");
2314 }
2315
2316 return output;
2317 };
2318 return func;
2319 }
2320
2321 Manager::FunctionPtr hasTrueTargets(const std::vector<std::string>& arguments)
2322 {
2323 if (arguments.size() != 1)
2324 B2FATAL("Wrong number of arguments for meta function hasTrueTargets");
2325
2326 std::string categoryName = arguments[0];
2327 auto func = [categoryName](const Particle * particle) -> double {
2328
2329 double output = Const::doubleNaN;
2330 auto* flavorTaggerInfo = particle->getRelatedTo<FlavorTaggerInfo>();
2331
2332 if (flavorTaggerInfo && flavorTaggerInfo->getUseModeFlavorTagger() == "Expert")
2333 {
2334 std::map<std::string, float> iHasTrueTargets = flavorTaggerInfo->getMethodMap("FBDT")->getHasTrueTarget();
2335 if (iHasTrueTargets.find(categoryName) != iHasTrueTargets.end()) output = iHasTrueTargets.at(categoryName);
2336 else if (iHasTrueTargets.size() != 0) B2FATAL("hasTrueTargets: Category with name " << categoryName
2337 << " not found. Check the official category names or if this category is included in the flavor tagger categories list.");
2338 }
2339
2340 return output;
2341 };
2342 return func;
2343 }
2344
2345 VARIABLE_GROUP("Flavor Tagger Expert Variables");
2346
2347 REGISTER_VARIABLE("pMissTag", momentumMissingTagSide, R"DOC(
2348[Expert] Calculates the missing momentum for a given particle on the tag side.
2349:noindex:
2350)DOC","GeV/c");
2351 REGISTER_METAVARIABLE("pMissTag(maskName)", momentumMissingTagSideWithMask,
2352 "[Expert] Calculates the missing momentum for a given particle on the tag side. The unit of the missing momentum is ``GeV/c`` ",
2353 Manager::VariableDataType::c_double);
2354 REGISTER_VARIABLE("cosTPTO" , cosTPTO, R"DOC(
2355[Expert] Returns cosine of angle between thrust axis of given particle and thrust axis of ROE.
2356:noindex:
2357)DOC");
2358 REGISTER_METAVARIABLE("cosTPTO(maskName)", cosTPTOWithMask,
2359 "[Expert] Returns cosine of angle between thrust axis of given particle and thrust axis of ROE.",
2360 Manager::VariableDataType::c_double);
2361 REGISTER_VARIABLE("lambdaFlavor", lambdaFlavor,
2362 "[Expert] Returns 1.0 if particle is ``Lambda0``, -1.0 in case of ``anti-Lambda0``, 0.0 otherwise.");
2363 REGISTER_VARIABLE("isLambda", isLambda, "[Expert] Returns 1.0 if particle is truth-matched to ``Lambda0``, 0.0 otherwise.");
2364 REGISTER_VARIABLE("lambdaZError", lambdaZError, "[Expert] Returns the variance of the z-component of the decay vertex.\n\n",":math:`\\text{cm}^2`");
2365 REGISTER_VARIABLE("momentumOfSecondDaughter", momentumOfSecondDaughter,
2366 "[Expert] Returns the momentum of second daughter if exists, 0. otherwise.\n\n","GeV/c");
2367 REGISTER_VARIABLE("momentumOfSecondDaughterCMS", momentumOfSecondDaughterCMS,
2368 "[Expert] Returns the momentum of the second daughter in the centre-of-mass system, 0. if this daughter doesn't exist.\n\n","GeV/c");
2369 REGISTER_VARIABLE("chargeTimesKaonLiklihood", chargeTimesKaonLiklihood,
2370 "[Expert] Returns ``q*(highest PID_Likelihood for Kaons)``, 0. otherwise.");
2371 REGISTER_VARIABLE("ptTracksRoe", transverseMomentumOfChargeTracksInRoe, R"DOC(
2372[Expert] Returns the transverse momentum of all charged tracks of the ROE related to the given particle, 0.0 if particle has no related ROE.
2373:noindex:
2374)DOC","GeV/c");
2375 REGISTER_METAVARIABLE("ptTracksRoe(maskName)", transverseMomentumOfChargeTracksInRoeWithMask,
2376 "[Expert] Returns the transverse momentum of all charged tracks of the ROE related to the given particle, 0.0 if particle has no related ROE. The unit of the momentum is ``GeV/c`` ",
2377 Manager::VariableDataType::c_double);
2378 REGISTER_METAVARIABLE("pt2TracksRoe(maskName)", transverseMomentumSquaredOfChargeTracksInRoeWithMask,
2379 "[Expert] Returns the transverse momentum squared of all charged tracks of the ROE related to the given particle, 0.0 if particle has no related ROE. The unit of the momentum squared is :math:`[\\text{GeV}/\\text{c}]^2` ",
2380 Manager::VariableDataType::c_double);
2381 REGISTER_VARIABLE("NumberOfKShortsInRoe", NumberOfKShortsInRoe,
2382 "[Expert] Returns the number of ``K_S0`` in the rest of event. The particle list ``K_S0:inRoe`` has to be filled beforehand.");
2383
2384 REGISTER_VARIABLE("isInElectronOrMuonCat", isInElectronOrMuonCat,
2385 "[Expert] Returns 1.0 if the particle has been selected as target in the Muon or Electron Category, 0.0 otherwise.");
2386
2387 REGISTER_VARIABLE("isMajorityInRestOfEventFromB0", isMajorityInRestOfEventFromB0, R"DOC(
2388[Eventbased][Expert] Checks if the majority of the tracks in the current RestOfEvent are from a ``B0``.
2389:noindex:
2390)DOC");
2391 REGISTER_METAVARIABLE("isMajorityInRestOfEventFromB0(maskName)", isMajorityInRestOfEventFromB0WithMask,
2392 "[Eventbased][Expert] Checks if the majority of the tracks in the current RestOfEvent are from a ``B0``.",
2393 Manager::VariableDataType::c_bool);
2394 REGISTER_VARIABLE("isMajorityInRestOfEventFromB0bar", isMajorityInRestOfEventFromB0bar, R"DOC(
2395[Eventbased][Expert] Check if the majority of the tracks in the current RestOfEvent are from a ``anti-B0``.
2396:noindex:
2397)DOC");
2398 REGISTER_METAVARIABLE("isMajorityInRestOfEventFromB0bar(maskName)", isMajorityInRestOfEventFromB0barWithMask,
2399 "[Eventbased][Expert] Check if the majority of the tracks in the current RestOfEvent are from a ``anti-B0``.",
2400 Manager::VariableDataType::c_bool);
2401 REGISTER_VARIABLE("hasRestOfEventTracks", hasRestOfEventTracks, R"DOC(
2402[Expert] Returns the number of tracks in the RestOfEvent related to the given Particle. -2 if the RestOfEvent is empty.
2403:noindex:
2404)DOC");
2405 REGISTER_METAVARIABLE("hasRestOfEventTracks(maskName)", hasRestOfEventTracksWithMask,
2406 "[Expert] Returns the number of tracks in the RestOfEvent related to the given Particle. -2 if the RestOfEvent is empty.",
2407 Manager::VariableDataType::c_bool);
2408
2409 REGISTER_VARIABLE("qrCombined", isRestOfEventB0Flavor, R"DOC(
2410[Eventbased][Expert] Returns -1 (1) if current RestOfEvent is related to a ``anti-B0`` (``B0``).
2411The ``MCError`` bit of Breco has to be 0, 1, 2, 16 or 1024.
2412The output of the variable is 0 otherwise.
2413If one particle in the RestOfEvent is found to belong to the reconstructed ``B0``, the output is -2(2) for a ``anti-B0`` (``B0``) on the reconstructed side.
2414)DOC");
2415 REGISTER_VARIABLE("ancestorHasWhichFlavor", ancestorHasWhichFlavor,
2416 "[Expert] Checks the decay chain of the given particle upwards up to the ``Upsilon(4S)`` resonance and outputs 0 (1) if an ancestor is found to be a ``anti-B0`` (``B0``), if not -2.");
2417 REGISTER_VARIABLE("B0mcErrors", B0mcErrors, "[Expert] Returns MC-matching flag, see :b2:var:`mcErrors` for the particle, e.g. ``B0`` .");
2418 REGISTER_VARIABLE("isRelatedRestOfEventMajorityB0Flavor", isRelatedRestOfEventMajorityB0Flavor, R"DOC(
2419[Expert] Returns 0 (1) if the majority of tracks and clusters of the RestOfEvent related to the given Particle are related to a ``anti-B0`` (``B0``).
2420:noindex:
2421)DOC");
2422 REGISTER_METAVARIABLE("isRelatedRestOfEventMajorityB0Flavor(maskName)", isRelatedRestOfEventMajorityB0FlavorWithMask,
2423 "[Expert] Returns 0 (1) if the majority of tracks and clusters of the RestOfEvent related to the given Particle are related to a ``anti-B0`` (``B0``).",
2424 Manager::VariableDataType::c_int);
2425 REGISTER_VARIABLE("isRestOfEventMajorityB0Flavor", isRestOfEventMajorityB0Flavor,
2426 "[Expert] Returns 0 (1) if the majority of tracks and clusters of the current RestOfEvent are related to a ``anti-B0`` (``B0``).");
2427 REGISTER_VARIABLE("mcFlavorOfOtherB", mcFlavorOfOtherB, R"DOC(
2428[Expert] Returns the MC flavor (+1 or -1) of the accompanying tag-side B meson if the given particle is a correctly truth-matched B candidate, 0 otherwise.
2429In other words, this variable checks the generated flavor of the other generated ``Upsilon(4S)`` daughter.
2430)DOC");
2431
2432
2433 REGISTER_METAVARIABLE("BtagToWBosonVariables(requestedVariable[, maskName])", BtagToWBosonVariables, R"DOC(
2434[Eventbased][Expert] Returns values of FlavorTagging-specific kinematical variables assuming a semileptonic decay with the given particle as target.
2435The input values of ``requestedVariable`` can be the following: recoilMass in GeV/c^2 , pMissCMS in ``GeV/c``, cosThetaMissCMS and EW90.
2436)DOC", Manager::VariableDataType::c_double);
2437 REGISTER_METAVARIABLE("KaonPionVariables(requestedVariable)" , KaonPionVariables , R"DOC(
2438[Expert] Returns values of FlavorTagging-specific kinematical variables for ``KaonPion`` category.
2439The input values of ``requestedVariable`` can be the following: cosKaonPion, HaveOpositeCharges.
2440)DOC", Manager::VariableDataType::c_double);
2441 REGISTER_METAVARIABLE("FSCVariables(requestedVariable)", FSCVariables, R"DOC(
2442[Eventbased][Expert] Returns values of FlavorTagging-specific kinematical variables for ``FastSlowCorrelated`` category.
2443The input values of ``requestedVariable`` can be the following: pFastCMS in ``GeV/c``, cosSlowFast, SlowFastHaveOpositeCharges, or cosTPTOFast.
2444)DOC", Manager::VariableDataType::c_double);
2445 REGISTER_METAVARIABLE("hasHighestProbInCat(particleListName, extraInfoName)", hasHighestProbInCat, R"DOC(
2446[Expert] Returns 1.0 if the given Particle is classified as target track, i.e. if it has the highest target track probability in particleListName.
2447The probability is accessed via ``extraInfoName``, which can have the following input values:
2448
2449* isRightTrack(Electron),
2450* isRightTrack(IntermediateElectron),
2451* isRightTrack(Muon),
2452* isRightTrack(IntermediateMuon),
2453* isRightTrack(KinLepton),
2454* isRightTrack(IntermediateKinLepton),
2455* isRightTrack(Kaon),
2456* isRightTrack(SlowPion),
2457* isRightTrack(FastHadron),
2458* isRightTrack(MaximumPstar),
2459* isRightTrack(Lambda),
2460* isRightCategory(Electron),
2461* isRightCategory(IntermediateElectron),
2462* isRightCategory(Muon),
2463* isRightCategory(IntermediateMuon),
2464* isRightCategory(KinLepton),
2465* isRightCategory(IntermediateKinLepton),
2466* isRightCategory(Kaon),
2467* isRightCategory(SlowPion),
2468* isRightCategory(FastHadron),
2469* isRightCategory(MaximumPstar),
2470* isRightCategory(Lambda),
2471* isRightCategory(KaonPion),
2472* isRightCategory(FSC).
2473
2474)DOC", Manager::VariableDataType::c_bool);
2475 REGISTER_METAVARIABLE("HighestProbInCat(particleListName, extraInfoName)", HighestProbInCat,
2476 "[Expert] Returns the highest target track probability value for the given category, for allowed input values for ``extraInfoName`` see :b2:var:`hasHighestProbInCat`.", Manager::VariableDataType::c_double);
2477
2478 REGISTER_METAVARIABLE("isRightTrack(particleName)", isRightTrack, R"DOC(
2479[Expert] Returns 1.0 if the given particle was really from a B-meson depending on category provided in ``particleName`` argument, 0.0 otherwise.
2480Allowed input values for ``particleName`` argument in this variable are the following:
2481
2482* Electron,
2483* IntermediateElectron,
2484* Muon,
2485* IntermediateMuon,
2486* KinLepton,
2487* IntermediateKinLepton,
2488* Kaon,
2489* SlowPion,
2490* FastHadron,
2491* Lambda,
2492* mcAssociated.
2493
2494)DOC", Manager::VariableDataType::c_int);
2495 REGISTER_METAVARIABLE("isRightCategory(particleName)", isRightCategory, R"DOC(
2496[Expert] Returns 1.0 if the class track by ``particleName`` category has the same flavor as the MC target track, 0.0 otherwise.
2497Allowed input values for ``particleName`` argument in this variable are the following:
2498
2499* Electron,
2500* IntermediateElectron,
2501* Muon,
2502* IntermediateMuon,
2503* KinLepton,
2504* IntermediateKinLepton
2505* Kaon,
2506* SlowPion,
2507* FastHadron,
2508* KaonPion,
2509* MaximumPstar,
2510* FSC,
2511* Lambda,
2512* mcAssociated.
2513
2514)DOC", Manager::VariableDataType::c_int);
2515 REGISTER_METAVARIABLE("QpOf(particleListName, outputExtraInfo, rankingExtraInfo)", QpOf, R"DOC(
2516[Eventbased][Expert] Returns the :math:`q*p` value for a given particle list provided as the 1st argument,
2517where :math:`p` is the probability of a category stored as extraInfo, provided as the 2nd argument,
2518allowed values are same as in :b2:var:`hasHighestProbInCat`.
2519The particle is selected after ranking according to a flavor tagging extraInfo, provided as the 3rd argument,
2520allowed values are same as in :b2:var:`hasHighestProbInCat`.
2521)DOC", Manager::VariableDataType::c_double);
2522 REGISTER_METAVARIABLE("weightedQpOf(particleListName, outputExtraInfo, rankingExtraInfo)", weightedQpOf, R"DOC(
2523[Eventbased][Expert] Returns the weighted :math:`q*p` value for a given particle list, provided as the 1st argument,
2524where :math:`p` is the probability of a category stored as extraInfo, provided in the 2nd argument,
2525allowed values are same as in :b2:var:`hasHighestProbInCat`.
2526The particles in the list are ranked according to a flavor tagging extraInfo, provided as the 3rd argument,
2527allowed values are same as in :b2:var:`hasHighestProbInCat`.
2528The values for the three top particles is combined into an effective (weighted) output.
2529)DOC", Manager::VariableDataType::c_double);
2530 REGISTER_METAVARIABLE("QpTrack(particleListName, outputExtraInfo, rankingExtraInfo)", QpTrack, R"DOC(
2531[Expert] Returns the :math:`q*p` value of the particle in a given particle list provided as the 1st argument that is originated from the same Track of given particle.
2532where :math:`p` is the probability of a category stored as extraInfo, provided as the 2nd argument,
2533allowed values are same as in :b2:var:`hasHighestProbInCat`.
2534The particle is selected after ranking according to a flavor tagging extraInfo, provided as the 3rd argument,
2535allowed values are same as in :b2:var:`hasHighestProbInCat`.
2536)DOC", Manager::VariableDataType::c_double);
2537
2538 REGISTER_METAVARIABLE("variableOfTarget(particleListName, inputVariable, rankingExtraInfo)", variableOfTarget, R"DOC(
2539[Eventbased][Expert] Returns the value of an input variable provided as the 2nd argument for a particle selected from the given list provided as the 1st argument.
2540The particles are ranked according to a flavor tagging extraInfo, provided as the 2nd argument,
2541allowed values are same as in :b2:var:`hasHighestProbInCat`.
2542)DOC", Manager::VariableDataType::c_double);
2543
2544 REGISTER_METAVARIABLE("hasTrueTarget(categoryName)", hasTrueTarget,
2545 "[Expert] Returns 1 if the given category has a target, 0 otherwise.", Manager::VariableDataType::c_double);
2546 REGISTER_METAVARIABLE("isTrueCategory(categoryName)", isTrueCategory,
2547 "[Expert] Returns 1 if the given category tags the B0 MC flavor correctly, 0 otherwise.", Manager::VariableDataType::c_double);
2548
2549 REGISTER_METAVARIABLE("qpCategory(categoryName)", qpCategory, R"DOC(
2550[Expert] Returns the output :math:`q` (charge of target track) times :math:`p` (probability that this is the right category) of the category with the given name.
2551The allowed categories are the official Flavor Tagger Category Names.
2552)DOC", Manager::VariableDataType::c_double);
2553 REGISTER_METAVARIABLE("isTrueFTCategory(categoryName)", isTrueFTCategory, R"DOC(
2554[Expert] Returns 1 if the target particle (checking the decay chain) of the category with the given name is found in the MC particles,
2555and if it provides the right flavor. The allowed categories are the official Flavor Tagger Category Names.
2556)DOC", Manager::VariableDataType::c_double);
2557 REGISTER_METAVARIABLE("hasTrueTargets(categoryName)", hasTrueTargets, R"DOC(
2558[Expert] Returns 1 if target particles (checking only the decay chain) of the category with the given name is found in the MC particles.
2559The allowed categories are the official Flavor Tagger Category Names.
2560)DOC", Manager::VariableDataType::c_double);
2561
2562 VARIABLE_GROUP("Flavor Tagger Analysis Variables")
2563
2564 REGISTER_METAVARIABLE("rBinBelle(combinerMethod)", rBinBelle, R"DOC(
2565Returns the corresponding :math:`r` (dilution) bin according to the Belle binning for the given ``combinerMethod``.
2566The available methods are 'FBDT' and 'FANN' (category-based combiners), and 'DNN' (DNN tagger output).
2567The return values and the corresponding dilution ranges are the following:
2568
2569* 0: :math:`0.000 < r < 0.100`;
2570* 1: :math:`0.100 < r < 0.250`;
2571* 2: :math:`0.250 < r < 0.500`;
2572* 3: :math:`0.500 < r < 0.625`;
2573* 4: :math:`0.625 < r < 0.750`;
2574* 5: :math:`0.750 < r < 0.875`;
2575* 6: :math:`0.875 < r < 1.000`.
2576
2577.. warning:: You have to run the Flavor Tagger for this variable to be meaningful.
2578.. seealso:: :ref:`FlavorTagger` and :func:`flavorTagger.flavorTagger`.
2579)DOC", Manager::VariableDataType::c_int);
2580 REGISTER_METAVARIABLE("qrOutput(combinerMethod)", qrOutput, R"DOC(
2581Returns the output of the flavorTagger, flavor tag :math:`q` times the dilution factor :math:`r`, for the given combiner method.
2582The available methods are 'FBDT' and 'FANN' (category-based combiners), and 'DNN' (DNN tagger output).
2583
2584.. warning:: You have to run the Flavor Tagger for this variable to be meaningful.
2585.. seealso:: :ref:`FlavorTagger` and :func:`flavorTagger.flavorTagger`.
2586)DOC", Manager::VariableDataType::c_double);
2587 REGISTER_METAVARIABLE("qOutput(combinerMethod)", qOutput, R"DOC(
2588Returns the flavor tag :math:`q` output of the flavorTagger for the given combinerMethod.
2589The available methods are 'FBDT' and 'FANN' (category-based combiners), and 'DNN' (DNN tagger output).
2590
2591.. warning:: You have to run the Flavor Tagger for this variable to be meaningful.
2592.. seealso:: :ref:`FlavorTagger` and :func:`flavorTagger.flavorTagger`.
2593)DOC", Manager::VariableDataType::c_double);
2594 REGISTER_VARIABLE("isRelatedRestOfEventB0Flavor", isRelatedRestOfEventB0Flavor, R"DOC(
2595Returns -1 (1) if the RestOfEvent related to the given particle is related to a ``anti-B0`` (``B0``).
2596The ``MCError`` bit of Breco has to be 0, 1, 2, 16 or 1024.
2597The output of the variable is 0 otherwise.
2598If one particle in the RestOfEvent is found to belong to the reconstructed ``B0``, the output is -2(2) for a ``anti-B0`` (``B0``) on the reconstructed side.
2599)DOC");
2600 }
2602}
double R
typedef autogenerated by FFTW
int getPDGCode() const
PDG code.
Definition Const.h:473
static const ParticleType Lambda
Lambda particle.
Definition Const.h:679
static const ChargedStable muon
muon particle
Definition Const.h:660
static const ChargedStable pion
charged pion particle
Definition Const.h:661
static const ParticleType antiLambda
Anti-Lambda particle.
Definition Const.h:680
static const double doubleNaN
quiet_NaN
Definition Const.h:703
static const ChargedStable kaon
charged kaon particle
Definition Const.h:662
static const ChargedStable electron
electron particle
Definition Const.h:659
@ c_nPhotons
CR is split into n photons (N1)
Definition ECLCluster.h:41
static ROOT::Math::PxPyPzMVector labToCms(const ROOT::Math::PxPyPzMVector &vec)
Transforms Lorentz vector into CM System.
static constexpr const char * c_defaultMaskName
Default mask name.
Definition RestOfEvent.h:58
static ROOT::Math::XYZVector calculateThrust(const std::vector< ROOT::Math::PxPyPzEVector > &momenta)
calculates the thrust axis
Definition Thrust.cc:71
Global list of available variables.
Definition Manager.h:100
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
std::string particleName(int pdgCode)
Returns the name of a particle with given pdg code.
Definition EvtPDLUtil.cc:20
Abstract base class for different kinds of events.
STL namespace.
static int getMCErrors(const Particle *particle, const MCParticle *mcParticle=nullptr)
Returns quality indicator of the match as a bit pattern where the individual bits indicate the the ty...