Jpp  16.0.1
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JEvtReweightMupage.cc
Go to the documentation of this file.
1 #include <iostream>
2 #include <iomanip>
3 #include <string>
4 #include <vector>
5 
8 
9 #include "Jeep/JPrint.hh"
10 #include "Jeep/JParser.hh"
11 #include "Jeep/JMessage.hh"
12 
13 #include "JLang/JToken.hh"
14 
15 #include "JAAnet/JHead.hh"
16 #include "JAAnet/JHeadToolkit.hh"
17 #include "JAAnet/JAAnetToolkit.hh"
20 
21 #include "JSupport/JMeta.hh"
22 #include "JSupport/JSupport.hh"
25 
26 #include "JGizmo/JGizmoToolkit.hh"
27 
28 
29 /**
30  * \file
31  * Program for reweighting mupage data according to a specifiable TFormula.
32  *
33  * The TFormula may take any number of parameters, but is limited to the physical variables\n
34  * defined in `JEvtWeightFactorMupage`. The indices of the physical variables are defined in\n
35  * the enum `JEvtWeightFactorMupage::variables`.
36  *
37  * \author bjung
38  */
39 int main(int argc, char **argv)
40 {
41  using namespace std;
42  using namespace JPP;
43 
44  typedef typename JTYPELIST<Head, Evt, JMeta>::typelist typelist;
45  typedef JToken<';'> JToken_t;
46 
47 
48  string inputFile;
49  string outputFile;
50 
51  string formula;
53 
54  int debug;
55 
56  try {
57 
58  JParser<> zap;
59 
60  zap['f'] = make_field(inputFile);
61  zap['o'] = make_field(outputFile);
62  zap['F'] = make_field(formula, "Reweighting formula, e.g.: \"[0] + [1]*x[0]\"") = "";
63  zap['@'] = make_field(parameters, "Parameter values, e.g.: \"p0 = 0.5; p1 = 0.2;\"") = JPARSER::initialised();
64  zap['d'] = make_field(debug) = 1;
65 
66  zap(argc, argv);
67  }
68  catch(const exception& error) {
69  FATAL(error.what() << endl);
70  }
71 
72 
73  // Create file scanner
74 
75  JEvtWeightFileScanner<> scanner(inputFile);
76 
77  const JHead& header = scanner.getHeader();
78 
79  if (!is_mupage(header)) {
80  FATAL("MC-header is incompatible with MUPAGE.");
81  }
82 
83 
84  // Define reweighting function
85 
86  JEvtWeightFactorMupage reweighter("reweighter", formula.c_str());
87 
88  if (reweighter.GetNdim() > JEvtWeightFactorMupage::NUMBER_OF_VARIABLES) {
89 
90  FATAL("Invalid formula " << formula << " (number of dimensions must be < " <<
91  JEvtWeightFactorMupage::NUMBER_OF_VARIABLES << "; see `JAANET::JEvtWeightFactorMupage`).");
92  }
93 
94  if (size_t(reweighter.GetNpar()) != parameters.size()) {
95 
96  FATAL("Not all parameters have been specified.");
97  }
98 
99 
100  // Read parameter settings
101 
102  try {
103 
104  for (vector<JToken_t>::const_iterator i = parameters.begin(); i != parameters.cend(); ++i) {
105 
106  const int index = getParameter(*i);
107 
108  reweighter.SetParameter(index, getValue(*i, 0));
109  }
110  }
111  catch (JLANG::JParseError& error) {
112  FATAL(error << endl);
113  }
114 
115 
116  // Set reweighting factors for mupage-files
117 
118  if (!scanner.setEvtWeightFactor(TRACK_TYPE_ANTIMUON, make_weightFactor(reweighter))) {
119  FATAL("Setting of reweighting factor was unsuccessful.");
120  }
121 
122 
123  // Store events with new weights
124 
125  JFileRecorder<typelist> out(outputFile.c_str());
126 
127  out.open();
128 
129  Head head;
130  copy(header, head);
131  out.put(head);
132 
133  out.put(JMeta(argc, argv));
134 
135  for (JMultipleFileScanner<JMeta> in(scanner.getFilelist()); in.hasNext(); ) {
136  out.put(*in.next());
137  }
138 
139  STATUS(RIGHT(15) << "event" << RIGHT(15) << "weight" << endl);
140 
141  while (scanner.hasNext()) {
142 
143  Evt* event = scanner.next();
144 
145  if (event != NULL) {
146 
147  double weight = scanner.getWeight(*event);
148 
149  STATUS(RIGHT (15) << scanner.getCounter() <<
150  SCIENTIFIC(15, 3) << weight << '\r'); DEBUG(endl);
151 
152  if (!isfinite(weight)) {
153 
154  WARNING("Non-finite reweighting factor " << weight <<
155  " for event " << scanner.getCounter() << "!");
156  }
157 
158  if (event->w.size() < RESCALED_WEIGHT_INDEX + 1) {
159  event->w.resize(RESCALED_WEIGHT_INDEX + 1);
160  }
161 
162  event->w.at(RESCALED_WEIGHT_INDEX) = weight;
163 
164  out.put(*event);
165 
166  } else {
167 
168  WARNING("Event " << scanner.getCounter() << " is empty; skip.");
169  }
170  }
171 
172  out.close();
173 
174  return 0;
175 }
bool is_mupage(const JHead &header)
Check for generator.
Definition: JHeadToolkit.hh:84
static const int RESCALED_WEIGHT_INDEX
Index of rescaled weight.
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:1500
#define WARNING(A)
Definition: JMessage.hh:65
double getValue(const JScale_t scale)
Get numerical value corresponding to scale.
Definition: JScale.hh:47
int main(int argc, char *argv[])
Definition: Main.cc:15
int getParameter(const std::string &text)
Get parameter number from text string.
ROOT TTree parameter settings of various packages.
bool setEvtWeightFactor(const int type, const JEvtWeightFactor_t &factor)
Set event-weight factor corresponding to a given PDG code.
const JHead & getHeader() const
Get header.
Definition: JHead.hh:1197
#define STATUS(A)
Definition: JMessage.hh:63
*fatal Wrong number of arguments esac JCookie sh typeset Z DETECTOR typeset Z SOURCE_RUN typeset Z TARGET_RUN set_variable PARAMETERS_FILE $WORKDIR parameters
Definition: diff-Tuna.sh:38
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:66
string outputFile
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.
I/O formatting auxiliaries.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1961
ROOT I/O of application specific meta data.
Definition of hit and track types and auxiliary methods for handling Monte Carlo data.
int debug
debug level
Definition: JSirene.cc:63
Implementation of reweighting factor for mupage events according to a specifiable ROOT TFormula...
General purpose messaging.
Monte Carlo run header.
Definition: JHead.hh:1164
double getWeight(const Evt &evt) const
Get weight of given event.
The Head class reflects the header of Monte-Carlo event files, which consists of keys (also referred ...
Definition: Head.hh:67
#define FATAL(A)
Definition: JMessage.hh:67
Scanning of objects from multiple files according a format that follows from the extension of each fi...
General purpose class for object reading from a list of file names.
Utility class to parse command line options.
Wrapper class around string.
Definition: JToken.hh:23
then usage $script< input_file >< detector_file > fi set_variable OUTPUT_DIR set_variable SELECTOR JDAQTimesliceL1 set_variable DEBUG case set_variable DEBUG
void copy(const Head &from, JHead &to)
Copy header from from to to.
Definition: JHead.cc:139
Exception for parsing value.
Definition: JException.hh:180
virtual void open(const char *file_name) override
Open file.
Template event-weighter-associated file scanner.
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:42
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:417