root icon indicating copy to clipboard operation
root copied to clipboard

Add Relativistic Voigt Function to TMath

Open JackHLindon opened this issue 2 years ago • 4 comments

This Pull request:

This follows on from a recent pull request which added a relativistic breit wigner to ROOT https://github.com/root-project/root/pull/10760, to add a relativistic voigt which is a natural extension.

Brief summary Files added: math/mathmore/inc/Math/VoigtRelativistic.h - Defines VoigtRelativistic class with two functions, evaluate which calculates the relativistic voigt, and dumpingFunction which calculates the relativisitic voigt's dumping function. math/mathmore/src/VoigtRelativistic.cxx - Implements the two functions defined in VoigtRelativistic.h . This implementation uses the redefinition of variables in https://arxiv.org/pdf/1711.09304.pdf which is required to write the integral in a computationally feasiable way. tutorials/math/Voigt.C - Add a tutorial which shows the use of the new VoigtRelativistic class two produce two plots, comparing relativistic voigt to non relativistic voigt, both as the voigt itself and their dumping functions. Files modified: math/mathmore/CMakeLists.txt math/mathmore/inc/Math/LinkDef.h - These two files needed to be modified to compile my changes adding in the relativisitc voigt.

Full description Implemented relativistic version of voigt (non-relativistic case already exists in TMath). As previously discussed in the ROOT forum https://root-forum.cern.ch/t/relativistic-breit-wigner-and-relativistic-voigt/48844 , the relativistic version has been added to mathmore rather than TMath, as it requires including a header file with an integrator which is CPU intensive. This has hence been added as it's own class in mathmore, with two functions:

static double evaluate(double x, double median, double sigma, double lg, double integrationRange=26.615717509251260); static double dumpingFunction(double median, double sigma, double lg, double integrationRange=26.615717509251260);

The default integrationRange was chosen as the integral has been wrote in a form where it is multiplied by exp(-t^2), and this value of integrationRange results in going to the point where exp(-t^2) reaches the minimum value stored by a double.

A relativistic voigt is a convolution of a relativistic breit wigner and a gaussian, doing this naively is computationally infeasible however a fairly simple redefinition of variables allows rewriting the integral in a form that is reasonable to compute. This is shown in https://arxiv.org/pdf/1711.09304.pdf equation 9 which has been used here. I am unsure how referencing works with ROOT. [Note: the function used here is not actually identical to equation 9, as the paper linked uses a different normalization than ROOT does for the non relativistic voigt which is already implemented in TMath. I have renormalized this equation to be consistent with the non relativistic case already in ROOT).

The two functions that have been defined are:

evaluate, simply gives the value of the relativisitic voigt at the point requested dumpingFunction, this is a commonly used function which is defined as the value of the relativistic voigt at its peak, divided by the value of the relativistic breit wigner at the peak that is part of the relativistic voigt. This is useful as it quantifies how much the voigt is smeared by the gaussian, which as a physical effect in particle physics is usually easily translated to how much sensitivity to a signal you lose due to detector effects.

A tutorial Voigt.C has been added in tutorials/math which produces plots comparing the non relativistic and relativistic case of voigt and the dumping function.

The output plots from tutorials/math/Voigt.C which produces a plot using the new VoigtRelativistic class which is added in this pull request is here Voigt DumpingFunction.

As mentioned this relativistic voigt is implemented in mathmore, however there is already a non relativistic voigt in TMath ( called Voigt https://root.cern.ch/doc/master/namespaceTMath.html ). I'm unsure if this is possible but it would be ideal if the documentation here could mention in the non relativistic voigt case that there is a relativistic version available as I think as is very few people even if they want a relativistic case would realise there is one now available (in fact I think many people would just assume the voigt in TMath is relativistic).

I have not added a test as I am unsure exactly how to and the non relativistic voigt that already exists does not have a test either to use as a template. I will try to learn how to add tests in the future and hopefully add one for this and the relativistic breit wigner case in the future.

Thank you for your help, Jack

Changes or fixes:

Checklist:

  • [x] tested changes locally
  • [x] updated the docs (if necessary)

JackHLindon avatar Jul 26 '22 00:07 JackHLindon

Can one of the admins verify this patch?

phsft-bot avatar Jul 26 '22 00:07 phsft-bot

Hi, I had to rebase to change the code format to pass the clang-tools code analysis check (is this new? I didn't have to do this in my previous pull request which had the same formatting as this PR which it is now unhappy with. Also weirdly one of the things it required was in LinkDef.h putting a space before and after "+" which isn't consistent with the other lines.

Also it required some very weird/bad formatting for calling the function "plotTwoTGraphs" in the tutorial added, I think it is trying to reduce the number of characters per line, but it does it in a quite poor way.

Also the command that the script "https://github.com/root-project/root/blob/master/.ci/format_script.sh" which does this format check says to run to rebase is incorrect. It says to do:

git rebase -i -x "git-clang-format-7 master && git commit -a --allow-empty --fixup=HEAD" --strategy-option=theirs origin/master
git rebase --autosquash -i master

But this does not run and complains that master does not exist. It should be

git rebase -i -x "git-clang-format-8 origin/master && git commit -a --allow-empty --fixup=HEAD" --strategy-option=theirs origin/master
git rebase --autosquash -i origin/master

[i.e. master-> origin/master] ).

This rebase added some spurious commit messages. When (/if) this pull request is accepted, could you please select "squash and merge", it's a lot simpler than me having to rebase and squashing manually in the terminal (which last time I tried I messed up so bad I ended up just having to delete my fork and start over), the title of the PR works as a commit message for the full thing "Add Relativistic Voigt Function to TMath"

Thanks in advance for any help, Jack

JackHLindon avatar Jul 26 '22 04:07 JackHLindon

Hi @lmoneta ,

Thank you for the help :). I'm not sure I exactly understand though, I thought I had implemented it as a class? I can't see much difference between this and the format of https://github.com/root-project/root/blob/master/math/mathmore/src/VavilovAccuratePdf.cxx (I assume this is the class you mean by VavilovPdf)

Sure I'm happy to add a test, though I'm not entirely sure exactly how, is following the format of the test for Vavilov ok? (e.g. https://github.com/root-project/root/blob/master/math/mathmore/test/VavilovTest.h https://github.com/root-project/root/blob/master/math/mathmore/test/VavilovTest.cxx )

This seems to just be a function that returns 0 if the values are as expected, is this what's wanted?

Cheers, Jack

JackHLindon avatar Aug 02 '22 02:08 JackHLindon

Hi, yes following the linked example of the VavilovAccuratePdf is fine. It is true you have made a class, but only with static functions. Like this is not needed to be a class, can be simple function in a namespace. The advantage of having as a class is that one can cache some information as class data members.

Then test of the Vavilov is also a good example, however it is better to use the google test framework, resulting in a much simpler test to write. You need just to use the macro defined by test to compare the function value with the reference one. You have many examples in ROOT for gtest, for example in hist/hist/test directory like test_TH1.cxx or test_TFormula.cxx, test_TF123_Moments.cxx. If you need any help for the test, please let me know

lmoneta avatar Aug 02 '22 08:08 lmoneta