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