Re: Wigner D-function in ROOT or C++?

From: Rene Brun <Rene.Brun_at_cern.ch>
Date: Tue, 04 Jul 2006 09:32:30 +0200


Art,

I have translated the function in CERNLIB to C++. See attachment. I just tested that it compiles and executes with some random input numbers. Could you test it and let me know if it is what you want?

Rene Brun

Arthur E. Snyder wrote:
> Hi again Rene,
>
> I see D-functions are not in mathlib yet though lots of other good stuff
> is. The so called beta term is available in cernlib (called RDJMNB) from
> the MATHLIB library. Is it possible to access this from ROOT?
>
> -Art S.
>
> A.E. Snyder, Group EC \!c*p?/
> SLAC Mail Stop #95 ((. .))
> Box 4349 |
> Stanford, Ca, USA, 94309 '\|/`
> e-mail:snyder_at_slac.stanford.edu o
> phone:650-926-2701 _
> http://www.slac.stanford.edu/~snyder BaBar
> FAX:650-926-2657 Collaboration
>
>
>
> On Mon, 3 Jul 2006, Arthur E. Snyder wrote:
>
>
>> Hi Rene,
>>
>> I'm not looking for Breit-Wigner but for Wigner's so called D-function
>> (D(alpha,beta,gamma)) which are functions of anlges used in partial wave
>> analysis based on the Jacob and Wick formalism.
>>
>> The article I refer to proposes somekind of mathlib to collect in C++ form
>> a lot of mathmatical functions available dispersed among a lot libraries
>> of functions like cernlib, clhep, napl, and so forth. The D-functions are
>> just mentioned in one of the tables .. google hunted them out .. with a
>> search for "Wigner D-function C++".
>>
>> Anyway I need the D-function in a form I can access from ROOT. I may just
>> have to code the few lower order ones I need myself, but I'd rather have a
>> nice (well debugged) set from some package or class library.
>>
>> -Art S.
>>
>> A.E. Snyder, Group EC \!c*p?/
>> SLAC Mail Stop #95 ((. .))
>> Box 4349 |
>> Stanford, Ca, USA, 94309 '\|/`
>> e-mail:snyder_at_slac.stanford.edu o
>> phone:650-926-2701 _
>> http://www.slac.stanford.edu/~snyder BaBar
>> FAX:650-926-2657 Collaboration
>>
>>
>>
>> On Mon, 3 Jul 2006, Rene Brun wrote:
>>
>>
>>> Art,
>>>
>>> I am curious to know the talk/article you are referring to ::)
>>> In ROOT we have the Breit-Wigner distributtion. See
>>> http://root.cern.ch/root/htmldoc/src/TMath.cxx.html#TMath:BreitWigner
>>> http://root.cern.ch/root/htmldoc/src/TRandom.cxx.html#TRandom:BreitWigner
>>>
>>> see also:
>>> http://root.cern.ch/root/roottalk/roottalk05/0954.html
>>>
>>> Rene Brun
>>>
>>> Snyder, Arthur E. wrote:
>>>
>>>> Does anybody know where to find a implentation of Wigner D-functions that I can use in ROOT?
>>>>
>>>> I found a talk/article by Rene where they are mentioned but the implentation seems to be in
>>>> the future.
>>>>
>>>> -Art S.
>>>>
>>>>
>>>> -----Original Message-----
>>>> From: owner-roottalk_at_pcroot.cern.ch on behalf of Fine, Valeri
>>>> Sent: Sun 7/2/2006 6:56 PM
>>>> To: nasim; roottalk_at_pcroot.cern.ch
>>>> Subject: RE: [ROOT] geting the values stored in the histograms.
>>>>
>>>>
>>>> You need either
>>>> http://root.cern.ch/root/htmldoc/TH1.html#TH1:GetBinLowEdge
>>>> or
>>>> http://root.cern.ch/root/htmldoc/TH1.html#TH1:GetBinCenter
>>>>
>>>> Valeri
>>>>
>>>>
>>>>> There are two ways:
>>>>>
>>>>> 1. Each concrete histogram class is derived from TArray class
>>>>> http://root.cern.ch/root/htmldoc/TH1F.html
>>>>> This means one can treat each histogram object as a plain array of
>>>>> "float" number and use it as such. (For example method "[int]} returns
>>>>> what you want)
>>>>>
>>>>> TH1F *histogram = new TH1F( some parameters);
>>>>> TH1F &bins = *histogram;
>>>>> for (int i=0; i< 100; i++)
>>>>> printf(" Bin %i contains the value %f\n", i, bins[i]);
>>>>>
>>>>> 2. One can use the methods like TH1::GetBin<XXX>
>>>>> , where XXX could be "Center", "Content", "Error" etc
>>>>> see doc please
>>>>> http://root.cern.ch/root/htmldoc/TH1.html#TH1:GetBinContent
>>>>>
>>>>>
>>>>>
>>>>
>>>>
>>>>> Hope this helps.
>>>>>
>>>>> ----
>>>>> Best regards
>>>>> Valeri
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>> -----Original Message-----
>>>>>> From: owner-roottalk_at_pcroot.cern.ch
>>>>>>
>>>>>>
>>>>> [mailto:owner-roottalk_at_pcroot.cern.ch] On
>>>>>
>>>>>
>>>>>> Behalf Of nasim
>>>>>> Sent: Sunday, July 02, 2006 7:43 PM
>>>>>> To: roottalk_at_pcroot.cern.ch
>>>>>> Subject: [ROOT] geting the values stored in the histograms.
>>>>>>
>>>>>>
>>>>>> Hi,
>>>>>> If we have one dimensional histogram with 100 entries.
>>>>>> Is it possible to extract the value of its x axis for each entry?
>>>>>>
>>>>>> Any help would be appreciated,
>>>>>>
>>>>>> Nasim
>>>>>>
>>>>>>
>>>>>> --
>>>>>>
>>>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>
>>

