From 9b0a4f7c63226e92fe68f9be9875d2fe488b76a2 Mon Sep 17 00:00:00 2001 From: JackHLindon Date: Tue, 14 Jun 2022 12:07:56 +0100 Subject: [PATCH 1/3] Add Relativistic Breit Wigner Function to TMath Implemented relativistic version of breit wigner (non-relativistic case already exists in TMath). Define BreitWignerRelativistic(Double_t x, Double_t median=0, Double_t gamma=1) in TMath.h and then calculate in TMath.cxx, a similar non relativistic function BreitWigner already exists in the same location which was used as a template. A tutorial BreitWigner.C has been added in tutorials/math which produces plots comparing the non relativistic and relativistic case --- math/mathcore/inc/TMath.h | 1 + math/mathcore/src/TMath.cxx | 18 ++++ math/mathcore/test/testTMath.cxx | 18 ++++ tutorials/math/BreitWigner.C | 173 +++++++++++++++++++++++++++++++ 4 files changed, 210 insertions(+) create mode 100644 tutorials/math/BreitWigner.C diff --git a/math/mathcore/inc/TMath.h b/math/mathcore/inc/TMath.h index 17ebef994a5b9..70ebdaadd68cb 100644 --- a/math/mathcore/inc/TMath.h +++ b/math/mathcore/inc/TMath.h @@ -482,6 +482,7 @@ struct Limits { Double_t Binomial(Int_t n,Int_t k); // Calculate the binomial coefficient n over k Double_t BinomialI(Double_t p, Int_t n, Int_t k); Double_t BreitWigner(Double_t x, Double_t mean=0, Double_t gamma=1); + Double_t BreitWignerRelativistic(Double_t x, Double_t median=0, Double_t gamma=1); Double_t CauchyDist(Double_t x, Double_t t=0, Double_t s=1); Double_t ChisquareQuantile(Double_t p, Double_t ndf); Double_t FDist(Double_t F, Double_t N, Double_t M); diff --git a/math/mathcore/src/TMath.cxx b/math/mathcore/src/TMath.cxx index a6a3cafb8ec6d..6d3d3f86cb3f5 100644 --- a/math/mathcore/src/TMath.cxx +++ b/math/mathcore/src/TMath.cxx @@ -445,6 +445,24 @@ Double_t TMath::BreitWigner(Double_t x, Double_t mean, Double_t gamma) return bw/(2*Pi()); } +//////////////////////////////////////////////////////////////////////////////// +/// Calculates a Relativistic Breit Wigner function with median and gamma. +// \f$ BW(E) = \frac{2\sqrt{2}}{\pi}\frac{M^{2}\gamma\sqrt{M^{2} + \gamma^{2}}}{\left(\sqrt{M^{2}+M\sqrt{M^{2} + \gamma^{2}}}\right)\left(\left(E^{2} - M^{2}\right)^{2} + M^{2}\gamma^{2}\right)} \f$ + +Double_t TMath::BreitWignerRelativistic(Double_t x, Double_t median, Double_t gamma) +{ + Double_t mm = median*median; + Double_t gg = gamma*gamma; + Double_t mg = median*gamma; + Double_t xxMinusmm = x*x - mm; + + Double_t y = sqrt(mm * (mm + gg)); + Double_t k = (0.90031631615710606*mg*y)/(sqrt(mm+y)); //2*sqrt(2)/pi = 0.90031631615710606 + + Double_t bw = k/(xxMinusmm*xxMinusmm + mg*mg); + return bw; +} + //////////////////////////////////////////////////////////////////////////////// /// Calculates a gaussian function with mean and sigma. /// If norm=kTRUE (default is kFALSE) the result is divided diff --git a/math/mathcore/test/testTMath.cxx b/math/mathcore/test/testTMath.cxx index ee0204b1325e9..6f96a240d68c9 100644 --- a/math/mathcore/test/testTMath.cxx +++ b/math/mathcore/test/testTMath.cxx @@ -141,6 +141,20 @@ void testPlane() << endl; } +void testBreitWignerRelativistic() +{ + Double_t median = 5000; + Double_t gamma = 100; + Int_t nPoints = 10; + Double_t xMinimum = 0; Double_t xMaximum = 10000; + Double_t xStepSize = (xMaximum-xMinimum)/nPoints; + + for (Int_t i=0;i<=nPoints;i++) { + Double_t currentX = xMinimum+i*xStepSize; + cout << "BreitWignerRelativistic(" << currentX << "," << median << "," << gamma << ") = " << BreitWignerRelativistic(currentX,median,gamma) << endl; + } +} + void testTMath() { cout << "Starting tests on TMath..." << endl; @@ -179,6 +193,10 @@ void testTMath() testPlane(); testPlane(); + + cout << "\nBreitWignerRelativistic tests: " << endl; + + testBreitWignerRelativistic(); } int main() diff --git a/tutorials/math/BreitWigner.C b/tutorials/math/BreitWigner.C new file mode 100644 index 0000000000000..32cfaf5bffd10 --- /dev/null +++ b/tutorials/math/BreitWigner.C @@ -0,0 +1,173 @@ +/// \file +/// \ingroup tutorial_math +/// \notebook +/// Tutorial illustrating how to create a plot comparing a Breit Wigner to a Relativistic Breit Wigner +/// +/// can be run with: +/// +/// ~~~{.cpp} +/// root[0] .x BreitWigner.C +/// ~~~ +/// +/// \macro_image +/// \macro_code +/// +/// \author Jack Lindon + +#include "TMath.h" + +#include +#include +#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 pdfTitle){ + + /////////////////////////////////////////////////////// + //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 plotTitleXPos = 0.23; Double_t plotTitleYPos = 0.25; + Double_t fontSize = 0.04; + Double_t lineWidth = 2; + bool setLimitPlotLogScale = true; + std::string xAxisTitle = "E [GeV]"; std::string yAxisTitle = "Events"; + 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 three 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); + + + ///////////////////////////////////// + // Draw everything + grAxis->Draw("AL"); + gr1->Draw("L same"); + gr2->Draw("L same"); + legend->Draw(); + tex_Title->Draw(); + tex_Title2->Draw(); + tex_Title3->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 BreitWigner(){ + + ///////////////////////////////////////////////////////// + // 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_nonRelBW[nPoints], y_relBW[nPoints]; + + + ////////////////////////////////// + // Define Breit-Wigner 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::max(); + Double_t yMaximum = TMath::BreitWignerRelativistic(median,median,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 Date: Fri, 22 Jul 2022 16:37:31 +0100 Subject: [PATCH 2/3] Add Relativistic Voigt Function to mathmore and a tutorial for this function. --- math/mathmore/CMakeLists.txt | 2 + math/mathmore/inc/Math/LinkDef.h | 2 + math/mathmore/inc/Math/VoigtRelativistic.h | 63 +++++++ math/mathmore/src/VoigtRelativistic.cxx | 85 +++++++++ tutorials/math/Voigt.C | 196 +++++++++++++++++++++ 5 files changed, 348 insertions(+) create mode 100644 math/mathmore/inc/Math/VoigtRelativistic.h create mode 100644 math/mathmore/src/VoigtRelativistic.cxx create mode 100644 tutorials/math/Voigt.C diff --git a/math/mathmore/CMakeLists.txt b/math/mathmore/CMakeLists.txt index 8200bde048358..e6f0ac4792912 100644 --- a/math/mathmore/CMakeLists.txt +++ b/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 diff --git a/math/mathmore/inc/Math/LinkDef.h b/math/mathmore/inc/Math/LinkDef.h index 98e423721beeb..36107e0c8598e 100644 --- a/math/mathmore/inc/Math/LinkDef.h +++ b/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__ diff --git a/math/mathmore/inc/Math/VoigtRelativistic.h b/math/mathmore/inc/Math/VoigtRelativistic.h new file mode 100644 index 0000000000000..cda9c2433dc08 --- /dev/null +++ b/math/mathmore/inc/Math/VoigtRelativistic.h @@ -0,0 +1,63 @@ +// @(#)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 diff --git a/math/mathmore/src/VoigtRelativistic.cxx b/math/mathmore/src/VoigtRelativistic.cxx new file mode 100644 index 0000000000000..c06f1d3423eae --- /dev/null +++ b/math/mathmore/src/VoigtRelativistic.cxx @@ -0,0 +1,85 @@ +// @(#)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 diff --git a/tutorials/math/Voigt.C b/tutorials/math/Voigt.C new file mode 100644 index 0000000000000..81e90d9079f2d --- /dev/null +++ b/tutorials/math/Voigt.C @@ -0,0 +1,196 @@ +/// \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 +#include +#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::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 Date: Fri, 22 Jul 2022 16:37:31 +0100 Subject: [PATCH 3/3] Add Relativistic Voigt Function to mathmore and a tutorial for this function. --- math/mathmore/CMakeLists.txt | 2 + math/mathmore/inc/Math/LinkDef.h | 2 + math/mathmore/inc/Math/VoigtRelativistic.h | 65 +++++++ math/mathmore/src/VoigtRelativistic.cxx | 90 +++++++++ tutorials/math/Voigt.C | 216 +++++++++++++++++++++ 5 files changed, 375 insertions(+) create mode 100644 math/mathmore/inc/Math/VoigtRelativistic.h create mode 100644 math/mathmore/src/VoigtRelativistic.cxx create mode 100644 tutorials/math/Voigt.C diff --git a/math/mathmore/CMakeLists.txt b/math/mathmore/CMakeLists.txt index 8200bde048358..e6f0ac4792912 100644 --- a/math/mathmore/CMakeLists.txt +++ b/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 diff --git a/math/mathmore/inc/Math/LinkDef.h b/math/mathmore/inc/Math/LinkDef.h index 98e423721beeb..347e7a06d725b 100644 --- a/math/mathmore/inc/Math/LinkDef.h +++ b/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__ diff --git a/math/mathmore/inc/Math/VoigtRelativistic.h b/math/mathmore/inc/Math/VoigtRelativistic.h new file mode 100644 index 0000000000000..9aacc8fb06319 --- /dev/null +++ b/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 diff --git a/math/mathmore/src/VoigtRelativistic.cxx b/math/mathmore/src/VoigtRelativistic.cxx new file mode 100644 index 0000000000000..05b842a2aab55 --- /dev/null +++ b/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 diff --git a/tutorials/math/Voigt.C b/tutorials/math/Voigt.C new file mode 100644 index 0000000000000..33e405f23b957 --- /dev/null +++ b/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 +#include +#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::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); +}