Belle II Software  release-05-01-25
EventVariables.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2015-2020 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Marko Staric, Anze Zupanc, Thomas Keck, Sam Cunliffe *
7  * for the EventKinematics variables: Ami Rostomyan, *
8  * Michel Villanueva *
9  * *
10  * This software is provided "as is" without any warranty. *
11  **************************************************************************/
12 
13 // Own include
14 #include <analysis/variables/EventVariables.h>
15 
16 #include <analysis/VariableManager/Manager.h>
17 
18 // framework - DataStore
19 #include <framework/datastore/StoreArray.h>
20 #include <framework/datastore/StoreObjPtr.h>
21 #include <framework/dataobjects/EventMetaData.h>
22 
23 // dataobjects
24 #include <analysis/dataobjects/Particle.h>
25 #include <analysis/dataobjects/EventKinematics.h>
26 
27 #include <mdst/dataobjects/MCParticle.h>
28 #include <mdst/dataobjects/Track.h>
29 #include <mdst/dataobjects/ECLCluster.h>
30 #include <mdst/dataobjects/KLMCluster.h>
31 
32 #include <framework/dataobjects/EventT0.h>
33 
34 // database
35 #include <framework/database/DBObjPtr.h>
36 #include <mdst/dbobjects/BeamSpot.h>
37 
38 #include <analysis/utility/PCmsLabTransform.h>
39 
40 #include <framework/logging/Logger.h>
41 
42 #include <TLorentzVector.h>
43 #include <TVector3.h>
44 
45 namespace Belle2 {
50  namespace Variable {
51 
52  // Event ------------------------------------------------
53  double eventType(const Particle*)
54  {
55  StoreArray<MCParticle> mcparticles;
56  return (mcparticles.getEntries()) > 0 ? 0 : 1;
57  }
58 
59  double isContinuumEvent(const Particle*)
60  {
61  return (isNotContinuumEvent(nullptr) == 1.0 ? 0.0 : 1.0);
62  }
63 
64  double isChargedBEvent(const Particle*)
65  {
66  StoreArray<MCParticle> mcParticles;
67  for (const auto& mcp : mcParticles) {
68  int pdg_no = mcp.getPDG();
69  if (abs(pdg_no) == 521) return 1.0;
70  }
71  return 0.0;
72  }
73 
74  double isUnmixedBEvent(const Particle*)
75  {
76  StoreArray<MCParticle> mcParticles;
77  std::vector<int> bPDGs;
78  for (const auto& mcp : mcParticles) {
79  int pdg_no = mcp.getPDG();
80  if (abs(pdg_no) == 511) bPDGs.push_back(pdg_no);
81  }
82  if (bPDGs.size() == 2) {
83  return bPDGs[0] * bPDGs[1] < 0;
84  }
85  return std::numeric_limits<float>::quiet_NaN();
86  }
87 
88 
89  double isNotContinuumEvent(const Particle*)
90  {
91  StoreArray<MCParticle> mcParticles;
92  for (const MCParticle& mcp : mcParticles) {
93  int pdg_no = mcp.getPDG();
94  if (mcp.getMother() == nullptr &&
95  ((pdg_no == 553) ||
96  (pdg_no == 100553) ||
97  (pdg_no == 200553) ||
98  (pdg_no == 300553) ||
99  (pdg_no == 9000553) ||
100  (pdg_no == 9010553)))
101  return 1.0;
102  }
103  return 0.0;
104  }
105 
106  double nMCParticles(const Particle*)
107  {
108  StoreArray<MCParticle> mcps;
109  return mcps.getEntries();
110  }
111 
112  double nTracks(const Particle*)
113  {
114  StoreArray<Track> tracks;
115  return tracks.getEntries();
116  }
117 
118  double nChargeZeroTrackFits(const Particle*)
119  {
120  StoreArray<TrackFitResult> tfrs;
121  int out = 0;
122  for (const auto& t : tfrs)
123  if (t.getChargeSign() == 0) out++;
124  return double(out);
125  }
126 
127  double belleECLEnergy(const Particle*)
128  {
129  StoreArray<ECLCluster> eclClusters;
130  double result = 0;
131  for (int i = 0; i < eclClusters.getEntries(); ++i) {
132  // sum only ECLClusters which have the N1 (n photons) hypothesis
133  if (!eclClusters[i]->hasHypothesis(ECLCluster::EHypothesisBit::c_nPhotons))
134  continue;
135 
136  result += eclClusters[i]->getEnergy(ECLCluster::EHypothesisBit::c_nPhotons);
137  }
138  return result;
139  }
140 
141  double nKLMClusters(const Particle*)
142  {
143  StoreArray<KLMCluster> klmClusters;
144  return klmClusters.getEntries();
145  }
146 
147  double expNum(const Particle*)
148  {
149  StoreObjPtr<EventMetaData> evtMetaData;
150  int exp_no = evtMetaData->getExperiment();
151  return exp_no;
152  }
153 
154  double productionIdentifier(const Particle*)
155  {
156  StoreObjPtr<EventMetaData> evtMetaData;
157  int eventProduction = evtMetaData->getProduction();
158  return eventProduction;
159  }
160 
161  double evtNum(const Particle*)
162  {
163  StoreObjPtr<EventMetaData> evtMetaData;
164  int evt_no = evtMetaData->getEvent();
165  return evt_no;
166  }
167 
168  double runNum(const Particle*)
169  {
170  StoreObjPtr<EventMetaData> evtMetaData;
171  int run_no = evtMetaData->getRun();
172  return run_no;
173  }
174 
175  // Beam Energies
176  double getCMSEnergy(const Particle*)
177  {
178  PCmsLabTransform T;
179  return T.getCMSEnergy();
180  }
181 
182  double getBeamPx(const Particle*)
183  {
184  PCmsLabTransform T;
185  return (T.getBeamFourMomentum()).Px();
186  }
187 
188  double getBeamPy(const Particle*)
189  {
190  PCmsLabTransform T;
191  return (T.getBeamFourMomentum()).Py();
192  }
193 
194  double getBeamPz(const Particle*)
195  {
196  PCmsLabTransform T;
197  return (T.getBeamFourMomentum()).Pz();
198  }
199 
200  double getBeamE(const Particle*)
201  {
202  PCmsLabTransform T;
203  return (T.getBeamFourMomentum()).E();
204  }
205 
206  double getIPX(const Particle*)
207  {
208  static DBObjPtr<BeamSpot> beamSpotDB;
209  return (beamSpotDB->getIPPosition()).X();
210  }
211 
212  double getIPY(const Particle*)
213  {
214  static DBObjPtr<BeamSpot> beamSpotDB;
215  return (beamSpotDB->getIPPosition()).Y();
216  }
217 
218  double getIPZ(const Particle*)
219  {
220  static DBObjPtr<BeamSpot> beamSpotDB;
221  return (beamSpotDB->getIPPosition()).Z();
222  }
223 
224  double ipCovMatrixElement(const Particle*, const std::vector<double>& element)
225  {
226  int elementI = int(std::lround(element[0]));
227  int elementJ = int(std::lround(element[1]));
228 
229  if (elementI < 0 || elementI > 3) {
230  B2WARNING("Requested IP covariance matrix element is out of boundaries [0 - 3]: i = " << elementI);
231  return std::numeric_limits<float>::quiet_NaN();
232  }
233  if (elementJ < 0 || elementJ > 3) {
234  B2WARNING("Requested particle's momentumVertex covariance matrix element is out of boundaries [0 - 3]: j = " << elementJ);
235  return std::numeric_limits<float>::quiet_NaN();
236  }
237 
238  static DBObjPtr<BeamSpot> beamSpotDB;
239  return beamSpotDB->getCovVertex()(elementI, elementJ);
240  }
241 
242  // Event kinematics -> missing momentum in lab and CMS, missing energy and mass2, visible energy
243  double missingMomentumOfEvent(const Particle*)
244  {
245  StoreObjPtr<EventKinematics> evtShape;
246  if (!evtShape) {
247  B2WARNING("Cannot find missing momentum information, did you forget to run EventKinematicsModule?");
248  return std::numeric_limits<float>::quiet_NaN();
249  }
250  double missing = evtShape->getMissingMomentum().Mag();
251  return missing;
252  }
253 
254  double missingMomentumOfEvent_Px(const Particle*)
255  {
256  StoreObjPtr<EventKinematics> evtShape;
257  if (!evtShape) {
258  B2WARNING("Cannot find missing momentum information, did you forget to run EventKinematicsModule?");
259  return std::numeric_limits<float>::quiet_NaN();
260  }
261  double missing = evtShape->getMissingMomentum().Px();
262  return missing;
263  }
264 
265  double missingMomentumOfEvent_Py(const Particle*)
266  {
267  StoreObjPtr<EventKinematics> evtShape;
268  if (!evtShape) {
269  B2WARNING("Cannot find missing momentum information, did you forget to run EventKinematicsModule?");
270  return std::numeric_limits<float>::quiet_NaN();
271  }
272  double missing = evtShape->getMissingMomentum().Py();
273  return missing;
274  }
275 
276  double missingMomentumOfEvent_Pz(const Particle*)
277  {
278  StoreObjPtr<EventKinematics> evtShape;
279  if (!evtShape) {
280  B2WARNING("Cannot find missing momentum information, did you forget to run EventKinematicsModule?");
281  return std::numeric_limits<float>::quiet_NaN();
282  }
283  double missing = evtShape->getMissingMomentum().Pz();
284  return missing;
285  }
286 
287  double missingMomentumOfEvent_theta(const Particle*)
288  {
289  StoreObjPtr<EventKinematics> evtShape;
290  if (!evtShape) {
291  B2WARNING("Cannot find missing momentum information, did you forget to run EventKinematicsModule?");
292  return std::numeric_limits<float>::quiet_NaN();
293  }
294  double missing = evtShape->getMissingMomentum().Theta();
295  return missing;
296  }
297 
298  double missingMomentumOfEventCMS(const Particle*)
299  {
300  StoreObjPtr<EventKinematics> evtShape;
301  if (!evtShape) {
302  B2WARNING("Cannot find missing momentum information, did you forget to run EventKinematicsModule?");
303  return std::numeric_limits<float>::quiet_NaN();
304  }
305  double missing = evtShape->getMissingMomentumCMS().Mag();
306  return missing;
307  }
308 
309  double missingMomentumOfEventCMS_Px(const Particle*)
310  {
311  StoreObjPtr<EventKinematics> evtShape;
312  if (!evtShape) {
313  B2WARNING("Cannot find missing momentum information, did you forget to run EventKinematicsModule?");
314  return std::numeric_limits<float>::quiet_NaN();
315  }
316  double missing = evtShape->getMissingMomentumCMS().Px();
317  return missing;
318  }
319 
320  double missingMomentumOfEventCMS_Py(const Particle*)
321  {
322  StoreObjPtr<EventKinematics> evtShape;
323  if (!evtShape) {
324  B2WARNING("Cannot find missing momentum information, did you forget to run EventKinematicsModule?");
325  return std::numeric_limits<float>::quiet_NaN();
326  }
327  double missing = evtShape->getMissingMomentumCMS().Py();
328  return missing;
329  }
330 
331  double missingMomentumOfEventCMS_Pz(const Particle*)
332  {
333  StoreObjPtr<EventKinematics> evtShape;
334  if (!evtShape) {
335  B2WARNING("Cannot find missing momentum information, did you forget to run EventKinematicsModule?");
336  return std::numeric_limits<float>::quiet_NaN();
337  }
338  double missing = evtShape->getMissingMomentumCMS().Pz();
339  return missing;
340  }
341 
342  double missingMomentumOfEventCMS_theta(const Particle*)
343  {
344  StoreObjPtr<EventKinematics> evtShape;
345  if (!evtShape) {
346  B2WARNING("Cannot find missing momentum information, did you forget to run EventKinematicsModule?");
347  return std::numeric_limits<float>::quiet_NaN();
348  }
349  double theta = evtShape->getMissingMomentumCMS().Theta();
350  return theta;
351  }
352 
353  double missingEnergyOfEventCMS(const Particle*)
354  {
355  StoreObjPtr<EventKinematics> evtShape;
356  if (!evtShape) {
357  B2WARNING("Cannot find missing momentum information, did you forget to run EventKinematicsModule?");
358  return std::numeric_limits<float>::quiet_NaN();
359  }
360  double missing = evtShape->getMissingEnergyCMS();
361  return missing;
362  }
363 
364  double missingMass2OfEvent(const Particle*)
365  {
366  StoreObjPtr<EventKinematics> evtShape;
367  if (!evtShape) {
368  B2WARNING("Cannot find missing momentum information, did you forget to run EventKinematicsModule?");
369  return std::numeric_limits<float>::quiet_NaN();
370  }
371  double missing = evtShape->getMissingMass2();
372  return missing;
373  }
374 
375  double visibleEnergyOfEventCMS(const Particle*)
376  {
377  StoreObjPtr<EventKinematics> evtShape;
378  if (!evtShape) {
379  B2WARNING("Cannot find missing momentum information, did you forget to run EventKinematicsModule?");
380  return std::numeric_limits<float>::quiet_NaN();
381  }
382  double visible = evtShape->getVisibleEnergyCMS();
383  return visible;
384  }
385 
386  double totalPhotonsEnergyOfEvent(const Particle*)
387  {
388  StoreObjPtr<EventKinematics> evtShape;
389  if (!evtShape) {
390  B2WARNING("Cannot find missing momentum information, did you forget to run EventKinematicsModule?");
391  return std::numeric_limits<float>::quiet_NaN();
392  }
393  double energyOfPhotons = evtShape->getTotalPhotonsEnergy();
394  return energyOfPhotons;
395  }
396 
397  double eventYearMonthDay(const Particle*)
398  {
399  StoreObjPtr<EventMetaData> evtMetaData;
400  if (!evtMetaData) {
401  return std::numeric_limits<float>::quiet_NaN();
402  }
403  std::time_t rawtime = trunc(evtMetaData->getTime() / 1e9);
404  auto tt = std::gmtime(&rawtime); // GMT
405  int y = tt->tm_year + 1900; // years since 1900
406  int m = tt->tm_mon + 1; // months since January
407  int d = tt->tm_mday; // day of the month
408  return (y * 1e4) + (m * 1e2) + d;
409  }
410 
411  double eventYear(const Particle*)
412  {
413  StoreObjPtr<EventMetaData> evtMetaData;
414  if (!evtMetaData) {
415  return std::numeric_limits<float>::quiet_NaN();
416  }
417  std::time_t rawtime = trunc(evtMetaData->getTime() / 1e9);
418  auto tt = std::gmtime(&rawtime);
419  return tt->tm_year + 1900;
420  }
421 
422  double eventTimeSeconds(const Particle*)
423  {
424  StoreObjPtr<EventMetaData> evtMetaData;
425 
426  if (!evtMetaData) {
427  return std::numeric_limits<float>::quiet_NaN();
428  }
429  double evtTime = trunc(evtMetaData->getTime() / 1e9);
430 
431  return evtTime;
432  }
433 
434  double eventTimeSecondsFractionRemainder(const Particle*)
435  {
436  StoreObjPtr<EventMetaData> evtMetaData;
437 
438  if (!evtMetaData) {
439  return std::numeric_limits<float>::quiet_NaN();
440  }
441  double evtTime = trunc(evtMetaData->getTime() / 1e9);
442 
443  double evtTimeFrac = (evtMetaData->getTime() - evtTime * 1e9) / 1e9;
444 
445  return evtTimeFrac;
446  }
447 
448  double eventT0(const Particle*)
449  {
450  StoreObjPtr<EventT0> evtT0;
451 
452  if (!evtT0) {
453  B2WARNING("StoreObjPtr<EventT0> does not exist, are you running over cDST data?");
454  return std::numeric_limits<float>::quiet_NaN();
455  }
456 
457  if (evtT0->hasEventT0()) {
458  return evtT0->getEventT0();
459  } else {
460  return std::numeric_limits<float>::quiet_NaN();
461  }
462  }
463 
464 
465  VARIABLE_GROUP("Event");
466 
467  REGISTER_VARIABLE("EventType", eventType, "[Eventbased] EventType (0 MC, 1 Data)");
468  REGISTER_VARIABLE("isContinuumEvent", isContinuumEvent,
469  "[Eventbased] true if event doesn't contain an Y(4S)");
470  REGISTER_VARIABLE("isNotContinuumEvent", isNotContinuumEvent,
471  "[Eventbased] 1.0 if event does contain an Y(4S) and therefore is not a continuum Event");
472 
473  REGISTER_VARIABLE("isChargedBEvent", isChargedBEvent,
474  "[Eventbased] true if event contains a charged B-meson");
475  REGISTER_VARIABLE("isUnmixedBEvent", isUnmixedBEvent,
476  R"DOC([Eventbased] true if event contains opposite flavor neutral B-mesons,
477 false in case of same flavor B-mesons and NaN if an event has no generated neutral B)DOC");
478 
479  REGISTER_VARIABLE("nTracks", nTracks,
480  "[Eventbased] number of tracks in the event");
481  REGISTER_VARIABLE("nChargeZeroTrackFits", nChargeZeroTrackFits,
482  "[Eventbased] number of track fits with a zero charge."
483  "Sometimes this can happen if background or non IP originating "
484  "tracks (for example) are fit from the IP. These tracks are "
485  "removed from particle lists but a large number charge zero "
486  "fits them may indicate problems with whole event constraints "
487  "or abnominally high beam backgrounds and/or noisy events.")
488  REGISTER_VARIABLE("belleECLEnergy", belleECLEnergy,
489  "[Eventbased] legacy total energy in ECL in the event as used in Belle 1 analyses. For Belle II "
490  "consider totalEnergyOfParticlesInList(gamma:all) instead");
491  REGISTER_VARIABLE("nKLMClusters", nKLMClusters,
492  "[Eventbased] number of KLM in the event");
493  REGISTER_VARIABLE("nMCParticles", nMCParticles,
494  "[Eventbased] number of MCParticles in the event");
495 
496  REGISTER_VARIABLE("expNum", expNum, "[Eventbased] experiment number");
497  REGISTER_VARIABLE("evtNum", evtNum, "[Eventbased] event number");
498  REGISTER_VARIABLE("runNum", runNum, "[Eventbased] run number");
499  REGISTER_VARIABLE("productionIdentifier", productionIdentifier, "[Eventbased] production identifier");
500 
501  REGISTER_VARIABLE("Ecms", getCMSEnergy, "[Eventbased] CMS energy");
502  REGISTER_VARIABLE("beamE", getBeamE, "[Eventbased] Beam energy (lab)");
503  REGISTER_VARIABLE("beamPx", getBeamPx, "[Eventbased] Beam momentum Px (lab)");
504  REGISTER_VARIABLE("beamPy", getBeamPy, "[Eventbased] Beam momentum Py (lab)");
505  REGISTER_VARIABLE("beamPz", getBeamPz, "[Eventbased] Beam momentum Pz (lab)");
506 
507  REGISTER_VARIABLE("IPX", getIPX, R"DOC(
508 [Eventbased] x coordinate of the measured interaction point.
509 
510 .. note:: For old data and uncalibrated MC files this will return 0.0.
511 
512 .. note:: You might hear tracking and calibration people refer to this as the ``BeamSpot``.
513 )DOC");
514  REGISTER_VARIABLE("IPY", getIPY, "[Eventbased] y coordinate of the measured interaction point");
515  REGISTER_VARIABLE("IPZ", getIPZ, "[Eventbased] z coordinate of the measured interaction point");
516  REGISTER_VARIABLE("IPCov(i,j)", ipCovMatrixElement, "[Eventbased] (i,j)-th element of the covariance matrix of the measured interaction point")
517 
518  REGISTER_VARIABLE("date", eventYearMonthDay,
519  "[Eventbased] Returns the date when the event was recorded, a number of the form YYYYMMDD (in UTC).\n"
520  " See also eventYear, provided for convenience."
521  " For more precise eventTime, see eventTimeSeconds and eventTimeSecondsFractionRemainder.");
522  REGISTER_VARIABLE("year", eventYear,
523  "[Eventbased] Returns the year when the event was recorded (in UTC).\n"
524  "For more precise eventTime, see eventTimeSeconds and eventTimeSecondsFractionRemainder.");
525  REGISTER_VARIABLE("eventTimeSeconds", eventTimeSeconds,
526  "[Eventbased] Time of the event in seconds (truncated down) since 1970/1/1 (Unix epoch).");
527  REGISTER_VARIABLE("eventTimeSecondsFractionRemainder", eventTimeSecondsFractionRemainder,
528  "[Eventbased] Remainder of the event time in fractions of a second.\n"
529  "Use eventTimeSeconds + eventTimeSecondsFractionRemainder to get the total event time in seconds.");
530 
531  VARIABLE_GROUP("EventKinematics");
532 
533  REGISTER_VARIABLE("missingMomentumOfEvent", missingMomentumOfEvent,
534  "[Eventbased] The magnitude of the missing momentum in lab obtained with EventKinematics module")
535  REGISTER_VARIABLE("missingMomentumOfEvent_Px", missingMomentumOfEvent_Px,
536  "[Eventbased] The x component of the missing momentum in lab obtained with EventKinematics module")
537  REGISTER_VARIABLE("missingMomentumOfEvent_Py", missingMomentumOfEvent_Py,
538  "[Eventbased] The y component of the missing momentum in lab obtained with EventKinematics module")
539  REGISTER_VARIABLE("missingMomentumOfEvent_Pz", missingMomentumOfEvent_Pz,
540  "[Eventbased] The z component of the missing momentum in lab obtained with EventKinematics module")
541  REGISTER_VARIABLE("missingMomentumOfEvent_theta", missingMomentumOfEvent_theta,
542  "[Eventbased] The theta angle of the missing momentum of the event in lab obtained with EventKinematics module")
543  REGISTER_VARIABLE("missingMomentumOfEventCMS", missingMomentumOfEventCMS,
544  "[Eventbased] The magnitude of the missing momentum in CMS obtained with EventKinematics module")
545  REGISTER_VARIABLE("missingMomentumOfEventCMS_Px", missingMomentumOfEventCMS_Px,
546  "[Eventbased] The x component of the missing momentum in CMS obtained with EventKinematics module")
547  REGISTER_VARIABLE("missingMomentumOfEventCMS_Py", missingMomentumOfEventCMS_Py,
548  "[Eventbased] The y component of the missing momentum in CMS obtained with EventKinematics module")
549  REGISTER_VARIABLE("missingMomentumOfEventCMS_Pz", missingMomentumOfEventCMS_Pz,
550  "[Eventbased] The z component of the missing momentum in CMS obtained with EventKinematics module")
551  REGISTER_VARIABLE("missingMomentumOfEventCMS_theta", missingMomentumOfEventCMS_theta,
552  "[Eventbased] The theta angle of the missing momentum in CMS obtained with EventKinematics module")
553  REGISTER_VARIABLE("missingEnergyOfEventCMS", missingEnergyOfEventCMS,
554  "[Eventbased] The missing energy in CMS obtained with EventKinematics module")
555  REGISTER_VARIABLE("missingMass2OfEvent", missingMass2OfEvent,
556  "[Eventbased] The missing mass squared obtained with EventKinematics module")
557  REGISTER_VARIABLE("visibleEnergyOfEventCMS", visibleEnergyOfEventCMS,
558  "[Eventbased] The visible energy in CMS obtained with EventKinematics module")
559  REGISTER_VARIABLE("totalPhotonsEnergyOfEvent", totalPhotonsEnergyOfEvent,
560  "[Eventbased] The energy in lab of all the photons obtained with EventKinematics module");
561 
562  VARIABLE_GROUP("Event (cDST only)");
563  REGISTER_VARIABLE("eventT0", eventT0,
564  "[Eventbased][Calibration] The Event t0, measured in ns, is the time of the event relative to the\n"
565  "trigger time. The event time can be measured by several sub-detectors including the CDC, ECL, and TOP.\n"
566  "This Event t0 variable is the final combined value of all the event time measurements.\n"
567  "(Currently only the CDC and ECL are used in this combination.)");
568  }
570 }
Belle2::ECLCluster::EHypothesisBit::c_nPhotons
@ c_nPhotons
CR is split into n photons (N1)
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19