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