StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
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(-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 //----------------------------------------------------------------------
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 
PhotosUtilities
Support functions.
Generated by   doxygen 1.8.5