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