Jpp 19.3.0-rc.2
the software that should make you happy
Loading...
Searching...
No Matches
JHondaFluxInterpolator.hh
Go to the documentation of this file.
1#ifndef __JAANET__JHONDAFLUXINTERPOLATOR__
2#define __JAANET__JHONDAFLUXINTERPOLATOR__
3
4#include <iostream>
5#include <string>
6#include <regex>
7
8#include "Jeep/JPrint.hh"
9#include "Jeep/JMessage.hh"
10
11#include "JLang/JClonable.hh"
12#include "JLang/JException.hh"
15
16#include "JTools/JRange.hh"
17#include "JTools/JPolint.hh"
18#include "JTools/JMapList.hh"
19#include "JTools/JCollection.hh"
24
26
27
28/**
29 * \author bjung
30 */
31
32namespace JAANET {}
33namespace JPP { using namespace JAANET; }
34
35namespace JAANET {
36
37 using JEEP::JMessage;
38
39 using JLANG::JClonable;
42
43 using JTOOLS::JRange;
44 using JTOOLS::JArray;
45
46 using JTOOLS::JMAPLIST;
47
51
52
53 static const char* const HONDA_ENERGY_KEYWORD = "Enu(GeV)";
54 static const char* const HONDA_COSINE_ZENITH_ANGLE_KEYWORD = "cosZ";
55 static const char* const HONDA_AZIMUTH_ANGLE_KEYWORD = "phi_Az";
56
57 static const char* const HONDA_MUON_NEUTRINO_FLUX_KEYWORD = "NuMu";
58 static const char* const HONDA_MUON_ANTINEUTRINO_FLUX_KEYWORD = "NuMubar";
59 static const char* const HONDA_ELECTRON_NEUTRINO_FLUX_KEYWORD = "NuE";
60 static const char* const HONDA_ELECTRON_ANTINEUTRINO_FLUX_KEYWORD = "NuEbar";
61
62 static const char* const HONDA_BINSPECS_SEPARATOR = "--";
63 static const char HONDA_BINSPECS_BEGIN = '[';
64 static const char HONDA_BINSPECS_END = ']';
65
66
67 /**
68 * Auxiliary data structure for reading Honda bin ranges.
69 */
71 JRange<double>
72 {
74
75
76 /**
77 * Default constructor.
78 */
81
82
83 /**
84 * Constructor.
85 *
86 * \param minimum minimum value
87 * \param maximum maximum value
88 */
89 JHondaBinRange(const double minimum,
90 const double maximum) :
91 JRange_t(minimum, maximum)
92 {}
93
94
95 /**
96 * Get central value of this range.
97 *
98 * \return central value
99 */
100 double getCentralValue() const
101 {
102 return (getLowerLimit() + getUpperLimit()) / 2.0;
103 }
104
105
106 /**
107 * Read bin range from input.
108 *
109 * \param in input stream
110 * \param object bin range
111 * \return input stream
112 */
113 friend inline std::istream& operator>>(std::istream& in, JHondaBinRange& object)
114 {
115 using namespace std;
116
117 double min = JRange_t::getMinimum();
118 double max = JRange_t::getMaximum();
119
120 string buffer;
121
122 if (in >> min >> buffer >> max &&
123 min <= max &&
124 buffer == HONDA_BINSPECS_SEPARATOR) {
125 static_cast<JRange_t&>(object) = JRange_t(min, max);
126 }
127
128 return in;
129 }
130 };
131
132
133 /**
134 * Auxiliary data structure for reading angular binning specifications of Honda flux table.
135 */
137 {
138 /**
139 * Default constructor.
140 */
142 coszRange(-1.0, 1.0),
143 phiRange ( 0.0, 360.0)
144 {}
145
146
147 /**
148 * Constructor.
149 *
150 * \param minCosz minimum cosine zenith angle
151 * \param maxCosz maximum cosine zenith angle
152 * \param minPhi minimum azimuth angle
153 * \param maxPhi maximum azimuth angle
154 */
155 JHondaAngularBinSpecs(const double minCosz,
156 const double maxCosz,
157 const double minPhi,
158 const double maxPhi) :
159 coszRange(minCosz, maxCosz),
160 phiRange (minPhi, maxPhi)
161 {}
162
163
164 /**
165 * Read bin specifications from input.
166 *
167 * \param in input stream
168 * \param binspecs bin specifications
169 * \return input stream
170 */
171 friend inline std::istream& operator>>(std::istream& in, JHondaAngularBinSpecs& binspecs)
172 {
173 static const JEquationParameters eqpars("=", ",", "./", "#");
174
175 JProperties properties(eqpars, 1);
176
177 properties[HONDA_COSINE_ZENITH_ANGLE_KEYWORD] = binspecs.coszRange;
178 properties[HONDA_AZIMUTH_ANGLE_KEYWORD] = binspecs.phiRange;
179
180 in >> properties;
181
182 return in >> properties;
183 }
184
185
186 JHondaBinRange coszRange; //!< Cosine zenith-angle range
187 JHondaBinRange phiRange; //!< Azimuthal angle range [deg]
188 };
189
190
191 /**
192 * Template definition for Honda flux table interpolator.
193 */
194 template<class JPhiFunction_t = JConstantFunction1D<double, JArray<4, double> >,
195 template<class, class, class> class JCoszFunctionalMap_t = JPolint1FunctionalMap,
196 template<class, class, class> class JEnergyFunctionalMap_t = JPolint1FunctionalMap>
198 public JMultiFunction <JPhiFunction_t, typename JMAPLIST<JEnergyFunctionalMap_t, JCoszFunctionalMap_t>::maplist>,
199 public JClonable<JFlux, JHondaFluxInterpolator <JPhiFunction_t, JCoszFunctionalMap_t, JEnergyFunctionalMap_t> >,
200 public JObjectStreamIO <JHondaFluxInterpolator <JPhiFunction_t, JCoszFunctionalMap_t, JEnergyFunctionalMap_t> >,
201 public JMessage <JHondaFluxInterpolator <JPhiFunction_t, JCoszFunctionalMap_t, JEnergyFunctionalMap_t> >
202 {
203 public:
204
206
209
211
216
218
222
224
225
226 /**
227 * Default constructor.
228 */
230 table (*this),
231 coszRange(0.0, -1.0),
232 phiRange (0.0, -1.0)
233 {}
234
235
236 /**
237 * Constructor.
238 *
239 * \param filename Honda flux table filename
240 */
241 JHondaFluxInterpolator(const char* filename) :
242 table (*this, filename),
243 coszRange(0.0, -1.0),
244 phiRange (0.0, -1.0)
245 {
246 load();
247 }
248
249
250 /**
251 * Copy constructor.
252 *
253 * \param interpolator Honda flux interpolator
254 */
256 table (*this, interpolator.table.c_str()),
257 coszRange(interpolator.coszRange)
258 {
259 this->getMultiFunction() = interpolator.getMultiFunction();
260
261 this->compile();
262 }
263
264
265 /**
266 * Check whether this interpolator is valid.
267 *
268 * \return true if valid; else false
269 */
270 bool is_valid() const override final
271 {
272 return (!table.empty() ?
273 this->size() > 0 && coszRange.getLength() > 0.0 && phiRange.getLength() > 0.0 :
274 true);
275 }
276
277
278 /**
279 * Get diffuse flux for given particle PDG-identifier, energy, cosine zenith angle and azimuth angle.
280 *
281 * \param type PDG particle type
282 * \param log10E logarithmic neutrino energy [GeV]
283 * \param costh cosine zenith angle
284 * \param phi azimuth angle (defined clockwise w.r.t. the North) [deg]
285 * \return diffuse flux [GeV^-1 * m^-2 * sr^-1 * s^-1]
286 */
287 double getFactor(const int type,
288 const double log10E,
289 const double costh,
290 const double phi) const
291 {
292 using namespace std;
293 using namespace JPP;
294
295 static const double epsilon_costh = 1.0e-3; //!< [--]
296 static const double epsilon_phi = 1.0; //!< [deg]
297
298 static const JRange_t mod(0.0, 360.0);
299
300 const multifunction_type& function = static_cast<const multifunction_type&>(*this);
301
302 const double _costh = min(max(costh, -1.0), 1.0);
303 const double _phi = mod.mod(180.0 - phi); //!< Honda azimuth angle is defined counterclockwise w.r.t. the South\n
304 //!< (c.f. arXiv:1102.2688), so a conversion is necessary.
305
306 result_type lnFluxes;
307
308 if (coszRange(_costh) && phiRange(_phi)) {
309
310 lnFluxes = function(log10E, _costh, _phi);
311
312 } else { // Perform linear extrapolation for edges of angular ranges
313
314 if (_costh > coszRange.getUpperLimit() && phiRange(_phi)) {
315
317 const abscissa_type x2 = coszRange.getUpperLimit() - epsilon_costh;
318
319 const result_type& y1 = function(log10E, x1, _phi);
320 const result_type& y2 = function(log10E, x2, _phi);
321
322 lnFluxes = y1 + ((y2 - y1) / (x2 - x1)) * (_costh - x1);
323
324 } else if (_costh < coszRange.getLowerLimit() && phiRange(_phi)) {
325
327 const abscissa_type x2 = coszRange.getLowerLimit() + epsilon_costh;
328
329 const result_type& y1 = function(log10E, x1, _phi);
330 const result_type& y2 = function(log10E, x2, _phi);
331
332 lnFluxes = y1 + ((y2 - y1) / (x2 - x1)) * (_costh - x1);
333
334 } else if (coszRange(_costh) && _phi > phiRange.getUpperLimit()) {
335
337 const abscissa_type x2 = phiRange.getUpperLimit() - epsilon_costh;
338
339 const result_type& y1 = function(log10E, _costh, x1);
340 const result_type& y2 = function(log10E, _costh, x2);
341
342 lnFluxes = y1 + ((y2 - y1) / (x2 - x1)) * (_phi - x1);
343
344 } else if (coszRange(_costh) && _phi < phiRange.getLowerLimit()) {
345
347 const abscissa_type x2 = phiRange.getLowerLimit() + epsilon_costh;
348
349 const result_type& y1 = function(log10E, _costh, x1);
350 const result_type& y2 = function(log10E, _costh, x2);
351
352 lnFluxes = y1 + ((y2 - y1) / (x2 - x1)) * (_phi - x1);
353
354 } else if (_costh > coszRange.getUpperLimit() && _phi > phiRange.getUpperLimit()) {
355
357 const abscissa_type c2 = coszRange.getUpperLimit() - epsilon_costh;
358
360 const abscissa_type p2 = phiRange .getUpperLimit() - epsilon_phi;
361
362 const result_type& y1 = function(log10E, c1, p1);
363 const result_type& y2 = function(log10E, c2, p2);
364
365 lnFluxes = y1 + ((y2 - y1) / (c2 - c1) / (p2 - p1)) * (_costh - p1) * (_phi - p1);
366
367 } else if (_costh > coszRange.getUpperLimit() && _phi < phiRange.getUpperLimit()) {
368
370 const abscissa_type c2 = coszRange.getUpperLimit() - epsilon_costh;
371
373 const abscissa_type p2 = phiRange .getLowerLimit() + epsilon_phi;
374
375 const result_type& y1 = function(log10E, c1, p1);
376 const result_type& y2 = function(log10E, c2, p2);
377
378 lnFluxes = y1 + ((y2 - y1) / (c2 - c1) / (p2 - p1)) * (_costh - p1) * (_phi - p1);
379
380 } else if (_costh < coszRange.getLowerLimit() && phi > phiRange.getUpperLimit()) {
381
383 const abscissa_type c2 = coszRange.getLowerLimit() + epsilon_costh;
384
386 const abscissa_type p2 = phiRange .getUpperLimit() - epsilon_phi;
387
388 const result_type& y1 = function(log10E, c1, p1);
389 const result_type& y2 = function(log10E, c2, p2);
390
391 lnFluxes = y1 + ((y2 - y1) / (c2 - c1) / (p2 - p1)) * (_costh - p1) * (_phi - p1);
392
393 } else {
394
396 const abscissa_type c2 = coszRange.getLowerLimit() + epsilon_costh;
397
399 const abscissa_type p2 = phiRange .getLowerLimit() + epsilon_phi;
400
401 const result_type& y1 = function(log10E, c1, p1);
402 const result_type& y2 = function(log10E, c2, p2);
403
404 lnFluxes = y1 + ((y2 - y1) / (c2 - c1) / (p2 - p1)) * (_costh - p1) * (_phi - p1);
405 }
406 }
407
408 switch(type) {
409 case TRACK_TYPE_NUMU:
410 return exp(lnFluxes[0]);
412 return exp(lnFluxes[1]);
413 case TRACK_TYPE_NUE:
414 return exp(lnFluxes[2]);
416 return exp(lnFluxes[3]);
417 default:
418 return 0.0;
419 }
420 }
421
422
423 /**
424 * Get flux of given event.
425 *
426 * \param type PDG particle type
427 * \param log10E logarithmic neutrino energy [GeV]
428 * \param costh cosine zenith angle
429 * \param phi azimuth angle (defined clockwise w.r.t. the North) [deg]
430 * \return diffuse flux [GeV^-1 * m^-2 * sr^-1 * s^-1]
431 */
432 double getFlux(const int type,
433 const double log10E,
434 const double costh,
435 const double phi) const
436 {
437 return getFactor(type, log10E, costh, phi);
438 }
439
440
441 /**
442 * Get diffuse flux for given particle PDG-identifier, energy, cosine zenith angle and azimuth angle.
443 *
444 * \param type PDG particle type
445 * \param log10E logarithmic neutrino energy [GeV]
446 * \param costh cosine zenith angle
447 * \param phi azimuth angle (defined clockwise w.r.t. the North) [deg]
448 * \return diffuse flux [GeV^-1 * m^-2 * sr^-1 * s^-1]
449 */
450 double operator()(const int type,
451 const double log10E,
452 const double costh,
453 const double phi) const
454 {
455 return getFactor(type, log10E, costh, phi);
456 }
457
458
459 /**
460 * Get flux of given event.
461 *
462 * \param evt event
463 * \return flux [GeV^-1 * m^-2 * sr^-1 * s^-1]
464 */
465 double getFactor(const Evt& evt) const override final
466 {
467 const Trk& neutrino = get_neutrino(evt);
468
469 const double log10E = log10(neutrino.E);
470 const double costh = neutrino.dir.z;
471
472 //!< The gSeaGen azimuth angle is defined clockwise w.r.t. the North (c.f. arXiv:2003.14040)
473 const double phi = atan(neutrino.dir.y / neutrino.dir.x) * 180.0 / M_PI;
474
475 return getFactor(neutrino.type, log10E, costh, phi);
476 }
477
478
479 /**
480 * Get flux of given event.
481 *
482 * \param event event
483 * \return flux [GeV^-1 * m^-2 * sr^-1 * s^-1]
484 */
485 double getFlux(const Evt& event) const
486 {
487 return getFactor(event);
488 }
489
490
491 /**
492 * Load Honda flux table.
493 */
494 void load()
495 {
496 using namespace std;
497 using namespace JPP;
498
499 if (!table.empty()) {
500
501 NOTICE("Loading Honda flux table from file " << table << "... " << flush);
502
504
505 this->check_validity();
506
507 NOTICE("OK" << endl);
508 }
509 }
510
511
512 /**
513 * Load Honda flux table.
514 *
515 * \param table Honda flux table filename
516 */
517 void load(const char* table)
518 {
519 this->table = table;
520
521 load();
522 }
523
524
525 /**
526 * Get properties of this class.
527 *
528 * \param eqpars equation parameters
529 */
531 {
532 return JHondaFluxInterpolatorHelper(*this, eqpars);
533 }
534
535
536 /**
537 * Get properties of this class.
538 *
539 * \param eqpars equation parameters
540 */
542 {
543 return JHondaFluxInterpolatorHelper(*this, eqpars);
544 }
545
546
547 /**
548 * Read Honda atmospheric neutrino flux interpolator from file.
549 *
550 * \param ifs input file stream
551 * \param interpolator Honda flux interpolator
552 * \return input file stream
553 */
554 friend std::ifstream& operator>>(std::ifstream& ifs, interpolator_type& interpolator)
555 {
556 using namespace std;
557
558 multimap_type& zmap = static_cast<multimap_type&>(interpolator);
559
560 JHondaAngularBinSpecs binspecs;
561
562 for (string line; getline(ifs, line); ) {
563
564 double energy = 0.0;
565
566 result_type values;
567
568 istringstream iss1(line);
569
570 const double cosz = binspecs.coszRange.getCentralValue();
571 const double phi = binspecs.phiRange .getCentralValue();
572
573 interpolator.coszRange.include(cosz);
574 interpolator.phiRange .include(phi);
575
576 if (iss1 >> energy >> values) {
577
578 const double log10E = log10(energy);
579
580 for (typename result_type::iterator i = values.begin(); i != values.end(); ++i) {
581 *i = log(*i);
582 }
583
584 zmap[log10E][cosz][phi] = values;
585
586 } else {
587
588 const size_t p1 = line.find (HONDA_BINSPECS_BEGIN);
589 const size_t p2 = line.rfind(HONDA_BINSPECS_END);
590
591 if (p1 != string::npos && p2 != string::npos) {
592
593 const string binspecstr = line.substr(p1 + 1, p2 - p1 - 1);
594
595 istringstream iss2(binspecstr);
596
597 iss2 >> binspecs;
598 }
599 }
600 }
601
602 interpolator.compile();
603
604 return ifs;
605 }
606
607
608 private:
609
610 /**
611 * Auxiliary class for I/O of Honda flux interpolator.
612 */
614 public JProperties
615 {
616 /**
617 * Constructor.
618 *
619 * \param interpolator Honda flux interpolator
620 * \param eqpars equation parameters
621 */
622 template<class JHondaFluxInterpolator_t>
623 JHondaFluxInterpolatorHelper(JHondaFluxInterpolator_t& interpolator,
624 const JEquationParameters& eqpars) :
625 JProperties(eqpars, 1)
626 {
627 (*this)[JEvtWeightFactor::getTypeKey()] = "Honda flux interpolator";
628
629 this->insert(gmake_property(interpolator.table));
630 }
631 };
632
633
635
636 JRange_t coszRange; //!< Zenith angle range
637 JRange_t phiRange; //!< Azimuth angle range [deg]
638
639 using JMessage<interpolator_type>::debug; //!< Debug level
640 };
641
642
643 /**
644 * Template specialisation for interpolator of azimuth-averaged Honda flux table.
645 */
646 template<template<class, class, class> class JCoszFunctionalMap_t,
647 template<class, class, class> class JEnergyFunctionalMap_t>
648 class JHondaFluxInterpolator<JConstantFunction1D<double, JArray<4, double> >, JCoszFunctionalMap_t, JEnergyFunctionalMap_t> final :
649 public JMultiFunction <JConstantFunction1D<double, JArray<4, double> >, typename JMAPLIST<JEnergyFunctionalMap_t, JCoszFunctionalMap_t>::maplist>,
650 public JClonable<JDiffuseFlux, JHondaFluxInterpolator <JConstantFunction1D<double, JArray<4, double> >, JCoszFunctionalMap_t, JEnergyFunctionalMap_t> >,
651 public JObjectStreamIO <JHondaFluxInterpolator <JConstantFunction1D<double, JArray<4, double> >, JCoszFunctionalMap_t, JEnergyFunctionalMap_t> >,
652 public JMessage <JHondaFluxInterpolator <JConstantFunction1D<double, JArray<4, double> >, JCoszFunctionalMap_t, JEnergyFunctionalMap_t> >
653 {
654 public:
655
657
659
662
664
669
671
675
677
678
679 /**
680 * Default constructor.
681 */
683 table (*this),
684 coszRange(0.0, -1.0)
685 {}
686
687
688 /**
689 * Constructor.
690 *
691 * \param filename Honda flux table filename
692 */
693 JHondaFluxInterpolator(const char* filename) :
694 table (*this, filename),
695 coszRange(0.0, -1.0)
696 {
697 load();
698 }
699
700
701 /**
702 * Copy constructor.
703 *
704 * \param interpolator Honda flux interpolator
705 */
707 table (*this, interpolator.table.c_str()),
708 coszRange(interpolator.coszRange)
709 {
710 this->getMultiFunction() = interpolator.getMultiFunction();
711
712 this->compile();
713 }
714
715
716 /**
717 * Check whether this interpolator is valid.
718 *
719 * \return true if valid; else false
720 */
721 bool is_valid() const override final
722 {
723 return (!table.empty() ? this->size() > 0 && coszRange.getLength() : true);
724 }
725
726
727 /**
728 * Get diffuse flux for given particle PDG-identifier, energy and zenith-angle.
729 *
730 * \param type PDG particle type
731 * \param log10E logarithmic neutrino energy [GeV]
732 * \param costh cosine zenith angle
733 * \return diffuse flux [GeV^-1 * m^-2 * sr^-1 * s^-1]
734 */
735 double dNdEdOmega(const int type,
736 const double log10E,
737 const double costh) const override final
738 {
739 using namespace std;
740 using namespace JPP;
741
742 static const double epsilon_costh = 1.0e-3;
743
744 const multifunction_type& function = static_cast<const multifunction_type&>(*this);
745
746 const double _costh = min(max(costh, -1.0), 1.0);
747
748 result_type lnFluxes;
749
750 if (coszRange(_costh)) {
751
752 lnFluxes = function(log10E, _costh);
753
754 } else { // Perform linear extrapolation for edges of cosine zenith angle range
755
756 if (_costh > coszRange.getUpperLimit()) {
757
759 const abscissa_type x2 = coszRange.getUpperLimit() - epsilon_costh;
760
761 const result_type& y1 = function(log10E, x1);
762 const result_type& y2 = function(log10E, x2);
763
764 lnFluxes = y1 + ((y2 - y1) / (x2 - x1)) * (_costh - x1);
765
766 } else {
767
769 const abscissa_type x2 = coszRange.getLowerLimit() + epsilon_costh;
770
771 const result_type& y1 = function(log10E, x1);
772 const result_type& y2 = function(log10E, x2);
773
774 lnFluxes = y1 + ((y2 - y1) / (x2 - x1)) * (_costh - x1);
775 }
776 }
777
778 switch(type) {
779 case TRACK_TYPE_NUMU:
780 return exp(lnFluxes[0]);
782 return exp(lnFluxes[1]);
783 case TRACK_TYPE_NUE:
784 return exp(lnFluxes[2]);
786 return exp(lnFluxes[3]);
787 default:
788 return 0.0;
789 }
790 }
791
792
793 /**
794 * Load Honda flux table.
795 */
796 void load()
797 {
798 using namespace std;
799 using namespace JPP;
800
801 if (!table.empty()) {
802
803 NOTICE("Loading Honda flux table from file " << table << "... " << flush);
804
806
807 this->check_validity();
808
809 NOTICE("OK" << endl);
810 }
811 }
812
813
814 /**
815 * Load Honda flux table.
816 *
817 * \param table Honda flux table filename
818 */
819 void load(const char* table)
820 {
821 this->table = table;
822
823 load();
824 }
825
826
827 /**
828 * Get properties of this class.
829 *
830 * \param eqpars equation parameters
831 */
833 {
834 return JHondaFluxInterpolatorHelper(*this, eqpars);
835 }
836
837
838 /**
839 * Get properties of this class.
840 *
841 * \param eqpars equation parameters
842 */
844 {
845 return JHondaFluxInterpolatorHelper(*this, eqpars);
846 }
847
848
849 /**
850 * Read Honda atmospheric neutrino flux interpolator from file.
851 *
852 * \param ifs input file stream
853 * \param interpolator Honda flux interpolator
854 * \return input file stream
855 */
856 friend std::ifstream& operator>>(std::ifstream& ifs, interpolator_type& interpolator)
857 {
858 using namespace std;
859
860 multimap_type& zmap = static_cast<multimap_type&>(interpolator);
861
862 JHondaAngularBinSpecs binspecs;
863
864 for (string line; getline(ifs, line); ) {
865
866 double energy = 0.0;
867
868 result_type values;
869
870 istringstream iss1(line);
871
872 const double cosz = binspecs.coszRange.getCentralValue();
873
874 interpolator.coszRange.include(cosz);
875
876 if (iss1 >> energy >> values) {
877
878 const double log10E = log10(energy);
879
880 for (typename result_type::iterator i = values.begin(); i != values.end(); ++i) {
881 *i = log(*i);
882 }
883
884 zmap[log10E][cosz] = values;
885
886 } else {
887
888 const size_t p1 = line.find (HONDA_BINSPECS_BEGIN);
889 const size_t p2 = line.rfind(HONDA_BINSPECS_END);
890
891 if (p1 != string::npos && p2 != string::npos) {
892
893 const string binspecstr = line.substr(p1 + 1, p2 - p1 - 1);
894
895 istringstream iss2(binspecstr);
896
897 iss2 >> binspecs;
898 }
899 }
900 }
901
902 interpolator.compile();
903
904 return ifs;
905 }
906
907
908 private:
909
910 /**
911 * Auxiliary class for I/O of Honda flux interpolator.
912 */
914 public JProperties
915 {
916 /**
917 * Constructor.
918 *
919 * \param interpolator Honda flux interpolator
920 * \param eqpars equation parameters
921 */
922 template<class JHondaFluxInterpolator_t>
923 JHondaFluxInterpolatorHelper(JHondaFluxInterpolator_t& interpolator,
924 const JEquationParameters& eqpars) :
925 JProperties(eqpars, 1)
926 {
927 (*this)[JEvtWeightFactor::getTypeKey()] = "Honda flux interpolator";
928
929 this->insert(gmake_property(interpolator.table));
930 }
931 };
932
933
935
936 JRange_t coszRange; //!< Cosine zenith angle range
937
938 using JMessage<interpolator_type>::debug; //!< Debug level
939 };
940
941
942 /**
943 * Alias template definition for 2-dimensional Honda flux interpolator.
944 */
945 template<template<class, class, class> class JCoszFunctionalMap_t = JPolint1FunctionalMap,
946 template<class, class, class> class JEnergyFunctionalMap_t = JPolint1FunctionalMap>
948 JCoszFunctionalMap_t,
949 JEnergyFunctionalMap_t>;
950}
951
952
953namespace JEEP {
954
955 /**
956 * JMessage template specialization for Honda atmospheric neutrino flux interpolators.
957 */
958 template<class JPhiFunction_t,
959 template<class, class, class> class JCoszFunctionalMap_t,
960 template<class, class, class> class JEnergyFunctionalMap_t>
961 struct JMessage<JAANET::JHondaFluxInterpolator<JPhiFunction_t, JCoszFunctionalMap_t, JEnergyFunctionalMap_t> >
962 {
963 static int debug;
964 };
965
966
967 /**
968 * Default verbosity for Honda atmospheric neutrino flux interpolators.
969 */
970 template<class JPhiFunction_t,
971 template<class, class, class> class JCoszFunctionalMap_t,
972 template<class, class, class> class JEnergyFunctionalMap_t>
974
975}
976
977#endif
General purpose class for a collection of sorted elements.
TCanvas * c1
Global variables to handle mouse events.
TPaveText * p1
Exceptions.
Various implementations of functional maps.
General purpose messaging.
#define NOTICE(A)
Definition JMessage.hh:64
int debug
debug level
Definition JSirene.cc:72
I/O formatting auxiliaries.
#define gmake_property(A)
macros to convert (template) parameter to JPropertiesElement object
Auxiliary class to define a range between two values.
double dNdEdOmega(const int type, const double log10E, const double costh) const override final
Get diffuse flux for given particle PDG-identifier, energy and zenith-angle.
JProperties getProperties(const JEquationParameters &eqpars=JEvtWeightFactor::getEquationParameters()) const override final
Get properties of this class.
friend std::ifstream & operator>>(std::ifstream &ifs, interpolator_type &interpolator)
Read Honda atmospheric neutrino flux interpolator from file.
JProperties getProperties(const JEquationParameters &eqpars=JEvtWeightFactor::getEquationParameters()) override final
Get properties of this class.
Template definition for Honda flux table interpolator.
multifunction_type::super_iterator super_iterator
bool is_valid() const override final
Check whether this interpolator is valid.
multifunction_type::multimap_type multimap_type
multifunction_type::function_type function_type
JMAPLIST< JEnergyFunctionalMap_t, JCoszFunctionalMap_t >::maplist JFunctionalMaplist_t
double getFlux(const Evt &event) const
Get flux of given event.
JHondaFluxInterpolator(const interpolator_type &interpolator)
Copy constructor.
double getFactor(const Evt &evt) const override final
Get flux of given event.
void load()
Load Honda flux table.
double getFactor(const int type, const double log10E, const double costh, const double phi) const
Get diffuse flux for given particle PDG-identifier, energy, cosine zenith angle and azimuth angle.
multifunction_type::result_type result_type
JProperties getProperties(const JEquationParameters &eqpars=JEvtWeightFactor::getEquationParameters()) override final
Get properties of this class.
double getFlux(const int type, const double log10E, const double costh, const double phi) const
Get flux of given event.
friend std::ifstream & operator>>(std::ifstream &ifs, interpolator_type &interpolator)
Read Honda atmospheric neutrino flux interpolator from file.
multifunction_type::abscissa_type abscissa_type
multifunction_type::super_const_iterator super_const_iterator
JRange_t phiRange
Azimuth angle range [deg].
JHondaFluxInterpolator< JPhiFunction_t, JCoszFunctionalMap_t, JEnergyFunctionalMap_t > interpolator_type
JProperties getProperties(const JEquationParameters &eqpars=JEvtWeightFactor::getEquationParameters()) const override final
Get properties of this class.
multifunction_type::value_type value_type
JHondaFluxInterpolator(const char *filename)
Constructor.
JRange_t coszRange
Zenith angle range.
void load(const char *table)
Load Honda flux table.
multifunction_type::argument_type argument_type
double operator()(const int type, const double log10E, const double costh, const double phi) const
Get diffuse flux for given particle PDG-identifier, energy, cosine zenith angle and azimuth angle.
JLoadProperty< interpolator_type, std::string > table
Table filename.
JMultiFunction< JPhiFunction_t, JFunctionalMaplist_t > multifunction_type
Utility class to parse parameter values.
Simple data structure to support I/O of equations (see class JLANG::JEquation).
Data structure for object properties which require reloading whenever the property is reread.
One dimensional array of template objects with fixed length.
Definition JArray.hh:43
Template implementation of function object in one dimension returning a constant value.
Multidimensional interpolation method.
multimap_type::super_iterator super_iterator
multimap_type::super_const_iterator super_const_iterator
multimap_type::result_type result_type
void compile()
Compilation.
void insert(const JMultiFunction< __JFunction_t, __JMaplist_t, __JDistance_t > &input)
Insert multidimensional input.
const JMultiFunction & getMultiFunction() const
Get multidimensional function.
function_type::value_type value_type
function_type::argument_type argument_type
multimap_type::abscissa_type abscissa_type
Multidimensional map.
Definition JMultiMap.hh:52
Range of values.
Definition JRange.hh:42
T getLength() const
Get length (difference between upper and lower limit).
Definition JRange.hh:289
range_type & include(argument_type x)
Include given value to range.
Definition JRange.hh:397
T getLowerLimit() const
Get lower limit.
Definition JRange.hh:202
T mod(argument_type x) const
Modulo value with respect to range.
Definition JRange.hh:365
static T getMaximum()
Get maximum possible value.
Definition JRange.hh:545
static T getMinimum()
Get minimum possible value.
Definition JRange.hh:534
T getUpperLimit() const
Get upper limit.
Definition JRange.hh:213
Extensions to Evt data format.
static const char *const HONDA_MUON_ANTINEUTRINO_FLUX_KEYWORD
static const char *const HONDA_MUON_NEUTRINO_FLUX_KEYWORD
static const char *const HONDA_ENERGY_KEYWORD
static const char *const HONDA_ELECTRON_ANTINEUTRINO_FLUX_KEYWORD
static const char *const HONDA_AZIMUTH_ANGLE_KEYWORD
static const char *const HONDA_ELECTRON_NEUTRINO_FLUX_KEYWORD
static const char *const HONDA_COSINE_ZENITH_ANGLE_KEYWORD
static const char HONDA_BINSPECS_END
static const char *const HONDA_BINSPECS_SEPARATOR
static const char HONDA_BINSPECS_BEGIN
const Trk & get_neutrino(const Evt &evt)
Get incoming neutrino.
@ TRACK_TYPE_ANTINUMU
General puprpose classes and methods.
@ notice_t
notice
Definition JMessage.hh:32
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition Evt.hh:21
static const char *const getTypeKey()
Get type keyword.
static JEquationParameters & getEquationParameters()
Get equation parameters.
Auxiliary data structure for reading angular binning specifications of Honda flux table.
JHondaBinRange coszRange
Cosine zenith-angle range.
JHondaBinRange phiRange
Azimuthal angle range [deg].
JHondaAngularBinSpecs(const double minCosz, const double maxCosz, const double minPhi, const double maxPhi)
Constructor.
friend std::istream & operator>>(std::istream &in, JHondaAngularBinSpecs &binspecs)
Read bin specifications from input.
JHondaAngularBinSpecs()
Default constructor.
Auxiliary data structure for reading Honda bin ranges.
friend std::istream & operator>>(std::istream &in, JHondaBinRange &object)
Read bin range from input.
JHondaBinRange()
Default constructor.
JHondaBinRange(const double minimum, const double maximum)
Constructor.
double getCentralValue() const
Get central value of this range.
Auxiliary class for I/O of Honda flux interpolator.
JHondaFluxInterpolatorHelper(JHondaFluxInterpolator_t &interpolator, const JEquationParameters &eqpars)
Constructor.
Auxiliary class for handling debug parameter within a class.
Definition JMessage.hh:44
static int debug
debug level (default is off).
Definition JMessage.hh:45
Template class for object cloning.
Definition JClonable.hh:59
Auxiliary base class for storing and loading a single object to and from an ASCII file,...
void load(const char *file_name)
Load from input file.
Auxiliary class for recursive map list generation.
Definition JMapList.hh:109
Map list.
Definition JMapList.hh:25
Type definition of a 1st degree polynomial interpolation based on a JMap implementation.
The Trk class represents a Monte Carlo (MC) particle as well as a reconstructed track/shower.
Definition Trk.hh:15
int type
MC: particle type in PDG encoding.
Definition Trk.hh:24
Vec dir
track direction
Definition Trk.hh:18
double E
Energy [GeV] (either MC truth or reconstructed)
Definition Trk.hh:20
double z
Definition Vec.hh:14
double x
Definition Vec.hh:14
double y
Definition Vec.hh:14