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> 83 sqrt_s(parameters.globalParameters().get<
double>(
"energy")) {
84 s12 = get<double>(parameters.get<
InputTag>(
"s12"));
85 s123 = get<double>(parameters.get<
InputTag>(
"s123"));
86 s1234 = get<double>(parameters.get<
InputTag>(
"s1234"));
88 m1 = parameters.get<
double>(
"m1", 0.);
90 m_p2 = get<LorentzVector>(parameters.get<
InputTag>(
"p2"));
91 m_p3 = get<LorentzVector>(parameters.get<
InputTag>(
"p3"));
92 m_p4 = get<LorentzVector>(parameters.get<
InputTag>(
"p4"));
95 virtual Status
work()
override {
99 const double sq_m1 =
SQ(m1);
100 const double m2 = m_p2->M();
101 const double sq_m2 =
SQ(m2);
102 const double m3 = m_p3->M();
103 const double sq_m3 =
SQ(m3);
104 const double m4 = m_p4->M();
105 const double sq_m4 =
SQ(m4);
108 if (*s1234 >=
SQ(sqrt_s) || *s12 + sq_m3 >= *s123 || *s123 + sq_m4 >= *s1234 || sq_m1 + sq_m2 >= *s12)
111 const double p2x = m_p2->Px();
112 const double p2y = m_p2->Py();
113 const double p2z = m_p2->Pz();
114 const double E2 = m_p2->E();
116 const double p3x = m_p3->Px();
117 const double p3y = m_p3->Py();
118 const double p3z = m_p3->Pz();
119 const double E3 = m_p3->E();
121 const double p4x = m_p4->Px();
122 const double p4y = m_p4->Py();
123 const double p4z = m_p4->Pz();
124 const double E4 = m_p4->E();
126 const double p2p3 = m_p2->Dot(*m_p3);
127 const double p2p4 = m_p2->Dot(*m_p4);
128 const double p3p4 = m_p3->Dot(*m_p4);
138 const double a11 = p2x;
139 const double a12 = p2y;
140 const double a13 = p2z;
142 const double a21 = p2x + p3x;
143 const double a22 = p2y + p3y;
144 const double a23 = p2z + p3z;
146 const double a31 = p2x + p3x + p4x;
147 const double a32 = p2y + p3y + p4y;
148 const double a33 = p2z + p3z + p4z;
150 const double b1 = E2;
151 const double c1 = 0.5 * (sq_m1 + sq_m2 - *s12);
152 const double b2 = E2 + E3;
153 const double c2 = 0.5 * (sq_m1 + sq_m2 + sq_m3 - *s123) + p2p3;
154 const double b3 = E2 + E3 + E4;
155 const double c3 = 0.5 * (sq_m1 + sq_m2 + sq_m3 + sq_m4 - *s1234) + p2p3 + p3p4 + p2p4;
157 const double inv_det = 1 / ((a13 * a22 - a12 * a23) * a31 - (a13 * a21 - a11 * a23) * a32 + (a12 * a21 - a11 * a22) * a33);
159 const double Ax = ((a13 * a22 - a12 * a23) * b3 - (a13 * a32 - a12 * a33) * b2 + (a23 * a32 - a22 * a33) * b1) * inv_det;
160 const double Bx = ((a13 * a22 - a12 * a23) * c3 - (a13 * a32 - a12 * a33) * c2 + (a23 * a32 - a22 * a33) * c1) * inv_det;
162 const double Ay = - ((a13 * a21 - a11 * a23) * b3 - (a13 * a31 - a11 * a33) * b2 + (a23 * a31 - a21 * a33) * b1) * inv_det;
163 const double By = - ((a13 * a21 - a11 * a23) * c3 - (a13 * a31 - a11 * a33) * c2 + (a23 * a31 - a21 * a33) * c1) * inv_det;
165 const double Az = ((a12 * a21 - a11 * a22) * b3 - (a12 * a31 - a11 * a32) * b2 + (a22 * a31 - a21 * a32) * b1) * inv_det;
166 const double Bz = ((a12 * a21 - a11 * a22) * c3 - (a12 * a31 - a11 * a32) * c2 + (a22 * a31 - a21 * a32) * c1) * inv_det;
169 std::vector<double> E1_sol;
170 bool foundSolution =
solveQuadratic(
SQ(Ax) +
SQ(Ay) +
SQ(Az) - 1, 2 * (Ax * Bx + Ay * By + Az * Bz),
SQ(Bx) +
SQ(By) +
SQ(Bz) + sq_m1, E1_sol);
176 for (
const double E1: E1_sol) {
181 const double p1x = Ax * E1 + Bx;
182 const double p1y = Ay * E1 + By;
183 const double p1z = Az * E1 + Bz;
185 LorentzVector p1(p1x, p1y, p1z, E1);
189 LOG(trace) <<
"[SecondaryBlockA] Throwing solution because of invalid mass. " <<
190 "Expected " << m1 <<
", got " << p1.M();
197 LOG(trace) <<
"[SecondaryBlockA] Throwing solution because of invalid invariant mass. " <<
198 "Expected " << *s1234 <<
", got " << (p1 + *m_p2 + *m_p3 + *m_p4).M2();
205 LOG(trace) <<
"[SecondaryBlockA] Throwing solution because of invalid invariant mass. " <<
206 "Expected " << *s123 <<
", got " << (p1 + *m_p2 + *m_p3).M2();
213 LOG(trace) <<
"[SecondaryBlockA] Throwing solution because of invalid invariant mass. " <<
214 "Expected " << *s12 <<
", got " << (p1 + *m_p2).M2();
220 const double jacobian = 1. / (128 * std::pow(M_PI, 3) * std::abs(
221 E4 * (p1z*p2y*p3x - p1y*p2z*p3x - p1z*p2x*p3y + p1x*p2z*p3y + p1y*p2x*p3z - p1x*p2y*p3z)
222 + E2 * (p1z*p3y*p4x - p1y*p3z*p4x - p1z*p3x*p4y + p1x*p3z*p4y + p1y*p3x*p4z - p1x*p3y*p4z)
223 + E1 * (- p2z*p3y*p4x + p2y*p3z*p4x + p2z*p3x*p4y - p2x*p3z*p4y - p2y*p3x*p4z + p2x*p3y*p4z)
224 + E3 * (- p1z*p2y*p4x + p1y*p2z*p4x + p1z*p2x*p4y - p1x*p2z*p4y - p1y*p2x*p4z + p1x*p2y*p4z)
227 Solution solution { { p1 }, jacobian,
true };
228 solutions->push_back(solution);
231 return (solutions->size() > 0) ? Status::OK : Status::NEXT;
248 std::shared_ptr<SolutionCollection> solutions = produce<SolutionCollection>(
"solutions");
260 .GlobalAttr(
"energy: double")
261 .Attr(
"m1: double=0");
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.
bool solveQuadratic(const double a, const double b, const double c, std::vector< double > &roots, bool verbose=false)
Finds the real solutions to .
virtual Status work() override
Main function.
Secondary Block A, describing
Parent class for all the modules.
A class encapsulating a lua table.
Module(PoolPtr pool, const std::string &name)
Constructor.