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