Jpp  18.6.0-rc.1
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JEvtWeightFactorGSeaGen.hh
Go to the documentation of this file.
1 #ifndef __JAANET__JEVTWEIGHTFACTORGSEAGEN__
2 #define __JAANET__JEVTWEIGHTFACTORGSEAGEN__
3 
5 
8 
9 #include "JLang/JManip.hh"
10 #include "JLang/JToken.hh"
11 #include "JLang/JException.hh"
12 
13 #include "JAAnet/JFlux.hh"
15 
16 #include "JGizmo/JGizmoToolkit.hh"
17 
18 #include "TFormula.h"
19 
20 
21 /**
22  * \author bjung
23  */
24 
25 namespace JAANET {}
26 namespace JPP { using namespace JAANET; }
27 
28 namespace JAANET {
29 
30  /**
31  * Implementation of reweighting factor for simulated neutrino interactions according to a specifiable ROOT TFormula.
32  *
33  * Note: The ROOT TFormula may assume any number of parameters, but should be restricted to\n
34  * the physical variables listed among `JEvtWeightFactorGSeaGen::variables`.\n
35  * These variables may be specified within the ROOT TFormula as 'x[<index>]',\n
36  * where <index> corresponds to the index of the desired physical variable within `JEvtWeightFactorGSeaGen::variables`.
37  */
39  public JFluxHelper,
40  public TFormula
41  {
42  /**
43  * Indices of reweighting variables for GSeaGen.
44  */
45  enum variables {
46  COSTH, //!< Cosine zenith angle
47  ENERGY, //!< Initial state energy
48  BJORKEN_X, //!< Björken-x (= fractional momentum carried by the struck nucleon)
49  INELASTICITY, //!< Inelasticity (= Björken-y)
50  INTERACTION_TYPE, //!< Interaction channel type
51  CURRENT_TYPE, //!< Weak current type (CC or NC)
52  XSEC, //!< Exclusive total cross-section of the interaction
53  XSEC_MEAN, //!< Average interaction cross-section per nucleon along neutrino path [m2]
54  XSEC_DIFFERENTIAL, //!< Differential cross-section of the interaction
55  XSEC_WATER, //!< Inclusive cross-section in water
56  INT_LENGTH_WATER, //!< Interaction length in pure water
57  COLUMN_DEPTH, //!< Column density [m.w.e]
58  P_EARTH, //!< Earth transmission probability
59  P_SCALE, //!< GENIE ineraction probability scale
60  TARGET_A, //!< Number of nucleons in the target
61  TARGET_Z, //!< Number of protons in the target
62  NUMBER_OF_VARIABLES //!< Number of reweighting variables; \n
63  //!< N.B.\ This enum value needs to be specified last!
64  };
65 
66 
67  /**
68  * Default constructor.
69  */
71  JFluxHelper(),
72  TFormula()
73  {}
74 
75 
76  /**
77  * Constructor.
78  *
79  * \param flux flux function
80  * \param name name
81  * \param formula formula
82  */
84  const char* name,
85  const char* formula) :
86  JFluxHelper(flux),
87  TFormula(name, formula)
88  {}
89 
90 
91  /**
92  * Get weighting factor for given event.
93  *
94  * \param evt event
95  * \return weighting factor
96  */
97  double operator()(const Evt& evt) const
98  {
99  using namespace std;
100  using namespace JPP;
101 
102  if (static_cast<const JFluxHelper&>(*this)) {
103 
104  Double_t vars[NUMBER_OF_VARIABLES] = { 0.0 };
105 
106  const Trk& neutrino = get_neutrino(evt);
107  const double flux = getFactor(evt);
108 
109  vars[COSTH] = neutrino.dir.z;
110  vars[ENERGY] = getE0(evt);
111  vars[BJORKEN_X] = evt.w2list[W2LIST_GSEAGEN_BX];
115  vars[XSEC] = evt.w2list[W2LIST_GSEAGEN_XSEC];
125 
126  return this->DoEval(&vars[0]) * flux;
127 
128  } else {
129 
130  THROW(JNullPointerException, "JEvtWeightFactorGSeaGen::operator(): Unspecified flux function.");
131  }
132  }
133 
134 
135  /**
136  * Stream input.
137  *
138  * \param in input stream
139  * \param object gSeaGen event-weight factor
140  * \return input stream
141  */
142  inline friend std::istream& operator>>(std::istream& in,
143  JEvtWeightFactorGSeaGen& object)
144  {
145  using namespace std;
146  using namespace JPP;
147 
148  typedef JToken<';'> JToken_t;
149 
150  JToken_t formula;
151  JToken_t equation;
152 
153  in >> formula;
154 
155  object.Clear();
156  object.Compile(formula.c_str());
157 
158  for (int count = 0; count < object.GetNpar() && in >> equation; ++count) {
159 
160  const int index = getParameter(equation);
161  const double value = getValue(equation, 0);
162 
163  object.SetParameter(index, value);
164  }
165 
166  return in;
167  }
168 
169 
170  /**
171  * Stream output.
172  *
173  * \param out output stream
174  * \param object gSeaGen event-weight factor
175  * \return output stream
176  */
177  inline friend std::ostream& operator<<(std::ostream& out,
178  const JEvtWeightFactorGSeaGen& object)
179  {
180  using namespace std;
181  using namespace JPP;
182 
183  out << object.GetExpFormula() << ';';
184 
185  for (int i = 0; i < object.GetNpar(); ++i) {
186  out << " p" << i << " = " << SCIENTIFIC(6,3) << object.GetParameter(i) << ';';
187  }
188 
189  return out;
190  }
191  };
192 }
193 
194 #endif
static const int W2LIST_GSEAGEN_BX
Bjorken x.
static const int W2LIST_GSEAGEN_P_SCALE
Interaction probability scale.
GENIE ineraction probability scale.
Exceptions.
double getValue(const JScale_t scale)
Get numerical value corresponding to scale.
Definition: JScale.hh:47
int getParameter(const std::string &text)
Get parameter number from text string.
double z
Definition: Vec.hh:14
double getE0(const Evt &evt)
Get initial state energy of a neutrino interaction.
static const int W2LIST_GSEAGEN_XSEC
exclusive total cross section of the interaction
friend std::ostream & operator<<(std::ostream &out, const JEvtWeightFactorGSeaGen &object)
Stream output.
Vec dir
track direction
Definition: Trk.hh:18
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:712
then usage $script< input file >[option] nPossible options count
Definition: JVolume1D.sh:31
variables
Indices of reweighting variables for GSeaGen.
friend std::istream & operator>>(std::istream &in, JEvtWeightFactorGSeaGen &object)
Stream input.
static const int W2LIST_GSEAGEN_TARGETZ
number of protons in the target
Exclusive total cross-section of the interaction.
Average interaction cross-section per nucleon along neutrino path [m2].
Number of reweighting variables; N.B. This enum value needs to be specified last! ...
static const int W2LIST_GSEAGEN_TARGETA
number of nuclons in the target
Differential cross-section of the interaction.
JEvtWeightFactorGSeaGen(const JFlux &flux, const char *name, const char *formula)
Constructor.
double getFactor(const Evt &evt) const
Get weight-factor of given event.
static const int W2LIST_GSEAGEN_ICHAN
Interaction channel.
Neutrino flux.
Definition: JHead.hh:906
Helper class for event-weight factor.
Björken-x (= fractional momentum carried by the struck nucleon)
Low-level interface for retrieving the flux corresponding to a given event.
Definition: JFlux.hh:21
then fatal The output file must have the wildcard in the name
Definition: JCanberra.sh:31
I/O manipulators.
JEvtWeightFactorGSeaGen()
Default constructor.
Implementation of reweighting factor for simulated neutrino interactions according to a specifiable R...
static const int W2LIST_GSEAGEN_WATER_INT_LEN
Interaction length in pure water in m.
static const int W2LIST_GSEAGEN_CC
Charged current interaction flag.
static const int W2LIST_GSEAGEN_BY
Bjorken y.
static const int W2LIST_GSEAGEN_P_EARTH
Transmission probability in the Earth (XSEC_MEAN and COLUMN_DEPTH used to compute PEarth) ...
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
static const int W2LIST_GSEAGEN_WATERXSEC
inclusive xsec in water
static const int W2LIST_GSEAGEN_COLUMN_DEPTH
Line integrated column density through the Earth for the neutrino direction.
static const int W2LIST_GSEAGEN_DXSEC
differential cross section of the interaction (dsigma/dxdy) extracted from genie
double operator()(const Evt &evt) const
Get weighting factor for given event.
static const int W2LIST_GSEAGEN_XSEC_MEAN
Average interaction cross-section per nucleon along the neutrino path throuh the Earth (in units of m...
const Trk & get_neutrino(const Evt &evt)
Get incoming neutrino.
The Trk class represents a Monte Carlo (MC) particle as well as a reconstructed track/shower.
Definition: Trk.hh:14
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:486
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
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition: Evt.hh:20