#include "TMath.h"
Double_t WignerD(Double_t aj, Double_t am, Double_t an, Double_t beta) {

   // Calculates the beta-term
   //                         d j mn (beta)
   // in the matrix element of the finite rotation operator
   // (Wigner's D-function), according to formula 4.3.1(3) in
   // D.A. Varshalovich, A.N. Moskalev, and V.K. Khersonskii,
   // Quantum Theory of Angular Momentum, World Scientific,
   // Singapore 1988.
   // CERNLIB DDJMNB function translated from Fortran to C++ by Rene Brun

   const Double_t f = 8.72664625997164788e-3;    const Double_t fcl[51] = { 0 ,0,

                              6.93147180559945309e-1 ,1.79175946922805500e00,
                              3.17805383034794562e00 ,4.78749174278204599e00,
                              6.57925121201010100e00 ,8.52516136106541430e00,
                              1.06046029027452502e01 ,1.28018274800814696e01,
                              1.51044125730755153e01 ,1.75023078458738858e01,
                              1.99872144956618861e01 ,2.25521638531234229e01,
                              2.51912211827386815e01 ,2.78992713838408916e01,
                              3.06718601060806728e01 ,3.35050734501368889e01,
                              3.63954452080330536e01 ,3.93398841871994940e01,
                              4.23356164607534850e01 ,4.53801388984769080e01,
                              4.84711813518352239e01 ,5.16066755677643736e01,
                              5.47847293981123192e01 ,5.80036052229805199e01,
                              6.12617017610020020e01 ,6.45575386270063311e01,
                              6.78897431371815350e01 ,7.12570389671680090e01,
                              7.46582363488301644e01 ,7.80922235533153106e01,
                              8.15579594561150372e01 ,8.50544670175815174e01,
                              8.85808275421976788e01 ,9.21361756036870925e01,
                              9.57196945421432025e01 ,9.93306124547874269e01,
                              1.02968198614513813e02 ,1.06631760260643459e02,
                              1.10320639714757395e02 ,1.14034211781461703e02,
                              1.17771881399745072e02 ,1.21533081515438634e02,
                              1.25317271149356895e02 ,1.29123933639127215e02,
                              1.32952575035616310e02 ,1.36802722637326368e02,
                              1.40673923648234259e02 ,1.44565743946344886e02,
                              1.48477766951773032e02};

   Int_t jpm = TMath::Nint(aj+am);
   Int_t jpn = TMath::Nint(aj+an);
   Int_t jmm = TMath::Nint(aj-am);

   Int_t jmn = TMath::Nint(aj-an);
   Int_t mpn = TMath::Nint(am+an);

   Double_t r = 0;
   if (jpm < 0 || jpn < 0 || jmm < 0 || jmn < 0

      || aj < 0 || aj > 25 || beta < 0 || beta > 360) {
       printf("WignerD: Illegal argument(s) aj=%g, am=%g, an=%g, beta=%g\n",aj,am,an,beta);
    } else if (beta == 0) {
       if (jpm == jpn) r = 1;
    } else if (beta == 180) {
       if (jpm == jmn) {
          r = 1;
          if (TMath::Abs(jpm)%2 == 1) r = -1;
       }
    } else if (beta == 360) {
       if (jpm == jpn) {
          r = 1;
          if (TMath::Abs(mpn)%2 == 1) r = -1;
       }
    } else {
       Double_t b  = f*beta;
       Double_t s  = TMath::Log(TMath::Sin(b));
       Double_t c  = TMath::Log(TMath::Abs(TMath::Cos(b)));
       Double_t rt = 0.5*(fcl[jpm]+fcl[jmm]+fcl[jpn]+fcl[jmn]);
       Int_t k0    = TMath::Max(0,mpn);
       Int_t kq    = k0+jpm;
       if (beta > 180) kq += mpn;
       Double_t q  = 1;
       if (kq%2 == 1) q = -1;
       kq = k0+k0;
       Double_t cx = kq-mpn;
       Double_t sx = jpm+jpn-kq;
       for (Int_t k=k0;k<=TMath::Min(jpm,jpn);k++) {
          r  += q*TMath::Exp(rt-fcl[k]-fcl[jpm-k]-fcl[jpn-k]-fcl[k-mpn]+ cx*c+sx*s);
          cx += 2;
          sx -= 2;
          q   = -q;
       }

   }
   return r;
} Received on Tue Jul 04 2006 - 09:32:42 MEST

This archive was generated by hypermail 2.2.0 : Mon Jan 01 2007 - 16:31:59 MET