Jpp
JFrodo.cc
Go to the documentation of this file.
1 #include <string>
2 #include <iostream>
3 #include <iomanip>
4 #include <vector>
5 
6 #include "TROOT.h"
7 #include "TFile.h"
8 #include "TH1D.h"
9 #include "TH2D.h"
10 #include "TF1.h"
11 #include "TMath.h"
12 #include "TFitResult.h"
13 
14 #include "JTools/JRange.hh"
15 #include "JDetector/JCalibration.hh"
16 #include "JDetector/JPMTAnalogueSignalProcessor.hh"
17 
18 #include "Jeep/JPrint.hh"
19 #include "Jeep/JParser.hh"
20 #include "Jeep/JMessage.hh"
21 
22 namespace {
23 
24  /**
25  * Name of histogram with fit results from JGandalf.
26  */
27  const char* h2_t = "ha";
28 
29  const char* gauss_t = "Gauss"; //!< Gaussian
30  const char* emg_t = "EMG"; //!< Exponentially Modified Gaussian
31 
32 
33  using namespace JPP;
34 
35  JPMTAnalogueSignalProcessor cpu;
36 
37  /**
38  * Fit function.
39  *
40  * \param x abscissa values
41  * \param parameters parameter values
42  * \return ordinate value
43  */
44  Double_t getRisetime(Double_t* x, Double_t* parameters)
45  {
46  const double t0 = parameters[0];
47  const double p1 = parameters[1];
48  const double t1 = parameters[2];
49 
50  const double npe = cpu.getNPE(x[0]);
51 
52  double ts = t0 + cpu.getRiseTime(npe) + t1 * (1.0 - exp(-npe*p1));
53 
54  return ts;
55  }
56 }
57 
58 
59 /**
60  * \file
61  * Fit time-slewing histogram in output of JGandalf.cc in calibration mode.
62  * \author mdejong
63  */
64 int main(int argc, char **argv)
65 {
66  using namespace std;
67  using namespace JPP;
68 
69  typedef JRange<double> JRange_t;
70 
71  vector<string> inputFile;
72  string outputFile;
73  string function;
74  JRange_t T_ns;
75  bool writeFits;
76  JPMTParameters parameters;
77  JRange_t X;
78  string option;
79  int debug;
80 
81  try {
82 
83  JProperties properties = parameters.getProperties();
84 
85  JParser<> zap("Program to fit time-slewing histogram in output of JGandalf.cc in calibration mode.");
86 
87  zap['f'] = make_field(inputFile, "input files (output from JGandalf -c).");
88  zap['o'] = make_field(outputFile, "output file.") = "frodo.root";
89  zap['F'] = make_field(function, "fit function") = gauss_t, emg_t;
90  zap['T'] = make_field(T_ns, "ROOT fit range around maximum [ns].") = JRange_t(-7.5,+7.5);
91  zap['w'] = make_field(writeFits, "write fit results.");
92  zap['P'] = make_field(properties) = JPARSER::initialised();
93  zap['x'] = make_field(X, "ROOT fit range for time-over threshold") = JRange_t(0.0, 256.0);
94  zap['O'] = make_field(option, "ROOT fit option, see TH1::Fit.") = "";
95  zap['d'] = make_field(debug) = 1;
96 
97  zap(argc, argv);
98  }
99  catch(const exception &error) {
100  FATAL(error.what() << endl);
101  }
102 
103 
104  cout.tie(&cerr);
105 
106  if (!T_ns.is_valid()) {
107  FATAL("Invalid fit range [ns] " << T_ns << endl);
108  }
109 
110  cpu.setPMTParameters(parameters);
111 
112  if (option.find('S') == string::npos) { option += 'S'; }
113  if (option.find('Q') == string::npos && debug < JEEP::debug_t) { option += 'Q'; }
114 
115  // rebin time-over-threshold
116 
117  const Double_t x[] = {
118  -0.5, 3.5, 6.5, 9.5, 12.5, 15.5, 18.5,
119  20.5, 21.5, 22.5, 23.5, 24.5, 25.5, 26.5, 27.5, 28.5, 29.5,
120  30.5, 31.5, 32.5, 33.5, 34.5, 35.5, 36.5, 37.5, 38.5, 39.5,
121  40.5, 42.5, 44.5, 46.5, 48.5,
122  50.5, 52.5, 54.5, 56.5, 58.5,
123  60.5, 65.5,
124  70.5, 75.5,
125  80.5, 85.5,
126  90.5, 95.5,
127  100.5, 120.5, 140.5, 160.5, 180.5, 200.5, 250.5
128  };
129 
130  const int N = sizeof(x)/sizeof(x[0]);
131 
132  TH2D* h2 = NULL;
133 
134  for (vector<string>::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {
135 
136  NOTICE("Processing " << *i << endl);
137 
138  TFile in(i->c_str(), "exist");
139 
140  if (!in.IsOpen()) {
141  FATAL("File " << *i << " not opened." << endl);
142  }
143 
144  TH2D* p = dynamic_cast<TH2D*>(in.Get(h2_t));
145 
146  if (p == NULL) {
147  FATAL("File " << *i << " has no histogram <" << h2_t << ">." << endl);
148  }
149 
150  if (h2 == NULL)
151  h2 = (TH2D*) p->Clone();
152  else
153  h2->Add(p);
154 
155  h2->SetDirectory(0);
156 
157  in.Close();
158  }
159 
160  if (h2 == NULL) {
161  FATAL("No histogram <" << h2_t << ">." << endl);
162  }
163 
164  // histogram fit results
165 
166  TH1D h0("h0", NULL, N-1, x);
167 
168  TFile out(outputFile.c_str(), "recreate");
169 
170  for (int i = 1; i <= h0.GetXaxis()->GetNbins(); ++i) {
171 
172  const Int_t imin = h2->GetXaxis()->FindBin(h0.GetXaxis()->GetBinLowEdge(i));
173  const Int_t imax = h2->GetXaxis()->FindBin(h0.GetXaxis()->GetBinUpEdge (i));
174 
175  TH1D h1(MAKE_CSTRING(i << ".1D"), NULL, h2->GetYaxis()->GetNbins(), h2->GetYaxis()->GetXmin(), h2->GetYaxis()->GetXmax());
176 
177  // start values
178 
179  Double_t ymax = 0.0;
180  Double_t x0 = 0.0; // [ns]
181  Double_t sigma = 3.5; // [ns]
182 
183  for (int iy = 1; iy <= h2->GetYaxis()->GetNbins(); ++iy) {
184 
185  Double_t x = h1.GetBinCenter (iy);
186  Double_t y = 0.0;
187 
188  for (int ix = imin; ix != imax; ++ix) {
189  y += h2->GetBinContent(ix,iy);
190  }
191 
192  h1.SetBinContent(iy, y);
193 
194  if (y > ymax) {
195  ymax = y;
196  x0 = x;
197  }
198  }
199 
200 
201  TF1* f1 = NULL;
202 
203  if (function == gauss_t) {
204 
205  f1 = new TF1(function.c_str(), "[0]*TMath::Gaus(x,[1],[2]) + [3]");
206 
207  f1->SetParameter(0, ymax);
208  f1->SetParameter(1, x0);
209  f1->SetParameter(2, sigma);
210  f1->SetParameter(3, 0.0);
211 
212  f1->SetParLimits(3, 0.0, ymax);
213 
214  } else if (function == emg_t) {
215 
216  f1 = new TF1(function.c_str(), "[0]*TMath::Exp(0.5*[3]*(2.0*[1]+[3]*[2]*[2]-2.0*x))*TMath::Erfc(([1]+[3]*[2]*[2]-x)/(TMath::Sqrt(2.0)*[2])) +[4]");
217 
218  f1->SetParameter(0, ymax);
219  f1->SetParameter(1, x0 - 0.5*sigma);
220  f1->SetParameter(2, sigma);
221  f1->SetParameter(3, 0.04);
222  f1->SetParameter(4, 0.0);
223 
224  f1->SetParLimits(4, 0.0, ymax);
225 
226  } else {
227 
228  FATAL("Unknown fit function " << function << endl);
229  }
230 
231  // fit
232 
233  Double_t xmin = x0 + T_ns.getLowerLimit();
234  Double_t xmax = x0 + T_ns.getUpperLimit();
235 
236  if (xmin < h1.GetXaxis()->GetXmin()) { xmin = h1.GetXaxis()->GetXmin(); }
237  if (xmax > h1.GetXaxis()->GetXmax()) { xmax = h1.GetXaxis()->GetXmax(); }
238 
239  TFitResultPtr result = h1.Fit(f1, option.c_str(), "same", xmin, xmax);
240 
241  if (debug >= notice_t || !result->IsValid()) {
242  cout << "Bin: "
243  << setw(4) << i << ' '
244  << FIXED(7,3) << f1->GetParameter(1) << " +/- "
245  << FIXED(7,3) << f1->GetParError (1) << ' '
246  << FIXED(7,3) << result->Chi2() << '/'
247  << result->Ndf() << ' '
248  << (result->IsValid() ? "" : "failed") << endl;
249  }
250 
251  if (result->IsValid()) {
252  h0.SetBinContent(i, f1->GetParameter(1));
253  h0.SetBinError (i, f1->GetParError (1));
254  }
255 
256  if (writeFits) {
257  h1.Write();
258  }
259 
260  delete f1;
261  }
262 
263  // Fit
264 
265  TF1 f1("f1", getRisetime, x[0], x[N-1], 3);
266 
267  f1.SetParameter(0, 0.0);
268  f1.SetParameter(1, 2.0e-2);
269  f1.FixParameter(2, -7.0);
270 
271  TFitResultPtr result = h0.Fit(&f1, "S", "same", X.constrain(x[0]), X.constrain(x[N-1]));
272 
273  if (debug >= notice_t || !result->IsValid()) {
274  cout << "Time-slewing correction" << endl;
275  cout << FIXED(7,3) << result->Chi2() << '/'
276  << result->Ndf() << ' '
277  << (result->IsValid() ? "" : "failed") << endl;
278 
279  for (int i = 0; i != f1.GetNpar(); ++i) {
280  cout << "parameter[" << i << "] = " << FIXED(8,4) << f1.GetParameter(i) << " +/- " << FIXED(8,4) << f1.GetParError (i) << endl;
281  }
282  }
283 
284  cout << "// Produced by JFrodo.cc; to be included on JGetRiseTime.hh." << endl;
285 
286  const double t0 = f1.Eval(TIME_OVER_THRESHOLD_NS);
287 
288  for (int i = 0; i != 256; ++i) {
289  cout << "this->push_back(" << FIXED(6,2) << f1.Eval((Double_t) i) - t0 << ");" << endl;
290  }
291 
292  h0.Write();
293 
294  out.Close();
295 }
JPMTAnalogueSignalProcessor.hh
FIXED
Auxiliary data structure for floating point format specification.
Definition: JPrint.hh:481
JMessage.hh
JPrint.hh
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:29
JEEP::notice_t
notice
Definition: JMessage.hh:32
JPARSER::initialised
Empty structure for specification of parser element that is initialised (i.e.
Definition: JParser.hh:63
std::vector
Definition: JSTDTypes.hh:12
main
int main(int argc, char **argv)
Definition: JFrodo.cc:64
JTOOLS::JRange< double >
JDETECTOR::JPMTAnalogueSignalProcessor::setPMTParameters
void setPMTParameters(const JPMTParameters &parameters)
Set PMT parameters.
Definition: JPMTAnalogueSignalProcessor.hh:319
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1493
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
MAKE_CSTRING
#define MAKE_CSTRING(A)
Make C-string.
Definition: JPrint.hh:708
JRange.hh
debug
int debug
debug level
Definition: JSirene.cc:59
JDETECTOR::JPMTParameters::getProperties
JProperties getProperties(const JEquationParameters &equation=JPMTParameters::getEquationParameters())
Get properties of this class.
Definition: JPMTParameters.hh:199
JTOOLS::result
return result
Definition: JPolint.hh:695
JDETECTOR::JPMTAnalogueSignalProcessor::getRiseTime
double getRiseTime(const double npe, const double th) const
Get time to pass from threshold to top of analogue pulse.
Definition: JPMTAnalogueSignalProcessor.hh:214
JCalibration.hh
JParser.hh
JDETECTOR::JPMTSignalProcessorInterface::getNPE
double getNPE(const double tot_ns, const double eps=1.0e-3) const
Get number of photo-electrons.
Definition: JPMTSignalProcessorInterface.hh:316
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1954
JEEP::JProperties
Utility class to parse parameter values.
Definition: JProperties.hh:496
std
Definition: jaanetDictionary.h:36
p1
TPaveText * p1
Definition: JDrawModule3D.cc:35
FATAL
#define FATAL(A)
Definition: JMessage.hh:67
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
JDETECTOR::JPMTParameters
Data structure for PMT parameters.
Definition: JPMTParameters.hh:29
JEEP::debug_t
debug
Definition: JMessage.hh:29
JDETECTOR::JPMTAnalogueSignalProcessor
PMT analogue signal processor.
Definition: JPMTAnalogueSignalProcessor.hh:46
Generated by   doxygen 1.8.16