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