11 #include <analysis/ContinuumSuppression/HarmonicMoments.h>
23 double pMag = p.Mag();
24 double cTheta = p.Dot(
m_axis) / pMag;
30 double cTheta2 = cTheta * cTheta;
31 double cTheta3 = cTheta2 * cTheta;
32 double cTheta4 = cTheta2 * cTheta2;
36 m_moment[2] += pMag * 0.5 * (3.*cTheta2 - 1);
37 m_moment[3] += pMag * 0.5 * (5.*cTheta3 - 3.*cTheta);
38 m_moment[4] += pMag * 0.125 * (35.*cTheta4 - 30.*cTheta2 + 3.);
49 double pMag = p.Mag();
50 double cTheta = p.Dot(
m_axis) / pMag;
56 double cTheta2 = cTheta * cTheta;
57 double cTheta3 = cTheta2 * cTheta;
58 double cTheta4 = cTheta2 * cTheta2;
59 double cTheta5 = cTheta4 * cTheta;
60 double cTheta6 = cTheta3 * cTheta3;
61 double cTheta7 = cTheta6 * cTheta;
62 double cTheta8 = cTheta4 * cTheta4;
66 m_moment[2] += pMag * 0.5 * (3.*cTheta2 - 1);
67 m_moment[3] += pMag * 0.5 * (5.*cTheta3 - 3.*cTheta);
68 m_moment[4] += pMag * 0.125 * (35.*cTheta4 - 30.*cTheta2 + 3.);
69 m_moment[5] += pMag * 0.125 * (63.*cTheta5 - 70 * cTheta3 + 15.*cTheta);
70 m_moment[6] += pMag * 0.0625 * (231.*cTheta6 - 315 * cTheta4 + 105 * cTheta2 - 5.);
71 m_moment[7] += pMag * 0.0625 * (429.*cTheta7 - 693.*cTheta5 + 315.*cTheta3 - 35.*cTheta);
72 m_moment[8] += pMag * 0.0078125 * (6435.*cTheta8 - 12012.*cTheta6 + 6930.*cTheta4 - 1260.*cTheta2 + 35.);