Jpp  18.0.1-rc.1
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions
JEvtReweightGSeaGen.cc File Reference

Program for reweighting gSeaGen data according to a specifiable TFormula. More...

#include <uuid/uuid.h>
#include <iostream>
#include <string>
#include <vector>
#include "km3net-dataformat/definitions/weightlist.hh"
#include "km3net-dataformat/offline/MultiHead.hh"
#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/Evt.hh"
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
#include "Jeep/JProperties.hh"
#include "JLang/JToken.hh"
#include "JAAnet/JHead.hh"
#include "JAAnet/JHeadToolkit.hh"
#include "JAAnet/JAAnetToolkit.hh"
#include "JAAnet/JParticleTypes.hh"
#include "JAAnet/JEvtWeightFactorGSeaGen.hh"
#include "JAAnet/JEvtWeightFactorFunction.hh"
#include "JAAnet/JAtmosphericNeutrinoFlux.hh"
#include "JAAnet/JFlux.hh"
#include "JOscProb/JOscParameters.hh"
#include "JOscProb/JOscProbInterpolator.hh"
#include "JSupport/JMeta.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JFileRecorder.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JEvtWeightFileScanner.hh"
#include "JGizmo/JGizmoToolkit.hh"

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Detailed Description

Program for reweighting gSeaGen data according to a specifiable TFormula.

One TFormula can be specified for each neutrino type. The TFormula may take any number of parameters, but is limited to the physical variables
defined in JEvtWeightFactorGSeaGen. The indices of the physical variables are defined in
the enum JEvtWeightFactorGSeaGen::variables.
Parameter settings can be specified behind the TFormula separated by semicolons, e.g.:
-#anumu="[0]*(x[1]-[1]); p0=1; p1=2.5;"

Author
bjjung

Definition in file JEvtReweightGSeaGen.cc.

Function Documentation

int main ( int  argc,
char **  argv 
)

Definition at line 69 of file JEvtReweightGSeaGen.cc.

