Jpp  18.0.0
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions
JOpera1D.cc File Reference

Auxiliary program for operations on 1D histograms. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <cmath>
#include "TROOT.h"
#include "TFile.h"
#include "TKey.h"
#include "TH1.h"
#include "TF1.h"
#include "TProfile.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 1D 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 JOpera1D.cc.

Function Documentation

int main ( int  argc,
char **  argv 
)

Definition at line 59 of file JOpera1D.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::sqrt(),
87  JOpera::Replace();
88  zap['o'] = make_field(outputFile, "ROOT file with histogram(s)") = "opera.root";
89  zap['O'] = make_field(name, "histogram name:"
90  << "\n\t\"" << JOpera::SAME_AS_OPERATION() << "\" -> same as operation; or"
91  << "\n\t\"" << JOpera::SAME_AS_INPUT() << "\" -> same as input; else"
92  << "\n\t as specified") = JOpera::SAME_AS_OPERATION();
93  zap['d'] = make_field(debug) = 1;
94 
95  zap(argc, argv);
96  }
97  catch(const exception &error) {
98  FATAL(error.what() << endl);
99  }
100 
101 
102  vector<TObject*> listOfObjects;
103 
104  for (vector<JRootObjectID>::const_iterator input = inputFile.begin(); input != inputFile.end(); ++input) {
105 
106  DEBUG("Input: " << *input << endl);
107 
108  TDirectory* dir = getDirectory(*input);
109 
110  if (dir == NULL) {
111  ERROR("File: " << input->getFullFilename() << " not opened." << endl);
112  continue;
113  }
114 
115  const TRegexp regexp(input->getObjectName());
116 
117  TIter iter(dir->GetListOfKeys());
118 
119  for (TKey* key; (key = (TKey*) iter.Next()) != NULL; ) {
120 
121  const TString tag(key->GetName());
122 
123  DEBUG("Key: " << tag << " match = " << tag.Contains(regexp) << endl);
124 
125  // option match
126 
127  if (tag.Contains(regexp) && isTObject(key)) {
128 
129  TObject* p = key->ReadObj();
130  TProfile* q = dynamic_cast<TProfile*>(p);
131 
132  if (q != NULL) {
133  p = q->ProjectionX();
134  }
135 
136  if (dynamic_cast<TH1*>(p) == NULL) {
137  FATAL("Object " << p->GetName() << " not compatible with histogram operations." << endl);
138  }
139 
140  listOfObjects.push_back(p);
141  }
142  }
143  }
144 
145 
146  TH1* h3 = NULL;
147 
148  if (listOfObjects.size() == 0) {
149 
150  FATAL("Number of histograms " << listOfObjects.size() << " = 0." << endl);
151 
152  } else if (listOfObjects.size() == 1 && (opera == JOpera::Add() ||
153  opera == JOpera::Subtract() ||
154  opera == JOpera::Multiply() ||
155  opera == JOpera::Divide() ||
156  opera == JOpera::Replace())) {
157 
158  // operation between histogram contents and first associated function
159 
160  TH1* h1 = dynamic_cast<TH1*>(listOfObjects[0]);
161  TF1* f1 = (TF1*) h1->GetListOfFunctions()->First();
162 
163  if (f1 == NULL) {
164  FATAL(h1->GetName() << " has no associated function." << endl);
165  }
166 
167  NOTICE(h1->GetName() << ' ' << opera << ' ' << f1->GetName() << endl);
168 
169  if (name == JOpera::SAME_AS_OPERATION())
170  h3 = (TH1*) h1->Clone(opera.c_str());
171  else if (name == JOpera::SAME_AS_INPUT())
172  h3 = (TH1*) h1->Clone(h1->GetName());
173  else
174  h3 = (TH1*) h1->Clone(name.c_str());
175 
176  if (opera == JOpera::Add()) {
177 
178  h3->Add (f1, +1.0);
179 
180  } else if (opera == JOpera::Subtract()) {
181 
182  h3->Add (f1, -1.0);
183 
184  } else if (opera == JOpera::Multiply()) {
185 
186  h3->Multiply(f1);
187 
188  } else if (opera == JOpera::Divide()) {
189 
190  h3->Divide (f1);
191 
192  } else if (opera == JOpera::Replace()) {
193 
194  h3->Reset();
195  h3->Fill(0.0, 0.0); // ensure one entry for drawing after reset
196  h3->Add (f1);
197  }
198 
199  } else if (listOfObjects.size() == 2) {
200 
201  // operation between two histograms
202 
203  TH1* h1 = dynamic_cast<TH1*>(listOfObjects[0]);
204  TH1* h2 = dynamic_cast<TH1*>(listOfObjects[1]);
205 
206  NOTICE(h1->GetName() << ' ' << opera << ' ' << h2->GetName() << endl);
207 
208  if (name == JOpera::SAME_AS_OPERATION())
209  h3 = (TH1*) h1->Clone(opera.c_str());
210  else if (name == JOpera::SAME_AS_INPUT())
211  h3 = (TH1*) h1->Clone(h1->GetName());
212  else
213  h3 = (TH1*) h1->Clone(name.c_str());
214 
215  if (opera == JOpera::Add()) {
216 
217  h3->Add (h1, h2, +1.0, +1.0);
218 
219  } else if (opera == JOpera::Subtract()) {
220 
221  h3->Add (h1, h2, +1.0, -1.0);
222 
223  } else if (opera == JOpera::Multiply()) {
224 
225  h3->Multiply(h1, h2, +1.0, +1.0);
226 
227  } else if (opera == JOpera::Divide()) {
228 
229  h3->Divide (h1, h2, +1.0, +1.0);
230 
231  } else if (opera == JOpera::add() ||
232  opera == JOpera::subtract() ||
233  opera == JOpera::multiply() ||
234  opera == JOpera::divide()) {
235 
236  for (Int_t i = 1; i <= h1->GetNbinsX(); ++i) {
237 
238  const Double_t x = h1->GetBinCenter (i);
239  const Double_t w1 = h1->GetBinContent(i);
240 
241  const Int_t j = h2->FindBin(x);
242  const Double_t w2 = h2->GetBinContent(j);
243 
244  Double_t w3 = 0.0;
245 
246  if (opera == JOpera::add()) {
247 
248  w3 = w1 + w2;
249 
250  } else if (opera == JOpera::subtract()) {
251 
252  w3 = w1 - w2;
253 
254  } else if (opera == JOpera::multiply()) {
255 
256  w3 = w1 * w2;
257 
258  } else if (opera == JOpera::divide()) {
259 
260  if (w2 == 0.0) {
261  ERROR("Bin " << h2->GetName() << "[" << j << "] empty." << endl);
262  continue;
263  }
264 
265  w3 = w1 / w2;
266  }
267 
268  h3->SetBinContent(i, w3);
269  }
270 
271  } else if (opera == JOpera::efficiency() ||
272  opera == JOpera::stdev() ||
273  opera == JOpera::sqrt()) {
274 
275  for (Int_t i = 1; i <= h1->GetNbinsX(); ++i) {
276 
277  const Double_t y1 = h1->GetBinContent(i);
278  const Double_t y2 = h2->GetBinContent(i);
279 
280  Double_t w1 = h1->GetBinError(i);
281  Double_t w2 = h2->GetBinError(i);
282 
283  Double_t y3 = 0.0;
284  Double_t w3 = 0.0;
285 
286  if (opera == JOpera::efficiency()) {
287 
288  if (y2 == 0.0) {
289  ERROR("Bin " << h2->GetName() << "[" << i << "] empty." << endl);
290  continue;
291  }
292 
293  y3 = y1 / y2;
294 
295  if (y1 != 0.0 &&
296  y2 != 0.0) {
297 
298  w1 /= y1;
299  w2 /= y2;
300 
301  w3 = y3 * fabs(y3 - 1.0) * sqrt(w1*w1 + w2*w2);
302  }
303 
304  } else if (opera == JOpera::stdev()) {
305 
306  w3 = sqrt(w1*w1 + w2*w2);
307 
308  if (w3 != 0.0) {
309  y3 = (y1 - y2) / w3;
310  w3 = 0.0;
311  }
312 
313  } else if (opera == JOpera::sqrt()) {
314 
315  y3 = (y1+y2) * (y1-y2);
316 
317  if (y3 < 0.0)
318  y3 = -sqrt(-y3);
319  else
320  y3 = +sqrt(+y3);
321 
322  w3 = 0.0;
323  }
324 
325  h3->SetBinContent(i, y3);
326  h3->SetBinError (i, w3);
327  }
328  }
329 
330  } else if (opera == JOpera::Add() ||
331  opera == JOpera::Multiply()) {
332 
333  // operation between more than two histograms
334 
335  for (vector<TObject*>::iterator i = listOfObjects.begin(); i != listOfObjects.end(); ++i) {
336 
337  TH1* h1 = dynamic_cast<TH1*>(*i);
338 
339  NOTICE(h1->GetName() << ' ' << opera << endl);
340 
341  if (h3 == NULL) {
342 
343  if (name == JOpera::SAME_AS_OPERATION())
344  h3 = (TH1*) h1->Clone(opera.c_str());
345  else if (name == JOpera::SAME_AS_INPUT())
346  h3 = (TH1*) h1->Clone(h1->GetName());
347  else
348  h3 = (TH1*) h1->Clone(name.c_str());
349 
350  } else {
351 
352  if (opera == JOpera::Add()) {
353 
354  h3->Add (h3, h1, +1.0, +1.0);
355 
356  } else if (opera == JOpera::Multiply()) {
357 
358  h3->Multiply(h3, h1, +1.0, +1.0);
359  }
360  }
361  }
362 
363  } else {
364 
365  FATAL("Incompatible number of histograms " << listOfObjects.size() << " with option " << opera << endl);
366  }
367 
368  if (h3 != NULL) {
369 
370  TFile out(outputFile.c_str(), "recreate");
371 
372  h3->Write();
373 
374  out.Write();
375  out.Close();
376  }
377 }
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
JTOOLS::j
int j
Definition: JPolint.hh:703
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