Jpp  19.1.0
the software that should make you happy
JNeutrinoInteractionCategory.hh
Go to the documentation of this file.
1 #ifndef __JAANET__JNEUTRINOINTERACTIONCATEGORY__
2 #define __JAANET__JNEUTRINOINTERACTIONCATEGORY__
3 
4 #include <utility>
5 #include <iostream>
6 
8 
11 
12 #include "JLang/JClonable.hh"
13 #include "JLang/JParameter.hh"
14 #include "JLang/JPredicate.hh"
15 
16 #include "JTools/JRange.hh"
17 
18 #include "JAAnet/JHead.hh"
19 #include "JAAnet/JHeadToolkit.hh"
20 #include "JAAnet/JParticleTypes.hh"
21 #include "JAAnet/JGENIETypes.hh"
22 #include "JAAnet/JAAnetToolkit.hh"
23 #include "JAAnet/JEvtCategory.hh"
24 
25 
26 /**
27  * \file
28  *
29  * Classes and methods for defining neutrino interaction categories.
30  * \author bjung
31  */
32 namespace JAANET {}
33 namespace JPP { using namespace JAANET; }
34 
35 namespace JAANET {
36 
37  using JLANG::JClonable;
38  using JLANG::JParameter;
39 
40  using JTOOLS::JRange;
41 
42 
43  static const char* const NEUTRAL_CURRENT = "NC"; //!< Neutral current header string identifier
44  static const char* const CHARGED_CURRENT = "CC"; //!< Charged current header string identifier
45 
46 
47  static const JRange<double> DEFAULT_NEUTRINO_ENERGY_RANGE = JRange<double>(0, 1e11); //!< Default neutrino energy range [GeV]
48  static const JRange<double> DEFAULT_NEUTRINO_COSINE_ZENITH_ANGLE_RANGE = JRange<double>(-1.0, 1.0); //!< Default neutrino cosine zenith angle range [-]
49 
50 
51  /**
52  * Auxiliary function to check if given PDG code corresponds to a neutrino.
53  *
54  * \param type PDG type
55  * \return true if given PDG type corresponds to a neutrino; else false
56  */
57  inline bool is_neutrino_primary(const int type)
58  {
59  const int absType = abs(type);
60 
61  return (absType == TRACK_TYPE_NUE ||
62  absType == TRACK_TYPE_NUMU ||
63  absType == TRACK_TYPE_NUTAU);
64  }
65 
66 
67  /**
68  * Auxiliary function to retrieve the GENIE interaction type corresponding to a given header.
69  *
70  * \param header MC header
71  * \return GENIE interaction type
72  */
74  {
75  using namespace std;
76  using namespace JPP;
77 
78  if (is_gseagen(header)) {
79 
80  for (vector<JAANET::physics>::const_iterator i = header.physics.cbegin(); i != header.physics.cend(); ++i) {
81 
82  istringstream iss(i->buffer);
83 
84  for (string buffer; iss >> buffer ;) {
85 
86  if (buffer == string(NEUTRAL_CURRENT)) {
87  return JInteractionTypeGENIE_t::WEAK_NEUTRAL_CURRENT;
88  } else if (buffer == string(CHARGED_CURRENT)) {
89  return JInteractionTypeGENIE_t::WEAK_CHARGED_CURRENT;
90  }
91  }
92  }
93  }
94 
96  }
97 
98 
99  /**
100  * Auxiliary function to retrieve the target atomic number and mass number.
101  *
102  * \param header MC header
103  * \return pair of atomic number (Z) and mass number (A)
104  */
107  {
108  using namespace std;
109  using namespace JPP;
110 
111  JParameter<int> Z;
112  JParameter<int> A;
113 
114  if (!header.target.buffer.empty()) {
115 
116  istringstream iss(header.target.buffer);
117 
118  if (! (iss.peek() == (int) 'A' &&
119  iss.ignore() && iss >> A &&
120  iss.peek() == (int) 'Z' &&
121  iss.ignore() && iss >> Z &&
122  Z > 0 && A >= Z) ) {
123  THROW(JValueOutOfRange, "getTargetInfo(): Invalid target header-field: " << header.target.buffer << endl);
124  }
125  }
126 
127  return make_pair(Z,A);
128  }
129 
130 
131  /**
132  * Class for neutrino interaction categories.
133  */
135  public JClonable<JEvtCategory, JNeutrinoInteractionCategory>
136  {
138 
139 
140  /**
141  * Default constructor.
142  */
144  interactionType(),
145  scatteringType (),
146  targetZ (),
147  targetA (),
150  {}
151 
152 
153  /**
154  * Constructor.
155  *
156  * \param neutrinoType neutrino PDG type
157  */
158  JNeutrinoInteractionCategory(const int neutrinoType) :
159  interactionType(),
160  scatteringType (),
161  targetZ (),
162  targetA (),
165  {
166  this->primaryType = neutrinoType;
167 
168  check_validity();
169  }
170 
171 
172  /**
173  * Constructor.
174  *
175  * \param neutrinoType neutrino PDG type
176  * \param Erange neutrino energy range
177  * \param coszRange zenith-angle range
178  */
179  JNeutrinoInteractionCategory(const int neutrinoType,
180  const JRange_t& Erange,
181  const JRange_t& coszRange) :
182  interactionType(),
183  scatteringType (),
184  targetZ (),
185  targetA (),
186  Erange (Erange),
187  coszRange (coszRange)
188  {
189  this->primaryType = neutrinoType;
190 
191  check_validity();
192  }
193 
194 
195  /**
196  * Constructor.
197  *
198  * \param neutrinoType neutrino PDG type
199  * \param interactionType interaction type
200  * \param scatteringType scattering type
201  * \param targetZ target atomic number
202  * \param targetA target mass number
203  * \param Erange neutrino energy range
204  * \param coszRange zenith-angle range
205  */
206  JNeutrinoInteractionCategory(const int neutrinoType,
207  const JInteractionTypeGENIE_t interactionType,
208  const JScatteringTypeGENIE_t scatteringType,
209  const int targetZ,
210  const int targetA,
211  const JRange_t& Erange,
212  const JRange_t& coszRange) :
213  interactionType((int) interactionType),
214  scatteringType ((int) scatteringType),
215  targetZ (targetZ),
216  targetA (targetA),
217  Erange (Erange),
218  coszRange (coszRange)
219  {
220  this->primaryType = neutrinoType;
221 
222  check_validity();
223  }
224 
225 
226  /**
227  * Constructor.
228  *
229  * \param header MC header
230  * \param primaryType primary PDG type
231  */
233  const int primaryType) :
235  {
236  using namespace std;
237  using namespace JPP;
238 
239  if (is_gseagen(header)) {
240 
241  vector<JAANET::flux>::const_iterator i = find_if(header.flux.cbegin(), header.flux.cend(),
242  make_predicate(&JAANET::flux::type, primaryType));
243 
244  if (primaryType == header.primary.type || i != header.flux.cend()) {
245  this->primaryType = primaryType;
246  } else {
247  THROW(JValueOutOfRange, "JNeutrinoInteractionCategory::JNeutrinoInteractionCategory(): Given header does not match with given primary type " << primaryType << ".");
248  }
249 
250  const JInteractionTypeGENIE_t intType = getInteractionType(header);
251 
252  if (intType != JInteractionTypeGENIE_t::UNDEFINED) {
253  interactionType = (int) intType;
254  }
255 
256  const pair<JParameter<int>,
257  JParameter<int> > targetInfo = getTargetInfo(header);
258 
259  targetZ = targetInfo.first;
260  targetA = targetInfo.second;
261 
262  if (header.cut_nu.E.is_valid() &&
263  header.cut_nu.E != cut().E) {
264  Erange = header.cut_nu.E;
265  }
266 
267  if (header.cut_nu.cosT.is_valid() &&
268  header.cut_nu.cosT != cut().cosT) {
269  coszRange = header.cut_nu.cosT;
270  }
271  }
272 
273  check_validity();
274  }
275 
276 
277  /**
278  * Constructor.
279  *
280  * N.B.: In this method, the neutrino PDG type is taken from `JAANET::JHead::primary::type`.
281  *
282  * \param header MC header
283  */
285  JNeutrinoInteractionCategory(header, header.primary.type)
286  {}
287 
288 
289  /**
290  * Check if neutrino interaction categories is valid.
291  *
292  * \return true if neutrino interaction category is valid; else false
293  */
294  inline bool is_valid() const override final
295  {
296  if (!is_neutrino_primary(this->primaryType)) {
297  return false;
298  }
299 
300  if ( interactionType.isDefined() &&
301  !(interactionType.getValue() > (int) JInteractionTypeGENIE_t::INTERACTION_TYPES_START &&
302  interactionType.getValue() < (int) JInteractionTypeGENIE_t::INTERACTION_TYPES_END) ) {
303  return false;
304  }
305 
306  if ( scatteringType.isDefined() &&
307  !((scatteringType > (int) JScatteringTypeGENIE_t::NEUTRINO_SCATTERING_TYPES_START &&
308  scatteringType < (int) JScatteringTypeGENIE_t::NEUTRINO_SCATTERING_TYPES_END) ||
309  (scatteringType > (int) JScatteringTypeGENIE_t::DARK_MATTER_SCATTERING_TYPES_START &&
310  scatteringType < (int) JScatteringTypeGENIE_t::DARK_MATTER_SCATTERING_TYPES_END)) ) {
311  return false;
312  }
313 
314  if (targetZ.isDefined() && targetZ.getValue() <= 0) {
315  return false;
316  }
317 
318  if ( targetA.isDefined() &&
319  (!targetZ.isDefined() || targetA.getValue() < targetZ.getValue())) {
320  return false;
321  }
322 
323  return Erange.is_valid() && coszRange.is_valid();
324  }
325 
326 
327  /**
328  * Check whether this event category matches with the given MC header.
329  *
330  * \param header header
331  * \return true if matching; else false
332  */
333  bool match(const JHead& header) const override final
334  {
335  if (is_gseagen(header)) {
336 
337  const JNeutrinoInteractionCategory cat(header);
338 
339  return (this->primaryType == cat.primaryType &&
340  this->interactionType == cat.interactionType &&
341  this->scatteringType == cat.scatteringType &&
342  this->targetZ == cat.targetZ &&
343  this->targetA == cat.targetA &&
344  this->Erange .overlap(cat.Erange) &&
345  this->coszRange .overlap(cat.coszRange));
346 
347  } else {
348 
349  return false;
350  }
351  }
352 
353 
354  /**
355  * Check whether given event matches with this event category.
356  *
357  * \param event event
358  * \return true if matching; else false
359  */
360  bool match(const Evt& event) const override final
361  {
362  using namespace std;
363  using namespace JPP;
364 
365  const Trk& primary = get_primary(event);
366 
367  if (!is_neutrino(primary)) { return false; }
368 
369  if (interactionType.isDefined() &&
370  event.w2list.size() > W2LIST_GSEAGEN_CC &&
371  interactionType.getValue() != event.w2list[W2LIST_GSEAGEN_CC]) {
372  return false;
373  }
374 
375  if (scatteringType.isDefined() &&
376  event.w2list.size() > W2LIST_GSEAGEN_ICHAN &&
377  scatteringType.getValue() != event.w2list[W2LIST_GSEAGEN_ICHAN]) {
378  return false;
379  }
380 
381  if (targetZ.isDefined() &&
382  event.w2list.size() > W2LIST_GSEAGEN_TARGETZ &&
383  targetZ.getValue() != event.w2list[W2LIST_GSEAGEN_TARGETZ]) {
384  return false;
385  }
386 
387  if (targetA != 0 &&
388  event.w2list.size() > W2LIST_GSEAGEN_TARGETA &&
389  targetA.getValue() != event.w2list[W2LIST_GSEAGEN_TARGETA]) {
390  return false;
391  }
392 
393  return Erange(primary.E) && coszRange(-primary.dir.z);
394  }
395 
396 
397  /**
398  * Less-than method.
399  *
400  * Note: The energy and cosine zenith-angle ranges are compared based on the lower limit first.\n
401  * If the lower limits are equal, the upper limits are compared.
402  *
403  * \param category event category
404  * \return true if first this event category less than given event category; else false
405  */
406  inline bool less(const JEvtCategory& category) const override final
407  {
408 #define RETURN_IF_DIFFERENT(A, B) \
409  if (A != B) { return A < B; }
410 
411  RETURN_IF_DIFFERENT(this->primaryType, category.getPrimaryType());
412 
413  const JNeutrinoInteractionCategory* p = dynamic_cast<const JNeutrinoInteractionCategory*>(&category);
414 
415  if (p == NULL) { return false; }
416 
417  RETURN_IF_DIFFERENT(this->interactionType, p->interactionType);
418  RETURN_IF_DIFFERENT(this->scatteringType, p->scatteringType);
419  RETURN_IF_DIFFERENT(this->targetZ, p->targetZ);
420  RETURN_IF_DIFFERENT(this->targetA, p->targetA);
421 
422  if (p->Erange != this->Erange) {
423  return (this->Erange.getLowerLimit() == p->Erange.getLowerLimit() ?
424  this->Erange.getUpperLimit() < p->Erange.getUpperLimit() :
425  this->Erange.getLowerLimit() < p->Erange.getLowerLimit());
426  } else {
427  return (this->coszRange.getLowerLimit() == p->coszRange.getLowerLimit() ?
428  this->coszRange.getUpperLimit() < p->coszRange.getUpperLimit() :
429  this->coszRange.getLowerLimit() < p->coszRange.getLowerLimit());
430  }
431 #undef RETURN_IF_DIFFERENT
432  }
433 
434 
435  /**
436  * Get properties of this class.
437  *
438  * \param eqpars equation parameters
439  */
441  {
442  return JNeutrinoInteractionCategoryHelper(*this, eqpars);
443  }
444 
445 
446  /**
447  * Get properties of this class.
448  *
449  * \param eqpars equation parameters
450  */
452  {
453  return JNeutrinoInteractionCategoryHelper(*this, eqpars);
454  }
455 
456 
457  private:
458 
459  JParameter<int> interactionType; //!< interaction type
460  JParameter<int> scatteringType; //!< scattering type
461  JParameter<int> targetZ; //!< target atomic number
462  JParameter<int> targetA; //!< target mass number
463  JRange_t Erange; //!< neutrino energy range
464  JRange_t coszRange; //!< cosine zenith range
465 
466 
467  /**
468  * Auxiliary class for I/O of neutrino interaction category.
469  */
471  public JProperties
472  {
473  /**
474  * Constructor.
475  *
476  * \param cat neutrino interaction category
477  * \param eqpars equation parameters
478  */
479  template<class JNeutrinoInteractionCategory_t>
480  JNeutrinoInteractionCategoryHelper(JNeutrinoInteractionCategory_t& cat,
481  const JEquationParameters& eqpars) :
482  JProperties(eqpars, 1)
483  {
484  (*this)["neutrinoType"] = cat.primaryType;
485 
486  this->insert(gmake_property(cat.interactionType));
487  this->insert(gmake_property(cat.scatteringType));
488  this->insert(gmake_property(cat.targetZ));
489  this->insert(gmake_property(cat.targetA));
490  this->insert(gmake_property(cat.Erange));
491  this->insert(gmake_property(cat.coszRange));
492  }
493  };
494  };
495 }
496 
497 #endif
Definition of hit and track types and auxiliary methods for handling Monte Carlo data.
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:712
Definition of GENIE neutrino interaction types.
#define RETURN_IF_DIFFERENT(A, B)
Definition of particle types.
#define gmake_property(A)
macros to convert (template) parameter to JPropertiesElement object
Auxiliary class to define a range between two values.
Monte Carlo run header.
Definition: JHead.hh:1236
JAANET::cut_nu cut_nu
Definition: JHead.hh:1596
JAANET::primary primary
Definition: JHead.hh:1611
std::vector< JAANET::flux > flux
Definition: JHead.hh:1612
JAANET::target target
Definition: JHead.hh:1589
std::vector< JAANET::physics > physics
Definition: JHead.hh:1590
Utility class to parse parameter values.
Definition: JProperties.hh:501
Simple data structure to support I/O of equations (see class JLANG::JEquation).
Parameter class.
Definition: JParameter.hh:36
Exception for accessing a value in a collection that is outside of its range.
Definition: JException.hh:180
Range of values.
Definition: JRange.hh:42
bool is_valid() const
Check validity of range.
Definition: JRange.hh:311
T getLowerLimit() const
Get lower limit.
Definition: JRange.hh:202
T getUpperLimit() const
Get upper limit.
Definition: JRange.hh:213
Extensions to Evt data format.
JScatteringTypeGENIE_t
Enumeration of GENIE scattering types.
Definition: JGENIETypes.hh:39
const Trk & get_primary(const Evt &evt)
Get primary.
std::pair< JParameter< int >, JParameter< int > > getTargetInfo(const JHead &header)
Auxiliary function to retrieve the target atomic number and mass number.
static const char *const CHARGED_CURRENT
Charged current header string identifier.
JInteractionTypeGENIE_t getInteractionType(const JHead &header)
Auxiliary function to retrieve the GENIE interaction type corresponding to a given header.
bool is_neutrino(const Trk &track)
Test whether given track is a neutrino.
JInteractionTypeGENIE_t
Enumeration of GENIE interaction types.
Definition: JGENIETypes.hh:20
@ UNDEFINED
Unknown interaction type.
static const char *const NEUTRAL_CURRENT
Neutral current header string identifier.
bool is_gseagen(const JHead &header)
Check for generator.
Definition: JHeadToolkit.hh:61
static const JRange< double > DEFAULT_NEUTRINO_ENERGY_RANGE
Default neutrino energy range [GeV].
bool is_neutrino_primary(const int type)
Auxiliary function to check if given PDG code corresponds to a neutrino.
@ TRACK_TYPE_NUTAU
@ TRACK_TYPE_NUE
@ TRACK_TYPE_NUMU
static const JRange< double > DEFAULT_NEUTRINO_COSINE_ZENITH_ANGLE_RANGE
Default neutrino cosine zenith angle range [-].
JPredicate< JResult_t T::*, JComparison::eq > make_predicate(JResult_t T::*member, const JResult_t value)
Helper method to create predicate for data member.
Definition: JPredicate.hh:128
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
Low-level interface for event categories.
Definition: JEvtCategory.hh:38
static JEquationParameters & getEquationParameters()
Get equation parameters.
JNeutrinoInteractionCategoryHelper(JNeutrinoInteractionCategory_t &cat, const JEquationParameters &eqpars)
Constructor.
Class for neutrino interaction categories.
JNeutrinoInteractionCategory(const JHead &header, const int primaryType)
Constructor.
bool is_valid() const override final
Check if neutrino interaction categories is valid.
JProperties getProperties(const JEquationParameters &eqpars=JEvtCategory::getEquationParameters()) override final
Get properties of this class.
bool match(const Evt &event) const override final
Check whether given event matches with this event category.
JParameter< int > interactionType
interaction type
JNeutrinoInteractionCategory(const int neutrinoType, const JRange_t &Erange, const JRange_t &coszRange)
Constructor.
JNeutrinoInteractionCategory(const JHead &header)
Constructor.
JNeutrinoInteractionCategory(const int neutrinoType, const JInteractionTypeGENIE_t interactionType, const JScatteringTypeGENIE_t scatteringType, const int targetZ, const int targetA, const JRange_t &Erange, const JRange_t &coszRange)
Constructor.
bool less(const JEvtCategory &category) const override final
Less-than method.
JNeutrinoInteractionCategory(const int neutrinoType)
Constructor.
JParameter< int > targetZ
target atomic number
JProperties getProperties(const JEquationParameters &eqpars=JEvtCategory::getEquationParameters()) const override final
Get properties of this class.
bool match(const JHead &header) const override final
Check whether this event category matches with the given MC header.
JParameter< int > targetA
target mass number
JParameter< int > scatteringType
scattering type
std::string buffer
General purpose name.
Definition: JHead.hh:217
General purpose class of phase space generation.
Definition: JHead.hh:341
JRange_t cosT
Cosine zenith angle range
Definition: JHead.hh:420
JRange_t E
Energy range [GeV].
Definition: JHead.hh:419
int type
Type.
Definition: JHead.hh:940
Primary particle.
Definition: JHead.hh:1174
int type
Particle type.
Definition: JHead.hh:1204
Template class for object cloning.
Definition: JClonable.hh:59
The Trk class represents a Monte Carlo (MC) particle as well as a reconstructed track/shower.
Definition: Trk.hh:15
static const int W2LIST_GSEAGEN_CC
Charged current interaction flag.
static const int W2LIST_GSEAGEN_TARGETZ
number of protons in the target
static const int W2LIST_GSEAGEN_TARGETA
number of nuclons in the target
static const int W2LIST_GSEAGEN_ICHAN
Interaction channel.