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

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

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <cmath>
#include "TROOT.h"
#include "TFile.h"
#include "TKey.h"
#include "TH1.h"
#include "TF1.h"
#include "TProfile.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 1D 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 JOpera1D.cc.

Function Documentation

◆ main()

int main ( int argc,
char ** argv )

Definition at line 59 of file JOpera1D.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::sqrt(),
87 JOpera::Replace();
88 zap['o'] = make_field(outputFile, "ROOT file with histogram(s)") = "opera.root";
89 zap['O'] = make_field(name, "histogram name:"
90 << "\n\t\"" << JOpera::SAME_AS_OPERATION() << "\" -> same as operation; or"
91 << "\n\t\"" << JOpera::SAME_AS_INPUT() << "\" -> same as input; else"
92 << "\n\t as specified") = JOpera::SAME_AS_OPERATION();
93 zap['d'] = make_field(debug) = 1;
94
95 zap(argc, argv);
96 }
97 catch(const exception &error) {
98 FATAL(error.what() << endl);
99 }
100
101
102 vector<TObject*> listOfObjects;
103
104 for (vector<JRootObjectID>::const_iterator input = inputFile.begin(); input != inputFile.end(); ++input) {
105
106 DEBUG("Input: " << *input << endl);
107
108 TDirectory* dir = getDirectory(*input);
109
110 if (dir == NULL) {
111 ERROR("File: " << input->getFullFilename() << " not opened." << endl);
112 continue;
113 }
114
115 const TRegexp regexp(input->getObjectName());
116
117 TIter iter(dir->GetListOfKeys());
118
119 for (TKey* key; (key = (TKey*) iter.Next()) != NULL; ) {
120
121 const TString tag(key->GetName());
122
123 DEBUG("Key: " << tag << " match = " << tag.Contains(regexp) << endl);
124
125 // option match
126
127 if (tag.Contains(regexp) && isTObject(key)) {
128
129 TObject* p = key->ReadObj();
130 TProfile* q = dynamic_cast<TProfile*>(p);
131
132 if (q != NULL) {
133 p = q->ProjectionX();
134 }
135
136 if (dynamic_cast<TH1*>(p) == NULL) {
137 FATAL("Object " << p->GetName() << " not compatible with histogram operations." << endl);
138 }
139
140 listOfObjects.push_back(p);
141 }
142 }
143 }
144
145
146 TH1* h3 = NULL;
147
148 if (listOfObjects.size() == 0) {
149
150 FATAL("Number of histograms " << listOfObjects.size() << " = 0." << endl);
151
152 } else if (listOfObjects.size() == 1 && (opera == JOpera::Add() ||
153 opera == JOpera::Subtract() ||
154 opera == JOpera::Multiply() ||
155 opera == JOpera::Divide() ||
156 opera == JOpera::Replace())) {
157
158 // operation between histogram contents and first associated function
159
160 TH1* h1 = dynamic_cast<TH1*>(listOfObjects[0]);
161 TF1* f1 = (TF1*) h1->GetListOfFunctions()->First();
162
163 if (f1 == NULL) {
164 FATAL(h1->GetName() << " has no associated function." << endl);
165 }
166
167 NOTICE(h1->GetName() << ' ' << opera << ' ' << f1->GetName() << endl);
168
169 if (name == JOpera::SAME_AS_OPERATION())
170 h3 = (TH1*) h1->Clone(opera.c_str());
171 else if (name == JOpera::SAME_AS_INPUT())
172 h3 = (TH1*) h1->Clone(h1->GetName());
173 else
174 h3 = (TH1*) h1->Clone(name.c_str());
175
176 if (opera == JOpera::Add()) {
177
178 h3->Add (f1, +1.0);
179
180 } else if (opera == JOpera::Subtract()) {
181
182 h3->Add (f1, -1.0);
183
184 } else if (opera == JOpera::Multiply()) {
185
186 h3->Multiply(f1);
187
188 } else if (opera == JOpera::Divide()) {
189
190 h3->Divide (f1);
191
192 } else if (opera == JOpera::Replace()) {
193
194 h3->Reset();
195 h3->Fill(0.0, 0.0); // ensure one entry for drawing after reset
196 h3->Add (f1);
197 }
198
199 } else if (listOfObjects.size() == 2) {
200
201 // operation between two histograms
202
203 TH1* h1 = dynamic_cast<TH1*>(listOfObjects[0]);
204 TH1* h2 = dynamic_cast<TH1*>(listOfObjects[1]);
205
206 NOTICE(h1->GetName() << ' ' << opera << ' ' << h2->GetName() << endl);
207
208 if (name == JOpera::SAME_AS_OPERATION())
209 h3 = (TH1*) h1->Clone(opera.c_str());
210 else if (name == JOpera::SAME_AS_INPUT())
211 h3 = (TH1*) h1->Clone(h1->GetName());
212 else
213 h3 = (TH1*) h1->Clone(name.c_str());
214
215 if (opera == JOpera::Add()) {
216
217 h3->Add (h1, h2, +1.0, +1.0);
218
219 } else if (opera == JOpera::Subtract()) {
220
221 h3->Add (h1, h2, +1.0, -1.0);
222
223 } else if (opera == JOpera::Multiply()) {
224
225 h3->Multiply(h1, h2, +1.0, +1.0);
226
227 } else if (opera == JOpera::Divide()) {
228
229 h3->Divide (h1, h2, +1.0, +1.0);
230
231 } else if (opera == JOpera::add() ||
232 opera == JOpera::subtract() ||
233 opera == JOpera::multiply() ||
234 opera == JOpera::divide()) {
235
236 for (Int_t i = 1; i <= h1->GetNbinsX(); ++i) {
237
238 const Double_t x = h1->GetBinCenter (i);
239 const Double_t w1 = h1->GetBinContent(i);
240
241 const Int_t j = h2->FindBin(x);
242 const Double_t w2 = h2->GetBinContent(j);
243
244 Double_t w3 = 0.0;
245
246 if (opera == JOpera::add()) {
247
248 w3 = w1 + w2;
249
250 } else if (opera == JOpera::subtract()) {
251
252 w3 = w1 - w2;
253
254 } else if (opera == JOpera::multiply()) {
255
256 w3 = w1 * w2;
257
258 } else if (opera == JOpera::divide()) {
259
260 if (w2 == 0.0) {
261 ERROR("Bin " << h2->GetName() << "[" << j << "] empty." << endl);
262 continue;
263 }
264
265 w3 = w1 / w2;
266 }
267
268 h3->SetBinContent(i, w3);
269 }
270
271 } else if (opera == JOpera::efficiency() ||
272 opera == JOpera::stdev() ||
273 opera == JOpera::sqrt()) {
274
275 for (Int_t i = 1; i <= h1->GetNbinsX(); ++i) {
276
277 const Double_t y1 = h1->GetBinContent(i);
278 const Double_t y2 = h2->GetBinContent(i);
279
280 Double_t w1 = h1->GetBinError(i);
281 Double_t w2 = h2->GetBinError(i);
282
283 Double_t y3 = 0.0;
284 Double_t w3 = 0.0;
285
286 if (opera == JOpera::efficiency()) {
287
288 if (y2 == 0.0) {
289 ERROR("Bin " << h2->GetName() << "[" << i << "] empty." << endl);
290 continue;
291 }
292
293 y3 = y1 / y2;
294
295 if (y1 != 0.0 &&
296 y2 != 0.0) {
297
298 w1 /= y1;
299 w2 /= y2;
300
301 w3 = y3 * fabs(y3 - 1.0) * sqrt(w1*w1 + w2*w2);
302 }
303
304 } else if (opera == JOpera::stdev()) {
305
306 w3 = sqrt(w1*w1 + w2*w2);
307
308 if (w3 != 0.0) {
309 y3 = (y1 - y2) / w3;
310 w3 = 0.0;
311 }
312
313 } else if (opera == JOpera::sqrt()) {
314
315 y3 = (y1+y2) * (y1-y2);
316
317 if (y3 < 0.0)
318 y3 = -sqrt(-y3);
319 else
320 y3 = +sqrt(+y3);
321
322 w3 = 0.0;
323 }
324
325 h3->SetBinContent(i, y3);
326 h3->SetBinError (i, w3);
327 }
328 }
329
330 } else if (opera == JOpera::Add() ||
331 opera == JOpera::Multiply()) {
332
333 // operation between more than two histograms
334
335 for (vector<TObject*>::iterator i = listOfObjects.begin(); i != listOfObjects.end(); ++i) {
336
337 TH1* h1 = dynamic_cast<TH1*>(*i);
338
339 NOTICE(h1->GetName() << ' ' << opera << endl);
340
341 if (h3 == NULL) {
342
343 if (name == JOpera::SAME_AS_OPERATION())
344 h3 = (TH1*) h1->Clone(opera.c_str());
345 else if (name == JOpera::SAME_AS_INPUT())
346 h3 = (TH1*) h1->Clone(h1->GetName());
347 else
348 h3 = (TH1*) h1->Clone(name.c_str());
349
350 } else {
351
352 if (opera == JOpera::Add()) {
353
354 h3->Add (h3, h1, +1.0, +1.0);
355
356 } else if (opera == JOpera::Multiply()) {
357
358 h3->Multiply(h3, h1, +1.0, +1.0);
359 }
360 }
361 }
362
363 } else {
364
365 FATAL("Incompatible number of histograms " << listOfObjects.size() << " with option " << opera << endl);
366 }
367
368 if (h3 != NULL) {
369
370 TFile out(outputFile.c_str(), "recreate");
371
372 h3->Write();
373
374 out.Write();
375 out.Close();
376 }
377}
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
JTOOLS::j
int j
Definition JPolint.hh:801
Generated by doxygen 1.10.0