source: trunk/MagicSoft/Simulation/Detector/Reflector/reflector.cxx@ 339

Last change on this file since 339 was 339, checked in by blanch, 25 years ago
Previous version modified by blanch are wrong. Now the change in the theta distribution is right.
File size: 95.8 KB
Line 
1//!/////////////////////////////////////////////////////////////////////
2//
3// reflector
4//
5// @file reflector.cxx
6// @title Simulation of the Reflector
7// @desc Program for the simulation of CT1 and MAGIC reflectors
8// @author J C Gonzalez
9// @email gonzalez@mppmu.mpg.de
10// @date Thu May 7 16:24:22 1998
11//
12//----------------------------------------------------------------------
13//
14// Created: Thu May 7 16:24:22 1998
15// Author: Jose Carlos Gonzalez
16// Purpose: Program for reflector simulation
17// Notes:
18//
19//----------------------------------------------------------------------
20//
21// $RCSfile: reflector.cxx,v $
22// $Revision: 1.6 $
23// $Author: blanch $
24// $Date: 2000-01-14 16:39:50 $
25//
26////////////////////////////////////////////////////////////////////////
27// @tableofcontents @coverpage
28
29//= Below you can find a sort of, incomplete for sure, documentation
30//= concerning this program. Please, do not hesiate to contact the
31//= author in case of problems.
32
33//=---------------------------------------------------------------------
34//!@section Details of the simulation of the detectors MAGIC and CT1.
35
36//=---------------------------------------------------------------------
37//!@subsection Overview of the simulation.
38
39/*!@"
40
41 This short section gives an idea about how the simulation
42 developes itself. The details on each step will be explained
43 later. At the beginning we have the Generated output file.
44
45 @itemize
46
47 @- B@Atmospheric Absorption@: A detailed description of the
48 atmospheric absorption is provided. This description takes into
49 account the effects of the Rayleigh scattering, the aerosol (Mie)
50 absorption, and the effect of the Ozone. It works for the wavelength
51 range @$290 \div 600@$ nm.
52
53 @- B@Reflection in the mirrors@: Each remaining photon is traced
54 till the mirrors, and then reflected till the camera plane. The
55 steps are the following:
56
57 @enumerate
58
59 @- Calculation of the mirror where the photon is
60 hitting.
61
62 @- We apply the E@reflectivity@ effect, using a continuous random
63 number in the range [0:1] : the photon can be lost.
64
65 @- If the photon is hitting in the E@black spot@ of the center,
66 the photon is also lost.
67
68 @- We take the E@mathematicaly exact@ normal of the mirror, and
69 apply the E@axis-deviation@ effect: the normal changes a little
70 bit.
71
72 @- We apply the E@imperfections@ in the mirrors (measurements are
73 needed).
74
75 @- Finally, with the displaced normal, we calculate the position of
76 the reflected photon in the camera plane. Here, even with an exact,
77 mathematical ray-tracing the image is E@blurred@, due to the
78 intrinsic behaviour of the telescope optics.
79
80 @endenumerate
81
82 At this step we have a B@RAW-Monte-Carlo@ output file.
83
84 @- B@Collection in the PMTs/IPCs@: The collection in the active
85 devices follows the next steps:
86
87 @enumerate
88
89 @- We apply a certain probability of transmision of the photon in
90 the E@plexiglas@, dependent on the incident angle of the photon in
91 the camera (measurements are needed).
92
93 @- We apply the transmitance factor of the E@light guides@.
94
95 \item Then we apply the E@Quantum efficiency@ of the PMT,
96 depending on the incident angle.
97
98 @endenumerate
99
100 @- B@Correction of the obtained signal@ In every pixel we have a
101 certain number of Cherenkov photons, what we take as a B@mean value@
102 of the measured charge in that pixel. Together to that value, we
103 take into account the E@Night Sky Background@ (NSB): we have for
104 this a mean value measured, and with this mean value we calculate
105 the contribution of the sky for each pixel.
106
107 At this level we will have a B@RAW-Data@ output file.
108
109 @enditemize
110
111 @"*/
112
113//=---------------------------------------------------------------------
114//!@subsection Simulation of the detector.
115
116//=---------------------------------------------------------------------
117//!@subsubsection Coordinate systems.
118
119/*!@"
120
121 The different angles used in the simulation of the detector are
122 shown in the next figures. The telescope (CT) is pointing to the
123 direction @$(\theta,\phi)@$, while the particle (background) is coming
124 from the direction @$(\theta_p,\phi_p)@$. This direction is calculated
125 using a pair of random generated angles @$\theta_s@$ and @$\phi_s@$, and
126 spherical trigonometry.
127
128 @F
129 \begin{figure}[htbp]
130 \begin{center}
131 \includegraphics[height=0.45\textheight]{coordinates.eps}
132 \caption{Coordinate system used in the reflector simulation}
133 \label{fig:coord}
134 \end{center}
135 \end{figure}
136 @F
137
138 @F
139 \stepcounter{figure}
140 \begin{figure}[htbp]
141 \begin{center}
142 \includegraphics[height=0.45\textheight]{reflector.eps}
143 \caption{Angles used in the reflector simulation}
144 \label{fig:reflector}
145 \end{center}
146 \end{figure}
147 @F
148
149 The angles are defined in the next table. Using this definitions, we
150 can start to explain how the reflector program works.
151
152 @F
153 \begin{table}
154 \begin{tabular}[hbt]{cl}
155 $\theta$ & Zenith angle of the position where the CT points to\\
156 $\phi$ & Horizontal angle from the North of the position\\
157 & where the CT points to\\
158 $\theta_s$ & Angular distance between the axis of the\\
159 & incoming shower and the direction where the CT is \\
160 & pointing to.\\
161 $\phi_s$ & Polar angle with axis equal to the axis of the CT.\\
162 $\theta_p$ & Zenith angle of the incoming shower.\\
163 $\phi_p$ & Horizontal angle (from the North) of the\\
164 & direction of the incoming shower.\\
165 $\theta_{CT}$ & Zenith angle of the dish ($= \theta$).\\
166 $\phi_{CT}$ & Horizontal angle of the dish ($= \phi$).\\
167 $\theta_m$ & Angle of the mirror axis relative to the CT axis.\\
168 $\phi_m$ & Polar angle of the position of the mirror
169 (arbitrary origin).\\
170 \end{tabular}
171 \caption{Definition of the different angles in the reflector program.}
172 \end{table}
173 @F
174
175 At the CORSIKA level, and also at the reflector simulation level, we
176 calculate the angles @$\theta_p@$ and @$\phi_p@$, which give the
177 direction of the incoming shower. To calculate them we use the
178 following.
179
180 @F
181 \begin{eqnarray}
182 \cos\theta_p &=& \cos\theta \cos\theta_s
183 + \sin\theta \sin\theta_s \cos\phi_s
184 \\
185 \sin\phi_p &=& \frac{\sin\theta_s}{\sin\theta_p} \sin\phi_s
186 \end{eqnarray}
187 @F
188
189 To finish with the definitions, we just write down some of
190 the parameters we will use in the reflexion algorithm, in
191 the next table.
192
193 @F
194 \stepcounter{table}
195 \begin{table}
196 \begin{tabular}[hbt]{cl}
197 $\mathbf{x}=(x,y,z)$ &
198 Position of the Cherenkov photon on the ground\\
199 & in the global system OXYZ.\\
200
201 $\mathbf{r}=(u,v,w)$ &
202 Vector of the Cherenkov photon's trayectory\\
203 & in the global system OXYZ.\\
204
205 $\mathbf{x_{CT}}=(x_{CT},y_{CT},z_{CT})$ &
206 Position of the Cherenkov photon on the ground\\
207 & in the system OXYZ$_{CT}$ of the CT.\\
208
209 $\mathbf{r_{CT}}=(u_{CT},v_{CT},w_{CT})$ &
210 Vector of the Cherenkov photon's trayectory\\
211 & in the system OXYZ$_{CT}$ of the CT.\\
212
213 $\mathbf{x_{m}}=(x_{m},y_{m},z_{m})$ &
214 Position of the Cherenkov photon on the ground\\
215 & in the system OXYZ$_{m}$ of the mirror.\\
216
217 $\mathbf{r_{m}}=(u_{m},v_{m},w_{m})$ &
218 Vector of the Cherenkov photon's trayectory\\
219 & in the system OXYZ$_{m}$ of the mirror.\\
220
221 $\mathbf{x_{cut}}=(x_{cut},y_{cut},z_{cut})$ &
222 Cut of the trayectory of the photon with the mirror,\\
223 & in the system OXYZ$_{m}$ of the mirror.\\
224
225 $\mathbf{r_\perp}=(u_\perp,v_\perp,w_\perp)$ &
226 Vector perpendicular to the mirror in the point $\mathbf{x_{cut}}$\\
227 & in the system OXYZ$_{m}$ of the mirror.\\
228
229 $\mathbf{r_r}=(u_r,v_r,w_r)$ &
230 Vector of the reflected photon, \\
231 & in the system OXYZ$_{m}$ of the mirror.\\
232
233 $\mathbf{r_r^{CT}}=(u_r^{CT},v_r^{CT},w_r^{CT})$ &
234 Vector of the reflected photon, \\
235 & in the system OXYZ$_{CT}$ of the CT.\\
236
237 $\mathbf{x_{cam}}=(x_{cam},y_{cam},z_{cam})$ &
238 Position of the photon in the camera plane,\\
239 & in the system OXYZ$_{CT}$ of the CT.\\
240
241 \end{tabular}
242 \caption{Different vectors and points used in the reflector
243 program.}
244 \end{table}
245 @F
246
247 @"*/
248
249//=---------------------------------------------------------------------
250//!@subsubsection Algorithm of the reflexion.
251
252/*!@"
253
254 We start calculating the matrices for changing of system of
255 coordinates, @$\Omega(\theta,\phi)@$ and
256 @$\Omega^{-1}(\theta,\phi)@$. They are defined to be:
257
258 @F
259 \begin{equation}
260 \begin{split}
261 \Omega(\theta,\phi) &= R_y(-\theta) R_z(-\phi) \\
262 \Omega^{-1}(\theta,\phi) &= R_z(\phi) R_y(\theta)
263 \end{split}
264 \end{equation}
265 @F
266
267 where the matrices @$R_x@$, @$R_y@$ and @$R_z@$ are rotations of
268 angle @$\alpha@$ around each axis OX, OY and OZ.
269
270 @F
271 \begin{eqnarray}
272 R_x(\alpha) =
273 \begin{bmatrix}
274 1 & 0 & 0 \\
275 0 & \cos\alpha & -\sin\alpha \\
276 0 & \sin\alpha & \cos\alpha \\
277 \end{bmatrix}
278 \label{eq:rotx}
279 \\
280 R_y(\alpha) =
281 \begin{bmatrix}
282 \cos\alpha & 0 & \sin\alpha \\
283 0 & 1 & 0 \\
284 -\sin\alpha & 0 & \cos\alpha \\
285 \end{bmatrix}
286 \label{eq:roty}
287 \\
288 R_z(\alpha) =
289 \begin{bmatrix}
290 \cos\alpha & -\sin\alpha & 0\\
291 \sin\alpha & \cos\alpha & 0\\
292 0 & 0 & 1 \\
293 \end{bmatrix}
294 \label{eq:rotz}
295 \end{eqnarray}
296 @F
297
298 With this, the matrices @$\Omega(\theta,\phi)@$ and
299 @$\Omega^{-1}(\theta,\phi)@$ result to be:
300
301 @F
302 \begin{eqnarray}
303 \Omega(\theta,\phi) =
304 \begin{bmatrix}
305 \cos\phi \cos\theta&\sin\phi \cos\theta&-\sin\theta \\
306 -\sin\phi & \cos\phi & 0 \\
307 \cos\phi \sin\theta&\sin\phi \sin\theta& \cos\theta \\
308 \end{bmatrix}
309 \\
310 \Omega^{-1}(\theta,\phi) =
311 \begin{bmatrix}
312 \cos\phi \cos\theta&-\sin\phi&\cos\phi \sin\theta\\
313 \sin\phi \cos\theta& \cos\phi&\sin\phi \sin\theta\\
314 -\sin\theta & 0 & \cos\theta \\
315 \end{bmatrix}
316 \end{eqnarray}
317 @F
318
319 and each change of system of coordinates uses these
320 matrices. In order to calculate the new coordinates X' of a
321 point X in the new system O', and viceversa, we apply the
322 following equations:
323
324 @F
325 \begin{equation}
326 \begin{split}
327 X' &= \Omega(\theta,\phi) X\\
328 X &= \Omega^{-1}(\theta,\phi) X'\\
329 \end{split}
330 \end{equation}
331 @F
332
333 The steps taken in the reflector simulation are the following:
334
335 @enumerate
336
337 @- We look for the mirror @$m@$ which is closer to the
338 trayectory of the photon (presumably, the photon will hit
339 it).
340
341 @- We move to the system @$\mbox{OXYZ}_{CT}@$ of the CT.
342
343 @- We move again, now to the system @$\mbox{OXYZ}_{m}@$ of the
344 mirror @$m@$.
345
346 @- Then the point where the trayectory of the photon intersects with
347 the mirror is calculated.
348
349 @- In this point, the normal (perpendicular) to the mirror is
350 obtained.
351
352 @- The reflected trayectory is calculated with the normal in the
353 mirror and the original trayectory.
354
355 @- We move back to the system @$\mbox{OXYZ}_{CT}@$, and calculate
356 there the point where the new trayectory of the photon hits the
357 camera plane.
358
359 @endenumerate
360
361 In between the previous points the effects of the atmospheric
362 absorption and those mentioned in the next section are included. At
363 the end, the position of the photon and the incident angle in the
364 camera plane is saved.
365
366 @"*/
367
368//=---------------------------------------------------------------------
369//!@subsection Effects to be included in the simulation of the reflector.
370
371/*!@"
372
373 @itemize
374
375 @- B@Incidence in the mirrors@:
376 For each single photon we trace the point where it hits a
377 given mirror. We can use either hexagonal mirrors or
378 circular mirrors with the same surface of the hexagonal ones.
379
380 @- B@Reflectivity of the mirrors@:
381 The reflectivity of the mirrors is a function @$R(\lambda)@$
382 of the wavelength @$\lambda@$ of the photon.
383
384 @- B@Black spots@:
385 In the center of each mirror there is a ``bad'' region,
386 where no reflection can be done at all. The radii of these
387 regions vary, but they can be @$r \simeq 2 \mathrm{\,cm}@$.
388
389 @- B@Imperfections@:
390 The surface of each mirror is locally not-perfect. We need
391 measurements to include this effect.
392
393 @- B@Axis--deviation@:
394 The axis of each mirror can be deviated from the
395 ``teorical'' point in the center of the camera. (@$\sigma
396 \simeq 2.5 \mathrm{\,mm}@$).
397
398 @- B@Blurring@:
399 Pure effect of optics: each mirror is ``seen'' by the camera
400 with a different angle.
401
402 @enditemize
403
404 @"*/
405
406//=---------------------------------------------------------------------
407//!@subsection Effects to be included in the simulation of the camera.
408
409/*!@"
410
411 @itemize
412
413 @- B@Plexiglas@:
414 In front of the light guides there is a ``plexiglas'' layer,
415 with a given transmitance (dependent on the incidence angle).
416
417 @- B@Light Guides@:
418 The looses in the transmitance in the light guides are
419 around @$15-20@$\%.
420
421 @- B@Quantum Efficiency@:
422 The QE of each PMT depends on the impact point of the photon
423 in the window of the PMT and the angle of incidence in this
424 point.
425
426 @- B@Excess noise factor@:
427 It produces a smearing of the signal with a sigma
428 @$\sqrt{2N}@$ instead of the typical @$\sqrt{N}@$.
429
430 @- B@Pick-up noise?@:
431
432 @- B@First dynode collection efficiency?@:
433
434 @enditemize
435
436 @"*/
437
438//=---------------------------------------------------------------------
439// @T \newpage
440
441//=---------------------------------------------------------------------
442//!@section Source code of |reflector.cxx|.
443
444/*!@"
445
446 In this section we show the (commented) code of the program
447 for the simulation of the reflector, in the current version.
448
449 @"*/
450
451//=---------------------------------------------------------------------
452//!@subsection Preliminary notes in the code.
453
454//!@{
455//=================================================================
456//
457// NOTE: Most of the equations used in this program can be found
458// in Bronstein, Semendjajew, Musiol and Muehlig, "Taschenbuch
459// der Mathematik", Verlag Harri Deutsch, Thun und Frankfurt
460// am Main, 1997. Some others were "deduced".
461//
462//
463// Matrices for rotation of coordinate systems
464//---------------------------------------------
465//
466// Omega(theta,phi) = Ry(theta) Rz(phi)
467// Omega-1(theta,phi) = Rz(-phi) Ry(-theta)
468//
469// X` = Omega(theta,phi) X
470// X = Omega-1(theta,phi) X`
471//
472// Omega(theta,phi) =
473// / \
474// | cos(phi)cos(theta) sin(phi)cos(theta) -sin(theta) |
475// | -sin(phi) cos(phi) 0 |
476// | cos(phi)sin(theta) sin(phi)sin(theta) cos(theta) |
477// \ /
478//
479// OmegaI(theta,phi) =
480// / \
481// | cos(phi)cos(theta) -sin(phi) cos(phi)sin(theta) |
482// | sin(phi)cos(theta) cos(phi) sin(phi)sin(theta) |
483// | -sin(theta) 0 cos(theta) |
484// \ /
485//
486//=================================================================
487//!@}
488
489//=---------------------------------------------------------------------
490//!@subsection Includes and Global variables definition.
491
492/*!@"
493
494 All the defines are located in the file |reflector.h|.
495
496 @"*/
497
498//!@{
499#include "reflector.h"
500//!@}
501
502//=---------------------------------------------------------------------
503//!@subsubsection Definition of global variables.
504
505/*!@"
506
507 We first define the matrices to change of coordinates system. The
508 definition of the matrices was already seen in the previous section.
509 All matrices are in the form @$m[i][j]@$ with i the running index
510 for rows, and j the one for columns.
511
512 @"*/
513
514//!@{
515
516//@: matrices to change to the system where the optical axis is OZ
517static float OmegaCT[3][3];
518
519//@: matrices to change to the system where the optical axis is OZ (inverse)
520static float OmegaICT[3][3];
521
522//@: matrices to change the system of coordinates
523static float Omega[3][3];
524
525//@: matrices to change the system of coordinates (inverse)
526static float OmegaI[3][3];
527
528//!@}
529
530/*!@"
531
532 Now we define some global variables with data about the telescope,
533 such as E@focal distance@, E@number of pixels/mirrors@, E@size of
534 the camera@, and so on.
535
536 @"*/
537
538/*!@"
539
540 Depending on the telescope we are using (CT1 or MAGIC), the
541 information stored in the definition file is different. The
542 variable |ct_Type| has the value 0 when we use CT1, and 1 when we
543 use MAGIC.
544
545 @"*/
546
547//!@{
548
549//@: Type of telescope (CT1: 0; MAGIC: 1)
550static int ct_Type;
551
552//!@}
553
554/*!@"
555
556 And this is the information about the whole telescope.
557
558@"*/
559
560//!@{
561// parameters of the CT (from the CT definition file)
562
563//@: Focal distances [cm]
564static float *ct_Focal;
565
566//@: Mean Focal distances [cm]
567static float ct_Focal_mean;
568
569//@: STDev. Focal distances [cm]
570static float ct_Focal_std;
571
572//@: Mean Point Spread function [cm]
573static float ct_PSpread_mean;
574
575//@: STDev. Point Spread function [cm]
576static float ct_PSpread_std;
577
578//@: STDev. Adjustmente deviation [cm]
579static float ct_Adjustment_std;
580
581//@: Radius of the Black Spot in mirror [cm]
582static float ct_BlackSpot_rad;
583
584//@: Radius of one mirror [cm]
585static float ct_RMirror;
586
587//@: Camera width [cm]
588static float ct_CameraWidth;
589
590//@: Pixel width [cm]
591static float ct_PixelWidth;
592
593//@: Number of mirrors
594static int ct_NMirrors = 0;
595
596//@: Number of pixels
597static int ct_NPixels;
598
599//!@}
600
601/*!@"
602
603 The following double-pointer is a 2-dimensional table with
604 information about each mirror in the dish. The routine
605 |read_ct_file()| will read this information from the file with the
606 name given by the user in the E@parameters file@, in the command
607 |ct_file|. The information stored in this file (and in this table)
608 depends on the type of telescope we are using.
609
610 @"*/
611
612//!@{
613
614//@: Pointer to a table with the following info.:
615static float **ct_data;
616
617/*
618 * TYPE=0 (CT1)
619 * i s rho theta x y z thetan phin xn yn zn
620 *
621 * i : number of the mirror
622 * s : arc length [cm]
623 * rho : polar rho of the position of the center of the mirror [cm]
624 * theta : polar angle of the position of the center of the mirror [cm]
625 * x : x coordinate of the center of the mirror [cm]
626 * y : y coordinate of the center of the mirror [cm]
627 * z : z coordinate of the center of the mirror [cm]
628 * thetan : polar theta angle of the direction where the mirror points to
629 * phin : polar phi angle of the direction where the mirror points to
630 * xn : xn coordinate of the normal vector in the center (normalized)
631 * yn : yn coordinate of the normal vector in the center (normalized)
632 * zn : zn coordinate of the normal vector in the center (normalized)
633 *
634 * TYPE=1 (MAGIC)
635 * i f sx sy x y z thetan phin
636 *
637 * i : number of the mirror
638 * f : focal distance of that mirror
639 * sx : curvilinear coordinate of mirror's center in X[cm]
640 * sy : curvilinear coordinate of mirror's center in X[cm]
641 * x : x coordinate of the center of the mirror [cm]
642 * y : y coordinate of the center of the mirror [cm]
643 * z : z coordinate of the center of the mirror [cm]
644 * thetan : polar theta angle of the direction where the mirror points to
645 * phin : polar phi angle of the direction where the mirror points to
646 * xn : xn coordinate of the normal vector in the center (normalized)
647 * yn : yn coordinate of the normal vector in the center (normalized)
648 * zn : zn coordinate of the normal vector in the center (normalized)
649 */
650
651//!@}
652
653/*!@"
654
655 Here we define the global variables where the parameters of a
656 parallel beam of light will be stored.
657
658 @"*/
659
660//!@{
661
662static int pb_ParallelBeam = FALSE;
663static int pb_ParallelBeamPM = FALSE;
664static float pb_Theta, pb_Phi, pb_X, pb_Y;
665static float pb_LengthX, pb_LengthY, pb_NX, pb_NY;
666static float pb_Scale;
667static float pb_Height;
668static char pb_Filename[40];
669
670//!@}
671
672/*!@"
673
674 We can analyze the data stored in the directories provided via the
675 |data_paths| entry in the parameters file in one go, or in blocks of
676 size |Block|, if this variable is a positive integer. This can be
677 done through a shell script. The mechanism is the following:
678
679 @itemize
680
681 @- First, we read the parameters file, and wait for a file called
682 |__DOIT| (written by the script when a new block of |Block| files
683 is available in a loop.
684
685 @- If this file exists, we remove it and create another file with
686 the name |__WORKING| (in order to let the script to know that we
687 are working), and process a block of |Block| files. After that, we
688 will rename the file |__WORKING| to |__READY|.
689
690 @- When the script sees a new block of |Block| files, it will write
691 again a file called |__DOIT|, unless there exist a file called
692 |__WORKING| (which means that the |reflector| is still working in
693 the previous block).
694
695 @- Each time a block is processed, is customary to the user to
696 delete or not that block of files (or to save it to tape, ...).
697
698 @- At the end, a file called |__DONE| is written telling us the
699 number of blocks processed.
700
701 @enditemize
702
703 @"*/
704
705//!@{
706
707static int Block=0;
708int WORKING=FALSE;
709int nfilesBlock=0;
710ofstream obfile;
711long posdir;
712
713//!@}
714
715/*!@"
716
717 We can read also the data from the standard input, and we could
718 write data into the standard output.
719
720 @"*/
721
722//!@{
723
724//@: boolean variable to control whether we read from STDIN or not
725static int Data_From_STDIN = FALSE;
726
727//@: boolean variable to control whether we write to STDOUT or not
728static int Data_To_STDOUT = FALSE;
729
730//!@}
731
732/*!@"
733
734 Next, we have two tables with data got from two files,
735 |reflectivity.dat| and |axisdev.dat|, which contain information
736 about the reflectivity of the mirrors, as a function of the
737 wavelength @$\lambda@$, and about the (simulated) deviation of the
738 mirrors' axes, with respect to the mathematical exact axes.
739
740 @"*/
741
742//!@{
743
744// table of reflectivity for each mirror
745
746//@: table with datapoints (wavelength,reflec.)
747static float **Reflectivity;
748
749//@: number of datapoints
750static int nReflectivity;
751
752//@: table with deviations of the mirrors' normals
753static float **AxisDeviation;
754
755//!@}
756
757/*!@"
758
759 We still define a table into where normal random numbers will be
760 stored by the routine |rnormal(double *r, int n)|.
761
762 @"*/
763
764//!@{
765
766//@: table of normal random numbers
767static double NormalRandomNumbers[500];
768
769//!@}
770
771/*!@"
772
773 This is a flag to change the E@verbosity@ of the output
774
775 @"*/
776
777//!@{
778
779//@: flag to change the verbosity
780int verbose = VERBOSE_DEFAULT;
781
782//!@}
783
784/*!@"
785
786 This option makes the telescope to point to a random position
787 relative to the shower axis, with a maximum angular separation of
788 |RANDOM_POINTING_MAX_SEPARATION| (defined in |reflector.h|.
789
790 @"*/
791
792//!@{
793
794//@: random pointing for the CT?
795int Random_Pointing = FALSE;
796
797//@: maximum random pointing distance
798float Random_Pointing_MaxDist; // [radians]
799
800//!@}
801
802//=---------------------------------------------------------------------
803//!@subsection Main program.
804
805/*!@"
806
807 In the main program we basically follow the items expressed in the
808 first section. Please, have a look into previous sections. I tried
809 to make a code already well commented. Most things should be clear.
810 Nevertheless, I'll try to comment again everything in following
811 releases.
812
813 @"*/
814
815//!@{
816
817//++++++++++++++++++++++++++++++++++++++++
818// MAIN PROGRAM
819//----------------------------------------
820
821int
822main(int argc, char **argv)
823{
824
825 DIR *directory; //@< directory of data
826 struct dirent *de; //@< directory entry
827
828 char pathname[256]; //@< name of directory
829 char cername[256]; //@< name of Cherenkov file
830 char staname[256]; //@< name of Statistics file
831 char outname[256]; //@< name of the output file
832 char parname[256]; //@< name of parameters file
833
834 ifstream cerfile; //@< stream for Cherenkov file
835 ifstream pmfile; //@< stream for pixmap file
836 ofstream outputfile; //@< stream for output file
837 ofstream pbfile; //@< stream for parallel beam file
838
839 COREventHeader evth; //@< Event Header class (from CORSIKA)
840 CORParticle photon; //@< Particle (photon) class (from CORSIKA)
841 CORStatfile stat; //@< Statistics file class (from CORSIKA)
842
843 MCEventHeader mcevth; //@< Event Header class
844 MCCphoton cphoton; //@< Particle (photon) class
845
846 float thetaCT, phiCT, xiCT; //@< parameters of a given shower
847 float coreD, coreX, coreY; //@< core position and distance
848 float u, v, w; //@< cosenos directores
849
850 float r[3]; //@< trayectory
851 float x[3]; //@< position of the photon in the ground
852
853 float rCT[3]; //@< trayectory in the system of the CT
854 float xCT[3]; //@< position of the photon in the ground (CT)
855 float rm[3]; //@< trayectory in the system of a mirror
856 float xmm[3]; //@< intermediate values
857 float xm[3]; //@< position of photon in ground
858 float xcut[3]; //@< location of the cut sphere-trayectory
859 float xcutCT[3]; //@< location of the cut sphere-trayectory (CT)
860
861 float rnor[3], rnorm; //@< normal in that point
862
863 float rrefl[3]; //@< reflected vector, from the mirror
864 float rreflCT[3]; //@< reflected vector, from the CT
865 float xcam[3]; //@< where the photon hits the camera plane
866
867 float calpha; //@< cos(alpha=angle incident/normal)
868 float phi; //@< angle between photon and camera plane
869
870 float a, b, c, t; //@< intermediate variables
871
872 float d; //@< minimum distance trayectory-mirror center
873
874 float wl; //@< wavelength of the Cphoton
875 float reflec; //@< reflectivity for a Cphoton
876
877 int ncer; //@< number of the shower;
878 int np; //@< number of path
879 int i, k; //@< simple counters
880 int i_mirror=-1; //@< number of a given mirror
881
882 int nCphotons; //@< number of a Cphotons written on output file
883
884 float TR; //@< atmospheric transmitance
885 int nbeforeTR, nafterTR; //@< number of Cph. before and after transmission
886 int num_cer_files; //@< number of cer files;
887 int max_num_cer_files; //@< maximum number of cer files to read
888 int f_ev, l_ev; //@< first and last events of a range
889
890 float lE, uE; //@< lower and upper edges of energy range allowed
891
892 float distmirr, distmirr2; //@< distances used in MAGIC reflexion routine
893 float sx, sy;
894 float t1, t2;
895 float pbx, pby, pbt, pbp;
896
897 char line[200];
898
899 int ch, errflg = 0; //@< used by getopt
900
901 /*!@'
902
903 We start in the main program. First we (could) make some
904 presentation, and follows the reading of the parameters file (now
905 from the STDIN), the reading of the CT parameters file, and the
906 creation of the output file, where the processed data will be
907 stored.
908
909 */
910
911 //++
912 // START
913 //--
914
915 // make unbuffered output
916
917 cout.setf ( ios::stdio );
918
919 // parse command line options (see reflector.h)
920
921 parname[0] = '\0';
922
923 optarg = NULL;
924 while ( !errflg && ((ch = getopt(argc, argv, COMMAND_LINE_OPTIONS)) != -1) )
925 switch (ch) {
926 case 'f':
927 strcpy(parname, optarg);
928 break;
929 case 'h':
930 usage();
931 break;
932 default :
933 errflg++;
934 }
935
936 // show help if error
937
938 if ( errflg>0 )
939 usage();
940
941 // make some sort of presentation
942 // present();
943
944 // read parameters from the stdin
945
946 if ( strlen(parname) < 1 )
947 readparam(NULL);
948 else
949 readparam(parname);
950
951 // blocking-mode?
952
953 Block = get_block();
954
955 // reading data from STDIN?
956
957 Data_From_STDIN = get_data_from_stdin();
958
959 // writing data to STDOUT?
960
961 Data_To_STDOUT = get_data_to_stdout();
962
963 // make it verbosely?
964
965 verbose = get_verbose();
966
967 /*!@'
968
969 @#### Parallel beam of light.
970
971 In order to test the behaviour of the program, we can feed it with
972 an artificial data file. This data file will be generated by the
973 program itself and it will consist in a grid of photons, forming
974 all together a parallel beam of light with certain
975 characteristics. These characteristics will be given by the user
976 in the |parameters file|, in the following form.
977
978 parallel\_beam |<theta> <phi> <Xcen> <Ycen> <lenX> <lenY> <nX> <nY> <H>|
979
980 or
981
982 parallel\_beam\_pm |<filename> <scale> <H>|
983
984 In the first form, the first two values give the arrival direction
985 of the beam (zenith angle and azimuth) in degrees, and the last
986 six values specify the geometry of the beam, in cm, and the height
987 of production, also in cm. In the second form, we use a pixmap
988 file in PPM ASCII format, and scale it to create (for 0 degrees of
989 zenith angle) a parallel beam of light.
990
991 */
992
993 //++
994 // PARALLEL BEAM OF LIGHT
995 //--
996
997 // are we going to use a parallel beam of light?
998
999 pb_ParallelBeam = is_parallel_beam();
1000 pb_ParallelBeamPM = is_parallel_beam_pm();
1001
1002 if ( pb_ParallelBeam == TRUE ) {
1003
1004 if ( pb_ParallelBeamPM == FALSE ) {
1005
1006 // get parameters of the parallel beam of light
1007
1008 get_parallel_beam(&pb_Theta, &pb_Phi,
1009 &pb_X, &pb_Y,
1010 &pb_LengthX, &pb_LengthY,
1011 &pb_NX, &pb_NY, &pb_Height);
1012
1013 // generate artificial file(s)
1014
1015 // a. cerfile
1016
1017 sprintf( cername, "%s/cer999999", get_path_name(0) );
1018
1019 pbfile.open( cername, ios::out );
1020
1021 evth.fill ( 1., 1., 1000., 1., 1.e6, 0., 0., 0., pb_Theta, pb_Phi,
1022 pb_X, pb_Y);
1023
1024 evth.write( pbfile );
1025
1026 for ( pbx = (pb_X - pb_LengthX/2.0);
1027 pbx <= (pb_X + pb_LengthX/2.0);
1028 pbx += (pb_LengthX/pb_NX) )
1029
1030 for ( pby = (pb_Y - pb_LengthY/2.0);
1031 pby <= (pb_Y + pb_LengthY/2.0);
1032 pby += (pb_LengthY/pb_NY) ) {
1033
1034 pbt = -atan2(sqrt(SQR(pbx)+SQR(pby)),pb_Height);
1035 pbp = atan2(pby, pbx);
1036 photon.fill(450.0,
1037 pbx, pby,
1038 -sin(pbt) * cos(pbp),
1039 -sin(pbt) * sin(pbp),
1040 0., pb_Height);
1041
1042 photon.write( pbfile );
1043
1044 }
1045
1046 photon.fill( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 );
1047
1048 photon.write( pbfile );
1049
1050 pbfile.close();
1051
1052 // b. statfile
1053
1054 sprintf( staname, "%s/sta999999", get_path_name(0) );
1055
1056 pbfile.open( staname );
1057
1058 stat.write( pbfile );
1059
1060 pbfile.close();
1061
1062 } else {
1063
1064 // get parameters of the parallel beam of light from PIXMAP
1065
1066 strcpy(pb_Filename, get_parallel_beam_pm(&pb_Scale,&pb_Height));
1067
1068 pmfile.open( pb_Filename, ios::out );
1069 pmfile.getline( line, 300 );
1070
1071 do {
1072 pmfile.getline( line, 300 );
1073 } while (line[0] == '#');
1074
1075 sscanf( line, "%f %f", &pb_NX, &pb_NY );
1076 pb_LengthX = pb_NX * pb_Scale;
1077 pb_LengthY = pb_NY * pb_Scale;
1078
1079 // generate artificial file(s)
1080
1081 // a. cerfile
1082
1083 sprintf( cername, "%s/cer999999", get_path_name(0) );
1084
1085 pbfile.open( cername, ios::out );
1086
1087 evth.fill ( 1., 1., 1000., 1., 1.e6, 0., 0., 0., 0., 0., 0., 0.);
1088
1089 evth.write( pbfile );
1090
1091 for ( pbx = (pb_X - pb_LengthX/2.0);
1092 pbx <= (pb_X + pb_LengthX/2.0);
1093 pbx += (pb_LengthX/pb_NX) )
1094
1095 for ( pby = (pb_Y - pb_LengthY/2.0);
1096 pby <= (pb_Y + pb_LengthY/2.0);
1097 pby += (pb_LengthY/pb_NY) ) {
1098
1099 pmfile >> i;
1100
1101 if ( i < 1 )
1102 continue;
1103
1104 pbt = -atan2(sqrt(SQR(pbx)+SQR(pby)),pb_Height);
1105 pbp = atan2(pby, pbx);
1106 photon.fill(450.0,
1107 pbx, pby,
1108 -sin(pbt) * cos(pbp),
1109 -sin(pbt) * sin(pbp),
1110 0., pb_Height);
1111
1112 photon.write( pbfile );
1113
1114 }
1115
1116 photon.fill( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 );
1117
1118 photon.write( pbfile );
1119
1120 pbfile.close();
1121
1122 // b. statfile
1123
1124 sprintf( staname, "%s/sta999999", get_path_name(0) );
1125
1126 pbfile.open( staname );
1127
1128 stat.write( pbfile );
1129
1130 pbfile.close();
1131
1132 pmfile.close();
1133
1134 }
1135
1136 } // YES parallel beam
1137
1138 // read parameters from the stdin
1139
1140 read_ct_file();
1141
1142 // set all random numbers seeds
1143
1144 setall( get_seeds(0), get_seeds(1) );
1145
1146 // get name of the output file
1147
1148 strcpy( outname, get_output_filename() );
1149
1150 // open((re)create) output file
1151
1152 if ( verbose >= VERBOSE_MINIMAL )
1153 log( SIGNATURE, "Openning/creating output file %s\n", outname );
1154
1155 outputfile.open( outname );
1156
1157 if ( outputfile.bad() )
1158 error( SIGNATURE, "Cannot open output file: %s\n", outname );
1159
1160 // save signature to the output file
1161
1162 if ( Data_To_STDOUT )
1163 cout.write( SIGNATURE, sizeof(SIGNATURE) );
1164 else
1165 outputfile.write( SIGNATURE, sizeof(SIGNATURE) );
1166
1167 // generate a sort of log information
1168
1169 if ( verbose >= VERBOSE_MINIMAL ) {
1170 log( SIGNATURE, "Atmospheric model: %s\n", get_atm_mod() );
1171 log( SIGNATURE, "Number of paths: %d\n", get_num_of_paths() );
1172 log( SIGNATURE, "Output filename: %s\n", outname );
1173 log( SIGNATURE, "CT def. filename: %s\n", get_ct_filename() );
1174 }
1175
1176 // get energy range
1177
1178 get_energy_cuts( &lE, &uE );
1179
1180 /*!@'
1181
1182 After reading the parameters file (from STDIN), and the CT
1183 definition file, we start the analysis. Basically, there are three
1184 loops. The outermost loop is a loop in the E@data directories@ we
1185 introduced in the parameters file; the middle loop follows each
1186 |cer*| file in the directory; and the innermost loop gets each
1187 Cherenkov from these files.
1188
1189 */
1190
1191 // for each path of data files
1192
1193 for (np=0; np<get_num_of_paths(); np++) {
1194
1195 strcpy(pathname, get_path_name(np));
1196
1197 // open directory
1198
1199 if ( verbose >= VERBOSE_MINIMAL )
1200 log( SIGNATURE, "Openning directory %s\n", pathname );
1201
1202 directory = opendir(pathname);
1203
1204 if ( directory == NULL )
1205 error( SIGNATURE,
1206 "Cannot open directory %s\n", pathname );
1207
1208 // write start-of-run mark
1209
1210 if ( Data_To_STDOUT )
1211 cout.write( FLAG_START_OF_RUN, SIZE_OF_FLAGS );
1212 else
1213 outputfile.write( FLAG_START_OF_RUN, SIZE_OF_FLAGS );
1214
1215 // trace the number of events (cer* files)
1216
1217 num_cer_files = 0;
1218
1219 // get maximum number, first and last event
1220
1221 max_num_cer_files = get_max_events();
1222
1223 get_range_events(&f_ev, &l_ev);
1224
1225 // if reading from disk, save position of handler in directory
1226
1227 posdir = telldir( directory );
1228
1229 if ( Data_From_STDIN ) {
1230
1231 if ( verbose >= VERBOSE_MINIMAL )
1232 log(SIGNATURE, " * * * READING DATA FROM STDIN * * *\n");
1233
1234 get_stdin_files( 0, lE, uE, TRUE);
1235
1236 }
1237
1238 //--- LOOP ON FILES ---
1239 // for each directory entry (files),
1240 // and(or) while still data has to be read
1241 // this condition is:
1242 // while ( ( there are files AND
1243 // maximum number is not reached )
1244 // OR
1245 // ( reading from stdin AND
1246 // maximum number is not reached ) )
1247
1248 while ( (
1249 ((de = readdir( directory )) != NULL) &&
1250 (num_cer_files < max_num_cer_files)
1251 )
1252 ||
1253 (
1254 ( Data_From_STDIN ) &&
1255 (num_cer_files < max_num_cer_files)
1256 ) ) {
1257
1258 // if Block > 0, then we wait till a file __DOIT is present
1259
1260 if ( Block > 0 ) {
1261
1262 if ( ! WORKING ) {
1263
1264 while ( ! WORKING ) {
1265 // sleep 10 minutes
1266 sleep( 60 );
1267
1268 // look for the file __DOIT
1269 obfile.open("__DOIT");
1270 if ( ! obfile.bad() ) {
1271 WORKING = TRUE;
1272 obfile.close();
1273 }
1274 }
1275
1276 obfile.open("__WORKING");
1277 remove("__DOIT");
1278
1279 rewinddir( directory );
1280 seekdir( directory, posdir );
1281
1282 nfilesBlock = 0;
1283
1284 continue;
1285
1286 } else {
1287
1288 if ( nfilesBlock == Block ) {
1289
1290 posdir = telldir( directory );
1291 obfile.close();
1292 obfile.open("__READY");
1293 obfile << nfilesBlock << " files/block, "
1294 << ncer << " showers.\n";
1295 obfile.close();
1296 remove("__WORKING");
1297 WORKING = FALSE;
1298
1299 continue;
1300 }
1301
1302 }
1303
1304 } // if Block
1305
1306 if ( Data_From_STDIN ) {
1307
1308 if ( verbose >= VERBOSE_MINIMAL )
1309 log(SIGNATURE, "Try to get data from STDIN\n");
1310
1311 if ( ! get_stdin_files(num_cer_files) ) {
1312 if ( verbose >= VERBOSE_MINIMAL )
1313 log(SIGNATURE, "get_stdin_files returned FALSE.\n");
1314 if (cin.eof()) {
1315 num_cer_files = max_num_cer_files + 1;
1316 continue;
1317 }
1318 }
1319
1320 ++num_cer_files;
1321 ncer = num_cer_files;
1322
1323 } else {
1324
1325 // skip removed files
1326
1327 if ( de->d_ino == 0)
1328 continue;
1329
1330 // keep only cer* files
1331
1332 if ( strstr(de->d_name, "cer") != de->d_name )
1333 continue;
1334
1335 if ( Block > 0 ) {
1336 ++nfilesBlock;
1337 obfile << nfilesBlock;
1338 }
1339
1340 // get cer* file number (number of the shower)
1341
1342 // increments the pointer in 3 to calculate the number
1343
1344 ncer = atoi(de->d_name + 3);
1345
1346 if ( (ncer < f_ev) || (ncer > l_ev) )
1347 continue;
1348
1349 // it is a real cer* file, and we want to read it
1350
1351 ++num_cer_files;
1352
1353 if ( verbose >= VERBOSE_NORMAL )
1354 log( SIGNATURE, "cerfile: %s\n", de->d_name );
1355
1356 } // data from disk
1357
1358 // Here we initialize the counters of Cherenkov photons at
1359 // different levels of the analysis.
1360
1361 nCphotons = nbeforeTR = nafterTR = 0;
1362
1363 // full names of the input files cer* and sta*
1364
1365 sprintf(cername, "%s/cer%06d", pathname, ncer);
1366 sprintf(staname, "%s/sta%06d", pathname, ncer);
1367
1368 // try to open cer* file
1369
1370 if ( verbose >= VERBOSE_MAXIMAL )
1371 log( SIGNATURE, "Opening %s\n", cername );
1372
1373 cerfile.open( cername );
1374
1375 if ( cerfile.bad() )
1376 error( SIGNATURE,
1377 "Cannot open input file: %s\n", cername );
1378
1379 /*!@'
1380
1381 Each shower has associated three files: E@ Cherenkov-photons
1382 file @ (|cer*|), E@ Particles file @ (|dat*|) and E@
1383 Statistics file @ (|sta*|). First we read the data from this
1384 last file, which is explained below, and then we read the E@
1385 Event-Header @ block from the |cer*| file.
1386
1387 */
1388
1389 // try to open the statistics file
1390
1391 if ( verbose >= VERBOSE_MAXIMAL )
1392 log( SIGNATURE, "Opening %s\n", staname );
1393
1394 stat.openfile ( staname );
1395
1396 // reading data from this sta* file
1397
1398 if ( verbose >= VERBOSE_MAXIMAL )
1399 log( SIGNATURE, "Reading %s\n", staname );
1400
1401 stat.read();
1402
1403 // read the Event Header from the cer* file
1404
1405 evth.read( cerfile );
1406
1407 if ( verbose >= VERBOSE_MAXIMAL ) {
1408 log( SIGNATURE, "EVENT HEADER:\n");
1409 evth.print();
1410 }
1411
1412 // check energy
1413
1414 if ( (evth.get_energy() < lE) ||
1415 (evth.get_energy() > uE) ) {
1416
1417 // close files
1418 cerfile.close();
1419 stat.closefile();
1420
1421 // skip it
1422 if ( verbose >= VERBOSE_MAXIMAL )
1423 log( SIGNATURE, "Bad energy.\n");
1424
1425 continue;
1426 }
1427
1428 // write start-of-event mark
1429
1430 if ( Data_To_STDOUT )
1431 cout.write( FLAG_START_OF_EVENT, SIZE_OF_FLAGS );
1432 else
1433 outputfile.write( FLAG_START_OF_EVENT, SIZE_OF_FLAGS );
1434
1435 // save information from the cer* to the mc-evth
1436
1437 mcevth.transport( &evth );
1438
1439 // save time interval in the mc-evth
1440
1441 mcevth.put_times ( stat.get_tfirst(), stat.get_tlast() );
1442
1443 // get direction where the CT is pointing to
1444 // (or, better, from when the shower is coming from)
1445
1446 // first, put the values from the shower
1447 thetaCT = evth.get_theta();
1448 phiCT = evth.get_phi();
1449
1450 // then, may be we have to change them
1451
1452 // do we want random pointing (around shower axis) ?
1453 if ( get_random_pointing( &Random_Pointing_MaxDist ) == TRUE ) {
1454
1455 // we do, then get a random position
1456 xiCT = get_new_ct_pointing(thetaCT, phiCT,
1457 Random_Pointing_MaxDist,
1458 &thetaCT, &phiCT);
1459
1460 xiCT = xiCT; // thi is to avoid "was set but not used" messages
1461
1462 // save deviations in mc-evth
1463 mcevth.put_deviations (thetaCT - evth.get_theta(),
1464 phiCT - evth.get_phi());
1465
1466 if ( verbose >= VERBOSE_MAXIMAL )
1467 log( SIGNATURE, "NEW POINTING: (%f,%f) -> (%f,%f)\n",
1468 evth.get_theta(), evth.get_phi(), thetaCT, phiCT);
1469
1470 } else {
1471
1472 mcevth.put_deviations (0.0, 0.0);
1473
1474 if ( get_fixed_target( &thetaCT, &phiCT ) == FALSE ) {
1475 thetaCT = evth.get_theta();
1476 phiCT = evth.get_phi();
1477 }
1478
1479 }
1480
1481 // save event-header to the output file
1482
1483 if ( Data_To_STDOUT )
1484 cout.write ( (char *)&mcevth, mcevth.mysize() );
1485 else
1486 mcevth.write( outputfile );
1487
1488 if ( verbose >= VERBOSE_MAXIMAL )
1489 mcevth.print();
1490
1491 // get core distance and position
1492
1493 coreD = evth.get_core(&coreX, &coreY);
1494 coreD=coreD; // to avoid "was set but not used" messages
1495
1496 // calculate matrices for changing of coordinate system from-to CT
1497
1498 makeOmega( thetaCT, phiCT );
1499 makeOmegaI( thetaCT, phiCT );
1500
1501 memcpy( OmegaCT, Omega, 9*sizeof(float) );
1502 memcpy( OmegaICT, OmegaI, 9*sizeof(float) );
1503
1504 // read each and every single photon (particle) from cer* file
1505
1506 while ( ! (cerfile.eof() || cerfile.bad()) ) {
1507
1508 // read photon data from the file
1509
1510 photon.read( cerfile );
1511
1512 wl = photon.get_wl();
1513
1514 if ( wl < 1.0 )
1515 break;
1516
1517 // get director cosines in z-axis
1518
1519 w = r[2] = photon.get_w();
1520
1521 /*!@'
1522
1523 @#### Atmospheric attenuation.
1524
1525 The routine |atm()| calls the proper attenuation routine,
1526 and returns the transmission coefficient for that
1527 whavelength, height and zenith angle.
1528
1529 */
1530
1531 //++
1532 // FILTER No. 1: ATMOSPHERE
1533 //--
1534
1535 ++nbeforeTR;
1536
1537 // get transmitance of the atmosphere
1538
1539 TR = atm( photon.get_wl(), photon.get_h(), acos(w) );
1540
1541 // passes the atmosphere?
1542
1543 // if random > TR then goes to the TOP of the loop again
1544
1545 if ( RandomNumber > TR )
1546 continue;
1547
1548 ++nafterTR;
1549
1550 /*!@'
1551
1552 @#### Reflectivity of the mirrors.
1553
1554 We make a 3rd. order interpolation using Lagrange
1555 polynomials, in order to calculate the reflectivity of the
1556 mirror for that wavelength. Further developments will
1557 include also a incidence-angle dependence (which is not very
1558 important).
1559
1560 */
1561
1562 //++
1563 // FILTER No. 2: REFLECTIVITY R(lambda)
1564 //--
1565
1566 // find data point to be used in Lagrange interpolation (-> k)
1567
1568 FindLagrange(Reflectivity,k,wl);
1569
1570 // if random > reflectivity then goes to the TOP of the loop again
1571
1572 reflec = Lagrange(Reflectivity,k,wl);
1573
1574 if ( RandomNumber > reflec )
1575 continue;
1576
1577 /*!@'
1578
1579 @#### Reflexion in the mirrors.
1580
1581 The reflexion is made now is a somehow dirty way: if the
1582 value of the variable |ct\_Type| is 0, the telescope is CT1,
1583 and if the value is 1, the telescope is MAGIC. Due to the
1584 completly different geometries of the dishes and the
1585 mirrors, now we support two fully separated and different
1586 algorithms.
1587
1588 */
1589
1590 //++
1591 // REFLEXION
1592 //--
1593
1594 // get director cosines x,y and save original trayectory
1595
1596 u = r[0] = photon.get_u();
1597 v = r[1] = photon.get_v();
1598 u=u, v=v; // to avoid "was set but not used" messages
1599
1600 // get position in the ground
1601
1602 x[0] = photon.get_x() - coreX;
1603 x[1] = photon.get_y() - coreY;
1604 x[2] = 0.0;
1605
1606 if ( verbose >= VERBOSE_MAXIMAL ) {
1607 cout << '@' << '1'
1608 << ' ' << ncer
1609 << ' ' << t
1610 << ' ' << x[0]
1611 << ' ' << x[1]
1612 << ' ' << x[2]
1613 << ' ' << r[0]
1614 << ' ' << r[1]
1615 << ' ' << r[2]
1616 << ' ' << thetaCT
1617 << ' ' << phiCT
1618 << endl;
1619 }
1620
1621 // change to the system of the CT
1622
1623 applyMxV( OmegaCT, x, xCT );
1624 applyMxV( OmegaCT, r, rCT );
1625
1626 // before moving to the system of the mirror, for MAGIC,
1627 // first we look whether the photon hits a mirror or not
1628
1629 // calculate the intersection of the trayectory of the photon
1630 // with the GLOBAL DISH !!!
1631 // we reproduce the calculation of the coefficients of the
1632 // second order polynomial in z (=xCT[2]), made with
1633 // Mathematica
1634
1635 /*
1636 * In[1]:= parab:=z-(x^2+y^2)/(4F)
1637 * par1=parab /. {x->x0+u/w(z-z0),y->y0+v/w(z-z0)}
1638 *
1639 * Out[1]=
1640 * u (z - z0) 2 v (z - z0) 2
1641 * (x0 + ----------) + (y0 + ----------)
1642 * w w
1643 * z - ---------------------------------------
1644 * 4 F
1645 *
1646 * In[2]:= CoefficientList[ExpandAll[par1*4F*w^2],z]
1647 *
1648 * Out[2]=
1649 * 2 2 2 2
1650 * {-(w x0 ) - w y0 + 2 u w x0 z0 +
1651 *
1652 * 2 2 2 2
1653 * 2 v w y0 z0 - u z0 - v z0 ,
1654 *
1655 * 2 2
1656 * 4 F w - 2 u w x0 - 2 v w y0 + 2 u z0 +
1657 *
1658 * 2 2 2
1659 * 2 v z0, -u - v }
1660 */
1661
1662 // the z coordinate is calculated
1663
1664 a = - SQR(rCT[0]) - SQR(rCT[1]);
1665 b = 4.0*ct_Focal_mean*SQR(rCT[2])
1666 - 2.0*rCT[0]*rCT[2]*xCT[0] - 2.0*rCT[1]*rCT[2]*xCT[1]
1667 + 2.0*SQR(rCT[0])*xCT[2] + 2.0*SQR(rCT[1])*xCT[2];
1668 c = 2*rCT[0]*rCT[2]*x[0]*x[2] + 2*rCT[1]*rCT[2]*x[1]*x[2]
1669 - SQR(rCT[2])*SQR(x[0]) - SQR(rCT[2])*SQR(x[1])
1670 - SQR(rCT[0])*SQR(x[2]) - SQR(rCT[1])*SQR(x[2]);
1671
1672 if ( fabs(a) < 1.e-6 ) {
1673
1674 // only one value
1675
1676 xcut[2] = -c / b;
1677
1678 } else {
1679
1680 d = sqrt( b*b - 4.0*a*c );
1681
1682 // two possible values for z
1683
1684 t1 = (-b+d) / (2.0*a);
1685 t2 = (-b-d) / (2.0*a);
1686
1687 // z must be the minimum of t1 and t2
1688
1689 xcut[2] = (t1 < t2) ? t1 : t2;
1690
1691 }
1692
1693 // xcut[] is NOW the cut between the GLOBAL dish of MAGIC and
1694 // the trayectory of the photon
1695
1696 xcut[0] = xCT[0] + rCT[0]/rCT[2]*(xcut[2]-xCT[2]);
1697 xcut[1] = xCT[1] + rCT[1]/rCT[2]*(xcut[2]-xCT[2]);
1698
1699 if ( ct_Type == 1 ) {
1700
1701 // ++ >>>>>>>>>> MAGIC <<<<<<<<<<
1702
1703 // convert to Curvilinear distance over the parabolic dish
1704
1705 sx = Lin2Curv( xcut[0] );
1706 sy = Lin2Curv( xcut[1] );
1707
1708 // is it outside the dish?
1709
1710 if ((fabs(sx) > 850.0) || (fabs(sy) > 850.0)) {
1711 //cout << "CONDICION 1 !" << endl << flush;
1712 //cout << '1';
1713 continue;
1714 }
1715
1716 // -- >>>>>>>>>> MAGIC <<<<<<<<<<
1717
1718 } // endif (ct_Type == 1)
1719
1720 // calculate the mirror to be used
1721
1722 distmirr = 1000000.;
1723
1724 for (i=0; i<ct_NMirrors && distmirr>=ct_RMirror; ++i) {
1725 distmirr2 = sqrt(SQR(ct_data[i][CT_X] - xcut[0]) +
1726 SQR(ct_data[i][CT_Y] - xcut[1]) +
1727 SQR(ct_data[i][CT_Z] - xcut[2]));
1728 if (distmirr2 < distmirr) {
1729 i_mirror = i;
1730 distmirr = distmirr2;
1731 }
1732 }
1733
1734 // the mirror to use is i_mirror (calculated several lines above)
1735 // check whether the photon is outside the nearest (this) mirror
1736
1737 if ( ct_Type == 0 ) {
1738
1739 // ++ >>>>>>>>>> CT1 <<<<<<<<<<
1740 if (distmirr > ct_RMirror) {
1741 //cout << "CONDICION 2 !" << endl << flush;
1742 //cout << '2';
1743 continue;
1744 }
1745 // -- >>>>>>>>>> CT1 <<<<<<<<<<
1746
1747 } else {
1748
1749 // ++ >>>>>>>>>> MAGIC <<<<<<<<<<
1750 if ((fabs(ct_data[i_mirror][CT_SX] - sx) > ct_RMirror) ||
1751 (fabs(ct_data[i_mirror][CT_SY] - sy) > ct_RMirror)) {
1752 //cout << "CONDICION 2 !" << endl << flush;
1753 //cout << '2';
1754 continue;
1755 }
1756 // -- >>>>>>>>>> MAGIC <<<<<<<<<<
1757
1758 }
1759
1760 // calculate matrices for the mirror
1761
1762 makeOmega (-ct_data[i_mirror][CT_THETAN],
1763 ct_data[i_mirror][CT_PHIN]);
1764 makeOmegaI(-ct_data[i_mirror][CT_THETAN],
1765 ct_data[i_mirror][CT_PHIN]);
1766
1767 // change to the system of the mirror
1768
1769 // first translation...
1770
1771 xmm[0] = xCT[0] - ct_data[i_mirror][CT_X];
1772 xmm[1] = xCT[1] - ct_data[i_mirror][CT_Y];
1773 xmm[2] = xCT[2] - ct_data[i_mirror][CT_Z];
1774
1775 // ...then rotation
1776
1777 applyMxV( Omega, xmm, xm );
1778 applyMxV( Omega, rCT, rm );
1779
1780 // the vector rCT should be normalized, and
1781 // so the vector rm remains normalized as well, but, anyhow...
1782
1783 rnorm = NORM( rm );
1784 rm[0] /= rnorm;
1785 rm[1] /= rnorm;
1786 rm[2] /= rnorm;
1787
1788 // calculate the intersection of the trayectory of the photon
1789 // with the mirror
1790 // we reproduce the calculation of the coefficients of the
1791 // second order polynomial in z (=xm[2]), made with
1792 // Mathematica
1793
1794 /*
1795 * In[1]:= esfera:=x^2+y^2+(z-R)^2-R^2;
1796 * recta:={x->x0+u/w(z-z0),y->y0+v/w(z-z0)}
1797 *
1798 * In[2]:= esfera
1799 *
1800 * 2 2 2 2
1801 * Out[2]= -R + x + y + (-R + z)
1802 *
1803 * In[3]:= recta
1804 *
1805 * u (z - z0) v (z - z0)
1806 * Out[3]= {x -> x0 + ----------, y -> y0 + ----------}
1807 * w w
1808 *
1809 * In[4]:= esf=esfera /. recta
1810 *
1811 * 2 2 u (z - z0) 2 v (z - z0) 2
1812 * Out[4]= -R + (-R + z) + (x0 + ----------) + (y0 + ----------)
1813 * w w
1814 *
1815 * In[5]:= coefs=CoefficientList[ExpandAll[esf],z]
1816 *
1817 * 2 2 2 2
1818 * 2 2 2 u x0 z0 2 v y0 z0 u z0 v z0
1819 * Out[5]= {x0 + y0 - --------- - --------- + ------ + ------,
1820 * w w 2 2
1821 * w w
1822 *
1823 * 2 2 2 2
1824 * 2 u x0 2 v y0 2 u z0 2 v z0 u v
1825 * > -2 R + ------ + ------ - ------- - -------, 1 + -- + --}
1826 * w w 2 2 2 2
1827 * w w w w
1828 * In[6]:= Simplify[ExpandAll[coefs*w^2]]
1829 *
1830 * 2 2 2 2 2 2
1831 * Out[6]= {w (x0 + y0 ) - 2 w (u x0 + v y0) z0 + (u + v ) z0 ,
1832 *
1833 * 2 2 2 2 2
1834 * > -2 (R w - u w x0 + u z0 + v (-(w y0) + v z0)), u + v + w }
1835 *
1836 */
1837
1838 // the z coordinate is calculated, using the coefficients
1839 // shown above
1840
1841 a = SQR(rm[0]) + SQR(rm[1]) + SQR(rm[2]);
1842 if ( ct_Type == 0 ) {
1843 // CT1
1844 b = -2*(2.*ct_Focal_mean*SQR(rm[2])
1845 - rm[0]*rm[2]*xm[0]
1846 + SQR(rm[0])*xm[2]
1847 + rm[1]*(-(rm[2]*xm[1]) + rm[1]*xm[2]));
1848 } else {
1849 // MAGIC
1850 b = -2*(2.*ct_data[i_mirror][CT_FOCAL]*SQR(rm[2])
1851 - rm[0]*rm[2]*xm[0]
1852 + SQR(rm[0])*xm[2]
1853 + rm[1]*(-(rm[2]*xm[1]) + rm[1]*xm[2]));
1854 }
1855 c = (SQR(rm[2])*(SQR(xm[0]) + SQR(xm[1]))
1856 - 2*rm[2]*(rm[0]*xm[0] + rm[1]*xm[1])*xm[2] +
1857 (SQR(rm[0]) + SQR(rm[1]))*SQR(xm[2]));
1858
1859 d = sqrt( b*b - 4.0*a*c );
1860
1861 // two possible values for z
1862
1863 t1 = (-b+d) / (2.0*a);
1864 t2 = (-b-d) / (2.0*a);
1865
1866 // z must be the minimum of t1 and t2
1867
1868 xcut[2] = (t1 < t2) ? t1 : t2;
1869 xcut[0] = xm[0] + rm[0]/rm[2]*(xcut[2]-xm[2]);
1870 xcut[1] = xm[1] + rm[1]/rm[2]*(xcut[2]-xm[2]);
1871
1872 //++
1873 // BLACK SPOTS: If the photon hits the black spot, it's lost
1874 //--
1875
1876 if ( sqrt(SQR(xcut[0]) + SQR(xcut[1])) < ct_BlackSpot_rad ) {
1877 //cout << "CONDITION 3!\n" << flush;
1878 //cout << '3';
1879 continue;
1880 }
1881
1882 // if still we have the photon, continue with the reflexion
1883
1884 // calculate normal vector in this point
1885 // (and normalize, with the sign changed)
1886
1887 rnor[0] = 2.0*xcut[0];
1888 rnor[1] = 2.0*xcut[1];
1889 rnor[2] = 2.0*(xcut[2] - 2.0*ct_Focal[i_mirror]);
1890
1891 rnorm = -NORM( rnor );
1892 rnor[0] /= rnorm;
1893 rnor[1] /= rnorm;
1894 rnor[2] /= rnorm;
1895
1896 // now, both "normal" vector and original trayectory are
1897 // normalized
1898 // just project the original vector in the normal, and
1899 // take it as the "mean" position of the original and
1900 // the "reflected" vector
1901 // from this, we can calculate the "reflected" vector
1902 // calpha = cos(angle(rnor,rm))
1903
1904 calpha = fabs(rnor[0]*rm[0] + rnor[1]*rm[1] + rnor[2]*rm[2]);
1905
1906 // finally!!! we have the reflected trayectory of the photon
1907
1908 rrefl[0] = 2.0*rnor[0]*calpha - rm[0];
1909 rrefl[1] = 2.0*rnor[1]*calpha - rm[1];
1910 rrefl[2] = 2.0*rnor[2]*calpha - rm[2];
1911
1912 rnorm = NORM( rrefl );
1913 rrefl[0] /= rnorm;
1914 rrefl[1] /= rnorm;
1915 rrefl[2] /= rnorm;
1916
1917 // let's go back to the coordinate system of the CT
1918
1919 // first rotation...
1920 applyMxV( OmegaI, xcut, xcutCT);
1921 applyMxV( OmegaI, rrefl, rreflCT);
1922
1923 // ...then translation
1924 xcutCT[0] += ct_data[i_mirror][CT_X];
1925 xcutCT[1] += ct_data[i_mirror][CT_Y];
1926 xcutCT[2] += ct_data[i_mirror][CT_Z];
1927
1928 // calculate intersection of this trayectory and the camera plane
1929 // in the system of the CT, this plane is z = ct_Focal
1930
1931 t = (ct_Focal_mean - xcutCT[2]) / rreflCT[2];
1932
1933 xcam[0] = xcutCT[0] + rreflCT[0]*t;
1934 xcam[1] = xcutCT[1] + rreflCT[1]*t;
1935 xcam[2] = xcutCT[2] + rreflCT[2]*t;
1936
1937 //++
1938 // AXIS DEVIATION: We introduce it here just as a first order
1939 // correction, by modifying the position of the reflected photon.
1940 //--
1941
1942 xcam[0] += AxisDeviation[0][i_mirror];
1943 xcam[1] += AxisDeviation[1][i_mirror];
1944
1945 //++
1946 // SMEARING: We apply the point spread function for the mirrors
1947 //--
1948
1949 // get two N(0;1) random numbers
1950
1951 rnormal( NormalRandomNumbers, 2 );
1952
1953 // modify the Cphoton position in the camera
1954
1955 xcam[0] += NormalRandomNumbers[0] * ct_PSpread_mean;
1956 xcam[1] += NormalRandomNumbers[1] * ct_PSpread_mean;
1957
1958 // check whether the photon goes out of the camera
1959
1960 if ( (SQR(xcam[0])+SQR(xcam[1])) > SQR(ct_CameraWidth) ) {
1961 continue;
1962 }
1963
1964 //++
1965 // ANGLE OF INCIDENCE
1966 //--
1967
1968 // calculate angle of incidence between tray. and camera plane
1969 // the camera plane is
1970 // 0 y + 0 y + z - ct_Focal = 0 => (A,B,C,D) = (0,0,1,-ct_Focal)
1971 // from Table 3.20 "Tasch. der Math."
1972
1973 phi = asin(rreflCT[2]);
1974
1975 //++
1976 // TIMING
1977 //--
1978
1979 // calculate the new time of the photon (in the camera)
1980
1981 // initial time since first interaction
1982
1983 t = photon.get_t();
1984
1985 // substract path from the mirror till the ground, 'cos
1986 // the photon actually hit the mirror!!
1987
1988 t = t + ((( xm[2] > 0. ) ? -1.0 : +1.0) *
1989 sqrt( SQR(xm[0] - xcut[0]) +
1990 SQR(xm[1] - xcut[1]) +
1991 SQR(xm[2] - xcut[2]) ) / Speed_of_Light_air_cmns);
1992 // add path from the mirror till the camera
1993
1994 t = t + sqrt( SQR(xcutCT[0] - xcam[0]) +
1995 SQR(xcutCT[1] - xcam[1]) +
1996 SQR(xcutCT[2] - xcam[2]) ) / Speed_of_Light_air_cmns;
1997
1998 // save data in the photon variable
1999
2000 cphoton.fill( photon.get_id(),
2001 xcam[0], xcam[1],
2002 photon.get_u(), photon.get_v(),
2003 photon.get_h(), t, phi);
2004
2005 // show it
2006
2007 if ( verbose >= VERBOSE_MAXIMAL ) {
2008 cout << '@' << '2'
2009 << ' ' << ncer
2010 << ' ' << xCT[0]
2011 << ' ' << xCT[1]
2012 << ' ' << xCT[2]
2013 << ' ' << rCT[0]
2014 << ' ' << rCT[1]
2015 << ' ' << rCT[2]
2016 << ' ' << xcut[0]
2017 << ' ' << xcut[1]
2018 << ' ' << xcut[2]
2019 << endl;
2020 cout << '@' << '3'
2021 << ' ' << ncer
2022 << ' ' << sx
2023 << ' ' << sy
2024 << ' ' << i_mirror
2025 << ' ' << ct_data[i_mirror][CT_SX]
2026 << ' ' << ct_data[i_mirror][CT_SY]
2027 << ' ' << ct_data[i_mirror][CT_SX] - sx
2028 << ' ' << ct_data[i_mirror][CT_SY] - sy
2029 << endl;
2030 if ( pb_ParallelBeam == FALSE ) {
2031 cout << '@' << '4'
2032 << ' ' << ncer
2033 << ' ' << xcam[0]
2034 << ' ' << xcam[1]
2035 << ' ' << xcam[2]
2036 << ' ' << phi
2037 << ' ' << photon.get_t()-stat.get_tfirst()
2038 << ' ' << t-stat.get_tfirst()
2039 << endl << flush;
2040 } else {
2041 cout << '@' << '4'
2042 << ' ' << ncer
2043 << ' ' << xcam[0]
2044 << ' ' << xcam[1]
2045 << ' ' << xcam[2]
2046 << ' ' << phi
2047 << endl << flush;
2048 }
2049 }
2050
2051 // write the photon data to the output file
2052
2053 if ( Data_To_STDOUT )
2054 cout.write ( (char *)&cphoton, cphoton.mysize() );
2055 else
2056 cphoton.write( outputfile );
2057
2058 // increase number of Cphotons written
2059
2060 ++nCphotons;
2061
2062 } // while still there are photons left
2063
2064 if ( verbose >= VERBOSE_NORMAL )
2065 cout << nCphotons << ' '
2066 << nbeforeTR << ' '
2067 << nafterTR << ' '
2068 << endl << flush;
2069
2070 // write end-of-event mark
2071
2072 if ( Data_To_STDOUT )
2073 cout.write( FLAG_END_OF_EVENT, SIZE_OF_FLAGS );
2074 else
2075 outputfile.write( FLAG_END_OF_EVENT, SIZE_OF_FLAGS );
2076
2077 // close files
2078
2079 cerfile.close();
2080 stat.closefile();
2081
2082 // show how many photons were written
2083
2084 if ( verbose >= VERBOSE_MAXIMAL )
2085 log( SIGNATURE, "%d C-photons written.\n", nCphotons );
2086
2087 // if we are reading data from STDIN, delete the last used files
2088
2089 if ( Data_From_STDIN ) {
2090
2091 unlink( cername );
2092 unlink( staname );
2093
2094 }
2095
2096 } // while there are still data left
2097
2098 // write end-of-run mark
2099
2100 if ( Data_To_STDOUT )
2101 cout.write( FLAG_END_OF_RUN, SIZE_OF_FLAGS );
2102 else
2103 outputfile.write( FLAG_END_OF_RUN, SIZE_OF_FLAGS );
2104
2105 // close directory
2106
2107 closedir( directory );
2108
2109 } // for each data directory
2110
2111 // write end-of-file mark
2112
2113 if ( Data_To_STDOUT )
2114 cout.write( FLAG_END_OF_FILE, SIZE_OF_FLAGS );
2115 else
2116 outputfile.write( FLAG_END_OF_FILE, SIZE_OF_FLAGS );
2117
2118 // close output file
2119
2120 if ( verbose >= VERBOSE_MINIMAL )
2121 log( SIGNATURE, "Closing output file %s\n", outname );
2122
2123 if ( ! Data_To_STDOUT )
2124 outputfile.close();
2125
2126 // clean everything
2127
2128 if ( verbose >= VERBOSE_MINIMAL )
2129 log( SIGNATURE, "Cleanning . . .\n" );
2130
2131 clean();
2132
2133 // program finished
2134
2135 if ( verbose >= VERBOSE_MINIMAL )
2136 log( SIGNATURE, "Done.\n");
2137
2138 return ( 0 );
2139}
2140//!@}
2141
2142// @T \newpage
2143
2144//=---------------------------------------------------------------------
2145// @subsection Functions definition
2146
2147//!---------------------------------------------------------------------
2148// @name present
2149//
2150// @desc Make some presentation
2151//
2152// @date Sat Jun 27 05:58:56 MET DST 1998
2153//----------------------------------------------------------------------
2154// @function
2155
2156//!@{
2157void
2158present(void)
2159{
2160 cout << "##################################################\n"
2161 << SIGNATURE << '\n' << '\n'
2162 << "Simulation of the reflector\n"
2163 << "J C Gonzalez, May 1998\n"
2164 << "##################################################\n\n";
2165}
2166//!@}
2167
2168
2169//!---------------------------------------------------------------------
2170// @name usage
2171//
2172// @desc show help
2173//
2174// @date Tue Dec 15 16:23:30 MET 1998
2175//----------------------------------------------------------------------
2176// @function
2177
2178//!@{
2179void
2180usage(void)
2181{
2182 present();
2183 cout << "\nusage ::\n\n"
2184 << "\t reflector "
2185 << " [ -@ paramfile ] "
2186 << " [ -h ] "
2187 << "\n\n or \n\n"
2188 << "\t reflector < paramfile"
2189 << "\n\n";
2190 exit(0);
2191}
2192//!@}
2193
2194
2195//!---------------------------------------------------------------------
2196// @name clean
2197//
2198// @desc basically, frees memory
2199//
2200// @date Sat Jun 27 05:58:56 MET DST 1998
2201//----------------------------------------------------------------------
2202// @function
2203
2204//!@{
2205void
2206clean(void)
2207{
2208 int i;
2209
2210 // delete focals table
2211
2212 delete [] ct_Focal;
2213 ct_Focal = NULL;
2214
2215 // delete reflectivity table
2216
2217 for (i=0; i<2; i++) {
2218 delete [] Reflectivity[i];
2219 }
2220
2221 delete [] Reflectivity;
2222 Reflectivity = NULL;
2223
2224 // delete mirrors' data table
2225
2226 for (i=0; i<ct_NMirrors; i++) {
2227 delete [] ct_data[i];
2228 }
2229
2230 delete [] ct_data;
2231 ct_data = NULL;
2232}
2233//!@}
2234
2235
2236//!---------------------------------------------------------------------
2237// @name log
2238//
2239// @desc function to send log information
2240//
2241// @var funct Name of the caller function
2242// @var fmt Format to be used (message)
2243// @var ... Other information to be shown
2244//
2245// @date Sat Jun 27 05:58:56 MET DST 1998
2246//----------------------------------------------------------------------
2247// @function
2248
2249//!@{
2250void
2251log(const char *funct, char *fmt, ...)
2252{
2253 va_list args;
2254
2255 // Display the name of the function that called error
2256 printf("[%s]: ", funct);
2257
2258 // Display the remainder of the message
2259 va_start(args, fmt);
2260 vprintf(fmt, args);
2261 va_end(args);
2262}
2263//!@}
2264
2265
2266//!---------------------------------------------------------------------
2267// @name logerr
2268//
2269// @desc function to send log information to stderr
2270//
2271// @var funct Name of the caller function
2272// @var fmt Format to be used (message)
2273// @var ... Other information to be shown
2274//
2275// @date Sat Jun 27 05:58:56 MET DST 1998
2276//----------------------------------------------------------------------
2277// @function
2278
2279//!@{
2280void
2281logerr(const char *funct, char *fmt, ...)
2282{
2283 va_list args;
2284
2285 // Display the name of the function that called error
2286 printf("[%s]: ", funct);
2287
2288 // Display the remainder of the message
2289 va_start(args, fmt);
2290 vfprintf(stderr, fmt, args);
2291 va_end(args);
2292}
2293//!@}
2294
2295
2296//!---------------------------------------------------------------------
2297// @name error
2298//
2299// @desc function to send an error message, and abort the program
2300//
2301// @var funct Name of the caller function
2302// @var fmt Format to be used (message)
2303// @var ... Other information to be shown
2304//
2305// @date Sat Jun 27 05:58:56 MET DST 1998
2306//----------------------------------------------------------------------
2307// @function
2308
2309//!@{
2310void
2311error(const char *funct, char *fmt, ...)
2312{
2313 va_list args;
2314
2315 // Display the name of the function that called error
2316 fprintf(stderr, "ERROR in %s: ", funct);
2317
2318 // Display the remainder of the message
2319 va_start(args, fmt);
2320 vfprintf(stderr, fmt, args);
2321 va_end(args);
2322
2323 abort();
2324}
2325//!@}
2326
2327
2328//!---------------------------------------------------------------------
2329// @name makeOmega
2330//
2331// @desc function to calculate the matrix Omega(theta,phi)
2332//
2333// @var theta Angle theta of the transformation
2334// @var phi Angle phi of the transformation
2335//
2336// @date Sat Jun 27 05:58:56 MET DST 1998
2337//----------------------------------------------------------------------
2338// @function
2339
2340//!@{
2341void
2342makeOmega (float theta, float phi)
2343{
2344 static float ct, st, cp, sp;
2345
2346 // shortcuts for cosine and sine of theta and phi
2347 ct = cos(theta);
2348 st = sin(theta);
2349 cp = cos(phi);
2350 sp = sin(phi);
2351
2352 // save values in the array (see top of file)
2353 Omega[0][0] = cp*ct;
2354 Omega[0][1] = sp*ct;
2355 Omega[0][2] = -st;
2356
2357 Omega[1][0] = -sp;
2358 Omega[1][1] = cp;
2359 Omega[1][2] = 0;
2360
2361 Omega[2][0] = cp*st;
2362 Omega[2][1] = sp*st;
2363 Omega[2][2] = ct;
2364}
2365//!@}
2366
2367
2368//!---------------------------------------------------------------------
2369// @name makeOmegaI
2370//
2371// @desc function to calculate the matrix Omega-1(theta,phi)
2372//
2373// @var theta Angle theta of the transformation
2374// @var phi Angle phi of the transformation
2375//
2376// @date Sat Jun 27 05:58:56 MET DST 1998
2377//----------------------------------------------------------------------
2378// @function
2379
2380//!@{
2381void
2382makeOmegaI(float theta, float phi)
2383{
2384 static float ct, st, cp, sp;
2385
2386 // shortcuts for cosine and sine of theta and phi
2387 ct = cos(theta);
2388 st = sin(theta);
2389 cp = cos(phi);
2390 sp = sin(phi);
2391
2392 // save values in the array (see top of file)
2393 OmegaI[0][0] = cp*ct;
2394 OmegaI[0][1] = -sp;
2395 OmegaI[0][2] = cp*st;
2396
2397 OmegaI[1][0] = sp*ct;
2398 OmegaI[1][1] = cp;
2399 OmegaI[1][2] = sp*st;
2400
2401 OmegaI[2][0] = -st;
2402 OmegaI[2][1] = 0;
2403 OmegaI[2][2] = ct;
2404}
2405//!@}
2406
2407
2408//!---------------------------------------------------------------------
2409// @name applyMxv
2410//
2411// @desc returns the vector v' such that v' = M x v
2412//
2413// @var M matrix of the transformation
2414// @var v vector to be multiplied
2415// @var vi resulting vector
2416//
2417// @date Sat Jun 27 05:58:56 MET DST 1998
2418//----------------------------------------------------------------------
2419// @function
2420
2421//!@{
2422void
2423applyMxV(float M[3][3], float *V, float *Vp)
2424{
2425 Vp[0] = (M[0][0] * V[0] +
2426 M[0][1] * V[1] +
2427 M[0][2] * V[2]);
2428 Vp[1] = (M[1][0] * V[0] +
2429 M[1][1] * V[1] +
2430 M[1][2] * V[2]);
2431 Vp[2] = (M[2][0] * V[0] +
2432 M[2][1] * V[1] +
2433 M[2][2] * V[2]);
2434}
2435//!@}
2436
2437
2438//!---------------------------------------------------------------------
2439// @name read_ct_file
2440//
2441// @desc read CT definition file
2442//
2443// @date Sat Jun 27 05:58:56 MET DST 1998
2444//----------------------------------------------------------------------
2445// @function
2446
2447//!@{
2448void
2449read_ct_file(void)
2450{
2451 char line[LINE_MAX_LENGTH]; // line to get from the ctin
2452 char token[ITEM_MAX_LENGTH]; // a single token
2453 char fmirr[40]; // name of the binary version of mirrors data
2454 int i; // dummy counters
2455
2456 if ( verbose >= VERBOSE_MINIMAL )
2457 log( "read_ct_file", "start.\n" );
2458
2459 ifstream ctin ( get_ct_filename() );
2460
2461 if ( ctin.bad() )
2462 error( "read_ct_file",
2463 "Cannot open CT def. file: %s\n", get_ct_filename() );
2464
2465 // loop till the "end" directive is reached
2466
2467 while (!ctin.eof()) {
2468
2469 // get line from stdin
2470
2471 ctin.getline(line, LINE_MAX_LENGTH);
2472
2473 // look for each item at the beginning of the line
2474
2475 for (i=0; i<=define_mirrors; i++)
2476 if (strstr(line, CT_ITEM_NAMES[i]) == line)
2477 break;
2478
2479 // if it is not a valid line, just ignore it
2480
2481 if (i == define_mirrors+1)
2482 continue;
2483
2484 // case block for each directive
2485
2486 switch ( i ) {
2487
2488 case type: // <type of telescope> (0:CT1 ¦ 1:MAGIC)
2489
2490 // get focal distance
2491
2492 sscanf(line, "%s %d", token, &ct_Type);
2493
2494 if ( verbose >= VERBOSE_MINIMAL )
2495 log( "read_ct_file", "<Type of Telescope>: %s\n",
2496 ((ct_Type==0) ? "CT1" : "MAGIC") );
2497
2498 break;
2499
2500 case focal_distance: // <focal distance> [cm]
2501
2502 // get focal distance
2503
2504 sscanf(line, "%s %f", token, &ct_Focal_mean);
2505
2506 if ( verbose >= VERBOSE_MINIMAL )
2507 log( "read_ct_file", "<Focal distance>: %f cm\n", ct_Focal_mean );
2508
2509 break;
2510
2511 case focal_std: // s(focal distance) [cm]
2512
2513 // get focal distance
2514
2515 sscanf(line, "%s %f", token, &ct_Focal_std);
2516
2517 if ( verbose >= VERBOSE_MINIMAL )
2518 log( "read_ct_file", "s(Focal distance): %f cm\n", ct_Focal_std );
2519
2520 break;
2521
2522 case point_spread: // <point spread> [cm]
2523
2524 // get point spread
2525
2526 sscanf(line, "%s %f", token, &ct_PSpread_mean);
2527
2528 if ( verbose >= VERBOSE_MINIMAL )
2529 log( "read_ct_file", "<Point spread>: %f cm\n", ct_PSpread_mean );
2530
2531 break;
2532
2533 case point_std: // s(point spread) [cm]
2534
2535 // get point spread
2536
2537 sscanf(line, "%s %f", token, &ct_PSpread_std);
2538
2539 if ( verbose >= VERBOSE_MINIMAL )
2540 log( "read_ct_file", "s(Point spread): %f cm\n", ct_PSpread_std );
2541
2542 break;
2543
2544 case adjustment_dev: // s(adjustment_dev) [cm]
2545
2546 // get point spread
2547
2548 sscanf(line, "%s %f", token, &ct_Adjustment_std);
2549
2550 if ( verbose >= VERBOSE_MINIMAL )
2551 log( "read_ct_file", "s(Adjustment): %f cm\n", ct_Adjustment_std );
2552
2553 break;
2554
2555 case black_spot: // radius of the black spot in the center [cm]
2556
2557 // get black spot radius
2558
2559 sscanf(line, "%s %f", token, &ct_BlackSpot_rad);
2560
2561 if ( verbose >= VERBOSE_MINIMAL )
2562 log( "read_ct_file", "Radius of the black spots: %f cm\n",
2563 ct_BlackSpot_rad);
2564
2565 break;
2566
2567 case r_mirror: // radius of the mirrors [cm]
2568
2569 // get radius of mirror
2570
2571 sscanf(line, "%s %f", token, &ct_RMirror);
2572
2573 if ( verbose >= VERBOSE_MINIMAL )
2574 log( "read_ct_file", "Radii of the mirrors: %f cm\n", ct_RMirror );
2575
2576 break;
2577
2578 case n_mirrors: // number of mirrors
2579
2580 // get the name of the output_file from the line
2581
2582 sscanf(line, "%s %d", token, &ct_NMirrors);
2583
2584 if ( verbose >= VERBOSE_MINIMAL )
2585 log( "read_ct_file", "Number of mirrors: %d\n", ct_NMirrors );
2586
2587 break;
2588
2589 case camera_width: // camera width [cm]
2590
2591 // get the name of the ct_file from the line
2592
2593 sscanf(line, "%s %f", token, &ct_CameraWidth);
2594
2595 if ( verbose >= VERBOSE_MINIMAL )
2596 log( "read_ct_file", "Camera width: %f cm\n", ct_CameraWidth );
2597
2598 break;
2599
2600 case n_pixels: // number of pixels
2601
2602 // get the name of the output_file from the line
2603
2604 sscanf(line, "%s %d", token, &ct_NPixels);
2605
2606 if ( verbose >= VERBOSE_MINIMAL )
2607 log( "read_ct_file", "Number of pixels: %d\n", ct_NPixels );
2608
2609 break;
2610
2611 case pixel_width: // pixel width [cm]
2612
2613 // get the name of the ct_file from the line
2614
2615 sscanf(line, "%s %f", token, &ct_PixelWidth);
2616
2617 if ( verbose >= VERBOSE_MINIMAL )
2618 log( "read_ct_file", "Pixel width: %f cm\n", ct_PixelWidth );
2619
2620 break;
2621
2622 case define_mirrors: // read table with the parameters of the mirrors
2623
2624 if ( verbose >= VERBOSE_MINIMAL )
2625 log( "read_ct_file", "Table of mirrors data:\n" );
2626
2627 // check whether the number of mirrors was already set
2628
2629 if ( ct_NMirrors == 0 )
2630 error( "read_ct_file", "NMirrors was not set.\n" );
2631
2632 // allocate memory for paths list
2633
2634 if ( verbose >= VERBOSE_MINIMAL )
2635 log( "read_ct_file", "Allocating memory for ct_data\n" );
2636
2637 ct_data = new float*[ct_NMirrors];
2638
2639 for (i=0; i<ct_NMirrors; i++)
2640 ct_data[i] = new float[CT_NDATA];
2641
2642 // read data
2643
2644 // look for binary version of mirror data (faster to load)
2645
2646 sprintf( fmirr, "%s.mirr", get_ct_filename() );
2647
2648 ifstream ctmirr;
2649 ctmirr.open( fmirr );
2650
2651 if ( ctmirr.bad() ) {
2652
2653 // the file does not exist
2654
2655 ctmirr.close(); // close not existing file
2656
2657 ofstream ctmirr_write ( fmirr ); // open file to save the data
2658
2659 if ( verbose >= VERBOSE_MINIMAL )
2660 log( "read_ct_file", "Reading mirrors data...\n" );
2661
2662 // read ASCII data
2663
2664 for (i=0; i<ct_NMirrors; i++) {
2665 // for (j=0; j<CT_NDATA; j++)
2666
2667 ctin.getline(line, LINE_MAX_LENGTH);
2668 sscanf(line, "%f %f %f %f %f %f %f %f %f %f %f %f",
2669 &ct_data[i][0],
2670 &ct_data[i][1],
2671 &ct_data[i][2],
2672 &ct_data[i][3],
2673 &ct_data[i][4],
2674 &ct_data[i][5],
2675 &ct_data[i][6],
2676 &ct_data[i][7],
2677 &ct_data[i][8],
2678 &ct_data[i][9],
2679 &ct_data[i][10],
2680 &ct_data[i][11]);
2681
2682 cout << '[' << i << ']' << flush;
2683 }
2684 cout << endl << flush;
2685
2686 // save binary data and close output file
2687
2688 for (i=0; i<ct_NMirrors; i++)
2689 ctmirr_write.write( (char*)ct_data[i], CT_NDATA*sizeof(float) );
2690
2691 ctmirr_write.close();
2692
2693 } else {
2694
2695 // Excellent!
2696 // read binary data and close input file
2697
2698 for (i=0; i<ct_NMirrors; i++)
2699 ctmirr.read( (char*)ct_data[i], CT_NDATA*sizeof(float) );
2700
2701 ctmirr.close();
2702
2703 }
2704
2705 break;
2706
2707 } // switch ( i )
2708
2709 } // while (! is_end)
2710
2711 // read reflectivity file
2712
2713 read_reflectivity();
2714
2715 // read axis deviations file
2716
2717 read_axisdev();
2718
2719 // read focals file
2720
2721 read_focals();
2722
2723 // end
2724
2725 if ( verbose >= VERBOSE_MINIMAL )
2726 log( "read_ct_file", "done.\n" );
2727 return;
2728}
2729//!@}
2730
2731
2732//!---------------------------------------------------------------------
2733// @name dist_r_P
2734//
2735// @desc distance straight line r - point P
2736//
2737// @var a coord. X of a fixed point of the straight line
2738// @var b coord. Y of a fixed point of the straight line
2739// @var c coord. Z of a fixed point of the straight line
2740// @var l component X of a vector of the straight line
2741// @var m component Y of a vector of the straight line
2742// @var n component Z of a vector of the straight line
2743// @var x coord. X of the point P
2744// @var y coord. Y of the point P
2745// @var z coord. Z of the point P
2746//
2747// @return Returns the distance d(r,P) between the line and the point P
2748//
2749// @date Sat Jun 27 05:58:56 MET DST 1998
2750//----------------------------------------------------------------------
2751// @function
2752
2753//!@{
2754float
2755dist_r_P(float a, float b, float c,
2756 float l, float m, float n,
2757 float x, float y, float z)
2758{
2759 return (
2760 sqrt((SQR((a-x)*m-(b-y)*l) +
2761 SQR((b-y)*n-(c-z)*m) +
2762 SQR((c-z)*l-(a-x)*n))/
2763 (SQR(l)+SQR(m)+SQR(n))
2764 )
2765 );
2766}
2767//!@}
2768
2769
2770//!---------------------------------------------------------------------
2771// @name read_reflectivity
2772//
2773// @desc read reflectivity data for the mirrors
2774//
2775// @date Sat Jun 27 05:58:56 MET DST 1998
2776//----------------------------------------------------------------------
2777// @function
2778
2779//!@{
2780void
2781read_reflectivity(void)
2782{
2783 ifstream reflfile;
2784 ofstream reflfile_write;
2785 int i, errors=0;
2786 char line[LINE_MAX_LENGTH];
2787
2788 // try to open the file
2789 reflfile.open( REFLECTIVITY_FILE );
2790
2791 // if it is wrong or does not exist:
2792
2793 while ( reflfile.bad() ) {
2794
2795 ++errors;
2796
2797 if ( verbose >= VERBOSE_MINIMAL )
2798 log( "read_reflectivity", "Cannot open reflectivity file: %s\n",
2799 REFLECTIVITY_FILE );
2800
2801 if ( errors > 5 )
2802 error( "read_reflectivity", "Exiting.");
2803
2804 // try to re-generate the file
2805
2806 if ( verbose >= VERBOSE_MINIMAL )
2807 log( "read_reflectivity", "Generating file: %s\n",
2808 REFLECTIVITY_FILE );
2809
2810 // try to open the file
2811 reflfile_write.open( REFLECTIVITY_FILE );
2812
2813 // write short comment in the beginning
2814
2815 reflfile_write << "# reflectivity file" << endl
2816 << "# reflectivity for each mirror in the frame" << endl
2817 << "# J C Gonzalez, 1998" << endl
2818 << "#" << endl;
2819
2820 // write data
2821
2822 // write number of data-points
2823
2824 reflfile_write << "# number of datapoints" << endl
2825 << 35 << endl;
2826
2827 // write data-points
2828
2829 reflfile_write << "# datapoints" << endl;
2830
2831 for ( i=0; i<35; ++i )
2832 reflfile_write << 270.0+i*10.0 << ' ' << .9 << endl;
2833
2834 // close it and try to re-open
2835
2836 reflfile_write.close();
2837 reflfile.close();
2838 reflfile.open( REFLECTIVITY_FILE );
2839
2840 }
2841
2842 if ( verbose >= VERBOSE_MINIMAL ) {
2843 if ( errors > 0 )
2844 log( "read_reflectivity", "Reading reflectivity file after %d fail%c\n",
2845 errors, ((errors>1)?'s':' ') );
2846 else
2847 log( "read_reflectivity", "Reading reflectivity file . . .\n" );
2848 }
2849
2850 // we set this to -1, in order to read in the first place the
2851 // number of data-points for the Reflectivity curve
2852 i = -1;
2853
2854 // scan the whole file
2855
2856 while ( ! reflfile.eof() ) {
2857
2858 // get line from the file
2859
2860 reflfile.getline(line, LINE_MAX_LENGTH);
2861
2862 // skip if comment
2863
2864 if ( *line == '#' )
2865 continue;
2866
2867 // get the value
2868
2869 if ( i < 0 ) { // if it is the first number we read
2870
2871 // read the number of datapoints
2872
2873 sscanf(line, "%d", &nReflectivity);
2874
2875 // allocate memory for Reflectivity table
2876
2877 Reflectivity = new float * [2];
2878 Reflectivity[0] = new float[nReflectivity];
2879 Reflectivity[1] = new float[nReflectivity];
2880
2881 } else { // if not, then it is a datapoint
2882
2883 // get the datapoint (wavelength, reflectivity)
2884
2885 sscanf(line, "%f %f", &Reflectivity[0][i], &Reflectivity[1][i]);
2886
2887 }
2888
2889 ++i;
2890
2891 }
2892
2893 // close it
2894
2895 reflfile.close();
2896
2897}
2898//!@}
2899
2900
2901//!---------------------------------------------------------------------
2902// @name read_axisdev
2903//
2904// @desc read axis deviations data for the mirrors
2905//
2906// @date Sat Jun 27 05:58:56 MET DST 1998
2907//----------------------------------------------------------------------
2908// @function
2909
2910//!@{
2911void
2912read_axisdev(void)
2913{
2914 ifstream axisdevfile;
2915 ofstream axisdevfile_write;
2916 int i, errors=0;
2917 char line[LINE_MAX_LENGTH];
2918
2919 // try to open the file
2920 axisdevfile.open( AXISDEVIATION_FILE );
2921
2922 // if it is wrong or does not exist:
2923
2924 while ( axisdevfile.bad() ) {
2925
2926 ++errors;
2927
2928 if ( verbose >= VERBOSE_MINIMAL )
2929 log( "read_axisdev", "Cannot open axis deviation file: %s\n",
2930 AXISDEVIATION_FILE );
2931
2932 if ( errors > 5 )
2933 error( "read_axisdev", "Exiting.");
2934
2935 // try to re-generate the file
2936
2937 if ( verbose >= VERBOSE_MINIMAL )
2938 log( "read_axisdev", "Generating file: %s\n",
2939 AXISDEVIATION_FILE );
2940
2941 // try to open the file
2942 axisdevfile_write.open( AXISDEVIATION_FILE );
2943
2944 // write short comment in the beginning
2945
2946 axisdevfile_write << "# axis deviation file" << endl
2947 << "# axis deviation for each mirror" << endl
2948 << "# J C Gonzalez, 1998" << endl
2949 << "#" << endl;
2950
2951 // write data
2952
2953 for ( i=0; i<ct_NMirrors; ++i )
2954 axisdevfile_write << gennor(0.0,ct_Adjustment_std)
2955 << ' '
2956 << gennor(0.0,ct_Adjustment_std)
2957 << endl;
2958
2959 // close it and try to re-open
2960
2961 axisdevfile_write.close();
2962 axisdevfile.close();
2963 axisdevfile.open( AXISDEVIATION_FILE );
2964
2965 }
2966
2967 if ( verbose >= VERBOSE_MINIMAL ) {
2968 if ( errors > 0 )
2969 log( "read_axisdev", "Reading axis deviation file after %d fail%c\n",
2970 errors, ((errors>1)?'s':' ') );
2971 else
2972 log( "read_axisdev", "Reading axis deviation file . . .\n" );
2973 }
2974
2975 //allocate memory for AxisDeviation table
2976
2977 AxisDeviation = new float * [2];
2978 AxisDeviation[0] = new float[ct_NMirrors];
2979 AxisDeviation[1] = new float[ct_NMirrors];
2980
2981 // scan the whole file
2982
2983 i = 0;
2984
2985 while ( (! axisdevfile.eof()) && (i<ct_NMirrors) ) {
2986
2987 // get line from the file
2988
2989 axisdevfile.getline(line, LINE_MAX_LENGTH);
2990
2991 // skip if comment
2992
2993 if ( *line == '#' )
2994 continue;
2995
2996 // get the value (dx, dy)
2997
2998 sscanf(line, "%f %f", &AxisDeviation[0][i], &AxisDeviation[1][i]);
2999
3000 ++i;
3001
3002 }
3003
3004 // did i do it for all the mirrors ?
3005
3006 if (i < ct_NMirrors) {
3007 error("read_axis_dev", "%s%s (%d < %d)\n",
3008 "Number of pairs in axisdev.dat file is",
3009 "\n\t\tsmaller than number of mirrors", i, ct_NMirrors);
3010 }
3011
3012 // close it
3013
3014 axisdevfile.close();
3015
3016}
3017//!@}
3018
3019
3020//!---------------------------------------------------------------------
3021// @name read_focals
3022//
3023// @desc read focals data for the mirrors
3024//
3025// @date Sat Jun 27 05:58:56 MET DST 1998
3026//----------------------------------------------------------------------
3027// @function
3028
3029//!@{
3030void
3031read_focals(void)
3032{
3033 ifstream focalfile;
3034 ofstream focalfile_write;
3035 int i, errors=0;
3036 char line[LINE_MAX_LENGTH];
3037
3038 /*!@'
3039
3040 Here we read the focals of each mirror, but only for CT1. For
3041 MAGIC we have already the Focals in the definition file. So, we
3042 copy these values to the |ct\_Focals| vector.
3043
3044 */
3045
3046 // >>>>>>>>>> MAGIC <<<<<<<<<<
3047 if ( ct_Type == 1) {
3048
3049 // allocate memory for the focals table
3050 ct_Focal = new float[ct_NMirrors];
3051
3052 for (i=0; i<ct_NMirrors; ++i)
3053 ct_Focal[i] = ct_data[i][CT_FOCAL];
3054
3055 return;
3056
3057 }
3058 // ! >>>>>>>>>> MAGIC <<<<<<<<<<
3059
3060 // try to open the file
3061 focalfile.open( FOCALS_FILE );
3062
3063 // if it is wrong or does not exist:
3064
3065 while ( focalfile.bad() ) {
3066
3067 ++errors;
3068
3069 if ( verbose >= VERBOSE_MINIMAL )
3070 log( "read_focals", "Cannot open focals file: %s\n",
3071 FOCALS_FILE );
3072
3073 if ( errors > 5 )
3074 error( "read_focals", "Exiting.");
3075
3076 // try to re-generate the file
3077
3078 if ( verbose >= VERBOSE_MINIMAL )
3079 log( "read_focals", "Generating file: %s\n",
3080 FOCALS_FILE );
3081
3082 // try to open the file
3083 focalfile_write.open( FOCALS_FILE );
3084
3085 // write short comment in the beginning
3086
3087 focalfile_write << "# focals file" << endl
3088 << "# focals for each mirror in the frame" << endl
3089 << "# J C Gonzalez, 1998" << endl;
3090
3091 // write data
3092
3093 rnormal( NormalRandomNumbers, ct_NMirrors+1 );
3094
3095 for (i=0; i<ct_NMirrors; i++)
3096 focalfile_write << NormalRandomNumbers[i] * ct_Focal_std + ct_Focal_mean
3097 << endl;
3098
3099 // close it and try to re-open
3100
3101 focalfile_write.close();
3102 focalfile.close();
3103 focalfile.open( FOCALS_FILE );
3104
3105 }
3106
3107 // allocate memory for the focals table
3108
3109 ct_Focal = new float[ct_NMirrors];
3110
3111 if ( verbose >= VERBOSE_MINIMAL ) {
3112
3113 if ( errors > 0 )
3114 log( "read_focals", "Reading focals file after %d fail%c\n",
3115 errors, ((errors>1) ? 's' : ' ') );
3116
3117 else
3118
3119 log( "read_focals", "Reading focals file . . .\n" );
3120 }
3121
3122 // scan the whole file
3123
3124 i = 0;
3125
3126 while ( (! focalfile.eof()) && (i<ct_NMirrors) ) {
3127
3128 // get line from the file
3129
3130 focalfile.getline(line, LINE_MAX_LENGTH);
3131
3132 // skip if comment
3133
3134 if ( *line == '#' )
3135 continue;
3136
3137 // get the value
3138
3139 sscanf(line, "%f", &ct_Focal[i]);
3140
3141 ++i;
3142
3143 }
3144
3145 // did i do it for all the mirrors ?
3146
3147 if (i < ct_NMirrors) {
3148 error("read_focals", "%s%s (%d < %d)\n",
3149 "Number of data in focals.dat file is",
3150 "\n\t\tsmaller than number of mirrors", i, ct_NMirrors);
3151 }
3152
3153 // close it
3154
3155 focalfile.close();
3156
3157}
3158//!@}
3159
3160
3161//!---------------------------------------------------------------------
3162// @name rnormal
3163//
3164// @desc returns n(=2k) normaly distributed numbers
3165//
3166// @var *r pointer to a vector where we write the numbers
3167// @var n how many numbers do we generate
3168//
3169// @date Sat Jun 27 05:58:56 MET DST 1998
3170//----------------------------------------------------------------------
3171// @function
3172
3173//!@{
3174void
3175rnormal(double *r, int n)
3176{
3177
3178 double z1, z2;
3179 int i;
3180
3181 for (i=0; i<n; i+=2) {
3182
3183 z1 = RandomNumber;
3184 z2 = RandomNumber;
3185
3186 r[i] = sqrt(-2.0*log(z1)) * cos(2.0*M_PI*z2);
3187 r[i+1] = sqrt(-2.0*log(z1)) * sin(2.0*M_PI*z2);
3188
3189 }
3190
3191}
3192//!@}
3193
3194
3195//!---------------------------------------------------------------------
3196// @name isA
3197//
3198// @desc returns TRUE(FALSE), if the flag is(is not) the given
3199//
3200// @var s1 String to be compared with s2
3201// @var s2 String in the form of a FLAG
3202//
3203// @return TRUE: both strings match; FALSE: otherwise
3204//
3205// @date Wed Jul 8 15:25:39 MET DST 1998
3206//----------------------------------------------------------------------
3207// @function
3208
3209//!@{
3210int
3211isA( char * s1, const char * flag )
3212{
3213 return ( (strncmp((char *)s1, flag, 8)==0) ? TRUE : FALSE );
3214}
3215//!@}
3216
3217
3218//!---------------------------------------------------------------------
3219// @name Curv2Lin
3220//
3221// @desc Curvilinear to Linear (Euclidean) distance
3222//
3223// @var s Curvilinear distance over the parabolic shape
3224//
3225// @return Radial distance from the axis of the paraboloid
3226//
3227// @date Wed Jul 8 15:25:39 MET DST 1998
3228//----------------------------------------------------------------------
3229// @function
3230
3231//!@{
3232float
3233Curv2Lin(float s)
3234{
3235 float x;
3236 short i;
3237
3238 x = s;
3239 for (i = 0; i < 4; i++)
3240 x = (s / 100.) / (1. + (float) 0.000144175317185 * x * x);
3241
3242 return (x*100.);
3243}
3244//!@}
3245
3246
3247//!---------------------------------------------------------------------
3248// @name Lin2Curv
3249//
3250// @desc Linear (Euclidean) to Curvilinear distance
3251//
3252// @var x Radial distance from the axis of the paraboloid
3253//
3254// @return Curvilinear distance over the parabolic shape
3255//
3256// @date Wed Jul 8 15:25:39 MET DST 1998
3257//----------------------------------------------------------------------
3258// @function
3259
3260//!@{
3261float
3262Lin2Curv(float x)
3263{
3264 x /= 100.;
3265 return ((x + (float) 0.000144175317185 * x * x * x)*100.);
3266}
3267//!@}
3268
3269
3270//!---------------------------------------------------------------------
3271// @name get_stdin_files
3272//
3273// @desc get data from stdin and create cerXXXXXX and staXXXXXX
3274//
3275// @var ncerf number of next Cherenkov file cerXXXXXX to be created
3276// @var El Lower limit for requested energy
3277// @var Eu Upper limit for requested energy
3278// @var flag Boolean flag: TRUE: initialize; FALSE: normal
3279//
3280// @return TRUE: succesful; FALSE: failure
3281//
3282// @date Tue Dec 15 15:03:40 MET 1998
3283//----------------------------------------------------------------------
3284// @function
3285
3286//!@{
3287int
3288get_stdin_files(int ncerf, float El, float Eu, int flag)
3289{
3290 static float *buffer; // buffer for read-out of STDIN
3291 static COREventHeader *evth;// Event Header class (from CORSIKA)
3292 static float Elow, Eup; // range for Energy
3293 static int reject;
3294
3295 char *chkbuffer;
3296
3297 char cername[256]; // output filename
3298 ofstream cerfile; // output file (stream)
3299 char staname[256]; // output filename sta
3300 ofstream stafile; // output file sta (stream)
3301 int bytes;
3302 float energy;
3303
3304 /*!@'
3305
3306 The main features of this function, as well as the algorithm,
3307 are taken from |corfilter|.
3308
3309 */
3310
3311 // if it's the first time, allocate buffers
3312
3313 if ( flag == TRUE ) {
3314
3315 if (verbose)
3316 log("get_stdin_files", "Allocating memory for buffer...\n");
3317
3318 buffer = new float [BUFFER_LENGTH];
3319 evth = (COREventHeader*)buffer;
3320
3321 Elow = El;
3322 Eup = Eu;
3323
3324 return ( TRUE );
3325 }
3326
3327 // main loop
3328
3329 while ( ! cin.eof() ) {
3330
3331 // read buffer from STDIN
3332
3333 cin.read( (char*) buffer, BUFFER_LENGTH * sizeof(float) );
3334 bytes = cin.gcount();
3335
3336 if (bytes < (BUFFER_LENGTH * sizeof(float))) {
3337 if (verbose)
3338 logerr("get_stdin_files",
3339 "Can read only %d bytes, instead of required %d.\n",
3340 bytes, BUFFER_LENGTH * sizeof(float));
3341 }
3342
3343 if (bytes == 0)
3344 return ( FALSE );
3345
3346 // check whether EVTH is there
3347
3348 chkbuffer = (char*)buffer;
3349 if (strstr(chkbuffer,EVTH) == 0) {
3350 if (verbose)
3351 logerr("get_stdin_files", "EVTH not found in this block\n");
3352 continue;
3353 } else {
3354 if (verbose)
3355 logerr("get_stdin_files", "EVTH was FOUND in this block!!\n");
3356 }
3357
3358
3359 // show log
3360 if (verbose)
3361 logerr("get_stdin_files", " (E=%8.1f) Reading event %d... ",
3362 (float)evth->Etotal, (int)evth->EvtNumber);
3363
3364 // CRITERIA
3365 energy = evth->Etotal;
3366
3367 if ( (energy<Elow) || (energy>Eup) ) {
3368 if (verbose)
3369 log("REJECTED", "\n");
3370 reject = TRUE;
3371 } else {
3372 if (verbose)
3373 log("ACCEPTED", "\n");
3374 reject = FALSE;
3375 }
3376
3377 if ( !reject ) {
3378
3379 // update no. of cer. file
3380 sprintf(cername, "%s/cer%06d", TMP_STDIN_DIR, ncerf+1);
3381
3382 cerfile.open( cername );
3383
3384 if ( verbose )
3385 log("get_stdin_files", " Writing event %d, will be %s\n",
3386 (int)evth->EvtNumber, cername);
3387
3388 // write it
3389 cerfile.write( (char*)buffer, BUFFER_LENGTH * sizeof(float) );
3390
3391 }
3392
3393 // read blocks till new EVTH is found
3394 // (now must be at the beginning of a block)
3395
3396 if (verbose)
3397 log("get_stdin_files", " Reading photons...\n");
3398
3399 bytes = 0;
3400 while ( TRUE ) {
3401
3402 // read photon
3403 cin.read( (char*)buffer, 7 * sizeof(float) );
3404
3405 // the last photon is reached when the EVTH is read, but
3406 // only when bytes % SIZE_OF_BLOCK == 0
3407 // SIZE_OF_BLOCK = 22932
3408 // with this, we only execute strstr once each 819 photons
3409
3410 // are we in a new block?
3411
3412 if ((bytes % SIZE_OF_BLOCK) == 0) {
3413 // is this the last photon?
3414 if ( strstr(chkbuffer, EVTH) == chkbuffer )
3415 break;
3416 }
3417
3418 // if not, write it to the cer-file
3419 if ( !reject )
3420 cerfile.write( (char*)buffer, 7 * sizeof(float) );
3421
3422 bytes += 7 * sizeof(float);
3423
3424 }
3425
3426 if (verbose)
3427 log("get_stdin_files", " Finished: %d bytes\n", bytes);
3428
3429 if ( !reject ) {
3430 // close cer-file
3431 cerfile.close();
3432 }
3433
3434 // at this point, we must be reading sta-file
3435
3436 // read the rest of the sta-file
3437 if (verbose)
3438 log("get_stdin_files", " Reading sta-data...\n");
3439 cin.read( (char*)buffer + 7*sizeof(float), 5956 - 7*sizeof(float));
3440
3441 if ( !reject ) {
3442
3443 // log
3444 if (verbose)
3445 log("get_stdin_files", " Saving sta-file...\n");
3446
3447 // open file
3448 sprintf(staname, "%s/sta%06d", TMP_STDIN_DIR, ncerf+1);
3449 stafile.open( staname );
3450
3451 // write it
3452 stafile.write( (char*)buffer, 5956 );
3453
3454 // close sta-file
3455 stafile.close();
3456
3457 // get out from here, return to main program
3458 return( TRUE );
3459
3460 }
3461
3462 }
3463
3464 return( TRUE );
3465
3466}
3467//!@}
3468
3469
3470//!---------------------------------------------------------------------
3471// @name get_new_ct_pointing
3472//
3473// @desc returns new random position of CT around a given one
3474//
3475// @var theta Theta (ZA) of the shower
3476// @var phi Phi (AZ) of the shower
3477// @var range Maximum allowed distance between CT and shower axis
3478// @var newtheta Resulting Theta of the CT
3479// @var newphi Resulting Phi of the CT
3480//
3481// @return Angular distance between orig. and new directions
3482//
3483// @date Sat Jun 27 05:58:56 MET DST 1998
3484//----------------------------------------------------------------------
3485// @function
3486
3487//!@{
3488float
3489get_new_ct_pointing(float theta, float phi, float range,
3490 float *newtheta, float *newphi)
3491{
3492 float distance;
3493 float it, ip,Nit;
3494 float sin_theta, cos_theta;
3495 float sin_newtheta, cos_newtheta;
3496 float sin_iphi, cos_iphi;
3497 float iphi;
3498 float range_aux;
3499 int i;
3500
3501 //Random position of the CT. Theta distribution
3502 //follows sin(theta) and ohi distribution is flat
3503
3504 range_aux=1.0-sin(M_PI/2-range);
3505 Nit = RandomNumber * range_aux;
3506 ip = RandomNumber * 2.0 * M_PI;
3507 it = M_PI/2-asin((1.0-Nit));
3508
3509 sin_theta = sin(theta);
3510 cos_theta = cos(theta);
3511
3512 cos_newtheta = cos_theta*cos(it) + sin_theta*sin(it)*cos(ip);
3513 *newtheta = acos( cos_newtheta );
3514 sin_newtheta = sin( *newtheta );
3515
3516 if(theta!=0){
3517 sin_iphi = sin(it)*sin(ip)/sin_newtheta;
3518 cos_iphi = ((cos(it)-cos_newtheta*cos_theta)/(sin_newtheta*sin_theta));
3519
3520 iphi = atan2( sin_iphi, cos_iphi );
3521 }
3522 else iphi=ip;
3523
3524 *newphi = phi + iphi;
3525
3526 return( it );
3527}
3528//!@}
3529
3530
3531//=------------------------------------------------------------
3532//!@subsection Log of this file.
3533
3534//!@{
3535//
3536// Revision 1.3 2000/01/03 12:41:26 harald
3537// There was a small mistake with +/- signs in the prevois versions. This
3538// line was founded by Jose Carlos. His remark is following:
3539//
3540// Look for the comment I include here in this block: you will see that
3541// now the following expression has a +- sign, where before it was only
3542// + (or -, I don't remember, but only one)
3543//
3544// Revision 1.2 1999/11/01 11:05:54 harald
3545// Small changes to comile the reflector program under linux.
3546// (Different use of NULL on DECalphas-osf1 and on linux)
3547//
3548// Revision 1.1.1.1 1999/10/29 07:00:33 harald
3549// This is the startpoint for the futher development of the Reflector program
3550// of Jose Carlos. For all developments use this CVS-controlled directory.
3551//
3552// Revision 1.18 1999/10/05 11:11:12 gonzalez
3553// Sep.1999
3554//
3555// Revision 1.17 1999/03/24 16:33:01 gonzalez
3556// REFLECTOR 1.1: Release
3557//
3558// Revision 1.16 1999/01/21 16:03:44 gonzalez
3559// Only small modifications
3560//
3561// Revision 1.15 1999/01/19 18:07:16 gonzalez
3562// Bugs in STDIN-STDOUT version corrected.
3563//
3564// Revision 1.14 1999/01/14 17:35:47 gonzalez
3565// Both reading from STDIN (data_from_stdin) and
3566// writing to STDOUT (data_to_STDOUT) working.
3567//
3568// Revision 1.13 1999/01/13 12:41:12 gonzalez
3569// THIS IS A WORKING (last?) VERSION
3570//
3571// Revision 1.12 1998/12/17 16:26:09 gonzalez
3572// *************************************************
3573// *************************************************
3574// WARNING:: Version 1.11 is completly wrong!!
3575// *************************************************
3576// *************************************************
3577//
3578// Revision 1.10 1998/12/15 10:47:22 gonzalez
3579// RELEASE-1.0
3580//
3581// Revision 1.9 1998/11/25 16:30:47 gonzalez
3582// Commit after inclusion of 'Blocking'
3583//
3584//!@}
3585
3586//=EOF
Note: See TracBrowser for help on using the repository browser.