TEvtGen/Photos/PhotosUtilities.cxx

   1 #include "PhotosUtilities.h"
   2 #include <cstdlib>
   3 #include <cstdio>
   4 using std::max;
   5
   6 namespace Photospp
   7 {
   8
   9 namespace PhotosUtilities
  10 {
  11
  12
  13 void fill_val(int beg, int end, double* array, double value)
  14 {
  15   for (int i = beg; i < end; i++)
  16     array[i] = value;
  17 }
  18
  19
  20 //----------------------------------------------------------------------
  21 //
  22 //    PHOEPS:   PHOeps vector product (normalized to unity)
  23 //
  24 //    Purpose:  calculates vector product, then normalizes its length.
  25 //              used to generate orthogonal vectors, i.e. to
  26 //              generate polarimetric vectors for photons.
  27 //
  28 //    Input Parameters:  VEC1,VEC2 - input 4-vectors
  29 //
  30 //    Output Parameters: EPS - normalized 4-vector, orthogonal to
  31 //                             VEC1 and VEC2
  32 //
  33 //    Author(s):  Z. Was, P.Golonka               Created at:  19/01/05
  34 //                                                Last Update: 10/06/13
  35 //
  36 //----------------------------------------------------------------------
  37
  38 void PHOEPS(double vec1[4], double vec2[4], double eps[4]){
  39   double xn;
  40   int j=1;  // convention of indices of Riemann space must be preserved.
  41
  42   eps[1-j]=vec1[2-j]*vec2[3-j] - vec1[3-j]*vec2[2-j];
  43   eps[2-j]=vec1[3-j]*vec2[1-j] - vec1[1-j]*vec2[3-j];
  44   eps[3-j]=vec1[1-j]*vec2[2-j] - vec1[2-j]*vec2[1-j];
  45   eps[4-j]=0.0;
  46
  47   xn=sqrt( eps[1-j]*eps[1-j] + eps[2-j]*eps[2-j] + eps[3-j]*eps[3-j]);
  48
  49   eps[1-j]=eps[1-j]/xn;
  50   eps[2-j]=eps[2-j]/xn;
  51   eps[3-j]=eps[3-j]/xn;
  52
  53 }
  54
  55
  56
  57 //----------------------------------------------------------------------
  58 //
  59 //    PHOTOS:   PHOton radiation  in decays function for SPIn determina-
  60 //              tion
  61 //
  62 //    Purpose:  Calculate  the spin  of particle  with  code IDHEP.  The
  63 //              code  of the particle  is  defined  by the Particle Data
  64 //              Group in Phys. Lett. B204 (1988) 1.
  65 //
  66 //    Input Parameter:   IDHEP
  67 //
  68 //    Output Parameter:  Funtion  value = spin  of  particle  with  code
  69 //                       IDHEP
  70 //
  71 //    Author(s):  E. Barberio and B. van Eijk     Created at:  29/11/89
  72 //                                                Last update: 10/06/13
  73 //
  74 //----------------------------------------------------------------------
  75 double PHOSPI(int idhep){
  76   static double SPIN[100] = { 0 };
  77   static int j=0;
  78   //--
  79   //--   Array 'SPIN' contains the spin  of  the first 100 particles accor-
  80   //--   ding to the PDG particle code...
  81
  82   if(j==0) // initialization
  83     {
  84       j=1;
  85       fill_val(0 ,  8, SPIN, 0.5);
  86       fill_val(8 ,  9, SPIN, 1.0);
  87       fill_val(9 , 10, SPIN, 0.0);
  88       fill_val(10, 18, SPIN, 0.5);
  89       fill_val(18, 20, SPIN, 0.0);
  90       fill_val(20, 24, SPIN, 1.0);
  91       fill_val(24,100, SPIN, 0.0);
  92     }
  93
  94   int idabs=abs(idhep);
  95   //--
  96   //--   Spin of quark, lepton, boson etc....
  97   if (idabs-1<100) return SPIN[idabs-1];
  98
  99   //--   ...other particles, however...
 100   double xx=((idabs % 10)-1.0)/2.0;
 101   //--
 102   //--   ...K_short and K_long are special !!
 103   xx=max(xx,0.0);
 104   return xx;
 105 }
 106
 107 //----------------------------------------------------------------------
 108 //
 109 //    PHOTOS:   PHOton radiation in decays CHArge determination
 110 //
 111 //    Purpose:  Calculate the charge  of particle  with code IDHEP.  The
 112 //              code  of the  particle  is  defined by the Particle Data
 113 //              Group in Phys. Lett. B204 (1988) 1.
 114 //
 115 //    Input Parameter:   IDHEP
 116 //
 117 //    Output Parameter:  Funtion value = charge  of  particle  with code
 118 //                       IDHEP
 119 //
 120 //    Author(s):  E. Barberio and B. van Eijk     Created at:  29/11/89
 121 //                                                Last update: 11/06/13
 122 //
 123 //----------------------------------------------------------------------
 124 double PHOCHA(int idhep){
 125   static double CHARGE[101] = { 0 };
 126   static int j=0;
 127   //--
 128   //--   Array 'SPIN' contains the spin  of  the first 100 particles accor-
 129   //--   ding to the PDG particle code...
 130
 131   if(j==0) // initialization
 132     {
 133       j=1;
 134       fill_val(0 ,  1, CHARGE, 0.0         );
 135       fill_val(1 ,  2, CHARGE,-0.3333333333);
 136       fill_val(2 ,  3, CHARGE, 0.6666666667);
 137       fill_val(3 ,  4, CHARGE,-0.3333333333);
 138       fill_val(4 ,  5, CHARGE, 0.6666666667);
 139       fill_val(5 ,  6, CHARGE,-0.3333333333);
 140       fill_val(6 ,  7, CHARGE, 0.6666666667);
 141       fill_val(7 ,  8, CHARGE,-0.3333333333);
 142       fill_val(8 ,  9, CHARGE, 0.6666666667);
 143       fill_val(9 , 11, CHARGE, 0.0         );
 144       fill_val(11 ,12, CHARGE,-1.0         );
 145       fill_val(12 ,13, CHARGE, 0.0         );
 146       fill_val(13 ,14, CHARGE,-1.0         );
 147       fill_val(14, 15, CHARGE, 0.0         );
 148       fill_val(15 ,16, CHARGE,-1.0         );
 149       fill_val(16, 17, CHARGE, 0.0         );
 150       fill_val(17 ,18, CHARGE,-1.0         );
 151       fill_val(18, 24, CHARGE, 0.0         );
 152       fill_val(24, 25, CHARGE, 1.0         );
 153       fill_val(25, 37, CHARGE, 0.0         );
 154       fill_val(37, 38, CHARGE, 1.0         );
 155       fill_val(38,101, CHARGE, 0.0         );
 156     }
 157
 158   int idabs=abs(idhep);
 159   double phoch=0.0;
 160
 161   //--
 162   //--   Charge of quark, lepton, boson etc....
 163   if (idabs<=100) phoch=CHARGE[idabs];
 164   else {
 165     int Q3= idabs/1000 % 10;
 166     int Q2= idabs/100  % 10;
 167     int Q1= idabs/10   % 10;
 168     if (Q3==0){
 169       //--
 170       //-- ...meson...
 171       if(Q2 % 2==0) phoch=CHARGE[Q2]-CHARGE[Q1];
 172       else          phoch=CHARGE[Q1]-CHARGE[Q2];
 173     }
 174     else{
 175       //--
 176       //--   ...diquarks or baryon.
 177       phoch=CHARGE[Q1]+CHARGE[Q2]+CHARGE[Q3];
 178     }
 179   }
 180   //--
 181   //--   Find the sign of the charge...
 182   if (idhep<0.0) phoch=-phoch;
 183   if (phoch*phoch<0.000001) phoch=0.0;
 184
 185   return phoch;
 186 }
 187
 188
 189
 190
 191 //----------------------------------------------------------------------
 192 //
 193 //    PHOTOS:   PHOton radiation in decays calculation of TRIangle fie
 194 //
 195 //    Purpose:  Calculation of triangle function for phase space.
 196 //
 197 //    Input Parameters:  A, B, C (Virtual) particle masses.
 198 //
 199 //    Output Parameter:  Function value =
 200 //                       SQRT(LAMBDA(A**2,B**2,C**2))/(2*A)
 201 //
 202 //    Author(s):  B. van Eijk                     Created at:  15/11/89
 203 //                                                Last Update: 12/06/13
 204 //
 205 //----------------------------------------------------------------------
 206 double PHOTRI(double A,double B,double C){
 207   double DA,DB,DC,DAPB,DAMB,DTRIAN;
 208   DA=A;
 209   DB=B;
 210   DC=C;
 211   DAPB=DA+DB;
 212   DAMB=DA-DB;
 213   DTRIAN=sqrt((DAMB-DC)*(DAPB+DC)*(DAMB+DC)*(DAPB-DC));
 214   return DTRIAN/(DA+DA);
 215 }
 216 //----------------------------------------------------------------------
 217 //
 218 //    PHOTOS:   PHOton radiation in decays calculation of ANgle '1'
 219 //
 220 //    Purpose:  Calculate angle from X and Y
 221 //
 222 //    Input Parameters:  X, Y
 223 //
 224 //    Output Parameter:  Function value
 225 //
 226 //    Author(s):  S. Jadach                       Created at:  01/01/89
 227 //                B. van Eijk                     Last Update: 12/06/13
 228 //
 229 //----------------------------------------------------------------------
 230 double PHOAN1(double X,double Y){
 231
 232   double phoan1 = 0.0;
 233
 234   static double PI=3.14159265358979324, TWOPI=6.28318530717958648;
 235
 236   if (fabs(Y)<fabs(X)){
 237     phoan1=atan(fabs(Y/X));
 238     if (X<0.0) phoan1=PI-phoan1;
 239   }
 240   else phoan1=acos(X/sqrt(X*X+Y*Y));
 241   //
 242   if (Y<0.0) phoan1=TWOPI-phoan1;
 243   return phoan1;
 244
 245 }
 246
 247 //----------------------------------------------------------------------
 248 //
 249 //    PHOTOS:   PHOton radiation in decays calculation of ANgle '2'
 250 //
 251 //    Purpose:  Calculate angle from X and Y
 252 //
 253 //    Input Parameters:  X, Y
 254 //
 255 //    Output Parameter:  Function value
 256 //
 257 //    Author(s):  S. Jadach                       Created at:  01/01/89
 258 //                B. van Eijk                     Last Update: 12/06/13
 259 //
 260 //----------------------------------------------------------------------
 261 double PHOAN2(double X,double Y){
 262
 263   double phoan2 = 0.0;
 264
 265   static double PI=3.14159265358979324; //, TWOPI=6.28318530717958648;
 266
 267   if (fabs(Y)<fabs(X)){
 268     phoan2=atan(fabs(Y/X));
 269     if (X<0.0) phoan2=PI-phoan2;
 270   }
 271   else phoan2=acos(X/sqrt(X*X+Y*Y));
 272   return phoan2;
 273 }
 274
 275 //----------------------------------------------------------------------
 276 //
 277 //    PHOTOS:   PHOton radiation in decays ROtation routine '2'
 278 //
 279 //    Purpose:  Rotate  x and z components  of vector PVEC  around angle
 280 //              'ANGLE'.
 281 //
 282 //    Input Parameters:  ANGLE, PVEC
 283 //
 284 //    Output Parameter:  PVEC
 285 //
 286 //    Author(s):  S. Jadach                       Created at:  01/01/89
 287 //                B. van Eijk                     Last Update: 12/06/13
 288 //
 289 //----------------------------------------------------------------------
 290 void PHORO2(double ANGLE,double PVEC[4]){
 291   int j=1;  // convention of indices of Riemann space must be preserved.
 292
 293   double CS,SN;
 294   CS= cos(ANGLE)*PVEC[1-j]+sin(ANGLE)*PVEC[3-j];
 295   SN=-sin(ANGLE)*PVEC[1-j]+cos(ANGLE)*PVEC[3-j];
 296   PVEC[1-j]=CS;
 297   PVEC[3-j]=SN;
 298 }
 299
 300 //----------------------------------------------------------------------
 301 //
 302 //    PHOTOS:   PHOton radiation in decays ROtation routine '3'
 303 //
 304 //    Purpose:  Rotate  x and y components  of vector PVEC  around angle
 305 //              'ANGLE'.
 306 //
 307 //    Input Parameters:  ANGLE, PVEC
 308 //
 309 //    Output Parameter:  PVEC
 310 //
 311 //    Author(s):  S. Jadach     RO                 Created at:  01/01/89
 312 //                B. van Eijk                     Last Update: 12/06/13
 313 //
 314 //----------------------------------------------------------------------
 315 void PHORO3(double ANGLE,double PVEC[4]){
 316   int j=1;  // convention of indices of Riemann space must be preserved.
 317   double CS,SN;
 318   CS=cos(ANGLE)*PVEC[1-j]-sin(ANGLE)*PVEC[2-j];
 319   SN=sin(ANGLE)*PVEC[1-j]+cos(ANGLE)*PVEC[2-j];
 320   PVEC[1-j]=CS;
 321   PVEC[2-j]=SN;
 322 }
 323
 324 //----------------------------------------------------------------------
 325 //
 326 //
 327 //    PHOB:     PHotosBoost
 328 //
 329 //    Purpose:  Boosts VEC to (MODE=1)  rest frame of PBOOS1;
 330 //              or back (MODE=1)
 331 //
 332 //    Input Parameters:   MODE,PBOOS1,VEC
 333 //
 334 //    Output Parameters:  VEC
 335 //
 336 //    Author(s):                                  Created at:  08/12/05
 337 //                Z. Was                          Last Update: 13/06/13
 338 //
 339 //----------------------------------------------------------------------
 340
 341 void PHOB(int MODE,double PBOOS1[4],double vec[4]){
 342   double BET1[3],GAM1,PB;
 343   static int j0=1;
 344   int J;
 345
 346
 347   PB=sqrt(PBOOS1[4-j0]*PBOOS1[4-j0]-PBOOS1[3-j0]*PBOOS1[3-j0]-PBOOS1[2-j0]*PBOOS1[2-j0]-PBOOS1[1-j0]*PBOOS1[1-j0]);
 348   for( J=1; J<4;J++){
 349     if (MODE==1) BET1[J-j0]=-PBOOS1[J-j0]/PB;
 350     else BET1[J-j0]= PBOOS1[J-j0]/PB;
 351   }
 352
 353   GAM1=PBOOS1[4-j0]/PB;
 354
 355   //--
 356   //--   Boost vector
 357
 358   PB=BET1[1-j0]*vec[1-j0]+BET1[2-j0]*vec[2-j0]+BET1[3-j0]*vec[3-j0];
 359
 360   for( J=1; J<4;J++) vec[J-j0]=vec[J-j0]+BET1[J-j0]*(vec[4-j0]+PB/(GAM1+1.0));
 361   vec[4-j0]=GAM1*vec[4-j0]+PB;
 362   //--
 363 }
 364
 365
 366 //     *******************************
 367 // Boost along arbitrary axis (as implemented by Ronald Kleiss).
 368 // The method is described in book of Bjorken and Drell
 369 // p boosted into r  from actual frame to rest frame of q
 370 // forth (mode = 1) or back (mode = -1).
 371 // q must be a timelike, p may be arbitrary.
 372 void bostdq(int mode,double qq[4],double pp[4],double r[4]){
 373   double q[4],p[4],amq,fac;
 374   static int i=1;
 375   int k;
 376
 377   for(k=1;k<=4;k++){
 378     p[k-i]=pp[k-i];
 379     q[k-i]=qq[k-i];
 380   }
 381   amq =sqrt(q[4-i]*q[4-i]-q[1-i]*q[1-i]-q[2-i]*q[2-i]-q[3-i]*q[3-i]);
 382
 383   if    (mode == -1){
 384     r[4-i] = (p[1-i]*q[1-i]+p[2-i]*q[2-i]+p[3-i]*q[3-i]+p[4-i]*q[4-i])/amq;
 385     fac  = (r[4-i]+p[4-i])/(q[4-i]+amq);
 386   }
 387   else if(mode ==  1){
 388     r[4-i] =(-p[1-i]*q[1-i]-p[2-i]*q[2-i]-p[3-i]*q[3-i]+p[4-i]*q[4-i])/amq;
 389     fac  =-(r[4-i]+p[4-i])/(q[4-i]+amq);
 390   }
 391   else{
 392     cout << " ++++++++ wrong mode in boostdq " << endl;
 393     exit(0);
 394   }
 395   r[1-i]=p[1-i]+fac*q[1-i];
 396   r[2-i]=p[2-i]+fac*q[2-i];
 397   r[3-i]=p[3-i]+fac*q[3-i];
 398 }
 399
 400
 401 //----------------------------------------------------------------------
 402 //
 403 //    PHOTOS:   PHOton radiation in decays BOost routine '3'
 404 //
 405 //    Purpose:  Boost  vector PVEC  along z-axis where ANGLE = EXP(ETA),
 406 //              ETA is the hyperbolic velocity.
 407 //
 408 //    Input Parameters:  ANGLE, PVEC
 409 //
 410 //    Output Parameter:  PVEC
 411 //
 412 //    Author(s):  S. Jadach                       Created at:  01/01/89
 413 //                B. van Eijk                     Last Update: 12/06/13
 414 //
 415 //----------------------------------------------------------------------
 416 void PHOBO3(double ANGLE,double PVEC[4]){
 417   int j=1;  // convention of indices of Riemann space must be preserved.
 418   double QPL,QMI;
 419   QPL=(PVEC[4-j]+PVEC[3-j])*ANGLE;
 420   QMI=(PVEC[4-j]-PVEC[3-j])/ANGLE;
 421   PVEC[3-j]=(QPL-QMI)/2.0;
 422   PVEC[4-j]=(QPL+QMI)/2.0;
 423 }
 424
 425 } // namespace PhotosUtilities
 426
 427 } // namespace Photospp
 428