Belle II Software development
FEIVariables.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9// Own header.
10#include <analysis/variables/FEIVariables.h>
11
12// include VariableManager
13#include <analysis/VariableManager/Manager.h>
14
15#include <analysis/dataobjects/Particle.h>
16#include <mdst/dataobjects/MCParticle.h>
17#include <framework/datastore/StoreArray.h>
18#include <framework/logging/Logger.h>
19#include <map>
20#include <cmath>
21#include <algorithm>
22#include <iterator>
23#include <limits>
24
25namespace Belle2 {
30 namespace Variable {
31
32 int mostcommonBTagIndex(const Particle* part)
33 {
34 std::map <int, int> tag_candidates;
35 const std::vector<const Particle*>& fsp_tag = part->getFinalStateDaughters();
36 for (const Particle* fsp : fsp_tag) {
37 const MCParticle* mc_fsp = fsp->getMCParticle();
38 if (mc_fsp) {
39 int tag_index = finddescendant(mc_fsp);
40 if (tag_candidates.find(tag_index) == tag_candidates.end()) {
41 tag_candidates[tag_index] = 1;
42 } else {
43 tag_candidates[tag_index]++;
44 }
45 }
46 }
47 int tag_index = -1;
48 int tag_N = 0;
49 // cppcheck-suppress unassignedVariable
50 for (const auto& [key, value] : tag_candidates) {
51 if (value > tag_N) {
52 tag_index = key;
53 tag_N = value;
54 }
55 }
56 return tag_index;
57 }
58 int finddescendant(const MCParticle* mcpart)
59 {
60 const MCParticle* i_mcpart = mcpart;
61 while (i_mcpart) {
62 auto* i_mcpart_mother = i_mcpart->getMother();
63 if (i_mcpart_mother) {
64 std::vector<int> B_PDG = {511, 521};
65 auto result = std::find(std::begin(B_PDG), std::end(B_PDG), abs(i_mcpart_mother->getPDG()));
66 if (result != std::end(B_PDG)) {
67 return i_mcpart_mother->getArrayIndex();
68 }
69 i_mcpart = i_mcpart_mother;
70 } else {
71 return -1;
72 }
73 }
74 return -2;
75 }
76 std::vector<int> truthFSPTag(int BTag_index)
77 {
78 StoreArray<MCParticle> MC_Particle_list;
79 std::vector<int> fsp_truth_index;
80 for (const MCParticle& iMCParticle : MC_Particle_list) {
81 if ((BTag_index == finddescendant(&iMCParticle)) && (iMCParticle.hasStatus(MCParticle::c_StableInGenerator) == true)
82 && (iMCParticle.hasStatus(MCParticle::c_IsISRPhoton) == false) && iMCParticle.hasStatus(MCParticle::c_IsFSRPhoton) == false
83 && (iMCParticle.hasStatus(MCParticle::c_IsPHOTOSPhoton) == false)) {
84 fsp_truth_index.push_back(iMCParticle.getArrayIndex());
85 }
86 }
87 return fsp_truth_index;
88 }
89 double percentageMissingParticlesBTag(const Particle* part)
90 {
91 int index = mostcommonBTagIndex(part);
92 if (index == -1) {
93 return -1;
94 } else {
95 std::vector<int> fsp_FEI_index;
96 const std::vector<const Particle*>& fsp_tag = part->getFinalStateDaughters();
97 for (const Particle* fsp : fsp_tag) {
98 const MCParticle* mc_fsp = fsp->getMCParticle();
99 if (mc_fsp) {
100 int tag_index = finddescendant(mc_fsp);
101 if (tag_index == index) {
102 fsp_FEI_index.push_back(mc_fsp->getArrayIndex());
103 }
104 }
105 }
106 std::vector<int> diff;
107 std::vector<int> fsp_truth_index = truthFSPTag(index);
108 std::sort(fsp_truth_index.begin(), fsp_truth_index.end());
109 std::sort(fsp_FEI_index.begin(), fsp_FEI_index.end());
110 std::set_difference(fsp_truth_index.begin(), fsp_truth_index.end(), fsp_FEI_index.begin(), fsp_FEI_index.end(), std::inserter(diff,
111 diff.begin()));
112 return double(diff.size()) / double(fsp_truth_index.size());
113 }
114 }
115 double percentageWrongParticlesBTag(const Particle* part)
116 {
117 int index = mostcommonBTagIndex(part);
118 if (index == -1) {
119 return -1;
120 } else {
121 const std::vector<const Particle*>& fsp_tag = part->getFinalStateDaughters();
122 int wrong_FEI = 0;
123 for (const Particle* fsp : fsp_tag) {
124 const MCParticle* mc_fsp = fsp->getMCParticle();
125 if (mc_fsp) {
126 int tag_index = finddescendant(mc_fsp);
127 if (tag_index != index) {
128 wrong_FEI ++;
129 }
130 }
131 }
132 return double(wrong_FEI) / double(truthFSPTag(index).size());
133 }
134 }
135 VARIABLE_GROUP("FEIVariables");
136 REGISTER_VARIABLE("mostcommonBTagIndex", mostcommonBTagIndex,
137 "By giving e.g. a FEI B meson candidate the B meson index on generator level is determined, where most reconstructed particles can be assigned to. If no B meson found on generator level -1 is returned.");
138 REGISTER_VARIABLE("percentageMissingParticlesBTag", percentageMissingParticlesBTag,
139 "Get the percentage of missing particles by using the mostcommonBTagIndex. So the number of particles not reconstructed by e.g. the FEI are determined and divided by the number of generated particles using the given index of the B meson. If no B meson found on generator level -1 is returned.");
140 REGISTER_VARIABLE("percentageWrongParticlesBTag", percentageWrongParticlesBTag,
141 "Get the percentage of wrong particles by using the mostcommonBTagIndex. In this context wrong means that the reconstructed particles originated from the other B meson. The absolute number is divided by the total number of generated FSP from the given B meson index. If no B meson found on generator level -1 is returned.");
142 }
144}
@ c_IsFSRPhoton
bit 7: Particle is from finial state radiation
Definition: MCParticle.h:61
@ c_IsPHOTOSPhoton
bit 8: Particle is an radiative photon from PHOTOS
Definition: MCParticle.h:63
@ c_StableInGenerator
bit 1: Particle is stable, i.e., not decaying in the generator.
Definition: MCParticle.h:49
@ c_IsISRPhoton
bit 6: Particle is from initial state radiation
Definition: MCParticle.h:59
Abstract base class for different kinds of events.