Jpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions
JFitToT.cc File Reference

Auxiliary program to fit time-over-threshold distributions. More...

#include <string>
#include <iostream>
#include <iomanip>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TF1.h"
#include "TFitResult.h"
#include "TMath.h"
#include "JLang/JLangToolkit.hh"
#include "JPhysics/JConstants.hh"
#include "km3net-dataformat/online/JDAQ.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JPMTParametersMap.hh"
#include "JDetector/JPMTAnalogueSignalProcessor.hh"
#include "JROOT/JRootToolkit.hh"
#include "JTools/JRange.hh"
#include "JSupport/JMeta.hh"
#include "JCalibrate/JCalibrateToT.hh"
#include "JCalibrate/JFitToT.hh"
#include "Jeep/JPrint.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 to fit time-over-threshold distributions.

The fitted parameters are the PMT gain and gain spread.

Author
mkarel

Definition in file JFitToT.cc.

Function Documentation

int main ( int  argc,
char **  argv 
)

< Minimal histogram weight in fit range

< Fit range left of time-over-threshold peak

< Fit range right of time-over-threshold peak

< Minimal right-to-left gradient value in time-over-threshold peak-search, normalized with respect to the total number of histogram entries

Definition at line 50 of file JFitToT.cc.

51 {
52  using namespace std;
53  using namespace JPP;
54  using namespace KM3NETDAQ;
55 
56  typedef JRange<double> JRange_t;
57 
58  string inputFile;
59  string outputFile;
60  string detectorFile;
61  string pmtFile;
62  bool writeFits;
63 
64  double Wmin = 5e3; //!< Minimal histogram weight in fit range
65  double LeftMargin = 10.0; //!< Fit range left of time-over-threshold peak
66  double RightMargin = 10.0; //!< Fit range right of time-over-threshold peak
67  double gradientThreshold = -0.005; //!< Minimal right-to-left gradient value in time-over-threshold peak-search, normalized with respect to the total number of histogram entries
68 
69  JRange_t fitRange;
70  string option;
71 
72  string regexp;
73  int debug;
74 
75  try {
76 
77  JProperties properties;
78 
79  properties.insert(gmake_property(Wmin));
80  properties.insert(gmake_property(LeftMargin));
81  properties.insert(gmake_property(RightMargin));
82 
83  properties.insert(gmake_property(gradientThreshold));
84 
85 
86  JParser<> zap("Auxiliary program to fit time-over-threshold distributions.");
87 
88  zap['f'] = make_field(inputFile, "input file (output from JCalibrateToT).");
89  zap['o'] = make_field(outputFile, "output file.") = "fit.root";
90  zap['a'] = make_field(detectorFile, "detector file.");
91  zap['P'] = make_field(pmtFile, "specify PMT file name that can be given to JTriggerEfficiency.") = "";
92  zap['w'] = make_field(writeFits, "write fit results.");
93  zap['O'] = make_field(option, "ROOT fit options, see TH1::Fit") = "";
94  zap['@'] = make_field(properties) = JPARSER::initialised();
95  zap['x'] = make_field(fitRange, "ROOT fit range (time-over-threshold).") = JRange_t(0.0, 35.0);
96  zap['R'] = make_field(regexp, "regular expression for histogram name.") = MAKE_CSTRING(WILDCARD << _2SToT);
97  zap['d'] = make_field(debug, "debug.") = 1;
98 
99  zap(argc, argv);
100  }
101  catch(const exception &error) {
102  FATAL(error.what() << endl);
103  }
104 
105 
106  //----------------------------------------------------------
107  // load detector file and PMT parameters
108  //----------------------------------------------------------
109 
110  JDetector detector;
111 
112  try {
113  load(detectorFile, detector);
114  }
115  catch(const JException& error) {
116  FATAL(error);
117  }
118 
119 
120  JPMTParametersMap parameters;
121 
122  if (pmtFile != "") {
123  try {
124  parameters.load(pmtFile.c_str());
125  }
126  catch(const JException& error) {}
127  }
128 
129  parameters.comment.add(JMeta(argc, argv));
130 
131  // Set ROOT fit options
132  if (option.find('R') == string::npos) { option += 'R'; }
133  if (option.find('S') == string::npos) { option += 'S'; }
134  if (option.find('Q') == string::npos && debug < JEEP::debug_t) { option += 'Q'; }
135 
136 
137  TFile* in = TFile::Open(inputFile.c_str(), "exist");
138 
139  if (in == NULL || !in->IsOpen()) {
140  FATAL("File: " << inputFile << " not opened." << endl);
141  }
142 
143 
144  TH2D gain("gain", NULL,
145  100, 0.0, 2.0,
146  100, 0.0, 1.0);
147  TH1D chi2("chi2", NULL, 100, 0.0, 5.0);
148 
149 
150  TFile out(outputFile.c_str(), "recreate");
151 
152  //----------------------------------------------------------
153  // loop over modules in detector file to fit histogram
154  //----------------------------------------------------------
155 
156  for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
157 
158  NOTICE("Module " << setw(10) << module->getID() << ' ' << module->getLocation() << endl);
159 
160  if (module->empty()) {
161  continue;
162  }
163 
164  //----------------------------------------------------------
165  // get histogram for this module
166  //----------------------------------------------------------
167 
168  TH2D* h2s = (TH2D*) in->Get(replace(regexp, WILDCARD, MAKE_STRING(module->getID())).c_str());
169 
170  if (h2s == NULL) {
171 
172  WARNING("No histogram " << module->getID() << "; skip fit." << endl);
173 
174  continue;
175  }
176 
177  const double ny = h2s->GetYaxis()->GetNbins();
178  const double ymin = h2s->GetYaxis()->GetXmin();
179  const double ymax = h2s->GetYaxis()->GetXmax();
180 
181  TH1D chi2_1d (MAKE_CSTRING(module->getID() << ".1chi2"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS-0.5);
182  TH1D gain_1d (MAKE_CSTRING(module->getID() << ".1gain"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS-0.5);
183  TH1D gainspread_1d(MAKE_CSTRING(module->getID() << ".1gainspread"), NULL, NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS-0.5);
184 
185  for (int ix = 1; ix <= h2s->GetXaxis()->GetNbins(); ++ix) {
186 
187  const JPMTIdentifier id(module->getID(), ix-1);
188 
189  parameters[id].setPMTParameters(parameters.getDefaultPMTParameters());
190 
191  const string name = MAKE_STRING(module->getID() << '.' << FILL(2,'0') << id.getPMTAddress() << FITTOT_SUFFIX);
192  TH1D h1(name.c_str(), NULL, ny, ymin, ymax);
193 
194  h1.Sumw2();
195 
196  for (int iy = h2s->GetNbinsY(); iy >= 1; --iy) {
197 
198  const double w = h1.GetBinWidth (iy);
199  const double y = h2s->GetBinContent(ix, iy);
200 
201  h1.SetBinContent(iy, y/w);
202  h1.SetBinError (iy, sqrt(y)/w);
203  }
204 
205  const double weight = h1.Integral(h1.FindBin(fitRange.getLowerLimit()),
206  h1.FindBin(fitRange.getUpperLimit()));
207 
208  if (weight <= Wmin) {
209 
210  WARNING("Insufficient histogram entries for PMT " << id <<
211  "(" << weight << " < " << Wmin << "); skip fit." << endl);
212 
213  if (writeFits) {
214  out << h1;
215  }
216 
217  continue;
218  }
219 
220  //----------------------------------------------------------
221  // find the mode of the time-over-threshold distribution
222  //----------------------------------------------------------
223 
224  double ToTmode = JDETECTOR::TIME_OVER_THRESHOLD_NS;
225  double max = 0.0;
226 
227  for (int i = h1.GetBinCenter(ny-1);
228  i > h1.GetBinCenter(h1.FindBin(parameters[id].mean_ns + parameters[id].sigma_ns));
229  --i) {
230 
231  const double x = h1.GetBinCenter (i);
232  const double y = h1.GetBinContent(i);
233 
234  const double gradient = ( (h1.GetBinContent(i-1) - h1.GetBinContent(i+1)) /
235  (h1.GetBinCenter (i+1) - h1.GetBinCenter (i-1)) );
236 
237  if (y > max) {
238 
239  ToTmode = x;
240  max = y;
241 
242  } else if (gradient < -fabs(gradientThreshold) * h1.Integral()) {
243 
244  break;
245  }
246  }
247 
248  //----------------------------------------------------------
249  // set fit range
250  //----------------------------------------------------------
251 
252  if (!fitRange.is_valid()) {
253  fitRange.setRange(ToTmode - LeftMargin,
254  ToTmode + RightMargin);
255  }
256 
257  fitRange.setRange(std::max(h1.GetXaxis()->GetXmin(), fitRange.getLowerLimit()),
258  std::min(h1.GetXaxis()->GetXmax(), fitRange.getUpperLimit()));
259 
260  //----------------------------------------------------------
261  // perform minimum chi-square fit
262  //----------------------------------------------------------
263 
264  JFitToT fit(parameters[id], fitRange);
265 
266  const double initGain = fit.getCPU().getNPE(ToTmode);
267  const double initGainSpread = initGain / 3.0;
268 
269  if (initGain > FITTOT_GAIN_MAX || initGain < FITTOT_GAIN_MIN) { // Invalid initial gain value
270 
271  WARNING("Invalid initial gain for PMT " << id << "; set default values." << endl);
272 
273  parameters[id].gain = (initGain > FITTOT_GAIN_MAX ? FITTOT_GAIN_MAX : FITTOT_GAIN_MIN);
274  parameters[id].gainSpread = (initGain > FITTOT_GAIN_MAX ? FITTOT_GAINSPREAD_MAX : FITTOT_GAINSPREAD_MIN);
275 
276  if (writeFits) {
277  out << h1;
278  }
279 
280  continue;
281  }
282 
283  fit.gain = initGain;
284  fit.gainSpread = initGainSpread;
285 
286  const TFitResultPtr result = fit(h1, option);
287 
288  if (int(result) < 0 || !result->IsValid()) { // Fit failed
289 
290  WARNING("Fit failure for PMT " << id << "; set initialization values." << endl);
291 
292  h1.GetListOfFunctions()->Clear();
293 
294  parameters[id].gain = initGain;
295  parameters[id].gainSpread = initGainSpread;
296 
297  if (writeFits) {
298  out << h1;
299  }
300 
301  continue;
302  }
303 
304  if ((ymin < fitRange.getLowerLimit() || ymax > fitRange.getUpperLimit()) &&
305  (option.find("0") == string::npos && option.find("N") == string::npos)) {
306 
307  // add function with full range
308 
309  TF1* f1 = new TF1(MAKE_CSTRING(FITTOT_FNAME << "_fullRange"),
310  &fit,
311  &JFitToT::getValue,
312  ymin,
313  ymax,
314  JFitToT::getNumberOfModelParameters());
315  f1->SetParameters(fit.GetParameters());
316  f1->SetLineStyle(kDotted);
317  f1->SetNpx(1000);
318 
319  h1.GetListOfFunctions()->Add(f1);
320  }
321 
322  // store fit results
323 
324  const double reducedChi2 = result->Chi2() / result->Ndf();
325 
326  const JFitToTParameters values(fit.GetParameters());
327  const JFitToTParameters errors(fit.GetParErrors());
328 
329  chi2_1d. SetBinContent(ix, reducedChi2);
330  gain_1d. SetBinContent(ix, values.gain);
331  gain_1d. SetBinError (ix, errors.gain);
332  gainspread_1d.SetBinContent(ix, values.gainSpread);
333  gainspread_1d.SetBinError (ix, errors.gainSpread);
334 
335  gain.Fill(fit.gain, fit.gainSpread);
336  chi2.Fill(TMath::Log10(reducedChi2));
337 
338  NOTICE("PMT " << id << " chi2 / NDF " << result->Chi2() << ' ' << result->Ndf() << endl);
339 
340  parameters[id].gain = fit.gain;
341  parameters[id].gainSpread = fit.gainSpread;
342 
343  if (writeFits) {
344  out << h1;
345  }
346  }
347 
348  if (writeFits) {
349  out << chi2_1d << gain_1d << gainspread_1d;
350  }
351  }
352 
353  if (pmtFile != "") {
354  parameters.store(pmtFile.c_str());
355  }
356 
357  out << gain << chi2;
358 
359  out.Write();
360  out.Close();
361 }
JSUPPORT::JMeta
Auxiliary class for ROOT I/O of application specific meta data.
Definition: JMeta.hh:70
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1500
JLANG::JException
General exception.
Definition: JException.hh:23
JTOOLS::w
data_type w[N+1][M+1]
Definition: JPolint.hh:741
WARNING
#define WARNING(A)
Definition: JMessage.hh:65
JEEP::getValue
double getValue(const JScale_t scale)
Get numerical value corresponding to scale.
Definition: JScale.hh:47
JEEP::debug_t
debug
Definition: JMessage.hh:29
replace
TString replace(const TString &target, const TRegexp &regexp, const T &replacement)
Replace regular expression in input by given replacement.
Definition: JPrintResult.cc:63
gmake_property
#define gmake_property(A)
macro to convert (template) parameter to JPropertiesElement object
Definition: JProperties.hh:1291
JCALIBRATE::FITTOT_GAIN_MAX
static const double FITTOT_GAIN_MAX
Default maximal gain.
Definition: JFitToT.hh:44
JDETECTOR::JDetector
Detector data structure.
Definition: JDetector.hh:80
JDETECTOR::TIME_OVER_THRESHOLD_NS
const double TIME_OVER_THRESHOLD_NS
Specification for time-over-threshold corresponding to a one photo-electron pulse.
Definition: JDetector/JCalibration.hh:30
JCALIBRATE::FITTOT_GAIN_MIN
static const double FITTOT_GAIN_MIN
Default minimal gain.
Definition: JFitToT.hh:43
JEEP::JProperties
Utility class to parse parameter values.
Definition: JProperties.hh:496
parameters
*fatal Wrong number of arguments esac JCookie sh typeset Z DETECTOR typeset Z SOURCE_RUN typeset Z TARGET_RUN set_variable PARAMETERS_FILE $WORKDIR parameters
Definition: diff-Tuna.sh:38
MAKE_CSTRING
#define MAKE_CSTRING(A)
Make C-string.
Definition: JPrint.hh:151
h1
then for HISTOGRAM in h0 h1
Definition: JMatrixNZ.sh:69
JPARSER::initialised
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:66
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
JCALIBRATE::JFitToT
Parametrisation of time-over-threshold distribution.
Definition: JFitToT.hh:179
JCALIBRATE::JFitToTParameters
Fit parameters for two-fold coincidence rate due to K40.
Definition: JFitToT.hh:53
MAKE_STRING
#define MAKE_STRING(A)
Make string.
Definition: JPrint.hh:142
JCALIBRATE::FITTOT_SUFFIX
static const std::string FITTOT_SUFFIX
Definition: JFitToT.hh:36
JAANET::detector
Detector file.
Definition: JHead.hh:196
JCALIBRATE::FITTOT_GAINSPREAD_MAX
static const double FITTOT_GAINSPREAD_MAX
Default maximal gain spread.
Definition: JFitToT.hh:47
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1961
JTOOLS::result
return result
Definition: JPolint.hh:727
NOTICE
#define NOTICE(A)
Definition: JMessage.hh:64
JDETECTOR::JPMTParametersMap
Auxiliary class for map of PMT parameters.
Definition: JPMTParametersMap.hh:89
debug
int debug
debug level
Definition: JSirene.cc:63
FILL
Auxiliary data structure for sequence of same character.
Definition: JManip.hh:327
FATAL
#define FATAL(A)
Definition: JMessage.hh:67
name
then echo n User name
Definition: JCookie.sh:45
gain
then $JPP_DIR examples JDetector JToT o $OUTPUT_FILE n N $NPE P gain
Definition: JToT.sh:45
JDETECTOR::load
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Definition: JDetectorToolkit.hh:445
JCALIBRATE::_2SToT
static const char *const _2SToT
Histogram naming.
Definition: JCalibrateToT.hh:16
JCALIBRATE::FITTOT_GAINSPREAD_MIN
static const double FITTOT_GAINSPREAD_MIN
Default minimal gain spread.
Definition: JFitToT.hh:46
values
possible values
Definition: detector-XY:run.sh:41
detector
do set_variable DETECTOR_TXT $WORKDIR detector
Definition: detector-XY:fit1d.sh:37
KM3NETDAQ::NUMBER_OF_PMTS
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
Definition: JDAQ.hh:26
in
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable INPUT_FILE $argv[2] eval JPrintDetector a $DETECTOR O IDENTIFIER eval JPrintDetector a $DETECTOR O SUMMARY source JAcoustics sh $DETECTOR_ID typeset A TRIPODS get_tripods $WORKDIR tripod txt TRIPODS for EMITTER in
Definition: JCanberra.sh:36
JCALIBRATE::FITTOT_FNAME
static const std::string FITTOT_FNAME
Definition: JFitToT.hh:37
JTRIGGER::weight
std::vector< double > weight
Definition: JAlgorithm.hh:428
Generated by   doxygen 1.8.5