Jpp
Functions
JOpera2D.cc File Reference
#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>:

Definition in file JOpera2D.cc.

Function Documentation

◆ main()

int main ( int  argc,
char **  argv 
)

Definition at line 52 of file JOpera2D.cc.

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