KinTools.h

Go to the documentation of this file.

// --------------------------------------------------------------------------
// KinTools.cpp - Helper tools for DalitzGUI
// --------------------------------------------------------------------------
//
// Original author: Klaus Goetzen - GSI Darmstadt
// Last modified  : 2014/08/15
// --------------------------------------------------------------------------
#pragma once

#include <iostream>
#include "TComplex.h"
#include "TColor.h"
#include "TArrayD.h"
#include "TStyle.h"
#include "CRes.h"
#include "DalitzGlobals.h"

using namespace DalitzGuiGlobals;

//**************************************
//**************************************
//**************************************

namespace KinTools {

inline Double_t breakup(Double_t M, Double_t m1, Double_t m2)
{
  Double_t sum12 = m1 + m2;
  Double_t dif12 = m1 - m2;

  Double_t q = 0.;

  if (M > sum12)
    q = sqrt((M * M - sum12 * sum12) * (M * M - dif12 * dif12)) / (2. * M);

  return q;
}

inline Double_t BarrierP(Double_t q, Double_t q0, Int_t L, Double_t d = 5.07)
{
  Double_t result = 1.0;

  Double_t z0 = q0 * q0 * d * d;
  Double_t z = q * q * d * d;

  switch (L) {
  case 1: result = sqrt((1. + z0) / (1. + z)); break;

  case 2: result = sqrt(((z0 - 3.) * (z0 - 3.) + 9. * z0) / ((z - 3.) * (z - 3.) + 9. * z)); break;

  default: break;
  }

  return result;
}

inline Double_t Barrier(Double_t q, Int_t L, Double_t d = 5.07)
{
  Double_t result = 1.0;

  Double_t z = q * q * d * d;

  switch (L) {
  case 1: result = sqrt(2. * z / (1. + z)); break;

  case 2: result = sqrt((13. * z * z) / ((z - 3.) * (z - 3.) + 9. * z)); break;

  default: break;
  }

  return result;
}

inline Double_t BWGamma(Double_t m, Double_t m0, Double_t Gamma0, Double_t m1, Double_t m2, Int_t L, Double_t d = 5.07)
{
  Double_t q0 = breakup(m0, m1, m2);
  Double_t q = breakup(m, m1, m2);

  if (0. == q)
    return 0.;

  Double_t phsp = q / q0;

  switch (L) {
  case 1: phsp = phsp * phsp * phsp; break;
  case 2: phsp = phsp * phsp * phsp * phsp * phsp; break;
  default: break;
  }

  Double_t bar = BarrierP(q, q0, L, d);

  Double_t G = Gamma0 * phsp * (m0 / m) * bar * bar;

  return G;
}

inline TComplex BreitWigner(Double_t m, Double_t m0, Double_t Gamma0, Double_t m1, Double_t m2, Int_t L, Double_t d = 5.07)
{
  TComplex I(0., 1.);
  TComplex G(BWGamma(m, m0, Gamma0, m1, m2, L, d));

  TComplex BW = G / (m0 * m0 - m * m - I * m0 * G);

  return BW;
}

inline Double_t lambda(Double_t x, Double_t y, Double_t z)
{
  return x * x + y * y + z * z - 2. * x * y - 2. * y * z - 2. * z * x;
}

inline Double_t simin(Int_t i, Double_t m0, Double_t m1, Double_t m2, Double_t m3)
{
  Double_t result = 0.;

  switch (i) {
  case 1: result = (m2 + m3) * (m2 + m3); break;
  case 2: result = (m1 + m3) * (m1 + m3); break;
  case 3: result = (m1 + m2) * (m1 + m2); break;
  }
  return result;
}

inline Double_t simax(Int_t i, Double_t m0, Double_t m1, Double_t m2, Double_t m3)
{
  Double_t result = 0.;

  switch (i) {
  case 1: result = (m0 - m1) * (m0 - m1); break;
  case 2: result = (m0 - m2) * (m0 - m2); break;
  case 3: result = (m0 - m3) * (m0 - m3); break;
  }
  return result;
}

inline Double_t s2min(Double_t s1, Double_t m0, Double_t m1, Double_t m2, Double_t m3)
{
  Double_t s = m0 * m0;
  Double_t lamterm = sqrt(lambda(s1, s, m1 * m1)) * sqrt(lambda(s1, m2 * m2, m3 * m3));

  Double_t result = m1 * m1 + m3 * m3 + ((s - s1 - m1 * m1) * (s1 - m2 * m2 + m3 * m3) - lamterm) / (2. * s1);

  return result;
}

inline Double_t s2max(Double_t s1, Double_t m0, Double_t m1, Double_t m2, Double_t m3)
{
  Double_t s = m0 * m0;
  Double_t lamterm = sqrt(lambda(s1, s, m1 * m1)) * sqrt(lambda(s1, m2 * m2, m3 * m3));

  Double_t result = m1 * m1 + m3 * m3 + ((s - s1 - m1 * m1) * (s1 - m2 * m2 + m3 * m3) + lamterm) / (2. * s1);

  return result;
}

inline Double_t s3min(Double_t s1, Double_t m0, Double_t m1, Double_t m2, Double_t m3)
{
  Double_t s = m0 * m0;
  Double_t lamterm = sqrt(lambda(s1, s, m1 * m1)) * sqrt(lambda(s1, m3 * m3, m1 * m1));

  Double_t result = m1 * m1 + m2 * m2 + ((s - s1 - m1 * m1) * (s1 - m1 * m1 + m2 * m2) - lamterm) / (2. * s1);

  return result;
}

inline Double_t s3max(Double_t s1, Double_t m0, Double_t m1, Double_t m2, Double_t m3)
{
  Double_t s = m0 * m0;
  Double_t lamterm = sqrt(lambda(s1, s, m1 * m1)) * sqrt(lambda(s1, m3 * m3, m1 * m1));

  Double_t result = m1 * m1 + m2 * m2 + ((s - s1 - m1 * m1) * (s1 - m1 * m1 + m3 * m3) + lamterm) / (2. * s1);

  return result;
}

inline Double_t s1min(Double_t s2, Double_t m0, Double_t m1, Double_t m2, Double_t m3)
{
  Double_t s = m0 * m0;
  Double_t lamterm = sqrt(lambda(s2, s, m2 * m2)) * sqrt(lambda(s2, m3 * m3, m1 * m1));

  Double_t result = m2 * m2 + m3 * m3 + ((s - s2 - m2 * m2) * (s2 - m1 * m1 + m3 * m3) - lamterm) / (2. * s2);

  return result;
}

inline Double_t s1max(Double_t s2, Double_t m0, Double_t m1, Double_t m2, Double_t m3)
{
  Double_t s = m0 * m0;
  Double_t lamterm = sqrt(lambda(s2, s, m2 * m2)) * sqrt(lambda(s2, m1 * m1, m3 * m3));

  Double_t result = m2 * m2 + m3 * m3 + ((s - s2 - m2 * m2) * (s2 - m1 * m1 + m3 * m3) + lamterm) / (2. * s2);

  return result;
}

inline void palette2(Int_t mode = 0)
{
  const int ncol = 99;
  static Int_t col0[99], col1[99], col2[99], col3[99];
  static Bool_t initialized = kFALSE;

  if (!initialized) {
    {
      Double_t blu[3] = {1, 0.2, 0};
      Double_t red[3] = {1, 0.2, 0};
      Double_t grn[3] = {1, 0.2, 0};
      Double_t stp[3] = {0, 0.5, 1};

      Int_t FI = TColor::CreateGradientColorTable(3, stp, red, grn, blu, ncol);
      for (int i = 0; i < ncol; i++)
        col0[i] = FI + i;
    }
    {
      Double_t blu2[6] = {1, 1, 1, 0, 0, 0};
      Double_t red2[6] = {1, 0, 0, 0, 1, 1};
      Double_t grn2[6] = {1, 0, 1, 1, 1, 0};
      Double_t stp2[6] = {0.0, 0.12, 0.34, 0.56, 0.78, 1.0};

      Int_t FI = TColor::CreateGradientColorTable(6, stp2, red2, grn2, blu2, ncol);
      for (int i = 0; i < ncol; i++)
        col1[i] = FI + i;
    }
    {
      Double_t blu3[5] = {0, 0, 1, 0, 0};
      Double_t red3[5] = {0, 1, 1, 0, 0};
      Double_t grn3[5] = {0, 0, 1, 1, 0};
      Double_t stp3[5] = {0., 0.25, 0.5, 0.75, 1.0};

      Int_t FI = TColor::CreateGradientColorTable(5, stp3, red3, grn3, blu3, ncol);
      for (int i = 0; i < ncol; i++)
        col2[i] = FI + i;
    }
    {
      Double_t blu4[6] = {1, 0, 0, 0, 0.06, 0};
      Double_t red4[6] = {1, 1, 1, 0.5, 0.55, 0};
      Double_t grn4[6] = {1, 1, 0.5, 0.25, 0.10, 0};
      Double_t stp4[6] = {0., 0.2, 0.4, 0.60, 0.75, 1.0};

      Int_t FI = TColor::CreateGradientColorTable(6, stp4, red4, grn4, blu4, ncol);
      for (int i = 0; i < ncol; i++)
        col3[i] = FI + i;
    }
    initialized = kTRUE;
  }

  switch (mode) {
  case 0: gStyle->SetPalette(ncol, col0); break;
  case 1: gStyle->SetPalette(ncol, col1); break;
  case 2: gStyle->SetPalette(ncol, col2); break;
  case 3: gStyle->SetPalette(56, 0); break;
  }
}

inline Double_t Z_Ralt(Double_t s, Double_t smin, Double_t smax, Int_t J, double m_R = 0., double m_ab = 0., double m_c = 0.)
{
  Double_t res = 1.0;

  Double_t x = (s - smin) / (smax - smin) * 2. - 1.0;

  double relcorr = 1.0, xi2 = 1.;

  if (m_R > 1e-6) {
    xi2 = pow((m_R * m_R + m_ab * m_ab - m_c * m_c) / (2 * m_ab * m_R), 2) - 1;
  }

  switch (J) {
  case 1:
    res = 3. * x * x;
    if (m_R > 1e-6)
      relcorr = sqrt(1. + xi2);
    break;
  case 2:
    res = 5. * (x * x - 1. / 3.) * (x * x - 1. / 3.) * 9. / 4.;
    if (m_R > 1e-6)
      relcorr = sqrt(xi2 + 1.5);
    break;
  }

  return relcorr * sqrt(res); // *** new: return now sqrt(angle) to cope for computation of the intensity = amp*amp
}

inline Double_t Z_R(Double_t s, Double_t sab, Double_t sac, Double_t sbc, int i, Double_t mR, Int_t J, double fma, double fmb, double fmc)
{
  double res = 1.;

  if (J == 0)
    return res;

  double mR2 = 2 * mR;

  double mA2 = fma * fma, mB2 = fmb * fmb, mC2 = fmc * fmc;

  switch (J) {
  case 1:
    res = sac - sbc + (s - mC2) * (mA2 - mB2) / mR2; //(mAC*mAC-mBC*mBC+((mD*mD-mC*mC)*(mB*mB-mA*mA)/(mdenom*mdenom)))
    break;
  case 2:
    res = pow(sbc - sac + ((s - mC2) * (mA2 - mB2) / mR), 2) - 1. / 3. * (sab - 2 * s - 2 * mC2 + pow((s - mC2) / mR, 2)) * (sab - 2 * mA2 - 2 * mB2 + pow((mA2 - mB2) / mR, 2));
    break;
  }

  return res;
}

inline TComplex getAmp(CRes *r, Double_t m, Double_t m1, Double_t m2, Double_t qR, Double_t qM)
{
  TComplex bw1 = BreitWigner(m, r->GetM0(), r->GetG0(), m1, m2, r->GetJ());
  TComplex c = r->GetCoeff();
  c *= 1. / sqrt(r->GetG0()); // weight taking into account width (give narrow resonances more weight)
  Double_t bar1 = Barrier(qR, r->GetJ());
  Double_t bar2 = BarrierP(qM, breakup(r->GetM0(), m1, m2), r->GetJ());
  return bar1 * bar2 * c * bw1;
}

inline TComplex getAmpEvt(CRes *r, double sab, double sac, double sbc, double ma, double mb, double mc)
{
  double mAB2 = sab;
  // double mBC2 = sbc;
  // double mAC2 = sac;
  double mA2 = ma * ma;
  double mB2 = mb * mb;
  double mC2 = mc * mc;
  double mD2 = sab + sac + sbc - mA2 - mB2 - mC2;

  double mR = r->GetM0();
  double mR2 = mR * mR;
  double gammaR = r->GetG0();
  // double mdenom = mR;

  double pAB = sqrt((((mAB2 - mA2 - mB2) * (mAB2 - mA2 - mB2) / 4.0) - mA2 * mB2) / (mAB2));
  double pR = sqrt((((mR2 - mA2 - mB2) * (mR2 - mA2 - mB2) / 4.0) - mA2 * mB2) / (mR2));
  double pD = (((mD2 - mR2 - mC2) * (mD2 - mR2 - mC2) / 4.0) - mR2 * mC2) / (mD2);
  if (pD > 0) {
    pD = sqrt(pD);
  } else {
    pD = 0;
  }
  double pDAB = sqrt((((mD2 - mAB2 - mC2) * (mD2 - mAB2 - mC2) / 4.0) - mAB2 * mC2) / (mD2));

  double fR = 1;
  double fD = 1;
  int power = 0;
  int _spin = r->GetJ();
  switch (_spin) {
  case 0:
    fR = 1.0;
    fD = 1.0;
    power = 1;
    break;
  case 1:
    fR = sqrt(1.0 + 1.5 * 1.5 * pR * pR) / sqrt(1.0 + 1.5 * 1.5 * pAB * pAB);
    fD = sqrt(1.0 + 5.0 * 5.0 * pD * pD) / sqrt(1.0 + 5.0 * 5.0 * pDAB * pDAB);
    power = 3;
    break;
  case 2:
    fR = sqrt((9 + 3 * pow((1.5 * pR), 2) + pow((1.5 * pR), 4)) / (9 + 3 * pow((1.5 * pAB), 2) + pow((1.5 * pAB), 4)));
    fD = sqrt((9 + 3 * pow((5.0 * pD), 2) + pow((5.0 * pD), 4)) / (9 + 3 * pow((5.0 * pDAB), 2) + pow((5.0 * pDAB), 4)));
    power = 5;
    break;
  }

  double gammaAB = gammaR * pow(pAB / pR, power) * (mR / sqrt(mAB2)) * fR * fR;

  TComplex ampl = sqrt(r->GetG0()) * r->GetCoeff() * fR * fD / (mR2 - mAB2 - TComplex(0.0, mR * gammaAB));

  return ampl;
}

} // namespace KinTools