19 #include <momemta/ParameterSet.h> 20 #include <momemta/Module.h> 21 #include <momemta/Solution.h> 23 #include <momemta/InputTag.h> 24 #include <momemta/Types.h> 26 #include <Math/GenVector/VectorUtil.h> 81 sqrt_s(parameters.globalParameters().get<
double>(
"energy")) {
82 s12 = get<double>(parameters.get<
InputTag>(
"s12"));
83 s123 = get<double>(parameters.get<
InputTag>(
"s123"));
85 m_p1 = get<LorentzVector>(parameters.get<
InputTag>(
"p1"));
86 m_p2 = get<LorentzVector>(parameters.get<
InputTag>(
"p2"));
87 m_p3 = get<LorentzVector>(parameters.get<
InputTag>(
"p3"));
90 virtual Status
work()
override {
94 const double m1 = m_p1->M();
95 const double sq_m1 =
SQ(m1);
96 const double m3 = m_p3->M();
97 const double sq_m3 =
SQ(m3);
100 if (*s123 >=
SQ(sqrt_s) || *s12 + sq_m3 >= *s123 || sq_m1 >= *s12)
103 const double p3 = m_p3->P();
104 const double E3 = m_p3->E();
105 const double sq_E3 =
SQ(E3);
107 const double c12 = ROOT::Math::VectorUtil::CosTheta(*m_p1, *m_p2);
108 const double c13 = ROOT::Math::VectorUtil::CosTheta(*m_p1, *m_p3);
109 const double c23 = ROOT::Math::VectorUtil::CosTheta(*m_p2, *m_p3);
111 double X = p3 * c23 - E3;
112 double Y = *s123 - *s12 - sq_m3;
114 std::vector<double> abs_p1, abs_p2;
115 solve2Quads(
SQ(X),
SQ(p3 * c13) - sq_E3, 2 * p3 * c13 * X, X * Y, p3 * c13 * Y, 0.25 *
SQ(Y) - sq_E3 * sq_m1,
116 2 * X / E3, 0, 2 * (p3 * c13 / E3 - c12), Y / E3, 0, sq_m1 - *s12,
120 for (std::size_t i = 0; i < abs_p1.size(); i++) {
122 if (abs_p1[i] <= 0 || abs_p2[i] <= 0)
125 const double sin_theta_1 = std::sin(m_p1->Theta());
126 const double sin_theta_2 = std::sin(m_p2->Theta());
127 const double E1 = std::sqrt(
SQ(abs_p1[i]) + sq_m1);
129 LorentzVector p1_sol, p2_sol;
131 abs_p1[i] * std::cos(m_p1->Phi()) * sin_theta_1,
132 abs_p1[i] * std::sin(m_p1->Phi()) * sin_theta_1,
133 abs_p1[i] * std::cos(m_p1->Theta()),
136 abs_p2[i] * std::cos(m_p2->Phi()) * sin_theta_2,
137 abs_p2[i] * std::sin(m_p2->Phi()) * sin_theta_2,
138 abs_p2[i] * std::cos(m_p2->Theta()),
143 LOG(trace) <<
"[SecondaryBlockE] Throwing solution because of invalid mass. " <<
144 "Expected " << m1 <<
", got " << p1_sol.M();
151 LOG(trace) <<
"[SecondaryBlockE] Throwing solution because of invalid mass. " <<
152 "Expected 0., got " << p2_sol.M();
159 LOG(trace) <<
"[SecondaryBlockE] Throwing solution because of invalid invariant mass. " <<
160 "Expected " << *s123 <<
", got " << (p1_sol + p2_sol + *m_p3).M2();
167 LOG(trace) <<
"[SecondaryBlockE] Throwing solution because of invalid invariant mass. " <<
168 "Expected " << *s12 <<
", got " << (p1_sol + p2_sol).M2();
174 const double jacobian = abs_p2[i] *
SQ(abs_p1[i]) * sin_theta_1 * sin_theta_2 /
175 (1024 * std::pow(M_PI, 6) * std::abs(
176 abs_p2[i] * (abs_p1[i] - E1 * c12) * X + (E3 * abs_p1[i] - E1 * p3 * c13) * (E1 - abs_p1[i] * c12) )
179 Solution solution { { p1_sol, p2_sol }, jacobian,
true };
180 solutions->push_back(solution);
183 return (solutions->size() > 0) ? Status::OK : Status::NEXT;
198 std::shared_ptr<SolutionCollection> solutions = produce<SolutionCollection>(
"solutions");
208 .GlobalAttr(
"energy: double");
bool ApproxComparison(double value, double expected)
Compare two doubles and return true if they are approximatively equal.
Generic solution structure representing a set of particles, along with its jacobian.
virtual Status work() override
Main function.
Parent class for all the modules.
Secondary Block E, describing
A class encapsulating a lua table.
Module(PoolPtr pool, const std::string &name)
Constructor.
bool solve2Quads(const double a20, const double a02, const double a11, const double a10, const double a01, const double a00, const double b20, const double b02, const double b11, const double b10, const double b01, const double b00, std::vector< double > &E1, std::vector< double > &E2, bool verbose=false)
Solve a system of two quadratic equations.