diff --git a/include/Rivet/Makefile.am b/include/Rivet/Makefile.am --- a/include/Rivet/Makefile.am +++ b/include/Rivet/Makefile.am @@ -1,345 +1,343 @@ ## Internal headers - not to be installed nobase_dist_noinst_HEADERS = ## Public headers - to be installed nobase_pkginclude_HEADERS = ## Rivet interface nobase_pkginclude_HEADERS += \ Rivet.hh \ Run.hh \ Event.hh \ ParticleBase.hh \ Particle.fhh Particle.hh \ Jet.fhh Jet.hh \ Projection.fhh Projection.hh \ ProjectionApplier.hh \ ProjectionHandler.hh \ Analysis.hh \ AnalysisHandler.hh \ AnalysisInfo.hh \ AnalysisBuilder.hh \ AnalysisLoader.hh ## Build config stuff nobase_pkginclude_HEADERS += \ Config/RivetCommon.hh ## Projections nobase_pkginclude_HEADERS += \ Projections/AxesDefinition.hh \ Projections/Beam.hh \ Projections/BeamThrust.hh \ Projections/CentralEtHCM.hh \ Projections/ChargedFinalState.hh \ Projections/ChargedLeptons.hh \ Projections/ConstLossyFinalState.hh \ Projections/DirectFinalState.hh \ Projections/DISFinalState.hh \ Projections/DISKinematics.hh \ Projections/DISLepton.hh \ Projections/DressedLeptons.hh \ Projections/FastJets.hh \ Projections/FinalPartons.hh \ Projections/FinalState.hh \ - Projections/FoxWolframMoments.hh \ Projections/FParameter.hh \ Projections/HadronicFinalState.hh \ Projections/HeavyHadrons.hh \ Projections/Hemispheres.hh \ Projections/IdentifiedFinalState.hh \ Projections/IndirectFinalState.hh \ Projections/InitialQuarks.hh \ Projections/InvMassFinalState.hh \ Projections/JetAlg.hh \ Projections/JetShape.hh \ Projections/LeadingParticlesFinalState.hh \ Projections/LossyFinalState.hh \ Projections/MergedFinalState.hh \ Projections/MissingMomentum.hh \ Projections/NeutralFinalState.hh \ Projections/NonHadronicFinalState.hh \ Projections/NonPromptFinalState.hh \ Projections/ParisiTensor.hh \ Projections/ParticleFinder.hh \ Projections/PartonicTops.hh \ Projections/PrimaryHadrons.hh \ Projections/PromptFinalState.hh \ Projections/SmearedParticles.hh \ Projections/SmearedJets.hh \ Projections/SmearedMET.hh \ Projections/Sphericity.hh \ Projections/Spherocity.hh \ Projections/TauFinder.hh \ Projections/Thrust.hh \ Projections/TriggerCDFRun0Run1.hh \ Projections/TriggerCDFRun2.hh \ Projections/TriggerUA5.hh \ Projections/UnstableFinalState.hh \ Projections/VetoedFinalState.hh \ Projections/VisibleFinalState.hh \ Projections/WFinder.hh \ Projections/ZFinder.hh ## Meta-projection convenience headers nobase_pkginclude_HEADERS += \ Projections/FinalStates.hh \ Projections/Smearing.hh ## Analysis base class headers # TODO: Move to Rivet/AnalysisTools header dir? nobase_pkginclude_HEADERS += \ Analyses/MC_ParticleAnalysis.hh \ Analyses/MC_JetAnalysis.hh \ Analyses/MC_JetSplittings.hh ## Tools nobase_pkginclude_HEADERS += \ Tools/BeamConstraint.hh \ Tools/BinnedHistogram.hh \ Tools/CentralityBinner.hh \ Tools/Cmp.fhh \ Tools/Cmp.hh \ Tools/Cutflow.hh \ Tools/Cuts.fhh \ Tools/Cuts.hh \ Tools/Exceptions.hh \ Tools/JetUtils.hh \ Tools/Logging.hh \ Tools/Random.hh \ Tools/ParticleBaseUtils.hh \ Tools/ParticleIdUtils.hh \ Tools/ParticleUtils.hh \ Tools/ParticleName.hh \ Tools/PrettyPrint.hh \ Tools/RivetPaths.hh \ Tools/RivetSTL.hh \ Tools/RivetFastJet.hh \ Tools/RivetHepMC.hh \ Tools/RivetYODA.hh \ Tools/RivetMT2.hh \ Tools/SmearingFunctions.hh \ Tools/MomentumSmearingFunctions.hh \ Tools/ParticleSmearingFunctions.hh \ Tools/JetSmearingFunctions.hh \ Tools/TypeTraits.hh \ Tools/Utils.hh \ Tools/Nsubjettiness/AxesDefinition.hh \ Tools/Nsubjettiness/ExtraRecombiners.hh \ Tools/Nsubjettiness/MeasureDefinition.hh \ Tools/Nsubjettiness/Njettiness.hh \ Tools/Nsubjettiness/NjettinessPlugin.hh \ Tools/Nsubjettiness/Nsubjettiness.hh \ Tools/Nsubjettiness/TauComponents.hh \ Tools/Nsubjettiness/XConePlugin.hh nobase_dist_noinst_HEADERS += \ Tools/osdir.hh \ Math/eigen3/COPYING.GPL \ Math/eigen3/README.md ## Maths nobase_pkginclude_HEADERS += \ Math/Matrices.hh \ Math/Vector3.hh \ Math/VectorN.hh \ Math/MatrixN.hh \ - Math/MatrixDiag.hh \ Math/MathHeader.hh \ Math/Vectors.hh \ Math/LorentzTrans.hh \ Math/Matrix3.hh \ Math/MathUtils.hh \ Math/Vector4.hh \ Math/Math.hh \ Math/Units.hh \ Math/Constants.hh \ Math/eigen3/Cholesky \ Math/eigen3/Core \ Math/eigen3/Dense \ Math/eigen3/Eigenvalues \ Math/eigen3/Geometry \ Math/eigen3/Householder \ Math/eigen3/Jacobi \ Math/eigen3/LU \ Math/eigen3/QR \ Math/eigen3/src/Cholesky/LDLT.h \ Math/eigen3/src/Cholesky/LLT.h \ Math/eigen3/src/Core/arch/AltiVec/Complex.h \ Math/eigen3/src/Core/arch/AltiVec/MathFunctions.h \ Math/eigen3/src/Core/arch/AltiVec/PacketMath.h \ Math/eigen3/src/Core/arch/AVX/Complex.h \ Math/eigen3/src/Core/arch/AVX/MathFunctions.h \ Math/eigen3/src/Core/arch/AVX/PacketMath.h \ Math/eigen3/src/Core/arch/AVX/TypeCasting.h \ Math/eigen3/src/Core/arch/AVX512/MathFunctions.h \ Math/eigen3/src/Core/arch/AVX512/PacketMath.h \ Math/eigen3/src/Core/arch/CUDA/Complex.h \ Math/eigen3/src/Core/arch/CUDA/Half.h \ Math/eigen3/src/Core/arch/CUDA/MathFunctions.h \ Math/eigen3/src/Core/arch/CUDA/PacketMath.h \ Math/eigen3/src/Core/arch/CUDA/PacketMathHalf.h \ Math/eigen3/src/Core/arch/CUDA/TypeCasting.h \ Math/eigen3/src/Core/arch/Default/Settings.h \ Math/eigen3/src/Core/arch/NEON/Complex.h \ Math/eigen3/src/Core/arch/NEON/MathFunctions.h \ Math/eigen3/src/Core/arch/NEON/PacketMath.h \ Math/eigen3/src/Core/arch/SSE/Complex.h \ Math/eigen3/src/Core/arch/SSE/MathFunctions.h \ Math/eigen3/src/Core/arch/SSE/PacketMath.h \ Math/eigen3/src/Core/arch/SSE/TypeCasting.h \ Math/eigen3/src/Core/arch/ZVector/Complex.h \ Math/eigen3/src/Core/arch/ZVector/MathFunctions.h \ Math/eigen3/src/Core/arch/ZVector/PacketMath.h \ Math/eigen3/src/Core/Array.h \ Math/eigen3/src/Core/ArrayBase.h \ Math/eigen3/src/Core/ArrayWrapper.h \ Math/eigen3/src/Core/Assign.h \ Math/eigen3/src/Core/AssignEvaluator.h \ Math/eigen3/src/Core/BandMatrix.h \ Math/eigen3/src/Core/Block.h \ Math/eigen3/src/Core/BooleanRedux.h \ Math/eigen3/src/Core/CommaInitializer.h \ Math/eigen3/src/Core/ConditionEstimator.h \ Math/eigen3/src/Core/CoreEvaluators.h \ Math/eigen3/src/Core/CoreIterators.h \ Math/eigen3/src/Core/CwiseBinaryOp.h \ Math/eigen3/src/Core/CwiseNullaryOp.h \ Math/eigen3/src/Core/CwiseTernaryOp.h \ Math/eigen3/src/Core/CwiseUnaryOp.h \ Math/eigen3/src/Core/CwiseUnaryView.h \ Math/eigen3/src/Core/DenseBase.h \ Math/eigen3/src/Core/DenseCoeffsBase.h \ Math/eigen3/src/Core/DenseStorage.h \ Math/eigen3/src/Core/Diagonal.h \ Math/eigen3/src/Core/DiagonalMatrix.h \ Math/eigen3/src/Core/DiagonalProduct.h \ Math/eigen3/src/Core/Dot.h \ Math/eigen3/src/Core/EigenBase.h \ Math/eigen3/src/Core/functors/AssignmentFunctors.h \ Math/eigen3/src/Core/functors/BinaryFunctors.h \ Math/eigen3/src/Core/functors/NullaryFunctors.h \ Math/eigen3/src/Core/functors/StlFunctors.h \ Math/eigen3/src/Core/functors/TernaryFunctors.h \ Math/eigen3/src/Core/functors/UnaryFunctors.h \ Math/eigen3/src/Core/Fuzzy.h \ Math/eigen3/src/Core/GeneralProduct.h \ Math/eigen3/src/Core/GenericPacketMath.h \ Math/eigen3/src/Core/GlobalFunctions.h \ Math/eigen3/src/Core/Inverse.h \ Math/eigen3/src/Core/IO.h \ Math/eigen3/src/Core/Map.h \ Math/eigen3/src/Core/MapBase.h \ Math/eigen3/src/Core/MathFunctions.h \ Math/eigen3/src/Core/MathFunctionsImpl.h \ Math/eigen3/src/Core/Matrix.h \ Math/eigen3/src/Core/MatrixBase.h \ Math/eigen3/src/Core/NestByValue.h \ Math/eigen3/src/Core/NoAlias.h \ Math/eigen3/src/Core/NumTraits.h \ Math/eigen3/src/Core/PermutationMatrix.h \ Math/eigen3/src/Core/PlainObjectBase.h \ Math/eigen3/src/Core/Product.h \ Math/eigen3/src/Core/ProductEvaluators.h \ Math/eigen3/src/Core/products/GeneralBlockPanelKernel.h \ Math/eigen3/src/Core/products/GeneralMatrixMatrix.h \ Math/eigen3/src/Core/products/GeneralMatrixMatrixTriangular.h \ Math/eigen3/src/Core/products/GeneralMatrixVector.h \ Math/eigen3/src/Core/products/Parallelizer.h \ Math/eigen3/src/Core/products/SelfadjointMatrixMatrix.h \ Math/eigen3/src/Core/products/SelfadjointMatrixVector.h \ Math/eigen3/src/Core/products/SelfadjointProduct.h \ Math/eigen3/src/Core/products/SelfadjointRank2Update.h \ Math/eigen3/src/Core/products/TriangularMatrixMatrix.h \ Math/eigen3/src/Core/products/TriangularMatrixVector.h \ Math/eigen3/src/Core/products/TriangularSolverMatrix.h \ Math/eigen3/src/Core/products/TriangularSolverVector.h \ Math/eigen3/src/Core/Random.h \ Math/eigen3/src/Core/Redux.h \ Math/eigen3/src/Core/Ref.h \ Math/eigen3/src/Core/Replicate.h \ Math/eigen3/src/Core/ReturnByValue.h \ Math/eigen3/src/Core/Reverse.h \ Math/eigen3/src/Core/Select.h \ Math/eigen3/src/Core/SelfAdjointView.h \ Math/eigen3/src/Core/SelfCwiseBinaryOp.h \ Math/eigen3/src/Core/Solve.h \ Math/eigen3/src/Core/SolverBase.h \ Math/eigen3/src/Core/SolveTriangular.h \ Math/eigen3/src/Core/StableNorm.h \ Math/eigen3/src/Core/Stride.h \ Math/eigen3/src/Core/Swap.h \ Math/eigen3/src/Core/Transpose.h \ Math/eigen3/src/Core/Transpositions.h \ Math/eigen3/src/Core/TriangularMatrix.h \ Math/eigen3/src/Core/util/BlasUtil.h \ Math/eigen3/src/Core/util/Constants.h \ Math/eigen3/src/Core/util/DisableStupidWarnings.h \ Math/eigen3/src/Core/util/ForwardDeclarations.h \ Math/eigen3/src/Core/util/Macros.h \ Math/eigen3/src/Core/util/Memory.h \ Math/eigen3/src/Core/util/Meta.h \ Math/eigen3/src/Core/util/MKL_support.h \ Math/eigen3/src/Core/util/NonMPL2.h \ Math/eigen3/src/Core/util/ReenableStupidWarnings.h \ Math/eigen3/src/Core/util/StaticAssert.h \ Math/eigen3/src/Core/util/XprHelper.h \ Math/eigen3/src/Core/VectorBlock.h \ Math/eigen3/src/Core/VectorwiseOp.h \ Math/eigen3/src/Core/Visitor.h \ Math/eigen3/src/Eigenvalues/ComplexEigenSolver.h \ Math/eigen3/src/Eigenvalues/ComplexSchur.h \ Math/eigen3/src/Eigenvalues/EigenSolver.h \ Math/eigen3/src/Eigenvalues/GeneralizedEigenSolver.h \ Math/eigen3/src/Eigenvalues/GeneralizedSelfAdjointEigenSolver.h \ Math/eigen3/src/Eigenvalues/HessenbergDecomposition.h \ Math/eigen3/src/Eigenvalues/MatrixBaseEigenvalues.h \ Math/eigen3/src/Eigenvalues/RealQZ.h \ Math/eigen3/src/Eigenvalues/RealSchur.h \ Math/eigen3/src/Eigenvalues/SelfAdjointEigenSolver.h \ Math/eigen3/src/Eigenvalues/Tridiagonalization.h \ Math/eigen3/src/Geometry/AlignedBox.h \ Math/eigen3/src/Geometry/AngleAxis.h \ Math/eigen3/src/Geometry/arch/Geometry_SSE.h \ Math/eigen3/src/Geometry/EulerAngles.h \ Math/eigen3/src/Geometry/Homogeneous.h \ Math/eigen3/src/Geometry/Hyperplane.h \ Math/eigen3/src/Geometry/OrthoMethods.h \ Math/eigen3/src/Geometry/ParametrizedLine.h \ Math/eigen3/src/Geometry/Quaternion.h \ Math/eigen3/src/Geometry/Rotation2D.h \ Math/eigen3/src/Geometry/RotationBase.h \ Math/eigen3/src/Geometry/Scaling.h \ Math/eigen3/src/Geometry/Transform.h \ Math/eigen3/src/Geometry/Translation.h \ Math/eigen3/src/Geometry/Umeyama.h \ Math/eigen3/src/Householder/BlockHouseholder.h \ Math/eigen3/src/Householder/Householder.h \ Math/eigen3/src/Householder/HouseholderSequence.h \ Math/eigen3/src/Jacobi/Jacobi.h \ Math/eigen3/src/LU/arch/Inverse_SSE.h \ Math/eigen3/src/LU/Determinant.h \ Math/eigen3/src/LU/FullPivLU.h \ Math/eigen3/src/LU/InverseImpl.h \ Math/eigen3/src/LU/PartialPivLU.h \ Math/eigen3/src/misc/Image.h \ Math/eigen3/src/misc/Kernel.h \ Math/eigen3/src/misc/RealSvd2x2.h \ Math/eigen3/src/plugins/ArrayCwiseBinaryOps.h \ Math/eigen3/src/plugins/ArrayCwiseUnaryOps.h \ Math/eigen3/src/plugins/BlockMethods.h \ Math/eigen3/src/plugins/CommonCwiseBinaryOps.h \ Math/eigen3/src/plugins/CommonCwiseUnaryOps.h \ Math/eigen3/src/plugins/MatrixCwiseBinaryOps.h \ Math/eigen3/src/plugins/MatrixCwiseUnaryOps.h \ Math/eigen3/src/QR/ColPivHouseholderQR.h \ Math/eigen3/src/QR/CompleteOrthogonalDecomposition.h \ Math/eigen3/src/QR/FullPivHouseholderQR.h \ Math/eigen3/src/QR/HouseholderQR.h \ Math/eigen3/src/SVD/BDCSVD.h \ Math/eigen3/src/SVD/JacobiSVD.h \ Math/eigen3/src/SVD/SVDBase.h \ Math/eigen3/src/SVD/UpperBidiagonalization.h \ Math/eigen3/SVD diff --git a/include/Rivet/Math/MatrixDiag.hh b/include/Rivet/Math/MatrixDiag.hh deleted file mode 100644 --- a/include/Rivet/Math/MatrixDiag.hh +++ /dev/null @@ -1,149 +0,0 @@ -#ifndef RIVET_MATH_MATRIXDIAG -#define RIVET_MATH_MATRIXDIAG - -#include "Rivet/Math/MathHeader.hh" -#include "Rivet/Math/MatrixN.hh" - -#include "gsl/gsl_vector.h" -#include "gsl/gsl_matrix.h" -#include "gsl/gsl_eigen.h" - -namespace Rivet { - - -template -class EigenSystem; -template -EigenSystem diagonalize(const Matrix& m); - - -/// Handy object containing results of a diagonalization. -template -class EigenSystem { - template - friend EigenSystem diagonalize(const Matrix&); - -public: - - typedef pair > EigenPair; - typedef vector EigenPairs; - - Vector getDiagVector() const { - assert(_eigenPairs.size() == N); - Vector ret; - for (size_t i = 0; i < N; ++i) { - ret.set(i, _eigenPairs[i].first); - } - return ret; - } - - Matrix getDiagMatrix() const { - return Matrix::mkDiag(getDiagVector()); - } - - EigenPairs getEigenPairs() const { - return _eigenPairs; - } - - vector getEigenValues() const { - assert(_eigenPairs.size() == N); - vector ret; - for (size_t i = 0; i < N; ++i) { - ret.push_back(_eigenPairs[i].first); - } - return ret; - } - - vector > getEigenVectors() const { - assert(_eigenPairs.size() == N); - vector > ret; - for (size_t i = 0; i < N; ++i) { - ret.push_back(_eigenPairs[i].second); - } - return ret; - } - -private: - - EigenPairs _eigenPairs; - -}; - - -/// Comparison functor for "eigen-pairs". -template -struct EigenPairCmp { - bool operator()(const typename EigenSystem::EigenPair& a, - const typename EigenSystem::EigenPair& b) const - { - return a.first < b.first; - } -}; - - -/// Diagonalize an NxN matrix, returning a collection of pairs of -/// eigenvalues and eigenvectors, ordered decreasing in eigenvalue. -template -EigenSystem diagonalize(const Matrix& m) { - EigenSystem esys; - - // Make a GSL matrix. - gsl_matrix* A = gsl_matrix_alloc(N, N); - for (size_t i = 0; i < N; ++i) { - for (size_t j = 0; j < N; ++j) { - gsl_matrix_set(A, i, j, m.get(i, j)); - } - } - - // Use GSL diagonalization. - gsl_matrix* vecs = gsl_matrix_alloc(N, N); - gsl_vector* vals = gsl_vector_alloc(N); - gsl_eigen_symmv_workspace* workspace = gsl_eigen_symmv_alloc(N); - gsl_eigen_symmv(A, vals, vecs, workspace); - gsl_eigen_symmv_sort(vals, vecs, GSL_EIGEN_SORT_VAL_DESC); - - // Build the vector of "eigen-pairs". - typename EigenSystem::EigenPairs eigensolns; - for (size_t i = 0; i < N; ++i) { - typename EigenSystem::EigenPair ep; - ep.first = gsl_vector_get(vals, i); - Vector ev; - for (size_t j = 0; j < N; ++j) { - ev.set(j, gsl_matrix_get(vecs, j, i)); - } - ep.second = ev; - eigensolns.push_back(ep); - } - - // Free GSL memory. - gsl_eigen_symmv_free(workspace); - gsl_matrix_free(A); - gsl_matrix_free(vecs); - gsl_vector_free(vals); - - // Populate the returned object. - esys._eigenPairs = eigensolns; - return esys; -} - - -template -inline const string toString(const typename EigenSystem::EigenPair& e) { - ostringstream ss; - //for (typename EigenSystem::EigenPairs::const_iterator i = e.begin(); i != e.end(); ++i) { - ss << e->first << " -> " << e->second; - // if (i+1 != e.end()) ss << '\n'; - //} - return ss.str(); -} - -template -inline ostream& operator<<(std::ostream& out, const typename EigenSystem::EigenPair& e) { - out << toString(e); - return out; -} - - -} - -#endif diff --git a/include/Rivet/Projections/FoxWolframMoments.hh b/include/Rivet/Projections/FoxWolframMoments.hh deleted file mode 100644 --- a/include/Rivet/Projections/FoxWolframMoments.hh +++ /dev/null @@ -1,72 +0,0 @@ -// -*- C++ -*- -#ifndef RIVET_FoxWolframMoments_HH -#define RIVET_FoxWolframMoments_HH - -#include "Rivet/Config/RivetCommon.hh" -#include "Rivet/Projection.hh" -#include "Rivet/Projections/FinalState.hh" -#include "Rivet/Projections/VetoedFinalState.hh" -#include "Rivet/Projections/VisibleFinalState.hh" -#include "Rivet/Particle.hh" -#include "Rivet/Event.hh" - -#include - -#define MAXMOMENT 5 - -namespace Rivet { - - - /// @brief Calculate Fox-Wolfram moments - class FoxWolframMoments : public Projection { - public: - - /// Constructor. - FoxWolframMoments(const FinalState& fsp) { - setName("FoxWolframMoments"); - declare(fsp, "FS"); - - /// @todo Let the user supply any projection they like? - VisibleFinalState vfs(fsp); - declare(vfs, "VFS"); - - // Initialize moments vector - for (int i = 0; i < MAXMOMENT ; ++i) { - _fwmoments.push_back(0.0); - } - } - - /// Clone on the heap. - DEFAULT_RIVET_PROJ_CLONE(FoxWolframMoments); - - - /// The projected Fox-Wolfram Moment of order l - double getFoxWolframMoment(unsigned int l) const { - if (l < MAXMOMENT) { - return _fwmoments[l]; - } - /// @todo What?!? - return -666.0; - } - - - protected: - - /// Apply the projection to the event. - void project(const Event& e); - - /// Compare projections. - CmpState compare(const Projection& p) const; - - - private: - - vector _fwmoments; - - }; - - -} - - -#endif diff --git a/src/Projections/FoxWolframMoments.cc b/src/Projections/FoxWolframMoments.cc deleted file mode 100644 --- a/src/Projections/FoxWolframMoments.cc +++ /dev/null @@ -1,62 +0,0 @@ -// -*- C++ -*- -#include "Rivet/Projections/FoxWolframMoments.hh" - -namespace Rivet { - - CmpState FoxWolframMoments::compare(const Projection& p) const { - return mkNamedPCmp(p, "FS"); - } - - - void FoxWolframMoments::project(const Event& e) { - // Project into final state and get total visible momentum - const FinalState& fs = applyProjection(e, "VFS"); - - // remember: # pairs = N! / ( r! * (N-r)! ) - // N.B.: Autocorrelations are included! Treat them separately as diagonal elements. - // see: http://cepa.fnal.gov/psm/simulation/mcgen/lund/pythia_manual/pythia6.3/pythia6301/node215.html - - double sumEnergy = 0.0; - for (Particles::const_iterator pi = fs.particles().begin(); pi != fs.particles().end(); ++pi) { - sumEnergy += pi->momentum().E(); - const FourMomentum pi_4 = pi->momentum(); - for (Particles::const_iterator pj = pi+1; pj != fs.particles().end(); ++pj) { - const FourMomentum pj_4 = pj->momentum(); - - // Calculate x_ij = cos(theta_ij) - double x_ij = 1.0; - if ( pi != pj ) { - double denom = pi_4.vector3().mod() * pj_4.vector3().mod(); - x_ij = pi_4.vector3().dot( pj_4.vector3() ) / denom; - } - - //const double core = fabs( pi_4 * pj_4 ); // / sumet2 ; - const double core = pi_4.vector3().mod() * pi_4.vector3().mod(); - - for (int order = 0; order < MAXMOMENT; ++order) { - // enter a factor 2.0 because ij = ji. Use symmetry to speed up! - _fwmoments[order] += 2.0 * core * gsl_sf_legendre_Pl( order, x_ij ) ; - } - } // end loop over p_j - - // Now add autocorrelations - // Obviously cos(theta_ij) = 1.0 - // Note that P_l(1) == 1 for each l - for (int order = 0; order < MAXMOMENT; ++order) { - _fwmoments[order] += fabs( pi_4 * pi_4 ); - } - } // end loop over p_i - - MSG_DEBUG("sumEnergy = " << sumEnergy); - - for (int order = 0; order < MAXMOMENT; ++order) { - _fwmoments[order] /= (sumEnergy*sumEnergy); - } - - // Normalize to H0 - for (int order = 1; order < MAXMOMENT; ++order) { - _fwmoments[order] /= _fwmoments[0]; - } - } - -}