Skip to content
Snippets Groups Projects
Commit 48fb8744 authored by Tamas Gal's avatar Tamas Gal :speech_balloon:
Browse files

Cleanup the PDF wrapper

parent eaa5bc12
Branches 1-addition-to-readme
Tags v3.2.3
No related merge requests found
Pipeline #10145 passed
Stage: init
Stage: build
Stage: test
Hide whitespace changes
Inline Side-by-side
Showing
with 1 addition and 271 deletions
src/pdf.cc
+ 1
− 271
View file @ 48fb8744
#include <pybind11/pybind11.h> #include <pybind11/pybind11.h>
#include "JLang/JException.hh" #include "JPhysics/JPDF_t.hh"
#include "JTools/JCollection.hh"
#include "JTools/JMap.hh"
#include "JTools/JGridMap.hh"
#include "JTools/JMapList.hh"
#include "JTools/JSpline.hh"
#include "JTools/JPolint.hh"
#include "JTools/JElement.hh"
#include "JTools/JResult.hh"
#include "JPhysics/JPDFTable.hh"
#include "JPhysics/JPDFToolkit.hh"
#include "JPhysics/JPDFTypes.hh"
#include "JMath/JZero.hh"
/**
* Auxiliary data structure for muon PDF.
*/
struct JMuonPDF_t {
typedef JPP::JSplineFunction1D<JPP::JSplineElement2S<double, double>,
JPP::JCollection,
JPP::JResultPDF<double> > JFunction1D_t;
typedef JPP::JMAPLIST<JPP::JPolint1FunctionalMap,
JPP::JPolint0FunctionalGridMap,
JPP::JPolint0FunctionalGridMap>::maplist JPDFMaplist_t;
typedef JPP::JPDFTable<JFunction1D_t, JPDFMaplist_t> JPDF_t;
typedef JFunction1D_t::result_type result_type;
/**
* Constructor.
*
* The PDF file descriptor should contain the wild card character JPHYSICS::WILD_CARD.\n
* The <tt>TTS</tt> corresponds to the additional time smearing applied to the PDFs.
*
* \param fileDescriptor PDF file descriptor
* \param TTS TTS [ns]
* \param numberOfPoints number of points for Gauss-Hermite integration of TTS
* \param epsilon precision for Gauss-Hermite integration of TTS
*/
JMuonPDF_t(const std::string& fileDescriptor,
const double TTS,
const int numberOfPoints = 25,
const double epsilon = 1.0e-10)
{
using namespace std;
using namespace JPP;
const JPDFType_t pdf_t[] = { DIRECT_LIGHT_FROM_MUON,
SCATTERED_LIGHT_FROM_MUON,
DIRECT_LIGHT_FROM_EMSHOWERS,
SCATTERED_LIGHT_FROM_EMSHOWERS,
DIRECT_LIGHT_FROM_DELTARAYS,
SCATTERED_LIGHT_FROM_DELTARAYS };
const int N = sizeof(pdf_t) / sizeof(pdf_t[0]);
JPDF_t pdf[N];
const JPDF_t::JSupervisor supervisor(new JPDF_t::JDefaultResult(zero));
for (int i = 0; i != N; ++i) {
const string file_name = getFilename(fileDescriptor, pdf_t[i]);
cout << "loading input from file " << file_name << "... " << flush;
pdf[i].load(file_name.c_str());
pdf[i].setExceptionHandler(supervisor);
cout << "OK" << endl;
}
// Add PDFs
cout << "adding PDFs... " << flush;
pdfA = pdf[1]; pdfA.add(pdf[0]);
pdfB = pdf[3]; pdfB.add(pdf[2]);
pdfC = pdf[5]; pdfC.add(pdf[4]);
cout << "OK" << endl;
if (TTS > 0.0) {
cout << "bluring PDFs... " << flush;
pdfA.blur(TTS, numberOfPoints, epsilon);
pdfB.blur(TTS, numberOfPoints, epsilon);
pdfC.blur(TTS, numberOfPoints, epsilon);
cout << "OK" << endl;
} else if (TTS < 0.0) {
THROW(JValueOutOfRange, "Illegal value of TTS [ns]: " << TTS);
}
}
/**
* Get PDF.
*
* The orientation of the PMT should be defined according this <a href="https://common.pages.km3net.de/jpp/JPDF.PDF">documentation</a>.\n
* In this, the zenith and azimuth angles are limited to \f[\left[0, \pi\right]\f].
*
* \param E muon energy at minimum distance of approach [GeV]
* \param R minimum distance of approach [m]
* \param theta PMT zenith angle [rad]
* \param phi PMT azimuth angle [rad]
* \param t1 arrival time relative to Cherenkov hypothesis [ns]
* \return hypothesis value
*/
result_type calculate(const double E,
const double R,
const double theta,
const double phi,
const double t1) const
{
using namespace JPP;
result_type h1 = (pdfA(R, theta, phi, t1) +
pdfB(R, theta, phi, t1) * E +
pdfC(R, theta, phi, t1) * getDeltaRaysFromMuon(E));
// safety measures
if (h1.f <= 0.0) {
h1.f = 0.0;
h1.fp = 0.0;
}
if (h1.v <= 0.0) {
h1.v = 0.0;
}
return h1;
}
JPDF_t pdfA; //!< PDF for minimum ionisong particle
JPDF_t pdfB; //!< PDF for average energy losses
JPDF_t pdfC; //!< PDF for delta-rays
};
/**
* Auxiliary data structure for shower PDF.
*/
struct JShowerPDF_t {
typedef JPP::JSplineFunction1D<JPP::JSplineElement2S<double, double>,
JPP::JCollection,
JPP::JResultPDF<double> > JFunction1D_t;
typedef JPP::JMAPLIST<JPP::JPolint1FunctionalMap,
JPP::JPolint1FunctionalMap,
JPP::JPolint0FunctionalGridMap,
JPP::JPolint0FunctionalGridMap>::maplist JPDFMaplist_t;
typedef JPP::JPDFTable<JFunction1D_t, JPDFMaplist_t> JPDF_t;
typedef JFunction1D_t::result_type result_type;
/**
* Constructor.
*
* The PDF file descriptor should contain the wild card character JPHYSICS::WILD_CARD.\n
* The <tt>TTS</tt> corresponds to the additional time smearing applied to the PDFs.
*
* \param fileDescriptor PDF file descriptor
* \param TTS TTS [ns]
* \param numberOfPoints number of points for Gauss-Hermite integration of TTS
* \param epsilon precision for Gauss-Hermite integration of TTS
*/
JShowerPDF_t(const std::string& fileDescriptor,
const double TTS,
const int numberOfPoints = 25,
const double epsilon = 1.0e-10)
{
using namespace std;
using namespace JPP;
const JPDFType_t pdf_t[] = { SCATTERED_LIGHT_FROM_EMSHOWER,
DIRECT_LIGHT_FROM_EMSHOWER };
const int N = sizeof(pdf_t) / sizeof(pdf_t[0]);
const JPDF_t::JSupervisor supervisor(new JPDF_t::JDefaultResult(zero));
for (int i = 0; i != N; ++i) {
const string file_name = getFilename(fileDescriptor, pdf_t[i]);
cout << "loading input from file " << file_name << "... " << flush;
JPDF_t pdf;
pdf.load(file_name.c_str());
pdf.setExceptionHandler(supervisor);
if (pdfA.empty())
pdfA = pdf;
else
pdfA.add(pdf);
cout << "OK" << endl;
}
if (TTS > 0.0) {
cout << "bluring PDFs... " << flush;
pdfA.blur(TTS, numberOfPoints, epsilon);
cout << "OK" << endl;
} else if (TTS < 0.0) {
THROW(JValueOutOfRange, "Illegal value of TTS [ns]: " << TTS);
}
}
/**
* Get PDF.
*
* The orientation of the PMT should be defined according this <a href="https://common.pages.km3net.de/jpp/JPDF.PDF">documentation</a>.\n
* In this, the zenith and azimuth angles are limited to \f[\left[0, \pi\right]\f].
*
* \param E shower energy at minimum distance of approach [GeV]
* \param D distance [m]
* \param cd cosine emission angle
* \param theta PMT zenith angle [rad]
* \param phi PMT azimuth angle [rad]
* \param t1 arrival time relative to Cherenkov hypothesis [ns]
* \return hypothesis value
*/
result_type calculate(const double E,
const double D,
const double cd,
const double theta,
const double phi,
const double t1) const
{
using namespace JPP;
result_type h1 = pdfA(D, cd, theta, phi, t1) * E;
// safety measures
if (h1.f <= 0.0) {
h1.f = 0.0;
h1.fp = 0.0;
}
if (h1.v <= 0.0) {
h1.v = 0.0;
}
return h1;
}
JPDF_t pdfA; //!< PDF for shower
};
/**
Python bindings.
*/
namespace py = pybind11; namespace py = pybind11;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment