Jpp 19.3.0-rc.1
the software that should make you happy
Loading...
Searching...
No Matches
JOpera2D.cc
Go to the documentation of this file.
1
2#include <string>
3#include <iostream>
4#include <iomanip>
5#include <vector>
6#include <cmath>
7
8#include "TROOT.h"
9#include "TFile.h"
10#include "TKey.h"
11#include "TH2.h"
12#include "TF2.h"
13#include "TProfile2D.h"
14#include "TString.h"
15#include "TRegexp.h"
16
17#include "JGizmo/JRootObjectID.hh"
18#include "JGizmo/JGizmoToolkit.hh"
19
20#include "Jeep/JParser.hh"
21#include "Jeep/JMessage.hh"
22
23
24/**
25 * \file
26 * Auxiliary program for operations on 2D histograms.\n
27 * The option <tt>-f</tt> corresponds to <tt><file name>:<object name></tt>.
28 *
29 * Possible operations (option <tt>-u <operation></tt>:
30 * - Add\n
31 * ROOT addition of histograms;
32 * - add\n
33 * Addition of histograms by look up of corresponding bin in second histogram;
34 * - Subtract\n
35 * ROOT subtraction of histograms;
36 * - subtract\n
37 * Subtraction of histograms by look up of corresponding bin in second histogram;
38 * - Multiply\n
39 * ROOT multiplication of histograms;
40 * - multiply\n
41 * Multiplication of histograms by look up of corresponding bin in second histogram;
42 * - Divide\n
43 * ROOT division of histograms;
44 * - divide\n
45 * Division of histograms by look up of corresponding bin in second histogram;
46 * - efficiency\n
47 * Like "Divide", but with correction for errors;
48 * - stdev\n
49 * Standard deviation;
50 * - Replace\n
51 * Replace histogram contents by function value;
52 *
53 * If only one histogram is specified, the operation applies to the histogram and the first associated function.\n
54 * Only the quoted ROOT operations and the additional operation Replace are then allowed.\n
55 * The option Replace will replace the histogram contents by the function value.
56 *
57 * \author mdejong
58 */
59int main(int argc, char **argv)
60{
61 using namespace std;
62 using namespace JPP;
63
64 vector<JRootObjectID> inputFile;
65 string outputFile;
66 string opera;
67 string name;
68 int debug;
69
70 try {
71
72 JParser<> zap("Auxiliary program for histogram operations.");
73
74 zap['f'] = make_field(inputFile, "<input file>:<object name>");
75 zap['u'] = make_field(opera) =
76 JOpera::Add(),
77 JOpera::add(),
78 JOpera::Subtract(),
79 JOpera::subtract(),
80 JOpera::Multiply(),
81 JOpera::multiply(),
82 JOpera::Divide(),
83 JOpera::divide(),
84 JOpera::efficiency(),
85 JOpera::stdev(),
86 JOpera::Replace();
87 zap['o'] = make_field(outputFile, "ROOT file with histogram(s)") = "opera.root";
88 zap['O'] = make_field(name, "histogram name:"
89 << "\n\t\"" << JOpera::SAME_AS_OPERATION() << "\" -> same as operation; or"
90 << "\n\t\"" << JOpera::SAME_AS_INPUT() << "\" -> same as input; else"
91 << "\n\t as specified") = JOpera::SAME_AS_OPERATION();
92 zap['d'] = make_field(debug) = 1;
93
94 zap(argc, argv);
95 }
96 catch(const exception &error) {
97 FATAL(error.what() << endl);
98 }
99
100
101 vector<TObject*> listOfObjects;
102
103 for (vector<JRootObjectID>::const_iterator input = inputFile.begin(); input != inputFile.end(); ++input) {
104
105 DEBUG("Input: " << *input << endl);
106
107 TDirectory* dir = getDirectory(*input);
108
109 if (dir == NULL) {
110 ERROR("File: " << input->getFullFilename() << " not opened." << endl);
111 continue;
112 }
113
114 const TRegexp regexp(input->getObjectName());
115
116 TIter iter(dir->GetListOfKeys());
117
118 for (TKey* key; (key = (TKey*) iter.Next()) != NULL; ) {
119
120 const TString tag(key->GetName());
121
122 DEBUG("Key: " << tag << " match = " << tag.Contains(regexp) << endl);
123
124 // option match
125
126 if (tag.Contains(regexp) && isTObject(key)) {
127
128 TObject* p = key->ReadObj();
129
130 if (dynamic_cast<TH2*>(p) == NULL) {
131 FATAL("Object " << p->GetName() << " not compatible with histogram operations." << endl);
132 }
133
134 listOfObjects.push_back(p);
135 }
136 }
137 }
138
139
140 TH2* h3 = NULL;
141
142 if (listOfObjects.size() == 0) {
143
144 FATAL("Number of histograms " << listOfObjects.size() << " = 0." << endl);
145
146 } else if (listOfObjects.size() == 1 && (opera == JOpera::Add() ||
147 opera == JOpera::Subtract() ||
148 opera == JOpera::Multiply() ||
149 opera == JOpera::Divide() ||
150 opera == JOpera::Replace())) {
151
152 // operation between histogram contents and first associated function
153
154 TH2* h1 = dynamic_cast<TH2*>(listOfObjects[0]);
155 TF2* f1 = dynamic_cast<TF2*>(h1->GetListOfFunctions()->First());
156
157 if (f1 == NULL) {
158 FATAL(h1->GetName() << " has no associated function." << endl);
159 }
160
161 NOTICE(h1->GetName() << ' ' << opera << ' ' << f1->GetName() << endl);
162
163 if (name == JOpera::SAME_AS_OPERATION())
164 h3 = (TH2*) h1->Clone(opera.c_str());
165 else if (name == JOpera::SAME_AS_INPUT())
166 h3 = (TH2*) h1->Clone(h1->GetName());
167 else
168 h3 = (TH2*) h1->Clone(name.c_str());
169
170 if (opera == JOpera::Add()) {
171
172 h3->Add (f1, +1.0);
173
174 } else if (opera == JOpera::Subtract()) {
175
176 h3->Add (f1, -1.0);
177
178 } else if (opera == JOpera::Multiply()) {
179
180 h3->Multiply(f1);
181
182 } else if (opera == JOpera::Divide()) {
183
184 h3->Divide (f1);
185
186 } else if (opera == JOpera::Replace()) {
187
188 h3->Reset();
189 h3->Fill(0.0, 0.0, 0.0); // ensure one entry for drawing after reset
190 h3->Add (f1);
191 }
192
193 } else if (listOfObjects.size() == 2) {
194
195 TH2* h1 = dynamic_cast<TH2*>(listOfObjects[0]);
196 TH2* h2 = dynamic_cast<TH2*>(listOfObjects[1]);
197
198 NOTICE(h1->GetName() << ' ' << opera << ' ' << h2->GetName() << endl);
199
200 if (name == JOpera::SAME_AS_OPERATION())
201 h3 = (TH2*) h1->Clone(opera.c_str());
202 else if (name == JOpera::SAME_AS_INPUT())
203 h3 = (TH2*) h1->Clone(h1->GetName());
204 else
205 h3 = (TH2*) h1->Clone(name.c_str());
206
207 if (opera == JOpera::Add()) {
208
209 h3->Add (h1, h2, +1.0, +1.0);
210
211 } else if (opera == JOpera::Subtract()) {
212
213 h3->Add (h1, h2, +1.0, -1.0);
214
215 } else if (opera == JOpera::Multiply()) {
216
217 h3->Multiply(h1, h2, +1.0, +1.0);
218
219 } else if (opera == JOpera::Divide()) {
220
221 h3->Divide (h1, h2, +1.0, +1.0);
222
223 } else if (opera == JOpera::add() ||
224 opera == JOpera::subtract() ||
225 opera == JOpera::multiply() ||
226 opera == JOpera::divide()) {
227
228 for (Int_t i1 = 1; i1 <= h1->GetNbinsX(); ++i1) {
229 for (Int_t i2 = 1; i2 <= h1->GetNbinsY(); ++i2) {
230
231 const Double_t x = h1->GetXaxis()->GetBinCenter(i1);
232 const Double_t y = h1->GetYaxis()->GetBinCenter(i2);
233
234 const Int_t j1 = h2->GetXaxis()->FindBin(x);
235 const Int_t j2 = h2->GetYaxis()->FindBin(y);
236
237 const Double_t w1 = h1->GetBinContent(i1,i2);
238 const Double_t w2 = h2->GetBinContent(j1,j2);
239
240 Double_t w3 = 0.0;
241
242 if (opera == JOpera::add()) {
243
244 w3 = w1 + w2;
245
246 } else if (opera == JOpera::subtract()) {
247
248 w3 = w1 - w2;
249
250 } else if (opera == JOpera::multiply()) {
251
252 w3 = w1 * w2;
253
254 } else if (opera == JOpera::divide()) {
255
256 if (w2 == 0.0) {
257 ERROR("Bin " << h2->GetName() << "[" << j1 << "," << j2 << "] empty." << endl);
258 continue;
259 }
260
261 w3 = w1 / w2;
262 }
263
264 h3->SetBinContent(i1, i2, w3);
265 }
266 }
267
268 } else if (opera == JOpera::efficiency() ||
269 opera == JOpera::stdev() ||
270 opera == JOpera::sqrt()) {
271
272 for (Int_t i1 = 1; i1 <= h1->GetNbinsX(); ++i1) {
273 for (Int_t i2 = 1; i2 <= h1->GetNbinsY(); ++i2) {
274
275 const Double_t y1 = h1->GetBinContent(i1,i2);
276 const Double_t y2 = h2->GetBinContent(i1,i2);
277
278 Double_t w1 = h1->GetBinError(i1,i2);
279 Double_t w2 = h2->GetBinError(i1,i2);
280
281 Double_t y3 = 0.0;
282 Double_t w3 = 0.0;
283
284 if (opera == JOpera::efficiency()) {
285
286 if (y2 == 0.0) {
287 ERROR("Bin " << h2->GetName() << "[" << i1 << "," << i2 << "] empty." << endl);
288 continue;
289 }
290
291 y3 = y1 / y2;
292
293 if (y1 != 0.0 &&
294 y2 != 0.0) {
295
296 w1 /= y1;
297 w2 /= y2;
298
299 w3 = y3 * fabs(y3 - 1.0) * sqrt(w1*w1 + w2*w2);
300 }
301
302 } else if (opera == JOpera::stdev()) {
303
304 w3 = sqrt(w1*w1 + w2*w2);
305
306 if (w3 != 0.0) {
307 y3 = (y1 - y2) / w3;
308 w3 = 0.0;
309 }
310
311 } else if (opera == JOpera::sqrt()) {
312
313 y3 = (y1+y2) * (y1-y2);
314
315 if (y3 < 0.0)
316 y3 = -sqrt(-y3);
317 else
318 y3 = +sqrt(+y3);
319
320 w3 = 0.0;
321 }
322
323 h3->SetBinContent(i1, i2, y3);
324 h3->SetBinError (i1, i2, w3);
325 }
326 }
327 }
328
329 } else if (opera == JOpera::Add() ||
330 opera == JOpera::Multiply()) {
331
332 for (vector<TObject*>::iterator i = listOfObjects.begin(); i != listOfObjects.end(); ++i) {
333
334 TH2* h1 = dynamic_cast<TH2*>(*i);
335
336 NOTICE(h1->GetName() << ' ' << opera << endl);
337
338 if (h3 == NULL) {
339
340 if (name == JOpera::SAME_AS_OPERATION())
341 h3 = (TH2*) h1->Clone(opera.c_str());
342 else if (name == JOpera::SAME_AS_INPUT())
343 h3 = (TH2*) h1->Clone(h1->GetName());
344 else
345 h3 = (TH2*) h1->Clone(name.c_str());
346
347 } else {
348
349 if (opera == JOpera::Add()) {
350
351 h3->Add (h3, h1, +1.0, +1.0);
352
353 } else if (opera == JOpera::Multiply()) {
354
355 h3->Multiply(h3, h1, +1.0, +1.0);
356 }
357 }
358 }
359
360 } else {
361
362 FATAL("Incompatible number of histograms " << listOfObjects.size() << " with option " << opera << endl);
363 }
364
365 if (h3 != NULL) {
366
367 TFile out(outputFile.c_str(), "recreate");
368
369 h3->Write();
370
371 out.Write();
372 out.Close();
373 }
374}
outputFile
string outputFile
Definition JDAQTimesliceSelector.cc:37
JMessage.hh
General purpose messaging.
DEBUG
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
ERROR
#define ERROR(A)
Definition JMessage.hh:66
NOTICE
#define NOTICE(A)
Definition JMessage.hh:64
FATAL
#define FATAL(A)
Definition JMessage.hh:67
debug
int debug
debug level
Definition JSirene.cc:72
main
int main(int argc, char **argv)
Definition JOpera2D.cc:59
JParser.hh
Utility class to parse command line options.
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
JPARSER::JParser
Utility class to parse command line options.
Definition JParser.hh:1698
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition JAAnetToolkit.hh:43
Generated by doxygen 1.12.0