Page MenuHomeHEPForge
Files F7878290

No OneTemporary
Actions

View Options

diff --git a/MatrixElement/Matchbox/Utility/AmplitudeCache.h b/MatrixElement/Matchbox/Utility/AmplitudeCache.h
--- a/MatrixElement/Matchbox/Utility/AmplitudeCache.h
+++ b/MatrixElement/Matchbox/Utility/AmplitudeCache.h
@@ -1,346 +1,346 @@
// -*- C++ -*-
//
// AmplitudeCache.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2012 The Herwig Collaboration
//
// Herwig is licenced under version 2 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 {
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief Caching for amplitudes using spinor helicity techniques.
*
*/
-template<class AmplitudeKey>
+template<typename AmplitudeKey>
class AmplitudeCache {
typedef map<AmplitudeKey,pair<bool,Complex> > AmplitudeCacheMap;
typedef map<AmplitudeKey,pair<bool,LorentzVector<Complex> > > CurrentCacheMap;
/**
+ * Maximum N we can handle, SYM_N is storage size for a symmetric matrix of N * N elements
+ */
+ enum { MAX_N = 6, SYM_N = MAX_N*(MAX_N+1)/2 };
+
+ /**
* The number of points
*/
int theNPoints;
/**
* The energy scale to obtain dimensionless
* quantities.
*/
mutable Energy theScale;
/**
* Masses indexed by partons
*/
- mutable vector<double> theMasses;
+ mutable array<double,MAX_N> theMasses;
/**
* Momenta indexed by partons
*/
- mutable vector<LorentzMomentum> theMomenta;
+ mutable array<LorentzMomentum,MAX_N> theMomenta;
/**
* Crossing signs indexed by partons
*/
- mutable vector<int> theCrossingSigns;
+ mutable array<int,MAX_N> theCrossingSigns;
/**
* Plus spinors indexed by partons
*/
- mutable vector<PlusSpinor> thePlusSpinors;
+ mutable array<PlusSpinor,MAX_N> thePlusSpinors;
/**
* Plus conjugate spinors indexed by partons
*/
- mutable vector<PlusConjugateSpinor> thePlusConjugateSpinors;
+ mutable array<PlusConjugateSpinor,MAX_N> thePlusConjugateSpinors;
/**
* Invariants indexed by partons
*/
- mutable vector<vector<double> > theInvariants;
+ mutable array<double,SYM_N> theInvariants;
/**
* Flag products to be recalculated
*/
- mutable vector<vector<bool> > getInvariant;
+ mutable array<bool,SYM_N> getInvariant;
/**
* Spinor products indexed by partons
*/
- mutable vector<vector<Complex> > thePlusProducts;
+ mutable array<Complex,SYM_N> thePlusProducts;
/**
* Flag products to be recalculated
*/
- mutable vector<vector<bool> > getPlusProduct;
+ mutable array<bool,SYM_N> getPlusProduct;
/**
* Spinor currents indexed by partons
*/
- mutable vector<vector<LorentzVector<Complex> > > thePlusCurrents;
+ mutable array<LorentzVector<Complex>,SYM_N> thePlusCurrents;
/**
* Flag currents to be recalculated
*/
- mutable vector<vector<bool> > getPlusCurrent;
+ mutable array<bool,SYM_N> getPlusCurrent;
/**
* Cache intermediate amplitudes by index
*/
mutable AmplitudeCacheMap theCachedAmplitudes;
/**
* The last query for a cached amplitude
*/
mutable typename AmplitudeCacheMap::iterator theLastAmplitude;
/**
* Cache intermediate currents by index
*/
mutable CurrentCacheMap theCachedCurrents;
/**
* The last query for a cached current
*/
mutable typename CurrentCacheMap::iterator theLastCurrent;
/**
+ * Helper function to index symmetric arrays, assumes i <= j.
+ * Usual indexing function (N*i + j) corrected by triangle number for i-th row.
+ */
+ inline size_t idx(size_t i, size_t j) const {
+ return MAX_N * i - (i + 1) * i / 2 + j;
+ }
+
+ /**
* Helper to reset flags
*/
struct boolResetter {
- void operator()(vector<bool>::reference flag) const {
- flag = true;
- }
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;
}
};
- /**
- * Helper to reset flags
- */
- struct boolVectorResetter {
- void operator()(vector<bool>& flags) const {
- for_each(flags.begin(),flags.end(),boolResetter());
- }
- };
-
public:
/**
* Constructor
*/
- AmplitudeCache()
- : theNPoints(0) {}
+ AmplitudeCache() : theNPoints(0) {}
/**
* Prepare for n-point amplitude
*/
void nPoints(int n);
/**
* Return the number of points
*/
int nPoints() const {
return theNPoints;
}
/**
* Set the energy scale to obtain dimensionless
* quantities and flag all quantities to be recalculated.
*/
void amplitudeScale(Energy s) const;
/**
* Set the momentum for the k'th parton
* and its associated mass.
*/
void momentum(int k, const LorentzMomentum& p,
bool getSpinors = true,
Energy mass = ZERO) const;
/**
* Reset flags
*/
void reset() const;
public:
/**
* Return the momentum for the k'th parton
*/
LorentzVector<double> momentum(int k) const { return theMomenta[k]/theScale; }
/**
* Get the energy scale to obtain dimensionless
* quantities and flag all quantities to be recalculated.
*/
Energy amplitudeScale() const { return theScale; }
/**
* Return the mass associated to the k'th parton
*/
double mass(int k) const { return theMasses[k]; }
/**
* Return the crossing sign for the
* i'th parton
*/
int crossingSign(int i) const { return theCrossingSigns[i]; }
/**
* Return the crossing sign for the
* i'th and j'th parton
*/
double crossingSign(int i, int j) const { return theCrossingSigns[i]*theCrossingSigns[j]; }
/**
* Return (ij)
*/
double invariant(int i, int j) const {
- if ( i== j )
- return 0.;
- if ( i > j )
- swap(i,j);
- if ( getInvariant[i][j] ) {
- getInvariant[i][j] = false;
- theInvariants[i][j] = 2.*(momentum(i)*momentum(j));
+ 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[i][j];
+ return theInvariants[idx(i,j)];
}
/**
* Return <ij>
*/
Complex plusProduct(int i, int j) const {
if ( i== j )
return 0.;
bool swapij = (i > j);
if ( swapij )
swap(i,j);
- if ( getPlusProduct[i][j] ) {
- getPlusProduct[i][j] = false;
- thePlusProducts[i][j] =
+ if ( getPlusProduct[idx(i,j)] ) {
+ getPlusProduct[idx(i,j)] = false;
+ thePlusProducts[idx(i,j)] =
PlusSpinorProduct(thePlusConjugateSpinors[i],
thePlusSpinors[j]).eval() / theScale;
}
- return swapij ? -thePlusProducts[i][j] : thePlusProducts[i][j];
+ return swapij ? -thePlusProducts[idx(i,j)] : thePlusProducts[idx(i,j)];
}
/**
* Return [ij]
*/
Complex minusProduct(int i, int j) const {
if ( i== j )
return 0.;
return -crossingSign(i,j)*conj(plusProduct(i,j));
}
/**
* Return <i|\gamma^\mu|j]
*/
LorentzVector<Complex> plusCurrent(int i, int j) const {
bool swapij = (i > j);
if ( swapij )
swap(i,j);
- if ( getPlusCurrent[i][j] ) {
- getPlusCurrent[i][j] = false;
+ if ( getPlusCurrent[idx(i,j)] ) {
+ getPlusCurrent[idx(i,j)] = false;
if ( i != j ) {
- thePlusCurrents[i][j] =
+ thePlusCurrents[idx(i,j)] =
PlusSpinorCurrent(thePlusConjugateSpinors[i],
MinusSpinor(theMomenta[j])).eval() / theScale;
} else {
- thePlusCurrents[i][j] = 2.*momentum(i);
+ thePlusCurrents[idx(i,j)] = 2.*momentum(i);
}
}
- return swapij ? crossingSign(i,j)*thePlusCurrents[i][j].conjugate() : thePlusCurrents[i][j];
+ return swapij ? crossingSign(i,j)*thePlusCurrents[idx(i,j)].conjugate() : thePlusCurrents[idx(i,j)];
}
/**
* Return [i|\gamma^\mu|j>
*/
LorentzVector<Complex> minusCurrent(int i, int j) const {
return plusCurrent(j,i);
}
public:
/**
* Return true, if the given amplitude
* needs to be recalculated.
*/
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;
}
/**
* Cache an amplitude
*/
void cacheAmplitude(Complex amp) const {
theLastAmplitude->second = make_pair(false,amp);
}
/**
* Return a cached amplitude
*/
const Complex& cachedAmplitude() const {
return theLastAmplitude->second.second;
}
/**
* Return true, if the given current
* needs to be recalculated.
*/
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;
}
/**
* Cache an current
*/
void cacheCurrent(const LorentzVector<Complex>& curr) const {
theLastCurrent->second = make_pair(false,curr);
}
/**
* Return a cached current
*/
const LorentzVector<Complex>& cachedCurrent() const {
return theLastCurrent->second.second;
}
};
}
}
#include "AmplitudeCache.tcc"
#endif // HERWIG_AmplitudeCache_H
diff --git a/MatrixElement/Matchbox/Utility/AmplitudeCache.tcc b/MatrixElement/Matchbox/Utility/AmplitudeCache.tcc
--- a/MatrixElement/Matchbox/Utility/AmplitudeCache.tcc
+++ b/MatrixElement/Matchbox/Utility/AmplitudeCache.tcc
@@ -1,73 +1,60 @@
// -*- C++ -*-
//
// AmplitudeCache.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2012 The Herwig Collaboration
//
// Herwig is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
namespace Herwig {
namespace SpinorHelicity {
-template<class AmplitudeKey>
+template<typename AmplitudeKey>
void AmplitudeCache<AmplitudeKey>::nPoints(int n) {
+ assert( n <= MAX_N );
+
theNPoints = n;
- theMasses.clear();
- theMomenta.clear();
- theCrossingSigns.clear();
- thePlusSpinors.clear();
- thePlusConjugateSpinors.clear();
- theInvariants.clear();
- thePlusProducts.clear();
- thePlusCurrents.clear();
- getInvariant.clear();
- getPlusProduct.clear();
- getPlusCurrent.clear();
-
- theMasses.resize(n);
- theMomenta.resize(n);
- theCrossingSigns.resize(n);
- thePlusSpinors.resize(n);
- thePlusConjugateSpinors.resize(n);
- theInvariants.resize(n,vector<double>(n));
- thePlusProducts.resize(n,vector<Complex>(n));
- thePlusCurrents.resize(n,vector<LorentzVector<Complex> >(n));
- getInvariant.resize(n,vector<bool>(n));
- getPlusProduct.resize(n,vector<bool>(n));
- getPlusCurrent.resize(n,vector<bool>(n));
+ theMasses.fill({});
+ theMomenta.fill({});
+ theCrossingSigns.fill({});
+ thePlusSpinors.fill(PlusSpinor());
+ thePlusConjugateSpinors.fill(PlusConjugateSpinor());
+ theInvariants.fill({});
+ thePlusProducts.fill({});
+ thePlusCurrents.fill({});
reset();
}
-template<class AmplitudeKey>
+template<typename AmplitudeKey>
void AmplitudeCache<AmplitudeKey>::amplitudeScale(Energy s) const {
theScale = s;
reset();
}
-template<class AmplitudeKey>
+template<typename AmplitudeKey>
void AmplitudeCache<AmplitudeKey>::momentum(int k, const LorentzMomentum& p,
bool getSpinors,
Energy mass) const {
theMasses[k] = mass/theScale;
theMomenta[k] = p;
if ( getSpinors ) {
theCrossingSigns[k] = p.t() > ZERO ? 1 : -1;
thePlusSpinors[k] = PlusSpinor(p);
thePlusConjugateSpinors[k] = PlusConjugateSpinor(p);
}
}
-template<class AmplitudeKey>
+template<typename AmplitudeKey>
void AmplitudeCache<AmplitudeKey>::reset() const {
- for_each(getInvariant.begin(),getInvariant.end(),boolVectorResetter());
- for_each(getPlusProduct.begin(),getPlusProduct.end(),boolVectorResetter());
- for_each(getPlusCurrent.begin(),getPlusCurrent.end(),boolVectorResetter());
+ getInvariant.fill(true);
+ getPlusProduct.fill(true);
+ getPlusCurrent.fill(true);
for_each(theCachedAmplitudes.begin(),theCachedAmplitudes.end(),boolResetter());
for_each(theCachedCurrents.begin(),theCachedCurrents.end(),boolResetter());
}
}}

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 5:38 PM (1 d, 12 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3797781
Default Alt Text
(12 KB)

Event Timeline