Jpp  master_rocky-40-g5f0272dcd
the software that should make you happy
JCheckK40.cc
Go to the documentation of this file.
1 #include <string>
2 #include <iostream>
3 #include <iomanip>
4 #include <limits>
5 #include <vector>
6 #include <map>
7 
8 #include "TROOT.h"
9 #include "TFile.h"
10 #include "TH1D.h"
11 #include "TH2D.h"
12 #include "TF1.h"
13 #include "TFitResult.h"
14 
15 #include "JROOT/JMinimizer.hh"
16 
18 #include "JDAQ/JDAQHeaderIO.hh"
19 #include "JROOT/JRootFileReader.hh"
20 
21 #include "JDB/JDB.hh"
22 #include "JDB/JSelector.hh"
24 #include "JDB/JDBToolkit.hh"
25 #include "JDB/JPMTHVRunSettings.hh"
26 #include "JDB/JLocation_t.hh"
27 #include "JDB/JSupport.hh"
28 
29 #include "JDetector/JDetector.hh"
31 
33 #include "JROOT/JManager.hh"
34 #include "JTools/JRange.hh"
35 #include "JTools/JQuantile.hh"
36 
37 #include "Jeep/JProperties.hh"
38 #include "Jeep/JPrint.hh"
39 #include "Jeep/JParser.hh"
40 #include "Jeep/JMessage.hh"
41 
42 
43 namespace {
44 
45  using namespace std;
46  using namespace JPP;
47 
48  bool usePMTID; // option for old data for which correction for PMT cable swaps does not apply
49 
50  /**
51  * Setup corresponding to an output file from JCalibrateK40.
52  */
53  struct JSetup {
54 
55  /**
56  * Constructor.
57  *
58  * \param file_name file name
59  */
60  JSetup(const string& file_name) :
61  file_name(file_name)
62  {
63  fp = TFile::Open(file_name.c_str(), "read");
64 
65  int counter = 0;
66 
67  for (JRootFileReader<JDAQHeader> in(file_name.c_str()); in.hasNext(); ) {
68 
69  const JDAQHeader* p = in.next();
70 
71  if (counter == 0)
72  header = *p;
73  else
74  THROW(JException, "Multiple headers in file " << file_name);
75  }
76 
77  ResultSet& rs = getResultSet(getTable <JPMTHVRunSettings>(),
78  getSelector<JPMTHVRunSettings>(getDetector(header.getDetectorID()), header.getRunNumber()));
79 
80  for (JPMTHVRunSettings parameters; rs >> parameters; ) {
81 
82  JLocation_t location;
83 
84  if (usePMTID)
85  location = JLocation_t(parameters.DUID, parameters.FLOORID, parameters.PMTINTID);
86  else
87  location = JLocation_t(parameters.DUID, parameters.FLOORID, parameters.CABLEPOS);
88 
89  HV[location] = parameters.HV_VALUE;
90  }
91 
92  rs.Close();
93  }
94 
95  /**
96  * Get histogram from file.
97  *
98  * \param key histogram name
99  * \return pointer to histogram
100  */
101  TH2D* get(const string& key) const
102  {
103  return (TH2D*) fp->Get(key.c_str());
104  }
105 
106  string file_name;
107  JDAQHeader header;
109 
110  private:
111  TFile* fp;
112  };
113 }
114 
115 
116 /**
117  * \file
118  *
119  * Auxiliary program to check t0's.
120  * \author mdejong
121  */
122 int main(int argc, char **argv)
123 {
124  using namespace std;
125  using namespace JPP;
126  using namespace KM3NETDAQ;
127 
128  typedef JRange<size_t> JRange_t;
129 
130  JServer server;
131  string usr;
132  string pwd;
133  string cookie;
134  string detectorFile;
135  pair<string, string> inputFile;
136  string outputFile;
137  double precision;
138  double Wmin = 100.0;
139  JRange_t X;
140  string option;
141  int debug;
142 
143  try {
144 
145  JProperties properties;
146 
147  properties.insert(gmake_property(Wmin));
148 
149  JParser<> zap("Auxiliary program to check t0's.");
150 
151  zap['s'] = make_field(server) = getServernames();
152  zap['u'] = make_field(usr) = "";
153  zap['!'] = make_field(pwd) = "";
154  zap['C'] = make_field(cookie) = "";
155  zap['a'] = make_field(detectorFile);
156  zap['f'] = make_field(inputFile, "pair of input files (output of JCalibrateK40)");
157  zap['o'] = make_field(outputFile) = "k40.root";
158  zap['e'] = make_field(precision, "precision for HV comparison") = 0.5;
159  zap['@'] = make_field(properties) = JPARSER::initialised();
160  zap['x'] = make_field(X, "ROOT fit range (PMT pairs).") = JRange_t(300, numeric_limits<size_t>::max());
161  zap['O'] = make_field(option, "ROOT fit option, see TH1::Fit.") = "";
162  zap['U'] = make_field(usePMTID);
163  zap['d'] = make_field(debug) = 2;
164 
165  zap(argc, argv);
166  }
167  catch(const exception &error) {
168  FATAL(error.what() << endl);
169  }
170 
171 
172  try {
173  JDB::reset(usr, pwd, cookie);
174  }
175  catch(const exception& error) {
176  FATAL(error.what() << endl);
177  }
178 
179  JSetup setups[] = {
180  JSetup(inputFile.first),
181  JSetup(inputFile.second)
182  };
183 
184  for (int i = 0; i != 2; ++i) {
185  DEBUG(setw(32) << setups[i].file_name << ' ' << setups[i].header.getDetectorID() << ' ' << setups[i].header.getRunNumber() << endl);
186  }
187 
189 
190  try {
191  load(detectorFile, detector);
192  }
193  catch(const JException& error) {
194  FATAL(error);
195  }
196 
197  if (option.find('S') == string::npos) { option += 'S'; }
198  if (option.find('Q') == string::npos && debug < JEEP::debug_t) { option += 'Q'; }
199 
200 
201  JManager<int, TH1D> H1(new TH1D("[%]", NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS - 0.5));
202 
203 
204  TF1 f1("f1", "[0]*TMath::Gaus(x,[1],[2]) + [3]");
205 
206  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
207 
208  TH2D* h2[] = {
209  setups[0].get(MAKE_CSTRING(module->getID() << _2S)),
210  setups[1].get(MAKE_CSTRING(module->getID() << _2S))
211  };
212 
213  DEBUG("Module " << setw(10) << module->getID() << ' ' << (h2[0] != NULL) << (h2[0] != NULL) << endl);
214 
215  if (h2[0] == NULL ||
216  h2[1] == NULL) {
217  continue;
218  }
219 
220  JCombinatorics combinatorics;
221 
222  combinatorics.configure(module->size());
223 
224  combinatorics.sort(JPairwiseComparator(*module));
225 
226  vector<JQuantile> Q(module->size());
227 
228  for (size_t ip = max(X.getLowerLimit(), size_t(0)); ip != min(X.getUpperLimit(), combinatorics.getNumberOfPairs()); ++ip) {
229 
230  const JCombinatorics::pair_type pair = combinatorics.getPair(ip);
231 
232  const JLocation_t location_1(module->getString(), module->getFloor(), pair.first);
233  const JLocation_t location_2(module->getString(), module->getFloor(), pair.second);
234 
235  const bool hv_1 = (fabs(setups[0].HV[location_1] - setups[1].HV[location_1]) < precision);
236  const bool hv_2 = (fabs(setups[0].HV[location_2] - setups[1].HV[location_2]) < precision);
237 
238  double t1[] = {
239  numeric_limits<double>::max(),
240  numeric_limits<double>::max()
241  };
242 
243  const Int_t ix = ip + 1;
244 
245  for (int i = 0; i != 2; ++i) {
246 
247  TH1D h1("__py", NULL, h2[i]->GetYaxis()->GetNbins(), h2[i]->GetYaxis()->GetXmin(), h2[i]->GetYaxis()->GetXmax());
248 
249  // start values
250 
251  Double_t ymin = numeric_limits<double>::max();
252  Double_t ymax = numeric_limits<double>::lowest();
253  Double_t mean = 0.0;
254  Double_t sigma = 4.5;
255  Double_t W = 0.0;
256 
257  for (int iy = 1; iy <= h1.GetNbinsX(); ++iy) {
258 
259  const Double_t x = h1.GetBinCenter(iy);
260  const Double_t y = h2[i]->GetBinContent(ix,iy);
261 
262  h1.SetBinContent(iy, y);
263  h1.SetBinError (iy, sqrt(y));
264 
265  if (y > ymax) {
266  mean = x;
267  ymax = y;
268  }
269 
270  if (y < ymin) {
271  ymin = y;
272  }
273 
274  W += y;
275  }
276 
277  if (W >= Wmin) {
278 
279  f1.SetParameter(0, ymax);
280  f1.SetParameter(1, mean);
281  f1.SetParameter(2, sigma);
282  f1.SetParameter(3, ymin);
283 
284  for (Int_t i = 0; i != f1.GetNpar(); ++i) {
285  f1.SetParError(i, 0.0);
286  }
287 
288  TFitResultPtr result = h1.Fit(&f1, option.c_str(), "same");
289 
290  if (result.Get() == NULL) {
291  FATAL("Invalid TFitResultPtr " << h1.GetName() << endl);
292  }
293 
294  if (debug >= debug_t || !result->IsValid()) {
295  cout << "Histogram slice: "
296  << setw(3) << ix << ' '
297  << FIXED(7,3) << f1.GetParameter(1) << " +/- "
298  << FIXED(7,3) << f1.GetParError(1) << ' '
299  << FIXED(7,3) << result->Chi2() << '/'
300  << result->Ndf() << ' '
301  << (result->IsValid() ? "" : "failed") << endl;
302  }
303 
304  t1[i] = f1.GetParameter(1);
305  }
306  }
307 
308  if (t1[0] != numeric_limits<double>::max() &&
309  t1[1] != numeric_limits<double>::max()) {
310 
311  if (hv_1 != hv_2) {
312 
314 
315  if (hv_1) {
316  p2 = JCombinatorics::pair_type(pair.second, pair.first);
317  }
318 
319  if (hv_2) {
320  p2 = JCombinatorics::pair_type(pair.first, pair.second);
321  }
322 
323  if (debug >= debug_t) {
324  cout << setw(10) << module->getID() << "." << FILL(2,'0') << p2.first << FILL() << ' ';
325  cout << "(" << FILL(2,'0') << p2.second << FILL() << ")" << ' ';
326  cout << FIXED(6,2) << (combinatorics.getSign(p2) * (t1[1] - t1[0])) << endl;
327  }
328 
329  Q[p2.first].put(combinatorics.getSign(p2) * (t1[1] - t1[0]));
330  }
331  }
332  }
333 
334  for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
335  if (Q[i].getCount() >= 2) {
336  H1[module->getID()]->SetBinContent(i+1, Q[i].getMean());
337  H1[module->getID()]->SetBinError (i+1, Q[i].getSTDev());
338  }
339  }
340  }
341 
342  if (outputFile != "") {
343  H1.Write(outputFile.c_str());
344  }
345 }
double getMean(vector< double > &v)
get mean of vector content
int main(int argc, char **argv)
Definition: JCheckK40.cc:122
string outputFile
KM3NeT DAQ constants, bit handling, etc.
ROOT TTree parameter settings.
Data structure for detector geometry and calibration.
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:712
Dynamic ROOT object management.
General purpose messaging.
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
#define FATAL(A)
Definition: JMessage.hh:67
int debug
debug level
Definition: JSirene.cc:69
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:2142
I/O formatting auxiliaries.
#define MAKE_CSTRING(A)
Make C-string.
Definition: JPrint.hh:72
Utility class to parse parameter values.
#define gmake_property(A)
macros to convert (template) parameter to JPropertiesElement object
Auxiliary class to define a range between two values.
Detector data structure.
Definition: JDetector.hh:96
Utility class to parse parameter values.
Definition: JProperties.hh:501
General exception.
Definition: JException.hh:24
Utility class to parse command line options.
Definition: JParser.hh:1698
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys.
Definition: JManager.hh:47
void Write(TDirectory &out, const bool wm=false)
Write objects to file.
Definition: JManager.hh:304
ROOT file reader.
Auxiliary class to convert pair of indices to unique index and back.
void configure(const int numberOfIndices)
Configure.
static int getSign(const int first, const int second)
Sign of pair of indices.
void sort(JComparator_t comparator)
Sort address pairs.
size_t getNumberOfPairs() const
Get number of pairs.
const pair_type & getPair(const int index) const
Get pair of indices for given index.
Range of values.
Definition: JRange.hh:42
T getLowerLimit() const
Get lower limit.
Definition: JRange.hh:202
T getUpperLimit() const
Get upper limit.
Definition: JRange.hh:213
const double sigma[]
Definition: JQuadrature.cc:74
const JPolynome f1(1.0, 2.0, 3.0)
Function.
static const char *const _2S
Name extension for 2D counts.
JCombinatorics::pair_type pair_type
ResultSet & getResultSet(const std::string &query)
Get result set.
Definition: JDB.hh:437
JDetectorsHelper & getDetector()
Auxiliary function for helper object initialisation.
Definition: JDBToolkit.hh:378
std::vector< JServer > getServernames()
Get list of names of available database servers.
Definition: JDB.hh:107
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
@ debug_t
debug
Definition: JMessage.hh:29
size_t getCount(const array_type< T > &buffer, const JCompare_t &compare)
Count number of unique values.
Definition: JVectorize.hh:261
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
void reset(T &value)
Reset value.
KM3NeT DAQ data structures and auxiliaries.
Definition: DataQueue.cc:39
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
Definition: JDAQ.hh:26
Definition: JSTDTypes.hh:14
Auxiliary data structure for sequence of same character.
Definition: JManip.hh:330
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:448
Type definition of range.
Definition: JHead.hh:43
Detector file.
Definition: JHead.hh:227
Auxiliary data structure for setup of complete system.
Definition: JSydney.cc:108
Auxiliary class to sort pairs of PMT addresses within optical module.
Auxiliary data structure for location of product in detector.
Definition: JLocation_t.hh:26
Wrapper class for server name.
Definition: JDB.hh:53
Template definition for getting table specific selector.
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:68
Data structure for a pair of indices.