Second emitted photon in γ cascade not sampled correctly

[op3]

When sampling from 0⁺ → 2⁺ → 2⁺ → 0⁺, the produced angular correlation is not valid. The direction of the second photon is (exclusively?) determined by the direction of the first photon.

Code to reproduce the problem:

#include <array>
#include <iostream>
#include <vector>

using std::array;
using std::cout;
using std::vector;

#include "AngularCorrelation.hh"
#include "CascadeRejectionSampler.hh"

int main() {

  auto M1 = Transition(magnetic, 2, electric, 4, 0.);
  auto E2 = Transition(electric, 4, magnetic, 6, 0.);
  auto S0p = State(0, positive);
  auto S2p = State(4, positive);

  vector<AngularCorrelation> cascade{
      AngularCorrelation(S0p, {{E2, S2p}, {M1, S2p}}),
      AngularCorrelation(S2p, {{M1, S2p}, {E2, S0p}})
  };
  CascadeRejectionSampler cas_rej_sam(cascade, 0, {0., 0., 0.}, false);

  cout << "array([\n";
  for (size_t i = 0; i < 1000000; i++) {
    vector<array<double, 2>> theta_phi = cas_rej_sam();
    auto [theta1, phi1] = theta_phi[1];
    if (0.49 * M_PI < theta1 && theta1 < 0.51 * M_PI && 0.49 * M_PI < phi1 &&
        phi1 < 0.51 * M_PI) {
      cout << "  [" << theta_phi[0][0] << ", " << theta_phi[0][1] << "],\n";
    }
  }
  cout << "])\n";
}

It can be seen that the first photon is always sent to ϑ = 0.25 π, ϕ = 0.5 π if the second photon is sent to ϑ = 0.5 π, ϕ = 0.5 π, which is certainly not what the angular distribution looks like, which seems to be calculated correctly.

[op3]

I updated the code-sample to work with my latest changes in #5:

#include <array>
#include <iostream>
#include <memory>
#include <vector>

#include "alpaca/AngCorrRejectionSampler.hh"
#include "alpaca/AngularCorrelation.hh"
#include "alpaca/CascadeSampler.hh"
#include "alpaca/EulerAngleRotation.hh"
#include "alpaca/ReferenceFrameSampler.hh"

using namespace alpaca;

int main() {

  auto M1 = Transition::M1();
  auto E2 = Transition::E2();
  auto S0p = State::EvenPlus(0);
  auto S2p = State::EvenPlus(2);

  auto corr_0p_2p_2p = AngularCorrelation(S0p, {{E2, S2p}, {M1, S2p}});
  auto corr_2p_2p_0p = AngularCorrelation(S2p, {{M1, S2p}, {E2, S0p}});

  std::vector<std::shared_ptr<ReferenceFrameSampler>> ac_rej_sam{
      std::make_shared<AngCorrRejectionSampler>(corr_0p_2p_2p, 0),
      std::make_shared<AngCorrRejectionSampler>(corr_2p_2p_0p, 1)};
  CascadeSampler cas_sam(ac_rej_sam);

  std::cout << "array([\n";
  for (size_t i = 0; i < 1000000; i++) {
    std::vector<std::array<double, 3>> Phi_Theta_Psis = cas_sam();
    std::vector<std::array<double, 2>> theta_phi;
    theta_phi.reserve(Phi_Theta_Psis.size());
    std::transform(Phi_Theta_Psis.cbegin(), Phi_Theta_Psis.cend(),
                   theta_phi.begin(),
                   [](const std::array<double, 3> &Phi_Theta_Psi) {
                     return euler_angle_transform::to_spherical(Phi_Theta_Psi);
                   });
    auto [theta1, phi1] = theta_phi[1];
    if (0.49 * M_PI < theta1 && theta1 < 0.51 * M_PI && 0.49 * M_PI < phi1 &&
        phi1 < 0.51 * M_PI) {
      std::cout << "  [" << theta_phi[0][0] << ", " << theta_phi[0][1]
                << "],\n";
    }
  }
  std::cout << "])\n";
}

The problem is not as obvious anymore, but the resulting sample still does not seem to follow the correct distribution (if my own analytical calculations are correct, which might not be the case.

[op3]

I just had a look at a few triple-cascades. The first emitted photon has the correct angular distribution, but the second one does not.

Also, the domain of θ is [0, 2π] for the first photon and [-1/2π, 3/2π] for the second photon.