herwig is hosted by Hepforge, IPPP Durham
Herwig  7.2.1
AmplitudeCache.h
1 // -*- C++ -*-
2 //
3 // AmplitudeCache.h is a part of Herwig - A multi-purpose Monte Carlo event generator
4 // Copyright (C) 2002-2019 The Herwig Collaboration
5 //
6 // Herwig is licenced under version 3 of the GPL, see COPYING for details.
7 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
8 //
9 #ifndef HERWIG_AmplitudeCache_H
10 #define HERWIG_AmplitudeCache_H
11 
12 #include "Herwig/MatrixElement/Matchbox/Utility/SpinorHelicity.h"
14 #include <array>
15 
16 namespace Herwig {
17 
18 using namespace ThePEG;
19 using std::array;
20 
21 namespace SpinorHelicity {
22 
30 template<typename AmplitudeKey>
32 
33  typedef map<AmplitudeKey,pair<bool,Complex> > AmplitudeCacheMap;
34  typedef map<AmplitudeKey,pair<bool,LorentzVector<Complex> > > CurrentCacheMap;
35 
39  enum { MAX_N = 7, SYM_N = MAX_N*(MAX_N+1)/2 };
40 
45 
50  mutable Energy theScale;
51 
55  mutable array<double,MAX_N> theMasses;
56 
60  mutable array<LorentzMomentum,MAX_N> theMomenta;
61 
65  mutable array<int,MAX_N> theCrossingSigns;
66 
70  mutable array<PlusSpinor,MAX_N> thePlusSpinors;
71 
75  mutable array<PlusConjugateSpinor,MAX_N> thePlusConjugateSpinors;
76 
80  mutable array<double,SYM_N> theInvariants;
81 
85  mutable array<bool,SYM_N> getInvariant;
86 
90  mutable array<Complex,SYM_N> thePlusProducts;
91 
95  mutable array<bool,SYM_N> getPlusProduct;
96 
100  mutable array<LorentzVector<Complex>,SYM_N> thePlusCurrents;
101 
105  mutable array<bool,SYM_N> getPlusCurrent;
106 
110  mutable AmplitudeCacheMap theCachedAmplitudes;
111 
115  mutable typename AmplitudeCacheMap::iterator theLastAmplitude;
116 
120  mutable CurrentCacheMap theCachedCurrents;
121 
125  mutable typename CurrentCacheMap::iterator theLastCurrent;
126 
131  inline size_t idx(size_t i, size_t j) const {
132  return MAX_N * i - (i + 1) * i / 2 + j;
133  }
134 
138  struct boolResetter {
139  void operator()(pair<const AmplitudeKey,pair<bool,Complex> >& flag) const {
140  flag.second.first = true;
141  }
142  void operator()(pair<const AmplitudeKey,pair<bool,LorentzVector<Complex> > >& flag) const {
143  flag.second.first = true;
144  }
145  };
146 
147 public:
148 
152  AmplitudeCache() : theNPoints(0) {}
153 
157  void nPoints(int n);
158 
162  int nPoints() const {
163  return theNPoints;
164  }
165 
170  void amplitudeScale(Energy s) const;
171 
176  void momentum(int k, const LorentzMomentum& p,
177  bool getSpinors = true,
178  Energy mass = ZERO) const;
179 
183  void reset() const;
184 
185 public:
186 
190  LorentzVector<double> momentum(int k) const { return theMomenta[k]/theScale; }
191 
196  Energy amplitudeScale() const { return theScale; }
197 
201  double mass(int k) const { return theMasses[k]; }
202 
207  int crossingSign(int i) const { return theCrossingSigns[i]; }
208 
213  double crossingSign(int i, int j) const { return theCrossingSigns[i]*theCrossingSigns[j]; }
214 
218  double invariant(int i, int j) const {
219  if ( i == j ) return 0.;
220  if ( i > j ) swap(i,j);
221  if ( getInvariant[idx(i,j)] ) {
222  getInvariant[idx(i,j)] = false;
223  theInvariants[idx(i,j)] = 2.*(momentum(i)*momentum(j));
224  }
225  return theInvariants[idx(i,j)];
226  }
227 
231  Complex plusProduct(int i, int j) const {
232  if ( i== j )
233  return 0.;
234  bool swapij = (i > j);
235  if ( swapij )
236  swap(i,j);
237  if ( getPlusProduct[idx(i,j)] ) {
238  getPlusProduct[idx(i,j)] = false;
239  thePlusProducts[idx(i,j)] =
240  Complex(PlusSpinorProduct(thePlusConjugateSpinors[i],
241  thePlusSpinors[j]).eval() / theScale);
242  }
243  return swapij ? -thePlusProducts[idx(i,j)] : thePlusProducts[idx(i,j)];
244  }
245 
249  Complex minusProduct(int i, int j) const {
250  if ( i== j )
251  return 0.;
252  return -crossingSign(i,j)*conj(plusProduct(i,j));
253  }
254 
258  LorentzVector<Complex> plusCurrent(int i, int j) const {
259  bool swapij = (i > j);
260  if ( swapij )
261  swap(i,j);
262  if ( getPlusCurrent[idx(i,j)] ) {
263  getPlusCurrent[idx(i,j)] = false;
264  if ( i != j ) {
265  thePlusCurrents[idx(i,j)] =
266  PlusSpinorCurrent(thePlusConjugateSpinors[i],
267  MinusSpinor(theMomenta[j])).eval() / theScale;
268  } else {
269  thePlusCurrents[idx(i,j)] = 2.*momentum(i);
270  }
271  }
272  return swapij ? crossingSign(i,j)*thePlusCurrents[idx(i,j)].conjugate() : thePlusCurrents[idx(i,j)];
273  }
274 
278  LorentzVector<Complex> minusCurrent(int i, int j) const {
279  return plusCurrent(j,i);
280  }
281 
282 public:
283 
288  bool getAmplitude(const AmplitudeKey& key) const {
289  static Complex czero;
290  if ( ( theLastAmplitude = theCachedAmplitudes.find(key) )
291  == theCachedAmplitudes.end() ) {
292  theLastAmplitude = theCachedAmplitudes.insert(make_pair(key,make_pair(true,czero))).first;
293  }
294  return theLastAmplitude->second.first;
295  }
296 
300  void cacheAmplitude(Complex amp) const {
301  theLastAmplitude->second = make_pair(false,amp);
302  }
303 
307  const Complex& cachedAmplitude() const {
308  return theLastAmplitude->second.second;
309  }
310 
315  bool getCurrent(const AmplitudeKey& key) const {
316  static LorentzVector<Complex> czero;
317  if ( ( theLastCurrent = theCachedCurrents.find(key) )
318  == theCachedCurrents.end() ) {
319  theLastCurrent = theCachedCurrents.insert(make_pair(key,make_pair(true,czero))).first;
320  }
321  return theLastCurrent->second.first;
322  }
323 
327  void cacheCurrent(const LorentzVector<Complex>& curr) const {
328  theLastCurrent->second = make_pair(false,curr);
329  }
330 
335  return theLastCurrent->second.second;
336  }
337 
338 };
339 
340 }
341 
342 }
343 
344 #include "AmplitudeCache.tcc"
345 
346 #endif // HERWIG_AmplitudeCache_H
int crossingSign(int i) const
Return the crossing sign for the i&#39;th parton.
bool getAmplitude(const AmplitudeKey &key) const
Return true, if the given amplitude needs to be recalculated.
array< PlusSpinor, MAX_N > thePlusSpinors
Plus spinors indexed by partons.
std::complex< double > Complex
CurrentCacheMap::iterator theLastCurrent
The last query for a cached current.
LorentzVector< double > momentum(int k) const
Return the momentum for the k&#39;th parton.
size_t idx(size_t i, size_t j) const
Helper function to index symmetric arrays, assumes i <= j.
array< Complex, SYM_N > thePlusProducts
Spinor products indexed by partons.
array< bool, SYM_N > getPlusProduct
Flag products to be recalculated.
double mass(int k) const
Return the mass associated to the k&#39;th parton.
ResultType eval() const
Return result.
int theNPoints
The number of points.
Base class for Weyl spinors.
const Complex & cachedAmplitude() const
Return a cached amplitude.
CurrentCacheMap theCachedCurrents
Cache intermediate currents by index.
double crossingSign(int i, int j) const
Return the crossing sign for the i&#39;th and j&#39;th parton.
array< double, MAX_N > theMasses
Masses indexed by partons.
array< PlusConjugateSpinor, MAX_N > thePlusConjugateSpinors
Plus conjugate spinors indexed by partons.
Energy amplitudeScale() const
Get the energy scale to obtain dimensionless quantities and flag all quantities to be recalculated...
array< bool, SYM_N > getPlusCurrent
Flag currents to be recalculated.
array< bool, SYM_N > getInvariant
Flag products to be recalculated.
Complex minusProduct(int i, int j) const
Return [ij].
const LorentzVector< Complex > & cachedCurrent() const
Return a cached current.
Complex plusProduct(int i, int j) const
Return <ij>
AmplitudeCacheMap theCachedAmplitudes
Cache intermediate amplitudes by index.
double invariant(int i, int j) const
Return (ij)
array< LorentzVector< Complex >, SYM_N > thePlusCurrents
Spinor currents indexed by partons.
void cacheAmplitude(Complex amp) const
Cache an amplitude.
array< LorentzMomentum, MAX_N > theMomenta
Momenta indexed by partons.
array< int, MAX_N > theCrossingSigns
Crossing signs indexed by partons.
-*- C++ -*-
array< double, SYM_N > theInvariants
Invariants indexed by partons.
LorentzVector< Complex > plusCurrent(int i, int j) const
Return <i|^|j].
Caching for amplitudes using spinor helicity techniques.
AmplitudeCacheMap::iterator theLastAmplitude
The last query for a cached amplitude.
int nPoints() const
Return the number of points.
void cacheCurrent(const LorentzVector< Complex > &curr) const
Cache an current.
bool getCurrent(const AmplitudeKey &key) const
Return true, if the given current needs to be recalculated.
Energy theScale
The energy scale to obtain dimensionless quantities.
constexpr ZeroUnit ZERO
LorentzVector< Complex > minusCurrent(int i, int j) const
Return [i|^|j>