19 #include <momemta/Logging.h> 20 #include <momemta/Module.h> 21 #include <momemta/ParameterSet.h> 22 #include <momemta/Types.h> 42 m_reco_input = get<LorentzVector>(parameters.get<
InputTag>(
"reco_particle"));
44 std::string file_path = parameters.get<std::string>(
"file");
45 std::string th2_name = parameters.get<std::string>(
"th2_name");
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);
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);
56 TAxis* yAxis = m_th2->GetYaxis();
57 m_deltaMin = yAxis->GetXmin();
58 m_deltaMax = yAxis->GetXmax();
59 m_deltaRange = m_deltaMax - m_deltaMin;
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();
69 m_fallBackEgenMax = xAxis->GetBinCenter(xAxis->GetNbins());
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.";
81 std::unique_ptr<TH2> m_th2;
83 double m_deltaMin, m_deltaMax, m_deltaRange;
84 double m_EgenMin, m_EgenMax;
85 double m_fallBackEgenMax;
91 class file_not_found_error:
public std::runtime_error{
92 using std::runtime_error::runtime_error;
95 class th2_not_found_error:
public std::runtime_error{
96 using std::runtime_error::runtime_error;
147 m_ps_point = get<double>(parameters.get<
InputTag>(
"ps_point"));
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;
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()),
167 const int bin = m_th2->FindFixBin(std::min(gen_E, m_fallBackEgenMax), delta);
170 *TF_times_jacobian = m_th2->GetBinContent(bin) * range *
dP_over_dE(*output);
181 std::shared_ptr<LorentzVector> output = produce<LorentzVector>(
"output");
182 std::shared_ptr<double> TF_times_jacobian = produce<double>(
"TF_times_jacobian");
188 inline double GetDeltaRange(
const double Erec,
const double Mrec)
const {
189 return GetDeltaMax(Erec, Mrec) - m_deltaMin;
192 inline double GetDeltaMax(
const double Erec,
const double Mrec)
const {
193 if (Erec <= Mrec || m_deltaMax <= Mrec)
195 return std::min(m_deltaMax, Erec - m_EgenMin) - Mrec;
240 m_gen_input = get<LorentzVector>(parameters.get<
InputTag>(
"gen_particle"));
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) {
253 const int bin = m_th2->FindFixBin(std::min(gen_E, m_fallBackEgenMax), delta);
256 *TF_value = m_th2->GetBinContent(bin);
267 std::shared_ptr<double> TF_value = produce<double>(
"TF");
271 .Input(
"reco_particle")
274 .Output(
"TF_times_jacobian")
276 .Attr(
"th2_name:string")
277 .Attr(
"min_E:double=0");
280 .Input(
"reco_particle")
281 .Input(
"gen_particle")
284 .Attr(
"th2_name:string")
285 .Attr(
"min_E:double=0");
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.
Parent class for all the modules.
A class encapsulating a lua table.
double dP_over_dE(const T &v)
Compute .
Evaluate a transfer function on energy described by a 2D histogram retrieved from a ROOT file...
Module(PoolPtr pool, const std::string &name)
Constructor.
virtual Status work() override
Main function.