70 {
71  using namespace std;
72  using namespace JPP;
73 
74  typedef typename JTYPELIST<JAAnetTypes_t, JMetaTypes_t>::typelist typelist;
76 
77 
78  JMultipleFileScanner_t inputFiles;
80 
81  string oscProbTable;
82  JOscParameters oscParameters;
83 
85 
86  int debug;
87 
88  try {
89 
90  JProperties formulaMap(JEquationParameters("=","\n\r","./","#"));
91 
92  for (map<int, string>::const_iterator i = nuLabels.cbegin(); i != nuLabels.cend(); ++i) {
93  formulaMap[i->second] = formulae[i->first];
94  }
95 
96  JParser<> zap;
97 
98  zap['f'] = make_field(inputFiles);
99  zap['o'] = make_field(outputFile);
100  zap['P'] = make_field(oscProbTable);
101  zap['@'] = make_field(oscParameters) = JPARSER::initialised();
102  zap['#'] = make_field(formulaMap) = JPARSER::initialised();
103  zap['d'] = make_field(debug) = 1;
104 
105  zap(argc, argv);
106  }
107  catch(const exception& error) {
108  FATAL(error.what() << endl);
109  }
110 
111 
112  // Create atmospheric neutrino flux function
113 
114  const JFluxFunction_t fluxFunction = make_atmosphericNeutrinoFluxFunction(oscProbTable.c_str(), oscParameters);
115 
116 
117  // Create file scanner
118 
119  JEvtWeightFileScanner<> scanner(inputFiles);
120 
121  JHead& header = scanner.getHeader();
122  header.createUUID();
123 
124  if (!is_gseagen(header)) {
125  FATAL("MC-header is incompatible with gSeaGen.");
126  }
127 
128 
129  // Define reweighting functions
130 
131  for (map<int, JEvtWeightFactorGSeaGen>::iterator i = formulae.begin(); i != formulae.end(); ++i) {
132 
133  const TString expression = i->second.GetExpFormula();
134 
135  if (!expression.IsNull() && i->second.IsValid() && i->second.GetNdim() <= JEvtWeightFactorGSeaGen::NUMBER_OF_VARIABLES) {
136 
137  i->second.configure(fluxFunction);
138 
139  scanner.setEvtWeightFactor(i->first, make_weightFactor(i->second));
140 
141  } else if (expression.IsNull()) {
142 
143  continue;
144 
145  } else {
146 
147  FATAL("Invalid formula " << expression << " (number of dimensions must be < " <<
148  JEvtWeightFactorGSeaGen::NUMBER_OF_VARIABLES << "; see `JAANET::JEvtWeightFactorGSeaGen`).");
149  }
150  }
151 
152 
153  // Store events with new weights
154 
155  outputFile.open();
156 
157  Head head;
158  copy(header, head);
159  outputFile.put(head);
160 
161  outputFile.put(JMeta(argc, argv));
162 
163  STATUS(RIGHT(15) << "event" << RIGHT(15) << "weight" << endl);
164 
165  while (scanner.hasNext()) {
166 
167  const Evt* event = scanner.next();
168 
169  if (event != NULL) {
170 
171  const double weight = scanner.getWeight(*event);
172 
173  STATUS(RIGHT (15) << scanner.getCounter() <<
174  SCIENTIFIC(15, 3) << weight << '\r'); DEBUG(endl);
175 
176  if (!isfinite(weight)) {
177 
178  WARNING("Non-finite reweighting factor " << weight <<
179  " for event " << scanner.getCounter() << "!");
180  }
181 
182  Evt out = *event;
183 
184  uuid_copy(out.header_uuid, header.UUID.uuid);
185 
186  if (out.w.size() <= WEIGHTLIST_RESCALED_EVENT_RATE) {
187  out.w.resize(WEIGHTLIST_RESCALED_EVENT_RATE + 1, 0.0);
188  }
189 
190  out.w.at(WEIGHTLIST_NORMALISATION) = scanner.getNormalisation(*event);
191  out.w.at(WEIGHTLIST_RESCALED_EVENT_RATE) = weight;
192 
193  outputFile.put(out);
194 
195  } else {
196 
197  WARNING("Event " << scanner.getCounter() << " is empty; skip.");
198  }
199  }
200 
202 
203  io >> outputFile;
204 
205  outputFile.close();
206 
207  return 0;
208 }
Auxiliary class for ROOT I/O of application specific meta data.
Definition: JMeta.hh:70
Object writing to file.
Utility class to parse command line options.
Definition: JParser.hh:1514
static const int WEIGHTLIST_RESCALED_EVENT_RATE
Rescaled event rate [s-1].
Definition: weightlist.hh:17
bool is_gseagen(const JHead &header)
Check for generator.
Definition: JHeadToolkit.hh:61
const JHead & getHeader() const
Get header.
Definition: JHead.hh:1256
#define STATUS(A)
Definition: JMessage.hh:63
Utility class to parse parameter values.
Definition: JProperties.hh:496
Implementation of event-weight factor interface.
Data structure for single set of oscillation parameters.
Simple data structure to support I/O of equations (see class JLANG::JEquation).
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
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:83
string outputFile
void createUUID()
Create UUID if not already set.
Definition: JHead.hh:1287
Type list.
Definition: JTypeList.hh:22
JEvtWeightFactorFunction< JFunction_t, JEvtWeightFactor_t > make_weightFactor(const JFunction_t &function)
Auxiliary method for creating an interface to an event-weight factor.
uuid_t header_uuid
UUID of header containing the event-weight information.
Definition: Evt.hh:35
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1989
Monte Carlo run header.
Definition: JHead.hh:1221
The Head class reflects the header of Monte-Carlo event files, which consists of keys (also referred ...
Definition: Head.hh:65
#define FATAL(A)
Definition: JMessage.hh:67
Auxiliary base class for list of file names.
General purpose class for object reading from a list of file names.
void copy(const Head &from, JHead &to)
Copy header from from to to.
Definition: JHead.cc:162
static const int WEIGHTLIST_NORMALISATION
Event rate normalisation.
Definition: weightlist.hh:16
Template event-weighter-associated file scanner.
then echo WARNING
Definition: JTuneHV.sh:91
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:484
JEvtWeightFactorFunction< JAtmosphericNeutrinoFlux, JFlux > make_atmosphericNeutrinoFluxFunction(const string &oscProbTableFile, const JOscParameters &oscParameters)
Auxiliary method for creating an interface to an atmospheric neutrino flux function using an oscillat...
int debug
debug level
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition: Evt.hh:20
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62