source: trunk/MagicSoft/Simulation/Detector/ReflectorII/ph2cph.c@ 801

Last change on this file since 801 was 787, checked in by magicsol, 24 years ago
A new parameter has been ontroduced to call FindLagrange to avoid acces to no allocated memory.
File size: 19.2 KB
Line 
1#include <stdio.h>
2#include <math.h>
3#include "diag.h"
4#include "init.h"
5#include "lagrange.h"
6/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
7
8/* random numbers */
9#define RandomNumber ranf()
10
11/* Speed of Light in vacuum, in m/s */
12#define Speed_of_Light_vacuum 299792458.0f
13#define Speed_of_Light_air (Speed_of_Light_vacuum / 1.000293f)
14
15/* Speed of Light in vacuum, in cm/ns */
16#define Speed_of_Light_vacuum_cmns (Speed_of_Light_vacuum / 1.0e7f)
17#define Speed_of_Light_air_cmns (Speed_of_Light_air / 1.0e7f)
18
19/* Macros */
20#define SQR(A) ((A)*(A))
21#define NORM(A) ((float) sqrt((SQR(A[0]))+(SQR(A[1]))+(SQR(A[2]))))
22
23/* Function declarations */
24extern float ranf(void);
25void rnormal(double *r, int n);
26
27/* Static definitions */
28float OmegaCT[3][3];
29float OmegaICT[3][3];
30float Omega[3][3];
31float OmegaI[3][3];
32
33static double NormalRandomNumbers[500];
34/*
35 From photons on ground, i.e. observation level,
36 to photons on focal plane, i.e. chamber !
37
38 Mirror reflectivity
39 Mirror reflection
40 Photon position on chamber
41 Position smearing
42 Timing
43
44 Returned values:
45
46 0 OK photon reached the chamber
47
48 1 Photon lost due to mirror reflectivity
49 2 Photon lost because out of mirror
50 3 Photon lost due to black spot
51 4 Photon lost because reflected out of chamber
52
53*/
54
55int ph2cph(photon *ph, cphoton *cph)
56{
57
58 float u, v, w; /* photon director cosines */
59
60 float r[3]; /* photon trajectory */
61 float x[3]; /* position of the photon on ground */
62
63 float rCT[3]; /* photon trajectory in the system of the CT */
64 float xCT[3]; /* photon position on ground (CT) */
65 float rm[3]; /* photon trajectory in the system of a mirror */
66 float xmm[3]; /* intermediate values */
67 float xm[3]; /* photon position on ground */
68 float xcut[3]; /* location of the cut sphere-trajectory */
69 float xcutCT[3]; /* location of the cut sphere-trajectory (CT) */
70
71 float rnor[3], rnorm; /* normal in that point */
72
73 float rrefl[3]; /* reflected vector, from the mirror */
74 float rreflCT[3]; /* reflected vector, from the CT */
75 float xcam[3]; /* where the photon hits the camera plane */
76
77 float calpha; /* cos(alpha=angle incident/normal) */
78 float phi; /* angle between photon and camera plane */
79
80 float a, b, c, t; /* intermediate variables */
81
82 float d; /* minimum distance trajectory-mirror center */
83
84 float wl; /* photon wavelength */
85 float reflec; /* reflectivity for a photon */
86
87 float h; /* photon production height */
88
89 int i, k; /* simple counters */
90 int i_mirror=-1; /* number of a given mirror */
91
92
93 float distmirr, distmirr2; /* distances used in MAGIC reflection routine */
94 float sx, sy;
95 float t1, t2;
96
97 void makeOmega(float theta, float phi);
98 void makeOmegaI(float theta, float phi);
99 void applyMxV(float M[3][3], float *V, float *Vp);
100 float Lin2Curv(float x);
101
102 /* begin code */
103
104 /* get photon wawelength */
105
106 wl = ph->w;
107
108 /* get position on ground */
109
110 x[0] = ph->x;
111 x[1] = ph->y;
112 x[2] = 0.0; /* ground => obs. level => z=0 */
113
114 /* get director cosines x,y on ground */
115
116 r[0] = ph->u;
117 r[1] = ph->v;
118 r[2] = (float) sqrt(1.0 - r[0]*r[0] - r[1]*r[1]);
119
120 /* get photon time and production height */
121
122 h = ph->h;
123
124 /*!@'
125
126 @#### Reflectivity of the mirrors.
127
128 We make a 3rd. order interpolation using Lagrange
129 polynomials, in order to calculate the reflectivity of the
130 mirror for that wavelength. Further developments will
131 include also a incidence-angle dependence (which is not very
132 important).
133
134 */
135
136 /* ++
137 FILTER: REFLECTIVITY R(lambda)
138 -- */
139
140 /* find data point to be used in Lagrange interpolation (-> k) */
141
142 FindLagrange(Reflectivity,k,wl,nReflectivity);
143
144 /* if random > reflectivity then goes to the TOP of the loop again */
145
146 reflec = Lagrange(Reflectivity,k,wl);
147
148 if ( RandomNumber > reflec ) return 1;
149
150
151 /*!@'
152
153 @#### Reflection on mirrors.
154 We calculate reflected photon direction
155
156 */
157
158 /* ++
159 REFLECTION
160 -- */
161
162 Debug("@1 %f %f %f %f %f %f\n", x[0], x[1], x[2], r[0], r[1], r[2]);
163
164 /* change to the system of the CT */
165
166 applyMxV( OmegaCT, x, xCT );
167 applyMxV( OmegaCT, r, rCT );
168
169 /*
170 before moving to the system of the mirror
171 we look whether the photon hits a mirror or not
172
173 calculate the intersection of the trajectory of the photon
174 with the GLOBAL DISH !!!
175 we reproduce the calculation of the coefficients of the
176 second order polynomial in z (=xCT[2]), made with
177 Mathematica
178 */
179
180 /*
181 * In[1]:= parab:=z-(x^2+y^2)/(4F)
182 * par1=parab /. {x->x0+u/w(z-z0),y->y0+v/w(z-z0)}
183 *
184 * Out[1]=
185 * u (z - z0) 2 v (z - z0) 2
186 * (x0 + ----------) + (y0 + ----------)
187 * w w
188 * z - ---------------------------------------
189 * 4 F
190 *
191 * In[2]:= CoefficientList[ExpandAll[par1*4F*w^2],z]
192 *
193 * Out[2]=
194 * 2 2 2 2
195 * {-(w x0 ) - w y0 + 2 u w x0 z0 +
196 *
197 * 2 2 2 2
198 * 2 v w y0 z0 - u z0 - v z0 ,
199 *
200 * 2 2
201 * 4 F w - 2 u w x0 - 2 v w y0 + 2 u z0 +
202 *
203 * 2 2 2
204 * 2 v z0, -u - v }
205 */
206
207 /* the z coordinate is calculated */
208
209 a = - SQR(rCT[0]) - SQR(rCT[1]);
210
211 b = (float) (4.0*ct_Focal_mean*SQR(rCT[2])
212 - 2.0*rCT[0]*rCT[2]*xCT[0] - 2.0*rCT[1]*rCT[2]*xCT[1]
213 + 2.0*SQR(rCT[0])*xCT[2] + 2.0*SQR(rCT[1])*xCT[2]);
214
215 c = 2*rCT[0]*rCT[2]*x[0]*x[2] + 2*rCT[1]*rCT[2]*x[1]*x[2]
216 - SQR(rCT[2])*SQR(x[0]) - SQR(rCT[2])*SQR(x[1])
217 - SQR(rCT[0])*SQR(x[2]) - SQR(rCT[1])*SQR(x[2]);
218
219 if ( fabs(a) < 1.e-6 ) {
220
221 /* only one value */
222
223 xcut[2] = -c / b;
224
225 } else {
226
227 d = (float) sqrt( b*b - 4.0*a*c );
228
229 /* two possible values for z */
230
231 t1 = (float) ((-b+d) / (2.0*a));
232 t2 = (float) ((-b-d) / (2.0*a));
233
234 /* z must be the minimum of t1 and t2 */
235
236 xcut[2] = (t1 < t2) ? t1 : t2;
237
238 }
239
240 /*
241 xcut[] is NOW the cut between the GLOBAL dish of MAGIC and
242 the trajectory of the photon
243 */
244
245 xcut[0] = xCT[0] + rCT[0]/rCT[2]*(xcut[2]-xCT[2]);
246 xcut[1] = xCT[1] + rCT[1]/rCT[2]*(xcut[2]-xCT[2]);
247
248 /* convert to Curvilinear distance over the parabolic dish */
249
250 sx = Lin2Curv( xcut[0] );
251 sy = Lin2Curv( xcut[1] );
252
253 /* is it outside the dish? */
254
255 if ((fabs(sx) > 850.0) || (fabs(sy) > 850.0)) {
256 /*
257 cout << "CONDITION 1 !" << endl << flush;
258 cout << '1';
259 */
260 return 2;
261 }
262
263 /* calculate the mirror to be used */
264
265 distmirr = 1000000.0f;
266
267 for (i=0; i<ct_NMirrors && distmirr>=ct_RMirror; ++i) {
268 distmirr2 = (float) sqrt(SQR(ct_data[i].x - xcut[0]) +
269 SQR(ct_data[i].y - xcut[1]) +
270 SQR(ct_data[i].z - xcut[2]));
271 if (distmirr2 < distmirr) {
272 i_mirror = i;
273 distmirr = distmirr2;
274 }
275 }
276
277 /*
278 the mirror to use is i_mirror (calculated several lines above)
279 check whether the photon is outside the nearest (this) mirror
280 */
281
282 if ((fabs(ct_data[i_mirror].sx - sx) > ct_RMirror) ||
283 (fabs(ct_data[i_mirror].sy - sy) > ct_RMirror)) {
284 /*
285 cout << "CONDITION 2 !" << endl << flush;
286 cout << '2';
287 */
288 return 2;
289 }
290
291 /* calculate matrices for the mirror */
292
293 makeOmega (-ct_data[i_mirror].theta,
294 ct_data[i_mirror].phi);
295 makeOmegaI(-ct_data[i_mirror].theta,
296 ct_data[i_mirror].phi);
297
298 /* change to the system of the mirror */
299
300 /* first translation... */
301
302 xmm[0] = xCT[0] - ct_data[i_mirror].x;
303 xmm[1] = xCT[1] - ct_data[i_mirror].y;
304 xmm[2] = xCT[2] - ct_data[i_mirror].z;
305
306 /* ...then rotation */
307
308 applyMxV( Omega, xmm, xm );
309 applyMxV( Omega, rCT, rm );
310
311 /*
312 the vector rCT should be normalized, and
313 so the vector rm remains normalized as well, but, anyhow...
314 */
315
316 rnorm = NORM( rm );
317 rm[0] /= rnorm;
318 rm[1] /= rnorm;
319 rm[2] /= rnorm;
320
321 /*
322 calculate the intersection of the trajectory of the photon
323 with the mirror
324 we reproduce the calculation of the coefficients of the
325 second order polynomial in z (=xm[2]), made with
326 Mathematica
327 */
328
329 /*
330 * In[1]:= esfera:=x^2+y^2+(z-R)^2-R^2;
331 * recta:={x->x0+u/w(z-z0),y->y0+v/w(z-z0)}
332 *
333 * In[2]:= esfera
334 *
335 * 2 2 2 2
336 * Out[2]= -R + x + y + (-R + z)
337 *
338 * In[3]:= recta
339 *
340 * u (z - z0) v (z - z0)
341 * Out[3]= {x -> x0 + ----------, y -> y0 + ----------}
342 * w w
343 *
344 * In[4]:= esf=esfera /. recta
345 *
346 * 2 2 u (z - z0) 2 v (z - z0) 2
347 * Out[4]= -R + (-R + z) + (x0 + ----------) + (y0 + ----------)
348 * w w
349 *
350 * In[5]:= coefs=CoefficientList[ExpandAll[esf],z]
351 *
352 * 2 2 2 2
353 * 2 2 2 u x0 z0 2 v y0 z0 u z0 v z0
354 * Out[5]= {x0 + y0 - --------- - --------- + ------ + ------,
355 * w w 2 2
356 * w w
357 *
358 * 2 2 2 2
359 * 2 u x0 2 v y0 2 u z0 2 v z0 u v
360 * > -2 R + ------ + ------ - ------- - -------, 1 + -- + --}
361 * w w 2 2 2 2
362 * w w w w
363 * In[6]:= Simplify[ExpandAll[coefs*w^2]]
364 *
365 * 2 2 2 2 2 2
366 * Out[6]= {w (x0 + y0 ) - 2 w (u x0 + v y0) z0 + (u + v ) z0 ,
367 *
368 * 2 2 2 2 2
369 * > -2 (R w - u w x0 + u z0 + v (-(w y0) + v z0)), u + v + w }
370 *
371 */
372
373 /*
374 the z coordinate is calculated, using the coefficient
375 shown above
376 */
377
378 a = SQR(rm[0]) + SQR(rm[1]) + SQR(rm[2]);
379
380 b = (float) (-2*(2.*ct_data[i_mirror].f*SQR(rm[2])
381 - rm[0]*rm[2]*xm[0]
382 + SQR(rm[0])*xm[2]
383 + rm[1]*(-(rm[2]*xm[1]) + rm[1]*xm[2])));
384
385 c = (SQR(rm[2])*(SQR(xm[0]) + SQR(xm[1]))
386 - 2*rm[2]*(rm[0]*xm[0] + rm[1]*xm[1])*xm[2]
387 + (SQR(rm[0]) + SQR(rm[1]))*SQR(xm[2]));
388
389 d = (float) sqrt( b*b - 4.0*a*c );
390
391 /* two possible values for z */
392
393 t1 = (float) ((-b+d) / (2.0*a));
394 t2 = (float) ((-b-d) / (2.0*a));
395
396 /* z must be the minimum of t1 and t2 */
397
398 xcut[2] = (t1 < t2) ? t1 : t2;
399 xcut[0] = xm[0] + rm[0]/rm[2]*(xcut[2]-xm[2]);
400 xcut[1] = xm[1] + rm[1]/rm[2]*(xcut[2]-xm[2]);
401
402 /*
403 ++
404 BLACK SPOTS: If the photon hits the black spot, it's lost
405 --
406 */
407
408 if ( sqrt(SQR(xcut[0]) + SQR(xcut[1])) < ct_BlackSpot_rad ) {
409 /*
410 cout << "CONDITION 3!\n" << flush;
411 cout << '3';
412 */
413 return 3;
414 }
415
416 /*
417 if we still have the photon, we continue with the reflexion;
418 we calculate normal vector in this point
419 (and normalize, with the sign changed)
420 */
421
422 rnor[0] = 2.0f*xcut[0];
423 rnor[1] = 2.0f*xcut[1];
424 rnor[2] = (float) (2.0*(xcut[2] - 2.0*ct_Focal[i_mirror]));
425
426 rnorm = -NORM( rnor );
427 rnor[0] /= rnorm;
428 rnor[1] /= rnorm;
429 rnor[2] /= rnorm;
430
431 /*
432 now, both "normal" vector and original trajectory are
433 normalized
434 just project the original vector in the normal, and
435 take it as the "mean" position of the original and
436 the "reflected" vector
437 from this, we can calculate the "reflected" vector
438 calpha = cos(angle(rnor,rm))
439 */
440
441 calpha = (float) fabs(rnor[0]*rm[0] + rnor[1]*rm[1] + rnor[2]*rm[2]);
442
443 /* finally!!! we have the reflected trajectory of the photon */
444
445 rrefl[0] = (float) (2.0*rnor[0]*calpha - rm[0]);
446 rrefl[1] = (float) (2.0*rnor[1]*calpha - rm[1]);
447 rrefl[2] = (float) (2.0*rnor[2]*calpha - rm[2]);
448
449 rnorm = NORM( rrefl );
450 rrefl[0] /= rnorm;
451 rrefl[1] /= rnorm;
452 rrefl[2] /= rnorm;
453
454 /* let's go back to the coordinate system of the CT */
455
456 /* first rotation... */
457
458 applyMxV( OmegaI, xcut, xcutCT);
459 applyMxV( OmegaI, rrefl, rreflCT);
460
461 /* ...then translation */
462
463 xcutCT[0] += ct_data[i_mirror].x;
464 xcutCT[1] += ct_data[i_mirror].y;
465 xcutCT[2] += ct_data[i_mirror].z;
466
467 /*
468 calculate intersection of this trajectory and the camera plane
469 in the system of the CT, this plane is z = ct_Focal
470 */
471
472 t = (ct_Focal_mean - xcutCT[2]) / rreflCT[2];
473
474 xcam[0] = xcutCT[0] + rreflCT[0]*t;
475 xcam[1] = xcutCT[1] + rreflCT[1]*t;
476 xcam[2] = xcutCT[2] + rreflCT[2]*t;
477
478 /*
479 ++
480 AXIS DEVIATION: We introduce it here just as a first order
481 correction, by modifying the position of the reflected photon.
482 --
483 */
484
485 xcam[0] += AxisDeviation[0][i_mirror];
486 xcam[1] += AxisDeviation[1][i_mirror];
487
488 /*
489 ++
490 SMEARING: We apply the point spread function for the mirrors
491 --
492 */
493
494 /* get two N(0;1) random numbers */
495
496 rnormal( NormalRandomNumbers, 2 );
497
498 /* modify the Cphoton position in the camera */
499
500 xcam[0] += (float) (NormalRandomNumbers[0] * ct_PSpread_mean);
501 xcam[1] += (float) (NormalRandomNumbers[1] * ct_PSpread_mean);
502
503 /* check whether the photon goes out of the camera */
504
505 if ( (SQR(xcam[0])+SQR(xcam[1])) > SQR(ct_CameraWidth) ) {
506 return 4;
507 }
508
509 /*
510 ++
511 ANGLE OF INCIDENCE
512 --
513
514 calculate angle of incidence between tray. and camera plane
515 the camera plane is
516 0 y + 0 y + z - ct_Focal = 0 => (A,B,C,D) = (0,0,1,-ct_Focal)
517 from Table 3.20 "Tasch. der Math."
518 */
519
520 phi = (float) asin(rreflCT[2]);
521
522 /*
523 ++
524 TIMING
525 --
526 */
527
528 /* calculate the new time of the photon (in the camera) */
529
530 t = ph->t;
531
532 /*
533 substract path from the mirror till the ground, 'cos
534 the photon actually hit the mirror!!
535 */
536
537 t = (float) (t + ((( xm[2] > 0. ) ? -1.0 : +1.0) *
538 sqrt( SQR(xm[0] - xcut[0]) +
539 SQR(xm[1] - xcut[1]) +
540 SQR(xm[2] - xcut[2]) ) / Speed_of_Light_air_cmns));
541
542 /* add path from the mirror till the camera */
543
544 t = (float) (t + sqrt( SQR(xcutCT[0] - xcam[0]) +
545 SQR(xcutCT[1] - xcam[1]) +
546 SQR(xcutCT[2] - xcam[2]) ) / Speed_of_Light_air_cmns);
547
548 /* show it */
549
550 Debug("@2 %f %f %f %f %f %f %f %f %f\n"
551 "@3 %f %f %d %f %f %f %f\n"
552 "@4 %f %f %f %f\n",
553 xCT[0], xCT[1], xCT[2], rCT[0], rCT[1], rCT[2],
554 xcut[0], xcut[1], xcut[2],
555 sx, sy, i_mirror, ct_data[i_mirror].sx, ct_data[i_mirror].sy,
556 ct_data[i_mirror].sx - sx, ct_data[i_mirror].sy - sy,
557 xcam[0], xcam[1], xcam[2], phi);
558
559 /* Output */
560
561/* cph->w = wl; */
562 cph->x = xcam[0];
563 cph->y = xcam[1];
564 cph->u = r[0];
565 cph->v = r[1];
566 cph->t = t;
567 cph->h = h;
568 cph->phi = phi;
569
570 return 0;
571
572} /* end of ph2cph */
573
574
575/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
576
577!---------------------------------------------------------------------
578 @name makeOmega
579
580 @desc function to calculate the matrix Omega(theta,phi)
581
582 @var theta Angle theta of the transformation
583 @var phi Angle phi of the transformation
584
585 @date Sat Jun 27 05:58:56 MET DST 1998
586----------------------------------------------------------------------
587 @function
588*/
589
590void
591makeOmega (float theta, float phi)
592{
593 static float ct, st, cp, sp;
594
595 /* shortcuts for cosine and sine of theta and phi */
596 ct = (float) cos(theta);
597 st = (float) sin(theta);
598 cp = (float) cos(phi);
599 sp = (float) sin(phi);
600
601 /* save values in the array (see top of file) */
602 Omega[0][0] = cp*ct;
603 Omega[0][1] = sp*ct;
604 Omega[0][2] = -st;
605
606 Omega[1][0] = -sp;
607 Omega[1][1] = cp;
608 Omega[1][2] = 0;
609
610 Omega[2][0] = cp*st;
611 Omega[2][1] = sp*st;
612 Omega[2][2] = ct;
613}
614
615
616/*
617!---------------------------------------------------------------------
618 @name makeOmegaI
619
620 @desc function to calculate the matrix Omega-1(theta,phi)
621
622 @var theta Angle theta of the transformation
623 @var phi Angle phi of the transformation
624
625 @date Sat Jun 27 05:58:56 MET DST 1998
626----------------------------------------------------------------------
627 @function
628*/
629
630void
631makeOmegaI(float theta, float phi)
632{
633 static float ct, st, cp, sp;
634
635 /* shortcuts for cosine and sine of theta and phi */
636 ct = (float) cos(theta);
637 st = (float) sin(theta);
638 cp = (float) cos(phi);
639 sp = (float) sin(phi);
640
641 /* save values in the array (see top of file) */
642 OmegaI[0][0] = cp*ct;
643 OmegaI[0][1] = -sp;
644 OmegaI[0][2] = cp*st;
645
646 OmegaI[1][0] = sp*ct;
647 OmegaI[1][1] = cp;
648 OmegaI[1][2] = sp*st;
649
650 OmegaI[2][0] = -st;
651 OmegaI[2][1] = 0;
652 OmegaI[2][2] = ct;
653}
654
655
656/*
657!---------------------------------------------------------------------
658 @name applyMxv
659
660 @desc returns the vector v' such that v' = M x v
661
662 @var M matrix of the transformation
663 @var v vector to be multiplied
664 @var vi resulting vector
665
666 @date Sat Jun 27 05:58:56 MET DST 1998
667----------------------------------------------------------------------
668 @function
669*/
670
671void
672applyMxV(float M[3][3], float *V, float *Vp)
673{
674 Vp[0] = (M[0][0] * V[0] +
675 M[0][1] * V[1] +
676 M[0][2] * V[2]);
677 Vp[1] = (M[1][0] * V[0] +
678 M[1][1] * V[1] +
679 M[1][2] * V[2]);
680 Vp[2] = (M[2][0] * V[0] +
681 M[2][1] * V[1] +
682 M[2][2] * V[2]);
683}
684
685/*
686!---------------------------------------------------------------------
687 @name Lin2Curv
688
689 @desc Linear (Euclidean) to Curvilinear distance
690
691 @var x Radial distance from the axis of the paraboloid
692
693 @return Curvilinear distance over the parabolic shape
694
695 @date Wed Jul 8 15:25:39 MET DST 1998
696----------------------------------------------------------------------
697 @function
698*/
699
700float
701Lin2Curv(float x)
702{
703 x /= 100.f;
704 return ((x + 0.000144175317185f * x * x * x)*100.f);
705}
706
707/*!---------------------------------------------------------------------
708// @name rnormal
709//
710// @desc returns n(=2k) normaly distributed numbers
711//
712// @var *r pointer to a vector where we write the numbers
713// @var n how many numbers do we generate
714//
715// @date Sat Jun 27 05:58:56 MET DST 1998
716//----------------------------------------------------------------------
717// @function */
718
719void rnormal(double *r, int n)
720{
721
722 double z1, z2;
723 int i;
724
725 for (i=0; i<n; i+=2) {
726
727 z1 = RandomNumber;
728 z2 = RandomNumber;
729
730 r[i] = sqrt(-2.0*log(z1)) * cos(2.0*M_PI*z2);
731 r[i+1] = sqrt(-2.0*log(z1)) * sin(2.0*M_PI*z2);
732
733 }
734
735}
Note: See TracBrowser for help on using the repository browser.