clebschgordan.h

Go to the documentation of this file.

/*
 *    Copyright (c) 2015-2018
 *      SMASH Team
 *
 *    GNU General Public License (GPLv3 or later)
 */

#ifndef SRC_INCLUDE_CLEBSCHGORDAN_H_
#define SRC_INCLUDE_CLEBSCHGORDAN_H_

#include <algorithm>

#include "particletype.h"

namespace smash {

double clebsch_gordan(const int j_a, const int j_b, const int j_c,
                      const int m_a, const int m_b, const int m_c);

inline double isospin_clebsch_gordan_sqr_2to1(const ParticleType &p_a,
                                              const ParticleType &p_b,
                                              const ParticleType &Res) {
  const double cg =
      clebsch_gordan(p_a.isospin(), p_b.isospin(), Res.isospin(),
                     p_a.isospin3(), p_b.isospin3(), Res.isospin3());
  return cg * cg;
}

double isospin_clebsch_gordan_sqr_3to1(const ParticleType &p_a,
                                       const ParticleType &p_b,
                                       const ParticleType &p_c,
                                       const ParticleType &Res);

double isospin_clebsch_gordan_sqr_2to2(const ParticleType &p_a,
                                       const ParticleType &p_b,
                                       const ParticleType &p_c,
                                       const ParticleType &p_d,
                                       const int I = -1);

class I_tot_range {
 private:
  int I_min_;
  int I_max_;

 public:
  I_tot_range(const ParticleType &p_a, const ParticleType &p_b) {
    // Compute total isospin range with given particles.
    const int I_z_abs = std::abs(p_a.isospin3() + p_b.isospin3());
    I_max_ = p_a.isospin() + p_b.isospin();
    I_min_ = std::max(std::abs(p_a.isospin() - p_b.isospin()), I_z_abs);
  }

  I_tot_range(const ParticleType &p_a, const ParticleType &p_b,
              const ParticleType &p_c, const ParticleType &p_d) {
    // Compute total isospin range with given initial and final particles.
    const int I_z = p_a.isospin3() + p_b.isospin3();
    if (I_z != p_c.isospin3() + p_d.isospin3()) {
      /* This reaction is forbidden by isospin conservation.
       * Set impossible values to make sure an empty range is returned. */
      I_min_ = 1;
      I_max_ = 0;
      return;
    }
    I_max_ =
        std::min(p_a.isospin() + p_b.isospin(), p_c.isospin() + p_d.isospin());
    I_min_ = std::max(std::abs(p_a.isospin() - p_b.isospin()),
                      std::abs(p_c.isospin() - p_d.isospin()));
    I_min_ = std::max(I_min_, std::abs(I_z));
  }

  class iterator : public std::iterator<std::forward_iterator_tag, int> {
   private:
    int c_;
    I_tot_range &parent_;

   public:
    iterator(int start, I_tot_range &parent) : c_(start), parent_(parent) {}
    int operator*() { return c_; }
    const iterator *operator++() {
      c_ -= 2;
      return this;
    }
    iterator operator++(int) {
      c_ -= 2;
      return iterator(c_ + 2, parent_);
    }
    bool operator==(const iterator &other) { return c_ == other.c_; }
    bool operator!=(const iterator &other) { return c_ != other.c_; }
  };

  iterator begin() { return iterator(I_max_, *this); }
  iterator end() {
    if (I_min_ > I_max_) {
      return begin();
    }
    return iterator(I_min_ - 2, *this);
  }
};

}  // namespace smash

#endif  // SRC_INCLUDE_CLEBSCHGORDAN_H_