Belle II Software  release-05-01-25
PIDPriors.cc
1 /****************************************************************
2  * This test file for the PID prior classes *
3  ****************************************************************/
4 
5 #include <analysis/dbobjects/PIDPriorsTable.h>
6 #include <analysis/dbobjects/PIDPriors.h>
7 
8 #include <framework/gearbox/Const.h>
9 #include <framework/utilities/TestHelpers.h>
10 #include <TRandom3.h>
11 #include <TH2F.h>
12 
13 #include <gtest/gtest.h>
14 
15 using namespace std;
16 
17 namespace Belle2 {
23  class PIDPriorsTest : public ::testing::Test {
24  protected:
25  };
26 
28  TEST_F(PIDPriorsTest, PIDPriorsTableTest)
29  {
30 
31  std::vector<float> edgesX = {0., 0.2, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.7, 1.9, 2.2, 2.5, 3., 3.5, 4., 5., 6., 7., 10.};
32  std::vector<float> edgesY = { -1., -0.95, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 1.};
33 
34  PIDPriorsTable PIDTable = PIDPriorsTable();
35  PIDTable.setBinEdges(edgesX, edgesY);
36 
37  PIDTable.printPrior();
38 
39  PIDTable.setPriorValue(1.1, 0.54, 0.1);
40  PIDTable.setPriorValue(9.5, -0.67, 0.3);
41 
42  std::cout << "after setting (1.1, 0.54) amd (9.5, -0.67)" << std::endl;
43  PIDTable.printPrior();
44 
45  EXPECT_B2WARNING(PIDTable.setPriorValue(11.0, -10.01, 3.));
46  EXPECT_B2WARNING(PIDTable.setPriorValue(4.23, -0.01, 4.));
47 
48  std::cout << " ----- Setting the priors from outside ------ " << std::endl;
49 
50  std::vector<float> prior;
51 
52  for (int j = 0; j < (int)edgesY.size() - 1; j++) {
53  for (int i = 0; i < (int)edgesX.size() - 1; i++) {
54  prior.push_back(0.001 * (i + j * (edgesX.size() - 1)));
55  }
56  }
57 
58 
59  PIDTable.setPriorsTable(prior);
60 
61  PIDTable.printPrior();
62 
63  EXPECT_FLOAT_EQ(PIDTable.getPriorValue(0.1, -0.98), 0.001 * (0 + 0 * (edgesX.size() - 1)));
64  EXPECT_FLOAT_EQ(PIDTable.getPriorValue(0.71, -0.65), 0.001 * (5 + 4 * (edgesX.size() - 1)));
65  EXPECT_FLOAT_EQ(PIDTable.getPriorValue(0.71, -1.65), 0.);
66  EXPECT_FLOAT_EQ(PIDTable.getPriorValue(11.71, 0.876), 0.);
67  EXPECT_B2WARNING(PIDTable.getPriorValue(0.71, -1.65));
68  EXPECT_B2WARNING(PIDTable.getPriorValue(11.71, 0.876));
69  EXPECT_FLOAT_EQ(0., PIDTable.getPriorValue(11.71, 0.876));
70 
71 
72 
73  // now let's try to construct a table with no or on edge only
74  std::vector<float> badEdges1 = {0.};
75  std::vector<float> badEdges2 = {};
76 
77 
78  PIDPriorsTable PIDTable2 = PIDPriorsTable();
79 
80  EXPECT_B2WARNING(PIDTable2.setBinEdges(edgesX, badEdges1));
81  EXPECT_B2WARNING(PIDTable2.setBinEdges(edgesX, badEdges2));
82 
83 
84  PIDTable2.setBinEdges(edgesX, badEdges2);
85 
86  PIDTable2.printPrior();
87 
88  PIDTable2.setPriorValue(1.1, 0.54, 0.3);
89  PIDTable2.setPriorValue(9.5, -0.67, 0.3);
90 
91  std::cout << "after setting (1.1, 0.54) and (9.5, -0.67)" << std::endl;
92  PIDTable2.printPrior();
93 
94  EXPECT_B2WARNING(PIDTable2.setPriorValue(11.0, -10.01, 0.3));
95  EXPECT_FLOAT_EQ(PIDTable2.getPriorValue(9.516, 0.), 0.3);
96 
97  }
98 
99 
101  TEST_F(PIDPriorsTest, PIDPriorTest)
102  {
103 
104  // Creates 6 prior objects in different formats (array, PIDPriorsTable, TH2F) to test different setters
105 
106  std::vector<float> edgesX = {0., 1, 2, 3, 4, 5, 6, 7, 8, 9, 10.};
107  std::vector<float> edgesY = { -1., -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.};
108 
109  std::vector<float> prior1;
110  std::vector<float> prior2;
111  std::vector<float> prior3;
112  std::vector<float> prior4;
113  std::vector<float> prior5;
114  std::vector<float> prior6;
115 
116  TH2F* counts1 = new TH2F("counts1", "counts1", 10, 0, 10., 20, -1, 1);
117  TH2F* counts2 = new TH2F("counts2", "counts2", 10, 0, 10., 20, -1, 1);
118  TH2F* counts3 = new TH2F("counts3", "counts3", 10, 0, 10., 20, -1, 1);
119  TH2F* counts4 = new TH2F("counts4", "counts4", 10, 0, 10., 20, -1, 1);
120  TH2F* counts5 = new TH2F("counts5", "counts5", 10, 0, 10., 20, -1, 1);
121  TH2F* counts6 = new TH2F("counts6", "counts6", 10, 0, 10., 20, -1, 1);
122 
123  TH2F* norms = new TH2F("norms", "norms", 10, 0, 10., 20, -1, 1);
124 
125  TRandom3* rn = new TRandom3();
126 
127 
128 
129  for (int i = 0; i < (int)edgesX.size() - 1; i++) {
130  for (int j = 0; j < (int)edgesY.size() - 1; j++) {
131 
132  float n1 = (float)rn->Poisson(100);
133  float n2 = (float)rn->Poisson(100);
134  float n3 = (float)rn->Poisson(700);
135  float n4 = (float)rn->Poisson(200);
136  float n5 = (float)rn->Poisson(100);
137  float n6 = (float)rn->Poisson(10);
138 
139  float n = n1 + n2 + n3 + n4 + n5 + n6;
140 
141  prior1.push_back(n1 / n);
142  prior2.push_back(n2 / n);
143  prior3.push_back(n3 / n);
144  prior4.push_back(n4 / n);
145  prior5.push_back(n5 / n);
146  prior6.push_back(n6 / n);
147 
148  counts1->SetBinContent(i + 1, j + 1, n1);
149  counts2->SetBinContent(i + 1, j + 1, n2);
150  counts3->SetBinContent(i + 1, j + 1, n3);
151  counts4->SetBinContent(i + 1, j + 1, n4);
152  counts5->SetBinContent(i + 1, j + 1, n5);
153  counts6->SetBinContent(i + 1, j + 1, n6);
154  norms->SetBinContent(i + 1, j + 1, n);
155 
156  }
157  }
158 
159 
160  PIDPriorsTable PIDTable1 = PIDPriorsTable();
161  PIDTable1.setBinEdges(edgesX, edgesY);
162  PIDTable1.setPriorsTable(prior1);
163 
164  PIDPriorsTable PIDTable2 = PIDPriorsTable();
165  PIDTable2.setBinEdges(edgesX, edgesY);
166  PIDTable2.setPriorsTable(prior2);
167 
168  PIDPriorsTable PIDTable3 = PIDPriorsTable();
169  PIDTable3.setBinEdges(edgesX, edgesY);
170  PIDTable3.setPriorsTable(prior3);
171 
172  PIDPriorsTable PIDTable4 = PIDPriorsTable();
173  PIDTable4.setBinEdges(edgesX, edgesY);
174  PIDTable4.setPriorsTable(prior4);
175 
176  PIDPriorsTable PIDTable5 = PIDPriorsTable();
177  PIDTable5.setBinEdges(edgesX, edgesY);
178  PIDTable5.setPriorsTable(prior5);
179 
180  PIDPriorsTable PIDTable6 = PIDPriorsTable();
181  PIDTable6.setBinEdges(edgesX, edgesY);
182  PIDTable6.setPriorsTable(prior6);
183 
184  Const::ChargedStable part_e = Const::ChargedStable(11);
185  Const::ChargedStable part_mu = Const::ChargedStable(13);
186  Const::ChargedStable part_pi = Const::ChargedStable(211);
187  Const::ChargedStable part_K = Const::ChargedStable(321);
188  Const::ChargedStable part_p = Const::ChargedStable(2212);
189  Const::ChargedStable part_d = Const::ChargedStable(1000010020);
190 
191 
192  std::cout << "setting the tables " << std::endl;
193  // setter with the priortable
194  PIDPriors priors = PIDPriors();
195  priors.setPriors(part_e, PIDTable1);
196  priors.setPriors(part_mu, PIDTable2);
197  priors.setPriors(part_pi, PIDTable3);
198  priors.setPriors(part_K, PIDTable4);
199  priors.setPriors(part_p, PIDTable5);
200  priors.setPriors(part_d, PIDTable6);
201 
202  // setter with the count histograms
203  PIDPriors priors2 = PIDPriors();
204  priors2.setPriors(part_e, counts1, norms);
205  priors2.setPriors(part_mu, counts2, norms);
206  priors2.setPriors(part_pi, counts3, norms);
207  priors2.setPriors(part_K, counts4, norms);
208  priors2.setPriors(part_p, counts5, norms);
209  priors2.setPriors(part_d, counts6, norms);
210 
211  counts1->Divide(norms);
212  counts2->Divide(norms);
213  counts3->Divide(norms);
214  counts4->Divide(norms);
215  counts5->Divide(norms);
216  counts6->Divide(norms);
217 
218  std::cout << "division done" << std::endl;
219  // setter with priors stored in a TH2F
220  std::cout << "setting the TH2 " << std::endl;
221  PIDPriors priors3 = PIDPriors();
222  priors3.setPriors(part_e, counts1);
223  priors3.setPriors(part_mu, counts2);
224  priors3.setPriors(part_pi, counts3);
225  priors3.setPriors(part_K, counts4);
226  priors3.setPriors(part_p, counts5);
227  priors3.setPriors(part_d, counts6);
228 
229 
230  std::vector<Const::ChargedStable> codes = {part_e, part_mu, part_pi, part_K, part_p, part_d};
231  for (auto code : codes) {
232  for (int i = 0; i < 3; i++) {
233  double p = rn->Uniform(0., 10);
234  double c = rn->Uniform(-1., 1);
235  EXPECT_FLOAT_EQ(priors.getPriorValue(code, p, c), priors2.getPriorValue(code, p, c));
236  EXPECT_FLOAT_EQ(priors.getPriorValue(code, p, c), priors3.getPriorValue(code, p, c));
237  }
238  }
239 
240  EXPECT_FLOAT_EQ(0., priors.getPriorValue(part_pi, 11.1, 0.));
241  EXPECT_FLOAT_EQ(0., priors.getPriorValue(part_pi, -1, 0.5));
242  EXPECT_FLOAT_EQ(0., priors.getPriorValue(part_pi, 0.98, -1.001));
243 
244  EXPECT_FLOAT_EQ(0., priors2.getPriorValue(part_pi, 11.1, 0.));
245  EXPECT_FLOAT_EQ(0., priors2.getPriorValue(part_pi, -1, 0.5));
246  EXPECT_FLOAT_EQ(0., priors2.getPriorValue(part_pi, 0.98, -1.001));
247 
248  EXPECT_FLOAT_EQ(0., priors3.getPriorValue(part_pi, 11.1, 0.));
249  EXPECT_FLOAT_EQ(0., priors3.getPriorValue(part_pi, -1, 0.5));
250  EXPECT_FLOAT_EQ(0., priors3.getPriorValue(part_pi, 0.98, -1.001));
251 
252  }
253 
255 } // namespace
Belle2::PIDPriorsTest
Definition: PIDPriors.cc:23
Belle2::PIDPriors::setPriors
void setPriors(const Const::ChargedStable &particle, const PIDPriorsTable &table)
Sets the prior table for a particle species from a PIDPriorsTable object.
Definition: PIDPriors.h:54
Belle2::PIDPriorsTable::setBinEdges
void setBinEdges(const std::vector< float > &binEdgesX, const std::vector< float > &binEdgesY)
Sets the bin edges arrays.
Definition: PIDPriorsTable.cc:18
Belle2::PIDPriorsTable::setPriorsTable
void setPriorsTable(const std::vector< float > &priors)
Sets the priors table from a 2D std::vector.
Definition: PIDPriorsTable.h:54
Belle2::PIDPriorsTable::printPrior
void printPrior() const
Prints the content of the table and the axes.
Definition: PIDPriorsTable.cc:100
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::PIDPriors::getPriorValue
float getPriorValue(const Const::ChargedStable &particle, float x, float y) const
Returns the prior probability associated to a particle with defined species and parameters.
Definition: PIDPriors.h:179
Belle2::PIDPriorsTable::setPriorValue
void setPriorValue(float x, float y, float value)
Sets the prior value for a single bin.
Definition: PIDPriorsTable.cc:53
Belle2::PIDPriors
A database class to hold the prior probability for the particle identification.
Definition: PIDPriors.h:33
Belle2::rn
TString rn()
Get random string.
Definition: tools.h:41
Belle2::Const::ChargedStable
Provides a type-safe way to pass members of the chargedStableSet set.
Definition: Const.h:465
Belle2::TEST_F
TEST_F(PIDPriorsTest, PIDPriorTest)
Test the PIDPriors dbobject.
Definition: PIDPriors.cc:101
Belle2::PIDPriorsTable::getPriorValue
float getPriorValue(float x, float y) const
Returns the prior probability for a given value of (x,y).
Definition: PIDPriorsTable.cc:80
Belle2::PIDPriorsTable
This class holds the prior distribution for a single particle species.
Definition: PIDPriorsTable.h:37