Matchbox/Utility/AmplitudeCache.h

0001 // -*- C++ -*-
0002 //
0003 // AmplitudeCache.h is a part of Herwig - A multi-purpose Monte Carlo event generator
0004 // Copyright (C) 2002-2019 The Herwig Collaboration
0005 //
0006 // Herwig is licenced under version 3 of the GPL, see COPYING for details.
0007 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
0008 //
0009 #ifndef HERWIG_AmplitudeCache_H
0010 #define HERWIG_AmplitudeCache_H
0011
0012 #include "Herwig/MatrixElement/Matchbox/Utility/SpinorHelicity.h"
0013 #include "ThePEG/Config/algorithm.h"
0014 #include <array>
0015
0016 namespace Herwig {
0017
0018 using namespace ThePEG;
0019 using std::array;
0020
0021 namespace SpinorHelicity {
0022
0023 /**
0024  * \ingroup Matchbox
0025  * \author Simon Platzer
0026  *
0027  * \brief Caching for amplitudes using spinor helicity techniques.
0028  *
0029  */
0030 template<typename AmplitudeKey>
0031 class AmplitudeCache {
0032
0033   typedef map<AmplitudeKey,pair<bool,Complex> > AmplitudeCacheMap;
0034   typedef map<AmplitudeKey,pair<bool,LorentzVector<Complex> > > CurrentCacheMap;
0035
0036   /**
0037    * Maximum N we can handle, SYM_N is storage size for a symmetric matrix of N * N elements
0038    */
0039   enum { MAX_N = 7, SYM_N = MAX_N*(MAX_N+1)/2 };
0040
0041   /**
0042    * The number of points
0043    */
0044   int theNPoints;
0045
0046   /**
0047    * The energy scale to obtain dimensionless
0048    * quantities.
0049    */
0050   mutable Energy theScale;
0051
0052   /**
0053    * Masses indexed by partons
0054    */
0055   mutable array<double,MAX_N> theMasses;
0056
0057   /**
0058    * Momenta indexed by partons
0059    */
0060   mutable array<LorentzMomentum,MAX_N> theMomenta;
0061
0062   /**
0063    * Crossing signs indexed by partons
0064    */
0065   mutable array<int,MAX_N> theCrossingSigns;
0066
0067   /**
0068    * Plus spinors indexed by partons
0069    */
0070   mutable array<PlusSpinor,MAX_N> thePlusSpinors;
0071
0072   /**
0073    * Plus conjugate spinors indexed by partons
0074    */
0075   mutable array<PlusConjugateSpinor,MAX_N> thePlusConjugateSpinors;
0076
0077   /**
0078    * Invariants indexed by partons
0079    */
0080   mutable array<double,SYM_N> theInvariants;
0081
0082   /**
0083    * Flag products to be recalculated
0084    */
0085   mutable array<bool,SYM_N> getInvariant;
0086
0087   /**
0088    * Spinor products indexed by partons
0089    */
0090   mutable array<Complex,SYM_N> thePlusProducts;
0091
0092   /**
0093    * Flag products to be recalculated
0094    */
0095   mutable array<bool,SYM_N> getPlusProduct;
0096
0097   /**
0098    * Spinor currents indexed by partons
0099    */
0100   mutable array<LorentzVector<Complex>,SYM_N> thePlusCurrents;
0101
0102   /**
0103    * Flag currents to be recalculated
0104    */
0105   mutable array<bool,SYM_N> getPlusCurrent;
0106
0107   /**
0108    * Cache intermediate amplitudes by index
0109    */
0110   mutable AmplitudeCacheMap theCachedAmplitudes;
0111
0112   /**
0113    * The last query for a cached amplitude
0114    */
0115   mutable typename AmplitudeCacheMap::iterator theLastAmplitude;
0116
0117   /**
0118    * Cache intermediate currents by index
0119    */
0120   mutable CurrentCacheMap theCachedCurrents;
0121
0122   /**
0123    * The last query for a cached current
0124    */
0125   mutable typename CurrentCacheMap::iterator theLastCurrent;
0126
0127   /**
0128    * Helper function to index symmetric arrays, assumes i <= j.
0129    * Usual indexing function (N*i + j) corrected by triangle number for i-th row.
0130    */
0131   inline size_t idx(size_t i, size_t j) const {
0132     return MAX_N * i - (i + 1) * i / 2 + j;
0133   }
0134
0135   /**
0136    * Helper to reset flags
0137    */
0138   struct boolResetter {
0139     void operator()(pair<const AmplitudeKey,pair<bool,Complex> >& flag) const {
0140       flag.second.first = true;
0141     }
0142     void operator()(pair<const AmplitudeKey,pair<bool,LorentzVector<Complex> > >& flag) const {
0143       flag.second.first = true;
0144     }
0145   };
0146
0147 public:
0148
0149   /**
0150    * Constructor
0151    */
0152   AmplitudeCache() : theNPoints(0) {}
0153
0154   /**
0155    * Prepare for n-point amplitude
0156    */
0157   void nPoints(int n);
0158
0159   /**
0160    * Return the number of points
0161    */
0162   int nPoints() const {
0163     return theNPoints;
0164   }
0165
0166   /**
0167    * Set the energy scale to obtain dimensionless
0168    * quantities and flag all quantities to be recalculated.
0169    */
0170   void amplitudeScale(Energy s) const;
0171
0172   /**
0173    * Set the momentum for the k'th parton
0174    * and its associated mass.
0175    */
0176   void momentum(int k, const LorentzMomentum& p,
0177         bool getSpinors = true,
0178         Energy mass = ZERO) const;
0179
0180   /**
0181    * Reset flags
0182    */
0183   void reset() const;
0184
0185 public:
0186
0187   /**
0188    * Return the momentum for the k'th parton
0189    */
0190   LorentzVector<double> momentum(int k) const { return theMomenta[k]/theScale; }
0191
0192   /**
0193    * Get the energy scale to obtain dimensionless
0194    * quantities and flag all quantities to be recalculated.
0195    */
0196   Energy amplitudeScale() const { return theScale; }
0197
0198   /**
0199    * Return the mass associated to the k'th parton
0200    */
0201   double mass(int k) const { return theMasses[k]; }
0202
0203   /**
0204    * Return the crossing sign for the
0205    * i'th parton
0206    */
0207   int crossingSign(int i) const { return theCrossingSigns[i]; }
0208
0209   /**
0210    * Return the crossing sign for the
0211    * i'th and j'th parton
0212    */
0213   double crossingSign(int i, int j) const { return theCrossingSigns[i]*theCrossingSigns[j]; }
0214
0215   /**
0216    * Return (ij)
0217    */
0218   double invariant(int i, int j) const {
0219     if ( i == j ) return 0.;
0220     if ( i > j  ) swap(i,j);
0221     if ( getInvariant[idx(i,j)] ) {
0222       getInvariant[idx(i,j)] = false;
0223       theInvariants[idx(i,j)] = 2.*(momentum(i)*momentum(j));
0224     }
0225     return theInvariants[idx(i,j)];
0226   }
0227
0228   /**
0229    * Return <ij>
0230    */
0231   Complex plusProduct(int i, int j) const {
0232     if ( i== j )
0233       return 0.;
0234     bool swapij = (i > j);
0235     if ( swapij )
0236       swap(i,j);
0237     if ( getPlusProduct[idx(i,j)] ) {
0238       getPlusProduct[idx(i,j)] = false;
0239       thePlusProducts[idx(i,j)] =
0240     Complex(PlusSpinorProduct(thePlusConjugateSpinors[i],
0241                   thePlusSpinors[j]).eval() / theScale);
0242     }
0243     return swapij ? -thePlusProducts[idx(i,j)] : thePlusProducts[idx(i,j)];
0244   }
0245
0246   /**
0247    * Return [ij]
0248    */
0249   Complex minusProduct(int i, int j) const {
0250     if ( i== j )
0251       return 0.;
0252     return -crossingSign(i,j)*conj(plusProduct(i,j));
0253   }
0254
0255   /**
0256    * Return <i|\gamma^\mu|j]
0257    */
0258   LorentzVector<Complex> plusCurrent(int i, int j) const {
0259     bool swapij = (i > j);
0260     if ( swapij )
0261       swap(i,j);
0262     if ( getPlusCurrent[idx(i,j)] ) {
0263       getPlusCurrent[idx(i,j)] = false;
0264       if ( i != j ) {
0265     thePlusCurrents[idx(i,j)] =
0266       PlusSpinorCurrent(thePlusConjugateSpinors[i],
0267                 MinusSpinor(theMomenta[j])).eval() / theScale;
0268       } else {
0269     thePlusCurrents[idx(i,j)] = 2.*momentum(i);
0270       }
0271     }
0272     return swapij ? crossingSign(i,j)*thePlusCurrents[idx(i,j)].conjugate() : thePlusCurrents[idx(i,j)];
0273   }
0274
0275   /**
0276    * Return [i|\gamma^\mu|j>
0277    */
0278   LorentzVector<Complex> minusCurrent(int i, int j) const {
0279     return plusCurrent(j,i);
0280   }
0281
0282 public:
0283
0284   /**
0285    * Return true, if the given amplitude
0286    * needs to be recalculated.
0287    */
0288   bool getAmplitude(const AmplitudeKey& key) const {
0289     static Complex czero;
0290     if ( ( theLastAmplitude = theCachedAmplitudes.find(key) )
0291      == theCachedAmplitudes.end() ) {
0292       theLastAmplitude = theCachedAmplitudes.insert(make_pair(key,make_pair(true,czero))).first;
0293     }
0294     return theLastAmplitude->second.first;
0295   }
0296
0297   /**
0298    * Cache an amplitude
0299    */
0300   void cacheAmplitude(Complex amp) const {
0301     theLastAmplitude->second = make_pair(false,amp);
0302   }
0303
0304   /**
0305    * Return a cached amplitude
0306    */
0307   const Complex& cachedAmplitude() const {
0308     return theLastAmplitude->second.second;
0309   }
0310
0311   /**
0312    * Return true, if the given current
0313    * needs to be recalculated.
0314    */
0315   bool getCurrent(const AmplitudeKey& key) const {
0316     static LorentzVector<Complex> czero;
0317     if ( ( theLastCurrent = theCachedCurrents.find(key) )
0318      == theCachedCurrents.end() ) {
0319       theLastCurrent = theCachedCurrents.insert(make_pair(key,make_pair(true,czero))).first;
0320     }
0321     return theLastCurrent->second.first;
0322   }
0323
0324   /**
0325    * Cache an current
0326    */
0327   void cacheCurrent(const LorentzVector<Complex>& curr) const {
0328     theLastCurrent->second = make_pair(false,curr);
0329   }
0330
0331   /**
0332    * Return a cached current
0333    */
0334   const LorentzVector<Complex>& cachedCurrent() const {
0335     return theLastCurrent->second.second;
0336   }
0337
0338 };
0339
0340 }
0341
0342 }
0343
0344 #include "AmplitudeCache.tcc"
0345
0346 #endif // HERWIG_AmplitudeCache_H