48 int main(
int argc,
char **argv)
53 typedef JAbstractHistogram<double> JHistogram_t;
62 JHistogram_t histogram;
71 JParser<> zap(
"Auxiliary program to draw PDF of Cherenkov light from muon.");
76 zap[
'e'] =
make_field(epsilon,
"precision for integration") = 1.0e-10;
79 zap[
'E'] =
make_field(E,
"muon energy [GeV]") = 1.0;
80 zap[
'R'] =
make_field(R,
"distance of approach [m]");
81 zap[
'D'] =
make_field(dir,
"(theta, phi) of PMT [rad]");
83 zap[
'H'] =
make_field(histogram,
"histogram binning") = JHistogram_t();
88 catch(
const exception &error) {
89 FATAL(error.what() << endl);
109 for (
double dt; cin >> dt; ) {
113 cout << setw(2) << *F <<
' '
115 <<
FIXED(5,1) << R <<
' '
116 <<
FIXED(5,2) << dir.getTheta() <<
' '
117 <<
FIXED(5,2) << dir.getPhi() <<
' '
118 <<
FIXED(5,1) << dt <<
' '
119 <<
SCIENTIFIC(9,3) << pdf.getLightFromMuon(*F, E, R, dir.getTheta(), dir.getPhi(), dt) * E << endl;
131 const double t0 = 0.0;
133 if (!histogram.is_valid()) {
137 histogram = JHistogram_t(t0 - 20.0, t0 + 50.0);
139 histogram.setBinWidth(0.1);
143 histogram = JHistogram_t(t0 - 20.0, t0 + 500.0);
145 histogram.setBinWidth(0.5);
149 TH1D h0(
"h0", NULL, histogram.getNumberOfBins(), histogram.getLowerLimit(), histogram.getUpperLimit());
151 JSplineFunction1S_t f1;
153 for (
int i = 1; i <= h0.GetNbinsX(); ++i) {
155 const double dt = h0.GetBinCenter(i) - t0;
160 value += pdf.getLightFromMuon(*F, E, R, dir.getTheta(), dir.getPhi(), dt);
163 h0.SetBinContent(i, value);
168 f1.setExceptionHandler(
new JSplineFunction1S_t::JDefaultResult(
JMATH::zero));
173 JQuantiles quantiles(f1);
175 const double t1 = quantiles.getX();
176 const double y = f1(t1 + T).v - f1(t1 - T).v;
178 DEBUG(
"E " << E << endl);
179 DEBUG(
"R " << R << endl);
180 DEBUG(
"theta " << dir.getTheta() << endl);
181 DEBUG(
"phi " << dir.getPhi() << endl);
182 DEBUG(
"int " << quantiles.getIntegral() << endl);
183 DEBUG(
"t1 " << t1 << endl);
184 DEBUG(
"max " << quantiles.getY() << endl);
185 DEBUG(
"FWHM " << quantiles.getFWHM() << endl);
186 DEBUG(
"int[] " << y << endl);