12 #include "TGraphErrors.h"
13 #include "TGraphSmooth.h"
26 const char*
const h0_1 =
"h0";
35 int main(
int argc,
char **argv)
40 typedef JRange<double> JRange_t;
53 JParser<> zap(
"Auxiliary program to model transition time distribution generator from data.");
55 zap[
'f'] =
make_field(inputFile,
"input file (output from JPulsar)");
57 zap[
'B'] =
make_field(bass,
"see TGraphSmooth") = 0.00;
58 zap[
'S'] =
make_field(span,
"see TGraphSmooth") = 0.01;
59 zap[
'P'] =
make_field(pdf,
"PDF output file; for inclusion in JPMTTransitTimeProbability.hh") =
"";
60 zap[
'C'] =
make_field(
cdf,
"CDF output file; for inclusion in JPMTTransitTimeGenerator.hh") =
"";
61 zap[
'T'] =
make_field(T_ns,
"time range for PDF and CDF") = JRange_t(-20.0, +100.0);
66 catch(
const exception &error) {
67 FATAL(error.what() << endl);
73 TFile
in(inputFile.c_str(),
"exist");
76 FATAL(
"File " << inputFile <<
" not opened." << endl);
79 h1 =
dynamic_cast<TH1D*
>(
in.Get(h0_1));
82 FATAL(
"No histogram <" << h0_1 <<
">." << endl);
88 TF1
f1(
"f1",
"[0]*exp(-0.5*(x-[1])*(x-[1])/([2]*[2])) + [3]");
97 for (Int_t
i = 1;
i <= h1->GetNbinsX(); ++
i) {
99 const Double_t
x = h1->GetBinCenter (
i);
100 const Double_t
y = h1->GetBinContent(
i);
108 f1.SetParameter(0, ymax);
109 f1.SetParameter(1, x0);
110 f1.SetParameter(2,
sigma);
111 f1.SetParameter(3, ymin);
115 h1->Fit(&f1,
"LL",
"same", x0 - 5 *
sigma, x0 + 5 *
sigma);
119 ymax = f1.GetParameter(0);
120 x0 = f1.GetParameter(1);
121 sigma = f1.GetParameter(2);
122 ymin = f1.GetParameter(3);
134 for (Int_t
i = 1;
i <= h1->GetNbinsX(); ++
i) {
136 const Double_t
x = h1->GetBinCenter (
i);
137 const Double_t
y = h1->GetBinContent(
i);
143 if (x < x0 - 5.0 *
sigma)
145 else if (x > x0 + 5.0 *
sigma)
165 NOTICE(
"Number of pulses: " << yy << endl);
166 NOTICE(
"Pre-pulses: " << yl / yy << endl);
167 NOTICE(
"Delayed pulses: " << yr / yy << endl);
179 for (Int_t
i = 1;
i <= h1->GetNbinsX(); ++
i) {
181 const Double_t
x = h1->GetBinCenter (
i) - x0;
182 const Double_t
y = h1->GetBinContent(
i) - y0;
188 EY.push_back(sqrt(y));
192 const int N =
X.size();
195 FATAL(
"Not enough points." << endl);
199 TGraphErrors g0(
N,
X.data(),
Y.data(), EX.data(), EY.data());
210 for (
int i = 0;
i !=
N; ++
i) {
212 const Double_t
x = g1.GetX()[
i];
213 const Double_t
f = f1.Eval(x + x0);
222 TGraphSmooth gs(
"gs");
224 TGraph* gout = gs.SmoothSuper(&g1,
"", bass, span,
false);
229 for (
int i = 0;
i !=
N; ++
i) {
230 g1.GetY()[
i] = gout->GetY()[
i];
236 for (
int i = 0;
i !=
N; ++
i) {
238 const Double_t
x = g1.GetX()[
i];
239 const Double_t
f = f1.Eval(x + x0);
254 for (
int i = 0;
i !=
N; ++
i) {
258 for (
int i = 0;
i !=
N; ++
i) {
262 ofstream out(pdf.c_str());
264 out <<
"\t// produced by JLegolas.cc" << endl;
266 const Double_t
xmin = T_ns.getLowerLimit();
267 const Double_t
xmax = T_ns.getUpperLimit();
268 const Double_t dx = 0.25;
270 for (Double_t
x = xmin;
x <= xmax + 0.5*dx;
x += dx) {
272 Double_t
y = g2.Eval(
x);
280 << showpos <<
FIXED(7,2) <<
x
282 << noshowpos <<
FIXED(8,6) << y
295 for (
int i = 0;
i+1 !=
N; ++
i) {
296 g2.GetY()[
i+1] += g2.GetY()[
i];
300 const Double_t
xmin = g2.GetX()[ 0 ];
301 const Double_t
xmax = g2.GetX()[
N-1];
302 const Double_t dx = (xmax -
xmin) /
N;
304 const Double_t ymin = g2.GetY()[ 0 ];
305 const Double_t ymax = g2.GetY()[
N-1];
307 for (
int i = 0;
i !=
N; ++
i) {
309 const Double_t
x = g2.GetX()[
i];
310 const Double_t
y = g2.GetY()[
i];
312 g2.GetX()[
i] = (y - ymin) / ymax;
313 g2.GetY()[
i] = x + 0.5 * dx;
317 ofstream out(
cdf.c_str());
319 out <<
"\t// produced by JLegolas.cc" << endl;
321 const Double_t
xmin = 0.0;
322 const Double_t
xmax = 1.0;
323 const Double_t dx = 0.001;
325 for (Double_t
x = xmin;
x <= xmax + 0.5*dx;
x += dx) {
329 << noshowpos <<
FIXED(6,4) <<
x
331 << showpos <<
FIXED(9,5) << g2.Eval(
x)
343 const Double_t
xmin =
X[ 0 ];
344 const Double_t
xmax =
X[
N-1];
345 const Double_t dx = (xmax -
xmin) /
N;
347 TH1D h2(
"pmt", NULL,
N, xmin - 0.5*dx, xmax + 0.5*dx);
355 for (
int i = 0;
i !=
N; ++
i) {
361 for (
int i = 0;
i !=
N; ++
i) {
363 h2.SetBinContent(
i+1,
Y [
i]/W);
364 h2.SetBinError (
i+1, EY[
i]/W);
Utility class to parse command line options.
int main(int argc, char *argv[])
o $QUALITY_ROOT d $DEBUG!CHECK_EXIT_CODE JPlot1D f
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
Auxiliary data structure for floating point format specification.
const JPolynome f1(1.0, 2.0, 3.0)
Function.
I/O formatting auxiliaries.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
General purpose messaging.
then for FUNCTION in pdf npe cdf
then usage $script< input file >[option[primary[working directory]]] nWhere option can be N
Auxiliary class to define a range between two values.
Utility class to parse command line options.
no fit printf nominal n $STRING awk v X
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 JAcoustics sh $DETECTOR_ID source JAcousticsToolkit sh CHECK_EXIT_CODE typeset A EMITTERS get_tripods $WORKDIR tripod txt EMITTERS get_transmitters $WORKDIR transmitter txt EMITTERS for EMITTER in
Double_t g1(const Double_t x)
Function.