25 #include <analysis/ContinuumSuppression/KsfwMoments.h>
37 inline double legendre(
const double z,
const int i)
42 case 2:
return 1.5 * z * z - 0.5;
43 case 3:
return z * (2.5 * z * z - 1.5);
44 case 4:
return (4.375 * z * z * z * z - 3.75 * z * z + 0.375);
53 std::vector<std::pair<TVector3, int>> p3_cms_q_sigA,
54 std::vector<std::pair<TVector3, int>> p3_cms_q_sigB,
55 std::vector<std::pair<TVector3, int>> p3_cms_q_roe,
56 const TLorentzVector& p_cms_missA,
57 const TLorentzVector& p_cms_missB,
67 m_mm2[0] = p_cms_missA.E() > 0
69 : -p_cms_missA.E() * p_cms_missA.E() - p_cms_missA.Vect().Mag2();
70 m_mm2[1] = p_cms_missB.E() > 0
72 : -p_cms_missB.E() * p_cms_missB.E() - p_cms_missB.Vect().Mag2();
77 std::vector<std::pair<TVector3, int>>::iterator pqi, pqj;
80 for (
int i = 0; i < 3; i++) {
81 for (
int k = 0; k < 5; k++) {
87 for (pqi = p3_cms_q_sigA.begin(); pqi != p3_cms_q_sigA.end(); ++pqi) {
88 const double pi_mag((pqi->first).Mag());
89 for (pqj = p3_cms_q_roe.begin(); pqj != p3_cms_q_roe.end(); ++pqj) {
90 const double pj_mag((pqj->first).Mag());
91 const double ij_cos((pqi->first) * (pqj->first) / pi_mag / pj_mag);
92 const int c_or_n(0 == (pqj->second) ? 1 : 0);
93 for (
int k = 0; k < 5; k++) {
94 m_Hso[0][c_or_n][k] += (k % 2)
95 ? (pqi->second) * (pqj->second) * pj_mag *
legendre(ij_cos, k)
99 const double p_miss_mag(p_cms_missA.Rho());
100 const double i_miss_cos((pqi->first)*p_cms_missA.Vect() / pi_mag / p_miss_mag);
101 for (
int k = 0; k < 5; k++) {
102 m_Hso[0][2][k] += (k % 2) ? 0 : p_miss_mag *
legendre(i_miss_cos, k);
107 for (pqi = p3_cms_q_sigB.begin(); pqi != p3_cms_q_sigB.end(); ++pqi) {
108 const double pi_mag((pqi->first).Mag());
109 for (pqj = p3_cms_q_roe.begin(); pqj != p3_cms_q_roe.end(); ++pqj) {
110 const double pj_mag((pqj->first).Mag());
111 const double ij_cos((pqi->first) * (pqj->first) / pi_mag / pj_mag);
112 const int c_or_n(0 == (pqj->second) ? 1 : 0);
113 for (
int k = 0; k < 5; k++) {
114 m_Hso[1][c_or_n][k] += (k % 2)
115 ? (pqi->second) * (pqj->second) * pj_mag *
legendre(ij_cos, k)
119 const double p_miss_mag(p_cms_missB.Rho());
120 const double i_miss_cos((pqi->first)*p_cms_missB.Vect() / pi_mag / p_miss_mag);
121 for (
int k = 0; k < 5; k++) {
122 m_Hso[1][2][k] += (k % 2) ? 0 : p_miss_mag *
legendre(i_miss_cos, k);
127 std::vector<std::pair<TVector3, int>> p3_cms_q_roeA(p3_cms_q_roe), p3_cms_q_roeB(p3_cms_q_roe);
128 p3_cms_q_roeA.emplace_back(p_cms_missA.Vect(), 0);
129 p3_cms_q_roeB.emplace_back(p_cms_missB.Vect(), 0);
132 for (
int k = 0; k < 5; k++) {
135 for (pqi = p3_cms_q_roeA.begin(); pqi != p3_cms_q_roeA.end(); ++pqi) {
136 const double pi_mag((pqi->first).Mag());
137 for (pqj = p3_cms_q_roeA.begin(); pqj != pqi; ++pqj) {
138 const double pj_mag((pqj->first).Mag());
139 const double ij_cos((pqi->first) * (pqj->first) / pi_mag / pj_mag);
140 for (
int k = 0; k < 5; k++) {
141 m_Hoo[0][k] += (k % 2)
142 ? (pqi->second) * (pqj->second) * pi_mag * pj_mag *
legendre(ij_cos, k)
143 : pi_mag * pj_mag *
legendre(ij_cos, k);
147 for (pqi = p3_cms_q_roeB.begin(); pqi != p3_cms_q_roeB.end(); ++pqi) {
148 const double pi_mag((pqi->first).Mag());
149 for (pqj = p3_cms_q_roeB.begin(); pqj != pqi; ++pqj) {
150 const double pj_mag((pqj->first).Mag());
151 const double ij_cos((pqi->first) * (pqj->first) / pi_mag / pj_mag);
152 for (
int k = 0; k < 5; k++) {
153 m_Hoo[1][k] += (k % 2)
154 ? (pqi->second) * (pqj->second) * pi_mag * pj_mag *
legendre(ij_cos, k)
155 : pi_mag * pj_mag *
legendre(ij_cos, k);
161 for (
int k = 0; k < 5; k++) {
162 for (
int j = 0; j < ((k % 2) ? 1 : 3); j++) {
163 m_Hso[0][j][k] /= Hso0_max;
164 m_Hso[1][j][k] /= Hso0_max;
166 m_Hoo[0][k] /= (Hso0_max * Hso0_max);
167 m_Hoo[1][k] /= (Hso0_max * Hso0_max);