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