packages/tpc2019/TpcPrototypeDefs.cc

0001 #include "TpcPrototypeDefs.h"
0002
0003 #include <TAttMarker.h>
0004 #include <TCanvas.h>
0005 #include <TF1.h>
0006 #include <TGraph.h>
0007 #include <TPaveText.h>
0008 #include <TStyle.h>
0009
0010 #include <cmath>
0011 #include <iostream>
0012 #include <limits>
0013 #include <memory>
0014 #include <string>
0015
0016 using namespace std;
0017
0018 namespace TpcPrototypeDefs
0019 {
0020 //! TPC v1 FEE test stand decoder
0021 namespace FEEv2
0022 {
0023 SampleFit_PowerLawDoubleExp_PDFMaker::SampleFit_PowerLawDoubleExp_PDFMaker()
0024 {
0025   gStyle->SetOptFit(1111);
0026
0027   m_canvas = new TCanvas("SampleFit_PowerLawDoubleExp_PDFMaker", "SampleFit_PowerLawDoubleExp_PDFMaker");
0028   m_pavedtext = new TPaveText(.05, .1, .95, .8);
0029
0030   m_pavedtext->AddText("SampleFit_PowerLawDoubleExp Fit output");
0031   m_pavedtext->AddText("A double-component power-law exponential fit of time-dependent ADC pulses.");
0032   m_pavedtext->AddText("Magenta curve is the sum of the two component, the red and blue curves.");
0033   m_pavedtext->AddText("Red dot denote the max points");
0034   m_pavedtext->Draw();
0035
0036   m_canvas->Print("SampleFit_PowerLawDoubleExp.pdf(");  //open multiplage PDF
0037 }
0038 SampleFit_PowerLawDoubleExp_PDFMaker::~SampleFit_PowerLawDoubleExp_PDFMaker()
0039 {
0040   if (m_pavedtext) delete m_pavedtext;
0041   if (m_canvas) delete m_canvas;
0042
0043   m_canvas = new TCanvas("SampleFit_PowerLawDoubleExp_PDFMaker", "SampleFit_PowerLawDoubleExp_PDFMaker");
0044   m_pavedtext = new TPaveText(.05, .1, .95, .8);
0045
0046   m_pavedtext->AddText("SampleFit_PowerLawDoubleExp Fit output");
0047   m_pavedtext->AddText("End of pages");
0048   m_pavedtext->Draw();
0049
0050   m_canvas->Print("SampleFit_PowerLawDoubleExp.pdf)");  //close multiplage PDF
0051 }
0052
0053 void SampleFit_PowerLawDoubleExp_PDFMaker::MakeSectionPage(const string &title)
0054 {
0055   if (m_pavedtext) delete m_pavedtext;
0056   if (m_canvas) delete m_canvas;
0057
0058   m_canvas = new TCanvas("SampleFit_PowerLawDoubleExp_PDFMaker", "SampleFit_PowerLawDoubleExp_PDFMaker");
0059
0060   m_pavedtext = new TPaveText(.05, .1, .95, .8);
0061
0062   m_pavedtext->AddText(title.c_str());
0063   m_pavedtext->Draw();
0064
0065   m_canvas->Print("SampleFit_PowerLawDoubleExp.pdf");
0066 }
0067
0068 bool SampleFit_PowerLawDoubleExp(        //
0069     const std::vector<double> &samples,  //
0070     double &peak,                        //
0071     double &peak_sample,                 //
0072     double &pedestal,                    //
0073     std::map<int, double> &parameters_io,
0074     const int verbosity)
0075 {
0076   static const int n_parameter = 7;
0077
0078   // inital guesses
0079   int peakPos = 0.;
0080
0081   //  assert(samples.size() == n_samples);
0082   const int n_samples = samples.size();
0083
0084   TGraph gpulse(n_samples);
0085   for (int i = 0; i < n_samples; i++)
0086   {
0087     (gpulse.GetX())[i] = i;
0088
0089     (gpulse.GetY())[i] = samples[i];
0090   }
0091
0092   //Saturation correction - Abhisek
0093   //  for (int ipoint = 0; ipoint < gpulse.GetN(); ipoint++)
0094   //    if ((gpulse.GetY())[ipoint] >= ((1 << 10) - 10)  // drop point if touching max or low limit on ADCs
0095   //        or (not isnormal((gpulse.GetY())[ipoint])))
0096   //    {
0097   //      gpulse.RemovePoint(ipoint);
0098   //      ipoint--;
0099   //    }
0100
0101   pedestal = gpulse.GetY()[0];  //(double) PEDESTAL;
0102   double peakval = pedestal;
0103   const double risetime = 1.5;
0104
0105   for (int iSample = 0; iSample < n_samples - risetime * 3; iSample++)
0106   {
0107     if (abs(gpulse.GetY()[iSample] - pedestal) > abs(peakval - pedestal))
0108     {
0109       peakval = gpulse.GetY()[iSample];
0110       peakPos = iSample;
0111     }
0112   }
0113   peakval -= pedestal;
0114
0115   if (verbosity)
0116   {
0117     cout << "SampleFit_PowerLawDoubleExp - "
0118          << "pedestal = " << pedestal << ", "
0119          << "peakval = " << peakval << ", "
0120          << "peakPos = " << peakPos << endl;
0121   }
0122
0123   // build default value
0124   struct default_values_t
0125   {
0126     default_values_t(double default_value, double min_value, double max_value)
0127       : def(default_value)
0128       , min(min_value)
0129       , max(max_value)
0130     {
0131     }
0132     double def;
0133     double min;
0134     double max;
0135   };
0136
0137   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()));
0138
0139   default_values[0] = default_values_t(peakval * .7, peakval * -1.5, peakval * 1.5);
0140   default_values[1] = default_values_t(peakPos - risetime, peakPos - 3 * risetime, peakPos + risetime);
0141   default_values[2] = default_values_t(5., 1, 10.);
0142   default_values[3] = default_values_t(risetime, risetime * .2, risetime * 10);
0143   default_values[4] = default_values_t(pedestal, pedestal - abs(peakval), pedestal + abs(peakval));
0144   //  default_values[5] = default_values_t(0.3, 0, 1);
0145   //  default_values[6] = default_values_t(5, risetime * .2, risetime * 10);
0146   default_values[5] = default_values_t(0, 0, 0);  // disable 2nd component
0147   default_values[6] = default_values_t(risetime, risetime, risetime);
0148
0149   // fit function
0150   static TF1 fits("f_SignalShape_PowerLawDoubleExp", SignalShape_PowerLawDoubleExp, 0., n_samples, n_parameter);
0151   fits.SetParNames("Amplitude", "Sample Start", "Power", "Peak Time 1", "Pedestal", "Amplitude ratio", "Peak Time 2");
0152
0153   for (int i = 0; i < n_parameter; ++i)
0154   {
0155     if (parameters_io.find(i) == parameters_io.end())
0156     {
0157       fits.SetParameter(i, default_values[i].def);
0158
0159       if (default_values[i].min < default_values[i].max)
0160       {
0161         fits.SetParLimits(i, default_values[i].min, default_values[i].max);
0162       }
0163       else
0164       {
0165         fits.FixParameter(i, default_values[i].def);
0166       }
0167
0168       if (verbosity)
0169       {
0170         cout << "SampleFit_PowerLawDoubleExp - parameter [" << i << "]: "
0171              << "default value = " << default_values[i].def
0172              << ", min value = " << default_values[i].min
0173              << ", max value = " << default_values[i].max << endl;
0174       }
0175     }
0176     else
0177     {
0178 //      fits.SetParLimits(i, parameters_io[i], parameters_io[i]);
0179       fits.SetParameter(i, parameters_io[i]);
0180       fits.FixParameter(i, parameters_io[i]);
0181
0182       if (verbosity)
0183       {
0184         cout << "SampleFit_PowerLawDoubleExp - parameter [" << i << "]: fixed to " << parameters_io[i] << endl;
0185       }
0186     }
0187   }
0188
0189   if (verbosity <= 1)
0190     gpulse.Fit(&fits, "QRN0W", "goff", 0., (double) n_samples);
0191   else
0192     gpulse.Fit(&fits, "RN0VW+", "goff", 0., (double) n_samples);
0193
0194   // store results
0195   pedestal = fits.GetParameter(4);
0196
0197   const double peakpos1 = fits.GetParameter(3);
0198   const double peakpos2 = fits.GetParameter(6);
0199   double max_peakpos = fits.GetParameter(1) + (peakpos1 > peakpos2 ? peakpos1 : peakpos2);
0200   if (max_peakpos > n_samples - 1) max_peakpos = n_samples - 1;
0201
0202   if (fits.GetParameter(0) > 0)
0203     peak_sample = fits.GetMaximumX(fits.GetParameter(1), max_peakpos);
0204   else
0205     peak_sample = fits.GetMinimumX(fits.GetParameter(1), max_peakpos);
0206
0207   peak = fits.Eval(peak_sample) - pedestal;
0208
0209   if (verbosity)
0210   {
0211     static int id = 0;
0212     ++id;
0213
0214     string c_name(string("SampleFit_PowerLawDoubleExp_") + to_string(id));
0215
0216     TCanvas *canvas = new TCanvas(
0217         c_name.c_str(), c_name.c_str());
0218     canvas->Update();
0219
0220     TGraph *g_plot = static_cast<TGraph *>(gpulse.DrawClone("ap*l"));
0221     g_plot->SetTitle((string("ADC data and fit #") + to_string(id) + string(";Sample number;ADC value")).c_str());
0222
0223     fits.SetLineColor(kMagenta);
0224     fits.DrawClone("same");
0225     fits.Print();
0226
0227     TF1 f1("f_SignalShape_PowerLawExp1", SignalShape_PowerLawExp, 0., n_samples, 5);
0228     f1.SetParameters(
0229         fits.GetParameter(0) * (1 - fits.GetParameter(5)) / pow(fits.GetParameter(3), fits.GetParameter(2)) * exp(fits.GetParameter(2)),
0230         fits.GetParameter(1),
0231         fits.GetParameter(2),
0232         fits.GetParameter(2) / fits.GetParameter(3),
0233         fits.GetParameter(4));
0234     f1.SetLineColor(kBlue);
0235     f1.DrawClone("same");
0236
0237     TF1 f2("f_SignalShape_PowerLawExp2", SignalShape_PowerLawExp, 0., n_samples, 5);
0238     f2.SetParameters(
0239         fits.GetParameter(0) * fits.GetParameter(5) / pow(fits.GetParameter(6), fits.GetParameter(2)) * exp(fits.GetParameter(2)),
0240         fits.GetParameter(1),
0241         fits.GetParameter(2),
0242         fits.GetParameter(2) / fits.GetParameter(6),
0243         fits.GetParameter(4));
0244     f2.SetLineColor(kRed);
0245     f2.DrawClone("same");
0246
0247     TGraph g_max(1);
0248
0249     g_max.GetX()[0] = peak_sample;
0250     g_max.GetY()[0] = peak + pedestal;
0251
0252     g_max.SetMarkerStyle(kFullCircle);
0253     g_max.SetMarkerSize(2);
0254     g_max.SetMarkerColor(kRed);
0255
0256     g_max.DrawClone("p");
0257
0258     canvas->Update();
0259
0260     //    if (id == 1)
0261     //    {
0262     //      canvas->Print("SampleFit_PowerLawDoubleExp.pdf(");
0263     //    }
0264     canvas->Print("SampleFit_PowerLawDoubleExp.pdf");
0265   }
0266
0267   for (int i = 0; i < n_parameter; ++i)
0268   {
0269     parameters_io[i] = fits.GetParameter(i);
0270   }
0271
0272   if (verbosity)
0273   {
0274     cout << "SampleFit_PowerLawDoubleExp - "
0275          << "peak_sample = " << peak_sample << ", "
0276          << "max_peakpos = " << max_peakpos << ", "
0277          << "fits.GetParameter(1) = " << fits.GetParameter(1) << ", "
0278          << "peak = " << peak << ", "
0279          << "pedestal = " << pedestal << endl;
0280   }
0281
0282   return true;
0283 }
0284
0285 double
0286 SignalShape_PowerLawExp(double *x, double *par)
0287 {
0288   double pedestal = par[4];
0289   //                        + ((x[0] - 1.5 * par[1]) > 0) * par[5];  // quick fix on exting tails on the signal function
0290   if (x[0] < par[1])
0291     return pedestal;
0292   //double  signal = (-1)*par[0]*pow((x[0]-par[1]),par[2])*exp(-(x[0]-par[1])*par[3]);
0293   double signal = par[0] * pow((x[0] - par[1]), par[2]) * exp(-(x[0] - par[1]) * par[3]);
0294   return pedestal + signal;
0295 }
0296
0297 double
0298 SignalShape_PowerLawDoubleExp(double *x, double *par)
0299 {
0300   double pedestal = par[4];
0301   //                        + ((x[0] - 1.5 * par[1]) > 0) * par[5];  // quick fix on exting tails on the signal function
0302   if (x[0] < par[1])
0303     return pedestal;
0304   //double  signal = (-1)*par[0]*pow((x[0]-par[1]),par[2])*exp(-(x[0]-par[1])*par[3]);
0305   //  peak / pow(fits.GetParameter(2) / fits.GetParameter(3), fits.GetParameter(2)) * exp(fits.GetParameter(2)) = fits.GetParameter(0);  // exact peak height is (p0*Power(p2/p3,p2))/Power(E,p2)
0306   //  fits.GetParameter(2) / peak_shift =  fits.GetParameter(3);  // signal peak time
0307
0308   double signal =                                                                                         //
0309       par[0]                                                                                              //
0310       * pow((x[0] - par[1]), par[2])                                                                      //
0311       * (((1. - par[5]) / pow(par[3], par[2]) * exp(par[2])) * exp(-(x[0] - par[1]) * (par[2] / par[3]))  //
0312          + (par[5] / pow(par[6], par[2]) * exp(par[2])) * exp(-(x[0] - par[1]) * (par[2] / par[6]))       //
0313         );
0314   return pedestal + signal;
0315 }
0316
0317
0318 }  // namespace FEEv2
0319
0320 }  // namespace TpcPrototypeDefs