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