Belle II Software  release-06-00-14
photosC.cc
1 #include "Photos.h"
2 #include "forW-MEc.h"
3 #include "forZ-MEc.h"
4 #include "pairs.h"
5 #include "Log.h"
6 #include <cstdio>
7 #include <cmath>
8 #include <iostream>
9 #include "photosC.h"
10 #include "HEPEVT_struct.h"
11 #include "PhotosUtilities.h"
12 using std::cout;
13 using std::endl;
14 using std::max;
15 using namespace Photospp;
16 using namespace PhotosUtilities;
17 
18 namespace Photospp {
19 
20 // Instantiating structs declared in photosC.h
21 
22  struct HEPEVT hep;
23  struct HEPEVT pho;
24 
25  struct PHOCOP phocop;
26  struct PHNUM phnum;
27  struct PHOKEY phokey;
28  struct PHOSTA phosta;
29  struct PHOLUN pholun;
30  struct PHOPHS phophs;
31  struct TOFROM tofrom;
32  struct PHOPRO phopro;
33  struct PHOREST phorest;
34  struct PHWT phwt;
35  struct PHOCORWT phocorwt;
36  struct PHOMOM phomom;
37  struct PHOCMS phocms;
38  struct PHOEXP phoexp;
39  struct PHLUPY phlupy;
40  struct DARKR darkr;
45  bool F(int m, int i)
46  {
47  return Photos::IPHQRK_setQarknoEmission(0, i) && (i <= 41 || i > 100)
48  && i != 21
49  && i != 2101 && i != 3101 && i != 3201
50  && i != 1103 && i != 2103 && i != 2203
51  && i != 3103 && i != 3203 && i != 3303;
52  }
53 
54 
55 // --- can be used with VARIANT A. For B use PHINT1 or 2 --------------
56 //----------------------------------------------------------------------
57 //
58 // PHINT: PHotos universal INTerference correction weight
59 //
60 // Purpose: calculates correction weight as expressed by
61 // formula (17) from CPC 79 (1994), 291.
62 //
63 // Input Parameters: Common /PHOEVT/, with photon added.
64 //
65 // Output Parameters: correction weight
66 //
67 // Author(s): Z. Was, P.Golonka Created at: 19/01/05
68 // Last Update: 23/06/13
69 //
70 //----------------------------------------------------------------------
71 
72  double PHINT(int IDUM)
73  {
74 
75  double PHINT2;
76  double EPS1[4], EPS2[4], PH[4], PL[4];
77  static int i = 1;
78  int K, L;
79  // DOUBLE PRECISION EMU,MCHREN,BETA,COSTHG,MPASQR,XPH, XC1, XC2
80  double XNUM1, XNUM2, XDENO, XC1, XC2;
81 
82  // REAL*8 PHOCHA
83  //--
84 
85  // Calculate polarimetric vector: ph, eps1, eps2 are orthogonal
86 
87  for (K = 1; K <= 4; K++) {
88  PH[K - i] = pho.phep[pho.nhep - i][K - i];
89  EPS2[K - i] = 1.0;
90  }
91 
92 
93  PHOEPS(PH, EPS2, EPS1);
94  PHOEPS(PH, EPS1, EPS2);
95 
96 
97  XNUM1 = 0.0;
98  XNUM2 = 0.0;
99  XDENO = 0.0;
100 
101  for (K = pho.jdahep[1 - i][1 - i]; K <= pho.nhep - i; K++) {
102 
103  // momenta of charged particle in PL
104 
105  for (L = 1; L <= 4; L++) PL[L - i] = pho.phep[K - i][L - i];
106 
107  // scalar products: epsilon*p/k*p
108 
109  XC1 = - PHOCHA(pho.idhep[K - i]) *
110  (PL[1 - i] * EPS1[1 - i] + PL[2 - i] * EPS1[2 - i] + PL[3 - i] * EPS1[3 - i]) /
111  (PH[4 - i] * PL[4 - i] - PH[1 - i] * PL[1 - i] - PH[2 - i] * PL[2 - i] - PH[3 - i] * PL[3 - i]);
112 
113  XC2 = - PHOCHA(pho.idhep[K - i]) *
114  (PL[1 - i] * EPS2[1 - i] + PL[2 - i] * EPS2[2 - i] + PL[3 - i] * EPS2[3 - i]) /
115  (PH[4 - i] * PL[4 - i] - PH[1 - i] * PL[1 - i] - PH[2 - i] * PL[2 - i] - PH[3 - i] * PL[3 - i]);
116 
117 
118  // accumulate the currents
119  XNUM1 = XNUM1 + XC1;
120  XNUM2 = XNUM2 + XC2;
121 
122  XDENO = XDENO + XC1 * XC1 + XC2 * XC2;
123  }
124 
125  PHINT2 = (XNUM1 * XNUM1 + XNUM2 * XNUM2) / XDENO;
126  return PHINT2;
127 
128  }
129 
130 
131 
132 //----------------------------------------------------------------------
133 //
134 // PHINT: PHotos INTerference (Old version kept for tests only.
135 //
136 // Purpose: Calculates interference between emission of photons from
137 // different possible chaged daughters stored in
138 // the HEP common /PHOEVT/.
139 //
140 // Input Parameter: commons /PHOEVT/ /PHOMOM/ /PHOPHS/
141 //
142 //
143 // Output Parameters:
144 //
145 //
146 // Author(s): Z. Was, Created at: 10/08/93
147 // Last Update: 15/03/99
148 //
149 //----------------------------------------------------------------------
150 
151  double PHINT1(int IDUM)
152  {
153 
154  double PHINT;
155 
156  /*
157  DOUBLE PRECISION MCHSQR,MNESQR
158  REAL*8 PNEUTR
159  COMMON/PHOMOM/MCHSQR,MNESQR,PNEUTR(5)
160  DOUBLE PRECISION COSTHG,SINTHG
161  REAL*8 XPHMAX,XPHOTO
162  COMMON/PHOPHS/XPHMAX,XPHOTO,COSTHG,SINTHG
163 
164  */
165  double MPASQR, XX, BETA;
166  bool IFINT;
167  int K, IDENT;
168  double& COSTHG = phophs.costhg;
169  double& XPHOTO = phophs.xphoto;
170  //double *PNEUTR = phomom.pneutr; // unused variable
171  double& MCHSQR = phomom.mchsqr;
172  double& MNESQR = phomom.mnesqr;
173 
174  static int i = 1;
175  IDENT = pho.nhep;
176  //
177  for (K = pho.jdahep[1 - i][2 - i]; K >= pho.jdahep[1 - i][1 - i]; K--) {
178  if (pho.idhep[K - i] != 22) {
179  IDENT = K;
180  break;
181  }
182  }
183 
184  // check if there is a photon
185  IFINT = pho.nhep > IDENT;
186  // check if it is two body + gammas reaction
187  IFINT = IFINT && (IDENT - pho.jdahep[1 - i][1 - i]) == 1;
188  // check if two body was particle antiparticle
189  IFINT = IFINT && pho.idhep[pho.jdahep[1 - i][1 - i] - i] == -pho.idhep[IDENT - i];
190  // check if particles were charged
191  IFINT = IFINT && PHOCHA(pho.idhep[IDENT - i]) != 0;
192  // calculates interference weight contribution
193  if (IFINT) {
194  MPASQR = pho.phep[1 - i][5 - i] * pho.phep[1 - i][5 - i];
195  XX = 4.0 * MCHSQR / MPASQR * (1.0 - XPHOTO) / (1.0 - XPHOTO + (MCHSQR - MNESQR) / MPASQR) / (1.0 - XPHOTO +
196  (MCHSQR - MNESQR) / MPASQR);
197  BETA = sqrt(1.0 - XX);
198  PHINT = 2.0 / (1.0 + COSTHG * COSTHG * BETA * BETA);
199  } else {
200  PHINT = 1.0;
201  }
202 
203  return PHINT;
204  }
205 
206 
207 //----------------------------------------------------------------------
208 //
209 // PHINT: PHotos INTerference
210 //
211 // Purpose: Calculates interference between emission of photons from
212 // different possible chaged daughters stored in
213 // the HEP common /PHOEVT/.
214 //
215 // Input Parameter: commons /PHOEVT/ /PHOMOM/ /PHOPHS/
216 //
217 //
218 // Output Parameters:
219 //
220 //
221 // Author(s): Z. Was, Created at: 10/08/93
222 // Last Update:
223 //
224 //----------------------------------------------------------------------
225 
226  double PHINT2(int IDUM)
227  {
228 
229 
230  /*
231  DOUBLE PRECISION MCHSQR,MNESQR
232  REAL*8 PNEUTR
233  COMMON/PHOMOM/MCHSQR,MNESQR,PNEUTR(5)
234  DOUBLE PRECISION COSTHG,SINTHG
235  REAL*8 XPHMAX,XPHOTO
236  COMMON/PHOPHS/XPHMAX,XPHOTO,COSTHG,SINTHG
237  */
238  double& COSTHG = phophs.costhg;
239  double& XPHOTO = phophs.xphoto;
240  double& MCHSQR = phomom.mchsqr;
241  double& MNESQR = phomom.mnesqr;
242 
243 
244  double MPASQR, XX, BETA, pq1[4], pq2[4], pphot[4];
245  double SS, PP2, PP, E1, E2, q1, q2, costhe, PHINT;
246  bool IFINT;
247  int K, k, IDENT;
248  static int i = 1;
249  IDENT = pho.nhep;
250  //
251  for (K = pho.jdahep[1 - i][2 - i]; K >= pho.jdahep[1 - i][1 - i]; K--) {
252  if (pho.idhep[K - i] != 22) {
253  IDENT = K;
254  break;
255  }
256  }
257 
258  // check if there is a photon
259  IFINT = pho.nhep > IDENT;
260  // check if it is two body + gammas reaction
261  IFINT = IFINT && (IDENT - pho.jdahep[1 - i][1 - i]) == 1;
262  // check if two body was particle antiparticle (we improve on it !
263  // IFINT= IFINT.AND.pho.idhep(JDAPHO(1,1)).EQ.-pho.idhep(IDENT)
264  // check if particles were charged
265  IFINT = IFINT && fabs(PHOCHA(pho.idhep[IDENT - i])) > 0.01;
266  // check if they have both charge
267  IFINT = IFINT && fabs(PHOCHA(pho.idhep[pho.jdahep[1 - i][1 - i] - i])) > 0.01;
268  // calculates interference weight contribution
269  if (IFINT) {
270  MPASQR = pho.phep[1 - i][5 - i] * pho.phep[1 - i][5 - i];
271  XX = 4.0 * MCHSQR / MPASQR * (1.0 - XPHOTO) / pow(1. - XPHOTO + (MCHSQR - MNESQR) / MPASQR, 2);
272  BETA = sqrt(1.0 - XX);
273  PHINT = 2.0 / (1.0 + COSTHG * COSTHG * BETA * BETA);
274  SS = MPASQR * (1.0 - XPHOTO);
275  PP2 = ((SS - MCHSQR - MNESQR) * (SS - MCHSQR - MNESQR) - 4 * MCHSQR * MNESQR) / SS / 4;
276  PP = sqrt(PP2);
277  E1 = sqrt(PP2 + MCHSQR);
278  E2 = sqrt(PP2 + MNESQR);
279  PHINT = (E1 + E2) * (E1 + E2) / ((E2 + COSTHG * PP) * (E2 + COSTHG * PP) + (E1 - COSTHG * PP) * (E1 - COSTHG * PP));
280  // return PHINT;
281  //
282  q1 = PHOCHA(pho.idhep[pho.jdahep[1 - i][1 - i] - i]);
283  q2 = PHOCHA(pho.idhep[IDENT - i]);
284  for (k = 1; k <= 4; k++) {
285  pq1[k - i] = pho.phep[pho.jdahep[1 - i][1 - i] - i][k - i];
286  pq2[k - i] = pho.phep[pho.jdahep[1 - i][1 - i] + 1 - i][k - i];
287  pphot[k - i] = pho.phep[pho.nhep - i][k - i];
288  }
289  costhe = (pphot[1 - i] * pq1[1 - i] + pphot[2 - i] * pq1[2 - i] + pphot[3 - i] * pq1[3 - i]);
290  costhe = costhe / sqrt(pq1[1 - i] * pq1[1 - i] + pq1[2 - i] * pq1[2 - i] + pq1[3 - i] * pq1[3 - i]);
291  costhe = costhe / sqrt(pphot[1 - i] * pphot[1 - i] + pphot[2 - i] * pphot[2 - i] + pphot[3 - i] * pphot[3 - i]);
292  //
293  // --- this IF checks whether JDAPHO(1,1) was MCH or MNE.
294  // --- COSTHG angle (and in-generation variables) may be better choice
295  // --- than costhe. note that in the formulae below amplitudes were
296  // --- multiplied by (E2+COSTHG*PP)*(E1-COSTHG*PP).
297  if (COSTHG * costhe > 0) {
298 
299  PHINT = pow(q1 * (E2 + COSTHG * PP) - q2 * (E1 - COSTHG * PP),
300  2) / (q1 * q1 * (E2 + COSTHG * PP) * (E2 + COSTHG * PP) + q2 * q2 * (E1 - COSTHG * PP) * (E1 - COSTHG * PP));
301  } else {
302 
303  PHINT = pow(q1 * (E1 - COSTHG * PP) - q2 * (E2 + COSTHG * PP),
304  2) / (q1 * q1 * (E1 - COSTHG * PP) * (E1 - COSTHG * PP) + q2 * q2 * (E2 + COSTHG * PP) * (E2 + COSTHG * PP));
305  }
306  } else {
307  PHINT = 1.0;
308  }
309  return PHINT;
310  }
311 
312 
313 //*****************************************************************
314 //*****************************************************************
315 //*****************************************************************
316 // beginning of the class of methods reading from PH_HEPEVT
317 //*****************************************************************
318 //*****************************************************************
319 //*****************************************************************
320 
321 
322 //----------------------------------------------------------------------
323 //
324 // PHOTOS: PHOton radiation in decays event DuMP routine
325 //
326 // Purpose: Print event record.
327 //
328 // Input Parameters: Common /PH_HEPEVT/
329 //
330 // Output Parameters: None
331 //
332 // Author(s): B. van Eijk Created at: 05/06/90
333 // Last Update: 20/06/13
334 //
335 //----------------------------------------------------------------------
336  void PHODMP()
337  {
338 
339  double SUMVEC[5];
340  int I, J;
341  static int i = 1;
342  const char eq80[81] = "================================================================================";
343  const char X29[30] = " ";
344  const char X23[24 ] = " ";
345  const char X1[2] = " ";
346  const char X2[3] = " ";
347  const char X3[4] = " ";
348  const char X4[5] = " ";
349  const char X6[7] = " ";
350  const char X7[8] = " ";
351  FILE* PHLUN = stdout;
352 
353  for (I = 0; I < 5; I++) SUMVEC[I] = 0.0;
354  //--
355  //-- Print event number...
356  fprintf(PHLUN, "%s", eq80);
357  fprintf(PHLUN, "%s Event No.: %10i\n", X29, hep.nevhep);
358  fprintf(PHLUN, "%s Particle Parameters\n", X6);
359  fprintf(PHLUN, "%s Nr %s Type %s Parent(s) %s Daughter(s) %s Px %s Py %s Pz %s E %s Inv. M.\n", X1, X3, X3, X2, X6, X7, X7, X7, X4);
360  for (I = 1; I <= hep.nhep; I++) {
361  //--
362  //-- For 'stable particle' calculate vector momentum sum
363  if (hep.jdahep[I - i][1 - i] == 0) {
364  for (J = 1; J <= 4; J++) {
365  SUMVEC[J - i] = SUMVEC[J - i] + hep.phep[I - i][J - i];
366  }
367  if (hep.jmohep[I - i][2 - i] == 0) {
368  fprintf(PHLUN, "%4i %7i %s %4i %s Stable %9.2f %9.2f %9.2f %9.2f %9.2f\n" , I, hep.idhep[I - i], X3, hep.jmohep[I - i][1 - i], X7,
369  hep.phep[I - i][1 - i], hep.phep[I - i][2 - i], hep.phep[I - i][3 - i], hep.phep[I - i][4 - i], hep.phep[I - i][5 - i]);
370  } else {
371  fprintf(PHLUN, "%4i %7i %4i - %4i %s Stable %9.2f %9.2f %9.2f %9.2f %9.2f\n", I, hep.idhep[I - i], hep.jmohep[I - i][1 - i],
372  hep.jmohep[I - i][2 - i], X4, hep.phep[I - i][1 - i], hep.phep[I - i][2 - i], hep.phep[I - i][3 - i], hep.phep[I - i][4 - i],
373  hep.phep[I - i][5 - i]);
374  }
375  } else {
376  if (hep.jmohep[I - i][2 - i] == 0) {
377  fprintf(PHLUN, "%4i %7i %s %4i %s %4i - %4i %9.2f %9.2f %9.2f %9.2f %9.2f\n" , I, hep.idhep[I - i], X3, hep.jmohep[I - i][1 - i],
378  X2, hep.jdahep[I - i][1 - i], hep.jdahep[I - i][2 - i], hep.phep[I - i][1 - i], hep.phep[I - i][2 - i], hep.phep[I - i][3 - i],
379  hep.phep[I - i][4 - i], hep.phep[I - i][5 - i]);
380  } else {
381  fprintf(PHLUN, "%4i %7i %4i - %4i %4i - %4i %9.2f %9.2f %9.2f %9.2f %9.2f\n", I, hep.idhep[I - i], hep.jmohep[I - i][1 - i],
382  hep.jmohep[I - i][2 - i], hep.jdahep[I - i][1 - i], hep.jdahep[I - i][2 - i], hep.phep[I - i][1 - i], hep.phep[I - i][2 - i],
383  hep.phep[I - i][3 - i], hep.phep[I - i][4 - i], hep.phep[I - i][5 - i]);
384  }
385  }
386  }
387  SUMVEC[5 - i] = sqrt(SUMVEC[4 - i] * SUMVEC[4 - i] - SUMVEC[1 - i] * SUMVEC[1 - i] - SUMVEC[2 - i] * SUMVEC[2 - i] - SUMVEC[3 - i] *
388  SUMVEC[3 - i]);
389  fprintf(PHLUN, "%s Vector Sum: %9.2f %9.2f %9.2f %9.2f %9.2f\n", X23, SUMVEC[1 - i], SUMVEC[2 - i], SUMVEC[3 - i], SUMVEC[4 - i],
390  SUMVEC[5 - i]);
391 
392 
393 
394 
395 // 9030 FORMAT(1H ,I4,I7,3X,I4,9X,'Stable',2X,5F9.2)
396 //"%4i %7i %s %4i %s Stable %s %9.2f %9.2f %9.2f %9.2f %9.2f " X3,9X,X2
397 
398  // 9050 FORMAT(1H ,I4,I7,3X,I4,6X,I4,' - ',I4,5F9.2)
399  //"%4i %7i %s %4i %s %4i - %4i %9.2f %9.2f %9.2f %9.2f %9.2f " X3,X6
400 
401 
402 
403 
404  //"%4i %7i %4i - %4i %s Stable %s %9.2f %9.2f %9.2f %9.2f %9.2f " X5,X2
405 
406 
407 //9060 FORMAT(1H ,I4,I7,I4,' - ',I4,2X,I4,' - ',I4,5F9.2)
408  //"%4i %7i %4i - %4i %s %4i - %4i %9.2f %9.2f %9.2f %9.2f %9.2f " X2,
409  }
410 
411 
412 
413 //----------------------------------------------------------------------
414 //
415 // PHLUPAB: debugging tool
416 //
417 // Purpose: NONE, eventually may printout content of the
418 // /PH_HEPEVT/ common
419 //
420 // Input Parameters: Common /PH_HEPEVT/ and /PHNUM/
421 // latter may have number of the event.
422 //
423 // Output Parameters: None
424 //
425 // Author(s): Z. Was Created at: 30/05/93
426 // Last Update: 20/06/13
427 //
428 //----------------------------------------------------------------------
429 
430  void PHLUPAB(int IPOINT)
431  {
432  char name[12] = "/PH_HEPEVT/";
433  int I, J;
434  static int IPOIN0 = -5;
435  static int i = 1;
436  double SUM[5];
437  FILE* PHLUN = stdout;
438 
439  if (IPOIN0 < 0) {
440  IPOIN0 = 400000; // ! maximal no-print point
441  phlupy.ipoin = IPOIN0;
442  phlupy.ipoinm = 400001; // ! minimal no-print point
443  }
444 
445  if (IPOINT <= phlupy.ipoinm || IPOINT >= phlupy.ipoin) return;
446  if ((int)phnum.iev < 1000) {
447  for (I = 1; I <= 5; I++) SUM[I - i] = 0.0;
448 
449  fprintf(PHLUN, "EVENT NR= %i WE ARE TESTING %s at IPOINT=%i \n", (int)phnum.iev, name, IPOINT);
450  fprintf(PHLUN, " ID p_x p_y p_z E m ID-MO_DA1 ID-MO_DA2\n");
451  I = 1;
452  fprintf(PHLUN, "%4i %14.9f %14.9f %14.9f %14.9f %14.9f %9i %9i\n", hep.idhep[I - i], hep.phep[I - i][1 - i], hep.phep[I - i][2 - i],
453  hep.phep[I - i][3 - i], hep.phep[I - i][4 - i], hep.phep[I - i][5 - i], hep.jdahep[I - i][1 - i], hep.jdahep[I - i][2 - i]);
454  I = 2;
455  fprintf(PHLUN, "%4i %14.9f %14.9f %14.9f %14.9f %14.9f %9i %9i\n", hep.idhep[I - i], hep.phep[I - i][1 - i], hep.phep[I - i][2 - i],
456  hep.phep[I - i][3 - i], hep.phep[I - i][4 - i], hep.phep[I - i][5 - i], hep.jdahep[I - i][1 - i], hep.jdahep[I - i][2 - i]);
457  fprintf(PHLUN, " \n");
458  for (I = 3; I <= hep.nhep; I++) {
459  fprintf(PHLUN, "%4i %14.9f %14.9f %14.9f %14.9f %14.9f %9i %9i\n", hep.idhep[I - i], hep.phep[I - i][1 - i], hep.phep[I - i][2 - i],
460  hep.phep[I - i][3 - i], hep.phep[I - i][4 - i], hep.phep[I - i][5 - i], hep.jmohep[I - i][1 - i], hep.jmohep[I - i][2 - i]);
461  for (J = 1; J <= 4; J++) SUM[J - i] = SUM[J - i] + hep.phep[I - i][J - i];
462  }
463 
464 
465  SUM[5 - i] = sqrt(fabs(SUM[4 - i] * SUM[4 - i] - SUM[1 - i] * SUM[1 - i] - SUM[2 - i] * SUM[2 - i] - SUM[3 - i] * SUM[3 - i]));
466  fprintf(PHLUN, " SUM %14.9f %14.9f %14.9f %14.9f %14.9f\n", SUM[1 - i], SUM[2 - i], SUM[3 - i], SUM[4 - i], SUM[5 - i]);
467 
468  }
469 
470 
471  // 10 FORMAT(1X,' ID ','p_x ','p_y ','p_z ',
472  //$ 'E ','m ',
473  //$ 'ID-MO_DA1','ID-MO DA2' )
474  // 20 FORMAT(1X,I4,5(F14.9),2I9)
475  //"%i4 %14.9f %14.9f %14.9f %14.9f %i9 i9"
476  // 30 FORMAT(1X,' SUM',5(F14.9))
477  }
478 
479 
480 
481 
482 
483 
484 
485 
486 
487 //----------------------------------------------------------------------
488 //
489 // PHLUPA: debugging tool
490 //
491 // Purpose: NONE, eventually may printout content of the
492 // /PHOEVT/ common
493 //
494 // Input Parameters: Common /PHOEVT/ and /PHNUM/
495 // latter may have number of the event.
496 //
497 // Output Parameters: None
498 //
499 // Author(s): Z. Was Created at: 30/05/93
500 // Last Update: 21/06/13
501 //
502 //----------------------------------------------------------------------
503 
504  void PHLUPA(int IPOINT)
505  {
506  char name[9] = "/PHOEVT/";
507  int I, J;
508  static int IPOIN0 = -5;
509  static int i = 1;
510  double SUM[5];
511  FILE* PHLUN = stdout;
512 
513  if (IPOIN0 < 0) {
514  IPOIN0 = 400000; // ! maximal no-print point
515  phlupy.ipoin = IPOIN0;
516  phlupy.ipoinm = 400001; // ! minimal no-print point
517  }
518 
519  if (IPOINT <= phlupy.ipoinm || IPOINT >= phlupy.ipoin) return;
520  if ((int)phnum.iev < 1000) {
521  for (I = 1; I <= 5; I++) SUM[I - i] = 0.0;
522 
523  fprintf(PHLUN, "EVENT NR= %i WE ARE TESTING %s at IPOINT=%i \n", (int)phnum.iev, name, IPOINT);
524  fprintf(PHLUN, " ID p_x p_y p_z E m ID-MO_DA1 ID-MO_DA2\n");
525  I = 1;
526  fprintf(PHLUN, "%4i %14.9f %14.9f %14.9f %14.9f %14.9f %9i %9i\n", pho.idhep[I - i], pho.phep[I - i][1 - i], pho.phep[I - i][2 - i],
527  pho.phep[I - i][3 - i], pho.phep[I - i][4 - i], pho.phep[I - i][5 - i], pho.jdahep[I - i][1 - i], pho.jdahep[I - i][2 - i]);
528  I = 2;
529  fprintf(PHLUN, "%4i %14.9f %14.9f %14.9f %14.9f %14.9f %9i %9i\n", pho.idhep[I - i], pho.phep[I - i][1 - i], pho.phep[I - i][2 - i],
530  pho.phep[I - i][3 - i], pho.phep[I - i][4 - i], pho.phep[I - i][5 - i], pho.jdahep[I - i][1 - i], pho.jdahep[I - i][2 - i]);
531  fprintf(PHLUN, " \n");
532  for (I = 3; I <= pho.nhep; I++) {
533  fprintf(PHLUN, "%4i %14.9f %14.9f %14.9f %14.9f %14.9f %9i %9i\n", pho.idhep[I - i], pho.phep[I - i][1 - i], pho.phep[I - i][2 - i],
534  pho.phep[I - i][3 - i], pho.phep[I - i][4 - i], pho.phep[I - i][5 - i], pho.jmohep[I - i][1 - i], pho.jmohep[I - i][2 - i]);
535  for (J = 1; J <= 4; J++) SUM[J - i] = SUM[J - i] + pho.phep[I - i][J - i];
536  }
537 
538 
539  SUM[5 - i] = sqrt(fabs(SUM[4 - i] * SUM[4 - i] - SUM[1 - i] * SUM[1 - i] - SUM[2 - i] * SUM[2 - i] - SUM[3 - i] * SUM[3 - i]));
540  fprintf(PHLUN, " SUM %14.9f %14.9f %14.9f %14.9f %14.9f\n", SUM[1 - i], SUM[2 - i], SUM[3 - i], SUM[4 - i], SUM[5 - i]);
541 
542  }
543 
544 
545  // 10 FORMAT(1X,' ID ','p_x ','p_y ','p_z ',
546  //$ 'E ','m ',
547  //$ 'ID-MO_DA1','ID-MO DA2' )
548  // 20 FORMAT(1X,I4,5(F14.9),2I9)
549  //"%4i %14.9f %14.9f %14.9f %14.9f %9i %9i"
550  // 30 FORMAT(1X,' SUM',5(F14.9))
551  }
552 
553 
554  void PHOtoRF()
555  {
556 
557 
558  // COMMON /PH_TOFROM/ QQ[4],XM,th1,fi1
559  double PP[4], RR[4];
560 
561  int K, L;
562  static int i = 1;
563 
564  for (K = 1; K <= 4; K++) {
565  tofrom.QQ[K - i] = 0.0;
566  }
567  for (L = hep.jdahep[hep.jmohep[hep.nhep - i][1 - i] - i][1 - i]; L <= hep.jdahep[hep.jmohep[hep.nhep - i][1 - i] - i][2 - i]; L++) {
568  for (K = 1; K <= 4; K++) {
569  tofrom.QQ[K - i] = tofrom.QQ[K - i] + hep.phep[L - i][K - i];
570  }
571  }
572  tofrom.XM = tofrom.QQ[4 - i] * tofrom.QQ[4 - i] - tofrom.QQ[3 - i] * tofrom.QQ[3 - i] - tofrom.QQ[2 - i] * tofrom.QQ[2 - i] -
573  tofrom.QQ[1 - i] * tofrom.QQ[1 - i];
574  if (tofrom.XM > 0.0) tofrom.XM = sqrt(tofrom.XM);
575  if (tofrom.XM <= 0.0) return;
576 
577  for (L = 1; L <= hep.nhep; L++) {
578  for (K = 1; K <= 4; K++) {
579  PP[K - i] = hep.phep[L - i][K - i];
580  }
581  bostdq(1, tofrom.QQ, PP, RR);
582  for (K = 1; K <= 4; K++) {
583  hep.phep[L - i][K - i] = RR[K - i];
584  }
585  }
586 
587  tofrom.fi1 = 0.0;
588  tofrom.th1 = 0.0;
589  if (fabs(hep.phep[1 - i][1 - i]) + fabs(hep.phep[1 - i][2 - i]) > 0.0) tofrom.fi1 = PHOAN1(hep.phep[1 - i][1 - i],
590  hep.phep[1 - i][2 - i]);
591  if (fabs(hep.phep[1 - i][1 - i]) + fabs(hep.phep[1 - i][2 - i]) + fabs(hep.phep[1 - i][3 - i]) > 0.0)
592  tofrom.th1 = PHOAN2(hep.phep[1 - i][3 - i],
593  sqrt(hep.phep[1 - i][1 - i] * hep.phep[1 - i][1 - i] + hep.phep[1 - i][2 - i] * hep.phep[1 - i][2 - i]));
594 
595  for (L = 1; L <= hep.nhep; L++) {
596  for (K = 1; K <= 4; K++) {
597  RR[K - i] = hep.phep[L - i][K - i];
598  }
599 
600  PHORO3(-tofrom.fi1, RR);
601  PHORO2(-tofrom.th1, RR);
602  for (K = 1; K <= 4; K++) {
603  hep.phep[L - i][K - i] = RR[K - i];
604  }
605  }
606 
607  return;
608  }
609 
610  void PHOtoLAB()
611  {
612 
613  // // REAL*8 QQ(4),XM,th1,fi1
614  // COMMON /PH_TOFROM/ QQ,XM,th1,fi1
615  double PP[4], RR[4];
616  int K, L;
617  static int i = 1;
618 
619  if (tofrom.XM <= 0.0) return;
620 
621 
622  for (L = 1; L <= hep.nhep; L++) {
623  for (K = 1; K <= 4; K++) {
624  PP[K - i] = hep.phep[L - i][K - i];
625  }
626 
627  PHORO2(tofrom.th1, PP);
628  PHORO3(tofrom.fi1, PP);
629  bostdq(-1, tofrom.QQ, PP, RR);
630 
631  for (K = 1; K <= 4; K++) {
632  hep.phep[L - i][K - i] = RR[K - i];
633  }
634  }
635  return;
636  }
637  void PHOcorrectDARK(int IDaction)
638  {
639  int K, L;
640  static int i = 1;
641  // if(darkr.ifspecial==1)
642  {
643  for (L = 1; L <= hep.nhep; L++) {
644  if (hep.idhep[L - i] == darkr.IDspecial) {
645  //PHODMP();
646  hep.jdahep[L - i][1 - i] = L - i + 1 + 1;
647  hep.jdahep[L - i][2 - i] = L - i + 2 + 1;
648  hep.jmohep[L - i + 1][1 - i] = L - i + 1;
649  hep.jmohep[L - i + 2][1 - i] = L - i + 1;
650 
651  hep.jmohep[L - i - 2][2 - i] = 0;
652  hep.jmohep[L - i - 1][2 - i] = 0;
653  hep.jmohep[L - i][2 - i] = 0;
654 
655  hep.jmohep[L - i + 1][2 - i] = 0;
656  hep.jmohep[L - i + 2][2 - i] = 0;
657  // cout << "adres prosze= " << hep.jmohep[L-i][1-i] << endl;
658  hep.jdahep[ hep.jmohep[L - i][1 - i] - i ][2 - i] = hep.jdahep[ hep.jmohep[L - i][1 - i] - i ][2 - i] - 2;
659  //PHODMP();
660  }
661  }
662  }
663  }
664 
665 
666 
667 
668 
669 
670 // 2) GENERAL INTERFACE:
671 // PHOTOS_GET
672 // PHOTOS_MAKE
673 
674 
675 // COMMONS:
676 // NAME USED IN SECT. # OF OC// Comment
677 // PHOQED 1) 2) 3 Flags whether emisson to be gen.
678 // PHOLUN 1) 4) 6 Output device number
679 // PHOCOP 1) 3) 4 photon coupling & min energy
680 // PHPICO 1) 3) 4) 5 PI & 2*PI
681 // PHSEED 1) 4) 3 RN seed
682 // PHOSTA 1) 4) 3 Status information
683 // PHOKEY 1) 2) 3) 7 Keys for nonstandard application
684 // PHOVER 1) 1 Version info for outside
685 // HEPEVT 2) 2 PDG common
686 // PH_HEPEVT2) 8 PDG common internal
687 // PHOEVT 2) 3) 10 PDG branch
688 // PHOIF 2) 3) 2 emission flags for PDG branch
689 // PHOMOM 3) 5 param of char-neutr system
690 // PHOPHS 3) 5 photon momentum parameters
691 // PHOPRO 3) 4 var. for photon rep. (in branch)
692 // PHOCMS 2) 3 parameters of boost to branch CMS
693 // PHNUM 4) 1 event number from outside
694 //----------------------------------------------------------------------
695 
696 
697 //----------------------------------------------------------------------
698 //
699 // PHOTOS_MAKE: General search routine
700 //
701 // Purpose: Search through the /PH_HEPEVT/ standard HEP common, sta-
702 // rting from the IPPAR-th particle. Whenevr branching
703 // point is found routine PHTYPE(IP) is called.
704 // Finally if calls on PHTYPE(IP) modified entries, common
705 // /PH_HEPEVT/ is ordered.
706 //
707 // Input Parameter: IPPAR: Pointer to decaying particle in
708 // /PH_HEPEVT/ and the common itself,
709 //
710 // Output Parameters: Common /PH_HEPEVT/, either with or without
711 // new particles added.
712 //
713 // Author(s): Z. Was, B. van Eijk Created at: 26/11/89
714 // Last Update: 30/08/93
715 //
716 //----------------------------------------------------------------------
717 
718  void PHOTOS_MAKE_C(int IPARR)
719  {
720  static int i = 1;
721  int IPPAR, I, J, NLAST, MOTHER;
722  //--
723  PHLUPAB(3);
724 
725  // write(*,*) 'at poczatek'
726  // PHODMP();
727  IPPAR = abs(IPARR);
728  //-- Store pointers for cascade treatement...
729  NLAST = hep.nhep;
730 
731 
732  //--
733  //-- Check decay multiplicity and minimum of correctness..
734  if ((hep.jdahep[IPPAR - i][1 - i] == 0) || (hep.jmohep[hep.jdahep[IPPAR - i][1 - i] - i][1 - i] != IPPAR)) return;
735 
736  PHOtoRF();
737 
738  // write(*,*) 'at przygotowany'
739  // PHODMP();
740 
741  //--
742  //-- single branch mode
743  //-- IPPAR is original position where the program was called
744 
745  //-- let-s do generation
746  PHTYPE(IPPAR);
747 
748 
749  //-- rearrange /PH_HEPEVT/ for added particles.
750  //-- at present this may be not needed as information
751  //-- is set at HepMC level.
752  if (hep.nhep > NLAST) {
753  for (I = NLAST + 1; I <= hep.nhep; I++) {
754  //--
755  //-- Photon mother and vertex...
756  MOTHER = hep.jmohep[I - i][1 - i];
757  hep.jdahep[MOTHER - i][2 - i] = I;
758  for (J = 1; J <= 4; J++) {
759  hep.vhep[I - i][J - i] = hep.vhep[I - 1 - i][J - i];
760  }
761  }
762  }
763  // write(*,*) 'at po dzialaniu '
764  // PHODMP();
765  PHOtoLAB();
766  PHOcorrectDARK(1);
767  // write(*,*) 'at koniec'
768  // PHODMP();
769  return;
770  }
771 
772 
773 
774 
775 //----------------------------------------------------------------------
776 //
777 // PHCORK: corrects kinmatics of subbranch needed if host program
778 // produces events with the shaky momentum conservation
779 //
780 // Input Parameters: Common /PHOEVT/, MODCOR
781 // MODCOR >0 type of action
782 // =1 no action
783 // =2 corrects energy from mass
784 // =3 corrects mass from energy
785 // =4 corrects energy from mass for
786 // particles up to .4 GeV mass,
787 // for heavier ones corrects mass,
788 // =5 most complete correct also of mother
789 // often necessary for exponentiation.
790 // =0 execution mode
791 //
792 // Output Parameters: corrected /PHOEVT/
793 //
794 // Author(s): P.Golonka, Z. Was Created at: 01/02/99
795 // Modified : 07/07/13
796 //----------------------------------------------------------------------
797 
798  void PHCORK(int MODCOR)
799  {
800 
801  double M, P2, PX, PY, PZ, E, EN, XMS;
802  int I, K;
803  FILE* PHLUN = stdout;
804 
805 
806  static int MODOP = 0;
807  static int IPRINT = 0;
808  static double MCUT = 0.4;
809  static int i = 1;
810 
811  if (MODCOR != 0) {
812  // INITIALIZATION
813  MODOP = MODCOR;
814 
815  fprintf(PHLUN, "Message from PHCORK(MODCOR):: initialization\n");
816  if (MODOP == 1) fprintf(PHLUN, "MODOP=1 -- no corrections on event: DEFAULT\n");
817  else if (MODOP == 2) fprintf(PHLUN, "MODOP=2 -- corrects Energy from mass\n");
818  else if (MODOP == 3) fprintf(PHLUN, "MODOP=3 -- corrects mass from Energy\n");
819  else if (MODOP == 4) {
820  fprintf(PHLUN, "MODOP=4 -- corrects Energy from mass to Mcut\n");
821  fprintf(PHLUN, " and mass from energy above Mcut\n");
822  fprintf(PHLUN, " Mcut=%6.3f GeV", MCUT);
823  } else if (MODOP == 5) fprintf(PHLUN, "MODOP=5 -- corrects Energy from mass+flow\n");
824 
825  else {
826  fprintf(PHLUN, "PHCORK wrong MODCOR=%4i\n", MODCOR);
827  exit(-1);
828  }
829  return;
830  }
831 
832  if (MODOP == 0 && MODCOR == 0) {
833  fprintf(PHLUN, "PHCORK lack of initialization\n");
834  exit(-1);
835  }
836 
837  // execution mode
838  // ==============
839  // ==============
840 
841 
842  PX = 0.0;
843  PY = 0.0;
844  PZ = 0.0;
845  E = 0.0;
846 
847  if (MODOP == 1) {
848  // -----------------------
849  // In this case we do nothing
850  return;
851  } else if (MODOP == 2) {
852  // -----------------------
853  // lets loop thru all daughters and correct their energies
854  // according to E^2=p^2+m^2
855 
856  for (I = 3; I <= pho.nhep; I++) {
857 
858  PX = PX + pho.phep[I - i][1 - i];
859  PY = PY + pho.phep[I - i][2 - i];
860  PZ = PZ + pho.phep[I - i][3 - i];
861 
862  P2 = pho.phep[I - i][1 - i] * pho.phep[I - i][1 - i] + pho.phep[I - i][2 - i] * pho.phep[I - i][2 - i] + pho.phep[I - i][3 - i] *
863  pho.phep[I - i][3 - i];
864 
865  EN = sqrt(pho.phep[I - i][5 - i] * pho.phep[I - i][5 - i] + P2);
866 
867  if (IPRINT == 1)fprintf(PHLUN, "CORRECTING ENERGY OF %6i: %14.9f => %14.9f\n", I, pho.phep[I - i][4 - i], EN);
868 
869  pho.phep[I - i][4 - i] = EN;
870  E = E + pho.phep[I - i][4 - i];
871 
872  }
873  }
874 
875  else if (MODOP == 5) {
876  // -----------------------
877  //C lets loop thru all daughters and correct their energies
878  //C according to E^2=p^2+m^2
879 
880  for (I = 3; I <= pho.nhep; I++) {
881  PX = PX + pho.phep[I - i][1 - i];
882  PY = PY + pho.phep[I - i][2 - i];
883  PZ = PZ + pho.phep[I - i][3 - i];
884 
885  P2 = pho.phep[I - i][1 - i] * pho.phep[I - i][1 - i] + pho.phep[I - i][2 - i] * pho.phep[I - i][2 - i] + pho.phep[I - i][3 - i] *
886  pho.phep[I - i][3 - i];
887 
888  EN = sqrt(pho.phep[I - i][5 - i] * pho.phep[I - i][5 - i] + P2);
889 
890  if (IPRINT == 1)fprintf(PHLUN, "CORRECTING ENERGY OF %6i: %14.9f => %14.9f\n", I, pho.phep[I - i][4 - i], EN);
891 
892  pho.phep[I - i][4 - i] = EN;
893  E = E + pho.phep[I - i][4 - i];
894 
895  }
896  for (K = 1; K <= 4; K++) {
897  pho.phep[1 - i][K - i] = 0.0;
898  for (I = 3; I <= pho.nhep; I++) {
899  pho.phep[1 - i][K - i] = pho.phep[1 - i][K - i] + pho.phep[I - i][K - i];
900  }
901  }
902  XMS = sqrt(pho.phep[1 - i][4 - i] * pho.phep[1 - i][4 - i] - pho.phep[1 - i][3 - i] * pho.phep[1 - i][3 - i] - pho.phep[1 - i][2 -
903  i] * pho.phep[1 - i][2 - i] - pho.phep[1 - i][1 - i] * pho.phep[1 - i][1 - i]);
904  pho.phep[1 - i][5 - i] = XMS;
905  } else if (MODOP == 3) {
906  // -----------------------
907 
908  // lets loop thru all daughters and correct their masses
909  // according to E^2=p^2+m^2
910 
911  for (I = 3; I <= pho.nhep; I++) {
912 
913  PX = PX + pho.phep[I - i][1 - i];
914  PY = PY + pho.phep[I - i][2 - i];
915  PZ = PZ + pho.phep[I - i][3 - i];
916  E = E + pho.phep[I - i][4 - i];
917 
918  P2 = pho.phep[I - i][1 - i] * pho.phep[I - i][1 - i] + pho.phep[I - i][2 - i] * pho.phep[I - i][2 - i] + pho.phep[I - i][3 - i] *
919  pho.phep[I - i][3 - i];
920 
921  M = sqrt(fabs(pho.phep[I - i][4 - i] * pho.phep[I - i][4 - i] - P2));
922 
923  if (IPRINT == 1) fprintf(PHLUN, "CORRECTING MASS OF %6i: %14.9f => %14.9f\n", I, pho.phep[I - i][5 - i], M);
924 
925  pho.phep[I - i][5 - i] = M;
926 
927  }
928 
929  } else if (MODOP == 4) {
930  // -----------------------
931 
932  // lets loop thru all daughters and correct their masses
933  // or energies according to E^2=p^2+m^2
934 
935  for (I = 3; I <= pho.nhep; I++) {
936 
937  PX = PX + pho.phep[I - i][1 - i];
938  PY = PY + pho.phep[I - i][2 - i];
939  PZ = PZ + pho.phep[I - i][3 - i];
940  P2 = pho.phep[I - i][1 - i] * pho.phep[I - i][1 - i] + pho.phep[I - i][2 - i] * pho.phep[I - i][2 - i] + pho.phep[I - i][3 - i] *
941  pho.phep[I - i][3 - i];
942  M = sqrt(fabs(pho.phep[I - i][4 - i] * pho.phep[I - i][4 - i] - P2));
943 
944 
945  if (M > MCUT) {
946  if (IPRINT == 1) fprintf(PHLUN, "CORRECTING MASS OF %6i: %14.9f => %14.9f\n", I, pho.phep[I - i][5 - i], M);
947  pho.phep[I - i][5 - i] = M;
948  E = E + pho.phep[I - i][4 - i];
949  } else {
950 
951  EN = sqrt(pho.phep[I - i][5 - i] * pho.phep[I - i][5 - i] + P2);
952  if (IPRINT == 1) fprintf(PHLUN, "CORRECTING ENERGY OF %6i: %14.9f =>% 14.9f\n", I , pho.phep[I - i][4 - i], EN);
953 
954  pho.phep[I - i][4 - i] = EN;
955  E = E + pho.phep[I - i][4 - i];
956  }
957 
958 
959  }
960  }
961 
962  // -----
963 
964  if (IPRINT == 1) {
965  fprintf(PHLUN, "CORRECTING MOTHER");
966  fprintf(PHLUN, "PX:%14.9f =>%14.9f", pho.phep[1 - i][1 - i], PX - pho.phep[2 - i][1 - i]);
967  fprintf(PHLUN, "PY:%14.9f =>%14.9f", pho.phep[1 - i][2 - i], PY - pho.phep[2 - i][2 - i]);
968  fprintf(PHLUN, "PZ:%14.9f =>%14.9f", pho.phep[1 - i][3 - i], PZ - pho.phep[2 - i][3 - i]);
969  fprintf(PHLUN, " E:%14.9f =>%14.9f", pho.phep[1 - i][4 - i], E - pho.phep[2 - i][4 - i]);
970  }
971 
972  pho.phep[1 - i][1 - i] = PX - pho.phep[2 - i][1 - i];
973  pho.phep[1 - i][2 - i] = PY - pho.phep[2 - i][2 - i];
974  pho.phep[1 - i][3 - i] = PZ - pho.phep[2 - i][3 - i];
975  pho.phep[1 - i][4 - i] = E - pho.phep[2 - i][4 - i];
976 
977 
978  P2 = pho.phep[1 - i][1 - i] * pho.phep[1 - i][1 - i] + pho.phep[1 - i][2 - i] * pho.phep[1 - i][2 - i] + pho.phep[1 - i][3 - i] *
979  pho.phep[1 - i][3 - i];
980  if (pho.phep[1 - i][4 - i]*pho.phep[1 - i][4 - i] > P2) {
981  M = sqrt(pho.phep[1 - i][4 - i] * pho.phep[1 - i][4 - i] - P2);
982  if (IPRINT == 1)fprintf(PHLUN, " M: %14.9f => %14.9f\n", pho.phep[1 - i][5 - i], M);
983  pho.phep[1 - i][5 - i] = M;
984  }
985 
986  PHLUPA(25);
987 
988  }
989 
990 
991 
992 
993 
994 
995 //----------------------------------------------------------------------
996 //
997 // PHOTOS: PHOton radiation in decays DOing of KINematics
998 //
999 // Purpose: Starting from the charged particle energy/momentum,
1000 // PNEUTR, photon energy fraction and photon angle with
1001 // respect to the axis formed by charged particle energy/
1002 // momentum vector and PNEUTR, scale the energy/momentum,
1003 // keeping the original direction of the neutral system in
1004 // the lab. frame untouched.
1005 //
1006 // Input Parameters: IP: Pointer to decaying particle in
1007 // /PHOEVT/ and the common itself
1008 // NCHARB: pointer to the charged radiating
1009 // daughter in /PHOEVT/.
1010 // NEUDAU: pointer to the first neutral daughter
1011 // Output Parameters: Common /PHOEVT/, with photon added.
1012 //
1013 // Author(s): Z. Was, B. van Eijk Created at: 26/11/89
1014 // Last Update: 27/05/93
1015 //
1016 //----------------------------------------------------------------------
1017 
1018  void PHODO(int IP, int NCHARB, int NEUDAU)
1019  {
1020  static int i = 1;
1021  double QNEW, QOLD, EPHOTO, PMAVIR;
1022  double GNEUT, DATA;
1023  double CCOSTH, SSINTH, PVEC[4], PARNE;
1024  double TH3, FI3, TH4, FI4, FI5, ANGLE;
1025  int I, J, FIRST, LAST;
1026  double& COSTHG = phophs.costhg;
1027  double& SINTHG = phophs.sinthg;
1028  double& XPHOTO = phophs.xphoto;
1029  double* PNEUTR = phomom.pneutr;
1030  double& MNESQR = phomom.mnesqr;
1031 
1032  //--
1033  EPHOTO = XPHOTO * pho.phep[IP - i][5 - i] / 2.0;
1034  PMAVIR = sqrt(pho.phep[IP - i][5 - i] * (pho.phep[IP - i][5 - i] - 2.0 * EPHOTO));
1035  //--
1036  //-- Reconstruct kinematics of charged particle and neutral system
1037  phorest.fi1 = PHOAN1(PNEUTR[1 - i], PNEUTR[2 - i]);
1038  //--
1039  //-- Choose axis along z of PNEUTR, calculate angle between x and y
1040  //-- components and z and x-y plane and perform Lorentz transform...
1041  phorest.th1 = PHOAN2(PNEUTR[3 - i], sqrt(PNEUTR[1 - i] * PNEUTR[1 - i] + PNEUTR[2 - i] * PNEUTR[2 - i]));
1042  PHORO3(-phorest.fi1, PNEUTR);
1043  PHORO2(-phorest.th1, PNEUTR);
1044  //--
1045  //-- Take away photon energy from charged particle and PNEUTR ! Thus
1046  //-- the onshell charged particle decays into virtual charged particle
1047  //-- and photon. The virtual charged particle mass becomes:
1048  //-- SQRT(pho.phep[5,IP)*(pho.phep[5,IP)-2*EPHOTO)). Construct new PNEUTR mo-
1049  //-- mentum in the rest frame of the parent:
1050  //-- 1) Scaling parameters...
1051  QNEW = PHOTRI(PMAVIR, PNEUTR[5 - i], pho.phep[NCHARB - i][5 - i]);
1052  QOLD = PNEUTR[3 - i];
1053  GNEUT = (QNEW * QNEW + QOLD * QOLD + MNESQR) / (QNEW * QOLD + sqrt((QNEW * QNEW + MNESQR) * (QOLD * QOLD + MNESQR)));
1054  if (GNEUT < 1.0) {
1055  DATA = 0.0;
1056  PHOERR(4, "PHOKIN", DATA);
1057  }
1058  PARNE = GNEUT - sqrt(max(GNEUT * GNEUT - 1.0, 0.0));
1059  //--
1060  //-- 2) ...reductive boost...
1061  PHOBO3(PARNE, PNEUTR);
1062  //--
1063  //-- ...calculate photon energy in the reduced system...
1064  pho.nhep = pho.nhep + 1;
1065  pho.isthep[pho.nhep - i] = 1;
1066  pho.idhep[pho.nhep - i] = 22;
1067  //-- Photon mother and daughter pointers !
1068  pho.jmohep[pho.nhep - i][1 - i] = IP;
1069  pho.jmohep[pho.nhep - i][2 - i] = 0;
1070  pho.jdahep[pho.nhep - i][1 - i] = 0;
1071  pho.jdahep[pho.nhep - i][2 - i] = 0;
1072  pho.phep[pho.nhep - i][4 - i] = EPHOTO * pho.phep[IP - i][5 - i] / PMAVIR;
1073  //--
1074  //-- ...and photon momenta
1075  CCOSTH = -COSTHG;
1076  SSINTH = SINTHG;
1077  TH3 = PHOAN2(CCOSTH, SSINTH);
1078  FI3 = TWOPI * Photos::randomDouble();
1079  pho.phep[pho.nhep - i][1 - i] = pho.phep[pho.nhep - i][4 - i] * SINTHG * cos(FI3);
1080  pho.phep[pho.nhep - i][2 - i] = pho.phep[pho.nhep - i][4 - i] * SINTHG * sin(FI3);
1081  //--
1082  //-- Minus sign because axis opposite direction of charged particle !
1083  pho.phep[pho.nhep - i][3 - i] = -pho.phep[pho.nhep - i][4 - i] * COSTHG;
1084  pho.phep[pho.nhep - i][5 - i] = 0.0;
1085  //--
1086  //-- Rotate in order to get photon along z-axis
1087  PHORO3(-FI3, PNEUTR);
1088  PHORO3(-FI3, pho.phep[pho.nhep - i]);
1089  PHORO2(-TH3, PNEUTR);
1090  PHORO2(-TH3, pho.phep[pho.nhep - i]);
1091  ANGLE = EPHOTO / pho.phep[pho.nhep - i][4 - i];
1092  //--
1093  //-- Boost to the rest frame of decaying particle
1094  PHOBO3(ANGLE, PNEUTR);
1095  PHOBO3(ANGLE, pho.phep[pho.nhep - i]);
1096  //--
1097  //-- Back in the parent rest frame but PNEUTR not yet oriented !
1098  FI4 = PHOAN1(PNEUTR[1 - i], PNEUTR[2 - i]);
1099  TH4 = PHOAN2(PNEUTR[3 - i], sqrt(PNEUTR[1 - i] * PNEUTR[1 - i] + PNEUTR[2 - i] * PNEUTR[2 - i]));
1100  PHORO3(FI4, PNEUTR);
1101  PHORO3(FI4, pho.phep[pho.nhep - i]);
1102  //--
1103  for (I = 2; I <= 4; I++) PVEC[I - i] = 0.0;
1104  PVEC[1 - i] = 1.0;
1105 
1106  PHORO3(-FI3, PVEC);
1107  PHORO2(-TH3, PVEC);
1108  PHOBO3(ANGLE, PVEC);
1109  PHORO3(FI4, PVEC);
1110  PHORO2(-TH4, PNEUTR);
1111  PHORO2(-TH4, pho.phep[pho.nhep - i]);
1112  PHORO2(-TH4, PVEC);
1113  FI5 = PHOAN1(PVEC[1 - i], PVEC[2 - i]);
1114  //--
1115  //-- Charged particle restores original direction
1116  PHORO3(-FI5, PNEUTR);
1117  PHORO3(-FI5, pho.phep[pho.nhep - i]);
1118  PHORO2(phorest.th1, PNEUTR);
1119  PHORO2(phorest.th1, pho.phep[pho.nhep - i]);
1120  PHORO3(phorest.fi1, PNEUTR);
1121  PHORO3(phorest.fi1, pho.phep[pho.nhep - i]);
1122  //-- See whether neutral system has multiplicity larger than 1...
1123 
1124  if ((pho.jdahep[IP - i][2 - i] - pho.jdahep[IP - i][1 - i]) > 1) {
1125  //-- Find pointers to components of 'neutral' system
1126  //--
1127  FIRST = NEUDAU;
1128  LAST = pho.jdahep[IP - i][2 - i];
1129  for (I = FIRST; I <= LAST; I++) {
1130  if (I != NCHARB && (pho.jmohep[I - i][1 - i] == IP)) {
1131  //--
1132  //-- Reconstruct kinematics...
1133  PHORO3(-phorest.fi1, pho.phep[I - i]);
1134  PHORO2(-phorest.th1, pho.phep[I - i]);
1135  //--
1136  //-- ...reductive boost
1137  PHOBO3(PARNE, pho.phep[I - i]);
1138  //--
1139  //-- Rotate in order to get photon along z-axis
1140  PHORO3(-FI3, pho.phep[I - i]);
1141  PHORO2(-TH3, pho.phep[I - i]);
1142  //--
1143  //-- Boost to the rest frame of decaying particle
1144  PHOBO3(ANGLE, pho.phep[I - i]);
1145  //--
1146  //-- Back in the parent rest-frame but PNEUTR not yet oriented.
1147  PHORO3(FI4, pho.phep[I - i]);
1148  PHORO2(-TH4, pho.phep[I - i]);
1149  //--
1150  //-- Charged particle restores original direction
1151  PHORO3(-FI5, pho.phep[I - i]);
1152  PHORO2(phorest.th1, pho.phep[I - i]);
1153  PHORO3(phorest.fi1, pho.phep[I - i]);
1154  }
1155  }
1156  } else {
1157  //--
1158  // ...only one 'neutral' particle in addition to photon!
1159  for (J = 1; J <= 4; J++) pho.phep[NEUDAU - i][J - i] = PNEUTR[J - i];
1160  }
1161  //--
1162  //-- All 'neutrals' treated, fill /PHOEVT/ for charged particle...
1163  for (J = 1; J <= 3; J++) pho.phep[NCHARB - i][J - i] = -(pho.phep[pho.nhep - i][J - i] + PNEUTR[J - i]);
1164  pho.phep[NCHARB - i][4 - i] = pho.phep[IP - i][5 - i] - (pho.phep[pho.nhep - i][4 - i] + PNEUTR[4 - i]);
1165  //--
1166  }
1167 
1168 
1169 //----------------------------------------------------------------------
1170 //
1171 // PHOTOS: PHOtos Boson W correction weight
1172 //
1173 // Purpose: calculates correction weight due to amplitudes of
1174 // emission from W boson.
1175 //
1176 //
1177 //
1178 //
1179 //
1180 // Input Parameters: Common /PHOEVT/, with photon added.
1181 // wt to be corrected
1182 //
1183 //
1184 //
1185 // Output Parameters: wt
1186 //
1187 // Author(s): G. Nanava, Z. Was Created at: 13/03/03
1188 // Last Update: 08/07/13
1189 //
1190 //----------------------------------------------------------------------
1191 
1192  void PHOBW(double* WT)
1193  {
1194  static int i = 1;
1195  int I;
1196  double EMU, MCHREN, BETA, COSTHG, MPASQR, XPH;
1197  //--
1198  if (abs(pho.idhep[1 - i]) == 24 &&
1199  abs(pho.idhep[pho.jdahep[1 - i][1 - i] - i]) >= 11 &&
1200  abs(pho.idhep[pho.jdahep[1 - i][1 - i] - i]) <= 16 &&
1201  abs(pho.idhep[pho.jdahep[1 - i][1 - i] + 1 - i]) >= 11 &&
1202  abs(pho.idhep[pho.jdahep[1 - i][1 - i] + 1 - i]) <= 16) {
1203 
1204  if (
1205  abs(pho.idhep[pho.jdahep[1 - i][1 - i] - i]) == 11 ||
1206  abs(pho.idhep[pho.jdahep[1 - i][1 - i] - i]) == 13 ||
1207  abs(pho.idhep[pho.jdahep[1 - i][1 - i] - i]) == 15) {
1208  I = pho.jdahep[1 - i][1 - i];
1209  } else {
1210  I = pho.jdahep[1 - i][1 - i] + 1;
1211  }
1212 
1213  EMU = pho.phep[I - i][4 - i];
1214  MCHREN = fabs(pow(pho.phep[I - i][4 - i], 2) - pow(pho.phep[I - i][3 - i], 2)
1215  - pow(pho.phep[I - i][2 - i], 2) - pow(pho.phep[I - i][1 - i], 2));
1216  BETA = sqrt(1.0 - MCHREN / pho.phep[I - i][4 - i] / pho.phep[I - i][4 - i]);
1217  COSTHG = (pho.phep[I - i][3 - i] * pho.phep[pho.nhep - i][3 - i] + pho.phep[I - i][2 - i] * pho.phep[pho.nhep - i][2 - i]
1218  + pho.phep[I - i][1 - i] * pho.phep[pho.nhep - i][1 - i]) /
1219  sqrt(pho.phep[I - i][3 - i] * pho.phep[I - i][3 - i] + pho.phep[I - i][2 - i] * pho.phep[I - i][2 - i] + pho.phep[I - i][1 - i] *
1220  pho.phep[I - i][1 - i]) /
1221  sqrt(pho.phep[pho.nhep - i][3 - i] * pho.phep[pho.nhep - i][3 - i] + pho.phep[pho.nhep - i][2 - i] * pho.phep[pho.nhep - i][2 - i] +
1222  pho.phep[pho.nhep - i][1 - i] * pho.phep[pho.nhep - i][1 - i]);
1223  MPASQR = pho.phep[1 - i][4 - i] * pho.phep[1 - i][4 - i];
1224  XPH = pho.phep[pho.nhep - i][4 - i];
1225  *WT = (*WT) * (1 - 8 * EMU * XPH * (1 - COSTHG * BETA) *
1226  (MCHREN + 2 * XPH * sqrt(MPASQR)) /
1227  (MPASQR * MPASQR) / (1 - MCHREN / MPASQR) / (4 - MCHREN / MPASQR));
1228  }
1229  // write(*,*) pho.idhep[1),pho.idhep[pho.jdahep[1,1)),pho.idhep[pho.jdahep[1,1)+1)
1230  // write(*,*) emu,xph,costhg,beta,mpasqr,mchren
1231 
1232  }
1233 
1234 
1235 
1236 //----------------------------------------------------------------------
1237 //
1238 // PHOTOS: PHOton radiation in decays control FACtor
1239 //
1240 // Purpose: This is the control function for the photon spectrum and
1241 // final weighting. It is called from PHOENE for genera-
1242 // ting the raw photon energy spectrum (MODE=0) and in PHO-
1243 // COR to scale the final weight (MODE=1). The factor con-
1244 // sists of 3 terms. Addition of the factor FF which mul-
1245 // tiplies PHOFAC for MODE=0 and divides PHOFAC for MODE=1,
1246 // does not affect the results for the MC generation. An
1247 // appropriate choice for FF can speed up the calculation.
1248 // Note that a too small value of FF may cause weight over-
1249 // flow in PHOCOR and will generate a warning, halting the
1250 // execution. PRX should be included for repeated calls
1251 // for the same event, allowing more particles to radiate
1252 // photons. At the first call IREP=0, for more than 1
1253 // charged decay products, IREP >= 1. Thus, PRSOFT (no
1254 // photon radiation probability in the previous calls)
1255 // appropriately scales the strength of the bremsstrahlung.
1256 //
1257 // Input Parameters: MODE, PROBH, XF
1258 //
1259 // Output Parameter: Function value
1260 //
1261 // Author(s): S. Jadach, Z. Was Created at: 01/01/89
1262 // B. van Eijk, P.Golonka Last Update: 09/07/13
1263 //
1264 //----------------------------------------------------------------------
1265 
1266  double PHOFAC(int MODE)
1267  {
1268  static double FF = 0.0, PRX = 0.0;
1269 
1270  if (phokey.iexp) return 1.0; // In case of exponentiation this routine is useles
1271 
1272  if (MODE == -1) {
1273  PRX = 1.0;
1274  FF = 1.0;
1275  phopro.probh = 0.0;
1276  } else if (MODE == 0) {
1277  if (phopro.irep == 0) PRX = 1.0;
1278  PRX = PRX / (1.0 - phopro.probh);
1279  FF = 1.0;
1280  //--
1281  //-- Following options are not considered for the time being...
1282  //-- (1) Good choice, but does not save very much time:
1283  //-- FF=(1.0-sqrt(phopro.xf)/2.0)/(1.0+sqrt(phopro.xf)/2.0)
1284  //-- (2) Taken from the blue, but works without weight overflows...
1285  //-- FF=(1.0-phopro.xf/(1-pow((1-sqrt(phopro.xf)),2)))*(1+(1-sqrt(phopro.xf))/sqrt(1-phopro.xf))/2.0
1286  return FF * PRX;
1287  } else {
1288  return 1.0 / FF;
1289  }
1290 
1291  return NAN;
1292  }
1293 
1294 
1295 
1296 // ######
1297 // replace with,
1298 // ######
1299 
1300 //----------------------------------------------------------------------
1301 //
1302 // PHOTOS: PHOton radiation in decays CORrection weight from
1303 // matrix elements This version for spin 1/2 is verified for
1304 // W decay only
1305 // Purpose: Calculate photon angle. The reshaping functions will
1306 // have to depend on the spin S of the charged particle.
1307 // We define: ME = 2 * S + 1 !
1308 // THIS IS POSSIBLY ALWAYS BETTER THAN PHOCOR()
1309 //
1310 // Input Parameters: MPASQR: Parent mass squared,
1311 // MCHREN: Renormalised mass of charged system,
1312 // ME: 2 * spin + 1 determines matrix element
1313 //
1314 // Output Parameter: Function value.
1315 //
1316 // Author(s): Z. Was, B. van Eijk, G. Nanava Created at: 26/11/89
1317 // Last Update: 01/11/12
1318 //
1319 //----------------------------------------------------------------------
1320 
1321  double PHOCORN(double MPASQR, double MCHREN, int ME)
1322  {
1323  double wt1, wt2, wt3;
1324  double beta0, beta1, XX, YY, DATA;
1325  double S1, PHOCOR;
1326  double& COSTHG = phophs.costhg;
1327  double& XPHMAX = phophs.xphmax;
1328  double& XPHOTO = phophs.xphoto;
1329  double& MCHSQR = phomom.mchsqr;
1330  double& MNESQR = phomom.mnesqr;
1331 
1332 
1333 
1334  //--
1335  //-- Shaping (modified by ZW)...
1336  XX = 4.0 * MCHSQR / MPASQR * (1.0 - XPHOTO) / pow(1.0 - XPHOTO + (MCHSQR - MNESQR) / MPASQR, 2);
1337  if (ME == 1) {
1338  S1 = MPASQR * (1.0 - XPHOTO);
1339  beta0 = 2 * PHOTRI(1.0, sqrt(MCHSQR / MPASQR), sqrt(MNESQR / MPASQR));
1340  beta1 = 2 * PHOTRI(1.0, sqrt(MCHSQR / S1), sqrt(MNESQR / S1));
1341  wt1 = (1.0 - COSTHG * sqrt(1.0 - MCHREN))
1342  / ((1.0 + pow(1.0 - XPHOTO / XPHMAX, 2)) / 2.0) * XPHOTO; // de-presampler
1343 
1344  wt2 = beta1 / beta0 * XPHOTO; //phase space jacobians
1345  wt3 = beta1 * beta1 * (1.0 - COSTHG * COSTHG) * (1.0 - XPHOTO) / XPHOTO / XPHOTO
1346  / pow(1.0 + MCHSQR / S1 - MNESQR / S1 - beta1 * COSTHG, 2) / 2.0; // matrix element
1347  } else if (ME == 2) {
1348  S1 = MPASQR * (1.0 - XPHOTO);
1349  beta0 = 2 * PHOTRI(1.0, sqrt(MCHSQR / MPASQR), sqrt(MNESQR / MPASQR));
1350  beta1 = 2 * PHOTRI(1.0, sqrt(MCHSQR / S1), sqrt(MNESQR / S1));
1351  wt1 = (1.0 - COSTHG * sqrt(1.0 - MCHREN))
1352  / ((1.0 + pow(1.0 - XPHOTO / XPHMAX, 2)) / 2.0) * XPHOTO; // de-presampler
1353 
1354  wt2 = beta1 / beta0 * XPHOTO; // phase space jacobians
1355 
1356  wt3 = beta1 * beta1 * (1.0 - COSTHG * COSTHG) * (1.0 - XPHOTO) / XPHOTO / XPHOTO // matrix element
1357  / pow(1.0 + MCHSQR / S1 - MNESQR / S1 - beta1 * COSTHG, 2) / 2.0 ;
1358  wt3 = wt3 * (1 - XPHOTO / XPHMAX + 0.5 * pow(XPHOTO / XPHMAX, 2)) / (1 - XPHOTO / XPHMAX);
1359  // print*,"wt3=",wt3
1360  phocorwt.phocorwt3 = wt3;
1361  phocorwt.phocorwt2 = wt2;
1362  phocorwt.phocorwt1 = wt1;
1363 
1364  // YY=0.5D0*(1.D0-XPHOTO/XPHMAX+1.D0/(1.D0-XPHOTO/XPHMAX))
1365  // phwt.beta=SQRT(1.D0-XX)
1366  // wt1=(1.D0-COSTHG*SQRT(1.D0-MCHREN))/(1.D0-COSTHG*phwt.beta)
1367  // wt2=(1.D0-XX/YY/(1.D0-phwt.beta**2*COSTHG**2))*(1.D0+COSTHG*phwt.beta)/2.D0
1368  // wt3=1.D0
1369  } else if ((ME == 3) || (ME == 4) || (ME == 5)) {
1370  YY = 1.0;
1371  phwt.beta = sqrt(1.0 - XX);
1372  wt1 = (1.0 - COSTHG * sqrt(1.0 - MCHREN)) / (1.0 - COSTHG * phwt.beta);
1373  wt2 = (1.0 - XX / YY / (1.0 - phwt.beta * phwt.beta * COSTHG * COSTHG)) * (1.0 + COSTHG * phwt.beta) / 2.0;
1374  wt3 = (1.0 + pow(1.0 - XPHOTO / XPHMAX, 2) - pow(XPHOTO / XPHMAX, 3)) /
1375  (1.0 + pow(1.0 - XPHOTO / XPHMAX, 2));
1376  } else {
1377  DATA = (ME - 1.0) / 2.0;
1378  PHOERR(6, "PHOCORN", DATA);
1379  YY = 1.0;
1380  phwt.beta = sqrt(1.0 - XX);
1381  wt1 = (1.0 - COSTHG * sqrt(1.0 - MCHREN)) / (1.0 - COSTHG * phwt.beta);
1382  wt2 = (1.0 - XX / YY / (1.0 - phwt.beta * phwt.beta * COSTHG * COSTHG)) * (1.0 + COSTHG * phwt.beta) / 2.0;
1383  wt3 = 1.0;
1384  }
1385  wt2 = wt2 * PHOFAC(1);
1386  PHOCOR = wt1 * wt2 * wt3;
1387 
1388  phopro.corwt = PHOCOR;
1389  if (PHOCOR > 1.0) {
1390  DATA = PHOCOR;
1391  PHOERR(3, "PHOCOR", DATA);
1392  }
1393  return PHOCOR;
1394  }
1395 
1396 
1397 
1398 
1399 
1400 //----------------------------------------------------------------------
1401 //
1402 // PHOTOS: PHOton radiation in decays CORrection weight from
1403 // matrix elements
1404 //
1405 // Purpose: Calculate photon angle. The reshaping functions will
1406 // have to depend on the spin S of the charged particle.
1407 // We define: ME = 2 * S + 1 !
1408 //
1409 // Input Parameters: MPASQR: Parent mass squared,
1410 // MCHREN: Renormalised mass of charged system,
1411 // ME: 2 * spin + 1 determines matrix element
1412 //
1413 // Output Parameter: Function value.
1414 //
1415 // Author(s): Z. Was, B. van Eijk Created at: 26/11/89
1416 // Last Update: 21/03/93
1417 //
1418 //----------------------------------------------------------------------
1419 
1420  double PHOCOR(double MPASQR, double MCHREN, int ME)
1421  {
1422  double XX, YY, DATA;
1423  double PHOC;
1424  int IscaNLO;
1425  double& COSTHG = phophs.costhg;
1426  double& XPHMAX = phophs.xphmax;
1427  double& XPHOTO = phophs.xphoto;
1428  double& MCHSQR = phomom.mchsqr;
1429  double& MNESQR = phomom.mnesqr;
1430 
1431 
1432  //--
1433  //-- Shaping (modified by ZW)...
1434  XX = 4.0 * MCHSQR / MPASQR * (1.0 - XPHOTO) / pow((1.0 - XPHOTO + (MCHSQR - MNESQR) / MPASQR), 2);
1435  if (ME == 1) {
1436  YY = 1.0;
1437  phwt.wt3 = (1.0 - XPHOTO / XPHMAX) / ((1.0 + pow((1.0 - XPHOTO / XPHMAX), 2)) / 2.0);
1438  } else if (ME == 2) {
1439  YY = 0.5 * (1.0 - XPHOTO / XPHMAX + 1.0 / (1.0 - XPHOTO / XPHMAX));
1440  phwt.wt3 = 1.0;
1441  } else if ((ME == 3) || (ME == 4) || (ME == 5)) {
1442  YY = 1.0;
1443  phwt.wt3 = (1.0 + pow(1.0 - XPHOTO / XPHMAX, 2) - pow(XPHOTO / XPHMAX, 3)) /
1444  (1.0 + pow(1.0 - XPHOTO / XPHMAX, 2));
1445  } else {
1446  DATA = (ME - 1.0) / 2.0;
1447  PHOERR(6, "PHOCOR", DATA);
1448  YY = 1.0;
1449  phwt.wt3 = 1.0;
1450  }
1451 
1452 
1453  phwt.beta = sqrt(1.0 - XX);
1454  phwt.wt1 = (1.0 - COSTHG * sqrt(1.0 - MCHREN)) / (1.0 - COSTHG * phwt.beta);
1455  phwt.wt2 = (1.0 - XX / YY / (1.0 - phwt.beta * phwt.beta * COSTHG * COSTHG)) * (1.0 + COSTHG * phwt.beta) / 2.0;
1456 
1457 
1458  IscaNLO = Photos::meCorrectionWtForScalar;
1459  if (ME == 1 && IscaNLO == 1) { // this switch NLO in scalar decays.
1460  // overrules default calculation.
1461  // Need tests including basic ones
1462  PHOC = PHOCORN(MPASQR, MCHREN, ME);
1463  phwt.wt1 = 1.0;
1464  phwt.wt2 = 1.0;
1465  phwt.wt3 = PHOC;
1466  } else {
1467  phwt.wt2 = phwt.wt2 * PHOFAC(1);
1468  }
1469  PHOC = phwt.wt1 * phwt.wt2 * phwt.wt3;
1470 
1471  phopro.corwt = PHOC;
1472  if (PHOC > 1.0) {
1473  DATA = PHOC;
1474  PHOERR(3, "PHOCOR", DATA);
1475  }
1476  return PHOC;
1477  }
1478 
1479 
1480 //----------------------------------------------------------------------
1481 //
1482 // PHOTWO: PHOtos but TWO mothers allowed
1483 //
1484 // Purpose: Combines two mothers into one in /PHOEVT/
1485 // necessary eg in case of g g (q qbar) --> t tbar
1486 //
1487 // Input Parameters: Common /PHOEVT/ (/PHOCMS/)
1488 //
1489 // Output Parameters: Common /PHOEVT/, (stored mothers)
1490 //
1491 // Author(s): Z. Was Created at: 5/08/93
1492 // Last Update:10/08/93
1493 //
1494 //----------------------------------------------------------------------
1495 
1496  void PHOTWO(int MODE)
1497  {
1498 
1499  int I;
1500  static int i = 1;
1501  double MPASQR;
1502  bool IFRAD;
1503  // logical IFRAD is used to tag cases when two mothers may be
1504  // merged to the sole one.
1505  // So far used in case:
1506  // 1) of t tbar production
1507  //
1508  // t tbar case
1509  if (MODE == 0) {
1510  IFRAD = (pho.idhep[1 - i] == 21) && (pho.idhep[2 - i] == 21);
1511  IFRAD = IFRAD || (pho.idhep[1 - i] == -pho.idhep[2 - i] && abs(pho.idhep[1 - i]) <= 6);
1512  IFRAD = IFRAD && (abs(pho.idhep[3 - i]) == 6) && (abs(pho.idhep[4 - i]) == 6);
1513  MPASQR = pow(pho.phep[1 - i][4 - i] + pho.phep[2 - i][4 - i], 2) - pow(pho.phep[1 - i][3 - i] + pho.phep[2 - i][3 - i], 2)
1514  - pow(pho.phep[1 - i][2 - i] + pho.phep[2 - i][2 - i], 2) - pow(pho.phep[1 - i][1 - i] + pho.phep[2 - i][1 - i], 2);
1515  IFRAD = IFRAD && (MPASQR > 0.0);
1516  if (IFRAD) {
1517  //.....combining first and second mother
1518  for (I = 1; I <= 4; I++) {
1519  pho.phep[1 - i][I - i] = pho.phep[1 - i][I - i] + pho.phep[2 - i][I - i];
1520  }
1521  pho.phep[1 - i][5 - i] = sqrt(MPASQR);
1522  //.....removing second mother,
1523  for (I = 1; I <= 5; I++) {
1524  pho.phep[2 - i][I - i] = 0.0;
1525  }
1526  }
1527  } else {
1528  // boosting of the mothers to the reaction frame not implemented yet.
1529  // to do it in mode 0 original mothers have to be stored in new comon (?)
1530  // and in mode 1 boosted to cms.
1531  }
1532  }
1533 
1534 
1535 
1536 //----------------------------------------------------------------------
1537 //
1538 // PHOTOS: PHOtos CDE-s
1539 //
1540 // Purpose: Keep definitions for PHOTOS QED correction Monte Carlo.
1541 //
1542 // Input Parameters: None
1543 //
1544 // Output Parameters: None
1545 //
1546 // Author(s): Z. Was, B. van Eijk Created at: 29/11/89
1547 // Last Update: 10/08/93
1548 //
1549 // =========================================================
1550 // General Structure Information: =
1551 // =========================================================
1552 //: ROUTINES:
1553 // 1) INITIALIZATION (all in C++ now)
1554 // 2) GENERAL INTERFACE:
1555 // PHOBOS
1556 // PHOIN
1557 // PHOTWO (specific interface
1558 // PHOOUT
1559 // PHOCHK
1560 // PHTYPE (specific interface
1561 // PHOMAK (specific interface
1562 // 3) QED PHOTON GENERATION:
1563 // PHINT
1564 // PHOBW
1565 // PHOPRE
1566 // PHOOMA
1567 // PHOENE
1568 // PHOCOR
1569 // PHOFAC
1570 // PHODO
1571 // 4) UTILITIES:
1572 // PHOTRI
1573 // PHOAN1
1574 // PHOAN2
1575 // PHOBO3
1576 // PHORO2
1577 // PHORO3
1578 // PHOCHA
1579 // PHOSPI
1580 // PHOERR
1581 // PHOREP
1582 // PHLUPA
1583 // PHCORK
1584 // IPHQRK
1585 // IPHEKL
1586 // COMMONS:
1587 // NAME USED IN SECT. # OF OC// Comment
1588 // PHOQED 1) 2) 3 Flags whether emisson to be gen.
1589 // PHOLUN 1) 4) 6 Output device number
1590 // PHOCOP 1) 3) 4 photon coupling & min energy
1591 // PHPICO 1) 3) 4) 5 PI & 2*PI
1592 // PHOSTA 1) 4) 3 Status information
1593 // PHOKEY 1) 2) 3) 7 Keys for nonstandard application
1594 // PHOVER 1) 1 Version info for outside
1595 // HEPEVT 2) 2 PDG common
1596 // PH_HEPEVT2) 8 PDG common internal
1597 // PHOEVT 2) 3) 10 PDG branch
1598 // PHOIF 2) 3) 2 emission flags for PDG branch
1599 // PHOMOM 3) 5 param of char-neutr system
1600 // PHOPHS 3) 5 photon momentum parameters
1601 // PHOPRO 3) 4 var. for photon rep. (in branch)
1602 // PHOCMS 2) 3 parameters of boost to branch CMS
1603 // PHNUM 4) 1 event number from outside
1604 //----------------------------------------------------------------------
1605 
1606 
1607 //----------------------------------------------------------------------
1608 //
1609 // PHOIN: PHOtos INput
1610 //
1611 // Purpose: copies IP branch of the common /PH_HEPEVT/ into /PHOEVT/
1612 // moves branch into its CMS system.
1613 //
1614 // Input Parameters: IP: pointer of particle starting branch
1615 // to be copied
1616 // BOOST: Flag whether boost to CMS was or was
1617 // . replace stri not performed.
1618 //
1619 // Output Parameters: Commons: /PHOEVT/, /PHOCMS/
1620 //
1621 // Author(s): Z. Was Created at: 24/05/93
1622 // Last Update: 16/11/93
1623 //
1624 //----------------------------------------------------------------------
1625  void PHOIN(int IP, bool* BOOST, int* NHEP0)
1626  {
1627  int FIRST, LAST, I, LL, IP2, J, NA;
1628  double PB;
1629  static int i = 1;
1630  int& nhep0 = *NHEP0;
1631  // double &BET[3]=BET;
1632  double& GAM = phocms.gam;
1633  double* BET = phocms.bet;
1634 
1635  //--
1636  // let-s calculate size of the little common entry
1637  FIRST = hep.jdahep[IP - i][1 - i];
1638  LAST = hep.jdahep[IP - i][2 - i];
1639  pho.nhep = 3 + LAST - FIRST + hep.nhep - nhep0;
1640  pho.nevhep = pho.nhep;
1641 
1642  // let-s take in decaying particle
1643  pho.idhep[1 - i] = hep.idhep[IP - i];
1644  pho.jdahep[1 - i][1 - i] = 3;
1645  pho.jdahep[1 - i][2 - i] = 3 + LAST - FIRST;
1646  for (I = 1; I <= 5; I++) pho.phep[1 - i][I - i] = hep.phep[IP - i][I - i];
1647 
1648  // let-s take in eventual second mother
1649  IP2 = hep.jmohep[hep.jdahep[IP - i][1 - i] - i][2 - i];
1650  if ((IP2 != 0) && (IP2 != IP)) {
1651  pho.idhep[2 - i] = hep.idhep[IP2 - i];
1652  pho.jdahep[2 - i][1 - i] = 3;
1653  pho.jdahep[2 - i][2 - i] = 3 + LAST - FIRST;
1654  for (I = 1; I <= 5; I++)
1655  pho.phep[2 - i][I - i] = hep.phep[IP2 - i][I - i];
1656  } else {
1657  pho.idhep[2 - i] = 0;
1658  for (I = 1; I <= 5; I++) pho.phep[2 - i][I - i] = 0.0;
1659  }
1660 
1661  // let-s take in daughters
1662  for (LL = 0; LL <= LAST - FIRST; LL++) {
1663  pho.idhep[3 + LL - i] = hep.idhep[FIRST + LL - i];
1664  pho.jmohep[3 + LL - i][1 - i] = hep.jmohep[FIRST + LL - i][1 - i];
1665  if (hep.jmohep[FIRST + LL - i][1 - i] == IP) pho.jmohep[3 + LL - i][1 - i] = 1;
1666  for (I = 1; I <= 5; I++) pho.phep[3 + LL - i][I - i] = hep.phep[FIRST + LL - i][I - i];
1667 
1668  }
1669  if (hep.nhep > nhep0) {
1670  // let-s take in illegitimate daughters
1671  NA = 3 + LAST - FIRST;
1672  for (LL = 1; LL <= hep.nhep - nhep0; LL++) {
1673  pho.idhep[NA + LL - i] = hep.idhep[nhep0 + LL - i];
1674  pho.jmohep[NA + LL - i][1 - i] = hep.jmohep[nhep0 + LL - i][1 - i];
1675  if (hep.jmohep[nhep0 + LL - i][1 - i] == IP) pho.jmohep[NA + LL - i][1 - i] = 1;
1676  for (I = 1; I <= 5; I++) pho.phep[NA + LL - i][I - i] = hep.phep[nhep0 + LL - i][I - i];
1677 
1678  }
1679  //-- there is hep.nhep-nhep0 daugters more.
1680  pho.jdahep[1 - i][2 - i] = 3 + LAST - FIRST + hep.nhep - nhep0;
1681  }
1682  if (pho.idhep[pho.nhep - i] == 22) PHLUPA(100001);
1683  // if (pho.idhep[pho.nhep-i]==22) exit(-1);
1684  PHCORK(0);
1685  if (pho.idhep[pho.nhep - i] == 22) PHLUPA(100002);
1686 
1687  // special case of t tbar production process
1688  if (phokey.iftop) PHOTWO(0);
1689  *BOOST = false;
1690 
1691  //-- Check whether parent is in its rest frame...
1692  // ZBW ND 27.07.2009:
1693  // bug reported by Vladimir Savinov localized and fixed.
1694  // protection against rounding error was back-firing if soft
1695  // momentum of mother was physical. Consequence was that PHCORK was
1696  // messing up masses of final state particles in vertex of the decay.
1697  // Only configurations with previously generated photons of energy fraction
1698  // smaller than 0.0001 were affected. Effect was numerically insignificant.
1699 
1700  // IF ( (ABS(pho.phep[4,1)-pho.phep[5,1)).GT.pho.phep[5,1)*1.D-8)
1701  // $ .AND.(pho.phep[5,1).NE.0)) THEN
1702 
1703  if ((fabs(pho.phep[1 - i][1 - i] + fabs(pho.phep[1 - i][2 - i]) + fabs(pho.phep[1 - i][3 - i])) >
1704  pho.phep[1 - i][5 - i] * 1.E-8) && (pho.phep[1 - i][5 - i] != 0)) {
1705 
1706  *BOOST = true;
1707  //PHOERR(404,"PHOIN",1.0); // we need to improve this warning: program should never
1708  // enter this place
1709  // may be exit(-1);
1710  //--
1711  //-- Boost daughter particles to rest frame of parent...
1712  //-- Resultant neutral system already calculated in rest frame !
1713  for (J = 1; J <= 3; J++) BET[J - i] = -pho.phep[1 - i][J - i] / pho.phep[1 - i][5 - i];
1714  GAM = pho.phep[1 - i][4 - i] / pho.phep[1 - i][5 - i];
1715  for (I = pho.jdahep[1 - i][1 - i]; I <= pho.jdahep[1 - i][2 - i]; I++) {
1716  PB = BET[1 - i] * pho.phep[I - i][1 - i] + BET[2 - i] * pho.phep[I - i][2 - i] + BET[3 - i] * pho.phep[I - i][3 - i];
1717  for (J = 1; J <= 3;
1718  J++) pho.phep[I - i][J - i] = pho.phep[I - i][J - i] + BET[J - i] * (pho.phep[I - i][4 - i] + PB / (GAM + 1.0));
1719  pho.phep[I - i][4 - i] = GAM * pho.phep[I - i][4 - i] + PB;
1720  }
1721  //-- Finally boost mother as well
1722  I = 1;
1723  PB = BET[1 - i] * pho.phep[I - i][1 - i] + BET[2 - i] * pho.phep[I - i][2 - i] + BET[3 - i] * pho.phep[I - i][3 - i];
1724  for (J = 1; J <= 3; J++) pho.phep[I - i][J - i] = pho.phep[I - i][J - i] + BET[J - i] * (pho.phep[I - i][4 - i] + PB / (GAM + 1.0));
1725 
1726  pho.phep[I - i][4 - i] = GAM * pho.phep[I - i][4 - i] + PB;
1727  }
1728 
1729 
1730  // special case of t tbar production process
1731  if (phokey.iftop) PHOTWO(1);
1732  PHLUPA(2);
1733  if (pho.idhep[pho.nhep - i] == 22) PHLUPA(10000);
1734  //if (pho.idhep[pho.nhep-1-i]==22) exit(-1); // this is probably form very old times ...
1735  return;
1736  }
1737 
1738 
1739 //----------------------------------------------------------------------
1740 //
1741 // PHOOUT: PHOtos OUTput
1742 //
1743 // Purpose: copies back IP branch of the common /PH_HEPEVT/ from
1744 // /PHOEVT/ moves branch back from its CMS system.
1745 //
1746 // Input Parameters: IP: pointer of particle starting branch
1747 // to be given back.
1748 // BOOST: Flag whether boost to CMS was or was
1749 // . not performed.
1750 //
1751 // Output Parameters: Common /PHOEVT/,
1752 //
1753 // Author(s): Z. Was Created at: 24/05/93
1754 // Last Update:
1755 //
1756 //----------------------------------------------------------------------
1757  void PHOOUT(int IP, bool BOOST, int nhep0)
1758  {
1759  int LL, FIRST, LAST, I;
1760  int NN, J, K, NA;
1761  double PB;
1762  double& GAM = phocms.gam;
1763  double* BET = phocms.bet;
1764 
1765  static int i = 1;
1766  if (pho.nhep == pho.nevhep) return;
1767  //-- When parent was not in its rest-frame, boost back...
1768  PHLUPA(10);
1769  if (BOOST) {
1770  //PHOERR(404,"PHOOUT",1.0); // we need to improve this warning: program should never
1771  // enter this place
1772 
1773  double phocms_check = fabs(1 - GAM) + fabs(BET[1 - i]) + fabs(BET[2 - i]) + fabs(BET[3 - i]);
1774  if (phocms_check > 0.001) {
1775  Log::Error() << "Msg. from PHOOUT: possible problems with boosting due to the rounding errors." << endl
1776  << "Boost parameters: (" << GAM << ","
1777  << BET[1 - i] << "," << BET[2 - i] << "," << BET[3 - i] << ")" << endl
1778  << "should be equal to: (1,0,0,0) up to at least several digits." << endl;
1779  } else {
1780  Log::Warning() << "Msg. from PHOOUT: possible problems with boosting due to the rounding errors." << endl
1781  << "Boost parameters: (" << GAM << ","
1782  << BET[1 - i] << "," << BET[2 - i] << "," << BET[3 - i] << ")" << endl
1783  << "should be equal to: (1,0,0,0) up to at least several digits." << endl;
1784  }
1785 
1786  for (J = pho.jdahep[1 - i][1 - i]; J <= pho.jdahep[1 - i][2 - i]; J++) {
1787  PB = -BET[1 - i] * pho.phep[J - i][1 - i] - BET[2 - i] * pho.phep[J - i][2 - i] - BET[3 - i] * pho.phep[J - i][3 - i];
1788  for (K = 1; K <= 3; K++) pho.phep[J - i][K - i] = pho.phep[J - i][K - i] - BET[K - i] * (pho.phep[J - i][4 - i] + PB / (GAM + 1.0));
1789  pho.phep[J - i][4 - i] = GAM * pho.phep[J - i][4 - i] + PB;
1790  }
1791 
1792  //-- ...boost photon, or whatever else has shown up
1793  for (NN = pho.nevhep + 1; NN <= pho.nhep; NN++) {
1794  PB = -BET[1 - i] * pho.phep[NN - i][1 - i] - BET[2 - i] * pho.phep[NN - i][2 - i] - BET[3 - i] * pho.phep[NN - i][3 - i];
1795  for (K = 1; K <= 3;
1796  K++) pho.phep[NN - i][K - i] = pho.phep[NN - i][K - i] - BET[K - i] * (pho.phep[NN - i][4 - i] + PB / (GAM + 1.0));
1797  pho.phep[NN - i][4 - i] = GAM * pho.phep[NN][4 - i] + PB;
1798  }
1799  }
1800  PHCORK(0); // we have to use it because it clears input
1801  // for grandaughters modified in C++
1802  FIRST = hep.jdahep[IP - i][1 - i];
1803  LAST = hep.jdahep[IP - i][2 - i];
1804  // let-s take in original daughters
1805  for (LL = 0; LL <= LAST - FIRST; LL++) {
1806  hep.idhep[FIRST + LL - i] = pho.idhep[3 + LL - i];
1807  for (I = 1; I <= 5; I++) hep.phep[FIRST + LL - i][I - i] = pho.phep[3 + LL - i][I - i];
1808  }
1809 
1810  // let-s take newcomers to the end of HEPEVT.
1811  NA = 3 + LAST - FIRST;
1812  for (LL = 1; LL <= pho.nhep - NA; LL++) {
1813  hep.idhep[nhep0 + LL - i] = pho.idhep[NA + LL - i];
1814  hep.isthep[nhep0 + LL - i] = pho.isthep[NA + LL - i];
1815  hep.jmohep[nhep0 + LL - i][1 - i] = IP;
1816  hep.jmohep[nhep0 + LL - i][2 - i] = hep.jmohep[hep.jdahep[IP - i][1 - i] - i][2 - i];
1817  hep.jdahep[nhep0 + LL - i][1 - i] = 0;
1818  hep.jdahep[nhep0 + LL - i][2 - i] = 0;
1819  for (I = 1; I <= 5; I++) hep.phep[nhep0 + LL - i][I - i] = pho.phep[NA + LL - i][I - i];
1820  }
1821  hep.nhep = hep.nhep + pho.nhep - pho.nevhep;
1822  PHLUPA(20);
1823  return;
1824  }
1825 
1826 //----------------------------------------------------------------------
1827 //
1828 // PHOCHK: checking branch.
1829 //
1830 // Purpose: checks whether particles in the common block /PHOEVT/
1831 // can be served by PHOMAK.
1832 // JFIRST is the position in /PH_HEPEVT/ (!) of the first
1833 // daughter of sub-branch under action.
1834 //
1835 //
1836 // Author(s): Z. Was Created at: 22/10/92
1837 // Last Update: 11/12/00
1838 //
1839 //----------------------------------------------------------------------
1840 // ********************
1841 
1842  void PHOCHK(int JFIRST)
1843  {
1844 
1845  int IDABS, NLAST, I;
1846  bool IFRAD;
1847  int IDENT, K;
1848  static int i = 1, IPPAR = 1;
1849 
1850  NLAST = pho.nhep;
1851  //
1852 
1853  for (I = IPPAR; I <= NLAST; I++) {
1854  IDABS = abs(pho.idhep[I - i]);
1855  // possibly call on PHZODE is a dead (to be omitted) code.
1856  pho.qedrad[I - i] = pho.qedrad[I - i] && F(0, IDABS) && F(0, abs(pho.idhep[1 - i]))
1857  && (pho.idhep[2 - i] == 0);
1858 
1859  if (I > 2) pho.qedrad[I - i] = pho.qedrad[I - i] && hep.qedrad[JFIRST + I - IPPAR - 2 - i];
1860  }
1861 
1862  //--
1863  // now we go to special cases, where pho.qedrad[I) will be overwritten
1864  //--
1865  IDENT = pho.nhep;
1866  if (phokey.iftop) {
1867  // special case of top pair production
1868  for (K = pho.jdahep[1 - i][2 - i]; K >= pho.jdahep[1 - i][1 - i]; K--) {
1869  if (pho.idhep[K - i] != 22) {
1870  IDENT = K;
1871  break; // from loop over K
1872  }
1873  }
1874 
1875  IFRAD = ((pho.idhep[1 - i] == 21) && (pho.idhep[2 - i] == 21))
1876  || ((abs(pho.idhep[1 - i]) <= 6) && (pho.idhep[2 - i] == (-pho.idhep[1 - i])));
1877  IFRAD = IFRAD
1878  && (abs(pho.idhep[3 - i]) == 6) && (pho.idhep[4 - i] == (-pho.idhep[3 - i]))
1879  && (IDENT == 4);
1880  if (IFRAD) {
1881  for (I = IPPAR; I <= NLAST; I++) {
1882  pho.qedrad[I - i] = true;
1883  if (I > 2) pho.qedrad[I - i] = pho.qedrad[I - i] && hep.qedrad[JFIRST + I - IPPAR - 2 - i];
1884  }
1885  }
1886  }
1887  //--
1888  //--
1889  if (phokey.iftop) {
1890  // special case of top decay
1891  for (K = pho.jdahep[1 - i][2 - i]; K >= pho.jdahep[1 - i][1 - i]; K--) {
1892  if (pho.idhep[K - i] != 22) {
1893  IDENT = K;
1894  break;
1895  }
1896  }
1897  IFRAD = ((abs(pho.idhep[1 - i]) == 6) && (pho.idhep[2 - i] == 0));
1898  IFRAD = IFRAD
1899  && (((abs(pho.idhep[3 - i]) == 24) && (abs(pho.idhep[4 - i]) == 5))
1900  || ((abs(pho.idhep[3 - i]) == 5) && (abs(pho.idhep[4 - i]) == 24)))
1901  && (IDENT == 4);
1902 
1903  if (IFRAD) {
1904  for (I = IPPAR; I <= NLAST; I++) {
1905  pho.qedrad[I - i] = true;
1906  if (I > 2) pho.qedrad[I - i] = (pho.qedrad[I - i] && hep.qedrad[JFIRST + I - IPPAR - 2 - i]);
1907  }
1908  }
1909  }
1910  //--
1911  //--
1912  return;
1913  }
1914 
1915 
1916 
1917 //----------------------------------------------------------------------
1918 //
1919 // PHOTOS: PHOton radiation in decays calculation of photon ENErgy
1920 // fraction
1921 //
1922 // Purpose: Subroutine returns photon energy fraction (in (parent
1923 // mass)/2 units) for the decay bremsstrahlung.
1924 //
1925 // Input Parameters: MPASQR: Mass of decaying system squared,
1926 // XPHCUT: Minimum energy fraction of photon,
1927 // XPHMAX: Maximum energy fraction of photon.
1928 //
1929 // Output Parameter: MCHREN: Renormalised mass squared,
1930 // BETA: Beta factor due to renormalisation,
1931 // XPHOTO: Photon energy fraction,
1932 // XF: Correction factor for PHOFA//
1933 //
1934 // Author(s): S. Jadach, Z. Was Created at: 01/01/89
1935 // B. van Eijk, P.Golonka Last Update: 11/07/13
1936 //
1937 //----------------------------------------------------------------------
1938 
1939  void PHOENE(double MPASQR, double* pMCHREN, double* pBETA, double* pBIGLOG, int IDENT)
1940  {
1941  double DATA;
1942  double PRSOFT, PRHARD;
1943  double PRKILL, RRR;
1944  int K, IDME;
1945  double PRSUM;
1946  static int i = 1;
1947  double& MCHREN = *pMCHREN;
1948  double& BETA = *pBETA;
1949  double& BIGLOG = *pBIGLOG;
1950  int& NCHAN = phoexp.nchan;
1951  double& XPHMAX = phophs.xphmax;
1952  double& XPHOTO = phophs.xphoto;
1953  double& MCHSQR = phomom.mchsqr;
1954  double& MNESQR = phomom.mnesqr;
1955 
1956  //--
1957  if (XPHMAX <= phocop.xphcut) {
1958  BETA = PHOFAC(-1); // to zero counter, here beta is dummy
1959  XPHOTO = 0.0;
1960  return;
1961  }
1962  //-- Probabilities for hard and soft bremstrahlung...
1963  MCHREN = 4.0 * MCHSQR / MPASQR / pow(1.0 + MCHSQR / MPASQR, 2);
1964  BETA = sqrt(1.0 - MCHREN);
1965 
1966 #ifdef VARIANTB
1967  // ----------- VARIANT B ------------------
1968  // we replace 1D0/BETA*BIGLOG with (1.0/BETA*BIGLOG+2*phokey.fint)
1969  // for integral of new crude
1970  BIGLOG = log(MPASQR / MCHSQR * (1.0 + BETA) * (1.0 + BETA) / 4.0 *
1971  pow(1.0 + MCHSQR / MPASQR, 2));
1972  PRHARD = phocop.alpha / PI * (1.0 / BETA * BIGLOG + 2 * phokey.fint)
1973  * (log(XPHMAX / phocop.xphcut) - .75 + phocop.xphcut / XPHMAX - .25 * phocop.xphcut * phocop.xphcut / XPHMAX / XPHMAX);
1974  PRHARD = PRHARD * PHOCHA(IDENT) * PHOCHA(IDENT) * phokey.fsec;
1975  // ----------- END OF VARIANT B ------------------
1976 #else
1977  // ----------- VARIANT A ------------------
1978  BIGLOG = log(MPASQR / MCHSQR * (1.0 + BETA) * (1.0 + BETA) / 4.0 *
1979  pow(1.0 + MCHSQR / MPASQR, 2));
1980  PRHARD = phocop.alpha / PI * (1.0 / BETA * BIGLOG) *
1981  (log(XPHMAX / phocop.xphcut) - .75 + phocop.xphcut / XPHMAX - .25 * phocop.xphcut * phocop.xphcut / XPHMAX / XPHMAX);
1982  PRHARD = PRHARD * PHOCHA(IDENT) * PHOCHA(IDENT) * phokey.fsec * phokey.fint;
1983  //me_channel_(&IDME);
1984  IDME = HEPEVT_struct::ME_channel;
1985  // write(*,*) 'KANALIK IDME=',IDME
1986  if (IDME == 0) {
1987  // do nothing
1988  }
1989 
1990  else if (IDME == 1) {
1991  PRHARD = PRHARD / (1.0 + 0.75 * phocop.alpha / PI); // NLO
1992  } else if (IDME == 2) {
1993  // work on virtual crrections in W decay to be done.
1994  } else {
1995  cout << "problem with ME_CHANNEL IDME= " << IDME << endl;
1996  exit(-1);
1997  }
1998 
1999  //----------- END OF VARIANT A ------------------
2000 #endif
2001  if (phopro.irep == 0) phopro.probh = 0.0;
2002  PRKILL = 0.0;
2003  if (phokey.iexp) { // IEXP
2004  NCHAN = NCHAN + 1;
2005  if (phoexp.expini) { // EXPINI
2006  phoexp.pro[NCHAN - i] = PRHARD + 0.25 * (1.0 + phokey.fint); // we store hard photon emission prob
2007  //for leg NCHAN
2008  PRHARD = 0.0; // to kill emission at initialization call
2009  phopro.probh = PRHARD;
2010  } else { // EXPINI
2011  PRSUM = 0.0;
2012  for (K = NCHAN; K <= phoexp.NX; K++) PRSUM = PRSUM + phoexp.pro[K - i];
2013  PRHARD = PRHARD / PRSUM; // note that PRHARD may be smaller than
2014  //phoexp.pro[NCHAN) because it is calculated
2015  // for kinematical configuartion as is
2016  // (with effects of previous photons)
2017  PRKILL = phoexp.pro[NCHAN - i] / PRSUM - PRHARD;
2018 
2019  } // EXPINI
2020  PRSOFT = 1.0 - PRHARD;
2021  } else { // IEXP
2022  PRHARD = PRHARD * PHOFAC(0); // PHOFAC is used to control eikonal
2023  // formfactors for non exp version only
2024  // here PHOFAC(0)=1 at least now.
2025  phopro.probh = PRHARD;
2026  } // IEXP
2027  PRSOFT = 1.0 - PRHARD;
2028  //--
2029  //-- Check on kinematical bounds
2030  if (phokey.iexp) {
2031  if (PRSOFT < -5.0E-8) {
2032  DATA = PRSOFT;
2033  PHOERR(2, "PHOENE", DATA);
2034  }
2035  } else {
2036  if (PRSOFT < 0.1) {
2037  DATA = PRSOFT;
2038  PHOERR(2, "PHOENE", DATA);
2039  }
2040  }
2041 
2042  RRR = Photos::randomDouble();
2043  if (RRR < PRSOFT) {
2044  //--
2045  //-- No photon... (ie. photon too soft)
2046  XPHOTO = 0.0;
2047  if (RRR < PRKILL) XPHOTO = -5.0; //No photon...no further trials
2048  } else {
2049  //--
2050  //-- Hard photon... (ie. photon hard enough).
2051  //-- Calculate Altarelli-Parisi Kernel
2052  do {
2053  XPHOTO = exp(Photos::randomDouble() * log(phocop.xphcut / XPHMAX));
2054  XPHOTO = XPHOTO * XPHMAX;
2055  } while (Photos::randomDouble() > ((1.0 + pow(1.0 - XPHOTO / XPHMAX, 2)) / 2.0));
2056  }
2057 
2058  //--
2059  //-- Calculate parameter for PHOFAC function
2060  phopro.xf = 4.0 * MCHSQR * MPASQR / pow(MPASQR + MCHSQR - MNESQR, 2);
2061  return;
2062  }
2063 
2064 
2065 //----------------------------------------------------------------------
2066 //
2067 // PHOTOS: Photon radiation in decays
2068 //
2069 // Purpose: Order (alpha) radiative corrections are generated in
2070 // the decay of the IPPAR-th particle in the HEP-like
2071 // common /PHOEVT/. Photon radiation takes place from one
2072 // of the charged daughters of the decaying particle IPPAR
2073 // WT is calculated, eventual rejection will be performed
2074 // later after inclusion of interference weight.
2075 //
2076 // Input Parameter: IPPAR: Pointer to decaying particle in
2077 // /PHOEVT/ and the common itself,
2078 //
2079 // Output Parameters: Common /PHOEVT/, either with or without a
2080 // photon(s) added.
2081 // WT weight of the configuration
2082 //
2083 // Author(s): Z. Was, B. van Eijk Created at: 26/11/89
2084 // Last Update: 12/07/13
2085 //
2086 //----------------------------------------------------------------------
2087 
2088  void PHOPRE(int IPARR, double* pWT, int* pNEUDAU, int* pNCHARB)
2089  {
2090  int CHAPOI[pho.nmxhep];
2091  double MINMAS, MPASQR, MCHREN;
2092  double EPS, DEL1, DEL2, DATA, BIGLOG;
2093  double MASSUM;
2094  int IP, IPPAR, I, J, ME, NCHARG, NEUPOI, NLAST;
2095  int IDABS;
2096  double WGT;
2097  int IDME;
2098  double a, b;
2099  double& WT = *pWT;
2100  int& NEUDAU = *pNEUDAU;
2101  int& NCHARB = *pNCHARB;
2102  double& COSTHG = phophs.costhg;
2103  double& SINTHG = phophs.sinthg;
2104  double& XPHOTO = phophs.xphoto;
2105  double& XPHMAX = phophs.xphmax;
2106  double* PNEUTR = phomom.pneutr;
2107  double& MCHSQR = phomom.mchsqr;
2108  double& MNESQR = phomom.mnesqr;
2109 
2110  static int i = 1;
2111 
2112  //--
2113  IPPAR = IPARR;
2114  //-- Store pointers for cascade treatement...
2115  IP = IPPAR;
2116  NLAST = pho.nhep;
2117 
2118  //--
2119  //-- Check decay multiplicity..
2120  if (pho.jdahep[IP - i][1 - i] == 0) return;
2121 
2122  //--
2123  //-- Loop over daughters, determine charge multiplicity
2124 
2125  NCHARG = 0;
2126  phopro.irep = 0;
2127  MINMAS = 0.0;
2128  MASSUM = 0.0;
2129  for (I = pho.jdahep[IP - i][1 - i]; I <= pho.jdahep[IP - i][2 - i]; I++) {
2130  //--
2131  //--
2132  //-- Exclude marked particles, quarks and gluons etc...
2133  IDABS = abs(pho.idhep[I - i]);
2134  if (pho.qedrad[I - pho.jdahep[IP - i][1 - i] + 3 - i]) {
2135  if (PHOCHA(pho.idhep[I - i]) != 0) {
2136  NCHARG = NCHARG + 1;
2137  if (NCHARG > pho.nmxhep) {
2138  DATA = NCHARG;
2139  PHOERR(1, "PHOTOS", DATA);
2140  }
2141  CHAPOI[NCHARG - i] = I;
2142  }
2143  MINMAS = MINMAS + pho.phep[I - i][5 - i] * pho.phep[I - i][5 - i];
2144  }
2145  MASSUM = MASSUM + pho.phep[I - i][5 - i];
2146  }
2147 
2148  if (NCHARG != 0) {
2149  //--
2150  //-- Check that sum of daughter masses does not exceed parent mass
2151  if ((pho.phep[IP - i][5 - i] - MASSUM) / pho.phep[IP - i][5 - i] > 2.0 * phocop.xphcut) {
2152  //--
2153 label30:
2154 
2155 // do{
2156 
2157  for (J = 1; J <= 3; J++) PNEUTR[J - i] = -pho.phep[CHAPOI[NCHARG - i] - i][J - i];
2158  PNEUTR[4 - i] = pho.phep[IP - i][5 - i] - pho.phep[CHAPOI[NCHARG - i] - i][4 - i];
2159  //--
2160  //-- Calculate invariant mass of 'neutral' etc. systems
2161  MPASQR = pho.phep[IP - i][5 - i] * pho.phep[IP - i][5 - i];
2162  MCHSQR = pow(pho.phep[CHAPOI[NCHARG - i] - i][5 - i], 2);
2163  if ((pho.jdahep[IP - i][2 - i] - pho.jdahep[IP - i][1 - i]) == 1) {
2164  NEUPOI = pho.jdahep[IP - i][1 - i];
2165  if (NEUPOI == CHAPOI[NCHARG - i]) NEUPOI = pho.jdahep[IP - i][2 - i];
2166  MNESQR = pho.phep[NEUPOI - i][5 - i] * pho.phep[NEUPOI - i][5 - i];
2167  PNEUTR[5 - i] = pho.phep[NEUPOI - i][5 - i];
2168  } else {
2169  MNESQR = pow(PNEUTR[4 - i], 2) - pow(PNEUTR[1 - i], 2) - pow(PNEUTR[2 - i], 2) - pow(PNEUTR[3 - i], 2);
2170  MNESQR = max(MNESQR, MINMAS - MCHSQR);
2171  PNEUTR[5 - i] = sqrt(MNESQR);
2172  }
2173 
2174  //--
2175  //-- Determine kinematical limit...
2176  XPHMAX = (MPASQR - pow(PNEUTR[5 - i] + pho.phep[CHAPOI[NCHARG - i] - i][5 - i], 2)) / MPASQR;
2177 
2178  //--
2179  //-- Photon energy fraction...
2180  PHOENE(MPASQR, &MCHREN, &phwt.beta, &BIGLOG, pho.idhep[CHAPOI[NCHARG - i] - i]);
2181  //--
2182 
2183  if (XPHOTO < -4.0) {
2184  NCHARG = 0; // we really stop trials
2185  XPHOTO = 0.0; // in this case !!
2186  //-- Energy fraction not too large (very seldom) ? Define angle.
2187  } else if ((XPHOTO < phocop.xphcut) || (XPHOTO > XPHMAX)) {
2188  //--
2189  //-- No radiation was accepted, check for more daughters that may ra-
2190  //-- diate and correct radiation probability...
2191  NCHARG = NCHARG - 1;
2192  if (NCHARG > 0) phopro.irep = phopro.irep + 1;
2193  if (NCHARG > 0) goto label30;
2194  } else {
2195  //--
2196  //-- Angle is generated in the frame defined by charged vector and
2197  //-- PNEUTR, distribution is taken in the infrared limit...
2198  EPS = MCHREN / (1.0 + phwt.beta);
2199  //--
2200  //-- Calculate sin(theta) and cos(theta) from interval variables
2201  DEL1 = (2.0 - EPS) * pow(EPS / (2.0 - EPS), Photos::randomDouble());
2202  DEL2 = 2.0 - DEL1;
2203 
2204 #ifdef VARIANTB
2205  // ----------- VARIANT B ------------------
2206  // corrections for more efiicient interference correction,
2207  // instead of doubling crude distribution, we add flat parallel channel
2208  if (Photos::randomDouble() < BIGLOG / phwt.beta / (BIGLOG / phwt.beta + 2 * phokey.fint)) {
2209  COSTHG = (1.0 - DEL1) / phwt.beta;
2210  SINTHG = sqrt(DEL1 * DEL2 - MCHREN) / phwt.beta;
2211  } else {
2212  COSTHG = -1.0 + 2 * Photos::randomDouble();
2213  SINTHG = sqrt(1.0 - COSTHG * COSTHG);
2214  }
2215 
2216  if (phokey.fint > 1.0) {
2217 
2218  WGT = 1.0 / (1.0 - phwt.beta * COSTHG);
2219  WGT = WGT / (WGT + phokey.fint);
2220  // WGT=1.0 // ??
2221  } else {
2222  WGT = 1.0;
2223  }
2224  //
2225  // ----------- END OF VARIANT B ------------------
2226 #else
2227  // ----------- VARIANT A ------------------
2228  COSTHG = (1.0 - DEL1) / phwt.beta;
2229  SINTHG = sqrt(DEL1 * DEL2 - MCHREN) / phwt.beta;
2230  WGT = 1.0;
2231  // ----------- END OF VARIANT A ------------------
2232 #endif
2233  //--
2234  //-- Determine spin of particle and construct code for matrix element
2235  ME = (int)(2.0 * PHOSPI(pho.idhep[CHAPOI[NCHARG - i] - i]) + 1.0);
2236  //--
2237  //-- Weighting procedure with 'exact' matrix element, reconstruct kine-
2238  //-- matics for photon, neutral and charged system and update /PHOEVT/.
2239  //-- Find pointer to the first component of 'neutral' system
2240  for (I = pho.jdahep[IP - i][1 - i]; I <= pho.jdahep[IP - i][2 - i]; I++) {
2241  if (I != CHAPOI[NCHARG - i]) {
2242  NEUDAU = I;
2243  goto label51; //break; // to 51
2244  }
2245  }
2246  //--
2247  //-- Pointer not found...
2248  DATA = NCHARG;
2249  PHOERR(5, "PHOKIN", DATA);
2250 label51:
2251 
2252  NCHARB = CHAPOI[NCHARG - i];
2253  NCHARB = NCHARB - pho.jdahep[IP - i][1 - i] + 3;
2254  NEUDAU = NEUDAU - pho.jdahep[IP - i][1 - i] + 3;
2255 
2256  IDME = HEPEVT_struct::ME_channel;
2257  // two options introduced temporarily.
2258  // In future always PHOCOR-->PHOCORN
2259  // Tests and adjustment of wts for Znlo needed.
2260  // otherwise simple change. PHOCORN implements
2261  // exact ME for scalar to 2 scalar decays.
2262  if (IDME == 2) {
2263  b = PHOCORN(MPASQR, MCHREN, ME);
2264  WT = b * WGT;
2265  WT = WT / (1 - XPHOTO / XPHMAX + 0.5 * pow(XPHOTO / XPHMAX, 2)) * (1 - XPHOTO / XPHMAX) / 2; // factor to go to WnloWT
2266  } else if (IDME == 1) {
2267 
2268  a = PHOCOR(MPASQR, MCHREN, ME);
2269  b = PHOCORN(MPASQR, MCHREN, ME);
2270  WT = b * WGT ;
2271  WT = WT * phwt.wt1 * phwt.wt2 * phwt.wt3 / phocorwt.phocorwt1 / phocorwt.phocorwt2 / phocorwt.phocorwt3; // factor to go to ZnloWT
2272  // write(*,*) ' -----------'
2273  // write(*,*) phwt.wt1,' ',phwt.wt2,' ',phwt.wt3
2274  // write(*,*) phocorwt.phocorwt1,' ',phocorwt.phocorwt2,' ',phocorwt.phocorwt3
2275  } else {
2276  a = PHOCOR(MPASQR, MCHREN, ME);
2277  WT = a * WGT;
2278 // WT=b*WGT; // /(1-XPHOTO/XPHMAX+0.5*pow(XPHOTO/XPHMAX,2))*(1-XPHOTO/XPHMAX)/2;
2279  }
2280 
2281 
2282 
2283  }
2284  } else {
2285  DATA = pho.phep[IP - i][5 - i] - MASSUM;
2286  PHOERR(10, "PHOTOS", DATA);
2287  }
2288  }
2289 
2290  //--
2291  return;
2292  }
2293 
2294 
2295 //----------------------------------------------------------------------
2296 //
2297 // PHOMAK: PHOtos MAKe
2298 //
2299 // Purpose: Single or double bremstrahlung radiative corrections
2300 // are generated in the decay of the IPPAR-th particle in
2301 // the HEP common /PH_HEPEVT/. Example of the use of
2302 // general tools.
2303 //
2304 // Input Parameter: IPPAR: Pointer to decaying particle in
2305 // /PH_HEPEVT/ and the common itself
2306 //
2307 // Output Parameters: Common /PH_HEPEVT/, either with or without
2308 // particles added.
2309 //
2310 // Author(s): Z. Was, Created at: 26/05/93
2311 // Last Update: 29/01/05
2312 //
2313 //----------------------------------------------------------------------
2314 
2315  void PHOMAK(int IPPAR, int NHEP0)
2316  {
2317 
2318  double DATA;
2319  int IP, NCHARG, IDME;
2320  int IDUM;
2321  int NCHARB, NEUDAU;
2322  double RN, WT;
2323  bool BOOST;
2324  static int i = 1;
2325  //--
2326  IP = IPPAR;
2327  IDUM = 1;
2328  NCHARG = 0;
2329  //--
2330  PHOIN(IP, &BOOST, &NHEP0);
2331  PHOCHK(hep.jdahep[IP - i][1 - i]);
2332  WT = 0.0;
2333  PHOPRE(1, &WT, &NEUDAU, &NCHARB);
2334 
2335  if (WT == 0.0) return;
2336  RN = Photos::randomDouble();
2337  // PHODO is caling randomDouble(), thus change of series if it is moved before if
2338  PHODO(1, NCHARB, NEUDAU);
2339 
2340 #ifdef VARIANTB
2341  // we eliminate divisions /phokey.fint in variant B. ???
2342 #endif
2343  // get ID of channel dependent ME, ID=0 means no
2344 
2345  IDME = HEPEVT_struct::ME_channel;
2346  // corrections for matrix elements
2347  // controlled by IDME
2348  // write(*,*) 'KANALIK IDME=',IDME
2349 
2350  if (IDME == 0) { // default
2351 
2352  if (phokey.interf) WT = WT * PHINT(IDUM);
2353  if (phokey.ifw) PHOBW(&WT); // extra weight for leptonic W decay
2354  } else if (IDME == 2) { // ME weight for leptonic W decay
2355 
2356  PhotosMEforW::PHOBWnlo(&WT);
2357  WT = WT * 2.0;
2358  } else if (IDME == 1) { // ME weight for leptonic Z decay
2359 
2360  WT = WT * PhotosMEforZ::phwtnlo();
2361  } else {
2362  cout << "problem with ME_CHANNEL IDME= " << IDME << endl;
2363  exit(-1);
2364  }
2365 
2366 #ifndef VARIANTB
2367  WT = WT / phokey.fint; // FINT must be in variant A
2368 #endif
2369 
2370  DATA = WT;
2371  if (WT > 1.0) PHOERR(3, "WT_INT", DATA);
2372  // weighting
2373  if (RN <= WT) {
2374  PHOOUT(IP, BOOST, NHEP0);
2375  }
2376  return;
2377  }
2378 
2379 //----------------------------------------------------------------------
2380 //
2381 // PHTYPE: Central manadgement routine.
2382 //
2383 // Purpose: defines what kind of the
2384 // actions will be performed at point ID.
2385 //
2386 // Input Parameters: ID: pointer of particle starting branch
2387 // in /PH_HEPEVT/ to be treated.
2388 //
2389 // Output Parameters: Common /PH_HEPEVT/.
2390 //
2391 // Author(s): Z. Was Created at: 24/05/93
2392 // P. Golonka Last Update: 27/06/04
2393 //
2394 //----------------------------------------------------------------------
2395  void PHTYPE(int ID)
2396  {
2397 
2398  int K;
2399  double PRSUM, ESU;
2400  int NHEP0;
2401  bool IPAIR;
2402  bool IPHOT;
2403  double RN, SUM;
2404  bool IFOUR;
2405  int& NCHAN = phoexp.nchan;
2406 
2407  static int i = 1;
2408 
2409 
2410  //--
2411  IFOUR = phokey.itre; // we can make internal choice whether
2412  // we want 3 or four photons at most.
2413  IPAIR = false;
2414  IPAIR = Photos::IfPair;
2415  IPHOT = true;
2416  IPHOT = Photos::IfPhot;
2417  //-- Check decay multiplicity..
2418  if (hep.jdahep[ID - i][1 - i] == 0) return;
2419  // if (hep.jdahep[ID-i][1-i]==hep.jdahep[ID-i][2-i]) return;
2420  //--
2421  NHEP0 = hep.nhep;
2422 
2423  // initialization of pho.qedrad for new event.
2424  // some of (old style and doubling) restrictions introduced with PHOCHK,
2425  // also new pairs have emissions blocked with pho.qedrad[]
2426  // most of the restrictions are introduced prior decay vertex is copied
2427  // to struct pho.
2428 
2429  // Establish size for the struct pho: number of daughters + 2 places for mothers (no grandmothers)
2430  // This solution is `hep.nhep hardy'. Use of hep.nhep was perilous
2431  // if decaying particle (ID-i) was the first in the event. That was the case of EvtGen
2432  // interface. We adopt to such non-standard HepMC fill.
2433  // NOTE: here 'max' is used as a safety for future changes to hep or pho content.
2434  // TP ZW (26.09.15): Thanks to Michal Kreps and John Back
2435 
2436  int pho_size = max(NHEP0, (hep.jdahep[ID - i][2 - i] - hep.jdahep[ID - i][1 - i] + 1) + 2);
2437 
2438  for (int I = 0; I < pho_size; ++I) {
2439  pho.qedrad[I] = true;
2440  }
2441 
2442  double elMass = 0.000511;
2443  double muMass = 0.1057;
2444  double STRENG = 0.5;
2445 
2446  if (IPAIR) {
2447 
2448  switch (Photos::momentumUnit) {
2449  case Photos::GEV:
2450  elMass = 0.000511;
2451  muMass = 0.1057;
2452  break;
2453  case Photos::MEV:
2454  elMass = 0.511;
2455  muMass = 105.7;
2456  break;
2457  default:
2458  Log::Error() << "GEV or MEV unit must be set for pair emission" << endl;
2459  break;
2460  };
2461  PHOPAR(ID, NHEP0, 11, elMass, &STRENG);
2462  PHOPAR(ID, NHEP0, 13, muMass, &STRENG);
2463  }
2464  //--
2465  if (IPHOT) {
2466  if (phokey.iexp) {
2467  phoexp.expini = true; // Initialization/cleaning
2468  for (NCHAN = 1; NCHAN <= phoexp.NX; NCHAN++)
2469  phoexp.pro[NCHAN - i] = 0.0;
2470  NCHAN = 0;
2471 
2472  phokey.fsec = 1.0;
2473  PHOMAK(ID, NHEP0); // Initialization/crude formfactors into
2474  // phoexp.pro[NCHAN)
2475  phoexp.expini = false;
2476  RN = Photos::randomDouble();
2477  PRSUM = 0.0;
2478  for (K = 1; K <= phoexp.NX; K++)PRSUM = PRSUM + phoexp.pro[K - i];
2479 
2480  ESU = exp(-PRSUM);
2481  // exponent for crude Poissonian multiplicity
2482  // distribution, will be later overwritten
2483  // to give probability for k
2484  SUM = ESU;
2485  // distribuant for the crude Poissonian
2486  // at first for k=0
2487  for (K = 1; K <= 100; K++) { // hard coded max (photon) multiplicity is 100
2488  if (RN < SUM) break;
2489  ESU = ESU * PRSUM / K; // we get at K ESU=EXP(-PRSUM)*PRSUM**K/K!
2490  SUM = SUM + ESU; // thus we get distribuant at K.
2491  NCHAN = 0;
2492  PHOMAK(ID, NHEP0); // LOOPING
2493  if (SUM > 1.0 - phokey.expeps) break;
2494  }
2495 
2496  } else if (IFOUR) {
2497  //-- quatro photon emission
2498  phokey.fsec = 1.0;
2499  RN = Photos::randomDouble();
2500  if (RN >= 23.0 / 24.0) {
2501  PHOMAK(ID, NHEP0);
2502  PHOMAK(ID, NHEP0);
2503  PHOMAK(ID, NHEP0);
2504  PHOMAK(ID, NHEP0);
2505  } else if (RN >= 17.0 / 24.0) {
2506  PHOMAK(ID, NHEP0);
2507  PHOMAK(ID, NHEP0);
2508  } else if (RN >= 9.0 / 24.0) {
2509  PHOMAK(ID, NHEP0);
2510  } else {
2511  }
2512  } else if (phokey.itre) {
2513  //-- triple photon emission
2514  phokey.fsec = 1.0;
2515  RN = Photos::randomDouble();
2516  if (RN >= 5.0 / 6.0) {
2517  PHOMAK(ID, NHEP0);
2518  PHOMAK(ID, NHEP0);
2519  PHOMAK(ID, NHEP0);
2520  } else if (RN >= 2.0 / 6.0) {
2521  PHOMAK(ID, NHEP0);
2522  }
2523  } else if (phokey.isec) {
2524  //-- double photon emission
2525  phokey.fsec = 1.0;
2526  RN = Photos::randomDouble();
2527  if (RN >= 0.5) {
2528  PHOMAK(ID, NHEP0);
2529  PHOMAK(ID, NHEP0);
2530  }
2531  } else {
2532  //-- single photon emission
2533  phokey.fsec = 1.0;
2534  PHOMAK(ID, NHEP0);
2535  }
2536  }
2537  //--
2538  //-- lepton anti-lepton pair(s)
2539  // we prepare to migrate half of tries to before photons accordingly to LL
2540  // pho.qedrad is not yet used by PHOPAR
2541  if (IPAIR) {
2542  PHOPAR(ID, NHEP0, 11, elMass, &STRENG);
2543  PHOPAR(ID, NHEP0, 13, muMass, &STRENG);
2544  }
2545 
2546 // Fill Photos::EventNo in user main program to provide
2547 // debug input for specific events, e.g.:
2548 // if(Photos::EventNo==1331094) printf("PHOTOS: event no: %10i finished\n",Photos::EventNo);
2549  }
2550 
2551  /*----------------------------------------------------------------------
2552 
2553  PHOTOS: Photon radiation in decays
2554 
2555  Purpose: e+e- pairs are generated in
2556  the decay of the IPPAR-th particle in the HEP-like
2557  common /PHOEVT/. Radiation takes place from one
2558  of the charged daughters of the decaying particle IPPAR
2559 
2560 
2561 
2562  Input Parameter: IPPAR: Pointer to decaying particle in
2563  /PHOEVT/ and the common itself,
2564  NHEP0 length of the /HEPEVT/ entry
2565  before starting any activity on this
2566  IPPAR decay.
2567  Output Parameters: Common /HEPEVT/, either with or without a
2568  e+e-(s) added.
2569 
2570 
2571  Author(s): Z. Was, Created at: 01/06/93
2572  Last Update:
2573 
2574  ----------------------------------------------------------------------*/
2575  void PHOPAR(int IPARR, int NHEP0, int idlep, double masslep, double* pSTRENG)
2576  {
2577  double PCHAR[4], PNEU[4], PELE[4], PPOZ[4], BUF[4];
2578  int IP, IPPAR, NLAST;
2579  bool BOOST, JESLI;
2580  static int i = 1;
2581  IPPAR = IPARR;
2582 
2583  double& STRENG = *pSTRENG;
2584 
2585  // Store pointers for cascade treatment...
2586  IP = 0;
2587  NLAST = pho.nhep;
2588  // Check decay multiplicity..
2589  PHOIN(IPPAR, &BOOST, &NHEP0);
2590  PHOCHK(pho.jdahep[IP][0]); // should be loop over all mothers?
2591  PHLUPA(100);
2592 
2593  if (pho.jdahep[IP][0] == 0) return;
2594  if (pho.jdahep[IP][0] == pho.jdahep[IP][1]) return;
2595 
2596  // Loop over charged daughters
2597  for (int I = pho.jdahep[IP][0] - i; I <= pho.jdahep[IP][1] - i; ++I) {
2598 
2599 
2600  // Skip this particle if it has no charge
2601  if (PHOCHA(pho.idhep[I]) == 0) continue;
2602 
2603  int IDABS = abs(pho.idhep[I]);
2604  // at the moment the following re-definition make not much sense as constraints
2605  // were already checked before for photons tries.
2606  // we have to come back to this when we will have pairs emitted before photons.
2607 
2608  pho.qedrad[I] = pho.qedrad[I] && F(0, IDABS) && F(0, abs(pho.idhep[1 - i]))
2609  && (pho.idhep[2 - i] == 0);
2610 
2611  if (!pho.qedrad[I]) continue; //
2612 
2613 
2614  // Set 3-vectors
2615  for (int J = 0; J < 3; ++J) {
2616  PCHAR[J] = pho.phep[I][J];
2617  PNEU [J] = -pho.phep[I][J];
2618  }
2619 
2620  // Set energy
2621  PNEU[3] = pho.phep[IP][3] - pho.phep[I][3];
2622  PCHAR[3] = pho.phep[I][3];
2623  // Set mass
2624  double AMCH = pho.phep[I][4];
2625  //here we attempt generating pair from PCHAR. One of the charged
2626  //decay products; that is why algorithm works in a loop.
2627  //PNEU is four vector of all decay products except PCHAR
2628  //we do not care on rare cases when two pairs could be generated
2629  //we assume it is negligibly rare and fourth order in alpha anyway
2630  //TRYPAR should take as an input electron mass.
2631  //then it can be used for muons.
2632  // printf ("wrotki %10.7f\n",STRENG);
2633  /*
2634  double PCH[4]={0};
2635  double PNEu[4]={0};
2636  double CC1;
2637  double CC2;
2638 
2639  for(int K = 0; K<4; ++K) {
2640  PCH[K]=PCHAR[K];
2641  PNEu[K]=PNEU[K];
2642  }
2643  */
2644  // printf ("idlep,pdgidid= %10i %10i\n",idlep,pho.idhep[I]);
2645 
2646  // arrangements for the case when emitted lept6ons have
2647  // the same flavour as emitters
2648  bool sameflav = abs(idlep) == abs(pho.idhep[I]);
2649  int idsign = 1;
2650  if (pho.idhep[I] < 0) idsign = -1; // this is to ensure
2651  //that first lepton has the same PDGID as emitter
2652 
2653  trypar(&JESLI, &STRENG, AMCH, masslep, PCHAR, PNEU, PELE, PPOZ, &sameflav);
2654  // printf ("rowerek %10.7f\n",STRENG);
2655  //emitted pair four momenta are stored in PELE PPOZ
2656  //then JESLI=.true.
2657  /*
2658  if (JESLI) {
2659  // printf ("PCHAR %10.7f %10.7f %10.7f %10.7f\n",PCHAR[0],PCHAR[1],PCHAR[2],PCHAR[3]);
2660  //printf ("PNEU %10.7f %10.7f %10.7f %10.7f\n",PNEU[0],PNEU[1],PNEU[2],PNEU[3]);
2661 
2662  // printf ("PNEu %10.7f %10.7f %10.7f %10.7f\n",PNEu[0],PNEu[1],PNEu[2],PNEu[3]);
2663 
2664  printf ("PELE %10.7f %10.7f %10.7f %10.7f\n",PELE[0],PELE[1],PELE[2],PELE[3]);
2665  printf ("PPOZ %10.7f %10.7f %10.7f %10.7f\n",PPOZ[0],PPOZ[1],PPOZ[2],PPOZ[3]);
2666  printf ("-----------------\n");
2667  printf ("PCH %10.7f %10.7f %10.7f %10.7f\n",PCH[0],PCH[1],PCH[2],PCH[3]);
2668  CC1=(PELE[0]*PCHAR[0]+PELE[1]*PCHAR[1]+PELE[2]*PCHAR[2])/sqrt(PELE[0]*PELE[0]+PELE[1]*PELE[1]+PELE[2]*PELE[2])/sqrt(PCHAR[0]*PCHAR[0]+PCHAR[1]*PCHAR[1]+PCHAR[2]*PCHAR[2]);
2669  CC2=(PPOZ[0]*PCHAR[0]+PPOZ[1]*PCHAR[1]+PPOZ[2]*PCHAR[2])/sqrt(PPOZ[0]*PPOZ[0]+PPOZ[1]*PPOZ[1]+PPOZ[2]*PPOZ[2])/sqrt(PCHAR[0]*PCHAR[0]+PCHAR[1]*PCHAR[1]+PCHAR[2]*PCHAR[2]);
2670 
2671  printf ("-=================-\n");
2672 
2673  }
2674  */
2675  // If JESLI = true, we modify old particles of the vertex
2676  if (JESLI) {
2677  PHLUPA(1010);
2678 
2679  // we have to correct 4-momenta
2680  // of all decay products
2681  // we use PARTRA for that
2682  // PELE PPOZ are in right frame
2683  for (int J = pho.jdahep[IP][0] - i; J <= pho.jdahep[IP][1] - i; ++J) {
2684  for (int K = 0; K < 4; ++K) {
2685  BUF[K] = pho.phep[J][K];
2686  }
2687  if (J == I) partra(1, BUF);
2688  else partra(-1, BUF);
2689 
2690  for (int K = 0; K < 4; ++K) {
2691  pho.phep[J][K] = BUF[K];
2692  }
2693  /*
2694  if (J == I){
2695  printf ("PCHar %10.7f %10.7f %10.7f %10.7f\n",pho.phep[J][0],pho.phep[J][1],pho.phep[J][2],pho.phep[J][3]);
2696  printf ("c1= %10.7f\n",CC1);
2697  printf ("c2= %10.7f\n",CC2);
2698  printf ("-=#####################################====-\n");
2699  }
2700  */
2701  }
2702 
2703  PHLUPA(1011);
2704 
2705  if (darkr.ifspecial == 1) {
2706  // virtual adding to vertex
2707  pho.nhep = pho.nhep + 1;
2708  pho.isthep[pho.nhep - i] = 2;
2709  pho.idhep [pho.nhep - i] = darkr.IDspecial;
2710  pho.jmohep[pho.nhep - i][0] = IP;
2711  pho.jmohep[pho.nhep - i][1] = 0;
2712  pho.jdahep[pho.nhep - i][0] = 0;
2713  pho.jdahep[pho.nhep - i][1] = 0;
2714  pho.qedrad[pho.nhep - i] = false;
2715 
2716  pho.phep[pho.nhep - i][4] = sqrt(-(PELE[0] + PPOZ[0]) * (PELE[0] + PPOZ[0])
2717  - (PELE[1] + PPOZ[1]) * (PELE[1] + PPOZ[1])
2718  - (PELE[2] + PPOZ[2]) * (PELE[2] + PPOZ[2])
2719  + (PELE[3] + PPOZ[3]) * (PELE[3] + PPOZ[3]));
2720  }
2721  //double life=darkr.SpecialLife;
2722  double life = darkr.SpecialLife * (-log(Photos::randomDouble()));
2723 
2724  // here was missing if(darkr.ifspecial==1) zbw:16.09.2021
2725  if (darkr.ifspecial == 1) {
2726  for (int K = 1; K <= 4; ++K) {
2727  pho.phep[pho.nhep - i][K - i] = PELE[K - i] + PPOZ[K - i];
2728  pho.vhep[pho.nhep - i][K - i] = pho.vhep[IP][K - i]
2729  + (PELE[K - i] + PPOZ[K - i]) / pho.phep[pho.nhep - i][4] * life;
2730  }
2731  }
2732 
2733  // electron: adding to vertex
2734  pho.nhep = pho.nhep + 1;
2735  pho.isthep[pho.nhep - i] = 1;
2736  pho.idhep [pho.nhep - i] = idlep * idsign;
2737  pho.jmohep[pho.nhep - i][0] = IP;
2738  pho.jmohep[pho.nhep - i][1] = 0;
2739  pho.jdahep[pho.nhep - i][0] = 0;
2740  pho.jdahep[pho.nhep - i][1] = 0;
2741  pho.qedrad[pho.nhep - i] = false;
2742 
2743 
2744  for (int K = 1; K <= 4; ++K) {
2745  pho.phep[pho.nhep - i][K - i] = PELE[K - i];
2746  if (darkr.ifspecial == 1)
2747  pho.vhep[pho.nhep - i][K - i] = pho.vhep[pho.nhep - i - 1][K - i];
2748  }
2749 
2750  pho.phep[pho.nhep - i][4] = masslep;
2751 
2752  // positron: adding
2753  pho.nhep = pho.nhep + 1;
2754  pho.isthep[pho.nhep - i] = 1;
2755  pho.idhep [pho.nhep - i] = -idlep * idsign;
2756  pho.jmohep[pho.nhep - i][0] = IP;
2757  pho.jmohep[pho.nhep - i][1] = 0;
2758  pho.jdahep[pho.nhep - i][0] = 0;
2759  pho.jdahep[pho.nhep - i][1] = 0;
2760  pho.qedrad[pho.nhep - i] = false;
2761 
2762  for (int K = 1; K <= 4; ++K) {
2763  pho.phep[pho.nhep - i][K - i] = PPOZ[K - i];
2764  if (darkr.ifspecial == 1)
2765  pho.vhep[pho.nhep - i][K - i] = pho.vhep[pho.nhep - i - 1][K - i];
2766  }
2767 
2768  // for mc-test with KORALW, mumu from mu mu emissions: BEGIN
2769  /*
2770  double RRX[2];
2771  for( int k=0;k<=1;k++) RRX[k]=Photos::randomDouble();
2772 
2773  for(int KK=0;KK<=pho.nhep-i;KK++){
2774  for(int KJ=KK+1;KJ<=pho.nhep-i;KJ++){
2775  // 1 <-> 3
2776  if(RRX[0]>.5&&pho.idhep[KK]==13&&pho.idhep[KJ]==13){
2777  for( int k=0;k<=3;k++){
2778  double stored=pho.phep[KK][k];
2779  pho.phep[KK][k]=pho.phep[KJ][k];
2780  pho.phep[KJ][k]=stored;
2781  }
2782  }
2783  // 2 <-> 4
2784 
2785  if(RRX[1]>.5&&pho.idhep[KK]==-13&&pho.idhep[KJ]==-13){
2786  for( int k=0;k<=3;k++){
2787  double stored=pho.phep[KK][k];
2788  pho.phep[KK][k]=pho.phep[KJ][k];
2789  pho.phep[KJ][k]=stored;
2790 
2791  }
2792  }
2793 
2794  }
2795  }
2796 
2797  // for mc-test with KORALW, mumu from mu mu emissions: END
2798  */
2799 
2800  pho.phep[pho.nhep - i][4] = masslep;
2801  PHCORK(0);
2802  // write in
2803  PHLUPA(1012);
2804  PHOOUT(IPPAR, BOOST, NHEP0);
2805  PHOIN(IPPAR, &BOOST, &NHEP0);
2806  PHLUPA(1013);
2807  } // end of if (JESLI)
2808  } // end of loop over charged particles
2809  }
2810 
2811 
2812 } // namespace Photospp
2813 
PhotosUtilities
Support functions.
Photospp::HEPEVT_struct::ME_channel
static int ME_channel
Number of channel to be used - flag for fortran routines.
Definition: HEPEVT_struct.h:36
Photospp::PhotosMEforW::PHOBWnlo
static void PHOBWnlo(double *WT)
Definition: forW-MEc.cc:918
Photospp::Photos::IfPhot
static bool IfPhot
Flag for generating emission of photons.
Definition: Photos.h:219
Photospp::Photos::meCorrectionWtForScalar
static bool meCorrectionWtForScalar
Flag for complete effects of matrix element (in scalars decays)
Definition: Photos.h:204
Photospp::Photos::randomDouble
static double(* randomDouble)()
Pointer to random generator function.
Definition: Photos.h:226
Photospp::Photos::IfPair
static bool IfPair
Flag for generating emission of pairs.
Definition: Photos.h:216
Photospp::HEPEVT
Definition of the PHOEVT common block.
Definition: photosC.h:11