Jpp  18.0.1-rc.1
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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

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 }
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1514
name
then echo Enter input within $TIMEOUT_S seconds echo n User name
Definition: JCookie.sh:42
TObject
Definition: JRoot.hh:19
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
GAUSS_LEGENDRE::f1
const JPolynome f1(1.0, 2.0, 3.0)
Function.
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1989
NOTICE
#define NOTICE(A)
Definition: JMessage.hh:64
ERROR
#define ERROR(A)
Definition: JMessage.hh:66
FATAL
#define FATAL(A)
Definition: JMessage.hh:67
JGIZMO::isTObject
bool isTObject(const TKey *key)
Check if given key corresponds to a TObject.
Definition: JGizmoToolkit.hh:1517
JGIZMO::getDirectory
TDirectory * getDirectory(const JRootObjectID &id)
Get TDirectory pointer.
Definition: JGizmoToolkit.hh:137
debug
int debug
debug level
Definition: archive-put-wiki-detectors.sh:92
DEBUG
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
Generated by   doxygen 1.8.5