forZ-MEc.cxx
1#include "forZ-MEc.h"
2#include "Photos.h"
3#include "PhotosUtilities.h"
4#include "photosC.h"
5#include <cmath>
6#include <cstdio>
7#include <cstdlib>
8#include <iostream>
9using std::cout;
10using std::endl;
11using namespace Photospp;
12using namespace PhotosUtilities;
13
14namespace Photospp
15{
16
17 double (*PhotosMEforZ::currentVakPol)(double[4], double[4], double[4], double[4], double[4], int, int) = default_VakPol;
18
19// from photosC.cxx
20
21void PHODMP();
22double PHINT(int idumm);
23// ----------------------------------------------------------------------
24// PROVIDES ELECTRIC CHARGE AND WEAK IZOSPIN OF A FAMILY FERMION
25// IDFERM=1,2,3,4 DENOTES NEUTRINO, LEPTON, UP AND DOWN QUARK
26// NEGATIVE IDFERM=-1,-2,-3,-4, DENOTES ANTIPARTICLE
27// IHELIC=+1,-1 DENOTES RIGHT AND LEFT HANDEDNES ( CHIRALITY)
28// SIZO3 IS THIRD PROJECTION OF WEAK IZOSPIN (PLUS MINUS HALF)
29// AND CHARGE IS ELECTRIC CHARGE IN UNITS OF ELECTRON CHARGE
30// KOLOR IS A QCD COLOUR, 1 FOR LEPTON, 3 FOR QUARKS
31//
32// called by : EVENTE, EVENTM, FUNTIH, .....
33// ----------------------------------------------------------------------
34
35void PhotosMEforZ::GIVIZO(int IDFERM,int IHELIC,double *SIZO3,double *CHARGE,int *KOLOR) {
36 //
37 int IH, IDTYPE, IC, LEPQUA, IUPDOW;
38 if (IDFERM==0 || abs(IDFERM)>4 || abs(IHELIC)!=1){
39 cout << "STOP IN GIVIZO: WRONG PARAMS" << endl;
40 exit(-1);
41 }
42
43 IH =IHELIC;
44 IDTYPE =abs(IDFERM);
45 IC =IDFERM/IDTYPE;
46 LEPQUA=(int)(IDTYPE*0.4999999);
47 IUPDOW=IDTYPE-2*LEPQUA-1;
48 *CHARGE =(-IUPDOW+2.0/3.0*LEPQUA)*IC;
49 *SIZO3 =0.25*(IC-IH)*(1-2*IUPDOW);
50 *KOLOR=1+2*LEPQUA;
51 //** NOTE THAT CONVENTIONALY Z0 COUPLING IS
52 //** XOUPZ=(SIZO3-CHARGE*SWSQ)/SQRT(SWSQ*(1-SWSQ))
53 return;
54}
55
56
57////////////////////////////////////////////////////////////////////////////
58/// //
59/// This routine provides unsophisticated Born differential cross section //
60/// at the crude x-section level, with Z and gamma s-chanel exchange. //
61///////////////////////////////////////////////////////////////////////////
62double PhotosMEforZ::PHBORNM(double svar,double costhe,double T3e,double qe,double T3f,double qf,int NCf){
63
64 double s,Sw2,MZ,MZ2,GammZ,AlfInv,GFermi; // t,MW,MW2,
65 double Ve,Ae,thresh; // sum,deno,
66 double xe,yf,xf,ye,ff0,ff1,amx2,amfin,Vf,Af;
67 double ReChiZ,SqChiZ,RaZ; //,RaW,ReChiW,SqChiW;
68 double Born; //, BornS;
69 // int KeyZet,HadMin,KFbeam;
70 // int i,ke,KFfin,kf,IsGenerated,iKF;
71 int KeyWidFix;
72
73 AlfInv= 137.0359895;
74 GFermi=1.16639e-5;
75
76 //--------------------------------------------------------------------
77 s = svar;
78 //------------------------------
79 // EW paratemetrs taken from BornV
80 MZ=91.187;
81 GammZ=2.50072032;
82 Sw2=.22276773;
83 //------------------------------
84 // Z and gamma couplings to beams (electrons)
85 // Z and gamma couplings to final fermions
86 // Loop over all flavours defined in m_xpar(400+i)
87
88
89 //------ incoming fermion
90 Ve= 2*T3e -4*qe*Sw2;
91 Ae= 2*T3e;
92 //------ final fermion couplings
93 amfin = 0.000511; // m_xpar(kf+6)
94 Vf = 2*T3f -4*qf*Sw2;
95 Af = 2*T3f;
96 if(fabs(costhe) > 1.0){
97 cout << "+++++STOP in PHBORN: costhe>0 =" << costhe << endl;
98 exit(-1);
99 }
100 MZ2 = MZ*MZ;
101 RaZ = (GFermi *MZ2 *AlfInv )/( sqrt(2.0) *8.0 *PI); //
102 RaZ = 1/(16.0*Sw2*(1.0-Sw2));
103 KeyWidFix = 1; // fixed width
104 KeyWidFix = 0; // variable width
105 if( KeyWidFix == 0 ){
106 ReChiZ=(s-MZ2)*s/((s-MZ2)*(s-MZ2)+(GammZ*s/MZ)*(GammZ*s/MZ)) *RaZ; // variable width
107 SqChiZ= s*s/((s-MZ2)*(s-MZ2)+(GammZ*s/MZ)*(GammZ*s/MZ)) *RaZ*RaZ; // variable width
108 }
109 else{
110 ReChiZ=(s-MZ2)*s/((s-MZ2)*(s-MZ2)+(GammZ*MZ)*(GammZ*MZ)) *RaZ; // fixed width
111 SqChiZ= s*s/((s-MZ2)*(s-MZ2)+(GammZ*MZ)*(GammZ*MZ)) *RaZ*RaZ; // fixed width
112 }
113 xe= Ve*Ve +Ae*Ae;
114 xf= Vf*Vf +Af*Af;
115 ye= 2*Ve*Ae;
116 yf= 2*Vf*Af;
117 ff0= qe*qe*qf*qf +2*ReChiZ*qe*qf*Ve*Vf +SqChiZ*xe*xf;
118 ff1= +2*ReChiZ*qe*qf*Ae*Af +SqChiZ*ye*yf;
119 Born = (1.0+ costhe*costhe)*ff0 +2.0*costhe*ff1;
120 // Colour factor
121 Born = NCf*Born;
122 // Crude method of correcting threshold, cos(theta) depencence incorrect!!!
123 if( svar < 4.0*amfin*amfin){
124 thresh=0.0;
125 }
126 else if(svar < 16.0*amfin*amfin){
127 amx2=4.0*amfin*amfin/svar;
128 thresh=sqrt(1.0-amx2)*(1.0+amx2/2.0);
129 }
130 else{
131 thresh=1.0;
132 }
133
134 Born= Born*thresh;
135 return Born;
136}
137
138
139// ----------------------------------------------------------------------
140// THIS ROUTINE CALCULATES BORN ASYMMETRY.
141// IT EXPLOITS THE FACT THAT BORN X. SECTION = A + B*C + D*C**2
142//
143// called by : EVENTM
144// ----------------------------------------------------------------------
145//
146double PhotosMEforZ::AFBCALC(double SVAR,int IDEE,int IDFF){
147 int KOLOR,KOLOR1;
148 double T3e,qe,T3f,qf,A,B;
149 GIVIZO(IDEE,-1,&T3e,&qe,&KOLOR);
150 GIVIZO(IDFF,-1,&T3f,&qf,&KOLOR1);
151
152 A=PHBORNM(SVAR,0.5,T3e,qe,T3f,qf,KOLOR*KOLOR1);
153 B=PHBORNM(SVAR,-0.5,T3e,qe,T3f,qf,KOLOR*KOLOR1);
154 return (A-B)/(A+B)*5.0/2.0 *3.0/8.0;
155}
156
157
158int PhotosMEforZ::GETIDEE(int IDE){
159
160 int IDEE;
161 IDEE=-555;
162 if((IDE==11) || (IDE== 13) || (IDE== 15)){
163 IDEE=2;
164 }
165 else if((IDE==-11) || (IDE==-13) || (IDE==-15)){
166 IDEE=-2;
167 }
168 else if((IDE== 12) || (IDE== 14) || (IDE== 16)){
169 IDEE=1;
170 }
171 else if((IDE==-12) || (IDE==-14) || (IDE==-16)){
172 IDEE=-1;
173 }
174 else if((IDE== 1) || (IDE== 3) || (IDE== 5)){
175 IDEE=4;
176 }
177 else if((IDE== -1) || (IDE== -3) || (IDE== -5)){
178 IDEE=-4;
179 }
180 else if((IDE== 2) || (IDE== 4) || (IDE== 6)){
181 IDEE=3;
182 }
183 else if((IDE==- 2) || (IDE== -4) || (IDE== -6)){
184 IDEE=-3;
185 }
186 if(IDEE==-555) {cout << " ERROR IN GETIDEE of PHOTS Z-ME: I3= &4i"<<IDEE<<endl;}
187 return IDEE;
188}
189
190
191
192
193//----------------------------------------------------------------------
194//
195// PHASYZ: PHotosASYmmetry of Z
196//
197// Purpose: Calculates born level asymmetry for Z
198// between distributions (1-C)**2 and (1+C)**2
199// At present dummy, requrires effective Z and gamma
200// Couplings and also spin polarization states
201// For initial and final states.
202// To be correct this function need to be tuned
203// to host generator. Axes orientation polarisation
204// conventions etc etc.
205// Modularity of PHOTOS would break.
206//
207// Input Parameters: SVAR
208//
209// Output Parameters: Function value
210//
211// Author(s): Z. Was Created at: 10/12/05
212// Last Update: 19/06/13
213//
214//----------------------------------------------------------------------
215double PhotosMEforZ::PHASYZ(double SVAR,int IDE, int IDF){
216
217 double AFB;
218 int IDEE,IDFF;
219
220 IDEE=abs(GETIDEE(IDE));
221 IDFF=abs(GETIDEE(IDF));
222 AFB= -AFBCALC(SVAR,IDEE,IDFF);
223 // AFB=0
224 return 4.0/3.0*AFB;
225 // write(*,*) 'IDE=',IDE,' IDF=',IDF,' SVAR=',SVAR,'AFB=',AFB
226}
227
228//----------------------------------------------------------------------
229//
230// PHWTNLO: PHotosWTatNLO
231//
232// Purpose: calculates instead of interference weight
233// complete NLO weight for vector boson decays
234// of pure vector (or pseudovector) couplings
235// Proper orientation of beams required.
236// This is not standard in PHOTOS.
237// At NLO more tuning than in standard is needed.
238//
239//
240//
241// Input Parameters: as in function declaration
242//
243// Output Parameters: Function value
244//
245// Author(s): Z. Was Created at: 08/12/05
246// Last Update: 20/06/13
247//
248//----------------------------------------------------------------------
249double PhotosMEforZ::Zphwtnlo(double svar,double xk,int IDHEP3,int IREP,double qp[4],double qm[4],double ph[4],double pp[4],double pm[4],double COSTHG,double BETA,double th1,int IDE,int IDF){
250 double C,s,xkaM,xkaP,t,u,t1,u1,BT,BU;
251 double waga,wagan2;
252 static int i=1;
253 int IBREM;
254
255
256 // IBREM is spurious but it numbers branches of MUSTRAAL
257 IBREM=1;
258 if (IREP==1) IBREM=-1;
259
260 // we calculate C and S, note that TH1 exists in MUSTRAAL as well.
261
262 C=cos(th1); // this parameter is calculated outside of the class
263
264 // from off line application we had:
265 if(IBREM==-1) C=-C;
266 // ... we may want to re-check it.
267 s=sqrt(1.0-C*C);
268
269 if (IBREM==1){
270 xkaM=(qp[4-i]*ph[4-i]-qp[3-i]*ph[3-i]-qp[2-i]*ph[2-i]-qp[1-i]*ph[1-i])/xk;
271 xkaP=(qm[4-i]*ph[4-i]-qm[3-i]*ph[3-i]-qm[2-i]*ph[2-i]-qm[1-i]*ph[1-i])/xk;
272 }
273 else{
274 xkaP=(qp[4-i]*ph[4-i]-qp[3-i]*ph[3-i]-qp[2-i]*ph[2-i]-qp[1-i]*ph[1-i])/xk;
275 xkaM=(qm[4-i]*ph[4-i]-qm[3-i]*ph[3-i]-qm[2-i]*ph[2-i]-qm[1-i]*ph[1-i])/xk;
276 }
277
278 // XK=2*PHEP(4,nhep)/PHEP(4,1)/xphmax ! it is not used becuse here
279 // ! order of emissions is meaningless
280 //
281 // DELTA=2*PHEP(5,4)**2/svar/(1+(1-XK)**2)*(xKAP/xKAM+xKAM/xKAP)
282 // waga=svar/4./xkap
283 // waga=waga*(1.D0-COSTHG*BETA) ! sprawdzone 1= svar/xKAp/4 * (1.D0-COSTHG*BETA)
284 // waga=waga*(1-delta) /wt2 ! sprawdzone ze to jest =2/(1.D0+COSTHG*BETA)
285 // ! czyli ubija de-interferencje
286
287
288 // this is true only for intermediate resonances with afb=0!
289 t =2*(qp[4-i]*pp[4-i]-qp[3-i]*pp[3-i]-qp[2-i]*pp[2-i]-qp[1-i]*pp[1-i]);
290 u =2*(qm[4-i]*pp[4-i]-qm[3-i]*pp[3-i]-qm[2-i]*pp[2-i]-qm[1-i]*pp[1-i]);
291 u1=2*(qp[4-i]*pm[4-i]-qp[3-i]*pm[3-i]-qp[2-i]*pm[2-i]-qp[1-i]*pm[1-i]);
292 t1=2*(qm[4-i]*pm[4-i]-qm[3-i]*pm[3-i]-qm[2-i]*pm[2-i]-qm[1-i]*pm[1-i]);
293
294 // basically irrelevant lines ...
295 t =t - (qp[4-i]*qp[4-i]-qp[3-i]*qp[3-i]-qp[2-i]*qp[2-i]-qp[1-i]*qp[1-i]);
296 u =u - (qm[4-i]*qm[4-i]-qm[3-i]*qm[3-i]-qm[2-i]*qm[2-i]-qm[1-i]*qm[1-i]);
297 u1=u1- (qp[4-i]*qp[4-i]-qp[3-i]*qp[3-i]-qp[2-i]*qp[2-i]-qp[1-i]*qp[1-i]);
298 t1=t1- (qm[4-i]*qm[4-i]-qm[3-i]*qm[3-i]-qm[2-i]*qm[2-i]-qm[1-i]*qm[1-i]);
299
300 // we adjust to what is f-st,s-nd beam flavour
301 if (IDE*IDHEP3>0){
302 BT=1.0+PHASYZ(svar,IDE,IDF);
303 BU=1.0-PHASYZ(svar,IDE,IDF);
304 }
305 else{
306 BT=1.0-PHASYZ(svar,IDE,IDF);
307 BU=1.0+PHASYZ(svar,IDE,IDF);
308 }
309 wagan2=2*(BT*t*t+BU*u*u+BT*t1*t1+BU*u1*u1)
310 /(1+(1-xk)*(1-xk))* 2.0/(BT*(1-C)*(1-C)+BU*(1+C)*(1+C))/svar/svar;
311
312 wagan2=wagan2*VakPol(qp,qm,ph,pp,pm,IDE,IDF); // default VakWpol returns 1. Hook for the user function
313 //! waga=waga*wagan2
314 //! waga=waga*(1-delta) /wt2 ! sprawdzone ze to jest =2/(1.D0+COSTHG*BETA)
315 waga=2/(1.0+COSTHG*BETA)*wagan2;
316 //! % * svar/4./xkap*(1.D0-COSTHG*BETA)*sqrt(1.0-xk)
317
318
319 if(wagan2<=3.8) return waga;
320
321 //
322 // exceptional case wagan2>3.8
323 // it should correspond to extremely high bremssthahlung in multiphot conf.
324 //
325 FILE *PHLUN = stdout;
326
327
328 // fprintf(PHLUN,"") 'phwtnlo= ',phwtnlo
329 // fprintf(PHLUN,"") 'idhepy= ',IDHEP[1-i],IDHEP[2-i],IDHEP[3-i],IDHEP[4-i],IDHEP(5)
330 fprintf(PHLUN," IDE= %i IDF= %i",IDE,IDF);
331 fprintf(PHLUN,"bt,bu,bt+bu= %f %f %f",BT,BU,BT+BU);
332 PHODMP(); // we will activate this once PHODMP(); is re-written
333
334 fprintf(PHLUN," ");
335 fprintf(PHLUN,"%i %i <-- IREP,IBREM", IREP,IBREM);
336 //! fprintf(PHLUN,"%f %f %f %f %f") 'pneutr= ',phomom.pneutr[0],phomom.pneutr[1],phomom.pneutr[2],phomom.pneutr[3],phomom.pneutr[4];
337 fprintf(PHLUN,"%f %f %f %f qp = ",qp[0],qp[1],qp[2],qp[3]);
338 fprintf(PHLUN,"%f %f %f %f qm = ",qm[0],qm[1],qm[2],qm[3]);
339 fprintf(PHLUN," ");
340 fprintf(PHLUN,"%f %f %f %f ph = ",ph[0],ph[1],ph[2],ph[3]);
341 // fprintf(PHLUN,"") 'p1= ',PHEP(1,1),PHEP(2,1),PHEP(3,1),PHEP(4,1)
342 // fprintf(PHLUN,"") 'p2= ',PHEP(1,2),PHEP(2,2),PHEP(3,2),PHEP(4,2)
343 // fprintf(PHLUN,"") 'p3= ',PHEP(1,3),PHEP(2,3),PHEP(3,3),PHEP(4,3)
344 // fprintf(PHLUN,"") 'p4= ',PHEP(1,4),PHEP(2,4),PHEP(3,4),PHEP(4,4)
345 // fprintf(PHLUN,"") 'p5= ',PHEP(1,5),PHEP(2,5),PHEP(3,5),PHEP(4,5)
346
347 fprintf(PHLUN," c= %f theta= %f",C,th1);
348 // fprintf(PHLUN,"") 'photos waga daje ... IBREM=',IBREM,' waga=',waga
349 // fprintf(PHLUN,"") 'xk,COSTHG,c',xk,COSTHG,c
350 // fprintf(PHLUN,"") svar/4./xkap*(1.D0-COSTHG*BETA),
351 // $ (1-delta) /wt2 *(1.D0+COSTHG*BETA)/2, wagan2
352 // fprintf(PHLUN,"") ' delta, wt2,beta', delta, wt2,beta
353 fprintf(PHLUN," - ");
354 fprintf(PHLUN,"t,u = %f %f",t,u);
355 fprintf(PHLUN,"t1,u1 = %f %f",t1,u1);
356 fprintf(PHLUN,"sredniaki = %f %f",svar*(1-C)/2,svar*(1+C)/2);
357 // ! fprintf(PHLUN,"") 'xk= %f c= %f COSTHG= %f' ,xk,c,COSTHG
358 fprintf(PHLUN,"PHASYZ(svar)=',%f,' svar= %f',' waga= %f",PHASYZ(svar,IDE,IDF),svar,waga);
359 fprintf(PHLUN," - ");
360 fprintf(PHLUN,"BT-part= %f BU-part= %f",
361 2*(BT*t*t+BT*t1*t1)
362 /(1+(1-xk)*(1-xk))* 2.0/(BT*(1-C)*(1-C))/svar/svar,
363 2*(BU*u*u+BU*u1*u1)
364 /(1+(1-xk)*(1-xk))* 2.0/(BU*(1+C)*(1+C))/svar/svar);
365 fprintf(PHLUN,"BT-part*BU-part= %f wagan2= %f",
366 2*(BT*t*t+BT*t1*t1)
367 /(1+(1-xk)*(1-xk))* 2.0/(BT*(1-C)*(1-C))/svar/svar
368 *2*(BU*u*u+BU*u1*u1)
369 /(1+(1-xk)*(1-xk))* 2.0/(BU*(1+C)*(1+C))/svar/svar, wagan2);
370
371 fprintf(PHLUN,"wagan2= %f",wagan2);
372 fprintf(PHLUN," ################### ");
373
374 // wagan2=wagan2*VakPol(qp,qm,ph,pp,pm,IDE,IDF); // default VakWpol returns 1. Hook for the user function
375 wagan2=3.8; // ! overwrite
376 waga=2/(1.0+COSTHG*BETA)*wagan2 ;
377 // % * svar/4./xkap*(1.D0-COSTHG*BETA)*sqrt(1.0-xk)
378
379
380
381
382 return waga;
383
384}
385
386double PhotosMEforZ::VakPol(double qp[4],double qm[4],double ph[4],double pp[4],double pm[4],int IDE,int IDF)
387{
388 return PhotosMEforZ::currentVakPol(qp,qm,ph,pp,pm,IDE,IDF);
389}
390
391double PhotosMEforZ::default_VakPol(double qp[4],double qm[4],double ph[4],double pp[4],double pm[4],int IDE,int IDF)
392{
393 return 1; // that is default.
394}
395
396void PhotosMEforZ::set_VakPol( double (*fun)(double[4], double[4], double[4], double[4], double[4], int, int) )
397{
398 if (fun==NULL) PhotosMEforZ::currentVakPol = default_VakPol;
399 else PhotosMEforZ::currentVakPol = fun;
400}
401
402//----------------------------------------------------------------------
403//
404// PHWTNLO: PHotosWTatNLO
405//
406// Purpose: calculates instead of interference weight
407// complete NLO weight for vector boson decays
408// of pure vector (or pseudovector) couplings
409// Proper orientation of beams required.
410// Uses Zphwtnlo encapsulating actual matrix element
411// At NLO more tuning on kinematical conf.
412// than in standard is needed.
413// Works with KORALZ and KKM//
414// Note some commented out commons from MUSTAAL, KORALZ
415//
416// Input Parameters: Common /PHOEVT/ /PHOPS/ /PHOREST/ /PHOPRO/
417//
418// Output Parameters: Function value
419//
420// Author(s): Z. Was Created at: 08/12/05
421// Last Update: 23/06/13
422//
423//----------------------------------------------------------------------
424
425double PhotosMEforZ::phwtnlo(){
426 // fi3 orientation of photon, fi1,th1 orientation of neutral
427
428 // COMMON/PHOPHS/XPHMAX,XPHOTO,COSTHG,SINTHG
429
430 // COMMON /PHOREST/ FI3,fi1,th1
431 // COMMON /PHWT/ BETA,WT1,WT2,WT3
432 // COMMON/PHOPRO/PROBH,CORWT,XF,IREP
433 // COMMON/PHOMOM/MCHSQR,MNESQR,PNEUTR(5)
434 // static double PI=3.141592653589793238462643;
435 static int i=1;
436 int K,L,IDHEP3,IDUM=0;
437 int IDE,IDF;
438 double QP[4],QM[4],PH[4],QQ[4],PP[4],PM[4],QQS[4];
439 double XK,ENE,svar;
440
441 // REAL*8 s,c,svar,xkaM,xkaP,xk,phwtnlo,xdumm,PHINT
442 // REAL*8 ENE,a,t,u,t1,u1,wagan2,waga,PHASYZ,BT,BU,ENEB
443 // INTEGER IBREM,K,L,IREP,IDUM,IDHEP3
444 // integer icont,ide,idf
445 // REAL*8 delta
446
447/////////////////////
448// phlupa(299500);
449
450
451/////////////////////
452// phlupa(299500);
453
454 XK=2.0*pho.phep[pho.nhep-i][4-i]/pho.phep[1-i][4-i];
455
456// XK=2.0*pho.phep[pho.nhep-i][4-i]/pho.phep[1-i][4-i]/phophs.xphmax; // it is not used becuse here
457 //order of emissions is meaningless
458 if(pho.nhep<=4) XK=0.0;
459 // the mother must be Z or gamma* !!!!
460
461 if (XK>1.0e-10 &&(pho.idhep[1-i]==22 || pho.idhep[1-i]==23)){
462
463 // write(*,*) 'nhep=',nhep
464 // DO K=1,3 ENDDO
465 // IF (K.EQ.1) IBREM= 1
466 // IF (K.EQ.2) IBREM=-1
467 // ICONT=ICONT+1
468 // IBREM=IBREX ! that will be input parameter.
469 // IBREM=IBREY ! that IS now input parameter.
470
471 // We initialize twice 4-vectors, here and again later after boost
472 // must be the same way. Important is how the reduction procedure will work.
473 // It seems at present that the beams must be translated to be back to back.
474 // this may be done after initialising, thus on 4-vectors.
475
476 for( K=1;K<5;K++){
477 PP[K-i]=pho.phep[1-i][K-i];
478 PM[K-i]=pho.phep[2-i][K-i];
479 QP[K-i]=pho.phep[3-i][K-i];
480 QM[K-i]=pho.phep[4-i][K-i];
481 PH[K-i]=pho.phep[pho.nhep-i][K-i];
482 QQ[K-i]=0.0;
483 QQS[K-i]=QP[K-i]+QM[K-i];
484 }
485
486
487 PP[4-i]=(pho.phep[1-i][4-i]+pho.phep[2-i][4-i])/2.0;
488 PM[4-i]=(pho.phep[1-i][4-i]+pho.phep[2-i][4-i])/2.0;
489 PP[3-i]= PP[4-i];
490 PM[3-i]=-PP[4-i];
491
492 for(L=5;L<=pho.nhep-1;L++){
493 for( K=1;K<5;K++){
494 QQ [K-i]=QQ [K-i]+ pho.phep[L-i][K-i];
495 QQS[K-i]=QQS[K-i]+ pho.phep[L-i][K-i];
496 }
497 }
498
499 // go to the restframe of 3
500 PHOB(1,QQS,QP);
501 PHOB(1,QQS,QM);
502 PHOB(1,QQS,QQ);
503 ENE=(QP[4-i]+QM[4-i]+QQ[4-i])/2;
504
505 // preserve direction of emitting particle and wipeout QQ
506 if (phopro.irep==1){
507 double a=sqrt(ENE*ENE-pho.phep[3-i][5-i]*pho.phep[3-i][5-i])/sqrt(QM[4-i]*QM[4-i]-pho.phep[3-i][5-i]*pho.phep[3-i][5-i]);
508 QM[1-i]= QM[1-i]*a;
509 QM[2-i]= QM[2-i]*a;
510 QM[3-i]= QM[3-i]*a;
511 QP[1-i]=-QM[1-i];
512 QP[2-i]=-QM[2-i];
513 QP[3-i]=-QM[3-i];
514 }
515 else{
516 double a=sqrt(ENE*ENE-pho.phep[3-i][5-i]*pho.phep[3-i][5-i])/sqrt(QP[4-i]*QP[4-i]-pho.phep[3-i][5-i]*pho.phep[3-i][5-i]);
517 QP[1-i]= QP[1-i]*a;
518 QP[2-i]= QP[2-i]*a;
519 QP[3-i]= QP[3-i]*a;
520 QM[1-i]=-QP[1-i];
521 QM[2-i]=-QP[2-i];
522 QM[3-i]=-QP[3-i];
523 }
524 QP[4-i]=ENE;
525 QM[4-i]=ENE;
526 // go back to reaction frame (QQ eliminated)
527 PHOB(-1,QQS,QP);
528 PHOB(-1,QQS,QM);
529 PHOB(-1,QQS,QQ);
530
531 svar=pho.phep[1-i][4-i]*pho.phep[1-i][4-i];
532
533 IDE=hep.idhep[1-i];
534 IDF=hep.idhep[4-i];
535 if(abs(hep.idhep[4-i])==abs(hep.idhep[3-i])) IDF=hep.idhep[3-i];
536
537 IDHEP3=pho.idhep[3-i];
538 return Zphwtnlo(svar,XK,IDHEP3,phopro.irep,QP,QM,PH,PP,PM,phophs.costhg,phwt.beta,phorest.th1,IDE,IDF);
539 }
540 else{
541 // in other cases we just use default setups.
542 return PHINT(IDUM);
543 }
544}
545
546} // namespace Photospp
547
Support functions.
static double Zphwtnlo(double svar, double xk, int IDHEP3, int IREP, double qp[4], double qm[4], double ph[4], double pp[4], double pm[4], double COSTHG, double BETA, double th1, int IDE, int IDF)
Definition: forZ-MEc.cxx:249
static double PHBORNM(double svar, double costhe, double T3e, double qe, double T3f, double qf, int Ncf)
Definition: forZ-MEc.cxx:62