34int main(
int argc,
char **argv)
53 JParser<> zap(
"Auxiliary program to determine PDF of L1 hit.");
57 zap[
'E'] =
make_field(E,
"muon energy [GeV]") = 1.0;
58 zap[
'R'] =
make_field(R,
"distance of approach [m]");
59 zap[
'D'] =
make_field(dir,
"(theta, phi) of PMT [rad]");
60 zap[
'T'] =
make_field(T_ns,
"time window [ns]") = 10.0;
61 zap[
'B'] =
make_field(R_Hz,
"background rate [Hz]") = 0.0;
67 catch(
const exception &error) {
68 FATAL(error.what() << endl);
77 JFunction1D_t::JSupervisor supervisor(
new JFunction1D_t::JDefaultResult(zero));
79 const JPDFType_t pdf_t[] = { DIRECT_LIGHT_FROM_MUON,
80 SCATTERED_LIGHT_FROM_MUON,
81 DIRECT_LIGHT_FROM_DELTARAYS,
82 SCATTERED_LIGHT_FROM_DELTARAYS,
83 DIRECT_LIGHT_FROM_EMSHOWERS,
84 SCATTERED_LIGHT_FROM_EMSHOWERS };
86 const int N =
sizeof(pdf_t) /
sizeof(pdf_t[0]);
90 for (
int i = 0; i != N; ++i) {
94 const string file_name = getFilename(inputFile, pdf_t[i]);
96 NOTICE(
"loading input from file " << file_name <<
"... " << flush);
98 pdf[i].load(file_name.c_str());
106 pdf[i].setExceptionHandler(supervisor);
117 JModule module = getModule<JKM3NeT_t>(1);
119 module.rotate(JRotation3D(dir));
130 const double t2 = t1 + T_ns;
132 JFunction1D_t::result_type y1;
133 JFunction1D_t::result_type y2;
135 for (
int i = 0; i != N; ++i) {
137 JFunction1D_t::result_type
p1;
138 JFunction1D_t::result_type p2;
142 for (
const auto& pmt : module) {
143 p1 += pdf[i](R, pmt.getTheta(), pmt.getPhi(), t1);
144 p2 += pdf[i](R, pmt.getTheta(), pmt.getPhi(), t2);
147 if (is_deltarays(pdf_t[i])) {
148 p1 *= getDeltaRaysFromMuon(E);
149 p2 *= getDeltaRaysFromMuon(E);
150 }
else if (is_bremsstrahlung(pdf_t[i])) {
165 zsp[t1] = y1.f * (1.0 - exp(y1.v - y2.v));
170 for (
int ix = 1; ix <= h0.GetNbinsX(); ++ix) {
172 const double t1 = h0.GetBinCenter(ix);
176 const JFunction1D_t::result_type result = zsp(t1);
178 const double W = exp(-result.v) * result.f / (1.0 - exp(-result.V));
180 h0.SetBinContent(ix, W);
182 catch(
const exception&) {}