Add Relativistic Voigt Function to TMath

math/mathmore/CMakeLists.txt

@@ -42,6 +42,7 @@ ROOT_STANDARD_LIBRARY_PACKAGE(MathMore
     Math/VavilovAccuratePdf.h
     Math/VavilovAccurateQuantile.h
     Math/VavilovFast.h
+    Math/VoigtRelativistic.h
 SOURCES
     src/ChebyshevApprox.cxx
     src/Derivator.cxx
@@ -77,6 +78,7 @@ SOURCES
     src/VavilovAccurateQuantile.cxx
     src/VavilovFast.cxx
     src/cblas.cxx
+    src/VoigtRelativistic.cxx
   LINKDEF
     Math/LinkDef.h
   DEPENDENCIES

math/mathmore/inc/Math/LinkDef.h

@@ -109,4 +109,6 @@
 #pragma link C++ class ROOT::Math::GSLMultiRootFinder+;
 #pragma link C++ typedef ROOT::Math::MultiRootFinder;
 
+#pragma link C++ class ROOT::Math::VoigtRelativistic + ;
+
 #endif //__CINT__

math/mathmore/inc/Math/VoigtRelativistic.h

@@ -0,0 +1,65 @@
+// @(#)root/mathmore:$Id$
+// Author: J. Lindon Wed Jun 15 02:35:26 2022
+
+/**********************************************************************
+ *                                                                    *
+ * Copyright (c) 2022  LCG ROOT Math Team, CERN/PH-SFT                *
+ *                                                                    *
+ * This library is free software; you can redistribute it and/or      *
+ * modify it under the terms of the GNU General Public License        *
+ * as published by the Free Software Foundation; either version 2     *
+ * of the License, or (at your option) any later version.             *
+ *                                                                    *
+ * This library is distributed in the hope that it will be useful,    *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of     *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU   *
+ * General Public License for more details.                           *
+ *                                                                    *
+ * You should have received a copy of the GNU General Public License  *
+ * along with this library (see file COPYING); if not, write          *
+ * to the Free Software Foundation, Inc., 59 Temple Place, Suite      *
+ * 330, Boston, MA 02111-1307 USA, or contact the author.             *
+ *                                                                    *
+ **********************************************************************/
+
+// Header file for class VoigtRelativistic
+
+#ifndef ROOT_Math_VoigtRelativistic
+#define ROOT_Math_VoigtRelativistic
+
+//////////////////////////////////////////////////////////////////////////////
+//                                                                          //
+// VoigtRelativistic                                                        //
+//                                                                          //
+// Calculates the relativistic voigt function                               //
+// (convolution of a relativistic breit wigner and a Gaussian)              //
+// Implemented as equation 9 of https://doi.org/10.1016/j.jmaa.2018.03.065  //
+// With normalization changed to match TMath::Voigt (non relativistic voigt)//
+//                                                                          //
+//////////////////////////////////////////////////////////////////////////////
+
+namespace ROOT {
+namespace Math {
+
+class VoigtRelativistic {
+public:
+   // The relativistic voigt function
+   static double
+   evaluate(double x, double median, double sigma, double lg,
+            double integrationRange =
+               26.615717509251260); // t=26.615717509251260 gives exp(-t^2)=minimum value stored by double.
+   static double
+   dumpingFunction(double median, double sigma, double lg,
+                   double integrationRange =
+                      26.615717509251260); // t=26.615717509251260 gives exp(-t^2)=minimum value stored by double.
+
+   // Include an empty virtual desctructor to eliminate compiler warnings
+   virtual ~VoigtRelativistic() {}
+
+protected:
+};
+
+} // namespace Math
+} // namespace ROOT
+
+#endif

math/mathmore/src/VoigtRelativistic.cxx

@@ -0,0 +1,90 @@
+// @(#)root/mathmore:$Id$
+// Author: J. Lindon Wed Jun 15 02:35:26 2022
+
+/**********************************************************************
+ *                                                                    *
+ * Copyright (c) 2022  LCG ROOT Math Team, CERN/PH-SFT                *
+ *                                                                    *
+ *                                                                    *
+ **********************************************************************/
+
+// Implementation file for class VoigtRelativistic
+
+#include "Math/VoigtRelativistic.h"
+#include "Math/Integrator.h"
+#include "Math/Functor.h"
+#include "TMath.h"
+
+namespace ROOT {
+namespace Math {
+
+////////////////////////////////////////////////////////////////////////////////
+// Computation of a relativistic voigt function (normalised), the convolution of
+// a guassian and relativistic breit wigner function.
+//
+// \f$ V(E;\mu,\sigma,\Gamma) = \left(\frac{\sqrt{2}}{2\pi^{\frac{3}{2}}}\right)\left(\frac{\Gamma\mu^{2}
+// \sqrt{\mu^{2}+\Gamma^{2} }}{\sigma^{4}\sqrt{\mu^{2}+\mu\sqrt{\mu^{2}+\Gamma^{2} }}}\right) \int_{-\infty}^{+\infty}
+// \frac{e^{-t^{2}}}{\left(\frac{1}{\sqrt{2}}\left(E-\mu\right)\sigma-t\right)^{2}\left(\frac{1}{\sqrt{2}}\left(E+\mu\right)\sigma-t\right)^{2}+\frac{\Gamma^{2}\mu^{2}}{4\sigma^{4}}}
+// \f$ Where all varialbles are real.
+//
+// E is the independent variable (typically the energy)
+// \f$\mu\f$ is the median (typically the pole mass of the resonance)
+// \f$\sigma\f$ is the width of the gaussian
+// \f$\Gamma\f$ is the width of the relativistic breit wigner
+// for this function to be exact the integrationRange must be taken to infinity
+// however, it converges rapidly and the integrationRange by default is taken
+// to the maximum precision allowed by this method for a double.
+////////////////////////////////////////////////////////////////////////////////
+
+double VoigtRelativistic::evaluate(double x, double median, double sigma, double lg, double integrationRange)
+{
+
+   double inverse_sSqrt2 = 1 / (1.41421356237309504 * sigma); // 1.41421356237309504=sqrt(2)
+   double ss = sigma * sigma;
+   double mm = median * median;
+
+   double u1 = (x - median) * inverse_sSqrt2;
+   double u2 = (x + median) * inverse_sSqrt2;
+   double a = (lg * median) / (2 * ss);
+   double aa = a * a;
+
+   double y = median * sqrt(mm + lg * lg);
+   double k = (0.25397454373696387 * a * y) / (ss * sqrt(mm + y)); // 0.25397454373696387=sqrt(2)/(pi^(3/2))
+
+   auto integrand = [&](double t) {
+      double u1Minust = (u1 - t);
+      double u2Minust = (u2 - t);
+      return ((exp(-t * t)) / (u1Minust * u1Minust * u2Minust * u2Minust + (aa)));
+   };
+
+   ROOT::Math::Functor1D f(integrand);
+   ROOT::Math::Integrator integrator(ROOT::Math::IntegrationOneDim::kDEFAULT);
+   integrator.SetFunction(f);
+   double voigt = k * integrator.Integral(-integrationRange, integrationRange);
+
+   return voigt;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+// Computation of a relativistic voigt function's dumping function, the ratio of
+// the peak of the voigt to the peak of the breit wigner that makes it up
+//
+// \f$ D(\sigma;\Gamma,\mu) = \frac{V(\mu;\mu,\sigma,\Gamma)}{BW(\mu;\mu,\Gamma)} \f$
+// Where V(\mu;\mu,\sigma,\Gamma) is the value of the relativistic voigt at the median
+// and   BW(\mu;\mu,\Gamma) is the value of the relativistic breit wigner at the median
+// where all varialbles are real.
+// \f$\mu\f$ is the median (typically the pole mass of the resonance).
+// \f$\sigma\f$ is the width of the gaussian
+// \f$\Gamma\f$ is the width of the relativistic breit wigner
+// for this function to be exact the integrationRange must be taken to infinity
+// however, it converges rapidly and the integrationRange by default is taken
+// to the maximum precision allowed by this method for a double.
+////////////////////////////////////////////////////////////////////////////////
+double VoigtRelativistic::dumpingFunction(double median, double sigma, double lg, double integrationRange)
+{
+   return VoigtRelativistic::evaluate(median, median, sigma, lg, integrationRange) /
+          TMath::BreitWignerRelativistic(median, median, lg);
+}
+
+} // namespace Math
+} // namespace ROOT

tutorials/math/Voigt.C

@@ -0,0 +1,216 @@
+/// \file
+/// \ingroup tutorial_math
+/// \notebook
+/// Tutorial illustrating how to create a plot comparing a Voigt to a Relativistic Voigt
+///
+/// can be run with:
+///
+/// ~~~{.cpp}
+///  root[0] .x Voigt.C
+/// ~~~
+///
+/// \macro_image
+/// \macro_code
+///
+/// \author Jack Lindon
+
+#include "TMath.h"
+#include "Math/VoigtRelativistic.h"
+
+#include <limits>
+#include <string>
+#include "TAxis.h"
+#include "TGraph.h"
+#include "TCanvas.h"
+#include "TLatex.h"
+#include "TLegend.h"
+#include "TStyle.h" //For gStyle to remove stat box.
+
+void plotTwoTGraphs(Double_t x[], Double_t y1[], Double_t y2[], const Int_t nPoints, Double_t lowerXLimit,
+                    Double_t upperXLimit, Double_t lowerYLimit, Double_t upperYLimit, std::string legend1,
+                    std::string legend2, std::string plotTitle1, std::string plotTitle2, std::string plotTitle3,
+                    std::string plotTitle4, std::string pdfTitle, std::string xAxisTitle = "E [GeV]",
+                    std::string yAxisTitle = "Events", bool setLimitPlotLogScale = true, Double_t plotTitleXPos = 0.23,
+                    Double_t plotTitleYPos = 0.25)
+{
+
+   ///////////////////////////////////////////////////////
+   // Define variables for plot aesthetics and positioning
+   Double_t legendXPos = 0.63;
+   Double_t legendYPos = 0.85;
+   Double_t legendXWidth = 0.29;
+   Double_t legendYHeight = 0.1;
+   Double_t fontSize = 0.04;
+   Double_t lineWidth = 2;
+   Double_t xAxisTitleOffset = 1;
+   Double_t yAxisTitleOffset = 1.3;
+   gStyle->SetOptStat(0);
+
+   ///////////////////////////////////////////////////////
+   // Initialize TGraphs
+   TGraph *gr1 = new TGraph(nPoints, x, y1);
+   TGraph *gr2 = new TGraph(nPoints, x, y2);
+   gr1->SetLineWidth(lineWidth);
+   gr2->SetLineWidth(lineWidth);
+   gr1->SetLineColor(kBlack);
+   gr2->SetLineColor(kBlue);
+
+   /////////////////////////////////////////////////////////
+   // Initialize canvas
+   TCanvas *c1 = new TCanvas("c1", "transparent pad", 200, 10, 600, 600);
+   c1->SetLogy(setLimitPlotLogScale);
+   c1->SetTicks(1, 1);
+   c1->SetRightMargin(0.02);
+   c1->SetTopMargin(0.02);
+
+   ///////////////////////////////////////////////////////
+   // Make just a basic invisible TGraph just for the axes
+   const Double_t axis_x[2] = {lowerXLimit, upperXLimit};
+   const Double_t axis_y[2] = {lowerYLimit, upperYLimit};
+   TGraph *grAxis = new TGraph(2, axis_x, axis_y);
+   grAxis->SetTitle("");
+   grAxis->GetYaxis()->SetTitle(yAxisTitle.c_str());
+   grAxis->GetXaxis()->SetTitle(xAxisTitle.c_str());
+   grAxis->GetXaxis()->SetRangeUser(lowerXLimit, upperXLimit);
+   grAxis->GetYaxis()->SetRangeUser(lowerYLimit, upperYLimit);
+   grAxis->GetXaxis()->SetLabelSize(fontSize);
+   grAxis->GetYaxis()->SetLabelSize(fontSize);
+   grAxis->GetXaxis()->SetTitleSize(fontSize);
+   grAxis->GetYaxis()->SetTitleSize(fontSize);
+   grAxis->GetXaxis()->SetTitleOffset(xAxisTitleOffset);
+   grAxis->GetYaxis()->SetTitleOffset(yAxisTitleOffset);
+   grAxis->SetLineWidth(0); // So invisible
+
+   ///////////////////////////////////////////////////////////
+   // Make legend and set aesthetics
+   auto legend = new TLegend(legendXPos, legendYPos, legendXPos + legendXWidth, legendYPos + legendYHeight);
+   legend->SetFillStyle(0);
+   legend->SetBorderSize(0);
+   legend->SetTextSize(fontSize);
+   legend->AddEntry(gr1, legend1.c_str(), "L");
+   legend->AddEntry(gr2, legend2.c_str(), "L");
+
+   /////////////////////////////////////////////////////////////
+   // Add plot title to plot. Make in four lines so not crowded.
+   // Shift each line down by shiftY
+   float shiftY{0.037};
+   TLatex *tex_Title = new TLatex(plotTitleXPos, plotTitleYPos - 0 * shiftY, plotTitle1.c_str());
+   tex_Title->SetNDC();
+   tex_Title->SetTextFont(42);
+   tex_Title->SetTextSize(fontSize);
+   tex_Title->SetLineWidth(lineWidth);
+   TLatex *tex_Title2 = new TLatex(plotTitleXPos, plotTitleYPos - 1 * shiftY, plotTitle2.c_str());
+   tex_Title2->SetNDC();
+   tex_Title2->SetTextFont(42);
+   tex_Title2->SetTextSize(fontSize);
+   tex_Title2->SetLineWidth(lineWidth);
+   TLatex *tex_Title3 = new TLatex(plotTitleXPos, plotTitleYPos - 2 * shiftY, plotTitle3.c_str());
+   tex_Title3->SetNDC();
+   tex_Title3->SetTextFont(42);
+   tex_Title3->SetTextSize(fontSize);
+   tex_Title3->SetLineWidth(lineWidth);
+   TLatex *tex_Title4 = new TLatex(plotTitleXPos, plotTitleYPos - 3 * shiftY, plotTitle4.c_str());
+   tex_Title4->SetNDC();
+   tex_Title4->SetTextFont(42);
+   tex_Title4->SetTextSize(fontSize);
+   tex_Title4->SetLineWidth(lineWidth);
+
+   /////////////////////////////////////
+   // Draw everything
+   grAxis->Draw("AL");
+   gr1->Draw("L same");
+   gr2->Draw("L same");
+   legend->Draw();
+   tex_Title->Draw();
+   tex_Title2->Draw();
+   tex_Title3->Draw();
+   tex_Title4->Draw();
+   c1->RedrawAxis(); // Be sure to redraw axis AFTER plotting TGraphs otherwise TGraphs will be on top of tick marks and
+                     // axis borders.
+
+   gPad->Print(pdfTitle.c_str());
+}
+
+void Voigt()
+{
+
+   /////////////////////////////////////////////////////////
+   // Define x axis limits and steps for each plotted point
+   const Int_t nPoints = 1000;
+   Double_t xMinimum = 0;
+   Double_t xMaximum = 13000;
+   Double_t xStepSize = (xMaximum - xMinimum) / nPoints;
+
+   ///////////////////////////////////////////////////////
+   // Define arrays of (x,y) points.
+   Double_t x[nPoints];
+   Double_t y_nonRelVoigt[nPoints], y_relVoigt[nPoints];
+   Double_t y_nonRelVoigtDumpingFunction[nPoints], y_relVoigtDumpingFunction[nPoints];
+
+   //////////////////////////////////
+   // Define Voigt parameters
+   Double_t width = 1350;
+   Double_t sigma = 269.7899;
+   Double_t median = 9000;
+
+   ///////////////////////////////////////////////////
+   // Loop over x axis range, filling in (x,y) points,
+   // and finding y minimums and maximums for axis limit.
+   Double_t yMinimum = std::numeric_limits<Double_t>::max();
+   Double_t yMaximum = ROOT::Math::VoigtRelativistic::evaluate(
+      median, median, sigma, width); // y maximum is at x=median (and non relativistic = relativistic at median so
+                                     // choice of function does not matter).
+   for (Int_t i = 0; i < nPoints; i++) {
+      Double_t currentX = xMinimum + i * xStepSize;
+      x[i] = currentX;
+      y_nonRelVoigt[i] = TMath::Voigt(currentX - median, sigma, width);
+      y_relVoigt[i] = ROOT::Math::VoigtRelativistic::evaluate(currentX, median, sigma, width);
+
+      if (i != 0) { // calculate the voigt dumping functions with varying sigma. Skip sigma=0 as voigt is undefined
+                    // here.
+         y_nonRelVoigtDumpingFunction[i] =
+            TMath::Voigt(median - median, currentX, width) / TMath::BreitWigner(median, median, width);
+         y_relVoigtDumpingFunction[i] = ROOT::Math::VoigtRelativistic::dumpingFunction(median, currentX, width);
+         ;
+      }
+
+      if (y_nonRelVoigt[i] < yMinimum) {
+         yMinimum = y_nonRelVoigt[i];
+      }
+      if (y_relVoigt[i] < yMinimum) {
+         yMinimum = y_relVoigt[i];
+      }
+   }
+
+   plotTwoTGraphs(x, y_nonRelVoigt, y_relVoigt, nPoints, xMinimum, xMaximum // xAxis limits
+                  ,
+                  yMinimum / 4, yMaximum * 4 // yAxis limits, expand for aesthetics.
+                  ,
+                  "NonRel Voigt", "Rel Voigt" // Legend entries
+                  ,
+                  "Comparing Voigt", "M = " + std::to_string(int(round(median))) + " GeV",
+                  "#Gamma = " + std::to_string(int(round(width))) + " GeV",
+                  "#sigma = " + std::to_string(int(round(sigma))) + " GeV" // Plot Title entry four lines)
+                  ,
+                  "Voigt_M" + std::to_string(int(round(median))) + "_Gamma" + std::to_string(int(round(width))) +
+                     "_sigma" + std::to_string(int(round(sigma))) + ".pdf)" // PDF file title.
+   );
+
+   plotTwoTGraphs(x, y_nonRelVoigtDumpingFunction, y_relVoigtDumpingFunction, nPoints, 200, xMaximum // xAxis limits
+                  ,
+                  0, 1.2 // yAxis limits, expand for aesthetics.
+                  ,
+                  "NonRel Voigt DF", "Rel Voigt DF" // Legend entries
+                  ,
+                  "Voigt Dumping Functions", "M = " + std::to_string(int(round(median))) + " GeV",
+                  "#Gamma = " + std::to_string(int(round(width))) + " GeV", "" // Plot Title entry three lines)
+                  ,
+                  "VoigtDF_M" + std::to_string(int(round(median))) + "_Gamma" + std::to_string(int(round(width))) +
+                     ".pdf)" // PDF file title.
+                  ,
+                  "#sigma [GeV]", "DF" // axis titles
+                  ,
+                  false // no log
+                  ,
+                  0.23, 0.7);
+}