SecondaryBlockE.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 
77 class SecondaryBlockE: public Module {
78  public:
79 
80  SecondaryBlockE(PoolPtr pool, const ParameterSet& parameters): Module(pool, parameters.getModuleName()),
81  sqrt_s(parameters.globalParameters().get<double>("energy")) {
82  s12 = get<double>(parameters.get<InputTag>("s12"));
83  s123 = get<double>(parameters.get<InputTag>("s123"));
84 
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"));
88  };
89 
90  virtual Status work() override {
91 
92  solutions->clear();
93 
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);
98 
99  // Don't spend time on unphysical part of phase-space
100  if (*s123 >= SQ(sqrt_s) || *s12 + sq_m3 >= *s123 || sq_m1 >= *s12)
101  return Status::NEXT;
102 
103  const double p3 = m_p3->P();
104  const double E3 = m_p3->E();
105  const double sq_E3 = SQ(E3);
106 
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);
110 
111  double X = p3 * c23 - E3;
112  double Y = *s123 - *s12 - sq_m3;
113 
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,
117  abs_p2, abs_p1);
118 
119  // Use now the obtained |p1| and |p2| solutions to build p1 and p2 (m2=0!)
120  for (std::size_t i = 0; i < abs_p1.size(); i++) {
121  // Skip unphysical solutions
122  if (abs_p1[i] <= 0 || abs_p2[i] <= 0)
123  continue;
124 
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);
128 
129  LorentzVector p1_sol, p2_sol;
130  p1_sol.SetPxPyPzE(
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()),
134  E1);
135  p2_sol.SetPxPyPzE(
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()),
139  abs_p2[i]);
140 
141  if (!ApproxComparison(p1_sol.M() / p1_sol.E(), m1 / p1_sol.E())) {
142 #ifndef NDEBUG
143  LOG(trace) << "[SecondaryBlockE] Throwing solution because of invalid mass. " <<
144  "Expected " << m1 << ", got " << p1_sol.M();
145 #endif
146  continue;
147  }
148 
149  if (!ApproxComparison(p2_sol.M() / p2_sol.E(), 0.)) {
150 #ifndef NDEBUG
151  LOG(trace) << "[SecondaryBlockE] Throwing solution because of invalid mass. " <<
152  "Expected 0., got " << p2_sol.M();
153 #endif
154  continue;
155  }
156 
157  if (!ApproxComparison((p1_sol + p2_sol + *m_p3).M2(), *s123)) {
158 #ifndef NDEBUE
159  LOG(trace) << "[SecondaryBlockE] Throwing solution because of invalid invariant mass. " <<
160  "Expected " << *s123 << ", got " << (p1_sol + p2_sol + *m_p3).M2();
161 #endif
162  continue;
163  }
164 
165  if (!ApproxComparison((p1_sol + p2_sol).M2(), *s12)) {
166 #ifndef NDEBUG
167  LOG(trace) << "[SecondaryBlockE] Throwing solution because of invalid invariant mass. " <<
168  "Expected " << *s12 << ", got " << (p1_sol + p2_sol).M2();
169 #endif
170  continue;
171  }
172 
173  // Compute jacobian
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) )
177  );
178 
179  Solution solution { { p1_sol, p2_sol }, jacobian, true };
180  solutions->push_back(solution);
181  }
182 
183  return (solutions->size() > 0) ? Status::OK : Status::NEXT;
184 
185  }
186 
187  private:
188  double sqrt_s;
189 
190  // Inputs
191  Value<double> s12;
192  Value<double> s123;
196 
197  // Output
198  std::shared_ptr<SolutionCollection> solutions = produce<SolutionCollection>("solutions");
199 };
200 
201 REGISTER_MODULE(SecondaryBlockE)
202  .Input("s12")
203  .Input("s123")
204  .Input("p1")
205  .Input("p2")
206  .Input("p3")
207  .Output("solutions")
208  .GlobalAttr("energy: double");
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.
virtual Status work() override
Main function.
An identifier of a module&#39;s output.
Definition: InputTag_fwd.h:37
Parent class for all the modules.
Definition: Module.h:37
Secondary Block E, describing
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
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.
Definition: Math.cc:192