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