Jpp - the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JMultiParticleFlux.cc
Go to the documentation of this file.
1 #include <iostream>
2 #include <iomanip>
3 #include <string>
4 #include <set>
5 
8 
9 #include "JLang/JException.hh"
10 
11 #include "JAAnet/JAAnetToolkit.hh"
12 #include "JAAnet/JFluxFunction.hh"
14 
16 #include "JSupport/JLimit.hh"
18 
19 #include "Jeep/JPrint.hh"
20 #include "Jeep/JParser.hh"
21 #include "Jeep/JMessage.hh"
22 #include "Jeep/JProperties.hh"
23 
24 
25 namespace {
26 
29 
30 
31  /**
32  * Set flux of given event to zero.
33  *
34  * \param evt event
35  * \return zero [GeV * m^-2 * sr^-1 * s^-1]
36  */
37  inline double zeroFlux(const Evt& evt)
38  {
39  return 0.0;
40  }
41 
42 
43  /**
44  * Example function object for constant flux.
45  */
46  struct JFlatFlux
47  {
48  /**
49  * Default constructor.
50  */
51  JFlatFlux() :
52  flux(0.0)
53  {}
54 
55 
56  /**
57  * Constructor.
58  *
59  * \param flux [GeV * m^-2 * sr^-1 * s^-1]
60  */
61  JFlatFlux(const double flux) :
62  flux(flux)
63  {}
64 
65 
66  /**
67  * Get flux of given event.
68  *
69  * \param evt event
70  * \return flux [GeV * m^-2 * sr^-1 * s^-1]
71  */
72  double operator()(const Evt& evt) const
73  {
74  return flux;
75  }
76 
77 
78  /**
79  * Stream input.
80  *
81  * \param in input stream
82  * \param object uniform flux
83  * \return input stream
84  */
85  friend inline std::istream& operator>>(std::istream& in, JFlatFlux& object)
86  {
87  in >> object.flux;
88 
89  return in;
90  }
91 
92 
93  /**
94  * Write power-law neutrino fluxes to output stream.
95  *
96  * \param out output stream
97  * \param object uniform
98  * \return output stream
99  */
100  friend inline std::ostream& operator<<(std::ostream& out, const JFlatFlux& object)
101  {
102  const JFormat format(out);
103 
104  out << FIXED(5,3) << object.flux;
105 
106  return out;
107  }
108 
109  double flux; //!< flux [GeV * m^-2 * sr^-1 * s^-1]
110  };
111 
112 
113  /**
114  * Example function object for computing power-law flux.
115  */
116  struct JPowerLawFlux
117  {
118  /**
119  * Default constructor.
120  */
121  JPowerLawFlux() :
122  normalisation(1.0),
123  spectralIndex(0.0)
124  {}
125 
126 
127  /**
128  * Constructor.
129  *
130  * \param normalisation normalisation
131  * \param spectralIndex spectral index
132  */
133  JPowerLawFlux(const double normalisation,
134  const double spectralIndex) :
135  normalisation(normalisation),
136  spectralIndex(spectralIndex)
137  {}
138 
139 
140  /**
141  * Get flux of given event.
142  *
143  * \param evt event
144  * \return flux [GeV * m^-2 * sr^-1 * s^-1]
145  */
146  double operator()(const Evt& evt) const
147  {
148  using namespace std;
149  using namespace JPP;
150 
151  if (has_neutrino(evt)) {
152 
153  const Trk& neutrino = get_neutrino(evt);
154 
155  return normalisation * pow(neutrino.E, -spectralIndex);
156 
157  } else {
158 
159  THROW(JNullPointerException, "JPowerLawFlux::getFlux(): No neutrino for event " << evt.id << '.' << endl);
160  }
161  }
162 
163 
164  /**
165  * Stream input.
166  *
167  * \param in input stream
168  * \param object power-law flux
169  * \return input stream
170  */
171  inline friend std::istream& operator>>(std::istream& in,
172  JPowerLawFlux& object)
173  {
174  return in >> object.normalisation
175  >> object.spectralIndex;
176  }
177 
178 
179  /**
180  * Write power-law parameters to output stream.
181  *
182  * \param out output stream
183  * \param object power-law flux
184  * \return output stream
185  */
186  inline friend std::ostream& operator<<(std::ostream& out,
187  const JPowerLawFlux& object)
188  {
189  const JFormat format(out);
190 
191  out << FIXED(5,3) << object.normalisation << ' '
192  << FIXED(5,3) << object.spectralIndex;
193 
194  return out;
195  }
196 
197 
198  double normalisation; //!< normalisation [GeV * m^-2 * sr^-1 * s^-1]
199  double spectralIndex; //!< spectral index <= 0
200  };
201 }
202 
203 
204 /**
205  * \file
206  * Example program for scanning event-weights of Monte Carlo files\n
207  * containing a given set of primaries one primary.
208  *
209  * The list of possible options for the flux includes:
210  * <pre>
211  * -@ zero <type>
212  * -@ flat <type> <value>
213  * -@ powerlaw <type> <normalisation> <spectral index>
214  * </pre>
215  *
216  * \author bjung
217  */
218 int main(int argc, char **argv)
219 {
220  using namespace std;
221  using namespace JPP;
222 
223  JMultipleFileScanner_t inputFiles;
224  JLimit numberOfEvents;
225 
226  vector<int> zeroFluxes;
227  map<int, JFlatFlux> flatFluxes;
228  map<int, JPowerLawFlux> powerlawFluxes;
229 
230  int debug;
231 
232  try {
233 
234  JProperties fluxMaps;
235 
236  fluxMaps["zero"] = zeroFluxes;
237  fluxMaps["flat"] = flatFluxes;
238  fluxMaps["powerlaw"] = powerlawFluxes;
239 
240  JParser<> zap;
241 
242  zap['f'] = make_field(inputFiles);
243  zap['@'] = make_field(fluxMaps) = JPARSER::initialised();
244  zap['n'] = make_field(numberOfEvents) = JLimit::max();
245  zap['d'] = make_field(debug) = 1;
246 
247  zap(argc, argv);
248  }
249  catch(const exception& error) {
250  FATAL(error.what() << endl);
251  }
252 
253  // Sanity checks
254 
255  for (map<int, JPowerLawFlux>::const_iterator i = powerlawFluxes.cbegin(); i != powerlawFluxes.cend(); ++i) {
256  if (abs(i->first) != TRACK_TYPE_NUE &&
257  abs(i->first) != TRACK_TYPE_NUMU &&
258  abs(i->first) != TRACK_TYPE_NUTAU) {
259  FATAL("Particle type " << i->first << " does not correspond to a neutrino." << endl);
260  }
261  }
262 
263  // Create multi-particle flux function
264 
265  JMultiParticleFlux multiFlux;
266 
267  for (vector<int>::const_iterator i = zeroFluxes.cbegin(); i != zeroFluxes.cend(); ++i) {
268  multiFlux.insert(*i, make_fluxFunction(zeroFlux));
269  }
270 
271  for (map<int, JFlatFlux>::const_iterator i = flatFluxes.cbegin(); i != flatFluxes.cend(); ++i) {
272  multiFlux.insert(i->first, make_fluxFunction(i->second));
273  }
274 
275  for (map<int, JPowerLawFlux>::const_iterator i = powerlawFluxes.cbegin(); i != powerlawFluxes.cend(); ++i) {
276  multiFlux.insert(i->first, make_fluxFunction(i->second));
277  }
278 
279 
280  // Set event weighter
281 
282  JWeightFileScannerSet<> scanners(inputFiles, numberOfEvents);
283 
284  size_t n = scanners.setFlux(multiFlux);
285 
286  if (n == 0) {
287  WARNING("No file found containing all given primaries; Flux function not set." << endl);
288  }
289 
290 
291  // Scan events
292 
293  for (JWeightFileScannerSet<>::iterator scanner = scanners.begin(); scanner != scanners.end(); ++scanner) {
294 
295  if (scanner->simul.size() > 0) {
296  STATUS("Scanning " << scanner->simul[0].program << " files..." << endl);
297  }
298 
299  STATUS(LEFT(15) << "event" << RIGHT(15) << "weight" << endl);
300 
301  while (scanner->hasNext()) {
302 
303  const Evt* event = scanner->next();
304  const double weight = scanner->getWeight(*event);
305 
306  STATUS(LEFT (15) << scanner->getCounter() <<
307  SCIENTIFIC(15, 3) << weight << '\r'); DEBUG(endl);
308  }
309  }
310 
311  return 0;
312 }
Utility class to parse command line options.
Definition: JParser.hh:1500
Implementation of flux function for multiple particle types.
#define WARNING(A)
Definition: JMessage.hh:65
Exceptions.
int main(int argc, char *argv[])
Definition: Main.cc:15
#define STATUS(A)
Definition: JMessage.hh:63
bool has_neutrino(const Evt &evt)
Test whether given event has an incoming neutrino.
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:670
Utility class to parse parameter values.
Definition: JProperties.hh:496
void insert(const int type, const JFlux &flux)
Insert pair of particle code and flux function.
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:66
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:446
double E
Energy [GeV] (either MC truth or reconstructed)
Definition: Trk.hh:19
Auxiliary class for organising Monte Carlo file scanners.
Utility class to parse parameter values.
std::vector< value_type >::iterator iterator
Auxiliary class to temporarily define format specifications.
Definition: JManip.hh:632
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
I/O formatting auxiliaries.
Exception for null pointer operation.
Definition: JException.hh:216
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1961
Neutrino flux.
Definition: JHead.hh:839
size_t setFlux(const int type, const JFlux &function)
Set flux function for all MC-files corresponding to a given PDG code.
Definition of hit and track types and auxiliary methods for handling Monte Carlo data.
T pow(const T &x, const double y)
Power .
Definition: JMath.hh:98
int debug
debug level
Definition: JSirene.cc:63
General purpose messaging.
#define FATAL(A)
Definition: JMessage.hh:67
Scanning of objects from multiple files according a format that follows from the extension of each fi...
JFluxFunction< JFunction_t > make_fluxFunction(const JFunction_t &function)
Auxiliary method for creating flux function.
Auxiliary base class for list of file names.
std::istream & operator>>(std::istream &in, JAANET::JHead &header)
Read header from input.
Definition: JHead.hh:1618
Utility class to parse command line options.
alias put_queue eval echo n
Definition: qlib.csh:19
std::ostream & operator<<(std::ostream &stream, const CLBCommonHeader &header)
int id
offline event identifier
Definition: Evt.hh:21
Exception for accessing a value in a collection that is outside of its range.
Definition: JException.hh:162
const Trk & get_neutrino(const Evt &evt)
Get incoming neutrino.
Auxiliaries for defining the range of iterations of objects.
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable INPUT_FILE $argv[2] eval JPrintDetector a $DETECTOR O IDENTIFIER eval JPrintDetector a $DETECTOR O SUMMARY source JAcoustics sh $DETECTOR_ID CHECK_EXIT_CODE typeset A TRIPODS get_tripods $WORKDIR tripod txt TRIPODS for EMITTER in
Definition: JCanberra.sh:38
The Trk class represents a Monte Carlo (MC) particle as well as a reconstructed track/shower.
Definition: Trk.hh:13
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:484
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition: Evt.hh:19
std::vector< double > weight
Definition: JAlgorithm.hh:428
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62