Belle II Software  release-06-02-00
pairs.cc
1 #include "Photos.h"
2 #include "pairs.h"
3 
4 #include <cmath>
5 #include <stdio.h>
6 #include <stdlib.h>
7 using namespace Photospp;
8 
9 namespace Photospp {
10 
11 
12 //
13  inline double xlam(double A, double B, double C) {return sqrt((A - B - C) * (A - B - C) - 4.0 * B * C);}
14 
15  inline double max(double a, double b)
16  {
17  return (a > b) ? a : b;
18  }
19  //
20  //extern "C" void varran_( double RRR[], int *N);
21  double angfi(double X, double Y)
22  {
23  double THE;
24  //const double PI=3.141592653589793238462643;
25 
26  if (X == 0.0 && Y == 0.0) return 0.0;
27  if (fabs(Y) < fabs(X)) {
28  THE = atan(fabs(Y / X));
29  if (X <= 0.0) THE = PI - THE;
30  } else {
31  THE = acos(X / sqrt(X * X + Y * Y));
32  }
33  if (Y < 0.0) THE = 2 * PI - THE;
34  return THE;
35  }
36 
37  double angxy(double X, double Y)
38  {
39  double THE;
40  //const double PI=3.141592653589793238462643;
41 
42  if (X == 0.0 && Y == 0.0) return 0.0;
43 
44  if (fabs(Y) < fabs(X)) {
45  THE = atan(fabs(Y / X));
46  if (X <= 0.0) THE = PI - THE;
47  } else {
48  THE = acos(X / sqrt(X * X + Y * Y));
49  }
50  return THE;
51  }
52 
53  void bostd3(double EXE, double PVEC[4], double QVEC[4])
54  {
55  // ----------------------------------------------------------------------
56  // BOOST ALONG Z AXIS, EXE=EXP(ETA), ETA= HIPERBOLIC VELOCITY.
57  //
58  // USED BY : KORALZ RADKOR
60  int j = 1; // convention of indices of Riemann space must be preserved.
61  double RPL, RMI, QPL, QMI;
62  double RVEC[4];
63 
64 
65  RVEC[1 - j] = PVEC[1 - j];
66  RVEC[2 - j] = PVEC[2 - j];
67  RVEC[3 - j] = PVEC[3 - j];
68  RVEC[4 - j] = PVEC[4 - j];
69  RPL = RVEC[4 - j] + RVEC[3 - j];
70  RMI = RVEC[4 - j] - RVEC[3 - j];
71  QPL = RPL * EXE;
72  QMI = RMI / EXE;
73  QVEC[1 - j] = RVEC[1 - j];
74  QVEC[2 - j] = RVEC[2 - j];
75  QVEC[3 - j] = (QPL - QMI) / 2;
76  QVEC[4 - j] = (QPL + QMI) / 2;
77  }
78 
79 // after investigations PHORO3 of PhotosUtilities.cxx will be used instead
80 // but it must be checked first if it works
81 
82  void rotod3(double ANGLE, double PVEC[4], double QVEC[4])
83  {
84 
85 
86  int j = 1; // convention of indices of Riemann space must be preserved.
87  double CS, SN;
88  // printf ("%5.2f\n",cos(ANGLE));
89  CS = cos(ANGLE) * PVEC[1 - j] - sin(ANGLE) * PVEC[2 - j];
90  SN = sin(ANGLE) * PVEC[1 - j] + cos(ANGLE) * PVEC[2 - j];
91 
92  QVEC[1 - j] = CS;
93  QVEC[2 - j] = SN;
94  QVEC[3 - j] = PVEC[3 - j];
95  QVEC[4 - j] = PVEC[4 - j];
96  }
97 
98 
99 
100  void rotod2(double PHI, double PVEC[4], double QVEC[4])
101  {
102 
103  double RVEC[4];
104  int j = 1; // convention of indices of Riemann space must be preserved.
105  double CS, SN;
106 
107  CS = cos(PHI);
108  SN = sin(PHI);
109 
110  RVEC[1 - j] = PVEC[1 - j];
111  RVEC[2 - j] = PVEC[2 - j];
112  RVEC[3 - j] = PVEC[3 - j];
113  RVEC[4 - j] = PVEC[4 - j];
114 
115  QVEC[1 - j] = CS * RVEC[1 - j] + SN * RVEC[3 - j];
116  QVEC[2 - j] = RVEC[2 - j];
117  QVEC[3 - j] = -SN * RVEC[1 - j] + CS * RVEC[3 - j];
118  QVEC[4 - j] = RVEC[4 - j];
119  // printf ("%15.12f %15.12f %15.12f %15.12f \n",QVEC[0],QVEC[1],QVEC[2],QVEC[3]);
120  // exit(-1);
121  }
122 
123  void lortra(int KEY, double PRM, double PNEUTR[4], double PNU[4], double PAA[4], double PP[4], double PE[4])
124  {
125  // ---------------------------------------------------------------------
126  // THIS ROUTINE PERFORMS LORENTZ TRANSFORMATION ON MANY 4-VECTORS
127  // KEY =1 BOOST ALONG 3RD AXIS
128  // =2 ROTATION AROUND 2ND AXIS
129  // =3 ROTATION AROUND 3RD AXIS
130  // PRM TRANSFORMATION PARAMETER - ANGLE OR EXP(HIPERANGLE).
131  //
132  // called by : RADCOR
133  // ---------------------------------------------------------------------
134  if (KEY == 1) {
135  bostd3(PRM, PNEUTR, PNEUTR);
136  bostd3(PRM, PNU , PNU);
137  bostd3(PRM, PAA , PAA);
138  bostd3(PRM, PE , PE);
139  bostd3(PRM, PP , PP);
140  } else if (KEY == 2) {
141  rotod2(PRM, PNEUTR, PNEUTR);
142  rotod2(PRM, PNU , PNU);
143  rotod2(PRM, PAA , PAA);
144  rotod2(PRM, PE , PE);
145  rotod2(PRM, PP , PP);
146  } else if (KEY == 3) {
147  rotod3(PRM, PNEUTR, PNEUTR);
148  rotod3(PRM, PNU , PNU);
149  rotod3(PRM, PAA , PAA);
150  rotod3(PRM, PE , PE);
151  rotod3(PRM, PP , PP);
152  } else {
153  printf(" STOP IN LOTRA. WRONG KEYTRA");
154  exit(-1);
155  }
156  }
157 
158  double amast(double VEC[4])
159  {
160  int i = 1; // convention of indices of Riemann space must be preserved
161  double ama = VEC[4 - i] * VEC[4 - i] - VEC[1 - i] * VEC[1 - i] - VEC[2 - i] * VEC[2 - i] - VEC[3 - i] * VEC[3 - i];
162  ama = sqrt(fabs(ama));
163  return ama;
164  }
165 
166  void spaj(int KUDA, double P2[4], double Q2[4], double PP[4], double PE[4])
167  {
168  // **********************
169  // THIS PRINTS OUT FOUR MOMENTA OF PHOTONS
170  // ON OUTPUT UNIT NOUT
171 
172  double SUM[4];
173  const int KLUCZ = 0;
174  if (KLUCZ == 0) return;
175 
176  printf(" %10i =====================SPAJ==================== \n", KUDA);
177  printf(" P2 %18.13f %18.13f %18.13f %18.13f \n", P2[0], P2[1], P2[2], P2[3]);
178  printf(" Q2 %18.13f %18.13f %18.13f %18.13f \n", Q2[0], Q2[1], Q2[2], Q2[3]);
179  printf(" PE %18.13f %18.13f %18.13f %18.13f \n", PE[0], PE[1], PE[2], PE[3]);
180  printf(" PP %18.13f %18.13f %18.13f %18.13f \n", PP[0], PP[1], PP[2], PP[3]);
181 
182  for (int k = 0; k <= 3; k++) SUM[k] = P2[k] + Q2[k] + PE[k] + PP[k];
183 
184  printf("SUM %18.13f %18.13f %18.13f %18.13f \n", SUM[0], SUM[1], SUM[2], SUM[3]);
185  }
186 
187 //extern "C" {
188  struct PARKIN {
189  double fi0; // FI0
190  double fi1; // FI1
191  double fi2; // FI2
192  double fi3; // FI3
193  double fi4; // FI4
194  double fi5; // FI5
195  double th0; // TH0
196  double th1; // TH1
197  double th3; // TH3
198  double th4; // TH4
199  double parneu; // PARNEU
200  double parch; // PARCH
201  double bpar; // BPAR
202  double bsta; // BSTA
203  double bstb; // BSTB
204  } parkin;
205 //}
206 
207 //struct PARKIN parkin;
208 
209  void partra(int IBRAN, double PHOT[4])
210  {
211 
212 
213 
214  rotod3(-parkin.fi0, PHOT, PHOT);
215  rotod2(-parkin.th0, PHOT, PHOT);
216  bostd3(parkin.bsta, PHOT, PHOT);
217  rotod3(-parkin.fi1, PHOT, PHOT);
218  rotod2(-parkin.th1, PHOT, PHOT);
219  rotod3(-parkin.fi2, PHOT, PHOT);
220 
221  if (IBRAN == -1) {
222  bostd3(parkin.parneu, PHOT, PHOT);
223  } else {
224  bostd3(parkin.parch, PHOT, PHOT);
225  }
226 
227  rotod3(-parkin.fi3, PHOT, PHOT);
228  rotod2(-parkin.th3, PHOT, PHOT);
229  bostd3(parkin.bpar, PHOT, PHOT);
230  rotod3(parkin.fi4, PHOT, PHOT);
231  rotod2(-parkin.th4, PHOT, PHOT);
232  rotod3(-parkin.fi5, PHOT, PHOT);
233  rotod3(parkin.fi2, PHOT, PHOT);
234  rotod2(parkin.th1, PHOT, PHOT);
235  rotod3(parkin.fi1, PHOT, PHOT);
236  bostd3(parkin.bstb, PHOT, PHOT);
237  rotod2(parkin.th0, PHOT, PHOT);
238  rotod3(parkin.fi0, PHOT, PHOT);
239 
240  }
241 
242 
243  void trypar(bool* JESLI, double* pSTRENG, double AMCH, double AMEL, double PA[4], double PB[4], double PE[4], double PP[4],
244  bool* sameflav)
245  {
246  double& STRENG = *pSTRENG;
247  // COMMON /PARKIN/
248  double& FI0 = parkin.fi0;
249  double& FI1 = parkin.fi1;
250  double& FI2 = parkin.fi2;
251  double& FI3 = parkin.fi3;
252  double& FI4 = parkin.fi4;
253  double& FI5 = parkin.fi5;
254  double& TH0 = parkin.th0;
255  double& TH1 = parkin.th1;
256  double& TH3 = parkin.th3;
257  double& TH4 = parkin.th4;
258  double& PARNEU = parkin.parneu;
259  double& PARCH = parkin.parch;
260  double& BPAR = parkin.bpar;
261  double& BSTA = parkin.bsta;
262  double& BSTB = parkin.bstb;
263 
264  double PNEUTR[4], PAA[4], PHOT[4], PSUM[4];
265  double VEC[4];
266  double RRR[8];
267  bool JESLIK;
268  //const double PI=3.141592653589793238462643;
269  const double ALFINV = 137.01;
270  const int j = 1; // convention of indices of Riemann space must be preserved.
271 
272  PA[4 - j] = max(PA[4 - j], sqrt(PA[1 - j] * PA[1 - j] + PA[2 - j] * PA[2 - j] + PA[3 - j] * PA[3 - j]));
273  PB[4 - j] = max(PB[4 - j], sqrt(PB[1 - j] * PB[1 - j] + PB[2 - j] * PB[2 - j] + PB[3 - j] * PB[3 - j]));
274 
275  // 4-MOMENTUM OF THE NEUTRAL SYSTEM
276  for (int k = 0; k <= 3; k++) {
277  PE[k] = 0.0;
278  PP[k] = 0.0;
279  PSUM[k] = PA[k] + PB[k];
280  PAA[k] = PA[k];
281  PNEUTR[k] = PB[k];
282  }
283  if ((PAA[4 - j] + PNEUTR[4 - j]) < 0.01) {
284  printf(" too small energy to emit %10.7f\n", PAA[4 - j] + PNEUTR[4 - j]);
285  *JESLI = false;
286  return;
287  }
288 
289  // MASSES OF THE NEUTRAL AND CHARGED SYSTEMS AND OVERALL MASS
290  // FIRST WE HAVE TO GO TO THE RESTFRAME TO GET RID OF INSTABILITIES
291  // FROM BHLUMI OR ANYTHING ELSE
292  // THIRD AXIS ALONG PNEUTR
293  double X1 = PSUM[1 - j];
294  double X2 = PSUM[2 - j];
295  FI0 = angfi(X1, X2);
296  X1 = PSUM[3 - j];
297  X2 = sqrt(PSUM[1 - j] * PSUM[1 - j] + PSUM[2 - j] * PSUM[2 - j]);
298  TH0 = angxy(X1, X2) ;
299  spaj(-2, PNEUTR, PAA, PP, PE);
300  lortra(3, -FI0, PNEUTR, VEC, PAA, PP, PE);
301  lortra(2, -TH0, PNEUTR, VEC, PAA, PP, PE);
302  rotod3(-FI0, PSUM, PSUM);
303  rotod2(-TH0, PSUM, PSUM);
304  BSTA = (PSUM[4 - j] - PSUM[3 - j]) / sqrt(PSUM[4 - j] * PSUM[4 - j] - PSUM[3 - j] * PSUM[3 - j]);
305  BSTB = (PSUM[4 - j] + PSUM[3 - j]) / sqrt(PSUM[4 - j] * PSUM[4 - j] - PSUM[3 - j] * PSUM[3 - j]);
306  lortra(1, BSTA, PNEUTR, VEC, PAA, PP, PE);
307  spaj(-1, PNEUTR, PAA, PP, PE);
308  double AMNE = amast(PNEUTR);
309  AMCH = amast(PAA); // to be improved. May be dangerous because of rounding error
310  if (AMCH < 0.0) AMCH = AMEL;
311  if (AMNE < 0.0) AMNE = 0.0;
312  double AMTO = PAA[4 - j] + PNEUTR[4 - j];
313 
314 
315  for (int k = 0; k <= 7; k++) RRR[k] = Photos::randomDouble();
316 
317  if (STRENG == 0.0) {*JESLI = false; return;}
318 
319  double PRHARD;
320  PRHARD = STRENG // NOTE: logs from phase space presamplers not MEs
321  * 0.5 * (1.0 / PI / ALFINV) * (1.0 / PI / ALFINV) // normalization of triple log 1/36 from
322  // journals.aps.org/prd/pdf/10.1103/PhysRevD.49.1178
323  // must come from rejection
324  // 0.5 is because it is for 1-leg only
325  // STRENG=0,5 because of calls before and after photons
326  // other logs should come from rejection
327  * 2 * log(AMTO / AMEL / 2.0) // virtuality
328  *log(AMTO / AMEL / 2.0) // soft
329  *log((AMTO * AMTO + 2 * AMCH * AMCH) / 2.0 / AMCH / AMCH); // collinear
330  // ZBW-2021
331  // artificial ad hoc increase of probability for e+e-/mu+mu- pair appearance.
332  // Should be calculated from MXX GXX etc. but now it remains a shadow of QED.
333  if (darkr.ifspecial == 1) {
334  if (AMEL < 0.001) PRHARD = PRHARD * darkr.NormFact;
335  else PRHARD = PRHARD * darkr.NormFmu; // for muons we need even more.
336  // PRHARD= PRHARD*darkr.NormFact;
337  //if(AMEL>0.001) PRHARD=PRHARD*darkr.NormFmu;
338  }
339 // printf(" PRHARD/amel= %18.13f %18.13f \n",PRHARD, AMEL);
340  // ZBW-2021 end
341 
342  double FREJECT = 2.; // to make room for interference second pair posiblty.
343  PRHARD = PRHARD * FREJECT;
344  // PRHARD=PRHARD*50; // to increase number of pairs in test of mu mu from mu
345  // fror mumuee set *15
346  // enforces hard pairs to be generated 'always'
347  // for the sake of tests with high statistics, also for flat phase space.
348  // PRHARD=0.99* STRENG*2;
349  // STRENG=0.0;
350  if (PRHARD > 1.0) {
351  printf(" stop from Photos pairs.cxx PRHARD= %18.13f \n", PRHARD);
352  exit(0);
353  }
354 // delta is for tests with PhysRevD.49.1178, default is AMTO*2 no restriction on pair phase space
355  double delta = AMTO * 2; //5;//.125; //AMTO*2; //.125; //AMTO*2; ;0.25;
356 
357 
358  if (RRR[7 - j] > PRHARD) { // compensate crude probablilities; for pairs from consecutive sources
359  STRENG = STRENG / (1.0 - PRHARD);
360  *JESLI = false;
361  return;
362  } else STRENG = 0.0;
363 
364 
365 
366 
367 
368  //
369 
370  //virtuality of lepton pair
371  // ZBW-2021
372  // mass and width of intermediate resoinance.
373  double XMP, ALP1, ALP2, ALP;
374  if (darkr.ifspecial == 1) {
375  ALP1 = atan((4 * AMEL * AMEL - darkr.MXX * darkr.MXX) / darkr.MXX / darkr.GXX);
376  ALP2 = atan((AMTO * AMTO - darkr.MXX) / darkr.MXX / darkr.GXX);
377  ALP = ALP1 + RRR[1 - j] * (ALP2 - ALP1);
378  XMP = sqrt(darkr.MXX * darkr.MXX + darkr.MXX * darkr.GXX * tan(ALP));
379  // ZBW-2021 end
380  } else {
381  XMP = 2.0 * AMEL * exp(RRR[1 - j] * log(AMTO / 2.0 / AMEL));
382  // XMP=2.0*AMEL*2.0*AMEL+RRR[1-j]*(AMTO-2.0*AMEL)*(AMTO-2.0*AMEL); XMP=sqrt(XMP); // option of no presampler
383  }
384 
385  // energy of lepton pair replace virtuality of CHAR+NEU system in phase space parametrization
386  double XPmin = 2.0 * AMEL;
387  if (darkr.ifspecial == 1) {
388  XPmin = darkr.MXX - 5 * darkr.GXX;
389  if (XPmin < 2.0 * AMEL) XPmin = 2.0 * AMEL;
390  // ZBW-2021 end
391  }
392 
393  double XPdelta = AMTO - XPmin;
394  double XP = XPmin * exp(RRR[2 - j] * log((XPdelta + XPmin) / XPmin));
395  // XP= XPmin +RRR[2-j]*XPdelta; // option of no presampler
396  double XMK2 = (AMTO * AMTO + XMP * XMP) - 2.0 * AMTO * XP;
397 
398  // angles of lepton pair
399  double eps = 4.0 * AMCH * AMCH / AMTO / AMTO;
400  if (darkr.ifspecial == 1) {
401  eps = (darkr.MXX - 5 * darkr.GXX) * (darkr.MXX - 5 * darkr.GXX) / AMTO / AMTO;
402  if (eps < 4.0 * AMCH * AMCH / AMTO / AMTO) eps = 4.0 * AMCH * AMCH / AMTO / AMTO;
403  // ZBW-2021 end
404  }
405  double C1 = 1.0 + eps - eps * exp(RRR[3 - j] * log((2 + eps) / eps));
406  // C1=1.0-2.0*RRR[3-j]; // option of no presampler
407  double FIX1 = 2.0 * PI * RRR[4 - j];
408 
409  // angles of lepton in restframe of lepton pair
410  double C2 = 1.0 - 2.0 * RRR[5 - j];
411  double FIX2 = 2.0 * PI * RRR[6 - j];
412 
413 
414 
415  // histograming .......................
416  JESLIK = (XP < ((AMTO * AMTO + XMP * XMP - (AMCH + AMNE) * (AMCH + AMNE)) / 2.0 / AMTO));
417  double WTA = 0.0;
418  WTA = WTA * 5;
419  if (JESLIK) WTA = 1.0;
420  // GMONIT( 0,101 ,WTA,1D0,0D0)
421  JESLIK = (XMP < (AMTO - AMNE - AMCH));
422  WTA = 0.0;
423  if (JESLIK) WTA = 1.0;
424  // GMONIT( 0,102 ,WTA,1D0,0D0)
425  JESLIK = (XMP < (AMTO - AMNE - AMCH)) && (XP > XMP);
426 
427  WTA = 0.0;
428  if (JESLIK) WTA = 1.0;
429  // GMONIT( 0,103 ,WTA,1D0,0D0)
430  JESLIK =
431  (XMP < (AMTO - AMNE - AMCH)) &&
432  (XP > XMP) &&
433  (XP < ((AMTO * AMTO + XMP * XMP - (AMCH + AMNE) * (AMCH + AMNE)) / 2.0 / AMTO));
434  WTA = 0.0;
435  if (JESLIK) WTA = 1.0;
436  // GMONIT( 0,104 ,WTA,1D0,0D0)
437  // end of histograming ................
438 
439  // phase space limits rejection variable
440  *JESLI = (RRR[7 - j] < PRHARD) &&
441  (XMP < (AMTO - AMNE - AMCH)) &&
442  (XP > XMP) &&
443  (XP < ((AMTO * AMTO + XMP * XMP - (AMCH + AMNE) * (AMCH + AMNE)) / 2.0 / AMTO));
444 
445 
446 // rejection for phase space restriction: for tests with PhysRevD.49.1178
447  *JESLI = *JESLI && XP < delta;
448  if (!*JESLI) return;
449 
450  // for events in phase: jacobians weights etc.
451 
452  // virtuality of added lepton pair
453  double F = (AMTO * AMTO - 4.0 * AMEL * AMEL) // flat phase space
454  / (AMTO * AMTO - 4.0 * AMEL * AMEL) * XMP * XMP; // use this when presampler is on (log moved to PRHARD)
455  // ZBW-2021
456  if (darkr.ifspecial == 1) F = (ALP2 - ALP1) * ((XMP * XMP - darkr.MXX * darkr.MXX) * (XMP * XMP - darkr.MXX * darkr.MXX) +
457  (darkr.MXX * darkr.GXX) * (darkr.MXX * darkr.GXX)) / (darkr.MXX * darkr.GXX);
458  // ZBW-2021 end
459  // Energy of added lepton pair represented by virtuality of CH+N pair
460  double G = 2 * AMTO * XPdelta // flat phase space
461  / (2 * AMTO * XPdelta) * 2 * AMTO * XP; // use this when presampler is on (log moved to PRHARD)
462 
463 
464  // scattering angle of emitted lepton pair (also flat factors for other angles)
465  double H = 2.0 // flat phase space
466  / 2.0 * (1.0 + eps - C1) / 2.0; // use this when presampler is on (log moved to PRHARD)
467 
468  double H1 = 2.0 // for other generation arm: char neutr replaced
469  / 2.0 * (1.0 + eps - C1);
470  double H2 = 2.0
471  / 2.0 * (1.0 + eps + C1);
472  H = 1. / (0.5 / H1 + 0.5 / H2);
473 
474  //*2*PI*4*PI /2/PI/4/PI; // other angles normalization of transformation to random numbers.
475 
476  double XPMAX = (AMTO * AMTO + XMP * XMP - (AMCH + AMNE) * (AMCH + AMNE)) / 2.0 / AMTO;
477 
478  double YOT3 = F * G * H; // jacobians for phase space variables
479  double YOT2 = // lambda factors:
480  xlam(1.0, AMEL * AMEL / XMP / XMP, AMEL * AMEL / XMP / XMP) *
481  xlam(1.0, XMK2 / AMTO / AMTO, XMP * XMP / AMTO / AMTO) *
482  xlam(1.0, AMCH * AMCH / XMK2, AMNE * AMNE / XMK2)
483  / xlam(1.0, AMCH * AMCH / AMTO / AMTO, AMNE * AMNE / AMTO / AMTO);
484  // if(darkr.ifspecial==1) YOT2=YOT2/xlam(1.0,XMK2/AMTO/AMTO,XMP*XMP/AMTO/AMTO);
485 
486  //C histograming .......................
487  // GMONIT( 0,105 ,WT ,1D0,0D0)
488  // GMONIT( 0,106 ,YOT1,1D0,0D0)
489  // GMONIT( 0,107 ,YOT2,1D0,0D0)
490  // GMONIT( 0,108 ,YOT3,1D0,0D0)
491  // GMONIT( 0,109 ,YOT4,1D0,0D0)
492  // end of histograming ................
493 
494 
495  //
496  //
497  // FRAGMENTATION COMES !!
498  //
499  // THIRD AXIS ALONG PNEUTR
500  X1 = PNEUTR[1 - j];
501  X2 = PNEUTR[2 - j];
502  FI1 = angfi(X1, X2);
503  X1 = PNEUTR[3 - j];
504  X2 = sqrt(PNEUTR[1 - j] * PNEUTR[1 - j] + PNEUTR[2 - j] * PNEUTR[2 - j]) ;
505  TH1 = angxy(X1, X2);
506  spaj(0, PNEUTR, PAA, PP, PE);
507  lortra(3, -FI1, PNEUTR, VEC, PAA, PP, PE);
508  lortra(2, -TH1, PNEUTR, VEC, PAA, PP, PE);
509  VEC[4 - j] = 0.0;
510  VEC[3 - j] = 0.0;
511  VEC[2 - j] = 0.0;
512  VEC[1 - j] = 1.0;
513  FI2 = angfi(VEC[1 - j], VEC[2 - j]);
514  lortra(3, -FI2, PNEUTR, VEC, PAA, PP, PE);
515  spaj(1, PNEUTR, PAA, PP, PE);
516 
517  // STEALING FROM PAA AND PNEUTR ENERGY FOR THE pair
518  // ====================================================
519  // NEW MOMENTUM OF PAA AND PNEUTR (IN THEIR VERY REST FRAME)
520  // 1) PARAMETERS.....
521  double AMCH2 = AMCH * AMCH;
522  double AMNE2 = AMNE * AMNE;
523  double AMTOST = XMK2;
524  double QNEW = xlam(AMTOST, AMNE2, AMCH2) / sqrt(AMTOST) / 2.0;
525  double QOLD = PNEUTR[3 - j];
526  double GCHAR = (QNEW * QNEW + QOLD * QOLD + AMCH * AMCH) /
527  (QNEW * QOLD + sqrt((QNEW * QNEW + AMCH * AMCH) * (QOLD * QOLD + AMCH * AMCH)));
528  double GNEU = (QNEW * QNEW + QOLD * QOLD + AMNE * AMNE) /
529  (QNEW * QOLD + sqrt((QNEW * QNEW + AMNE * AMNE) * (QOLD * QOLD + AMNE * AMNE)));
530 
531  // GCHAR=(QOLD**2-QNEW**2)/(
532  // & QNEW*SQRT(QOLD**2+AMCH2)+QOLD*SQRT(QNEW**2+AMCH2)
533  // & )
534  // GCHAR=SQRT(1D0+GCHAR**2)
535  // GNEU=(QOLD**2-QNEW**2)/(
536  // & QNEW*SQRT(QOLD**2+AMNE2)+QOLD*SQRT(QNEW**2+AMNE2)
537  // & )
538  // GNEU=SQRT(1D0+GNEU**2)
539  if (GNEU < 1. || GCHAR < 1.) {
540  printf(" PHOTOS TRYPAR GBOOST LT 1., LIMIT OF PHASE SPACE %18.13f %18.13f %18.13f %18.13f %18.13f %18.13f %18.13f %18.13f \n"
541  , GNEU, GCHAR, QNEW, QOLD, AMTO, AMTOST, AMNE, AMCH);
542  return;
543  }
544  PARCH = GCHAR + sqrt(GCHAR * GCHAR - 1.0);
545  PARNEU = GNEU - sqrt(GNEU * GNEU - 1.0);
546 
547  // 2) REDUCTIEV BOOSTS
548  bostd3(PARNEU, VEC , VEC);
549  bostd3(PARNEU, PNEUTR, PNEUTR);
550  bostd3(PARCH, PAA , PAA);
551  spaj(2, PNEUTR, PAA, PP, PE);
552 
553  // TIME FOR THE PHOTON that is electron pair
554  double PMOD = xlam(XMP * XMP, AMEL * AMEL, AMEL * AMEL) / XMP / 2.0;
555  double S2 = sqrt(1.0 - C2 * C2);
556  PP[4 - j] = XMP / 2.0;
557  PP[3 - j] = PMOD * C2;
558  PP[2 - j] = PMOD * S2 * cos(FIX2);
559  PP[1 - j] = PMOD * S2 * sin(FIX2);
560  PE[4 - j] = PP[4 - j];
561  PE[3 - j] = -PP[3 - j];
562  PE[2 - j] = -PP[2 - j];
563  PE[1 - j] = -PP[1 - j];
564  // PHOTON ENERGY and momentum IN THE REDUCED SYSTEM
565  double PENE = (AMTO * AMTO - XMP * XMP - XMK2) / 2.0 / sqrt(XMK2);
566  double PPED = sqrt(PENE * PENE - XMP * XMP);
567  FI3 = FIX1;
568  double COSTHG = C1;
569  double SINTHG = sqrt(1.0 - C1 * C1);
570  X1 = -COSTHG;
571  X2 = SINTHG;
572  TH3 = angxy(X1, X2);
573  PHOT[1 - j] = PMOD * SINTHG * cos(FI3);
574  PHOT[2 - j] = PMOD * SINTHG * sin(FI3);
575  // MINUS BECAUSE AXIS OPPOSITE TO PAA
576  PHOT[3 - j] = -PMOD * COSTHG;
577  PHOT[4 - j] = PENE;
578  // ROTATE TO PUT PHOTON ALONG THIRD AXIS
579  X1 = PHOT[1 - j];
580  X2 = PHOT[2 - j];
581  lortra(3, -FI3, PNEUTR, VEC, PAA, PP, PE);
582  rotod3(-FI3, PHOT, PHOT);
583  lortra(2, -TH3, PNEUTR, VEC, PAA, PP, PE);
584  rotod2(-TH3, PHOT, PHOT);
585  spaj(21, PNEUTR, PAA, PP, PE);
586  // ... now get the pair !
587  double PAIRB = PENE / XMP + PPED / XMP;
588  bostd3(PAIRB, PE, PE);
589  bostd3(PAIRB, PP, PP);
590  spaj(3, PNEUTR, PAA, PP, PE);
591  double GAMM = (PNEUTR[4 - j] + PAA[4 - j] + PP[4 - j] + PE[4 - j]) / AMTO;
592 
593  // TP and ZW: 25.II.2015: fix for cases when mother is very close
594  // to its rest frame and pair is generated after photon emission.
595  // Then GAMM can be slightly less than 1.0 due to rounding error
596  if (GAMM < 1.0) {
597  if (GAMM > 0.9999999) GAMM = 1.0;
598  else {
599  printf("Photos::partra: GAMM = %20.18e\n", GAMM);
600  printf(" BELOW 0.9999999 THRESHOLD!\n");
601  GAMM = 1.0;
602  }
603  }
604 
605  BPAR = GAMM - sqrt(GAMM * GAMM - 1.0);
606  lortra(1, BPAR, PNEUTR, VEC, PAA, PP, PE);
607  bostd3(BPAR, PHOT, PHOT);
608  spaj(4, PNEUTR, PAA, PP, PE);
609  // BACK IN THE TAU REST FRAME BUT PNEUTR NOT YET ORIENTED.
610  X1 = PNEUTR[1 - j];
611  X2 = PNEUTR[2 - j];
612  FI4 = angfi(X1, X2);
613  X1 = PNEUTR[3 - j];
614  X2 = sqrt(PNEUTR[1 - j] * PNEUTR[1 - j] + PNEUTR[2 - j] * PNEUTR[2 - j]);
615  TH4 = angxy(X1, X2);
616  lortra(3, FI4, PNEUTR, VEC, PAA, PP, PE);
617  rotod3(FI4, PHOT, PHOT);
618  lortra(2, -TH4, PNEUTR, VEC, PAA, PP, PE);
619  rotod2(-TH4, PHOT, PHOT);
620  X1 = VEC[1 - j];
621  X2 = VEC[2 - j];
622  FI5 = angfi(X1, X2);
623  lortra(3, -FI5, PNEUTR, VEC, PAA, PP, PE);
624  rotod3(-FI5, PHOT, PHOT);
625  // PAA RESTORES ORIGINAL DIRECTION
626  lortra(3, FI2, PNEUTR, VEC, PAA, PP, PE);
627  rotod3(FI2, PHOT, PHOT);
628  lortra(2, TH1, PNEUTR, VEC, PAA, PP, PE);
629  rotod2(TH1, PHOT, PHOT);
630  lortra(3, FI1, PNEUTR, VEC, PAA, PP, PE);
631  rotod3(FI1, PHOT, PHOT);
632  spaj(10, PNEUTR, PAA, PP, PE);
633  lortra(1, BSTB, PNEUTR, VEC, PAA, PP, PE);
634  lortra(2, TH0, PNEUTR, VEC, PAA, PP, PE);
635  lortra(3, FI0, PNEUTR, VEC, PAA, PP, PE);
636  spaj(11, PNEUTR, PAA, PP, PE);
637 
638 
639 // matrix element as formula 1 from journals.aps.org/prd/pdf/10.1103/PhysRevD.49.1178
640 
641  double pq = PAA[3] * PP[3] - PAA[2] * PP[2] - PAA[1] * PP[1] - PAA[0] * PP[0];
642  pq = pq + PAA[3] * PE[3] - PAA[2] * PE[2] - PAA[1] * PE[1] - PAA[0] * PE[0];
643 
644  double ppq = PNEUTR[3] * PP[3] - PNEUTR[2] * PP[2] - PNEUTR[1] * PP[1] - PNEUTR[0] * PP[0];
645  ppq = ppq + PNEUTR[3] * PE[3] - PNEUTR[2] * PE[2] - PNEUTR[1] * PE[1] - PNEUTR[0] * PE[0];
646  double pq1 = PAA[3] * PP[3] - PAA[2] * PP[2] - PAA[1] * PP[1] - PAA[0] * PP[0];
647  double pq2 = PAA[3] * PE[3] - PAA[2] * PE[2] - PAA[1] * PE[1] - PAA[0] * PE[0];
648 
649  double ppq1 = PNEUTR[3] * PP[3] - PNEUTR[2] * PP[2] - PNEUTR[1] * PP[1] - PNEUTR[0] * PP[0];
650  double ppq2 = PNEUTR[3] * PE[3] - PNEUTR[2] * PE[2] - PNEUTR[1] * PE[1] - PNEUTR[0] * PE[0];
651 
652  double ppp = PNEUTR[3] * PAA[3] - PNEUTR[2] * PAA[2] - PNEUTR[1] * PAA[1] - PNEUTR[0] * PAA[0];
653  double mneutr2 = PNEUTR[3] * PNEUTR[3] - PNEUTR[2] * PNEUTR[2] - PNEUTR[1] * PNEUTR[1] - PNEUTR[0] * PNEUTR[0];
654  double maa2 = PAA[3] * PAA[3] - PAA[2] * PAA[2] - PAA[1] * PAA[1] - PAA[0] * PAA[0];
655 
656  double YOT1 = 1. / 2. / XMP / XMP / XMP / XMP *
657  (4 * (pq1 / pq - ppq1 / ppq) * (pq2 / pq - ppq2 / ppq)
658  - XMP * XMP * (AMCH2 / pq / pq + AMNE2 / ppq / ppq - ppp / pq / ppq - ppp / pq / ppq));
659  // ZBW-2021
660 
661  if (darkr.ifspecial == 1 && darkr.iboson == 1) {
662  YOT1 = YOT1 * XMP * XMP * XMP * XMP / ((darkr.MXX * darkr.MXX - XMP * XMP) * (darkr.MXX * darkr.MXX - XMP * XMP) + darkr.GXX *
663  darkr.GXX * darkr.MXX * darkr.MXX);
664  YOT1 = YOT1 * darkr.GXX * darkr.MXX ; // factor of total rate should be elsewhere
665  // YOT1=YOT1* (AMTO*AMTO-4*AMCH*AMCH)/(AMTO*AMTO); //oct 21 11:00
666  // YOT1=YOT1* XMK2/(XMK2-4*AMCH*AMCH);//* XMK2/(AMTO*AMTO);//oct 21 11:00
667 
668  YOT1 = YOT1 * AMTO / sqrt(XMK2); //oct 21 10:55
669  YOT1 = YOT1 * AMTO / sqrt(XMK2); //oct 21 10:55
670  YOT1 = YOT1 * sqrt((AMTO * AMTO - XMK2 + 2 * AMCH * AMCH) / (AMTO * AMTO - XMK2)); //oct 21 10:55
671  } else if (darkr.ifspecial == 1) {
672  double mcr = 0.5 * darkr.MXX * darkr.MXX;
673  mcr = 0.5 * XMP * XMP;
674  YOT1 = 1. / 2. / XMP / XMP / XMP / XMP *
675  (4 * (1. / (pq + mcr) - 1. / (ppq + mcr)) * (1. / (pq + mcr) - 1. / (ppq + mcr)) * (XMP * XMP + 0 * AMCH * AMCH)
676  - 4 * XMP * XMP / AMTO / AMTO * (AMCH2 / pq / pq + AMNE2 / ppq / ppq - ppp / pq / ppq - ppp / pq / ppq));
677 
678  YOT1 = YOT1 * XMP * XMP * XMP * XMP / ((darkr.MXX * darkr.MXX - XMP * XMP) * (darkr.MXX * darkr.MXX - XMP * XMP) + darkr.GXX *
679  darkr.GXX * darkr.MXX * darkr.MXX);
680  YOT1 = YOT1 * darkr.GXX * darkr.MXX ; // factor of total rate should be elsewhere
681  YOT1 = YOT1 * (AMTO * AMTO - 4 * AMCH * AMCH) / (AMTO * AMTO);
682  YOT1 = YOT1 * XMK2 / (XMK2 - 4 * AMCH * AMCH); //* XMK2/(AMTO*AMTO);
683  YOT1 = YOT1 * AMTO / sqrt(XMK2);
684  YOT1 = YOT1 * AMTO / sqrt(XMK2);
685  YOT1 = YOT1 * sqrt((AMTO * AMTO - XMK2 + 2 * AMCH * AMCH) / (AMTO * AMTO - XMK2)); //sep 5 13:10
686 
687  }
688 // ZBW-2021 end
689 
690 
691  if (*sameflav) {
692 // we interchange: PAA <--> pp
693  for (int k = 0; k <= 3; k++) {
694  double stored = PAA[k];
695  PAA[k] = PE[k];
696  PE[k] = stored;
697  }
698 
699  pq = PAA[3] * PP[3] - PAA[2] * PP[2] - PAA[1] * PP[1] - PAA[0] * PP[0];
700  pq = pq + PAA[3] * PE[3] - PAA[2] * PE[2] - PAA[1] * PE[1] - PAA[0] * PE[0];
701 
702  ppq = PNEUTR[3] * PP[3] - PNEUTR[2] * PP[2] - PNEUTR[1] * PP[1] - PNEUTR[0] * PP[0];
703  ppq = ppq + PNEUTR[3] * PE[3] - PNEUTR[2] * PE[2] - PNEUTR[1] * PE[1] - PNEUTR[0] * PE[0];
704  pq1 = PAA[3] * PP[3] - PAA[2] * PP[2] - PAA[1] * PP[1] - PAA[0] * PP[0];
705  pq2 = PAA[3] * PE[3] - PAA[2] * PE[2] - PAA[1] * PE[1] - PAA[0] * PE[0];
706 
707  ppq1 = PNEUTR[3] * PP[3] - PNEUTR[2] * PP[2] - PNEUTR[1] * PP[1] - PNEUTR[0] * PP[0];
708  ppq2 = PNEUTR[3] * PE[3] - PNEUTR[2] * PE[2] - PNEUTR[1] * PE[1] - PNEUTR[0] * PE[0];
709 
710  ppp = PNEUTR[3] * PAA[3] - PNEUTR[2] * PAA[2] - PNEUTR[1] * PAA[1] - PNEUTR[0] * PAA[0];
711 
712  XMP = -(PP[0] + PE[0]) * (PP[0] + PE[0]) - (PP[1] + PE[1]) * (PP[1] + PE[1])
713  - (PP[2] + PE[2]) * (PP[2] + PE[2]) + (PP[3] + PE[3]) * (PP[3] + PE[3]);
714  XMP = sqrt(fabs(XMP));
715 
716 
717  double YOT1p = 1. / 2. / XMP / XMP / XMP / XMP *
718  (4 * (pq1 / pq - ppq1 / ppq) * (pq2 / pq - ppq2 / ppq)
719  - XMP * XMP * (AMCH2 / pq / pq + AMNE2 / ppq / ppq - ppp / pq / ppq - ppp / pq / ppq));
720  // *(1-XP/XPMAX+0.5*(XP/XPMAX)*(XP/XPMAX)); // A-P kernel divide by (1-XP/XPMAX)?
721  double wtint = 0.; // not yet installed
722  wtint = 1; //(YOT1+YOT1p+wtint)/(YOT1+YOT1p);
723 // ZBW-2021
724 //if(darkr.ifspecial==1) YOT1=YOT1*XMP*XMP*XMP*XMP/((MXX*MXX-XMP*XMP)*(MXX*MXX-XMP*XMP)+GXX*GXX*MXX*MXX);
725 // ZBW-2021 end
726  YOT1 = YOT1 * wtint;
727 
728 // we interchange: PAA <--> pp back into place
729  for (int k = 0; k <= 3; k++) {
730  double stored = PAA[k];
731  PAA[k] = PE[k];
732  PE[k] = stored;
733  }
734  } // end sameflav
735 
736  double WT = YOT1 * YOT2 * YOT3;
737 
738  WT = WT / 8 / FREJECT; // origin must be understood
739  if (darkr.ifspecial == 1 && darkr.ifforce == 1 && AMEL < 0.001) {
740  darkr.Fel = max(darkr.Fel, WT);
741  darkr.Fel = std::min(darkr.Fel, 1.0);
742  WT = WT / darkr.Fel;
743  }
744  if (darkr.ifspecial == 1 && darkr.ifforce == 1 && AMEL > 0.001) {
745  darkr.Fmu = max(darkr.Fmu, WT);
746  darkr.Fmu = std::min(darkr.Fmu, 1.0);
747  WT = WT / darkr.Fmu;
748  }
749  // printf (" from Photos pairs.cxx Fel/Fmu= %15.8f %15.8f %5i \n",darkr.Fel,darkr.Fmu,darkr.ifforce);
750  // if(WT>1.0){
751  // printf (" from Photos pairs.cxx WT= %15.8f \n",WT);
752  // printf (" from Photos pairs.cxx WT= %15.8f %15.8f %15.8f\n",YOT1,YOT2,YOT3);
753  // if(WT>20.0){
754  // printf ("XMP,MXX,GXX = %15.8f %15.8f %15.8f\n",XMP,darkr.MXX,darkr.GXX);
755  // printf ("AMEL AMTO = %15.8f %15.8f\n",AMEL,AMTO);
756  // printf ("ALP1 ALP2 ALP = %15.8f %15.8f %15.8f\n",ALP1,ALP2,ALP);
757  // printf ("YOT1,YOT2,YOT3 = %15.8f %15.8f %15.8f\n",YOT1,YOT2,YOT3);
758  // }
759  // }
760  // printf(" WT/amel/frej,yot1/yot2/yot3 = %15.10f %15.10f %15.10f %15.10f %15.10f %15.10f \n",WT, AMEL,FREJECT,YOT1,YOT2,YOT3);
761  if (RRR[8 - j] > WT) {
762  *JESLI = false;
763  return;
764  }
765 
766  }
767 
768 } // namespace Photospp
769 
Photospp::Photos::randomDouble
static double(* randomDouble)()
Pointer to random generator function.
Definition: Photos.h:226