Belle II Software light-2406-ragdoll
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
35using namespace Belle2;
36using namespace Belle2::Variable;
37
38namespace {
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 ROOT::Math::XYZVector position(d.X(), d.Y(), generator.Uniform(-1, 1));
87 ROOT::Math::XYZVector 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);
95 e1->setHypothesis(ECLCluster::EHypothesisBit::c_nPhotons);
96 e1->setClusterId(1);
97 e1->setTime(1);
98 e1->setDeltaTime99(0.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);
106 e2->setHypothesis(ECLCluster::EHypothesisBit::c_nPhotons);
107 e2->setClusterId(2);
108 e2->setTime(2);
109 e2->setDeltaTime99(0.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);
116 e3->setHypothesis(ECLCluster::EHypothesisBit::c_nPhotons);
117 e3->addHypothesis(ECLCluster::EHypothesisBit::c_neutralHadron);
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);
129 e4->setHypothesis(ECLCluster::EHypothesisBit::c_nPhotons);
130 e4->setClusterId(4);
131 t1->addRelationTo(e4);
132 e4->setIsTrack(true);
133
134 ECLCluster* e5 = eclclusters.appendNew(ECLCluster());
135 e5->setEnergy(0.3);
136 e5->setHypothesis(ECLCluster::EHypothesisBit::c_nPhotons);
137 e5->setClusterId(5);
138 t2->addRelationTo(e5);
139 e5->setIsTrack(true);
140
141 ECLCluster* e6 = eclclusters.appendNew(ECLCluster());
142 e6->setEnergy(0.2);
143 e6->setHypothesis(ECLCluster::EHypothesisBit::c_nPhotons);
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.0276554);
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.6035342);
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.7839828);
205 EXPECT_FLOAT_EQ(std::get<double>(b2bClusterPhi->function(gammalist->getParticle(2))), -1.3155545);
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.0442522);
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.2625608);
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());
470 eclst0->setHypothesis(ECLCluster::EHypothesisBit::c_nPhotons);
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());
477 eclst1->setHypothesis(ECLCluster::EHypothesisBit::c_nPhotons);
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 gammalist->setEditable(true);
574
575 // make the photons from clusters
576 for (int i = 0; i < eclclusters.getEntries(); ++i) {
577 if (!eclclusters[i]->isTrack()) {
578 const Particle* p = particles.appendNew(Particle(eclclusters[i]));
579 gammalist->addParticle(p);
580 }
581 }
582 gammalist->setEditable(false);
583
584 // initialize the pion lists
585 pionlist.create();
586 pionlist->initialize(211, pionlist.getName());
587 piminuslist.create();
588 piminuslist->initialize(-211, piminuslist.getName());
589 piminuslist->bindAntiParticleList(*(pionlist));
590 pionlist->setEditable(true);
591
592 // fill the pion list with the tracks
593 for (int i = 0; i < tracks.getEntries(); ++i) {
594 const Particle* p = particles.appendNew(Particle(tracks[i], Const::pion));
595 pionlist->addParticle(p);
596 }
597 pionlist->setEditable(false);
598
599 // check overlap without any arguments
600 const Manager::Var* photonHasOverlapNoArgs = Manager::Instance().getVariable("photonHasOverlap()");
601 // the track list e-:all is missing so NaN is returned
602 EXPECT_TRUE(std::isnan(std::get<double>(photonHasOverlapNoArgs->function(particles[0]))));
603
604 // check overlap without any requirement on other photons
605 const Manager::Var* photonHasOverlapAll = Manager::Instance().getVariable("photonHasOverlap(, gamma:all, pi+:all)");
606 // cluster 3 and cluster 1 share the connected region
607 EXPECT_TRUE(std::get<double>(photonHasOverlapAll->function(particles[0])));
608 // photonHasOverlap is designed for photons, so calling it for a pion returns NaN
609 EXPECT_TRUE(std::isnan(std::get<double>(photonHasOverlapAll->function(particles[3]))));
610
611 // check overlap with photons in barrel
612 const Manager::Var* photonHasOverlapBarrel = Manager::Instance().getVariable("photonHasOverlap(clusterReg==2, gamma:all, pi+:all)");
613 // cluster 3 is in the forward end cap so it doesn't matter that it has the same connected region like cluster 1
614 EXPECT_FALSE(std::get<double>(photonHasOverlapBarrel->function(particles[0])));
615 }
616
617 class KLMVariableTest : public ::testing::Test {
618 protected:
620 void SetUp() override
621 {
622 // setup the DataStore
624
625 // particles (to be filled)
626 StoreArray<Particle> particles;
627 particles.registerInDataStore();
628
629 // mock up mdst objects
630 StoreArray<Track> tracks;
631 tracks.registerInDataStore();
633 trackFits.registerInDataStore();
634 StoreArray<KLMCluster> klmClusters;
635 klmClusters.registerInDataStore();
636
637 // tracks can be matched to clusters
638 tracks.registerRelationTo(klmClusters);
639
640 // we're done setting up the datastore
642 }
643
645 void TearDown() override
646 {
648 }
649 };
650
651 TEST_F(KLMVariableTest, WholeEventClosure)
652 {
653 // we need the Particles, Tracks, TrackFitResults and KLMClusters StoreArrays
654 StoreArray<Particle> particles;
655 StoreArray<Track> tracks;
657 StoreArray<KLMCluster> klmClusters;
658
659 // create a KLong (clusters) and muon (tracks) lists
660 StoreObjPtr<ParticleList> kLongList("K0_L:testKLong");
661 StoreObjPtr<ParticleList> muonsList("mu-:testMuons");
662 StoreObjPtr<ParticleList> amuonsList("mu+:testMuons");
663
664 // register the lists in the datastore
666 kLongList.registerInDataStore(DataStore::c_DontWriteOut);
667 muonsList.registerInDataStore(DataStore::c_DontWriteOut);
668 amuonsList.registerInDataStore(DataStore::c_DontWriteOut);
670
671 // initialise the lists
672 kLongList.create();
673 kLongList->initialize(Const::Klong.getPDGCode(), kLongList.getName());
674 muonsList.create();
675 muonsList->initialize(Const::muon.getPDGCode(), muonsList.getName());
676 amuonsList.create();
677 amuonsList->initialize(-Const::muon.getPDGCode(), amuonsList.getName());
678 amuonsList->bindAntiParticleList(*(muonsList));
679
680 // add some tracks
681 const Track* t1 = tracks.appendNew(Track());
682 const Track* t2 = tracks.appendNew(Track());
683 const Track* t3 = tracks.appendNew(Track());
684 tracks.appendNew(Track());
685 tracks.appendNew(Track());
686
687 // mock up some TrackFits for them (all muons)
688 TRandom3 generator;
689 TMatrixDSym cov6(6);
690 auto CDCValue = static_cast<unsigned long long int>(0x300000000000000);
691
692 for (int i = 0; i < tracks.getEntries(); ++i) {
693 int charge = (i % 2 == 0) ? +1 : -1;
694 ROOT::Math::Cartesian2D d(generator.Uniform(-1, 1), generator.Uniform(-1, 1));
695 ROOT::Math::Cartesian2D pt(generator.Uniform(-1, 1), generator.Uniform(-1, 1));
696 d.SetXY(d.X(), -(d.X()*pt.X()) / pt.Y());
697 ROOT::Math::XYZVector position(d.X(), d.Y(), generator.Uniform(-1, 1));
698 ROOT::Math::XYZVector momentum(pt.X(), pt.Y(), generator.Uniform(-1, 1));
699 trackFits.appendNew(position, momentum, cov6, charge, Const::muon, 0.5, 1.5, CDCValue, 16777215, 0);
700 tracks[i]->setTrackFitResultIndex(Const::muon, i);
701 }
702
703 // add some clusters
704 KLMCluster* klm1 = klmClusters.appendNew(KLMCluster());
705 klm1->setTime(1.1);
706 klm1->setClusterPosition(1.1, 1.1, 1.0);
707 klm1->setLayers(1);
708 klm1->setInnermostLayer(1);
709 klm1->setMomentumMag(1.0);
710 KLMCluster* klm2 = klmClusters.appendNew(KLMCluster());
711 klm2->setTime(1.2);
712 klm2->setClusterPosition(1.2, 1.2, 2.0);
713 klm2->setLayers(2);
714 klm2->setInnermostLayer(2);
715 klm2->setMomentumMag(1.0);
716 KLMCluster* klm3 = klmClusters.appendNew(KLMCluster());
717 klm3->setTime(1.3);
718 klm3->setClusterPosition(1.3, 1.3, 3.0);
719 klm3->setLayers(3);
720 klm3->setInnermostLayer(3);
721 klm3->setMomentumMag(1.0);
722
723 // and add clusters related to the tracks
724 // case 1: 1 track --> 1 cluster
725 KLMCluster* klm4 = klmClusters.appendNew(KLMCluster());
726 klm4->setTime(1.4);
727 klm4->setClusterPosition(-1.4, -1.4, 1.0);
728 klm4->setLayers(4);
729 klm4->setInnermostLayer(4);
730 klm4->setMomentumMag(1.0);
731 t1->addRelationTo(klm4);
732
733 // case 2: 2 tracks --> 1 cluster
734 KLMCluster* klm5 = klmClusters.appendNew(KLMCluster());
735 klm5->setTime(1.5);
736 klm5->setClusterPosition(-1.5, -1.5, 1.0);
737 klm5->setLayers(5);
738 klm5->setInnermostLayer(5);
739 klm5->setMomentumMag(1.0);
740 t2->addRelationTo(klm5);
741 t3->addRelationTo(klm5);
742
743 // case 3: 1 track --> 2 clusters
744 // not possible
745
746 // make the KLong from clusters (and sum up the total KLM momentum magnitude)
747 double klmMomentum = 0.0;
748 for (int i = 0; i < klmClusters.getEntries(); ++i) {
749 klmMomentum += klmClusters[i]->getMomentumMag();
750 if (!klmClusters[i]->getAssociatedTrackFlag()) {
751 const Particle* p = particles.appendNew(Particle(klmClusters[i]));
752 kLongList->addParticle(p);
753 }
754 }
755
756 // make the muons from tracks
757 for (int i = 0; i < tracks.getEntries(); ++i) {
758 const Particle* p = particles.appendNew(Particle(tracks[i], Const::muon));
759 muonsList->addParticle(p);
760 }
761
762 // grab variables
763 const Manager::Var* vClusterP = Manager::Instance().getVariable("klmClusterMomentum");
764 const Manager::Var* vClNTrack = Manager::Instance().getVariable("nKLMClusterTrackMatches");
765
766 // calculate the total KLM momentum from the KLong list --> VM
767 double totalKLongMomentum = 0.0;
768 for (size_t i = 0; i < kLongList->getListSize(); ++i)
769 totalKLongMomentum += std::get<double>(vClusterP->function(kLongList->getParticle(i)));
770
771 // calculate the total KLM momentum from muon-matched list --> VM
772 double totalMuonMomentum = 0.0;
773 for (size_t i = 0; i < muonsList->getListSize(); ++i) { // includes antiparticles
774 double muonMomentum = std::get<double>(vClusterP->function(muonsList->getParticle(i)));
775 double nOtherCl = std::get<double>(vClNTrack->function(muonsList->getParticle(i)));
776 if (nOtherCl > 0)
777 totalMuonMomentum += muonMomentum / nOtherCl;
778 }
779
780 EXPECT_FLOAT_EQ(5.0, klmMomentum);
781 EXPECT_FLOAT_EQ(totalKLongMomentum + totalMuonMomentum, klmMomentum);
782 }
783
784 TEST_F(KLMVariableTest, TrackToKLMClusterMatchingTest)
785 {
786 StoreArray<Particle> particles;
787 StoreArray<Track> tracks;
789 StoreArray<KLMCluster> klmClusters;
790
791 // add a TrackFitResult
792 ROOT::Math::XYZVector position(1.0, 0, 0);
793 ROOT::Math::XYZVector momentum(0, 1.0, 0);
794 TMatrixDSym cov6(6);
795 const int charge = 1;
796 const float pValue = 0.5;
797 const float bField = 1.5;
798 auto CDCValue = static_cast<unsigned long long int>(0x300000000000000);
799 trackFits.appendNew(position, momentum, cov6, charge, Const::muon, pValue, bField, CDCValue, 16777215, 0);
800
801 // add one Track
802 Track myTrack;
804 Track* muonTrack = tracks.appendNew(myTrack);
805
806 // add two KLMClusters
807 KLMCluster* klm1 = klmClusters.appendNew(KLMCluster());
808 klm1->setTime(1.1);
809 klm1->setClusterPosition(1.1, 1.1, 1.0);
810 klm1->setLayers(1);
811 klm1->setInnermostLayer(1);
812 klm1->setMomentumMag(1.0);
813 KLMCluster* klm2 = klmClusters.appendNew(KLMCluster());
814 klm2->setTime(1.2);
815 klm2->setClusterPosition(1.2, 1.2, 2.0);
816 klm2->setLayers(2);
817 klm2->setInnermostLayer(2);
818 klm2->setMomentumMag(1.0);
819
820 // and add a weighted relationship between the track and both clusters
821 // only the relation with klm1 must be returned
822 // in reconstruction we set a relation to only one cluster (if any),
823 // so here we test that getKLMCluster() returns the first cluster
824 // stored in the RelationVector
825 float distance1 = 11.1;
826 muonTrack->addRelationTo(klm1, 1. / distance1);
827 float distance2 = 2.2;
828 muonTrack->addRelationTo(klm2, 1. / distance2);
829
830 // add a Particle
831 const Particle* muon = particles.appendNew(Particle(muonTrack, Const::muon));
832
833 // grab variables
834 const Manager::Var* vTrNClusters = Manager::Instance().getVariable("nMatchedKLMClusters");
835 const Manager::Var* vClusterInnermostLayer = Manager::Instance().getVariable("klmClusterInnermostLayer");
836 const Manager::Var* vClusterTrackDistance = Manager::Instance().getVariable("klmClusterTrackDistance");
837
838 EXPECT_POSITIVE(std::get<double>(vTrNClusters->function(muon)));
839 EXPECT_FLOAT_EQ(1.0, std::get<double>(vClusterInnermostLayer->function(muon)));
840 EXPECT_FLOAT_EQ(distance1, std::get<double>(vClusterTrackDistance->function(muon)));
841
842 // add a Pion - no clusters matched here
843 trackFits.appendNew(position, momentum, cov6, charge, Const::pion, pValue, bField, CDCValue, 16777215, 0);
844 Track mySecondTrack;
845 mySecondTrack.setTrackFitResultIndex(Const::pion, 0);
846 Track* pionTrack = tracks.appendNew(mySecondTrack);
847 const Particle* pion = particles.appendNew(Particle(pionTrack, Const::pion));
848
849 EXPECT_FLOAT_EQ(0.0, std::get<double>(vTrNClusters->function(pion)));
850 }
851}
static const ChargedStable muon
muon particle
Definition: Const.h:660
static const ChargedStable pion
charged pion particle
Definition: Const.h:661
static const ParticleType Klong
K^0_L particle.
Definition: Const.h:678
@ 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:54
void setInitializeActive(bool active)
Setter for m_initializeActive.
Definition: DataStore.cc:94
void reset(EDurability durability)
Frees memory occupied by data store items and removes all objects from the map.
Definition: DataStore.cc:86
ECL cluster data.
Definition: ECLCluster.h:27
void 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
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:75
EFlavorType getFlavorType() const
Returns flavor type of the decay (for FS particles: flavor type of particle)
Definition: Particle.h:469
int getPDGCode(void) const
Returns PDG code.
Definition: Particle.h:454
@ c_Unflavored
Is its own antiparticle or we don't know whether it is a particle/antiparticle.
Definition: Particle.h:95
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:96
Class that bundles various TrackFitResults.
Definition: Track.h:25
void setTrackFitResultIndex(const Const::ChargedStable &chargedStable, short index)
Set an index (for positive values) or unavailability-code (index = -1) for a specific mass hypothesis...
Definition: Track.h:174
const Var * getVariable(std::string name)
Get the variable belonging to the given key.
Definition: Manager.cc:57
static Manager & Instance()
get singleton instance.
Definition: Manager.cc:25
static Gearbox & getInstance()
Return reference to the Gearbox instance.
Definition: Gearbox.cc:81
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:24
A variable returning a floating-point value for a given Particle.
Definition: Manager.h:146
FunctionPtr function
Pointer to function.
Definition: Manager.h:147