PhotosUtilities.cxx
1#include "PhotosUtilities.h"
2#include <cstdlib>
3#include <cstdio>
4using std::max;
5
6namespace Photospp
7{
8
9namespace PhotosUtilities
10{
11
12
13void 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
38void 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//----------------------------------------------------------------------
75double 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//----------------------------------------------------------------------
124double 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//----------------------------------------------------------------------
206double 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//----------------------------------------------------------------------
230double 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//----------------------------------------------------------------------
261double 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//----------------------------------------------------------------------
290void 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//----------------------------------------------------------------------
315void 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
341void 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.
372void 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(-1);
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//----------------------------------------------------------------------
416void 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
Support functions.