Belle II Software development
ROEVariables.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// Own header.
10#include <analysis/variables/ROEVariables.h>
11
12#include <analysis/variables/Variables.h>
13
14// framework - DataStore
15#include <framework/datastore/StoreArray.h>
16#include <framework/datastore/StoreObjPtr.h>
17
18// dataobjects
19#include <analysis/dataobjects/Particle.h>
20#include <analysis/dataobjects/ParticleList.h>
21#include <mdst/dataobjects/MCParticle.h>
22#include <mdst/dataobjects/ECLCluster.h>
23
24// framework aux
25#include <framework/logging/Logger.h>
26#include <framework/utilities/Conversion.h>
27#include <framework/gearbox/Const.h>
28
29// utility
30#include <analysis/utility/PCmsLabTransform.h>
31#include <analysis/utility/ReferenceFrame.h>
32
33#include <TRandom.h>
34#include <TMath.h>
35#include <Math/AxisAngle.h>
36
37#include <iostream>
38
39using namespace std;
40
41namespace Belle2 {
46 namespace Variable {
47
48 double isInRestOfEvent(const Particle* particle)
49 {
50
51 StoreObjPtr<RestOfEvent> roeobjptr;
52 if (not roeobjptr.isValid())
53 return 0;
54
55 const RestOfEvent* roe = &(*roeobjptr);
56
57 return isInThisRestOfEvent(particle, roe);
58 }
59
60 double isCloneOfSignalSide(const Particle* particle)
61 {
62
63 StoreObjPtr<RestOfEvent> roe;
64 if (not roe.isValid()) {
65 B2WARNING("Please use isCloneOfSignalSide variable in for_each ROE loop!");
66 return Const::doubleNaN;
67 }
68 auto* particleMC = particle->getMCParticle();
69 if (!particleMC) {
70 return 0.0;
71 }
72 auto* signal = roe->getRelatedFrom<Particle>();
73 auto signalFSPs = signal->getFinalStateDaughters();
74 for (auto* daughter : signalFSPs) {
75 auto* daughterMC = daughter->getMCParticle();
76 if (daughterMC == particleMC) {
77 return 1.0;
78 }
79 }
80 return 0.0;
81 }
82
83 double hasAncestorFromSignalSide(const Particle* particle)
84 {
85 StoreObjPtr<RestOfEvent> roe;
86 if (!roe.isValid()) {
87 B2WARNING("Please use hasAncestorFromSignalSide variable in for_each ROE loop!");
88 return Const::doubleNaN;
89 }
90 auto* particleMC = particle->getMCParticle();
91 if (!particleMC) {
92 return 0.0;
93 }
94 auto* signalReco = roe->getRelatedFrom<Particle>();
95 auto* signalMC = signalReco->getMCParticle();
96 MCParticle* ancestorMC = particleMC->getMother();
97 while (ancestorMC) {
98 if (ancestorMC == signalMC) {
99 return 1.0;
100 }
101 ancestorMC = ancestorMC->getMother();
102 }
103 return 0.0;
104 }
105
106 Manager::FunctionPtr currentROEIsInList(const std::vector<std::string>& arguments)
107 {
108 if (arguments.size() != 1)
109 B2FATAL("Wrong number of arguments (1 required) for meta function currentROEIsInList");
110
111 std::string listName = arguments[0];
112
113 auto func = [listName](const Particle*) -> double {
114
115 StoreObjPtr<ParticleList> particleList(listName);
116 if (!(particleList.isValid()))
117 {
118 B2FATAL("Invalid Listname " << listName << " given to currentROEIsInList!");
119 }
120 StoreObjPtr<RestOfEvent> roe("RestOfEvent");
121
122 if (not roe.isValid())
123 return 0;
124
125 auto* particle = roe->getRelatedFrom<Particle>();
126 if (particle == nullptr)
127 {
128 B2ERROR("Relation between particle and ROE doesn't exist! currentROEIsInList() variable has to be called from ROE loop");
129 return Const::doubleNaN;
130 }
131 return particleList->contains(particle) ? 1 : 0;
132
133 };
134 return func;
135 }
136
137 Manager::FunctionPtr particleRelatedToCurrentROE(const std::vector<std::string>& arguments)
138 {
139 if (arguments.size() != 1)
140 B2FATAL("Wrong number of arguments (1 required) for meta function particleRelatedToCurrentROE");
141
142 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
143 auto func = [var](const Particle*) -> double {
144
145 StoreObjPtr<RestOfEvent> roe("RestOfEvent");
146
147 if (not roe.isValid())
148 return Const::doubleNaN;
149
150 auto* particle = roe->getRelatedFrom<Particle>();
151 if (particle == nullptr)
152 {
153 B2ERROR("Relation between particle and ROE doesn't exist! particleRelatedToCurrentROE() variable has to be called from ROE loop");
154 return Const::doubleNaN;
155 }
156 if (std::holds_alternative<double>(var->function(particle)))
157 {
158 return std::get<double>(var->function(particle));
159 } else if (std::holds_alternative<int>(var->function(particle)))
160 {
161 return std::get<int>(var->function(particle));
162 } else if (std::holds_alternative<bool>(var->function(particle)))
163 {
164 return std::get<bool>(var->function(particle));
165 } else return Const::doubleNaN;
166
167 };
168 return func;
169 }
170
171 Manager::FunctionPtr useROERecoilFrame(const std::vector<std::string>& arguments)
172 {
173 if (arguments.size() == 1) {
174 const Variable::Manager::Var* var = Manager::Instance().getVariable(arguments[0]);
175 auto func = [var](const Particle * particle) -> double {
176 // Here we prioritize old variable behaviour first:
177 const RestOfEvent* roe = particle->getRelatedTo<RestOfEvent>();
178 // if related ROE not found, get the StoreArray pointer
179 if (roe == nullptr)
180 {
181 StoreObjPtr<RestOfEvent> roeObjPtr("RestOfEvent");
182 if (roeObjPtr.isValid()) {
183 roe = &*roeObjPtr;
184 }
185 }
186 if (roe == nullptr)
187 {
188 B2ERROR("Neither relation between particle and ROE doesn't exist nor ROE object has not been found!");
189 return Const::doubleNaN;
190 }
191 PCmsLabTransform T;
192 ROOT::Math::PxPyPzEVector pRecoil = T.getBeamFourMomentum() - roe->get4Vector();
193 Particle tmp(pRecoil, 0);
194 UseReferenceFrame<RestFrame> frame(&tmp);
195 if (std::holds_alternative<double>(var->function(particle)))
196 {
197 return std::get<double>(var->function(particle));
198 } else if (std::holds_alternative<int>(var->function(particle)))
199 {
200 return std::get<int>(var->function(particle));
201 } else if (std::holds_alternative<bool>(var->function(particle)))
202 {
203 return std::get<bool>(var->function(particle));
204 } else return Const::doubleNaN;
205 };
206 return func;
207 } else {
208 B2WARNING("Wrong number of arguments for meta function useROERecoilFrame");
209 return nullptr;
210 }
211 }
212
213 // only the helper function
214 int nRemainingTracksInROE(const Particle* particle, const std::string& maskName)
215 {
216 StoreObjPtr<RestOfEvent> roe("RestOfEvent");
217 if (not roe.isValid())
218 return 0;
219 int n_roe_tracks = roe->getNTracks(maskName);
220 int n_par_tracks = 0;
221 const auto& daughters = particle->getFinalStateDaughters();
222 for (const auto& daughter : daughters) {
223 if (daughter->getParticleSource() == Particle::EParticleSourceObject::c_Track && roe->hasParticle(daughter, maskName)) {
224 n_par_tracks++;
225 }
226 }
227 return n_roe_tracks - n_par_tracks;
228 }
229
230 Manager::FunctionPtr nROE_RemainingTracksWithMask(const std::vector<std::string>& arguments)
231 {
232 std::string maskName = RestOfEvent::c_defaultMaskName;
233
234 if (arguments.size() == 1)
235 maskName = arguments[0];
236 else if (arguments.size() > 1)
237 B2FATAL("At most 1 argument (name of mask) accepted for meta function nROETracks");
238
239 auto func = [maskName](const Particle * particle) -> int {
240 return nRemainingTracksInROE(particle, maskName);
241 };
242 return func;
243 }
244
245 int nROE_RemainingTracks(const Particle* particle)
246 {
247 return nRemainingTracksInROE(particle);
248 }
249
250 double nROE_KLMClusters(const Particle* particle)
251 {
252 // Get related ROE object
253 const RestOfEvent* roe = getRelatedROEObject(particle);
254
255 if (!roe) {
256 B2ERROR("Relation between particle and ROE doesn't exist!");
257 return Const::doubleNaN;
258 }
259
260 return roe->getNKLMClusters();
261 }
262
263 double ROE_MC_E(const Particle* particle)
264 {
265 const MCParticle* mcp = particle->getMCParticle();
266
267 if (!mcp)
268 return Const::doubleNaN;
269
270 PCmsLabTransform T;
271 ROOT::Math::PxPyPzEVector boostvec = T.getBeamFourMomentum();
272 auto mcroe4vector = boostvec - mcp->get4Vector();
273 const auto& frame = ReferenceFrame::GetCurrent();
274 auto frameMCRoe4Vector = frame.getMomentum(mcroe4vector);
275 return frameMCRoe4Vector.energy();
276 }
277
278 double ROE_MC_P(const Particle* particle)
279 {
280 const MCParticle* mcp = particle->getMCParticle();
281
282 if (!mcp)
283 return Const::doubleNaN;
284
285 PCmsLabTransform T;
286 ROOT::Math::PxPyPzEVector boostvec = T.getBeamFourMomentum();
287 auto mcroe4vector = boostvec - mcp->get4Vector();
288 const auto& frame = ReferenceFrame::GetCurrent();
289 auto frameMCRoe4Vector = frame.getMomentum(mcroe4vector);
290 return frameMCRoe4Vector.P();
291 }
292
293 double ROE_MC_Px(const Particle* particle)
294 {
295 const MCParticle* mcp = particle->getMCParticle();
296
297 if (!mcp)
298 return Const::doubleNaN;
299
300 PCmsLabTransform T;
301 ROOT::Math::PxPyPzEVector boostvec = T.getBeamFourMomentum();
302 auto mcroe4vector = boostvec - mcp->get4Vector();
303 const auto& frame = ReferenceFrame::GetCurrent();
304 auto frameMCRoe4Vector = frame.getMomentum(mcroe4vector);
305
306 return frameMCRoe4Vector.px();
307 }
308
309 double ROE_MC_Py(const Particle* particle)
310 {
311 const MCParticle* mcp = particle->getMCParticle();
312
313 if (!mcp)
314 return Const::doubleNaN;
315
316 PCmsLabTransform T;
317 ROOT::Math::PxPyPzEVector boostvec = T.getBeamFourMomentum();
318 auto mcroe4vector = boostvec - mcp->get4Vector();
319 const auto& frame = ReferenceFrame::GetCurrent();
320 auto frameMCRoe4Vector = frame.getMomentum(mcroe4vector);
321
322 return frameMCRoe4Vector.py();
323 }
324
325 double ROE_MC_Pz(const Particle* particle)
326 {
327 const MCParticle* mcp = particle->getMCParticle();
328
329 if (!mcp)
330 return Const::doubleNaN;
331
332 PCmsLabTransform T;
333 ROOT::Math::PxPyPzEVector boostvec = T.getBeamFourMomentum();
334 auto mcroe4vector = boostvec - mcp->get4Vector();
335 const auto& frame = ReferenceFrame::GetCurrent();
336 auto frameMCRoe4Vector = frame.getMomentum(mcroe4vector);
337
338 return frameMCRoe4Vector.pz();
339 }
340
341 double ROE_MC_Pt(const Particle* particle)
342 {
343 const MCParticle* mcp = particle->getMCParticle();
344
345 if (!mcp)
346 return Const::doubleNaN;
347
348 PCmsLabTransform T;
349 ROOT::Math::PxPyPzEVector boostvec = T.getBeamFourMomentum();
350 auto mcroe4vector = boostvec - mcp->get4Vector();
351 const auto& frame = ReferenceFrame::GetCurrent();
352 auto frameMCRoe4Vector = frame.getMomentum(mcroe4vector);
353
354 return frameMCRoe4Vector.pt();
355 }
356
357 double ROE_MC_PTheta(const Particle* particle)
358 {
359 const MCParticle* mcp = particle->getMCParticle();
360
361 if (!mcp)
362 return Const::doubleNaN;
363
364 PCmsLabTransform T;
365 ROOT::Math::PxPyPzEVector boostvec = T.getBeamFourMomentum();
366 auto mcroe4vector = boostvec - mcp->get4Vector();
367 const auto& frame = ReferenceFrame::GetCurrent();
368 auto frameMCRoe4Vector = frame.getMomentum(mcroe4vector);
369
370 return frameMCRoe4Vector.Theta();
371 }
372
373 double ROE_MC_M(const Particle* particle)
374 {
375 const MCParticle* mcp = particle->getMCParticle();
376
377 if (!mcp)
378 return Const::doubleNaN;
379
380 PCmsLabTransform T;
381 ROOT::Math::PxPyPzEVector boostvec = T.getBeamFourMomentum();
382
383 return (boostvec - mcp->get4Vector()).M();
384 }
385
386 Manager::FunctionPtr ROE_MC_MissingFlags(const std::vector<std::string>& arguments)
387 {
388 std::string maskName = RestOfEvent::c_defaultMaskName;
389
390 if (arguments.size() == 1)
391 maskName = arguments[0];
392 else if (arguments.size() > 1)
393 B2FATAL("At most 1 argument (name of mask) accepted for meta function nROETracks");
394
395 auto func = [maskName](const Particle * particle) -> double {
396
397 StoreArray<Particle> particles;
398
399 //Get MC Particle of the B meson
400 const MCParticle* mcParticle = particle->getMCParticle();
401
402 if (!mcParticle)
403 return Const::doubleNaN;
404
405 // Get Mother
406 const MCParticle* mcMother = mcParticle->getMother();
407
408 if (!mcMother)
409 return Const::doubleNaN;
410
411 // Get daughters
412 std::vector<MCParticle*> mcDaughters = mcMother->getDaughters();
413
414 if (mcDaughters.size() != 2)
415 return Const::doubleNaN;
416
417 // Get the companion B meson
418 MCParticle* mcROE = nullptr;
419 if (mcDaughters[0]->getArrayIndex() == mcParticle->getArrayIndex())
420 mcROE = mcDaughters[1];
421 else
422 mcROE = mcDaughters[0];
423
424 // Get related ROE object
425 const RestOfEvent* roe = getRelatedROEObject(particle);
426
427 std::set<const MCParticle*> mcROEObjects;
428
429 auto roeParticles = roe->getParticles(maskName);
430 for (auto* roeParticle : roeParticles)
431 {
432 auto* mcroeParticle = roeParticle->getMCParticle();
433 if (mcroeParticle != nullptr) {
434 mcROEObjects.insert(mcroeParticle);
435 }
436 }
437 int flags = 0;
438 checkMCParticleMissingFlags(mcROE, mcROEObjects, flags);
439
440 return flags;
441 };
442 return func;
443 }
444
445 Manager::FunctionPtr nROE_Tracks(const std::vector<std::string>& arguments)
446 {
447 std::string maskName = RestOfEvent::c_defaultMaskName;
448
449 if (arguments.size() == 1)
450 maskName = arguments[0];
451 else if (arguments.size() > 1)
452 B2FATAL("At most 1 argument (name of mask) accepted for meta function nROETracks");
453
454 auto func = [maskName](const Particle * particle) -> double {
455
456 // Get related ROE object
457 const RestOfEvent* roe = getRelatedROEObject(particle);
458
459 if (!roe)
460 {
461 B2ERROR("Relation between particle and ROE doesn't exist!");
462 return Const::doubleNaN;
463 }
464
465 return roe->getNTracks(maskName);
466 };
467 return func;
468 }
469
470 Manager::FunctionPtr nROE_ECLClusters(const std::vector<std::string>& arguments)
471 {
472 std::string maskName = RestOfEvent::c_defaultMaskName;
473
474 if (arguments.size() == 1)
475 maskName = arguments[0];
476 else if (arguments.size() > 1)
477 B2FATAL("At most 1 argument (name of mask) accepted for meta function nROEECLClusters");
478
479 auto func = [maskName](const Particle * particle) -> double {
480
481 // Get related ROE object
482 const RestOfEvent* roe = getRelatedROEObject(particle);
483
484 if (!roe)
485 {
486 B2ERROR("Relation between particle and ROE doesn't exist!");
487 return Const::doubleNaN;
488 }
489
490 return roe->getNECLClusters(maskName);
491 };
492 return func;
493 }
494
495 Manager::FunctionPtr nROE_NeutralECLClusters(const std::vector<std::string>& arguments)
496 {
497 std::string maskName = RestOfEvent::c_defaultMaskName;
498
499 if (arguments.size() == 1)
500 maskName = arguments[0];
501 else if (arguments.size() > 1)
502 B2FATAL("At most 1 argument (name of mask) accepted for meta function nROENeutralECLClusters");
503
504 auto func = [maskName](const Particle * particle) -> double {
505
506 // Get related ROE object
507 const RestOfEvent* roe = getRelatedROEObject(particle);
508
509 if (!roe)
510 {
511 B2ERROR("Relation between particle and ROE doesn't exist!");
512 return Const::doubleNaN;
513 }
514
515 return roe->getPhotons(maskName).size();
516 };
517 return func;
518 }
519
520 Manager::FunctionPtr nROE_Photons(const std::vector<std::string>& arguments)
521 {
522 std::string maskName = RestOfEvent::c_defaultMaskName;
523
524 if (arguments.size() == 1) {
525 maskName = arguments[0];
526 } else if (arguments.size() > 1) {
527 B2FATAL("At most 1 argument (name of mask) accepted for meta function nROE_Photons");
528 }
529 auto func = [maskName](const Particle * particle) -> double {
530
531 // Get related ROE object
532 const RestOfEvent* roe = getRelatedROEObject(particle);
533
534 if (!roe)
535 {
536 B2ERROR("Relation between particle and ROE doesn't exist!");
537 return Const::doubleNaN;
538 }
539
540 // Get unused ECLClusters in ROE
541 auto roeClusters = roe->getPhotons(maskName);
542 int nPhotons = 0;
543
544 // Select ECLClusters with photon hypothesis
545 for (auto& roeCluster : roeClusters)
546 if (roeCluster->getECLClusterEHypothesisBit() == ECLCluster::EHypothesisBit::c_nPhotons)
547 nPhotons++;
548 return nPhotons;
549 };
550 return func;
551 }
552
553 Manager::FunctionPtr nROE_NeutralHadrons(const std::vector<std::string>& arguments)
554 {
555 std::string maskName = RestOfEvent::c_defaultMaskName;
556
557 if (arguments.size() == 1) {
558 maskName = arguments[0];
559 } else if (arguments.size() > 1) {
560 B2FATAL("At most 1 argument (name of mask) accepted for meta function nROE_NeutralHadrons");
561 }
562 auto func = [maskName](const Particle * particle) -> double {
563
564 // Get related ROE object
565 const RestOfEvent* roe = getRelatedROEObject(particle);
566
567 if (!roe)
568 {
569 B2ERROR("Relation between particle and ROE doesn't exist!");
570 return Const::doubleNaN;
571 }
572
573 return roe->getHadrons(maskName).size();
574 };
575 return func;
576 }
577
578 Manager::FunctionPtr nROE_ChargedParticles(const std::vector<std::string>& arguments)
579 {
580 std::string maskName = RestOfEvent::c_defaultMaskName;
581 int pdgCode = 0;
582 if (arguments.size() == 1) {
583 maskName = arguments[0];
584 } else if (arguments.size() == 2) {
585 maskName = arguments[0];
586 try {
587 pdgCode = Belle2::convertString<int>(arguments[1]);
588 } catch (std::invalid_argument&) {
589 B2ERROR("First argument of nROE_ChargedParticles must be a PDG code");
590 return nullptr;
591 }
592 } else if (arguments.size() > 2) {
593 B2FATAL("Wrong number of arguments (2 optional) for meta function nROE_ChargedParticles");
594 }
595 auto func = [maskName, pdgCode](const Particle * particle) -> double {
596
597 // Get related ROE object
598 const RestOfEvent* roe = getRelatedROEObject(particle);
599
600 if (!roe)
601 {
602 B2ERROR("Relation between particle and ROE doesn't exist!");
603 return Const::doubleNaN;
604 }
605
606 return roe->getChargedParticles(maskName, abs(pdgCode)).size();
607 };
608 return func;
609 }
610
611 Manager::FunctionPtr nROE_Composites(const std::vector<std::string>& arguments)
612 {
613 std::string maskName = RestOfEvent::c_defaultMaskName;
614
615 if (arguments.size() == 1) {
616 maskName = arguments[0];
617 } else if (arguments.size() > 1) {
618 B2FATAL("At most 1 argument (name of mask) accepted for meta function nROE_Composites");
619 }
620 auto func = [maskName](const Particle * particle) -> double {
621
622 // Get related ROE object
623 const RestOfEvent* roe = getRelatedROEObject(particle);
624
625 if (!roe)
626 {
627 B2ERROR("Relation between particle and ROE doesn't exist!");
628 return Const::doubleNaN;
629 }
630 int result = 0;
631 auto particles = roe->getParticles(maskName, false);
632
633 for (auto roeParticle : particles)
634 {
635 if (roeParticle->getParticleSource() == Particle::c_Composite or
636 roeParticle->getParticleSource() == Particle::c_V0) {
637 result++;
638 }
639 }
640 return result;
641 };
642 return func;
643 }
644
645 Manager::FunctionPtr nROE_ParticlesInList(const std::vector<std::string>& arguments)
646 {
647 std::string pListName;
648 std::string maskName = RestOfEvent::c_defaultMaskName;
649
650 if (arguments.size() == 1) {
651 pListName = arguments[0];
652 } else if (arguments.size() == 2) {
653 pListName = arguments[0];
654 maskName = arguments[1];
655 } else
656 B2FATAL("Wrong number of arguments (1 or 2 required) for meta function nROE_ParticlesInList");
657
658 auto func = [pListName, maskName](const Particle * particle) -> double {
659
660 // Get related ROE object
661 const RestOfEvent* roe = getRelatedROEObject(particle);
662
663 if (!roe)
664 {
665 B2ERROR("Relation between particle and ROE doesn't exist!");
666 return Const::doubleNaN;
667 }
668
669 int nPart = 0;
670
671 // Get particle list
672 StoreObjPtr<ParticleList> pList(pListName);
673 if (!pList.isValid())
674 B2FATAL("ParticleList " << pListName << " could not be found or is not valid!");
675
676 for (unsigned int i = 0; i < pList->getListSize(); i++)
677 {
678 const Particle* part = pList->getParticle(i);
679 if (isInThisRestOfEvent(part, roe, maskName))
680 ++nPart;
681 }
682
683 return nPart;
684 };
685 return func;
686 }
687
688 Manager::FunctionPtr ROE_Charge(const std::vector<std::string>& arguments)
689 {
690 std::string maskName = RestOfEvent::c_defaultMaskName;
691
692 if (arguments.size() == 1)
693 maskName = arguments[0];
694 else if (arguments.size() > 1)
695 B2FATAL("At most 1 argument (name of mask) accepted for meta function ROECharge");
696
697 auto func = [maskName](const Particle * particle) -> double {
698
699 // Get related ROE object
700 const RestOfEvent* roe = getRelatedROEObject(particle);
701
702 if (!roe)
703 {
704 B2ERROR("Relation between particle and ROE doesn't exist!");
705 return Const::doubleNaN;
706 }
707
708 // Get tracks in ROE
709 auto roeParticles = roe->getParticles(maskName);
710 int roeCharge = 0;
711
712 for (auto* roeParticle : roeParticles)
713 {
714 roeCharge += roeParticle->getCharge();
715 }
716
717 return roeCharge;
718 };
719 return func;
720 }
721
722 Manager::FunctionPtr ROE_ExtraEnergy(const std::vector<std::string>& arguments)
723 {
724 std::string maskName = RestOfEvent::c_defaultMaskName;
725
726 if (arguments.size() == 1)
727 maskName = arguments[0];
728 else if (arguments.size() > 1)
729 B2FATAL("At most 1 argument (name of mask) accepted for meta function extraEnergy");
730
731 auto func = [maskName](const Particle * particle) -> double {
732
733 // Get related ROE object
734 const RestOfEvent* roe = getRelatedROEObject(particle);
735
736 if (!roe)
737 {
738 B2ERROR("Relation between particle and ROE doesn't exist!");
739 return Const::doubleNaN;
740 }
741
742 double extraE = 0.0;
743
744 auto roeClusters = roe->getPhotons(maskName);
745
746 for (auto& roeCluster : roeClusters)
747 if (roeCluster->getECLClusterEHypothesisBit() == ECLCluster::EHypothesisBit::c_nPhotons)
748 extraE += roeCluster->getECLClusterEnergy();
749
750 auto roeChargedParticles = roe->getChargedParticles(maskName);
751
752 for (auto& roeChargedParticle : roeChargedParticles)
753 {
754 if (roeChargedParticle->getECLCluster())
755 extraE += roeChargedParticle->getECLClusterEnergy();
756 }
757
758 return extraE;
759 };
760 return func;
761 }
762
763 Manager::FunctionPtr ROE_NeutralExtraEnergy(const std::vector<std::string>& arguments)
764 {
765 std::string maskName = RestOfEvent::c_defaultMaskName;
766
767 if (arguments.size() == 1)
768 maskName = arguments[0];
769 else if (arguments.size() > 1)
770 B2FATAL("At most 1 argument (name of mask) accepted for meta function extraEnergy");
771
772 auto func = [maskName](const Particle * particle) -> double {
773
774 // Get related ROE object
775 const RestOfEvent* roe = getRelatedROEObject(particle);
776
777 if (!roe)
778 {
779 B2ERROR("Relation between particle and ROE doesn't exist!");
780 return Const::doubleNaN;
781 }
782 auto roephotons = roe->getPhotons(maskName);
783 ROOT::Math::PxPyPzEVector total4vector;
784 for (auto* photon : roephotons)
785 {
786 total4vector += photon->get4Vector();
787 }
788 const auto& frame = ReferenceFrame::GetCurrent();
789 auto frameRoe4Vector = frame.getMomentum(total4vector);
790 return frameRoe4Vector.energy();
791 };
792 return func;
793 }
794
795 Manager::FunctionPtr ROE_E(const std::vector<std::string>& arguments)
796 {
797 std::string maskName = RestOfEvent::c_defaultMaskName;
798
799 if (arguments.size() == 1)
800 maskName = arguments[0];
801 else if (arguments.size() > 1)
802 B2FATAL("At most 1 argument (name of mask) accepted for meta function ROE_E");
803
804 auto func = [maskName](const Particle * particle) -> double {
805 const RestOfEvent* roe = particle->getRelatedTo<RestOfEvent>();
806 if (!roe)
807 {
808 B2ERROR("Relation between particle and ROE doesn't exist!");
809 return Const::doubleNaN;
810 }
811 const auto& frame = ReferenceFrame::GetCurrent();
812 auto frameRoe4Vector = frame.getMomentum(roe->get4Vector(maskName));
813 return frameRoe4Vector.energy();
814 };
815 return func;
816 }
817
818 Manager::FunctionPtr ROE_M(const std::vector<std::string>& arguments)
819 {
820 std::string maskName = RestOfEvent::c_defaultMaskName;
821
822 if (arguments.size() == 1)
823 maskName = arguments[0];
824 else if (arguments.size() > 1)
825 B2FATAL("At most 1 argument (name of mask) accepted for meta function ROE_M");
826
827 auto func = [maskName](const Particle * particle) -> double {
828
829 // Get related ROE object
830 const RestOfEvent* roe = getRelatedROEObject(particle);
831
832 if (!roe)
833 {
834 B2ERROR("Relation between particle and ROE doesn't exist!");
835 return Const::doubleNaN;
836 }
837
838 return roe->get4Vector(maskName).M();
839 };
840 return func;
841 }
842
843 Manager::FunctionPtr ROE_P(const std::vector<std::string>& arguments)
844 {
845 std::string maskName = RestOfEvent::c_defaultMaskName;
846
847 if (arguments.size() == 1)
848 maskName = arguments[0];
849 else if (arguments.size() > 1)
850 B2FATAL("At most 1 argument (name of mask) accepted for meta function ROE_P");
851
852 auto func = [maskName](const Particle * particle) -> double {
853
854 // Get related ROE object
855 const RestOfEvent* roe = getRelatedROEObject(particle);
856
857 if (!roe)
858 {
859 B2ERROR("Relation between particle and ROE doesn't exist!");
860 return Const::doubleNaN;
861 }
862
863 const auto& frame = ReferenceFrame::GetCurrent();
864 auto frameRoe4Vector = frame.getMomentum(roe->get4Vector(maskName));
865 return frameRoe4Vector.P();
866 };
867 return func;
868 }
869
870 Manager::FunctionPtr ROE_Px(const std::vector<std::string>& arguments)
871 {
872 std::string maskName = RestOfEvent::c_defaultMaskName;
873
874 if (arguments.size() == 1)
875 maskName = arguments[0];
876 else if (arguments.size() > 1)
877 B2FATAL("At most 1 argument (name of mask) accepted for meta function ROE_Px");
878
879 auto func = [maskName](const Particle * particle) -> double {
880
881 // Get related ROE object
882 const RestOfEvent* roe = getRelatedROEObject(particle);
883
884 if (!roe)
885 {
886 B2ERROR("Relation between particle and ROE doesn't exist!");
887 return Const::doubleNaN;
888 }
889
890 const auto& frame = ReferenceFrame::GetCurrent();
891 auto frameRoe4Vector = frame.getMomentum(roe->get4Vector(maskName));
892 return frameRoe4Vector.px();
893 };
894 return func;
895 }
896
897 Manager::FunctionPtr ROE_Py(const std::vector<std::string>& arguments)
898 {
899 std::string maskName = RestOfEvent::c_defaultMaskName;
900
901 if (arguments.size() == 1)
902 maskName = arguments[0];
903 else if (arguments.size() > 1)
904 B2FATAL("At most 1 argument (name of mask) accepted for meta function ROE_Py");
905
906 auto func = [maskName](const Particle * particle) -> double {
907
908 // Get related ROE object
909 const RestOfEvent* roe = getRelatedROEObject(particle);
910
911 if (!roe)
912 {
913 B2ERROR("Relation between particle and ROE doesn't exist!");
914 return Const::doubleNaN;
915 }
916
917 const auto& frame = ReferenceFrame::GetCurrent();
918 auto frameRoe4Vector = frame.getMomentum(roe->get4Vector(maskName));
919 return frameRoe4Vector.py();
920 };
921 return func;
922 }
923
924 Manager::FunctionPtr ROE_Pt(const std::vector<std::string>& arguments)
925 {
926 std::string maskName = RestOfEvent::c_defaultMaskName;
927
928 if (arguments.size() == 1)
929 maskName = arguments[0];
930 else if (arguments.size() > 1)
931 B2FATAL("At most 1 argument (name of mask) accepted for meta function ROE_Pt");
932
933 auto func = [maskName](const Particle * particle) -> double {
934
935 // Get related ROE object
936 const RestOfEvent* roe = getRelatedROEObject(particle);
937
938 if (!roe)
939 {
940 B2ERROR("Relation between particle and ROE doesn't exist!");
941 return Const::doubleNaN;
942 }
943
944 const auto& frame = ReferenceFrame::GetCurrent();
945 auto frameRoe4Vector = frame.getMomentum(roe->get4Vector(maskName));
946 return frameRoe4Vector.pt();
947 };
948 return func;
949 }
950
951 Manager::FunctionPtr ROE_Pz(const std::vector<std::string>& arguments)
952 {
953 std::string maskName = RestOfEvent::c_defaultMaskName;
954
955 if (arguments.size() == 1)
956 maskName = arguments[0];
957 else if (arguments.size() > 1)
958 B2FATAL("At most 1 argument (name of mask) accepted for meta function ROE_Pz");
959
960 auto func = [maskName](const Particle * particle) -> double {
961
962 // Get related ROE object
963 const RestOfEvent* roe = getRelatedROEObject(particle);
964
965 if (!roe)
966 {
967 B2ERROR("Relation between particle and ROE doesn't exist!");
968 return Const::doubleNaN;
969 }
970
971 const auto& frame = ReferenceFrame::GetCurrent();
972 auto frameRoe4Vector = frame.getMomentum(roe->get4Vector(maskName));
973 return frameRoe4Vector.pz();
974 };
975 return func;
976 }
977
978 Manager::FunctionPtr ROE_PTheta(const std::vector<std::string>& arguments)
979 {
980 std::string maskName = RestOfEvent::c_defaultMaskName;
981
982 if (arguments.size() == 1)
983 maskName = arguments[0];
984 else if (arguments.size() > 1)
985 B2FATAL("At most 1 argument (name of mask) accepted for meta function ROE_PTheta");
986
987 auto func = [maskName](const Particle * particle) -> double {
988
989 // Get related ROE object
990 const RestOfEvent* roe = getRelatedROEObject(particle);
991
992 if (!roe)
993 {
994 B2ERROR("Relation between particle and ROE doesn't exist!");
995 return Const::doubleNaN;
996 }
997
998 const auto& frame = ReferenceFrame::GetCurrent();
999 auto frameRoe4Vector = frame.getMomentum(roe->get4Vector(maskName));
1000 return frameRoe4Vector.Theta();
1001 };
1002 return func;
1003 }
1004
1005 Manager::FunctionPtr ROE_DeltaE(const std::vector<std::string>& arguments)
1006 {
1007 std::string maskName = RestOfEvent::c_defaultMaskName;
1008
1009 if (arguments.size() == 1)
1010 maskName = arguments[0];
1011 else if (arguments.size() > 1)
1012 B2FATAL("At most 1 argument (name of mask) accepted for meta function ROE_deltae");
1013
1014 auto func = [maskName](const Particle * particle) -> double {
1015
1016 // Get related ROE object
1017 const RestOfEvent* roe = getRelatedROEObject(particle);
1018
1019 if (!roe)
1020 {
1021 B2ERROR("Relation between particle and ROE doesn't exist!");
1022 return Const::doubleNaN;
1023 }
1024
1025 PCmsLabTransform T;
1026 ROOT::Math::PxPyPzEVector vec = T.rotateLabToCms() * roe->get4Vector(maskName);
1027 return vec.E() - T.getCMSEnergy() / 2;
1028 };
1029 return func;
1030 }
1031
1032 Manager::FunctionPtr ROE_Mbc(const std::vector<std::string>& arguments)
1033 {
1034 std::string maskName = RestOfEvent::c_defaultMaskName;
1035
1036 if (arguments.size() == 1)
1037 maskName = arguments[0];
1038 else if (arguments.size() > 1)
1039 B2FATAL("At most 1 argument (name of mask) accepted for meta function ROE_mbc");
1040
1041 auto func = [maskName](const Particle * particle) -> double {
1042
1043 // Get related ROE object
1044 const RestOfEvent* roe = getRelatedROEObject(particle);
1045
1046 if (!roe)
1047 {
1048 B2ERROR("Relation between particle and ROE doesn't exist!");
1049 return Const::doubleNaN;
1050 }
1051
1052 PCmsLabTransform T;
1053 ROOT::Math::PxPyPzEVector vec = T.rotateLabToCms() * roe->get4Vector(maskName);
1054
1055 double E = T.getCMSEnergy() / 2;
1056 double m2 = E * E - vec.P2();
1057 double mbc = m2 > 0 ? sqrt(m2) : 0;
1058
1059 return mbc;
1060 };
1061 return func;
1062 }
1063
1064 Manager::FunctionPtr bssMassDifference(const std::vector<std::string>& arguments)
1065 {
1066 std::string maskName = RestOfEvent::c_defaultMaskName;
1067
1068 if (arguments.size() == 1)
1069 maskName = arguments[0];
1070 else if (arguments.size() > 1)
1071 B2FATAL("At most 1 argument (name of mask) accepted for meta function bssMassDifference");
1072
1073 auto func = [maskName](const Particle * particle) -> double {
1074
1075 // Get related ROE object
1076 ROOT::Math::PxPyPzEVector neutrino4vec = missing4Vector(particle->getDaughter(0), maskName, "6");
1077 ROOT::Math::PxPyPzEVector sig4vec = particle->getDaughter(0)->get4Vector();
1078
1079 ROOT::Math::PxPyPzEVector bsMom = neutrino4vec + sig4vec;
1080 ROOT::Math::PxPyPzEVector bssMom = bsMom + particle->getDaughter(1)->get4Vector();
1081
1082 return bssMom.M() - bsMom.M();
1083 };
1084 return func;
1085 }
1086
1087 Manager::FunctionPtr WE_DeltaE(const std::vector<std::string>& arguments)
1088 {
1089 std::string maskName;
1090 std::string opt;
1091
1092 if (arguments.size() == 1) {
1094 opt = arguments[0];
1095 } else if (arguments.size() == 2) {
1096 maskName = arguments[0];
1097 opt = arguments[1];
1098 } else
1099 B2FATAL("Wrong number of arguments (2 required) for meta function correctedB_deltae");
1100
1101 auto func = [maskName, opt](const Particle * particle) -> double {
1102
1103 PCmsLabTransform T;
1104 ROOT::Math::PxPyPzEVector boostvec = T.getBeamFourMomentum();
1105 ROOT::Math::PxPyPzEVector sig4vec = T.rotateLabToCms() * particle->get4Vector();
1106 ROOT::Math::PxPyPzEVector sig4vecLAB = particle->get4Vector();
1107 ROOT::Math::PxPyPzEVector neutrino4vec = missing4Vector(particle, maskName, "1");
1108 ROOT::Math::PxPyPzEVector neutrino4vecLAB = missing4Vector(particle, maskName, "6");
1109
1110 double deltaE = Const::doubleNaN;
1111
1112 // Definition 0: CMS
1113 if (opt == "0")
1114 {
1115 double totalSigEnergy = (sig4vec + neutrino4vec).energy();
1116 double E = T.getCMSEnergy() / 2;
1117 deltaE = totalSigEnergy - E;
1118 }
1119
1120 // Definition 1: LAB
1121 else if (opt == "1")
1122 {
1123 double Ecms = T.getCMSEnergy();
1124 double s = Ecms * Ecms;
1125 deltaE = ((sig4vecLAB + neutrino4vecLAB).Dot(boostvec) - s / 2.0) / sqrt(s);
1126 }
1127
1128 else
1129 B2FATAL("Option for correctedB_deltae variable should only be 0/1 (CMS/LAB)");
1130
1131 return deltaE;
1132 };
1133 return func;
1134 }
1135
1136 Manager::FunctionPtr WE_Mbc(const std::vector<std::string>& arguments)
1137 {
1138 std::string maskName;
1139 std::string opt;
1140
1141 if (arguments.size() == 1) {
1143 opt = arguments[0];
1144 } else if (arguments.size() == 2) {
1145 maskName = arguments[0];
1146 opt = arguments[1];
1147 } else
1148 B2FATAL("Wrong number of arguments (2 required) for meta function correctedB_mbc");
1149
1150 auto func = [maskName, opt](const Particle * particle) -> double {
1151
1152 PCmsLabTransform T;
1153 ROOT::Math::PxPyPzEVector boostvec = T.getBeamFourMomentum();
1154 ROOT::Math::PxPyPzEVector sig4vec = T.rotateLabToCms() * particle->get4Vector();
1155 ROOT::Math::PxPyPzEVector sig4vecLAB = particle->get4Vector();
1156 ROOT::Math::PxPyPzEVector neutrino4vec;
1157
1158 double mbc = Const::doubleNaN;
1159
1160 // Definition 0: CMS
1161 if (opt == "0")
1162 {
1163 neutrino4vec = missing4Vector(particle, maskName, "1");
1164 B2Vector3D bmom = (sig4vec + neutrino4vec).Vect();
1165 double E = T.getCMSEnergy() / 2;
1166 double m2 = E * E - bmom.Mag2();
1167 mbc = m2 > 0 ? sqrt(m2) : 0;
1168 }
1169
1170 // Definition 1: LAB
1171 else if (opt == "1")
1172 {
1173 neutrino4vec = missing4Vector(particle, maskName, "6");
1174 B2Vector3D bmom = (sig4vecLAB + neutrino4vec).Vect();
1175 double Ecms = T.getCMSEnergy();
1176 double s = Ecms * Ecms;
1177 double m2 = pow((s / 2.0 + bmom * B2Vector3D(boostvec.Vect())) / boostvec.energy(), 2.0) - bmom.Mag2();
1178 mbc = m2 > 0 ? sqrt(m2) : 0;
1179 }
1180
1181 // Definition 2: CMS with factor alpha (so that dE == 0)
1182 else if (opt == "2")
1183 {
1184 neutrino4vec = missing4Vector(particle, maskName, "7");
1185 B2Vector3D bmom = (sig4vec + neutrino4vec).Vect();
1186 double E = T.getCMSEnergy() / 2;
1187 double m2 = E * E - bmom.Mag2();
1188 mbc = m2 > 0 ? sqrt(m2) : 0;
1189 }
1190
1191 else
1192 B2FATAL("Option for correctedB_mbc variable should only be 0/1/2 (CMS/LAB/CMS with factor)");
1193
1194 return mbc;
1195 };
1196 return func;
1197 }
1198
1199 Manager::FunctionPtr WE_MissM2(const std::vector<std::string>& arguments)
1200 {
1201 std::string maskName;
1202 std::string opt;
1203
1204 if (arguments.size() == 1) {
1206 opt = arguments[0];
1207 } else if (arguments.size() == 2) {
1208 maskName = arguments[0];
1209 opt = arguments[1];
1210 } else
1211 B2FATAL("Wrong number of arguments (2 required) for meta function missM2");
1212
1213 auto func = [maskName, opt](const Particle * particle) -> double {
1214
1215 return missing4Vector(particle, maskName, opt).M2();
1216 };
1217 return func;
1218 }
1219
1220 Manager::FunctionPtr WE_MissPTheta(const std::vector<std::string>& arguments)
1221 {
1222 std::string maskName;
1223 std::string opt;
1224
1225 if (arguments.size() == 1) {
1227 opt = arguments[0];
1228 } else if (arguments.size() == 2) {
1229 maskName = arguments[0];
1230 opt = arguments[1];
1231 } else
1232 B2FATAL("Wrong number of arguments (2 required) for meta function WE_MissPTheta");
1233
1234 auto func = [maskName, opt](const Particle * particle) -> double {
1235
1236 // Get related ROE object
1237 const RestOfEvent* roe = getRelatedROEObject(particle);
1238
1239 if (!roe)
1240 {
1241 B2ERROR("Relation between particle and ROE doesn't exist!");
1242 return Const::doubleNaN;
1243 }
1244
1245 return missing4Vector(particle, maskName, opt).Theta();
1246 };
1247 return func;
1248 }
1249
1250 Manager::FunctionPtr WE_MissP(const std::vector<std::string>& arguments)
1251 {
1252 std::string maskName;
1253 std::string opt;
1254
1255 if (arguments.size() == 1) {
1257 opt = arguments[0];
1258 } else if (arguments.size() == 2) {
1259 maskName = arguments[0];
1260 opt = arguments[1];
1261 } else
1262 B2FATAL("Wrong number of arguments (2 required) for meta function WE_MissP");
1263
1264 auto func = [maskName, opt](const Particle * particle) -> double {
1265
1266 // Get related ROE object
1267 const RestOfEvent* roe = getRelatedROEObject(particle);
1268
1269 if (!roe)
1270 {
1271 B2ERROR("Relation between particle and ROE doesn't exist!");
1272 return Const::doubleNaN;
1273 }
1274
1275 return missing4Vector(particle, maskName, opt).P();
1276 };
1277 return func;
1278 }
1279
1280 Manager::FunctionPtr WE_MissPx(const std::vector<std::string>& arguments)
1281 {
1282 std::string maskName;
1283 std::string opt;
1284
1285 if (arguments.size() == 1) {
1287 opt = arguments[0];
1288 } else if (arguments.size() == 2) {
1289 maskName = arguments[0];
1290 opt = arguments[1];
1291 } else
1292 B2FATAL("Wrong number of arguments (2 required) for meta function WE_MissPx");
1293
1294 auto func = [maskName, opt](const Particle * particle) -> double {
1295
1296 // Get related ROE object
1297 const RestOfEvent* roe = getRelatedROEObject(particle);
1298
1299 if (!roe)
1300 {
1301 B2ERROR("Relation between particle and ROE doesn't exist!");
1302 return Const::doubleNaN;
1303 }
1304
1305 return missing4Vector(particle, maskName, opt).Px();
1306 };
1307 return func;
1308 }
1309
1310 Manager::FunctionPtr WE_MissPy(const std::vector<std::string>& arguments)
1311 {
1312 std::string maskName;
1313 std::string opt;
1314
1315 if (arguments.size() == 1) {
1317 opt = arguments[0];
1318 } else if (arguments.size() == 2) {
1319 maskName = arguments[0];
1320 opt = arguments[1];
1321 } else
1322 B2FATAL("Wrong number of arguments (2 required) for meta function WE_MissPy");
1323
1324 auto func = [maskName, opt](const Particle * particle) -> double {
1325
1326 // Get related ROE object
1327 const RestOfEvent* roe = getRelatedROEObject(particle);
1328
1329 if (!roe)
1330 {
1331 B2ERROR("Relation between particle and ROE doesn't exist!");
1332 return Const::doubleNaN;
1333 }
1334
1335 return missing4Vector(particle, maskName, opt).Py();
1336 };
1337 return func;
1338 }
1339
1340 Manager::FunctionPtr WE_MissPz(const std::vector<std::string>& arguments)
1341 {
1342 std::string maskName;
1343 std::string opt;
1344
1345 if (arguments.size() == 1) {
1347 opt = arguments[0];
1348 } else if (arguments.size() == 2) {
1349 maskName = arguments[0];
1350 opt = arguments[1];
1351 } else
1352 B2FATAL("Wrong number of arguments (2 required) for meta function WE_MissPz");
1353
1354 auto func = [maskName, opt](const Particle * particle) -> double {
1355
1356 // Get related ROE object
1357 const RestOfEvent* roe = getRelatedROEObject(particle);
1358
1359 if (!roe)
1360 {
1361 B2ERROR("Relation between particle and ROE doesn't exist!");
1362 return Const::doubleNaN;
1363 }
1364
1365 return missing4Vector(particle, maskName, opt).Pz();
1366 };
1367 return func;
1368 }
1369
1370 Manager::FunctionPtr WE_MissE(const std::vector<std::string>& arguments)
1371 {
1372 std::string maskName;
1373 std::string opt;
1374
1375 if (arguments.size() == 1) {
1377 opt = arguments[0];
1378 } else if (arguments.size() == 2) {
1379 maskName = arguments[0];
1380 opt = arguments[1];
1381 } else
1382 B2FATAL("Wrong number of arguments (2 required) for meta function WE_MissE");
1383
1384 auto func = [maskName, opt](const Particle * particle) -> double {
1385
1386 // Get related ROE object
1387 const RestOfEvent* roe = getRelatedROEObject(particle);
1388
1389 if (!roe)
1390 {
1391 B2ERROR("Relation between particle and ROE doesn't exist!");
1392 return Const::doubleNaN;
1393 }
1394
1395 return missing4Vector(particle, maskName, opt).energy();
1396 };
1397 return func;
1398 }
1399
1400 Manager::FunctionPtr WE_xiZ(const std::vector<std::string>& arguments)
1401 {
1402 std::string maskName = RestOfEvent::c_defaultMaskName;
1403
1404 if (arguments.size() == 1)
1405 maskName = arguments[0];
1406 else if (arguments.size() > 1)
1407 B2FATAL("At most 1 argument (name of mask) accepted for meta function xiZ");
1408
1409 auto func = [maskName](const Particle * particle) -> double {
1410
1411 // Get related ROE object
1412 const RestOfEvent* roe = getRelatedROEObject(particle);
1413
1414 if (!roe)
1415 {
1416 B2ERROR("Relation between particle and ROE doesn't exist!");
1417 return Const::doubleNaN;
1418 }
1419
1420 double pz = 0;
1421 double energy = 0;
1422
1423 // Get all Tracks on reconstructed side
1424 std::vector<const Particle*> recTrackParticles = particle->getFinalStateDaughters();
1425
1426 // Loop the reconstructed side
1427 for (auto& recTrackParticle : recTrackParticles)
1428 {
1429 pz += recTrackParticle->getPz();
1430 energy += recTrackParticle->getEnergy();
1431 }
1432
1433 // Loop the ROE side
1434 auto roeParticles = roe->getChargedParticles(maskName);
1435 for (auto* roeParticle : roeParticles)
1436 {
1437 pz += roeParticle->getPz();
1438 energy += roeParticle->getEnergy();
1439 }
1440
1441 return pz / energy;
1442 };
1443 return func;
1444 }
1445
1446 Manager::FunctionPtr WE_MissM2OverMissE(const std::vector<std::string>& arguments)
1447 {
1448 std::string maskName = RestOfEvent::c_defaultMaskName;
1449
1450 if (arguments.size() == 1)
1451 maskName = arguments[0];
1452 else if (arguments.size() > 1)
1453 B2FATAL("At most 1 argument (name of mask) accepted for meta function WE_MissM2OverMissE");
1454
1455 auto func = [maskName](const Particle * particle) -> double {
1456
1457 // Get related ROE object
1458 const RestOfEvent* roe = getRelatedROEObject(particle);
1459
1460 if (!roe)
1461 {
1462 B2ERROR("Relation between particle and ROE doesn't exist!");
1463 return Const::doubleNaN;
1464 }
1465
1466 return missing4Vector(particle, maskName, "5").M2() / (2.0 * missing4Vector(particle, maskName, "5").energy());
1467 };
1468 return func;
1469 }
1470
1471 Manager::FunctionPtr WE_q2lnuSimple(const std::vector<std::string>& arguments)
1472 {
1473 std::string maskName(RestOfEvent::c_defaultMaskName);
1474 std::string option("1");
1475
1476 if (arguments.size() == 1) {
1477 maskName = arguments[0];
1478 } else if (arguments.size() == 2) {
1479 maskName = arguments[0];
1480 option = arguments[1];
1481 } else if (arguments.size() > 2) {
1482 B2FATAL("Too many arguments. At most two arguments are allowed for meta function q2lnuSimple(maskname,option)");
1483 }
1484
1485 auto func = [maskName, option](const Particle * particle) -> double {
1486
1487 // Get related ROE object
1488 const RestOfEvent* roe = getRelatedROEObject(particle);
1489
1490 if (!roe)
1491 {
1492 B2ERROR("Relation between particle and ROE doesn't exist!");
1493 return Const::doubleNaN;
1494 }
1495
1496 int n = particle->getNDaughters();
1497
1498 if (n < 1)
1499 return Const::doubleNaN;
1500
1501 // Assumes lepton is the last particle in the reco decay chain!
1502 PCmsLabTransform T;
1503 const Particle* lep = particle->getDaughter(n - 1);
1504 ROOT::Math::PxPyPzEVector lep4vec = T.rotateLabToCms() * lep->get4Vector();
1505 ROOT::Math::PxPyPzEVector nu4vec = missing4Vector(particle, maskName, option);
1506
1507 return (lep4vec + nu4vec).M2();
1508 };
1509 return func;
1510 }
1511
1512 Manager::FunctionPtr WE_q2lnu(const std::vector<std::string>& arguments)
1513 {
1514 std::string maskName(RestOfEvent::c_defaultMaskName);
1515 std::string option("7");
1516
1517 if (arguments.size() == 1) {
1518 maskName = arguments[0];
1519 } else if (arguments.size() == 2) {
1520 maskName = arguments[0];
1521 option = arguments[1];
1522 } else if (arguments.size() > 2) {
1523 B2FATAL("Too many arguments. At most two arguments are allowed for meta function q2lnu(maskname, option)");
1524 }
1525
1526 auto func = [maskName, option](const Particle * particle) -> double {
1527
1528 // Get related ROE object
1529 const RestOfEvent* roe = getRelatedROEObject(particle);
1530
1531 if (!roe)
1532 {
1533 B2ERROR("Relation between particle and ROE doesn't exist!");
1534 return Const::doubleNaN;
1535 }
1536
1537 int n = particle->getNDaughters();
1538
1539 if (n < 1)
1540 return Const::doubleNaN;
1541
1542 PCmsLabTransform T;
1543 const Particle* lep = particle->getDaughter(n - 1);
1544 ROOT::Math::PxPyPzEVector lep_cm = T.rotateLabToCms() * lep->get4Vector();
1545
1546 ROOT::Math::PxPyPzEVector Y_cm = T.rotateLabToCms() * particle->get4Vector();
1547 ROOT::Math::PxPyPzEVector neu_cm = missing4Vector(particle, maskName, option);
1548
1549 double e_beam = T.getCMSEnergy() / 2.0;
1550
1551 //just to make the formula simpler
1552 double bmass = particle->getPDGMass();
1553 double pB2 = e_beam * e_beam - bmass * bmass;
1554
1555 //angle between the Y and the neutrino, from the Mbc=M_B constraint
1556 double cos_angle_nu = (pB2 - Y_cm.P2() - neu_cm.P2()) / (2.0 * Y_cm.P() * neu_cm.P());
1557 if (abs(cos_angle_nu) > 1)
1558 {
1559 return (lep_cm + neu_cm).M2();
1560 }
1561
1562 double angle_nu = TMath::ACos(cos_angle_nu);
1563 //first get one random neutrino, on the allowed cone for the constraint
1564 ROOT::Math::PtEtaPhiEVector rotated_neu(-Y_cm); //first get reverse Y
1565 rotated_neu.SetE(Y_cm.E());
1566
1567 double nu_eta = -log(tan((rotated_neu.Theta() + (TMath::Pi() - angle_nu)) / 2.));
1568 rotated_neu.SetEta(nu_eta);
1569 rotated_neu.SetPt(neu_cm.pt());
1570 rotated_neu.SetE(neu_cm.E());
1571
1572 B2Vector3D Yneu_norm = Y_cm.Vect().Cross(neu_cm.Vect());
1573 B2Vector3D Yrot_norm = Y_cm.Vect().Cross(rotated_neu.Vect());
1574 //angle between the two crossproducts -> angle between the two vectors perpendicular to the Y-p_inc and Y-B planes -> angle between the planes
1575 //this angle needs to come out as zero
1576
1577 double rot_angle = Yneu_norm.Angle(Yrot_norm);
1578
1579 ROOT::Math::PxPyPzEVector rotated_neu2(rotated_neu);
1580 //unfortunately don't -and probably can't- know in which direction to rotate without trying
1581 //so create a copy of the vector, and later choose the correct one
1582 //However, rotation by 180 degrees is never needed, direction of the cross-product vector assures that when after rotation
1583 //the B-vector is in the plane, it always is on the side closer to pcm_lv_inc.
1584 //rotate around Y into the Y-neutrino-plane
1585 ROOT::Math::AxisAngle rotation(Y_cm.Vect(), rot_angle);
1586 rotation(rotated_neu);
1587 rotation.Invert();
1588 rotation(rotated_neu2);
1589
1590 double dot1 = rotated_neu.Vect().Dot(Yneu_norm);
1591 double dot2 = rotated_neu2.Vect().Dot(Yneu_norm);
1592
1593 if (abs(dot2) < abs(dot1)) rotated_neu = rotated_neu2;
1594
1595 return (lep_cm + rotated_neu).M2();
1596 };
1597 return func;
1598 }
1599
1600 Manager::FunctionPtr WE_cosThetaEll(const std::vector<std::string>& arguments)
1601 {
1602 std::string maskName = RestOfEvent::c_defaultMaskName;
1603
1604 if (arguments.size() == 1)
1605 maskName = arguments[0];
1606 else if (arguments.size() > 1)
1607 B2FATAL("At most 1 argument (name of mask) accepted for meta function cosThetaEll");
1608
1609 auto func = [maskName](const Particle * particle) -> double {
1610
1611 ROOT::Math::PxPyPzEVector pNu = missing4Vector(particle, maskName, "6");
1612
1613 ROOT::Math::PxPyPzEVector pLep;
1614 // TODO: avoid hardocoded values
1615 for (unsigned i = 0; i < particle->getNDaughters(); i++)
1616 {
1617 int absPDG = abs(particle->getDaughter(i)->getPDGCode());
1618 if (absPDG == Const::electron.getPDGCode() || absPDG == Const::muon.getPDGCode() || absPDG == 15) {
1619 pLep = particle->getDaughter(i)->get4Vector();
1620 break;
1621 }
1622 }
1623
1624 ROOT::Math::PxPyPzEVector pW = pNu + pLep;
1625 ROOT::Math::PxPyPzEVector pB = particle->get4Vector() + pNu;
1626
1627 // boost lepton and B momentum to W frame
1628 B2Vector3D boost2W = pW.BoostToCM();
1629 pLep = ROOT::Math::Boost(boost2W) * pLep;
1630 pB = ROOT::Math::Boost(boost2W) * pB;
1631
1632 B2Vector3D lep3Vector = pLep.Vect();
1633 B2Vector3D B3Vector = pB.Vect();
1634 double numerator = lep3Vector.Dot(B3Vector);
1635 double denominator = (lep3Vector.Mag()) * (B3Vector.Mag());
1636
1637 return numerator / denominator;
1638 };
1639 return func;
1640 }
1641
1642 Manager::FunctionPtr passesROEMask(const std::vector<std::string>& arguments)
1643 {
1644 std::string maskName = RestOfEvent::c_defaultMaskName;
1645
1646 if (arguments.size() == 1)
1647 maskName = arguments[0];
1648 else if (arguments.size() > 1)
1649 B2FATAL("At most 1 argument (name of mask) accepted for meta function passesROEMask");
1650
1651 auto func = [maskName](const Particle * particle) -> double {
1652
1653 double result = 0;
1654
1655 StoreObjPtr<RestOfEvent> roe("RestOfEvent");
1656 if (not roe.isValid())
1657 return Const::doubleNaN;
1658
1659 if (roe->hasParticle(particle, maskName))
1660 {
1661 return 1.0;
1662 }
1663
1664 return result;
1665 };
1666 return func;
1667 }
1668
1669 double printROE(const Particle* particle)
1670 {
1671 const RestOfEvent* roe = getRelatedROEObject(particle);
1672
1673 if (!roe) {
1674 B2ERROR("Relation between particle and ROE doesn't exist!");
1675 } else roe->print();
1676 return 0;
1677 }
1678
1679 double hasCorrectROECombination(const Particle* particle)
1680 {
1681 unsigned nDaughters = particle->getNDaughters();
1682 if (nDaughters < 2) {
1683 B2ERROR("The particle must have at least two daughters.");
1684 return Const::doubleNaN;
1685 }
1686
1687 for (unsigned i = 0; i < particle->getNDaughters(); i++) {
1688
1689 // Find a daughter that is loaded from a ROE object
1690 auto daughter = particle->getDaughter(i);
1691 auto roe = daughter->getRelatedFrom<RestOfEvent>();
1692 if (!roe)
1693 continue;
1694
1695 auto sourceParticle = roe->getRelatedFrom<Particle>();
1696 for (unsigned j = 0; j < particle->getNDaughters(); j++) {
1697 if (i == j) continue;
1698 const auto anotherDaughter = particle->getDaughter(j);
1699
1700 if (anotherDaughter == sourceParticle)
1701 return 1.0;
1702 }
1703 return 0.0;
1704 }
1705
1706 B2ERROR("There is no daughter particle loaded from the ROE object.");
1707 return Const::doubleNaN;
1708 }
1709
1710 Manager::FunctionPtr pi0Prob(const std::vector<std::string>& arguments)
1711 {
1712 if (arguments.size() != 1)
1713 B2ERROR("Wrong number of arguments (1 required) for pi0Prob");
1714
1715 std::string mode;
1716 mode = arguments[0];
1717
1718 if (mode != "standard" and mode != "tight" and mode != "cluster" and mode != "both")
1719 B2ERROR("the given argument is not supported in pi0Prob!");
1720
1721 auto func = [mode](const Particle * particle) -> double {
1722 if (mode == "standard")
1723 {
1724 if (particle->hasExtraInfo("Pi0ProbOrigin")) {
1725 return particle->getExtraInfo("Pi0ProbOrigin");
1726 } else {
1727 B2WARNING("Pi0ProbOrigin is not registerted in extraInfo! \n"
1728 "the function writePi0EtaVeto has to be executed to register this extraInfo.");
1729 return Const::doubleNaN;
1730 }
1731 } else if (mode == "tight")
1732 {
1733 if (particle->hasExtraInfo("Pi0ProbTightEnergyThreshold")) {
1734 return particle->getExtraInfo("Pi0ProbTightEnergyThreshold");
1735 } else {
1736 B2WARNING("Pi0ProbTightEnergyThreshold is not registerted in extraInfo! \n"
1737 "the function writePi0EtaVeto has to be executed to register this extraInfo.");
1738 return Const::doubleNaN;
1739 }
1740 } else if (mode == "cluster")
1741 {
1742 if (particle->hasExtraInfo("Pi0ProbLargeClusterSize")) {
1743 return particle->getExtraInfo("Pi0ProbLargeClusterSize");
1744 } else {
1745 B2WARNING("Pi0ProbLargeClusterSize is not registerted in extraInfo! \n"
1746 "the function writePi0EtaVeto has to be executed to register this extraInfo.");
1747 return Const::doubleNaN;
1748 }
1749 } else if (mode == "both")
1750 {
1751 if (particle->hasExtraInfo("Pi0ProbTightEnergyThresholdAndLargeClusterSize")) {
1752 return particle->getExtraInfo("Pi0ProbTightEnergyThresholdAndLargeClusterSize");
1753 } else {
1754 B2WARNING("Pi0ProbTightEnergyThresholdAndLargeClusterSize is not registerted in extraInfo! \n"
1755 "the function writePi0EtaVeto has to be executed to register this extraInfo.");
1756 return Const::doubleNaN;
1757 }
1758 } else
1759 {
1760 return Const::doubleNaN;
1761 }
1762 };
1763 return func;
1764 }
1765
1766 Manager::FunctionPtr etaProb(const std::vector<std::string>& arguments)
1767 {
1768 if (arguments.size() != 1)
1769 B2ERROR("Wrong number of arguments (1 required) for etaProb");
1770
1771 std::string mode;
1772 mode = arguments[0];
1773
1774 if (mode != "standard" and mode != "tight" and mode != "cluster" and mode != "both")
1775 B2ERROR("the given argument is not supported in etaProb!");
1776
1777 auto func = [mode](const Particle * particle) -> double {
1778 if (mode == "standard")
1779 {
1780 if (particle->hasExtraInfo("EtaProbOrigin")) {
1781 return particle->getExtraInfo("EtaProbOrigin");
1782 } else {
1783 B2WARNING("EtaProbOrigin is not registerted in extraInfo! \n"
1784 "the function writePi0EtaVeto has to be executed to register this extraInfo.");
1785 return Const::doubleNaN;
1786 }
1787 } else if (mode == "tight")
1788 {
1789 if (particle->hasExtraInfo("EtaProbTightEnergyThreshold")) {
1790 return particle->getExtraInfo("EtaProbTightEnergyThreshold");
1791 } else {
1792 B2WARNING("EtaProbTightEnergyThreshold is not registerted in extraInfo! \n"
1793 "the function writePi0EtaVeto has to be executed to register this extraInfo.");
1794 return Const::doubleNaN;
1795 }
1796 } else if (mode == "cluster")
1797 {
1798 if (particle->hasExtraInfo("EtaProbLargeClusterSize")) {
1799 return particle->getExtraInfo("EtaProbLargeClusterSize");
1800 } else {
1801 B2WARNING("EtaProbLargeClusterSize is not registerted in extraInfo! \n"
1802 "the function writePi0EtaVeto has to be executed to register this extraInfo.");
1803 return Const::doubleNaN;
1804 }
1805 } else if (mode == "both")
1806 {
1807 if (particle->hasExtraInfo("EtaProbTightEnergyThresholdAndLargeClusterSize")) {
1808 return particle->getExtraInfo("EtaProbTightEnergyThresholdAndLargeClusterSize");
1809 } else {
1810 B2WARNING("EtaProbTightEnergyThresholdAndLargeClusterSize is not registerted in extraInfo! \n"
1811 "the function writePi0EtaVeto has to be executed to register this extraInfo.");
1812 return Const::doubleNaN;
1813 }
1814 } else
1815 {
1816 return Const::doubleNaN;
1817 }
1818 };
1819 return func;
1820 }
1821
1822 // ------------------------------------------------------------------------------
1823 // Below are some functions for ease of usage, they are not a part of variables
1824 // ------------------------------------------------------------------------------
1825
1826 ROOT::Math::PxPyPzEVector missing4Vector(const Particle* particle, const std::string& maskName, const std::string& opt)
1827 {
1828 // Get related ROE object
1829 const RestOfEvent* roe = getRelatedROEObject(particle);
1830
1831 if (!roe) {
1832 B2ERROR("Relation between particle and ROE doesn't exist!");
1833 ROOT::Math::PxPyPzEVector empty;
1834 return empty;
1835 }
1836
1837 PCmsLabTransform T;
1838 ROOT::Math::PxPyPzEVector boostvec = T.getBeamFourMomentum();
1839
1840 ROOT::Math::PxPyPzEVector rec4vecLAB = particle->get4Vector();
1841 ROOT::Math::PxPyPzEVector roe4vecLAB = roe->get4Vector(maskName);
1842
1843 ROOT::Math::PxPyPzEVector rec4vec = T.rotateLabToCms() * rec4vecLAB;
1844 ROOT::Math::PxPyPzEVector roe4vec = T.rotateLabToCms() * roe4vecLAB;
1845
1846 ROOT::Math::PxPyPzEVector miss4vec;
1847 double E_beam_cms = T.getCMSEnergy() / 2.0;
1848
1849 // Definition 0: CMS, use energy and momentum of tracks and clusters
1850 if (opt == "0") {
1851 miss4vec = - (rec4vec + roe4vec);
1852 miss4vec.SetE(2 * E_beam_cms - (rec4vec.E() + roe4vec.E()));
1853 }
1854
1855 // Definition 1: CMS, same as 0, fix Emiss = pmiss
1856 else if (opt == "1") {
1857 miss4vec = - (rec4vec + roe4vec);
1858 miss4vec.SetE(miss4vec.P());
1859 }
1860
1861 // Definition 2: CMS, same as 0, fix Eroe = Ecms/2
1862 else if (opt == "2") {
1863 miss4vec = - (rec4vec + roe4vec);
1864 miss4vec.SetE(E_beam_cms - rec4vec.E());
1865 }
1866
1867 // Definition 3: CMS, use only energy and momentum of signal side
1868 else if (opt == "3") {
1869 miss4vec = - rec4vec;
1870 miss4vec.SetE(E_beam_cms - rec4vec.E());
1871 }
1872
1873 // Definition 4: CMS, same as 3, update with direction of ROE momentum
1874 else if (opt == "4") {
1875 B2Vector3D pB = - roe4vec.Vect();
1876 pB.SetMag(0.340);
1877 pB -= rec4vec.Vect();
1878 miss4vec.SetPxPyPzE(pB.X(), pB.Y(), pB.Z(), E_beam_cms - rec4vec.E());
1879 }
1880
1881 // Definition 5: LAB, use energy and momentum of tracks and clusters from whole event
1882 else if (opt == "5") {
1883 miss4vec = boostvec - (rec4vecLAB + roe4vecLAB);
1884 }
1885
1886 // Definition 6: LAB, same as 5, fix Emiss = pmiss
1887 else if (opt == "6") {
1888 miss4vec = boostvec - (rec4vecLAB + roe4vecLAB);
1889 miss4vec.SetE(miss4vec.P());
1890 }
1891
1892 // Definition 7: CMS, correct pmiss 3-momentum vector with factor alpha so that dE = 0 (used for Mbc calculation)
1893 else if (opt == "7") {
1894 miss4vec = - (rec4vec + roe4vec);
1895 miss4vec.SetE(miss4vec.P());
1896 double factorAlpha = (E_beam_cms - rec4vec.E()) / miss4vec.E();
1897 miss4vec *= factorAlpha;
1898 miss4vec.SetE(miss4vec.P());
1899 }
1900
1901 return miss4vec;
1902 }
1903
1904 void checkMCParticleMissingFlags(const MCParticle* mcp, std::set<const MCParticle*> mcROEObjects, int& missingFlags)
1905 {
1906 std::vector<MCParticle*> daughters = mcp->getDaughters();
1907 for (auto& daughter : daughters) {
1908
1909 if (!daughter->hasStatus(MCParticle::c_PrimaryParticle))
1910 continue;
1911
1912 if (mcROEObjects.find(daughter) == mcROEObjects.end()) {
1913
1914 int pdg = abs(daughter->getPDG());
1915
1916 // photon
1917 if (pdg == Const::photon.getPDGCode() and (missingFlags & 1) == 0)
1918 missingFlags += 1;
1919
1920 // electrons
1921 else if (pdg == Const::electron.getPDGCode() and (missingFlags & 2) == 0)
1922 missingFlags += 2;
1923
1924 // muons
1925 else if (pdg == Const::muon.getPDGCode() and (missingFlags & 4) == 0)
1926 missingFlags += 4;
1927
1928 // pions
1929 else if (pdg == Const::pion.getPDGCode() and (missingFlags & 8) == 0)
1930 missingFlags += 8;
1931
1932 // kaons
1933 else if (pdg == Const::kaon.getPDGCode() and (missingFlags & 16) == 0)
1934 missingFlags += 16;
1935
1936 // protons
1937 else if (pdg == Const::proton.getPDGCode() and (missingFlags & 32) == 0)
1938 missingFlags += 32;
1939
1940 // neutrons
1941 else if (pdg == Const::neutron.getPDGCode() and (missingFlags & 64) == 0)
1942 missingFlags += 64;
1943
1944 // kshort
1945 else if (pdg == Const::Kshort.getPDGCode() and ((missingFlags & 128) == 0 or (missingFlags & 256) == 0)) {
1946 std::vector<MCParticle*> ksDaug = daughter->getDaughters();
1947 if (ksDaug.size() == 2) {
1948 // K_S0 -> pi+ pi-
1949 if (abs(ksDaug[0]->getPDG()) == Const::pion.getPDGCode() and abs(ksDaug[1]->getPDG()) == Const::pion.getPDGCode()
1950 and (missingFlags & 128) == 0) {
1951 if (mcROEObjects.find(ksDaug[0]) == mcROEObjects.end() or mcROEObjects.find(ksDaug[1]) == mcROEObjects.end())
1952 missingFlags += 128;
1953 }
1954 // K_S0 -> pi0 pi0
1955 else if (abs(ksDaug[0]->getPDG()) == Const::pi0.getPDGCode() and abs(ksDaug[1]->getPDG()) == Const::pi0.getPDGCode()
1956 and (missingFlags & 256) == 0) {
1957 std::vector<MCParticle*> pi0Daug0 = ksDaug[0]->getDaughters();
1958 std::vector<MCParticle*> pi0Daug1 = ksDaug[1]->getDaughters();
1959 if (mcROEObjects.find(pi0Daug0[0]) == mcROEObjects.end() or
1960 mcROEObjects.find(pi0Daug0[1]) == mcROEObjects.end() or
1961 mcROEObjects.find(pi0Daug1[0]) == mcROEObjects.end() or
1962 mcROEObjects.find(pi0Daug1[1]) == mcROEObjects.end())
1963 missingFlags += 256;
1964 }
1965 }
1966 }
1967
1968 // klong
1969 else if (pdg == Const::Klong.getPDGCode() and (missingFlags & 512) == 0)
1970 missingFlags += 512;
1971
1972 // neutrinos, which are not in the Const::
1973 else if ((pdg == 12 or pdg == 14 or pdg == 16) and (missingFlags & 1024) == 0)
1974 missingFlags += 1024;
1975 }
1976 checkMCParticleMissingFlags(daughter, mcROEObjects, missingFlags);
1977 }
1978 }
1979
1980 double isInThisRestOfEvent(const Particle* particle, const RestOfEvent* roe, const std::string& maskName)
1981 {
1982 if (particle->getParticleSource() == Particle::c_Composite or
1983 particle->getParticleSource() == Particle::c_V0) {
1984 std::vector<const Particle*> fspDaug = particle->getFinalStateDaughters();
1985 for (auto& i : fspDaug) {
1986 if (isInThisRestOfEvent(i, roe, maskName) == 0)
1987 return 0;
1988 }
1989 return 1;
1990 }
1991 return roe->hasParticle(particle, maskName);
1992 }
1993
1994 const RestOfEvent* getRelatedROEObject(const Particle* particle, bool returnHostOnly)
1995 {
1996 // Get related ROE object
1997 const RestOfEvent* roe = particle->getRelatedTo<RestOfEvent>();
1998 if (!roe && !returnHostOnly) {
1999 roe = particle->getRelatedTo<RestOfEvent>("NestedRestOfEvents");
2000
2001 }
2002 return roe;
2003 }
2004
2005 VARIABLE_GROUP("Rest Of Event");
2006
2007 REGISTER_METAVARIABLE("useROERecoilFrame(variable)", useROERecoilFrame,
2008 "Returns the value of the variable using the rest frame of the ROE recoil as current reference frame.\n"
2009 "Can be used inside for_each loop or outside of it if the particle has associated Rest of Event.\n"
2010 "E.g. ``useROERecoilFrame(E)`` returns the energy of a particle in the ROE recoil frame.", Manager::VariableDataType::c_double);
2011
2012 REGISTER_VARIABLE("isInRestOfEvent", isInRestOfEvent,
2013 "Returns 1 if a track, ecl or klmCluster associated to particle is in the current RestOfEvent object, 0 otherwise."
2014 "One can use this variable only in a for_each loop over the RestOfEvent StoreArray.");
2015
2016 REGISTER_VARIABLE("isCloneOfSignalSide", isCloneOfSignalSide,
2017 "Returns 1 if a particle is a clone of signal side final state particles, 0 otherwise. "
2018 "Requires generator information and truth-matching. "
2019 "One can use this variable only in a ``for_each`` loop over the RestOfEvent StoreArray.");
2020
2021 REGISTER_VARIABLE("hasAncestorFromSignalSide", hasAncestorFromSignalSide,
2022 "Returns 1 if a particle has ancestor from signal side, 0 otherwise. "
2023 "Requires generator information and truth-matching. "
2024 "One can use this variable only in a ``for_each`` loop over the RestOfEvent StoreArray.");
2025
2026 REGISTER_METAVARIABLE("currentROEIsInList(particleList)", currentROEIsInList,
2027 "[Eventbased] Returns 1 the associated particle of the current ROE is contained in the given list or its charge-conjugated."
2028 "Useful to restrict the for_each loop over ROEs to ROEs of a certain ParticleList.", Manager::VariableDataType::c_double);
2029
2030 REGISTER_VARIABLE("nROE_RemainingTracks", nROE_RemainingTracks,
2031 "Returns number of tracks in ROE - number of tracks of given particle"
2032 "One can use this variable only in a for_each loop over the RestOfEvent StoreArray.");
2033
2034 REGISTER_METAVARIABLE("nROE_RemainingTracks(maskName)", nROE_RemainingTracksWithMask,
2035 "Returns number of remaining tracks between the ROE (specified via a mask) and the given particle. For the given particle only tracks are counted which are in the RoE."
2036 "One can use this variable only in a for_each loop over the RestOfEvent StoreArray."
2037 "Is required for the specific FEI. :noindex:", Manager::VariableDataType::c_int);
2038 // nROE_RemainingTracks is overloaded (two C++ functions sharing one
2039 // variable name) so one of the two needs to be made the indexed
2040 // variable in sphinx
2041
2042 REGISTER_VARIABLE("nROE_KLMClusters", nROE_KLMClusters,
2043 "Returns number of all remaining KLM clusters in the related RestOfEvent object.");
2044
2045 REGISTER_METAVARIABLE("nROE_Charged(maskName, PDGcode = 0)", nROE_ChargedParticles,
2046 "Returns number of all charged particles in the related RestOfEvent object. First optional argument is ROE mask name. "
2047 "Second argument is a PDG code to count only one charged particle species, independently of charge. "
2048 "For example: ``nROE_Charged(cleanMask, 321)`` will output number of kaons in Rest Of Event with ``cleanMask``. "
2049 "PDG code 0 is used to count all charged particles", Manager::VariableDataType::c_double);
2050
2051 REGISTER_METAVARIABLE("nROE_Photons(maskName)", nROE_Photons,
2052 "Returns number of all photons in the related RestOfEvent object, accepts 1 optional argument of ROE mask name. ",
2053 Manager::VariableDataType::c_double);
2054
2055 REGISTER_METAVARIABLE("nROE_NeutralHadrons(maskName)", nROE_NeutralHadrons,
2056 "Returns number of all neutral hadrons in the related RestOfEvent object, accepts 1 optional argument of ROE mask name. ",
2057 Manager::VariableDataType::c_double);
2058
2059 REGISTER_METAVARIABLE("particleRelatedToCurrentROE(var)", particleRelatedToCurrentROE,
2060 "[Eventbased] Returns variable applied to the particle which is related to the current RestOfEvent object"
2061 "One can use this variable only in a for_each loop over the RestOfEvent StoreArray.", Manager::VariableDataType::c_double);
2062
2063 REGISTER_VARIABLE("roeMC_E", ROE_MC_E,
2064 "Returns true energy of unused tracks and clusters in ROE, can be used with ``use***Frame()`` function.\n\n", "GeV");
2065
2066 REGISTER_VARIABLE("roeMC_M", ROE_MC_M,
2067 "Returns true invariant mass of unused tracks and clusters in ROE\n\n", "GeV/:math:`\\text{c}^2`");
2068
2069 REGISTER_VARIABLE("roeMC_P", ROE_MC_P,
2070 "Returns true momentum of unused tracks and clusters in ROE, can be used with ``use***Frame()`` function.\n\n", "GeV/c");
2071
2072 REGISTER_VARIABLE("roeMC_Px", ROE_MC_Px,
2073 "Returns x component of true momentum of unused tracks and clusters in ROE, can be used with ``use***Frame()`` function.\n\n",
2074 "GeV/c");
2075
2076 REGISTER_VARIABLE("roeMC_Py", ROE_MC_Py,
2077 "Returns y component of true momentum of unused tracks and clusters in ROE, can be used with ``use***Frame()`` function.\n\n",
2078 "GeV/c");
2079
2080 REGISTER_VARIABLE("roeMC_Pz", ROE_MC_Pz,
2081 "Returns z component of true momentum of unused tracks and clusters in ROE, can be used with ``use***Frame()`` function.\n\n",
2082 "GeV/c");
2083
2084 REGISTER_VARIABLE("roeMC_Pt", ROE_MC_Pt,
2085 "Returns transverse component of true momentum of unused tracks and clusters in ROE, can be used with ``use***Frame()`` function.\n\n",
2086 "GeV/c");
2087
2088 REGISTER_VARIABLE("roeMC_PTheta", ROE_MC_PTheta,
2089 "Returns polar angle of true momentum of unused tracks and clusters in ROE, can be used with ``use***Frame()`` function.\n\n",
2090 "rad");
2091
2092 REGISTER_METAVARIABLE("roeMC_MissFlags(maskName)", ROE_MC_MissingFlags,
2093 "Returns flags corresponding to missing particles on ROE side.", Manager::VariableDataType::c_double);
2094
2095 REGISTER_METAVARIABLE("nROE_Tracks(maskName)", nROE_Tracks,
2096 "Returns number of tracks in the related RestOfEvent object that pass the selection criteria.",
2097 Manager::VariableDataType::c_double);
2098
2099 REGISTER_METAVARIABLE("nROE_ECLClusters(maskName)", nROE_ECLClusters,
2100 "Returns number of ECL clusters in the related RestOfEvent object that pass the selection criteria.",
2101 Manager::VariableDataType::c_double);
2102
2103 REGISTER_METAVARIABLE("nROE_NeutralECLClusters(maskName)", nROE_NeutralECLClusters,
2104 "Returns number of neutral ECL clusters in the related RestOfEvent object that pass the selection criteria.",
2105 Manager::VariableDataType::c_double);
2106
2107 REGISTER_METAVARIABLE("nROE_Composites(maskName)", nROE_Composites,
2108 "Returns number of composite particles or V0s in the related RestOfEvent object that pass the selection criteria.",
2109 Manager::VariableDataType::c_double);
2110
2111 REGISTER_METAVARIABLE("nROE_ParticlesInList(pListName[, maskName])", nROE_ParticlesInList,
2112 "Returns the number of particles in ROE from the given particle list. If a mask name is provided the selection criteria are applied.\n"
2113 "Use of variable aliases is advised.", Manager::VariableDataType::c_double);
2114
2115 REGISTER_METAVARIABLE("roeCharge(maskName)", ROE_Charge,
2116 "Returns total charge of the related RestOfEvent object. The unit of the charge is ``e`` ", Manager::VariableDataType::c_double);
2117
2118 REGISTER_METAVARIABLE("roeEextra(maskName)", ROE_ExtraEnergy,
2119 "Returns extra energy from ECLClusters in the calorimeter that is not associated to the given Particle. The unit of the energy is ``GeV`` ",
2120 Manager::VariableDataType::c_double);
2121
2122 REGISTER_METAVARIABLE("roeNeextra(maskName)", ROE_NeutralExtraEnergy,
2123 "Returns extra energy from neutral ECLClusters in the calorimeter that is not associated to the given Particle, can be used with ``use***Frame()`` function. The unit of the energy is ``GeV`` ",
2124 Manager::VariableDataType::c_double);
2125
2126 REGISTER_METAVARIABLE("roeE(maskName)", ROE_E,
2127 "Returns energy of unused tracks and clusters in ROE, can be used with ``use***Frame()`` function. The unit of the energy is ``GeV`` ",
2128 Manager::VariableDataType::c_double);
2129
2130 REGISTER_METAVARIABLE("roeM(maskName)", ROE_M,
2131 "Returns invariant mass of unused tracks and clusters in ROE. The unit of the invariant mass is :math:`\\text{GeV/c}^2`",
2132 Manager::VariableDataType::c_double);
2133
2134 REGISTER_METAVARIABLE("roeP(maskName)", ROE_P,
2135 "Returns momentum of unused tracks and clusters in ROE, can be used with ``use***Frame()`` function. The unit of the momentum is ``GeV/c`` ",
2136 Manager::VariableDataType::c_double);
2137
2138 REGISTER_METAVARIABLE("roePt(maskName)", ROE_Pt,
2139 "Returns transverse component of momentum of unused tracks and clusters in ROE, can be used with ``use***Frame()`` function. The unit of the momentum is ``GeV/c`` ",
2140 Manager::VariableDataType::c_double);
2141
2142 REGISTER_METAVARIABLE("roePx(maskName)", ROE_Px,
2143 "Returns x component of momentum of unused tracks and clusters in ROE, can be used with ``use***Frame()`` function. The unit of the momentum is ``GeV/c`` ",
2144 Manager::VariableDataType::c_double);
2145
2146 REGISTER_METAVARIABLE("roePy(maskName)", ROE_Py,
2147 "Returns y component of momentum of unused tracks and clusters in ROE, can be used with ``use***Frame()`` function. The unit of the momentum is ``GeV/c`` ",
2148 Manager::VariableDataType::c_double);
2149
2150 REGISTER_METAVARIABLE("roePz(maskName)", ROE_Pz,
2151 "Returns z component of momentum of unused tracks and clusters in ROE, can be used with ``use***Frame()`` function. The unit of the momentum is ``GeV/c`` ",
2152 Manager::VariableDataType::c_double);
2153
2154 REGISTER_METAVARIABLE("roePTheta(maskName)", ROE_PTheta,
2155 "Returns theta angle of momentum of unused tracks and clusters in ROE, can be used with ``use***Frame()`` function. The unit of the angle is ``rad`` ",
2156 Manager::VariableDataType::c_double);
2157
2158 REGISTER_METAVARIABLE("roeDeltae(maskName)", ROE_DeltaE,
2159 "Returns energy difference of the related RestOfEvent object with respect to :math:`E_\\mathrm{cms}/2`. The unit of the energy is ``GeV`` ",
2160 Manager::VariableDataType::c_double);
2161
2162 REGISTER_METAVARIABLE("roeMbc(maskName)", ROE_Mbc,
2163 "Returns beam constrained mass of the related RestOfEvent object with respect to :math:`E_\\mathrm{cms}/2`. The unit of the beam constrained mass is :math:`\\text{GeV/c}^2`.",
2164 Manager::VariableDataType::c_double);
2165
2166 REGISTER_METAVARIABLE("weDeltae(maskName, opt)", WE_DeltaE,
2167 "Returns the energy difference of the B meson, corrected with the missing neutrino momentum (reconstructed side + neutrino) with respect to :math:`E_\\mathrm{cms}/2`. The unit of the energy is ``GeV`` ",
2168 Manager::VariableDataType::c_double);
2169
2170 REGISTER_METAVARIABLE("weMbc(maskName, opt)", WE_Mbc,
2171 "Returns beam constrained mass of B meson, corrected with the missing neutrino momentum (reconstructed side + neutrino) with respect to :math:`E_\\mathrm{cms}/2`. The unit of the beam constrained mass is :math:`\\text{GeV/c}^2`.",
2172 Manager::VariableDataType::c_double);
2173
2174 REGISTER_METAVARIABLE("weMissM2(maskName, opt)", WE_MissM2,
2175 "Returns the invariant mass squared of the missing momentum (see :b2:var:`weMissE` possible options). The unit of the invariant mass squared is :math:`[\\text{GeV}/\\text{c}^2]^2`.",
2176 Manager::VariableDataType::c_double);
2177
2178 REGISTER_METAVARIABLE("weMissPTheta(maskName, opt)", WE_MissPTheta,
2179 "Returns the polar angle of the missing momentum (see possible :b2:var:`weMissE` options). The unit of the polar angle is ``rad`` ",
2180 Manager::VariableDataType::c_double);
2181
2182 REGISTER_METAVARIABLE("weMissP(maskName, opt)", WE_MissP,
2183 "Returns the magnitude of the missing momentum (see possible :b2:var:`weMissE` options). The unit of the magnitude of missing momentum is ``GeV/c`` ",
2184 Manager::VariableDataType::c_double);
2185
2186 REGISTER_METAVARIABLE("weMissPx(maskName, opt)", WE_MissPx,
2187 "Returns the x component of the missing momentum (see :b2:var:`weMissE` possible options). The unit of the missing momentum is ``GeV/c`` ",
2188 Manager::VariableDataType::c_double);
2189
2190 REGISTER_METAVARIABLE("weMissPy(maskName, opt)", WE_MissPy,
2191 "Returns the y component of the missing momentum (see :b2:var:`weMissE` possible options). The unit of the missing momentum is ``GeV/c`` ",
2192 Manager::VariableDataType::c_double);
2193
2194 REGISTER_METAVARIABLE("weMissPz(maskName, opt)", WE_MissPz,
2195 "Returns the z component of the missing momentum (see :b2:var:`weMissE` possible options). The unit of the missing momentum is ``GeV/c`` ",
2196 Manager::VariableDataType::c_double);
2197
2198 REGISTER_METAVARIABLE("weMissE(maskName, opt)", WE_MissE,
2199 R"DOC(Returns the energy of the missing momentum. The unit of the Energy is ``GeV`` . Possible options ``opt`` are the following:
2200
2201- ``0``: CMS, use energy and momentum of charged particles and photons
2202- ``1``: CMS, same as ``0``, fix :math:`E_\mathrm{miss} = p_\mathrm{miss}`
2203- ``2``: CMS, same as ``0``, fix :math:`E_\mathrm{roe} = E_\mathrm{cms}/2`
2204- ``3``: CMS, use only energy and momentum of signal side
2205- ``4``: CMS, same as ``3``, update with direction of ROE momentum
2206- ``5``: LAB, use energy and momentum of charged particles and photons from whole event
2207- ``6``: LAB, same as ``5``, fix :math:`E_\mathrm{miss} = p_\mathrm{miss}``
2208- ``7``: CMS, correct pmiss 3-momentum vector with factor alpha so that :math:`d_E = 0`` (used for :math:`M_\mathrm{bc}` calculation).)DOC",
2209 Manager::VariableDataType::c_double);
2210
2211 REGISTER_METAVARIABLE("weXiZ(maskName)", WE_xiZ,
2212 "Returns Xi_z in event (for Bhabha suppression and two-photon scattering). The unit of this variable is ``1/c`` ",
2213 Manager::VariableDataType::c_double);
2214
2215 REGISTER_METAVARIABLE("bssMassDifference(maskName)", bssMassDifference,
2216 "Bs* - Bs mass difference. The unit of the mass is :math:`\\text{GeV/c}^2`.", Manager::VariableDataType::c_double);
2217
2218 REGISTER_METAVARIABLE("weCosThetaEll(maskName)", WE_cosThetaEll, R"DOC(
2219
2220Returns the cosine of the angle between :math:`M` and lepton in :math:`W` rest frame in the decays of the type:
2221:math:`M \to h_1 ... h_n \ell`, where W 4-momentum is given as
2222
2223.. math::
2224 p_W = p_\ell + p_\nu.
2225
2226The neutrino momentum is calculated from ROE taking into account the specified mask, and setting
2227
2228.. math::
2229 E_{\nu} = |p_{miss}|.
2230
2231)DOC", Manager::VariableDataType::c_double);
2232
2233 REGISTER_METAVARIABLE("weQ2lnuSimple(maskName,option)", WE_q2lnuSimple,
2234 "Returns the momentum transfer squared, :math:`q^2`, calculated in CMS as :math:`q^2 = (p_l + p_\\nu)^2`, \n"
2235 "where :math:`B \\to H_1\\dots H_n \\ell \\nu_\\ell`. Lepton is assumed to be the last reconstructed daughter. \n"
2236 "By default, option is set to ``1`` (see :b2:var:`weMissE`). Unless you know what you are doing, keep this default value. The unit of the momentum transfer squared is :math:`[\\text{GeV}/\\text{c}]^2`.", Manager::VariableDataType::c_double);
2237
2238 REGISTER_METAVARIABLE("weQ2lnu(maskName,option)", WE_q2lnu,
2239 "Returns the momentum transfer squared, :math:`q^2`, calculated in CMS as :math:`q^2 = (p_l + p_\\nu)^2`, \n"
2240 "where :math:`B \\to H_1\\dots H_n \\ell \\nu_\\ell`. Lepton is assumed to be the last reconstructed daughter. \n"
2241 "This calculation uses constraints from dE = 0 and Mbc = Mb to correct the neutrino direction. \n"
2242 "By default, option is set to ``7`` (see :b2:var:`weMissE`). Unless you know what you are doing, keep this default value. The unit of the momentum transfer squared is :math:`[\\text{GeV}/\\text{c}]^2`.", Manager::VariableDataType::c_double);
2243
2244 REGISTER_METAVARIABLE("weMissM2OverMissE(maskName)", WE_MissM2OverMissE,
2245 "Returns missing mass squared over missing energy. The unit of the missing mass squared is :math:`\\text{GeV/c}^4`.", Manager::VariableDataType::c_double);
2246
2247 REGISTER_METAVARIABLE("passesROEMask(maskName)", passesROEMask,
2248 "Returns boolean value if a particle passes a certain mask or not. Only to be used in for_each path, otherwise returns quiet NaN.", Manager::VariableDataType::c_double);
2249
2250 REGISTER_VARIABLE("printROE", printROE,
2251 "For debugging, prints indices of all particles in the ROE and all masks. Returns 0.");
2252
2253 REGISTER_VARIABLE("hasCorrectROECombination", hasCorrectROECombination,
2254 "Returns 1 if there is correct combination of daughter particles between the particle that is the basis of the ROE and the particle loaded from the ROE. "
2255 "Returns 0 if there is not correct combination. "
2256 "If there is no daughter particle loaded from the ROE, returns quiet NaN.");
2257
2258 REGISTER_METAVARIABLE("pi0Prob(mode)", pi0Prob,
2259 "Returns pi0 probability, where mode is used to specify the selection criteria for soft photon. \n"
2260 "The following strings are available. \n\n"
2261 "- ``standard``: loose energy cut and no clusterNHits cut are applied to soft photon \n"
2262 "- ``tight``: tight energy cut and no clusterNHits cut are applied to soft photon \n"
2263 "- ``cluster``: loose energy cut and clusterNHits cut are applied to soft photon \n"
2264 "- ``both``: tight energy cut and clusterNHits cut are applied to soft photon \n\n"
2265 "You can find more details in `writePi0EtaVeto` function in modularAnalysis.py.", Manager::VariableDataType::c_double);
2266
2267 REGISTER_METAVARIABLE("etaProb(mode)", etaProb,
2268 "Returns eta probability, where mode is used to specify the selection criteria for soft photon. \n"
2269 "The following strings are available. \n\n"
2270 "- ``standard``: loose energy cut and no clusterNHits cut are applied to soft photon \n"
2271 "- ``tight``: tight energy cut and no clusterNHits cut are applied to soft photon \n"
2272 "- ``cluster``: loose energy cut and clusterNHits cut are applied to soft photon \n"
2273 "- ``both``: tight energy cut and clusterNHits cut are applied to soft photon \n\n"
2274 "You can find more details in `writePi0EtaVeto` function in modularAnalysis.py.", Manager::VariableDataType::c_double);
2275
2276 }
2278}
R E
internal precision of FFTW codelets
void SetMag(DataType mag)
Set magnitude keeping theta and phi constant.
Definition: B2Vector3.h:182
B2Vector3< DataType > Cross(const B2Vector3< DataType > &p) const
Cross product.
Definition: B2Vector3.h:296
DataType Mag2() const
The magnitude squared (rho^2 in spherical coordinate system).
Definition: B2Vector3.h:157
DataType Angle(const B2Vector3< DataType > &q) const
The angle w.r.t.
Definition: B2Vector3.h:302
static const ParticleType neutron
neutron particle
Definition: Const.h:675
static const ParticleType pi0
neutral pion particle
Definition: Const.h:674
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
static const ChargedStable proton
proton particle
Definition: Const.h:663
static const ParticleType Kshort
K^0_S particle.
Definition: Const.h:677
static const double doubleNaN
quiet_NaN
Definition: Const.h:703
static const ChargedStable kaon
charged kaon particle
Definition: Const.h:662
static const ParticleType photon
photon particle
Definition: Const.h:673
static const ChargedStable electron
electron particle
Definition: Const.h:659
@ c_nPhotons
CR is split into n photons (N1)
@ c_PrimaryParticle
bit 0: Particle is primary particle.
Definition: MCParticle.h:47
static const ReferenceFrame & GetCurrent()
Get current rest frame.
static constexpr const char * c_defaultMaskName
Default mask name.
Definition: RestOfEvent.h:60
std::function< VarVariant(const Particle *)> FunctionPtr
functions stored take a const Particle* and return VarVariant.
Definition: Manager.h:113
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
B2Vector3< double > B2Vector3D
typedef for common usage with double
Definition: B2Vector3.h:516
double sqrt(double a)
sqrt for double
Definition: beamHelpers.h:28
Abstract base class for different kinds of events.
STL namespace.