8 #include <TVirtualFitter.h>
23 SampleFit_PowerLawDoubleExp_PDFMaker::SampleFit_PowerLawDoubleExp_PDFMaker()
25 gStyle->SetOptFit(1111);
27 m_canvas =
new TCanvas(
"SampleFit_PowerLawDoubleExp_PDFMaker",
"SampleFit_PowerLawDoubleExp_PDFMaker");
28 m_pavedtext =
new TPaveText(.05, .1, .95, .8);
30 m_pavedtext->AddText(
"SampleFit_PowerLawDoubleExp Fit output");
31 m_pavedtext->AddText(
"A double-component power-law exponential fit of time-dependent ADC pulses.");
32 m_pavedtext->AddText(
"Magenta curve is the sum of the two component, the red and blue curves.");
33 m_pavedtext->AddText(
"Red dot denote the max points");
36 m_canvas->Print(
"SampleFit_PowerLawDoubleExp.pdf(");
38 SampleFit_PowerLawDoubleExp_PDFMaker::~SampleFit_PowerLawDoubleExp_PDFMaker()
40 if (m_pavedtext)
delete m_pavedtext;
41 if (m_canvas)
delete m_canvas;
43 m_canvas =
new TCanvas(
"SampleFit_PowerLawDoubleExp_PDFMaker",
"SampleFit_PowerLawDoubleExp_PDFMaker");
44 m_pavedtext =
new TPaveText(.05, .1, .95, .8);
46 m_pavedtext->AddText(
"SampleFit_PowerLawDoubleExp Fit output");
47 m_pavedtext->AddText(
"End of pages");
50 m_canvas->Print(
"SampleFit_PowerLawDoubleExp.pdf)");
53 void SampleFit_PowerLawDoubleExp_PDFMaker::MakeSectionPage(
const string &
title)
55 if (m_pavedtext)
delete m_pavedtext;
56 if (m_canvas)
delete m_canvas;
58 m_canvas =
new TCanvas(
"SampleFit_PowerLawDoubleExp_PDFMaker",
"SampleFit_PowerLawDoubleExp_PDFMaker");
60 m_pavedtext =
new TPaveText(.05, .1, .95, .8);
62 m_pavedtext->AddText(title.c_str());
65 m_canvas->Print(
"SampleFit_PowerLawDoubleExp.pdf");
69 const std::vector<double> &samples,
73 std::map<int, double> ¶meters_io,
76 static const int n_parameter = 7;
82 const int n_samples = samples.size();
84 TGraph gpulse(n_samples);
85 for (
int i = 0; i < n_samples; i++)
87 (gpulse.GetX())[i] = i;
89 (gpulse.GetY())[i] = samples[i];
101 pedestal = gpulse.GetY()[0];
102 double peakval = pedestal;
103 const double risetime = 1.5;
105 for (
int iSample = 0; iSample < n_samples - risetime * 3; iSample++)
107 if (
abs(gpulse.GetY()[iSample] - pedestal) >
abs(peakval - pedestal))
109 peakval = gpulse.GetY()[iSample];
117 cout <<
"SampleFit_PowerLawDoubleExp - "
118 <<
"pedestal = " << pedestal <<
", "
119 <<
"peakval = " << peakval <<
", "
120 <<
"peakPos = " << peakPos << endl;
124 struct default_values_t
126 default_values_t(
double default_value,
double min_value,
double max_value)
137 vector<default_values_t> default_values(n_parameter, default_values_t(numeric_limits<double>::signaling_NaN(), numeric_limits<double>::signaling_NaN(), numeric_limits<double>::signaling_NaN()));
139 default_values[0] = default_values_t(peakval * .7, peakval * -1.5, peakval * 1.5);
140 default_values[1] = default_values_t(peakPos - risetime, peakPos - 3 * risetime, peakPos + risetime);
141 default_values[2] = default_values_t(5., 1, 10.);
142 default_values[3] = default_values_t(risetime, risetime * .2, risetime * 10);
143 default_values[4] = default_values_t(pedestal, pedestal -
abs(peakval), pedestal +
abs(peakval));
146 default_values[5] = default_values_t(0, 0, 0);
147 default_values[6] = default_values_t(risetime, risetime, risetime);
151 fits.SetParNames(
"Amplitude",
"Sample Start",
"Power",
"Peak Time 1",
"Pedestal",
"Amplitude ratio",
"Peak Time 2");
153 for (
int i = 0; i < n_parameter; ++i)
155 if (parameters_io.find(i) == parameters_io.end())
157 fits.SetParameter(i, default_values[i].def);
159 if (default_values[i].
min < default_values[i].
max)
161 fits.SetParLimits(i, default_values[i].
min, default_values[i].max);
165 fits.FixParameter(i, default_values[i].def);
170 cout <<
"SampleFit_PowerLawDoubleExp - parameter [" << i <<
"]: "
171 <<
"default value = " << default_values[i].def
172 <<
", min value = " << default_values[i].min
173 <<
", max value = " << default_values[i].max << endl;
179 fits.SetParameter(i, parameters_io[i]);
180 fits.FixParameter(i, parameters_io[i]);
184 cout <<
"SampleFit_PowerLawDoubleExp - parameter [" << i <<
"]: fixed to " << parameters_io[i] << endl;
190 gpulse.Fit(&fits,
"QRN0W",
"goff", 0., (
double) n_samples);
192 gpulse.Fit(&fits,
"RN0VW+",
"goff", 0., (
double) n_samples);
195 pedestal = fits.GetParameter(4);
197 const double peakpos1 = fits.GetParameter(3);
198 const double peakpos2 = fits.GetParameter(6);
199 double max_peakpos = fits.GetParameter(1) + (peakpos1 > peakpos2 ? peakpos1 : peakpos2);
200 if (max_peakpos > n_samples - 1) max_peakpos = n_samples - 1;
202 if (fits.GetParameter(0) > 0)
203 peak_sample = fits.GetMaximumX(fits.GetParameter(1), max_peakpos);
205 peak_sample = fits.GetMinimumX(fits.GetParameter(1), max_peakpos);
207 peak = fits.Eval(peak_sample) - pedestal;
214 string c_name(
string(
"SampleFit_PowerLawDoubleExp_") +
to_string(
id));
216 TCanvas *canvas =
new TCanvas(
217 c_name.c_str(), c_name.c_str());
220 TGraph *g_plot =
static_cast<TGraph *
>(gpulse.DrawClone(
"ap*l"));
221 g_plot->SetTitle((
string(
"ADC data and fit #") +
to_string(
id) +
string(
";Sample number;ADC value")).c_str());
223 fits.SetLineColor(kMagenta);
224 fits.DrawClone(
"same");
229 fits.GetParameter(0) * (1 - fits.GetParameter(5)) / pow(fits.GetParameter(3), fits.GetParameter(2)) * exp(fits.GetParameter(2)),
230 fits.GetParameter(1),
231 fits.GetParameter(2),
232 fits.GetParameter(2) / fits.GetParameter(3),
233 fits.GetParameter(4));
234 f1.SetLineColor(kBlue);
235 f1.DrawClone(
"same");
239 fits.GetParameter(0) * fits.GetParameter(5) / pow(fits.GetParameter(6), fits.GetParameter(2)) * exp(fits.GetParameter(2)),
240 fits.GetParameter(1),
241 fits.GetParameter(2),
242 fits.GetParameter(2) / fits.GetParameter(6),
243 fits.GetParameter(4));
244 f2.SetLineColor(kRed);
245 f2.DrawClone(
"same");
249 g_max.GetX()[0] = peak_sample;
250 g_max.GetY()[0] = peak + pedestal;
252 g_max.SetMarkerStyle(kFullCircle);
253 g_max.SetMarkerSize(2);
254 g_max.SetMarkerColor(kRed);
256 static_cast<TGraph *
>(g_max.DrawClone(
"p"));
264 canvas->Print(
"SampleFit_PowerLawDoubleExp.pdf");
267 for (
int i = 0; i < n_parameter; ++i)
269 parameters_io[i] = fits.GetParameter(i);
274 cout <<
"SampleFit_PowerLawDoubleExp - "
275 <<
"peak_sample = " << peak_sample <<
", "
276 <<
"max_peakpos = " << max_peakpos <<
", "
277 <<
"fits.GetParameter(1) = " << fits.GetParameter(1) <<
", "
278 <<
"peak = " << peak <<
", "
279 <<
"pedestal = " << pedestal << endl;
288 double pedestal = par[4];
293 double signal = par[0] * pow((x[0] - par[1]), par[2]) * exp(-(x[0] - par[1]) * par[3]);
294 return pedestal + signal;
300 double pedestal = par[4];
310 * pow((x[0] - par[1]), par[2])
311 * (((1. - par[5]) / pow(par[3], par[2]) * exp(par[2])) * exp(-(x[0] - par[1]) * (par[2] / par[3]))
312 + (par[5] / pow(par[6], par[2]) * exp(par[2])) * exp(-(x[0] - par[1]) * (par[2] / par[6]))
314 return pedestal + signal;
319 static const int pad_map_to_xy[512][2] = {
833 static const int sampa_chan_to_pad[] = {
1093 if (fee_channel >= 256)
1096 const int pad_number = sampa_chan_to_pad[fee_channel];
1097 const int pad_x = pad_map_to_xy[pad_number][0];
1098 const int pad_y = pad_map_to_xy[pad_number][1];
1100 return make_pair(pad_x, pad_y);