Line data Source code
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 0 : for (int i = beg; i < end; i++)
16 0 : array[i] = value;
17 0 : }
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 0 : eps[1-j]=vec1[2-j]*vec2[3-j] - vec1[3-j]*vec2[2-j];
43 0 : eps[2-j]=vec1[3-j]*vec2[1-j] - vec1[1-j]*vec2[3-j];
44 0 : eps[3-j]=vec1[1-j]*vec2[2-j] - vec1[2-j]*vec2[1-j];
45 0 : eps[4-j]=0.0;
46 :
47 0 : xn=sqrt( eps[1-j]*eps[1-j] + eps[2-j]*eps[2-j] + eps[3-j]*eps[3-j]);
48 :
49 0 : eps[1-j]=eps[1-j]/xn;
50 0 : eps[2-j]=eps[2-j]/xn;
51 0 : eps[3-j]=eps[3-j]/xn;
52 :
53 0 : }
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 0 : if(j==0) // initialization
83 : {
84 0 : j=1;
85 0 : fill_val(0 , 8, SPIN, 0.5);
86 0 : fill_val(8 , 9, SPIN, 1.0);
87 0 : fill_val(9 , 10, SPIN, 0.0);
88 0 : fill_val(10, 18, SPIN, 0.5);
89 0 : fill_val(18, 20, SPIN, 0.0);
90 0 : fill_val(20, 24, SPIN, 1.0);
91 0 : fill_val(24,100, SPIN, 0.0);
92 0 : }
93 :
94 0 : int idabs=abs(idhep);
95 : //--
96 : //-- Spin of quark, lepton, boson etc....
97 0 : if (idabs-1<100) return SPIN[idabs-1];
98 :
99 : //-- ...other particles, however...
100 0 : double xx=((idabs % 10)-1.0)/2.0;
101 : //--
102 : //-- ...K_short and K_long are special !!
103 0 : xx=max(xx,0.0);
104 : return xx;
105 0 : }
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 0 : if(j==0) // initialization
132 : {
133 0 : j=1;
134 0 : fill_val(0 , 1, CHARGE, 0.0 );
135 0 : fill_val(1 , 2, CHARGE,-0.3333333333);
136 0 : fill_val(2 , 3, CHARGE, 0.6666666667);
137 0 : fill_val(3 , 4, CHARGE,-0.3333333333);
138 0 : fill_val(4 , 5, CHARGE, 0.6666666667);
139 0 : fill_val(5 , 6, CHARGE,-0.3333333333);
140 0 : fill_val(6 , 7, CHARGE, 0.6666666667);
141 0 : fill_val(7 , 8, CHARGE,-0.3333333333);
142 0 : fill_val(8 , 9, CHARGE, 0.6666666667);
143 0 : fill_val(9 , 11, CHARGE, 0.0 );
144 0 : fill_val(11 ,12, CHARGE,-1.0 );
145 0 : fill_val(12 ,13, CHARGE, 0.0 );
146 0 : fill_val(13 ,14, CHARGE,-1.0 );
147 0 : fill_val(14, 15, CHARGE, 0.0 );
148 0 : fill_val(15 ,16, CHARGE,-1.0 );
149 0 : fill_val(16, 17, CHARGE, 0.0 );
150 0 : fill_val(17 ,18, CHARGE,-1.0 );
151 0 : fill_val(18, 24, CHARGE, 0.0 );
152 0 : fill_val(24, 25, CHARGE, 1.0 );
153 0 : fill_val(25, 37, CHARGE, 0.0 );
154 0 : fill_val(37, 38, CHARGE, 1.0 );
155 0 : fill_val(38,101, CHARGE, 0.0 );
156 0 : }
157 :
158 0 : int idabs=abs(idhep);
159 : double phoch=0.0;
160 :
161 : //--
162 : //-- Charge of quark, lepton, boson etc....
163 0 : if (idabs<=100) phoch=CHARGE[idabs];
164 : else {
165 0 : int Q3= idabs/1000 % 10;
166 0 : int Q2= idabs/100 % 10;
167 0 : int Q1= idabs/10 % 10;
168 0 : if (Q3==0){
169 : //--
170 : //-- ...meson...
171 0 : if(Q2 % 2==0) phoch=CHARGE[Q2]-CHARGE[Q1];
172 0 : else phoch=CHARGE[Q1]-CHARGE[Q2];
173 : }
174 : else{
175 : //--
176 : //-- ...diquarks or baryon.
177 0 : phoch=CHARGE[Q1]+CHARGE[Q2]+CHARGE[Q3];
178 : }
179 : }
180 : //--
181 : //-- Find the sign of the charge...
182 0 : if (idhep<0.0) phoch=-phoch;
183 0 : if (phoch*phoch<0.000001) phoch=0.0;
184 :
185 0 : 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 0 : DAPB=DA+DB;
212 0 : DAMB=DA-DB;
213 0 : DTRIAN=sqrt((DAMB-DC)*(DAPB+DC)*(DAMB+DC)*(DAPB-DC));
214 0 : 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 0 : if (fabs(Y)<fabs(X)){
237 0 : phoan1=atan(fabs(Y/X));
238 0 : if (X<0.0) phoan1=PI-phoan1;
239 : }
240 0 : else phoan1=acos(X/sqrt(X*X+Y*Y));
241 : //
242 0 : if (Y<0.0) phoan1=TWOPI-phoan1;
243 0 : 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 0 : if (fabs(Y)<fabs(X)){
268 0 : phoan2=atan(fabs(Y/X));
269 0 : if (X<0.0) phoan2=PI-phoan2;
270 : }
271 0 : else phoan2=acos(X/sqrt(X*X+Y*Y));
272 0 : 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 0 : CS= cos(ANGLE)*PVEC[1-j]+sin(ANGLE)*PVEC[3-j];
295 0 : SN=-sin(ANGLE)*PVEC[1-j]+cos(ANGLE)*PVEC[3-j];
296 0 : PVEC[1-j]=CS;
297 0 : PVEC[3-j]=SN;
298 0 : }
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 0 : CS=cos(ANGLE)*PVEC[1-j]-sin(ANGLE)*PVEC[2-j];
319 0 : SN=sin(ANGLE)*PVEC[1-j]+cos(ANGLE)*PVEC[2-j];
320 0 : PVEC[1-j]=CS;
321 0 : PVEC[2-j]=SN;
322 0 : }
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 0 : double BET1[3],GAM1,PB;
343 : static int j0=1;
344 : int J;
345 :
346 :
347 0 : 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 0 : for( J=1; J<4;J++){
349 0 : if (MODE==1) BET1[J-j0]=-PBOOS1[J-j0]/PB;
350 0 : else BET1[J-j0]= PBOOS1[J-j0]/PB;
351 : }
352 :
353 0 : GAM1=PBOOS1[4-j0]/PB;
354 :
355 : //--
356 : //-- Boost vector
357 :
358 0 : PB=BET1[1-j0]*vec[1-j0]+BET1[2-j0]*vec[2-j0]+BET1[3-j0]*vec[3-j0];
359 :
360 0 : for( J=1; J<4;J++) vec[J-j0]=vec[J-j0]+BET1[J-j0]*(vec[4-j0]+PB/(GAM1+1.0));
361 0 : vec[4-j0]=GAM1*vec[4-j0]+PB;
362 : //--
363 0 : }
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 0 : double q[4],p[4],amq,fac;
374 : static int i=1;
375 : int k;
376 :
377 0 : for(k=1;k<=4;k++){
378 0 : p[k-i]=pp[k-i];
379 0 : q[k-i]=qq[k-i];
380 : }
381 0 : 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 0 : if (mode == -1){
384 0 : 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 0 : fac = (r[4-i]+p[4-i])/(q[4-i]+amq);
386 0 : }
387 0 : else if(mode == 1){
388 0 : 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 0 : fac =-(r[4-i]+p[4-i])/(q[4-i]+amq);
390 : }
391 : else{
392 0 : cout << " ++++++++ wrong mode in boostdq " << endl;
393 0 : exit(-1);
394 : }
395 0 : r[1-i]=p[1-i]+fac*q[1-i];
396 0 : r[2-i]=p[2-i]+fac*q[2-i];
397 0 : r[3-i]=p[3-i]+fac*q[3-i];
398 0 : }
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 0 : QPL=(PVEC[4-j]+PVEC[3-j])*ANGLE;
420 0 : QMI=(PVEC[4-j]-PVEC[3-j])/ANGLE;
421 0 : PVEC[3-j]=(QPL-QMI)/2.0;
422 0 : PVEC[4-j]=(QPL+QMI)/2.0;
423 0 : }
424 :
425 : } // namespace PhotosUtilities
426 :
427 : } // namespace Photospp
428 :
|
