Jpp 20.0.0-27-g39925593c-D
the software that should make you happy
Loading...
Searching...
No Matches
Functions
JOpera2D.cc File Reference

Auxiliary program for operations on 2D histograms. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <cmath>
#include "TROOT.h"
#include "TFile.h"
#include "TKey.h"
#include "TH2.h"
#include "TF2.h"
#include "TProfile2D.h"
#include "TString.h"
#include "TRegexp.h"
#include "JGizmo/JRootObjectID.hh"
#include "JGizmo/JGizmoToolkit.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Detailed Description

Auxiliary program for operations on 2D histograms.


The option -f corresponds to <file name>:<object name>.

Possible operations (option -u <operation>:

If only one histogram is specified, the operation applies to the histogram and the first associated function.
Only the quoted ROOT operations and the additional operation Replace are then allowed.
The option Replace will replace the histogram contents by the function value.

Author
mdejong

Definition in file JOpera2D.cc.

Function Documentation

◆ main()

int main ( int argc,
char ** argv )

Definition at line 59 of file JOpera2D.cc.

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
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
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
JOSCPROB::JOscProbInterpolator
Template definition of a multi-dimensional oscillation probability interpolation table.
Definition JOscProbInterpolator.hh:71
JPARSER::JParser
Utility class to parse command line options.
Definition JParser.hh:1698
JGIZMO::getDirectory
TDirectory * getDirectory(const JRootObjectID &id)
Get TDirectory pointer.
Definition JGizmoToolkit.hh:144
JGIZMO::isTObject
bool isTObject(const TKey *key)
Check if given key corresponds to a TObject.
Definition JGizmoToolkit.hh:1652
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition JAAnetToolkit.hh:45
Generated by doxygen 1.10.0