AngDist.cc

Go to the documentation of this file.

 
/* C++ */
#include <cmath>
#include <complex>
#include <iostream>
#include <vector>
 
/* CERN ROOT */
#include <Math/SpecFunc.h> /* legendre polynomials, 3j and 6j symbols */
 
#include "AngDist.hh"
#include "WignerD.hh"
 
namespace GamR {
namespace AngDist {

double Fk(double k, double jf, double lambda, double lambdaPrime, double ji)
{
   return pow((-1), -1 + ji + jf) * pow((2 * lambda + 1) * (2 * lambdaPrime + 1) * (2 * ji + 1) * (2 * k + 1), 0.5) *
          ROOT::Math::wigner_3j((int)(2 * lambda), (int)(2 * lambdaPrime), (int)(2 * k),
                                2, -2, 0) *
          ROOT::Math::wigner_6j((int)(2 * ji), (int)(2 * ji), (int)(2 * k),
                                (int)(2 * lambdaPrime), (int)(2 * lambda), (int)(2 * jf));
}

double Fgen(double k, double k2, double k1, double jf, double lambda, double lambdaPrime, double ji)
{
   return pow(-1, lambdaPrime + k2 + k1 + 1) *
          sqrt((2 * ji + 1) * (2 * jf + 1) * (2 * lambda + 1) * (2 * lambdaPrime + 1) * (2 * k1 + 1) * (2 * k2 + 1) *
               (2 * k + 1)) *
          ROOT::Math::wigner_3j((int)(2 * lambda), (int)(2 * lambdaPrime), (int)(2 * k), 2, -2, 0) *
          ROOT::Math::wigner_9j((int)(2 * jf), (int)(2 * lambda), (int)(2 * ji), (int)(2 * jf), (int)(2 * lambdaPrime),
                                (int)(2 * ji), (int)(2 * k2), (int)(2 * k), (int)(2 * k1));
}

double Ak(double k, double ji, double lambda, double lambdaPrime, double jf, double delta)
{
   return 1. / (1. + pow(delta, 2)) *
          (Fk(k, jf, lambda, lambda, ji) + 2 * delta * Fk(k, jf, lambda, lambdaPrime, ji) +
           pow(delta, 2) * Fk(k, jf, lambdaPrime, lambdaPrime, ji));
}

double Agen(double k, double k1, double k2, double jf, double lambda, double lambdaPrime, double ji, double delta)
{
   return 1. / (1. + pow(delta, 2)) *
          (Fgen(k, k2, k1, jf, lambda, lambda, ji) + 2 * delta * Fgen(k, k2, k1, jf, lambda, lambdaPrime, ji) +
           pow(delta, 2) * Fgen(k, k2, k1, jf, lambdaPrime, lambdaPrime, ji));
}

double W(double theta, double ji, double lambda, double lambdaPrime, double jf, double delta, StatTensor *rho,
         SolidAttenuation *Qk, double kmax /*= 4*/)
{
   std::complex<double> wtheta = 0;
   for (int ik = 0; ik <= int(kmax); ++ik) {
      double k = (double)ik;
      wtheta = wtheta + Ak(k, ji, lambda, lambdaPrime, jf, delta) * rho->GetBk(k) * Qk->Get()[ik] *
                           ROOT::Math::legendre(ik, cos(theta));
   }
   if (abs(wtheta.imag()) >= 1e-15)
      std::cout << "Warning: imaginary parts of statistical tensor have not cancelled" << std::endl;
   return wtheta.real();
}

 
  std::function<std::complex<double>(double, double, double, double, double) >
  TransFrame(std::complex<double> (&Gk)(double, double, double, double, double), double alpha, double beta, double gamma)
  {
    std::function<std::complex<double>(double, double, double, double, double)> newGk;
        
    newGk = [alpha, beta, gamma, &Gk](double q, double qprime, double k, double kprime, double t) {
      //std::cout << "Frame transformation" << std::endl;
      std::complex<double> sum = 0;
      for (int iQ = 0; iQ < 2*k + 1; ++iQ) {
        for (int iQprime = 0; iQprime < 2*kprime+1; ++iQprime) {
          double Q = (double)iQ - k;
          double Qprime = (double)iQprime - kprime;
          sum += std::conj(GamR::AngDist::BigD(k, q, Q, alpha, beta, gamma))*
            std::conj(Gk(Q, Qprime, k, kprime, t))*
            GamR::AngDist::BigD(kprime, qprime, Qprime, alpha, beta, gamma);
 
          // if (Q==Qprime && k==kprime) {
          //   std::cout << k << "   " << q << "  " << Q << "  " << std::conj(GamR::AngDist::BigD(k, q, Q, alpha, beta, gamma)) << std::endl;
          //   std::cout << kprime << "  " << qprime << "  " << Qprime << "  " <<  GamR::AngDist::BigD(kprime, qprime, Qprime, alpha, beta, gamma) << std::endl;
          //}
        }
      }
      //printf("Gkk = (%15.14f, %15.14f)\n", sum.real(), -sum.imag());
      return std::conj(sum);
    };
 
    return newGk;
  }
  
 
std::vector<double> W(std::vector<double> theta, double ji, double lambda, double lambdaPrime, double jf, double delta,
                      StatTensor *rho, SolidAttenuation *Qk, double kmax /*= 4*/)
{
   /* equation 2.43 */
   std::vector<std::complex<double>> Aks;
   std::vector<double> W;
   for (int ik = 0; ik <= int(kmax); ++ik) {
      double k = (double)ik;
      Aks.push_back(Ak(k, ji, lambda, lambdaPrime, jf, delta) * rho->GetBk(k));
   }
   for (int it = 0; it < (int)theta.size(); ++it) {
      std::complex<double> wtheta = 0;
      for (int ik = 0; ik <= int(kmax); ++ik) {
         double k = (double)ik;
         wtheta = wtheta + Aks[ik] * Qk->Get()[ik] * ROOT::Math::legendre(k, cos(theta[it]));
      }
      if (wtheta.imag() != 0)
         std::cout << "Warning: imaginary parts of statistical tensor have not "
                      "cancelled"
                   << std::endl;
 
      W.push_back(wtheta.real());
   }
   return W;
}
 
} /* namespace AngDist */
} /* namespace GamR */