BlockF.cc
1 /*
2  * MoMEMta: a modular implementation of the Matrix Element Method
3  * Copyright (C) 2016 Universite catholique de Louvain (UCL), Belgium
4  *
5  * This program is free software: you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation, either version 3 of the License, or
8  * (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program. If not, see <http://www.gnu.org/licenses/>.
17  */
18 
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>
24 #include <momemta/Math.h>
25 
95 class BlockF: public Module {
96  public:
97 
98  BlockF(PoolPtr pool, const ParameterSet& parameters): Module(pool, parameters.getModuleName()) {
99 
100  m_ps_point1 = get<double>(parameters.get<InputTag>("q1"));
101  m_ps_point2 = get<double>(parameters.get<InputTag>("q2"));
102 
103  m1 = parameters.get<double>("m1", 0.);
104  m2 = parameters.get<double>("m2", 0.);
105 
106  sqrt_s = parameters.globalParameters().get<double>("energy");
107 
108  s13 = get<double>(parameters.get<InputTag>("s13"));
109  s24 = get<double>(parameters.get<InputTag>("s24"));
110 
111  p3 = get<LorentzVector>(parameters.get<InputTag>("p3"));
112  p4 = get<LorentzVector>(parameters.get<InputTag>("p4"));
113 
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));
118  }
119  };
120 
121  virtual Status work() override {
122 
123  solutions->clear();
124 
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();
129 
130  // Don't spend time on unphysical part of phase-space
131  if (sq_m1 + sq_m3 >= *s13 || sq_m2 + sq_m4 >= *s24 || *s13 + *s24 > SQ(sqrt_s))
132  return Status::NEXT;
133 
134  const double p3x = p3->Px();
135  const double p3y = p3->Py();
136  const double p3z = p3->Pz();
137  const double E3 = p3->E();
138 
139  const double p4x = p4->Px();
140  const double p4y = p4->Py();
141  const double p4z = p4->pz();
142  const double E4 = p4->E();
143 
144  // Total visible momentum
145  LorentzVector pb = *p3 + *p4;
146  for (size_t i = 0; i < m_branches.size(); i++) {
147  pb += *m_branches[i];
148  }
149 
150  const double Eb = pb.E();
151  const double pbx = pb.Px();
152  const double pby = pb.Py();
153  const double pbz = pb.Pz();
154 
155  const double q1 = *m_ps_point1;
156  const double q2 = *m_ps_point2;
157 
158  const double Etot = sqrt_s * (q1 + q2) / 2 - Eb;
159  const double ptotz = sqrt_s * (q1 - q2) / 2 - pbz;
160 
161  const double X = 0.5 * (- sq_m1 - sq_m3 + *s13);
162  const double Y = 0.5 * (- sq_m2 - sq_m4 + *s24);
163 
164 
165  /* Solve the linear system:
166  * p1x + p2x = - pbx
167  * p1y + p2y = - pby
168  * p1z + p2z = ptotz
169  * p3x p1x + p3y p1y + p3z p1z = - X + E3 E1
170  * p4x p2x + p4z p2z = - Y + E4 E2 -p4y p2y
171  *
172  * The solutions are parameterised by E2 and p2y:
173  * p1x = A1x E2 + B1x p2y + C1x
174  * p1y = - p2y - pby
175  * p1z = A1z E2 + B1z p2y + C1z
176  * p2x = - A1x E2 - B1x p2y - (C1x + pbx)
177  * p2z = - A1z E2 - B1z p2y - (C1z - ptotz)
178  *
179  * where one has used that E1 = Etot - E2
180  */
181 
182  const double den = p3z * p4x - p3x * p4z;
183 
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;
187 
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;
191 
192 
193  /* Now, insert those expressions into the mass-shell conditions for p1 and p2:
194  * (Etot - E2)^2 - p1x^2 - p1y^2 - p1z^2 - m1^2 = 0 (1)
195  * E2^2 - p2x^2 - p2y^2 - p2z^2 - m2^2 = 0 (2)
196  *
197  * Using the above, (1) is written as:
198  * a20 E2^2 + a02 p2y^2 + a11 E2 p2y + a10 E2 + a01 p2y + a00 = 0
199  *
200  * From (2)-(1), is it possible to write E2 = a p2y + b. This is then inserted into (1),
201  * yielding a quadratic equation in p2y only.
202  */
203 
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;
207 
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);
214 
215  std::vector<double> p2y_sol;
216  const bool foundSolution = solveQuadratic(a02 + SQ(a) * a20 + a * a11,
217  2 * a * b * a20 + b * a11 + a01 + a * a10,
218  SQ(b) * a20 + b * a10 + a00,
219  p2y_sol);
220 
221  if (!foundSolution)
222  return Status::NEXT;
223 
224  for (const double p2y: p2y_sol) {
225  const double E2 = a * p2y + b;
226 
227  if (E2 <= 0)
228  continue;
229 
230  const double E1 = Etot - E2;
231 
232  if (E1 <= 0)
233  continue;
234 
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);
239 
240  const double p2x = - p1x - pbx;
241  const double p2z = - p1z + ptotz;
242  const LorentzVector p2(p2x, p2y, p2z, E2);
243 
244  // Check if solutions are physical
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.;
248 
249  if (q1Pz > sqrt_s / 2 || q2Pz > sqrt_s / 2)
250  continue;
251 
252  if (!ApproxComparison(tot.Pt(), 0.)) {
253 #ifndef NDEBUG
254  LOG(trace) << "[BlockF] Throwing solution because total Pt is incorrect. "
255  << "Expected " << 0. << ", got " << tot.Pt();
256 #endif
257  continue;
258  }
259 
260  if (!ApproxComparison(p1.M() / p1.E(), m1 / p1.E())) {
261 #ifndef NDEBUG
262  LOG(trace) << "[BlockF] Throwing solution because p1 has an invalid mass. " <<
263  "Expected " << m1 << ", got " << p1.M();
264 #endif
265  continue;
266  }
267 
268  if (!ApproxComparison(p2.M() / p2.E(), m2 / p2.E())) {
269 #ifndef NDEBUG
270  LOG(trace) << "[BlockF] Throwing solution because p2 has an invalid mass. " <<
271  "Expected " << m2 << ", got " << p2.M();
272 #endif
273  continue;
274  }
275 
276  if (!ApproxComparison((p1 + *p3).M2(), *s13)) {
277 #ifndef NDEBUG
278  LOG(trace) << "[BlockF] Throwing solution because of invalid invariant mass. " <<
279  "Expected " << *s13 << ", got " << (p1 + *p3).M2();
280 #endif
281  continue;
282  }
283 
284  if (!ApproxComparison((p2 + *p4).M2(), *s24)) {
285 #ifndef NDEBUG
286  LOG(trace) << "[BlockF] Throwing solution because of invalid invariant mass. " <<
287  "Expected " << *s24 << ", got " << (p2 + *p4).M2();
288 #endif
289  continue;
290  }
291 
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)));
293 
294  Solution s { {p1, p2}, jacobian, true };
295  solutions->push_back(s);
296  }
297 
298  return solutions->size() > 0 ? Status::OK : Status::NEXT;
299  }
300 
301  private:
302  double sqrt_s;
303 
304  // Inputs
305  Value<double> s13;
306  Value<double> s24;
307  Value<double> m_ps_point1;
308  Value<double> m_ps_point2;
309  double m1;
310  double m2;
311  Value<LorentzVector> p3, p4;
312  std::vector<Value<LorentzVector>> m_branches;
313 
314  // Outputs
315  std::shared_ptr<SolutionCollection> solutions = produce<SolutionCollection>("solutions");
316 };
317 
318 REGISTER_MODULE(BlockF)
319  .Input("q1")
320  .Input("q2")
321  .Input("s13")
322  .Input("s24")
323  .Input("p3")
324  .Input("p4")
325  .OptionalInputs("branches")
326  .Output("solutions")
327  .GlobalAttr("energy:double")
328  .Attr("m1:double=0")
329  .Attr("m2:double=0");
bool ApproxComparison(double value, double expected)
Compare two doubles and return true if they are approximatively equal.
Definition: Math.cc:409
Generic solution structure representing a set of particles, along with its jacobian.
Definition: Solution.h:28
Mathematical functions.
bool solveQuadratic(const double a, const double b, const double c, std::vector< double > &roots, bool verbose=false)
Finds the real solutions to .
Definition: Math.cc:26
Final (main) Block F, describing
Definition: BlockF.cc:95
An identifier of a module&#39;s output.
Definition: InputTag_fwd.h:37
Parent class for all the modules.
Definition: Module.h:37
A class encapsulating a lua table.
Definition: ParameterSet.h:82
#define SQ(x)
Compute .
Definition: Math.h:25
Module(PoolPtr pool, const std::string &name)
Constructor.
Definition: Module.h:61
virtual Status work() override
Main function.
Definition: BlockF.cc:121