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

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

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <cmath>
#include "TROOT.h"
#include "TFile.h"
#include "TKey.h"
#include "TH3.h"
#include "TF3.h"
#include "TProfile3D.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 3D 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 JOpera3D.cc.

Function Documentation

◆ main()

int main ( int argc,
char ** argv )

Definition at line 59 of file JOpera3D.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<TH3*>(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 TH3* 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 TH3* h1 = dynamic_cast<TH3*>(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 = (TH3*) h1->Clone(opera.c_str());
165 else if (name == JOpera::SAME_AS_INPUT())
166 h3 = (TH3*) h1->Clone(h1->GetName());
167 else
168 h3 = (TH3*) 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 TH3* h1 = dynamic_cast<TH3*>(listOfObjects[0]);
196 TH3* h2 = dynamic_cast<TH3*>(listOfObjects[1]);
197
198 NOTICE(h1->GetName() << ' ' << opera << ' ' << h2->GetName() << endl);
199
200 if (name == JOpera::SAME_AS_OPERATION())
201 h3 = (TH3*) h1->Clone(opera.c_str());
202 else if (name == JOpera::SAME_AS_INPUT())
203 h3 = (TH3*) h1->Clone(h1->GetName());
204 else
205 h3 = (TH3*) 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 for (Int_t i3 = 1; i3 <= h1->GetNbinsZ(); ++i3) {
231
232 const Double_t x = h1->GetXaxis()->GetBinCenter(i1);
233 const Double_t y = h1->GetYaxis()->GetBinCenter(i2);
234 const Double_t z = h1->GetZaxis()->GetBinCenter(i3);
235
236 const Int_t j1 = h2->GetXaxis()->FindBin(x);
237 const Int_t j2 = h2->GetYaxis()->FindBin(y);
238 const Int_t j3 = h2->GetYaxis()->FindBin(z);
239
240 const Double_t w1 = h1->GetBinContent(i1,i2,i3);
241 const Double_t w2 = h2->GetBinContent(j1,j2,j3);
242
243 Double_t w3 = 0.0;
244
245 if (opera == JOpera::add()) {
246
247 w3 = w1 + w2;
248
249 } else if (opera == JOpera::subtract()) {
250
251 w3 = w1 - w2;
252
253 } else if (opera == JOpera::multiply()) {
254
255 w3 = w1 * w2;
256
257 } else if (opera == JOpera::divide()) {
258
259 if (w2 == 0.0) {
260 ERROR("Bin " << h2->GetName() << "[" << j1 << "," << j2 << "] empty." << endl);
261 continue;
262 }
263
264 w3 = w1 / w2;
265 }
266
267 h3->SetBinContent(i1, i2, i3, w3);
268 }
269 }
270 }
271
272 } else if (opera == JOpera::efficiency() ||
273 opera == JOpera::stdev() ||
274 opera == JOpera::sqrt()) {
275
276 for (Int_t i1 = 1; i1 <= h1->GetNbinsX(); ++i1) {
277 for (Int_t i2 = 1; i2 <= h1->GetNbinsY(); ++i2) {
278 for (Int_t i3 = 1; i3 <= h1->GetNbinsZ(); ++i3) {
279
280 const Double_t y1 = h1->GetBinContent(i1,i2,i3);
281 const Double_t y2 = h2->GetBinContent(i1,i2,i3);
282
283 Double_t w1 = h1->GetBinError(i1,i2,i3);
284 Double_t w2 = h2->GetBinError(i1,i2,i3);
285
286 Double_t y3 = 0.0;
287 Double_t w3 = 0.0;
288
289 if (opera == JOpera::efficiency()) {
290
291 if (y2 == 0.0) {
292 ERROR("Bin " << h2->GetName() << "[" << i1 << "," << i2 << "] empty." << endl);
293 continue;
294 }
295
296 y3 = y1 / y2;
297
298 if (y1 != 0.0 &&
299 y2 != 0.0) {
300
301 w1 /= y1;
302 w2 /= y2;
303
304 w3 = y3 * fabs(y3 - 1.0) * sqrt(w1*w1 + w2*w2);
305 }
306
307 } else if (opera == JOpera::stdev()) {
308
309 w3 = sqrt(w1*w1 + w2*w2);
310
311 if (w3 != 0.0) {
312 y3 = (y1 - y2) / w3;
313 w3 = 0.0;
314 }
315
316 } else if (opera == JOpera::sqrt()) {
317
318 y3 = (y1+y2) * (y1-y2);
319
320 if (y3 < 0.0)
321 y3 = -sqrt(-y3);
322 else
323 y3 = +sqrt(+y3);
324
325 w3 = 0.0;
326 }
327
328 h3->SetBinContent(i1, i2, i3, y3);
329 h3->SetBinError (i1, i2, i3, w3);
330 }
331 }
332 }
333 }
334
335 } else if (opera == JOpera::Add() ||
336 opera == JOpera::Multiply()) {
337
338 for (vector<TObject*>::iterator i = listOfObjects.begin(); i != listOfObjects.end(); ++i) {
339
340 TH3* h1 = dynamic_cast<TH3*>(*i);
341
342 NOTICE(h1->GetName() << ' ' << opera << endl);
343
344 if (h3 == NULL) {
345
346 if (name == JOpera::SAME_AS_OPERATION())
347 h3 = (TH3*) h1->Clone(opera.c_str());
348 else if (name == JOpera::SAME_AS_INPUT())
349 h3 = (TH3*) h1->Clone(h1->GetName());
350 else
351 h3 = (TH3*) h1->Clone(name.c_str());
352
353 } else {
354
355 if (opera == JOpera::Add()) {
356
357 h3->Add (h3, h1, +1.0, +1.0);
358
359 } else if (opera == JOpera::Multiply()) {
360
361 h3->Multiply(h3, h1, +1.0, +1.0);
362 }
363 }
364 }
365
366 } else {
367
368 FATAL("Incompatible number of histograms " << listOfObjects.size() << " with option " << opera << endl);
369 }
370
371 if (h3 != NULL) {
372
373 TFile out(outputFile.c_str(), "recreate");
374
375 h3->Write();
376
377 out.Write();
378 out.Close();
379 }
380}
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