TEvtGen/Photos/photosC.cxx

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