GaussianTransferFunctionOnPt.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 
20 #include <momemta/Logging.h>
21 #include <momemta/Module.h>
22 #include <momemta/ParameterSet.h>
23 #include <momemta/Types.h>
24 #include <momemta/Math.h>
25 
26 #include <Math/DistFunc.h>
27 
36  public:
37 
38  GaussianTransferFunctionOnPtBase(PoolPtr pool, const ParameterSet& parameters): Module(pool, parameters.getModuleName()) {
39  m_reco_input = get<LorentzVector>(parameters.get<InputTag>("reco_particle"));
40 
41  m_sigma = parameters.get<double>("sigma", 0.10);
42  m_sigma_range = parameters.get<double>("sigma_range", 5);
43  m_min_Pt = parameters.get<double>("min_Pt", 0.);
44  }
45 
46  protected:
47  double m_min_Pt;
48  double m_sigma;
49  double m_sigma_range;
50 
51  // Input
52  Value<LorentzVector> m_reco_input;
53 };
54 
96  public:
97  GaussianTransferFunctionOnPt(PoolPtr pool, const ParameterSet& parameters): GaussianTransferFunctionOnPtBase(pool, parameters) {
98  m_ps_point = get<double>(parameters.get<InputTag>("ps_point"));
99  }
100 
101  virtual Status work() override {
102  // Estimate the width over which to integrate using the width of the TF at Pt_rec ...
103  const double sigma_Pt_rec = m_reco_input->Pt() * m_sigma;
104 
105  const double cosh_eta = std::cosh(m_reco_input->Eta());
106  double range_min = std::max(m_min_Pt, m_reco_input->Pt() - (m_sigma_range * sigma_Pt_rec));
107  double range_max = m_reco_input->Pt() + (m_sigma_range * sigma_Pt_rec);
108  double range = (range_max - range_min);
109 
110  double gen_Pt = range_min + range * (*m_ps_point);
111 
112  // To change the particle's Pt without changing its direction and mass:
113  const double gen_E = std::sqrt(SQ(m_reco_input->M()) + SQ(cosh_eta * gen_Pt));
114 
115  output->SetCoordinates(
116  gen_Pt * std::cos(m_reco_input->Phi()),
117  gen_Pt * std::sin(m_reco_input->Phi()),
118  gen_Pt * std::sinh(m_reco_input->Eta()),
119  gen_E);
120 
121  // ... but compute the width of the TF at Pt_gen!
122  const double sigma_Pt_gen = gen_Pt * m_sigma;
123 
124  // Compute TF*jacobian, where the jacobian includes the transformation of [0,1]->[range_min,range_max] and d|P|/dPt = cosh(eta)
125  *TF_times_jacobian = ROOT::Math::normal_pdf(gen_Pt, sigma_Pt_gen, m_reco_input->Pt()) * range * cosh_eta;
126 
127  return Status::OK;
128  }
129 
130  private:
131  // Input
132  Value<double> m_ps_point;
133 
134  // Outputs
135  std::shared_ptr<LorentzVector> output = produce<LorentzVector>("output");
136  std::shared_ptr<double> TF_times_jacobian = produce<double>("TF_times_jacobian");
137 
138 };
139 
174  public:
175  GaussianTransferFunctionOnPtEvaluator(PoolPtr pool, const ParameterSet& parameters): GaussianTransferFunctionOnPtBase(pool, parameters) {
176  m_gen_input = get<LorentzVector>(parameters.get<InputTag>("gen_particle"));
177  }
178 
179  virtual Status work() override {
180  // Compute TF value
181  *TF_value = ROOT::Math::normal_pdf(m_gen_input->Pt(), m_gen_input->Pt() * m_sigma, m_reco_input->Pt());
182 
183  return Status::OK;
184  }
185 
186  private:
187  // Input
188  Value<LorentzVector> m_gen_input;
189 
190  // Outputs
191  std::shared_ptr<double> TF_value = produce<double>("TF");
192 
193 };
194 
195 REGISTER_MODULE(GaussianTransferFunctionOnPt)
196  .Input("ps_point")
197  .Input("reco_particle")
198  .Output("output")
199  .Output("TF_times_jacobian")
200  .Attr("sigma:double=0.10")
201  .Attr("sigma_range:double=5")
202  .Attr("min_Pt:double=0");
203 
205  .Input("gen_particle")
206  .Input("reco_particle")
207  .Output("TF")
208  .Attr("sigma:double=0.10")
209  .Attr("sigma_range:double=5")
210  .Attr("min_Pt:double=0");
Mathematical functions.
Helper class for Gaussian transfer function modules.
Integrate over a transfer function on Pt described by a Gaussian distribution.
virtual Status work() override
Main function.
An identifier of a module&#39;s output.
Definition: InputTag_fwd.h:37
virtual Status work() override
Main function.
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
Evaluate a transfer function on Pt described by a Gaussian distribution.