Jpp  18.5.2
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JOpera2D.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 "TH2.h"
12 #include "TF2.h"
13 #include "TProfile2D.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 2D histograms.\n
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  */
59 int 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::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 }
Utility class to parse command line options.
Definition: JParser.hh:1514
int main(int argc, char *argv[])
Definition: Main.cc:15
then echo Enter input within $TIMEOUT_S seconds echo n User name
Definition: JCookie.sh:42
Definition: JRoot.hh:19
string outputFile
const JPolynome f1(1.0, 2.0, 3.0)
Function.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1989
#define NOTICE(A)
Definition: JMessage.hh:64
#define ERROR(A)
Definition: JMessage.hh:66
General purpose messaging.
#define FATAL(A)
Definition: JMessage.hh:67
Utility class to parse command line options.
bool isTObject(const TKey *key)
Check if given key corresponds to a TObject.
TDirectory * getDirectory(const JRootObjectID &id)
Get TDirectory pointer.
int debug
debug level
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62