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"
12using std::cout;
13using std::endl;
14using std::max;
15using namespace Photospp;
16using namespace PhotosUtilities;
17
18namespace Photospp
19{
20
21// Instantiating structs declared in photosC.h
22
23struct HEPEVT hep;
24struct HEPEVT pho;
25
26struct PHOCOP phocop;
27struct PHNUM phnum;
28struct PHOKEY phokey;
29struct PHOSTA phosta;
30struct PHOLUN pholun;
31struct PHOPHS phophs;
32struct TOFROM tofrom;
33struct PHOPRO phopro;
34struct PHOREST phorest;
35struct PHWT phwt;
36struct PHOCORWT phocorwt;
37struct PHOMOM phomom;
38struct PHOCMS phocms;
39struct PHOEXP phoexp;
40struct PHLUPY phlupy;
41
42/** Logical function used deep inside algorithm to check if emitted
43 particles are to emit. For mother it blocks the vertex,
44 but for daughters individually: bad sisters will not prevent electron to emit.
45 top quark has further exception method. */
46bool 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
73double 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++){ //! or jdahep[1-i][2-i]
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
151double 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
225double 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//----------------------------------------------------------------------
334void 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
421void 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
491void 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
537void 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
589void 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
668void 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
746void 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
964void 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
1138void 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
1210double 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
1266double 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
1367double 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
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
1446void 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//----------------------------------------------------------------------
1575void 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//----------------------------------------------------------------------
1706void 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
1790void 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
1886void 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);
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.25*(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
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
2040void 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
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
2274void 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;
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
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
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//----------------------------------------------------------------------
2356void 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;
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;
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;
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;
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 ----------------------------------------------------------------------*/
2543void 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 /*
2625if (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
static void PHOBWnlo(double *WT)
Definition forW-MEc.cxx:914
static bool meCorrectionWtForScalar
static double(* randomDouble)()