Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-06 08:17:32 

0001 #include "CaloWaveformFitting.h"
0002 
0003 #include <TF1.h>
0004 #include <TFile.h>
0005 #include <TProfile.h>
0006 #include <TSpline.h>
0007 
0008 #include <Fit/BinData.h>
0009 #include <Fit/Chi2FCN.h>
0010 #include <Fit/Fitter.h>
0011 #include <Fit/UnBinData.h>
0012 #include <HFitInterface.h>
0013 #include <Math/WrappedMultiTF1.h>
0014 #include <Math/WrappedTF1.h>
0015 #include <ROOT/TThreadExecutor.hxx>
0016 #include <ROOT/TThreadedObject.hxx>
0017 
0018 #include <pthread.h>
0019 #include <algorithm>
0020 #include <iostream>
0021 #include <limits>
0022 #include <string>
0023 
0024 static ROOT::TThreadExecutor *t = new ROOT::TThreadExecutor(1);// NOLINT(misc-use-anonymous-namespace)
0025 double CaloWaveformFitting::template_function(double *x, double *par)
0026 {
0027   Double_t v1 = (par[0] * h_template->Interpolate(x[0] - par[1])) + par[2];
0028   return v1;
0029 }
0030 
0031 CaloWaveformFitting::~CaloWaveformFitting()
0032 {
0033   delete h_template;
0034 }
0035 
0036 void CaloWaveformFitting::initialize_processing(const std::string &templatefile)
0037 {
0038   TFile *fin = TFile::Open(templatefile.c_str());
0039   assert(fin);
0040   assert(fin->IsOpen());
0041   h_template = dynamic_cast<TProfile *>(fin->Get("waveform_template"));
0042   h_template->SetDirectory(nullptr);
0043   fin->Close();
0044   delete fin;
0045   m_peakTimeTemp = h_template->GetBinCenter(h_template->GetMaximumBin());
0046   t = new ROOT::TThreadExecutor(_nthreads);
0047 }
0048 
0049 std::vector<std::vector<float>> CaloWaveformFitting::process_waveform(std::vector<std::vector<float>> waveformvector)
0050 {
0051   int size1 = waveformvector.size();
0052   std::vector<std::vector<float>> fitresults;
0053   for (int i = 0; i < size1; i++)
0054   {
0055     waveformvector.at(i).push_back(i);
0056   }
0057   fitresults = calo_processing_templatefit(waveformvector);
0058   return fitresults;
0059 }
0060 
0061 std::vector<std::vector<float>> CaloWaveformFitting::calo_processing_templatefit(std::vector<std::vector<float>> chnlvector)
0062 {
0063   auto func = [&](std::vector<float> &v)
0064   {
0065     int size1 = v.size() - 1;
0066     if (size1 == _nzerosuppresssamples)
0067     {
0068       v.push_back(v.at(1) - v.at(0));                        // returns peak sample - pedestal sample
0069       v.push_back(std::numeric_limits<float>::quiet_NaN());  // set time to qnan for ZS
0070       v.push_back(v.at(0));
0071       if (v.at(0) != 0 && v.at(1) == 0)  // check if post-sample is 0, if so set high chi2
0072       {
0073         v.push_back(1000000);
0074       }
0075       else
0076       {
0077         v.push_back(std::numeric_limits<float>::quiet_NaN());
0078       }
0079       v.push_back(0);
0080     }
0081     else
0082     {
0083       float maxheight = 0;
0084       int maxbin = 0;
0085       for (int i = 0; i < size1; i++)
0086       {
0087         if (v.at(i) > maxheight)
0088         {
0089           maxheight = v.at(i);
0090           maxbin = i;
0091         }
0092       }
0093       float pedestal = 1500;
0094       if (maxbin > 4)
0095       {
0096         pedestal = 0.5 * (v.at(maxbin - 4) + v.at(maxbin - 5));
0097       }
0098       else if (maxbin > 3)
0099       {
0100         pedestal = (v.at(maxbin - 4));
0101       }
0102       else
0103       {
0104         pedestal = 0.5 * (v.at(size1 - 3) + v.at(size1 - 2));
0105       }
0106 
0107       if ((_bdosoftwarezerosuppression && v.at(6) - v.at(0) < _nsoftwarezerosuppression) || (_maxsoftwarezerosuppression && maxheight - pedestal < _nsoftwarezerosuppression))
0108       {
0109         v.push_back(v.at(6) - v.at(0));
0110         v.push_back(std::numeric_limits<float>::quiet_NaN());
0111         v.push_back(v.at(0));
0112         if (v.at(0) != 0 && v.at(1) == 0)  // check if post-sample is 0, if so set high chi2
0113         {
0114           v.push_back(1000000);
0115         }
0116         else
0117         {
0118           v.push_back(std::numeric_limits<float>::quiet_NaN());
0119         }
0120         v.push_back(0);
0121       }
0122       else
0123       {
0124         auto *h = new TH1F(std::string("h_" + std::to_string((int) round(v.at(size1)))).c_str(), "", size1, -0.5, size1 - 0.5);
0125 
0126         int ndata = 0;
0127         for (int i = 0; i < size1; ++i)
0128         {
0129           if ((v.at(i) == 16383) && _handleSaturation)
0130           {
0131             continue;
0132           }
0133 
0134           h->SetBinContent(i + 1, v.at(i));
0135           h->SetBinError(i + 1, 1);
0136           ndata++;
0137         }
0138         // if too many are saturated don't do the saturation recovery need enough ndf
0139         if (ndata > (size1 - 4))
0140         {
0141           ndata = size1;
0142           for (int i = 0; i < size1; ++i)
0143           {
0144             h->SetBinContent(i + 1, v.at(i));
0145             h->SetBinError(i + 1, 1);
0146           }
0147         }
0148 
0149         auto *f = new TF1(std::string("f_" + std::to_string((int) round(v.at(size1)))).c_str(), this, &CaloWaveformFitting::template_function, 0, 31, 3, "CaloWaveformFitting", "template_function");
0150         ROOT::Math::WrappedMultiTF1 *fitFunction = new ROOT::Math::WrappedMultiTF1(*f, 3);
0151         ROOT::Fit::BinData data(v.size() - 1, 1);
0152         ROOT::Fit::FillData(data, h);
0153         ROOT::Fit::Chi2Function *EPChi2 = new ROOT::Fit::Chi2Function(data, *fitFunction);
0154         ROOT::Fit::Fitter *fitter = new ROOT::Fit::Fitter();
0155         fitter->Config().MinimizerOptions().SetMinimizerType("GSLMultiFit");
0156         fitter->Config().MinimizerOptions().SetPrintLevel(-1);
0157         double params[] = {static_cast<double>(maxheight - pedestal), static_cast<double>(maxbin - m_peakTimeTemp), static_cast<double>(pedestal)};
0158         // double params[] = {static_cast<double>(maxheight - pedestal), 0, static_cast<double>(pedestal)};
0159         fitter->Config().SetParamsSettings(3, params);
0160         fitter->Config().ParSettings(1).SetLimits(-1 * m_peakTimeTemp, size1 - m_peakTimeTemp);  // set lim on time par
0161         if (m_setTimeLim)
0162         {
0163           fitter->Config().ParSettings(1).SetLimits(m_timeLim_low, m_timeLim_high);
0164         }
0165         fitter->FitFCN(*EPChi2, nullptr, data.Size(), true);
0166         ROOT::Fit::FitResult fitres = fitter->Result();
0167         // get the result status
0168         /*
0169         bool validfit = fitres.IsValid();
0170         if(!validfit)
0171         {
0172           std::cout<<"invalid fit"<<std::endl;
0173           for (int i = 0; i < size1; ++i)
0174         {
0175           std::cout<<v.at(i)<<std::endl;
0176         }
0177         }
0178         */
0179         double chi2min = fitres.MinFcnValue();
0180         // chi2min /= size1 - 3;  // divide by the number of dof
0181         chi2min /= ndata - 3;  // divide by the number of dof
0182         if (chi2min > _chi2threshold && (f->GetParameter(2) < _bfr_highpedestalthreshold || pedestal < _bfr_highpedestalthreshold) && (f->GetParameter(2) > _bfr_lowpedestalthreshold || pedestal > _bfr_lowpedestalthreshold) && _dobitfliprecovery)
0183         {
0184           std::vector<float> rv;  // temporary recovered waveform
0185           rv.reserve(size1);
0186           for (int i = 0; i < size1; i++)
0187           {
0188             rv.push_back(v.at(i));
0189           }
0190           unsigned int bits[3] = {8192, 4096, 2048};
0191           for (auto bit : bits)
0192           {
0193             for (int i = 0; i < size1; i++)
0194             {
0195               if (((unsigned int) rv.at(i) & bit) && ((unsigned int) rv.at(i) % bit > _bfr_lowpedestalthreshold))
0196               {
0197                 rv.at(i) = rv.at(i) - bit;
0198               }
0199             }
0200           }
0201           for (int i = 0; i < size1; i++)
0202           {
0203             h->SetBinContent(i + 1, rv.at(i));
0204             h->SetBinError(i + 1, 1);
0205           }
0206 
0207           maxheight = 0;
0208           maxbin = 0;
0209           for (int i = 0; i < size1; i++)
0210           {
0211             if (rv.at(i) > maxheight)
0212             {
0213               maxheight = rv.at(i);
0214               maxbin = i;
0215             }
0216           }
0217           if (maxbin > 4)
0218           {
0219             pedestal = 0.5 * (rv.at(maxbin - 4) + rv.at(maxbin - 5));
0220           }
0221           else if (maxbin > 3)
0222           {
0223             pedestal = (rv.at(maxbin - 4));
0224           }
0225           else
0226           {
0227             pedestal = 0.5 * (rv.at(size1 - 3) + rv.at(size1 - 2));
0228           }
0229 
0230           auto *recover_f = new TF1(std::string("recover_f_" + std::to_string((int) round(v.at(size1)))).c_str(), this, &CaloWaveformFitting::template_function, 0, 31, 3, "CaloWaveformFitting", "template_function");
0231           ROOT::Math::WrappedMultiTF1 *recoverFitFunction = new ROOT::Math::WrappedMultiTF1(*recover_f, 3);
0232           ROOT::Fit::BinData recoverData(rv.size() - 1, 1);
0233           ROOT::Fit::FillData(recoverData, h);
0234           ROOT::Fit::Chi2Function *recoverEPChi2 = new ROOT::Fit::Chi2Function(recoverData, *recoverFitFunction);
0235           ROOT::Fit::Fitter *recoverFitter = new ROOT::Fit::Fitter();
0236           recoverFitter->Config().MinimizerOptions().SetMinimizerType("GSLMultiFit");
0237           double recover_params[] = {static_cast<double>(maxheight - pedestal), 0, static_cast<double>(pedestal)};
0238           recoverFitter->Config().SetParamsSettings(3, recover_params);
0239           recoverFitter->Config().ParSettings(1).SetLimits(-1 * m_peakTimeTemp, size1 - m_peakTimeTemp);  // set lim on time par
0240           recoverFitter->FitFCN(*recoverEPChi2, nullptr, recoverData.Size(), true);
0241           ROOT::Fit::FitResult recover_fitres = recoverFitter->Result();
0242           double recover_chi2min = recover_fitres.MinFcnValue();
0243           recover_chi2min /= size1 - 3;  // divide by the number of dof
0244           if (recover_chi2min < _chi2lowthreshold && recover_f->GetParameter(2) < _bfr_highpedestalthreshold && recover_f->GetParameter(2) > _bfr_lowpedestalthreshold)
0245           {
0246             for (int i = 0; i < size1; i++)
0247             {
0248               v.at(i) = rv.at(i);
0249             }
0250             for (int i = 0; i < 3; i++)
0251             {
0252               v.push_back(recover_f->GetParameter(i));
0253             }
0254             v.push_back(recover_chi2min);
0255             v.push_back(1);
0256           }
0257           else
0258           {
0259             for (int i = 0; i < 3; i++)
0260             {
0261               v.push_back(f->GetParameter(i));
0262             }
0263             v.push_back(chi2min);
0264             v.push_back(0);
0265           }
0266           recover_f->Delete();
0267           delete recoverFitFunction;
0268           delete recoverFitter;
0269           delete recoverEPChi2;
0270         }
0271         else
0272         {
0273           for (int i = 0; i < 3; i++)
0274           {
0275             v.push_back(f->GetParameter(i));
0276           }
0277           v.push_back(chi2min);
0278           v.push_back(0);
0279         }
0280         h->Delete();
0281         f->Delete();
0282         delete fitFunction;
0283         delete fitter;
0284         delete EPChi2;
0285       }
0286     }
0287   };
0288 
0289   t->Foreach(func, chnlvector);
0290   int size3 = chnlvector.size();
0291   std::vector<std::vector<float>> fit_params;
0292   std::vector<float> fit_params_tmp;
0293   for (int i = 0; i < size3; i++)
0294   {
0295     const std::vector<float> &tv = chnlvector.at(i);
0296     int size2 = tv.size();
0297     for (int q = 5; q > 0; q--)
0298     {
0299       fit_params_tmp.push_back(tv.at(size2 - q));
0300     }
0301     fit_params.push_back(fit_params_tmp);
0302     fit_params_tmp.clear();
0303   }
0304   chnlvector.clear();
0305   return fit_params;
0306 }
0307 
0308 void CaloWaveformFitting::FastMax(float x0, float x1, float x2, float y0, float y1, float y2, float &xmax, float &ymax)
0309 {
0310   int n = 3;
0311   double xp[3] = {x0, x1, x2};
0312   double yp[3] = {y0, y1, y2};
0313   TSpline3 *sp = new TSpline3("", xp, yp, n, "b2e2", 0, 0);
0314   double X;
0315   double Y;
0316   double B;
0317   double C;
0318   double D;
0319   ymax = y1;
0320   xmax = x1;
0321   if (y0 > ymax)
0322   {
0323     ymax = y0;
0324     xmax = x0;
0325   }
0326   if (y2 > ymax)
0327   {
0328     ymax = y2;
0329     xmax = x2;
0330   }
0331   for (int i = 0; i <= 1; i++)
0332   {
0333     sp->GetCoeff(i, X, Y, B, C, D);
0334     if (D == 0)
0335     {
0336       if (C < 0)
0337       {
0338         // TSpline is a quadratic equation
0339 
0340         float root = (-B / (2 * C)) + X;
0341         if (root >= xp[i] && root <= xp[i + 1])
0342         {
0343           float yvalue = sp->Eval(root);
0344           if (yvalue > ymax)
0345           {
0346             ymax = yvalue;
0347             xmax = root;
0348           }
0349         }
0350       }
0351     }
0352     else
0353     {
0354       // find x when derivative = 0
0355       float root = ((-2 * C + sqrt((4 * C * C) - (12 * B * D))) / (6 * D)) + X;
0356       if (root >= xp[i] && root <= xp[i + 1])
0357       {
0358         float yvalue = sp->Eval(root);
0359         if (yvalue > ymax)
0360         {
0361           ymax = yvalue;
0362           xmax = root;
0363         }
0364       }
0365       root = (-2 * C - sqrt((4 * C * C) - (12 * B * D))) / (6 * D) + X;
0366       if (root >= xp[i] && root <= xp[i + 1])
0367       {
0368         float yvalue = sp->Eval(root);
0369         if (yvalue > ymax)
0370         {
0371           ymax = yvalue;
0372           xmax = root;
0373         }
0374       }
0375     }
0376   }
0377   delete sp;
0378   return;
0379 }
0380 std::vector<std::vector<float>> CaloWaveformFitting::calo_processing_fast(const std::vector<std::vector<float>> &chnlvector)
0381 {
0382   std::vector<std::vector<float>> fit_values;
0383   int nchnls = chnlvector.size();
0384   for (int m = 0; m < nchnls; m++)
0385   {
0386     const std::vector<float> &v = chnlvector.at(m);
0387     int nsamples = v.size();
0388 
0389     double maxy = v.at(0);
0390     float amp = 0;
0391     float time = 0;
0392     float ped = 0;
0393     float chi2 = std::numeric_limits<float>::quiet_NaN();
0394     if (nsamples == 2)
0395     {
0396       amp = v.at(1);
0397       time = std::numeric_limits<float>::quiet_NaN();
0398       ped = v.at(0);
0399       if (v.at(0) != 0 && v.at(1) == 0)  // check if post-sample is 0, if so set high chi2
0400       {
0401         chi2 = 1000000;
0402       }
0403     }
0404     else if (nsamples >= 3)
0405     {
0406       int maxx = 0;
0407       for (int i = 0; i < nsamples; i++)
0408       {
0409         if (i < 3)
0410         {
0411           ped += v.at(i);
0412         }
0413         if (v.at(i) > maxy)
0414         {
0415           maxy = v.at(i);
0416           maxx = i;
0417         }
0418       }
0419       ped /= 3;
0420       // if maxx <=5 nsample >=10 use the last two sample for pedestal(for HCal TP)
0421       if (maxx <= 5 && nsamples >= 10)
0422       {
0423         ped = 0.5 * (v.at(nsamples - 2) + v.at(nsamples - 1));
0424       }
0425       if (maxx == 0 || maxx == nsamples - 1)
0426       {
0427         amp = maxy;
0428         time = maxx;
0429       }
0430       else
0431       {
0432         FastMax(maxx - 1, maxx, maxx + 1, v.at(maxx - 1), v.at(maxx), v.at(maxx + 1), time, amp);
0433       }
0434     }
0435     amp -= ped;
0436     std::vector<float> val = {amp, time, ped, chi2, 0};
0437     fit_values.push_back(val);
0438     val.clear();
0439   }
0440   return fit_values;
0441 }
0442 
0443 std::vector<std::vector<float>> CaloWaveformFitting::calo_processing_nyquist(const std::vector<std::vector<float>> &chnlvector)
0444 {
0445   std::vector<std::vector<float>> fit_values;
0446   int nchnls = chnlvector.size();
0447   for (int m = 0; m < nchnls; m++)
0448   {
0449     std::vector<float> v = chnlvector.at(m);
0450     int nsamples = (int) v.size();
0451 
0452     if (nsamples == 2)
0453     {
0454       float chi2 = std::numeric_limits<float>::quiet_NaN();
0455       if (v.at(0) != 0 && v.at(1) == 0)  // check if post-sample is 0, if so set high chi2
0456       {
0457         chi2 = 1000000;
0458       }
0459       fit_values.push_back({v.at(1) - v.at(0), std::numeric_limits<float>::quiet_NaN(), v.at(0), chi2, 0});
0460       continue;
0461     }
0462 
0463     std::vector<float> result = NyquistInterpolation(v);
0464     fit_values.push_back(result);
0465   }
0466   return fit_values;
0467 }
0468 // mabye I can find a way to make it thread safe
0469 std::vector<float> CaloWaveformFitting::NyquistInterpolation(std::vector<float> &vec_signal_samples)
0470 {
0471   // int N = (int) vec_signal_samples.size();
0472   auto max_elem_iter = std::max_element(vec_signal_samples.begin(), vec_signal_samples.end());
0473   int maxx = std::distance(vec_signal_samples.begin(), max_elem_iter);
0474   float max = *max_elem_iter;
0475 
0476   float maxpos = maxx;
0477   float steplength = 0.5;
0478 
0479   while (steplength > 0.001)
0480   {
0481     // use 1.5 instead of 1 to avoid the floating point error...
0482     float starttime = maxpos - (1 * steplength);
0483     float endtime = maxpos + (1.5 * steplength);
0484 
0485     for (float i = starttime; i < endtime; i += steplength)
0486     {
0487       float yval = max;
0488       if (i != maxpos)
0489       {
0490         yval = psinc(i, vec_signal_samples);
0491       }
0492       if (yval > max)
0493       {
0494         max = yval;
0495         maxpos = i;
0496       }
0497     }
0498     steplength /= 2;
0499   }
0500 
0501   float pedestal = 0;
0502 
0503   if (maxpos > 5)
0504   {
0505     for (int i = 0; i < 3; i++)
0506     {
0507       pedestal += vec_signal_samples[i];
0508     }
0509     pedestal = pedestal / 3;
0510   }
0511   else if (maxpos > 4)
0512   {
0513     pedestal = (vec_signal_samples[0] + vec_signal_samples[1]) / 2;
0514   }
0515   // need more consideration for what is the most effieicnt
0516   else
0517   {
0518     pedestal = max;
0519     for (float i = maxpos - 5; i < maxpos; i += 0.1)
0520     {
0521       float yval = psinc(i, vec_signal_samples);
0522       pedestal = std::min(yval, pedestal);
0523     }
0524   }
0525   // calculate chi2 using the tempalte
0526   float chi2 = 0;
0527   double par[3] = {max - pedestal, maxpos - m_peakTimeTemp, pedestal};
0528   for (int i = 0; i < (int) vec_signal_samples.size(); i++)
0529   {
0530     double xval[1] = {(double) i};
0531     float diff = vec_signal_samples[i] - template_function(xval, par);
0532     chi2 += diff * diff;
0533   }
0534   std::vector<float> val = {max - pedestal, maxpos, pedestal, chi2, 0};
0535   return val;
0536 }
0537 
0538 // for odd N
0539 double CaloWaveformFitting::Dkernelodd(double x, int N)
0540 {
0541   double sum = 0;
0542   for (int k = 0; k < (N + 1) / 2; k++)
0543   {
0544     sum += 2 * std::cos(2 * M_PI * k * x / N);
0545   }
0546   sum -= 1;
0547   sum = sum / N;
0548   return sum;
0549 }
0550 
0551 // for even N
0552 double CaloWaveformFitting::Dkernel(double x, int N)
0553 {
0554   double sum = 0;
0555   for (int k = 0; k < N / 2; k++)
0556   {
0557     sum += 2 * std::cos(2 * M_PI * k * x / N);
0558   }
0559   sum -= 1;
0560   sum += std::cos(M_PI * x);
0561   sum = sum / N;
0562   return sum;
0563 }
0564 
0565 float CaloWaveformFitting::stablepsinc(float time, std::vector<float> &vec_signal_samples)
0566 {
0567   int N = (int) vec_signal_samples.size();
0568   float sum = 0;
0569   if (N % 2 == 0)
0570   {
0571     for (int n = 0; n < N; n++)
0572     {
0573       sum += vec_signal_samples[n] * Dkernel(time - n, N);
0574     }
0575   }
0576   else
0577   {
0578     for (int n = 0; n < N; n++)
0579     {
0580       sum += vec_signal_samples[n] * Dkernelodd(time - n, N);
0581     }
0582   }
0583   return sum;
0584 }
0585 
0586 float CaloWaveformFitting::psinc(float time, std::vector<float> &vec_signal_samples)
0587 {
0588   int N = (int) vec_signal_samples.size();
0589 
0590   if (std::abs(std::round(time) - time) < 1e-6)
0591   {
0592     if (time < 0 || time >= N)
0593     {
0594       return stablepsinc(time, vec_signal_samples);
0595     }
0596 
0597     return vec_signal_samples.at(std::round(time));
0598   }
0599 
0600   float sum = 0;
0601   if (N % 2 == 0)
0602   {
0603     for (int n = 0; n < N; n++)
0604     {
0605       double piu = M_PI * (time - n);
0606       double piuN = piu / N;
0607       sum += vec_signal_samples[n] * std::sin(piu) / (std::tan(piuN)) / N;
0608     }
0609   }
0610   else
0611   {
0612     for (int n = 0; n < N; n++)
0613     {
0614       double piu = M_PI * (time - n);
0615       double piuN = piu / N;
0616       sum += vec_signal_samples[n] * std::sin(piu) / (std::sin(piuN)) / N;
0617     }
0618   }
0619 
0620   return sum;
0621 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!