Jpp test-rotations-old
the software that should make you happy
Loading...
Searching...
No Matches
JOpera1D.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 "TH1.h"
12#include "TF1.h"
13#include "TProfile.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 1D histograms.
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::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
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 JOpera1D.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