GooFit  v2.1.3
SpinFactors.cu
Go to the documentation of this file.
1 /*
2 04/05/2016 Christoph Hasse
3 The Spinfactors are an adaptation from the MINT implementation, by Jonas Rademacker.
4 
5 DISCLAIMER:
6 This code is not sufficently tested yet and still under heavy development!
7 */
8 #include <goofit/PDFs/physics/SpinFactors.h>
9 #include <goofit/PDFs/physics/SpinHelper.h>
10 
11 namespace GooFit {
12 
13 #define ZEMACH 1
14 
15 __device__ fptype DtoPP1_PtoSP2_StoP3P4(fptype *Vecs, unsigned int *indices) { return 1.0; }
16 
17 __device__ fptype ONE(fptype *Vecs, unsigned int *indices) { return 1.0; }
18 
19 __device__ fptype FF_12_34_L1(fptype *Vecs, unsigned int *indices) {
20  fptype mother_radius = functorConstants[indices[1]];
21  unsigned int p1 = indices[2];
22  unsigned int p2 = indices[3];
23  unsigned int p3 = indices[4];
24  unsigned int p4 = indices[5];
25  gpuLVec P1(Vecs[0 + 4 * p1], Vecs[1 + 4 * p1], Vecs[2 + 4 * p1], Vecs[3 + 4 * p1]);
26  gpuLVec P2(Vecs[0 + 4 * p2], Vecs[1 + 4 * p2], Vecs[2 + 4 * p2], Vecs[3 + 4 * p2]);
27  gpuLVec P3(Vecs[0 + 4 * p3], Vecs[1 + 4 * p3], Vecs[2 + 4 * p3], Vecs[3 + 4 * p3]);
28  gpuLVec P4(Vecs[0 + 4 * p4], Vecs[1 + 4 * p4], Vecs[2 + 4 * p4], Vecs[3 + 4 * p4]);
29 
30  fptype m1 = (P1 + P2).Mag2();
31  fptype m2 = (P3 + P4).Mag2();
32  fptype s = (P1 + P2 + P3 + P4).Mag2();
33  fptype q2 = s / 4. - (m1 + m2) / 2. + (m1 - m2) * (m1 - m2) / (4 * s);
34  fptype z2 = q2 * mother_radius * mother_radius;
35  fptype ff = 1.0 / (1 + z2);
36  // printf("%.5g, %.5g, %.5g, %.5g\n",s,m1,m2,sqrt(ff) );
37  return sqrt(ff);
38 }
39 
40 __device__ fptype FF_12_34_L2(fptype *Vecs, unsigned int *indices) {
41  fptype mother_radius = functorConstants[indices[1]];
42  unsigned int p1 = indices[2];
43  unsigned int p2 = indices[3];
44  unsigned int p3 = indices[4];
45  unsigned int p4 = indices[5];
46  gpuLVec P1(Vecs[0 + 4 * p1], Vecs[1 + 4 * p1], Vecs[2 + 4 * p1], Vecs[3 + 4 * p1]);
47  gpuLVec P2(Vecs[0 + 4 * p2], Vecs[1 + 4 * p2], Vecs[2 + 4 * p2], Vecs[3 + 4 * p2]);
48  gpuLVec P3(Vecs[0 + 4 * p3], Vecs[1 + 4 * p3], Vecs[2 + 4 * p3], Vecs[3 + 4 * p3]);
49  gpuLVec P4(Vecs[0 + 4 * p4], Vecs[1 + 4 * p4], Vecs[2 + 4 * p4], Vecs[3 + 4 * p4]);
50 
51  fptype m1 = (P1 + P2).Mag2();
52  fptype m2 = (P3 + P4).Mag2();
53  fptype s = (P1 + P2 + P3 + P4).Mag2();
54  fptype q2 = s / 4. - (m1 + m2) / 2. + (m1 - m2) * (m1 - m2) / (4 * s);
55  fptype z2 = q2 * mother_radius * mother_radius;
56  fptype ff = 1.0 / (z2 * z2 + 3 * z2 + 9);
57 
58  return sqrt(ff);
59 }
60 
61 __device__ fptype FF_123_4_L1(fptype *Vecs, unsigned int *indices) {
62  fptype mother_radius = functorConstants[indices[1]];
63  unsigned int p1 = indices[2];
64  unsigned int p2 = indices[3];
65  unsigned int p3 = indices[4];
66  unsigned int p4 = indices[5];
67  gpuLVec P1(Vecs[0 + 4 * p1], Vecs[1 + 4 * p1], Vecs[2 + 4 * p1], Vecs[3 + 4 * p1]);
68  gpuLVec P2(Vecs[0 + 4 * p2], Vecs[1 + 4 * p2], Vecs[2 + 4 * p2], Vecs[3 + 4 * p2]);
69  gpuLVec P3(Vecs[0 + 4 * p3], Vecs[1 + 4 * p3], Vecs[2 + 4 * p3], Vecs[3 + 4 * p3]);
70  gpuLVec P4(Vecs[0 + 4 * p4], Vecs[1 + 4 * p4], Vecs[2 + 4 * p4], Vecs[3 + 4 * p4]);
71 
72  fptype m1 = (P1 + P2 + P3).Mag2();
73  fptype m2 = P4.Mag2();
74  fptype s = (P1 + P2 + P3 + P4).Mag2();
75  fptype q2 = s / 4. - (m1 + m2) / 2. + (m1 - m2) * (m1 - m2) / (4 * s);
76  fptype z2 = q2 * mother_radius * mother_radius;
77  fptype ff = 1.0 / (1 + z2);
78 
79  return sqrt(ff);
80 }
81 
82 __device__ fptype FF_123_4_L2(fptype *Vecs, unsigned int *indices) {
83  fptype mother_radius = functorConstants[indices[1]];
84  unsigned int p1 = indices[2];
85  unsigned int p2 = indices[3];
86  unsigned int p3 = indices[4];
87  unsigned int p4 = indices[5];
88  gpuLVec P1(Vecs[0 + 4 * p1], Vecs[1 + 4 * p1], Vecs[2 + 4 * p1], Vecs[3 + 4 * p1]);
89  gpuLVec P2(Vecs[0 + 4 * p2], Vecs[1 + 4 * p2], Vecs[2 + 4 * p2], Vecs[3 + 4 * p2]);
90  gpuLVec P3(Vecs[0 + 4 * p3], Vecs[1 + 4 * p3], Vecs[2 + 4 * p3], Vecs[3 + 4 * p3]);
91  gpuLVec P4(Vecs[0 + 4 * p4], Vecs[1 + 4 * p4], Vecs[2 + 4 * p4], Vecs[3 + 4 * p4]);
92 
93  fptype m1 = (P1 + P2 + P3).Mag2();
94  fptype m2 = P4.Mag2();
95  fptype s = (P1 + P2 + P3 + P4).Mag2();
96  fptype q2 = s / 4. - (m1 + m2) / 2. + (m1 - m2) * (m1 - m2) / (4 * s);
97  fptype z2 = q2 * mother_radius * mother_radius;
98  fptype ff = 1.0 / (z2 * z2 + 3 * z2 + 9);
99 
100  return sqrt(ff);
101 }
102 
103 __device__ fptype DtoPP1_PtoVP2_VtoP3P4(fptype *Vecs, unsigned int *indices) {
104  unsigned int p1 = indices[2];
105  unsigned int p2 = indices[3];
106  unsigned int p3 = indices[4];
107  unsigned int p4 = indices[5];
108  gpuLVec P1(Vecs[0 + 4 * p1], Vecs[1 + 4 * p1], Vecs[2 + 4 * p1], Vecs[3 + 4 * p1]);
109  gpuLVec P2(Vecs[0 + 4 * p2], Vecs[1 + 4 * p2], Vecs[2 + 4 * p2], Vecs[3 + 4 * p2]);
110  gpuLVec P3(Vecs[0 + 4 * p3], Vecs[1 + 4 * p3], Vecs[2 + 4 * p3], Vecs[3 + 4 * p3]);
111  gpuLVec P4(Vecs[0 + 4 * p4], Vecs[1 + 4 * p4], Vecs[2 + 4 * p4], Vecs[3 + 4 * p4]);
112 
113  gpuLVec pV = P3 + P4;
114  gpuLVec qV = P3 - P4;
115 
116 #ifdef ZEMACH
117  gpuLVec pP = pV + P2;
118  gpuLVec qP = pV - P2;
119 
120  ZTspin1 LP(qP, pP, pP.M());
121  ZTspin1 LV(qV, pV, pV.M());
122 
123  return LP.Dot(LV);
124 #else
125 
126  return (P2.Dot(qV) - P2.Dot(pV) * pV.Dot(qV)) / (pV.M() * pV.M());
127 #endif
128 }
129 
130 __device__ fptype DtoV1V2_V1toP1P2_V2toP3P4_S(fptype *Vecs, unsigned int *indices) {
131  unsigned int p1 = indices[2];
132  unsigned int p2 = indices[3];
133  unsigned int p3 = indices[4];
134  unsigned int p4 = indices[5];
135  gpuLVec P1(Vecs[0 + 4 * p1], Vecs[1 + 4 * p1], Vecs[2 + 4 * p1], Vecs[3 + 4 * p1]);
136  gpuLVec P2(Vecs[0 + 4 * p2], Vecs[1 + 4 * p2], Vecs[2 + 4 * p2], Vecs[3 + 4 * p2]);
137  gpuLVec P3(Vecs[0 + 4 * p3], Vecs[1 + 4 * p3], Vecs[2 + 4 * p3], Vecs[3 + 4 * p3]);
138  gpuLVec P4(Vecs[0 + 4 * p4], Vecs[1 + 4 * p4], Vecs[2 + 4 * p4], Vecs[3 + 4 * p4]);
139 
140  // printf("vec%i %.5g, %.5g, %.5g, %.5g\n",0, P1.getX(), P1.getY(), P1.getZ(),P1.getE());
141  // printf("vec%i %.5g, %.5g, %.5g, %.5g\n",1, P2.getX(), P2.getY(), P2.getZ(),P2.getE());
142  // printf("vec%i %.5g, %.5g, %.5g, %.5g\n",2, P3.getX(), P3.getY(), P3.getZ(),P3.getE());
143  // printf("vec%i %.5g, %.5g, %.5g, %.5g\n",3, P4.getX(), P4.getY(), P4.getZ(),P4.getE());
144 
145  gpuLVec pV1 = P1 + P2;
146  gpuLVec qV1 = P1 - P2;
147  gpuLVec pV2 = P3 + P4;
148  gpuLVec qV2 = P3 - P4;
149 
150  fptype MV1 = sqrt(pV1.Dot(pV1));
151  fptype MV2 = sqrt(pV2.Dot(pV2));
152 
153  fptype returnVal
154  = (qV1.Dot(qV2) - qV1.Dot(pV1) * pV1.Dot(qV2) / (MV1 * MV1) - qV1.Dot(pV2) * pV2.Dot(qV2) / (MV2 * MV2)
155  + qV1.Dot(pV1) * pV1.Dot(pV2) * pV2.Dot(qV2) / (MV1 * MV1 * MV2 * MV2));
156  // printf("s1 %.5g; %i,%i,%i,%i\n",returnVal, indices[2], indices[3], indices[4], indices[5]);
157  return returnVal;
158 }
159 
160 __device__ fptype DtoV1V2_V1toP1P2_V2toP3P4_P(fptype *Vecs, unsigned int *indices) {
161  unsigned int p1 = indices[2];
162  unsigned int p2 = indices[3];
163  unsigned int p3 = indices[4];
164  unsigned int p4 = indices[5];
165  gpuLVec P1(Vecs[0 + 4 * p1], Vecs[1 + 4 * p1], Vecs[2 + 4 * p1], Vecs[3 + 4 * p1]);
166  gpuLVec P2(Vecs[0 + 4 * p2], Vecs[1 + 4 * p2], Vecs[2 + 4 * p2], Vecs[3 + 4 * p2]);
167  gpuLVec P3(Vecs[0 + 4 * p3], Vecs[1 + 4 * p3], Vecs[2 + 4 * p3], Vecs[3 + 4 * p3]);
168  gpuLVec P4(Vecs[0 + 4 * p4], Vecs[1 + 4 * p4], Vecs[2 + 4 * p4], Vecs[3 + 4 * p4]);
169 
170  gpuLVec pV1 = P1 + P2;
171  gpuLVec qV1 = P2 - P1;
172  gpuLVec pV2 = P3 + P4;
173  gpuLVec qV2 = P4 - P3;
174 
175  gpuLVec pD = pV1 + pV2;
176  gpuLVec qD = pV1 - pV2;
177 
178 #ifdef ZEMACH
179  fptype MV1 = sqrt(pV1.Dot(pV1));
180  fptype MV2 = sqrt(pV2.Dot(pV2));
181 
182  ZTspin1 LD(qD, pD, pD.M());
183  ZTspin1 LV1(qV1, pV1, MV1);
184  ZTspin1 LV2(qV2, pV2, MV2);
185 
186  return -LeviCivita(pD, LD, LV1).Dot(LV2); // minus gives the same result as MINT3
187 
188 #else
189  return LeviCivita(pD, qD, qV1, qV2);
190 #endif
191 }
192 
193 __device__ fptype DtoV1V2_V1toP1P2_V2toP3P4_D(fptype *Vecs, unsigned int *indices) {
194  unsigned int p1 = indices[2];
195  unsigned int p2 = indices[3];
196  unsigned int p3 = indices[4];
197  unsigned int p4 = indices[5];
198  gpuLVec P1(Vecs[0 + 4 * p1], Vecs[1 + 4 * p1], Vecs[2 + 4 * p1], Vecs[3 + 4 * p1]);
199  gpuLVec P2(Vecs[0 + 4 * p2], Vecs[1 + 4 * p2], Vecs[2 + 4 * p2], Vecs[3 + 4 * p2]);
200  gpuLVec P3(Vecs[0 + 4 * p3], Vecs[1 + 4 * p3], Vecs[2 + 4 * p3], Vecs[3 + 4 * p3]);
201  gpuLVec P4(Vecs[0 + 4 * p4], Vecs[1 + 4 * p4], Vecs[2 + 4 * p4], Vecs[3 + 4 * p4]);
202 
203  gpuLVec pV1 = P1 + P2;
204  gpuLVec qV1 = P1 - P2;
205  gpuLVec pV2 = P3 + P4;
206  gpuLVec qV2 = P3 - P4;
207 
208 #ifdef ZEMACH
209  gpuLVec pD = pV1 + pV2;
210  gpuLVec qD = pV1 - pV2;
211  double mD = pD.M();
212 
213  ZTspin1 tV1(qV1, pV1, pV1.M());
214  ZTspin1 tV2(qV2, pV2, pV2.M());
215  ZTspin2 tD(qD, pD, mD);
216  double returnVal = tV1.Contract(tD.Contract(tV2));
217  // printf("%f, %f, %f, %f,",P1.GetX(), P1.GetY(), P1.GetZ(), P1.GetE() );
218  // printf("%f, %f, %f, %f,",P2.GetX(), P2.GetY(), P2.GetZ(), P2.GetE() );
219  // printf("%f, %f, %f, %f,",P3.GetX(), P3.GetY(), P3.GetZ(), P3.GetE() );
220  // printf("%f, %f, %f, %f,",P4.GetX(), P4.GetY(), P4.GetZ(), P4.GetE() );
221  // printf("%f\n",returnVal );
222  return returnVal;
223 #else
224 
225  fptype MV1 = sqrt(pV1.Dot(pV1));
226  fptype MV2 = sqrt(pV2.Dot(pV2));
227  fptype returnVal = (qV1.Dot(pV2) - qV1.Dot(pV1) * pV1.Dot(pV2) / (MV1 * MV1))
228  * (qV2.Dot(pV1) - qV2.Dot(pV2) * pV2.Dot(pV1) / (MV2 * MV2));
229 
230  return returnVal;
231 #endif
232 }
233 
234 __device__ fptype DtoV1P1_V1toV2P2_V2toP3P4(fptype *Vecs, unsigned int *indices) {
235  unsigned int p1 = indices[2];
236  unsigned int p2 = indices[3];
237  unsigned int p3 = indices[4];
238  unsigned int p4 = indices[5];
239  gpuLVec P1(Vecs[0 + 4 * p1], Vecs[1 + 4 * p1], Vecs[2 + 4 * p1], Vecs[3 + 4 * p1]);
240  gpuLVec P2(Vecs[0 + 4 * p2], Vecs[1 + 4 * p2], Vecs[2 + 4 * p2], Vecs[3 + 4 * p2]);
241  gpuLVec P3(Vecs[0 + 4 * p3], Vecs[1 + 4 * p3], Vecs[2 + 4 * p3], Vecs[3 + 4 * p3]);
242  gpuLVec P4(Vecs[0 + 4 * p4], Vecs[1 + 4 * p4], Vecs[2 + 4 * p4], Vecs[3 + 4 * p4]);
243 
244  gpuLVec pV1 = P2 + P3 + P4;
245  gpuLVec pV2 = P3 + P4;
246  gpuLVec qV1 = (P3 + P4) - P2;
247  gpuLVec qV2 = P3 - P4;
248 
249 #ifdef ZEMACH
250  double MV1 = pV1.M();
251  double MV2 = pV2.M();
252 
253  gpuLVec pD = pV1 + P1;
254  gpuLVec qD = pV1 - P1;
255 
256  ZTspin1 LD(qD, pD, pD.M());
257  ZTspin1 LV1(qV1, pV1, MV1);
258  ZTspin1 LV2(qV2, pV2, MV2);
259  SpinSumV PV1(pV1, MV1);
260 
261  gpuLVec tmp = PV1.Dot(LD);
262 
263  return LeviCivita(pV1, LV1, LV2).Dot(tmp);
264 #else
265 
266  // printf("%f, %f, %f, %f\n",P1.getX(), P1.getY(), P1.getZ(), P1.getE() );
267  // printf("%f, %f, %f, %f\n",P2.getX(), P2.getY(), P2.getZ(), P2.getE() );
268  // printf("%f, %f, %f, %f\n",P3.getX(), P3.getY(), P3.getZ(), P3.getE() );
269  // printf("%f, %f, %f, %f\n",P4.getX(), P4.getY(), P4.getZ(), P4.getE() );
270 
271  fptype returnVal = LeviCivita(pV1, qV1, P1, qV2);
272  return returnVal;
273 #endif
274 }
275 
276 __device__ fptype DtoVS_VtoP1P2_StoP3P4(fptype *Vecs, unsigned int *indices) {
277  unsigned int p1 = indices[2];
278  unsigned int p2 = indices[3];
279  unsigned int p3 = indices[4];
280  unsigned int p4 = indices[5];
281  gpuLVec P1(Vecs[0 + 4 * p1], Vecs[1 + 4 * p1], Vecs[2 + 4 * p1], Vecs[3 + 4 * p1]);
282  gpuLVec P2(Vecs[0 + 4 * p2], Vecs[1 + 4 * p2], Vecs[2 + 4 * p2], Vecs[3 + 4 * p2]);
283  gpuLVec P3(Vecs[0 + 4 * p3], Vecs[1 + 4 * p3], Vecs[2 + 4 * p3], Vecs[3 + 4 * p3]);
284  gpuLVec P4(Vecs[0 + 4 * p4], Vecs[1 + 4 * p4], Vecs[2 + 4 * p4], Vecs[3 + 4 * p4]);
285 
286  gpuLVec pS = P3 + P4;
287  gpuLVec pV = P1 + P2;
288  gpuLVec qV = P1 - P2;
289 
290 #ifdef ZEMACH
291  gpuLVec pD = pV + pS;
292  gpuLVec qD = pV - pS;
293  ZTspin1 LD(qD, pD, pD.M());
294  ZTspin1 LV(qV, pV, pV.M());
295  return (LD.Dot(LV));
296 #else
297 
298  // printf("%f, %f, %f, %f\n",P1.getX(), P1.getY(), P1.getZ(), P1.getE() );
299  // printf("%f, %f, %f, %f\n",P2.getX(), P2.getY(), P2.getZ(), P2.getE() );
300  // printf("%f, %f, %f, %f\n",P3.getX(), P3.getY(), P3.getZ(), P3.getE() );
301  // printf("%f, %f, %f, %f\n",P4.getX(), P4.getY(), P4.getZ(), P4.getE() );
302 
303  fptype MV = sqrt(pV.Dot(pV));
304 
305  fptype returnVal = (pS.Dot(qV) - pS.Dot(pV) * pV.Dot(qV) / (MV * MV));
306  return returnVal;
307 #endif
308 }
309 
310 __device__ fptype DtoAP1_AtoSP2_StoP3P4(fptype *Vecs, unsigned int *indices) {
311  unsigned int p1 = indices[2];
312  unsigned int p2 = indices[3];
313  unsigned int p3 = indices[4];
314  unsigned int p4 = indices[5];
315  gpuLVec P1(Vecs[0 + 4 * p1], Vecs[1 + 4 * p1], Vecs[2 + 4 * p1], Vecs[3 + 4 * p1]);
316  gpuLVec P2(Vecs[0 + 4 * p2], Vecs[1 + 4 * p2], Vecs[2 + 4 * p2], Vecs[3 + 4 * p2]);
317  gpuLVec P3(Vecs[0 + 4 * p3], Vecs[1 + 4 * p3], Vecs[2 + 4 * p3], Vecs[3 + 4 * p3]);
318  gpuLVec P4(Vecs[0 + 4 * p4], Vecs[1 + 4 * p4], Vecs[2 + 4 * p4], Vecs[3 + 4 * p4]);
319 
320  gpuLVec pS = P3 + P4;
321  gpuLVec pA = P2 + pS;
322  gpuLVec qA = pS - P2;
323  gpuLVec pD = pA + P1;
324  gpuLVec qD = pA - P1;
325 
326 #ifdef ZEMACH
327  ZTspin1 LD(qD, pD, pD.M());
328  ZTspin1 LA(qA, pA, pA.M());
329  return (LD.Dot(LA));
330 #else
331 
332  fptype MA = sqrt(pA.Dot(pA));
333  fptype returnVal = (P1.Dot(qA) - P1.Dot(pA) * pA.Dot(qA) / (MA * MA));
334  return returnVal;
335 #endif
336 }
337 
338 __device__ fptype DtoAP1_AtoVP2_VtoP3P4(fptype *Vecs, unsigned int *indices) {
339  unsigned int p1 = indices[2];
340  unsigned int p2 = indices[3];
341  unsigned int p3 = indices[4];
342  unsigned int p4 = indices[5];
343  gpuLVec P1(Vecs[0 + 4 * p1], Vecs[1 + 4 * p1], Vecs[2 + 4 * p1], Vecs[3 + 4 * p1]);
344  gpuLVec P2(Vecs[0 + 4 * p2], Vecs[1 + 4 * p2], Vecs[2 + 4 * p2], Vecs[3 + 4 * p2]);
345  gpuLVec P3(Vecs[0 + 4 * p3], Vecs[1 + 4 * p3], Vecs[2 + 4 * p3], Vecs[3 + 4 * p3]);
346  gpuLVec P4(Vecs[0 + 4 * p4], Vecs[1 + 4 * p4], Vecs[2 + 4 * p4], Vecs[3 + 4 * p4]);
347 
348  gpuLVec pV = P3 + P4;
349  gpuLVec qV = P3 - P4;
350  gpuLVec pA = P2 + pV;
351  gpuLVec pD = P1 + pA;
352  gpuLVec qD = pA - P1;
353 
354 #ifdef ZEMACH
355  ZTspin1 LB(qD, pD, pD.M());
356  ZTspin1 LV(qV, pV, pV.M());
357  SpinSumV PA(pA, pA.M());
358 
359  gpuLVec tmp = PA.Dot(LV);
360  return -(LB.Dot(tmp)); // minus to be equal to MINT3
361 #else
362  gpuLVec p0 = P1;
363 
364  fptype MA = sqrt(pA.Dot(pA));
365  fptype MV = sqrt(pV.Dot(pV));
366  fptype returnVal = P1.Dot(qV) - p0.Dot(pA) * pA.Dot(qV) / (MA * MA) - p0.Dot(pV) * pV.Dot(qV) / (MV * MV)
367  + p0.Dot(pA) * pA.Dot(pV) * pV.Dot(qV) / (MA * MA * MV * MV);
368  // printf("spin %.7g\n",returnVal );
369  return returnVal;
370 #endif
371 }
372 
373 __device__ fptype DtoAP1_AtoVP2Dwave_VtoP3P4(fptype *Vecs, unsigned int *indices) {
374  unsigned int p1 = indices[2];
375  unsigned int p2 = indices[3];
376  unsigned int p3 = indices[4];
377  unsigned int p4 = indices[5];
378  gpuLVec P1(Vecs[0 + 4 * p1], Vecs[1 + 4 * p1], Vecs[2 + 4 * p1], Vecs[3 + 4 * p1]);
379  gpuLVec P2(Vecs[0 + 4 * p2], Vecs[1 + 4 * p2], Vecs[2 + 4 * p2], Vecs[3 + 4 * p2]);
380  gpuLVec P3(Vecs[0 + 4 * p3], Vecs[1 + 4 * p3], Vecs[2 + 4 * p3], Vecs[3 + 4 * p3]);
381  gpuLVec P4(Vecs[0 + 4 * p4], Vecs[1 + 4 * p4], Vecs[2 + 4 * p4], Vecs[3 + 4 * p4]);
382 
383  gpuLVec pV = P3 + P4;
384  gpuLVec qV = P3 - P4;
385  gpuLVec pA = P2 + pV;
386  gpuLVec qA = P2 - pV;
387  gpuLVec pD = P1 + pA;
388  gpuLVec qD = pA - P1;
389 
390 #ifdef ZEMACH
391  ZTspin1 LD(qD, pD, pD.M());
392  ZTspin1 LV(qV, pV, pV.M());
393  ZTspin2 LA(qA, pA, pA.M());
394  gpuLVec tmp = LA.Contract(LV);
395 
396  return (LD.Dot(tmp));
397 #else
398 
399  fptype MA = sqrt(pA.Dot(pA));
400  fptype MV = sqrt(pV.Dot(pV));
401 
402  fptype returnVal
403  = P1.Dot((qA - qA.Dot(pA) * pA * (1. / (MA * MA)))) * (qV - qV.Dot(pV) * pV * (1. / (MV * MV))).Dot(pA);
404 
405  // printf("spin %.7g\n",returnVal );
406  return returnVal;
407 #endif
408 }
409 
410 __device__ fptype DtoTP1_TtoVP2_VtoP3P4(fptype *Vecs, unsigned int *indices) {
411  unsigned int p1 = indices[2];
412  unsigned int p2 = indices[3];
413  unsigned int p3 = indices[4];
414  unsigned int p4 = indices[5];
415  gpuLVec P1(Vecs[0 + 4 * p1], Vecs[1 + 4 * p1], Vecs[2 + 4 * p1], Vecs[3 + 4 * p1]);
416  gpuLVec P2(Vecs[0 + 4 * p2], Vecs[1 + 4 * p2], Vecs[2 + 4 * p2], Vecs[3 + 4 * p2]);
417  gpuLVec P3(Vecs[0 + 4 * p3], Vecs[1 + 4 * p3], Vecs[2 + 4 * p3], Vecs[3 + 4 * p3]);
418  gpuLVec P4(Vecs[0 + 4 * p4], Vecs[1 + 4 * p4], Vecs[2 + 4 * p4], Vecs[3 + 4 * p4]);
419 
420  gpuLVec pV = P3 + P4;
421  gpuLVec qV = P3 - P4;
422  gpuLVec pT = P2 + pV;
423  gpuLVec qT = pV - P2;
424  gpuLVec pD = P1 + pT;
425  gpuLVec qD = pT - P1;
426 
427  ZTspin2 t2T(qT, pT, pT.M());
428  ZTspin1 tV(qV, pV, pV.M());
429 
430 #ifdef ZEMACH
431  ZTspin1 tD(qD, pD, pD.M());
432  ZTspin1 t1T(qT, pT, pT.M());
433 
434  SpinSumV P1T(pT, pT.M());
435  SpinSumT P2T(pT, pT.M());
436 
437  gpuLVec tmp = LeviCivita(t1T, pT, P1T.Dot(tV));
438  return P2T.Sandwich(tD, tD, tmp, t1T) + 1. / 3. * tD.Dot(tD) / (pT.Dot(pT)) * P2T.Sandwich(pD, pD, tmp, t1T);
439 #else
440 
441  fptype returnVal = LeviCivita(pD, qD, DT, tV);
442 
443  // printf("spin %.7g\n",returnVal );
444  return returnVal;
445 #endif
446 }
447 
448 __device__ spin_function_ptr ptr_to_DtoPP1_PtoSP2_StoP3P4 = DtoPP1_PtoSP2_StoP3P4;
449 __device__ spin_function_ptr ptr_to_DtoPP1_PtoVP2_VtoP3P4 = DtoPP1_PtoVP2_VtoP3P4;
450 __device__ spin_function_ptr ptr_to_DtoV1V2_V1toP1P2_V2toP3P4_S = DtoV1V2_V1toP1P2_V2toP3P4_S;
451 __device__ spin_function_ptr ptr_to_DtoV1V2_V1toP1P2_V2toP3P4_P = DtoV1V2_V1toP1P2_V2toP3P4_P;
452 __device__ spin_function_ptr ptr_to_DtoV1V2_V1toP1P2_V2toP3P4_D = DtoV1V2_V1toP1P2_V2toP3P4_D;
453 __device__ spin_function_ptr ptr_to_DtoVS_VtoP1P2_StoP3P4 = DtoVS_VtoP1P2_StoP3P4;
454 __device__ spin_function_ptr ptr_to_DtoV1P1_V1toV2P2_V2toP3P4 = DtoV1P1_V1toV2P2_V2toP3P4;
455 __device__ spin_function_ptr ptr_to_DtoAP1_AtoSP2_StoP3P4 = DtoAP1_AtoSP2_StoP3P4;
456 __device__ spin_function_ptr ptr_to_DtoAP1_AtoVP2_VtoP3P4 = DtoAP1_AtoVP2_VtoP3P4;
457 __device__ spin_function_ptr ptr_to_DtoAP1_AtoVP2Dwave_VtoP3P4 = DtoAP1_AtoVP2Dwave_VtoP3P4;
458 __device__ spin_function_ptr ptr_to_DtoTP1_TtoVP2_VtoP3P4 = DtoTP1_TtoVP2_VtoP3P4;
459 
460 __device__ spin_function_ptr ptr_to_FF_12_34_L1 = FF_12_34_L1;
461 __device__ spin_function_ptr ptr_to_FF_12_34_L2 = FF_12_34_L2;
462 __device__ spin_function_ptr ptr_to_FF_123_4_L1 = FF_123_4_L1;
463 __device__ spin_function_ptr ptr_to_FF_123_4_L2 = FF_123_4_L2;
464 __device__ spin_function_ptr ptr_to_ONE = ONE;
465 
466 SpinFactor::SpinFactor(
467  std::string name, SF_4Body SF, unsigned int P0, unsigned int P1, unsigned int P2, unsigned int P3)
468  : GooPdf(name)
469  , _SF(SF)
470  , _P0(P0)
471  , _P1(P1)
472  , _P2(P2)
473  , _P3(P3) {
474  std::vector<unsigned int> pindices;
475  pindices.push_back(0); // dummy for index to constants.
476  pindices.push_back(P0);
477  pindices.push_back(P1);
478  pindices.push_back(P2);
479  pindices.push_back(P3);
480 
481  switch(SF) {
482  case SF_4Body::DtoPP1_PtoSP2_StoP3P4:
483  GET_FUNCTION_ADDR(ptr_to_DtoPP1_PtoSP2_StoP3P4);
484  break;
485 
486  case SF_4Body::DtoPP1_PtoVP2_VtoP3P4:
487  GET_FUNCTION_ADDR(ptr_to_DtoPP1_PtoVP2_VtoP3P4);
488  break;
489 
490  case SF_4Body::DtoV1V2_V1toP1P2_V2toP3P4_S:
491  GET_FUNCTION_ADDR(ptr_to_DtoV1V2_V1toP1P2_V2toP3P4_S);
492  break;
493 
494  case SF_4Body::DtoV1V2_V1toP1P2_V2toP3P4_P:
495  GET_FUNCTION_ADDR(ptr_to_DtoV1V2_V1toP1P2_V2toP3P4_P);
496  break;
497 
498  case SF_4Body::DtoV1V2_V1toP1P2_V2toP3P4_D:
499  GET_FUNCTION_ADDR(ptr_to_DtoV1V2_V1toP1P2_V2toP3P4_D);
500  break;
501 
502  case SF_4Body::DtoAP1_AtoVP2_VtoP3P4:
503  GET_FUNCTION_ADDR(ptr_to_DtoAP1_AtoVP2_VtoP3P4);
504  break;
505 
506  case SF_4Body::DtoAP1_AtoVP2Dwave_VtoP3P4:
507  GET_FUNCTION_ADDR(ptr_to_DtoAP1_AtoVP2Dwave_VtoP3P4);
508  break;
509 
510  case SF_4Body::DtoVS_VtoP1P2_StoP3P4:
511  GET_FUNCTION_ADDR(ptr_to_DtoVS_VtoP1P2_StoP3P4);
512  break;
513 
514  case SF_4Body::DtoV1P1_V1toV2P2_V2toP3P4:
515  GET_FUNCTION_ADDR(ptr_to_DtoV1P1_V1toV2P2_V2toP3P4);
516  break;
517 
518  case SF_4Body::DtoAP1_AtoSP2_StoP3P4:
519  GET_FUNCTION_ADDR(ptr_to_DtoAP1_AtoSP2_StoP3P4);
520  break;
521 
522  case SF_4Body::DtoTP1_TtoVP2_VtoP3P4:
523  GET_FUNCTION_ADDR(ptr_to_DtoTP1_TtoVP2_VtoP3P4);
524  break;
525 
526  case SF_4Body::FF_12_34_L1:
527  GET_FUNCTION_ADDR(ptr_to_FF_12_34_L1);
528  break;
529 
530  case SF_4Body::FF_12_34_L2:
531  GET_FUNCTION_ADDR(ptr_to_FF_12_34_L2);
532  break;
533 
534  case SF_4Body::FF_123_4_L1:
535  GET_FUNCTION_ADDR(ptr_to_FF_123_4_L1);
536  break;
537 
538  case SF_4Body::FF_123_4_L2:
539  GET_FUNCTION_ADDR(ptr_to_FF_123_4_L2);
540  break;
541 
542  case SF_4Body::ONE:
543  GET_FUNCTION_ADDR(ptr_to_ONE);
544  break;
545 
546  default:
547  std::cout << "No Spinfunction implemented for that kind." << std::endl;
548  exit(0);
549  break;
550  }
551 
552  initialize(pindices);
553 }
554 
555 } // namespace GooFit
Generated on Sat Apr 21 2018 20:13:18 for GooFit by   doxygen 1.8.13