Belle II Software  release-05-02-19
variables.cc
1 
2 #include <analysis/variables/Variables.h>
3 #include <analysis/variables/BasicParticleInformation.h>
4 #include <analysis/variables/VertexVariables.h>
5 #include <analysis/variables/PIDVariables.h>
6 #include <analysis/variables/TrackVariables.h>
7 
8 #include <analysis/VariableManager/Manager.h>
9 #include <analysis/VariableManager/Utility.h>
10 
11 #include <analysis/dataobjects/Particle.h>
12 #include <analysis/dataobjects/ParticleExtraInfoMap.h>
13 #include <analysis/dataobjects/ParticleList.h>
14 #include <analysis/dataobjects/EventExtraInfo.h>
15 #include <analysis/dataobjects/RestOfEvent.h>
16 #include <analysis/utility/ReferenceFrame.h>
17 
18 #include <framework/datastore/StoreArray.h>
19 #include <framework/datastore/StoreObjPtr.h>
20 #include <framework/utilities/TestHelpers.h>
21 #include <framework/logging/Logger.h>
22 #include <framework/gearbox/Gearbox.h>
23 
24 #include <mdst/dataobjects/MCParticle.h>
25 #include <mdst/dataobjects/MCParticleGraph.h>
26 #include <mdst/dataobjects/PIDLikelihood.h>
27 #include <mdst/dataobjects/Track.h>
28 #include <mdst/dataobjects/V0.h>
29 #include <mdst/dataobjects/ECLCluster.h>
30 #include <mdst/dataobjects/KLMCluster.h>
31 
32 #include <gtest/gtest.h>
33 
34 #include <TMatrixFSym.h>
35 #include <TRandom3.h>
36 #include <TLorentzVector.h>
37 #include <TMath.h>
38 #include <utility>
39 
40 using namespace std;
41 using namespace Belle2;
42 using namespace Belle2::Variable;
43 
44 namespace {
45 
47  TEST(KinematicVariableTest, Variable)
48  {
49 
50  // Connect gearbox for CMS variables
51 
52  Gearbox& gearbox = Gearbox::getInstance();
53  gearbox.setBackends({std::string("file:")});
54  gearbox.close();
55  gearbox.open("geometry/Belle2.xml", false);
56 
57  {
58  Particle p({ 0.1 , -0.4, 0.8, 1.0 }, 11);
59 
60  TMatrixFSym error(7);
61  error.Zero();
62  error(0, 0) = 0.05;
63  error(1, 1) = 0.2;
64  error(2, 2) = 0.4;
65  error(0, 1) = -0.1;
66  error(0, 2) = 0.9;
67  p.setMomentumVertexErrorMatrix(error);
68 
69  EXPECT_FLOAT_EQ(0.9, particleP(&p));
70  EXPECT_FLOAT_EQ(1.0, particleE(&p));
71  EXPECT_FLOAT_EQ(0.1, particlePx(&p));
72  EXPECT_FLOAT_EQ(-0.4, particlePy(&p));
73  EXPECT_FLOAT_EQ(0.8, particlePz(&p));
74  EXPECT_FLOAT_EQ(0.412310562, particlePt(&p));
75  EXPECT_FLOAT_EQ(0.8 / 0.9, particleCosTheta(&p));
76  EXPECT_FLOAT_EQ(-1.325817664, particlePhi(&p));
77 
78  EXPECT_FLOAT_EQ(0.737446378, particlePErr(&p));
79  EXPECT_FLOAT_EQ(sqrt(0.05), particlePxErr(&p));
80  EXPECT_FLOAT_EQ(sqrt(0.2), particlePyErr(&p));
81  EXPECT_FLOAT_EQ(sqrt(0.4), particlePzErr(&p));
82  EXPECT_FLOAT_EQ(0.488093530, particlePtErr(&p));
83  EXPECT_FLOAT_EQ(0.156402664, particleCosThetaErr(&p));
84  EXPECT_FLOAT_EQ(0.263066820, particlePhiErr(&p));
85 
86 
87  {
88  UseReferenceFrame<CMSFrame> dummy;
89  EXPECT_FLOAT_EQ(0.68176979, particleP(&p));
90  EXPECT_FLOAT_EQ(0.80920333, particleE(&p));
91  EXPECT_FLOAT_EQ(0.061728548, particlePx(&p));
92  EXPECT_FLOAT_EQ(-0.40000001, particlePy(&p));
93  EXPECT_FLOAT_EQ(0.54863429, particlePz(&p));
94  EXPECT_FLOAT_EQ(0.404735, particlePt(&p));
95  EXPECT_FLOAT_EQ(0.80472076, particleCosTheta(&p));
96  EXPECT_FLOAT_EQ(-1.4176828, particlePhi(&p));
97 
98  EXPECT_FLOAT_EQ(sqrt(0.2), particlePyErr(&p));
99  }
100 
101  {
102  UseReferenceFrame<RestFrame> dummy(&p);
103  EXPECT_ALL_NEAR(particleP(&p), 0.0, 1e-9);
104  EXPECT_FLOAT_EQ(0.4358899, particleE(&p));
105  EXPECT_ALL_NEAR(0.0, particlePx(&p), 1e-9);
106  EXPECT_ALL_NEAR(0.0, particlePy(&p), 1e-9);
107  EXPECT_ALL_NEAR(0.0, particlePz(&p), 1e-9);
108  EXPECT_ALL_NEAR(0.0, particlePt(&p), 1e-9);
109 
110  }
111 
112  {
113  UseReferenceFrame<LabFrame> dummy;
114  EXPECT_FLOAT_EQ(0.9, particleP(&p));
115  EXPECT_FLOAT_EQ(1.0, particleE(&p));
116  EXPECT_FLOAT_EQ(0.1, particlePx(&p));
117  EXPECT_FLOAT_EQ(-0.4, particlePy(&p));
118  EXPECT_FLOAT_EQ(0.8, particlePz(&p));
119  EXPECT_FLOAT_EQ(0.412310562, particlePt(&p));
120  EXPECT_FLOAT_EQ(0.8 / 0.9, particleCosTheta(&p));
121  EXPECT_FLOAT_EQ(-1.325817664, particlePhi(&p));
122 
123  EXPECT_FLOAT_EQ(0.737446378, particlePErr(&p));
124  EXPECT_FLOAT_EQ(sqrt(0.05), particlePxErr(&p));
125  EXPECT_FLOAT_EQ(sqrt(0.2), particlePyErr(&p));
126  EXPECT_FLOAT_EQ(sqrt(0.4), particlePzErr(&p));
127  EXPECT_FLOAT_EQ(0.488093530, particlePtErr(&p));
128  EXPECT_FLOAT_EQ(0.156402664, particleCosThetaErr(&p));
129  EXPECT_FLOAT_EQ(0.263066820, particlePhiErr(&p));
130  }
131 
132  {
133  UseReferenceFrame<RotationFrame> dummy(TVector3(1, 0, 0), TVector3(0, 1, 0), TVector3(0, 0, 1));
134  EXPECT_FLOAT_EQ(0.9, particleP(&p));
135  EXPECT_FLOAT_EQ(1.0, particleE(&p));
136  EXPECT_FLOAT_EQ(0.1, particlePx(&p));
137  EXPECT_FLOAT_EQ(-0.4, particlePy(&p));
138  EXPECT_FLOAT_EQ(0.8, particlePz(&p));
139  EXPECT_FLOAT_EQ(0.412310562, particlePt(&p));
140  EXPECT_FLOAT_EQ(0.8 / 0.9, particleCosTheta(&p));
141  EXPECT_FLOAT_EQ(-1.325817664, particlePhi(&p));
142 
143  EXPECT_FLOAT_EQ(0.737446378, particlePErr(&p));
144  EXPECT_FLOAT_EQ(sqrt(0.05), particlePxErr(&p));
145  EXPECT_FLOAT_EQ(sqrt(0.2), particlePyErr(&p));
146  EXPECT_FLOAT_EQ(sqrt(0.4), particlePzErr(&p));
147  EXPECT_FLOAT_EQ(0.488093530, particlePtErr(&p));
148  EXPECT_FLOAT_EQ(0.156402664, particleCosThetaErr(&p));
149  EXPECT_FLOAT_EQ(0.263066820, particlePhiErr(&p));
150 
151  const auto& frame = ReferenceFrame::GetCurrent();
152  EXPECT_FLOAT_EQ(-0.1, frame.getMomentumErrorMatrix(&p)(0, 1));
153  EXPECT_FLOAT_EQ(0.9, frame.getMomentumErrorMatrix(&p)(0, 2));
154  }
155 
156  {
157  UseReferenceFrame<RotationFrame> dummy(TVector3(1, 0, 0), TVector3(0, 0, -1), TVector3(0, 1, 0));
158  EXPECT_FLOAT_EQ(0.9, particleP(&p));
159  EXPECT_FLOAT_EQ(1.0, particleE(&p));
160  EXPECT_FLOAT_EQ(0.1, particlePx(&p));
161  EXPECT_FLOAT_EQ(-0.8, particlePy(&p));
162  EXPECT_FLOAT_EQ(-0.4, particlePz(&p));
163 
164  EXPECT_FLOAT_EQ(0.737446378, particlePErr(&p));
165  EXPECT_FLOAT_EQ(sqrt(0.05), particlePxErr(&p));
166  EXPECT_FLOAT_EQ(sqrt(0.4), particlePyErr(&p));
167  EXPECT_FLOAT_EQ(sqrt(0.2), particlePzErr(&p));
168 
169  const auto& frame = ReferenceFrame::GetCurrent();
170  EXPECT_FLOAT_EQ(-0.9, frame.getMomentumErrorMatrix(&p)(0, 1));
171  EXPECT_FLOAT_EQ(-0.1, frame.getMomentumErrorMatrix(&p)(0, 2));
172  }
173 
174  {
175  UseReferenceFrame<CMSRotationFrame> dummy(TVector3(1, 0, 0), TVector3(0, 1, 0), TVector3(0, 0, 1));
176  EXPECT_FLOAT_EQ(0.68176979, particleP(&p));
177  EXPECT_FLOAT_EQ(0.80920333, particleE(&p));
178  EXPECT_FLOAT_EQ(0.061728548, particlePx(&p));
179  EXPECT_FLOAT_EQ(-0.40000001, particlePy(&p));
180  EXPECT_FLOAT_EQ(0.54863429, particlePz(&p));
181  EXPECT_FLOAT_EQ(0.404735, particlePt(&p));
182  EXPECT_FLOAT_EQ(0.80472076, particleCosTheta(&p));
183  EXPECT_FLOAT_EQ(-1.4176828, particlePhi(&p));
184 
185  EXPECT_FLOAT_EQ(sqrt(0.2), particlePyErr(&p));
186  }
187 
188  {
189  Particle pinv({ -0.1 , 0.4, -0.8, 1.0 }, 11);
190  UseReferenceFrame<RestFrame> dummy(&pinv);
191  Particle p2({ 0.0 , 0.0, 0.0, 0.4358899}, 11);
192  EXPECT_FLOAT_EQ(0.9, particleP(&p2));
193  EXPECT_FLOAT_EQ(1.0, particleE(&p2));
194  EXPECT_FLOAT_EQ(0.1, particlePx(&p2));
195  EXPECT_FLOAT_EQ(-0.4, particlePy(&p2));
196  EXPECT_FLOAT_EQ(0.8, particlePz(&p2));
197  EXPECT_FLOAT_EQ(0.412310562, particlePt(&p2));
198  EXPECT_FLOAT_EQ(0.8 / 0.9, particleCosTheta(&p2));
199  EXPECT_FLOAT_EQ(-1.325817664, particlePhi(&p2));
200  }
201  }
202 
203  {
204  Particle p({ 0.0 , 0.0, 0.0, 0.0 }, 11);
205  EXPECT_FLOAT_EQ(0.0, particleP(&p));
206  EXPECT_FLOAT_EQ(0.0, particleE(&p));
207  EXPECT_FLOAT_EQ(0.0, particlePx(&p));
208  EXPECT_FLOAT_EQ(0.0, particlePy(&p));
209  EXPECT_FLOAT_EQ(0.0, particlePz(&p));
210  EXPECT_FLOAT_EQ(0.0, particlePt(&p));
211  EXPECT_FLOAT_EQ(1.0, particleCosTheta(&p));
212  EXPECT_FLOAT_EQ(0.0, particlePhi(&p));
213 
214  UseReferenceFrame<CMSFrame> dummy;
215  EXPECT_FLOAT_EQ(0.0, particleP(&p));
216  EXPECT_FLOAT_EQ(0.0, particleE(&p));
217  EXPECT_FLOAT_EQ(0.0, particlePx(&p));
218  EXPECT_FLOAT_EQ(0.0, particlePy(&p));
219  EXPECT_FLOAT_EQ(0.0, particlePz(&p));
220  EXPECT_FLOAT_EQ(0.0, particlePt(&p));
221  EXPECT_FLOAT_EQ(1.0, particleCosTheta(&p));
222  EXPECT_FLOAT_EQ(0.0, particlePhi(&p));
223  }
224 
225  {
226  DataStore::Instance().setInitializeActive(true);
227  StoreArray<Particle> particles;
228  particles.registerInDataStore();
229  DataStore::Instance().setInitializeActive(false);
230  PCmsLabTransform T;
231  TLorentzVector vec0 = {0.0, 0.0, 0.0, T.getCMSEnergy()};
232  TLorentzVector vec1 = {0.0, +0.332174566, 0.0, T.getCMSEnergy() / 2.};
233  TLorentzVector vec2 = {0.0, -0.332174566, 0.0, T.getCMSEnergy() / 2.};
234  Particle* p0 = particles.appendNew(Particle(T.rotateCmsToLab() * vec0, 22));
235  Particle* p1 = particles.appendNew(Particle(T.rotateCmsToLab() * vec1, 22, Particle::c_Unflavored, Particle::c_Undefined, 1));
236  Particle* p2 = particles.appendNew(Particle(T.rotateCmsToLab() * vec2, 22, Particle::c_Unflavored, Particle::c_Undefined, 2));
237 
238  p0->appendDaughter(p1->getArrayIndex());
239  p0->appendDaughter(p2->getArrayIndex());
240 
241  EXPECT_ALL_NEAR(m2RecoilSignalSide(p0), 0.0, 1e-7);
242  }
243 
244 
245  }
246 
247 
248  TEST(VertexVariableTest, Variable)
249  {
250 
251  // Connect gearbox for CMS variables
252 
253  Gearbox& gearbox = Gearbox::getInstance();
254  gearbox.setBackends({std::string("file:")});
255  gearbox.close();
256  gearbox.open("geometry/Belle2.xml", false);
257 
258  Particle p({ 0.1 , -0.4, 0.8, 1.0 }, 11);
259  p.setPValue(0.5);
260  p.setVertex(TVector3(1.0, 2.0, 2.0));
261 
262  EXPECT_FLOAT_EQ(1.0, particleDX(&p));
263  EXPECT_FLOAT_EQ(2.0, particleDY(&p));
264  EXPECT_FLOAT_EQ(2.0, particleDZ(&p));
265  EXPECT_FLOAT_EQ(std::sqrt(5.0), particleDRho(&p));
266  EXPECT_FLOAT_EQ(3.0, particleDistance(&p));
267  EXPECT_FLOAT_EQ(0.5, particlePvalue(&p));
268 
269  {
270  UseReferenceFrame<CMSFrame> dummy;
271  EXPECT_FLOAT_EQ(1.0382183, particleDX(&p));
272  EXPECT_FLOAT_EQ(2.0, particleDY(&p));
273  EXPECT_FLOAT_EQ(2.2510159, particleDZ(&p));
274  EXPECT_FLOAT_EQ(std::sqrt(2.0 * 2.0 + 1.0382183 * 1.0382183), particleDRho(&p));
275  EXPECT_FLOAT_EQ(3.185117, particleDistance(&p));
276  EXPECT_FLOAT_EQ(0.5, particlePvalue(&p));
277  }
278 
279  {
280  Particle p2({ 0.1 , -0.4, 0.8, 1.0 }, 11);
281  p2.setPValue(0.5);
282  p2.setVertex(TVector3(1.0, 2.0, 2.0));
283 
284  UseReferenceFrame<RestFrame> dummy(&p2);
285  EXPECT_FLOAT_EQ(0.0, particleDX(&p));
286  EXPECT_FLOAT_EQ(0.0, particleDY(&p));
287  EXPECT_FLOAT_EQ(0.0, particleDZ(&p));
288  EXPECT_FLOAT_EQ(0.0, particleDRho(&p));
289  EXPECT_FLOAT_EQ(0.0, particleDistance(&p));
290  EXPECT_FLOAT_EQ(0.5, particlePvalue(&p));
291  }
292 
293  /* Test with a distance between mother and daughter vertex. One
294  * has to calculate the result by hand to test the code....
295 
296  {
297  Particle p2({ 0.0 , 1.0, 0.0, 1.0 }, 11);
298  p2.setPValue(0.5);
299  p2.setVertex(TVector3(1.0, 0.0, 2.0));
300 
301  UseReferenceFrame<RestFrame> dummy(&p2);
302  EXPECT_FLOAT_EQ(0.0, particleDX(&p));
303  EXPECT_FLOAT_EQ(2.0, particleDY(&p));
304  EXPECT_FLOAT_EQ(0.0, particleDZ(&p));
305  EXPECT_FLOAT_EQ(2.0, particleDRho(&p));
306  EXPECT_FLOAT_EQ(2.0, particleDistance(&p));
307  EXPECT_FLOAT_EQ(0.5, particlePvalue(&p));
308  }
309  */
310 
311  }
312 
313  TEST(TrackVariablesTest, Variable)
314  {
315  DataStore::Instance().setInitializeActive(true);
316  StoreArray<TrackFitResult> myResults;
317  StoreArray<Track> myTracks;
318  StoreArray<V0> myV0s;
319  StoreArray<Particle> myParticles;
320  myResults.registerInDataStore();
321  myTracks.registerInDataStore();
322  myV0s.registerInDataStore();
323  myParticles.registerInDataStore();
324  DataStore::Instance().setInitializeActive(false);
325 
326  TRandom3 generator;
327 
328  const float pValue = 0.5;
329  const float bField = 1.5;
330  const int charge = 1;
331  TMatrixDSym cov6(6);
332 
333  // Generate a random put orthogonal pair of vectors in the r-phi plane
334  TVector2 d(generator.Uniform(-1, 1), generator.Uniform(-1, 1));
335  TVector2 pt(generator.Uniform(-1, 1), generator.Uniform(-1, 1));
336  d.Set(d.X(), -(d.X()*pt.Px()) / pt.Py());
337 
338  // Add a random z component
339  TVector3 position(d.X(), d.Y(), generator.Uniform(-1, 1));
340  TVector3 momentum(pt.Px(), pt.Py(), generator.Uniform(-1, 1));
341 
342  auto CDCValue = static_cast<unsigned long long int>(0x300000000000000);
343 
344  myResults.appendNew(position, momentum, cov6, charge, Const::electron, pValue, bField, CDCValue, 16777215, 0);
345  Track mytrack;
346  mytrack.setTrackFitResultIndex(Const::electron, 0);
347  Track* savedTrack = myTracks.appendNew(mytrack);
348 
349  Particle* part = myParticles.appendNew(savedTrack, Const::ChargedStable(11));
350 
351  const Manager::Var* vIsFromECL = Manager::Instance().getVariable("isFromECL");
352  const Manager::Var* vIsFromKLM = Manager::Instance().getVariable("isFromKLM");
353  const Manager::Var* vIsFromTrack = Manager::Instance().getVariable("isFromTrack");
354  const Manager::Var* vIsFromV0 = Manager::Instance().getVariable("isFromV0");
355 
356  EXPECT_TRUE(vIsFromTrack->function(part));
357  EXPECT_FALSE(vIsFromECL->function(part));
358  EXPECT_FALSE(vIsFromKLM->function(part));
359  EXPECT_FALSE(vIsFromV0->function(part));
360  EXPECT_FLOAT_EQ(0.5, trackPValue(part));
361  EXPECT_FLOAT_EQ(position.Z(), trackZ0(part));
362  EXPECT_FLOAT_EQ(sqrt(pow(position.X(), 2) + pow(position.Y(), 2)), trackD0(part));
363  EXPECT_FLOAT_EQ(3, trackNCDCHits(part));
364  EXPECT_FLOAT_EQ(24, trackNSVDHits(part));
365  EXPECT_FLOAT_EQ(12, trackNPXDHits(part));
366 
367  //-----------------------------------------------------------------------
368  // now add another track and mock up a V0 and a V0-based particle
369  myResults.appendNew(position, momentum, cov6, charge * -1,
370  Const::electron, pValue, bField, CDCValue, 16777215, 0);
371  Track secondTrack;
372  secondTrack.setTrackFitResultIndex(Const::electron, 1);
373  Track* savedTrack2 = myTracks.appendNew(secondTrack);
374  myParticles.appendNew(savedTrack2, Const::ChargedStable(11));
375  myV0s.appendNew(V0(std::pair(savedTrack, myResults[0]), std::pair(savedTrack2, myResults[1])));
376  const TLorentzVector v0Momentum(momentum * 2, (momentum * 2).Mag());
377  auto v0particle = myParticles.appendNew(v0Momentum, 22,
378  Particle::c_Unflavored, Particle::c_V0, 0);
379  v0particle->appendDaughter(0, false);
380  v0particle->appendDaughter(1, false);
381  //-----------------------------------------------------------------------
382 
383  EXPECT_FALSE(vIsFromTrack->function(v0particle));
384  EXPECT_FALSE(vIsFromECL->function(v0particle));
385  EXPECT_FALSE(vIsFromKLM->function(v0particle));
386  EXPECT_TRUE(vIsFromV0->function(v0particle));
387 
388  const Manager::Var* vNDaughters = Manager::Instance().getVariable("nDaughters");
389  EXPECT_FLOAT_EQ(vNDaughters->function(v0particle), 2);
390  }
391 
392  class MCTruthVariablesTest : public ::testing::Test {
393  protected:
394  virtual void SetUp()
395  {
396  // datastore things
397  DataStore::Instance().reset();
398  DataStore::Instance().setInitializeActive(true);
399 
400  // needed to mock up
401  StoreArray<ECLCluster> clusters;
402  StoreArray<MCParticle> mcparticles;
403  StoreArray<Track> tracks;
404  StoreArray<TrackFitResult> trackfits;
405  StoreArray<Particle> particles;
406 
407  // register the arrays
408  clusters.registerInDataStore();
409  mcparticles.registerInDataStore();
410  tracks.registerInDataStore();
411  trackfits.registerInDataStore();
412  particles.registerInDataStore();
413 
414  // register the relations for mock up mcmatching
415  clusters.registerRelationTo(mcparticles);
416  tracks.registerRelationTo(mcparticles);
417  particles.registerRelationTo(mcparticles);
418 
419  // register the relation for mock up track <--> cluster matching
420  //clusters.registerRelationTo(tracks);
421  tracks.registerRelationTo(clusters);
422 
423  // end datastore things
424  DataStore::Instance().setInitializeActive(false);
425 
426  /* mock up an electron (track with a cluster AND a track-cluster match)
427  * and a photon (cluster, no track) and MCParticles for both
428  *
429  * this assumes that everything (tracking, clustering, track-cluster
430  * matching *and* mcmatching all worked)
431  *
432  * this can be extended to pions, kaons, etc but leave it simple for now
433  */
434 
435  // create the true underlying mcparticles
436  auto* true_photon = mcparticles.appendNew(MCParticle());
437  true_photon->setPDG(22);
438  auto* true_electron = mcparticles.appendNew(MCParticle());
439  true_electron->setPDG(11);
440  auto* true_pion = mcparticles.appendNew(MCParticle());
441  true_pion->setPDG(-211);
442 
443  // create the reco clusters
444  auto* cl0 = clusters.appendNew(ECLCluster());
445  cl0->setEnergy(1.0);
446  cl0->setHypothesis(ECLCluster::EHypothesisBit::c_nPhotons);
447  cl0->setClusterId(0);
448 
449  auto* cl1 = clusters.appendNew(ECLCluster());
450  cl1->setEnergy(0.5);
451  cl1->setHypothesis(ECLCluster::EHypothesisBit::c_nPhotons);
452  cl1->setClusterId(1);
453 
454  // create a reco track (one has to also mock up a track fit result)
455  TMatrixDSym cov(6);
456  trackfits.appendNew(
457  TVector3(), TVector3(), cov, -1, Const::electron, 0.5, 1.5,
458  static_cast<unsigned long long int>(0x300000000000000), 16777215, 0);
459  auto* electron_tr = tracks.appendNew(Track());
460  electron_tr->setTrackFitResultIndex(Const::electron, 0);
461  electron_tr->addRelationTo(cl1); // a track <--> cluster match
462 
463  TMatrixDSym cov1(6);
464  trackfits.appendNew(
465  TVector3(), TVector3(), cov1, -1, Const::pion, 0.51, 1.5,
466  static_cast<unsigned long long int>(0x300000000000000), 16777215, 0);
467  auto* pion_tr = tracks.appendNew(Track());
468  pion_tr->setTrackFitResultIndex(Const::pion, 0);
469  pion_tr->addRelationTo(cl1); // a track <--> cluster match
470 
471  // now set mcmatch relations
472  cl0->addRelationTo(true_photon, 12.3);
473  cl0->addRelationTo(true_electron, 2.3);
474  cl1->addRelationTo(true_electron, 45.6);
475  cl1->addRelationTo(true_photon, 5.6);
476  cl1->addRelationTo(true_pion, 15.6);
477 
478  electron_tr->addRelationTo(true_electron);
479  pion_tr->addRelationTo(true_pion);
480 
481  // create belle2::Particles from the mdst objects
482  const auto* photon = particles.appendNew(Particle(cl0));
483  const auto* electron = particles.appendNew(Particle(electron_tr, Const::electron));
484  const auto* pion = particles.appendNew(Particle(pion_tr, Const::pion));
485  const auto* misid_photon = particles.appendNew(Particle(cl1));
486 
487  // now set mcmatch relations
488  photon->addRelationTo(true_photon);
489  electron->addRelationTo(true_electron);
490  pion->addRelationTo(true_pion);
491  misid_photon->addRelationTo(true_electron); // assume MC matching caught this
492  }
493 
494  virtual void TearDown()
495  {
496  DataStore::Instance().reset();
497  }
498  };
499 
500  TEST_F(MCTruthVariablesTest, ECLMCMatchWeightVariable)
501  {
502  StoreArray<Particle> particles;
503  const auto* photon = particles[0];
504  const auto* electron = particles[1];
505  const auto* pion = particles[2];
506 
507  const auto* weight = Manager::Instance().getVariable("clusterMCMatchWeight");
508  EXPECT_FLOAT_EQ(weight->function(photon), 12.3);
509  EXPECT_FLOAT_EQ(weight->function(electron), 45.6);
510  EXPECT_FLOAT_EQ(weight->function(pion), 15.6);
511  }
512 
513  TEST_F(MCTruthVariablesTest, ECLBestMCMatchVariables)
514  {
515  StoreArray<Particle> particles;
516  const auto* photon = particles[0];
517  const auto* electron = particles[1];
518  const auto* pion = particles[2];
519  const auto* misid_photon = particles[3];
520 
521 
522  const auto* pdgcode = Manager::Instance().getVariable("clusterBestMCPDG");
523  EXPECT_EQ(pdgcode->function(photon), 22);
524  EXPECT_EQ(pdgcode->function(electron), 11);
525  EXPECT_EQ(pdgcode->function(pion), 11);
526  EXPECT_EQ(pdgcode->function(misid_photon), 11);
527 
528  const auto* weight = Manager::Instance().getVariable("clusterBestMCMatchWeight");
529  EXPECT_FLOAT_EQ(weight->function(photon), 12.3);
530  EXPECT_FLOAT_EQ(weight->function(electron), 45.6);
531  EXPECT_FLOAT_EQ(weight->function(pion), 45.6);
532  EXPECT_FLOAT_EQ(weight->function(misid_photon), 45.6);
533  }
534 
535  class ROEVariablesTest : public ::testing::Test {
536  protected:
538  void SetUp() override
539  {
540 
541  StoreObjPtr<ParticleList> pi0ParticleList("pi0:vartest");
542  DataStore::Instance().setInitializeActive(true);
543  pi0ParticleList.registerInDataStore(DataStore::c_DontWriteOut);
544  StoreArray<ECLCluster> myECLClusters;
545  StoreArray<KLMCluster> myKLMClusters;
546  StoreArray<TrackFitResult> myTFRs;
547  StoreArray<Track> myTracks;
548  StoreArray<Particle> myParticles;
549  StoreArray<RestOfEvent> myROEs;
550  StoreArray<PIDLikelihood> myPIDLikelihoods;
551  myECLClusters.registerInDataStore();
552  myKLMClusters.registerInDataStore();
553  myTFRs.registerInDataStore();
554  myTracks.registerInDataStore();
555  myParticles.registerInDataStore();
556  myROEs.registerInDataStore();
557  myPIDLikelihoods.registerInDataStore();
558  myParticles.registerRelationTo(myROEs);
559  myTracks.registerRelationTo(myPIDLikelihoods);
560  DataStore::Instance().setInitializeActive(false);
561  }
562 
564  void TearDown() override
565  {
566  DataStore::Instance().reset();
567  }
568  };
569 //
570 // TODO: redo all ROE variable tests
571 //
572 
573  TEST_F(ROEVariablesTest, Variable)
574  {
575  Gearbox& gearbox = Gearbox::getInstance();
576  gearbox.setBackends({std::string("file:")});
577  gearbox.close();
578  gearbox.open("geometry/Belle2.xml", false);
579  StoreObjPtr<ParticleList> pi0ParticleList("pi0:vartest");
580  StoreArray<ECLCluster> myECLClusters;
581  StoreArray<KLMCluster> myKLMClusters;
582  StoreArray<TrackFitResult> myTFRs;
583  StoreArray<Track> myTracks;
584  StoreArray<Particle> myParticles;
585  StoreArray<RestOfEvent> myROEs;
586  StoreArray<PIDLikelihood> myPIDLikelihoods;
587 
588  pi0ParticleList.create();
589  pi0ParticleList->initialize(111, "pi0:vartest");
590 
591  // Neutral ECLCluster on reconstructed side
592  ECLCluster myECL;
593  myECL.setIsTrack(false);
594  float eclREC = 0.5;
595  myECL.setEnergy(eclREC);
596  myECL.setHypothesis(ECLCluster::EHypothesisBit::c_nPhotons);
597  ECLCluster* savedECL = myECLClusters.appendNew(myECL);
598 
599  // Particle on reconstructed side from ECLCluster
600  Particle p(savedECL);
601  Particle* part = myParticles.appendNew(p);
602 
603  // Create ECLCluster on ROE side
604  ECLCluster myROEECL;
605  myROEECL.setIsTrack(false);
606  float eclROE = 1.0;
607  myROEECL.setEnergy(eclROE);
608  myROEECL.setHypothesis(ECLCluster::EHypothesisBit::c_nPhotons);
609  ECLCluster* savedROEECL = myECLClusters.appendNew(myROEECL);
610  Particle* roeECLParticle = myParticles.appendNew(savedROEECL);
611  // Create KLMCluster on ROE side
612  KLMCluster myROEKLM;
613  KLMCluster* savedROEKLM = myKLMClusters.appendNew(myROEKLM);
614  Particle* roeKLMParticle = myParticles.appendNew(savedROEKLM);
615 
616  // Create Track on ROE side
617  // - create TFR
618 
619  const float pValue = 0.5;
620  const float bField = 1.5;
621  const int charge = 1;
622  TMatrixDSym cov6(6);
623 
624  TVector3 position(1.0, 0, 0);
625  TVector3 momentum(0, 1.0, 0);
626 
627  auto CDCValue = static_cast<unsigned long long int>(0x300000000000000);
628 
629  myTFRs.appendNew(position, momentum, cov6, charge, Const::muon, pValue, bField, CDCValue, 16777215, 0);
630 
631  // - create Track
632  Track myROETrack;
633  myROETrack.setTrackFitResultIndex(Const::muon, 0);
634  Track* savedROETrack = myTracks.appendNew(myROETrack);
635  // - create PID information, add relation
636  PIDLikelihood myPID;
637  myPID.setLogLikelihood(Const::TOP, Const::muon, 0.15);
638  myPID.setLogLikelihood(Const::ARICH, Const::muon, 0.152);
639  myPID.setLogLikelihood(Const::ECL, Const::muon, 0.154);
640  myPID.setLogLikelihood(Const::CDC, Const::muon, 0.156);
641  myPID.setLogLikelihood(Const::SVD, Const::muon, 0.158);
642  myPID.setLogLikelihood(Const::TOP, Const::pion, 0.5);
643  myPID.setLogLikelihood(Const::ARICH, Const::pion, 0.52);
644  myPID.setLogLikelihood(Const::ECL, Const::pion, 0.54);
645  myPID.setLogLikelihood(Const::CDC, Const::pion, 0.56);
646  myPID.setLogLikelihood(Const::SVD, Const::pion, 0.58);
647  PIDLikelihood* savedPID = myPIDLikelihoods.appendNew(myPID);
648 
649  savedROETrack->addRelationTo(savedPID);
650  Particle* roeTrackParticle = myParticles.appendNew(savedROETrack, Const::muon);
651 
652  // Create ROE object, append tracks, clusters, add relation to particle
653  //TODO: make particles
654  RestOfEvent roe;
655  vector<const Particle*> roeParticlesToAdd;
656  roeParticlesToAdd.push_back(roeTrackParticle);
657  roeParticlesToAdd.push_back(roeECLParticle);
658  roeParticlesToAdd.push_back(roeKLMParticle);
659  roe.addParticles(roeParticlesToAdd);
660  RestOfEvent* savedROE = myROEs.appendNew(roe);
661  /*
662  std::map<std::string, std::map<unsigned int, bool>> tMasks;
663  std::map<std::string, std::map<unsigned int, bool>> cMasks;
664  std::map<std::string, std::vector<double>> fracs;
665 
666  std::map<unsigned int, bool> tMask1;
667  std::map<unsigned int, bool> tMask2;
668  tMask1[savedROETrack->getArrayIndex()] = true;
669  tMask2[savedROETrack->getArrayIndex()] = false;
670 
671  std::map<unsigned int, bool> cMask1;
672  std::map<unsigned int, bool> cMask2;
673  cMask1[savedROEECL->getArrayIndex()] = true;
674  cMask2[savedROEECL->getArrayIndex()] = false;
675 
676  std::vector<double> frac1 = {0, 0, 1, 0, 0, 0};
677  std::vector<double> frac2 = {1, 1, 1, 1, 1, 1};
678 
679  tMasks["mask1"] = tMask1;
680  tMasks["mask2"] = tMask2;
681 
682  cMasks["mask1"] = cMask1;
683  cMasks["mask2"] = cMask2;
684 
685  fracs["mask1"] = frac1;
686  fracs["mask2"] = frac2;
687 
688  savedROE->appendTrackMasks(tMasks);
689  savedROE->appendECLClusterMasks(cMasks);
690  savedROE->appendChargedStableFractionsSet(fracs);
691  */
692  savedROE->initializeMask("mask1", "test");
693  std::shared_ptr<Variable::Cut> trackSelection = std::shared_ptr<Variable::Cut>(Variable::Cut::compile("p > 2"));
694  std::shared_ptr<Variable::Cut> eclSelection = std::shared_ptr<Variable::Cut>(Variable::Cut::compile("p > 2"));
695  savedROE->updateMaskWithCuts("mask1");
696  savedROE->initializeMask("mask2", "test");
697  savedROE->updateMaskWithCuts("mask2", trackSelection, eclSelection);
698  part->addRelationTo(savedROE);
699 
700  // ROE variables
701  PCmsLabTransform T;
702  float E0 = T.getCMSEnergy() / 2;
703  B2INFO("E0 is " << E0);
704  //*/
705  TLorentzVector pTrack_ROE_Lab(momentum, TMath::Sqrt(Const::muon.getMass()*Const::muon.getMass() + 1.0 /*momentum.Mag2()*/));
706  pTrack_ROE_Lab = roeTrackParticle->get4Vector();
707  TLorentzVector pECL_ROE_Lab(0, 0, eclROE, eclROE);
708  TLorentzVector pECL_REC_Lab(0, 0, eclREC, eclREC);
709 
710  TLorentzVector rec4vec;
711  rec4vec.SetE(pECL_REC_Lab.E());
712  rec4vec.SetVect(pECL_REC_Lab.Vect());
713 
714  TLorentzVector roe4vec;
715  roe4vec.SetE(pTrack_ROE_Lab.E() + pECL_ROE_Lab.E());
716  roe4vec.SetVect(pTrack_ROE_Lab.Vect() + pECL_ROE_Lab.Vect());
717 
718  TLorentzVector rec4vecCMS = T.rotateLabToCms() * rec4vec;
719  TLorentzVector roe4vecCMS = T.rotateLabToCms() * roe4vec;
720 
721  TVector3 pB = - roe4vecCMS.Vect();
722  pB.SetMag(0.340);
723 
724  TLorentzVector m4v0;
725  m4v0.SetE(2 * E0 - (rec4vecCMS.E() + roe4vecCMS.E()));
726  m4v0.SetVect(- (rec4vecCMS.Vect() + roe4vecCMS.Vect()));
727 
728  TLorentzVector m4v1;
729  m4v1.SetE(E0 - rec4vecCMS.E());
730  m4v1.SetVect(- (rec4vecCMS.Vect() + roe4vecCMS.Vect()));
731 
732  TLorentzVector m4v2;
733  m4v2.SetE(E0 - rec4vecCMS.E());
734  m4v2.SetVect(- rec4vecCMS.Vect());
735 
736  TLorentzVector m4v3;
737  m4v3.SetE(E0 - rec4vecCMS.E());
738  m4v3.SetVect(pB - rec4vecCMS.Vect());
739 
740  TLorentzVector neutrino4vecCMS;
741  neutrino4vecCMS.SetVect(- (roe4vecCMS.Vect() + rec4vecCMS.Vect()));
742  neutrino4vecCMS.SetE(neutrino4vecCMS.Vect().Mag());
743 
744  TLorentzVector corrRec4vecCMS = rec4vecCMS + neutrino4vecCMS;
745  B2INFO("roe4vecCMS.E() = " << roe4vecCMS.E());
746  // TESTS FOR ROE STRUCTURE
747  //EXPECT_B2FATAL(savedROE->getTrackMask("noSuchMask"));
748  //EXPECT_B2FATAL(savedROE->getECLClusterMask("noSuchMask"));
749  //double fArray[6];
750  //EXPECT_B2FATAL(savedROE->fillFractions(fArray, "noSuchMask"));
751  EXPECT_B2FATAL(savedROE->updateMaskWithCuts("noSuchMask"));
752  EXPECT_B2FATAL(savedROE->updateMaskWithV0("noSuchMask", part));
753  EXPECT_B2FATAL(savedROE->hasParticle(part, "noSuchMask"));
754 
755  // TESTS FOR ROE VARIABLES
756 
757  const Manager::Var* var = Manager::Instance().getVariable("nROE_Charged(mask1)");
758  ASSERT_NE(var, nullptr);
759  EXPECT_FLOAT_EQ(var->function(part), 1.0);
760 
761  var = Manager::Instance().getVariable("nROE_Charged(mask2)");
762  ASSERT_NE(var, nullptr);
763  EXPECT_FLOAT_EQ(var->function(part), 0.0);
764 
765  var = Manager::Instance().getVariable("nROE_Charged(mask1, 13)");
766  ASSERT_NE(var, nullptr);
767  EXPECT_FLOAT_EQ(var->function(part), 1.0);
768 
769  var = Manager::Instance().getVariable("nROE_Charged(mask1, 211)");
770  ASSERT_NE(var, nullptr);
771  EXPECT_FLOAT_EQ(var->function(part), 0.0);
772 
773  var = Manager::Instance().getVariable("nROE_Photons(mask1)");
774  ASSERT_NE(var, nullptr);
775  EXPECT_FLOAT_EQ(var->function(part), 1.0);
776 
777  var = Manager::Instance().getVariable("nROE_Photons(mask2)");
778  ASSERT_NE(var, nullptr);
779  EXPECT_FLOAT_EQ(var->function(part), 0.0);
780 
781  var = Manager::Instance().getVariable("nROE_NeutralHadrons(mask1)");
782  ASSERT_NE(var, nullptr);
783  EXPECT_FLOAT_EQ(var->function(part), 1.0);
784 
785  var = Manager::Instance().getVariable("nROE_Tracks(mask1)");
786  ASSERT_NE(var, nullptr);
787  EXPECT_FLOAT_EQ(var->function(part), 1.0);
788 
789  var = Manager::Instance().getVariable("nROE_Tracks(mask2)");
790  ASSERT_NE(var, nullptr);
791  EXPECT_FLOAT_EQ(var->function(part), 0.0);
792 
793  var = Manager::Instance().getVariable("nROE_ECLClusters(mask1)");
794  ASSERT_NE(var, nullptr);
795  EXPECT_FLOAT_EQ(var->function(part), 1.0);
796 
797  var = Manager::Instance().getVariable("nROE_ECLClusters(mask2)");
798  ASSERT_NE(var, nullptr);
799  EXPECT_FLOAT_EQ(var->function(part), 0.0);
800 
801  var = Manager::Instance().getVariable("nROE_NeutralECLClusters(mask1)");
802  ASSERT_NE(var, nullptr);
803  EXPECT_FLOAT_EQ(var->function(part), 1.0);
804 
805  var = Manager::Instance().getVariable("nROE_NeutralECLClusters(mask2)");
806  ASSERT_NE(var, nullptr);
807  EXPECT_FLOAT_EQ(var->function(part), 0.0);
808 
809  var = Manager::Instance().getVariable("nROE_ParticlesInList(pi0:vartest)");
810  ASSERT_NE(var, nullptr);
811  EXPECT_FLOAT_EQ(var->function(part), 0.0);
812 
813  var = Manager::Instance().getVariable("roeCharge(mask1)");
814  ASSERT_NE(var, nullptr);
815  EXPECT_FLOAT_EQ(var->function(part), 1.0);
816 
817  var = Manager::Instance().getVariable("roeCharge(mask2)");
818  ASSERT_NE(var, nullptr);
819  EXPECT_FLOAT_EQ(var->function(part), 0.0);
820 
821  var = Manager::Instance().getVariable("roeEextra(mask1)");
822  ASSERT_NE(var, nullptr);
823  EXPECT_FLOAT_EQ(var->function(part), savedROEECL->getEnergy(ECLCluster::EHypothesisBit::c_nPhotons));
824 
825  var = Manager::Instance().getVariable("roeEextra(mask2)");
826  ASSERT_NE(var, nullptr);
827  EXPECT_FLOAT_EQ(var->function(part), 0.0);
828 
829  var = Manager::Instance().getVariable("roeDeltae(mask1)");
830  ASSERT_NE(var, nullptr);
831  EXPECT_FLOAT_EQ(var->function(part), roe4vecCMS.E() - E0);
832 
833  var = Manager::Instance().getVariable("roeDeltae(mask2)");
834  ASSERT_NE(var, nullptr);
835  EXPECT_FLOAT_EQ(var->function(part), -E0);
836 
837  var = Manager::Instance().getVariable("roeMbc(mask1)");
838  ASSERT_NE(var, nullptr);
839  EXPECT_FLOAT_EQ(var->function(part), TMath::Sqrt(E0 * E0 - roe4vecCMS.Vect().Mag2()));
840 
841  var = Manager::Instance().getVariable("roeMbc(mask2)");
842  ASSERT_NE(var, nullptr);
843  EXPECT_FLOAT_EQ(var->function(part), E0);
844 
845  var = Manager::Instance().getVariable("weDeltae(mask1,0)");
846  ASSERT_NE(var, nullptr);
847  EXPECT_FLOAT_EQ(var->function(part), corrRec4vecCMS.E() - E0);
848 
849  var = Manager::Instance().getVariable("weDeltae(mask2,0)");
850  ASSERT_NE(var, nullptr);
851  EXPECT_FLOAT_EQ(var->function(part), rec4vecCMS.E() + rec4vecCMS.Vect().Mag() - E0);
852 
853  var = Manager::Instance().getVariable("weMbc(mask1,0)");
854  ASSERT_NE(var, nullptr);
855  EXPECT_FLOAT_EQ(var->function(part), TMath::Sqrt(E0 * E0 - corrRec4vecCMS.Vect().Mag2()));
856 
857  var = Manager::Instance().getVariable("weMbc(mask2,0)");
858  ASSERT_NE(var, nullptr);
859  EXPECT_FLOAT_EQ(var->function(part), E0);
860 
861  var = Manager::Instance().getVariable("weMissM2(mask1,0)");
862  ASSERT_NE(var, nullptr);
863  EXPECT_FLOAT_EQ(var->function(part), m4v0.Mag2());
864 
865  var = Manager::Instance().getVariable("weMissM2(mask2,0)");
866  ASSERT_NE(var, nullptr);
867  EXPECT_FLOAT_EQ(var->function(part), (2 * E0 - rec4vecCMS.E()) * (2 * E0 - rec4vecCMS.E()) - rec4vecCMS.Vect().Mag2());
868 
869  }
870 
871 
872  class EventVariableTest : public ::testing::Test {
873  protected:
875  void SetUp() override
876  {
877  DataStore::Instance().setInitializeActive(true);
878  StoreArray<Particle>().registerInDataStore();
879  StoreArray<MCParticle>().registerInDataStore();
880  DataStore::Instance().setInitializeActive(false);
881 
882  }
883 
885  void TearDown() override
886  {
887  DataStore::Instance().reset();
888  }
889  };
890 
891  TEST_F(EventVariableTest, ExperimentRunEventDateAndTime)
892  {
893  const Manager::Var* exp = Manager::Instance().getVariable("expNum");
894  const Manager::Var* run = Manager::Instance().getVariable("runNum");
895  const Manager::Var* evt = Manager::Instance().getVariable("evtNum");
896  const Manager::Var* date = Manager::Instance().getVariable("date");
897  const Manager::Var* year = Manager::Instance().getVariable("year");
898  const Manager::Var* time = Manager::Instance().getVariable("eventTimeSeconds");
899 
900  // there is no EventMetaData so expect nan
901  EXPECT_FALSE(date->function(nullptr) == date->function(nullptr));
902  EXPECT_FALSE(year->function(nullptr) == year->function(nullptr));
903  EXPECT_FALSE(time->function(nullptr) == time->function(nullptr));
904 
905  DataStore::Instance().setInitializeActive(true);
906  StoreObjPtr<EventMetaData> evtMetaData;
907  evtMetaData.registerInDataStore();
908  DataStore::Instance().setInitializeActive(false);
909  evtMetaData.create();
910  evtMetaData->setExperiment(1337);
911  evtMetaData->setRun(12345);
912  evtMetaData->setEvent(54321);
913  evtMetaData->setTime(1288569600e9);
914  // 01/11/2010 is the date TDR was uploaded to arXiv ... experiment's birthday?
915 
916 
917  // -
918  EXPECT_FLOAT_EQ(exp->function(NULL), 1337.);
919  EXPECT_FLOAT_EQ(run->function(NULL), 12345.);
920  EXPECT_FLOAT_EQ(evt->function(NULL), 54321.);
921  EXPECT_FLOAT_EQ(date->function(NULL), 20101101.);
922  EXPECT_FLOAT_EQ(year->function(NULL), 2010.);
923  EXPECT_FLOAT_EQ(time->function(NULL), 1288569600);
924  }
925 
926  TEST_F(EventVariableTest, TestGlobalCounters)
927  {
928  StoreArray<MCParticle> mcParticles; // empty
929  const Manager::Var* var = Manager::Instance().getVariable("nMCParticles");
930  EXPECT_FLOAT_EQ(var->function(NULL), 0.0);
931 
932  for (unsigned i = 0; i < 10; ++i)
933  mcParticles.appendNew();
934 
935  EXPECT_FLOAT_EQ(var->function(NULL), 10.0);
936 
937  // TODO: add other counters nTracks etc in here
938  }
939 
940  TEST_F(EventVariableTest, TestIfContinuumEvent_ForContinuumEvent)
941  {
942  DataStore::Instance().setInitializeActive(true);
943  StoreArray<MCParticle> mcParticles;
944  StoreArray<Particle> particles;
945  particles.registerRelationTo(mcParticles);
946  DataStore::Instance().setInitializeActive(false);
947 
948  auto* mcParticle = mcParticles.appendNew();
949  mcParticle->setPDG(11);
950  mcParticle->setStatus(MCParticle::c_PrimaryParticle);
951  auto* p1 = particles.appendNew(TLorentzVector({ 0.0 , -0.4, 0.8, 1.0}), 11);
952  p1->addRelationTo(mcParticle);
953 
954  mcParticle = mcParticles.appendNew();
955  mcParticle->setPDG(-11);
956  mcParticle->setStatus(MCParticle::c_PrimaryParticle);
957  auto* p2 = particles.appendNew(TLorentzVector({ 0.0 , -0.4, 0.8, 1.0}), 11);
958  p2->addRelationTo(mcParticle);
959 
960  const Manager::Var* var = Manager::Instance().getVariable("isContinuumEvent");
961  ASSERT_NE(var, nullptr);
962  EXPECT_FLOAT_EQ(var->function(p1), 1.0);
963  EXPECT_FLOAT_EQ(var->function(p2), 1.0);
964  const Manager::Var* varN = Manager::Instance().getVariable("isNotContinuumEvent");
965  ASSERT_NE(varN, nullptr);
966  EXPECT_FLOAT_EQ(varN->function(p1), 0.0);
967  EXPECT_FLOAT_EQ(varN->function(p2), 0.0);
968  }
969 
970  TEST_F(EventVariableTest, TestIfContinuumEvent_ForUpsilon4SEvent)
971  {
972  DataStore::Instance().setInitializeActive(true);
973  StoreArray<MCParticle> mcParticles2;
974  StoreArray<Particle> particles2;
975  particles2.registerRelationTo(mcParticles2);
976  DataStore::Instance().setInitializeActive(false);
977 
978  auto* mcParticle = mcParticles2.appendNew();
979  mcParticle->setPDG(22);
980  mcParticle->setStatus(MCParticle::c_PrimaryParticle);
981  auto* p3 = particles2.appendNew(TLorentzVector({ 0.0 , -0.4, 0.8, 1.0}), 11);
982  p3->addRelationTo(mcParticle);
983 
984  mcParticle = mcParticles2.appendNew();
985  mcParticle->setPDG(300553);
986  mcParticle->setStatus(MCParticle::c_PrimaryParticle);
987  auto* p4 = particles2.appendNew(TLorentzVector({ 0.0 , -0.4, 0.8, 1.0}), 300553);
988  p4->addRelationTo(mcParticle);
989 
990  const Manager::Var* var2 = Manager::Instance().getVariable("isContinuumEvent");
991  ASSERT_NE(var2, nullptr);
992  EXPECT_FLOAT_EQ(var2->function(p3), 0.0);
993  EXPECT_FLOAT_EQ(var2->function(p4), 0.0);
994  const Manager::Var* var2N = Manager::Instance().getVariable("isNotContinuumEvent");
995  ASSERT_NE(var2N, nullptr);
996  EXPECT_FLOAT_EQ(var2N->function(p3), 1.0);
997  EXPECT_FLOAT_EQ(var2N->function(p4), 1.0);
998  }
999 
1000  TEST_F(EventVariableTest, TestIfContinuumEvent_ForWrongReconstructedUpsilon4SEvent)
1001  {
1002  DataStore::Instance().setInitializeActive(true);
1003  StoreArray<MCParticle> mcParticles3;
1004  StoreArray<Particle> particles3;
1005  particles3.registerRelationTo(mcParticles3);
1006  DataStore::Instance().setInitializeActive(false);
1007 
1008  auto* mcParticle = mcParticles3.appendNew();
1009  mcParticle->setPDG(22);
1010  mcParticle->setStatus(MCParticle::c_PrimaryParticle);
1011  auto* p5 = particles3.appendNew(TLorentzVector({ 0.0 , -0.4, 0.8, 1.0}), 11);
1012  p5->addRelationTo(mcParticle);
1013 
1014  mcParticle = mcParticles3.appendNew();
1015  mcParticle->setPDG(300553);
1016  mcParticle->setStatus(MCParticle::c_PrimaryParticle);
1017  auto* p6 = particles3.appendNew(TLorentzVector({ 0.0 , -0.4, 0.8, 1.0}), 15);
1018  p6->addRelationTo(mcParticle);
1019 
1020  const Manager::Var* var3 = Manager::Instance().getVariable("isContinuumEvent");
1021  ASSERT_NE(var3, nullptr);
1022  EXPECT_FLOAT_EQ(var3->function(p5), 0.0);
1023  EXPECT_FLOAT_EQ(var3->function(p6), 0.0);
1024  const Manager::Var* var3N = Manager::Instance().getVariable("isNotContinuumEvent");
1025  ASSERT_NE(var3N, nullptr);
1026  EXPECT_FLOAT_EQ(var3N->function(p5), 1.0);
1027  EXPECT_FLOAT_EQ(var3N->function(p6), 1.0);
1028  }
1029 
1030 
1031  class MetaVariableTest : public ::testing::Test {
1032  protected:
1034  void SetUp() override
1035  {
1036  DataStore::Instance().setInitializeActive(true);
1037  StoreObjPtr<ParticleExtraInfoMap>().registerInDataStore();
1038  StoreObjPtr<EventExtraInfo>().registerInDataStore();
1039  StoreArray<Particle>().registerInDataStore();
1040  StoreArray<MCParticle>().registerInDataStore();
1041  DataStore::Instance().setInitializeActive(false);
1042  }
1043 
1045  void TearDown() override
1046  {
1047  DataStore::Instance().reset();
1048  }
1049  };
1050 
1051  TEST_F(MetaVariableTest, countDaughters)
1052  {
1053  TLorentzVector momentum;
1054  const int nDaughters = 6;
1055  StoreArray<Particle> particles;
1056  std::vector<int> daughterIndices;
1057  for (int i = 0; i < nDaughters; i++) {
1058  Particle d(TLorentzVector(1, 0, 0, 3.0), (i % 2) ? 211 : -211);
1059  momentum += d.get4Vector();
1060  Particle* newDaughters = particles.appendNew(d);
1061  daughterIndices.push_back(newDaughters->getArrayIndex());
1062  }
1063  const Particle* p = particles.appendNew(momentum, 411, Particle::c_Unflavored, daughterIndices);
1064 
1065  const Manager::Var* var = Manager::Instance().getVariable("countDaughters(charge > 0)");
1066  ASSERT_NE(var, nullptr);
1067  EXPECT_DOUBLE_EQ(var->function(p), 3.0);
1068 
1069  var = Manager::Instance().getVariable("countDaughters(abs(charge) > 0)");
1070  ASSERT_NE(var, nullptr);
1071  EXPECT_DOUBLE_EQ(var->function(p), 6.0);
1072 
1073  }
1074 
1075  TEST_F(MetaVariableTest, useRestFrame)
1076  {
1077  Gearbox& gearbox = Gearbox::getInstance();
1078  gearbox.setBackends({std::string("file:")});
1079  gearbox.close();
1080  gearbox.open("geometry/Belle2.xml", false);
1081 
1082  Particle p({ 0.1 , -0.4, 0.8, 1.0 }, 11);
1083  p.setVertex(TVector3(1.0, 2.0, 2.0));
1084 
1085  const Manager::Var* var = Manager::Instance().getVariable("p");
1086  ASSERT_NE(var, nullptr);
1087  EXPECT_FLOAT_EQ(var->function(&p), 0.9);
1088 
1089  var = Manager::Instance().getVariable("E");
1090  ASSERT_NE(var, nullptr);
1091  EXPECT_FLOAT_EQ(var->function(&p), 1.0);
1092 
1093  var = Manager::Instance().getVariable("distance");
1094  ASSERT_NE(var, nullptr);
1095  EXPECT_FLOAT_EQ(var->function(&p), 3.0);
1096 
1097  var = Manager::Instance().getVariable("useRestFrame(p)");
1098  ASSERT_NE(var, nullptr);
1099  EXPECT_ALL_NEAR(var->function(&p), 0.0, 1e-9);
1100 
1101  var = Manager::Instance().getVariable("useRestFrame(E)");
1102  ASSERT_NE(var, nullptr);
1103  EXPECT_FLOAT_EQ(var->function(&p), 0.4358899);
1104 
1105  var = Manager::Instance().getVariable("useRestFrame(distance)");
1106  ASSERT_NE(var, nullptr);
1107  EXPECT_FLOAT_EQ(var->function(&p), 0.0);
1108  }
1109 
1110  TEST_F(MetaVariableTest, useLabFrame)
1111  {
1112  Particle p({ 0.1 , -0.4, 0.8, 1.0 }, 11);
1113  p.setVertex(TVector3(1.0, 2.0, 2.0));
1114 
1115  const Manager::Var* var = Manager::Instance().getVariable("p");
1116  ASSERT_NE(var, nullptr);
1117  EXPECT_FLOAT_EQ(var->function(&p), 0.9);
1118 
1119  var = Manager::Instance().getVariable("E");
1120  ASSERT_NE(var, nullptr);
1121  EXPECT_FLOAT_EQ(var->function(&p), 1.0);
1122 
1123  var = Manager::Instance().getVariable("distance");
1124  ASSERT_NE(var, nullptr);
1125  EXPECT_FLOAT_EQ(var->function(&p), 3.0);
1126 
1127  var = Manager::Instance().getVariable("useLabFrame(p)");
1128  ASSERT_NE(var, nullptr);
1129  EXPECT_FLOAT_EQ(var->function(&p), 0.9);
1130 
1131  var = Manager::Instance().getVariable("useLabFrame(E)");
1132  ASSERT_NE(var, nullptr);
1133  EXPECT_FLOAT_EQ(var->function(&p), 1.0);
1134 
1135  var = Manager::Instance().getVariable("useLabFrame(distance)");
1136  ASSERT_NE(var, nullptr);
1137  EXPECT_FLOAT_EQ(var->function(&p), 3.0);
1138  }
1139 
1140  TEST_F(MetaVariableTest, useCMSFrame)
1141  {
1142  Gearbox& gearbox = Gearbox::getInstance();
1143  gearbox.setBackends({std::string("file:")});
1144  gearbox.close();
1145  gearbox.open("geometry/Belle2.xml", false);
1146 
1147  Particle p({ 0.1 , -0.4, 0.8, 1.0 }, 11);
1148  p.setVertex(TVector3(1.0, 2.0, 2.0));
1149 
1150  const Manager::Var* var = Manager::Instance().getVariable("p");
1151  ASSERT_NE(var, nullptr);
1152  EXPECT_FLOAT_EQ(var->function(&p), 0.9);
1153 
1154  var = Manager::Instance().getVariable("E");
1155  ASSERT_NE(var, nullptr);
1156  EXPECT_FLOAT_EQ(var->function(&p), 1.0);
1157 
1158  var = Manager::Instance().getVariable("distance");
1159  ASSERT_NE(var, nullptr);
1160  EXPECT_FLOAT_EQ(var->function(&p), 3.0);
1161 
1162  var = Manager::Instance().getVariable("useCMSFrame(p)");
1163  ASSERT_NE(var, nullptr);
1164  EXPECT_FLOAT_EQ(var->function(&p), 0.68176979);
1165 
1166  var = Manager::Instance().getVariable("useCMSFrame(E)");
1167  ASSERT_NE(var, nullptr);
1168  EXPECT_FLOAT_EQ(var->function(&p), 0.80920333);
1169 
1170  var = Manager::Instance().getVariable("useCMSFrame(distance)");
1171  ASSERT_NE(var, nullptr);
1172  EXPECT_FLOAT_EQ(var->function(&p), 3.185117);
1173  }
1174 
1175  TEST_F(MetaVariableTest, useTagSideRecoilRestFrame)
1176  {
1177  DataStore::Instance().setInitializeActive(true);
1178  StoreArray<Particle> particles;
1179  particles.registerInDataStore();
1180  DataStore::Instance().setInitializeActive(false);
1181  PCmsLabTransform T;
1182  TLorentzVector vec0 = {0.0, 0.0, 0.0, T.getCMSEnergy()};
1183  TLorentzVector vec1 = {0.0, +0.332174566, 0.0, T.getCMSEnergy() / 2.};
1184  TLorentzVector vec2 = {0.0, -0.332174566, 0.0, T.getCMSEnergy() / 2.};
1185  Particle* p0 = particles.appendNew(Particle(T.rotateCmsToLab() * vec0, 300553));
1186  Particle* p1 = particles.appendNew(Particle(T.rotateCmsToLab() * vec1, 511, Particle::c_Unflavored, Particle::c_Undefined, 1));
1187  Particle* p2 = particles.appendNew(Particle(T.rotateCmsToLab() * vec2, -511, Particle::c_Unflavored, Particle::c_Undefined, 2));
1188 
1189  p0->appendDaughter(p1->getArrayIndex());
1190  p0->appendDaughter(p2->getArrayIndex());
1191 
1192  const Manager::Var* var = Manager::Instance().getVariable("useTagSideRecoilRestFrame(daughter(1, p), 0)");
1193  ASSERT_NE(var, nullptr);
1194  EXPECT_NEAR(var->function(p0), 0., 1e-6);
1195 
1196  var = Manager::Instance().getVariable("useTagSideRecoilRestFrame(daughter(1, px), 0)");
1197  ASSERT_NE(var, nullptr);
1198  EXPECT_NEAR(var->function(p0), 0., 1e-6);
1199 
1200  var = Manager::Instance().getVariable("useTagSideRecoilRestFrame(daughter(1, py), 0)");
1201  ASSERT_NE(var, nullptr);
1202  EXPECT_NEAR(var->function(p0), 0., 1e-6);
1203 
1204  var = Manager::Instance().getVariable("useTagSideRecoilRestFrame(daughter(1, pz), 0)");
1205  ASSERT_NE(var, nullptr);
1206  EXPECT_NEAR(var->function(p0), 0., 1e-6);
1207 
1208  var = Manager::Instance().getVariable("useTagSideRecoilRestFrame(daughter(1, E), 0)");
1209  ASSERT_NE(var, nullptr);
1210  EXPECT_NEAR(var->function(p0), p1->getMass(), 1e-6);
1211  }
1212 
1213 
1214  TEST_F(MetaVariableTest, extraInfo)
1215  {
1216  Particle p({ 0.1 , -0.4, 0.8, 1.0 }, 11);
1217  p.addExtraInfo("pi", 3.14);
1218 
1219  const Manager::Var* var = Manager::Instance().getVariable("extraInfo(pi)");
1220  ASSERT_NE(var, nullptr);
1221  EXPECT_FLOAT_EQ(var->function(&p), 3.14);
1222 
1223  // If nullptr is given, NaN is returned
1224  EXPECT_TRUE(std::isnan(var->function(nullptr)));
1225  }
1226 
1227  TEST_F(MetaVariableTest, eventExtraInfo)
1228  {
1229  StoreObjPtr<EventExtraInfo> eventExtraInfo;
1230  if (not eventExtraInfo.isValid())
1231  eventExtraInfo.create();
1232  eventExtraInfo->addExtraInfo("pi", 3.14);
1233  const Manager::Var* var = Manager::Instance().getVariable("eventExtraInfo(pi)");
1234  ASSERT_NE(var, nullptr);
1235  EXPECT_FLOAT_EQ(var->function(nullptr), 3.14);
1236  }
1237 
1238  TEST_F(MetaVariableTest, eventCached)
1239  {
1240  const Manager::Var* var = Manager::Instance().getVariable("eventCached(constant(3.14))");
1241  ASSERT_NE(var, nullptr);
1242  EXPECT_FLOAT_EQ(var->function(nullptr), 3.14);
1243  StoreObjPtr<EventExtraInfo> eventExtraInfo;
1244  EXPECT_TRUE(eventExtraInfo.isValid());
1245  EXPECT_TRUE(eventExtraInfo->hasExtraInfo("__constant__bo3__pt14__bc"));
1246  EXPECT_FLOAT_EQ(eventExtraInfo->getExtraInfo("__constant__bo3__pt14__bc"), 3.14);
1247  eventExtraInfo->addExtraInfo("__eventExtraInfo__bopi__bc", 3.14);
1248  var = Manager::Instance().getVariable("eventCached(eventExtraInfo(pi))");
1249  ASSERT_NE(var, nullptr);
1250  EXPECT_FLOAT_EQ(var->function(nullptr), 3.14);
1251  }
1252 
1253  TEST_F(MetaVariableTest, particleCached)
1254  {
1255  Particle p({ 0.1 , -0.4, 0.8, 2.0 }, 11);
1256  const Manager::Var* var = Manager::Instance().getVariable("particleCached(px)");
1257  ASSERT_NE(var, nullptr);
1258  EXPECT_FLOAT_EQ(var->function(&p), 0.1);
1259  EXPECT_TRUE(p.hasExtraInfo("__px"));
1260  EXPECT_FLOAT_EQ(p.getExtraInfo("__px"), 0.1);
1261  p.addExtraInfo("__py", -0.5); // NOT -0.4 because we want to see if the cache is used instead of py!
1262  var = Manager::Instance().getVariable("particleCached(py)");
1263  ASSERT_NE(var, nullptr);
1264  EXPECT_FLOAT_EQ(var->function(&p), -0.5);
1265  }
1266 
1267  TEST_F(MetaVariableTest, basicMathTest)
1268  {
1269  Particle p({ 0.1 , -0.4, 0.8, 2.0 }, 11);
1270 
1271  const Manager::Var* var = Manager::Instance().getVariable("abs(py)");
1272  ASSERT_NE(var, nullptr);
1273  EXPECT_FLOAT_EQ(var->function(&p), 0.4);
1274 
1275  var = Manager::Instance().getVariable("min(E, pz)");
1276  ASSERT_NE(var, nullptr);
1277  EXPECT_FLOAT_EQ(var->function(&p), 0.8);
1278 
1279  var = Manager::Instance().getVariable("max(E, pz)");
1280  ASSERT_NE(var, nullptr);
1281  EXPECT_FLOAT_EQ(var->function(&p), 2.0);
1282 
1283  var = Manager::Instance().getVariable("log10(px)");
1284  ASSERT_NE(var, nullptr);
1285  EXPECT_FLOAT_EQ(var->function(&p), -1.0);
1286  }
1287 
1288  TEST_F(MetaVariableTest, formula)
1289  {
1290  // see also unit tests in framework/formula_parser.cc
1291  //
1292  // keep particle-based tests here, and operator precidence tests (etc) in
1293  // framework with the parser itself
1294 
1295  Particle p({ 0.1 , -0.4, 0.8, 2.0 }, 11);
1296 
1297  const Manager::Var* var = Manager::Instance().getVariable("formula(px + py)");
1298  ASSERT_NE(var, nullptr);
1299  EXPECT_FLOAT_EQ(var->function(&p), -0.3);
1300 
1301  var = Manager::Instance().getVariable("formula(px - py)");
1302  ASSERT_NE(var, nullptr);
1303  EXPECT_FLOAT_EQ(var->function(&p), 0.5);
1304 
1305  var = Manager::Instance().getVariable("formula(px * py)");
1306  ASSERT_NE(var, nullptr);
1307  EXPECT_FLOAT_EQ(var->function(&p), -0.04);
1308 
1309  var = Manager::Instance().getVariable("formula(py / px)");
1310  ASSERT_NE(var, nullptr);
1311  EXPECT_FLOAT_EQ(var->function(&p), -4.0);
1312 
1313  var = Manager::Instance().getVariable("formula(px ^ E)");
1314  ASSERT_NE(var, nullptr);
1315  EXPECT_FLOAT_EQ(var->function(&p), 0.01);
1316 
1317  var = Manager::Instance().getVariable("formula(px * py + pz)");
1318  ASSERT_NE(var, nullptr);
1319  EXPECT_ALL_NEAR(var->function(&p), 0.76, 1e-6);
1320 
1321  var = Manager::Instance().getVariable("formula(pz + px * py)");
1322  ASSERT_NE(var, nullptr);
1323  EXPECT_ALL_NEAR(var->function(&p), 0.76, 1e-6);
1324 
1325  var = Manager::Instance().getVariable("formula(pt)");
1326  ASSERT_NE(var, nullptr);
1327  EXPECT_FLOAT_EQ(var->function(&p), 0.41231057);
1328  double pt = var->function(&p);
1329 
1330  var = Manager::Instance().getVariable("formula((px**2 + py**2)**(1/2))");
1331  ASSERT_NE(var, nullptr);
1332  EXPECT_FLOAT_EQ(var->function(&p), pt);
1333 
1334  var = Manager::Instance().getVariable("formula(charge)");
1335  ASSERT_NE(var, nullptr);
1336  EXPECT_FLOAT_EQ(var->function(&p), -1.0);
1337 
1338  var = Manager::Instance().getVariable("formula(charge**2)");
1339  ASSERT_NE(var, nullptr);
1340  EXPECT_FLOAT_EQ(var->function(&p), 1.0);
1341 
1342  var = Manager::Instance().getVariable("formula(charge^2)");
1343  ASSERT_NE(var, nullptr);
1344  EXPECT_FLOAT_EQ(var->function(&p), 1.0);
1345 
1346  var = Manager::Instance().getVariable("formula(PDG * charge)");
1347  ASSERT_NE(var, nullptr);
1348  EXPECT_FLOAT_EQ(var->function(&p), -11.0);
1349 
1350  var = Manager::Instance().getVariable("formula(PDG**2 * charge)");
1351  ASSERT_NE(var, nullptr);
1352  EXPECT_FLOAT_EQ(var->function(&p), -121.0);
1353 
1354  var = Manager::Instance().getVariable("formula(10.58 - (px + py + pz - E)**2)");
1355  ASSERT_NE(var, nullptr);
1356  EXPECT_FLOAT_EQ(var->function(&p), 8.33);
1357 
1358  var = Manager::Instance().getVariable("formula(-10.58 + (px + py + pz - E)**2)");
1359  ASSERT_NE(var, nullptr);
1360  EXPECT_FLOAT_EQ(var->function(&p), -8.33);
1361 
1362  var = Manager::Instance().getVariable("formula(-1.0 * PDG)");
1363  ASSERT_NE(var, nullptr);
1364  EXPECT_FLOAT_EQ(var->function(&p), -11);
1365  }
1366 
1367  TEST_F(MetaVariableTest, passesCut)
1368  {
1369  Particle p({ 0.1 , -0.4, 0.8, 2.0 }, 11);
1370  Particle p2({ 0.1 , -0.4, 0.8, 4.0 }, 11);
1371 
1372  const Manager::Var* var = Manager::Instance().getVariable("passesCut(E < 3)");
1373  ASSERT_NE(var, nullptr);
1374  EXPECT_FLOAT_EQ(var->function(&p), 1);
1375  EXPECT_FLOAT_EQ(var->function(&p2), 0);
1376  EXPECT_TRUE(std::isnan(var->function(nullptr)));
1377 
1378  }
1379 
1380  TEST_F(MetaVariableTest, unmask)
1381  {
1382  DataStore::Instance().setInitializeActive(true);
1383  StoreArray<MCParticle> mcParticles;
1384  StoreArray<Particle> particles;
1385  particles.registerInDataStore();
1386  mcParticles.registerInDataStore();
1387  particles.registerRelationTo(mcParticles);
1388  DataStore::Instance().setInitializeActive(false);
1389 
1390  // Create MC graph for B -> (muon -> electron + muon_neutrino) + anti_muon_neutrino
1391  MCParticleGraph mcGraph;
1392 
1393  MCParticleGraph::GraphParticle& graphParticleGrandMother = mcGraph.addParticle();
1394 
1395  MCParticleGraph::GraphParticle& graphParticleMother = mcGraph.addParticle();
1396  MCParticleGraph::GraphParticle& graphParticleAunt = mcGraph.addParticle();
1397 
1398  MCParticleGraph::GraphParticle& graphParticleDaughter1 = mcGraph.addParticle();
1399  MCParticleGraph::GraphParticle& graphParticleDaughter2 = mcGraph.addParticle();
1400 
1401  graphParticleGrandMother.setPDG(-521);
1402  graphParticleMother.setPDG(13);
1403  graphParticleAunt.setPDG(-14);
1404  graphParticleDaughter1.setPDG(11);
1405  graphParticleDaughter2.setPDG(14);
1406 
1407  graphParticleMother.comesFrom(graphParticleGrandMother);
1408  graphParticleAunt.comesFrom(graphParticleGrandMother);
1409  graphParticleDaughter1.comesFrom(graphParticleMother);
1410  graphParticleDaughter2.comesFrom(graphParticleMother);
1411  mcGraph.generateList();
1412 
1413 
1414  // Get MC Particles from StoreArray
1415  auto* mcGrandMother = mcParticles[0];
1416  mcGrandMother->setStatus(MCParticle::c_PrimaryParticle);
1417 
1418  auto* mcMother = mcParticles[1];
1419  mcMother->setStatus(MCParticle::c_PrimaryParticle);
1420 
1421  auto* mcAunt = mcParticles[2];
1422  mcAunt->setStatus(MCParticle::c_PrimaryParticle);
1423 
1424  auto* mcDaughter1 = mcParticles[3];
1425  mcDaughter1->setStatus(MCParticle::c_PrimaryParticle);
1426 
1427  auto* mcDaughter2 = mcParticles[4];
1428  mcDaughter2->setStatus(MCParticle::c_PrimaryParticle);
1429 
1430  auto* pGrandMother = particles.appendNew(TLorentzVector({ 0.0 , -0.4, 0.8, 1.0}), -521);
1431  pGrandMother->addRelationTo(mcGrandMother);
1432 
1433  auto* pMother = particles.appendNew(TLorentzVector({ 0.0 , -0.4, 0.8, 1.0}), 13);
1434  pMother->addRelationTo(mcMother);
1435 
1436 
1437  pMother->writeExtraInfo("mcErrors", 8);
1438  pGrandMother->writeExtraInfo("mcErrors", 8 | 16);
1439  const Manager::Var* var1 = Manager::Instance().getVariable("unmask(mcErrors, 8)");
1440  const Manager::Var* var2 = Manager::Instance().getVariable("unmask(mcErrors, 8, 16, 32, 64)");
1441  ASSERT_NE(var1, nullptr);
1442  EXPECT_FLOAT_EQ(var1->function(pMother), 0);
1443  EXPECT_FLOAT_EQ(var1->function(pGrandMother), 16);
1444  ASSERT_NE(var2, nullptr);
1445  EXPECT_FLOAT_EQ(var2->function(pMother), 0);
1446  EXPECT_FLOAT_EQ(var2->function(pGrandMother), 0);
1447 
1448 
1449  pMother->writeExtraInfo("mcErrors", 8 | 128);
1450  pGrandMother->writeExtraInfo("mcErrors", 8 | 16 | 512);
1451  ASSERT_NE(var1, nullptr);
1452  EXPECT_FLOAT_EQ(var1->function(pMother), 128);
1453  EXPECT_FLOAT_EQ(var1->function(pGrandMother), 16 | 512);
1454  ASSERT_NE(var2, nullptr);
1455  EXPECT_FLOAT_EQ(var2->function(pMother), 128);
1456  EXPECT_FLOAT_EQ(var2->function(pGrandMother), 512);
1457 
1458  // unmask variable needs at least two arguments
1459  EXPECT_B2FATAL(Manager::Instance().getVariable("unmask(mcErrors)"));
1460 
1461  // all but the first argument have to be integers
1462  EXPECT_B2FATAL(Manager::Instance().getVariable("unmask(mcErrors, NOTINT)"));
1463  }
1464 
1465  TEST_F(MetaVariableTest, conditionalVariableSelector)
1466  {
1467  Particle p({ 0.1, -0.4, 0.8, 2.0 }, 11);
1468 
1469  const Manager::Var* var = Manager::Instance().getVariable("conditionalVariableSelector(E>1, px, py)");
1470  ASSERT_NE(var, nullptr);
1471  EXPECT_FLOAT_EQ(var->function(&p), 0.1);
1472 
1473  var = Manager::Instance().getVariable("conditionalVariableSelector(E<1, px, py)");
1474  ASSERT_NE(var, nullptr);
1475  EXPECT_FLOAT_EQ(var->function(&p), -0.4);
1476 
1477  }
1478 
1479  TEST_F(MetaVariableTest, nCleanedTracks)
1480  {
1481  DataStore::Instance().setInitializeActive(true);
1482  StoreArray<TrackFitResult> track_fit_results;
1483  StoreArray<Track> tracks;
1484  track_fit_results.registerInDataStore();
1485  tracks.registerInDataStore();
1486  DataStore::Instance().setInitializeActive(false);
1487 
1488  Particle p({ 0.1 , -0.4, 0.8, 2.0 }, 11);
1489  Particle p2({ 0.1 , -0.4, 0.8, 4.0 }, 11);
1490 
1491  track_fit_results.appendNew(TVector3(0.1, 0.1, 0.1), TVector3(0.1, 0.0, 0.0),
1492  TMatrixDSym(6), 1, Const::pion, 0.01, 1.5, 0, 0, 0);
1493  track_fit_results.appendNew(TVector3(0.1, 0.1, 0.1), TVector3(0.15, 0.0, 0.0),
1494  TMatrixDSym(6), 1, Const::pion, 0.01, 1.5, 0, 0, 0);
1495  track_fit_results.appendNew(TVector3(0.1, 0.1, 0.1), TVector3(0.4, 0.0, 0.0),
1496  TMatrixDSym(6), 1, Const::pion, 0.01, 1.5, 0, 0, 0);
1497  track_fit_results.appendNew(TVector3(0.1, 0.1, 0.1), TVector3(0.6, 0.0, 0.0),
1498  TMatrixDSym(6), 1, Const::pion, 0.01, 1.5, 0, 0, 0);
1499 
1500  tracks.appendNew()->setTrackFitResultIndex(Const::pion, 0);
1501  tracks.appendNew()->setTrackFitResultIndex(Const::pion, 1);
1502  tracks.appendNew()->setTrackFitResultIndex(Const::pion, 2);
1503  tracks.appendNew()->setTrackFitResultIndex(Const::pion, 3);
1504 
1505  const Manager::Var* var1 = Manager::Instance().getVariable("nCleanedTracks(p > 0.5)");
1506  EXPECT_FLOAT_EQ(var1->function(nullptr), 1);
1507 
1508  const Manager::Var* var2 = Manager::Instance().getVariable("nCleanedTracks(p > 0.2)");
1509  EXPECT_FLOAT_EQ(var2->function(nullptr), 2);
1510 
1511  const Manager::Var* var3 = Manager::Instance().getVariable("nCleanedTracks()");
1512  EXPECT_FLOAT_EQ(var3->function(nullptr), 4);
1513 
1514 
1515  }
1516 
1517  TEST_F(MetaVariableTest, NumberOfMCParticlesInEvent)
1518  {
1519  Particle p({ 0.1 , -0.4, 0.8, 2.0 }, 11);
1520  Particle p2({ 0.1 , -0.4, 0.8, 4.0 }, 11);
1521 
1522  StoreArray<MCParticle> mcParticles;
1523  auto* mcParticle = mcParticles.appendNew();
1524  mcParticle->setPDG(11);
1525  mcParticle->setStatus(MCParticle::c_PrimaryParticle);
1526  mcParticle = mcParticles.appendNew();
1527  mcParticle->setPDG(22);
1528  mcParticle->setStatus(MCParticle::c_PrimaryParticle);
1529  mcParticle = mcParticles.appendNew();
1530  mcParticle->setPDG(-11);
1531  mcParticle->setStatus(MCParticle::c_PrimaryParticle);
1532  mcParticle = mcParticles.appendNew();
1533  mcParticle->setPDG(11);
1534 
1535 
1536  const Manager::Var* var = Manager::Instance().getVariable("NumberOfMCParticlesInEvent(11)");
1537  ASSERT_NE(var, nullptr);
1538  EXPECT_FLOAT_EQ(var->function(nullptr), 2);
1539 
1540  }
1541 
1542  TEST_F(MetaVariableTest, daughterInvM)
1543  {
1544  TLorentzVector momentum;
1545  const int nDaughters = 6;
1546  StoreArray<Particle> particles;
1547  std::vector<int> daughterIndices;
1548  for (int i = 0; i < nDaughters; i++) {
1549  Particle d(TLorentzVector(2, 2, 2, 4.0), (i % 2) ? 211 : -211);
1550  momentum += d.get4Vector();
1551  Particle* newDaughters = particles.appendNew(d);
1552  daughterIndices.push_back(newDaughters->getArrayIndex());
1553  }
1554  const Particle* p = particles.appendNew(momentum, 411, Particle::c_Unflavored, daughterIndices);
1555 
1556  const Manager::Var* var = Manager::Instance().getVariable("daughterInvM(6,5)");
1557  ASSERT_NE(var, nullptr);
1558  EXPECT_TRUE(std::isnan(var->function(p)));
1559 
1560  var = Manager::Instance().getVariable("daughterInvM(0, 1)");
1561  ASSERT_NE(var, nullptr);
1562  EXPECT_FLOAT_EQ(var->function(p), 4.0);
1563 
1564  var = Manager::Instance().getVariable("daughterInvM(0, 1, 2)");
1565  ASSERT_NE(var, nullptr);
1566  EXPECT_FLOAT_EQ(var->function(p), 6.0);
1567  }
1568 
1569  TEST_F(MetaVariableTest, daughter)
1570  {
1571  TLorentzVector momentum;
1572  const int nDaughters = 6;
1573  StoreArray<Particle> particles;
1574  std::vector<int> daughterIndices;
1575  for (int i = 0; i < nDaughters; i++) {
1576  Particle d(TLorentzVector(i * 1.0, 1, 1, 1), (i % 2) ? 211 : -211);
1577  momentum += d.get4Vector();
1578  Particle* newDaughters = particles.appendNew(d);
1579  daughterIndices.push_back(newDaughters->getArrayIndex());
1580  }
1581  const Particle* p = particles.appendNew(momentum, 411, Particle::c_Unflavored, daughterIndices);
1582 
1583  const Manager::Var* var = Manager::Instance().getVariable("daughter(6, px)");
1584  ASSERT_NE(var, nullptr);
1585  EXPECT_TRUE(std::isnan(var->function(p)));
1586 
1587  var = Manager::Instance().getVariable("daughter(0, px)");
1588  ASSERT_NE(var, nullptr);
1589  EXPECT_ALL_NEAR(var->function(p), 0.0, 1e-6);
1590 
1591  var = Manager::Instance().getVariable("daughter(1, px)");
1592  ASSERT_NE(var, nullptr);
1593  EXPECT_FLOAT_EQ(var->function(p), 1.0);
1594 
1595  var = Manager::Instance().getVariable("daughter(2, px)");
1596  ASSERT_NE(var, nullptr);
1597  EXPECT_FLOAT_EQ(var->function(p), 2.0);
1598  }
1599 
1600  TEST_F(MetaVariableTest, mcDaughter)
1601  {
1602  DataStore::Instance().setInitializeActive(true);
1603  StoreArray<MCParticle> mcParticles;
1604  StoreArray<Particle> particles;
1605  particles.registerInDataStore();
1606  mcParticles.registerInDataStore();
1607  particles.registerRelationTo(mcParticles);
1608  DataStore::Instance().setInitializeActive(false);
1609 
1610  // Create MC graph for B -> (muon -> electron + muon_neutrino) + anti_muon_neutrino
1611  MCParticleGraph mcGraph;
1612 
1613  MCParticleGraph::GraphParticle& graphParticleGrandMother = mcGraph.addParticle();
1614 
1615  MCParticleGraph::GraphParticle& graphParticleMother = mcGraph.addParticle();
1616  MCParticleGraph::GraphParticle& graphParticleAunt = mcGraph.addParticle();
1617 
1618  MCParticleGraph::GraphParticle& graphParticleDaughter1 = mcGraph.addParticle();
1619  MCParticleGraph::GraphParticle& graphParticleDaughter2 = mcGraph.addParticle();
1620 
1621  graphParticleGrandMother.setPDG(-521);
1622  graphParticleMother.setPDG(13);
1623  graphParticleAunt.setPDG(-14);
1624  graphParticleDaughter1.setPDG(11);
1625  graphParticleDaughter2.setPDG(14);
1626 
1627  graphParticleMother.comesFrom(graphParticleGrandMother);
1628  graphParticleAunt.comesFrom(graphParticleGrandMother);
1629  graphParticleDaughter1.comesFrom(graphParticleMother);
1630  graphParticleDaughter2.comesFrom(graphParticleMother);
1631  mcGraph.generateList();
1632 
1633  // Get MC Particles from StoreArray
1634  auto* mcGrandMother = mcParticles[0];
1635  mcGrandMother->setStatus(MCParticle::c_PrimaryParticle);
1636 
1637  auto* mcMother = mcParticles[1];
1638  mcMother->setStatus(MCParticle::c_PrimaryParticle);
1639 
1640  auto* mcAunt = mcParticles[2];
1641  mcAunt->setStatus(MCParticle::c_PrimaryParticle);
1642 
1643  auto* mcDaughter1 = mcParticles[3];
1644  mcDaughter1->setStatus(MCParticle::c_PrimaryParticle);
1645 
1646  auto* mcDaughter2 = mcParticles[4];
1647  mcDaughter2->setStatus(MCParticle::c_PrimaryParticle);
1648 
1649  auto* pGrandMother = particles.appendNew(TLorentzVector({ 0.0 , -0.4, 0.8, 1.0}), -521);
1650  pGrandMother->addRelationTo(mcGrandMother);
1651 
1652  auto* pMother = particles.appendNew(TLorentzVector({ 0.0 , -0.4, 0.8, 1.0}), 13);
1653  pMother->addRelationTo(mcMother);
1654 
1655  // Test for particle that has no MC match
1656  auto* p_noMC = particles.appendNew(TLorentzVector({ 0.0 , -0.4, 0.8, 1.0}), 13);
1657 
1658  // Test for particle that has MC match, but MC match has no daughter
1659  auto* p_noDaughter = particles.appendNew(TLorentzVector({ 0.0 , -0.4, 0.8, 1.0}), 11);
1660  p_noDaughter->addRelationTo(mcDaughter1);
1661 
1662  const Manager::Var* var = Manager::Instance().getVariable("mcDaughter(0, PDG)");
1663  ASSERT_NE(var, nullptr);
1664  EXPECT_FLOAT_EQ(var->function(pGrandMother), 13);
1665  EXPECT_FLOAT_EQ(var->function(pMother), 11);
1666  EXPECT_TRUE(std::isnan(var->function(p_noMC)));
1667  EXPECT_TRUE(std::isnan(var->function(p_noDaughter)));
1668  var = Manager::Instance().getVariable("mcDaughter(1, PDG)");
1669  EXPECT_FLOAT_EQ(var->function(pGrandMother), -14);
1670  EXPECT_FLOAT_EQ(var->function(pMother), 14);
1671  // Test for particle where mc daughter index is out of range of mc daughters
1672  var = Manager::Instance().getVariable("mcDaughter(2, PDG)");
1673  EXPECT_TRUE(std::isnan(var->function(pGrandMother)));
1674  EXPECT_TRUE(std::isnan(var->function(pMother)));
1675  // Test nested application of mcDaughter
1676  var = Manager::Instance().getVariable("mcDaughter(0, mcDaughter(0, PDG))");
1677  EXPECT_FLOAT_EQ(var->function(pGrandMother), 11);
1678  EXPECT_TRUE(std::isnan(var->function(pMother)));
1679  var = Manager::Instance().getVariable("mcDaughter(0, mcDaughter(1, PDG))");
1680  EXPECT_FLOAT_EQ(var->function(pGrandMother), 14);
1681  var = Manager::Instance().getVariable("mcDaughter(0, mcDaughter(2, PDG))");
1682  EXPECT_TRUE(std::isnan(var->function(pGrandMother)));
1683  var = Manager::Instance().getVariable("mcDaughter(1, mcDaughter(0, PDG))");
1684  EXPECT_TRUE(std::isnan(var->function(pGrandMother)));
1685  }
1686 
1687  TEST_F(MetaVariableTest, mcMother)
1688  {
1689  DataStore::Instance().setInitializeActive(true);
1690  StoreArray<MCParticle> mcParticles;
1691  StoreArray<Particle> particles;
1692  particles.registerInDataStore();
1693  mcParticles.registerInDataStore();
1694  particles.registerRelationTo(mcParticles);
1695  DataStore::Instance().setInitializeActive(false);
1696 
1697  // Create MC graph for B -> (muon -> electron + muon_neutrino) + anti_muon_neutrino
1698  MCParticleGraph mcGraph;
1699 
1700  MCParticleGraph::GraphParticle& graphParticleGrandMother = mcGraph.addParticle();
1701 
1702  MCParticleGraph::GraphParticle& graphParticleMother = mcGraph.addParticle();
1703  MCParticleGraph::GraphParticle& graphParticleAunt = mcGraph.addParticle();
1704 
1705  MCParticleGraph::GraphParticle& graphParticleDaughter1 = mcGraph.addParticle();
1706  MCParticleGraph::GraphParticle& graphParticleDaughter2 = mcGraph.addParticle();
1707 
1708  graphParticleGrandMother.setPDG(-521);
1709  graphParticleMother.setPDG(13);
1710  graphParticleAunt.setPDG(-14);
1711  graphParticleDaughter1.setPDG(11);
1712  graphParticleDaughter2.setPDG(14);
1713 
1714  graphParticleMother.comesFrom(graphParticleGrandMother);
1715  graphParticleAunt.comesFrom(graphParticleGrandMother);
1716  graphParticleDaughter1.comesFrom(graphParticleMother);
1717  graphParticleDaughter2.comesFrom(graphParticleMother);
1718 
1719  mcGraph.generateList();
1720 
1721  // Get MC Particles from StoreArray
1722  auto* mcGrandMother = mcParticles[0];
1723  mcGrandMother->setStatus(MCParticle::c_PrimaryParticle);
1724 
1725  auto* mcMother = mcParticles[1];
1726  mcMother->setStatus(MCParticle::c_PrimaryParticle);
1727 
1728  auto* mcAunt = mcParticles[2];
1729  mcAunt->setStatus(MCParticle::c_PrimaryParticle);
1730 
1731  auto* mcDaughter1 = mcParticles[3];
1732  mcDaughter1->setStatus(MCParticle::c_PrimaryParticle);
1733 
1734  auto* mcDaughter2 = mcParticles[4];
1735  mcDaughter2->setStatus(MCParticle::c_PrimaryParticle);
1736 
1737  auto* p1 = particles.appendNew(TLorentzVector({ 0.0 , -0.4, 0.8, 1.0}), 11);
1738  p1->addRelationTo(mcDaughter1);
1739 
1740  auto* p2 = particles.appendNew(TLorentzVector({ 0.0 , -0.4, 0.8, 1.0}), 14);
1741  p2->addRelationTo(mcDaughter2);
1742 
1743  auto* pMother = particles.appendNew(TLorentzVector({ 0.0 , -0.4, 0.8, 1.0}), 13);
1744  pMother->addRelationTo(mcMother);
1745 
1746  // For test of particle that has no MC match
1747  auto* p_noMC = particles.appendNew(TLorentzVector({ 0.0 , -0.4, 0.8, 1.0}), 11);
1748 
1749  // For test of particle that has MC match, but MC match has no mother
1750  auto* p_noMother = particles.appendNew(TLorentzVector({ 0.0 , -0.4, 0.8, 1.0}), -521);
1751  p_noMother->addRelationTo(mcGrandMother);
1752 
1753  const Manager::Var* var = Manager::Instance().getVariable("mcMother(PDG)");
1754  ASSERT_NE(var, nullptr);
1755  EXPECT_FLOAT_EQ(var->function(p1), 13);
1756  EXPECT_FLOAT_EQ(var->function(p2), 13);
1757  EXPECT_FLOAT_EQ(var->function(pMother), -521);
1758  EXPECT_TRUE(std::isnan(var->function(p_noMC)));
1759  EXPECT_TRUE(std::isnan(var->function(p_noMother)));
1760 
1761  // Test if nested calls of mcMother work correctly
1762  var = Manager::Instance().getVariable("mcMother(mcMother(PDG))");
1763  EXPECT_FLOAT_EQ(var->function(p1), -521);
1764  }
1765 
1766  TEST_F(MetaVariableTest, genParticle)
1767  {
1768  DataStore::Instance().setInitializeActive(true);
1769  StoreArray<MCParticle> mcParticles;
1770  StoreArray<Particle> particles;
1771  particles.registerInDataStore();
1772  mcParticles.registerInDataStore();
1773  particles.registerRelationTo(mcParticles);
1774  DataStore::Instance().setInitializeActive(false);
1775 
1776  // Create MC graph for Upsilon(4S) -> (B^- -> electron + anti_electron_neutrino) + B^+
1777  MCParticleGraph mcGraph;
1778 
1779  MCParticleGraph::GraphParticle& graphParticleGrandMother = mcGraph.addParticle();
1780 
1781  MCParticleGraph::GraphParticle& graphParticleMother = mcGraph.addParticle();
1782  MCParticleGraph::GraphParticle& graphParticleAunt = mcGraph.addParticle();
1783 
1784  MCParticleGraph::GraphParticle& graphParticleDaughter1 = mcGraph.addParticle();
1785  MCParticleGraph::GraphParticle& graphParticleDaughter2 = mcGraph.addParticle();
1786 
1787  graphParticleGrandMother.setPDG(300553);
1788  graphParticleMother.setPDG(-521);
1789  graphParticleAunt.setPDG(521);
1790  graphParticleDaughter1.setPDG(11);
1791  graphParticleDaughter2.setPDG(-12);
1792 
1793  graphParticleGrandMother.setMomentum(0.0, 0.0, 0.4);
1794  graphParticleMother.setMomentum(1.1, 1.3, 1.5);
1795 
1796  graphParticleMother.comesFrom(graphParticleGrandMother);
1797  graphParticleAunt.comesFrom(graphParticleGrandMother);
1798  graphParticleDaughter1.comesFrom(graphParticleMother);
1799  graphParticleDaughter2.comesFrom(graphParticleMother);
1800 
1801  mcGraph.generateList();
1802 
1803  // Get MC Particles from StoreArray
1804  auto* mcGrandMother = mcParticles[0];
1805  mcGrandMother->setStatus(MCParticle::c_PrimaryParticle);
1806 
1807  auto* mcMother = mcParticles[1];
1808  mcMother->setStatus(MCParticle::c_PrimaryParticle);
1809 
1810  auto* mcAunt = mcParticles[2];
1811  mcAunt->setStatus(MCParticle::c_PrimaryParticle);
1812 
1813  auto* mcDaughter1 = mcParticles[3];
1814  mcDaughter1->setStatus(MCParticle::c_PrimaryParticle);
1815 
1816  auto* mcDaughter2 = mcParticles[4];
1817  mcDaughter2->setStatus(MCParticle::c_PrimaryParticle);
1818 
1819  auto* p1 = particles.appendNew(TLorentzVector({ 0.0 , -0.4, 0.8, 1.0}), 11);
1820  p1->addRelationTo(mcDaughter1);
1821 
1822  // For test of particle that has no MC match
1823  auto* p_noMC = particles.appendNew(TLorentzVector({ 0.0 , -0.4, 0.8, 1.0}), 211);
1824 
1825  const Manager::Var* var = Manager::Instance().getVariable("genParticle(0, PDG)");
1826  ASSERT_NE(var, nullptr);
1827  EXPECT_FLOAT_EQ(var->function(p1), 300553);
1828  EXPECT_FLOAT_EQ(var->function(p_noMC), 300553);
1829 
1830  var = Manager::Instance().getVariable("genParticle(0, matchedMC(pz))");
1831  ASSERT_NE(var, nullptr);
1832  EXPECT_FLOAT_EQ(var->function(p1), 0.4);
1833  EXPECT_FLOAT_EQ(var->function(p_noMC), 0.4);
1834 
1835  var = Manager::Instance().getVariable("genParticle(0, mcDaughter(0, PDG))");
1836  ASSERT_NE(var, nullptr);
1837  EXPECT_FLOAT_EQ(var->function(p1), -521);
1838  EXPECT_FLOAT_EQ(var->function(p_noMC), -521);
1839 
1840  var = Manager::Instance().getVariable("genParticle(0, mcDaughter(0, matchedMC(px)))");
1841  ASSERT_NE(var, nullptr);
1842  EXPECT_FLOAT_EQ(var->function(p1), 1.1);
1843  EXPECT_FLOAT_EQ(var->function(p_noMC), 1.1);
1844 
1845  var = Manager::Instance().getVariable("genParticle(1, PDG)");
1846  ASSERT_NE(var, nullptr);
1847  EXPECT_FLOAT_EQ(var->function(p1), -521);
1848  EXPECT_FLOAT_EQ(var->function(p_noMC), -521);
1849 
1850  var = Manager::Instance().getVariable("genParticle(4, PDG)");
1851  ASSERT_NE(var, nullptr);
1852  EXPECT_FLOAT_EQ(var->function(p1), -12);
1853  EXPECT_FLOAT_EQ(var->function(p_noMC), -12);
1854 
1855  var = Manager::Instance().getVariable("genParticle(5, PDG)");
1856  ASSERT_NE(var, nullptr);
1857  EXPECT_TRUE(std::isnan(var->function(p1)));
1858  EXPECT_TRUE(std::isnan(var->function(p_noMC)));
1859  }
1860 
1861  TEST_F(MetaVariableTest, genUpsilon4S)
1862  {
1863  DataStore::Instance().setInitializeActive(true);
1864  StoreArray<MCParticle> mcParticles;
1865  StoreArray<Particle> particles;
1866  particles.registerInDataStore();
1867  mcParticles.registerInDataStore();
1868  particles.registerRelationTo(mcParticles);
1869  DataStore::Instance().setInitializeActive(false);
1870 
1871  // Create MC graph for Upsilon(4S) -> (B^- -> electron + anti_electron_neutrino) + B^+
1872  MCParticleGraph mcGraph;
1873 
1874  MCParticleGraph::GraphParticle& graphParticleGrandMother = mcGraph.addParticle();
1875 
1876  MCParticleGraph::GraphParticle& graphParticleMother = mcGraph.addParticle();
1877  MCParticleGraph::GraphParticle& graphParticleAunt = mcGraph.addParticle();
1878 
1879  MCParticleGraph::GraphParticle& graphParticleDaughter1 = mcGraph.addParticle();
1880  MCParticleGraph::GraphParticle& graphParticleDaughter2 = mcGraph.addParticle();
1881 
1882  graphParticleGrandMother.setPDG(300553);
1883  graphParticleMother.setPDG(-521);
1884  graphParticleAunt.setPDG(521);
1885  graphParticleDaughter1.setPDG(11);
1886  graphParticleDaughter2.setPDG(-12);
1887 
1888  graphParticleGrandMother.setMomentum(0.0, 0.0, 0.4);
1889  graphParticleMother.setMomentum(1.1, 1.3, 1.5);
1890 
1891  graphParticleMother.comesFrom(graphParticleGrandMother);
1892  graphParticleAunt.comesFrom(graphParticleGrandMother);
1893  graphParticleDaughter1.comesFrom(graphParticleMother);
1894  graphParticleDaughter2.comesFrom(graphParticleMother);
1895 
1896  mcGraph.generateList();
1897 
1898  // Get MC Particles from StoreArray
1899  auto* mcGrandMother = mcParticles[0];
1900  mcGrandMother->setStatus(MCParticle::c_PrimaryParticle);
1901 
1902  auto* mcMother = mcParticles[1];
1903  mcMother->setStatus(MCParticle::c_PrimaryParticle);
1904 
1905  auto* mcAunt = mcParticles[2];
1906  mcAunt->setStatus(MCParticle::c_PrimaryParticle);
1907 
1908  auto* mcDaughter1 = mcParticles[3];
1909  mcDaughter1->setStatus(MCParticle::c_PrimaryParticle);
1910 
1911  auto* mcDaughter2 = mcParticles[4];
1912  mcDaughter2->setStatus(MCParticle::c_PrimaryParticle);
1913 
1914  auto* p1 = particles.appendNew(TLorentzVector({ 0.0 , -0.4, 0.8, 1.0}), 11);
1915  p1->addRelationTo(mcDaughter1);
1916 
1917  // For test of particle that has no MC match
1918  auto* p_noMC = particles.appendNew(TLorentzVector({ 0.0 , -0.4, 0.8, 1.0}), 211);
1919 
1920  const Manager::Var* var = Manager::Instance().getVariable("genUpsilon4S(PDG)");
1921  ASSERT_NE(var, nullptr);
1922  EXPECT_FLOAT_EQ(var->function(p1), 300553);
1923  EXPECT_FLOAT_EQ(var->function(p_noMC), 300553);
1924 
1925  var = Manager::Instance().getVariable("genUpsilon4S(matchedMC(pz))");
1926  ASSERT_NE(var, nullptr);
1927  EXPECT_FLOAT_EQ(var->function(p1), 0.4);
1928  EXPECT_FLOAT_EQ(var->function(p_noMC), 0.4);
1929 
1930  var = Manager::Instance().getVariable("genUpsilon4S(mcDaughter(0, PDG))");
1931  ASSERT_NE(var, nullptr);
1932  EXPECT_FLOAT_EQ(var->function(p1), -521);
1933  EXPECT_FLOAT_EQ(var->function(p_noMC), -521);
1934 
1935  var = Manager::Instance().getVariable("genUpsilon4S(mcDaughter(0, matchedMC(px)))");
1936  ASSERT_NE(var, nullptr);
1937  EXPECT_FLOAT_EQ(var->function(p1), 1.1);
1938  EXPECT_FLOAT_EQ(var->function(p_noMC), 1.1);
1939 
1941  mcParticles.clear();
1942  particles.clear();
1943  MCParticleGraph mcGraph2;
1944 
1945  MCParticleGraph::GraphParticle& graphParticle1 = mcGraph2.addParticle();
1946  MCParticleGraph::GraphParticle& graphParticle2 = mcGraph2.addParticle();
1947 
1948  graphParticle1.setPDG(11);
1949  graphParticle2.setPDG(-11);
1950 
1951  graphParticle1.setMomentum(1.1, 1.3, 1.4);
1952  graphParticle1.setMomentum(-1.1, -1.3, 1.4);
1953 
1954  mcGraph2.generateList();
1955 
1956  auto* mcP1 = mcParticles[0];
1957  mcP1->setStatus(MCParticle::c_PrimaryParticle);
1958 
1959  auto* mcP2 = mcParticles[1];
1960  mcP2->setStatus(MCParticle::c_PrimaryParticle);
1961 
1962  auto* someParticle = particles.appendNew(TLorentzVector({ 0.0 , -0.4, 0.8, 1.0}), 11);
1963  someParticle->addRelationTo(mcP1);
1964 
1965  var = Manager::Instance().getVariable("genUpsilon4S(PDG)");
1966  ASSERT_NE(var, nullptr);
1967  EXPECT_TRUE(std::isnan(var->function(someParticle)));
1968  }
1969 
1970  TEST_F(MetaVariableTest, daughterProductOf)
1971  {
1972  TLorentzVector momentum;
1973  const int nDaughters = 4;
1974  StoreArray<Particle> particles;
1975  std::vector<int> daughterIndices;
1976  for (int i = 0; i < nDaughters; i++) {
1977  Particle d(TLorentzVector(1, 1, 1, i * 1.0 + 1.0), (i % 2) ? 211 : -211);
1978  momentum += d.get4Vector();
1979  Particle* newDaughters = particles.appendNew(d);
1980  daughterIndices.push_back(newDaughters->getArrayIndex());
1981  }
1982  const Particle* p = particles.appendNew(momentum, 411, Particle::c_Unflavored, daughterIndices);
1983 
1984  const Manager::Var* var = Manager::Instance().getVariable("daughterProductOf(E)");
1985  ASSERT_NE(var, nullptr);
1986  EXPECT_FLOAT_EQ(var->function(p), 24.0);
1987  }
1988 
1989  TEST_F(MetaVariableTest, daughterSumOf)
1990  {
1991  TLorentzVector momentum;
1992  const int nDaughters = 4;
1993  StoreArray<Particle> particles;
1994  std::vector<int> daughterIndices;
1995  for (int i = 0; i < nDaughters; i++) {
1996  Particle d(TLorentzVector(1, 1, 1, i * 1.0 + 1.0), (i % 2) ? 211 : -211);
1997  momentum += d.get4Vector();
1998  Particle* newDaughters = particles.appendNew(d);
1999  daughterIndices.push_back(newDaughters->getArrayIndex());
2000  }
2001  const Particle* p = particles.appendNew(momentum, 411, Particle::c_Unflavored, daughterIndices);
2002 
2003  const Manager::Var* var = Manager::Instance().getVariable("daughterSumOf(E)");
2004  ASSERT_NE(var, nullptr);
2005  EXPECT_FLOAT_EQ(var->function(p), 10.0);
2006 
2007  }
2008 
2009  TEST_F(MetaVariableTest, daughterLowest)
2010  {
2011  TLorentzVector momentum;
2012  const int nDaughters = 4;
2013  StoreArray<Particle> particles;
2014  std::vector<int> daughterIndices;
2015  for (int i = 0; i < nDaughters; i++) {
2016  Particle d(TLorentzVector(1, 1, 1, i * 1.0 + 1.0), (i % 2) ? 211 : -211);
2017  momentum += d.get4Vector();
2018  Particle* newDaughters = particles.appendNew(d);
2019  daughterIndices.push_back(newDaughters->getArrayIndex());
2020  }
2021  const Particle* p = particles.appendNew(momentum, 411, Particle::c_Unflavored, daughterIndices);
2022 
2023  const Manager::Var* var = Manager::Instance().getVariable("daughterLowest(E)");
2024  ASSERT_NE(var, nullptr);
2025  EXPECT_FLOAT_EQ(var->function(p), 1.0);
2026  }
2027 
2028  TEST_F(MetaVariableTest, daughterHighest)
2029  {
2030  TLorentzVector momentum;
2031  const int nDaughters = 4;
2032  StoreArray<Particle> particles;
2033  std::vector<int> daughterIndices;
2034  for (int i = 0; i < nDaughters; i++) {
2035  Particle d(TLorentzVector(1, 1, 1, i * 1.0 + 1.0), (i % 2) ? 211 : -211);
2036  momentum += d.get4Vector();
2037  Particle* newDaughters = particles.appendNew(d);
2038  daughterIndices.push_back(newDaughters->getArrayIndex());
2039  }
2040  const Particle* p = particles.appendNew(momentum, 411, Particle::c_Unflavored, daughterIndices);
2041 
2042  const Manager::Var* var = Manager::Instance().getVariable("daughterHighest(E)");
2043  ASSERT_NE(var, nullptr);
2044  EXPECT_FLOAT_EQ(var->function(p), 4.0);
2045  }
2046 
2047  TEST_F(MetaVariableTest, daughterDiffOf)
2048  {
2049  TLorentzVector momentum;
2050  const int nDaughters = 4;
2051  StoreArray<Particle> particles;
2052  std::vector<int> daughterIndices;
2053  for (int i = 0; i < nDaughters; i++) {
2054  Particle d(TLorentzVector(1, 1, 1, i * 1.0 + 1.0), (i % 2) ? -11 : 211);
2055  momentum += d.get4Vector();
2056  Particle* newDaughters = particles.appendNew(d);
2057  daughterIndices.push_back(newDaughters->getArrayIndex());
2058  }
2059  const Particle* p = particles.appendNew(momentum, 411, Particle::c_Unflavored, daughterIndices);
2060 
2061  const Manager::Var* var = Manager::Instance().getVariable("daughterDiffOf(0, 1, PDG)");
2062  ASSERT_NE(var, nullptr);
2063  EXPECT_FLOAT_EQ(var->function(p), -222);
2064 
2065  var = Manager::Instance().getVariable("daughterDiffOf(1, 0, PDG)");
2066  ASSERT_NE(var, nullptr);
2067  EXPECT_FLOAT_EQ(var->function(p), 222);
2068 
2069  var = Manager::Instance().getVariable("daughterDiffOf(0, 1, abs(PDG))");
2070  ASSERT_NE(var, nullptr);
2071  EXPECT_FLOAT_EQ(var->function(p), -200);
2072 
2073  var = Manager::Instance().getVariable("daughterDiffOf(1, 1, PDG)");
2074  ASSERT_NE(var, nullptr);
2075  EXPECT_FLOAT_EQ(var->function(p), 0);
2076 
2077  var = Manager::Instance().getVariable("daughterDiffOf(1, 3, abs(PDG))");
2078  ASSERT_NE(var, nullptr);
2079  EXPECT_FLOAT_EQ(var->function(p), 0);
2080 
2081  var = Manager::Instance().getVariable("daughterDiffOf(0, 2, PDG)");
2082  ASSERT_NE(var, nullptr);
2083  EXPECT_FLOAT_EQ(var->function(p), 0);
2084 
2085  EXPECT_B2FATAL(Manager::Instance().getVariable("daughterDiffOf(0, NOTINT, PDG)"));
2086  }
2087 
2088  TEST_F(MetaVariableTest, daughterClusterAngleInBetween)
2089  {
2090  // declare all the array we need
2091  StoreArray<Particle> particles, particles_noclst;
2092  std::vector<int> daughterIndices, daughterIndices_noclst;
2093 
2094  //proxy initialize where to declare the needed array
2095  DataStore::Instance().setInitializeActive(true);
2096  StoreArray<ECLCluster> eclclusters;
2097  eclclusters.registerInDataStore();
2098  particles.registerRelationTo(eclclusters);
2099  DataStore::Instance().setInitializeActive(false);
2100 
2101  // create two Lorentz vectors that are back to back in the CMS and boost them to the Lab frame
2102  const float px_CM = 2.;
2103  const float py_CM = 1.;
2104  const float pz_CM = 3.;
2105  float E_CM;
2106  E_CM = sqrt(pow(px_CM, 2) + pow(py_CM, 2) + pow(pz_CM, 2));
2107  TLorentzVector momentum, momentum_noclst;
2108  TLorentzVector dau0_4vec_CM(px_CM, py_CM, pz_CM, E_CM), dau1_4vec_CM(-px_CM, -py_CM, -pz_CM, E_CM);
2109  TLorentzVector dau0_4vec_Lab, dau1_4vec_Lab;
2110  dau0_4vec_Lab = PCmsLabTransform::cmsToLab(
2111  dau0_4vec_CM); //why is everybody using the extended method when there are the functions that do all the steps for us?
2112  dau1_4vec_Lab = PCmsLabTransform::cmsToLab(dau1_4vec_CM);
2113 
2114  // add the two photons (now in the Lab frame) as the two daughters of some particle and create the latter
2115  Particle dau0_noclst(dau0_4vec_Lab, 22);
2116  momentum += dau0_noclst.get4Vector();
2117  Particle* newDaughter0_noclst = particles.appendNew(dau0_noclst);
2118  daughterIndices_noclst.push_back(newDaughter0_noclst->getArrayIndex());
2119  Particle dau1_noclst(dau1_4vec_Lab, 22);
2120  momentum += dau1_noclst.get4Vector();
2121  Particle* newDaughter1_noclst = particles.appendNew(dau1_noclst);
2122  daughterIndices_noclst.push_back(newDaughter1_noclst->getArrayIndex());
2123  const Particle* par_noclst = particles.appendNew(momentum, 111, Particle::c_Unflavored, daughterIndices_noclst);
2124 
2125  // grab variables
2126  const Manager::Var* var = Manager::Instance().getVariable("daughterClusterAngleInBetween(0, 1)");
2127  const Manager::Var* varCMS = Manager::Instance().getVariable("useCMSFrame(daughterClusterAngleInBetween(0, 1))");
2128 
2129  // when no relations are set between the particles and the eclClusters, nan is expected to be returned
2130  ASSERT_NE(var, nullptr);
2131  EXPECT_TRUE(std::isnan(var->function(par_noclst)));
2132 
2133  // set relations between particles and eclClusters
2134  ECLCluster* eclst0 = eclclusters.appendNew(ECLCluster());
2135  eclst0->setEnergy(dau0_4vec_Lab.E());
2136  eclst0->setHypothesis(ECLCluster::EHypothesisBit::c_nPhotons);
2137  eclst0->setClusterId(1);
2138  eclst0->setTheta(dau0_4vec_Lab.Theta());
2139  eclst0->setPhi(dau0_4vec_Lab.Phi());
2140  eclst0->setR(148.4);
2141  ECLCluster* eclst1 = eclclusters.appendNew(ECLCluster());
2142  eclst1->setEnergy(dau1_4vec_Lab.E());
2143  eclst1->setHypothesis(ECLCluster::EHypothesisBit::c_nPhotons);
2144  eclst1->setClusterId(2);
2145  eclst1->setTheta(dau1_4vec_Lab.Theta());
2146  eclst1->setPhi(dau1_4vec_Lab.Phi());
2147  eclst1->setR(148.5);
2148 
2149  const Particle* newDaughter0 = particles.appendNew(Particle(eclclusters[0]));
2150  daughterIndices.push_back(newDaughter0->getArrayIndex());
2151  const Particle* newDaughter1 = particles.appendNew(Particle(eclclusters[1]));
2152  daughterIndices.push_back(newDaughter1->getArrayIndex());
2153 
2154  const Particle* par = particles.appendNew(momentum, 111, Particle::c_Unflavored, daughterIndices);
2155 
2156  //now we expect non-nan results
2157  EXPECT_FLOAT_EQ(var->function(par), 2.8638029);
2158  EXPECT_FLOAT_EQ(varCMS->function(par), M_PI);
2159  }
2160 
2161  TEST_F(MetaVariableTest, grandDaughterDiffOfs)
2162  {
2163  // declare all the array we need
2164  StoreArray<Particle> particles, particles_noclst;
2165  std::vector<int> daughterIndices0_noclst, daughterIndices1_noclst, daughterIndices2_noclst;
2166  std::vector<int> daughterIndices0, daughterIndices1, daughterIndices2;
2167 
2168  //proxy initialize where to declare the needed array
2169  DataStore::Instance().setInitializeActive(true);
2170  StoreArray<ECLCluster> eclclusters;
2171  eclclusters.registerInDataStore();
2172  particles.registerRelationTo(eclclusters);
2173  DataStore::Instance().setInitializeActive(false);
2174 
2175  // create two Lorentz vectors
2176  const float px_0 = 2.;
2177  const float py_0 = 1.;
2178  const float pz_0 = 3.;
2179  const float px_1 = 1.5;
2180  const float py_1 = 1.5;
2181  const float pz_1 = 2.5;
2182  float E_0, E_1;
2183  E_0 = sqrt(pow(px_0, 2) + pow(py_0, 2) + pow(pz_0, 2));
2184  E_1 = sqrt(pow(px_1, 2) + pow(py_1, 2) + pow(pz_1, 2));
2185  TLorentzVector momentum_0, momentum_1, momentum;
2186  TLorentzVector dau0_4vec(px_0, py_0, pz_0, E_0), dau1_4vec(px_1, py_1, pz_1, E_1);
2187 
2188  // add the two photons as the two daughters of some particle and create the latter
2189  // Particle dau0_noclst(dau0_4vec, 22);
2190  // momentum += dau0_noclst.get4Vector();
2191  // Particle* newDaughter0_noclst = particles.appendNew(dau0_noclst);
2192  // daughterIndices_noclst.push_back(newDaughter0_noclst->getArrayIndex());
2193  // Particle dau1_noclst(dau1_4vec, 22);
2194  // momentum += dau1_noclst.get4Vector();
2195  // Particle* newDaughter1_noclst = particles.appendNew(dau1_noclst);
2196  // daughterIndices_noclst.push_back(newDaughter1_noclst->getArrayIndex());
2197  // const Particle* par_noclst = particles.appendNew(momentum, 111, Particle::c_Unflavored, daughterIndices_noclst);
2198 
2199  Particle dau0_noclst(dau0_4vec, 22);
2200  momentum_0 = dau0_4vec;
2201  Particle* newDaughter0_noclst = particles.appendNew(dau0_noclst);
2202  daughterIndices0_noclst.push_back(newDaughter0_noclst->getArrayIndex());
2203  const Particle* par0_noclst = particles.appendNew(momentum_0, 111, Particle::c_Unflavored, daughterIndices0_noclst);
2204  Particle dau1_noclst(dau1_4vec, 22);
2205  momentum_1 = dau1_4vec;
2206  Particle* newDaughter1_noclst = particles.appendNew(dau1_noclst);
2207  daughterIndices1_noclst.push_back(newDaughter1_noclst->getArrayIndex());
2208  const Particle* par1_noclst = particles.appendNew(momentum_1, 111, Particle::c_Unflavored, daughterIndices1_noclst);
2209 
2210  momentum = momentum_0 + momentum_1;
2211  daughterIndices2_noclst.push_back(par0_noclst->getArrayIndex());
2212  daughterIndices2_noclst.push_back(par1_noclst->getArrayIndex());
2213  const Particle* parGranny_noclst = particles.appendNew(momentum, 111, Particle::c_Unflavored, daughterIndices2_noclst);
2214 
2215  // grab variables
2216  const Manager::Var* var_Theta = Manager::Instance().getVariable("grandDaughterDiffOf(0,1,0,0,theta)");
2217  const Manager::Var* var_ClusterTheta = Manager::Instance().getVariable("grandDaughterDiffOf(0,1,0,0,clusterTheta)");
2218  const Manager::Var* var_E = Manager::Instance().getVariable("grandDaughterDiffOf(0,1,0,0,E)");
2219  const Manager::Var* var_ClusterE = Manager::Instance().getVariable("grandDaughterDiffOf(0,1,0,0,clusterE)");
2220  const Manager::Var* var_E_wrongIndexes = Manager::Instance().getVariable("grandDaughterDiffOf(0,1,2,3,E)");
2221  const Manager::Var* var_ClusterE_wrongIndexes = Manager::Instance().getVariable("grandDaughterDiffOf(0,1,2,3,clusterE)");
2222 
2223  const Manager::Var* var_ClusterPhi = Manager::Instance().getVariable("grandDaughterDiffOfClusterPhi(0,1,0,0)");
2224  const Manager::Var* var_Phi = Manager::Instance().getVariable("grandDaughterDiffOfPhi(0,1,0,0)");
2225  const Manager::Var* var_ClusterPhi_wrongIndexes = Manager::Instance().getVariable("grandDaughterDiffOfClusterPhi(0,1,2,3)");
2226  const Manager::Var* var_Phi_wrongIndexes = Manager::Instance().getVariable("grandDaughterDiffOfPhi(0,1,2,3)");
2227 
2228  // when no relations are set between the particles and the eclClusters, nan is expected to be returned for the Cluster- vars
2229  // no problems are supposed to happen for non-Cluster- vars
2230  // also, we expect NaN when we pass wrong indexes
2231  ASSERT_NE(var_ClusterPhi, nullptr);
2232  EXPECT_TRUE(std::isnan(var_ClusterPhi->function(parGranny_noclst)));
2233  EXPECT_TRUE(std::isnan(var_ClusterTheta->function(parGranny_noclst)));
2234  EXPECT_TRUE(std::isnan(var_ClusterE->function(parGranny_noclst)));
2235  EXPECT_FLOAT_EQ(var_Phi->function(parGranny_noclst), 0.32175055);
2236  EXPECT_FLOAT_EQ(var_Theta->function(parGranny_noclst), 0.06311664);
2237  EXPECT_FLOAT_EQ(var_E->function(parGranny_noclst), -0.46293831);
2238  EXPECT_TRUE(std::isnan(var_ClusterPhi_wrongIndexes->function(parGranny_noclst)));
2239  EXPECT_TRUE(std::isnan(var_Phi_wrongIndexes->function(parGranny_noclst)));
2240  EXPECT_TRUE(std::isnan(var_ClusterE_wrongIndexes->function(parGranny_noclst)));
2241  EXPECT_TRUE(std::isnan(var_E_wrongIndexes->function(parGranny_noclst)));
2242 
2243  // set relations between particles and eclClusters
2244  ECLCluster* eclst0 = eclclusters.appendNew(ECLCluster());
2245  eclst0->setEnergy(dau0_4vec.E());
2246  eclst0->setHypothesis(ECLCluster::EHypothesisBit::c_nPhotons);
2247  eclst0->setClusterId(1);
2248  eclst0->setTheta(dau0_4vec.Theta());
2249  eclst0->setPhi(dau0_4vec.Phi());
2250  eclst0->setR(148.4);
2251  ECLCluster* eclst1 = eclclusters.appendNew(ECLCluster());
2252  eclst1->setEnergy(dau1_4vec.E());
2253  eclst1->setHypothesis(ECLCluster::EHypothesisBit::c_nPhotons);
2254  eclst1->setClusterId(2);
2255  eclst1->setTheta(dau1_4vec.Theta());
2256  eclst1->setPhi(dau1_4vec.Phi());
2257  eclst1->setR(148.5);
2258 
2259  const Particle* newDaughter0 = particles.appendNew(Particle(eclclusters[0]));
2260  daughterIndices0.push_back(newDaughter0->getArrayIndex());
2261  const Particle* par0 = particles.appendNew(momentum_0, 111, Particle::c_Unflavored, daughterIndices0);
2262 
2263  const Particle* newDaughter1 = particles.appendNew(Particle(eclclusters[1]));
2264  daughterIndices1.push_back(newDaughter1->getArrayIndex());
2265  const Particle* par1 = particles.appendNew(momentum_1, 111, Particle::c_Unflavored, daughterIndices1);
2266 
2267  daughterIndices2.push_back(par0->getArrayIndex());
2268  daughterIndices2.push_back(par1->getArrayIndex());
2269  const Particle* parGranny = particles.appendNew(momentum, 111, Particle::c_Unflavored, daughterIndices2);
2270  //const Particle* par = particles.appendNew(momentum, 111, Particle::c_Unflavored, daughterIndices);
2271 
2272  //now we expect non-nan results
2273  EXPECT_FLOAT_EQ(var_ClusterPhi->function(parGranny), 0.32175055);
2274  EXPECT_FLOAT_EQ(var_Phi->function(parGranny), 0.32175055);
2275  EXPECT_FLOAT_EQ(var_ClusterTheta->function(parGranny), 0.06311664);
2276  EXPECT_FLOAT_EQ(var_Theta->function(parGranny), 0.06311664);
2277  EXPECT_FLOAT_EQ(var_ClusterE->function(parGranny), -0.46293831);
2278  EXPECT_FLOAT_EQ(var_E->function(parGranny), -0.46293813);
2279  }
2280 
2281  TEST_F(MetaVariableTest, daughterNormDiffOf)
2282  {
2283  TLorentzVector momentum;
2284  const int nDaughters = 4;
2285  StoreArray<Particle> particles;
2286  std::vector<int> daughterIndices;
2287  for (int i = 0; i < nDaughters; i++) {
2288  Particle d(TLorentzVector(1, 1, 1, i * 1.0 + 1.0), (i % 2) ? -11 : 211);
2289  momentum += d.get4Vector();
2290  Particle* newDaughters = particles.appendNew(d);
2291  daughterIndices.push_back(newDaughters->getArrayIndex());
2292  }
2293  const Particle* p = particles.appendNew(momentum, 411, Particle::c_Unflavored, daughterIndices);
2294 
2295  const Manager::Var* var = Manager::Instance().getVariable("daughterNormDiffOf(0, 1, PDG)");
2296  ASSERT_NE(var, nullptr);
2297  EXPECT_FLOAT_EQ(var->function(p), -222 / 200.);
2298 
2299  var = Manager::Instance().getVariable("daughterNormDiffOf(1, 0, PDG)");
2300  ASSERT_NE(var, nullptr);
2301  EXPECT_FLOAT_EQ(var->function(p), 222 / 200.);
2302 
2303  var = Manager::Instance().getVariable("daughterNormDiffOf(0, 1, abs(PDG))");
2304  ASSERT_NE(var, nullptr);
2305  EXPECT_FLOAT_EQ(var->function(p), -200 / 222.);
2306 
2307  var = Manager::Instance().getVariable("daughterNormDiffOf(1, 1, PDG)");
2308  ASSERT_NE(var, nullptr);
2309  EXPECT_FLOAT_EQ(var->function(p), -0 / 22.);
2310 
2311  var = Manager::Instance().getVariable("daughterNormDiffOf(1, 3, abs(PDG))");
2312  ASSERT_NE(var, nullptr);
2313  EXPECT_FLOAT_EQ(var->function(p), 0 / 22.);
2314 
2315  var = Manager::Instance().getVariable("daughterNormDiffOf(0, 2, PDG)");
2316  ASSERT_NE(var, nullptr);
2317  EXPECT_FLOAT_EQ(var->function(p), 0 / 422.);
2318 
2319  }
2320 
2321  TEST_F(MetaVariableTest, daughterMotherDiffOf)
2322  {
2323  TLorentzVector momentum;
2324  const int nDaughters = 4;
2325  StoreArray<Particle> particles;
2326  std::vector<int> daughterIndices;
2327  for (int i = 0; i < nDaughters; i++) {
2328  Particle d(TLorentzVector(1, 1, 1, i * 1.0 + 1.0), (i % 2) ? -11 : 211);
2329  momentum += d.get4Vector();
2330  Particle* newDaughters = particles.appendNew(d);
2331  daughterIndices.push_back(newDaughters->getArrayIndex());
2332  }
2333  const Particle* p = particles.appendNew(momentum, 411, Particle::c_Unflavored, daughterIndices);
2334 
2335  const Manager::Var* var = Manager::Instance().getVariable("daughterMotherDiffOf(1, PDG)");
2336  ASSERT_NE(var, nullptr);
2337  EXPECT_FLOAT_EQ(var->function(p), 422);
2338 
2339  var = Manager::Instance().getVariable("daughterMotherDiffOf(1, abs(PDG))");
2340  ASSERT_NE(var, nullptr);
2341  EXPECT_FLOAT_EQ(var->function(p), 400);
2342 
2343  var = Manager::Instance().getVariable("daughterMotherDiffOf(0, PDG)");
2344  ASSERT_NE(var, nullptr);
2345  EXPECT_FLOAT_EQ(var->function(p), 200);
2346 
2347  }
2348 
2349  TEST_F(MetaVariableTest, daughterMotherNormDiffOf)
2350  {
2351  TLorentzVector momentum;
2352  const int nDaughters = 4;
2353  StoreArray<Particle> particles;
2354  std::vector<int> daughterIndices;
2355  for (int i = 0; i < nDaughters; i++) {
2356  Particle d(TLorentzVector(1, 1, 1, i * 1.0 + 1.0), (i % 2) ? -11 : 211);
2357  momentum += d.get4Vector();
2358  Particle* newDaughters = particles.appendNew(d);
2359  daughterIndices.push_back(newDaughters->getArrayIndex());
2360  }
2361  const Particle* p = particles.appendNew(momentum, 411, Particle::c_Unflavored, daughterIndices);
2362 
2363  const Manager::Var* var = Manager::Instance().getVariable("daughterMotherNormDiffOf(1, PDG)");
2364  ASSERT_NE(var, nullptr);
2365  EXPECT_FLOAT_EQ(var->function(p), 422 / 400.);
2366 
2367  var = Manager::Instance().getVariable("daughterMotherNormDiffOf(1, abs(PDG))");
2368  ASSERT_NE(var, nullptr);
2369  EXPECT_FLOAT_EQ(var->function(p), 400 / 422.);
2370 
2371  var = Manager::Instance().getVariable("daughterMotherNormDiffOf(0, PDG)");
2372  ASSERT_NE(var, nullptr);
2373  EXPECT_FLOAT_EQ(var->function(p), 200 / 622.);
2374 
2375  }
2376 
2377  TEST_F(MetaVariableTest, constant)
2378  {
2379 
2380  const Manager::Var* var = Manager::Instance().getVariable("constant(1)");
2381  ASSERT_NE(var, nullptr);
2382  EXPECT_FLOAT_EQ(var->function(nullptr), 1.0);
2383 
2384  var = Manager::Instance().getVariable("constant(0)");
2385  ASSERT_NE(var, nullptr);
2386  EXPECT_FLOAT_EQ(var->function(nullptr), 0.0);
2387 
2388  }
2389 
2390  TEST_F(MetaVariableTest, abs)
2391  {
2392  Particle p({ 0.1 , -0.4, 0.8, 2.0 }, 11);
2393  Particle p2({ -0.1 , -0.4, 0.8, 4.0 }, -11);
2394 
2395  const Manager::Var* var = Manager::Instance().getVariable("abs(px)");
2396  ASSERT_NE(var, nullptr);
2397  EXPECT_FLOAT_EQ(var->function(&p), 0.1);
2398  EXPECT_FLOAT_EQ(var->function(&p2), 0.1);
2399 
2400  }
2401 
2402  TEST_F(MetaVariableTest, sin)
2403  {
2404  Particle p({ 3.14159265359 / 2.0 , -0.4, 0.8, 1.0}, 11);
2405  Particle p2({ 0.0 , -0.4, 0.8, 1.0 }, -11);
2406 
2407  const Manager::Var* var = Manager::Instance().getVariable("sin(px)");
2408  ASSERT_NE(var, nullptr);
2409  EXPECT_FLOAT_EQ(var->function(&p), 1.0);
2410  EXPECT_ALL_NEAR(var->function(&p2), 0.0, 1e-6);
2411 
2412  }
2413 
2414  TEST_F(MetaVariableTest, cos)
2415  {
2416  Particle p({ 3.14159265359 / 2.0 , -0.4, 0.8, 1.0}, 11);
2417  Particle p2({ 0.0 , -0.4, 0.8, 1.0 }, -11);
2418 
2419  const Manager::Var* var = Manager::Instance().getVariable("cos(px)");
2420  ASSERT_NE(var, nullptr);
2421  EXPECT_ALL_NEAR(var->function(&p), 0.0, 1e-6);
2422  EXPECT_FLOAT_EQ(var->function(&p2), 1.0);
2423 
2424  }
2425 
2426  TEST_F(MetaVariableTest, NBDeltaIfMissingDeathTest)
2427  {
2428  //Variable got removed, test for absence
2429  EXPECT_B2FATAL(Manager::Instance().getVariable("NBDeltaIfMissing(TOP, 11)"));
2430  EXPECT_B2FATAL(Manager::Instance().getVariable("NBDeltaIfMissing(ARICH, 11)"));
2431  }
2432 
2433  TEST_F(MetaVariableTest, matchedMC)
2434  {
2435  DataStore::Instance().setInitializeActive(true);
2436  StoreArray<MCParticle> mcParticles;
2437  StoreArray<Particle> particles;
2438  particles.registerRelationTo(mcParticles);
2439  DataStore::Instance().setInitializeActive(false);
2440 
2441  auto* mcParticle = mcParticles.appendNew();
2442  mcParticle->setPDG(11);
2443  mcParticle->setStatus(MCParticle::c_PrimaryParticle);
2444  auto* p1 = particles.appendNew(TLorentzVector({ 0.0 , -0.4, 0.8, 1.0}), 11);
2445  p1->addRelationTo(mcParticle);
2446 
2447  mcParticle = mcParticles.appendNew();
2448  mcParticle->setPDG(-11);
2449  mcParticle->setStatus(MCParticle::c_PrimaryParticle);
2450  auto* p2 = particles.appendNew(TLorentzVector({ 0.0 , -0.4, 0.8, 1.0}), 11);
2451  p2->addRelationTo(mcParticle);
2452 
2453  mcParticle = mcParticles.appendNew();
2454  mcParticle->setPDG(22);
2455  mcParticle->setStatus(MCParticle::c_PrimaryParticle);
2456  auto* p3 = particles.appendNew(TLorentzVector({ 0.0 , -0.4, 0.8, 1.0}), 11);
2457  p3->addRelationTo(mcParticle);
2458 
2459  // Test if matchedMC also works for particle which already is an MCParticle.
2460  auto* p4 = particles.appendNew(mcParticle);
2461 
2462  const Manager::Var* var = Manager::Instance().getVariable("matchedMC(charge)");
2463  ASSERT_NE(var, nullptr);
2464  EXPECT_FLOAT_EQ(var->function(p1), -1);
2465  EXPECT_FLOAT_EQ(var->function(p2), 1);
2466  EXPECT_FLOAT_EQ(var->function(p3), 0);
2467  EXPECT_FLOAT_EQ(var->function(p4), 0);
2468  }
2469 
2470  TEST_F(MetaVariableTest, countInList)
2471  {
2472  StoreArray<Particle> particles;
2473  DataStore::EStoreFlags flags = DataStore::c_DontWriteOut;
2474 
2475  StoreObjPtr<ParticleList> outputList("pList1");
2476  DataStore::Instance().setInitializeActive(true);
2477  outputList.registerInDataStore(flags);
2478  DataStore::Instance().setInitializeActive(false);
2479  outputList.create();
2480  outputList->initialize(22, "pList1");
2481 
2482  particles.appendNew(Particle({0.5 , 0.4 , 0.5 , 0.8}, 22, Particle::c_Unflavored, Particle::c_Undefined, 2));
2483  particles.appendNew(Particle({0.5 , 0.2 , 0.7 , 0.9}, 22, Particle::c_Unflavored, Particle::c_Undefined, 3));
2484  particles.appendNew(Particle({0.4 , 0.2 , 0.7 , 0.9}, 22, Particle::c_Unflavored, Particle::c_Undefined, 4));
2485  particles.appendNew(Particle({0.5 , 0.4 , 0.8 , 1.1}, 22, Particle::c_Unflavored, Particle::c_Undefined, 5));
2486  particles.appendNew(Particle({0.3 , 0.3 , 0.4 , 0.6}, 22, Particle::c_Unflavored, Particle::c_Undefined, 6));
2487 
2488  outputList->addParticle(0, 22, Particle::c_Unflavored);
2489  outputList->addParticle(1, 22, Particle::c_Unflavored);
2490  outputList->addParticle(2, 22, Particle::c_Unflavored);
2491  outputList->addParticle(3, 22, Particle::c_Unflavored);
2492  outputList->addParticle(4, 22, Particle::c_Unflavored);
2493 
2494  const Manager::Var* var = Manager::Instance().getVariable("countInList(pList1, E < 0.85)");
2495  ASSERT_NE(var, nullptr);
2496  EXPECT_DOUBLE_EQ(var->function(nullptr), 2);
2497 
2498  var = Manager::Instance().getVariable("countInList(pList1)");
2499  ASSERT_NE(var, nullptr);
2500  EXPECT_DOUBLE_EQ(var->function(nullptr), 5);
2501 
2502  var = Manager::Instance().getVariable("countInList(pList1, E > 5)");
2503  ASSERT_NE(var, nullptr);
2504  EXPECT_DOUBLE_EQ(var->function(nullptr), 0);
2505 
2506  var = Manager::Instance().getVariable("countInList(pList1, E < 5)");
2507  ASSERT_NE(var, nullptr);
2508  EXPECT_DOUBLE_EQ(var->function(nullptr), 5);
2509  }
2510 
2511  TEST_F(MetaVariableTest, isInList)
2512  {
2513  // we need the particles StoreArray
2514  StoreArray<Particle> particles;
2515  DataStore::EStoreFlags flags = DataStore::c_DontWriteOut;
2516 
2517  // create a photon list for testing
2518  StoreObjPtr<ParticleList> gammalist("testGammaList");
2519  DataStore::Instance().setInitializeActive(true);
2520  gammalist.registerInDataStore(flags);
2521  DataStore::Instance().setInitializeActive(false);
2522  gammalist.create();
2523  gammalist->initialize(22, "testGammaList");
2524 
2525  // mock up two photons
2526  Particle goingin({0.5 , 0.4 , 0.5 , 0.8}, 22, Particle::c_Unflavored, Particle::c_Undefined, 0);
2527  Particle notgoingin({0.3 , 0.3 , 0.4 , 0.6}, 22, Particle::c_Unflavored, Particle::c_Undefined, 1);
2528  auto* inthelist = particles.appendNew(goingin);
2529  auto* notinthelist = particles.appendNew(notgoingin);
2530 
2531  // put the the zeroth one in the list the first on not in the list
2532  gammalist->addParticle(0, 22, Particle::c_Unflavored);
2533 
2534  // get the variables
2535  const Manager::Var* vnonsense = Manager::Instance().getVariable("isInList(NONEXISTANTLIST)");
2536  const Manager::Var* vsensible = Manager::Instance().getVariable("isInList(testGammaList)");
2537 
2538  // -
2539  EXPECT_B2FATAL(vnonsense->function(notinthelist));
2540  EXPECT_FLOAT_EQ(vsensible->function(inthelist), 1.0);
2541  EXPECT_FLOAT_EQ(vsensible->function(notinthelist), 0.0);
2542  }
2543 
2544  TEST_F(MetaVariableTest, sourceObjectIsInList)
2545  {
2546  // datastore things
2547  DataStore::Instance().reset();
2548  DataStore::Instance().setInitializeActive(true);
2549 
2550  // needed to mock up
2551  StoreArray<ECLCluster> clusters;
2552  StoreArray<Particle> particles;
2553  StoreObjPtr<ParticleList> gammalist("testGammaList");
2554 
2555  clusters.registerInDataStore();
2556  particles.registerInDataStore();
2557  DataStore::EStoreFlags flags = DataStore::c_DontWriteOut;
2558  gammalist.registerInDataStore(flags);
2559 
2560  // end datastore things
2561  DataStore::Instance().setInitializeActive(false);
2562 
2563  // of course we have to create the list...
2564  gammalist.create();
2565  gammalist->initialize(22, "testGammaList");
2566 
2567  // mock up two clusters from the ECL let's say they both came from true Klongs
2568  // but one looked a little bit photon-like
2569  auto* cl0 = clusters.appendNew(ECLCluster());
2570  cl0->setEnergy(1.0);
2571  cl0->setHypothesis(ECLCluster::EHypothesisBit::c_nPhotons);
2572  cl0->addHypothesis(ECLCluster::EHypothesisBit::c_neutralHadron);
2573  cl0->setClusterId(0);
2574  auto* cl1 = clusters.appendNew(ECLCluster());
2575  cl1->setEnergy(1.0);
2576  cl1->setHypothesis(ECLCluster::EHypothesisBit::c_neutralHadron);
2577  cl1->setClusterId(1);
2578 
2579  // create particles from the clusters
2580  Particle myphoton(cl0, Const::photon);
2581  Particle iscopiedin(cl0, Const::Klong);
2582  Particle notcopiedin(cl1, Const::Klong);
2583 
2584  // add the particle created from cluster zero to the gamma list
2585  auto* myphoton_ = particles.appendNew(myphoton);
2586  gammalist->addParticle(myphoton_);
2587 
2588  auto* iscopied = particles.appendNew(iscopiedin); // a clone of this guy is now in the gamma list
2589  auto* notcopied = particles.appendNew(notcopiedin);
2590 
2591  // get the variables
2592  const Manager::Var* vnonsense = Manager::Instance().getVariable("sourceObjectIsInList(NONEXISTANTLIST)");
2593  const Manager::Var* vsensible = Manager::Instance().getVariable("sourceObjectIsInList(testGammaList)");
2594 
2595  // -
2596  EXPECT_B2FATAL(vnonsense->function(iscopied));
2597  EXPECT_FLOAT_EQ(vsensible->function(iscopied), 1.0);
2598  EXPECT_FLOAT_EQ(vsensible->function(notcopied), 0.0);
2599 
2600  // now mock up some other type particles
2601  Particle composite({0.5 , 0.4 , 0.5 , 0.8}, 512, Particle::c_Unflavored, Particle::c_Composite, 0);
2602  Particle undefined({0.3 , 0.3 , 0.4 , 0.6}, 22, Particle::c_Unflavored, Particle::c_Undefined, 1);
2603  auto* composite_ = particles.appendNew(undefined);
2604  auto* undefined_ = particles.appendNew(composite);
2605  EXPECT_FLOAT_EQ(vsensible->function(composite_), -1.0);
2606  EXPECT_FLOAT_EQ(vsensible->function(undefined_), -1.0);
2607  }
2608 
2609  TEST_F(MetaVariableTest, mcParticleIsInMCList)
2610  {
2611  // datastore things
2612  DataStore::Instance().reset();
2613  DataStore::Instance().setInitializeActive(true);
2614 
2615  // needed to mock up
2616  StoreArray<MCParticle> mcparticles;
2617  StoreArray<Particle> particles;
2618  StoreObjPtr<ParticleList> list("testList");
2619  StoreObjPtr<ParticleList> anotherlist("supplimentaryList");
2620 
2621  mcparticles.registerInDataStore();
2622  particles.registerInDataStore();
2623  particles.registerRelationTo(mcparticles);
2624  DataStore::EStoreFlags flags = DataStore::c_DontWriteOut;
2625  list.registerInDataStore(flags);
2626  anotherlist.registerInDataStore(flags);
2627 
2628  DataStore::Instance().setInitializeActive(false);
2629  // end datastore setup
2630 
2631  list.create();
2632  list->initialize(22, "testList");
2633 
2634  anotherlist.create();
2635  anotherlist->initialize(22, "supplimentaryList");
2636 
2637  // MCParticles
2638  auto* mcphoton = mcparticles.appendNew();
2639  mcphoton->setPDG(22);
2640  mcphoton->setStatus(MCParticle::c_PrimaryParticle);
2641 
2642  auto* mcelectron = mcparticles.appendNew();
2643  mcelectron->setPDG(11);
2644  mcelectron->setStatus(MCParticle::c_PrimaryParticle);
2645 
2646  auto* mcanotherelectron = mcparticles.appendNew();
2647  mcanotherelectron->setPDG(22);
2648  mcanotherelectron->setStatus(MCParticle::c_PrimaryParticle);
2649 
2650  auto* mcyetanotherelectron = mcparticles.appendNew();
2651  mcyetanotherelectron->setPDG(22);
2652  mcyetanotherelectron->setStatus(MCParticle::c_PrimaryParticle);
2653 
2654  // particles
2655  auto* photon = particles.appendNew(TLorentzVector({ 0.0 , -0.4, 0.8, 1.0}), 22);
2656  photon->addRelationTo(mcphoton);
2657  list->addParticle(photon);
2658 
2659  auto* electron = particles.appendNew(TLorentzVector({ 0.0 , -0.4, 0.8, 1.0}), 22);
2660  electron->addRelationTo(mcelectron);
2661  list->addParticle(electron);
2662 
2663  auto* other = particles.appendNew(TLorentzVector({ 0.0 , -0.4, 0.8, 1.0}), 22);
2664  other->addRelationTo(mcanotherelectron);
2665 
2666  auto* yetanotherelectron = particles.appendNew(TLorentzVector({ 0.0 , -0.4, 0.8, 1.0}), 22);
2667  yetanotherelectron->addRelationTo(mcyetanotherelectron);
2668  anotherlist->addParticle(yetanotherelectron);
2669  // not in the list
2670 
2671  // get the variable
2672  const Manager::Var* vnonsense = Manager::Instance().getVariable("mcParticleIsInMCList(NONEXISTANTLIST)");
2673  const Manager::Var* vsensible = Manager::Instance().getVariable("mcParticleIsInMCList(testList)");
2674 
2675  // -
2676  EXPECT_B2FATAL(vnonsense->function(photon));
2677  EXPECT_FLOAT_EQ(vsensible->function(photon), 1.0);
2678  EXPECT_FLOAT_EQ(vsensible->function(electron), 1.0);
2679  EXPECT_FLOAT_EQ(vsensible->function(other), 0.0);
2680  EXPECT_FLOAT_EQ(vsensible->function(yetanotherelectron), 0.0);
2681 
2682  // now mock up some other type particles
2683  Particle composite({0.5 , 0.4 , 0.5 , 0.8}, 512, Particle::c_Unflavored, Particle::c_Composite, 0);
2684  Particle undefined({0.3 , 0.3 , 0.4 , 0.6}, 22, Particle::c_Unflavored, Particle::c_Undefined, 1);
2685  auto* composite_ = particles.appendNew(undefined);
2686  auto* undefined_ = particles.appendNew(composite);
2687  EXPECT_FLOAT_EQ(vsensible->function(composite_), 0.0);
2688  EXPECT_FLOAT_EQ(vsensible->function(undefined_), 0.0);
2689  }
2690 
2691  TEST_F(MetaVariableTest, mostB2BAndClosestParticles)
2692  {
2693  /* Mock up an event with a "photon" and an "electron" which are nearly back to
2694  * back, and second "photon" which is close-ish to the "electron".
2695  *
2696  * Other test of non-existent / empty lists and variables also included.
2697  */
2698 
2699  // Connect gearbox for CMS variables
2700  Gearbox& gearbox = Gearbox::getInstance();
2701  gearbox.setBackends({std::string("file:")});
2702  gearbox.close();
2703  gearbox.open("geometry/Belle2.xml", false);
2704 
2705  // we need the particles StoreArray
2706  StoreArray<Particle> particles;
2707  DataStore::EStoreFlags flags = DataStore::c_DontWriteOut;
2708 
2709  // create a photon list for testing
2710  StoreObjPtr<ParticleList> gammalist("testGammaList");
2711  StoreObjPtr<ParticleList> emptylist("testEmptyList");
2712  DataStore::Instance().setInitializeActive(true);
2713  gammalist.registerInDataStore(flags);
2714  emptylist.registerInDataStore(flags);
2715  DataStore::Instance().setInitializeActive(false);
2716  gammalist.create();
2717  gammalist->initialize(22, "testGammaList");
2718  emptylist.create();
2719  emptylist->initialize(22, "testEmptyList");
2720 
2721  // create some photons in an stdvector
2722  std::vector<Particle> gammavector = {
2723  Particle({ -1.0 , -1.0 , 0.8, 1.2}, // this should be the most b2b to our reference particle
2724  22, Particle::c_Unflavored, Particle::c_Undefined, 0),
2725  Particle({0.2 , 0.7 , 0.9, 3.4}, // should be the closest
2726  22, Particle::c_Unflavored, Particle::c_Undefined, 1),
2727  };
2728  // put the photons in the StoreArray
2729  for (const auto& g : gammavector)
2730  particles.appendNew(g);
2731 
2732  // put the photons in the test list
2733  for (size_t i = 0; i < gammavector.size(); i++)
2734  gammalist->addParticle(i, 22, Particle::c_Unflavored);
2735 
2736  // add the reference particle (electron) to the StoreArray
2737  const auto* electron = particles.appendNew(
2738  Particle({1.0 , 1.0 , 0.5, 0.8}, // somewhere in the +ve quarter of the detector
2739  11, Particle::c_Unflavored, Particle::c_Undefined, 2) // needs to be incremented if we add to gamma vector
2740  );
2741 
2742  {
2743  EXPECT_B2FATAL(Manager::Instance().getVariable("angleToClosestInList"));
2744  EXPECT_B2FATAL(Manager::Instance().getVariable("angleToClosestInList(A, B)"));
2745 
2746  const auto* nonexistant = Manager::Instance().getVariable("angleToClosestInList(NONEXISTANTLIST)");
2747  EXPECT_B2FATAL(nonexistant->function(electron));
2748 
2749  const auto* empty = Manager::Instance().getVariable("angleToClosestInList(testEmptyList)");
2750  EXPECT_TRUE(std::isnan(empty->function(electron)));
2751 
2752  const auto* closest = Manager::Instance().getVariable("angleToClosestInList(testGammaList)");
2753  EXPECT_FLOAT_EQ(closest->function(electron), 0.68014491);
2754 
2755  const auto* closestCMS = Manager::Instance().getVariable("useCMSFrame(angleToClosestInList(testGammaList))");
2756  EXPECT_FLOAT_EQ(closestCMS->function(electron), 0.72592634);
2757  }
2758 
2759  {
2760  EXPECT_B2FATAL(Manager::Instance().getVariable("closestInList"));
2761  EXPECT_B2FATAL(Manager::Instance().getVariable("closestInList(A, B, C)"));
2762 
2763  const auto* nonexistant = Manager::Instance().getVariable("closestInList(NONEXISTANTLIST, E)");
2764  EXPECT_B2FATAL(nonexistant->function(electron));
2765 
2766  const auto* empty = Manager::Instance().getVariable("closestInList(testEmptyList, E)");
2767  EXPECT_TRUE(std::isnan(empty->function(electron)));
2768 
2769  const auto* closest = Manager::Instance().getVariable("closestInList(testGammaList, E)");
2770  EXPECT_FLOAT_EQ(closest->function(electron), 3.4);
2771 
2772  const auto* closestCMS = Manager::Instance().getVariable("useCMSFrame(closestInList(testGammaList, E))");
2773  EXPECT_FLOAT_EQ(closestCMS->function(electron), 3.2732551); // the energy gets smeared because of boost
2774 
2775  const auto* closestCMSLabE = Manager::Instance().getVariable("useCMSFrame(closestInList(testGammaList, useLabFrame(E)))");
2776  EXPECT_FLOAT_EQ(closestCMSLabE->function(electron), 3.4); // aaand should be back to the lab frame value
2777  }
2778 
2779  {
2780  EXPECT_B2FATAL(Manager::Instance().getVariable("angleToMostB2BInList"));
2781  EXPECT_B2FATAL(Manager::Instance().getVariable("angleToMostB2BInList(A, B)"));
2782 
2783  const auto* nonexistant = Manager::Instance().getVariable("angleToMostB2BInList(NONEXISTANTLIST)");
2784  EXPECT_B2FATAL(nonexistant->function(electron));
2785 
2786  const auto* empty = Manager::Instance().getVariable("angleToMostB2BInList(testEmptyList)");
2787  EXPECT_TRUE(std::isnan(empty->function(electron)));
2788 
2789  const auto* mostB2B = Manager::Instance().getVariable("angleToMostB2BInList(testGammaList)");
2790  EXPECT_FLOAT_EQ(mostB2B->function(electron), 2.2869499);
2791 
2792  const auto* mostB2BCMS = Manager::Instance().getVariable("useCMSFrame(angleToMostB2BInList(testGammaList))");
2793  EXPECT_FLOAT_EQ(mostB2BCMS->function(electron), 2.6054888);
2794  }
2795 
2796  {
2797  EXPECT_B2FATAL(Manager::Instance().getVariable("mostB2BInList"));
2798  EXPECT_B2FATAL(Manager::Instance().getVariable("mostB2BInList(A, B, C)"));
2799 
2800  const auto* nonexistant = Manager::Instance().getVariable("mostB2BInList(NONEXISTANTLIST, E)");
2801  EXPECT_B2FATAL(nonexistant->function(electron));
2802 
2803  const auto* empty = Manager::Instance().getVariable("mostB2BInList(testEmptyList, E)");
2804  EXPECT_TRUE(std::isnan(empty->function(electron)));
2805 
2806  const auto* mostB2B = Manager::Instance().getVariable("mostB2BInList(testGammaList, E)");
2807  EXPECT_FLOAT_EQ(mostB2B->function(electron), 1.2);
2808 
2809  const auto* mostB2BCMS = Manager::Instance().getVariable("useCMSFrame(mostB2BInList(testGammaList, E))");
2810  EXPECT_FLOAT_EQ(mostB2BCMS->function(electron), 1.0647389); // the energy gets smeared because of boost
2811 
2812  const auto* mostB2BCMSLabE = Manager::Instance().getVariable("useCMSFrame(mostB2BInList(testGammaList, useLabFrame(E)))");
2813  EXPECT_FLOAT_EQ(mostB2BCMSLabE->function(electron), 1.2); // aaand should be back to the lab frame value
2814  }
2815  }
2816 
2817  TEST_F(MetaVariableTest, totalEnergyOfParticlesInList)
2818  {
2819  // we need the particles StoreArray
2820  StoreArray<Particle> particles;
2821  DataStore::EStoreFlags flags = DataStore::c_DontWriteOut;
2822 
2823  // create a photon list for testing
2824  StoreObjPtr<ParticleList> gammalist("testGammaList");
2825  DataStore::Instance().setInitializeActive(true);
2826  gammalist.registerInDataStore(flags);
2827  DataStore::Instance().setInitializeActive(false);
2828  gammalist.create();
2829  gammalist->initialize(22, "testGammaList");
2830 
2831  // create some photons in an stdvector
2832  std::vector<Particle> gammavector = {
2833  Particle({0.5 , 0.4 , 0.5 , 0.8}, 22, Particle::c_Unflavored, Particle::c_Undefined, 0),
2834  Particle({0.5 , 0.2 , 0.7 , 0.9}, 22, Particle::c_Unflavored, Particle::c_Undefined, 1),
2835  Particle({0.4 , 0.2 , 0.7 , 0.9}, 22, Particle::c_Unflavored, Particle::c_Undefined, 2),
2836  Particle({0.5 , 0.4 , 0.8 , 1.1}, 22, Particle::c_Unflavored, Particle::c_Undefined, 3),
2837  Particle({0.3 , 0.3 , 0.4 , 0.6}, 22, Particle::c_Unflavored, Particle::c_Undefined, 4)
2838  };
2839 
2840  // put the photons in the StoreArray
2841  for (const auto& g : gammavector)
2842  particles.appendNew(g);
2843 
2844  // put the photons in the test list
2845  for (size_t i = 0; i < gammavector.size(); i++)
2846  gammalist->addParticle(i, 22, Particle::c_Unflavored);
2847 
2848  // get their total energy
2849  const Manager::Var* vnonsense = Manager::Instance().getVariable(
2850  "totalEnergyOfParticlesInList(NONEXISTANTLIST)");
2851  const Manager::Var* vsensible = Manager::Instance().getVariable(
2852  "totalEnergyOfParticlesInList(testGammaList)");
2853 
2854  // -
2855  EXPECT_B2FATAL(vnonsense->function(nullptr));
2856  EXPECT_FLOAT_EQ(vsensible->function(nullptr), 4.3);
2857  }
2858  TEST_F(MetaVariableTest, totalPxOfParticlesInList)
2859  {
2860  // we need the particles StoreArray
2861  StoreArray<Particle> particles;
2862  DataStore::EStoreFlags flags = DataStore::c_DontWriteOut;
2863 
2864  // create a photon list for testing
2865  StoreObjPtr<ParticleList> gammalist("testGammaList");
2866  DataStore::Instance().setInitializeActive(true);
2867  gammalist.registerInDataStore(flags);
2868  DataStore::Instance().setInitializeActive(false);
2869  gammalist.create();
2870  gammalist->initialize(22, "testGammaList");
2871 
2872  // create some photons in an stdvector
2873  std::vector<Particle> gammavector = {
2874  Particle({0.5 , 0.4 , 0.5 , 0.8}, 22, Particle::c_Unflavored, Particle::c_Undefined, 0),
2875  Particle({0.5 , 0.2 , 0.7 , 0.9}, 22, Particle::c_Unflavored, Particle::c_Undefined, 1),
2876  Particle({0.4 , 0.2 , 0.7 , 0.9}, 22, Particle::c_Unflavored, Particle::c_Undefined, 2),
2877  Particle({0.5 , 0.4 , 0.8 , 1.1}, 22, Particle::c_Unflavored, Particle::c_Undefined, 3),
2878  Particle({0.3 , 0.3 , 0.4 , 0.6}, 22, Particle::c_Unflavored, Particle::c_Undefined, 4)
2879  };
2880 
2881  // put the photons in the StoreArray
2882  for (const auto& g : gammavector)
2883  particles.appendNew(g);
2884 
2885  // put the photons in the test list
2886  for (size_t i = 0; i < gammavector.size(); i++)
2887  gammalist->addParticle(i, 22, Particle::c_Unflavored);
2888 
2889  // get their total energy
2890  const Manager::Var* vnonsense = Manager::Instance().getVariable(
2891  "totalPxOfParticlesInList(NONEXISTANTLIST)");
2892  const Manager::Var* vsensible = Manager::Instance().getVariable(
2893  "totalPxOfParticlesInList(testGammaList)");
2894 
2895  // -
2896  EXPECT_B2FATAL(vnonsense->function(nullptr));
2897  EXPECT_FLOAT_EQ(vsensible->function(nullptr), 2.2);
2898  }
2899  TEST_F(MetaVariableTest, totalPyOfParticlesInList)
2900  {
2901  // we need the particles StoreArray
2902  StoreArray<Particle> particles;
2903  DataStore::EStoreFlags flags = DataStore::c_DontWriteOut;
2904 
2905  // create a photon list for testing
2906  StoreObjPtr<ParticleList> gammalist("testGammaList");
2907  DataStore::Instance().setInitializeActive(true);
2908  gammalist.registerInDataStore(flags);
2909  DataStore::Instance().setInitializeActive(false);
2910  gammalist.create();
2911  gammalist->initialize(22, "testGammaList");
2912 
2913  // create some photons in an stdvector
2914  std::vector<Particle> gammavector = {
2915  Particle({0.5 , 0.4 , 0.5 , 0.8}, 22, Particle::c_Unflavored, Particle::c_Undefined, 0),
2916  Particle({0.5 , 0.2 , 0.7 , 0.9}, 22, Particle::c_Unflavored, Particle::c_Undefined, 1),
2917  Particle({0.4 , 0.2 , 0.7 , 0.9}, 22, Particle::c_Unflavored, Particle::c_Undefined, 2),
2918  Particle({0.5 , 0.4 , 0.8 , 1.1}, 22, Particle::c_Unflavored, Particle::c_Undefined, 3),
2919  Particle({0.3 , 0.3 , 0.4 , 0.6}, 22, Particle::c_Unflavored, Particle::c_Undefined, 4)
2920  };
2921 
2922  // put the photons in the StoreArray
2923  for (const auto& g : gammavector)
2924  particles.appendNew(g);
2925 
2926  // put the photons in the test list
2927  for (size_t i = 0; i < gammavector.size(); i++)
2928  gammalist->addParticle(i, 22, Particle::c_Unflavored);
2929 
2930  // get their total energy
2931  const Manager::Var* vnonsense = Manager::Instance().getVariable(
2932  "totalPyOfParticlesInList(NONEXISTANTLIST)");
2933  const Manager::Var* vsensible = Manager::Instance().getVariable(
2934  "totalPyOfParticlesInList(testGammaList)");
2935 
2936  // -
2937  EXPECT_B2FATAL(vnonsense->function(nullptr));
2938  EXPECT_FLOAT_EQ(vsensible->function(nullptr), 1.5);
2939  }
2940  TEST_F(MetaVariableTest, totalPzOfParticlesInList)
2941  {
2942  // we need the particles StoreArray
2943  StoreArray<Particle> particles;
2944  DataStore::EStoreFlags flags = DataStore::c_DontWriteOut;
2945 
2946  // create a photon list for testing
2947  StoreObjPtr<ParticleList> gammalist("testGammaList");
2948  DataStore::Instance().setInitializeActive(true);
2949  gammalist.registerInDataStore(flags);
2950  DataStore::Instance().setInitializeActive(false);
2951  gammalist.create();
2952  gammalist->initialize(22, "testGammaList");
2953 
2954  // create some photons in an stdvector
2955  std::vector<Particle> gammavector = {
2956  Particle({0.5 , 0.4 , 0.5 , 0.8}, 22, Particle::c_Unflavored, Particle::c_Undefined, 0),
2957  Particle({0.5 , 0.2 , 0.7 , 0.9}, 22, Particle::c_Unflavored, Particle::c_Undefined, 1),
2958  Particle({0.4 , 0.2 , 0.7 , 0.9}, 22, Particle::c_Unflavored, Particle::c_Undefined, 2),
2959  Particle({0.5 , 0.4 , 0.8 , 1.1}, 22, Particle::c_Unflavored, Particle::c_Undefined, 3),
2960  Particle({0.3 , 0.3 , 0.4 , 0.6}, 22, Particle::c_Unflavored, Particle::c_Undefined, 4)
2961  };
2962 
2963  // put the photons in the StoreArray
2964  for (const auto& g : gammavector)
2965  particles.appendNew(g);
2966 
2967  // put the photons in the test list
2968  for (size_t i = 0; i < gammavector.size(); i++)
2969  gammalist->addParticle(i, 22, Particle::c_Unflavored);
2970 
2971  // get their total energy
2972  const Manager::Var* vnonsense = Manager::Instance().getVariable(
2973  "totalPzOfParticlesInList(NONEXISTANTLIST)");
2974  const Manager::Var* vsensible = Manager::Instance().getVariable(
2975  "totalPzOfParticlesInList(testGammaList)");
2976 
2977  // -
2978  EXPECT_B2FATAL(vnonsense->function(nullptr));
2979  EXPECT_FLOAT_EQ(vsensible->function(nullptr), 3.1);
2980  }
2981  TEST_F(MetaVariableTest, maxPtInList)
2982  {
2983  // we need the particles StoreArray
2984  StoreArray<Particle> particles;
2985  DataStore::EStoreFlags flags = DataStore::c_DontWriteOut;
2986 
2987  // create a photon list for testing
2988  StoreObjPtr<ParticleList> gammalist("testGammaList");
2989  DataStore::Instance().setInitializeActive(true);
2990  gammalist.registerInDataStore(flags);
2991  DataStore::Instance().setInitializeActive(false);
2992  gammalist.create();
2993  gammalist->initialize(22, "testGammaList");
2994 
2995  // create some photons in an stdvector
2996  std::vector<Particle> gammavector = {
2997  Particle({0.5 , 0.4 , 0.5 , 0.8}, 22, Particle::c_Unflavored, Particle::c_Undefined, 0),
2998  Particle({0.5 , 0.2 , 0.7 , 0.9}, 22, Particle::c_Unflavored, Particle::c_Undefined, 1),
2999  Particle({0.4 , 0.2 , 0.7 , 0.9}, 22, Particle::c_Unflavored, Particle::c_Undefined, 2),
3000  Particle({0.5 , 0.4 , 0.8 , 1.1}, 22, Particle::c_Unflavored, Particle::c_Undefined, 3),
3001  Particle({0.3 , 0.3 , 0.4 , 0.6}, 22, Particle::c_Unflavored, Particle::c_Undefined, 4)
3002  };
3003 
3004  // put the photons in the StoreArray
3005  for (const auto& g : gammavector)
3006  particles.appendNew(g);
3007 
3008  // put the photons in the test list
3009  for (size_t i = 0; i < gammavector.size(); i++)
3010  gammalist->addParticle(i, 22, Particle::c_Unflavored);
3011 
3012  // get their total energy
3013  const Manager::Var* vnonsense = Manager::Instance().getVariable(
3014  "maxPtInList(NONEXISTANTLIST)");
3015  const Manager::Var* vsensible = Manager::Instance().getVariable(
3016  "maxPtInList(testGammaList)");
3017 
3018  // -
3019  EXPECT_B2FATAL(vnonsense->function(nullptr));
3020  EXPECT_FLOAT_EQ(vsensible->function(nullptr), sqrt(0.5 * 0.5 + 0.4 * 0.4));
3021  }
3022 
3023 
3024  TEST_F(MetaVariableTest, numberOfNonOverlappingParticles)
3025  {
3026  StoreArray<Particle> particles;
3027  DataStore::EStoreFlags flags = DataStore::c_DontWriteOut;
3028 
3029  StoreObjPtr<ParticleList> outputList("pList1");
3030  DataStore::Instance().setInitializeActive(true);
3031  outputList.registerInDataStore(flags);
3032  DataStore::Instance().setInitializeActive(false);
3033  outputList.create();
3034  outputList->initialize(22, "pList1");
3035 
3036  auto* p1 = particles.appendNew(Particle({0.5 , 0.4 , 0.5 , 0.8}, 22, Particle::c_Unflavored, Particle::c_Undefined, 2));
3037  auto* p2 = particles.appendNew(Particle({0.5 , 0.2 , 0.7 , 0.9}, 22, Particle::c_Unflavored, Particle::c_Undefined, 3));
3038  auto* p3 = particles.appendNew(Particle({0.5 , 0.2 , 0.7 , 0.9}, 22, Particle::c_Unflavored, Particle::c_Undefined, 4));
3039 
3040  outputList->addParticle(0, 22, Particle::c_Unflavored);
3041  outputList->addParticle(1, 22, Particle::c_Unflavored);
3042 
3043  const Manager::Var* var = Manager::Instance().getVariable("numberOfNonOverlappingParticles(pList1)");
3044  ASSERT_NE(var, nullptr);
3045  EXPECT_DOUBLE_EQ(var->function(p1), 1);
3046  EXPECT_DOUBLE_EQ(var->function(p2), 1);
3047  EXPECT_DOUBLE_EQ(var->function(p3), 2);
3048 
3049  }
3050 
3051  TEST_F(MetaVariableTest, veto)
3052  {
3053  StoreArray<Particle> particles;
3054  DataStore::EStoreFlags flags = DataStore::c_DontWriteOut;
3055 
3056  const Particle* p = particles.appendNew(Particle({0.8 , 0.8 , 1.131370849898476039041351 , 1.6}, 22,
3057  Particle::c_Unflavored, Particle::c_Undefined, 1));
3058 
3059  StoreObjPtr<ParticleList> outputList("pList1");
3060  DataStore::Instance().setInitializeActive(true);
3061  outputList.registerInDataStore(flags);
3062  DataStore::Instance().setInitializeActive(false);
3063  outputList.create();
3064  outputList->initialize(22, "pList1");
3065 
3066  particles.appendNew(Particle({0.5 , 0.4953406774856531014212777 , 0.5609256753154148484773173 , 0.9}, 22,
3067  Particle::c_Unflavored, Particle::c_Undefined, 2)); //m=0.135
3068  particles.appendNew(Particle({0.5 , 0.2 , 0.72111 , 0.9}, 22, Particle::c_Unflavored, Particle::c_Undefined, 3)); //m=0.3582
3069  particles.appendNew(Particle({0.4 , 0.2 , 0.78102 , 0.9}, 22, Particle::c_Unflavored, Particle::c_Undefined, 4)); //m=0.3908
3070  particles.appendNew(Particle({0.5 , 0.4 , 0.89443 , 1.1}, 22, Particle::c_Unflavored, Particle::c_Undefined, 5)); //m=0.2369
3071  particles.appendNew(Particle({0.3 , 0.3 , 0.42426 , 0.6}, 22, Particle::c_Unflavored, Particle::c_Undefined, 6)); //m=0.0036
3072 
3073  outputList->addParticle(1, 22, Particle::c_Unflavored);
3074  outputList->addParticle(2, 22, Particle::c_Unflavored);
3075  outputList->addParticle(3, 22, Particle::c_Unflavored);
3076  outputList->addParticle(4, 22, Particle::c_Unflavored);
3077  outputList->addParticle(5, 22, Particle::c_Unflavored);
3078 
3079  StoreObjPtr<ParticleList> outputList2("pList2");
3080  DataStore::Instance().setInitializeActive(true);
3081  outputList2.registerInDataStore(flags);
3082  DataStore::Instance().setInitializeActive(false);
3083  outputList2.create();
3084  outputList2->initialize(22, "pList2");
3085 
3086  particles.appendNew(Particle({0.5 , -0.4 , 0.63246 , 0.9}, 22, Particle::c_Unflavored, Particle::c_Undefined, 7)); //m=1.1353
3087  particles.appendNew(Particle({0.5 , 0.2 , 0.72111 , 0.9}, 22, Particle::c_Unflavored, Particle::c_Undefined, 8)); //m=0.3582
3088  particles.appendNew(Particle({0.4 , 0.2 , 0.78102 , 0.9}, 22, Particle::c_Unflavored, Particle::c_Undefined, 9)); //m=0.3908
3089  particles.appendNew(Particle({0.5 , 0.4 , 0.89443 , 1.1}, 22, Particle::c_Unflavored, Particle::c_Undefined, 10)); //m=0.2369
3090  particles.appendNew(Particle({0.3 , 0.3 , 0.42426 , 0.6}, 22, Particle::c_Unflavored, Particle::c_Undefined, 11)); //m=0.0036
3091 
3092  outputList2->addParticle(6, 22, Particle::c_Unflavored);
3093  outputList2->addParticle(7, 22, Particle::c_Unflavored);
3094  outputList2->addParticle(8, 22, Particle::c_Unflavored);
3095  outputList2->addParticle(9, 22, Particle::c_Unflavored);
3096  outputList2->addParticle(10, 22, Particle::c_Unflavored);
3097 
3098  const Manager::Var* var = Manager::Instance().getVariable("veto(pList1, 0.130 < M < 0.140)");
3099  ASSERT_NE(var, nullptr);
3100  EXPECT_DOUBLE_EQ(var->function(p), 1);
3101 
3102  var = Manager::Instance().getVariable("veto(pList2, 0.130 < M < 0.140)");
3103  ASSERT_NE(var, nullptr);
3104  EXPECT_DOUBLE_EQ(var->function(p), 0);
3105 
3106  }
3107 
3108  TEST_F(MetaVariableTest, averageValueInList)
3109  {
3110  // we need the particles StoreArray
3111  StoreArray<Particle> particles;
3112  DataStore::EStoreFlags flags = DataStore::c_DontWriteOut;
3113 
3114  // create a photon list for testing
3115  StoreObjPtr<ParticleList> gammalist("testGammaList");
3116  DataStore::Instance().setInitializeActive(true);
3117  gammalist.registerInDataStore(flags);
3118  DataStore::Instance().setInitializeActive(false);
3119  gammalist.create();
3120  gammalist->initialize(22, "testGammaList");
3121 
3122  // create some photons in an stdvector
3123  std::vector<Particle> gammavector = {
3124  Particle({0.5 , 0.4 , 0.5 , 0.8}, 22, Particle::c_Unflavored, Particle::c_Undefined, 0),
3125  Particle({0.5 , 0.2 , 0.7 , 0.9}, 22, Particle::c_Unflavored, Particle::c_Undefined, 1),
3126  Particle({0.4 , 0.2 , 0.7 , 0.9}, 22, Particle::c_Unflavored, Particle::c_Undefined, 2),
3127  Particle({0.5 , 0.4 , 0.8 , 1.1}, 22, Particle::c_Unflavored, Particle::c_Undefined, 3),
3128  Particle({0.3 , 0.3 , 0.4 , 0.6}, 22, Particle::c_Unflavored, Particle::c_Undefined, 4)
3129  };
3130 
3131  // put the photons in the StoreArray
3132  for (const auto& g : gammavector)
3133  particles.appendNew(g);
3134 
3135  // put the photons in the test list
3136  for (size_t i = 0; i < gammavector.size(); i++)
3137  gammalist->addParticle(i, 22, Particle::c_Unflavored);
3138 
3139  // get the average px, py, pz, E of the gammas in the list
3140  const Manager::Var* vmeanpx = Manager::Instance().getVariable(
3141  "averageValueInList(testGammaList, px)");
3142  const Manager::Var* vmeanpy = Manager::Instance().getVariable(
3143  "averageValueInList(testGammaList, py)");
3144  const Manager::Var* vmeanpz = Manager::Instance().getVariable(
3145  "averageValueInList(testGammaList, pz)");
3146  const Manager::Var* vmeanE = Manager::Instance().getVariable(
3147  "averageValueInList(testGammaList, E)");
3148 
3149  EXPECT_FLOAT_EQ(vmeanpx->function(nullptr), 0.44);
3150  EXPECT_FLOAT_EQ(vmeanpy->function(nullptr), 0.3);
3151  EXPECT_FLOAT_EQ(vmeanpz->function(nullptr), 0.62);
3152  EXPECT_FLOAT_EQ(vmeanE->function(nullptr), 0.86);
3153 
3154  // wrong number of arguments (no variable provided)
3155  EXPECT_B2FATAL(Manager::Instance().getVariable("averageValueInList(testGammaList)"));
3156 
3157  // non-existing variable
3158  EXPECT_B2FATAL(Manager::Instance().getVariable("averageValueInList(testGammaList, NONEXISTANTVARIABLE)"));
3159 
3160  // non-existing list
3161  const Manager::Var* vnolist = Manager::Instance().getVariable(
3162  "averageValueInList(NONEXISTANTLIST, px)");
3163 
3164  EXPECT_B2FATAL(vnolist->function(nullptr));
3165  }
3166 
3167  TEST_F(MetaVariableTest, medianValueInList)
3168  {
3169  // we need the particles StoreArray
3170  StoreArray<Particle> particles;
3171  DataStore::EStoreFlags flags = DataStore::c_DontWriteOut;
3172 
3173  // create two photon lists for testing (one with odd and one with even number of particles)
3174  StoreObjPtr<ParticleList> oddgammalist("oddGammaList");
3175  DataStore::Instance().setInitializeActive(true);
3176  oddgammalist.registerInDataStore(flags);
3177  DataStore::Instance().setInitializeActive(false);
3178  oddgammalist.create();
3179  oddgammalist->initialize(22, "oddGammaList");
3180  StoreObjPtr<ParticleList> evengammalist("evenGammaList");
3181  DataStore::Instance().setInitializeActive(true);
3182  evengammalist.registerInDataStore(flags);
3183  DataStore::Instance().setInitializeActive(false);
3184  evengammalist.create();
3185  evengammalist->initialize(22, "evenGammaList");
3186 
3187  // create some photons in an stdvector
3188  std::vector<Particle> gammavector = {
3189  Particle({0.5 , 0.4 , 0.5 , 0.8}, 22, Particle::c_Unflavored, Particle::c_Undefined, 0),
3190  Particle({0.5 , 0.2 , 0.7 , 0.9}, 22, Particle::c_Unflavored, Particle::c_Undefined, 1),
3191  Particle({0.4 , 0.2 , 0.7 , 0.9}, 22, Particle::c_Unflavored, Particle::c_Undefined, 2),
3192  Particle({0.5 , 0.4 , 0.8 , 1.1}, 22, Particle::c_Unflavored, Particle::c_Undefined, 3),
3193  Particle({0.3 , 0.3 , 0.4 , 0.6}, 22, Particle::c_Unflavored, Particle::c_Undefined, 4)
3194  };
3195 
3196  // put the photons in the StoreArray
3197  for (const auto& g : gammavector)
3198  particles.appendNew(g);
3199 
3200  // put the photons in the test lists
3201  if (gammavector.size() % 2 == 0) {
3202  evengammalist->addParticle(0, 22, Particle::c_Unflavored);
3203  } else
3204  oddgammalist->addParticle(0, 22, Particle::c_Unflavored);
3205  for (size_t i = 1; i < gammavector.size(); i++) {
3206  oddgammalist->addParticle(i, 22, Particle::c_Unflavored);
3207  evengammalist->addParticle(i, 22, Particle::c_Unflavored);
3208  }
3209 
3210  // get the median px, py, pz, E of the gammas in the list with odd number of particles
3211  const Manager::Var* voddmedianpx = Manager::Instance().getVariable(
3212  "medianValueInList(oddGammaList, px)");
3213  const Manager::Var* voddmedianpy = Manager::Instance().getVariable(
3214  "medianValueInList(oddGammaList, py)");
3215  const Manager::Var* voddmedianpz = Manager::Instance().getVariable(
3216  "medianValueInList(oddGammaList, pz)");
3217  const Manager::Var* voddmedianE = Manager::Instance().getVariable(
3218  "medianValueInList(oddGammaList, E)");
3219 
3220  EXPECT_FLOAT_EQ(voddmedianpx->function(nullptr), 0.5);
3221  EXPECT_FLOAT_EQ(voddmedianpy->function(nullptr), 0.3);
3222  EXPECT_FLOAT_EQ(voddmedianpz->function(nullptr), 0.7);
3223  EXPECT_FLOAT_EQ(voddmedianE->function(nullptr), 0.9);
3224 
3225  // get the median px, py, pz, E of the gammas in the list with odd number of particles
3226  const Manager::Var* vevenmedianpx = Manager::Instance().getVariable(
3227  "medianValueInList(evenGammaList, px)");
3228  const Manager::Var* vevenmedianpy = Manager::Instance().getVariable(
3229  "medianValueInList(evenGammaList, py)");
3230  const Manager::Var* vevenmedianpz = Manager::Instance().getVariable(
3231  "medianValueInList(evenGammaList, pz)");
3232  const Manager::Var* vevenmedianE = Manager::Instance().getVariable(
3233  "medianValueInList(evenGammaList, E)");
3234 
3235  EXPECT_FLOAT_EQ(vevenmedianpx->function(nullptr), 0.45);
3236  EXPECT_FLOAT_EQ(vevenmedianpy->function(nullptr), 0.25);
3237  EXPECT_FLOAT_EQ(vevenmedianpz->function(nullptr), 0.7);
3238  EXPECT_FLOAT_EQ(vevenmedianE->function(nullptr), 0.9);
3239 
3240  // wrong number of arguments (no variable provided)
3241  EXPECT_B2FATAL(Manager::Instance().getVariable("medianValueInList(oddGammaList)"));
3242 
3243  // non-existing variable
3244  EXPECT_B2FATAL(Manager::Instance().getVariable("medianValueInList(oddGammaList, NONEXISTANTVARIABLE)"));
3245 
3246  // non-existing list
3247  const Manager::Var* vnolist = Manager::Instance().getVariable(
3248  "medianValueInList(NONEXISTANTLIST, px)");
3249 
3250  EXPECT_B2FATAL(vnolist->function(nullptr));
3251  }
3252 
3253  TEST_F(MetaVariableTest, pValueCombination)
3254  {
3255  TLorentzVector momentum;
3256  StoreArray<Particle> particles;
3257  std::vector<int> daughterIndices;
3258  Particle KS(TLorentzVector(1.164, 1.55200, 0, 2), 310, Particle::c_Unflavored, Particle::c_Composite, 0);
3259  KS.setPValue(0.1);
3260  momentum += KS.get4Vector();
3261  Particle* newDaughters = particles.appendNew(KS);
3262  daughterIndices.push_back(newDaughters->getArrayIndex());
3263  Particle Jpsi(TLorentzVector(-1, 1, 1, 3.548), 443, Particle::c_Unflavored, Particle::c_Composite, 1);
3264  Jpsi.setPValue(0.9);
3265  momentum += Jpsi.get4Vector();
3266  newDaughters = particles.appendNew(Jpsi);
3267  daughterIndices.push_back(newDaughters->getArrayIndex());
3268  Particle* B = particles.appendNew(momentum, 521, Particle::c_Flavored, daughterIndices);
3269  B->setPValue(0.5);
3270 
3271  const Manager::Var* singlePvalue = Manager::Instance().getVariable("pValueCombination(chiProb)");
3272  ASSERT_NE(singlePvalue, nullptr);
3273  EXPECT_FLOAT_EQ(singlePvalue->function(B), 0.5);
3274 
3275  const Manager::Var* twoPvalues = Manager::Instance().getVariable("pValueCombination(chiProb, daughter(0, chiProb))");
3276  ASSERT_NE(twoPvalues, nullptr);
3277  EXPECT_FLOAT_EQ(twoPvalues->function(B), 0.05 * (1 - log(0.05)));
3278 
3279  const Manager::Var* threePvalues =
3280  Manager::Instance().getVariable("pValueCombination(chiProb, daughter(0, chiProb), daughter(1, chiProb))");
3281  ASSERT_NE(threePvalues, nullptr);
3282  EXPECT_FLOAT_EQ(threePvalues->function(B), 0.045 * (1 - log(0.045) + 0.5 * log(0.045) * log(0.045)));
3283 
3284  // wrong number of arguments
3285  EXPECT_B2FATAL(Manager::Instance().getVariable("pValueCombination()"));
3286 
3287  // non-existing variable
3288  EXPECT_B2FATAL(Manager::Instance().getVariable("pValueCombination(chiProb, NONEXISTANTVARIABLE)"));
3289  }
3290 
3291 
3292  TEST_F(MetaVariableTest, daughterCombinationOneGeneration)
3293  {
3294  const int nDaughters = 5;
3295  TLorentzVector momentum(0, 0, 0, 0);
3296  StoreArray<Particle> particles;
3297  std::vector<int> daughterIndices;
3298  std::vector<TLorentzVector> daughterMomenta;
3299 
3300  for (int i = 0; i < nDaughters; i++) {
3301  TLorentzVector mom(1, i * 0.5, 1, i * 1.0 + 2.0);
3302  Particle d(mom, (i % 2) ? -11 : 211);
3303  Particle* newDaughters = particles.appendNew(d);
3304  daughterIndices.push_back(newDaughters->getArrayIndex());
3305  daughterMomenta.push_back(mom);
3306  momentum = momentum + mom;
3307  }
3308  const Particle* p = particles.appendNew(momentum, 411, Particle::c_Flavored, daughterIndices);
3309 
3310  // Test the invariant mass of several combinations
3311  const Manager::Var* var = Manager::Instance().getVariable("daughterCombination(M, 0,1,2)");
3312  double M_test = (daughterMomenta[0] + daughterMomenta[1] + daughterMomenta[2]).Mag();
3313  EXPECT_FLOAT_EQ(var->function(p), M_test);
3314 
3315  var = Manager::Instance().getVariable("daughterCombination(M, 0,4)");
3316  M_test = (daughterMomenta[0] + daughterMomenta[4]).Mag();
3317  EXPECT_FLOAT_EQ(var->function(p), M_test);
3318 
3319 
3320  // Try with a non-lorentz invariant quantity
3321  var = Manager::Instance().getVariable("daughterCombination(p, 1, 0, 4)");
3322  double p_test = (daughterMomenta[0] + daughterMomenta[1] + daughterMomenta[4]).Vect().Mag();
3323  EXPECT_FLOAT_EQ(var->function(p), p_test);
3324 
3325 
3326  // errors and bad stuff
3327  EXPECT_B2FATAL(Manager::Instance().getVariable("daughterCombination(aVeryNonExistingVariableSillyName, 1, 0, 4)"));
3328 
3329  var = Manager::Instance().getVariable("daughterCombination(M, 1, 0, 100)");
3330  EXPECT_B2WARNING(var->function(p));
3331  EXPECT_TRUE(std::isnan(var->function(p)));
3332 
3333 
3334  var = Manager::Instance().getVariable("daughterCombination(M, 1, -1)");
3335  EXPECT_B2WARNING(var->function(p));
3336  EXPECT_TRUE(std::isnan(var->function(p)));
3337 
3338 
3339  var = Manager::Instance().getVariable("daughterCombination(M, 1, 0:1:0:0:1)");
3340  EXPECT_B2WARNING(var->function(p));
3341  EXPECT_TRUE(std::isnan(var->function(p)));
3342 
3343  }
3344 
3345 
3346  TEST_F(MetaVariableTest, daughterCombinationTwoGenerations)
3347  {
3348  StoreArray<Particle> particles;
3349 
3350  // make a 1 -> 3 particle
3351 
3352  TLorentzVector momentum_1(0, 0, 0, 0);
3353  std::vector<TLorentzVector> daughterMomenta_1;
3354  std::vector<int> daughterIndices_1;
3355 
3356  for (int i = 0; i < 3; i++) {
3357  TLorentzVector mom(i * 0.2, 1, 1, i * 1.0 + 2.0);
3358  Particle d(mom, (i % 2) ? -11 : 211);
3359  Particle* newDaughters = particles.appendNew(d);
3360  daughterIndices_1.push_back(newDaughters->getArrayIndex());
3361  daughterMomenta_1.push_back(mom);
3362  momentum_1 = momentum_1 + mom;
3363  }
3364 
3365  const Particle* compositeDau_1 = particles.appendNew(momentum_1, 411, Particle::c_Flavored, daughterIndices_1);
3366 
3367 
3368  // make a 1 -> 2 particle
3369 
3370  TLorentzVector momentum_2(0, 0, 0, 0);
3371  std::vector<TLorentzVector> daughterMomenta_2;
3372  std::vector<int> daughterIndices_2;
3373 
3374  for (int i = 0; i < 2; i++) {
3375  TLorentzVector mom(1, 1, i * 0.3, i * 1.0 + 2.0);
3376  Particle d(mom, (i % 2) ? -11 : 211);
3377  Particle* newDaughters = particles.appendNew(d);
3378  daughterIndices_2.push_back(newDaughters->getArrayIndex());
3379  daughterMomenta_2.push_back(mom);
3380  momentum_2 = momentum_2 + mom;
3381  }
3382 
3383  const Particle* compositeDau_2 = particles.appendNew(momentum_2, 411, Particle::c_Flavored, daughterIndices_2);
3384 
3385 
3386  // make the composite particle
3387  std::vector<int> daughterIndices = {compositeDau_1->getArrayIndex(), compositeDau_2->getArrayIndex()};
3388  const Particle* p = particles.appendNew(momentum_2 + momentum_1, 111, Particle::c_Unflavored, daughterIndices);
3389 
3390 
3391  // Test the invariant mass of several combinations
3392  const Manager::Var* var = Manager::Instance().getVariable("daughterCombination(M, 0,1)");
3393  double M_test = (momentum_1 + momentum_2).Mag();
3394  EXPECT_FLOAT_EQ(var->function(p), M_test);
3395 
3396  // this should be the mass of the first daughter
3397  var = Manager::Instance().getVariable("daughterCombination(M, 0:0, 0:1, 0:2)");
3398  M_test = (momentum_1).Mag();
3399  EXPECT_FLOAT_EQ(var->function(p), M_test);
3400 
3401  // this should be a generic combinations
3402  var = Manager::Instance().getVariable("daughterCombination(M, 0:0, 0:1, 1:0)");
3403  M_test = (daughterMomenta_1[0] + daughterMomenta_1[1] + daughterMomenta_2[0]).Mag();
3404  EXPECT_FLOAT_EQ(var->function(p), M_test);
3405 
3406  }
3407 
3408 
3409  TEST_F(MetaVariableTest, useAlternativeDaughterHypothesis)
3410  {
3411  const int nDaughters = 5;
3412  StoreArray<Particle> particles;
3413 
3414  // Build a first Particle
3415  TLorentzVector momentum(0, 0, 0, 0);
3416  std::vector<int> daughterIndices;
3417  for (int i = 0; i < nDaughters; i++) {
3418  double px = i * 0.1;
3419  double py = i * 0.3;
3420  double pz = -i * 0.1 - 0.2;
3421 
3422  TLorentzVector mom(px, py, pz, 1);
3423  // all pions
3424  int pdgCode = Const::pion.getPDGCode();
3425  Particle d(mom, pdgCode);
3426  d.updateMass(pdgCode);
3427  mom.SetXYZM(px, py, pz, d.getMass());
3428 
3429  Particle* daughters = particles.appendNew(d);
3430  daughterIndices.push_back(daughters->getArrayIndex());
3431  momentum = momentum + mom;
3432  }
3433  const Particle* p = particles.appendNew(momentum, 411, Particle::c_Flavored, daughterIndices);
3434 
3435 
3436  // Build a second Particle with same momenta, but different mass hyp.
3437  TLorentzVector momentumAlt(0, 0, 0, 0);
3438  std::vector<int> daughterIndicesAlt;
3439  for (int i = 0; i < nDaughters; i++) {
3440  double px = i * 0.1;
3441  double py = i * 0.3;
3442  double pz = -i * 0.1 - 0.2;
3443 
3444  TLorentzVector mom(px, py, pz, 1);
3445  // all pions but the first two
3446  int pdgCode = Const::pion.getPDGCode();
3447  if (i == 0)
3448  pdgCode = Const::proton.getPDGCode(); // a proton
3449  if (i == 1)
3450  pdgCode = Const::kaon.getPDGCode(); // a K
3451  Particle d(mom, pdgCode);
3452  d.updateMass(pdgCode);
3453  mom.SetXYZM(px, py, pz, d.getMass());
3454 
3455  Particle* daughters = particles.appendNew(d);
3456  daughterIndicesAlt.push_back(daughters->getArrayIndex());
3457  momentumAlt = momentumAlt + mom;
3458  }
3459  const Particle* pAlt = particles.appendNew(momentumAlt, 411, Particle::c_Flavored, daughterIndicesAlt);
3460 
3461 
3462  // Test the invariant mass under the alternative hypothesis
3463  std::cout << "mass test" << std::endl;
3464  const Manager::Var* var = Manager::Instance().getVariable("useAlternativeDaughterHypothesis(M, 0:p+,1:K+)");
3465  const Manager::Var* varAlt = Manager::Instance().getVariable("M");
3466  EXPECT_FLOAT_EQ(var->function(p), varAlt->function(pAlt));
3467 
3468  // check it's really charge-insensitive...
3469  std::cout << "charge test" << std::endl;
3470  var = Manager::Instance().getVariable("useAlternativeDaughterHypothesis(M, 0:p+,1:K-)");
3471  EXPECT_FLOAT_EQ(var->function(p), varAlt->function(pAlt));
3472 
3473  // check the variable is not changing the 3-momentum
3474  std::cout << "momentum test" << std::endl;
3475  var = Manager::Instance().getVariable("useAlternativeDaughterHypothesis(p, 0:p+,1:K-)");
3476  varAlt = Manager::Instance().getVariable("p");
3477  EXPECT_FLOAT_EQ(var->function(p), varAlt->function(pAlt));
3478  EXPECT_FLOAT_EQ(var->function(p), varAlt->function(p));
3479  EXPECT_FLOAT_EQ(var->function(pAlt), varAlt->function(pAlt));
3480  }
3481 
3482 
3483 
3484 
3485  TEST_F(MetaVariableTest, daughterAngle)
3486  {
3487  StoreArray<Particle> particles;
3488 
3489  // make a 1 -> 3 particle
3490 
3491  TLorentzVector momentum_1(0, 0, 0, 0);
3492  std::vector<TLorentzVector> daughterMomenta_1;
3493  std::vector<int> daughterIndices_1;
3494 
3495  for (int i = 0; i < 3; i++) {
3496  TLorentzVector mom(i * 0.2, 1, 1, i * 1.0 + 2.0);
3497  Particle d(mom, (i % 2) ? -11 : 211);
3498  Particle* newDaughters = particles.appendNew(d);
3499  daughterIndices_1.push_back(newDaughters->getArrayIndex());
3500  daughterMomenta_1.push_back(mom);
3501  momentum_1 = momentum_1 + mom;
3502  }
3503 
3504  const Particle* compositeDau_1 = particles.appendNew(momentum_1, 411, Particle::c_Flavored, daughterIndices_1);
3505 
3506 
3507  // make a 1 -> 2 particle
3508 
3509  TLorentzVector momentum_2(0, 0, 0, 0);
3510  std::vector<TLorentzVector> daughterMomenta_2;
3511  std::vector<int> daughterIndices_2;
3512 
3513  for (int i = 0; i < 2; i++) {
3514  TLorentzVector mom(1, 1, i * 0.3, i * 1.0 + 2.0);
3515  Particle d(mom, (i % 2) ? -11 : 211);
3516  Particle* newDaughters = particles.appendNew(d);
3517  daughterIndices_2.push_back(newDaughters->getArrayIndex());
3518  daughterMomenta_2.push_back(mom);
3519  momentum_2 = momentum_2 + mom;
3520  }
3521 
3522  const Particle* compositeDau_2 = particles.appendNew(momentum_2, 411, Particle::c_Flavored, daughterIndices_2);
3523 
3524 
3525  // make the composite particle
3526  std::vector<int> daughterIndices = {compositeDau_1->getArrayIndex(), compositeDau_2->getArrayIndex()};
3527  const Particle* p = particles.appendNew(momentum_2 + momentum_1, 111, Particle::c_Unflavored, daughterIndices);
3528 
3529 
3530  // Test the invariant mass of several combinations
3531  const Manager::Var* var = Manager::Instance().getVariable("daughterAngle(0, 1)");
3532  double v_test = momentum_1.Vect().Angle(momentum_2.Vect());
3533  EXPECT_FLOAT_EQ(var->function(p), v_test);
3534 
3535  // this should be a generic combinations
3536  var = Manager::Instance().getVariable("daughterAngle(0:0, 1:0)");
3537  v_test = daughterMomenta_1[0].Vect().Angle(daughterMomenta_2[0].Vect());
3538  EXPECT_FLOAT_EQ(var->function(p), v_test);
3539 
3540  var = Manager::Instance().getVariable("daughterAngle( 1, -1)");
3541  EXPECT_B2WARNING(var->function(p));
3542  EXPECT_TRUE(std::isnan(var->function(p)));
3543 
3544  var = Manager::Instance().getVariable("daughterAngle(1, 0:1:0:0:1)");
3545  EXPECT_B2WARNING(var->function(p));
3546  EXPECT_TRUE(std::isnan(var->function(p)));
3547 
3548  }
3549 
3550 
3551  TEST_F(MetaVariableTest, varForFirstMCAncestorOfType)
3552  {
3553  DataStore::Instance().setInitializeActive(true);
3554  StoreArray<MCParticle> mcParticles;
3555  StoreArray<Particle> particles;
3556  particles.registerInDataStore();
3557  mcParticles.registerInDataStore();
3558  particles.registerRelationTo(mcParticles);
3559  StoreObjPtr<ParticleList> DList("D0:vartest");
3560  DList.registerInDataStore();
3561  DList.create();
3562  DList->initialize(421, "D0:vartest");
3563  DataStore::Instance().setInitializeActive(false);
3564  TLorentzVector momentum;
3565  TLorentzVector momentum_0;
3566  TLorentzVector momentum_1;
3567  std::vector<int> D_daughterIndices;
3568  std::vector<int> D_grandDaughterIndices_0;
3569  std::vector<int> D_grandDaughterIndices_1;
3570 
3571 
3572  // Create MC graph for D -> (K0s -> pi+ + pi-) (K0s -> pi+ + pi-)
3573  MCParticleGraph mcGraph;
3574 
3575  MCParticleGraph::GraphParticle& mcg_m = mcGraph.addParticle();
3576  MCParticleGraph::GraphParticle& mcg_d_0 = mcGraph.addParticle();
3577  MCParticleGraph::GraphParticle& mcg_d_1 = mcGraph.addParticle();
3578  MCParticleGraph::GraphParticle& mcg_gd_0_0 = mcGraph.addParticle();
3579  MCParticleGraph::GraphParticle& mcg_gd_0_1 = mcGraph.addParticle();
3580  MCParticleGraph::GraphParticle& mcg_gd_1_0 = mcGraph.addParticle();
3581  MCParticleGraph::GraphParticle& mcg_gd_1_1 = mcGraph.addParticle();
3582  MCParticleGraph::GraphParticle& mcg_not_child = mcGraph.addParticle();
3583 
3584  mcg_m.setPDG(421);
3585  mcg_m.set4Vector(TLorentzVector(7, 7, 7, 7));
3586  mcg_d_0.setPDG(-310);
3587  mcg_d_0.set4Vector(TLorentzVector(6, 6, 6, 6));
3588  mcg_d_1.setPDG(310);
3589  mcg_d_1.set4Vector(TLorentzVector(5, 5, 5, 5));
3590  mcg_gd_0_0.setPDG(211);
3591  mcg_gd_0_0.set4Vector(TLorentzVector(4, 4, 4, 4));
3592  mcg_gd_0_1.setPDG(-211);
3593  mcg_gd_0_1.set4Vector(TLorentzVector(3, 3, 3, 3));
3594  mcg_gd_1_0.setPDG(211);
3595  mcg_gd_1_0.set4Vector(TLorentzVector(2, 1, 2, 2));
3596  mcg_gd_1_1.setPDG(-211);
3597  mcg_gd_1_1.set4Vector(TLorentzVector(1, 1, 1, 1));
3598  mcg_not_child.setPDG(211);
3599  mcg_not_child.set4Vector(TLorentzVector(10, 10, 10, 10));
3600 
3601  mcg_d_0.comesFrom(mcg_m);
3602  mcg_d_1.comesFrom(mcg_m);
3603  mcg_gd_0_0.comesFrom(mcg_d_0);
3604  mcg_gd_0_1.comesFrom(mcg_d_0);
3605  mcg_gd_1_0.comesFrom(mcg_d_1);
3606  mcg_gd_1_1.comesFrom(mcg_d_1);
3607 
3608  mcGraph.generateList();
3609 
3610  // Get MC Particles from StoreArray
3611  auto* mc_not_child = mcParticles[0];
3612  auto* mc_m = mcParticles[1];
3613  auto* mc_d_0 = mcParticles[2];
3614  auto* mc_d_1 = mcParticles[3];
3615  auto* mc_gd_0_0 = mcParticles[4];
3616  auto* mc_gd_0_1 = mcParticles[5];
3617  auto* mc_gd_1_0 = mcParticles[6];
3618  auto* mc_gd_1_1 = mcParticles[7];
3619 
3620 
3621  mc_m->setStatus(MCParticle::c_PrimaryParticle);
3622  mc_d_0->setStatus(MCParticle::c_PrimaryParticle);
3623  mc_d_1->setStatus(MCParticle::c_PrimaryParticle);
3624  mc_gd_0_0->setStatus(MCParticle::c_PrimaryParticle);
3625  mc_gd_0_1->setStatus(MCParticle::c_PrimaryParticle);
3626  mc_gd_1_0->setStatus(MCParticle::c_PrimaryParticle);
3627  mc_gd_1_1->setStatus(MCParticle::c_PrimaryParticle);
3628  mc_not_child->setStatus(MCParticle::c_PrimaryParticle);
3629 
3630  // Creation of D decay: D->K0s(->pi pi) K0s(->pi pi)
3631 
3632  const Particle* D_gd_0_0 = particles.appendNew(TLorentzVector(0.0, 1, 1, 1), 211);
3633  const Particle* D_gd_0_1 = particles.appendNew(TLorentzVector(1.0, 1, 1, 1), -211);
3634  const Particle* D_gd_1_0 = particles.appendNew(TLorentzVector(2.0, 1, 1, 1), 211);
3635  const Particle* D_gd_1_1 = particles.appendNew(TLorentzVector(3.0, 1, 1, 1), -211);
3636 
3637  D_grandDaughterIndices_0.push_back(D_gd_0_0->getArrayIndex());
3638  D_grandDaughterIndices_0.push_back(D_gd_0_1->getArrayIndex());
3639  D_grandDaughterIndices_1.push_back(D_gd_1_0->getArrayIndex());
3640  D_grandDaughterIndices_1.push_back(D_gd_1_1->getArrayIndex());
3641  momentum_0 = D_gd_0_0->get4Vector() + D_gd_0_1->get4Vector();
3642  momentum_1 = D_gd_1_0->get4Vector() + D_gd_1_1->get4Vector();
3643 
3644 
3645  const Particle* D_d_0 = particles.appendNew(momentum_0, 310, Particle::c_Unflavored, D_grandDaughterIndices_0);
3646  const Particle* D_d_1 = particles.appendNew(momentum_1, 310, Particle::c_Unflavored, D_grandDaughterIndices_1);
3647 
3648 
3649  momentum = D_d_0->get4Vector() + D_d_1->get4Vector();
3650  D_daughterIndices.push_back(D_d_0->getArrayIndex());
3651  D_daughterIndices.push_back(D_d_1->getArrayIndex());
3652 
3653  const Particle* D_m = particles.appendNew(momentum, 421, Particle::c_Unflavored, D_daughterIndices);
3654  DList->addParticle(D_m);
3655 
3656  // Particle that is not an child
3657  const Particle* not_child = particles.appendNew(TLorentzVector(5.0, 1, 1, 1), 211);
3658 
3659  // Particle that is not an child and doesn't have MC particle
3660  const Particle* not_child_2 = particles.appendNew(TLorentzVector(6.0, 1, 1, 1), 211);
3661 
3662  // MC matching
3663  D_gd_0_0->addRelationTo(mc_gd_0_0);
3664  D_gd_0_1->addRelationTo(mc_gd_0_1);
3665  D_gd_1_0->addRelationTo(mc_gd_1_0);
3666  D_gd_1_1->addRelationTo(mc_gd_1_1);
3667  D_d_0->addRelationTo(mc_d_0);
3668  D_d_1->addRelationTo(mc_d_1);
3669  D_m->addRelationTo(mc_m);
3670  not_child->addRelationTo(mc_not_child);
3671 
3672  // All pions should have common D mother
3673  const Manager::Var* var_d = Manager::Instance().getVariable("varForFirstMCAncestorOfType(D0, mdstIndex)");
3674  ASSERT_NE(var_d, nullptr);
3675  EXPECT_TRUE(var_d->function(D_gd_0_0) >= 0);
3676  EXPECT_FLOAT_EQ(var_d->function(D_gd_0_0), var_d->function(D_gd_0_1));
3677  EXPECT_FLOAT_EQ(var_d->function(D_gd_1_0), var_d->function(D_gd_1_1));
3678  EXPECT_FLOAT_EQ(var_d->function(D_gd_0_0), var_d->function(D_gd_1_0));
3679  EXPECT_FLOAT_EQ(var_d->function(D_gd_0_1), var_d->function(D_gd_1_1));
3680  EXPECT_TRUE(std::isnan(var_d->function(not_child)));
3681  EXPECT_TRUE(std::isnan(var_d->function(not_child_2)));
3682 
3683 
3684  // // All but they have different K0s mothers
3685  const Manager::Var* var_310 = Manager::Instance().getVariable("varForFirstMCAncestorOfType(310, mdstIndex)");
3686  ASSERT_NE(var_310, nullptr);
3687  EXPECT_FLOAT_EQ(var_310->function(D_gd_0_0), var_310->function(D_gd_0_1));
3688  EXPECT_FLOAT_EQ(var_310->function(D_gd_1_0), var_310->function(D_gd_1_1));
3689  EXPECT_NE(var_310->function(D_gd_0_0), var_310->function(D_gd_1_0));
3690  EXPECT_NE(var_310->function(D_gd_0_1), var_310->function(D_gd_1_1));
3691  EXPECT_TRUE(std::isnan(var_310->function(not_child)));
3692  EXPECT_TRUE(std::isnan(var_310->function(not_child_2)));
3693  EXPECT_FLOAT_EQ(int(Manager::Instance().getVariable("varForFirstMCAncestorOfType(310, E)")->function(D_gd_0_0)), 10);
3694  }
3695 
3696  TEST_F(MetaVariableTest, isDescendantOfList)
3697  {
3698  DataStore::Instance().setInitializeActive(true);
3699  StoreObjPtr<ParticleList> DList("D0:vartest");
3700  DList.registerInDataStore();
3701  DList.create();
3702  DList->initialize(421, "D0:vartest");
3703  StoreObjPtr<ParticleList> BList("B:vartest");
3704  BList.registerInDataStore();
3705  BList.create();
3706  BList->initialize(521, "B:vartest");
3707  DataStore::Instance().setInitializeActive(false);
3708 
3709  TLorentzVector momentum;
3710  TLorentzVector momentum_0;
3711  TLorentzVector momentum_1;
3712  StoreArray<Particle> particles;
3713  std::vector<int> D_daughterIndices;
3714  std::vector<int> D_grandDaughterIndices_0;
3715  std::vector<int> D_grandDaughterIndices_1;
3716  std::vector<int> B_daughterIndices;
3717  std::vector<int> B_grandDaughterIndices;
3718  std::vector<int> B_grandGrandDaughterIndices;
3719 
3720  // Creation of D decay: D->K0s(->pi pi) K0s(->pi pi)
3721 
3722  const Particle* D_gd_0_0 = particles.appendNew(TLorentzVector(0.0, 1, 1, 1), 211, Particle::c_Flavored, Particle::c_Track, 0);
3723  const Particle* D_gd_0_1 = particles.appendNew(TLorentzVector(1.0, 1, 1, 1), -211, Particle::c_Flavored, Particle::c_Track, 1);
3724  const Particle* D_gd_1_0 = particles.appendNew(TLorentzVector(2.0, 1, 1, 1), 211, Particle::c_Flavored, Particle::c_Track, 2);
3725  const Particle* D_gd_1_1 = particles.appendNew(TLorentzVector(3.0, 1, 1, 1), -211, Particle::c_Flavored, Particle::c_Track, 3);
3726 
3727  D_grandDaughterIndices_0.push_back(D_gd_0_0->getArrayIndex());
3728  D_grandDaughterIndices_0.push_back(D_gd_0_1->getArrayIndex());
3729  D_grandDaughterIndices_1.push_back(D_gd_1_0->getArrayIndex());
3730  D_grandDaughterIndices_1.push_back(D_gd_1_1->getArrayIndex());
3731  momentum_0 = D_gd_0_0->get4Vector() + D_gd_0_1->get4Vector();
3732  momentum_1 = D_gd_1_0->get4Vector() + D_gd_1_1->get4Vector();
3733 
3734 
3735  const Particle* D_d_0 = particles.appendNew(momentum_0, 310, Particle::c_Unflavored, D_grandDaughterIndices_0);
3736  const Particle* D_d_1 = particles.appendNew(momentum_1, 310, Particle::c_Unflavored, D_grandDaughterIndices_1);
3737 
3738 
3739  momentum = D_d_0->get4Vector() + D_d_1->get4Vector();
3740  D_daughterIndices.push_back(D_d_0->getArrayIndex());
3741  D_daughterIndices.push_back(D_d_1->getArrayIndex());
3742 
3743  const Particle* D_m = particles.appendNew(momentum, 421, Particle::c_Unflavored, D_daughterIndices);
3744  DList->addParticle(D_m);
3745 
3746  // Creation of B decay B -> D(->K0s(->pi pi) pi) pi
3747 
3748  const Particle* B_d_1 = particles.appendNew(TLorentzVector(0.0, 1, 1, 1), 211, Particle::c_Flavored, Particle::c_Track, 4);
3749  const Particle* B_gd_0_1 = particles.appendNew(TLorentzVector(1.0, 1, 1, 1), -211, Particle::c_Flavored, Particle::c_Track, 5);
3750  const Particle* B_ggd_0_0_0 = particles.appendNew(TLorentzVector(2.0, 1, 1, 1), 211, Particle::c_Flavored, Particle::c_Track, 6);
3751  const Particle* B_ggd_0_0_1 = particles.appendNew(TLorentzVector(3.0, 1, 1, 1), -211, Particle::c_Flavored, Particle::c_Track, 7);
3752 
3753  B_grandGrandDaughterIndices.push_back(B_ggd_0_0_0->getArrayIndex());
3754  B_grandGrandDaughterIndices.push_back(B_ggd_0_0_1->getArrayIndex());
3755  momentum_0 = B_ggd_0_0_0->get4Vector() + B_ggd_0_0_1->get4Vector();
3756  const Particle* B_gd_0_0 = particles.appendNew(momentum_0, 310, Particle::c_Unflavored, B_grandGrandDaughterIndices);
3757 
3758  B_grandDaughterIndices.push_back(B_gd_0_0->getArrayIndex());
3759  B_grandDaughterIndices.push_back(B_gd_0_1->getArrayIndex());
3760  momentum_1 = B_gd_0_0->get4Vector() + B_gd_0_1->get4Vector();
3761  const Particle* B_d_0 = particles.appendNew(momentum_1, -411, Particle::c_Unflavored, B_grandDaughterIndices);
3762 
3763  B_daughterIndices.push_back(B_d_0->getArrayIndex());
3764  B_daughterIndices.push_back(B_d_1->getArrayIndex());
3765  momentum = B_d_0->get4Vector() + B_d_1->get4Vector();
3766  const Particle* B_m = particles.appendNew(momentum, 521, Particle::c_Unflavored, B_daughterIndices);
3767  BList->addParticle(B_m);
3768 
3769  // Particle that is not an child
3770  const Particle* not_child = particles.appendNew(TLorentzVector(5.0, 1, 1, 1), 211, Particle::c_Flavored, Particle::c_Track, 8);
3771 
3772 
3773  const Manager::Var* var_0 = Manager::Instance().getVariable("isDescendantOfList(D0:vartest)");
3774  ASSERT_NE(var_0, nullptr);
3775  EXPECT_FLOAT_EQ(var_0->function(D_gd_0_0), 1.);
3776  EXPECT_FLOAT_EQ(var_0->function(D_gd_0_1), 1.);
3777  EXPECT_FLOAT_EQ(var_0->function(D_gd_1_0), 1.);
3778  EXPECT_FLOAT_EQ(var_0->function(D_gd_1_1), 1.);
3779  EXPECT_FLOAT_EQ(var_0->function(D_d_0), 1.);
3780  EXPECT_FLOAT_EQ(var_0->function(D_d_1), 1.);
3781  EXPECT_FLOAT_EQ(var_0->function(B_ggd_0_0_0), 0.);
3782  EXPECT_FLOAT_EQ(var_0->function(B_ggd_0_0_1), 0.);
3783  EXPECT_FLOAT_EQ(var_0->function(B_gd_0_0), 0.);
3784  EXPECT_FLOAT_EQ(var_0->function(B_gd_0_1), 0.);
3785  EXPECT_FLOAT_EQ(var_0->function(B_d_0), 0.);
3786  EXPECT_FLOAT_EQ(var_0->function(B_d_1), 0.);
3787  EXPECT_FLOAT_EQ(var_0->function(not_child), 0.);
3788 
3789  const Manager::Var* var_0a = Manager::Instance().getVariable("isDaughterOfList(D0:vartest)");
3790  ASSERT_NE(var_0a, nullptr);
3791  EXPECT_FLOAT_EQ(var_0a->function(D_gd_0_0), 0.);
3792  EXPECT_FLOAT_EQ(var_0a->function(D_gd_0_1), 0.);
3793  EXPECT_FLOAT_EQ(var_0a->function(D_gd_1_0), 0.);
3794  EXPECT_FLOAT_EQ(var_0a->function(D_gd_1_1), 0.);
3795  EXPECT_FLOAT_EQ(var_0a->function(D_d_0), 1.);
3796  EXPECT_FLOAT_EQ(var_0a->function(D_d_1), 1.);
3797  EXPECT_FLOAT_EQ(var_0a->function(B_ggd_0_0_0), 0.);
3798  EXPECT_FLOAT_EQ(var_0a->function(B_ggd_0_0_1), 0.);
3799  EXPECT_FLOAT_EQ(var_0a->function(B_gd_0_0), 0.);
3800  EXPECT_FLOAT_EQ(var_0a->function(B_gd_0_1), 0.);
3801  EXPECT_FLOAT_EQ(var_0a->function(B_d_0), 0.);
3802  EXPECT_FLOAT_EQ(var_0a->function(B_d_1), 0.);
3803  EXPECT_FLOAT_EQ(var_0a->function(not_child), 0.);
3804 
3805  const Manager::Var* var_0b = Manager::Instance().getVariable("isGrandDaughterOfList(D0:vartest)");
3806  ASSERT_NE(var_0b, nullptr);
3807  EXPECT_FLOAT_EQ(var_0b->function(D_gd_0_0), 1.);
3808  EXPECT_FLOAT_EQ(var_0b->function(D_gd_0_1), 1.);
3809  EXPECT_FLOAT_EQ(var_0b->function(D_gd_1_0), 1.);
3810  EXPECT_FLOAT_EQ(var_0b->function(D_gd_1_1), 1.);
3811  EXPECT_FLOAT_EQ(var_0b->function(D_d_0), 0.);
3812  EXPECT_FLOAT_EQ(var_0b->function(D_d_1), 0.);
3813  EXPECT_FLOAT_EQ(var_0b->function(B_ggd_0_0_0), 0.);
3814  EXPECT_FLOAT_EQ(var_0b->function(B_ggd_0_0_1), 0.);
3815  EXPECT_FLOAT_EQ(var_0b->function(B_gd_0_0), 0.);
3816  EXPECT_FLOAT_EQ(var_0b->function(B_gd_0_1), 0.);
3817  EXPECT_FLOAT_EQ(var_0b->function(B_d_0), 0.);
3818  EXPECT_FLOAT_EQ(var_0b->function(B_d_1), 0.);
3819  EXPECT_FLOAT_EQ(var_0b->function(not_child), 0.);
3820 
3821  const Manager::Var* var_1 = Manager::Instance().getVariable("isDescendantOfList(D0:vartest, 1)");
3822  ASSERT_NE(var_1, nullptr);
3823  EXPECT_FLOAT_EQ(var_1->function(D_gd_0_0), 0.);
3824  EXPECT_FLOAT_EQ(var_1->function(D_gd_0_1), 0.);
3825  EXPECT_FLOAT_EQ(var_1->function(D_gd_1_0), 0.);
3826  EXPECT_FLOAT_EQ(var_1->function(D_gd_1_1), 0.);
3827  EXPECT_FLOAT_EQ(var_1->function(D_d_0), 1.);
3828  EXPECT_FLOAT_EQ(var_1->function(D_d_1), 1.);
3829  EXPECT_FLOAT_EQ(var_1->function(B_ggd_0_0_0), 0.);
3830  EXPECT_FLOAT_EQ(var_1->function(B_ggd_0_0_1), 0.);
3831  EXPECT_FLOAT_EQ(var_1->function(B_gd_0_0), 0.);
3832  EXPECT_FLOAT_EQ(var_1->function(B_gd_0_1), 0.);
3833  EXPECT_FLOAT_EQ(var_1->function(B_d_0), 0.);
3834  EXPECT_FLOAT_EQ(var_1->function(B_d_1), 0.);
3835  EXPECT_FLOAT_EQ(var_1->function(not_child), 0.);
3836 
3837  const Manager::Var* var_2 = Manager::Instance().getVariable("isDescendantOfList(D0:vartest, 2)");
3838  ASSERT_NE(var_2, nullptr);
3839  EXPECT_FLOAT_EQ(var_2->function(D_gd_0_0), 1.);
3840  EXPECT_FLOAT_EQ(var_2->function(D_gd_0_1), 1.);
3841  EXPECT_FLOAT_EQ(var_2->function(D_gd_1_0), 1.);
3842  EXPECT_FLOAT_EQ(var_2->function(D_gd_1_1), 1.);
3843  EXPECT_FLOAT_EQ(var_2->function(D_d_0), 0.);
3844  EXPECT_FLOAT_EQ(var_2->function(D_d_1), 0.);
3845  EXPECT_FLOAT_EQ(var_2->function(B_ggd_0_0_0), 0.);
3846  EXPECT_FLOAT_EQ(var_2->function(B_ggd_0_0_1), 0.);
3847  EXPECT_FLOAT_EQ(var_2->function(B_gd_0_0), 0.);
3848  EXPECT_FLOAT_EQ(var_2->function(B_gd_0_1), 0.);
3849  EXPECT_FLOAT_EQ(var_2->function(B_d_0), 0.);
3850  EXPECT_FLOAT_EQ(var_2->function(B_d_1), 0.);
3851  EXPECT_FLOAT_EQ(var_2->function(not_child), 0.);
3852 
3853  const Manager::Var* var_3 = Manager::Instance().getVariable("isDescendantOfList(D0:vartest, B:vartest)");
3854  ASSERT_NE(var_3, nullptr);
3855  EXPECT_FLOAT_EQ(var_3->function(D_gd_0_0), 1.);
3856  EXPECT_FLOAT_EQ(var_3->function(D_gd_0_1), 1.);
3857  EXPECT_FLOAT_EQ(var_3->function(D_gd_1_0), 1.);
3858  EXPECT_FLOAT_EQ(var_3->function(D_gd_1_1), 1.);
3859  EXPECT_FLOAT_EQ(var_3->function(D_d_0), 1.);
3860  EXPECT_FLOAT_EQ(var_3->function(D_d_1), 1.);
3861  EXPECT_FLOAT_EQ(var_3->function(B_ggd_0_0_0), 1.);
3862  EXPECT_FLOAT_EQ(var_3->function(B_ggd_0_0_1), 1.);
3863  EXPECT_FLOAT_EQ(var_3->function(B_gd_0_0), 1.);
3864  EXPECT_FLOAT_EQ(var_3->function(B_gd_0_1), 1.);
3865  EXPECT_FLOAT_EQ(var_3->function(B_d_0), 1.);
3866  EXPECT_FLOAT_EQ(var_3->function(B_d_1), 1.);
3867  EXPECT_FLOAT_EQ(var_3->function(not_child), 0.);
3868 
3869  const Manager::Var* var_4 = Manager::Instance().getVariable("isDescendantOfList(D0:vartest, B:vartest, -1)");
3870  ASSERT_NE(var_4, nullptr);
3871  EXPECT_FLOAT_EQ(var_4->function(D_gd_0_0), 1.);
3872  EXPECT_FLOAT_EQ(var_4->function(D_gd_0_1), 1.);
3873  EXPECT_FLOAT_EQ(var_4->function(D_gd_1_0), 1.);
3874  EXPECT_FLOAT_EQ(var_4->function(D_gd_1_1), 1.);
3875  EXPECT_FLOAT_EQ(var_4->function(D_d_0), 1.);
3876  EXPECT_FLOAT_EQ(var_4->function(D_d_1), 1.);
3877  EXPECT_FLOAT_EQ(var_4->function(B_ggd_0_0_0), 1.);
3878  EXPECT_FLOAT_EQ(var_4->function(B_ggd_0_0_1), 1.);
3879  EXPECT_FLOAT_EQ(var_4->function(B_gd_0_0), 1.);
3880  EXPECT_FLOAT_EQ(var_4->function(B_gd_0_1), 1.);
3881  EXPECT_FLOAT_EQ(var_4->function(B_d_0), 1.);
3882  EXPECT_FLOAT_EQ(var_4->function(B_d_1), 1.);
3883  EXPECT_FLOAT_EQ(var_4->function(not_child), 0.);
3884 
3885 
3886  const Manager::Var* var_5 = Manager::Instance().getVariable("isDescendantOfList(D0:vartest, B:vartest, 1)");
3887  ASSERT_NE(var_5, nullptr);
3888  EXPECT_FLOAT_EQ(var_5->function(D_gd_0_0), 0.);
3889  EXPECT_FLOAT_EQ(var_5->function(D_gd_0_1), 0.);
3890  EXPECT_FLOAT_EQ(var_5->function(D_gd_1_0), 0.);
3891  EXPECT_FLOAT_EQ(var_5->function(D_gd_1_1), 0.);
3892  EXPECT_FLOAT_EQ(var_5->function(D_d_0), 1.);
3893  EXPECT_FLOAT_EQ(var_5->function(D_d_1), 1.);
3894  EXPECT_FLOAT_EQ(var_5->function(B_ggd_0_0_0), 0.);
3895  EXPECT_FLOAT_EQ(var_5->function(B_ggd_0_0_1), 0.);
3896  EXPECT_FLOAT_EQ(var_5->function(B_gd_0_0), 0.);
3897  EXPECT_FLOAT_EQ(var_5->function(B_gd_0_1), 0.);
3898  EXPECT_FLOAT_EQ(var_5->function(B_d_0), 1.);
3899  EXPECT_FLOAT_EQ(var_5->function(B_d_1), 1.);
3900  EXPECT_FLOAT_EQ(var_5->function(not_child), 0.);
3901 
3902 
3903  const Manager::Var* var_6 = Manager::Instance().getVariable("isDescendantOfList(D0:vartest, B:vartest, 2)");
3904  ASSERT_NE(var_6, nullptr);
3905  EXPECT_FLOAT_EQ(var_6->function(D_gd_0_0), 1.);
3906  EXPECT_FLOAT_EQ(var_6->function(D_gd_0_1), 1.);
3907  EXPECT_FLOAT_EQ(var_6->function(D_gd_1_0), 1.);
3908  EXPECT_FLOAT_EQ(var_6->function(D_gd_1_1), 1.);
3909  EXPECT_FLOAT_EQ(var_6->function(D_d_0), 0.);
3910  EXPECT_FLOAT_EQ(var_6->function(D_d_1), 0.);
3911  EXPECT_FLOAT_EQ(var_6->function(B_ggd_0_0_0), 0.);
3912  EXPECT_FLOAT_EQ(var_6->function(B_ggd_0_0_1), 0.);
3913  EXPECT_FLOAT_EQ(var_6->function(B_gd_0_0), 1.);
3914  EXPECT_FLOAT_EQ(var_6->function(B_gd_0_1), 1.);
3915  EXPECT_FLOAT_EQ(var_6->function(B_d_0), 0.);
3916  EXPECT_FLOAT_EQ(var_6->function(B_d_1), 0.);
3917  EXPECT_FLOAT_EQ(var_6->function(not_child), 0.);
3918 
3919  const Manager::Var* var_7 = Manager::Instance().getVariable("isDescendantOfList(D0:vartest, B:vartest, 3)");
3920  ASSERT_NE(var_7, nullptr);
3921  EXPECT_FLOAT_EQ(var_7->function(D_gd_0_0), 0.);
3922  EXPECT_FLOAT_EQ(var_7->function(D_gd_0_1), 0.);
3923  EXPECT_FLOAT_EQ(var_7->function(D_gd_1_0), 0.);
3924  EXPECT_FLOAT_EQ(var_7->function(D_gd_1_1), 0.);
3925  EXPECT_FLOAT_EQ(var_7->function(D_d_0), 0.);
3926  EXPECT_FLOAT_EQ(var_7->function(D_d_1), 0.);
3927  EXPECT_FLOAT_EQ(var_7->function(B_ggd_0_0_0), 1.);
3928  EXPECT_FLOAT_EQ(var_7->function(B_ggd_0_0_1), 1.);
3929  EXPECT_FLOAT_EQ(var_7->function(B_gd_0_0), 0.);
3930  EXPECT_FLOAT_EQ(var_7->function(B_gd_0_1), 0.);
3931  EXPECT_FLOAT_EQ(var_7->function(B_d_0), 0.);
3932  EXPECT_FLOAT_EQ(var_7->function(B_d_1), 0.);
3933  EXPECT_FLOAT_EQ(var_7->function(not_child), 0.);
3934  }
3935 
3936 
3937  TEST_F(MetaVariableTest, isMCDescendantOfList)
3938  {
3939  DataStore::Instance().setInitializeActive(true);
3940  StoreArray<MCParticle> mcParticles;
3941  StoreArray<Particle> particles;
3942  particles.registerInDataStore();
3943  mcParticles.registerInDataStore();
3944  particles.registerRelationTo(mcParticles);
3945  StoreObjPtr<ParticleList> BList("B:vartest");
3946  BList.registerInDataStore();
3947  BList.create();
3948  BList->initialize(521, "B:vartest");
3949  StoreObjPtr<ParticleList> DList("D0:vartest");
3950  DList.registerInDataStore();
3951  DList.create();
3952  DList->initialize(421, "D0:vartest");
3953  DataStore::Instance().setInitializeActive(false);
3954  TLorentzVector momentum;
3955  TLorentzVector momentum_0;
3956  TLorentzVector momentum_1;
3957  std::vector<int> daughterIndices;
3958  std::vector<int> grandDaughterIndices;
3959  std::vector<int> grandGrandDaughterIndices;
3960  std::vector<int> D_daughterIndices;
3961  std::vector<int> D_grandDaughterIndices_0;
3962  std::vector<int> D_grandDaughterIndices_1;
3963 
3964 
3965  // Create MC graph for B+ -> (D -> (K0s -> pi+ + pi-) pi-) + pi+
3966  MCParticleGraph mcGraph;
3967 
3968  MCParticleGraph::GraphParticle& mcg_m = mcGraph.addParticle();
3969  MCParticleGraph::GraphParticle& mcg_d_0 = mcGraph.addParticle();
3970  MCParticleGraph::GraphParticle& mcg_d_1 = mcGraph.addParticle();
3971  MCParticleGraph::GraphParticle& mcg_gd_0_0 = mcGraph.addParticle();
3972  MCParticleGraph::GraphParticle& mcg_gd_0_1 = mcGraph.addParticle();
3973  MCParticleGraph::GraphParticle& mcg_ggd_0_0_0 = mcGraph.addParticle();
3974  MCParticleGraph::GraphParticle& mcg_ggd_0_0_1 = mcGraph.addParticle();
3975  MCParticleGraph::GraphParticle& mcg_not_child = mcGraph.addParticle();
3976 
3977  mcg_m.setPDG(521);
3978  mcg_d_0.setPDG(-411);
3979  mcg_d_1.setPDG(211);
3980  mcg_gd_0_0.setPDG(310);
3981  mcg_gd_0_1.setPDG(-211);
3982  mcg_ggd_0_0_0.setPDG(211);
3983  mcg_ggd_0_0_1.setPDG(-211);
3984  mcg_not_child.setPDG(211);
3985 
3986  mcg_d_0.comesFrom(mcg_m);
3987  mcg_d_1.comesFrom(mcg_m);
3988  mcg_gd_0_0.comesFrom(mcg_d_0);
3989  mcg_gd_0_1.comesFrom(mcg_d_0);
3990  mcg_ggd_0_0_0.comesFrom(mcg_gd_0_1);
3991  mcg_ggd_0_0_1.comesFrom(mcg_gd_0_1);
3992 
3993  mcGraph.generateList();
3994 
3995  // Get MC Particles from StoreArray
3996  auto* mc_m = mcParticles[0];
3997  auto* mc_d_0 = mcParticles[1];
3998  auto* mc_d_1 = mcParticles[2];
3999  auto* mc_gd_0_0 = mcParticles[3];
4000  auto* mc_gd_0_1 = mcParticles[4];
4001  auto* mc_ggd_0_0_0 = mcParticles[5];
4002  auto* mc_ggd_0_0_1 = mcParticles[6];
4003  auto* mc_not_child = mcParticles[7];
4004 
4005  mc_m->setStatus(MCParticle::c_PrimaryParticle);
4006  mc_d_0->setStatus(MCParticle::c_PrimaryParticle);
4007  mc_d_1->setStatus(MCParticle::c_PrimaryParticle);
4008  mc_gd_0_0->setStatus(MCParticle::c_PrimaryParticle);
4009  mc_gd_0_1->setStatus(MCParticle::c_PrimaryParticle);
4010  mc_ggd_0_0_0->setStatus(MCParticle::c_PrimaryParticle);
4011  mc_ggd_0_0_1->setStatus(MCParticle::c_PrimaryParticle);
4012  mc_not_child->setStatus(MCParticle::c_PrimaryParticle);
4013 
4014  // Creation of D decay: D->K0s(->pi pi) K0s(->pi pi) (not matched)
4015 
4016  const Particle* D_gd_0_0 = particles.appendNew(TLorentzVector(0.0, 1, 1, 1), 211);
4017  const Particle* D_gd_0_1 = particles.appendNew(TLorentzVector(1.0, 1, 1, 1), -211);
4018  const Particle* D_gd_1_0 = particles.appendNew(TLorentzVector(2.0, 1, 1, 1), 211);
4019  const Particle* D_gd_1_1 = particles.appendNew(TLorentzVector(3.0, 1, 1, 1), -211);
4020 
4021  D_grandDaughterIndices_0.push_back(D_gd_0_0->getArrayIndex());
4022  D_grandDaughterIndices_0.push_back(D_gd_0_1->getArrayIndex());
4023  D_grandDaughterIndices_1.push_back(D_gd_1_0->getArrayIndex());
4024  D_grandDaughterIndices_1.push_back(D_gd_1_1->getArrayIndex());
4025  momentum_0 = D_gd_0_0->get4Vector() + D_gd_0_1->get4Vector();
4026  momentum_1 = D_gd_1_0->get4Vector() + D_gd_1_1->get4Vector();
4027 
4028 
4029  const Particle* D_d_0 = particles.appendNew(momentum_0, 310, Particle::c_Unflavored, D_grandDaughterIndices_0);
4030  const Particle* D_d_1 = particles.appendNew(momentum_1, 310, Particle::c_Unflavored, D_grandDaughterIndices_1);
4031 
4032 
4033  momentum = D_d_0->get4Vector() + D_d_1->get4Vector();
4034  D_daughterIndices.push_back(D_d_0->getArrayIndex());
4035  D_daughterIndices.push_back(D_d_1->getArrayIndex());
4036 
4037  const Particle* D_m = particles.appendNew(momentum, 421, Particle::c_Unflavored, D_daughterIndices);
4038  DList->addParticle(D_m);
4039 
4040  // Creating B decay
4041  const Particle* d_1 = particles.appendNew(TLorentzVector(0.0, 1, 1, 1), 211);
4042  const Particle* gd_0_1 = particles.appendNew(TLorentzVector(1.0, 1, 1, 1), -211);
4043  const Particle* ggd_0_0_0 = particles.appendNew(TLorentzVector(2.0, 1, 1, 1), 211);
4044  const Particle* ggd_0_0_1 = particles.appendNew(TLorentzVector(3.0, 1, 1, 1), -211);
4045 
4046  grandGrandDaughterIndices.push_back(ggd_0_0_0->getArrayIndex());
4047  grandGrandDaughterIndices.push_back(ggd_0_0_1->getArrayIndex());
4048  momentum_0 = ggd_0_0_0->get4Vector() + ggd_0_0_1->get4Vector();
4049  const Particle* gd_0_0 = particles.appendNew(momentum_0, 310, Particle::c_Unflavored, grandGrandDaughterIndices);
4050 
4051  grandDaughterIndices.push_back(gd_0_0->getArrayIndex());
4052  grandDaughterIndices.push_back(gd_0_1->getArrayIndex());
4053  momentum_1 = gd_0_0->get4Vector() + gd_0_1->get4Vector();
4054  const Particle* d_0 = particles.appendNew(momentum_1, -411, Particle::c_Unflavored, grandDaughterIndices);
4055 
4056  daughterIndices.push_back(d_0->getArrayIndex());
4057  daughterIndices.push_back(d_1->getArrayIndex());
4058  momentum = d_0->get4Vector() + d_1->get4Vector();
4059  const Particle* m = particles.appendNew(momentum, 521, Particle::c_Unflavored, daughterIndices);
4060  BList->addParticle(m);
4061 
4062  // Particle that is not an child
4063  const Particle* not_child = particles.appendNew(TLorentzVector(5.0, 1, 1, 1), 211);
4064 
4065  // Particle that is not an child and doesn't have MC particle
4066  const Particle* not_child_2 = particles.appendNew(TLorentzVector(6.0, 1, 1, 1), 211);
4067 
4068  gd_0_0->addRelationTo(mc_gd_0_0);
4069  gd_0_1->addRelationTo(mc_gd_0_1);
4070  ggd_0_0_0->addRelationTo(mc_ggd_0_0_0);
4071  ggd_0_0_1->addRelationTo(mc_ggd_0_0_1);
4072  d_0->addRelationTo(mc_d_0);
4073  d_1->addRelationTo(mc_d_1);
4074  m->addRelationTo(mc_m);
4075  not_child->addRelationTo(mc_not_child);
4076 
4077  const Manager::Var* var_0 = Manager::Instance().getVariable("isMCDescendantOfList(B:vartest)");
4078  ASSERT_NE(var_0, nullptr);
4079  EXPECT_FLOAT_EQ(var_0->function(D_gd_0_0), 0.);
4080  EXPECT_FLOAT_EQ(var_0->function(D_gd_0_1), 0.);
4081  EXPECT_FLOAT_EQ(var_0->function(D_gd_1_0), 0.);
4082  EXPECT_FLOAT_EQ(var_0->function(D_gd_1_1), 0.);
4083  EXPECT_FLOAT_EQ(var_0->function(D_d_0), 0.);
4084  EXPECT_FLOAT_EQ(var_0->function(D_d_1), 0.);
4085  EXPECT_FLOAT_EQ(var_0->function(ggd_0_0_0), 1.);
4086  EXPECT_FLOAT_EQ(var_0->function(ggd_0_0_1), 1.);
4087  EXPECT_FLOAT_EQ(var_0->function(gd_0_0), 1.);
4088  EXPECT_FLOAT_EQ(var_0->function(gd_0_1), 1.);
4089  EXPECT_FLOAT_EQ(var_0->function(d_0), 1.);
4090  EXPECT_FLOAT_EQ(var_0->function(d_1), 1.);
4091  EXPECT_FLOAT_EQ(var_0->function(not_child), 0.);
4092  EXPECT_FLOAT_EQ(var_0->function(not_child_2), 0.);
4093 
4094  const Manager::Var* var_1 = Manager::Instance().getVariable("isMCDescendantOfList(B:vartest, D0:vartest)");
4095  ASSERT_NE(var_1, nullptr);
4096  EXPECT_FLOAT_EQ(var_1->function(D_gd_0_0), 0.);
4097  EXPECT_FLOAT_EQ(var_1->function(D_gd_0_1), 0.);
4098  EXPECT_FLOAT_EQ(var_1->function(D_gd_1_0), 0.);
4099  EXPECT_FLOAT_EQ(var_1->function(D_gd_1_1), 0.);
4100  EXPECT_FLOAT_EQ(var_1->function(D_d_0), 0.);
4101  EXPECT_FLOAT_EQ(var_1->function(D_d_1), 0.);
4102  EXPECT_FLOAT_EQ(var_1->function(ggd_0_0_0), 1.);
4103  EXPECT_FLOAT_EQ(var_1->function(ggd_0_0_1), 1.);
4104  EXPECT_FLOAT_EQ(var_1->function(gd_0_0), 1.);
4105  EXPECT_FLOAT_EQ(var_1->function(gd_0_1), 1.);
4106  EXPECT_FLOAT_EQ(var_1->function(d_0), 1.);
4107  EXPECT_FLOAT_EQ(var_1->function(d_1), 1.);
4108  EXPECT_FLOAT_EQ(var_1->function(not_child), 0.);
4109  EXPECT_FLOAT_EQ(var_1->function(not_child_2), 0.);
4110 
4111  const Manager::Var* var_2 = Manager::Instance().getVariable("isMCDescendantOfList(B:vartest, -1)");
4112  ASSERT_NE(var_2, nullptr);
4113  EXPECT_FLOAT_EQ(var_2->function(D_gd_0_0), 0.);
4114  EXPECT_FLOAT_EQ(var_2->function(D_gd_0_1), 0.);
4115  EXPECT_FLOAT_EQ(var_2->function(D_gd_1_0), 0.);
4116  EXPECT_FLOAT_EQ(var_2->function(D_gd_1_1), 0.);
4117  EXPECT_FLOAT_EQ(var_2->function(D_d_0), 0.);
4118  EXPECT_FLOAT_EQ(var_2->function(D_d_1), 0.);
4119  EXPECT_FLOAT_EQ(var_2->function(ggd_0_0_0), 1.);
4120  EXPECT_FLOAT_EQ(var_2->function(ggd_0_0_1), 1.);
4121  EXPECT_FLOAT_EQ(var_2->function(gd_0_0), 1.);
4122  EXPECT_FLOAT_EQ(var_2->function(gd_0_1), 1.);
4123  EXPECT_FLOAT_EQ(var_2->function(d_0), 1.);
4124  EXPECT_FLOAT_EQ(var_2->function(d_1), 1.);
4125  EXPECT_FLOAT_EQ(var_2->function(not_child), 0.);
4126  EXPECT_FLOAT_EQ(var_2->function(not_child_2), 0.);
4127 
4128  const Manager::Var* var_3 = Manager::Instance().getVariable("isMCDescendantOfList(B:vartest, 1)");
4129  ASSERT_NE(var_3, nullptr);
4130  EXPECT_FLOAT_EQ(var_3->function(D_gd_0_0), 0.);
4131  EXPECT_FLOAT_EQ(var_3->function(D_gd_0_1), 0.);
4132  EXPECT_FLOAT_EQ(var_3->function(D_gd_1_0), 0.);
4133  EXPECT_FLOAT_EQ(var_3->function(D_gd_1_1), 0.);
4134  EXPECT_FLOAT_EQ(var_3->function(D_d_0), 0.);
4135  EXPECT_FLOAT_EQ(var_3->function(D_d_1), 0.);
4136  EXPECT_FLOAT_EQ(var_3->function(ggd_0_0_0), 0.);
4137  EXPECT_FLOAT_EQ(var_3->function(ggd_0_0_1), 0.);
4138  EXPECT_FLOAT_EQ(var_3->function(gd_0_0), 0.);
4139  EXPECT_FLOAT_EQ(var_3->function(gd_0_1), 0.);
4140  EXPECT_FLOAT_EQ(var_3->function(d_0), 1.);
4141  EXPECT_FLOAT_EQ(var_3->function(d_1), 1.);
4142  EXPECT_FLOAT_EQ(var_3->function(not_child), 0.);
4143  EXPECT_FLOAT_EQ(var_3->function(not_child_2), 0.);
4144 
4145  const Manager::Var* var_4 = Manager::Instance().getVariable("isMCDescendantOfList(B:vartest, 2)");
4146  ASSERT_NE(var_4, nullptr);
4147  EXPECT_FLOAT_EQ(var_4->function(D_gd_0_0), 0.);
4148  EXPECT_FLOAT_EQ(var_4->function(D_gd_0_1), 0.);
4149  EXPECT_FLOAT_EQ(var_4->function(D_gd_1_0), 0.);
4150  EXPECT_FLOAT_EQ(var_4->function(D_gd_1_1), 0.);
4151  EXPECT_FLOAT_EQ(var_4->function(D_d_0), 0.);
4152  EXPECT_FLOAT_EQ(var_4->function(D_d_1), 0.);
4153  EXPECT_FLOAT_EQ(var_4->function(ggd_0_0_0), 0.);
4154  EXPECT_FLOAT_EQ(var_4->function(ggd_0_0_1), 0.);
4155  EXPECT_FLOAT_EQ(var_4->function(gd_0_0), 1.);
4156  EXPECT_FLOAT_EQ(var_4->function(gd_0_1), 1.);
4157  EXPECT_FLOAT_EQ(var_4->function(d_0), 0.);
4158  EXPECT_FLOAT_EQ(var_4->function(d_1), 0.);
4159  EXPECT_FLOAT_EQ(var_4->function(not_child), 0.);
4160  EXPECT_FLOAT_EQ(var_4->function(not_child_2), 0.);
4161 
4162 
4163  const Manager::Var* var_5 = Manager::Instance().getVariable("isMCDescendantOfList(B:vartest, 3)");
4164  ASSERT_NE(var_5, nullptr);
4165  EXPECT_FLOAT_EQ(var_5->function(D_gd_0_0), 0.);
4166  EXPECT_FLOAT_EQ(var_5->function(D_gd_0_1), 0.);
4167  EXPECT_FLOAT_EQ(var_5->function(D_gd_1_0), 0.);
4168  EXPECT_FLOAT_EQ(var_5->function(D_gd_1_1), 0.);
4169  EXPECT_FLOAT_EQ(var_5->function(D_d_0), 0.);
4170  EXPECT_FLOAT_EQ(var_5->function(D_d_1), 0.);
4171  EXPECT_FLOAT_EQ(var_5->function(ggd_0_0_0), 1.);
4172  EXPECT_FLOAT_EQ(var_5->function(ggd_0_0_1), 1.);
4173  EXPECT_FLOAT_EQ(var_5->function(gd_0_0), 0.);
4174  EXPECT_FLOAT_EQ(var_5->function(gd_0_1), 0.);
4175  EXPECT_FLOAT_EQ(var_5->function(d_0), 0.);
4176  EXPECT_FLOAT_EQ(var_5->function(d_1), 0.);
4177  EXPECT_FLOAT_EQ(var_5->function(not_child), 0.);
4178  EXPECT_FLOAT_EQ(var_5->function(not_child_2), 0.);
4179  }
4180 
4181 
4182 
4183 
4184 
4185  class PIDVariableTest : public ::testing::Test {
4186  protected:
4188  void SetUp() override
4189  {
4190  DataStore::Instance().setInitializeActive(true);
4191  StoreObjPtr<ParticleExtraInfoMap> peim;
4192  StoreArray<TrackFitResult> tfrs;
4193  StoreArray<MCParticle> mcparticles;
4194  StoreArray<PIDLikelihood> likelihood;
4195  StoreArray<Particle> particles;
4196  StoreArray<Track> tracks;
4197  peim.registerInDataStore();
4198  tfrs.registerInDataStore();
4199  mcparticles.registerInDataStore();
4200  likelihood.registerInDataStore();
4201  particles.registerInDataStore();
4202  tracks.registerInDataStore();
4203  particles.registerRelationTo(likelihood);
4204  tracks.registerRelationTo(likelihood);
4205  DataStore::Instance().setInitializeActive(false);
4206  }
4207 
4209  void TearDown() override
4210  {
4211  DataStore::Instance().reset();
4212  }
4213  };
4214 
4215  TEST_F(PIDVariableTest, LogLikelihood)
4216  {
4217  StoreArray<PIDLikelihood> likelihood;
4218  StoreArray<Particle> particles;
4219  StoreArray<Track> tracks;
4220  StoreArray<TrackFitResult> tfrs;
4221 
4222  // create tracks and trackFitResutls
4223  TRandom3 generator;
4224  const float pValue = 0.5;
4225  const float bField = 1.5;
4226  const int charge = 1;
4227  TMatrixDSym cov6(6);
4228  // Generate a random put orthogonal pair of vectors in the r-phi plane
4229  TVector2 d(generator.Uniform(-1, 1), generator.Uniform(-1, 1));
4230  TVector2 pt(generator.Uniform(-1, 1), generator.Uniform(-1, 1));
4231  d.Set(d.X(), -(d.X()*pt.Px()) / pt.Py());
4232  // Add a random z component
4233  TVector3 position(d.X(), d.Y(), generator.Uniform(-1, 1));
4234  TVector3 momentum(pt.Px(), pt.Py(), generator.Uniform(-1, 1));
4235 
4236  auto CDCValue = static_cast<unsigned long long int>(0x300000000000000);
4237  tfrs.appendNew(position, momentum, cov6, charge, Const::electron, pValue, bField, CDCValue, 16777215, 0);
4238  Track mytrack;
4239  mytrack.setTrackFitResultIndex(Const::electron, 0);
4240  Track* allTrack = tracks.appendNew(mytrack);
4241  Track* noSVDTrack = tracks.appendNew(mytrack);
4242  Track* noSVDNoTOPTrack = tracks.appendNew(mytrack);
4243  Track* noPIDTrack = tracks.appendNew(mytrack);
4244  Track* dEdxTrack = tracks.appendNew(mytrack);
4245 
4246  // Fill by hand likelihood values for all the detectors and hypothesis
4247  // This is clearly not a physical case, since a particle cannot leave good
4248  // signals in both TOP and ARICH
4249  auto* lAll = likelihood.appendNew();
4250  lAll->setLogLikelihood(Const::TOP, Const::electron, 0.18);
4251  lAll->setLogLikelihood(Const::ARICH, Const::electron, 0.16);
4252  lAll->setLogLikelihood(Const::ECL, Const::electron, 0.14);
4253  lAll->setLogLikelihood(Const::CDC, Const::electron, 0.12);
4254  lAll->setLogLikelihood(Const::SVD, Const::electron, 0.1);
4255  lAll->setLogLikelihood(Const::KLM, Const::electron, 0.01);
4256 
4257  lAll->setLogLikelihood(Const::TOP, Const::muon, 0.5);
4258  lAll->setLogLikelihood(Const::ARICH, Const::muon, 0.52);
4259  lAll->setLogLikelihood(Const::ECL, Const::muon, 0.54);
4260  lAll->setLogLikelihood(Const::CDC, Const::muon, 0.56);
4261  lAll->setLogLikelihood(Const::SVD, Const::muon, 0.58);
4262  lAll->setLogLikelihood(Const::KLM, Const::muon, 0.8);
4263 
4264  lAll->setLogLikelihood(Const::TOP, Const::pion, 0.2);
4265  lAll->setLogLikelihood(Const::ARICH, Const::pion, 0.22);
4266  lAll->setLogLikelihood(Const::ECL, Const::pion, 0.24);
4267  lAll->setLogLikelihood(Const::CDC, Const::pion, 0.26);
4268  lAll->setLogLikelihood(Const::SVD, Const::pion, 0.28);
4269  lAll->setLogLikelihood(Const::KLM, Const::pion, 0.2);
4270 
4271  lAll->setLogLikelihood(Const::TOP, Const::kaon, 0.3);
4272  lAll->setLogLikelihood(Const::ARICH, Const::kaon, 0.32);
4273  lAll->setLogLikelihood(Const::ECL, Const::kaon, 0.34);
4274  lAll->setLogLikelihood(Const::CDC, Const::kaon, 0.36);
4275  lAll->setLogLikelihood(Const::SVD, Const::kaon, 0.38);
4276  lAll->setLogLikelihood(Const::KLM, Const::kaon, 0.2);
4277 
4278  lAll->setLogLikelihood(Const::TOP, Const::proton, 0.4);
4279  lAll->setLogLikelihood(Const::ARICH, Const::proton, 0.42);
4280  lAll->setLogLikelihood(Const::ECL, Const::proton, 0.44);
4281  lAll->setLogLikelihood(Const::CDC, Const::proton, 0.46);
4282  lAll->setLogLikelihood(Const::SVD, Const::proton, 0.48);
4283  lAll->setLogLikelihood(Const::KLM, Const::proton, 0.02);
4284 
4285  lAll->setLogLikelihood(Const::TOP, Const::deuteron, 0.6);
4286  lAll->setLogLikelihood(Const::ARICH, Const::deuteron, 0.62);
4287  lAll->setLogLikelihood(Const::ECL, Const::deuteron, 0.64);
4288  lAll->setLogLikelihood(Const::CDC, Const::deuteron, 0.66);
4289  lAll->setLogLikelihood(Const::SVD, Const::deuteron, 0.68);
4290  lAll->setLogLikelihood(Const::KLM, Const::deuteron, 0.02);
4291 
4292  // Likelihoods for all detectors but SVD
4293  auto* lAllNoSVD = likelihood.appendNew();
4294  // Likelihoods for all detectors but SVD and TOP
4295  auto* lAllNoSVDNoTOP = likelihood.appendNew();
4296 
4297  for (unsigned int iDet(0); iDet < Const::PIDDetectors::c_size; iDet++) {
4298  const auto det = Const::PIDDetectors::c_set[iDet];
4299  for (const auto& hypo : Const::chargedStableSet) {
4300  if (det != Const::SVD) {
4301  lAllNoSVD->setLogLikelihood(det, hypo, lAll->getLogL(hypo, det));
4302  if (det != Const::TOP) {
4303  lAllNoSVDNoTOP->setLogLikelihood(det, hypo, lAll->getLogL(hypo, det));
4304  }
4305  }
4306  }
4307  }
4308 
4309  // Likelihoods for a dEdx only case
4310  auto* ldEdx = likelihood.appendNew();
4311  ldEdx->setLogLikelihood(Const::CDC, Const::electron, 0.12);
4312  ldEdx->setLogLikelihood(Const::SVD, Const::electron, 0.1);
4313 
4314  ldEdx->setLogLikelihood(Const::CDC, Const::pion, 0.26);
4315  ldEdx->setLogLikelihood(Const::SVD, Const::pion, 0.28);
4316 
4317  ldEdx->setLogLikelihood(Const::CDC, Const::kaon, 0.36);
4318  ldEdx->setLogLikelihood(Const::SVD, Const::kaon, 0.38);
4319 
4320  ldEdx->setLogLikelihood(Const::CDC, Const::proton, 0.46);
4321  ldEdx->setLogLikelihood(Const::SVD, Const::proton, 0.48);
4322 
4323  ldEdx->setLogLikelihood(Const::CDC, Const::muon, 0.56);
4324  ldEdx->setLogLikelihood(Const::SVD, Const::muon, 0.58);
4325 
4326  ldEdx->setLogLikelihood(Const::CDC, Const::deuteron, 0.66);
4327  ldEdx->setLogLikelihood(Const::SVD, Const::deuteron, 0.68);
4328 
4329  allTrack->addRelationTo(lAll);
4330  noSVDTrack->addRelationTo(lAllNoSVD);
4331  noSVDNoTOPTrack->addRelationTo(lAllNoSVDNoTOP);
4332  dEdxTrack->addRelationTo(ldEdx);
4333 
4334  // Table with the sum(LogL) for several cases
4335  // All dEdx AllNoSVD AllNoSVDNoTOP
4336  // e 0.71 0.22 0.61 0.43
4337  // mu 3.5 1.14 2.92 2.42
4338  // pi 1.4 0.54 1.12 0.92
4339  // k 1.9 0.74 1.52 1.22
4340  // p 2.22 0.94 1.74 1.34
4341  // d 3.22 1.34 2.54 1.94
4342 
4343  auto* particleAll = particles.appendNew(allTrack, Const::pion);
4344  auto* particleNoSVD = particles.appendNew(noSVDTrack, Const::pion);
4345  auto* particleNoSVDNoTOP = particles.appendNew(noSVDNoTOPTrack, Const::pion);
4346  auto* particledEdx = particles.appendNew(dEdxTrack, Const::pion);
4347  auto* particleNoID = particles.appendNew(noPIDTrack, Const::pion);
4348 
4349  double numsumexp = std::exp(0.71) + std::exp(3.5) + std::exp(1.4) + std::exp(1.9) + std::exp(2.22) + std::exp(3.22);
4350  double numsumexp_noSVD = std::exp(0.61) + std::exp(2.92) + std::exp(1.12) + std::exp(1.52) + std::exp(1.74) + std::exp(2.54);
4351  double numsumexp_noSVDNoTOP = std::exp(0.43) + std::exp(2.42) + std::exp(0.92) + std::exp(1.22) + std::exp(1.34) + std::exp(1.94);
4352 
4353  // Basic PID quantities. Currently just wrappers for global probability.
4354  EXPECT_FLOAT_EQ(electronID(particleNoSVD), std::exp(0.61) / numsumexp_noSVD);
4355  EXPECT_FLOAT_EQ(muonID(particleNoSVD), std::exp(2.92) / numsumexp_noSVD);
4356  EXPECT_FLOAT_EQ(pionID(particleNoSVD), std::exp(1.12) / numsumexp_noSVD);
4357  EXPECT_FLOAT_EQ(kaonID(particleNoSVD), std::exp(1.52) / numsumexp_noSVD);
4358  EXPECT_FLOAT_EQ(protonID(particleNoSVD), std::exp(1.74) / numsumexp_noSVD);
4359  EXPECT_FLOAT_EQ(deuteronID(particleNoSVD), std::exp(2.54) / numsumexp_noSVD);
4360 
4361  // smart PID that takes the hypothesis into account
4362  auto* particleElectronNoSVD = particles.appendNew(noSVDTrack, Const::electron);
4363  auto* particleMuonNoSVD = particles.appendNew(noSVDTrack, Const::muon);
4364  auto* particleKaonNoSVD = particles.appendNew(noSVDTrack, Const::kaon);
4365  auto* particleProtonNoSVD = particles.appendNew(noSVDTrack, Const::proton);
4366  auto* particleDeuteronNoSVD = particles.appendNew(noSVDTrack, Const::deuteron);
4367 
4368  EXPECT_FLOAT_EQ(particleID(particleNoSVD), std::exp(1.12) / numsumexp_noSVD); // there's already a pion
4369  EXPECT_FLOAT_EQ(particleID(particleElectronNoSVD), std::exp(0.61) / numsumexp_noSVD);
4370  EXPECT_FLOAT_EQ(particleID(particleMuonNoSVD), std::exp(2.92) / numsumexp_noSVD);
4371  EXPECT_FLOAT_EQ(particleID(particleKaonNoSVD), std::exp(1.52) / numsumexp_noSVD);
4372  EXPECT_FLOAT_EQ(particleID(particleProtonNoSVD), std::exp(1.74) / numsumexp_noSVD);
4373  EXPECT_FLOAT_EQ(particleID(particleDeuteronNoSVD), std::exp(2.54) / numsumexp_noSVD);
4374 
4375  // TEMP: Electron ID w/o the TOP.
4376  EXPECT_FLOAT_EQ(electronID_noTOP(particleNoSVDNoTOP), std::exp(0.43) / numsumexp_noSVDNoTOP);
4377 
4378  // Hadron ID vars w/ SVD info included. Only pi, K, p.
4379  double numsumexp_had = std::exp(1.4) + std::exp(1.9) + std::exp(2.22);
4380  EXPECT_FLOAT_EQ(pionID_SVD(particleAll), std::exp(1.4) / numsumexp_had);
4381  EXPECT_FLOAT_EQ(kaonID_SVD(particleAll), std::exp(1.9) / numsumexp_had);
4382  EXPECT_FLOAT_EQ(protonID_SVD(particleAll), std::exp(2.22) / numsumexp_had);
4383  std::vector<double> v_pi_K {211., 321.};
4384  std::vector<double> v_pi_p {211., 2212.};
4385  std::vector<double> v_K_p {321., 2212.};
4386  EXPECT_FLOAT_EQ(binaryPID_SVD(particleAll, v_pi_K), std::exp(1.4) / (std::exp(1.4) + std::exp(1.9)));
4387  EXPECT_FLOAT_EQ(binaryPID_SVD(particleAll, v_pi_p), std::exp(1.4) / (std::exp(1.4) + std::exp(2.22)));
4388  EXPECT_FLOAT_EQ(binaryPID_SVD(particleAll, v_K_p), std::exp(1.9) / (std::exp(1.9) + std::exp(2.22)));
4389 
4390  // Check what happens if no Likelihood is available
4391  EXPECT_TRUE(std::isnan(electronID(particleNoID)));
4392  EXPECT_TRUE(std::isnan(muonID(particleNoID)));
4393  EXPECT_TRUE(std::isnan(pionID(particleNoID)));
4394  EXPECT_TRUE(std::isnan(kaonID(particleNoID)));
4395  EXPECT_TRUE(std::isnan(protonID(particleNoID)));
4396  EXPECT_TRUE(std::isnan(deuteronID(particleNoID)));
4397 
4398  //expert stuff: LogL values
4399  EXPECT_FLOAT_EQ(Manager::Instance().getVariable("pidLogLikelihoodValueExpert(11, TOP)")->function(particleAll), 0.18);
4400  EXPECT_FLOAT_EQ(Manager::Instance().getVariable("pidLogLikelihoodValueExpert(11, ALL)")->function(particleAll), 0.71);
4401  EXPECT_FLOAT_EQ(Manager::Instance().getVariable("pidLogLikelihoodValueExpert(2212, TOP, CDC)")->function(particleAll), 0.86);
4402 
4403  // global probability
4404  EXPECT_FLOAT_EQ(Manager::Instance().getVariable("pidProbabilityExpert(1000010020, ALL)")->function(particleAll),
4405  std::exp(3.22) / numsumexp);
4406  EXPECT_FLOAT_EQ(Manager::Instance().getVariable("pidProbabilityExpert(2212, ALL)")->function(particleAll),
4407  std::exp(2.22) / numsumexp);
4408  EXPECT_FLOAT_EQ(Manager::Instance().getVariable("pidProbabilityExpert(211, ALL)")->function(particleAll),
4409  std::exp(1.4) / numsumexp);
4410  EXPECT_FLOAT_EQ(Manager::Instance().getVariable("pidProbabilityExpert(321, ALL)")->function(particleAll),
4411  std::exp(1.9) / numsumexp);
4412  EXPECT_FLOAT_EQ(Manager::Instance().getVariable("pidProbabilityExpert(13, ALL)")->function(particleAll),
4413  std::exp(3.5) / numsumexp);
4414  EXPECT_FLOAT_EQ(Manager::Instance().getVariable("pidProbabilityExpert(11, ALL)")->function(particleAll),
4415  std::exp(0.71) / numsumexp);
4416  EXPECT_FLOAT_EQ(Manager::Instance().getVariable("pidProbabilityExpert(211, ALL)")->function(particledEdx),
4417  std::exp(0.54) / (std::exp(0.22) + std::exp(1.14) + std::exp(0.54) + std::exp(0.74) + std::exp(0.94) + std::exp(1.34)));
4418  EXPECT_FLOAT_EQ(Manager::Instance().getVariable("pidProbabilityExpert(211, ALL)")->function(particledEdx),
4419  Manager::Instance().getVariable("pidProbabilityExpert(211, CDC, SVD)")->function(particleAll));
4420  EXPECT_FLOAT_EQ(Manager::Instance().getVariable("pidProbabilityExpert(211, CDC)")->function(particledEdx),
4421  Manager::Instance().getVariable("pidProbabilityExpert(211, CDC)")->function(particleAll));
4422  EXPECT_FLOAT_EQ(Manager::Instance().getVariable("pidProbabilityExpert(321, CDC)")->function(particleAll),
4423  std::exp(0.36) / (std::exp(0.12) + std::exp(0.26) + std::exp(0.36) + std::exp(0.46) + std::exp(0.56) + std::exp(0.66)));
4424 
4425  // binary probability
4426  EXPECT_FLOAT_EQ(Manager::Instance().getVariable("pidPairProbabilityExpert(321, 2212, ALL)")->function(particleAll),
4427  1.0 / (1.0 + std::exp(2.22 - 1.9)));
4428  EXPECT_FLOAT_EQ(Manager::Instance().getVariable("pidPairProbabilityExpert(321, 2212, ALL)")->function(particledEdx),
4429  1.0 / (1.0 + std::exp(0.94 - 0.74)));
4430  EXPECT_FLOAT_EQ(Manager::Instance().getVariable("pidPairProbabilityExpert(321, 2212, CDC, SVD)")->function(particleAll),
4431  1.0 / (1.0 + std::exp(0.94 - 0.74)));
4432 
4433  // No likelihood available
4434  EXPECT_TRUE(std::isnan(Manager::Instance().getVariable("pidPairProbabilityExpert(321, 2212, KLM)")->function(particledEdx)));
4435  EXPECT_TRUE(std::isnan(Manager::Instance().getVariable("pidLogLikelihoodValueExpert(11, TOP, CDC, SVD)")->function(particleNoID)));
4436  EXPECT_TRUE(std::isnan(Manager::Instance().getVariable("pidLogLikelihoodValueExpert(11, TOP)")->function(particledEdx)));
4437  EXPECT_TRUE(std::isnan(Manager::Instance().getVariable("pidPairProbabilityExpert(321, 2212, KLM)")->function(particledEdx)));
4438  EXPECT_TRUE(std::isnan(Manager::Instance().getVariable("pidPairProbabilityExpert(321, 2212, ECL, TOP, ARICH)")->function(
4439  particledEdx)));
4440  EXPECT_FALSE(std::isnan(Manager::Instance().getVariable("pidPairProbabilityExpert(321, 2212, ECL, TOP, ARICH, SVD)")->function(
4441  particledEdx)));
4442  //Mostlikely PDG tests:
4443  EXPECT_FLOAT_EQ(Manager::Instance().getVariable("pidMostLikelyPDG()")->function(particledEdx), 1.00001e+09);
4444  EXPECT_FLOAT_EQ(Manager::Instance().getVariable("pidMostLikelyPDG(0.5, 0.1, 0.1, 0.1, 0.1, 0.1)")->function(particledEdx), 11);
4445  EXPECT_FLOAT_EQ(Manager::Instance().getVariable("pidMostLikelyPDG(0.1, 0.5, 0.1, 0.1, 0.1, 0.1)")->function(particledEdx), 13);
4446  EXPECT_FLOAT_EQ(Manager::Instance().getVariable("pidMostLikelyPDG(0.1, 0.1, 0.5, 0.1, 0.1, 0.1)")->function(particledEdx), 211);
4447  EXPECT_FLOAT_EQ(Manager::Instance().getVariable("pidMostLikelyPDG(0.1, 0.1, 0.1, 0.5, 0.1, 0.1)")->function(particledEdx), 321);
4448  EXPECT_FLOAT_EQ(Manager::Instance().getVariable("pidMostLikelyPDG(0.1, 0.1, 0.1, 0.1, 0.5, 0.1)")->function(particledEdx), 2212);
4449  EXPECT_FLOAT_EQ(Manager::Instance().getVariable("pidMostLikelyPDG(0, 1., 0, 0, 0, 0)")->function(particledEdx), 13);
4450  }
4451 
4452  TEST_F(PIDVariableTest, MissingLikelihood)
4453  {
4454  StoreArray<PIDLikelihood> likelihood;
4455  StoreArray<Particle> particles;
4456  StoreArray<Track> tracks;
4457  StoreArray<TrackFitResult> tfrs;
4458 
4459  // create tracks and trackFitResutls
4460  TRandom3 generator;
4461  const float pValue = 0.5;
4462  const float bField = 1.5;
4463  const int charge = 1;
4464  TMatrixDSym cov6(6);
4465  // Generate a random put orthogonal pair of vectors in the r-phi plane
4466  TVector2 d(generator.Uniform(-1, 1), generator.Uniform(-1, 1));
4467  TVector2 pt(generator.Uniform(-1, 1), generator.Uniform(-1, 1));
4468  d.Set(d.X(), -(d.X()*pt.Px()) / pt.Py());
4469  // Add a random z component
4470  TVector3 position(d.X(), d.Y(), generator.Uniform(-1, 1));
4471  TVector3 momentum(pt.Px(), pt.Py(), generator.Uniform(-1, 1));
4472 
4473  auto CDCValue = static_cast<unsigned long long int>(0x300000000000000);
4474  tfrs.appendNew(position, momentum, cov6, charge, Const::electron, pValue, bField, CDCValue, 16777215, 0);
4475  Track mytrack;
4476  mytrack.setTrackFitResultIndex(Const::electron, 0);
4477  Track* savedTrack1 = tracks.appendNew(mytrack);
4478  Track* savedTrack2 = tracks.appendNew(mytrack);
4479  Track* savedTrack3 = tracks.appendNew(mytrack);
4480  Track* savedTrack4 = tracks.appendNew(mytrack);
4481 
4482  auto* l1 = likelihood.appendNew();
4483  l1->setLogLikelihood(Const::TOP, Const::electron, 0.18);
4484  l1->setLogLikelihood(Const::ECL, Const::electron, 0.14);
4485  savedTrack1->addRelationTo(l1);
4486 
4487  auto* electron = particles.appendNew(savedTrack1, Const::electron);
4488 
4489  auto* l2 = likelihood.appendNew();
4490  l2->setLogLikelihood(Const::TOP, Const::pion, 0.2);
4491  l2->setLogLikelihood(Const::ARICH, Const::pion, 0.22);
4492  l2->setLogLikelihood(Const::ECL, Const::pion, 0.24);
4493  l2->setLogLikelihood(Const::CDC, Const::pion, 0.26);
4494  l2->setLogLikelihood(Const::SVD, Const::pion, 0.28);
4495  savedTrack2->addRelationTo(l2);
4496 
4497  auto* pion = particles.appendNew(savedTrack2, Const::pion);
4498 
4499  auto* l3 = likelihood.appendNew();
4500  l3->setLogLikelihood(Const::TOP, Const::kaon, 0.3);
4501  l3->setLogLikelihood(Const::ARICH, Const::kaon, 0.32);
4502  savedTrack3->addRelationTo(l3);
4503 
4504  auto* kaon = particles.appendNew(savedTrack3, Const::kaon);
4505 
4506  auto* l4 = likelihood.appendNew();
4507  l4->setLogLikelihood(Const::ARICH, Const::proton, 0.42);
4508  l4->setLogLikelihood(Const::ECL, Const::proton, 0.44);
4509  l4->setLogLikelihood(Const::CDC, Const::proton, 0.46);
4510  l4->setLogLikelihood(Const::SVD, Const::proton, 0.48);
4511  savedTrack4->addRelationTo(l4);
4512 
4513  auto* proton = particles.appendNew(savedTrack4, Const::proton);
4514 
4515  const Manager::Var* varMissECL = Manager::Instance().getVariable("pidMissingProbabilityExpert(ECL)");
4516  const Manager::Var* varMissTOP = Manager::Instance().getVariable("pidMissingProbabilityExpert(TOP)");
4517  const Manager::Var* varMissARICH = Manager::Instance().getVariable("pidMissingProbabilityExpert(ARICH)");
4518 
4519 
4520  EXPECT_FLOAT_EQ(varMissTOP->function(electron), 0.0);
4521  EXPECT_FLOAT_EQ(varMissTOP->function(pion), 0.0);
4522  EXPECT_FLOAT_EQ(varMissTOP->function(kaon), 0.0);
4523  EXPECT_FLOAT_EQ(varMissTOP->function(proton), 1.0);
4524 
4525  EXPECT_FLOAT_EQ(varMissARICH->function(electron), 1.0);
4526  EXPECT_FLOAT_EQ(varMissARICH->function(pion), 0.0);
4527  EXPECT_FLOAT_EQ(varMissARICH->function(kaon), 0.0);
4528  EXPECT_FLOAT_EQ(varMissARICH->function(proton), 0.0);
4529 
4530  EXPECT_FLOAT_EQ(varMissECL->function(electron), 0.0);
4531  EXPECT_FLOAT_EQ(varMissECL->function(pion), 0.0);
4532  EXPECT_FLOAT_EQ(varMissECL->function(kaon), 1.0);
4533  EXPECT_FLOAT_EQ(varMissECL->function(proton), 0.0);
4534  }
4535 
4536  class FlightInfoTest : public ::testing::Test {
4537  protected:
4539  void SetUp() override
4540  {
4541  DataStore::Instance().setInitializeActive(true);
4542  StoreArray<Particle>().registerInDataStore();
4543  StoreArray<MCParticle>().registerInDataStore();
4544  StoreArray<MCParticle> mcParticles;
4545  StoreArray<Particle> particles;
4546  particles.registerRelationTo(mcParticles);
4547  StoreObjPtr<ParticleExtraInfoMap>().registerInDataStore();
4548  DataStore::Instance().setInitializeActive(false);
4549 
4550 
4551  // Insert MC particle logic here
4552  MCParticle mcKs;
4553  mcKs.setPDG(310);
4554  mcKs.setProductionVertex(1.0, 1.0, 0.0);
4555  mcKs.setDecayVertex(4.0, 5.0, 0.0);
4556  mcKs.setProductionTime(0);
4557  mcKs.setMassFromPDG();
4558  mcKs.setMomentum(1.164, 1.55200, 0);
4559  float decayTime = 5 * mcKs.getMass() / mcKs.getEnergy();
4560  mcKs.setDecayTime(decayTime);
4561  mcKs.setStatus(MCParticle::c_PrimaryParticle);
4562  MCParticle* newMCKs = mcParticles.appendNew(mcKs);
4563 
4564 
4565 
4566  MCParticle mcDp;
4567  mcDp.setPDG(411);
4568  mcDp.setDecayVertex(1.0, 1.0, 0.0);
4569  mcDp.setMassFromPDG();
4570  mcDp.setStatus(MCParticle::c_PrimaryParticle);
4571  MCParticle* newMCDp = mcParticles.appendNew(mcDp);
4572 
4573  // Insert Reco particle logic here
4574  TLorentzVector momentum;
4575  TMatrixFSym error(7);
4576  error.Zero();
4577  error(0, 0) = 0.05;
4578  error(1, 1) = 0.2;
4579  error(2, 2) = 0.4;
4580  error(3, 3) = 0.01;
4581  error(4, 4) = 0.04;
4582  error(5, 5) = 0.00875;
4583  error(6, 6) = 0.01;
4584  Particle pi(TLorentzVector(1.59607, 1.19705, 0, 2), 211);
4585  momentum += pi.get4Vector();
4586  Particle* newpi = particles.appendNew(pi);
4587 
4588 
4589  Particle Ks(TLorentzVector(1.164, 1.55200, 0, 2), 310, Particle::c_Unflavored, Particle::c_Composite, 0);
4590  Ks.setVertex(TVector3(4.0, 5.0, 0.0));
4591  Ks.setMomentumVertexErrorMatrix(error); // (order: px,py,pz,E,x,y,z)
4592  momentum += Ks.get4Vector();
4593  Ks.addExtraInfo("prodVertX", 1.0);
4594  Ks.addExtraInfo("prodVertY", 1.0);
4595  Ks.addExtraInfo("prodVertZ", 0.0);
4596  Ks.addExtraInfo("prodVertSxx", 0.04);
4597  Ks.addExtraInfo("prodVertSxy", 0.0);
4598  Ks.addExtraInfo("prodVertSxz", 0.0);
4599  Ks.addExtraInfo("prodVertSyx", 0.0);
4600  Ks.addExtraInfo("prodVertSyy", 0.00875);
4601  Ks.addExtraInfo("prodVertSyz", 0.0);
4602  Ks.addExtraInfo("prodVertSzx", 0.0);
4603  Ks.addExtraInfo("prodVertSzy", 0.0);
4604  Ks.addExtraInfo("prodVertSzz", 0.01);
4605  Particle* newKs = particles.appendNew(Ks);
4606  newKs->addRelationTo(newMCKs);
4607 
4608 
4609  Particle Dp(momentum, 411, Particle::c_Flavored, Particle::c_Composite, 0);
4610  Dp.appendDaughter(newpi);
4611  Dp.appendDaughter(newKs);
4612  TVector3 motherVtx(1.0, 1.0, 0.0);
4613  Dp.setVertex(motherVtx);
4614  Dp.setMomentumVertexErrorMatrix(error); // (order: px,py,pz,E,x,y,z)
4615  Dp.addExtraInfo("prodVertX", 0.0);
4616  Dp.addExtraInfo("prodVertY", 1.0);
4617  Dp.addExtraInfo("prodVertZ", -2.0);
4618  Dp.addExtraInfo("prodVertSxx", 0.04);
4619  Dp.addExtraInfo("prodVertSxy", 0.0);
4620  Dp.addExtraInfo("prodVertSxz", 0.0);
4621  Dp.addExtraInfo("prodVertSyx", 0.0);
4622  Dp.addExtraInfo("prodVertSyy", 0.01);
4623  Dp.addExtraInfo("prodVertSyz", 0.0);
4624  Dp.addExtraInfo("prodVertSzx", 0.0);
4625  Dp.addExtraInfo("prodVertSzy", 0.0);
4626  Dp.addExtraInfo("prodVertSzz", 0.1575);
4627  Particle* newDp = particles.appendNew(Dp);
4628  newDp->addRelationTo(newMCDp);
4629 
4630  }
4631 
4633  void TearDown() override
4634  {
4635  DataStore::Instance().reset();
4636  }
4637  };
4638  TEST_F(FlightInfoTest, flightDistance)
4639  {
4640  StoreArray<Particle> particles;
4641  const Particle* newKs = particles[1]; // Ks had flight distance of 5 cm
4642 
4643  const Manager::Var* var = Manager::Instance().getVariable("flightDistance");
4644  ASSERT_NE(var, nullptr);
4645  EXPECT_FLOAT_EQ(var->function(newKs), 5.0);
4646  }
4647  TEST_F(FlightInfoTest, flightDistanceErr)
4648  {
4649  StoreArray<Particle> particles;
4650  const Particle* newKs = particles[1]; // Ks had flight distance of 5 cm
4651 
4652  const Manager::Var* var = Manager::Instance().getVariable("flightDistanceErr");
4653  ASSERT_NE(var, nullptr);
4654  EXPECT_GT(var->function(newKs), 0.0);
4655  }
4656  TEST_F(FlightInfoTest, flightTime)
4657  {
4658  StoreArray<Particle> particles;
4659  const Particle* newKs = particles[1]; // Ks had flight time of 0.0427 us (t = d/c * m/p)
4660 
4661  const Manager::Var* var = Manager::Instance().getVariable("flightTime");
4662  ASSERT_NE(var, nullptr);
4663  EXPECT_FLOAT_EQ(var->function(newKs), 5.0 / Const::speedOfLight * newKs->getPDGMass() / newKs->getP());
4664  }
4665 
4666  TEST_F(FlightInfoTest, flightTimeErr)
4667  {
4668  StoreArray<Particle> particles;
4669  const Particle* newKs = particles[1]; // Ks should have positive flight distance uncertainty
4670 
4671  const Manager::Var* var = Manager::Instance().getVariable("flightTimeErr");
4672  ASSERT_NE(var, nullptr);
4673  EXPECT_GT(var->function(newKs), 0.0);
4674  }
4675 
4676 
4677  TEST_F(FlightInfoTest, flightDistanceOfDaughter)
4678  {
4679  StoreArray<Particle> particles;
4680  const Particle* newDp = particles[2]; // Get D+, its daughter Ks had flight distance of 5 cm
4681 
4682  const Manager::Var* var = Manager::Instance().getVariable("flightDistanceOfDaughter(1)");
4683  ASSERT_NE(var, nullptr);
4684  EXPECT_FLOAT_EQ(var->function(newDp), 5.0);
4685 
4686  var = Manager::Instance().getVariable("flightDistanceOfDaughter(3)");
4687  ASSERT_NE(var, nullptr);
4688  EXPECT_TRUE(std::isnan(var->function(newDp)));
4689  }
4690  TEST_F(FlightInfoTest, flightDistanceOfDaughterErr)
4691  {
4692  StoreArray<Particle> particles;
4693  const Particle* newDp = particles[2]; // Get D+, its daughter Ks should have positive flight distance uncertainty
4694 
4695  const Manager::Var* var = Manager::Instance().getVariable("flightDistanceOfDaughterErr(1)");
4696  ASSERT_NE(var, nullptr);
4697  EXPECT_GT(var->function(newDp), 0.0);
4698 
4699  var = Manager::Instance().getVariable("flightDistanceOfDaughterErr(3)");
4700  ASSERT_NE(var, nullptr);
4701  EXPECT_TRUE(std::isnan(var->function(newDp)));
4702  }
4703  TEST_F(FlightInfoTest, flightTimeOfDaughter)
4704  {
4705  StoreArray<Particle> particles;
4706  const Particle* newDp = particles[2]; // Get D+, its daughter Ks had flight time of 0.0427 us (t = d/c * m/p)
4707 
4708  const Manager::Var* var = Manager::Instance().getVariable("flightTimeOfDaughter(1)");
4709  ASSERT_NE(var, nullptr);
4710  const Particle* Ks = newDp->getDaughter(1);
4711 
4712  EXPECT_FLOAT_EQ(var->function(newDp), 5.0 / Const::speedOfLight * Ks->getPDGMass() / Ks->getP());
4713 
4714  var = Manager::Instance().getVariable("flightTimeOfDaughter(3)");
4715  ASSERT_NE(var, nullptr);
4716  EXPECT_TRUE(std::isnan(var->function(newDp)));
4717  }
4718  TEST_F(FlightInfoTest, flightTimeOfDaughterErr)
4719  {
4720  StoreArray<Particle> particles;
4721  const Particle* newDp = particles[2]; // Get D+, its daughter Ks should have positive flight time uncertainty
4722 
4723  const Manager::Var* var = Manager::Instance().getVariable("flightTimeOfDaughterErr(1)");
4724  ASSERT_NE(var, nullptr);
4725  EXPECT_GT(var->function(newDp), 0.0);
4726 
4727  var = Manager::Instance().getVariable("flightTimeOfDaughterErr(3)");
4728  ASSERT_NE(var, nullptr);
4729  EXPECT_TRUE(std::isnan(var->function(newDp)));
4730  }
4731  TEST_F(FlightInfoTest, mcFlightDistanceOfDaughter)
4732  {
4733  StoreArray<Particle> particles;
4734  const Particle* newDp = particles[2]; // Get D+, its daughter Ks had flight distance of 5 cm
4735 
4736  const Manager::Var* var = Manager::Instance().getVariable("mcFlightDistanceOfDaughter(1)");
4737  ASSERT_NE(var, nullptr);
4738 
4739  EXPECT_FLOAT_EQ(var->function(newDp), 5.0);
4740 
4741  var = Manager::Instance().getVariable("mcFlightDistanceOfDaughter(3)");
4742  ASSERT_NE(var, nullptr);
4743  EXPECT_TRUE(std::isnan(var->function(newDp)));
4744  }
4745  TEST_F(FlightInfoTest, mcFlightTimeOfDaughter)
4746  {
4747  StoreArray<Particle> particles;
4748  const Particle* newDp = particles[2]; // Get D+, its daughter Ks had flight time of 0.0427 us (t = d/c * m/p)
4749 
4750  const Manager::Var* var = Manager::Instance().getVariable("mcFlightTimeOfDaughter(1)");
4751  ASSERT_NE(var, nullptr);
4752  auto* Ks = newDp->getDaughter(1)->getRelatedTo<MCParticle>();
4753  // double p = Ks->getMomentum().Mag();
4754  // EXPECT_FLOAT_EQ(var->function(newDp), 5.0 / Const::speedOfLight * Ks->getMass() / p);
4755 
4756  EXPECT_FLOAT_EQ(var->function(newDp), Ks->getLifetime() / Ks->getEnergy()*Ks->getMass());
4757 
4758  var = Manager::Instance().getVariable("mcFlightTimeOfDaughter(3)");
4759  ASSERT_NE(var, nullptr);
4760  EXPECT_TRUE(std::isnan(var->function(newDp)));
4761  }
4762 
4763  TEST_F(FlightInfoTest, vertexDistance)
4764  {
4765  StoreArray<Particle> particles;
4766  const Particle* newKS = particles[1]; // Get KS, as it has both a production and decay vertex
4767 
4768  const Manager::Var* var = Manager::Instance().getVariable("vertexDistance");
4769  ASSERT_NE(var, nullptr);
4770  EXPECT_FLOAT_EQ(var->function(newKS), 5.0);
4771  }
4772 
4773  TEST_F(FlightInfoTest, vertexDistanceError)
4774  {
4775  StoreArray<Particle> particles;
4776  const Particle* newKS = particles[1]; // Get KS, as it has both a production and decay vertex
4777 
4778  const Manager::Var* var = Manager::Instance().getVariable("vertexDistanceErr");
4779  ASSERT_NE(var, nullptr);
4780  EXPECT_FLOAT_EQ(var->function(newKS), 0.2);
4781  }
4782 
4783  TEST_F(FlightInfoTest, vertexDistanceSignificance)
4784  {
4785  StoreArray<Particle> particles;
4786  const Particle* newKS = particles[1]; // Get KS, as it has both a production and decay vertex
4787 
4788  const Manager::Var* var = Manager::Instance().getVariable("vertexDistanceSignificance");
4789  ASSERT_NE(var, nullptr);
4790  EXPECT_FLOAT_EQ(var->function(newKS), 25);
4791  }
4792 
4793  TEST_F(FlightInfoTest, vertexDistanceOfDaughter)
4794  {
4795  StoreArray<Particle> particles;
4796  const Particle* newDp = particles[2]; // Get D+, its daughter KS has both a production and decay vertex
4797 
4798  const Manager::Var* var = Manager::Instance().getVariable("vertexDistanceOfDaughter(1, noIP)");
4799  ASSERT_NE(var, nullptr);
4800  EXPECT_FLOAT_EQ(var->function(newDp), 5.0);
4801 
4802  var = Manager::Instance().getVariable("vertexDistanceOfDaughter(1)");
4803  ASSERT_NE(var, nullptr);
4804  EXPECT_FLOAT_EQ(var->function(newDp), 6.0);
4805 
4806  var = Manager::Instance().getVariable("vertexDistanceOfDaughter(2)");
4807  ASSERT_NE(var, nullptr);
4808  EXPECT_TRUE(std::isnan(var->function(newDp)));
4809  }
4810 
4811  TEST_F(FlightInfoTest, vertexDistanceOfDaughterError)
4812  {
4813  StoreArray<Particle> particles;
4814  const Particle* newDp = particles[2]; // Get D+, its daughter KS has both a production and decay vertex
4815 
4816  const Manager::Var* var = Manager::Instance().getVariable("vertexDistanceOfDaughterErr(1, noIP)");
4817  ASSERT_NE(var, nullptr);
4818  EXPECT_FLOAT_EQ(var->function(newDp), 0.2);
4819 
4820  var = Manager::Instance().getVariable("vertexDistanceOfDaughterErr(1)");
4821  ASSERT_NE(var, nullptr);
4822  EXPECT_FLOAT_EQ(var->function(newDp), 0.25);
4823  }
4824 
4825  TEST_F(FlightInfoTest, vertexDistanceOfDaughterSignificance)
4826  {
4827  StoreArray<Particle> particles;
4828  const Particle* newDp = particles[2]; // Get D+, its daughter KS has both a production and decay vertex
4829 
4830  const Manager::Var* var = Manager::Instance().getVariable("vertexDistanceOfDaughterSignificance(1, noIP)");
4831  ASSERT_NE(var, nullptr);
4832  EXPECT_FLOAT_EQ(var->function(newDp), 25);
4833 
4834  var = Manager::Instance().getVariable("vertexDistanceOfDaughterSignificance(1)");
4835  ASSERT_NE(var, nullptr);
4836  EXPECT_FLOAT_EQ(var->function(newDp), 24);
4837  }
4838 
4839  class VertexVariablesTest : public ::testing::Test {
4840  protected:
4842  void SetUp() override
4843  {
4844  DataStore::Instance().setInitializeActive(true);
4845  StoreArray<Particle>().registerInDataStore();
4846  StoreArray<MCParticle>().registerInDataStore();
4847  StoreArray<MCParticle> mcParticles;
4848  StoreArray<Particle> particles;
4849  particles.registerRelationTo(mcParticles);
4850  StoreObjPtr<ParticleExtraInfoMap>().registerInDataStore();
4851  DataStore::Instance().setInitializeActive(false);
4852 
4853 
4854  // Insert MC particle logic here
4855  MCParticle mcKs;
4856  mcKs.setPDG(310);
4857  mcKs.setDecayVertex(4.0, 5.0, 0.0);
4858  mcKs.setProductionVertex(TVector3(1.0, 2.0, 3.0));
4859  mcKs.setMassFromPDG();
4860  mcKs.setMomentum(1.164, 1.55200, 0);
4861  mcKs.setStatus(MCParticle::c_PrimaryParticle);
4862  MCParticle* newMCKs = mcParticles.appendNew(mcKs);
4863 
4864  Particle Ks(TLorentzVector(1.164, 1.55200, 0, 2), 310);
4865  Ks.setVertex(TVector3(4.0, 5.0, 0.0));
4866  Ks.addExtraInfo("prodVertX", 1.0);
4867  Ks.addExtraInfo("prodVertY", 2.0);
4868  Ks.addExtraInfo("prodVertZ", 3.0);
4869  Ks.addExtraInfo("prodVertSxx", 0.1);
4870  Ks.addExtraInfo("prodVertSxy", 0.2);
4871  Ks.addExtraInfo("prodVertSxz", 0.3);
4872  Ks.addExtraInfo("prodVertSyx", 0.4);
4873  Ks.addExtraInfo("prodVertSyy", 0.5);
4874  Ks.addExtraInfo("prodVertSyz", 0.6);
4875  Ks.addExtraInfo("prodVertSzx", 0.7);
4876  Ks.addExtraInfo("prodVertSzy", 0.8);
4877  Ks.addExtraInfo("prodVertSzz", 0.9);
4878  Particle* newKs = particles.appendNew(Ks);
4879  newKs->addRelationTo(newMCKs);
4880  }
4881 
4883  void TearDown() override
4884  {
4885  DataStore::Instance().reset();
4886  }
4887  };
4888 
4889  // MC vertex tests
4890  TEST_F(VertexVariablesTest, mcDecayVertexX)
4891  {
4892  StoreArray<Particle> particles;
4893  const Particle* newKs = particles[0]; // Ks had truth decay x is 4.0
4894 
4895  const Manager::Var* var = Manager::Instance().getVariable("mcDecayVertexX");
4896  ASSERT_NE(var, nullptr);
4897  EXPECT_FLOAT_EQ(var->function(newKs), 4.0);
4898  }
4899 
4900  TEST_F(VertexVariablesTest, mcDecayVertexY)
4901  {
4902  StoreArray<Particle> particles;
4903  const Particle* newKs = particles[0]; // Ks had truth decay y is 5.0
4904 
4905  const Manager::Var* var = Manager::Instance().getVariable("mcDecayVertexY");
4906  ASSERT_NE(var, nullptr);
4907  EXPECT_FLOAT_EQ(var->function(newKs), 5.0);
4908  }
4909 
4910  TEST_F(VertexVariablesTest, mcDecayVertexZ)
4911  {
4912  StoreArray<Particle> particles;
4913  const Particle* newKs = particles[0]; // Ks had truth decay z is 0.0
4914 
4915  const Manager::Var* var = Manager::Instance().getVariable("mcDecayVertexZ");
4916  ASSERT_NE(var, nullptr);
4917  EXPECT_FLOAT_EQ(var->function(newKs), 0.0);
4918  }
4919 
4920 
4921  TEST_F(VertexVariablesTest, mcDecayVertexFromIPDistance)
4922  {
4923  StoreArray<Particle> particles;
4924  const Particle* newKs = particles[0]; // Ks had truth distance of sqrt(41)
4925 
4926  const Manager::Var* var = Manager::Instance().getVariable("mcDecayVertexFromIPDistance");
4927  ASSERT_NE(var, nullptr);
4928  EXPECT_FLOAT_EQ(var->function(newKs), sqrt(4.0 * 4.0 + 5.0 * 5.0));
4929  }
4930 
4931  TEST_F(VertexVariablesTest, mcDecayVertexRho)
4932  {
4933  StoreArray<Particle> particles;
4934  const Particle* newKs = particles[0]; // Ks had truth rho of sqrt(41)
4935 
4936  const Manager::Var* var = Manager::Instance().getVariable("mcDecayVertexRho");
4937  ASSERT_NE(var, nullptr);
4938  EXPECT_FLOAT_EQ(var->function(newKs), sqrt(4.0 * 4.0 + 5.0 * 5.0));
4939  }
4940 
4941  TEST_F(VertexVariablesTest, mcProductionVertexX)
4942  {
4943  StoreArray<Particle> particles;
4944  const Particle* newKs = particles[0]; // Ks had production vertex x of 1.0 cm
4945 
4946  const Manager::Var* var = Manager::Instance().getVariable("mcProductionVertexX");
4947  ASSERT_NE(var, nullptr);
4948  EXPECT_FLOAT_EQ(var->function(newKs), 1.0);
4949  }
4950 
4951  TEST_F(VertexVariablesTest, mcProductionVertexY)
4952  {
4953  StoreArray<Particle> particles;
4954  const Particle* newKs = particles[0]; // Ks had production vertex y of 2.0 cm
4955 
4956  const Manager::Var* var = Manager::Instance().getVariable("mcProductionVertexY");
4957  ASSERT_NE(var, nullptr);
4958  EXPECT_FLOAT_EQ(var->function(newKs), 2.0);
4959  }
4960 
4961  TEST_F(VertexVariablesTest, mcProductionVertexZ)
4962  {
4963  StoreArray<Particle> particles;
4964  const Particle* newKs = particles[0]; // Ks had production vertex z of 3.0 cm
4965 
4966  const Manager::Var* var = Manager::Instance().getVariable("mcProductionVertexZ");
4967  ASSERT_NE(var, nullptr);
4968  EXPECT_FLOAT_EQ(var->function(newKs), 3.0);
4969  }
4970 
4971  // Production position tests
4972 
4973  TEST_F(VertexVariablesTest, prodVertexX)
4974  {
4975  StoreArray<Particle> particles;
4976  const Particle* newKs = particles[0]; // Ks had production vertex x of 1.0 cm
4977 
4978  const Manager::Var* var = Manager::Instance().getVariable("prodVertexX");
4979  ASSERT_NE(var, nullptr);
4980  EXPECT_FLOAT_EQ(var->function(newKs), 1.0);
4981  }
4982  TEST_F(VertexVariablesTest, prodVertexY)
4983  {
4984  StoreArray<Particle> particles;
4985  const Particle* newKs = particles[0]; // Ks had production vertex y of 2.0 cm
4986 
4987  const Manager::Var* var = Manager::Instance().getVariable("prodVertexY");
4988  ASSERT_NE(var, nullptr);
4989  EXPECT_FLOAT_EQ(var->function(newKs), 2.0);
4990  }
4991  TEST_F(VertexVariablesTest, prodVertexZ)
4992  {
4993  StoreArray<Particle> particles;
4994  const Particle* newKs = particles[0]; // Ks had production vertex z of 3.0 cm
4995 
4996  const Manager::Var* var = Manager::Instance().getVariable("prodVertexZ");
4997  ASSERT_NE(var, nullptr);
4998  EXPECT_FLOAT_EQ(var->function(newKs), 3.0);
4999  }
5000 
5001  // Production Covariance tests
5002 
5003  TEST_F(VertexVariablesTest, prodVertexCov)
5004  {
5005  StoreArray<Particle> particles;
5006  const Particle* newKs = particles[0]; // Ks had production vertex covariance xx of .1 cm
5007 
5008  //const Manager::Var* var = Manager::Instance().getVariable("prodVertexCovXX");
5009  const Manager::Var* var = Manager::Instance().getVariable("prodVertexCov(0,0)");
5010  ASSERT_NE(var, nullptr);
5011  EXPECT_FLOAT_EQ(var->function(newKs), 0.1);
5012  var = Manager::Instance().getVariable("prodVertexCov(0,1)");
5013  EXPECT_FLOAT_EQ(var->function(newKs), 0.2);
5014  var = Manager::Instance().getVariable("prodVertexCov(0,2)");
5015  EXPECT_FLOAT_EQ(var->function(newKs), 0.3);
5016  var = Manager::Instance().getVariable("prodVertexCov(1,0)");
5017  EXPECT_FLOAT_EQ(var->function(newKs), 0.4);
5018  var = Manager::Instance().getVariable("prodVertexCov(1,1)");
5019  EXPECT_FLOAT_EQ(var->function(newKs), 0.5);
5020  var = Manager::Instance().getVariable("prodVertexCov(1,2)");
5021  EXPECT_FLOAT_EQ(var->function(newKs), 0.6);
5022  var = Manager::Instance().getVariable("prodVertexCov(2,0)");
5023  EXPECT_FLOAT_EQ(var->function(newKs), 0.7);
5024  var = Manager::Instance().getVariable("prodVertexCov(2,1)");
5025  EXPECT_FLOAT_EQ(var->function(newKs), 0.8);
5026  var = Manager::Instance().getVariable("prodVertexCov(2,2)");
5027  EXPECT_FLOAT_EQ(var->function(newKs), 0.9);
5028  var = Manager::Instance().getVariable("prodVertexXErr");
5029  ASSERT_NE(var, nullptr);
5030  EXPECT_FLOAT_EQ(var->function(newKs), sqrt(0.1));
5031  var = Manager::Instance().getVariable("prodVertexYErr");
5032  ASSERT_NE(var, nullptr);
5033  EXPECT_FLOAT_EQ(var->function(newKs), sqrt(0.5));
5034  var = Manager::Instance().getVariable("prodVertexZErr");
5035  ASSERT_NE(var, nullptr);
5036  EXPECT_FLOAT_EQ(var->function(newKs), sqrt(0.9));
5037  }
5038 
5039  // Tests of ContinuumSuppressionVariables
5040 
5041  TEST_F(MetaVariableTest, KSFWVariables)
5042  {
5043  // simple tests that do not require the ROE builder nor the CS builder
5044 
5045  // check that garbage input throws helpful B2FATAL
5046  EXPECT_B2FATAL(Manager::Instance().getVariable("KSFWVariables(NONSENSE)"));
5047 
5048  // check for NaN if we don't have a CS object for this particle
5049  StoreArray<Particle> myParticles;
5050  const Particle* particle_with_no_cs = myParticles.appendNew();
5051  const Manager::Var* var = Manager::Instance().getVariable("KSFWVariables(mm2)");
5052  EXPECT_TRUE(std::isnan(var->function(particle_with_no_cs)));
5053  }
5054 
5055  TEST_F(MetaVariableTest, CleoConeCS)
5056  {
5057  // simple tests that do not require the ROE builder nor the CS builder
5058 
5059  // check that garbage input throws helpful B2FATAL
5060  EXPECT_B2FATAL(Manager::Instance().getVariable("CleoConeCS(NONSENSE)"));
5061 
5062  // check that string other than ROE for second argument throws B2FATAL
5063  EXPECT_B2FATAL(Manager::Instance().getVariable("CleoConeCS(0, NOTROE)"));
5064 
5065  // check for NaN if we don't have a CS object for this particle
5066  StoreArray<Particle> myParticles;
5067  const Particle* particle_with_no_cs = myParticles.appendNew();
5068  const Manager::Var* var = Manager::Instance().getVariable("CleoConeCS(0)");
5069  EXPECT_TRUE(std::isnan(var->function(particle_with_no_cs)));
5070  }
5071 
5072  TEST_F(MetaVariableTest, TransformedNetworkOutput)
5073  {
5074  // check that garbage input throws helpful B2FATAL
5075  EXPECT_B2FATAL(Manager::Instance().getVariable("transformedNetworkOutput(NONSENSE)"));
5076 
5077  // check that helpful B2FATAL is thrown if second or third argument is not a double
5078  EXPECT_B2FATAL(Manager::Instance().getVariable("transformedNetworkOutput(NONEXISTANT, 0, NOTDOUBLE)"));
5079  EXPECT_B2FATAL(Manager::Instance().getVariable("transformedNetworkOutput(NONEXISTANT, NOTDOUBLE, 1)"));
5080 
5081  // check for NaN if network output variable does not exist (no matter whether particle is provided or not)
5082  StoreArray<Particle> myParticles;
5083  const Particle* particle = myParticles.appendNew();
5084  const Manager::Var* var = Manager::Instance().getVariable("transformedNetworkOutput(NONEXISTANT, 0, 1)");
5085  EXPECT_TRUE(std::isnan(var->function(particle)));
5086  StoreObjPtr<EventExtraInfo> eventExtraInfo;
5087  if (not eventExtraInfo.isValid())
5088  eventExtraInfo.create();
5089  var = Manager::Instance().getVariable("transformedNetworkOutput(NONEXISTANT, 0, 1)");
5090  EXPECT_TRUE(std::isnan(var->function(nullptr)));
5091  }
5092 }
Belle2::EvtPDLUtil::charge
double charge(int pdgCode)
Returns electric charge of a particle with given pdg code.
Definition: EvtPDLUtil.cc:46
Belle2::StoreArray::appendNew
T * appendNew()
Construct a new T object at the end of the array.
Definition: StoreArray.h:256
Belle2::MCParticleGraph::generateList
void generateList(const std::string &name="", int options=c_setNothing)
Generates the MCParticle list and stores it in the StoreArray with the given name.
Definition: MCParticleGraph.cc:211
Belle2::MCParticle::getEnergy
float getEnergy() const
Return particle energy in GeV.
Definition: MCParticle.h:158
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::RestOfEvent::updateMaskWithCuts
void updateMaskWithCuts(const std::string &name, const std::shared_ptr< Variable::Cut > &trackCut=nullptr, const std::shared_ptr< Variable::Cut > &eclCut=nullptr, const std::shared_ptr< Variable::Cut > &klmCut=nullptr, bool updateExisting=false)
Update mask with cuts.
Definition: RestOfEvent.cc:180
Belle2::V0
Object holding information for V0s.
Definition: V0.h:40
Belle2::ECLCluster
ECL cluster data.
Definition: ECLCluster.h:39
Belle2::MCParticleGraph
Class to build, validate and sort a particle decay chain.
Definition: MCParticleGraph.h:48
Belle2::Particle::addExtraInfo
void addExtraInfo(const std::string &name, float value)
Sets the user-defined data of given name to the given value.
Definition: Particle.cc:1224
Belle2::DataStore::EStoreFlags
EStoreFlags
Flags describing behaviours of objects etc.
Definition: DataStore.h:71
Belle2::PCmsLabTransform::getCMSEnergy
double getCMSEnergy() const
Returns CMS energy of e+e- (aka.
Definition: PCmsLabTransform.h:57
Belle2::PIDLikelihood::setLogLikelihood
void setLogLikelihood(Const::EDetector det, const Const::ChargedStable &part, float logl)
Set log likelihood for a given detector and particle.
Definition: PIDLikelihood.cc:35
Belle2::MCParticle::setDecayVertex
void setDecayVertex(const TVector3 &vertex)
Set decay vertex.
Definition: MCParticle.h:452
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::MCParticle::setStatus
void setStatus(unsigned short int status)
Set Status code for the particle.
Definition: MCParticle.h:354
Belle2::RestOfEvent
This is a general purpose class for collecting reconstructed MDST data objects that are not used in r...
Definition: RestOfEvent.h:59
Belle2::PIDLikelihood
Class to collect log likelihoods from TOP, ARICH, dEdx, ECL and KLM aimed for output to mdst includes...
Definition: PIDLikelihood.h:37
Belle2::RelationsInterface::getRelatedTo
TO * getRelatedTo(const std::string &name="", const std::string &namedRelation="") const
Get the object to which this object has a relation.
Definition: RelationsObject.h:250
Belle2::RestOfEvent::hasParticle
bool hasParticle(const Particle *particle, const std::string &maskName="") const
Check if ROE has StoreArray index of given to the list of unused tracks in the event.
Definition: RestOfEvent.cc:125
Belle2::ECLCluster::setEnergy
void setEnergy(double energy)
Set Corrected Energy (GeV).
Definition: ECLCluster.h:238
Belle2::Particle::getDaughter
const Particle * getDaughter(unsigned i) const
Returns a pointer to the i-th daughter particle.
Definition: Particle.cc:595
Belle2::Particle::getP
float getP() const
Returns momentum magnitude (same as getMomentumMagnitude but with shorter name)
Definition: Particle.h:493
Belle2::Variable::Manager::Var::function
FunctionPtr function
Pointer to function.
Definition: Manager.h:138
Belle2::ECLCluster::getEnergy
double getEnergy(const EHypothesisBit &hypothesis) const
Return Energy (GeV).
Definition: ECLCluster.cc:21
Belle2::PCmsLabTransform::rotateCmsToLab
const TLorentzRotation rotateCmsToLab() const
Returns Lorentz transformation from CMS to Lab.
Definition: PCmsLabTransform.h:84
Belle2::Particle::getEnergy
float getEnergy() const
Returns total energy.
Definition: Particle.h:455
Belle2::MCParticle::setProductionTime
void setProductionTime(float time)
Set production time.
Definition: MCParticle.h:392
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::Particle::getMass
float getMass() const
Returns invariant mass (= nominal for FS particles)
Definition: Particle.h:433
Belle2::StoreArray::clear
void clear() override
Delete all entries in this array.
Definition: StoreArray.h:217
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::Particle::appendDaughter
void appendDaughter(const Particle *daughter, const bool updateType=true)
Appends index of daughter to daughters index array.
Definition: Particle.cc:633
Belle2::MCParticleGraph::addParticle
GraphParticle & addParticle()
Add new particle to the graph.
Definition: MCParticleGraph.h:305
Belle2::TEST_F
TEST_F(GlobalLabelTest, LargeNumberOfTimeDependentParameters)
Test large number of time-dep params for registration and retrieval.
Definition: globalLabel.cc:65
Belle2::MCParticle::setProductionVertex
void setProductionVertex(const TVector3 &vertex)
Set production vertex position.
Definition: MCParticle.h:404
Belle2::UseReferenceFrame
A template class to apply the reference frame.
Definition: ReferenceFrame.h:349
Belle2::ECLCluster::setClusterId
void setClusterId(int clusterid)
Set cluster id.
Definition: ECLCluster.h:156
Belle2::MCParticle::getMass
float getMass() const
Return the particle mass in GeV.
Definition: MCParticle.h:146
Belle2::Particle::setVertex
void setVertex(const TVector3 &vertex)
Sets position (decay vertex)
Definition: Particle.h:291
Belle2::RelationsInterface::getArrayIndex
int getArrayIndex() const
Returns this object's array index (in StoreArray), or -1 if not found.
Definition: RelationsObject.h:387
Belle2::KLMCluster
KLM cluster data.
Definition: KLMCluster.h:38
Belle2::ECLCluster::setR
void setR(double r)
Set R (in cm).
Definition: ECLCluster.h:235
Belle2::TEST
TEST(TestgetDetectorRegion, TestgetDetectorRegion)
Test Constructors.
Definition: utilityFunctions.cc:18
Belle2::RestOfEvent::initializeMask
void initializeMask(const std::string &name, const std::string &origin="unknown")
Initialize new mask.
Definition: RestOfEvent.cc:135
Belle2::Particle
Class to store reconstructed particles.
Definition: Particle.h:77
Belle2::PCmsLabTransform
Class to hold Lorentz transformations from/to CMS and boost vector.
Definition: PCmsLabTransform.h:37
Belle2::PCmsLabTransform::rotateLabToCms
const TLorentzRotation rotateLabToCms() const
Returns Lorentz transformation from Lab to CMS.
Definition: PCmsLabTransform.h:74
Belle2::Const::ChargedStable
Provides a type-safe way to pass members of the chargedStableSet set.
Definition: Const.h:465
Belle2::MCParticle::set4Vector
void set4Vector(const TLorentzVector &p4)
Sets the 4Vector of particle.
Definition: MCParticle.h:446
Belle2::Particle::getPDGMass
float getPDGMass(void) const
Returns uncertaint on the invariant mass (requiers valid momentum error matrix)
Definition: Particle.cc:577
Belle2::Track
Class that bundles various TrackFitResults.
Definition: Track.h:35
Belle2::MCParticleGraph::GraphParticle::comesFrom
void comesFrom(GraphParticle &mother)
Tells the graph that this particle is a decay product of mother.
Definition: MCParticleGraph.h:116
Belle2::MCParticle
A Class to store the Monte Carlo particle information.
Definition: MCParticle.h:43
Belle2::StoreArray< Particle >
Belle2::Particle::get4Vector
TLorentzVector get4Vector() const
Returns Lorentz vector.
Definition: Particle.h:464
Belle2::RestOfEvent::addParticles
void addParticles(const std::vector< const Particle * > &particle)
Add StoreArray indices of given Particles to the list of unused particles in the event.
Definition: RestOfEvent.cc:27
Belle2::MCParticle::setDecayTime
void setDecayTime(float time)
Set decay time.
Definition: MCParticle.h:398
Belle2::RestOfEvent::updateMaskWithV0
void updateMaskWithV0(const std::string &name, const Particle *particleV0)
Update mask with composite particle.
Definition: RestOfEvent.cc:216
Belle2::MCParticle::setMomentum
void setMomentum(const TVector3 &momentum)
Set particle momentum.
Definition: MCParticle.h:425
Belle2::StoreArray::create
bool create(bool replace=false)
Creating StoreArrays is unnecessary, only used internally.
Definition: StoreArray.h:339
Belle2::MCParticleGraph::GraphParticle
Class to represent Particle data in graph.
Definition: MCParticleGraph.h:86
Belle2::ECLCluster::setTheta
void setTheta(double theta)
Set Theta of Shower (radian).
Definition: ECLCluster.h:229