BinnedTransferFunctionOnEnergy.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/Types.h>
23 #include <momemta/Math.h>
24 
25 
26 #include <TFile.h>
27 #include <TH2.h>
28 #include <TAxis.h>
29 
39  public:
40 
41  BinnedTransferFunctionOnEnergyBase(PoolPtr pool, const ParameterSet& parameters): Module(pool, parameters.getModuleName()) {
42  m_reco_input = get<LorentzVector>(parameters.get<InputTag>("reco_particle"));
43 
44  std::string file_path = parameters.get<std::string>("file");
45  std::string th2_name = parameters.get<std::string>("th2_name");
46 
47  std::unique_ptr<TFile> file(TFile::Open(file_path.c_str()));
48  if(!file->IsOpen() || file->IsZombie())
49  throw file_not_found_error("Could not open file " + file_path);
50 
51  m_th2 = std::unique_ptr<TH2>(static_cast<TH2*>(file->Get(th2_name.c_str())));
52  if(!m_th2->InheritsFrom("TH2") || !m_th2.get())
53  throw th2_not_found_error("Could not retrieve object " + th2_name + " deriving from class TH2 in file " + file_path + ".");
54  m_th2->SetDirectory(0);
55 
56  TAxis* yAxis = m_th2->GetYaxis();
57  m_deltaMin = yAxis->GetXmin();
58  m_deltaMax = yAxis->GetXmax();
59  m_deltaRange = m_deltaMax - m_deltaMin;
60 
61  TAxis* xAxis = m_th2->GetXaxis();
62  double E_cut = parameters.get<double>("min_E", 0.);
63  m_EgenMin = std::max(xAxis->GetXmin(), E_cut);
64  m_EgenMax = xAxis->GetXmax();
65 
66  // Since we assume the TF continues as a constant for E->infty,
67  // we need to be able to retrieve the X axis' last bin, to avoid
68  // fetching the TH2's overflow bin.
69  m_fallBackEgenMax = xAxis->GetBinCenter(xAxis->GetNbins());
70 
71  LOG(debug) << "Loaded TH2 " << th2_name << " from file " << file_path << ".";
72  LOG(debug) << "\tDelta range is " << m_deltaMin << " to " << m_deltaMax << ".";
73  LOG(debug) << "\tEnergy range is " << m_EgenMin << " to " << m_EgenMax << ".";
74  LOG(debug) << "\tWill use values at Egen = " << m_fallBackEgenMax << " for out-of-range values.";
75 
76  file->Close();
77  file.reset();
78  };
79 
80  protected:
81  std::unique_ptr<TH2> m_th2;
82 
83  double m_deltaMin, m_deltaMax, m_deltaRange;
84  double m_EgenMin, m_EgenMax;
85  double m_fallBackEgenMax;
86 
87  // Input
88  Value<LorentzVector> m_reco_input;
89 
90  private:
91  class file_not_found_error: public std::runtime_error{
92  using std::runtime_error::runtime_error;
93  };
94 
95  class th2_not_found_error: public std::runtime_error{
96  using std::runtime_error::runtime_error;
97  };
98 };
99 
144  public:
145 
146  BinnedTransferFunctionOnEnergy(PoolPtr pool, const ParameterSet& parameters): BinnedTransferFunctionOnEnergyBase(pool, parameters) {
147  m_ps_point = get<double>(parameters.get<InputTag>("ps_point"));
148  }
149 
150  virtual Status work() override {
151  const double rec_E = m_reco_input->E();
152  const double rec_M = m_reco_input->M();
153  const double range = GetDeltaRange(rec_E, rec_M);
154  const double gen_E = rec_E - GetDeltaMax(rec_E, rec_M) + range * (*m_ps_point);
155  const double delta = rec_E - gen_E;
156 
157  // To change the particle's energy without changing its direction and mass
158  // Forcing positive value of (gen_E - rec_M) due to numeric precision issue
159  double gen_pt = std::sqrt(std::max(gen_E - rec_M, 0.) * (gen_E + rec_M)) / std::cosh(m_reco_input->Eta());
160  output->SetCoordinates(
161  gen_pt * std::cos(m_reco_input->Phi()),
162  gen_pt * std::sin(m_reco_input->Phi()),
163  gen_pt * std::sinh(m_reco_input->Eta()),
164  gen_E);
165 
166  // The bin number is a ROOT "global bin number" using a 1D representation of the TH2
167  const int bin = m_th2->FindFixBin(std::min(gen_E, m_fallBackEgenMax), delta);
168 
169  // Compute TF*jacobian, where the jacobian includes the transformation of [0,1]->[range_min,range_max] and d|P|/dE
170  *TF_times_jacobian = m_th2->GetBinContent(bin) * range * dP_over_dE(*output);
171 
172  return Status::OK;
173  }
174 
175  private:
176 
177  // Input
178  Value<double> m_ps_point;
179 
180  // Outputs
181  std::shared_ptr<LorentzVector> output = produce<LorentzVector>("output");
182  std::shared_ptr<double> TF_times_jacobian = produce<double>("TF_times_jacobian");
183 
184  // We assume TF=0 for Egen < lower limit of the TH2, and
185  // TF=constant for Egen > upper limit of the TH2
186  // In the former case, the integrated range is adapted (shortened) to the range where TF!=0
187  // The range also takes into account the particle's given mass: Egen has to be larger than M
188  inline double GetDeltaRange(const double Erec, const double Mrec) const {
189  return GetDeltaMax(Erec, Mrec) - m_deltaMin;
190  }
191 
192  inline double GetDeltaMax(const double Erec, const double Mrec) const {
193  if (Erec <= Mrec || m_deltaMax <= Mrec)
194  return 0;
195  return std::min(m_deltaMax, Erec - m_EgenMin) - Mrec;
196  }
197 };
198 
237  public:
238 
239  BinnedTransferFunctionOnEnergyEvaluator(PoolPtr pool, const ParameterSet& parameters): BinnedTransferFunctionOnEnergyBase(pool, parameters) {
240  m_gen_input = get<LorentzVector>(parameters.get<InputTag>("gen_particle"));
241  }
242 
243  virtual Status work() override {
244  const double rec_E = m_reco_input->E();
245  const double gen_E = m_gen_input->E();
246  double delta = rec_E - gen_E;
247  if (delta <= m_deltaMin || delta >= m_deltaMax) {
248  *TF_value = 0;
249  return Status::OK;
250  }
251 
252  // The bin number is a ROOT "global bin number" using a 1D representation of the TH2
253  const int bin = m_th2->FindFixBin(std::min(gen_E, m_fallBackEgenMax), delta);
254 
255  // Retrieve TF value
256  *TF_value = m_th2->GetBinContent(bin);
257 
258  return Status::OK;
259  }
260 
261  private:
262 
263  // Input
264  Value<LorentzVector> m_gen_input;
265 
266  // Outputs
267  std::shared_ptr<double> TF_value = produce<double>("TF");
268 };
269 
270 REGISTER_MODULE(BinnedTransferFunctionOnEnergy)
271  .Input("reco_particle")
272  .Input("ps_point")
273  .Output("output")
274  .Output("TF_times_jacobian")
275  .Attr("file:string")
276  .Attr("th2_name:string")
277  .Attr("min_E:double=0");
278 
280  .Input("reco_particle")
281  .Input("gen_particle")
282  .Output("TF")
283  .Attr("file:string")
284  .Attr("th2_name:string")
285  .Attr("min_E:double=0");
Mathematical functions.
Helper class for binned transfer function modules.
Integrate over a transfer function on energy described by a 2D histogram retrieved from a ROOT file...
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
A class encapsulating a lua table.
Definition: ParameterSet.h:82
double dP_over_dE(const T &v)
Compute .
Definition: Math.h:53
Evaluate a transfer function on energy described by a 2D histogram retrieved from a ROOT file...
Module(PoolPtr pool, const std::string &name)
Constructor.
Definition: Module.h:61
virtual Status work() override
Main function.