Jpp
JFitToT.cc
Go to the documentation of this file.
1 #include <string>
2 #include <iostream>
3 #include <iomanip>
4 
5 #include "TROOT.h"
6 #include "TFile.h"
7 #include "TH1D.h"
8 #include "TH2D.h"
9 #include "TF1.h"
10 #include "TFitResult.h"
11 #include "TMath.h"
12 
13 #include "JLang/JLangToolkit.hh"
14 #include "JTools/JConstants.hh"
15 #include "JDAQ/JDAQ.hh"
16 #include "JDetector/JDetector.hh"
17 #include "JDetector/JDetectorToolkit.hh"
18 #include "JDetector/JPMTParametersMap.hh"
19 
20 #include "JROOT/JRootToolkit.hh"
21 #include "JTools/JRange.hh"
22 #include "JSupport/JMeta.hh"
23 
24 #include "JCalibrate/JCalibrateToT.hh"
25 #include "JCalibrate/JFitToT.hh"
26 
27 #include "Jeep/JPrint.hh"
28 #include "Jeep/JParser.hh"
29 #include "Jeep/JMessage.hh"
30 
31 namespace {
32 
33  /**
34  * Wild card for histogram naming.\n
35  * The wild card will replaced by the module identifier.
36  */
37  static const char* const WILDCARD = "%";
38 }
39 
40 
41 /**
42  * \file
43  *
44  * Auxiliary program to fit time-over-threshold distributions.
45  * The fitted parameters are the PMT gain and gain spread.
46  * \author mkarel
47  */
48 int main(int argc, char **argv)
49 {
50  using namespace std;
51  using namespace JPP;
52  using namespace KM3NETDAQ;
53 
54  typedef JRange<double> JRange_t;
55 
56  string inputFile;
57  string outputFile;
58  string detectorFile;
59  string pmtFile;
60  bool writeFits;
61  JRange_t X;
62  string option;
63  double WMin;
64  string regexp;
65  int debug;
66 
67  try {
68 
69  JParser<> zap("Auxiliary program to fit time-over-threshold distributions.");
70 
71  zap['f'] = make_field(inputFile, "input file (output from JCalibrateToT).");
72  zap['o'] = make_field(outputFile, "output file.") = "fit.root";
73  zap['a'] = make_field(detectorFile, "detector file.");
74  zap['P'] = make_field(pmtFile, "specify PMT file name that can be given to JTriggerEfficiency.") = "";
75  zap['w'] = make_field(writeFits, "write fit results.");
76  zap['x'] = make_field(X, "ROOT fit range (time-over-threshold).") = JRange_t(10.0, 35.0);
77  zap['O'] = make_field(option, "ROOT fit option, see TH1::Fit.") = "";
78  zap['W'] = make_field(WMin, "minimal histogram contents.") = 1000.0;
79  zap['R'] = make_field(regexp, "regular expression for histogram name.") = MAKE_CSTRING(WILDCARD << _2SToT);
80  zap['d'] = make_field(debug, "debug.") = 1;
81 
82  zap(argc, argv);
83  }
84  catch(const exception &error) {
85  FATAL(error.what() << endl);
86  }
87 
88 
89  //----------------------------------------------------------
90  // load detector file and PMT parameters
91  //----------------------------------------------------------
92 
93  JDetector detector;
94 
95  try {
96  load(detectorFile, detector);
97  }
98  catch(const JException& error) {
99  FATAL(error);
100  }
101 
102 
103  JPMTParametersMap parameters;
104 
105  if (pmtFile != "") {
106  try {
107  parameters.load(pmtFile.c_str());
108  }
109  catch(const JException& error) {}
110  }
111 
112  parameters.comment.add(JMeta(argc, argv));
113 
114 
115  if (option.find('R') == string::npos) { option += 'R'; }
116  if (option.find('S') == string::npos) { option += 'S'; }
117  if (option.find('Q') == string::npos && debug < JEEP::debug_t) { option += 'Q'; }
118 
119 
120  TFile* in = TFile::Open(inputFile.c_str(), "exist");
121 
122  if (in == NULL || !in->IsOpen()) {
123  FATAL("File: " << inputFile << " not opened." << endl);
124  }
125 
126 
127  TH2D gain("gain", NULL,
128  100, 0.0, 2.0,
129  100, 0.0, 1.0);
130  TH1D chi2("chi2", NULL, 100, 0.0, 5.0);
131 
132 
133  TFile out(outputFile.c_str(), "recreate");
134 
135  //----------------------------------------------------------
136  // loop over modules in detector file to fit histogram
137  //----------------------------------------------------------
138 
139  for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
140 
141  DEBUG("Module " << module->getID() << endl);
142 
143  //----------------------------------------------------------
144  // get histogram for this module
145  //----------------------------------------------------------
146 
147  TH2D* h2s = (TH2D*) in->Get(replace(regexp, WILDCARD, MAKE_STRING(module->getID())).c_str());
148 
149  if (h2s == NULL) {
150 
151  WARNING("No histogram " << module->getID() << "; skip fit." << endl);
152 
153  continue;
154  }
155 
156  const int ny = h2s->GetYaxis()->GetNbins();
157  const double ymin = h2s->GetYaxis()->GetXmin();
158  const double ymax = h2s->GetYaxis()->GetXmax();
159 
160  TH1D chi2_1d (MAKE_CSTRING(module->getID() << ".1chi2"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS-0.5);
161  TH1D gain_1d (MAKE_CSTRING(module->getID() << ".1gain"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS-0.5);
162  TH1D gainspread_1d(MAKE_CSTRING(module->getID() << ".1gainspread"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS-0.5);
163 
164  for (int ix = 1; ix <= h2s->GetXaxis()->GetNbins(); ++ix) {
165 
166  const JPMTIdentifier id(module->getID(), ix-1);
167 
168  TH1D h1(MAKE_CSTRING(module->getID() << "." << id.getPMTAddress() << ".1ToT"), NULL, ny, ymin, ymax);
169 
170  h1.Sumw2();
171 
172  Double_t W = 0.0;
173 
174  for (int iy = 1; iy <= h2s->GetNbinsY(); ++iy) {
175 
176  const double x = h1.GetBinCenter(iy);
177  const double z = h1.GetBinWidth (iy);
178  const double y = h2s->GetBinContent(ix, iy);
179 
180  if (X(x)) {
181  W += y;
182  }
183 
184  h1.SetBinContent(iy, y/z);
185  h1.SetBinError (iy, sqrt(y)/z);
186  }
187 
188  if (W <= WMin) {
189 
190  WARNING("Not enough entries for PMT " << id << ' ' << W << "; skip fit." << endl);
191 
192  continue;
193  }
194 
195  h1.Scale(1.0/h1.Integral());
196 
197 
198  //----------------------------------------------------------
199  // set-up fitting procedure
200  //----------------------------------------------------------
201 
202  JFitToT fit(parameters.getPMTParameters(id),
203  max(ymin, X.first),
204  min(ymax, X.second));
205 
206 
207  const TFitResultPtr result = fit(h1, option.c_str());
208 
209  chi2.Fill(TMath::Log10(result->Chi2() / result->Ndf()));
210 
211  if (debug >= JEEP::notice_t) {
212  cout << "PMT " << id << " chi2 / NDF " << result->Chi2() << ' ' << result->Ndf() << ' ' << (result->IsValid() ? "" : "failed") << endl;
213  }
214 
215  if (!result->IsValid()) {
216 
217  WARNING("Fit failure for PMT " << id << "; skip." << endl);
218 
219  continue;
220  }
221 
222  if (ymin < X.getLowerLimit() || ymax > X.getUpperLimit()) {
223 
224  // add function with full range
225 
226  TF1* p = new TF1("f1&",
227  JFitToT::getValue,
228  ymin,
229  ymax,
230  JFitToT::getNumberOfModelParameters());
231 
232  p->SetParameters(fit.GetParameters());
233  p->SetNpx(1000);
234 
235  h1.GetListOfFunctions()->Add(p);
236  }
237 
238 
239  // store fit results
240 
241  const JFitToTParameters values(fit.GetParameters());
242  const JFitToTParameters errors(fit.GetParErrors());
243 
244  chi2_1d. SetBinContent(ix, result->Chi2() / result->Ndf());
245  gain_1d. SetBinContent(ix, values.gain);
246  gain_1d. SetBinError (ix, errors.gain);
247  gainspread_1d.SetBinContent(ix, values.gainSpread);
248  gainspread_1d.SetBinError (ix, errors.gainSpread);
249 
250  gain.Fill(fit.gain, fit.gainSpread);
251 
252  if (fit.gain > 0.0) {
253 
254  parameters[id].gain = fit.gain;
255  parameters[id].gainSpread = fit.gainSpread;
256 
257  } else {
258 
259  WARNING("Error fit of PMT " << id << "; set default values." << endl);
260 
261  parameters[id].gain = 0.0;
262  parameters[id].gainSpread = parameters.getDefaultPMTParameters().gainSpread;
263  }
264 
265  if (writeFits) {
266  out << h1;
267  }
268  }
269 
270  if (writeFits) {
271  out << chi2_1d << gain_1d << gainspread_1d;
272  }
273  }
274 
275  if (pmtFile != "") {
276  parameters.store(pmtFile.c_str());
277  }
278 
279  out << gain << chi2;
280 
281  out.Write();
282  out.Close();
283 }
284 
285 
JPMTParametersMap.hh
JMeta.hh
JDAQ.hh
main
int main(int argc, char **argv)
Definition: JFitToT.cc:48
JMessage.hh
JCalibrateToT.hh
JPrint.hh
JEEP::notice_t
notice
Definition: JMessage.hh:32
JDETECTOR::load
void load(const JString &file_name, JDetector &detector)
Load detector from input file.
Definition: JDetectorToolkit.hh:456
KM3NETDAQ::NUMBER_OF_PMTS
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
Definition: JDAQ.hh:26
JCALIBRATE::_2SToT
static const char *const _2SToT
Histogram naming.
Definition: JCalibrateToT.hh:16
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
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition: JAAnetToolkit.hh:37
MAKE_CSTRING
#define MAKE_CSTRING(A)
Make C-string.
Definition: JPrint.hh:708
JRange.hh
WARNING
#define WARNING(A)
Definition: JMessage.hh:65
debug
int debug
debug level
Definition: JSirene.cc:59
JConstants.hh
JTOOLS::result
return result
Definition: JPolint.hh:695
JRootToolkit.hh
JLANG::replace
std::string replace(const std::string &input, const std::string &target, const std::string &replacement)
Replace tokens in string.
Definition: JLangToolkit.hh:124
JParser.hh
JDetectorToolkit.hh
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
KM3NETDAQ
KM3NeT DAQ data structures and auxiliaries.
Definition: DataQueue.cc:39
JEEP::getValue
double getValue(const JScale_t scale)
Get numerical value corresponding to scale.
Definition: JScale.hh:47
JDetector.hh
JFitToT.hh
FATAL
#define FATAL(A)
Definition: JMessage.hh:67
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
JEEP::debug_t
debug
Definition: JMessage.hh:29
JLangToolkit.hh
Generated by   doxygen 1.8.16