Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TpcPrototypeDefs.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file TpcPrototypeDefs.cc
1 #include "TpcPrototypeDefs.h"
2 
3 #include <TAttMarker.h>
4 #include <TCanvas.h>
5 #include <TF1.h>
6 #include <TGraph.h>
7 #include <TPaveText.h>
8 #include <TStyle.h>
9 
10 #include <cmath>
11 #include <iostream>
12 #include <limits>
13 #include <memory>
14 #include <string>
15 
16 using namespace std;
17 
18 namespace TpcPrototypeDefs
19 {
21 namespace FEEv2
22 {
23 SampleFit_PowerLawDoubleExp_PDFMaker::SampleFit_PowerLawDoubleExp_PDFMaker()
24 {
25  gStyle->SetOptFit(1111);
26 
27  m_canvas = new TCanvas("SampleFit_PowerLawDoubleExp_PDFMaker", "SampleFit_PowerLawDoubleExp_PDFMaker");
28  m_pavedtext = new TPaveText(.05, .1, .95, .8);
29 
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");
34  m_pavedtext->Draw();
35 
36  m_canvas->Print("SampleFit_PowerLawDoubleExp.pdf("); //open multiplage PDF
37 }
38 SampleFit_PowerLawDoubleExp_PDFMaker::~SampleFit_PowerLawDoubleExp_PDFMaker()
39 {
40  if (m_pavedtext) delete m_pavedtext;
41  if (m_canvas) delete m_canvas;
42 
43  m_canvas = new TCanvas("SampleFit_PowerLawDoubleExp_PDFMaker", "SampleFit_PowerLawDoubleExp_PDFMaker");
44  m_pavedtext = new TPaveText(.05, .1, .95, .8);
45 
46  m_pavedtext->AddText("SampleFit_PowerLawDoubleExp Fit output");
47  m_pavedtext->AddText("End of pages");
48  m_pavedtext->Draw();
49 
50  m_canvas->Print("SampleFit_PowerLawDoubleExp.pdf)"); //close multiplage PDF
51 }
52 
53 void SampleFit_PowerLawDoubleExp_PDFMaker::MakeSectionPage(const string &title)
54 {
55  if (m_pavedtext) delete m_pavedtext;
56  if (m_canvas) delete m_canvas;
57 
58  m_canvas = new TCanvas("SampleFit_PowerLawDoubleExp_PDFMaker", "SampleFit_PowerLawDoubleExp_PDFMaker");
59 
60  m_pavedtext = new TPaveText(.05, .1, .95, .8);
61 
62  m_pavedtext->AddText(title.c_str());
63  m_pavedtext->Draw();
64 
65  m_canvas->Print("SampleFit_PowerLawDoubleExp.pdf");
66 }
67 
68 bool SampleFit_PowerLawDoubleExp( //
69  const std::vector<double> &samples, //
70  double &peak, //
71  double &peak_sample, //
72  double &pedestal, //
73  std::map<int, double> &parameters_io,
74  const int verbosity)
75 {
76  static const int n_parameter = 7;
77 
78  // inital guesses
79  int peakPos = 0.;
80 
81  // assert(samples.size() == n_samples);
82  const int n_samples = samples.size();
83 
84  TGraph gpulse(n_samples);
85  for (int i = 0; i < n_samples; i++)
86  {
87  (gpulse.GetX())[i] = i;
88 
89  (gpulse.GetY())[i] = samples[i];
90  }
91 
92  //Saturation correction - Abhisek
93  // for (int ipoint = 0; ipoint < gpulse.GetN(); ipoint++)
94  // if ((gpulse.GetY())[ipoint] >= ((1 << 10) - 10) // drop point if touching max or low limit on ADCs
95  // or (not isnormal((gpulse.GetY())[ipoint])))
96  // {
97  // gpulse.RemovePoint(ipoint);
98  // ipoint--;
99  // }
100 
101  pedestal = gpulse.GetY()[0]; //(double) PEDESTAL;
102  double peakval = pedestal;
103  const double risetime = 1.5;
104 
105  for (int iSample = 0; iSample < n_samples - risetime * 3; iSample++)
106  {
107  if (abs(gpulse.GetY()[iSample] - pedestal) > abs(peakval - pedestal))
108  {
109  peakval = gpulse.GetY()[iSample];
110  peakPos = iSample;
111  }
112  }
113  peakval -= pedestal;
114 
115  if (verbosity)
116  {
117  cout << "SampleFit_PowerLawDoubleExp - "
118  << "pedestal = " << pedestal << ", "
119  << "peakval = " << peakval << ", "
120  << "peakPos = " << peakPos << endl;
121  }
122 
123  // build default value
124  struct default_values_t
125  {
126  default_values_t(double default_value, double min_value, double max_value)
127  : def(default_value)
128  , min(min_value)
129  , max(max_value)
130  {
131  }
132  double def;
133  double min;
134  double max;
135  };
136 
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()));
138 
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));
144  // default_values[5] = default_values_t(0.3, 0, 1);
145  // default_values[6] = default_values_t(5, risetime * .2, risetime * 10);
146  default_values[5] = default_values_t(0, 0, 0); // disable 2nd component
147  default_values[6] = default_values_t(risetime, risetime, risetime);
148 
149  // fit function
150  static TF1 fits("f_SignalShape_PowerLawDoubleExp", SignalShape_PowerLawDoubleExp, 0., n_samples, n_parameter);
151  fits.SetParNames("Amplitude", "Sample Start", "Power", "Peak Time 1", "Pedestal", "Amplitude ratio", "Peak Time 2");
152 
153  for (int i = 0; i < n_parameter; ++i)
154  {
155  if (parameters_io.find(i) == parameters_io.end())
156  {
157  fits.SetParameter(i, default_values[i].def);
158 
159  if (default_values[i].min < default_values[i].max)
160  {
161  fits.SetParLimits(i, default_values[i].min, default_values[i].max);
162  }
163  else
164  {
165  fits.FixParameter(i, default_values[i].def);
166  }
167 
168  if (verbosity)
169  {
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;
174  }
175  }
176  else
177  {
178 // fits.SetParLimits(i, parameters_io[i], parameters_io[i]);
179  fits.SetParameter(i, parameters_io[i]);
180  fits.FixParameter(i, parameters_io[i]);
181 
182  if (verbosity)
183  {
184  cout << "SampleFit_PowerLawDoubleExp - parameter [" << i << "]: fixed to " << parameters_io[i] << endl;
185  }
186  }
187  }
188 
189  if (verbosity <= 1)
190  gpulse.Fit(&fits, "QRN0W", "goff", 0., (double) n_samples);
191  else
192  gpulse.Fit(&fits, "RN0VW+", "goff", 0., (double) n_samples);
193 
194  // store results
195  pedestal = fits.GetParameter(4);
196 
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;
201 
202  if (fits.GetParameter(0) > 0)
203  peak_sample = fits.GetMaximumX(fits.GetParameter(1), max_peakpos);
204  else
205  peak_sample = fits.GetMinimumX(fits.GetParameter(1), max_peakpos);
206 
207  peak = fits.Eval(peak_sample) - pedestal;
208 
209  if (verbosity)
210  {
211  static int id = 0;
212  ++id;
213 
214  string c_name(string("SampleFit_PowerLawDoubleExp_") + to_string(id));
215 
216  TCanvas *canvas = new TCanvas(
217  c_name.c_str(), c_name.c_str());
218  canvas->Update();
219 
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());
222 
223  fits.SetLineColor(kMagenta);
224  fits.DrawClone("same");
225  fits.Print();
226 
227  TF1 f1("f_SignalShape_PowerLawExp1", SignalShape_PowerLawExp, 0., n_samples, 5);
228  f1.SetParameters(
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");
236 
237  TF1 f2("f_SignalShape_PowerLawExp2", SignalShape_PowerLawExp, 0., n_samples, 5);
238  f2.SetParameters(
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");
246 
247  TGraph g_max(1);
248 
249  g_max.GetX()[0] = peak_sample;
250  g_max.GetY()[0] = peak + pedestal;
251 
252  g_max.SetMarkerStyle(kFullCircle);
253  g_max.SetMarkerSize(2);
254  g_max.SetMarkerColor(kRed);
255 
256  g_max.DrawClone("p");
257 
258  canvas->Update();
259 
260  // if (id == 1)
261  // {
262  // canvas->Print("SampleFit_PowerLawDoubleExp.pdf(");
263  // }
264  canvas->Print("SampleFit_PowerLawDoubleExp.pdf");
265  }
266 
267  for (int i = 0; i < n_parameter; ++i)
268  {
269  parameters_io[i] = fits.GetParameter(i);
270  }
271 
272  if (verbosity)
273  {
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;
280  }
281 
282  return true;
283 }
284 
285 double
286 SignalShape_PowerLawExp(double *x, double *par)
287 {
288  double pedestal = par[4];
289  // + ((x[0] - 1.5 * par[1]) > 0) * par[5]; // quick fix on exting tails on the signal function
290  if (x[0] < par[1])
291  return pedestal;
292  //double signal = (-1)*par[0]*pow((x[0]-par[1]),par[2])*exp(-(x[0]-par[1])*par[3]);
293  double signal = par[0] * pow((x[0] - par[1]), par[2]) * exp(-(x[0] - par[1]) * par[3]);
294  return pedestal + signal;
295 }
296 
297 double
298 SignalShape_PowerLawDoubleExp(double *x, double *par)
299 {
300  double pedestal = par[4];
301  // + ((x[0] - 1.5 * par[1]) > 0) * par[5]; // quick fix on exting tails on the signal function
302  if (x[0] < par[1])
303  return pedestal;
304  //double signal = (-1)*par[0]*pow((x[0]-par[1]),par[2])*exp(-(x[0]-par[1])*par[3]);
305  // 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)
306  // fits.GetParameter(2) / peak_shift = fits.GetParameter(3); // signal peak time
307 
308  double signal = //
309  par[0] //
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])) //
313  );
314  return pedestal + signal;
315 }
316 
317 
318 } // namespace FEEv2
319 
320 } // namespace TpcPrototypeDefs