19 #include <momemta/Logging.h> 20 #include <momemta/Module.h> 21 #include <momemta/ParameterSet.h> 22 #include <momemta/Solution.h> 23 #include <momemta/Types.h> 100 m_ps_point1 = get<double>(parameters.get<
InputTag>(
"q1"));
101 m_ps_point2 = get<double>(parameters.get<
InputTag>(
"q2"));
103 m1 = parameters.get<
double>(
"m1", 0.);
104 m2 = parameters.get<
double>(
"m2", 0.);
106 sqrt_s = parameters.globalParameters().get<
double>(
"energy");
108 s13 = get<double>(parameters.get<
InputTag>(
"s13"));
109 s24 = get<double>(parameters.get<
InputTag>(
"s24"));
111 p3 = get<LorentzVector>(parameters.get<
InputTag>(
"p3"));
112 p4 = get<LorentzVector>(parameters.get<
InputTag>(
"p4"));
114 if (parameters.exists(
"branches")) {
115 auto branches_tags = parameters.get<std::vector<InputTag>>(
"branches");
116 for (
auto& t: branches_tags)
117 m_branches.push_back(get<LorentzVector>(t));
121 virtual Status
work()
override {
125 const double sq_m1 =
SQ(m1);
126 const double sq_m2 =
SQ(m2);
127 const double sq_m3 = p3->M2();
128 const double sq_m4 = p4->M2();
131 if (sq_m1 + sq_m3 >= *s13 || sq_m2 + sq_m4 >= *s24 || *s13 + *s24 >
SQ(sqrt_s))
134 const double p3x = p3->Px();
135 const double p3y = p3->Py();
136 const double p3z = p3->Pz();
137 const double E3 = p3->E();
139 const double p4x = p4->Px();
140 const double p4y = p4->Py();
141 const double p4z = p4->pz();
142 const double E4 = p4->E();
145 LorentzVector pb = *p3 + *p4;
146 for (
size_t i = 0; i < m_branches.size(); i++) {
147 pb += *m_branches[i];
150 const double Eb = pb.E();
151 const double pbx = pb.Px();
152 const double pby = pb.Py();
153 const double pbz = pb.Pz();
155 const double q1 = *m_ps_point1;
156 const double q2 = *m_ps_point2;
158 const double Etot = sqrt_s * (q1 + q2) / 2 - Eb;
159 const double ptotz = sqrt_s * (q1 - q2) / 2 - pbz;
161 const double X = 0.5 * (- sq_m1 - sq_m3 + *s13);
162 const double Y = 0.5 * (- sq_m2 - sq_m4 + *s24);
182 const double den = p3z * p4x - p3x * p4z;
184 const double A1x = - (E4 * p3z - E3 * p4z) / den;
185 const double B1x = (p3z * p4y - p3y * p4z) / den;
186 const double C1x = - (p4z * (E3 * Etot - p3z * ptotz + p3y * pby - X) - p3z * (Y - p4x * pbx)) / den;
188 const double A1z = (E4 * p3x - E3 * p4x) / den;
189 const double B1z = (p3y * p4x - p3x * p4y) / den;
190 const double C1z = (p4x * (E3 * Etot + p3y * pby + p3x * pbx - X) - p3x * (Y + p4z * ptotz)) / den;
204 const double fac = 2 * (A1x * pbx - A1z * ptotz - Etot);
205 const double a = - 2 * (B1x * pbx - B1z * ptotz - pby) / fac;
206 const double b = - (
SQ(Etot) + pow(C1x + pbx, 2) + pow(C1z - ptotz, 2) + sq_m2 -
SQ(C1x) -
SQ(pby) -
SQ(C1z) - sq_m1) / fac;
208 const double a20 = 1 -
SQ(A1x) -
SQ(A1z);
209 const double a02 = - (
SQ(B1x) +
SQ(B1z) + 1);
210 const double a11 = - 2 * (A1x * B1x + A1z * B1z);
211 const double a10 = - 2 * (A1x * C1x + A1z * C1z + Etot);
212 const double a01 = - 2 * (B1x * C1x + B1z * C1z + pby);
213 const double a00 =
SQ(Etot) - (
SQ(C1x) +
SQ(C1z) +
SQ(pby) + sq_m1);
215 std::vector<double> p2y_sol;
217 2 * a * b * a20 + b * a11 + a01 + a * a10,
218 SQ(b) * a20 + b * a10 + a00,
224 for (
const double p2y: p2y_sol) {
225 const double E2 = a * p2y + b;
230 const double E1 = Etot - E2;
235 const double p1x = A1x * E2 + B1x * p2y + C1x;
236 const double p1y = - p2y - pby;
237 const double p1z = A1z * E2 + B1z * p2y + C1z;
238 const LorentzVector p1(p1x, p1y, p1z, E1);
240 const double p2x = - p1x - pbx;
241 const double p2z = - p1z + ptotz;
242 const LorentzVector p2(p2x, p2y, p2z, E2);
245 const LorentzVector tot = p1 + p2 + pb;
246 const double q1Pz = std::abs(tot.Pz() + tot.E()) / 2.;
247 const double q2Pz = std::abs(tot.Pz() - tot.E()) / 2.;
249 if (q1Pz > sqrt_s / 2 || q2Pz > sqrt_s / 2)
254 LOG(trace) <<
"[BlockF] Throwing solution because total Pt is incorrect. " 255 <<
"Expected " << 0. <<
", got " << tot.Pt();
262 LOG(trace) <<
"[BlockF] Throwing solution because p1 has an invalid mass. " <<
263 "Expected " << m1 <<
", got " << p1.M();
270 LOG(trace) <<
"[BlockF] Throwing solution because p2 has an invalid mass. " <<
271 "Expected " << m2 <<
", got " << p2.M();
278 LOG(trace) <<
"[BlockF] Throwing solution because of invalid invariant mass. " <<
279 "Expected " << *s13 <<
", got " << (p1 + *p3).M2();
286 LOG(trace) <<
"[BlockF] Throwing solution because of invalid invariant mass. " <<
287 "Expected " << *s24 <<
", got " << (p2 + *p4).M2();
292 const double jacobian = 1. / (64 *
SQ(M_PI) * std::abs(E4*(p1z*p2y*p3x - p1y*p2z*p3x - p1z*p2x*p3y + p1x*p2z*p3y + p1y*p2x*p3z - p1x*p2y*p3z) + E2*p1z*p3y*p4x - E1*p2z*p3y*p4x - E2*p1y*p3z*p4x + E1*p2y*p3z*p4x - E2*p1z*p3x*p4y + E1*p2z*p3x*p4y + E2*p1x*p3z*p4y - E1*p2x*p3z*p4y + (E2*p1y*p3x - E1*p2y*p3x - E2*p1x*p3y + E1*p2x*p3y)*p4z + E3*(-(p1z*p2y*p4x) + p1y*p2z*p4x + p1z*p2x*p4y - p1x*p2z*p4y - p1y*p2x*p4z + p1x*p2y*p4z)));
294 Solution s { {p1, p2}, jacobian,
true };
295 solutions->push_back(s);
298 return solutions->size() > 0 ? Status::OK : Status::NEXT;
312 std::vector<Value<LorentzVector>> m_branches;
315 std::shared_ptr<SolutionCollection> solutions = produce<SolutionCollection>(
"solutions");
325 .OptionalInputs(
"branches")
327 .GlobalAttr(
"energy:double")
329 .Attr(
"m2: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 .
Final (main) Block F, describing
Parent class for all the modules.
A class encapsulating a lua table.
Module(PoolPtr pool, const std::string &name)
Constructor.
virtual Status work() override
Main function.