GooFit  v2.1.3
DalitzPlotHelpers.cu
Go to the documentation of this file.
1 /*
2 04/05/2016 Christoph Hasse
3 DISCLAIMER:
4 This code is not sufficently tested yet and still under heavy development!
5 
6 Helper functions
7 */
8 
9 #include <goofit/PDFs/physics/DalitzPlotHelpers.h>
10 
11 namespace GooFit {
12 
13 __host__ __device__ bool inDalitz(
14  const fptype &m12, const fptype &m13, const fptype &bigM, const fptype &dm1, const fptype &dm2, const fptype &dm3) {
15  fptype dm1pdm2 = dm1 + dm2;
16  fptype bigMmdm3 = bigM - dm3;
17 
18  bool m12less = (m12 < dm1pdm2 * dm1pdm2) ? false : true;
19  // if (m12 < dm1pdm2*dm1pdm2) return false; // This m12 cannot exist, it's less than the square of the (1,2)
20  // particle mass.
21  bool m12grea = (m12 > bigMmdm3 * bigMmdm3) ? false : true;
22  // if (m12 > bigMmdm3*bigMmdm3) return false; // This doesn't work either, there's no room for an at-rest 3
23  // daughter.
24 
25  fptype sqrtM12 = sqrt(m12);
26  fptype dm11 = dm1 * dm1;
27  fptype dm22 = dm2 * dm2;
28  fptype dm33 = dm3 * dm3;
29 
30  // Calculate energies of 1 and 3 particles in m12 rest frame.
31  // fptype e1star = 0.5 * (m12 - dm2*dm2 + dm1*dm1) / sqrt(m12);
32  fptype e1star = 0.5 * (m12 - dm22 + dm11) / sqrtM12;
33  // fptype e3star = 0.5 * (bigM*bigM - m12 - dm3*dm3) / sqrt(m12);
34  fptype e3star = 0.5 * (bigM * bigM - m12 - dm33) / sqrtM12;
35 
36  fptype rte1mdm11 = sqrt(e1star * e1star - dm11);
37  fptype rte3mdm33 = sqrt(e3star * e3star - dm33);
38 
39  // Bounds for m13 at this value of m12.
40  // fptype minimum = (e1star + e3star)*(e1star + e3star) - pow(sqrt(e1star1 - dm11) + sqrt(e3star*e3star - dm33), 2);
41  fptype minimum = (e1star + e3star) * (e1star + e3star) - (rte1mdm11 + rte3mdm33) * (rte1mdm11 + rte3mdm33);
42 
43  bool m13less = (m13 < minimum) ? false : true;
44  // if (m13 < minimum) return false;
45 
46  // fptype maximum = pow(e1star + e3star, 2) - pow(sqrt(e1star*e1star - dm1*dm1) - sqrt(e3star*e3star - dm3*dm3), 2);
47  fptype maximum = (e1star + e3star) * (e1star + e3star) - (rte1mdm11 - rte3mdm33) * (rte1mdm11 - rte3mdm33);
48  bool m13grea = (m13 > maximum) ? false : true;
49  // if (m13 > maximum) return false;
50 
51  return m12less && m12grea && m13less && m13grea;
52 }
53 
54 __device__ fptype Mass(const fptype *P0) {
55  return sqrt(-P0[0] * P0[0] - P0[1] * P0[1] - P0[2] * P0[2] + P0[3] * P0[3]);
56 }
57 __device__ fptype Mass(const fptype *P0, const fptype *P1) {
58  return sqrt(-((P0[0] + P1[0]) * (P0[0] + P1[0])) - ((P0[1] + P1[1]) * (P0[1] + P1[1]))
59  - ((P0[2] + P1[2]) * (P0[2] + P1[2])) + ((P0[3] + P1[3]) * (P0[3] + P1[3])));
60 }
61 __device__ fptype Mass(const fptype *P0, const fptype *P1, const fptype *P2) {
62  return sqrt(
63  -((P0[0] + P1[0] + P2[0]) * (P0[0] + P1[0] + P2[0])) - ((P0[1] + P1[1] + P2[1]) * (P0[1] + P1[1] + P2[1]))
64  - ((P0[2] + P1[2] + P2[2]) * (P0[2] + P1[2] + P2[2])) + ((P0[3] + P1[3] + P2[3]) * (P0[3] + P1[3] + P2[3])));
65 }
66 __device__ fptype VecDot(const fptype *P0, const fptype *P1) {
67  return (P0[0] * P1[0] + P0[1] + P1[1] + P0[2] + P1[2] + P0[3] + P1[3]);
68 }
69 
70 __device__ void get4Vecs(fptype *Vecs,
71  const unsigned int &constants,
72  const fptype &m12,
73  const fptype &m34,
74  const fptype &cos12,
75  const fptype &cos34,
76  const fptype &phi) {
77  fptype M = functorConstants[constants + 1];
78  fptype m1 = functorConstants[constants + 2];
79  fptype m2 = functorConstants[constants + 3];
80  fptype m3 = functorConstants[constants + 4];
81  fptype m4 = functorConstants[constants + 5];
82  // printf("g4v %f, %f, %f, %f, %f\n",M, m1, m2, m3, m4 );
83  fptype E1 = (m12 * m12 + m1 * m1 - m2 * m2) / (2 * m12);
84  fptype E2 = (m12 * m12 - m1 * m1 + m2 * m2) / (2 * m12);
85  fptype E3 = (m34 * m34 + m3 * m3 - m4 * m4) / (2 * m34);
86  fptype E4 = (m34 * m34 - m3 * m3 + m4 * m4) / (2 * m34);
87  fptype p1 = sqrt(E1 * E1 - m1 * m1);
88  fptype p3 = sqrt(E3 * E3 - m3 * m3);
89  fptype sin12 = sqrt(1 - cos12 * cos12);
90  fptype sin34 = sqrt(1 - cos34 * cos34);
91  fptype ED1 = (M * M + m12 * m12 - m34 * m34) / (2 * m12);
92  fptype PD1 = sqrt(ED1 * ED1 - M * M);
93  fptype beta1 = PD1 / ED1;
94  fptype gamma1 = 1.0 / sqrt(1 - beta1 * beta1);
95  fptype ED2 = (M * M - m12 * m12 + m34 * m34) / (2 * m34);
96  fptype PD2 = sqrt(ED2 * ED2 - M * M);
97  fptype beta2 = -PD2 / ED2;
98  fptype gamma2 = 1.0 / sqrt(1 - beta2 * beta2);
99  // printf("g4v %f, %f, %f, %f, %f\n",E1, m1, E2, p1, p3 );
100 
101  // set X-component
102  Vecs[0] = cos12 * p1;
103  Vecs[4] = -cos12 * p1;
104  Vecs[8] = -cos34 * p3;
105  Vecs[12] = cos34 * p3;
106 
107  // set Y-component
108  Vecs[1] = sin12 * p1;
109  Vecs[5] = -sin12 * p1;
110  Vecs[9] = -sin34 * p3;
111  Vecs[13] = sin34 * p3;
112 
113  // set Z-component
114  Vecs[2] = 0;
115  Vecs[6] = 0;
116  Vecs[10] = 0;
117  Vecs[14] = 0;
118 
119  // set E-component
120  Vecs[3] = E1;
121  Vecs[7] = E2;
122  Vecs[11] = E3;
123  Vecs[15] = E4;
124 
125  fptype tmpE = Vecs[3];
126  fptype tmpX = Vecs[0];
127  Vecs[3] = gamma1 * (tmpE + beta1 * tmpX);
128  Vecs[0] = gamma1 * (tmpX + beta1 * tmpE);
129 
130  tmpE = Vecs[7];
131  tmpX = Vecs[4];
132  Vecs[7] = gamma1 * (tmpE + beta1 * tmpX);
133  Vecs[4] = gamma1 * (tmpX + beta1 * tmpE);
134 
135  tmpE = Vecs[11];
136  tmpX = Vecs[8];
137  Vecs[11] = gamma2 * (tmpE + beta2 * tmpX);
138  Vecs[8] = gamma2 * (tmpX + beta2 * tmpE);
139 
140  tmpE = Vecs[15];
141  tmpX = Vecs[12];
142  Vecs[15] = gamma2 * (tmpE + beta2 * tmpX);
143  Vecs[12] = gamma2 * (tmpX + beta2 * tmpE);
144 
145  // rotation around X-axis of the first two vectors.
146  fptype cosphi = cos(phi);
147  fptype sinphi = sin(phi);
148 
149  // note that Z-component is zero thus rotation is as easy as this:
150  Vecs[2] = sinphi * Vecs[1];
151  Vecs[1] = cosphi * Vecs[1];
152 
153  Vecs[6] = sinphi * Vecs[5];
154  Vecs[5] = cosphi * Vecs[5];
155 }
156 
157 __device__ fptype getmass(const unsigned int &pair,
158  fptype &d1,
159  fptype &d2,
160  const fptype *vecs,
161  const fptype &m1,
162  const fptype &m2,
163  const fptype &m3,
164  const fptype &m4) {
165  const fptype *P1 = vecs;
166  const fptype *P2 = (vecs + 4);
167  const fptype *P3 = (vecs + 8);
168  const fptype *P4 = (vecs + 12);
169  fptype mpair = 0;
170 
171  switch(pair) {
172  case 2:
173  d1 = m1;
174  d2 = m3;
175  mpair = Mass(P1, P3);
176  break;
177 
178  case 3:
179  d1 = m1;
180  d2 = m4;
181  mpair = Mass(P1, P4);
182  break;
183 
184  case 4:
185  d1 = m2;
186  d2 = m3;
187  mpair = Mass(P2, P3);
188  break;
189 
190  case 5:
191  d1 = m2;
192  d2 = m4;
193  mpair = Mass(P2, P4);
194  break;
195 
196  case 6:
197  d1 = Mass(P1, P2);
198  d2 = m3;
199  mpair = Mass(P1, P2, P3);
200  break;
201 
202  case 7:
203  d1 = Mass(P1, P3);
204  d2 = m2;
205  mpair = Mass(P1, P2, P3);
206  break;
207 
208  case 8:
209  d1 = Mass(P2, P3);
210  d2 = m1;
211  mpair = Mass(P1, P2, P3);
212  break;
213 
214  case 9:
215  d1 = Mass(P1, P2);
216  d2 = m4;
217  mpair = Mass(P1, P2, P4);
218  break;
219 
220  case 10:
221  d1 = Mass(P1, P4);
222  d2 = m2;
223  mpair = Mass(P1, P2, P4);
224  break;
225 
226  case 11:
227  d1 = Mass(P2, P4);
228  d2 = m1;
229  mpair = Mass(P1, P2, P4);
230  break;
231 
232  case 12:
233  d1 = Mass(P1, P3);
234  d2 = m4;
235  mpair = Mass(P1, P3, P4);
236  break;
237 
238  case 13:
239  d1 = Mass(P1, P4);
240  d2 = m3;
241  mpair = Mass(P1, P3, P4);
242  break;
243 
244  case 14:
245  d1 = Mass(P3, P4);
246  d2 = m1;
247  mpair = Mass(P1, P3, P4);
248  break;
249 
250  case 15:
251  d1 = Mass(P2, P3);
252  d2 = m4;
253  mpair = Mass(P2, P3, P4);
254  break;
255 
256  case 16:
257  d1 = Mass(P2, P4);
258  d2 = m3;
259  mpair = Mass(P2, P3, P4);
260  break;
261 
262  case 17:
263  d1 = Mass(P3, P4);
264  d2 = m2;
265  mpair = Mass(P2, P3, P4);
266  break;
267  }
268 
269  return mpair;
270 }
271 } // namespace GooFit
Generated on Sat Apr 21 2018 20:13:18 for GooFit by   doxygen 1.8.13