SecondaryBlockA.cc
1 /*
2  * MoMEMta: a modular implementation of the Matrix Element Method
3  * Copyright (C) 2017 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/ParameterSet.h>
20 #include <momemta/Module.h>
21 #include <momemta/Solution.h>
22 #include <momemta/Math.h>
23 #include <momemta/InputTag.h>
24 #include <momemta/Types.h>
25 
26 #include <Math/GenVector/VectorUtil.h>
27 
79 class SecondaryBlockA: public Module {
80  public:
81 
82  SecondaryBlockA(PoolPtr pool, const ParameterSet& parameters): Module(pool, parameters.getModuleName()),
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"));
87 
88  m1 = parameters.get<double>("m1", 0.);
89 
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"));
93  };
94 
95  virtual Status work() override {
96 
97  solutions->clear();
98 
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);
106 
107  // Don't spend time on unphysical part of phase-space
108  if (*s1234 >= SQ(sqrt_s) || *s12 + sq_m3 >= *s123 || *s123 + sq_m4 >= *s1234 || sq_m1 + sq_m2 >= *s12)
109  return Status::NEXT;
110 
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();
115 
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();
120 
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();
125 
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);
129 
130  /* Analytically solve the system:
131  * ai1 p1x + ai2 p1y + ai3 p1z = bi E1 + ci, with i = 1...3
132  *
133  * This gives p1 as a function of E1:
134  * p1x = Ax E1 + Bx
135  * p1y = Ay E1 + By
136  * p1z = Az E1 + Bz
137  */
138  const double a11 = p2x;
139  const double a12 = p2y;
140  const double a13 = p2z;
141 
142  const double a21 = p2x + p3x;
143  const double a22 = p2y + p3y;
144  const double a23 = p2z + p3z;
145 
146  const double a31 = p2x + p3x + p4x;
147  const double a32 = p2y + p3y + p4y;
148  const double a33 = p2z + p3z + p4z;
149 
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;
156 
157  const double inv_det = 1 / ((a13 * a22 - a12 * a23) * a31 - (a13 * a21 - a11 * a23) * a32 + (a12 * a21 - a11 * a22) * a33);
158 
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;
161 
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;
164 
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;
167 
168  // Now the mass-shell condition for p1 gives a quadratic equation in E1 with up to two solutions
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);
171 
172  if (!foundSolution)
173  return Status::NEXT;
174 
175  // Use now the obtained solutions of E1 to build the full p1
176  for (const double E1: E1_sol) {
177  // Skip unphysical solutions
178  if (E1 <= 0)
179  continue;
180 
181  const double p1x = Ax * E1 + Bx;
182  const double p1y = Ay * E1 + By;
183  const double p1z = Az * E1 + Bz;
184 
185  LorentzVector p1(p1x, p1y, p1z, E1);
186 
187  if (!ApproxComparison(p1.M() / p1.E(), m1 / p1.E())) {
188 #ifndef NDEBUG
189  LOG(trace) << "[SecondaryBlockA] Throwing solution because of invalid mass. " <<
190  "Expected " << m1 << ", got " << p1.M();
191 #endif
192  continue;
193  }
194 
195  if (!ApproxComparison((p1 + *m_p2 + *m_p3 + *m_p4).M2(), *s1234)) {
196 #ifndef NDEBUG
197  LOG(trace) << "[SecondaryBlockA] Throwing solution because of invalid invariant mass. " <<
198  "Expected " << *s1234 << ", got " << (p1 + *m_p2 + *m_p3 + *m_p4).M2();
199 #endif
200  continue;
201  }
202 
203  if (!ApproxComparison((p1 + *m_p2 + *m_p3).M2(), *s123)) {
204 #ifndef NDEBUG
205  LOG(trace) << "[SecondaryBlockA] Throwing solution because of invalid invariant mass. " <<
206  "Expected " << *s123 << ", got " << (p1 + *m_p2 + *m_p3).M2();
207 #endif
208  continue;
209  }
210 
211  if (!ApproxComparison((p1 + *m_p2).M2(), *s12)) {
212 #ifndef NDEBUG
213  LOG(trace) << "[SecondaryBlockA] Throwing solution because of invalid invariant mass. " <<
214  "Expected " << *s12 << ", got " << (p1 + *m_p2).M2();
215 #endif
216  continue;
217  }
218 
219  // Compute jacobian
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)
225  ));
226 
227  Solution solution { { p1 }, jacobian, true };
228  solutions->push_back(solution);
229  }
230 
231  return (solutions->size() > 0) ? Status::OK : Status::NEXT;
232 
233  }
234 
235  private:
236  double sqrt_s;
237  double m1;
238 
239  // Inputs
240  Value<double> s12;
241  Value<double> s123;
242  Value<double> s1234;
246 
247  // Output
248  std::shared_ptr<SolutionCollection> solutions = produce<SolutionCollection>("solutions");
249 };
250 
251 REGISTER_MODULE(SecondaryBlockA)
252  .Input("s12")
253  .Input("s123")
254  .Input("s1234")
255  .Input("p1")
256  .Input("p2")
257  .Input("p3")
258  .Input("p4")
259  .Output("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.
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
virtual Status work() override
Main function.
Secondary Block A, describing
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