AmplitudeCache.h

// -*- C++ -*-
//
// AmplitudeCache.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_AmplitudeCache_H
#define HERWIG_AmplitudeCache_H
 
#include "Herwig/MatrixElement/Matchbox/Utility/SpinorHelicity.h"
#include "ThePEG/Config/algorithm.h"
#include <array>
 
namespace Herwig {
 
using namespace ThePEG;
using std::array;
 
namespace SpinorHelicity {
 
template<typename AmplitudeKey>
class AmplitudeCache {
 
  typedef map<AmplitudeKey,pair<bool,Complex> > AmplitudeCacheMap;
  typedef map<AmplitudeKey,pair<bool,LorentzVector<Complex> > > CurrentCacheMap;
 
  enum { MAX_N = 7, SYM_N = MAX_N*(MAX_N+1)/2 };
 
  int theNPoints;
 
  mutable Energy theScale;
 
  mutable array<double,MAX_N> theMasses;
 
  mutable array<LorentzMomentum,MAX_N> theMomenta;
 
  mutable array<int,MAX_N> theCrossingSigns;
 
  mutable array<PlusSpinor,MAX_N> thePlusSpinors;
 
  mutable array<PlusConjugateSpinor,MAX_N> thePlusConjugateSpinors;
 
  mutable array<double,SYM_N> theInvariants;
 
  mutable array<bool,SYM_N> getInvariant;
 
  mutable array<Complex,SYM_N> thePlusProducts;
 
  mutable array<bool,SYM_N> getPlusProduct;
 
  mutable array<LorentzVector<Complex>,SYM_N> thePlusCurrents;
 
  mutable array<bool,SYM_N> getPlusCurrent;
 
  mutable AmplitudeCacheMap theCachedAmplitudes;
 
  mutable typename AmplitudeCacheMap::iterator theLastAmplitude;
 
  mutable CurrentCacheMap theCachedCurrents;
 
  mutable typename CurrentCacheMap::iterator theLastCurrent;
 
  inline size_t idx(size_t i, size_t j) const {
    return MAX_N * i - (i + 1) * i / 2 + j;
  }
 
  struct boolResetter {
    void operator()(pair<const AmplitudeKey,pair<bool,Complex> >& flag) const {
      flag.second.first = true;
    }
    void operator()(pair<const AmplitudeKey,pair<bool,LorentzVector<Complex> > >& flag) const {
      flag.second.first = true;
    }
  };
 
public:
 
  AmplitudeCache() : theNPoints(0) {}
 
  void nPoints(int n);
 
  int nPoints() const {
    return theNPoints;
  }
 
  void amplitudeScale(Energy s) const;
 
  void momentum(int k, const LorentzMomentum& p,
                bool getSpinors = true,
                Energy mass = ZERO) const;
 
  void reset() const;
 
public:
 
  LorentzVector<double> momentum(int k) const { return theMomenta[k]/theScale; }
 
  Energy amplitudeScale() const { return theScale; }
 
  double mass(int k) const { return theMasses[k]; }
 
  int crossingSign(int i) const { return theCrossingSigns[i]; }
 
  double crossingSign(int i, int j) const { return theCrossingSigns[i]*theCrossingSigns[j]; }
 
  double invariant(int i, int j) const {
    if ( i == j ) return 0.;
    if ( i > j  ) swap(i,j);
    if ( getInvariant[idx(i,j)] ) {
      getInvariant[idx(i,j)] = false;
      theInvariants[idx(i,j)] = 2.*(momentum(i)*momentum(j));
    }
    return theInvariants[idx(i,j)];
  }
 
  Complex plusProduct(int i, int j) const {
    if ( i== j )
      return 0.;
    bool swapij = (i > j);
    if ( swapij )
      swap(i,j);
    if ( getPlusProduct[idx(i,j)] ) {
      getPlusProduct[idx(i,j)] = false;
      thePlusProducts[idx(i,j)] = 
        Complex(PlusSpinorProduct(thePlusConjugateSpinors[i],
                                  thePlusSpinors[j]).eval() / theScale);
    }
    return swapij ? -thePlusProducts[idx(i,j)] : thePlusProducts[idx(i,j)];
  }
 
  Complex minusProduct(int i, int j) const {
    if ( i== j )
      return 0.;
    return -crossingSign(i,j)*conj(plusProduct(i,j));
  }
 
  LorentzVector<Complex> plusCurrent(int i, int j) const {
    bool swapij = (i > j);
    if ( swapij )
      swap(i,j);
    if ( getPlusCurrent[idx(i,j)] ) {
      getPlusCurrent[idx(i,j)] = false;
      if ( i != j ) {
        thePlusCurrents[idx(i,j)] = 
          PlusSpinorCurrent(thePlusConjugateSpinors[i],
                            MinusSpinor(theMomenta[j])).eval() / theScale;
      } else {
        thePlusCurrents[idx(i,j)] = 2.*momentum(i);
      }
    }
    return swapij ? crossingSign(i,j)*thePlusCurrents[idx(i,j)].conjugate() : thePlusCurrents[idx(i,j)];
  }
 
  LorentzVector<Complex> minusCurrent(int i, int j) const {
    return plusCurrent(j,i);
  }
 
public:
 
  bool getAmplitude(const AmplitudeKey& key) const {
    static Complex czero;
    if ( ( theLastAmplitude = theCachedAmplitudes.find(key) )
         == theCachedAmplitudes.end() ) {
      theLastAmplitude = theCachedAmplitudes.insert(make_pair(key,make_pair(true,czero))).first;
    }
    return theLastAmplitude->second.first;
  }
 
  void cacheAmplitude(Complex amp) const {
    theLastAmplitude->second = make_pair(false,amp);
  }
 
  const Complex& cachedAmplitude() const {
    return theLastAmplitude->second.second;
  }
 
  bool getCurrent(const AmplitudeKey& key) const {
    static LorentzVector<Complex> czero;
    if ( ( theLastCurrent = theCachedCurrents.find(key) )
         == theCachedCurrents.end() ) {
      theLastCurrent = theCachedCurrents.insert(make_pair(key,make_pair(true,czero))).first;
    }
    return theLastCurrent->second.first;
  }
 
  void cacheCurrent(const LorentzVector<Complex>& curr) const {
    theLastCurrent->second = make_pair(false,curr);
  }
 
  const LorentzVector<Complex>& cachedCurrent() const {
    return theLastCurrent->second.second;
  }
 
};
 
}
 
}
 
#include "AmplitudeCache.tcc"
 
#endif // HERWIG_AmplitudeCache_H