eemmg-lib/src/kinem.h

Go to the documentation of this file.

/*
 * kinem.h - kinematics classes
 *
 * this file is part of PJFry library
 * Copyright 2011 Valery Yundin
 */
 
#ifndef QUL_KINEM_H
#define QUL_KINEM_H
 
#include "common.h"
#include <cstdlib>
// #include <cfloat>
 
template <int R>
class Kinem
{
  public:
    bool operator == (const Kinem& kinem) const;
 
    inline
    double mass(int i) const
    {
      return kdata[(i-1)*(i+2)/2];
    }
 
  protected:
    Kinem() {}
 
    void zero(); // initialize to zero
 
    typedef union { double d64; int64_t i64; } DI64;
//     static const double mdelta=1e-14;
//     static const int64_t idelta=1+mdelta/DBL_EPSILON; // needs cfloat
    /*
    Comparison tolerance for double
    idelta   fdelta(%)
        5  - 1.11e-15
       46  - 1.02e-14
      451  - 1.00e-13
     4504  - 1.00e-12
    45036  - 1.00e-11
    */
    static const uint64_t idelta=5; // unsigned !!!
 
    inline
    bool eq(const double& a, const double &b) const
    {
        const DI64 ia={a};
        const DI64 ib={b};
        const int64_t diff=ia.i64-ib.i64;
        return diff == 0LL || static_cast<uint64_t>(llabs(diff)) <= idelta;
    }
    inline
    bool neq(const double& a, const double &b) const
    {
        const DI64 ia={a};
        const DI64 ib={b};
        const int64_t diff=ia.i64-ib.i64;
        return diff != 0LL && static_cast<uint64_t>(llabs(diff)) > idelta;
    }
 
    static const int KLEN=R*(R+1)/2;
    enum Invar {im1=0,
                ip1, im2,
                ip2, ip3, im3,
                ip4, is12, is23, im4,
                ip5, is34, is45, is15, im5,
                ip6, /*   6-point inv    */ im6};
 
    double kdata[KLEN];
};
 
template <int R>
bool Kinem<R>::operator == (const Kinem<R>& kinem) const
{
  for (int i=0; i<KLEN; i++) {
    if (not eq(kdata[i], kinem.kdata[i])) return false;
  }
  return true;
}
 
template <int R>
void Kinem<R>::zero()
{
  for (int i=0; i<KLEN; i++) {
    kdata[i]=0;
  }
}
 
// 1-point kinematics
class Kinem1 : public Kinem<1>
{
  public:
    Kinem1() { zero(); }
    Kinem1(double xm1) {
      kdata[im1]=xm1;
    }
 
    inline double m1()  const { return kdata[im1];  }
};
 
// 2-point kinematics
class Kinem2 : public Kinem<2>
{
  public:
    Kinem2() { zero(); }
    Kinem2(double xp1,
           double xm1,  double xm2)
    {
      kdata[ip1]=xp1;
      kdata[im1]=xm1;
      kdata[im2]=xm2;
    }
 
    inline double p1()  const { return kdata[ip1];  }
    inline double m1()  const { return kdata[im1];  }
    inline double m2()  const { return kdata[im2];  }
};
 
// 3-point kinematics
class Kinem3 : public Kinem<3>
{
  public:
    Kinem3() { zero(); }
    Kinem3(double xp1,  double xp2,  double xp3,
           double xm1,  double xm2,  double xm3)
    {
      kdata[ip1]=xp1;
      kdata[ip2]=xp2;
      kdata[ip3]=xp3;
      kdata[im1]=xm1;
      kdata[im2]=xm2;
      kdata[im3]=xm3;
    }
 
    inline double p1()  const { return kdata[ip1];  }
    inline double p2()  const { return kdata[ip2];  }
    inline double p3()  const { return kdata[ip3];  }
    inline double m1()  const { return kdata[im1];  }
    inline double m2()  const { return kdata[im2];  }
    inline double m3()  const { return kdata[im3];  }
};
 
// 4-point kinematics
class Kinem4 : public Kinem<4>
{
  public:
    Kinem4() { zero(); }
    Kinem4(double xp1,  double xp2,  double xp3,  double xp4,
           double xs12, double xs23,
           double xm1,  double xm2,  double xm3,  double xm4)
    {
      kdata[ip1]=xp1;
      kdata[ip2]=xp2;
      kdata[ip3]=xp3;
      kdata[ip4]=xp4;
      kdata[is12]=xs12;
      kdata[is23]=xs23;
      kdata[im1]=xm1;
      kdata[im2]=xm2;
      kdata[im3]=xm3;
      kdata[im4]=xm4;
    }
 
    inline double p1()  const { return kdata[ip1];  }
    inline double p2()  const { return kdata[ip2];  }
    inline double p3()  const { return kdata[ip3];  }
    inline double p4()  const { return kdata[ip4];  }
    inline double m1()  const { return kdata[im1];  }
    inline double m2()  const { return kdata[im2];  }
    inline double m3()  const { return kdata[im3];  }
    inline double m4()  const { return kdata[im4];  }
    inline double s12() const { return kdata[is12]; }
    inline double s23() const { return kdata[is23]; }
};
 
// 5-point kinematics
class Kinem5 : public Kinem<5>
{
  public:
    Kinem5() { zero(); }
    Kinem5(double xp1,  double xp2,  double xp3,  double xp4,  double xp5,
           double xs12, double xs23, double xs34, double xs45, double xs15,
           double xm1,  double xm2,  double xm3,  double xm4,  double xm5)
    {
      kdata[ip1]=xp1;
      kdata[ip2]=xp2;
      kdata[ip3]=xp3;
      kdata[ip4]=xp4;
      kdata[ip5]=xp5;
      kdata[is12]=xs12;
      kdata[is23]=xs23;
      kdata[is34]=xs34;
      kdata[is45]=xs45;
      kdata[is15]=xs15;
      kdata[im1]=xm1;
      kdata[im2]=xm2;
      kdata[im3]=xm3;
      kdata[im4]=xm4;
      kdata[im5]=xm5;
    }
 
    inline double p1()  const { return kdata[ip1];  }
    inline double p2()  const { return kdata[ip2];  }
    inline double p3()  const { return kdata[ip3];  }
    inline double p4()  const { return kdata[ip4];  }
    inline double p5()  const { return kdata[ip5];  }
    inline double m1()  const { return kdata[im1];  }
    inline double m2()  const { return kdata[im2];  }
    inline double m3()  const { return kdata[im3];  }
    inline double m4()  const { return kdata[im4];  }
    inline double m5()  const { return kdata[im5];  }
    inline double s12() const { return kdata[is12]; }
    inline double s23() const { return kdata[is23]; }
    inline double s34() const { return kdata[is34]; }
    inline double s45() const { return kdata[is45]; }
    inline double s15() const { return kdata[is15]; }
};
 
#endif /* QUL_KINEM_H */