14 #include <analysis/variables/EventVariables.h>
16 #include <analysis/VariableManager/Manager.h>
19 #include <framework/datastore/StoreArray.h>
20 #include <framework/datastore/StoreObjPtr.h>
21 #include <framework/dataobjects/EventMetaData.h>
24 #include <analysis/dataobjects/Particle.h>
25 #include <analysis/dataobjects/EventKinematics.h>
27 #include <mdst/dataobjects/MCParticle.h>
28 #include <mdst/dataobjects/Track.h>
29 #include <mdst/dataobjects/ECLCluster.h>
30 #include <mdst/dataobjects/KLMCluster.h>
32 #include <framework/dataobjects/EventT0.h>
35 #include <framework/database/DBObjPtr.h>
36 #include <mdst/dbobjects/BeamSpot.h>
38 #include <analysis/utility/PCmsLabTransform.h>
40 #include <framework/logging/Logger.h>
42 #include <TLorentzVector.h>
53 double eventType(
const Particle*)
55 StoreArray<MCParticle> mcparticles;
56 return (mcparticles.getEntries()) > 0 ? 0 : 1;
59 double isContinuumEvent(
const Particle*)
61 return (isNotContinuumEvent(
nullptr) == 1.0 ? 0.0 : 1.0);
64 double isChargedBEvent(
const Particle*)
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;
74 double isUnmixedBEvent(
const Particle*)
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);
82 if (bPDGs.size() == 2) {
83 return bPDGs[0] * bPDGs[1] < 0;
85 return std::numeric_limits<float>::quiet_NaN();
89 double isNotContinuumEvent(
const Particle*)
91 StoreArray<MCParticle> mcParticles;
92 for (
const MCParticle& mcp : mcParticles) {
93 int pdg_no = mcp.getPDG();
94 if (mcp.getMother() ==
nullptr &&
99 (pdg_no == 9000553) ||
100 (pdg_no == 9010553)))
106 double nMCParticles(
const Particle*)
108 StoreArray<MCParticle> mcps;
109 return mcps.getEntries();
112 double nTracks(
const Particle*)
114 StoreArray<Track> tracks;
115 return tracks.getEntries();
118 double nChargeZeroTrackFits(
const Particle*)
120 StoreArray<TrackFitResult> tfrs;
122 for (
const auto& t : tfrs)
123 if (t.getChargeSign() == 0) out++;
127 double belleECLEnergy(
const Particle*)
129 StoreArray<ECLCluster> eclClusters;
131 for (
int i = 0; i < eclClusters.getEntries(); ++i) {
141 double nKLMClusters(
const Particle*)
143 StoreArray<KLMCluster> klmClusters;
144 return klmClusters.getEntries();
147 double expNum(
const Particle*)
149 StoreObjPtr<EventMetaData> evtMetaData;
150 int exp_no = evtMetaData->getExperiment();
154 double productionIdentifier(
const Particle*)
156 StoreObjPtr<EventMetaData> evtMetaData;
157 int eventProduction = evtMetaData->getProduction();
158 return eventProduction;
161 double evtNum(
const Particle*)
163 StoreObjPtr<EventMetaData> evtMetaData;
164 int evt_no = evtMetaData->getEvent();
168 double runNum(
const Particle*)
170 StoreObjPtr<EventMetaData> evtMetaData;
171 int run_no = evtMetaData->getRun();
176 double getCMSEnergy(
const Particle*)
179 return T.getCMSEnergy();
182 double getBeamPx(
const Particle*)
185 return (T.getBeamFourMomentum()).Px();
188 double getBeamPy(
const Particle*)
191 return (T.getBeamFourMomentum()).Py();
194 double getBeamPz(
const Particle*)
197 return (T.getBeamFourMomentum()).Pz();
200 double getBeamE(
const Particle*)
203 return (T.getBeamFourMomentum()).E();
206 double getIPX(
const Particle*)
208 static DBObjPtr<BeamSpot> beamSpotDB;
209 return (beamSpotDB->getIPPosition()).X();
212 double getIPY(
const Particle*)
214 static DBObjPtr<BeamSpot> beamSpotDB;
215 return (beamSpotDB->getIPPosition()).Y();
218 double getIPZ(
const Particle*)
220 static DBObjPtr<BeamSpot> beamSpotDB;
221 return (beamSpotDB->getIPPosition()).Z();
224 double ipCovMatrixElement(
const Particle*,
const std::vector<double>& element)
226 int elementI = int(std::lround(element[0]));
227 int elementJ = int(std::lround(element[1]));
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();
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();
238 static DBObjPtr<BeamSpot> beamSpotDB;
239 return beamSpotDB->getCovVertex()(elementI, elementJ);
243 double missingMomentumOfEvent(
const Particle*)
245 StoreObjPtr<EventKinematics> evtShape;
247 B2WARNING(
"Cannot find missing momentum information, did you forget to run EventKinematicsModule?");
248 return std::numeric_limits<float>::quiet_NaN();
250 double missing = evtShape->getMissingMomentum().Mag();
254 double missingMomentumOfEvent_Px(
const Particle*)
256 StoreObjPtr<EventKinematics> evtShape;
258 B2WARNING(
"Cannot find missing momentum information, did you forget to run EventKinematicsModule?");
259 return std::numeric_limits<float>::quiet_NaN();
261 double missing = evtShape->getMissingMomentum().Px();
265 double missingMomentumOfEvent_Py(
const Particle*)
267 StoreObjPtr<EventKinematics> evtShape;
269 B2WARNING(
"Cannot find missing momentum information, did you forget to run EventKinematicsModule?");
270 return std::numeric_limits<float>::quiet_NaN();
272 double missing = evtShape->getMissingMomentum().Py();
276 double missingMomentumOfEvent_Pz(
const Particle*)
278 StoreObjPtr<EventKinematics> evtShape;
280 B2WARNING(
"Cannot find missing momentum information, did you forget to run EventKinematicsModule?");
281 return std::numeric_limits<float>::quiet_NaN();
283 double missing = evtShape->getMissingMomentum().Pz();
287 double missingMomentumOfEvent_theta(
const Particle*)
289 StoreObjPtr<EventKinematics> evtShape;
291 B2WARNING(
"Cannot find missing momentum information, did you forget to run EventKinematicsModule?");
292 return std::numeric_limits<float>::quiet_NaN();
294 double missing = evtShape->getMissingMomentum().Theta();
298 double missingMomentumOfEventCMS(
const Particle*)
300 StoreObjPtr<EventKinematics> evtShape;
302 B2WARNING(
"Cannot find missing momentum information, did you forget to run EventKinematicsModule?");
303 return std::numeric_limits<float>::quiet_NaN();
305 double missing = evtShape->getMissingMomentumCMS().Mag();
309 double missingMomentumOfEventCMS_Px(
const Particle*)
311 StoreObjPtr<EventKinematics> evtShape;
313 B2WARNING(
"Cannot find missing momentum information, did you forget to run EventKinematicsModule?");
314 return std::numeric_limits<float>::quiet_NaN();
316 double missing = evtShape->getMissingMomentumCMS().Px();
320 double missingMomentumOfEventCMS_Py(
const Particle*)
322 StoreObjPtr<EventKinematics> evtShape;
324 B2WARNING(
"Cannot find missing momentum information, did you forget to run EventKinematicsModule?");
325 return std::numeric_limits<float>::quiet_NaN();
327 double missing = evtShape->getMissingMomentumCMS().Py();
331 double missingMomentumOfEventCMS_Pz(
const Particle*)
333 StoreObjPtr<EventKinematics> evtShape;
335 B2WARNING(
"Cannot find missing momentum information, did you forget to run EventKinematicsModule?");
336 return std::numeric_limits<float>::quiet_NaN();
338 double missing = evtShape->getMissingMomentumCMS().Pz();
342 double missingMomentumOfEventCMS_theta(
const Particle*)
344 StoreObjPtr<EventKinematics> evtShape;
346 B2WARNING(
"Cannot find missing momentum information, did you forget to run EventKinematicsModule?");
347 return std::numeric_limits<float>::quiet_NaN();
349 double theta = evtShape->getMissingMomentumCMS().Theta();
353 double missingEnergyOfEventCMS(
const Particle*)
355 StoreObjPtr<EventKinematics> evtShape;
357 B2WARNING(
"Cannot find missing momentum information, did you forget to run EventKinematicsModule?");
358 return std::numeric_limits<float>::quiet_NaN();
360 double missing = evtShape->getMissingEnergyCMS();
364 double missingMass2OfEvent(
const Particle*)
366 StoreObjPtr<EventKinematics> evtShape;
368 B2WARNING(
"Cannot find missing momentum information, did you forget to run EventKinematicsModule?");
369 return std::numeric_limits<float>::quiet_NaN();
371 double missing = evtShape->getMissingMass2();
375 double visibleEnergyOfEventCMS(
const Particle*)
377 StoreObjPtr<EventKinematics> evtShape;
379 B2WARNING(
"Cannot find missing momentum information, did you forget to run EventKinematicsModule?");
380 return std::numeric_limits<float>::quiet_NaN();
382 double visible = evtShape->getVisibleEnergyCMS();
386 double totalPhotonsEnergyOfEvent(
const Particle*)
388 StoreObjPtr<EventKinematics> evtShape;
390 B2WARNING(
"Cannot find missing momentum information, did you forget to run EventKinematicsModule?");
391 return std::numeric_limits<float>::quiet_NaN();
393 double energyOfPhotons = evtShape->getTotalPhotonsEnergy();
394 return energyOfPhotons;
397 double eventYearMonthDay(
const Particle*)
399 StoreObjPtr<EventMetaData> evtMetaData;
401 return std::numeric_limits<float>::quiet_NaN();
403 std::time_t rawtime = trunc(evtMetaData->getTime() / 1e9);
404 auto tt = std::gmtime(&rawtime);
405 int y = tt->tm_year + 1900;
406 int m = tt->tm_mon + 1;
408 return (y * 1e4) + (m * 1e2) + d;
411 double eventYear(
const Particle*)
413 StoreObjPtr<EventMetaData> evtMetaData;
415 return std::numeric_limits<float>::quiet_NaN();
417 std::time_t rawtime = trunc(evtMetaData->getTime() / 1e9);
418 auto tt = std::gmtime(&rawtime);
419 return tt->tm_year + 1900;
422 double eventTimeSeconds(
const Particle*)
424 StoreObjPtr<EventMetaData> evtMetaData;
427 return std::numeric_limits<float>::quiet_NaN();
429 double evtTime = trunc(evtMetaData->getTime() / 1e9);
434 double eventTimeSecondsFractionRemainder(
const Particle*)
436 StoreObjPtr<EventMetaData> evtMetaData;
439 return std::numeric_limits<float>::quiet_NaN();
441 double evtTime = trunc(evtMetaData->getTime() / 1e9);
443 double evtTimeFrac = (evtMetaData->getTime() - evtTime * 1e9) / 1e9;
448 double eventT0(
const Particle*)
450 StoreObjPtr<EventT0> evtT0;
453 B2WARNING(
"StoreObjPtr<EventT0> does not exist, are you running over cDST data?");
454 return std::numeric_limits<float>::quiet_NaN();
457 if (evtT0->hasEventT0()) {
458 return evtT0->getEventT0();
460 return std::numeric_limits<float>::quiet_NaN();
465 VARIABLE_GROUP(
"Event");
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");
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");
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");
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");
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)");
507 REGISTER_VARIABLE("IPX", getIPX, R"DOC(
508 [Eventbased] x coordinate of the measured interaction point.
510 .. note:: For old data and uncalibrated MC files this will return 0.0.
512 .. note:: You might hear tracking and calibration people refer to this as the ``BeamSpot``.
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")
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.");
531 VARIABLE_GROUP("EventKinematics");
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");
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.)");