Jpp  18.2.1-ARCA-DF-PATCH
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JStdev1D.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 "TProfile.h"
13 #include "TString.h"
14 #include "TRegexp.h"
15 
16 #include "JGizmo/JRootObjectID.hh"
17 #include "JGizmo/JGizmoToolkit.hh"
18 
19 #include "Jeep/JParser.hh"
20 #include "Jeep/JMessage.hh"
21 
22 
23 /**
24  * \file
25  * Auxiliary program to determine standard deviation of a set of 1D histograms.
26  * The option <tt>-f</tt> corresponds to <tt><file name>:<object name></tt>.
27  * \author mdejong
28  */
29 int main(int argc, char **argv)
30 {
31  using namespace std;
32  using namespace JPP;
33 
34  vector<JRootObjectID> inputFile;
35  string outputFile;
36  int debug;
37 
38  try {
39 
40  JParser<> zap("Auxiliary program to determine standard deviation of a set of 1D histograms.");
41 
42  zap['f'] = make_field(inputFile, "<input file>:<object name>");
43  zap['o'] = make_field(outputFile, "ROOT file with histogram") = "stdev.root";
44  zap['d'] = make_field(debug) = 1;
45 
46  zap(argc, argv);
47  }
48  catch(const exception &error) {
49  FATAL(error.what() << endl);
50  }
51 
52 
53  TH1* h3 = NULL;
54 
55  vector<TH1*> listOfHistograms;
56 
57  for (vector<JRootObjectID>::const_iterator input = inputFile.begin(); input != inputFile.end(); ++input) {
58 
59  DEBUG("Input: " << *input << endl);
60 
61  TDirectory* dir = getDirectory(*input);
62 
63  if (dir == NULL) {
64  ERROR("File: " << input->getFullFilename() << " not opened." << endl);
65  continue;
66  }
67 
68  const TRegexp regexp(input->getObjectName());
69 
70  TIter iter(dir->GetListOfKeys());
71 
72  for (TKey* key; (key = (TKey*) iter.Next()) != NULL; ) {
73 
74  const TString tag(key->GetName());
75 
76  DEBUG("Key: " << tag << " match = " << tag.Contains(regexp) << endl);
77 
78  // option match
79 
80  if (tag.Contains(regexp) && isTObject(key)) {
81 
82  TObject* p = key->ReadObj();
83 
84  if (!p->InheritsFrom("TH1")) {
85  FATAL("Object " << p->GetName() << " not compatible with histogram operations." << endl);
86  }
87 
88  if (h3 == NULL) {
89 
90  h3 = (TH1*) p->Clone("stdev");
91 
92  } else {
93 
94  TH1* h1 = (TH1*) p;
95 
96  if (h3->GetXaxis()->GetNbins() != h1->GetXaxis()->GetNbins() ||
97  h3->GetXaxis()->GetXmin () != h1->GetXaxis()->GetXmin () ||
98  h3->GetXaxis()->GetXmax () != h1->GetXaxis()->GetXmax ()) {
99  FATAL("Objects " << h3->GetName() << " and " << p->GetName() << " have different binning." << endl);
100  }
101  }
102 
103  listOfHistograms.push_back((TH1*) p);
104  }
105  }
106  }
107 
108  if (listOfHistograms.size() <= 1) {
109  FATAL("Not enough histograms " << listOfHistograms.size() << endl);
110  }
111 
112  for (Int_t i = 1; i <= h3->GetNbinsX(); ++i) {
113 
114  Double_t Y = 0.0;
115  Double_t Z = 0.0;
116 
117  for (vector<TH1*>::const_iterator p = listOfHistograms.begin(); p != listOfHistograms.end(); ++p) {
118  Y += (*p)->GetBinContent(i);
119  Z += (*p)->GetBinError(i) * (*p)->GetBinError(i);
120  }
121 
122  Z = sqrt(Z);
123  Y /= listOfHistograms.size();
124  Z /= listOfHistograms.size();
125 
126  Double_t y = 0.0;
127 
128  for (vector<TH1*>::const_iterator p = listOfHistograms.begin(); p != listOfHistograms.end(); ++p) {
129  y += ((*p)->GetBinContent(i) - Y)*((*p)->GetBinContent(i) - Y) / ((*p)->GetBinError(i)*(*p)->GetBinError(i) + Z*Z);
130  }
131 
132  y = sqrt(y / listOfHistograms.size());
133 
134  h3->SetBinContent(i, y);
135  h3->SetBinError (i, 0.0);
136  }
137 
138  if (h3 != NULL) {
139 
140  TFile out(outputFile.c_str(), "recreate");
141 
142  h3->Write();
143 
144  out.Write();
145  out.Close();
146  }
147 }
Utility class to parse command line options.
Definition: JParser.hh:1514
int main(int argc, char *argv[])
Definition: Main.cc:15
Definition: JRoot.hh:19
then fatal Wrong number of arguments fi set_variable STRING $argv[1] set_variable DETECTORXY_TXT $WORKDIR $DETECTORXY_TXT tail read X Y CHI2 RMS printf optimum n $X $Y $CHI2 $RMS awk v Y
string outputFile
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1989
#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.
do set_variable MODULE getModule a $WORKDIR detector_a datx L $STRING JEditDetector a $WORKDIR detector_a datx M $MODULE setz o $WORKDIR detector_a datx JEditDetector a $WORKDIR detector_b datx M $MODULE setz o $WORKDIR detector_b datx done echo Output stored at $WORKDIR detector_a datx and $WORKDIR tripod_a txt JDrawDetector2D a $WORKDIR detector_a datx a $WORKDIR detector_b datx L BL o detector $FORMAT $BATCH JDrawDetector2D T $WORKDIR tripod_a txt T $WORKDIR tripod_b txt L BL o tripod $FORMAT $BATCH JCompareDetector a $WORKDIR detector_a datx b $WORKDIR detector_b datx o $WORKDIR abc root &dev null for KEY in X Y Z
int debug
debug level
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62