Belle II Software development
PIDVariables.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/PIDVariables.h>
11
12#include <analysis/dataobjects/Particle.h>
13#include <analysis/utility/ReferenceFrame.h>
14#include <mdst/dataobjects/PIDLikelihood.h>
15#include <mdst/dataobjects/TrackFitResult.h>
16
17// framework aux
18#include <framework/logging/Logger.h>
19#include <framework/utilities/Conversion.h>
20#include <framework/gearbox/Const.h>
21
22// database
23//#include <framework/database/DBObjPtr.h>
24#include <analysis/dbobjects/PIDCalibrationWeight.h>
25#include <analysis/utility/PIDCalibrationWeightUtil.h>
26#include <analysis/utility/PIDNeuralNetwork.h>
27
28#include <boost/algorithm/string.hpp>
29
30#include <iostream>
31#include <cmath>
32
33using namespace std;
34
35namespace Belle2 {
40 namespace Variable {
41
42 //*************
43 // Utilities
44 //*************
45
46 // converts Belle numbering scheme for charged final state particles
47 // to Belle II ChargedStable
48 Const::ChargedStable hypothesisConversion(const int hypothesis)
49 {
50 switch (hypothesis) {
51 case 0:
52 return Const::electron;
53 case 1:
54 return Const::muon;
55 case 2:
56 return Const::pion;
57 case 3:
58 return Const::kaon;
59 case 4:
60 return Const::proton;
61 }
62
63 return Const::pion;
64 }
65
66
67 Const::PIDDetectorSet parseDetectors(const std::vector<std::string>& arguments)
68 {
70 for (std::string val : arguments) {
71 boost::to_lower(val);
72 if (val == "all") return Const::PIDDetectors::set();
73 else if (val == "svd") result += Const::SVD;
74 else if (val == "cdc") result += Const::CDC;
75 else if (val == "top") result += Const::TOP;
76 else if (val == "arich") result += Const::ARICH;
77 else if (val == "ecl") result += Const::ECL;
78 else if (val == "klm") result += Const::KLM;
79 else B2ERROR("Unknown detector component: " << val);
80 }
81 return result;
82 }
83
84 // Specialisation of valid detectors parser for ChargedBDT.
85 Const::PIDDetectorSet parseDetectorsChargedBDT(const std::vector<std::string>& arguments)
86 {
88 for (std::string val : arguments) {
89 boost::to_lower(val);
90 if (val == "all") return Const::PIDDetectors::set();
91 else if (val == "ecl") result += Const::ECL;
92 else B2ERROR("Invalid detector component: " << val << " for charged BDT.");
93 }
94 return result;
95 }
96
97 //*************
98 // Belle II
99 //*************
100
101 // a "smart" variable:
102 // finds the global probability based on the PDG code of the input particle
103 double particleID(const Particle* p)
104 {
105 int pdg = abs(p->getPDGCode());
106 if (pdg == Const::electron.getPDGCode()) return electronID(p);
107 else if (pdg == Const::muon.getPDGCode()) return muonID(p);
108 else if (pdg == Const::pion.getPDGCode()) return pionID(p);
109 else if (pdg == Const::kaon.getPDGCode()) return kaonID(p);
110 else if (pdg == Const::proton.getPDGCode()) return protonID(p);
111 else if (pdg == Const::deuteron.getPDGCode()) return deuteronID(p);
112 else return Const::doubleNaN;
113 }
114
115 bool isPIDAvailable(const Particle* part)
116 {
117 const auto* pid = part->getPIDLikelihood();
118 if (not pid) return false;
119 return pid->isAvailable();
120 }
121
122 Manager::FunctionPtr isPIDAvailableFrom(const std::vector<std::string>& arguments)
123 {
124 Const::PIDDetectorSet detectorSet = parseDetectors(arguments);
125 auto func = [detectorSet](const Particle * part) -> bool {
126 const auto* pid = part->getPIDLikelihood();
127 if (not pid) return false;
128 return pid->isAvailable(detectorSet);
129 };
130 return func;
131 }
132
133 Manager::FunctionPtr pidLogLikelihoodValueExpert(const std::vector<std::string>& arguments)
134 {
135 if (arguments.size() < 2) {
136 B2ERROR("Need at least two arguments to pidLogLikelihoodValueExpert");
137 return nullptr;
138 }
139 int pdgCode;
140 try {
141 pdgCode = Belle2::convertString<int>(arguments[0]);
142 } catch (std::invalid_argument& e) {
143 B2ERROR("First argument of pidLogLikelihoodValueExpert must be a PDG code");
144 return nullptr;
145 }
146 std::vector<std::string> detectors(arguments.begin() + 1, arguments.end());
147
148 Const::PIDDetectorSet detectorSet = parseDetectors(detectors);
149 auto hypType = Const::ChargedStable(abs(pdgCode));
150
151 auto func = [hypType, detectorSet](const Particle * part) -> double {
152 const PIDLikelihood* pid = part->getPIDLikelihood();
153 if (!pid)
154 return Const::doubleNaN;
155 // No information from any subdetector in the list
156 if (not pid->isAvailable(detectorSet)) return Const::doubleNaN;
157
158 return pid->getLogL(hypType, detectorSet);
159 };
160 return func;
161 }
162
163
164 Manager::FunctionPtr pidDeltaLogLikelihoodValueExpert(const std::vector<std::string>& arguments)
165 {
166 if (arguments.size() < 3) {
167 B2ERROR("Need at least three arguments to pidDeltaLogLikelihoodValueExpert");
168 return nullptr;
169 }
170 int pdgCodeHyp, pdgCodeTest;
171 try {
172 pdgCodeHyp = Belle2::convertString<int>(arguments[0]);
173 } catch (std::invalid_argument& e) {
174 B2ERROR("First argument of pidDeltaLogLikelihoodValueExpert must be a PDG code");
175 return nullptr;
176 }
177 try {
178 pdgCodeTest = Belle2::convertString<int>(arguments[1]);
179 } catch (std::invalid_argument& e) {
180 B2ERROR("Second argument of pidDeltaLogLikelihoodValueExpert must be a PDG code");
181 return nullptr;
182 }
183
184 std::vector<std::string> detectors(arguments.begin() + 2, arguments.end());
185 Const::PIDDetectorSet detectorSet = parseDetectors(detectors);
186 auto hypType = Const::ChargedStable(abs(pdgCodeHyp));
187 auto testType = Const::ChargedStable(abs(pdgCodeTest));
188
189 auto func = [hypType, testType, detectorSet](const Particle * part) -> double {
190 const PIDLikelihood* pid = part->getPIDLikelihood();
191 if (!pid) return Const::doubleNaN;
192 // No information from any subdetector in the list
193 if (not pid->isAvailable(detectorSet)) return Const::doubleNaN;
194
195 return pid->getDeltaLogL(hypType, testType, detectorSet);
196 };
197 return func;
198 }
199
200
201 Manager::FunctionPtr pidPairProbabilityExpert(const std::vector<std::string>& arguments)
202 {
203 if (arguments.size() < 3) {
204 B2ERROR("Need at least three arguments to pidPairProbabilityExpert");
205 return nullptr;
206 }
207 int pdgCodeHyp = 0, pdgCodeTest = 0;
208 try {
209 pdgCodeHyp = Belle2::convertString<int>(arguments[0]);
210 } catch (std::invalid_argument& e) {
211 B2ERROR("First argument of pidPairProbabilityExpert must be PDG code");
212 return nullptr;
213 }
214 try {
215 pdgCodeTest = Belle2::convertString<int>(arguments[1]);
216 } catch (std::invalid_argument& e) {
217 B2ERROR("Second argument of pidPairProbabilityExpert must be PDG code");
218 return nullptr;
219 }
220
221 std::vector<std::string> detectors(arguments.begin() + 2, arguments.end());
222
223 Const::PIDDetectorSet detectorSet = parseDetectors(detectors);
224 auto hypType = Const::ChargedStable(abs(pdgCodeHyp));
225 auto testType = Const::ChargedStable(abs(pdgCodeTest));
226 auto func = [hypType, testType, detectorSet](const Particle * part) -> double {
227 const PIDLikelihood* pid = part->getPIDLikelihood();
228 if (!pid) return Const::doubleNaN;
229 // No information from any subdetector in the list
230 if (not pid->isAvailable(detectorSet)) return Const::doubleNaN;
231
232 return pid->getProbability(hypType, testType, detectorSet);
233 };
234 return func;
235 }
236
237
238 Manager::FunctionPtr pidProbabilityExpert(const std::vector<std::string>& arguments)
239 {
240 if (arguments.size() < 2) {
241 B2ERROR("Need at least two arguments for pidProbabilityExpert");
242 return nullptr;
243 }
244 int pdgCodeHyp = 0;
245 try {
246 pdgCodeHyp = Belle2::convertString<int>(arguments[0]);
247 } catch (std::invalid_argument& e) {
248 B2ERROR("First argument of pidProbabilityExpert must be PDG code");
249 return nullptr;
250 }
251
252 std::vector<std::string> detectors(arguments.begin() + 1, arguments.end());
253 Const::PIDDetectorSet detectorSet = parseDetectors(detectors);
254 auto hypType = Const::ChargedStable(abs(pdgCodeHyp));
255
256 // Placeholder for the priors
257 const unsigned int n = Const::ChargedStable::c_SetSize;
258 double frac[n];
259 for (double& i : frac) i = 1.0; // flat priors
260
261 auto func = [hypType, frac, detectorSet](const Particle * part) -> double {
262 const PIDLikelihood* pid = part->getPIDLikelihood();
263 if (!pid) return Const::doubleNaN;
264 // No information from any subdetector in the list
265 if (not pid->isAvailable(detectorSet)) return Const::doubleNaN;
266
267 return pid->getProbability(hypType, frac, detectorSet);
268 };
269 return func;
270 }
271
272
273 Manager::FunctionPtr pidLogarithmicProbabilityExpert(const std::vector<std::string>& arguments)
274 {
275 if (arguments.size() < 2) {
276 B2ERROR("Need at least two arguments for pidLogarithmicProbabilityExpert");
277 return nullptr;
278 }
279 int pdgCodeHyp = 0;
280 try {
281 pdgCodeHyp = Belle2::convertString<int>(arguments[0]);
282 } catch (std::invalid_argument& e) {
283 B2ERROR("First argument of pidLogarithmicProbabilityExpert must be PDG code");
284 return nullptr;
285 }
286
287 std::vector<std::string> detectors(arguments.begin() + 1, arguments.end());
288 Const::PIDDetectorSet detectorSet = parseDetectors(detectors);
289 auto hypType = Const::ChargedStable(abs(pdgCodeHyp));
290
291 // Placeholder for the priors
292 const unsigned int n = Const::ChargedStable::c_SetSize;
293 double frac[n];
294 for (double& i : frac) i = 1.0; // flat priors
295
296 auto func = [hypType, frac, detectorSet](const Particle * part) -> double {
297 const auto* pid = part->getPIDLikelihood();
298 if (not pid) return Const::doubleNaN;
299 if (not pid->isAvailable(detectorSet)) return Const::doubleNaN;
300 return pid->getLogarithmicProbability(hypType, frac, detectorSet);
301 };
302 return func;
303 }
304
305
306 Manager::FunctionPtr pidMissingProbabilityExpert(const std::vector<std::string>& arguments)
307 {
308 if (arguments.size() < 1) {
309 B2ERROR("Need at least one argument to pidMissingProbabilityExpert");
310 return nullptr;
311 }
312
313 std::vector<std::string> detectors(arguments.begin(), arguments.end());
314 Const::PIDDetectorSet detectorSet = parseDetectors(detectors);
315
316 auto func = [detectorSet](const Particle * part) -> double {
317 const PIDLikelihood* pid = part->getPIDLikelihood();
318 if (!pid) return Const::doubleNaN;
319 if (not pid->areAllAvailable(detectorSet)) return 1; // to keep the same behaviour as in previous releases
320 else return 0;
321 };
322 return func;
323 }
324
325 Manager::FunctionPtr pidWeightedLogLikelihoodValueExpert(const std::vector<std::string>& arguments)
326 {
327 if (arguments.size() < 3) {
328 B2ERROR("Need at least three arguments to pidWeightedLogLikelihoodValueExpert");
329 return nullptr;
330 }
331 std::string matrixName = arguments[0];
332
333 int pdgCode;
334 try {
335 pdgCode = Belle2::convertString<int>(arguments[1]);
336 } catch (std::invalid_argument& e) {
337 B2ERROR("Second argument of pidWeightedLogLikelihoodValueExpert must be a PDG code");
338 return nullptr;
339 }
340 std::vector<std::string> detectors(arguments.begin() + 2, arguments.end());
341 Const::PIDDetectorSet detectorSet = parseDetectors(detectors);
342 auto hypType = Const::ChargedStable(abs(pdgCode));
343
344 auto func = [hypType, detectorSet, matrixName](const Particle * part) -> double {
345 PIDCalibrationWeightUtil weightMatrix(matrixName);
346 const PIDLikelihood* pid = part->getPIDLikelihood();
347 if (!pid)
348 return Const::doubleNaN;
349 // No information from any subdetector in the list
350 if (not pid->isAvailable(detectorSet)) return Const::doubleNaN;
351
352 const auto& frame = ReferenceFrame::GetCurrent();
353 auto mom = frame.getMomentum(part);
354 auto p = mom.P();
355 auto theta = mom.Theta();
356
357 double LogL = 0;
358 for (const Const::EDetector& detector : Const::PIDDetectorSet::set())
359 {
360 if (detectorSet.contains(detector))
361 LogL += pid->getLogL(hypType, detector) * weightMatrix.getWeight(hypType.getPDGCode(), detector, p, theta);
362 }
363 return LogL;
364 };
365
366 return func;
367 }
368
369 Manager::FunctionPtr pidWeightedProbabilityExpert(const std::vector<std::string>& arguments)
370 {
371 if (arguments.size() < 3) {
372 B2ERROR("Need at least three arguments for pidWeightedProbabilityExpert");
373 return nullptr;
374 }
375 std::string matrixName = arguments[0];
376
377 int pdgCodeHyp = 0;
378 try {
379 pdgCodeHyp = Belle2::convertString<int>(arguments[1]);
380 } catch (std::invalid_argument& e) {
381 B2ERROR("Second argument of pidWeightedProbabilityExpert must be PDG code");
382 return nullptr;
383 }
384
385 std::vector<std::string> detectors(arguments.begin() + 2, arguments.end());
386 Const::PIDDetectorSet detectorSet = parseDetectors(detectors);
387 auto hypType = Const::ChargedStable(abs(pdgCodeHyp));
388
389 auto func = [hypType, detectorSet, matrixName](const Particle * part) -> double {
390 PIDCalibrationWeightUtil weightMatrix(matrixName);
391 const PIDLikelihood* pid = part->getPIDLikelihood();
392 if (!pid) return Const::doubleNaN;
393 // No information from any subdetector in the list
394 if (not pid->isAvailable(detectorSet)) return Const::doubleNaN;
395
396 const auto& frame = ReferenceFrame::GetCurrent();
397 auto mom = frame.getMomentum(part);
398 auto p = mom.P();
399 auto theta = mom.Theta();
400
401 std::vector<double> LogL(Const::ChargedStable::c_SetSize);
402 double LogL_max = 0;
403 bool hasMax = false;
404 for (const auto& pdgIter : Const::chargedStableSet)
405 {
406 const int index_pdg = pdgIter.getIndex();
407
408 LogL[index_pdg] = 0;
409 for (const Const::EDetector& detector : Const::PIDDetectorSet::set()) {
410 if (detectorSet.contains(detector))
411 LogL[index_pdg] += pid->getLogL(pdgIter, detector) * weightMatrix.getWeight(pdgIter.getPDGCode(), detector, p, theta);
412 }
413
414 if (!hasMax || (LogL[index_pdg] > LogL_max)) {
415 LogL_max = LogL[index_pdg];
416 hasMax = true;
417 }
418 }
419
420 double norm = 0;
421 for (auto LogL_i : LogL)
422 norm += exp(LogL_i - LogL_max);
423
424 if (norm > 0)
425 return exp(LogL[hypType.getIndex()] - LogL_max) / norm;
426 else
427 return -1;
428 };
429 return func;
430 }
431
432
433 Manager::FunctionPtr pidNeuralNetworkValueExpert(const std::vector<std::string>& arguments)
434 {
435 if (arguments.size() == 0) {
436 B2ERROR("Need pdg code for pidNeuralNetworkValueExpert");
437 return nullptr;
438 }
439 if (arguments.size() > 2) {
440 B2ERROR("pidNeuralNetworkValueExpert expects at most two arguments, i.e. the pdg code and the pidNeuralNetworkName");
441 return nullptr;
442 }
443 int pdgCode;
444 try {
445 pdgCode = abs(Belle2::convertString<int>(arguments[0]));
446 } catch (std::invalid_argument& e) {
447 B2ERROR("First argument of pidNeuralNetworkValueExpert must be a PDG code");
448 return nullptr;
449 }
450
451 std::shared_ptr<PIDNeuralNetwork> neuralNetworkPtr;
452 if (arguments.size() == 2) {
453 std::string parameterSetName = arguments[1];
454 neuralNetworkPtr = std::make_shared<PIDNeuralNetwork>(parameterSetName);
455 } else {
456 neuralNetworkPtr = std::make_shared<PIDNeuralNetwork>();
457 }
458
459 neuralNetworkPtr->hasPdgCode(pdgCode, true); // raise exception if pdg code is not predicted
460
461 auto func = [neuralNetworkPtr, pdgCode](const Particle * part) -> double {
462 const auto& extraInfoName = neuralNetworkPtr->getExtraInfoName(pdgCode);
463 if (part->hasExtraInfo(extraInfoName))
464 return part->getExtraInfo(extraInfoName);
465
466 const PIDLikelihood* pid = part->getPIDLikelihood();
467 if (not pid) return Const::doubleNaN;
468 if (not pid->isAvailable()) return Const::doubleNaN;
469
470 // collect only those inputs needed for the neural network in the correct order
471 std::vector<float> inputsNN;
472 inputsNN.reserve(neuralNetworkPtr->getInputSize());
473 for (const auto& inputName : neuralNetworkPtr->getInputBasf2Names())
474 inputsNN.push_back(std::get<double>(Variable::Manager::Instance().getVariable(inputName)->function(part)));
475
476 const auto probabilities = neuralNetworkPtr->predict(inputsNN);
477
478 // store all probabilities for all hypotheses in extraInfo
479 for (const auto element : probabilities)
480 {
481 const auto [pdgCodeElement, probability] = element;
482 const_cast<Particle*>(part)->addExtraInfo(neuralNetworkPtr->getExtraInfoName(pdgCodeElement), probability);
483 }
484
485 return probabilities.at(pdgCode);
486 };
487
488 return func;
489 }
490
491
492 Manager::FunctionPtr pidWeightedPairProbabilityExpert(const std::vector<std::string>& arguments)
493 {
494 if (arguments.size() < 4) {
495 B2ERROR("Need at least four arguments to pidWeightedPairProbabilityExpert");
496 return nullptr;
497 }
498 std::string matrixName = arguments[0];
499
500 int pdgCodeHyp = 0, pdgCodeTest = 0;
501 try {
502 pdgCodeHyp = Belle2::convertString<int>(arguments[1]);
503 } catch (std::invalid_argument& e) {
504 B2ERROR("Second argument of pidWeightedPairProbabilityExpert must be PDG code");
505 return nullptr;
506 }
507 try {
508 pdgCodeTest = Belle2::convertString<int>(arguments[2]);
509 } catch (std::invalid_argument& e) {
510 B2ERROR("Third argument of pidWeightedPairProbabilityExpert must be PDG code");
511 return nullptr;
512 }
513
514 std::vector<std::string> detectors(arguments.begin() + 3, arguments.end());
515
516 Const::PIDDetectorSet detectorSet = parseDetectors(detectors);
517 auto hypType = Const::ChargedStable(abs(pdgCodeHyp));
518 auto testType = Const::ChargedStable(abs(pdgCodeTest));
519
520 auto func = [hypType, testType, detectorSet, matrixName](const Particle * part) -> double {
521 PIDCalibrationWeightUtil weightMatrix(matrixName);
522
523 const PIDLikelihood* pid = part->getPIDLikelihood();
524 if (!pid) return Const::doubleNaN;
525 // No information from any subdetector in the list
526 if (not pid->isAvailable(detectorSet)) return Const::doubleNaN;
527
528 const auto& frame = ReferenceFrame::GetCurrent();
529 auto mom = frame.getMomentum(part);
530 auto p = mom.P();
531 auto theta = mom.Theta();
532
533 double LogL_hypType(0), LogL_testType(0);
534 for (const Const::EDetector& detector : Const::PIDDetectorSet::set())
535 {
536 if (detectorSet.contains(detector)) {
537 LogL_hypType += pid->getLogL(hypType, detector) * weightMatrix.getWeight(hypType.getPDGCode(), detector, p, theta);
538 LogL_testType += pid->getLogL(testType, detector) * weightMatrix.getWeight(testType.getPDGCode(), detector, p, theta);
539 }
540 }
541
542 double deltaLogL = LogL_testType - LogL_hypType;
543 double res;
544 if (deltaLogL < 0)
545 {
546 double eLogL = exp(deltaLogL);
547 res = 1. / (1. + eLogL);
548 } else
549 {
550 double eLogL = exp(-deltaLogL);
551 res = eLogL / (1.0 + eLogL);
552 }
553
554 if (std::isfinite(res))
555 return res;
556
557 return 0;
558 };
559 return func;
560 }
561
562 double electronID(const Particle* part)
563 {
564 static Manager::FunctionPtr pidFunction =
565 pidProbabilityExpert({"11", "ALL"});
566 return std::get<double>(pidFunction(part));
567 }
568
569 double muonID(const Particle* part)
570 {
571 static Manager::FunctionPtr pidFunction =
572 pidProbabilityExpert({"13", "ALL"});
573 return std::get<double>(pidFunction(part));
574 }
575
576 double pionID(const Particle* part)
577 {
578 static Manager::FunctionPtr pidFunction =
579 pidProbabilityExpert({"211", "ALL"});
580 return std::get<double>(pidFunction(part));
581 }
582
583 double kaonID(const Particle* part)
584 {
585 static Manager::FunctionPtr pidFunction =
586 pidProbabilityExpert({"321", "ALL"});
587 return std::get<double>(pidFunction(part));
588 }
589
590 double protonID(const Particle* part)
591 {
592 static Manager::FunctionPtr pidFunction =
593 pidProbabilityExpert({"2212", "ALL"});
594 return std::get<double>(pidFunction(part));
595 }
596
597 double deuteronID(const Particle* part)
598 {
599 static Manager::FunctionPtr pidFunction =
600 pidProbabilityExpert({"1000010020", "ALL"});
601 return std::get<double>(pidFunction(part));
602 }
603
604 double binaryPID(const Particle* part, const std::vector<double>& arguments)
605 {
606 if (arguments.size() != 2) {
607 B2ERROR("The variable binaryPID needs exactly two arguments: the PDG codes of two hypotheses.");
608 return Const::doubleNaN;;
609 }
610 int pdgCodeHyp = std::abs(int(std::lround(arguments[0])));
611 int pdgCodeTest = std::abs(int(std::lround(arguments[1])));
612 return std::get<double>(Manager::Instance().getVariable("pidPairProbabilityExpert(" + std::to_string(
613 pdgCodeHyp) + ", " + std::to_string(
614 pdgCodeTest) + ", ALL)")->function(part));
615 }
616
617 double electronID_noSVD(const Particle* part)
618 {
619 // Excluding SVD for electron ID. This variable is temporary. BII-8760
620 static Manager::FunctionPtr pidFunction =
621 pidProbabilityExpert({"11", "CDC", "TOP", "ARICH", "ECL", "KLM"});
622 return std::get<double>(pidFunction(part));
623 }
624
625 double muonID_noSVD(const Particle* part)
626 {
627 // Excluding SVD for muon ID. This variable is temporary. BII-8760
628 static Manager::FunctionPtr pidFunction =
629 pidProbabilityExpert({"13", "CDC", "TOP", "ARICH", "ECL", "KLM"});
630 return std::get<double>(pidFunction(part));
631 }
632
633 double pionID_noSVD(const Particle* part)
634 {
635 // Excluding SVD for pion ID. This variable is temporary. BII-8760
636 static Manager::FunctionPtr pidFunction =
637 pidProbabilityExpert({"211", "CDC", "TOP", "ARICH", "ECL", "KLM"});
638 return std::get<double>(pidFunction(part));
639 }
640
641 double kaonID_noSVD(const Particle* part)
642 {
643 // Excluding SVD for kaon ID. This variable is temporary. BII-8760
644 static Manager::FunctionPtr pidFunction =
645 pidProbabilityExpert({"321", "CDC", "TOP", "ARICH", "ECL", "KLM"});
646 return std::get<double>(pidFunction(part));
647 }
648
649 double protonID_noSVD(const Particle* part)
650 {
651 // Excluding SVD for proton ID. This variable is temporary. BII-8760
652 static Manager::FunctionPtr pidFunction =
653 pidProbabilityExpert({"2212", "CDC", "TOP", "ARICH", "ECL", "KLM"});
654 return std::get<double>(pidFunction(part));
655 }
656
657 double deuteronID_noSVD(const Particle* part)
658 {
659 // Excluding SVD for deuteron ID. This variable is temporary. BII-8760
660 static Manager::FunctionPtr pidFunction =
661 pidProbabilityExpert({"1000010020", "CDC", "TOP", "ARICH", "ECL", "KLM"});
662 return std::get<double>(pidFunction(part));
663 }
664
665 double binaryPID_noSVD(const Particle* part, const std::vector<double>& arguments)
666 {
667 // Excluding SVD for binary ID. This variable is temporary. BII-8760
668 if (arguments.size() != 2) {
669 B2ERROR("The variable binaryPID_noSVD needs exactly two arguments: the PDG codes of two hypotheses.");
670 return Const::doubleNaN;;
671 }
672 int pdgCodeHyp = std::abs(int(std::lround(arguments[0])));
673 int pdgCodeTest = std::abs(int(std::lround(arguments[1])));
674 return std::get<double>(Manager::Instance().getVariable("pidPairProbabilityExpert(" + std::to_string(
675 pdgCodeHyp) + ", " + std::to_string(
676 pdgCodeTest) + ", CDC, TOP, ARICH, ECL, KLM)")->function(part));
677 }
678
679 double electronID_noTOP(const Particle* part)
680 {
681 // Excluding TOP for electron ID. This variable is temporary. BII-8444
682 static Manager::FunctionPtr pidFunction =
683 pidProbabilityExpert({"11", "SVD", "CDC", "ARICH", "ECL", "KLM"});
684 return std::get<double>(pidFunction(part));
685 }
686
687 double binaryElectronID_noTOP(const Particle* part, const std::vector<double>& arguments)
688 {
689 // Excluding TOP for electron ID. This is temporary. BII-8444
690 if (arguments.size() != 1) {
691 B2ERROR("The variable binaryElectronID_noTOP needs exactly one argument: the PDG code of the test hypothesis.");
692 return Const::doubleNaN;;
693 }
694
695 int pdgCodeHyp = Const::electron.getPDGCode();
696 int pdgCodeTest = std::abs(int(std::lround(arguments[0])));
697
698 const auto var = "pidPairProbabilityExpert(" + std::to_string(pdgCodeHyp) + ", " +
699 std::to_string(pdgCodeTest) + ", SVD, CDC, ARICH, ECL, KLM)";
700
701 return std::get<double>(Manager::Instance().getVariable(var)->function(part));
702 }
703
704 double electronID_noSVD_noTOP(const Particle* part)
705 {
706 // Excluding SVD and TOP for electron ID. This variable is temporary. BII-8444, BII-8760.
707 static Manager::FunctionPtr pidFunction =
708 pidProbabilityExpert({"11", "CDC", "ARICH", "ECL", "KLM"});
709 return std::get<double>(pidFunction(part));
710 }
711
712 double binaryElectronID_noSVD_noTOP(const Particle* part, const std::vector<double>& arguments)
713 {
714 // Excluding SVD and TOP for electron ID. This is temporary. BII-8444, BII-8760.
715 if (arguments.size() != 1) {
716 B2ERROR("The variable binaryElectronID_noSVD_noTOP needs exactly one argument: the PDG code of the test hypothesis.");
717 return Const::doubleNaN;;
718 }
719
720 int pdgCodeHyp = Const::electron.getPDGCode();
721 int pdgCodeTest = std::abs(int(std::lround(arguments[0])));
722
723 const auto var = "pidPairProbabilityExpert(" + std::to_string(pdgCodeHyp) + ", " +
724 std::to_string(pdgCodeTest) + ", CDC, ARICH, ECL, KLM)";
725
726 return std::get<double>(Manager::Instance().getVariable(var)->function(part));
727 }
728
729
730 double pionID_noARICHwoECL(const Particle* part)
731 {
732 // remove arich if no ecl cluster + identified as kaon in arich
733 const ECLCluster* cluster = part->getECLCluster();
734 if (!cluster) {
735 const PIDLikelihood* pid = part->getPIDLikelihood();
736 if (!pid) return Const::doubleNaN;
737 if (pid->getDeltaLogL(Const::kaon, Const::pion, Const::ARICH) > 0) {
738 static Manager::FunctionPtr pidFunction =
739 pidProbabilityExpert({"211", "SVD", "CDC", "TOP", "ECL", "KLM"});
740 return std::get<double>(pidFunction(part));
741 }
742 }
743 return pionID(part);
744 }
745
746
747 double kaonID_noARICHwoECL(const Particle* part)
748 {
749 // remove arich if no ecl cluster + identified as kaon in arich
750 const ECLCluster* cluster = part->getECLCluster();
751 if (!cluster) {
752 const PIDLikelihood* pid = part->getPIDLikelihood();
753 if (!pid) return Const::doubleNaN;
754 if (pid->getDeltaLogL(Const::kaon, Const::pion, Const::ARICH) > 0) {
755 static Manager::FunctionPtr pidFunction =
756 pidProbabilityExpert({"321", "SVD", "CDC", "TOP", "ECL", "KLM"});
757 return std::get<double>(pidFunction(part));
758 }
759 }
760 return kaonID(part);
761 }
762
763
764 double binaryPID_noARICHwoECL(const Particle* part, const std::vector<double>& arguments)
765 {
766 // Excluding ARICH for tracks without ECL cluster and identified as heavier of the two hypotheses from binary ID.
767 if (arguments.size() != 2) {
768 B2ERROR("The variable binaryPID_noARICHwoECL needs exactly two arguments: the PDG codes of two hypotheses.");
769 return Const::doubleNaN;;
770 }
771 int pdgCodeHyp = std::abs(int(std::lround(arguments[0])));
772 int pdgCodeTest = std::abs(int(std::lround(arguments[1])));
773 auto hypType = Const::ChargedStable(abs(pdgCodeHyp));
774 auto testType = Const::ChargedStable(abs(pdgCodeTest));
775
776 const ECLCluster* cluster = part->getECLCluster();
777 if (!cluster) {
778 const PIDLikelihood* pid = part->getPIDLikelihood();
779 if (!pid) return Const::doubleNaN;
780 double lkhdiff = pid->getDeltaLogL(hypType, testType, Const::ARICH);
781 if ((lkhdiff > 0 && pdgCodeHyp > pdgCodeTest) || (lkhdiff < 0 && pdgCodeHyp < pdgCodeTest)) {
782 std::string varName = "pidPairProbabilityExpert(" +
783 std::to_string(pdgCodeHyp) + ", " + std::to_string(pdgCodeTest) + ", SVD, CDC, TOP, ECL, KLM)";
784 return std::get<double>(Manager::Instance().getVariable(varName)->function(part));
785 }
786 }
787
788 return binaryPID(part, arguments);
789
790 }
791
792
793
794 double antineutronID(const Particle* particle)
795 {
796 if (particle->hasExtraInfo("nbarID")) {
797 return particle->getExtraInfo("nbarID");
798 } else {
799 if (particle->getPDGCode() == -Const::neutron.getPDGCode()) {
800 B2WARNING("The extraInfo nbarID is not registered! \n"
801 "Please use function getNbarIDMVA in modularAnalysis.");
802 }
803 return Const::doubleNaN;
804 }
805 }
806
807 Manager::FunctionPtr pidChargedBDTScore(const std::vector<std::string>& arguments)
808 {
809 if (arguments.size() != 2) {
810 B2ERROR("Need exactly two arguments for pidChargedBDTScore: pdgCodeHyp, detector");
811 return nullptr;
812 }
813
814 int hypPdgId;
815 try {
816 hypPdgId = Belle2::convertString<int>(arguments.at(0));
817 } catch (std::invalid_argument& e) {
818 B2ERROR("First argument of pidChargedBDTScore must be an integer (PDG code).");
819 return nullptr;
820 }
821 Const::ChargedStable hypType = Const::ChargedStable(hypPdgId);
822
823 std::vector<std::string> detectors(arguments.begin() + 1, arguments.end());
824 Const::PIDDetectorSet detectorSet = parseDetectorsChargedBDT(detectors);
825
826 auto func = [hypType, detectorSet](const Particle * part) -> double {
827 auto name = "pidChargedBDTScore_" + std::to_string(hypType.getPDGCode());
828 for (const Const::EDetector& detector : detectorSet)
829 {
830 name += "_" + std::to_string(detector);
831 }
832 return (part->hasExtraInfo(name)) ? part->getExtraInfo(name) : Const::doubleNaN;
833 };
834 return func;
835 }
836
837 Manager::FunctionPtr pidPairChargedBDTScore(const std::vector<std::string>& arguments)
838 {
839 if (arguments.size() != 3) {
840 B2ERROR("Need exactly three arguments for pidPairChargedBDTScore: pdgCodeHyp, pdgCodeTest, detector.");
841 return nullptr;
842 }
843
844 int hypPdgId, testPdgId;
845 try {
846 hypPdgId = Belle2::convertString<int>(arguments.at(0));
847 } catch (std::invalid_argument& e) {
848 B2ERROR("First argument of pidPairChargedBDTScore must be an integer (PDG code).");
849 return nullptr;
850 }
851 try {
852 testPdgId = Belle2::convertString<int>(arguments.at(1));
853 } catch (std::invalid_argument& e) {
854 B2ERROR("First argument of pidPairChargedBDTScore must be an integer (PDG code).");
855 return nullptr;
856 }
857 Const::ChargedStable hypType = Const::ChargedStable(hypPdgId);
858 Const::ChargedStable testType = Const::ChargedStable(testPdgId);
859
860 std::vector<std::string> detectors(arguments.begin() + 2, arguments.end());
861 Const::PIDDetectorSet detectorSet = parseDetectorsChargedBDT(detectors);
862
863 auto func = [hypType, testType, detectorSet](const Particle * part) -> double {
864 auto name = "pidPairChargedBDTScore_" + std::to_string(hypType.getPDGCode()) + "_VS_" + std::to_string(testType.getPDGCode());
865 for (const Const::EDetector& detector : detectorSet)
866 {
867 name += "_" + std::to_string(detector);
868 }
869 return (part->hasExtraInfo(name)) ? part->getExtraInfo(name) : Const::doubleNaN;
870 };
871 return func;
872 }
873
874 double mostLikelyPDG(const Particle* part, const std::vector<double>& arguments)
875 {
876 if (arguments.size() != 0 and arguments.size() != Const::ChargedStable::c_SetSize) {
877 B2ERROR("Need zero or exactly " << Const::ChargedStable::c_SetSize << " arguments for pidMostLikelyPDG");
878 return Const::doubleNaN;
879 }
881 if (arguments.size() == 0) {
882 for (unsigned int i = 0; i < Const::ChargedStable::c_SetSize; i++) prob[i] = 1. / Const::ChargedStable::c_SetSize;
883 } else {
884 copy(arguments.begin(), arguments.end(), prob);
885 }
886
887 auto* pid = part->getPIDLikelihood();
888 if (!pid) return Const::doubleNaN;
889 return pid->getMostLikely(prob).getPDGCode();
890 }
891
892 bool isMostLikely(const Particle* part, const std::vector<double>& arguments)
893 {
894 if (arguments.size() != 0 and arguments.size() != Const::ChargedStable::c_SetSize) {
895 B2ERROR("Need zero or exactly " << Const::ChargedStable::c_SetSize << " arguments for pidIsMostLikely");
896 return false;
897 }
898 return mostLikelyPDG(part, arguments) == abs(part->getPDGCode());
899 }
900
901 Manager::FunctionPtr weightedElectronID(const std::vector<std::string>& arguments)
902 {
903 std::string varName;
904 if (arguments.size() == 0) {
905 varName = "pidWeightedProbabilityExpert(PIDCalibrationWeightMatrix, 11, ALL)";
906 } else if (arguments.size() == 1) {
907 varName = "pidWeightedProbabilityExpert(" + arguments[0] + ", 11, ALL)";
908 } else {
909 B2ERROR("Need zero or one argument for weightedElectronID");
910 return nullptr;
911 }
912
913 const Variable::Manager::Var* var = Manager::Instance().getVariable(varName);
914 auto func = [var](const Particle * particle) -> double {
915 return std::get<double>(var->function(particle));
916 };
917 return func;
918 };
919
920 Manager::FunctionPtr weightedMuonID(const std::vector<std::string>& arguments)
921 {
922 std::string varName;
923 if (arguments.size() == 0) {
924 varName = "pidWeightedProbabilityExpert(PIDCalibrationWeightMatrix, 13, ALL)";
925 } else if (arguments.size() == 1) {
926 varName = "pidWeightedProbabilityExpert(" + arguments[0] + ", 13, ALL)";
927 } else {
928 B2ERROR("Need zero or one argument for weightedMuonID");
929 return nullptr;
930 }
931
932 const Variable::Manager::Var* var = Manager::Instance().getVariable(varName);
933 auto func = [var](const Particle * particle) -> double {
934 return std::get<double>(var->function(particle));
935 };
936 return func;
937 };
938
939 Manager::FunctionPtr weightedPionID(const std::vector<std::string>& arguments)
940 {
941 std::string varName;
942 if (arguments.size() == 0) {
943 varName = "pidWeightedProbabilityExpert(PIDCalibrationWeightMatrix, 211, ALL)";
944 } else if (arguments.size() == 1) {
945 varName = "pidWeightedProbabilityExpert(" + arguments[0] + ", 211, ALL)";
946 } else {
947 B2ERROR("Need zero or one argument for weightedPionID");
948 return nullptr;
949 }
950
951 const Variable::Manager::Var* var = Manager::Instance().getVariable(varName);
952 auto func = [var](const Particle * particle) -> double {
953 return std::get<double>(var->function(particle));
954 };
955 return func;
956 };
957
958 Manager::FunctionPtr weightedKaonID(const std::vector<std::string>& arguments)
959 {
960 std::string varName;
961 if (arguments.size() == 0) {
962 varName = "pidWeightedProbabilityExpert(PIDCalibrationWeightMatrix, 321, ALL)";
963 } else if (arguments.size() == 1) {
964 varName = "pidWeightedProbabilityExpert(" + arguments[0] + ", 321, ALL)";
965 } else {
966 B2ERROR("Need zero or one argument for weightedKaonID");
967 return nullptr;
968 }
969
970 const Variable::Manager::Var* var = Manager::Instance().getVariable(varName);
971 auto func = [var](const Particle * particle) -> double {
972 return std::get<double>(var->function(particle));
973 };
974 return func;
975 };
976
977 Manager::FunctionPtr weightedProtonID(const std::vector<std::string>& arguments)
978 {
979 std::string varName;
980 if (arguments.size() == 0) {
981 varName = "pidWeightedProbabilityExpert(PIDCalibrationWeightMatrix, 2212, ALL)";
982 } else if (arguments.size() == 1) {
983 varName = "pidWeightedProbabilityExpert(" + arguments[0] + ", 2212, ALL)";
984 } else {
985 B2ERROR("Need zero or one argument for weightedProtonID");
986 return nullptr;
987 }
988
989 const Variable::Manager::Var* var = Manager::Instance().getVariable(varName);
990 auto func = [var](const Particle * particle) -> double {
991 return std::get<double>(var->function(particle));
992 };
993 return func;
994 };
995
996 Manager::FunctionPtr weightedDeuteronID(const std::vector<std::string>& arguments)
997 {
998 std::string varName;
999 if (arguments.size() == 0) {
1000 varName = "pidWeightedProbabilityExpert(PIDCalibrationWeightMatrix, 1000010020, ALL)";
1001 } else if (arguments.size() == 1) {
1002 varName = "pidWeightedProbabilityExpert(" + arguments[0] + ", 1000010020, ALL)";
1003 } else {
1004 B2ERROR("Need zero or one argument for weightedDeuteronID");
1005 return nullptr;
1006 }
1007
1008 const Variable::Manager::Var* var = Manager::Instance().getVariable(varName);
1009 auto func = [var](const Particle * particle) -> double {
1010 return std::get<double>(var->function(particle));
1011 };
1012 return func;
1013 };
1014
1015
1016 double pionIDNN(const Particle* particle)
1017 {
1018 if (particle->hasExtraInfo("pidNeuralNetworkValueExpert(211)"))
1019 return particle->getExtraInfo("pidNeuralNetworkValueExpert(211)");
1020 static auto func = pidNeuralNetworkValueExpert({"211"});
1021 return std::get<double>(func(particle));
1022 }
1023
1024 double kaonIDNN(const Particle* particle)
1025 {
1026 if (particle->hasExtraInfo("pidNeuralNetworkValueExpert(321)"))
1027 return particle->getExtraInfo("pidNeuralNetworkValueExpert(321)");
1028 static auto func = pidNeuralNetworkValueExpert({"321"});
1029 return std::get<double>(func(particle));
1030 }
1031
1032 double klmMuonIDDNN(const Particle* part)
1033 {
1034 const PIDLikelihood* pid = part->getPIDLikelihood();
1035 if (!pid) return Const::doubleNaN;
1036 double klmMuonIDDNNvalue = pid->getPreOfficialLikelihood("klmMuonIDDNN");
1037 // klmMuonIDDNNvalue == -1.0 means there is no valid output from KLMMuonIDDNNExpertModule
1038 if (klmMuonIDDNNvalue < 0.0) return Const::doubleNaN;
1039 return klmMuonIDDNNvalue;
1040 }
1041
1042 //*************
1043 // B2BII
1044 //*************
1045
1046 double muIDBelle(const Particle* particle)
1047 {
1048 const PIDLikelihood* pid = particle->getPIDLikelihood();
1049 if (!pid) return 0.5; // Belle standard
1050
1051 if (pid->isAvailable(Const::KLM))
1052 return exp(pid->getLogL(Const::muon, Const::KLM));
1053 else
1054 return 0; // Belle standard
1055 }
1056
1057 double muIDBelleQuality(const Particle* particle)
1058 {
1059 const PIDLikelihood* pid = particle->getPIDLikelihood();
1060 if (!pid) return 0;// Belle standard
1061
1062 return pid->isAvailable(Const::KLM);
1063 }
1064
1065 double atcPIDBelle(const Particle* particle, const std::vector<double>& sigAndBkgHyp)
1066 {
1067 int sigHyp = int(std::lround(sigAndBkgHyp[0]));
1068 int bkgHyp = int(std::lround(sigAndBkgHyp[1]));
1069
1070 const PIDLikelihood* pid = particle->getPIDLikelihood();
1071 if (!pid) return 0.5; // Belle standard
1072
1073 // ACC = ARICH
1074 Const::PIDDetectorSet set = Const::ARICH;
1075 double acc_sig = exp(pid->getLogL(hypothesisConversion(sigHyp), set));
1076 double acc_bkg = exp(pid->getLogL(hypothesisConversion(bkgHyp), set));
1077 double acc = 0.5; // Belle standard
1078 if (acc_sig + acc_bkg > 0.0)
1079 acc = acc_sig / (acc_sig + acc_bkg);
1080
1081 // TOF = TOP
1082 set = Const::TOP;
1083 double tof_sig = exp(pid->getLogL(hypothesisConversion(sigHyp), set));
1084 double tof_bkg = exp(pid->getLogL(hypothesisConversion(bkgHyp), set));
1085 double tof = 0.5; // Belle standard
1086 double tof_all = tof_sig + tof_bkg;
1087 if (tof_all != 0) {
1088 tof = tof_sig / tof_all;
1089 if (tof < 0.001) tof = 0.001;
1090 if (tof > 0.999) tof = 0.999;
1091 }
1092
1093 // dE/dx = CDC
1094 set = Const::CDC;
1095 double cdc_sig = exp(pid->getLogL(hypothesisConversion(sigHyp), set));
1096 double cdc_bkg = exp(pid->getLogL(hypothesisConversion(bkgHyp), set));
1097 double cdc = 0.5; // Belle standard
1098 double cdc_all = cdc_sig + cdc_bkg;
1099 if (cdc_all != 0) {
1100 cdc = cdc_sig / cdc_all;
1101 if (cdc < 0.001) cdc = 0.001;
1102 if (cdc > 0.999) cdc = 0.999;
1103 }
1104
1105 // Combined
1106 double pid_sig = acc * tof * cdc;
1107 double pid_bkg = (1. - acc) * (1. - tof) * (1. - cdc);
1108
1109 return pid_sig / (pid_sig + pid_bkg);
1110 }
1111
1112
1113 double eIDBelle(const Particle* part)
1114 {
1115 const PIDLikelihood* pid = part->getPIDLikelihood();
1116 if (!pid) return 0.5; // Belle standard
1117
1118 Const::PIDDetectorSet set = Const::ECL;
1119 return pid->getProbability(Const::electron, Const::pion, set);
1120 }
1121
1122
1123 // PID variables to be used for analysis
1124 VARIABLE_GROUP("PID");
1125 REGISTER_VARIABLE("particleID", particleID,
1126 "the particle identification probability under the particle's own hypothesis, using info from all available detectors");
1127 REGISTER_VARIABLE("isPIDAvailable", isPIDAvailable,
1128 "True if PID is available (for at least one of the PID detectors");
1129 REGISTER_METAVARIABLE("isPIDAvailableFrom(detectorList)", isPIDAvailableFrom,
1130 "True if PID is available for at least one of the detectors in the list)",
1131 Manager::VariableDataType::c_bool);
1132 REGISTER_VARIABLE("electronID", electronID,
1133 "electron identification probability defined as :math:`\\mathcal{L}_e/(\\mathcal{L}_e+\\mathcal{L}_\\mu+\\mathcal{L}_\\pi+\\mathcal{L}_K+\\mathcal{L}_p+\\mathcal{L}_d)`, using info from all available detectors");
1134 REGISTER_VARIABLE("muonID", muonID,
1135 "muon identification probability defined as :math:`\\mathcal{L}_\\mu/(\\mathcal{L}_e+\\mathcal{L}_\\mu+\\mathcal{L}_\\pi+\\mathcal{L}_K+\\mathcal{L}_p+\\mathcal{L}_d)`, using info from all available detectors");
1136 REGISTER_VARIABLE("pionID", pionID,
1137 "pion identification probability defined as :math:`\\mathcal{L}_\\pi/(\\mathcal{L}_e+\\mathcal{L}_\\mu+\\mathcal{L}_\\pi+\\mathcal{L}_K+\\mathcal{L}_p+\\mathcal{L}_d)`, using info from all available detectors");
1138 REGISTER_VARIABLE("kaonID", kaonID,
1139 "kaon identification probability defined as :math:`\\mathcal{L}_K/(\\mathcal{L}_e+\\mathcal{L}_\\mu+\\mathcal{L}_\\pi+\\mathcal{L}_K+\\mathcal{L}_p+\\mathcal{L}_d)`, using info from all available detectors");
1140 REGISTER_VARIABLE("protonID", protonID,
1141 "proton identification probability defined as :math:`\\mathcal{L}_p/(\\mathcal{L}_e+\\mathcal{L}_\\mu+\\mathcal{L}_\\pi+\\mathcal{L}_K+\\mathcal{L}_p+\\mathcal{L}_d)`, using info from all available detectors");
1142 REGISTER_VARIABLE("deuteronID", deuteronID,
1143 "deuteron identification probability defined as :math:`\\mathcal{L}_d/(\\mathcal{L}_e+\\mathcal{L}_\\mu+\\mathcal{L}_\\pi+\\mathcal{L}_K+\\mathcal{L}_p+\\mathcal{L}_d)`, using info from all available detectors");
1144 REGISTER_METAVARIABLE("binaryPID(pdgCode1, pdgCode2)", binaryPID,
1145 "Returns the binary probability for the first provided mass hypothesis with respect to the second mass hypothesis using all detector components",
1146 Manager::VariableDataType::c_double);
1147 REGISTER_METAVARIABLE("pidChargedBDTScore(pdgCodeHyp, detector)", pidChargedBDTScore,
1148 "Returns the charged Pid BDT score for a certain mass hypothesis with respect to all other charged stable particle hypotheses. The second argument specifies which BDT training to use: based on 'ALL' PID detectors (NB: 'SVD' is currently excluded), or 'ECL' only. The choice depends on the ChargedPidMVAMulticlassModule's configuration.",
1149 Manager::VariableDataType::c_double);
1150 REGISTER_METAVARIABLE("pidPairChargedBDTScore(pdgCodeHyp, pdgCodeTest, detector)", pidPairChargedBDTScore,
1151 "Returns the charged Pid BDT score for a certain mass hypothesis with respect to an alternative hypothesis. The second argument specifies which BDT training to use: based on 'ALL' PID detectors (NB: 'SVD' is currently excluded), or 'ECL' only. The choice depends on the ChargedPidMVAModule's configuration.",
1152 Manager::VariableDataType::c_double);
1153 REGISTER_VARIABLE("nbarID", antineutronID, R"DOC(
1154Returns MVA classifier for antineutron PID.
1155
1156- 1 signal(antineutron) like
1157- 0 background like
1158- -1 invalid using this PID due to some ECL variables used unavailable
1159
1160This PID is only for antineutron. Neutron is also considered as background.
1161The variables used are `clusterPulseShapeDiscriminationMVA`, `clusterE`, `clusterLAT`, `clusterE1E9`, `clusterE9E21`,
1162`clusterAbsZernikeMoment40`, `clusterAbsZernikeMoment51`, `clusterZernikeMVA`.)DOC");
1163
1164 // Special temporary variables defined for users' convenience.
1165 REGISTER_VARIABLE("electronID_noSVD", electronID_noSVD,
1166 "**(SPECIAL (TEMP) variable)** electron identification probability defined as :math:`\\mathcal{L}_e/(\\mathcal{L}_e+\\mathcal{L}_\\mu+\\mathcal{L}_\\pi+\\mathcal{L}_K+\\mathcal{L}_p+\\mathcal{L}_d)`, using info from all available detectors *excluding the SVD*");
1167 REGISTER_VARIABLE("muonID_noSVD", muonID_noSVD,
1168 "**(SPECIAL (TEMP) variable)** muon identification probability defined as :math:`\\mathcal{L}_\\mu/(\\mathcal{L}_e+\\mathcal{L}_\\mu+\\mathcal{L}_\\pi+\\mathcal{L}_K+\\mathcal{L}_p+\\mathcal{L}_d)`, using info from all available detectors *excluding the SVD*");
1169 REGISTER_VARIABLE("pionID_noSVD", pionID_noSVD,
1170 "**(SPECIAL (TEMP) variable)** pion identification probability defined as :math:`\\mathcal{L}_\\pi/(\\mathcal{L}_e+\\mathcal{L}_\\mu+\\mathcal{L}_\\pi+\\mathcal{L}_K+\\mathcal{L}_p+\\mathcal{L}_d)`, using info from all available detectors *excluding the SVD*");
1171 REGISTER_VARIABLE("kaonID_noSVD", kaonID_noSVD,
1172 "**(SPECIAL (TEMP) variable)** kaon identification probability defined as :math:`\\mathcal{L}_K/(\\mathcal{L}_e+\\mathcal{L}_\\mu+\\mathcal{L}_\\pi+\\mathcal{L}_K+\\mathcal{L}_p+\\mathcal{L}_d)`, using info from all available detectors *excluding the SVD*");
1173 REGISTER_VARIABLE("protonID_noSVD", protonID_noSVD,
1174 "**(SPECIAL (TEMP) variable)** proton identification probability defined as :math:`\\mathcal{L}_p/(\\mathcal{L}_e+\\mathcal{L}_\\mu+\\mathcal{L}_\\pi+\\mathcal{L}_K+\\mathcal{L}_p+\\mathcal{L}_d)`, using info from all available detectors *excluding the SVD*");
1175 REGISTER_VARIABLE("deuteronID_noSVD", deuteronID_noSVD,
1176 "**(SPECIAL (TEMP) variable)** deuteron identification probability defined as :math:`\\mathcal{L}_d/(\\mathcal{L}_e+\\mathcal{L}_\\mu+\\mathcal{L}_\\pi+\\mathcal{L}_K+\\mathcal{L}_p+\\mathcal{L}_d)`, using info from all available detectors *excluding the SVD*");
1177 REGISTER_METAVARIABLE("binaryPID_noSVD(pdgCode1, pdgCode2)", binaryPID_noSVD,
1178 "Returns the binary probability for the first provided mass hypothesis with respect to the second mass hypothesis using all detector components, *excluding the SVD*.",
1179 Manager::VariableDataType::c_double);
1180 REGISTER_VARIABLE("electronID_noTOP", electronID_noTOP,
1181 "**(SPECIAL (TEMP) variable)** electron identification probability defined as :math:`\\mathcal{L}_e/(\\mathcal{L}_e+\\mathcal{L}_\\mu+\\mathcal{L}_\\pi+\\mathcal{L}_K+\\mathcal{L}_p+\\mathcal{L}_d)`, using info from all available detectors *excluding the TOP*. *NB:* this variable must be used in place of `electronID` when analysing data (MC) processed (simulated) in *release 6*");
1182 REGISTER_METAVARIABLE("binaryElectronID_noTOP(pdgCodeTest)", binaryElectronID_noTOP,
1183 "**(SPECIAL (TEMP) variable)** Returns the binary probability for the electron mass hypothesis with respect to another mass hypothesis using all detector components, *excluding the TOP*. *NB:* this variable must be used in place of `binaryPID` (``pdgCode1=11``) when analysing data (MC) processed (simulated) in **release 6**",
1184 Manager::VariableDataType::c_double);
1185 REGISTER_VARIABLE("electronID_noSVD_noTOP", electronID_noSVD_noTOP,
1186 "**(SPECIAL (TEMP) variable)** electron identification probability defined as :math:`\\mathcal{L}_e/(\\mathcal{L}_e+\\mathcal{L}_\\mu+\\mathcal{L}_\\pi+\\mathcal{L}_K+\\mathcal{L}_p+\\mathcal{L}_d)`, using info from all available detectors *excluding the SVD and the TOP*. *NB:* this variable must be used in place of `electronID` when analysing data (MC) processed (simulated) in *release 5*");
1187 REGISTER_METAVARIABLE("binaryElectronID_noSVD_noTOP(pdgCodeTest)", binaryElectronID_noSVD_noTOP,
1188 "**(SPECIAL (TEMP) variable)** Returns the binary probability for the electron mass hypothesis with respect to another mass hypothesis using all detector components, *excluding the SVD and the TOP*. *NB:* this variable must be used in place of `binaryPID` (``pdgCode1=11``) when analysing data (MC) processed (simulated) in **release 5**",
1189 Manager::VariableDataType::c_double);
1190 REGISTER_VARIABLE("pionID_noARICHwoECL", pionID_noARICHwoECL,
1191 "**(SPECIAL (TEMP) variable)** pion identification probability defined as :math:`\\mathcal{L}_\\pi/(\\mathcal{L}_e+\\mathcal{L}_\\mu+\\mathcal{L}_\\pi+\\mathcal{L}_K+\\mathcal{L}_p+\\mathcal{L}_d)`, using info from all available detectors but ARICH info excluded for tracks without associated ECL cluster");
1192 REGISTER_VARIABLE("kaonID_noARICHwoECL", kaonID_noARICHwoECL,
1193 "**(SPECIAL (TEMP) variable)** kaon identification probability defined as :math:`\\mathcal{L}_K/(\\mathcal{L}_e+\\mathcal{L}_\\mu+\\mathcal{L}_\\pi+\\mathcal{L}_K+\\mathcal{L}_p+\\mathcal{L}_d)`, using info from all available detectors but ARICH info excluded for tracks without associated ECL cluster");
1194 REGISTER_METAVARIABLE("binaryPID_noARICHwoECL(pdgCode1, pdgCode2)", binaryPID_noARICHwoECL,
1195 "Returns the binary probability for the first provided mass hypothesis with respect to the second mass hypothesis using all detector components, but ARICH info excluded for tracks without associated ECL cluster",
1196 Manager::VariableDataType::c_double);
1197
1198
1199 REGISTER_METAVARIABLE("weightedElectronID(weightMatrixName)", weightedElectronID,
1200 R"DOC(
1201weighted electron identification probability defined as :math:`\frac{\mathcal{\tilde{L}}_e}{\sum_{i=e,\mu,\pi,K,p,d} \mathcal{\tilde{L}}_i}`,
1202where :math:`\mathcal{\tilde{L}}_i` is defined as :math:`\log\mathcal{\tilde{L}}_i = \sum_{j={\mathrm{SVD, CDC, TOP, ARICH, ECL, KLM}}} \mathcal{w}_{ij}\log\mathcal{L}_{ij}`.
1203The :math:`\mathcal{L}_{ij}` is the original likelihood and :math:`\mathcal{w}_{ij}` is the PID calibration weight of i-th particle type and j-th detector.
1204One can provide the name of the weight matrix as the argument.
1205)DOC",
1206 Manager::VariableDataType::c_double);
1207 REGISTER_METAVARIABLE("weightedMuonID(weightMatrixName)", weightedMuonID,
1208 R"DOC(
1209weighted muon identification probability defined as :math:`\frac{\mathcal{\tilde{L}}_\mu}{\sum_{i=e,\mu,\pi,K,p,d} \mathcal{\tilde{L}}_i}`,
1210where :math:`\mathcal{\tilde{L}}_i` is defined as :math:`\log\mathcal{\tilde{L}}_i = \sum_{j={\mathrm{SVD, CDC, TOP, ARICH, ECL, KLM}}} \mathcal{w}_{ij}\log\mathcal{L}_{ij}`.
1211The :math:`\mathcal{L}_{ij}` is the original likelihood and :math:`\mathcal{w}_{ij}` is the PID calibration weight of i-th particle type and j-th detector.
1212One can provide the name of the weight matrix as the argument.
1213)DOC",
1214 Manager::VariableDataType::c_double);
1215 REGISTER_METAVARIABLE("weightedPionID(weightMatrixName)", weightedPionID,
1216 R"DOC(
1217weighted pion identification probability defined as :math:`\frac{\mathcal{\tilde{L}}_\pi}{\sum_{i=e,\mu,\pi,K,p,d} \mathcal{\tilde{L}}_i}`,
1218where :math:`\mathcal{\tilde{L}}_i` is defined as :math:`\log\mathcal{\tilde{L}}_i = \sum_{j={\mathrm{SVD, CDC, TOP, ARICH, ECL, KLM}}} \mathcal{w}_{ij}\log\mathcal{L}_{ij}`.
1219The :math:`\mathcal{L}_{ij}` is the original likelihood and :math:`\mathcal{w}_{ij}` is the PID calibration weight of i-th particle type and j-th detector.
1220One can provide the name of the weight matrix as the argument.
1221)DOC",
1222 Manager::VariableDataType::c_double);
1223 REGISTER_METAVARIABLE("weightedKaonID(weightMatrixName)", weightedKaonID,
1224 R"DOC(
1225weighted kaon identification probability defined as :math:`\frac{\mathcal{\tilde{L}}_K}{\sum_{i=e,\mu,\pi,K,p,d} \mathcal{\tilde{L}}_i}`,
1226where :math:`\mathcal{\tilde{L}}_i` is defined as :math:`\log\mathcal{\tilde{L}}_i = \sum_{j={\mathrm{SVD, CDC, TOP, ARICH, ECL, KLM}}} \mathcal{w}_{ij}\log\mathcal{L}_{ij}`.
1227The :math:`\mathcal{L}_{ij}` is the original likelihood and :math:`\mathcal{w}_{ij}` is the PID calibration weight of i-th particle type and j-th detector.
1228One can provide the name of the weight matrix as the argument.
1229)DOC",
1230 Manager::VariableDataType::c_double);
1231 REGISTER_METAVARIABLE("weightedProtonID(weightMatrixName)", weightedProtonID,
1232 R"DOC(
1233weighted proton identification probability defined as :math:`\frac{\mathcal{\tilde{L}}_p}{\sum_{i=e,\mu,\pi,K,p,d} \mathcal{\tilde{L}}_i}`,
1234where :math:`\mathcal{\tilde{L}}_i` is defined as :math:`\log\mathcal{\tilde{L}}_i = \sum_{j={\mathrm{SVD, CDC, TOP, ARICH, ECL, KLM}}} \mathcal{w}_{ij}\log\mathcal{L}_{ij}`.
1235The :math:`\mathcal{L}_{ij}` is the original likelihood and :math:`\mathcal{w}_{ij}` is the PID calibration weight of i-th particle type and j-th detector.
1236One can provide the name of the weight matrix as the argument.
1237)DOC",
1238 Manager::VariableDataType::c_double);
1239 REGISTER_METAVARIABLE("weightedDeuteronID(weightMatrixName)", weightedDeuteronID,
1240 R"DOC(
1241weighted deuteron identification probability defined as :math:`\frac{\mathcal{\tilde{L}}_d}{\sum_{i=e,\mu,\pi,K,p,d} \mathcal{\tilde{L}}_i}`,
1242where :math:`\mathcal{\tilde{L}}_i` is defined as :math:`\log\mathcal{\tilde{L}}_i = \sum_{j={\mathrm{SVD, CDC, TOP, ARICH, ECL, KLM}}} \mathcal{w}_{ij}\log\mathcal{L}_{ij}`.
1243The :math:`\mathcal{L}_{ij}` is the original likelihood and :math:`\mathcal{w}_{ij}` is the PID calibration weight of i-th particle type and j-th detector.
1244One can provide the name of the weight matrix as the argument.
1245)DOC",
1246 Manager::VariableDataType::c_double);
1247
1248
1249 REGISTER_VARIABLE("pionIDNN", pionIDNN,
1250 R"DOC(
1251pion identification probability as calculated from the PID neural network.
1252)DOC");
1253 REGISTER_VARIABLE("kaonIDNN", kaonIDNN,
1254 R"DOC(
1255kaon identification probability as calculated from the PID neural network.
1256)DOC");
1257
1258
1259 // Metafunctions for experts to access the basic PID quantities
1260 VARIABLE_GROUP("PID_expert");
1261 REGISTER_METAVARIABLE("pidLogLikelihoodValueExpert(pdgCode, detectorList)", pidLogLikelihoodValueExpert,
1262 "returns the log likelihood value of for a specific mass hypothesis and set of detectors.", Manager::VariableDataType::c_double);
1263 REGISTER_METAVARIABLE("pidDeltaLogLikelihoodValueExpert(pdgCode1, pdgCode2, detectorList)", pidDeltaLogLikelihoodValueExpert,
1264 "returns LogL(hyp1) - LogL(hyp2) (aka DLL) for two mass hypotheses and a set of detectors.", Manager::VariableDataType::c_double);
1265 REGISTER_METAVARIABLE("pidPairProbabilityExpert(pdgCodeHyp, pdgCodeTest, detectorList)", pidPairProbabilityExpert,
1266 "Pair (or binary) probability for the pdgCodeHyp mass hypothesis respect to the pdgCodeTest one, using an arbitrary set of detectors. :math:`\\mathcal{L}_{hyp}/(\\mathcal{L}_{test}+\\mathcal{L}_{hyp})`",
1267 Manager::VariableDataType::c_double);
1268 REGISTER_METAVARIABLE("pidProbabilityExpert(pdgCodeHyp, detectorList)", pidProbabilityExpert,
1269 "probability for the pdgCodeHyp mass hypothesis respect to all the other ones, using an arbitrary set of detectors :math:`\\mathcal{L}_{hyp}/(\\Sigma_{\\text{all~hyp}}\\mathcal{L}_{i})`. ",
1270 Manager::VariableDataType::c_double);
1271 REGISTER_METAVARIABLE("pidLogarithmicProbabilityExpert(pdgCodeHyp, detectorList)", pidLogarithmicProbabilityExpert,
1272 "logarithmic equivalent of pidProbability (p) defined as log(p/(1-p)), which gives a smooth peak-like distribution",
1273 Manager::VariableDataType::c_double);
1274 REGISTER_METAVARIABLE("pidMissingProbabilityExpert(detectorList)", pidMissingProbabilityExpert,
1275 "returns 1 if PID is missing for at least one of the detectors in the list, otherwise 0. ",
1276 Manager::VariableDataType::c_double);
1277 REGISTER_VARIABLE("pidMostLikelyPDG(ePrior=1/6, muPrior=1/6, piPrior=1/6, KPrior=1/6, pPrior=1/6, dPrior=1/6)", mostLikelyPDG,
1278 R"DOC(
1279Returns PDG code of the largest PID likelihood, or NaN if PID information is not available.
1280This function accepts either no arguments, or 6 floats as priors for the charged particle hypotheses
1281following the order shown in the metavariable's declaration. Flat priors are assumed as default.)DOC");
1282 REGISTER_VARIABLE("pidIsMostLikely(ePrior=1/6, muPrior=1/6, piPrior=1/6, KPrior=1/6, pPrior=1/6, dPrior=1/6)", isMostLikely, R"DOC(
1283Returns True if the largest PID likelihood of a given particle corresponds to its particle hypothesis.
1284This function accepts either no arguments, or 6 floats as priors for the charged particle hypotheses
1285following the order shown in the metavariable's declaration. Flat priors are assumed as default.)DOC");
1286
1287 REGISTER_METAVARIABLE("pidWeightedLogLikelihoodValueExpert(weightMatrixName, pdgCode, detectorList)",
1288 pidWeightedLogLikelihoodValueExpert,
1289 "returns the weighted log likelihood value of for a specific mass hypothesis and set of detectors, "
1290 ":math:`\\log\\mathcal{\\tilde{L}}_{hyp} = \\sum_{j\\in\\mathrm{detectorList}} \\mathcal{w}_{hyp,j}\\log\\mathcal{L}_{hyp,j}`. "
1291 "The :math:`\\mathcal{L}_{ij}` is the original likelihood and :math:`\\mathcal{w}_{ij}` is the PID calibration weight of i-th particle type and j-th detector.",
1292 Manager::VariableDataType::c_double);
1293 REGISTER_METAVARIABLE("pidWeightedPairProbabilityExpert(weightMatrixName, pdgCodeHyp, pdgCodeTest, detectorList)",
1294 pidWeightedPairProbabilityExpert,
1295 "Weighted pair (or binary) probability for the pdgCodeHyp mass hypothesis with respect to the pdgCodeTest one, using an arbitrary set of detectors, "
1296 ":math:`\\mathcal{\\tilde{L}}_{hyp}/(\\mathcal{\\tilde{L}}_{test}+\\mathcal{\\tilde{L}}_{hyp})` where :math:`\\mathcal{\\tilde{L}}_{i}` is defined as "
1297 ":math:`\\log\\mathcal{\\tilde{L}}_{i} = \\sum_{j\\in\\mathrm{detectorList}} \\mathcal{w}_{i,j}\\log\\mathcal{L}_{i,j}`. "
1298 "The :math:`\\mathcal{L}_{ij}` is the original likelihood and :math:`\\mathcal{w}_{ij}` is the PID calibration weight of i-th particle type and j-th detector.",
1299 Manager::VariableDataType::c_double);
1300 REGISTER_METAVARIABLE("pidWeightedProbabilityExpert(weightMatrixName, pdgCodeHyp, detectorList)",
1301 pidWeightedProbabilityExpert,
1302 "Weighted probability for the pdgCodeHyp mass hypothesis with respect to all the other ones, using an arbitrary set of detectors, "
1303 ":math:`\\mathcal{\\tilde{L}}_{hyp}/\\sum_{i=e,\\mu,\\pi,K,p,d} \\mathcal{\\tilde{L}}_i` where :math:`\\mathcal{\\tilde{L}}_{i}` is defined as "
1304 ":math:`\\log\\mathcal{\\tilde{L}}_{i} = \\sum_{j\\in\\mathrm{detectorList}} \\mathcal{w}_{i,j}\\log\\mathcal{L}_{i,j}`. "
1305 "The :math:`\\mathcal{L}_{ij}` is the original likelihood and :math:`\\mathcal{w}_{ij}` is the PID calibration weight of i-th particle type and j-th detector.",
1306 Manager::VariableDataType::c_double);
1307 REGISTER_METAVARIABLE("pidNeuralNetworkValueExpert(pdgCodeHyp, PIDNeuralNetworkName)",
1308 pidNeuralNetworkValueExpert,
1309 "Probability for the particle hypothesis pdgCodeHype calculated from a neural network, which uses high-level information as inputs, "
1310 "such as the likelihood from the 6 subdetectors for PID for all 6 hypotheses, "
1311 ":math:`\\mathcal{\\tilde{L}}_{hyp}^{det}`, or the track momentum and charge",
1312 Manager::VariableDataType::c_double);
1313 REGISTER_VARIABLE("klmMuonIDDNN", klmMuonIDDNN,
1314 "Muon probability calculated from Neural Network with KLM information (expert use only)");
1315
1316 // B2BII PID
1317 VARIABLE_GROUP("Belle PID variables");
1318 REGISTER_METAVARIABLE("atcPIDBelle(i,j)", atcPIDBelle, R"DOC(
1319[Legacy] Returns Belle's PID atc variable: ``atc_pid(3,1,5,i,j).prob()``.
1320Parameters i,j are signal and background hypothesis: (0 = electron, 1 = muon, 2 = pion, 3 = kaon, 4 = proton)
1321Returns 0.5 in case there is no likelihood found and a factor of 0.5 will appear in the product if any of the subdetectors don't report a likelihood (Belle behaviour).
1322
1323.. warning:: The behaviour is different from Belle II PID variables which typically return NaN in case of error.
1324 )DOC", Manager::VariableDataType::c_double);
1325 REGISTER_VARIABLE("muIDBelle", muIDBelle, R"DOC(
1326[Legacy] Returns Belle's PID ``Muon_likelihood()`` variable.
1327Returns 0.5 in case there is no likelihood found and returns zero if the muon likelihood is not usable (Belle behaviour).
1328
1329.. warning:: The behaviour is different from Belle II PID variables which typically return NaN in case of error.
1330 )DOC");
1331 REGISTER_VARIABLE("muIDBelleQuality", muIDBelleQuality, R"DOC(
1332[Legacy] Returns true if Belle's PID ``Muon_likelihood()`` is usable (reliable).
1333Returns zero/false if not usable or if there is no PID found.
1334 )DOC");
1335 REGISTER_VARIABLE("eIDBelle", eIDBelle, R"DOC(
1336[Legacy] Returns Belle's electron ID ``eid(3,-1,5).prob()`` variable.
1337Returns 0.5 in case there is no likelihood found (Belle behaviour).
1338
1339.. warning:: The behaviour is different from Belle II PID variables which typically return NaN in case of error.
1340 )DOC");
1341 }
1343}
static const unsigned int c_SetSize
Number of elements (for use in array bounds etc.)
Definition: Const.h:615
static DetectorSet set()
Accessor function for the set of valid detectors.
Definition: Const.h:374
int getPDGCode() const
PDG code.
Definition: Const.h:473
static DetectorSet set()
Accessor for the set of valid detector IDs.
Definition: Const.h:333
static const ParticleType neutron
neutron particle
Definition: Const.h:675
static const ChargedStable muon
muon particle
Definition: Const.h:660
RestrictedDetectorSet< PIDDetectors > PIDDetectorSet
Typedef for set of PID detectors.
Definition: Const.h:379
static const ParticleSet chargedStableSet
set of charged stable particles
Definition: Const.h:618
EDetector
Enum for identifying the detector components (detector and subdetector).
Definition: Const.h:42
static const ChargedStable pion
charged pion particle
Definition: Const.h:661
static const ChargedStable proton
proton particle
Definition: Const.h:663
static const double doubleNaN
quiet_NaN
Definition: Const.h:703
static const ChargedStable kaon
charged kaon particle
Definition: Const.h:662
static const ChargedStable electron
electron particle
Definition: Const.h:659
static const ChargedStable deuteron
deuteron particle
Definition: Const.h:664
static const ReferenceFrame & GetCurrent()
Get current rest frame.
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
Abstract base class for different kinds of events.
STL namespace.
FunctionPtr function
Pointer to function.
Definition: Manager.h:147