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