Jpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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

int main ( int  argc,
char **  argv 
)

Definition at line 57 of file JOpera2D.cc.

58 {
59  using namespace std;
60  using namespace JPP;
61 
62  vector<JRootObjectID> inputFile;
63  string outputFile;
64  string opera;
65  string name;
66  int debug;
67 
68  try {
69 
70  JParser<> zap("Auxiliary program for histogram operations.");
71 
72  zap['f'] = make_field(inputFile, "<input file>:<object name>");
73  zap['u'] = make_field(opera) =
74  JOpera::Add(),
75  JOpera::add(),
76  JOpera::Subtract(),
77  JOpera::subtract(),
78  JOpera::Multiply(),
79  JOpera::multiply(),
80  JOpera::Divide(),
81  JOpera::divide(),
82  JOpera::efficiency(),
83  JOpera::stdev(),
84  JOpera::Replace();
85  zap['o'] = make_field(outputFile, "ROOT file with histogram(s)") = "opera.root";
86  zap['O'] = make_field(name, MAKE_STRING("histogram name;"
87  << "\n\"" << JOpera::SAME_AS_OPERATION() << "\" -> same as operation; or"
88  << "\n\"" << JOpera::SAME_AS_INPUT() << "\" -> same as input; else"
89  << "\nas specified")) = JOpera::SAME_AS_OPERATION();
90  zap['d'] = make_field(debug) = 1;
91 
92  zap(argc, argv);
93  }
94  catch(const exception &error) {
95  FATAL(error.what() << endl);
96  }
97 
98 
99  vector<TObject*> listOfObjects;
100 
101  for (vector<JRootObjectID>::const_iterator input = inputFile.begin(); input != inputFile.end(); ++input) {
102 
103  DEBUG("Input: " << *input << endl);
104 
105  TDirectory* dir = getDirectory(*input);
106 
107  if (dir == NULL) {
108  ERROR("File: " << input->getFullFilename() << " not opened." << endl);
109  continue;
110  }
111 
112  const TRegexp regexp(input->getObjectName());
113 
114  TIter iter(dir->GetListOfKeys());
115 
116  for (TKey* key; (key = (TKey*) iter.Next()) != NULL; ) {
117 
118  const TString tag(key->GetName());
119 
120  DEBUG("Key: " << tag << " match = " << tag.Contains(regexp) << endl);
121 
122  // option match
123 
124  if (tag.Contains(regexp)) {
125 
126  TObject* p = key->ReadObj();
127 
128  if (dynamic_cast<TH2*>(p) == NULL) {
129  FATAL("Object " << p->GetName() << " not compatible with histogram operations." << endl);
130  }
131 
132  listOfObjects.push_back(p);
133  }
134  }
135  }
136 
137 
138  TH2* h3 = NULL;
139 
140  if (listOfObjects.size() == 0) {
141 
142  FATAL("Number of histograms " << listOfObjects.size() << " = 0." << endl);
143 
144  } else if (listOfObjects.size() == 1 && (opera == JOpera::Add() ||
145  opera == JOpera::Subtract() ||
146  opera == JOpera::Multiply() ||
147  opera == JOpera::Divide() ||
148  opera == JOpera::Replace())) {
149 
150  // operation between histogram contents and first associated function
151 
152  TH2* h1 = dynamic_cast<TH2*>(listOfObjects[0]);
153  TF2* f1 = dynamic_cast<TF2*>(h1->GetListOfFunctions()->First());
154 
155  if (f1 == NULL) {
156  FATAL(h1->GetName() << " has no associated function." << endl);
157  }
158 
159  NOTICE(h1->GetName() << ' ' << opera << ' ' << f1->GetName() << endl);
160 
161  if (name == JOpera::SAME_AS_OPERATION())
162  h3 = (TH2*) h1->Clone(opera.c_str());
163  else if (name == JOpera::SAME_AS_INPUT())
164  h3 = (TH2*) h1->Clone(h1->GetName());
165  else
166  h3 = (TH2*) h1->Clone(name.c_str());
167 
168  if (opera == JOpera::Add()) {
169 
170  h3->Add (f1, +1.0);
171 
172  } else if (opera == JOpera::Subtract()) {
173 
174  h3->Add (f1, -1.0);
175 
176  } else if (opera == JOpera::Multiply()) {
177 
178  h3->Multiply(f1);
179 
180  } else if (opera == JOpera::Divide()) {
181 
182  h3->Divide (f1);
183 
184  } else if (opera == JOpera::Replace()) {
185 
186  h3->Reset();
187  h3->Fill(0.0, 0.0, 0.0); // ensure one entry for drawing after reset
188  h3->Add (f1);
189  }
190 
191  } else if (listOfObjects.size() == 2) {
192 
193  TH2* h1 = dynamic_cast<TH2*>(listOfObjects[0]);
194  TH2* h2 = dynamic_cast<TH2*>(listOfObjects[1]);
195 
196  NOTICE(h1->GetName() << ' ' << opera << ' ' << h2->GetName() << endl);
197 
198  if (name == JOpera::SAME_AS_OPERATION())
199  h3 = (TH2*) h1->Clone(opera.c_str());
200  else if (name == JOpera::SAME_AS_INPUT())
201  h3 = (TH2*) h1->Clone(h1->GetName());
202  else
203  h3 = (TH2*) h1->Clone(name.c_str());
204 
205  if (opera == JOpera::Add()) {
206 
207  h3->Add (h1, h2, +1.0, +1.0);
208 
209  } else if (opera == JOpera::Subtract()) {
210 
211  h3->Add (h1, h2, +1.0, -1.0);
212 
213  } else if (opera == JOpera::Multiply()) {
214 
215  h3->Multiply(h1, h2, +1.0, +1.0);
216 
217  } else if (opera == JOpera::Divide()) {
218 
219  h3->Divide (h1, h2, +1.0, +1.0);
220 
221  } else if (opera == JOpera::add() ||
222  opera == JOpera::subtract() ||
223  opera == JOpera::multiply() ||
224  opera == JOpera::divide()) {
225 
226  for (Int_t i1 = 1; i1 <= h1->GetNbinsX(); ++i1) {
227  for (Int_t i2 = 1; i2 <= h1->GetNbinsY(); ++i2) {
228 
229  const Double_t x = h1->GetXaxis()->GetBinCenter(i1);
230  const Double_t y = h1->GetYaxis()->GetBinCenter(i2);
231 
232  const Int_t j1 = h2->GetXaxis()->FindBin(x);
233  const Int_t j2 = h2->GetYaxis()->FindBin(y);
234 
235  const Double_t w1 = h1->GetBinContent(i1,i2);
236  const Double_t w2 = h2->GetBinContent(j1,j2);
237 
238  Double_t w3 = 0.0;
239 
240  if (opera == JOpera::add()) {
241 
242  w3 = w1 + w2;
243 
244  } else if (opera == JOpera::subtract()) {
245 
246  w3 = w1 - w2;
247 
248  } else if (opera == JOpera::multiply()) {
249 
250  w3 = w1 * w2;
251 
252  } else if (opera == JOpera::divide()) {
253 
254  if (w2 == 0.0) {
255  ERROR("Bin " << h2->GetName() << "[" << j1 << "," << j2 << "] empty." << endl);
256  continue;
257  }
258 
259  w3 = w1 / w2;
260  }
261 
262  h3->SetBinContent(i1, i2, w3);
263  }
264  }
265 
266  } else if (opera == JOpera::efficiency() ||
267  opera == JOpera::stdev() ||
268  opera == JOpera::sqrt()) {
269 
270  for (Int_t i1 = 1; i1 <= h1->GetNbinsX(); ++i1) {
271  for (Int_t i2 = 1; i2 <= h1->GetNbinsY(); ++i2) {
272 
273  const Double_t y1 = h1->GetBinContent(i1,i2);
274  const Double_t y2 = h2->GetBinContent(i1,i2);
275 
276  Double_t w1 = h1->GetBinError(i1,i2);
277  Double_t w2 = h2->GetBinError(i1,i2);
278 
279  Double_t y3 = 0.0;
280  Double_t w3 = 0.0;
281 
282  if (opera == JOpera::efficiency()) {
283 
284  if (y2 == 0.0) {
285  ERROR("Bin " << h2->GetName() << "[" << i1 << "," << i2 << "] empty." << endl);
286  continue;
287  }
288 
289  y3 = y1 / y2;
290 
291  if (y1 != 0.0 &&
292  y2 != 0.0) {
293 
294  w1 /= y1;
295  w2 /= y2;
296 
297  w3 = y3 * fabs(y3 - 1.0) * sqrt(w1*w1 + w2*w2);
298  }
299 
300  } else if (opera == JOpera::stdev()) {
301 
302  w3 = sqrt(w1*w1 + w2*w2);
303 
304  if (w3 != 0.0) {
305  y3 = (y1 - y2) / w3;
306  w3 = 0.0;
307  }
308 
309  } else if (opera == JOpera::sqrt()) {
310 
311  y3 = (y1+y2) * (y1-y2);
312 
313  if (y3 < 0.0)
314  y3 = -sqrt(-y3);
315  else
316  y3 = +sqrt(+y3);
317 
318  w3 = 0.0;
319  }
320 
321  h3->SetBinContent(i1, i2, y3);
322  h3->SetBinError (i1, i2, w3);
323  }
324  }
325  }
326 
327  } else if (opera == JOpera::Add() ||
328  opera == JOpera::Multiply()) {
329 
330  for (vector<TObject*>::iterator i = listOfObjects.begin(); i != listOfObjects.end(); ++i) {
331 
332  TH2* h1 = dynamic_cast<TH2*>(*i);
333 
334  NOTICE(h1->GetName() << ' ' << opera << endl);
335 
336  if (h3 == NULL) {
337 
338  if (name == JOpera::SAME_AS_OPERATION())
339  h3 = (TH2*) h1->Clone(opera.c_str());
340  else if (name == JOpera::SAME_AS_INPUT())
341  h3 = (TH2*) h1->Clone(h1->GetName());
342  else
343  h3 = (TH2*) h1->Clone(name.c_str());
344 
345  } else {
346 
347  if (opera == JOpera::Add()) {
348 
349  h3->Add (h3, h1, +1.0, +1.0);
350 
351  } else if (opera == JOpera::Multiply()) {
352 
353  h3->Multiply(h3, h1, +1.0, +1.0);
354  }
355  }
356  }
357 
358  } else {
359 
360  FATAL("Incompatible number of histograms " << listOfObjects.size() << " with option " << opera << endl);
361  }
362 
363  if (h3 != NULL) {
364 
365  TFile out(outputFile.c_str(), "recreate");
366 
367  h3->Write();
368 
369  out.Write();
370  out.Close();
371  }
372 }
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1500
TObject
Definition: JRoot.hh:19
h1
then for HISTOGRAM in h0 h1
Definition: JMatrixNZ.sh:69
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
MAKE_STRING
#define MAKE_STRING(A)
Make string.
Definition: JPrint.hh:142
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1961
NOTICE
#define NOTICE(A)
Definition: JMessage.hh:64
ERROR
#define ERROR(A)
Definition: JMessage.hh:66
debug
int debug
debug level
Definition: JSirene.cc:63
FATAL
#define FATAL(A)
Definition: JMessage.hh:67
name
then echo n User name
Definition: JCookie.sh:45
JGIZMO::getDirectory
TDirectory * getDirectory(const JRootObjectID &id)
Get TDirectory pointer.
Definition: JGizmoToolkit.hh:129
DEBUG
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
Generated by   doxygen 1.8.5