Jpp  19.0.0
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JEvtWeightToolkit.hh
Go to the documentation of this file.
1 #ifndef __JAANET__JEVTWEIGHTTOOLKIT__
2 #define __JAANET__JEVTWEIGHTTOOLKIT__
3 
4 #include <map>
5 #include <vector>
6 #include <memory>
7 
14 
15 #include "JLang/JException.hh"
16 #include "JLang/JType.hh"
17 
18 #include "Jeep/JProperties.hh"
19 
21 
22 #include "JPhysics/JConstants.hh"
23 
24 #include "JAAnet/JHead.hh"
25 #include "JAAnet/JHeadToolkit.hh"
26 
27 #include "JAAnet/JEvtWeight.hh"
28 #include "JAAnet/JEvtWeightDAQ.hh"
35 
36 #include "JAAnet/JFlatFlux.hh"
37 #include "JAAnet/JPowerLawFlux.hh"
40 #include "JAAnet/JParticleTypes.hh"
41 #include "JAAnet/JAAnetToolkit.hh"
42 
43 
44 /**
45  * \author bjung, mdejong
46  */
47 
48 namespace JAANET {}
49 namespace JPP { using namespace JAANET; }
50 
51 namespace JAANET {
52 
53  using JLANG::JType;
55  using JEEP::JProperties;
56 
57 
58  /**
59  * Look-up table for event weighters.
60  */
61  struct JEvtWeighter :
62  public std::vector<std::shared_ptr<JEvtWeight> >
63  {
64  /**
65  * Constructor
66  */
68  {
69  this->emplace_back(new JEvtWeightGSeaGen());
70  this->emplace_back(new JEvtWeightKM3BUU());
71  this->emplace_back(new JEvtWeightCorsika());
72  this->emplace_back(new JEvtWeightMupage());
73  this->emplace_back(new JEvtWeightGenhen());
74  this->emplace_back(new JEvtWeightDAQ()); // must be after Monte Carlo weighters
75  this->emplace_back(new JEvtWeightMiscellaneous());
76  }
77 
78 
79  /**
80  * Get event weighter corresponding to given header.
81  *
82  * \param header header
83  * \return event weighter
84  */
85  const JEvtWeight& operator()(const JHead& header) const
86  {
87  using namespace JPP;
88 
89  for (const_iterator i = this->begin(); i != this->end(); ++i) {
90 
91  if ((*i)->check(header)) {
92  return *(*i);
93  }
94  }
95 
96  THROW(JValueOutOfRange, "JEvtWeighter::operator(): No event weighter found for given header.");
97  }
98  };
99 
100 
101  /**
102  * Auxiliary class for parsing a vector of neutrino PDG identifiers.
103  */
105  {
106  /**
107  * Default constructor.
108  */
110  identifiers()
111  {}
112 
113 
114  /**
115  * Add identifier.
116  *
117  * \param identifier PDG identifier
118  */
119  void put(const int identifier)
120  {
121  if (abs(identifier) == TRACK_TYPE_NUE ||
122  abs(identifier) == TRACK_TYPE_NUMU ||
123  abs(identifier) == TRACK_TYPE_NUTAU) {
124 
125  identifiers.push_back(identifier);
126 
127  } else {
128 
129  THROW(JValueOutOfRange, "JNeutrinoTypeCollection::put(): Invalid PDG identifier: " << identifier);
130  }
131  }
132 
133 
134  /**
135  * Get size of this collection.
136  *
137  * \return size of neutrino type collection.
138  */
139  size_t size() const
140  {
141  return identifiers.size();
142  }
143 
144 
145  /**
146  * Get constant iterator to the first element of the collection.
147  *
148  * \return constant iterator to the first element of the collection
149  */
151  {
152  return identifiers.cbegin();
153  }
154 
155 
156  /**
157  * Get constant iterator to the last element of the collection.
158  *
159  * \return constant iterator to the last element of the collection
160  */
162  {
163  return identifiers.cend();
164  }
165 
166 
167  /**
168  * Stream input.
169  *
170  * \param in input stream
171  * \param collection collection of neutrino PDG types
172  * \return input stream
173  */
174  friend inline std::istream& operator>>(std::istream& in, JNeutrinoTypeCollection& collection)
175  {
176  for (int identifier; in >> identifier; ) {
177  collection.put(identifier);
178  }
179 
180  return in;
181  }
182 
183 
184  /**
185  * Stream output.
186  *
187  * \param out output stream
188  * \param collection collection of neutrino PDG types
189  * \return output stream
190  */
191  friend inline std::ostream& operator<<(std::ostream& out, const JNeutrinoTypeCollection& collection)
192  {
193  using namespace std;
194 
195  const vector<int>& identifiers = collection.identifiers;
196 
197  for (vector<int>::const_iterator i = identifiers.cbegin(); i != identifiers.cend(); ++i) {
198  out << ' ' << *i;
199  }
200 
201  return out;
202  }
203 
204 
205  private:
206 
207  std::vector<int> identifiers; //!< Container for identifiers
208  };
209 
210 
211  /**
212  * Auxiliary class for parsing multiparticle fluxes.
213  */
214  struct JFluxMap :
215  public JOscProbHelper
216  {
217  /**
218  * Default constructor.
219  */
221  {}
222 
223 
224  /**
225  * Get multiparticle flux function.
226  *
227  * \return multiparticle flux function
228  */
230  {
231  using namespace std;
232  using namespace JPP;
233 
234  JFluxMultiParticle multiFlux;
235 
236  for (map<int, JFlatFlux>::const_iterator i = flatFluxes.cbegin(); i != flatFluxes.cend(); ++i) {
237  multiFlux.insert(i->first, make_fluxFunction(i->second));
238  }
239 
240  for (map<int, JPowerLawFlux>::const_iterator i = powerLawFluxes.cbegin(); i != powerLawFluxes.cend(); ++i) {
241  multiFlux.insert(i->first, make_fluxFunction(i->second));
242  }
243 
244  const JOscProbInterface& oscProbInterface = getOscProbInterface();
245 
246  const JAtmosphericNeutrinoFlux atmFlux(oscProbInterface);
247 
249 
250  for (vector<int>::const_iterator i = atmosphericFluxes.cbegin(); i != atmosphericFluxes.cend(); ++i) {
251  multiFlux.insert(*i, atmFluxFunction);
252  }
253 
254  return multiFlux;
255  }
256 
257 
258  /**
259  * Conversion operator.
260  *
261  * \return multiparticle flux function
262  */
263  operator JFluxMultiParticle() const
264  {
265  return getMultiParticleFlux();
266  }
267 
268 
269  /**
270  * Stream input.
271  *
272  * \param in input stream
273  * \param fluxMap flux map
274  * \return input stream
275  */
276  friend inline std::istream& operator>>(std::istream& in, JFluxMap& fluxMap)
277  {
278  return getProperties(fluxMap).read(in);
279  }
280 
281 
282  /**
283  * Stream output.
284  *
285  * \param out output stream
286  * \param fluxMap flux map
287  * \return output stream
288  */
289  friend inline std::ostream& operator<<(std::ostream& out, const JFluxMap& fluxMap)
290  {
291  return getProperties(fluxMap).write(out);
292  }
293 
294 
295  std::map<int, JFlatFlux> flatFluxes; //!< Uniform flux functions
296  std::map<int, JPowerLawFlux> powerLawFluxes; //!< Power-law flux functions
297 
298  JNeutrinoTypeCollection atmosphericFluxes; //!< Atmospheric neutrino flux functions
299 
300  private:
301  /**
302  * Get properties of this class.
303  *
304  * \param object flux map object
305  * \return properties
306  */
307  template<class JFluxMap_t>
308  static inline JProperties getProperties(JFluxMap_t& object)
309  {
310  JProperties properties;
311 
312  properties.insert(gmake_property(object.flatFluxes));
313  properties.insert(gmake_property(object.powerLawFluxes));
314  properties.insert(gmake_property(object.atmosphericFluxes));
315 
316  return properties;
317  }
318  };
319 
320 
321  /**
322  * Get volume of given event according given weighter.
323  *
324  * \param type type
325  * \param evt event
326  * \return volume [m^3 Gev sr]
327  */
328  inline double getVolume(const JType<JEvtWeightGSeaGen>& type, const Evt& evt)
329  {
333  }
334 
335 
336  /**
337  * Get volume of given event according given weighter.
338  *
339  * \param type type
340  * \param evt event
341  * \return volume [m^3 Gev sr]
342  */
343  inline double getVolume(const JType<JEvtWeightKM3BUU>& type, const Evt& evt)
344  {
348  }
349 
350 
351  /**
352  * Get volume of given event according given weighter.
353  *
354  * \param type type
355  * \param evt event
356  * \return volume [m^3 Gev sr]
357  */
358  inline double getVolume(const JType<JEvtWeightGenhen>& type, const Evt& evt)
359  {
360  using namespace JPHYSICS;
361 
362  const double l_int = 1.0e-6 * NUCLEON_MOLAR_MASS / (evt.w2list[W2LIST_GENHEN_SIG] * DENSITY_SEA_WATER * AVOGADRO);
363  const double year = 60*60*24*365; // [s]
364 
366  l_int /
367  year /
369  }
370 
371 
372  /**
373  * Get volume of given event according given weighter.
374  *
375  * The return value should be normalised to
376  * - number of generated events;
377  * - solid angle covered by the generation; and
378  * - bin width of the histogram.
379  *
380  * \param weighter weighter
381  * \param evt event
382  * \return volume [m^3 GeV sr]
383  */
384  inline double getVolume(const JEvtWeight& weighter, const Evt& evt)
385  {
386  if (dynamic_cast<const JEvtWeightGSeaGen*>(&weighter) != NULL) { return getVolume(JType<JEvtWeightGSeaGen>(), evt); }
387  if (dynamic_cast<const JEvtWeightGenhen*> (&weighter) != NULL) { return getVolume(JType<JEvtWeightGenhen> (), evt); }
388  if (dynamic_cast<const JEvtWeightKM3BUU*> (&weighter) != NULL) { return getVolume(JType<JEvtWeightKM3BUU> (), evt); }
389 
390  THROW(JCastException, "No valid weighter.");
391  }
392 
393  extern JEvtWeighter getEventWeighter; //!< Function object for mapping header to event weighter.
394 }
395 
396 #endif
friend std::istream & operator>>(std::istream &in, JFluxMap &fluxMap)
Stream input.
static JProperties getProperties(JFluxMap_t &object)
Get properties of this class.
JEvtWeighter getEventWeighter
Function object for mapping header to event weighter.
JNeutrinoTypeCollection()
Default constructor.
bool read(const JEquation &equation)
Read equation.
Definition: JProperties.hh:679
Exceptions.
void insert(const int type, const JEvtWeightFactor_t &factor)
Insert pair of particle code and event-weight factor.
Implementation of event weighing for DAQ data.
JFluxMap()
Default constructor.
friend std::istream & operator>>(std::istream &in, JNeutrinoTypeCollection &collection)
Stream input.
std::vector< int >::const_iterator cbegin() const noexcept
Get constant iterator to the first element of the collection.
JFluxMultiParticle getMultiParticleFlux() const
Get multiparticle flux function.
JProperties & getProperties(T &object, const JEquationParameters &parameters=JEquationParameters(), const int debug=1)
Get properties of a given object.
std::vector< int > identifiers
Container for identifiers.
double getVolume(const JType< JEvtWeightGSeaGen > &type, const Evt &evt)
Get volume of given event according given weighter.
Implementation of event weighting for Genhen data.
friend std::ostream & operator<<(std::ostream &out, const JNeutrinoTypeCollection &collection)
Stream output.
#define gmake_property(A)
macros to convert (template) parameter to JPropertiesElement object
static const int W2LIST_KM3BUU_WATER_INT_LEN
Interaction length in pure water in m.
static const int W2LIST_GENHEN_P_EARTH
Transmission probability in the Earth.
Low-level interface for oscillation probability calculators.
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:712
JEvtWeightFactorMultiParticle< JFlux > JFluxMultiParticle
Type-definition of multi-particle event-weight factor for fluxes.
Utility class to parse parameter values.
Definition: JProperties.hh:497
Implementation of event-weight factor interface.
static const double AVOGADRO
Avogadro&#39;s number.
Helper class for oscillation probabilities.
std::vector< double > w
MC: Weights w[0]=w1, w[1]=w2, w[2]=w3 (see e.g. Tag list or km3net-dataformat/definitions) ...
Definition: Evt.hh:42
Implementation of event weighing for miscellaneous data such as a merged offline file containing neut...
static const double DENSITY_SEA_WATER
Fixed environment values.
Auxiliary class for parsing a vector of neutrino PDG identifiers.
Auxiliary class for a type holder.
Definition: JType.hh:19
JEvtWeightFactorFunction< JFunction_t, JFlux > make_fluxFunction(const JFunction_t &flux)
Auxiliary method for creating an interface to a flux function.
std::map< int, JPowerLawFlux > powerLawFluxes
Power-law flux functions.
std::vector< int >::const_iterator cend() const noexcept
Get constant iterator to the last element of the collection.
Utility class to parse parameter values.
std::ostream & write(std::ostream &out) const
Write the current parameter values.
Definition: JProperties.hh:847
Abstract base class for event weighing.
Definition: JEvtWeight.hh:28
Monte Carlo run header.
Definition: JHead.hh:1234
Implementation of atmospheric neutrino flux using official KM3NeT atmospheric flux function...
Implementation of event-weight factor for multiple particle types.
friend std::ostream & operator<<(std::ostream &out, const JFluxMap &fluxMap)
Stream output.
Implementation of event weighting for GSeaGen data.
size_t size() const
Get size of this collection.
Definition of hit and track types and auxiliary methods for handling Monte Carlo data.
Physics constants.
static const int W2LIST_KM3BUU_P_EARTH
Transmission probability in the Earth.
Auxiliary class for parsing multiparticle fluxes.
Exception for cast operation.
Definition: JException.hh:250
JNeutrinoTypeCollection atmosphericFluxes
Atmospheric neutrino flux functions.
Implementation of event weighting for KM3BUU data.
static const int W2LIST_GSEAGEN_WATER_INT_LEN
Interaction length in pure water in m.
std::map< int, JFlatFlux > flatFluxes
Uniform flux functions.
Implementation of event weighing for MUPAGE data.
static const int W2LIST_GSEAGEN_P_EARTH
Transmission probability in the Earth (XSEC_MEAN and COLUMN_DEPTH used to compute PEarth) ...
Look-up table for event weighters.
static const double NUCLEON_MOLAR_MASS
nucleon molar mass [g/mol]
then fatal The output file must have the wildcard in the e g root fi eval JPrintDetector a $DETECTOR O IDENTIFIER eval JPrintDetector a $DETECTOR O SUMMARY JAcoustics sh $DETECTOR_ID source JAcousticsToolkit sh CHECK_EXIT_CODE typeset A EMITTERS get_tripods $WORKDIR tripod txt EMITTERS get_transmitters $WORKDIR transmitter txt EMITTERS for EMITTER in
Definition: JCanberra.sh:48
JEvtWeighter()
Constructor.
static const int WEIGHTLIST_DIFFERENTIAL_EVENT_RATE
Event rate per unit of flux (c.f. taglist document) [GeV m2 sr].
Definition: weightlist.hh:14
Definition of particle types.
Implementation of event weighting for Corsika data.
then set_variable DETECTOR set_variable OUTPUT_FILE set_variable DAQ_FILE set_variable PMT_FILE else fatal Wrong number of arguments fi JPrintTree f $DAQ_FILE type
Exception for accessing a value in a collection that is outside of its range.
Definition: JException.hh:178
static const int W2LIST_GENHEN_SIG
Cross section of the neutrion interaction.
void put(const int identifier)
Add identifier.
std::vector< double > w2list
MC: factors that make up w[1]=w2 (see e.g. Tag list or km3net-dataformat/definitions) ...
Definition: Evt.hh:43
const JEvtWeight & operator()(const JHead &header) const
Get event weighter corresponding to given header.
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition: Evt.hh:20