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