Belle II Software  release-05-01-25
KsfwMoments.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2010 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Mikihiko Nakao (KEK), Pablo Goldenzweig (KIT) *
7  * Original module writen by M. Nakao for Belle *
8  * Ported to Belle II by P. Goldenzweig *
9  * *
10  * This software is provided "as is" without any warranty. *
11  **************************************************************************/
12 
14 // //
15 // KsfwMoments.cc //
16 // //
17 // moment-calculation of the k_sfw //
18 // improved Super-Fox-Wolfram moments //
19 // to be used with rooksfw.{cc,h} //
20 // //
21 // M. Nakao //
22 // //
24 
25 #include <analysis/ContinuumSuppression/KsfwMoments.h>
26 
27 
28 namespace Belle2 {
37  inline double legendre(const double z, const int i)
38  {
39  switch (i) {
40  case 0: return 1.;
41  case 1: return z;
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);
45  default: return 0;
46  }
47  }
48 // ----------------------------------------------------------------------
49 // Constructor
50 // (copied from k_sfw.cc with minimum modification)
51 // ----------------------------------------------------------------------
52  KsfwMoments::KsfwMoments(double Hso0_max,
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,
58  double et[2]
59  )
60  {
61  // Private member needs to be initialized. Here it is initalized to -1 (illegal value).
62  m_uf = -1;
63 
64  m_et[0] = et[0];
65  m_et[1] = et[1];
66 
67  m_mm2[0] = p_cms_missA.E() > 0
68  ? p_cms_missA.Mag2()
69  : -p_cms_missA.E() * p_cms_missA.E() - p_cms_missA.Vect().Mag2();
70  m_mm2[1] = p_cms_missB.E() > 0
71  ? p_cms_missB.Mag2()
72  : -p_cms_missB.E() * p_cms_missB.E() - p_cms_missB.Vect().Mag2();
73 
74  //========================
75  // Calculate discriminants
76  //========================
77  std::vector<std::pair<TVector3, int>>::iterator pqi, pqj;
78 
79  // Calculate Hso components
80  for (int i = 0; i < 3; i++) {
81  for (int k = 0; k < 5; k++) {
82  m_Hso[0][i][k] = m_Hso[1][i][k] = 0;
83  }
84  }
85 
86  // Signal A (use_finalstate_for_sig == 0)
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); // 0: charged 1: neutral
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)
96  : pj_mag * legendre(ij_cos, k);
97  }
98  }
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);
103  }
104  }
105 
106  // Signal B (use_finalstate_for_sig == 1)
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); // 0: charged 1: neutral
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)
116  : pj_mag * legendre(ij_cos, k);
117  }
118  }
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);
123  }
124  }
125 
126  // Add missing to the lists
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);
130 
131  // Calculate Hoo components
132  for (int k = 0; k < 5; k++) {
133  m_Hoo[0][k] = m_Hoo[1][k] = 0;
134  }
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);
144  }
145  }
146  }
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);
156  }
157  }
158  }
159 
160  // Nomalize so that it does not dependent on delta_e
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;
165  }
166  m_Hoo[0][k] /= (Hso0_max * Hso0_max);
167  m_Hoo[1][k] /= (Hso0_max * Hso0_max);
168  }
169 
170  }
171 
173 } // Belle2 namespace
Belle2::KsfwMoments::m_mm2
double m_mm2[2]
Missing mass squared.
Definition: KsfwMoments.h:142
Belle2::KsfwMoments::m_uf
int m_uf
Flag that specifiies we are using the final state for signal.
Definition: KsfwMoments.h:138
Belle2::KsfwMoments::m_Hso
double m_Hso[2][3][5]
KSFW moments.
Definition: KsfwMoments.h:139
Belle2::legendre
double legendre(const double z, const int i)
Legendre polynomials.
Definition: KsfwMoments.cc:47
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::KsfwMoments::KsfwMoments
KsfwMoments()
Initialize KSFW moments, et, and mm2 to 0.
Definition: KsfwMoments.h:55
Belle2::KsfwMoments::m_et
double m_et[2]
Transverse energy.
Definition: KsfwMoments.h:141
Belle2::KsfwMoments::et
double et(int uf=-1) const
Returns calculated transverse energy.
Definition: KsfwMoments.h:107
Belle2::KsfwMoments::m_Hoo
double m_Hoo[2][5]
KSFW moments.
Definition: KsfwMoments.h:140