Belle II Software light-2406-ragdoll
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#include <analysis/VariableManager/Manager.h>
10
11#include <analysis/dataobjects/Particle.h>
12#include <analysis/dataobjects/ParticleExtraInfoMap.h>
13#include <analysis/dataobjects/RestOfEvent.h>
14#include <analysis/utility/MCMatching.h>
15
16#include <framework/datastore/StoreArray.h>
17#include <framework/datastore/StoreObjPtr.h>
18#include <framework/logging/Logger.h>
19#include <framework/gearbox/Gearbox.h>
20
21#include <mdst/dataobjects/MCParticle.h>
22#include <mdst/dataobjects/PIDLikelihood.h>
23#include <mdst/dataobjects/Track.h>
24#include <mdst/dataobjects/ECLCluster.h>
25#include <mdst/dataobjects/KLMCluster.h>
26
27#include <gtest/gtest.h>
28
29#include <TMatrixFSym.h>
30
31using namespace std;
32using namespace Belle2;
33using namespace Belle2::Variable;
34
35namespace {
36
37 class FlavorTaggingVariablesTest : public ::testing::Test {
38 protected:
40 void SetUp() override
41 {
43 StoreArray<ECLCluster> testsECLClusters;
44 StoreArray<KLMCluster> testsKLMClusters;
46 StoreArray<Track> testsTracks;
47 StoreArray<Particle> testsParticles;
49 StoreArray<MCParticle> testsMCParticles;
50 StoreObjPtr<RestOfEvent> roe("RestOfEvent");
52 StoreArray<PIDLikelihood> testsPIDLikelihoods;
53 testsECLClusters.registerInDataStore();
54 testsKLMClusters.registerInDataStore();
55 testsTFRs.registerInDataStore();
56 testsTracks.registerInDataStore();
57 testsParticles.registerInDataStore();
58 extraInfoMap.registerInDataStore();
59 testsMCParticles.registerInDataStore();
60 roe.registerInDataStore();
61 testsROEs.registerInDataStore();
62 testsPIDLikelihoods.registerInDataStore();
63 testsParticles.registerRelationTo(testsROEs);
64 testsParticles.registerRelationTo(testsMCParticles);
65 testsTracks.registerRelationTo(testsPIDLikelihoods);
66 testsTracks.registerRelationTo(testsECLClusters);
67 testsTracks.registerRelationTo(testsMCParticles);
68 testsECLClusters.registerRelationTo(testsTracks);
70 }
71
73 void TearDown() override
74 {
76 }
77 };
78
80 TEST_F(FlavorTaggingVariablesTest, VariablesRunningForEachROE)
81 {
82
88 Gearbox& gearbox = Gearbox::getInstance();
89 gearbox.setBackends({std::string("file:")});
90 gearbox.close();
91 gearbox.open("geometry/Belle2.xml", false);
92
93
94 StoreArray<ECLCluster> testsECLClusters;
96 StoreArray<Track> testsTracks;
97 StoreArray<Particle> testsParticles;
98
99
101 vector<const Particle*> roeNeutralParticles;
102
104 vector<vector<int>> roeNeutralECLClusterIds{{4, 1}, {6, 1}, {9, 1}, {12, 1}, {17, 1}};
105
107 vector<vector<double>> roeNeutralECLClusterProperties{{0.0487526, 0.493606, -2.65695, 239.246},
108 {0.485431, 0.768093, 0.905049, 199.805},
109 {0.436428, 0.99517, 1.6799, 167.271},
110 {0.0230045, 1.05471, -1.20024, 161.013},
111 {0.174332, 1.49141, 1.3059, 141.107}};
112
114 for (unsigned i = 0; i < roeNeutralECLClusterIds.size(); ++i) {
115
116 ECLCluster ROEECL;
117 ROEECL.setIsTrack(false);
118 ROEECL.setHypothesis(ECLCluster::EHypothesisBit::c_nPhotons);
119 ROEECL.setConnectedRegionId(roeNeutralECLClusterIds[i][0]);
120 ROEECL.setClusterId(roeNeutralECLClusterIds[i][1]);
121 ROEECL.setEnergy(roeNeutralECLClusterProperties[i][0]);
122 ROEECL.setTheta(roeNeutralECLClusterProperties[i][1]);
123 ROEECL.setPhi(roeNeutralECLClusterProperties[i][2]);
124 ROEECL.setR(roeNeutralECLClusterProperties[i][3]);
125 ECLCluster* savedROEECL = testsECLClusters.appendNew(ROEECL);
126 Particle* roeECLParticle = testsParticles.appendNew(savedROEECL);
127 roeNeutralParticles.push_back(roeECLParticle);
128
129 }
130
132 vector<const Particle*> roeChargedParticles;
133 const float bField = 1.5;
134 TMatrixDSym cov6(6);
135
137 vector<vector<float>> roeTFRProperties{{ -0.578196, 0.252548, 0.158642, -0.520087, -1.19071, 0.417451, 0.139206},
138 {0.00434622, -0.0191058, 0.231542, 0.664887, 0.15125, 0.463762, 0.0313268},
139 {0.0738595, -0.00294903, 0.242319, 0.0186749, 0.467721, 0.288627, 0.0415814},
140 {1.3887, -2.43016, 1.00615, -0.27093, -0.15482, 0.0639155, 0.0132331},
141 {0.0041096, -0.0152487, 0.264326, 0.294061, 0.079251, 0.0389014, 0.00200114},
142 { -0.00371803, -5.22544e-05, 0.210148, 0.0019418, -0.138163, -0.0287232, 0.000622186},
143 {1.02296, 0.608721, 0.32273, -0.0505521, 0.0849532, 0.0839057, 0},
144 {6.49062, 3.20227, 117.684, 0.0469579, -0.0951783, -0.0706371, 0.117562}};
145
147 vector<short int> roeTRFCharges{1, 1, 1, -1, -1, -1, 1, -1};
148
150 vector<uint64_t> roeTRFCDCValues{3098476543630901248, 3026418949592973312, 3170534137668829184, 3386706919782612992,
151 3242591731706757120, 2954361355555045376, 1080863910568919040, 504403158265495552};
152
154 vector<uint32_t> roeTRFVXDValues{5570560, 5592320, 5592320, 5570560, 5592320, 5264640, 328960, 0};
155
156
158 vector<vector<int>> roeChargedECLClusterIds{{14, 1}, {10, 1}, {7, 1}, {8, 1}, {15, 1}};
159
161 vector<vector<double>> roeChargedECLClusterProperties{{0.964336, 1.23481, -2.25428, 148.729},
162 {0.214864, 0.965066, -0.232973, 170.008},
163 {0.0148855, 0.914396, 1.01693, 175.861},
164 {0.524092, 0.956389, 0.854331, 171.378},
165 {0.230255, 1.33317, -1.45326, 144.849}};
166
167 unsigned int chargedECLCLusterCounter = 0;
168
170 for (unsigned i = 0; i < roeTRFCharges.size(); ++i) {
171
172 ROOT::Math::XYZVector position(roeTFRProperties[i][0], roeTFRProperties[i][1], roeTFRProperties[i][2]);
173 ROOT::Math::XYZVector momentum(roeTFRProperties[i][3], roeTFRProperties[i][4], roeTFRProperties[i][5]);
174
175 testsTFRs.appendNew(position, momentum, cov6, roeTRFCharges[i], Const::pion, roeTFRProperties[i][6], bField, roeTRFCDCValues[i],
176 roeTRFVXDValues[i], 0);
177
178 Track ROETrack;
180 Track* savedROETrack = testsTracks.appendNew(ROETrack);
181
183 if (i == 0 || i == 1 || i == 2 || i == 3) {
184 ECLCluster ROEChargedECL;
185 ROEChargedECL.setIsTrack(true);
186 ROEChargedECL.setHypothesis(ECLCluster::EHypothesisBit::c_nPhotons);
187 ROEChargedECL.setConnectedRegionId(roeChargedECLClusterIds[chargedECLCLusterCounter][0]);
188 ROEChargedECL.setClusterId(roeChargedECLClusterIds[chargedECLCLusterCounter][1]);
189 ROEChargedECL.setEnergy(roeChargedECLClusterProperties[chargedECLCLusterCounter][0]);
190 ROEChargedECL.setTheta(roeChargedECLClusterProperties[chargedECLCLusterCounter][1]);
191 ROEChargedECL.setPhi(roeChargedECLClusterProperties[chargedECLCLusterCounter][2]);
192 ROEChargedECL.setR(roeChargedECLClusterProperties[chargedECLCLusterCounter][3]);
193 ECLCluster* savedROEChargedECL = testsECLClusters.appendNew(ROEChargedECL);
194 savedROEChargedECL->addRelationTo(savedROETrack);
195 chargedECLCLusterCounter++;
196
197 if (i == 2) {
199 ECLCluster ROEChargedECL2;
200 ROEChargedECL2.setIsTrack(true);
201 ROEChargedECL2.setHypothesis(ECLCluster::EHypothesisBit::c_nPhotons);
202 ROEChargedECL2.setConnectedRegionId(roeChargedECLClusterIds[chargedECLCLusterCounter][0]);
203 ROEChargedECL2.setClusterId(roeChargedECLClusterIds[chargedECLCLusterCounter][1]);
204 ROEChargedECL2.setEnergy(roeChargedECLClusterProperties[chargedECLCLusterCounter][0]);
205 ROEChargedECL2.setTheta(roeChargedECLClusterProperties[chargedECLCLusterCounter][1]);
206 ROEChargedECL2.setPhi(roeChargedECLClusterProperties[chargedECLCLusterCounter][2]);
207 ROEChargedECL2.setR(roeChargedECLClusterProperties[chargedECLCLusterCounter][3]);
208 ECLCluster* savedROEChargedECL2 = testsECLClusters.appendNew(ROEChargedECL2);
209 savedROETrack->addRelationTo(savedROEChargedECL2);
210 chargedECLCLusterCounter++;
211
212 }
213 }
214 Particle* roeTrackParticle = testsParticles.appendNew(savedROETrack, Const::pion);
215 roeChargedParticles.push_back(roeTrackParticle);
216 }
217
219 StoreObjPtr<RestOfEvent> roe("RestOfEvent");
220 roe.create();
221 roe->addParticles(roeNeutralParticles);
222 roe->addParticles(roeChargedParticles);
223
224 B2INFO("Is the ROE valid? " << roe.isValid());
225 ASSERT_TRUE(roe.isValid());
226
227
229 ROOT::Math::PxPyPzEVector roe1FourVectorECLClusters = roe->get4VectorNeutralECLClusters();
230
231 B2INFO("The total four momentum of the neutral clusters in the first test ROE is = ("
232 << roe1FourVectorECLClusters.E() << ", "
233 << roe1FourVectorECLClusters.X() << ", "
234 << roe1FourVectorECLClusters.Y() << ", "
235 << roe1FourVectorECLClusters.Z() << ")");
236
237 EXPECT_NEAR(roe1FourVectorECLClusters.E(), 1.16795, 0.00005);
238 EXPECT_NEAR(roe1FourVectorECLClusters.X(), 0.20075, 0.00005);
239 EXPECT_NEAR(roe1FourVectorECLClusters.Y(), 0.76747, 0.00005);
240 EXPECT_NEAR(roe1FourVectorECLClusters.Z(), 0.65482, 0.00005);
241
243 const Manager::Var* var = Manager::Instance().getVariable("BtagToWBosonVariables(recoilMass)");
244
247 ASSERT_NE(var, nullptr);
248
251 vector<double> refsBtagToWBosonRecoilMass{3.2093, 4.1099, 4.3019, 4.4436, 4.4924, 4.5995, 4.6493, 4.5963};
252
255 for (unsigned i = 0; i < roeChargedParticles.size(); i++) {
256
257 double output = std::get<double>(var->function(roeChargedParticles[i]));
258
262 EXPECT_NEAR(output, refsBtagToWBosonRecoilMass[i], 0.0005);
263
264 }
265
268 var = Manager::Instance().getVariable("BtagToWBosonVariables(recoilMassSqrd)");
269 ASSERT_NE(var, nullptr);
270
271 vector<double> refsBtagToWBosonRecoilMassSqrd{10.300, 16.891, 18.506, 19.746, 20.182, 21.155, 21.616, 21.126};
272
273 for (unsigned i = 0; i < roeChargedParticles.size(); i++) {
274
275 double output = std::get<double>(var->function(roeChargedParticles[i]));
276
277 EXPECT_NEAR(output, refsBtagToWBosonRecoilMassSqrd[i], 0.0005);
278
279 }
280
281 var = Manager::Instance().getVariable("BtagToWBosonVariables(pMissCMS)");
282 ASSERT_NE(var, nullptr);
283
284 double refsBtagToWBosonPMissCMS = 0.5425695;
285
286 for (auto& roeChargedParticle : roeChargedParticles) {
287
288 double output = std::get<double>(var->function(roeChargedParticle));
289
290 EXPECT_NEAR(output, refsBtagToWBosonPMissCMS, 0.000005);
291
292 }
293
294 var = Manager::Instance().getVariable("BtagToWBosonVariables(cosThetaMissCMS)");
295 ASSERT_NE(var, nullptr);
296
297 vector<double> refsBtagToWBosonCosThetaMissCMS{0.0621, -0.6164, -0.2176, 0.3516, -0.1165, 0.4694, -0.0754, 0.6302};
298
299 for (unsigned i = 0; i < roeChargedParticles.size(); i++) {
300
301 double output = std::get<double>(var->function(roeChargedParticles[i]));
302
303 EXPECT_NEAR(output, refsBtagToWBosonCosThetaMissCMS[i], 0.0005);
304
305 }
306
307
308 var = Manager::Instance().getVariable("BtagToWBosonVariables(EW90)");
309 ASSERT_NE(var, nullptr);
310
311 vector<double> refsBtagToWBosonEW90{0.3020, 1.63517, 1.09619, 0.96434, 0.17433, 1.19459, 1.13867, 1.36892};
312
313 for (unsigned i = 0; i < roeChargedParticles.size(); i++) {
314
315 double output = std::get<double>(var->function(roeChargedParticles[i]));
316
317 EXPECT_NEAR(output, refsBtagToWBosonEW90[i], 0.0005);
318
319 }
320
321 }
322
323
325 TEST_F(FlavorTaggingVariablesTest, isSignalVariable)
326 {
327
335 StoreArray<Particle> testsParticles;
336 StoreArray<MCParticle> testsMCParticles;
337
339 MCParticle MCB0;
340 MCB0.setEnergy(5.68161);
341 MCB0.setPDG(511);
342 MCB0.setMassFromPDG();
343 MCB0.setMomentum(-0.12593, -0.143672, 2.09072);
344 MCParticle* savedMCB0 = testsMCParticles.appendNew(MCB0);
345
347 Particle B0({ -0.129174, -0.148899, 2.09292, 5.67644}, 511);
348 Particle* savedB0 = testsParticles.appendNew(B0);
349 savedB0->addRelationTo(savedMCB0);
350
352 const Manager::Var* var = Manager::Instance().getVariable("isSignal");
353 ASSERT_NE(var, nullptr);
354
357 savedB0->addExtraInfo(MCMatching::c_extraInfoMCErrors, 0);
358 double output1 = std::get<double>(var->function(savedB0));
359 ASSERT_EQ(output1, 1.0);
360
365 savedB0->setExtraInfo(MCMatching::c_extraInfoMCErrors, 0);
366 savedB0->setProperty(6);
367 double output2 = std::get<double>(var->function(savedB0));
368 ASSERT_EQ(output2, 1.0);
369
372 vector<int> notAcceptedMCErrorFlags{4, 8, 16, 32, 64, 128, 256, 512};
373
374 for (int notAcceptedMCErrorFlag : notAcceptedMCErrorFlags) {
375
376 savedB0->setExtraInfo(MCMatching::c_extraInfoMCErrors, notAcceptedMCErrorFlag);
377 double output = std::get<double>(var->function(savedB0));
378 ASSERT_EQ(output, 0);
379
380 }
381
382 }
383
384}
static const ChargedStable pion
charged pion particle
Definition: Const.h:661
static DataStore & Instance()
Instance of singleton Store.
Definition: DataStore.cc:54
void setInitializeActive(bool active)
Setter for m_initializeActive.
Definition: DataStore.cc:94
void reset(EDurability durability)
Frees memory occupied by data store items and removes all objects from the map.
Definition: DataStore.cc:86
ECL cluster data.
Definition: ECLCluster.h:27
void setTheta(double theta)
Set Theta of Shower (radian).
Definition: ECLCluster.h:217
void setConnectedRegionId(int crid)
Set connected region id.
Definition: ECLCluster.h:142
void setPhi(double phi)
Set Phi of Shower (radian).
Definition: ECLCluster.h:220
void setClusterId(int clusterid)
Set cluster id.
Definition: ECLCluster.h:145
void setHypothesis(EHypothesisBit hypothesis)
Set hypotheses.
Definition: ECLCluster.h:123
void setEnergy(double energy)
Set Corrected Energy (GeV).
Definition: ECLCluster.h:226
void setIsTrack(bool istrack)
Set m_isTrack true if the cluster matches with a track.
Definition: ECLCluster.h:104
void setR(double r)
Set R (in cm).
Definition: ECLCluster.h:223
Singleton class responsible for loading detector parameters from an XML file.
Definition: Gearbox.h:34
A Class to store the Monte Carlo particle information.
Definition: MCParticle.h:32
void setEnergy(float energy)
Set energy.
Definition: MCParticle.h:372
void setPDG(int pdg)
Set PDG code of the particle.
Definition: MCParticle.h:335
void setMomentum(const ROOT::Math::XYZVector &momentum)
Set particle momentum.
Definition: MCParticle.h:417
void setMassFromPDG()
Sets the mass for the particle from the particle's PDG code.
Definition: MCParticle.cc:28
Class to store reconstructed particles.
Definition: Particle.h:75
void addRelationTo(const RelationsInterface< BASE > *object, float weight=1.0, const std::string &namedRelation="") const
Add a relation from this object to another object (with caching).
bool registerInDataStore(DataStore::EStoreFlags storeFlags=DataStore::c_WriteOut)
Register the object/array in the DataStore.
Accessor to arrays stored in the data store.
Definition: StoreArray.h:113
T * appendNew()
Construct a new T object at the end of the array.
Definition: StoreArray.h:246
bool registerRelationTo(const StoreArray< TO > &toArray, DataStore::EDurability durability=DataStore::c_Event, DataStore::EStoreFlags storeFlags=DataStore::c_WriteOut, const std::string &namedRelation="") const
Register a relation to the given StoreArray.
Definition: StoreArray.h:140
Type-safe access to single objects in the data store.
Definition: StoreObjPtr.h:96
Class that bundles various TrackFitResults.
Definition: Track.h:25
void setTrackFitResultIndex(const Const::ChargedStable &chargedStable, short index)
Set an index (for positive values) or unavailability-code (index = -1) for a specific mass hypothesis...
Definition: Track.h:174
const Var * getVariable(std::string name)
Get the variable belonging to the given key.
Definition: Manager.cc:57
static Manager & Instance()
get singleton instance.
Definition: Manager.cc:25
static Gearbox & getInstance()
Return reference to the Gearbox instance.
Definition: Gearbox.cc:81
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24
STL namespace.
static const std::string c_extraInfoMCErrors
Name of extra-info field stored in Particle.
Definition: MCMatching.h:30
A variable returning a floating-point value for a given Particle.
Definition: Manager.h:146