ptmax-ptmin module

cpp/src/AnalysisVariables.cc

+/* Belle2 headers. */
+#include <analysis/dataobjects/Particle.h>
+#include <analysis/utility/PCmsLabTransform.h>
+#include <analysis/utility/ReferenceFrame.h>
+#include <analysis/VariableManager/Manager.h>
+#include <framework/core/Module.h>
+
+/* ROOT headers. */
+#include <Math/Vector4D.h>
+#include <TMath.h>
+
+/* C++ headers. */
+#include <algorithm>
+#include <cmath>
+
+namespace Belle2 {
+
+  namespace Variable {
+
+    double ptmax(const Particle* particle)
+    {
+      if (!particle)
+        return std::numeric_limits<double>::quiet_NaN();
+      const auto& frame{ReferenceFrame::GetCurrent()};
+      const ROOT::Math::PxPyPzEVector p4_mu0{frame.getMomentum(particle->getDaughter(0)->getDaughter(0))};
+      const ROOT::Math::PxPyPzEVector p4_mu1{frame.getMomentum(particle->getDaughter(0)->getDaughter(1))};
+      const ROOT::Math::PxPyPzEVector p4_rec{frame.getMomentum(particle)};
+      const double plong_max{
+        p4_mu0.P() > p4_mu1.P() ?
+        (p4_mu0.X() * p4_rec.X() + p4_mu0.Y() * p4_rec.Y() + p4_mu0.Z() * p4_rec.Z()) / p4_mu0.P() :
+        (p4_mu1.X() * p4_rec.X() + p4_mu1.Y() * p4_rec.Y() + p4_mu1.Z() * p4_rec.Z()) / p4_mu1.P()
+      };
+      return std::sqrt(p4_rec.P() * p4_rec.P() - plong_max * plong_max);
+    }
+
+    double plongmax(const Particle* particle)
+    {
+      if (!particle)
+        return std::numeric_limits<double>::quiet_NaN();
+      const auto& frame{ReferenceFrame::GetCurrent()};
+      const ROOT::Math::PxPyPzEVector p4_mu0{frame.getMomentum(particle->getDaughter(0)->getDaughter(0))};
+      const ROOT::Math::PxPyPzEVector p4_mu1{frame.getMomentum(particle->getDaughter(0)->getDaughter(1))};
+      const ROOT::Math::PxPyPzEVector p4_rec{frame.getMomentum(particle)};
+      return p4_mu0.P() > p4_mu1.P() ?
+             (p4_mu0.X() * p4_rec.X() + p4_mu0.Y() * p4_rec.Y() + p4_mu0.Z() * p4_rec.Z()) / p4_mu0.P() :
+             (p4_mu1.X() * p4_rec.X() + p4_mu1.Y() * p4_rec.Y() + p4_mu1.Z() * p4_rec.Z()) / p4_mu1.P();
+    }
+
+    double ptmin(const Particle* particle)
+    {
+      if (!particle)
+        return std::numeric_limits<double>::quiet_NaN();
+      const auto& frame{ReferenceFrame::GetCurrent()};
+      const ROOT::Math::PxPyPzEVector p4_mu0{frame.getMomentum(particle->getDaughter(0)->getDaughter(0))};
+      const ROOT::Math::PxPyPzEVector p4_mu1{frame.getMomentum(particle->getDaughter(0)->getDaughter(1))};
+      const ROOT::Math::PxPyPzEVector p4_rec{frame.getMomentum(particle)};
+      const double plong_min{
+        p4_mu0.P() < p4_mu1.P() ?
+        (p4_mu0.X() * p4_rec.X() + p4_mu0.Y() * p4_rec.Y() + p4_mu0.Z() * p4_rec.Z()) / p4_mu0.P() :
+        (p4_mu1.X() * p4_rec.X() + p4_mu1.Y() * p4_rec.Y() + p4_mu1.Z() * p4_rec.Z()) / p4_mu1.P()
+      };
+      return std::sqrt(p4_rec.P() * p4_rec.P() - plong_min * plong_min);
+    }
+
+    double plongmin(const Particle* particle)
+    {
+      if (!particle)
+        return std::numeric_limits<double>::quiet_NaN();
+      const auto& frame{ReferenceFrame::GetCurrent()};
+      const ROOT::Math::PxPyPzEVector p4_mu0{frame.getMomentum(particle->getDaughter(0)->getDaughter(0))};
+      const ROOT::Math::PxPyPzEVector p4_mu1{frame.getMomentum(particle->getDaughter(0)->getDaughter(1))};
+      const ROOT::Math::PxPyPzEVector p4_rec{frame.getMomentum(particle)};
+      return p4_mu0.P() < p4_mu1.P() ?
+             (p4_mu0.X() * p4_rec.X() + p4_mu0.Y() * p4_rec.Y() + p4_mu0.Z() * p4_rec.Z()) / p4_mu0.P() :
+             (p4_mu1.X() * p4_rec.X() + p4_mu1.Y() * p4_rec.Y() + p4_mu1.Z() * p4_rec.Z()) / p4_mu1.P();
+    }
+
+    double deltaPhiL1(const Particle* particle)
+    {
+      if (!particle)
+        return std::numeric_limits<double>::quiet_NaN();
+      const auto& frame{ReferenceFrame::GetCurrent()};
+      const double phi_mu0{(frame.getMomentum(particle->getDaughter(0)->getDaughter(0))).Phi()};
+      const double phi_mu1{(frame.getMomentum(particle->getDaughter(0)->getDaughter(1))).Phi()};
+      return std::min(std::max(phi_mu0, phi_mu1) - std::min(phi_mu0, phi_mu1),
+                      std::min(phi_mu0, phi_mu1) - std::max(phi_mu0, phi_mu1) + 2 * TMath::Pi());
+    }
+
+    VARIABLE_GROUP("Analysis-specific variables");
+    REGISTER_VARIABLE("ptmax", ptmax,
+                      "This is ptmax.");
+    REGISTER_VARIABLE("ptmin", ptmin,
+                      "This is ptmin.");
+    REGISTER_VARIABLE("plong_min", plongmin,
+                      "This is plong_min.");
+    REGISTER_VARIABLE("plong_max", plongmax,
+                      "This is plot_max.");
+    REGISTER_VARIABLE("deltaPhiL1", deltaPhiL1,
+                      "This is deltaPhiL1.");
+  }
+
+  class CustomVariablesModule : public Module {};
+  REG_MODULE(CustomVariables);
+}