Belle II Software  release-06-01-15
PhotosUtilities.cc
1 #include "PhotosUtilities.h"
2 #include <cstdlib>
3 #include <cstdio>
4 using std::max;
5 
6 namespace Photospp {
7 
8  namespace PhotosUtilities {
9 
10 
11  void fill_val(int beg, int end, double* array, double value)
12  {
13  for (int i = beg; i < end; i++)
14  array[i] = value;
15  }
16 
17 
18 //----------------------------------------------------------------------
19 //
20 // PHOEPS: PHOeps vector product (normalized to unity)
21 //
22 // Purpose: calculates vector product, then normalizes its length.
23 // used to generate orthogonal vectors, i.e. to
24 // generate polarimetric vectors for photons.
25 //
26 // Input Parameters: VEC1,VEC2 - input 4-vectors
27 //
28 // Output Parameters: EPS - normalized 4-vector, orthogonal to
29 // VEC1 and VEC2
30 //
31 // Author(s): Z. Was, P.Golonka Created at: 19/01/05
32 // Last Update: 10/06/13
33 //
34 //----------------------------------------------------------------------
35 
36  void PHOEPS(double vec1[4], double vec2[4], double eps[4])
37  {
38  double xn;
39  int j = 1; // convention of indices of Riemann space must be preserved.
40 
41  eps[1 - j] = vec1[2 - j] * vec2[3 - j] - vec1[3 - j] * vec2[2 - j];
42  eps[2 - j] = vec1[3 - j] * vec2[1 - j] - vec1[1 - j] * vec2[3 - j];
43  eps[3 - j] = vec1[1 - j] * vec2[2 - j] - vec1[2 - j] * vec2[1 - j];
44  eps[4 - j] = 0.0;
45 
46  xn = sqrt(eps[1 - j] * eps[1 - j] + eps[2 - j] * eps[2 - j] + eps[3 - j] * eps[3 - j]);
47 
48  eps[1 - j] = eps[1 - j] / xn;
49  eps[2 - j] = eps[2 - j] / xn;
50  eps[3 - j] = eps[3 - j] / xn;
51 
52  }
53 
54 
55 
56 //----------------------------------------------------------------------
57 //
58 // PHOTOS: PHOton radiation in decays function for SPIn determina-
59 // tion
60 //
61 // Purpose: Calculate the spin of particle with code IDHEP. The
62 // code of the particle is defined by the Particle Data
63 // Group in Phys. Lett. B204 (1988) 1.
64 //
65 // Input Parameter: IDHEP
66 //
67 // Output Parameter: Funtion value = spin of particle with code
68 // IDHEP
69 //
70 // Author(s): E. Barberio and B. van Eijk Created at: 29/11/89
71 // Last update: 10/06/13
72 //
73 //----------------------------------------------------------------------
74  double PHOSPI(int idhep)
75  {
76  static double SPIN[100] = { 0 };
77  static int j = 0;
78  //--
79  //-- Array 'SPIN' contains the spin of the first 100 particles accor-
80  //-- ding to the PDG particle code...
81 
82  if (j == 0) { // initialization
83  j = 1;
84  fill_val(0 , 8, SPIN, 0.5);
85  fill_val(8 , 9, SPIN, 1.0);
86  fill_val(9 , 10, SPIN, 0.0);
87  fill_val(10, 18, SPIN, 0.5);
88  fill_val(18, 20, SPIN, 0.0);
89  fill_val(20, 24, SPIN, 1.0);
90  fill_val(24, 100, SPIN, 0.0);
91  }
92 
93  int idabs = abs(idhep);
94  //--
95  //-- Spin of quark, lepton, boson etc....
96  if (idabs - 1 < 100) return SPIN[idabs - 1];
97 
98  //-- ...other particles, however...
99  double xx = ((idabs % 10) - 1.0) / 2.0;
100  //--
101  //-- ...K_short and K_long are special !!
102  xx = max(xx, 0.0);
103  return xx;
104  }
105 
106 //----------------------------------------------------------------------
107 //
108 // PHOTOS: PHOton radiation in decays CHArge determination
109 //
110 // Purpose: Calculate the charge of particle with code IDHEP. The
111 // code of the particle is defined by the Particle Data
112 // Group in Phys. Lett. B204 (1988) 1.
113 //
114 // Input Parameter: IDHEP
115 //
116 // Output Parameter: Funtion value = charge of particle with code
117 // IDHEP
118 //
119 // Author(s): E. Barberio and B. van Eijk Created at: 29/11/89
120 // Last update: 11/06/13
121 //
122 //----------------------------------------------------------------------
123  double PHOCHA(int idhep)
124  {
125  static double CHARGE[101] = { 0 };
126  static int j = 0;
127  //--
128  //-- Array 'SPIN' contains the spin of the first 100 particles accor-
129  //-- ding to the PDG particle code...
130 
131  if (j == 0) { // initialization
132  j = 1;
133  fill_val(0 , 1, CHARGE, 0.0);
134  fill_val(1 , 2, CHARGE, -0.3333333333);
135  fill_val(2 , 3, CHARGE, 0.6666666667);
136  fill_val(3 , 4, CHARGE, -0.3333333333);
137  fill_val(4 , 5, CHARGE, 0.6666666667);
138  fill_val(5 , 6, CHARGE, -0.3333333333);
139  fill_val(6 , 7, CHARGE, 0.6666666667);
140  fill_val(7 , 8, CHARGE, -0.3333333333);
141  fill_val(8 , 9, CHARGE, 0.6666666667);
142  fill_val(9 , 11, CHARGE, 0.0);
143  fill_val(11 , 12, CHARGE, -1.0);
144  fill_val(12 , 13, CHARGE, 0.0);
145  fill_val(13 , 14, CHARGE, -1.0);
146  fill_val(14, 15, CHARGE, 0.0);
147  fill_val(15 , 16, CHARGE, -1.0);
148  fill_val(16, 17, CHARGE, 0.0);
149  fill_val(17 , 18, CHARGE, -1.0);
150  fill_val(18, 24, CHARGE, 0.0);
151  fill_val(24, 25, CHARGE, 1.0);
152  fill_val(25, 37, CHARGE, 0.0);
153  fill_val(37, 38, CHARGE, 1.0);
154  fill_val(38, 101, CHARGE, 0.0);
155  }
156 
157  int idabs = abs(idhep);
158  double phoch = 0.0;
159 
160  //--
161  //-- Charge of quark, lepton, boson etc....
162  if (idabs <= 100) phoch = CHARGE[idabs];
163  else {
164  int Q3 = idabs / 1000 % 10;
165  int Q2 = idabs / 100 % 10;
166  int Q1 = idabs / 10 % 10;
167  if (Q3 == 0) {
168  //--
169  //-- ...meson...
170  if (Q2 % 2 == 0) phoch = CHARGE[Q2] - CHARGE[Q1];
171  else phoch = CHARGE[Q1] - CHARGE[Q2];
172  } else {
173  //--
174  //-- ...diquarks or baryon.
175  phoch = CHARGE[Q1] + CHARGE[Q2] + CHARGE[Q3];
176  }
177  }
178  //--
179  //-- Find the sign of the charge...
180  if (idhep < 0.0) phoch = -phoch;
181  if (phoch * phoch < 0.000001) phoch = 0.0;
182 
183  return phoch;
184  }
185 
186 
187 
188 
189 //----------------------------------------------------------------------
190 //
191 // PHOTOS: PHOton radiation in decays calculation of TRIangle fie
192 //
193 // Purpose: Calculation of triangle function for phase space.
194 //
195 // Input Parameters: A, B, C (Virtual) particle masses.
196 //
197 // Output Parameter: Function value =
198 // SQRT(LAMBDA(A**2,B**2,C**2))/(2*A)
199 //
200 // Author(s): B. van Eijk Created at: 15/11/89
201 // Last Update: 12/06/13
202 //
203 //----------------------------------------------------------------------
204  double PHOTRI(double A, double B, double C)
205  {
206  double DA, DB, DC, DAPB, DAMB, DTRIAN;
207  DA = A;
208  DB = B;
209  DC = C;
210  DAPB = DA + DB;
211  DAMB = DA - DB;
212  DTRIAN = sqrt((DAMB - DC) * (DAPB + DC) * (DAMB + DC) * (DAPB - DC));
213  return DTRIAN / (DA + DA);
214  }
215 //----------------------------------------------------------------------
216 //
217 // PHOTOS: PHOton radiation in decays calculation of ANgle '1'
218 //
219 // Purpose: Calculate angle from X and Y
220 //
221 // Input Parameters: X, Y
222 //
223 // Output Parameter: Function value
224 //
225 // Author(s): S. Jadach Created at: 01/01/89
226 // B. van Eijk Last Update: 12/06/13
227 //
228 //----------------------------------------------------------------------
229  double PHOAN1(double X, double Y)
230  {
231 
232  double phoan1 = 0.0;
233 
234  static double PI = 3.14159265358979324, TWOPI = 6.28318530717958648;
235 
236  if (fabs(Y) < fabs(X)) {
237  phoan1 = atan(fabs(Y / X));
238  if (X < 0.0) phoan1 = PI - phoan1;
239  } else phoan1 = acos(X / sqrt(X * X + Y * Y));
240  //
241  if (Y < 0.0) phoan1 = TWOPI - phoan1;
242  return phoan1;
243 
244  }
245 
246 //----------------------------------------------------------------------
247 //
248 // PHOTOS: PHOton radiation in decays calculation of ANgle '2'
249 //
250 // Purpose: Calculate angle from X and Y
251 //
252 // Input Parameters: X, Y
253 //
254 // Output Parameter: Function value
255 //
256 // Author(s): S. Jadach Created at: 01/01/89
257 // B. van Eijk Last Update: 12/06/13
258 //
259 //----------------------------------------------------------------------
260  double PHOAN2(double X, double Y)
261  {
262 
263  double phoan2 = 0.0;
264 
265  static double PI = 3.14159265358979324; //, TWOPI=6.28318530717958648;
266 
267  if (fabs(Y) < fabs(X)) {
268  phoan2 = atan(fabs(Y / X));
269  if (X < 0.0) phoan2 = PI - phoan2;
270  } else phoan2 = acos(X / sqrt(X * X + Y * Y));
271  return phoan2;
272  }
273 
274 //----------------------------------------------------------------------
275 //
276 // PHOTOS: PHOton radiation in decays ROtation routine '2'
277 //
278 // Purpose: Rotate x and z components of vector PVEC around angle
279 // 'ANGLE'.
280 //
281 // Input Parameters: ANGLE, PVEC
282 //
283 // Output Parameter: PVEC
284 //
285 // Author(s): S. Jadach Created at: 01/01/89
286 // B. van Eijk Last Update: 12/06/13
287 //
288 //----------------------------------------------------------------------
289  void PHORO2(double ANGLE, double PVEC[4])
290  {
291  int j = 1; // convention of indices of Riemann space must be preserved.
292 
293  double CS, SN;
294  CS = cos(ANGLE) * PVEC[1 - j] + sin(ANGLE) * PVEC[3 - j];
295  SN = -sin(ANGLE) * PVEC[1 - j] + cos(ANGLE) * PVEC[3 - j];
296  PVEC[1 - j] = CS;
297  PVEC[3 - j] = SN;
298  }
299 
300 //----------------------------------------------------------------------
301 //
302 // PHOTOS: PHOton radiation in decays ROtation routine '3'
303 //
304 // Purpose: Rotate x and y components of vector PVEC around angle
305 // 'ANGLE'.
306 //
307 // Input Parameters: ANGLE, PVEC
308 //
309 // Output Parameter: PVEC
310 //
311 // Author(s): S. Jadach RO Created at: 01/01/89
312 // B. van Eijk Last Update: 12/06/13
313 //
314 //----------------------------------------------------------------------
315  void PHORO3(double ANGLE, double PVEC[4])
316  {
317  int j = 1; // convention of indices of Riemann space must be preserved.
318  double CS, SN;
319  CS = cos(ANGLE) * PVEC[1 - j] - sin(ANGLE) * PVEC[2 - j];
320  SN = sin(ANGLE) * PVEC[1 - j] + cos(ANGLE) * PVEC[2 - j];
321  PVEC[1 - j] = CS;
322  PVEC[2 - j] = SN;
323  }
324 
325 //----------------------------------------------------------------------
326 //
327 //
328 // PHOB: PHotosBoost
329 //
330 // Purpose: Boosts VEC to (MODE=1) rest frame of PBOOS1;
331 // or back (MODE=1)
332 //
333 // Input Parameters: MODE,PBOOS1,VEC
334 //
335 // Output Parameters: VEC
336 //
337 // Author(s): Created at: 08/12/05
338 // Z. Was Last Update: 13/06/13
339 //
340 //----------------------------------------------------------------------
341 
342  void PHOB(int MODE, double PBOOS1[4], double vec[4])
343  {
344  double BET1[3], GAM1, PB;
345  static int j0 = 1;
346  int J;
347 
348 
349  PB = sqrt(PBOOS1[4 - j0] * PBOOS1[4 - j0] - PBOOS1[3 - j0] * PBOOS1[3 - j0] - PBOOS1[2 - j0] * PBOOS1[2 - j0] - PBOOS1[1 - j0] *
350  PBOOS1[1 - j0]);
351  for (J = 1; J < 4; J++) {
352  if (MODE == 1) BET1[J - j0] = -PBOOS1[J - j0] / PB;
353  else BET1[J - j0] = PBOOS1[J - j0] / PB;
354  }
355 
356  GAM1 = PBOOS1[4 - j0] / PB;
357 
358  //--
359  //-- Boost vector
360 
361  PB = BET1[1 - j0] * vec[1 - j0] + BET1[2 - j0] * vec[2 - j0] + BET1[3 - j0] * vec[3 - j0];
362 
363  for (J = 1; J < 4; J++) vec[J - j0] = vec[J - j0] + BET1[J - j0] * (vec[4 - j0] + PB / (GAM1 + 1.0));
364  vec[4 - j0] = GAM1 * vec[4 - j0] + PB;
365  //--
366  }
367 
368 
369 // *******************************
370 // Boost along arbitrary axis (as implemented by Ronald Kleiss).
371 // The method is described in book of Bjorken and Drell
372 // p boosted into r from actual frame to rest frame of q
373 // forth (mode = 1) or back (mode = -1).
374 // q must be a timelike, p may be arbitrary.
375  void bostdq(int mode, double qq[4], double pp[4], double r[4])
376  {
377  double q[4], p[4], amq, fac;
378  static int i = 1;
379  int k;
380 
381  for (k = 1; k <= 4; k++) {
382  p[k - i] = pp[k - i];
383  q[k - i] = qq[k - i];
384  }
385  amq = sqrt(q[4 - i] * q[4 - i] - q[1 - i] * q[1 - i] - q[2 - i] * q[2 - i] - q[3 - i] * q[3 - i]);
386 
387  if (mode == -1) {
388  r[4 - i] = (p[1 - i] * q[1 - i] + p[2 - i] * q[2 - i] + p[3 - i] * q[3 - i] + p[4 - i] * q[4 - i]) / amq;
389  fac = (r[4 - i] + p[4 - i]) / (q[4 - i] + amq);
390  } else if (mode == 1) {
391  r[4 - i] = (-p[1 - i] * q[1 - i] - p[2 - i] * q[2 - i] - p[3 - i] * q[3 - i] + p[4 - i] * q[4 - i]) / amq;
392  fac = -(r[4 - i] + p[4 - i]) / (q[4 - i] + amq);
393  } else {
394  cout << " ++++++++ wrong mode in boostdq " << endl;
395  exit(-1);
396  }
397  r[1 - i] = p[1 - i] + fac * q[1 - i];
398  r[2 - i] = p[2 - i] + fac * q[2 - i];
399  r[3 - i] = p[3 - i] + fac * q[3 - i];
400  }
401 
402 
403 //----------------------------------------------------------------------
404 //
405 // PHOTOS: PHOton radiation in decays BOost routine '3'
406 //
407 // Purpose: Boost vector PVEC along z-axis where ANGLE = EXP(ETA),
408 // ETA is the hyperbolic velocity.
409 //
410 // Input Parameters: ANGLE, PVEC
411 //
412 // Output Parameter: PVEC
413 //
414 // Author(s): S. Jadach Created at: 01/01/89
415 // B. van Eijk Last Update: 12/06/13
416 //
417 //----------------------------------------------------------------------
418  void PHOBO3(double ANGLE, double PVEC[4])
419  {
420  int j = 1; // convention of indices of Riemann space must be preserved.
421  double QPL, QMI;
422  QPL = (PVEC[4 - j] + PVEC[3 - j]) * ANGLE;
423  QMI = (PVEC[4 - j] - PVEC[3 - j]) / ANGLE;
424  PVEC[3 - j] = (QPL - QMI) / 2.0;
425  PVEC[4 - j] = (QPL + QMI) / 2.0;
426  }
427 
428  } // namespace PhotosUtilities
429 
430 } // namespace Photospp
431 
PhotosUtilities
Support functions.