Belle II Software  light-2205-abys
clusteringVariables.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 // Unit tests for all variables related to clustering and clustering
10 // subdetectors (KLM and ECL variables, track <--> cluster matching etc)
11 
12 #include <gtest/gtest.h>
13 #include <TRandom3.h>
14 #include <Math/Cartesian2D.h>
15 
16 #include <analysis/ParticleCombiner/ParticleCombiner.h>
17 
18 // VariableManager and particle(list)
19 #include <analysis/VariableManager/Manager.h>
20 #include <analysis/dataobjects/Particle.h>
21 #include <analysis/dataobjects/ParticleList.h>
22 
23 // mdst dataobjects
24 #include <mdst/dataobjects/ECLCluster.h>
25 #include <mdst/dataobjects/KLMCluster.h>
26 #include <mdst/dataobjects/Track.h>
27 
28 // framework - set up mock events
29 #include <framework/datastore/StoreArray.h>
30 #include <framework/datastore/StoreObjPtr.h>
31 #include <framework/utilities/TestHelpers.h>
32 #include <framework/gearbox/Const.h>
33 #include <framework/gearbox/Gearbox.h>
34 
35 using namespace Belle2;
36 using namespace Belle2::Variable;
37 
38 namespace {
39 
41  class ECLVariableTest : public ::testing::Test {
42  protected:
44  void SetUp() override
45  {
46  // setup the DataStore
48 
49  // particles (to be filled)
50  StoreArray<Particle> particles;
51  particles.registerInDataStore();
52 
53  // mock up mdst objects
54  StoreArray<Track> tracks;
55  tracks.registerInDataStore();
57  trackFits.registerInDataStore();
58  StoreArray<ECLCluster> eclclusters;
59  eclclusters.registerInDataStore();
60 
61  // tracks can be matched to clusters
62  tracks.registerRelationTo(eclclusters);
63 
64  // we're done setting up the datastore
66 
67  // add some tracks the zeroth one is not going to be matched
68  tracks.appendNew(Track());
69  const Track* t1 = tracks.appendNew(Track());
70  const Track* t2 = tracks.appendNew(Track());
71  const Track* t3 = tracks.appendNew(Track());
72  const Track* t4 = tracks.appendNew(Track());
73  tracks.appendNew(Track());
74  tracks.appendNew(Track());
75 
76  // mock up some TrackFits for them (all pions)
77  TRandom3 generator;
78  TMatrixDSym cov6(6);
79  auto CDCValue = static_cast<unsigned long long int>(0x300000000000000);
80 
81  for (int i = 0; i < tracks.getEntries(); ++i) {
82  int charge = (i % 2 == 0) ? +1 : -1;
83  ROOT::Math::Cartesian2D d(generator.Uniform(-1, 1), generator.Uniform(-1, 1));
84  ROOT::Math::Cartesian2D pt(generator.Uniform(-1, 1), generator.Uniform(-1, 1));
85  d.SetXY(d.X(), -(d.X()*pt.x()) / pt.y());
86  B2Vector3D position(d.X(), d.Y(), generator.Uniform(-1, 1));
87  B2Vector3D momentum(pt.x(), pt.y(), generator.Uniform(-1, 1));
88  trackFits.appendNew(position, momentum, cov6, charge, Const::pion, 0.5, 1.5, CDCValue, 16777215, 0);
89  tracks[i]->setTrackFitResultIndex(Const::pion, i);
90  }
91 
92  // add some ECL clusters
93  ECLCluster* e1 = eclclusters.appendNew(ECLCluster());
94  e1->setEnergy(0.3);
96  e1->setClusterId(1);
97  e1->setTime(1);
98  e1->setDeltaTime99(0.1);
99  e1->setConnectedRegionId(1);
100  // leave this guy with default theta and phi
101  ECLCluster* e2 = eclclusters.appendNew(ECLCluster());
102  e2->setEnergy(0.6);
103  e2->setTheta(1.0); // somewhere in the barrel
104  e2->setPhi(2.0);
105  e2->setR(148.5);
107  e2->setClusterId(2);
108  e2->setTime(2);
109  e2->setDeltaTime99(0.2);
110  e2->setConnectedRegionId(2);
111  ECLCluster* e3 = eclclusters.appendNew(ECLCluster());
112  e3->setEnergy(0.15);
113  e3->setTheta(0.2); // somewhere in the fwd endcap
114  e3->setPhi(1.5);
115  e3->setR(200.0);
118  // let's suppose this cluster could also be due to a neutral hadron. In
119  // this case, the c_neutralHadron hypothesis bit would hopefully also have
120  // been set by the reconstruction... arbitrarily choose cluster 3
121  e3->setClusterId(3);
122  e3->setTime(3);
123  e3->setDeltaTime99(0.3);
124  e3->setConnectedRegionId(1); // shares the connected region with cluster 1
125 
126  // aaand add clusters related to the tracks
127  ECLCluster* e4 = eclclusters.appendNew(ECLCluster());
128  e4->setEnergy(0.2);
130  e4->setClusterId(4);
131  t1->addRelationTo(e4);
132  e4->setIsTrack(true);
133 
134  ECLCluster* e5 = eclclusters.appendNew(ECLCluster());
135  e5->setEnergy(0.3);
137  e5->setClusterId(5);
138  t2->addRelationTo(e5);
139  e5->setIsTrack(true);
140 
141  ECLCluster* e6 = eclclusters.appendNew(ECLCluster());
142  e6->setEnergy(0.2);
144  e6->setClusterId(6);
145  t3->addRelationTo(e6);
146  t4->addRelationTo(e6);
147  // two tracks are related to this cluster this can happen due to real
148  // physics and we should be able to cope
149  e6->setIsTrack(true);
150 
151  }
152 
154  void TearDown() override
155  {
157  }
158  };
159 
160  TEST_F(ECLVariableTest, b2bKinematicsTest)
161  {
162  // we need the particles and ECLClusters arrays
163  StoreArray<Particle> particles;
164  StoreArray<ECLCluster> eclclusters;
165  StoreArray<Track> tracks;
166 
167  // connect gearbox for CMS boosting etc
168  Gearbox& gearbox = Gearbox::getInstance();
169  gearbox.setBackends({std::string("file:")});
170  gearbox.close();
171  gearbox.open("geometry/Belle2.xml", false);
172 
173  // register in the datastore
174  StoreObjPtr<ParticleList> gammalist("gamma:testGammaAllList");
176  gammalist.registerInDataStore(DataStore::c_DontWriteOut);
178 
179  // initialise the lists
180  gammalist.create();
181  gammalist->initialize(22, gammalist.getName());
182 
183  // make the photons from clusters
184  for (int i = 0; i < eclclusters.getEntries(); ++i) {
185  if (!eclclusters[i]->isTrack()) {
186  const Particle* p = particles.appendNew(Particle(eclclusters[i]));
187  gammalist->addParticle(p);
188  }
189  }
190 
191  // get the zeroth track in the array (is not associated to a cluster)
192  const Particle* noclustertrack = particles.appendNew(Particle(tracks[0], Const::pion));
193 
194  // grab variables for testing
195  const Manager::Var* b2bClusterTheta = Manager::Instance().getVariable("b2bClusterTheta");
196  const Manager::Var* b2bClusterPhi = Manager::Instance().getVariable("b2bClusterPhi");
197 
198  EXPECT_EQ(gammalist->getListSize(), 3);
199 
200  EXPECT_FLOAT_EQ(std::get<double>(b2bClusterTheta->function(gammalist->getParticle(0))), 3.0276606);
201  EXPECT_FLOAT_EQ(std::get<double>(b2bClusterPhi->function(gammalist->getParticle(0))), 0.0);
202  EXPECT_FLOAT_EQ(std::get<double>(b2bClusterTheta->function(gammalist->getParticle(1))), 1.6036042);
203  EXPECT_FLOAT_EQ(std::get<double>(b2bClusterPhi->function(gammalist->getParticle(1))), -1.0607308);
204  EXPECT_FLOAT_EQ(std::get<double>(b2bClusterTheta->function(gammalist->getParticle(2))), 2.7840068);
205  EXPECT_FLOAT_EQ(std::get<double>(b2bClusterPhi->function(gammalist->getParticle(2))), -1.3155469);
206 
207  // track (or anything without a cluster) should be nan
208  ASSERT_TRUE(std::isnan(std::get<double>(b2bClusterTheta->function(noclustertrack))));
209  ASSERT_TRUE(std::isnan(std::get<double>(b2bClusterPhi->function(noclustertrack))));
210 
211  // the "normal" (not cluster based) variables should be the same for photons
212  // (who have no track information)
213  const Manager::Var* b2bTheta = Manager::Instance().getVariable("b2bTheta");
214  const Manager::Var* b2bPhi = Manager::Instance().getVariable("b2bPhi");
215 
216  EXPECT_FLOAT_EQ(std::get<double>(b2bClusterTheta->function(gammalist->getParticle(0))),
217  std::get<double>(b2bTheta->function(gammalist->getParticle(0))));
218  EXPECT_FLOAT_EQ(std::get<double>(b2bClusterPhi->function(gammalist->getParticle(0))),
219  std::get<double>(b2bPhi->function(gammalist->getParticle(0))));
220  }
221 
222  TEST_F(ECLVariableTest, clusterKinematicsTest)
223  {
224  // we need the particles and ECLClusters arrays
225  StoreArray<Particle> particles;
226  StoreArray<ECLCluster> eclclusters;
227  StoreArray<Track> tracks;
228 
229  // connect gearbox for CMS boosting etc
230  Gearbox& gearbox = Gearbox::getInstance();
231  gearbox.setBackends({std::string("file:")});
232  gearbox.close();
233  gearbox.open("geometry/Belle2.xml", false);
234 
235  // register in the datastore
236  StoreObjPtr<ParticleList> gammalist("gamma:testGammaAllList");
238  gammalist.registerInDataStore(DataStore::c_DontWriteOut);
240 
241  // initialise the lists
242  gammalist.create();
243  gammalist->initialize(22, gammalist.getName());
244 
245  // make the photons from clusters
246  for (int i = 0; i < eclclusters.getEntries(); ++i) {
247  if (!eclclusters[i]->isTrack()) {
248  const Particle* p = particles.appendNew(Particle(eclclusters[i]));
249  gammalist->addParticle(p);
250  }
251  }
252 
253  // grab variables for testing
254  const Manager::Var* clusterPhi = Manager::Instance().getVariable("clusterPhi");
255  const Manager::Var* clusterPhiCMS = Manager::Instance().getVariable("useCMSFrame(clusterPhi)");
256  const Manager::Var* clusterTheta = Manager::Instance().getVariable("clusterTheta");
257  const Manager::Var* clusterThetaCMS = Manager::Instance().getVariable("useCMSFrame(clusterTheta)");
258 
259  EXPECT_FLOAT_EQ(std::get<double>(clusterPhi->function(gammalist->getParticle(1))), 2.0);
260  EXPECT_FLOAT_EQ(std::get<double>(clusterPhiCMS->function(gammalist->getParticle(1))), 2.042609);
261  EXPECT_FLOAT_EQ(std::get<double>(clusterTheta->function(gammalist->getParticle(1))), 1.0);
262  EXPECT_FLOAT_EQ(std::get<double>(clusterThetaCMS->function(gammalist->getParticle(1))), 1.2599005);
263 
264  // test cluster quantities directly (lab system only)
265  EXPECT_FLOAT_EQ(std::get<double>(clusterPhi->function(gammalist->getParticle(0))), eclclusters[0]->getPhi());
266  EXPECT_FLOAT_EQ(std::get<double>(clusterTheta->function(gammalist->getParticle(0))), eclclusters[0]->getTheta());
267  }
268 
269  TEST_F(ECLVariableTest, HypothesisVariables)
270  {
271  // we need the particles and ECLClusters arrays
272  StoreArray<Particle> particles;
273  StoreArray<ECLCluster> eclclusters;
274 
275  // register in the datastore
276  StoreObjPtr<ParticleList> gammalist("gamma");
278  gammalist.registerInDataStore(DataStore::c_DontWriteOut);
280 
281  // initialise the lists
282  gammalist.create();
283  gammalist->initialize(22, gammalist.getName());
284 
285  // make the photons from clusters
286  for (int i = 0; i < eclclusters.getEntries(); ++i)
287  if (!eclclusters[i]->isTrack()) {
288  const Particle* p = particles.appendNew(Particle(eclclusters[i]));
289  gammalist->addParticle(p);
290  }
291 
292  // grab variables for testing
293  const Manager::Var* vHasNPhotons = Manager::Instance().getVariable("clusterHasNPhotons");
294  const Manager::Var* vHasNeutHadr = Manager::Instance().getVariable("clusterHasNeutralHadron");
295 
296  // check that the hypotheses are correctly propagated to the VM.
297  for (size_t i = 0; i < gammalist->getListSize(); ++i) {
298  EXPECT_FLOAT_EQ(std::get<double>(vHasNPhotons->function(gammalist->getParticle(i))), 1.0);
299  if (i == 2) { // third cluster arbitrarily chosen to test the behaviour of dual hypothesis clusters
300  EXPECT_FLOAT_EQ(std::get<double>(vHasNeutHadr->function(gammalist->getParticle(i))), 1.0);
301  } else {
302  EXPECT_FLOAT_EQ(std::get<double>(vHasNeutHadr->function(gammalist->getParticle(i))), 0.0);
303  }
304  } // end loop over test list
305  }
306 
307  TEST_F(ECLVariableTest, IsFromECL)
308  {
309  StoreArray<Particle> particles;
310  StoreArray<ECLCluster> eclclusters;
311 
312  const Manager::Var* vIsFromECL = Manager::Instance().getVariable("isFromECL");
313  const Manager::Var* vIsFromKLM = Manager::Instance().getVariable("isFromKLM");
314  const Manager::Var* vIsFromTrack = Manager::Instance().getVariable("isFromTrack");
315  const Manager::Var* vIsFromV0 = Manager::Instance().getVariable("isFromV0");
316 
317  for (int i = 0; i < eclclusters.getEntries(); ++i)
318  if (!eclclusters[i]->isTrack()) {
319  const Particle* p = particles.appendNew(Particle(eclclusters[i]));
320  EXPECT_TRUE(std::get<bool>(vIsFromECL->function(p)));
321  EXPECT_FALSE(std::get<bool>(vIsFromKLM->function(p)));
322  EXPECT_FALSE(std::get<bool>(vIsFromTrack->function(p)));
323  EXPECT_FALSE(std::get<bool>(vIsFromV0->function(p)));
324  }
325  }
326 
327  TEST_F(ECLVariableTest, ECLThetaAndPhiId)
328  {
329  StoreArray<Particle> particles;
330  StoreArray<ECLCluster> clusters{};
331  StoreArray<ECLCluster> eclclusters;
332  // make a particle from cluster #1
333  const Particle* p = particles.appendNew(Particle(eclclusters[0]));
334 
335  // get the variables to test
336  const Manager::Var* clusterThetaID = Manager::Instance().getVariable("clusterThetaID");
337  const Manager::Var* clusterPhiID = Manager::Instance().getVariable("clusterPhiID");
338 
339  {
340  clusters[0]->setMaxECellId(1);
341  EXPECT_FLOAT_EQ(std::get<double>(clusterThetaID->function(p)), 0);
342  EXPECT_FLOAT_EQ(std::get<double>(clusterPhiID->function(p)), 0);
343  }
344  {
345  clusters[0]->setMaxECellId(6903);
346  EXPECT_FLOAT_EQ(std::get<double>(clusterThetaID->function(p)), 52);
347  EXPECT_FLOAT_EQ(std::get<double>(clusterPhiID->function(p)), 134);
348  }
349  {
350  clusters[0]->setMaxECellId(8457);
351  EXPECT_FLOAT_EQ(std::get<double>(clusterThetaID->function(p)), 65);
352  EXPECT_FLOAT_EQ(std::get<double>(clusterPhiID->function(p)), 8);
353  }
354  }
355 
356 
357  TEST_F(ECLVariableTest, WholeEventClosure)
358  {
359  // we need the particles, tracks, and ECLClusters StoreArrays
360  StoreArray<Particle> particles;
361  StoreArray<Track> tracks;
362  StoreArray<ECLCluster> eclclusters;
363 
364  // create a photon (clusters) and pion (tracks) lists
365  StoreObjPtr<ParticleList> gammalist("gamma:testGammaAllList");
366  StoreObjPtr<ParticleList> pionslist("pi+:testPionAllList");
367  StoreObjPtr<ParticleList> apionslist("pi-:testPionAllList");
368 
369  // register the lists in the datastore
371  gammalist.registerInDataStore(DataStore::c_DontWriteOut);
372  pionslist.registerInDataStore(DataStore::c_DontWriteOut);
373  apionslist.registerInDataStore(DataStore::c_DontWriteOut);
375 
376  // initialise the lists
377  gammalist.create();
378  gammalist->initialize(22, gammalist.getName());
379  pionslist.create();
380  pionslist->initialize(211, pionslist.getName());
381  apionslist.create();
382  apionslist->initialize(-211, apionslist.getName());
383  apionslist->bindAntiParticleList(*(pionslist));
384 
385  // make the photons from clusters (and sum up the total ecl energy)
386  double eclEnergy = 0.0;
387  for (int i = 0; i < eclclusters.getEntries(); ++i) {
388  eclEnergy += eclclusters[i]->getEnergy(ECLCluster::EHypothesisBit::c_nPhotons);
389  if (!eclclusters[i]->isTrack()) {
390  const Particle* p = particles.appendNew(Particle(eclclusters[i]));
391  gammalist->addParticle(p);
392  }
393  }
394 
395 
396  // make the pions from tracks
397  for (int i = 0; i < tracks.getEntries(); ++i) {
398  const Particle* p = particles.appendNew(Particle(tracks[i], Const::pion));
399  pionslist->addParticle(p);
400  }
401 
402  // grab variables
403  const Manager::Var* vClusterE = Manager::Instance().getVariable("clusterE");
404  const Manager::Var* vClNTrack = Manager::Instance().getVariable("nECLClusterTrackMatches");
405 
406  // calculate the total neutral energy from the particle list --> VM
407  double totalNeutralClusterE = 0.0;
408  for (size_t i = 0; i < gammalist->getListSize(); ++i)
409  totalNeutralClusterE += std::get<double>(vClusterE->function(gammalist->getParticle(i)));
410 
411  // calculate the total track-matched cluster energy from the particle list --> VM
412  double totalTrackClusterE = 0.0;
413  for (size_t i = 0; i < pionslist->getListSize(); ++i) { // includes antiparticles
414  double clusterE = std::get<double>(vClusterE->function(pionslist->getParticle(i)));
415  double nOtherCl = std::get<double>(vClNTrack->function(pionslist->getParticle(i)));
416  if (nOtherCl > 0)
417  totalTrackClusterE += clusterE / nOtherCl;
418  }
419 
420  EXPECT_FLOAT_EQ(totalNeutralClusterE + totalTrackClusterE, eclEnergy);
421  }
422 
423  TEST_F(ECLVariableTest, eclClusterOnlyInvariantMass)
424  {
425  // declare all the array we need
426  StoreArray<Particle> particles;
427  std::vector<int> daughterIndices, daughterIndices_noclst;
428 
429  //proxy initialize where to declare the needed array
431  StoreArray<ECLCluster> eclclusters_new;
432  eclclusters_new.registerInDataStore();
433  particles.registerRelationTo(eclclusters_new);
435 
436  // create two Lorentz vectors
437  const float px_0 = 2.;
438  const float py_0 = 1.;
439  const float pz_0 = 3.;
440  const float px_1 = 1.5;
441  const float py_1 = 1.5;
442  const float pz_1 = 2.5;
443  float E_0, E_1;
444  E_0 = sqrt(pow(px_0, 2) + pow(py_0, 2) + pow(pz_0, 2));
445  E_1 = sqrt(pow(px_1, 2) + pow(py_1, 2) + pow(pz_1, 2));
446  ROOT::Math::PxPyPzEVector momentum;
447  ROOT::Math::PxPyPzEVector dau0_4vec(px_0, py_0, pz_0, E_0), dau1_4vec(px_1, py_1, pz_1, E_1);
448 
449  // add the two photons as the two daughters of some particle and create the latter
450  Particle dau0_noclst(dau0_4vec, 22);
451  momentum += dau0_noclst.get4Vector();
452  Particle* newDaughter0_noclst = particles.appendNew(dau0_noclst);
453  daughterIndices_noclst.push_back(newDaughter0_noclst->getArrayIndex());
454  Particle dau1_noclst(dau1_4vec, 22);
455  momentum += dau1_noclst.get4Vector();
456  Particle* newDaughter1_noclst = particles.appendNew(dau1_noclst);
457  daughterIndices_noclst.push_back(newDaughter1_noclst->getArrayIndex());
458  const Particle* par_noclst = particles.appendNew(momentum, 111, Particle::c_Unflavored, daughterIndices_noclst);
459 
460  // grab variables
461  const Manager::Var* var = Manager::Instance().getVariable("eclClusterOnlyInvariantMass");
462 
463  // when no relations are set between the particles and the eclClusters, nan is expected to be returned
464  ASSERT_NE(var, nullptr);
465  EXPECT_TRUE(std::isnan(std::get<double>(var->function(par_noclst))));
466 
467  // set relations between particles and eclClusters
468  ECLCluster* eclst0 = eclclusters_new.appendNew(ECLCluster());
469  eclst0->setEnergy(dau0_4vec.E());
471  eclst0->setClusterId(1);
472  eclst0->setTheta(dau0_4vec.Theta());
473  eclst0->setPhi(dau0_4vec.Phi());
474  eclst0->setR(148.4);
475  ECLCluster* eclst1 = eclclusters_new.appendNew(ECLCluster());
476  eclst1->setEnergy(dau1_4vec.E());
478  eclst1->setClusterId(2);
479  eclst1->setTheta(dau1_4vec.Theta());
480  eclst1->setPhi(dau1_4vec.Phi());
481  eclst1->setR(148.5);
482 
483  // use these new-created clusters rather than the 6 default ones
484  const Particle* newDaughter0 = particles.appendNew(Particle(eclclusters_new[6]));
485  daughterIndices.push_back(newDaughter0->getArrayIndex());
486  const Particle* newDaughter1 = particles.appendNew(Particle(eclclusters_new[7]));
487  daughterIndices.push_back(newDaughter1->getArrayIndex());
488 
489  const Particle* par = particles.appendNew(momentum, 111, Particle::c_Unflavored, daughterIndices);
490 
491  //now we expect non-nan results
492  EXPECT_FLOAT_EQ(std::get<double>(var->function(par)), 0.73190731);
493  }
494 
495  TEST_F(ECLVariableTest, averageECLTimeQuantities)
496  {
497  // we need the particles, and ECLClusters StoreArrays
498  StoreArray<Particle> particles;
499  StoreArray<ECLCluster> eclclusters;
500 
501  // create a photon (clusters) and a pi0 list
502  StoreObjPtr<ParticleList> gammalist("gamma:testGammaAllList");
503  StoreObjPtr<ParticleList> pionslist("pi0:testPizAllList");
504 
505  // register the lists in the datastore
507  gammalist.registerInDataStore(DataStore::c_DontWriteOut);
508  pionslist.registerInDataStore(DataStore::c_DontWriteOut);
510 
511  // initialise the lists
512  gammalist.create();
513  gammalist->initialize(22, gammalist.getName());
514  pionslist.create();
515  pionslist->initialize(111, pionslist.getName());
516 
517  // make the photons from clusters
518  for (int i = 0; i < eclclusters.getEntries(); ++i) {
519  if (!eclclusters[i]->isTrack()) {
520  const Particle* p = particles.appendNew(Particle(eclclusters[i]));
521  gammalist->addParticle(p);
522  }
523  }
524 
525  // make the pi0s as combinations of photons
526  ParticleGenerator combiner_pi0("pi0:testPizAllList -> gamma:testGammaAllList gamma:testGammaAllList");
527  combiner_pi0.init();
528  while (combiner_pi0.loadNext()) {
529  const Particle& particle = combiner_pi0.getCurrentParticle();
530 
531  particles.appendNew(particle);
532  int iparticle = particles.getEntries() - 1;
533 
534  pionslist->addParticle(iparticle, particle.getPDGCode(), particle.getFlavorType());
535  }
536 
537  // grab variables
538  const Manager::Var* weightedAverageECLTime = Manager::Instance().getVariable("weightedAverageECLTime");
539  const Manager::Var* maxDist = Manager::Instance().getVariable("maxWeightedDistanceFromAverageECLTime");
540 
541  // particles without daughters should have NaN as weighted average
542  EXPECT_TRUE(std::isnan(std::get<double>(weightedAverageECLTime->function(gammalist->getParticle(0)))));
543 
544  // check that weighted average of first pi0 is correct
545  EXPECT_FLOAT_EQ(std::get<double>(weightedAverageECLTime->function(pionslist->getParticle(0))), 1.2);
546 
547  // check that maximal difference to weighted average in units of uncertainty is calculated correctly
548  EXPECT_FLOAT_EQ(std::get<double>(maxDist->function(pionslist->getParticle(0))), 4.0);
549  }
550 
551  TEST_F(ECLVariableTest, photonHasOverlap)
552  {
553  // declare StoreArrays of Particles and ECLClusters
554  StoreArray<Particle> particles;
555  StoreArray<ECLCluster> eclclusters;
556  StoreArray<Track> tracks;
557 
558  // create a photon and the pion lists
559  StoreObjPtr<ParticleList> gammalist("gamma:all");
560  StoreObjPtr<ParticleList> pionlist("pi+:all");
561  StoreObjPtr<ParticleList> piminuslist("pi-:all");
562 
563  // register the particle lists in the datastore
565  gammalist.registerInDataStore(DataStore::c_DontWriteOut);
566  pionlist.registerInDataStore(DataStore::c_DontWriteOut);
567  piminuslist.registerInDataStore(DataStore::c_DontWriteOut);
569 
570  // initialise the photon list
571  gammalist.create();
572  gammalist->initialize(22, gammalist.getName());
573 
574  // make the photons from clusters
575  for (int i = 0; i < eclclusters.getEntries(); ++i) {
576  if (!eclclusters[i]->isTrack()) {
577  const Particle* p = particles.appendNew(Particle(eclclusters[i]));
578  gammalist->addParticle(p);
579  }
580  }
581 
582  // initialize the pion lists
583  pionlist.create();
584  pionlist->initialize(211, pionlist.getName());
585  piminuslist.create();
586  piminuslist->initialize(-211, piminuslist.getName());
587  piminuslist->bindAntiParticleList(*(pionlist));
588 
589  // fill the pion list with the tracks
590  for (int i = 0; i < tracks.getEntries(); ++i) {
591  const Particle* p = particles.appendNew(Particle(tracks[i], Const::pion));
592  pionlist->addParticle(p);
593  }
594 
595  // check overlap without any arguments
596  const Manager::Var* photonHasOverlapNoArgs = Manager::Instance().getVariable("photonHasOverlap()");
597  // the track list e-:all is missing so NaN is returned
598  EXPECT_TRUE(std::isnan(std::get<double>(photonHasOverlapNoArgs->function(particles[0]))));
599 
600  // check overlap without any requirement on other photons
601  const Manager::Var* photonHasOverlapAll = Manager::Instance().getVariable("photonHasOverlap(, gamma:all, pi+:all)");
602  // cluster 3 and cluster 1 share the connected region
603  EXPECT_TRUE(std::get<double>(photonHasOverlapAll->function(particles[0])));
604  // photonHasOverlap is designed for photons, so calling it for a pion returns NaN
605  EXPECT_TRUE(std::isnan(std::get<double>(photonHasOverlapAll->function(particles[3]))));
606 
607  // check overlap with photons in barrel
608  const Manager::Var* photonHasOverlapBarrel = Manager::Instance().getVariable("photonHasOverlap(clusterReg==2, gamma:all, pi+:all)");
609  // cluster 3 is in the forward end cap so it doesn't matter that it has the same connected region like cluster 1
610  EXPECT_FALSE(std::get<double>(photonHasOverlapBarrel->function(particles[0])));
611  }
612 
613  class KLMVariableTest : public ::testing::Test {
614  protected:
616  void SetUp() override
617  {
618  // setup the DataStore
620 
621  // particles (to be filled)
622  StoreArray<Particle> particles;
623  particles.registerInDataStore();
624 
625  // mock up mdst objects
626  StoreArray<Track> tracks;
627  tracks.registerInDataStore();
628  StoreArray<TrackFitResult> trackFits;
629  trackFits.registerInDataStore();
630  StoreArray<KLMCluster> klmClusters;
631  klmClusters.registerInDataStore();
632 
633  // tracks can be matched to clusters
634  tracks.registerRelationTo(klmClusters);
635 
636  // we're done setting up the datastore
638  }
639 
641  void TearDown() override
642  {
644  }
645  };
646 
647  TEST_F(KLMVariableTest, WholeEventClosure)
648  {
649  // we need the Particles, Tracks, TrackFitResults and KLMClusters StoreArrays
650  StoreArray<Particle> particles;
651  StoreArray<Track> tracks;
652  StoreArray<TrackFitResult> trackFits;
653  StoreArray<KLMCluster> klmClusters;
654 
655  // create a KLong (clusters) and muon (tracks) lists
656  StoreObjPtr<ParticleList> kLongList("K0_L:testKLong");
657  StoreObjPtr<ParticleList> muonsList("mu-:testMuons");
658  StoreObjPtr<ParticleList> amuonsList("mu+:testMuons");
659 
660  // register the lists in the datastore
662  kLongList.registerInDataStore(DataStore::c_DontWriteOut);
663  muonsList.registerInDataStore(DataStore::c_DontWriteOut);
664  amuonsList.registerInDataStore(DataStore::c_DontWriteOut);
666 
667  // initialise the lists
668  kLongList.create();
669  kLongList->initialize(Const::Klong.getPDGCode(), kLongList.getName());
670  muonsList.create();
671  muonsList->initialize(Const::muon.getPDGCode(), muonsList.getName());
672  amuonsList.create();
673  amuonsList->initialize(-Const::muon.getPDGCode(), amuonsList.getName());
674  amuonsList->bindAntiParticleList(*(muonsList));
675 
676  // add some tracks
677  const Track* t1 = tracks.appendNew(Track());
678  const Track* t2 = tracks.appendNew(Track());
679  const Track* t3 = tracks.appendNew(Track());
680  tracks.appendNew(Track());
681  tracks.appendNew(Track());
682 
683  // mock up some TrackFits for them (all muons)
684  TRandom3 generator;
685  TMatrixDSym cov6(6);
686  auto CDCValue = static_cast<unsigned long long int>(0x300000000000000);
687 
688  for (int i = 0; i < tracks.getEntries(); ++i) {
689  int charge = (i % 2 == 0) ? +1 : -1;
690  ROOT::Math::Cartesian2D d(generator.Uniform(-1, 1), generator.Uniform(-1, 1));
691  ROOT::Math::Cartesian2D pt(generator.Uniform(-1, 1), generator.Uniform(-1, 1));
692  d.SetXY(d.X(), -(d.X()*pt.x()) / pt.y());
693  B2Vector3D position(d.X(), d.Y(), generator.Uniform(-1, 1));
694  B2Vector3D momentum(pt.x(), pt.y(), generator.Uniform(-1, 1));
695  trackFits.appendNew(position, momentum, cov6, charge, Const::muon, 0.5, 1.5, CDCValue, 16777215, 0);
696  tracks[i]->setTrackFitResultIndex(Const::muon, i);
697  }
698 
699  // add some clusters
700  KLMCluster* klm1 = klmClusters.appendNew(KLMCluster());
701  klm1->setTime(1.1);
702  klm1->setClusterPosition(1.1, 1.1, 1.0);
703  klm1->setLayers(1);
704  klm1->setInnermostLayer(1);
705  klm1->setMomentumMag(1.0);
706  KLMCluster* klm2 = klmClusters.appendNew(KLMCluster());
707  klm2->setTime(1.2);
708  klm2->setClusterPosition(1.2, 1.2, 2.0);
709  klm2->setLayers(2);
710  klm2->setInnermostLayer(2);
711  klm2->setMomentumMag(1.0);
712  KLMCluster* klm3 = klmClusters.appendNew(KLMCluster());
713  klm3->setTime(1.3);
714  klm3->setClusterPosition(1.3, 1.3, 3.0);
715  klm3->setLayers(3);
716  klm3->setInnermostLayer(3);
717  klm3->setMomentumMag(1.0);
718 
719  // and add clusters related to the tracks
720  // case 1: 1 track --> 1 cluster
721  KLMCluster* klm4 = klmClusters.appendNew(KLMCluster());
722  klm4->setTime(1.4);
723  klm4->setClusterPosition(-1.4, -1.4, 1.0);
724  klm4->setLayers(4);
725  klm4->setInnermostLayer(4);
726  klm4->setMomentumMag(1.0);
727  t1->addRelationTo(klm4);
728 
729  // case 2: 2 tracks --> 1 cluster
730  KLMCluster* klm5 = klmClusters.appendNew(KLMCluster());
731  klm5->setTime(1.5);
732  klm5->setClusterPosition(-1.5, -1.5, 1.0);
733  klm5->setLayers(5);
734  klm5->setInnermostLayer(5);
735  klm5->setMomentumMag(1.0);
736  t2->addRelationTo(klm5);
737  t3->addRelationTo(klm5);
738 
739  // case 3: 1 track --> 2 clusters
740  // not possible
741 
742  // make the KLong from clusters (and sum up the total KLM momentum magnitude)
743  double klmMomentum = 0.0;
744  for (int i = 0; i < klmClusters.getEntries(); ++i) {
745  klmMomentum += klmClusters[i]->getMomentumMag();
746  if (!klmClusters[i]->getAssociatedTrackFlag()) {
747  const Particle* p = particles.appendNew(Particle(klmClusters[i]));
748  kLongList->addParticle(p);
749  }
750  }
751 
752  // make the muons from tracks
753  for (int i = 0; i < tracks.getEntries(); ++i) {
754  const Particle* p = particles.appendNew(Particle(tracks[i], Const::muon));
755  muonsList->addParticle(p);
756  }
757 
758  // grab variables
759  const Manager::Var* vClusterP = Manager::Instance().getVariable("klmClusterMomentum");
760  const Manager::Var* vClNTrack = Manager::Instance().getVariable("nKLMClusterTrackMatches");
761 
762  // calculate the total KLM momentum from the KLong list --> VM
763  double totalKLongMomentum = 0.0;
764  for (size_t i = 0; i < kLongList->getListSize(); ++i)
765  totalKLongMomentum += std::get<double>(vClusterP->function(kLongList->getParticle(i)));
766 
767  // calculate the total KLM momentum from muon-matched list --> VM
768  double totalMuonMomentum = 0.0;
769  for (size_t i = 0; i < muonsList->getListSize(); ++i) { // includes antiparticles
770  double muonMomentum = std::get<double>(vClusterP->function(muonsList->getParticle(i)));
771  double nOtherCl = std::get<double>(vClNTrack->function(muonsList->getParticle(i)));
772  if (nOtherCl > 0)
773  totalMuonMomentum += muonMomentum / nOtherCl;
774  }
775 
776  EXPECT_FLOAT_EQ(5.0, klmMomentum);
777  EXPECT_FLOAT_EQ(totalKLongMomentum + totalMuonMomentum, klmMomentum);
778  }
779 
780  TEST_F(KLMVariableTest, TrackToKLMClusterMatchingTest)
781  {
782  StoreArray<Particle> particles;
783  StoreArray<Track> tracks;
784  StoreArray<TrackFitResult> trackFits;
785  StoreArray<KLMCluster> klmClusters;
786 
787  // add a TrackFitResult
788  B2Vector3D position(1.0, 0, 0);
789  B2Vector3D momentum(0, 1.0, 0);
790  TMatrixDSym cov6(6);
791  const int charge = 1;
792  const float pValue = 0.5;
793  const float bField = 1.5;
794  auto CDCValue = static_cast<unsigned long long int>(0x300000000000000);
795  trackFits.appendNew(position, momentum, cov6, charge, Const::muon, pValue, bField, CDCValue, 16777215, 0);
796 
797  // add one Track
798  Track myTrack;
800  Track* muonTrack = tracks.appendNew(myTrack);
801 
802  // add two KLMClusters
803  KLMCluster* klm1 = klmClusters.appendNew(KLMCluster());
804  klm1->setTime(1.1);
805  klm1->setClusterPosition(1.1, 1.1, 1.0);
806  klm1->setLayers(1);
807  klm1->setInnermostLayer(1);
808  klm1->setMomentumMag(1.0);
809  KLMCluster* klm2 = klmClusters.appendNew(KLMCluster());
810  klm2->setTime(1.2);
811  klm2->setClusterPosition(1.2, 1.2, 2.0);
812  klm2->setLayers(2);
813  klm2->setInnermostLayer(2);
814  klm2->setMomentumMag(1.0);
815 
816  // and add a weighted relationship between the track and both clusters
817  // only the relation with klm1 must be returned
818  // in reconstruction we set a relation to only one cluster (if any),
819  // so here we test that getKLMCluster() returns the first cluster
820  // stored in the RelationVector
821  float distance1 = 11.1;
822  muonTrack->addRelationTo(klm1, 1. / distance1);
823  float distance2 = 2.2;
824  muonTrack->addRelationTo(klm2, 1. / distance2);
825 
826  // add a Particle
827  const Particle* muon = particles.appendNew(Particle(muonTrack, Const::muon));
828 
829  // grab variables
830  const Manager::Var* vTrNClusters = Manager::Instance().getVariable("nMatchedKLMClusters");
831  const Manager::Var* vClusterInnermostLayer = Manager::Instance().getVariable("klmClusterInnermostLayer");
832  const Manager::Var* vClusterTrackDistance = Manager::Instance().getVariable("klmClusterTrackDistance");
833 
834  EXPECT_POSITIVE(std::get<double>(vTrNClusters->function(muon)));
835  EXPECT_FLOAT_EQ(1.0, std::get<double>(vClusterInnermostLayer->function(muon)));
836  EXPECT_FLOAT_EQ(distance1, std::get<double>(vClusterTrackDistance->function(muon)));
837 
838  // add a Pion - no clusters matched here
839  trackFits.appendNew(position, momentum, cov6, charge, Const::pion, pValue, bField, CDCValue, 16777215, 0);
840  Track mySecondTrack;
841  mySecondTrack.setTrackFitResultIndex(Const::pion, 0);
842  Track* pionTrack = tracks.appendNew(mySecondTrack);
843  const Particle* pion = particles.appendNew(Particle(pionTrack, Const::pion));
844 
845  EXPECT_FLOAT_EQ(0.0, std::get<double>(vTrNClusters->function(pion)));
846  }
847 }
static const ChargedStable muon
muon particle
Definition: Const.h:541
static const ChargedStable pion
charged pion particle
Definition: Const.h:542
static const ParticleType Klong
K^0_L particle.
Definition: Const.h:558
@ c_DontWriteOut
Object/array should be NOT saved by output modules.
Definition: DataStore.h:71
static DataStore & Instance()
Instance of singleton Store.
Definition: DataStore.cc:52
void setInitializeActive(bool active)
Setter for m_initializeActive.
Definition: DataStore.cc:92
void reset(EDurability durability)
Frees memory occupied by data store items and removes all objects from the map.
Definition: DataStore.cc:84
ECL cluster data.
Definition: ECLCluster.h:27
void addHypothesis(EHypothesisBit bitmask)
Add bitmask to current hypothesis.
Definition: ECLCluster.h:129
void setTheta(double theta)
Set Theta of Shower (radian).
Definition: ECLCluster.h:217
void setConnectedRegionId(int crid)
Set connected region id.
Definition: ECLCluster.h:142
void setPhi(double phi)
Set Phi of Shower (radian).
Definition: ECLCluster.h:220
void setTime(double time)
Set time information.
Definition: ECLCluster.h:211
void setDeltaTime99(double dtime99)
Set 99% time containment range.
Definition: ECLCluster.h:214
void setClusterId(int clusterid)
Set cluster id.
Definition: ECLCluster.h:145
void setHypothesis(EHypothesisBit hypothesis)
Set hypotheses.
Definition: ECLCluster.h:123
void setEnergy(double energy)
Set Corrected Energy (GeV).
Definition: ECLCluster.h:226
void setIsTrack(bool istrack)
Set m_isTrack true if the cluster matches with a track.
Definition: ECLCluster.h:104
@ c_nPhotons
CR is split into n photons (N1)
@ c_neutralHadron
CR is reconstructed as a neutral hadron (N2)
void setR(double r)
Set R (in cm).
Definition: ECLCluster.h:223
Singleton class responsible for loading detector parameters from an XML file.
Definition: Gearbox.h:34
KLM cluster data.
Definition: KLMCluster.h:28
void setLayers(int layers)
Set number of layers with hits.
Definition: KLMCluster.h:145
void setClusterPosition(float globalX, float globalY, float globalZ)
Set global position.
Definition: KLMCluster.h:161
void setTime(float time)
Set time.
Definition: KLMCluster.h:138
void setInnermostLayer(int innermostLayer)
Set number of the innermost layer with hits.
Definition: KLMCluster.h:152
void setMomentumMag(float momentumMag)
Set momentum magnitude.
Definition: KLMCluster.h:172
ParticleGenerator is a generator for all the particles combined from the given ParticleLists.
Class to store reconstructed particles.
Definition: Particle.h:74
EFlavorType getFlavorType() const
Returns flavor type of the decay (for FS particles: flavor type of particle)
Definition: Particle.h:418
int getPDGCode(void) const
Returns PDG code.
Definition: Particle.h:403
@ c_Unflavored
Is its own antiparticle or we don't know whether it is a particle/antiparticle.
Definition: Particle.h:93
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).
int getArrayIndex() const
Returns this object's array index (in StoreArray), or -1 if not found.
bool registerInDataStore(DataStore::EStoreFlags storeFlags=DataStore::c_WriteOut)
Register the object/array in the DataStore.
Accessor to arrays stored in the data store.
Definition: StoreArray.h:113
T * appendNew()
Construct a new T object at the end of the array.
Definition: StoreArray.h:246
int getEntries() const
Get the number of objects in the array.
Definition: StoreArray.h:216
bool registerRelationTo(const StoreArray< TO > &toArray, DataStore::EDurability durability=DataStore::c_Event, DataStore::EStoreFlags storeFlags=DataStore::c_WriteOut, const std::string &namedRelation="") const
Register a relation to the given StoreArray.
Definition: StoreArray.h:140
Type-safe access to single objects in the data store.
Definition: StoreObjPtr.h:95
Class that bundles various TrackFitResults.
Definition: Track.h:25
void setTrackFitResultIndex(const Const::ChargedStable &chargedStable, short index)
Set an index (for positive values) or unavailability-code (with negative values) for a specific mass ...
Definition: Track.h:94
const Var * getVariable(std::string name)
Get the variable belonging to the given key.
Definition: Manager.cc:57
static Manager & Instance()
get singleton instance.
Definition: Manager.cc:25
TEST_F(ChargedParticleIdentificatorTest, TestDBRep)
Test correct storage of weightfiles in the database representation inner structure.
static Gearbox & getInstance()
Return reference to the Gearbox instance.
Definition: Gearbox.cc:81
double charge(int pdgCode)
Returns electric charge of a particle with given pdg code.
Definition: EvtPDLUtil.cc:44
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:23
A variable returning a floating-point value for a given Particle.
Definition: Manager.h:146
FunctionPtr function
Pointer to function.
Definition: Manager.h:147