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