File indexing completed on 2025-08-05 08:12:21
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011 #include "PROTOTYPE4_FEM.h"
0012
0013 #include <TAttMarker.h> // for kFullCircle
0014 #include <TCanvas.h>
0015 #include <TF1.h>
0016 #include <TGraph.h>
0017
0018 #include <cassert>
0019 #include <cmath>
0020 #include <cstdlib> // for exit
0021 #include <iostream>
0022 #include <limits>
0023 #include <memory> // for allocator_traits<>::value_type
0024 #include <string>
0025
0026
0027 using namespace std;
0028
0029 int PROTOTYPE4_FEM::GetChannelNumber(const std::string &caloname, int i_column,
0030 int i_row)
0031 {
0032 if (caloname == "HCALIN")
0033 {
0034 assert(i_row >= 0);
0035 assert(i_row < NCH_IHCAL_ROWS);
0036 assert(i_column >= 0);
0037 assert(i_column < NCH_IHCAL_COLUMNS);
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051 static const int hbdchanIHC_col_0[4] = {67, 66, 65,
0052 64};
0053 static const int hbdchanIHC_col_1[4] = {71, 70, 69, 68};
0054 static const int hbdchanIHC_col_2[4] = {75, 74, 73, 72};
0055 static const int hbdchanIHC_col_3[4] = {79, 78, 77, 76};
0056
0057 if (i_column == 0)
0058 return hbdchanIHC_col_0[i_row];
0059 if (i_column == 1)
0060 return hbdchanIHC_col_1[i_row];
0061 if (i_column == 2)
0062 return hbdchanIHC_col_2[i_row];
0063 if (i_column == 3)
0064 return hbdchanIHC_col_3[i_row];
0065 }
0066 else if (caloname == "HCALOUT")
0067 {
0068 assert(i_row >= 0);
0069 assert(i_row < NCH_OHCAL_ROWS);
0070 assert(i_column >= 0);
0071 assert(i_column < NCH_OHCAL_COLUMNS);
0072
0073 static const int hbdchanOHC_col_0[4] = {116, 118, 112,
0074 114};
0075 static const int hbdchanOHC_col_1[4] = {124, 126, 120, 122};
0076 static const int hbdchanOHC_col_2[4] = {132, 134, 128, 130};
0077 static const int hbdchanOHC_col_3[4] = {140, 142, 136, 138};
0078
0079 if (i_column == 0)
0080 return hbdchanOHC_col_0[i_row];
0081 if (i_column == 1)
0082 return hbdchanOHC_col_1[i_row];
0083 if (i_column == 2)
0084 return hbdchanOHC_col_2[i_row];
0085 if (i_column == 3)
0086 return hbdchanOHC_col_3[i_row];
0087 }
0088 else if (caloname == "EMCAL")
0089 {
0090 assert(i_row >= 0);
0091 assert(i_row < NCH_EMCAL_ROWS);
0092 assert(i_column >= 0);
0093 assert(i_column < NCH_EMCAL_COLUMNS);
0094
0095
0096
0097
0098
0099
0100
0101
0102
0103
0104
0105
0106
0107
0108
0109
0110
0111
0112
0113
0114
0115
0116
0117
0118
0119
0120
0121
0122
0123
0124
0125
0126
0127
0128
0129
0130
0131
0132
0133
0134
0135
0136
0137
0138 const static int canmap[64] = {
0139 12, 13, 8, 9, 4, 5, 0, 1, 14, 15, 10, 11, 6, 7, 2, 3,
0140 28, 29, 24, 25, 20, 21, 16, 17, 30, 31, 26, 27, 22, 23, 18, 19,
0141 44, 45, 40, 41, 36, 37, 32, 33, 46, 47, 42, 43, 38, 39, 34, 35,
0142 60, 61, 56, 57, 52, 53, 48, 49, 62, 63, 58, 59, 54, 55, 50, 51};
0143
0144 const int linear_channel_ID =
0145 (NCH_EMCAL_ROWS - i_row - 1) * NCH_EMCAL_COLUMNS + i_column;
0146
0147 assert(linear_channel_ID >= 0);
0148 assert(linear_channel_ID < 64);
0149
0150 return canmap[linear_channel_ID];
0151 }
0152
0153 std::cout << "PROTOTYPE4_FEM::GetHBDCh - invalid input caloname " << caloname
0154 << " i_column " << i_column << " i_row " << i_row << std::endl;
0155 exit(1);
0156 return -9999;
0157 }
0158
0159 bool PROTOTYPE4_FEM::SampleFit_PowerLawExp(
0160 const std::vector<double> &samples,
0161 double &peak,
0162 double &peak_sample,
0163 double &pedestal,
0164 const int verbosity)
0165 {
0166 int peakPos = 0.;
0167
0168 assert(samples.size() == NSAMPLES);
0169
0170 TGraph gpulse(NSAMPLES);
0171 for (int i = 0; i < NSAMPLES; i++)
0172 {
0173 (gpulse.GetX())[i] = i;
0174
0175 (gpulse.GetY())[i] = samples[i];
0176 }
0177
0178 pedestal = gpulse.GetY()[0];
0179 double peakval = pedestal;
0180 const double risetime = 4;
0181
0182 for (int iSample = 0; iSample < NSAMPLES; iSample++)
0183 {
0184 if (abs(gpulse.GetY()[iSample] - pedestal) > abs(peakval - pedestal))
0185 {
0186 peakval = gpulse.GetY()[iSample];
0187 peakPos = iSample;
0188 }
0189 }
0190 peakval -= pedestal;
0191
0192
0193 TF1 fits("f_SignalShape_PowerLawExp", SignalShape_PowerLawExp, 0., NSAMPLES,
0194 5);
0195
0196 double par[6] = {0};
0197 par[0] = peakval;
0198 par[1] = peakPos - risetime;
0199 if (par[1] < 0.)
0200 par[1] = 0.;
0201 par[2] = 4.;
0202 par[3] = 1.5;
0203 par[4] = pedestal;
0204 par[5] = 0;
0205 fits.SetParameters(par);
0206 fits.SetParNames("Amplitude", "Sample Start", "Power", "Decay", "pedestal",
0207 "Baseline shift");
0208 fits.SetParLimits(0, peakval * 0.5, peakval * 10);
0209 fits.SetParLimits(1, 0, NSAMPLES);
0210 fits.SetParLimits(2, 0, 10.);
0211 fits.SetParLimits(3, 0, 10);
0212 fits.SetParLimits(4, pedestal - abs(peakval), pedestal + abs(peakval));
0213
0214
0215
0216
0217 for (int ipoint = 0; ipoint < gpulse.GetN(); ipoint++)
0218 if ((gpulse.GetY())[ipoint] <= 10 or
0219 (gpulse.GetY())[ipoint] >=
0220 ((1 << 14) - 10))
0221 {
0222 gpulse.RemovePoint(ipoint);
0223 ipoint--;
0224 }
0225
0226 if (verbosity <= 1)
0227 gpulse.Fit(&fits, "MQRN0W", "goff", 0., (double) NSAMPLES);
0228 else
0229 gpulse.Fit(&fits, "MRN0VW", "goff", 0., (double) NSAMPLES);
0230
0231 if (verbosity)
0232 {
0233 static int id = 0;
0234 ++id;
0235
0236 string c_name(string("PROTOTYPE4_FEM_SampleFit_PowerLawExp_") +
0237 to_string(id));
0238
0239 TCanvas *canvas = new TCanvas(c_name.c_str(), c_name.c_str());
0240 canvas->Update();
0241
0242 TGraph *g_plot = static_cast<TGraph *>(gpulse.DrawClone("ap*l"));
0243 g_plot->SetTitle("ADC data and fit;Sample number;ADC value");
0244 fits.DrawClone("same");
0245 fits.Print();
0246 canvas->Update();
0247
0248 }
0249
0250
0251 peak =
0252 (fits.GetParameter(0) *
0253 pow(fits.GetParameter(2) / fits.GetParameter(3), fits.GetParameter(2))) /
0254 exp(fits.GetParameter(
0255 2));
0256
0257
0258 peak_sample = fits.GetParameter(1) +
0259 fits.GetParameter(2) / fits.GetParameter(3);
0260
0261
0262
0263 pedestal = fits.GetParameter(4);
0264
0265 return true;
0266 }
0267
0268 double PROTOTYPE4_FEM::SignalShape_PowerLawExp(double *x, double *par)
0269 {
0270 double pedestal = par[4];
0271
0272
0273 if (x[0] < par[1])
0274 return pedestal;
0275
0276
0277 double signal =
0278 par[0] * pow((x[0] - par[1]), par[2]) * exp(-(x[0] - par[1]) * par[3]);
0279 return pedestal + signal;
0280 }
0281
0282 bool PROTOTYPE4_FEM::SampleFit_PowerLawDoubleExp(
0283 const std::vector<double> &samples,
0284 double &peak,
0285 double &peak_sample,
0286 double &pedestal,
0287 std::map<int, double> ¶meters_io, const int verbosity)
0288 {
0289 static const int n_parameter = 7;
0290
0291
0292 int peakPos = 0.;
0293
0294 assert(samples.size() == NSAMPLES);
0295
0296 TGraph gpulse(NSAMPLES);
0297 for (int i = 0; i < NSAMPLES; i++)
0298 {
0299 (gpulse.GetX())[i] = i;
0300
0301 (gpulse.GetY())[i] = samples[i];
0302 }
0303
0304
0305 for (int ipoint = 0; ipoint < gpulse.GetN(); ipoint++)
0306 if ((gpulse.GetY())[ipoint] <= 10 or
0307 (gpulse.GetY())[ipoint] >=
0308 ((1 << 14) - 10)
0309 or (not isnormal((gpulse.GetY())[ipoint])))
0310 {
0311 gpulse.RemovePoint(ipoint);
0312 ipoint--;
0313 }
0314
0315 pedestal = gpulse.GetY()[0];
0316 double peakval = pedestal;
0317 const double risetime = 2;
0318
0319 for (int iSample = 0; iSample < NSAMPLES - risetime * 3; iSample++)
0320 {
0321 if (abs(gpulse.GetY()[iSample] - pedestal) > abs(peakval - pedestal))
0322 {
0323 peakval = gpulse.GetY()[iSample];
0324 peakPos = iSample;
0325 }
0326 }
0327 peakval -= pedestal;
0328
0329 if (verbosity)
0330 {
0331 cout << "PROTOTYPE4_FEM::SampleFit_PowerLawDoubleExp - "
0332 << "pedestal = " << pedestal << ", "
0333 << "peakval = " << peakval << ", "
0334 << "peakPos = " << peakPos << endl;
0335 }
0336
0337
0338 struct default_values_t
0339 {
0340 default_values_t(double default_value, double min_value, double max_value)
0341 : def(default_value)
0342 , min(min_value)
0343 , max(max_value)
0344 {
0345 }
0346 double def;
0347 double min;
0348 double max;
0349 };
0350
0351 vector<default_values_t> default_values(
0352 n_parameter, default_values_t(numeric_limits<double>::signaling_NaN(),
0353 numeric_limits<double>::signaling_NaN(),
0354 numeric_limits<double>::signaling_NaN()));
0355
0356 default_values[0] =
0357 default_values_t(peakval * .7, peakval * -1.5, peakval * 1.5);
0358 default_values[1] = default_values_t(
0359 peakPos - risetime, peakPos - 3 * risetime, peakPos + risetime);
0360 default_values[2] = default_values_t(2., 1, 5.);
0361 default_values[3] = default_values_t(5, risetime * .5, risetime * 4);
0362 default_values[4] = default_values_t(pedestal, pedestal - abs(peakval),
0363 pedestal + abs(peakval));
0364 default_values[5] = default_values_t(.3, 0, 1);
0365 default_values[6] = default_values_t(5, risetime * .5, risetime * 4);
0366
0367
0368 TF1 fits("f_SignalShape_PowerLawDoubleExp", SignalShape_PowerLawDoubleExp, 0.,
0369 NSAMPLES, n_parameter);
0370 fits.SetParNames("Amplitude", "Sample Start", "Power", "Peak Time 1",
0371 "Pedestal", "Amplitude ratio", "Peak Time 2");
0372
0373 for (int i = 0; i < n_parameter; ++i)
0374 {
0375 if (parameters_io.find(i) == parameters_io.end())
0376 {
0377 fits.SetParameter(i, default_values[i].def);
0378 fits.SetParLimits(i, default_values[i].min, default_values[i].max);
0379
0380 if (verbosity)
0381 {
0382 cout << "PROTOTYPE4_FEM::SampleFit_PowerLawDoubleExp - parameter [" << i
0383 << "]: "
0384 << "default value = " << default_values[i].def
0385 << ", min value = " << default_values[i].min
0386 << ", max value = " << default_values[i].max << endl;
0387 }
0388 }
0389 else
0390 {
0391 fits.SetParLimits(i, parameters_io[i], parameters_io[i]);
0392 fits.SetParameter(i, parameters_io[i]);
0393
0394 if (verbosity)
0395 {
0396 cout << "PROTOTYPE4_FEM::SampleFit_PowerLawDoubleExp - parameter [" << i
0397 << "]: fixed to " << parameters_io[i] << endl;
0398 }
0399 }
0400 }
0401
0402 if (verbosity <= 1)
0403 gpulse.Fit(&fits, "QRN0W", "goff", 0., (double) NSAMPLES);
0404 else
0405 gpulse.Fit(&fits, "RN0VW", "goff", 0., (double) NSAMPLES);
0406
0407 if (verbosity)
0408 {
0409 static int id = 0;
0410 ++id;
0411
0412 string c_name(string("PROTOTYPE4_FEM_SampleFit_PowerLawDoubleExp_") +
0413 to_string(id));
0414
0415 TCanvas *canvas = new TCanvas(c_name.c_str(), c_name.c_str());
0416 canvas->Update();
0417
0418 TGraph *g_plot = static_cast<TGraph *>(gpulse.DrawClone("ap*l"));
0419 g_plot->SetTitle("ADC data and fit;Sample number;ADC value");
0420
0421 fits.SetLineColor(kMagenta);
0422 fits.DrawClone("same");
0423 fits.Print();
0424
0425 TF1 f1("f_SignalShape_PowerLawExp1", SignalShape_PowerLawExp, 0., NSAMPLES,
0426 5);
0427 f1.SetParameters(fits.GetParameter(0) * (1 - fits.GetParameter(5)) /
0428 pow(fits.GetParameter(3), fits.GetParameter(2)) *
0429 exp(fits.GetParameter(2)),
0430 fits.GetParameter(1), fits.GetParameter(2),
0431 fits.GetParameter(2) / fits.GetParameter(3),
0432 fits.GetParameter(4));
0433 f1.SetLineColor(kBlue);
0434 f1.DrawClone("same");
0435
0436 TF1 f2("f_SignalShape_PowerLawExp2", SignalShape_PowerLawExp, 0., NSAMPLES,
0437 5);
0438 f2.SetParameters(fits.GetParameter(0) * fits.GetParameter(5) /
0439 pow(fits.GetParameter(6), fits.GetParameter(2)) *
0440 exp(fits.GetParameter(2)),
0441 fits.GetParameter(1), fits.GetParameter(2),
0442 fits.GetParameter(2) / fits.GetParameter(6),
0443 fits.GetParameter(4));
0444 f2.SetLineColor(kRed);
0445 f2.DrawClone("same");
0446
0447 canvas->Update();
0448 }
0449
0450
0451 pedestal = fits.GetParameter(4);
0452
0453 const double peakpos1 = fits.GetParameter(3);
0454 const double peakpos2 = fits.GetParameter(6);
0455 double max_peakpos =
0456 fits.GetParameter(1) + (peakpos1 > peakpos2 ? peakpos1 : peakpos2);
0457 if (max_peakpos > NSAMPLES - 1)
0458 max_peakpos = NSAMPLES - 1;
0459
0460 if (fits.GetParameter(0) > 0)
0461 peak_sample = fits.GetMaximumX(fits.GetParameter(1), max_peakpos);
0462 else
0463 peak_sample = fits.GetMinimumX(fits.GetParameter(1), max_peakpos);
0464
0465 peak = fits.Eval(peak_sample) - pedestal;
0466
0467 if (verbosity)
0468 {
0469 TGraph g_max(NSAMPLES);
0470
0471 g_max.GetX()[0] = peak_sample;
0472 g_max.GetY()[0] = peak + pedestal;
0473
0474 g_max.SetMarkerStyle(kFullCircle);
0475 g_max.SetMarkerSize(2);
0476 g_max.SetMarkerColor(kRed);
0477
0478 g_max.DrawClone("p");
0479 }
0480
0481 for (int i = 0; i < n_parameter; ++i)
0482 {
0483 parameters_io[i] = fits.GetParameter(i);
0484 }
0485
0486 if (verbosity)
0487 {
0488 cout << "PROTOTYPE4_FEM::SampleFit_PowerLawDoubleExp - "
0489 << "peak_sample = " << peak_sample << ", "
0490 << "max_peakpos = " << max_peakpos << ", "
0491 << "fits.GetParameter(1) = " << fits.GetParameter(1) << ", "
0492 << "peak = " << peak << ", "
0493 << "pedestal = " << pedestal << endl;
0494 }
0495
0496 return true;
0497 }
0498
0499 double PROTOTYPE4_FEM::SignalShape_PowerLawDoubleExp(double *x, double *par)
0500 {
0501 double pedestal = par[4];
0502
0503
0504 if (x[0] < par[1])
0505 return pedestal;
0506
0507
0508
0509
0510
0511
0512
0513
0514 double signal =
0515 par[0]
0516 * pow((x[0] - par[1]), par[2])
0517 * (((1. - par[5]) / pow(par[3], par[2]) * exp(par[2])) *
0518 exp(-(x[0] - par[1]) * (par[2] / par[3]))
0519 + (par[5] / pow(par[6], par[2]) * exp(par[2])) *
0520 exp(-(x[0] - par[1]) * (par[2] / par[6]))
0521 );
0522 return pedestal + signal;
0523 }
0524
0525 bool PROTOTYPE4_FEM::SampleFit_PeakSample(
0526 const std::vector<double> &samples,
0527 double &peak,
0528 double &peak_sample,
0529 double &pedestal,
0530 const int verbosity)
0531 {
0532 int peakPos = 0.;
0533
0534 assert(samples.size() == NSAMPLES);
0535
0536 const static int N_pedestal = 3;
0537 pedestal = 0;
0538 for (int iSample = 0; iSample < N_pedestal; iSample++)
0539 {
0540 pedestal += samples[iSample];
0541 }
0542 pedestal /= N_pedestal;
0543
0544 double peakval = pedestal;
0545 for (int iSample = N_pedestal; iSample < NSAMPLES; iSample++)
0546 {
0547 if (abs(samples[iSample] - pedestal) > abs(peakval - pedestal))
0548 {
0549 peakval = samples[iSample];
0550 peakPos = iSample;
0551 }
0552 }
0553
0554 peak = peakval - pedestal;
0555 peak_sample = peakPos;
0556
0557 if (verbosity)
0558 {
0559 cout << "PROTOTYPE4_FEM::SampleFit_PeakSample - "
0560 << "pedestal = " << pedestal << ", "
0561 << "peakval = " << peakval << ", "
0562 << "peakPos = " << peakPos << endl;
0563 }
0564
0565 return true;
0566 }