Belle II Software development
ContinuumSuppressionVariables.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9#include <analysis/variables/ContinuumSuppressionVariables.h>
10
11#include <analysis/variables/ROEVariables.h>
12#include <analysis/VariableManager/Manager.h>
13#include <framework/dataobjects/EventExtraInfo.h>
14#include <analysis/dataobjects/Particle.h>
15#include <analysis/dataobjects/ContinuumSuppression.h>
16#include <analysis/utility/ReferenceFrame.h>
17#include <analysis/ClusterUtility/ClusterUtils.h>
18#include <analysis/ContinuumSuppression/FoxWolfram.h>
19
20#include <framework/logging/Logger.h>
21#include <framework/datastore/StoreArray.h>
22#include <framework/datastore/StoreObjPtr.h>
23#include <framework/utilities/Conversion.h>
24
25#include <mdst/dataobjects/PIDLikelihood.h>
26#include <mdst/dataobjects/Track.h>
27#include <mdst/dataobjects/ECLCluster.h>
28
29#include <cmath>
30
31
32namespace Belle2 {
37 namespace Variable {
38
39 Manager::FunctionPtr R2WithMask(const std::vector<std::string>& arguments)
40 {
41 if (arguments.size() != 1)
42 B2FATAL("An empty argument is not allowed for the variable R2."
43 "Either provide no argument or a valid mask name.");
44 std::string maskName = arguments[0];
45 auto func = [maskName](const Particle * particle) -> double {
46 const ContinuumSuppression* qq = particle->getRelatedTo<ContinuumSuppression>(maskName);
47 if (!qq)
48 return Const::doubleNaN;
49
50 return qq->getR2();
51 };
52 return func;
53 }
54
55 double R2(const Particle* particle)
56 {
57 RelationVector<ContinuumSuppression> continuumSuppressionRelations = particle->getRelationsTo<ContinuumSuppression>("ALL");
58 if (continuumSuppressionRelations.size() == 1) {
59 const ContinuumSuppression* qq = continuumSuppressionRelations[0];
60 return qq->getR2();
61 } else {
62 if (continuumSuppressionRelations.size() > 1) {
63 B2ERROR("The return value of R2 is ambiguous. Please provide the mask name as argument.");
64 }
65 return Const::doubleNaN;
66 }
67 }
68
69 Manager::FunctionPtr thrustBmWithMask(const std::vector<std::string>& arguments)
70 {
71 if (arguments.size() != 1)
72 B2FATAL("An empty argument is not allowed for the variable thrustBm."
73 "Either provide no argument or a valid mask name.");
74 std::string maskName = arguments[0];
75 auto func = [maskName](const Particle * particle) -> double {
76 const ContinuumSuppression* qq = particle->getRelatedTo<ContinuumSuppression>(maskName);
77 if (!qq)
78 return Const::doubleNaN;
79
80 return qq->getThrustBm();
81 };
82 return func;
83 }
84
85 double thrustBm(const Particle* particle)
86 {
87 RelationVector<ContinuumSuppression> continuumSuppressionRelations = particle->getRelationsTo<ContinuumSuppression>("ALL");
88 if (continuumSuppressionRelations.size() == 1) {
89 const ContinuumSuppression* qq = continuumSuppressionRelations[0];
90 return qq->getThrustBm();
91 } else {
92 if (continuumSuppressionRelations.size() > 1) {
93 B2ERROR("The return value of thrustBm is ambiguous. Please provide the mask name as argument.");
94 }
95 return Const::doubleNaN;
96 }
97 }
98
99 Manager::FunctionPtr thrustOmWithMask(const std::vector<std::string>& arguments)
100 {
101 if (arguments.size() != 1)
102 B2FATAL("An empty argument is not allowed for the variable thrustOm."
103 "Either provide no argument or a valid mask name.");
104 std::string maskName = arguments[0];
105 auto func = [maskName](const Particle * particle) -> double {
106 const ContinuumSuppression* qq = particle->getRelatedTo<ContinuumSuppression>(maskName);
107 if (!qq)
108 return Const::doubleNaN;
109
110 return qq->getThrustOm();
111 };
112 return func;
113 }
114
115 double thrustOm(const Particle* particle)
116 {
117 RelationVector<ContinuumSuppression> continuumSuppressionRelations = particle->getRelationsTo<ContinuumSuppression>("ALL");
118 if (continuumSuppressionRelations.size() == 1) {
119 const ContinuumSuppression* qq = continuumSuppressionRelations[0];
120 return qq->getThrustOm();
121 } else {
122 if (continuumSuppressionRelations.size() > 1) {
123 B2ERROR("The return value of thrustOm is ambiguous. Please provide the mask name as argument.");
124 }
125 return Const::doubleNaN;
126 }
127 }
128
129 Manager::FunctionPtr cosTBTOWithMask(const std::vector<std::string>& arguments)
130 {
131 if (arguments.size() != 1)
132 B2FATAL("An empty argument is not allowed for the variable cosTBTO."
133 "Either provide no argument or a valid mask name.");
134 std::string maskName = arguments[0];
135 auto func = [maskName](const Particle * particle) -> double {
136 const ContinuumSuppression* qq = particle->getRelatedTo<ContinuumSuppression>(maskName);
137 if (!qq)
138 return Const::doubleNaN;
139
140 return qq->getCosTBTO();
141 };
142 return func;
143 }
144
145 double cosTBTO(const Particle* particle)
146 {
147 RelationVector<ContinuumSuppression> continuumSuppressionRelations = particle->getRelationsTo<ContinuumSuppression>("ALL");
148 if (continuumSuppressionRelations.size() == 1) {
149 const ContinuumSuppression* qq = continuumSuppressionRelations[0];
150 return qq->getCosTBTO();
151 } else {
152 if (continuumSuppressionRelations.size() > 1) {
153 B2ERROR("The return value of cosTBTO is ambiguous. Please provide the mask name as argument.");
154 }
155 return Const::doubleNaN;
156 }
157 }
158
159 Manager::FunctionPtr cosTBzWithMask(const std::vector<std::string>& arguments)
160 {
161 if (arguments.size() != 1)
162 B2FATAL("An empty argument is not allowed for the variable cosTBz."
163 "Either provide no argument or a valid mask name.");
164 std::string maskName = arguments[0];
165 auto func = [maskName](const Particle * particle) -> double {
166 const ContinuumSuppression* qq = particle->getRelatedTo<ContinuumSuppression>(maskName);
167 if (!qq)
168 return Const::doubleNaN;
169
170 return qq->getCosTBz();
171 };
172 return func;
173 }
174
175 double cosTBz(const Particle* particle)
176 {
177 RelationVector<ContinuumSuppression> continuumSuppressionRelations = particle->getRelationsTo<ContinuumSuppression>("ALL");
178 if (continuumSuppressionRelations.size() == 1) {
179 const ContinuumSuppression* qq = continuumSuppressionRelations[0];
180 return qq->getCosTBz();
181 } else {
182 if (continuumSuppressionRelations.size() > 1) {
183 B2ERROR("The return value of cosTBz is ambiguous. Please provide the mask name as argument.");
184 }
185 return Const::doubleNaN;
186 }
187 }
188
189 Manager::FunctionPtr KSFWVariables(const std::vector<std::string>& arguments)
190 {
191 if (arguments.size() > 0 && arguments.size() < 4) {
192 bool useFS1 = false;
193 auto variableName = arguments[0];
194 std::string maskName = "";
195 if (arguments.size() >= 2) {
196 if (arguments[1] == "FS1") {
197 useFS1 = true;
198 } else {
199 maskName = arguments[1];
200 }
201 if (arguments.size() == 3) {
202 maskName = arguments[2];
203 if (maskName == "FS1") {
204 B2ERROR("It looks like you provided the arguments for KSFWVariables in the wrong order."
205 "If you want to use the KSFW moments calculated from the B final state particles, the second argument has to be 'FS1'."
206 "The third argument would then have to be the ROE mask name.");
207 }
208 }
209 }
210 int index = -1;
211
212 // all possible names
213 std::vector<std::string> names = {"mm2", "et",
214 "hso00", "hso01", "hso02", "hso03", "hso04",
215 "hso10", "hso12", "hso14",
216 "hso20", "hso22", "hso24",
217 "hoo0", "hoo1", "hoo2", "hoo3", "hoo4"
218 };
219
220 // find the index of the name
221 for (unsigned i = 0; i < names.size(); ++i) {
222 if (variableName == names[i])
223 index = i;
224 }
225
226 // throw helpful error if name provided was not in allowed list
227 if (index == -1) {
228 std::string allowed = "";
229 for (auto n : names)
230 allowed += n + ", ";
231 B2FATAL("Variable name provided: " << variableName << " is not one of the allowed options. Please choose from one of:" << allowed);
232 }
233
234 auto func = [index, useFS1, maskName](const Particle * particle) -> double {
235 RelationVector<ContinuumSuppression> continuumSuppressionRelations = particle->getRelationsTo<ContinuumSuppression>("ALL");
236 ContinuumSuppression* qq = nullptr;
237 if (maskName.empty())
238 {
239 if (continuumSuppressionRelations.size() == 1) {
240 qq = continuumSuppressionRelations[0];
241 } else if (continuumSuppressionRelations.size() > 1) {
242 B2ERROR("The return value of KSFWVariables is ambiguous. Please provide the mask name as argument.");
243 }
244 } else
245 {
246 qq = particle->getRelatedTo<ContinuumSuppression>(maskName);
247 }
248 if (!qq) return Const::doubleNaN;
249
250 // get the KSFW moments
251 std::vector<float> ksfw = qq->getKsfwFS0();
252 if (useFS1)
253 ksfw = qq->getKsfwFS1();
254
255 if (ksfw.size() == 0) B2FATAL("Could not find any KSFW moments");
256 return ksfw.at(index);
257 };
258 return func;
259 } else {
260 B2FATAL("Wrong number of arguments for meta function KSFWVariables. It only takes between one and three arguments."
261 " The first argument must be the variable. If you want to use the KSFW moments calculated from the B final state particles, set 'FS1' as second argument."
262 " You can also provide the ROE mask name as second or third argument.");
263 }
264 }
265
266 Manager::FunctionPtr CleoConesCS(const std::vector<std::string>& arguments)
267 {
268 if (arguments.size() > 0 && arguments.size() < 4) {
269
270 int coneNumber = 0;
271 try {
272 coneNumber = Belle2::convertString<int>(arguments[0]);
273 } catch (std::invalid_argument&) {
274 B2FATAL("The first argument of the CleoConeCS meta function must be an integer!");
275 }
276
277 bool useROE = false;
278 std::string maskName = "";
279 if (arguments.size() >= 2) {
280 if (arguments[1] == "ROE") {
281 useROE = true;
282 } else {
283 maskName = arguments[1];
284 }
285 if (arguments.size() == 3) {
286 maskName = arguments[2];
287 if (maskName == "ROE") {
288 B2ERROR("It looks like you provided the arguments for CleoConeCS in the wrong order."
289 "If you want to use the CleoCones calculated from all final state particles, the second argument has to be 'ROE'."
290 "The third argument would then have to be the ROE mask name.");
291 }
292 }
293 }
294
295 auto func = [coneNumber, useROE, maskName](const Particle * particle) -> double {
296 RelationVector<ContinuumSuppression> continuumSuppressionRelations = particle->getRelationsTo<ContinuumSuppression>("ALL");
297 ContinuumSuppression* qq = nullptr;
298 if (maskName.empty())
299 {
300 if (continuumSuppressionRelations.size() == 1) {
301 qq = continuumSuppressionRelations[0];
302 } else if (continuumSuppressionRelations.size() > 1) {
303 B2ERROR("The return value of CleoConeCS is ambiguous. Please provide the mask name as argument.");
304 }
305 } else
306 {
307 qq = particle->getRelatedTo<ContinuumSuppression>(maskName);
308 }
309 if (!qq)
310 return Const::doubleNaN;
311
312 std::vector<float> cleoCones = qq->getCleoConesALL();
313 if (useROE)
314 cleoCones = qq->getCleoConesROE();
315 return cleoCones.at(coneNumber - 1);
316 };
317 return func;
318 } else {
319 B2FATAL("Wrong number of arguments for CleoConeCS function. It only takes between one and three arguments."
320 "The first argument must be the cone number. If you want to use the CleoCones calculated from all final state particles, set 'ROE' as second argument."
321 "You can also provide the ROE mask name as second or third argument.");
322 }
323 }
324
325 Manager::FunctionPtr transformedNetworkOutput(const std::vector<std::string>& arguments)
326 {
327 if (arguments.size() == 3) {
328 double low = 0;
329 double high = 0;
330 try {
331 low = Belle2::convertString<double>(arguments[1]);
332 high = Belle2::convertString<double>(arguments[2]);
333 } catch (std::invalid_argument&) {
334 B2FATAL("Second and third argument of transformedNetworkOutput meta function must be doubles!");
335 }
336
337 auto extraInfoName = arguments[0];
338 auto func = [extraInfoName, low, high](const Particle * particle) -> double {
339 if (particle == nullptr)
340 {
341 StoreObjPtr<EventExtraInfo> eventExtraInfo;
342 if (eventExtraInfo->hasExtraInfo(extraInfoName)) {
343 return eventExtraInfo->getExtraInfo(extraInfoName);
344 } else {
345 return Const::doubleNaN;
346 }
347 }
348 if (particle->hasExtraInfo(extraInfoName))
349 {
350 return std::log(((particle->getExtraInfo(extraInfoName)) - low) / (high - (particle->getExtraInfo(extraInfoName))));
351 } else
352 {
353 return Const::doubleNaN;
354 }
355 };
356 return func;
357 } else {
358 B2FATAL("Wrong number of arguments for meta function transformedNetworkOutput");
359 }
360 }
361
362 Manager::FunctionPtr useBThrustFrame(const std::vector<std::string>& arguments)
363 {
364 if (arguments.size() == 2 || arguments.size() == 3) {
365 auto variableName = arguments[0];
366 std::string mode = arguments[1];
367
368 const bool modeisSignal = mode == "Signal";
369 const bool modeisAuto = mode == "Auto";
370
371 if (not modeisSignal and (mode != "ROE") and not modeisAuto)
372 B2FATAL("Second argument in useBThrustFrame can only be 'Signal', 'ROE' or 'Auto'. Your argument was " + mode);
373
374 const Variable::Manager::Var* var = Manager::Instance().getVariable(variableName);
375
376 std::string maskName = arguments.size() == 3 ? arguments[2] : "";
377
378 auto func = [var, modeisSignal, modeisAuto, maskName](const Particle * particle) -> double {
379 StoreObjPtr<RestOfEvent> roe("RestOfEvent");
380 const Particle* Bparticle = roe->getRelatedFrom<Particle>();
381 RelationVector<ContinuumSuppression> continuumSuppressionRelations = Bparticle->getRelationsTo<ContinuumSuppression>("ALL");
382 ContinuumSuppression* qq = nullptr;
383 if (maskName.empty())
384 {
385 if (continuumSuppressionRelations.size() == 1) {
386 qq = continuumSuppressionRelations[0];
387 } else if (continuumSuppressionRelations.size() > 1) {
388 B2ERROR("The return value of useBThrustFrame is ambiguous. Please provide the mask name as argument.");
389 }
390 } else
391 {
392 qq = Bparticle->getRelatedTo<ContinuumSuppression>(maskName);
393 }
394 if (!qq)
395 return Const::doubleNaN;
396
397 bool isinROE = isInRestOfEvent(particle);
398 ROOT::Math::XYZVector newZ;
399 if (modeisSignal or (modeisAuto and isinROE))
400 newZ = qq->getThrustB();
401 else
402 newZ = qq->getThrustO();
403
404 ROOT::Math::XYZVector newY(0, 0, 0);
405 if (newZ.Z() == 0 and newZ.Y() == 0)
406 newY.SetX(1);
407 else
408 {
409 newY.SetY(newZ.Z());
410 newY.SetZ(-newZ.Y());
411 }
412 ROOT::Math::XYZVector newX = newY.Cross(newZ);
413
414 UseReferenceFrame<CMSRotationFrame> signalframe(newX, newY, newZ);
415
416 return std::get<double>(var->function(particle));
417 };
418 return func;
419 } else {
420 B2FATAL("Wrong number of arguments for meta function useBThrustFrame. It only takes two or three arguments. The first argument must be the variable."
421 "The second can either be 'Signal', 'ROE', or 'Auto'."
422 "The third argument is optional (as long as the ContinuumSuppression was built only once) and can define a specific ROE mask name.");
423 }
424 }
425
426
427 VARIABLE_GROUP("Continuum Suppression");
428 REGISTER_METAVARIABLE("R2(maskname)", R2WithMask, R"DOC(
429Returns reduced Fox-Wolfram R2, defined as ratio of the i-th to the 0-th order Fox Wolfram moments.
430
431.. warning:: You have to run the Continuum Suppression builder module for this variable to be meaningful.
432.. seealso:: :ref:`analysis_continuumsuppression` and `buildContinuumSuppression`.
433)DOC", Manager::VariableDataType::c_double);
434 REGISTER_VARIABLE("R2", R2 , R"DOC(
435Returns reduced Fox-Wolfram R2, defined as ratio of the i-th to the 0-th order Fox Wolfram moments.
436
437.. warning:: You have to run the Continuum Suppression builder module for this variable to be meaningful.
438.. seealso:: :ref:`analysis_continuumsuppression` and `buildContinuumSuppression`.
439:noindex:
440)DOC");
441 REGISTER_METAVARIABLE("thrustBm(maskname)", thrustBmWithMask, R"DOC(
442Returns magnitude of the signal B thrust axis.
443
444.. warning:: You have to run the Continuum Suppression builder module for this variable to be meaningful.
445.. seealso:: :ref:`analysis_continuumsuppression` and `buildContinuumSuppression`.
446)DOC", Manager::VariableDataType::c_double);
447 REGISTER_VARIABLE("thrustBm", thrustBm, R"DOC(
448Returns magnitude of the signal B thrust axis.
449
450.. warning:: You have to run the Continuum Suppression builder module for this variable to be meaningful.
451.. seealso:: :ref:`analysis_continuumsuppression` and `buildContinuumSuppression`.
452:noindex:
453)DOC");
454 REGISTER_METAVARIABLE("thrustOm(maskname)", thrustOmWithMask, R"DOC(
455Returns magnitude of the ROE thrust axis.
456
457.. warning:: You have to run the Continuum Suppression builder module for this variable to be meaningful.
458.. seealso:: :ref:`analysis_continuumsuppression` and `buildContinuumSuppression`.
459)DOC", Manager::VariableDataType::c_double);
460 REGISTER_VARIABLE("thrustOm", thrustOm, R"DOC(
461Returns magnitude of the ROE thrust axis.
462
463.. warning:: You have to run the Continuum Suppression builder module for this variable to be meaningful.
464.. seealso:: :ref:`analysis_continuumsuppression` and `buildContinuumSuppression`.
465:noindex:
466)DOC");
467 REGISTER_METAVARIABLE("cosTBTO(maskname)", cosTBTOWithMask, R"DOC(
468Returns cosine of angle between thrust axis of the signal B and thrust axis of ROE.
469
470.. warning:: You have to run the Continuum Suppression builder module for this variable to be meaningful.
471.. seealso:: :ref:`analysis_continuumsuppression` and `buildContinuumSuppression`.
472)DOC", Manager::VariableDataType::c_double);
473 REGISTER_VARIABLE("cosTBTO", cosTBTO, R"DOC(
474Returns cosine of angle between thrust axis of the signal B and thrust axis of ROE.
475
476.. warning:: You have to run the Continuum Suppression builder module for this variable to be meaningful.
477.. seealso:: :ref:`analysis_continuumsuppression` and `buildContinuumSuppression`.
478:noindex:
479)DOC");
480 REGISTER_METAVARIABLE("cosTBz(maskname)", cosTBzWithMask, R"DOC(
481Returns cosine of angle between thrust axis of the signal B and z-axis.
482
483.. warning:: You have to run the Continuum Suppression builder module for this variable to be meaningful.
484.. seealso:: :ref:`analysis_continuumsuppression` and `buildContinuumSuppression`.
485)DOC", Manager::VariableDataType::c_double);
486 REGISTER_VARIABLE("cosTBz", cosTBz, R"DOC(
487Returns cosine of angle between thrust axis of the signal B and z-axis.
488
489.. warning:: You have to run the Continuum Suppression builder module for this variable to be meaningful.
490.. seealso:: :ref:`analysis_continuumsuppression` and `buildContinuumSuppression`.
491:noindex:
492)DOC");
493 REGISTER_METAVARIABLE("KSFWVariables(variable[, string, string])", KSFWVariables, R"DOC(
494Returns variable et in ``GeV/c``, mm2 in (GeV/c^2)^2, or one of the 16 KSFW moments.
495The second and third arguments are optional unless you have created multiple instances of the ContinuumSuppression with different ROE masks.
496In that case the desired ROE mask name must be provided as well.
497If the second argument is set to 'FS1', the KSFW moment is calculated from the B final state daughters.
498Otherwise, the KSFW moment is calculated from the B primary daughters.
499The ROE mask name is then either the second or the third argument and must not be called 'FS1'.
500Allowed input values for ``variable`` argument are the following:
501
502* mm2, et
503* hso00, hso01, hso02, hso03, hso04
504* hso10, hso12, hso14
505* hso20, hso22, hso24
506* hoo0, hoo1, hoo2, hoo3, hoo4.
507
508.. warning:: You have to run the Continuum Suppression builder module for this variable to be meaningful.
509.. seealso:: :ref:`analysis_continuumsuppression` and `buildContinuumSuppression`.
510)DOC", Manager::VariableDataType::c_double);
511
512 REGISTER_METAVARIABLE("CleoConeCS(integer[, string, string])", CleoConesCS, R"DOC(
513Returns i-th cleo cones from the continuum suppression. The allowed inputs for the ``integer`` argument are integers from *1* to *9*.
514The second and third arguments are optional unless you have created multiple instances of the ContinuumSuppression with different ROE masks.
515In that case the desired ROE mask name must be provided as well.
516If the second argument is set to 'ROE', the CleoCones are calculated only from ROE particles.
517Otherwise, the CleoCones are calculated from all final state particles.
518The ROE mask name is then either the second or the third argument and must not be called 'ROE'. The unit of the CleoConeCS is ``GeV/c``.
519
520.. warning:: You have to run the Continuum Suppression builder module for this variable to be meaningful.
521.. seealso:: :ref:`analysis_continuumsuppression` and `buildContinuumSuppression`.
522)DOC", Manager::VariableDataType::c_double);
523
524 REGISTER_METAVARIABLE("transformedNetworkOutput(name, low, high)", transformedNetworkOutput, R"DOC(
525Transforms the network output :math:`C \to C'` via: :math:`C'=\operatorname{log}((C-\mathrm{low})/(\mathrm{high}-C))`.
526The arguments of the metavariable are the following:
527
528* ``name`` is the `extraInfo` name, where the network output :math:`C` has been stored. If particle is not specified, event `extraInfo` is used instead;
529* ``low``, ``high`` are floating point numbers.
530
531Returns NaN, if the `extraInfo` has not been found.
532)DOC", Manager::VariableDataType::c_double);
533
534 REGISTER_METAVARIABLE("useBThrustFrame(variable, mode)", useBThrustFrame, R"DOC(
535Returns the variable with respect to rotated coordinates, in which z lies on the specified thrust axis.
536If mode is set to ``Signal`` it will use the thrust axis of the reconstructed B candidate, if mode is set to ROE it will use the ROE thrust axis.
537If mode is set to ``Auto`` the function use the thrust axis based on Rest Of Event (ROE) particles.
538Like :b2:var:`isInRestOfEvent`, one has to use this metavariable in ROE loop.
539
540.. warning:: You have to run the Continuum Suppression builder module for this variable to be meaningful.
541.. seealso:: :ref:`analysis_continuumsuppression` and `buildContinuumSuppression`.
542)DOC", Manager::VariableDataType::c_double);
543
544 }
546}
static const double doubleNaN
quiet_NaN
Definition: Const.h:703
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.