Jpp  debug
the software that should make you happy
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:69
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:2158
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1714
GAUSS_LEGENDRE::f1
const JPolynome f1(1.0, 2.0, 3.0)
Function.
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:1576
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition: JAAnetToolkit.hh:43
std
Definition: JSTDTypes.hh:14
TObject
Definition: JRoot.hh:19
Generated by doxygen 1.9.1