34int main(
int argc,
char **argv)
39 string fileDescriptor;
54 zap[
'O'] =
make_field(option) =
"KM3NeT",
"Antares";
55 zap[
'E'] =
make_field(E,
"muon/shower energy [GeV]");
56 zap[
'c'] =
make_field(cd,
"cosine emission angle");
57 zap[
'D'] =
make_field(dir,
"(theta, phi) of PMT [rad]");
58 zap[
'T'] =
make_field(TMax_ns,
"L1 coincidence time window [ns]") = 20.0;
65 catch(
const exception &error) {
66 FATAL(error.what() << endl);
70 const double epsilon = 1.0e-6;
78 if (option ==
"KM3NeT") {
112 }
else if (option ==
"Antares") {
120 FATAL(
"Fatal error at detector option " << option << endl);
136 TH1D h0m(
"L0m", NULL, 300, 1.0, 151.0);
137 TH1D h1m(
"L1m", NULL, 300, 1.0, 151.0);
139 TH1D h0s(
"L0s", NULL, 300, 1.0, 151.0);
140 TH1D h1s(
"L1s", NULL, 300, 1.0, 151.0);
153 const JPDFType_t type[] = { DIRECT_LIGHT_FROM_MUON,
154 SCATTERED_LIGHT_FROM_MUON,
155 DIRECT_LIGHT_FROM_EMSHOWERS,
156 SCATTERED_LIGHT_FROM_EMSHOWERS,
157 DIRECT_LIGHT_FROM_DELTARAYS,
158 SCATTERED_LIGHT_FROM_DELTARAYS };
160 const double TTS = 2.0;
162 const double epsilon = 1.0e-10;
164 const int NUMBER_OF_PDFS =
sizeof(type)/
sizeof(type[0]);
166 JPDF_t pdf[NUMBER_OF_PDFS];
168 const JPDF_t::JSupervisor supervisor(
new JPDF_t::JDefaultResult(
JMATH::zero));
170 for (
int i = 0; i != NUMBER_OF_PDFS; ++i) {
174 const string file_name = getFilename(fileDescriptor, type[i]);
176 NOTICE(
"loading PDF from file " << file_name <<
"... " << flush);
178 pdf[i].load(file_name.c_str());
182 pdf[i].setExceptionHandler(supervisor);
190 const double Tmin = -15.0;
191 const double Tmax = +250.0;
192 const double dt = 1.0;
194 for (
int ix = 1; ix <= h0m.GetNbinsX(); ++ix) {
196 const double R = h0m.GetBinCenter(ix);
200 for (vector<JAngle3D>::const_iterator pmt =
PMT.begin(); pmt !=
PMT.end(); ++pmt) {
202 const double theta = pi.
constrain(pmt->getTheta());
203 const double phi = pi.
constrain(fabs(pmt->getPhi()));
205 for (
int i = 0; i != NUMBER_OF_PDFS; ++i) {
207 double yi = pdf[i](R, theta, phi, Tmax).V;
209 if (is_bremsstrahlung(type[i])) {
211 }
else if (is_deltarays(type[i])) {
212 yi *= getDeltaRaysFromMuon(E);
218 h0m.SetBinContent(ix, 1.0 - TMath::PoissonI(0,V));
222 const int NUMBER_OF_PMTS =
PMT.size();
224 double Vi[NUMBER_OF_PMTS];
226 for (
int ix = 1; ix <= h1m.GetNbinsX(); ++ix) {
228 const double R = h1m.GetBinCenter(ix);
232 for (
double x = Tmin; x <= Tmax; x += dt) {
236 for (
int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
241 const double phi = pi.
constrain(fabs(
PMT[pmt].getPhi()));
243 for (
int i = 0; i != NUMBER_OF_PDFS; ++i) {
246 pdf[i](R, theta, phi, x).v,
247 pdf[i](R, theta, phi, x + TMax_ns).v
250 if (is_bremsstrahlung(type[i])) {
253 }
else if (is_deltarays(type[i])) {
254 yi[0] *= getDeltaRaysFromMuon(E);
255 yi[1] *= getDeltaRaysFromMuon(E);
258 Vi[pmt] += yi[1] - yi[0];
264 for (
int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
267 const double phi = pi.
constrain(fabs(
PMT[pmt].getPhi()));
271 for (
int i = 0; i != NUMBER_OF_PDFS; ++i) {
273 double yi = pdf[i](R, theta, phi, x).f;
275 if (is_bremsstrahlung(type[i])) {
277 }
else if (is_deltarays(type[i])) {
278 yi *= getDeltaRaysFromMuon(E);
285 Y += y * (1.0 - TMath::PoissonI(0, V - Vi[pmt])) * dt;
290 h1m.SetBinContent(ix, 1.0 - TMath::PoissonI(0,Y));
306 const JPDFType_t type[] = { DIRECT_LIGHT_FROM_EMSHOWER,
307 SCATTERED_LIGHT_FROM_EMSHOWER };
309 const double TTS = 2.0;
311 const double epsilon = 1.0e-10;
313 const int NUMBER_OF_PDFS =
sizeof(type)/
sizeof(type[0]);
315 JPDF_t pdf[NUMBER_OF_PDFS];
317 const JPDF_t::JSupervisor supervisor(
new JPDF_t::JDefaultResult(
JMATH::zero));
319 for (
int i = 0; i != NUMBER_OF_PDFS; ++i) {
323 const string file_name = getFilename(fileDescriptor, type[i]);
325 NOTICE(
"loading PDF from file " << file_name <<
"... " << flush);
327 pdf[i].load(file_name.c_str());
331 pdf[i].setExceptionHandler(supervisor);
339 const double Tmin = -15.0;
340 const double Tmax = +250.0;
341 const double dt = 1.0;
343 for (
int ix = 1; ix <= h0s.GetNbinsX(); ++ix) {
345 const double D = h0s.GetBinCenter(ix);
349 for (vector<JAngle3D>::const_iterator pmt =
PMT.begin(); pmt !=
PMT.end(); ++pmt) {
351 const double theta = pi.
constrain(pmt->getTheta());
352 const double phi = pi.
constrain(fabs(pmt->getPhi()));
354 for (
int i = 0; i != NUMBER_OF_PDFS; ++i) {
356 double yi = pdf[i](D, cd, theta, phi, Tmax).V;
364 h0s.SetBinContent(ix, 1.0 - TMath::PoissonI(0,V));
368 const int NUMBER_OF_PMTS =
PMT.size();
370 double Vi[NUMBER_OF_PMTS];
372 for (
int ix = 1; ix <= h1s.GetNbinsX(); ++ix) {
374 const double D = h1s.GetBinCenter(ix);
378 for (
double x = Tmin; x <= Tmax; x += dt) {
382 for (
int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
387 const double phi = pi.
constrain(fabs(
PMT[pmt].getPhi()));
389 for (
int i = 0; i != NUMBER_OF_PDFS; ++i) {
392 pdf[i](D, cd, theta, phi, x).v,
393 pdf[i](D, cd, theta, phi, x + TMax_ns).v
399 Vi[pmt] += yi[1] - yi[0];
405 for (
int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
408 const double phi = pi.
constrain(fabs(
PMT[pmt].getPhi()));
412 for (
int i = 0; i != NUMBER_OF_PDFS; ++i) {
414 double yi = pdf[i](D, cd, theta, phi, x).f;
422 Y += y * (1.0 - TMath::PoissonI(0, V - Vi[pmt])) * dt;
427 h1s.SetBinContent(ix, 1.0 - TMath::PoissonI(0,Y));