NewStructure/waveform_fitting/crystal_ball_fit.cc

0001 #include "waveform_fit_base.h"
0002 #include "crystal_ball_fit.h"
0003
0004 #include <map>
0005 #include <vector>
0006 #include <iostream>
0007
0008 #include <TF1.h>
0009 #include <TGraph.h>
0010 #include <TMath.h>
0011
0012 crystal_ball_fit::crystal_ball_fit() : waveform_fit_base{} {
0013     parameters[0] = 0;  // N
0014     parameters[1] = 0;  // alpha
0015     parameters[2] = 0;  // n
0016     parameters[3] = 0;  // mean
0017     parameters[4] = 0;  // sigma
0018     parameters[5] = 0;  // offset
0019
0020     crystal_ball_function = new TF1("crystal_ball_function", this, &crystal_ball_fit::crystal_ball, 0, 10, 6);
0021     linear_function = new TF1("constant", "[0]", 0, 10);
0022     crystal_ball_graph = new TGraph();
0023 }
0024
0025 crystal_ball_fit::~crystal_ball_fit() {
0026     delete crystal_ball_function;
0027     delete crystal_ball_graph;
0028 }
0029
0030 void crystal_ball_fit::fit() {
0031     // Fit the waveform with a crystal ball function
0032     // The crystal ball function is defined as:
0033     // f(x; alpha, n, mean, sigma) = N * exp(-0.5 * ((x - mean) / sigma)^2) for x < alpha
0034     //                                N * exp(-0.5 * ((x - mean) / sigma)^2) for x >= alpha
0035     //                                N * (n / abs(alpha)) * (n / abs(alpha)) * exp(-0.5 * alpha * alpha) * (n / abs(alpha) - abs(alpha) - (x - mean) / sigma)^(-n) for x < alpha
0036     //                                N * exp(-0.5 * ((x - mean) / sigma)^2) for x >= alpha
0037     // where N is the normalization factor, alpha is the transition point, n is the power, mean is the mean, and sigma is the standard deviation
0038     // The parameters are stored in the parameters map with the following keys:
0039     // 0: N
0040     // 1: alpha
0041     // 2: n
0042     // 3: mean
0043     // 4: sigma
0044
0045     // Check if the waveform has been set
0046     if (waveform.size() == 0) {
0047         std::cerr << "Waveform has not been set" << std::endl;
0048         return;
0049     }
0050     crystal_ball_function->SetRange(0, waveform.size());
0051
0052     // Fill the graph
0053     for (int i = 0; i < waveform.size(); i++) {
0054         crystal_ball_graph->SetPoint(i, i, waveform[i]);
0055     }
0056
0057     // Try a linear fit
0058     linear_function->SetParameter(0, parameters[5]);
0059     if (parameters.count(15) && parameters.count(25)) {   // If the parameter has a lower and upper limit
0060         linear_function->SetParLimits(0, parameters[15], parameters[25]);
0061     }
0062     crystal_ball_graph->Fit(linear_function, "Q");
0063     if (linear_function->GetChisquare() < 500) {
0064         E = 0;
0065         fit_ndf = linear_function->GetNDF();
0066         fit_chi2 = linear_function->GetChisquare();
0067         stale = false;
0068         return;
0069     }
0070
0071     // Set initial parameters
0072     for (int i = 0; i < 6; i++) {
0073         crystal_ball_function->SetParameter(i, parameters[i]);
0074         if (parameters.count(i + 10) && parameters.count(i + 20)) {   // If the parameter has a lower and upper limit
0075             crystal_ball_function->SetParLimits(i, parameters[i + 10], parameters[i + 20]);
0076         }
0077     }
0078     // Fit the graph
0079     crystal_ball_graph->Fit(crystal_ball_function, "Q");
0080
0081     // Get the results
0082     E = crystal_ball_function->GetParameter(4);
0083     fit_ndf = crystal_ball_function->GetNDF();
0084     fit_chi2 = crystal_ball_function->GetChisquare();
0085
0086     stale = false;
0087 }
0088
0089 int crystal_ball_fit::get_pedestal() {
0090     if (stale) {
0091         return -1;
0092     }
0093     return crystal_ball_function->GetParameter(5);
0094 }
0095
0096 double crystal_ball_fit::crystal_ball(double *inputs, double *parameters) {
0097 // Parameters
0098     // alpha: Where the gaussian transitions to the power law tail - fix?
0099     // n: The exponent of the power law tail - fix?
0100     // x_bar: The mean of the gaussian - free
0101     // sigma: The width of the gaussian - fix ?
0102     // N: The normalization of the gaussian - free
0103     // B baseline - fix?
0104
0105     double x = inputs[0];
0106
0107     double alpha = parameters[0];
0108     double n = parameters[1];
0109     double x_bar = parameters[2];
0110     double sigma = parameters[3];
0111     double N = parameters[4];
0112     double offset = parameters[5];
0113
0114     double A = pow(n / fabs(alpha), n) * exp(-0.5 * alpha * alpha);
0115     double B = n / fabs(alpha) - fabs(alpha);
0116     // std::cout << "A: " << A << std::endl;
0117
0118     // std::cout << "alpha: " << alpha << " n: " << n << " x_bar: " << x_bar << " sigma: " << sigma << " N: " << N << " B: " << B << " A: " << A << std::endl;
0119
0120     double ret_val;
0121     if ((x - x_bar) / sigma < alpha) {
0122         // std::cout << "path a" << std::endl;
0123         ret_val = exp(-0.5 * (x - x_bar) * (x - x_bar) / (sigma * sigma));
0124     } else {
0125         // std::cout << "path b" << std::endl;
0126         ret_val = A * pow(B + (x - x_bar) / sigma, -1 * n);
0127     }
0128     ret_val = N * ret_val + offset;
0129     // std::cout << "x: " << x << " y: " << ret_val << std::endl;
0130     return ret_val;
0131 }