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

Last change on this file since 347 was 347, checked in by harald, 25 years ago
smaller changes
File size: 95.9 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.8 $
23// $Author: harald $
24// $Date: 2000-01-28 08:59:49 $
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
1993 // add path from the mirror till the camera
1994
1995 t = t + sqrt( SQR(xcutCT[0] - xcam[0]) +
1996 SQR(xcutCT[1] - xcam[1]) +
1997 SQR(xcutCT[2] - xcam[2]) ) / Speed_of_Light_air_cmns;
1998
1999 // save data in the photon variable
2000
2001 cphoton.fill( photon.get_id(),
2002 xcam[0], xcam[1],
2003 photon.get_u(), photon.get_v(),
2004 photon.get_h(), t, phi);
2005
2006 // show it
2007
2008 if ( verbose >= VERBOSE_MAXIMAL ) {
2009 cout << '@' << '2'
2010 << ' ' << ncer
2011 << ' ' << xCT[0]
2012 << ' ' << xCT[1]
2013 << ' ' << xCT[2]
2014 << ' ' << rCT[0]
2015 << ' ' << rCT[1]
2016 << ' ' << rCT[2]
2017 << ' ' << xcut[0]
2018 << ' ' << xcut[1]
2019 << ' ' << xcut[2]
2020 << endl;
2021 cout << '@' << '3'
2022 << ' ' << ncer
2023 << ' ' << sx
2024 << ' ' << sy
2025 << ' ' << i_mirror
2026 << ' ' << ct_data[i_mirror][CT_SX]
2027 << ' ' << ct_data[i_mirror][CT_SY]
2028 << ' ' << ct_data[i_mirror][CT_SX] - sx
2029 << ' ' << ct_data[i_mirror][CT_SY] - sy
2030 << endl;
2031 if ( pb_ParallelBeam == FALSE ) {
2032 cout << '@' << '4'
2033 << ' ' << ncer
2034 << ' ' << xcam[0]
2035 << ' ' << xcam[1]
2036 << ' ' << xcam[2]
2037 << ' ' << phi
2038 << ' ' << photon.get_t()-stat.get_tfirst()
2039 << ' ' << t-stat.get_tfirst()
2040 << endl << flush;
2041 } else {
2042 cout << '@' << '4'
2043 << ' ' << ncer
2044 << ' ' << xcam[0]
2045 << ' ' << xcam[1]
2046 << ' ' << xcam[2]
2047 << ' ' << phi
2048 << endl << flush;
2049 }
2050 }
2051
2052 // write the photon data to the output file
2053
2054 if ( Data_To_STDOUT )
2055 cout.write ( (char *)&cphoton, cphoton.mysize() );
2056 else
2057 cphoton.write( outputfile );
2058
2059 // increase number of Cphotons written
2060
2061 ++nCphotons;
2062
2063 } // while still there are photons left
2064
2065 if ( verbose >= VERBOSE_NORMAL )
2066 cout << nCphotons << ' '
2067 << nbeforeTR << ' '
2068 << nafterTR << ' '
2069 << endl << flush;
2070
2071 // write end-of-event mark
2072
2073 if ( Data_To_STDOUT )
2074 cout.write( FLAG_END_OF_EVENT, SIZE_OF_FLAGS );
2075 else
2076 outputfile.write( FLAG_END_OF_EVENT, SIZE_OF_FLAGS );
2077
2078 // close files
2079
2080 cerfile.close();
2081 stat.closefile();
2082
2083 // show how many photons were written
2084
2085 if ( verbose >= VERBOSE_MAXIMAL )
2086 log( SIGNATURE, "%d C-photons written.\n", nCphotons );
2087
2088 // if we are reading data from STDIN, delete the last used files
2089
2090 if ( Data_From_STDIN ) {
2091
2092 unlink( cername );
2093 unlink( staname );
2094
2095 }
2096
2097 } // while there are still data left
2098
2099 // write end-of-run mark
2100
2101 if ( Data_To_STDOUT )
2102 cout.write( FLAG_END_OF_RUN, SIZE_OF_FLAGS );
2103 else
2104 outputfile.write( FLAG_END_OF_RUN, SIZE_OF_FLAGS );
2105
2106 // close directory
2107
2108 closedir( directory );
2109
2110 } // for each data directory
2111
2112 // write end-of-file mark
2113
2114 if ( Data_To_STDOUT )
2115 cout.write( FLAG_END_OF_FILE, SIZE_OF_FLAGS );
2116 else
2117 outputfile.write( FLAG_END_OF_FILE, SIZE_OF_FLAGS );
2118
2119 // close output file
2120
2121 if ( verbose >= VERBOSE_MINIMAL )
2122 log( SIGNATURE, "Closing output file %s\n", outname );
2123
2124 if ( ! Data_To_STDOUT )
2125 outputfile.close();
2126
2127 // clean everything
2128
2129 if ( verbose >= VERBOSE_MINIMAL )
2130 log( SIGNATURE, "Cleanning . . .\n" );
2131
2132 clean();
2133
2134 // program finished
2135
2136 if ( verbose >= VERBOSE_MINIMAL )
2137 log( SIGNATURE, "Done.\n");
2138
2139 return ( 0 );
2140}
2141//!@}
2142
2143// @T \newpage
2144
2145//=---------------------------------------------------------------------
2146// @subsection Functions definition
2147
2148//!---------------------------------------------------------------------
2149// @name present
2150//
2151// @desc Make some presentation
2152//
2153// @date Sat Jun 27 05:58:56 MET DST 1998
2154//----------------------------------------------------------------------
2155// @function
2156
2157//!@{
2158void
2159present(void)
2160{
2161 cout << "##################################################\n"
2162 << SIGNATURE << '\n' << '\n'
2163 << "Simulation of the reflector\n"
2164 << "J C Gonzalez, May 1998\n"
2165 << "##################################################\n\n";
2166}
2167//!@}
2168
2169
2170//!---------------------------------------------------------------------
2171// @name usage
2172//
2173// @desc show help
2174//
2175// @date Tue Dec 15 16:23:30 MET 1998
2176//----------------------------------------------------------------------
2177// @function
2178
2179//!@{
2180void
2181usage(void)
2182{
2183 present();
2184 cout << "\nusage ::\n\n"
2185 << "\t reflector "
2186 << " [ -@ paramfile ] "
2187 << " [ -h ] "
2188 << "\n\n or \n\n"
2189 << "\t reflector < paramfile"
2190 << "\n\n";
2191 exit(0);
2192}
2193//!@}
2194
2195
2196//!---------------------------------------------------------------------
2197// @name clean
2198//
2199// @desc basically, frees memory
2200//
2201// @date Sat Jun 27 05:58:56 MET DST 1998
2202//----------------------------------------------------------------------
2203// @function
2204
2205//!@{
2206void
2207clean(void)
2208{
2209 int i;
2210
2211 // delete focals table
2212
2213 delete [] ct_Focal;
2214 ct_Focal = NULL;
2215
2216 // delete reflectivity table
2217
2218 for (i=0; i<2; i++) {
2219 delete [] Reflectivity[i];
2220 }
2221
2222 delete [] Reflectivity;
2223 Reflectivity = NULL;
2224
2225 // delete mirrors' data table
2226
2227 for (i=0; i<ct_NMirrors; i++) {
2228 delete [] ct_data[i];
2229 }
2230
2231 delete [] ct_data;
2232 ct_data = NULL;
2233}
2234//!@}
2235
2236
2237//!---------------------------------------------------------------------
2238// @name log
2239//
2240// @desc function to send log information
2241//
2242// @var funct Name of the caller function
2243// @var fmt Format to be used (message)
2244// @var ... Other information to be shown
2245//
2246// @date Sat Jun 27 05:58:56 MET DST 1998
2247//----------------------------------------------------------------------
2248// @function
2249
2250//!@{
2251void
2252log(const char *funct, char *fmt, ...)
2253{
2254 va_list args;
2255
2256 // Display the name of the function that called error
2257 printf("[%s]: ", funct);
2258
2259 // Display the remainder of the message
2260 va_start(args, fmt);
2261 vprintf(fmt, args);
2262 va_end(args);
2263}
2264//!@}
2265
2266
2267//!---------------------------------------------------------------------
2268// @name logerr
2269//
2270// @desc function to send log information to stderr
2271//
2272// @var funct Name of the caller function
2273// @var fmt Format to be used (message)
2274// @var ... Other information to be shown
2275//
2276// @date Sat Jun 27 05:58:56 MET DST 1998
2277//----------------------------------------------------------------------
2278// @function
2279
2280//!@{
2281void
2282logerr(const char *funct, char *fmt, ...)
2283{
2284 va_list args;
2285
2286 // Display the name of the function that called error
2287 printf("[%s]: ", funct);
2288
2289 // Display the remainder of the message
2290 va_start(args, fmt);
2291 vfprintf(stderr, fmt, args);
2292 va_end(args);
2293}
2294//!@}
2295
2296
2297//!---------------------------------------------------------------------
2298// @name error
2299//
2300// @desc function to send an error message, and abort the program
2301//
2302// @var funct Name of the caller function
2303// @var fmt Format to be used (message)
2304// @var ... Other information to be shown
2305//
2306// @date Sat Jun 27 05:58:56 MET DST 1998
2307//----------------------------------------------------------------------
2308// @function
2309
2310//!@{
2311void
2312error(const char *funct, char *fmt, ...)
2313{
2314 va_list args;
2315
2316 // Display the name of the function that called error
2317 fprintf(stderr, "ERROR in %s: ", funct);
2318
2319 // Display the remainder of the message
2320 va_start(args, fmt);
2321 vfprintf(stderr, fmt, args);
2322 va_end(args);
2323
2324 abort();
2325}
2326//!@}
2327
2328
2329//!---------------------------------------------------------------------
2330// @name makeOmega
2331//
2332// @desc function to calculate the matrix Omega(theta,phi)
2333//
2334// @var theta Angle theta of the transformation
2335// @var phi Angle phi of the transformation
2336//
2337// @date Sat Jun 27 05:58:56 MET DST 1998
2338//----------------------------------------------------------------------
2339// @function
2340
2341//!@{
2342void
2343makeOmega (float theta, float phi)
2344{
2345 static float ct, st, cp, sp;
2346
2347 // shortcuts for cosine and sine of theta and phi
2348 ct = cos(theta);
2349 st = sin(theta);
2350 cp = cos(phi);
2351 sp = sin(phi);
2352
2353 // save values in the array (see top of file)
2354 Omega[0][0] = cp*ct;
2355 Omega[0][1] = sp*ct;
2356 Omega[0][2] = -st;
2357
2358 Omega[1][0] = -sp;
2359 Omega[1][1] = cp;
2360 Omega[1][2] = 0;
2361
2362 Omega[2][0] = cp*st;
2363 Omega[2][1] = sp*st;
2364 Omega[2][2] = ct;
2365}
2366//!@}
2367
2368
2369//!---------------------------------------------------------------------
2370// @name makeOmegaI
2371//
2372// @desc function to calculate the matrix Omega-1(theta,phi)
2373//
2374// @var theta Angle theta of the transformation
2375// @var phi Angle phi of the transformation
2376//
2377// @date Sat Jun 27 05:58:56 MET DST 1998
2378//----------------------------------------------------------------------
2379// @function
2380
2381//!@{
2382void
2383makeOmegaI(float theta, float phi)
2384{
2385 static float ct, st, cp, sp;
2386
2387 // shortcuts for cosine and sine of theta and phi
2388 ct = cos(theta);
2389 st = sin(theta);
2390 cp = cos(phi);
2391 sp = sin(phi);
2392
2393 // save values in the array (see top of file)
2394 OmegaI[0][0] = cp*ct;
2395 OmegaI[0][1] = -sp;
2396 OmegaI[0][2] = cp*st;
2397
2398 OmegaI[1][0] = sp*ct;
2399 OmegaI[1][1] = cp;
2400 OmegaI[1][2] = sp*st;
2401
2402 OmegaI[2][0] = -st;
2403 OmegaI[2][1] = 0;
2404 OmegaI[2][2] = ct;
2405}
2406//!@}
2407
2408
2409//!---------------------------------------------------------------------
2410// @name applyMxv
2411//
2412// @desc returns the vector v' such that v' = M x v
2413//
2414// @var M matrix of the transformation
2415// @var v vector to be multiplied
2416// @var vi resulting vector
2417//
2418// @date Sat Jun 27 05:58:56 MET DST 1998
2419//----------------------------------------------------------------------
2420// @function
2421
2422//!@{
2423void
2424applyMxV(float M[3][3], float *V, float *Vp)
2425{
2426 Vp[0] = (M[0][0] * V[0] +
2427 M[0][1] * V[1] +
2428 M[0][2] * V[2]);
2429 Vp[1] = (M[1][0] * V[0] +
2430 M[1][1] * V[1] +
2431 M[1][2] * V[2]);
2432 Vp[2] = (M[2][0] * V[0] +
2433 M[2][1] * V[1] +
2434 M[2][2] * V[2]);
2435}
2436//!@}
2437
2438
2439//!---------------------------------------------------------------------
2440// @name read_ct_file
2441//
2442// @desc read CT definition file
2443//
2444// @date Sat Jun 27 05:58:56 MET DST 1998
2445//----------------------------------------------------------------------
2446// @function
2447
2448//!@{
2449void
2450read_ct_file(void)
2451{
2452 char line[LINE_MAX_LENGTH]; // line to get from the ctin
2453 char token[ITEM_MAX_LENGTH]; // a single token
2454 char fmirr[40]; // name of the binary version of mirrors data
2455 int i; // dummy counters
2456
2457 if ( verbose >= VERBOSE_MINIMAL )
2458 log( "read_ct_file", "start.\n" );
2459
2460 ifstream ctin ( get_ct_filename() );
2461
2462 if ( ctin.bad() )
2463 error( "read_ct_file",
2464 "Cannot open CT def. file: %s\n", get_ct_filename() );
2465
2466 // loop till the "end" directive is reached
2467
2468 while (!ctin.eof()) {
2469
2470 // get line from stdin
2471
2472 ctin.getline(line, LINE_MAX_LENGTH);
2473
2474 // look for each item at the beginning of the line
2475
2476 for (i=0; i<=define_mirrors; i++)
2477 if (strstr(line, CT_ITEM_NAMES[i]) == line)
2478 break;
2479
2480 // if it is not a valid line, just ignore it
2481
2482 if (i == define_mirrors+1)
2483 continue;
2484
2485 // case block for each directive
2486
2487 switch ( i ) {
2488
2489 case type: // <type of telescope> (0:CT1 ¦ 1:MAGIC)
2490
2491 // get focal distance
2492
2493 sscanf(line, "%s %d", token, &ct_Type);
2494
2495 if ( verbose >= VERBOSE_MINIMAL )
2496 log( "read_ct_file", "<Type of Telescope>: %s\n",
2497 ((ct_Type==0) ? "CT1" : "MAGIC") );
2498
2499 break;
2500
2501 case focal_distance: // <focal distance> [cm]
2502
2503 // get focal distance
2504
2505 sscanf(line, "%s %f", token, &ct_Focal_mean);
2506
2507 if ( verbose >= VERBOSE_MINIMAL )
2508 log( "read_ct_file", "<Focal distance>: %f cm\n", ct_Focal_mean );
2509
2510 break;
2511
2512 case focal_std: // s(focal distance) [cm]
2513
2514 // get focal distance
2515
2516 sscanf(line, "%s %f", token, &ct_Focal_std);
2517
2518 if ( verbose >= VERBOSE_MINIMAL )
2519 log( "read_ct_file", "s(Focal distance): %f cm\n", ct_Focal_std );
2520
2521 break;
2522
2523 case point_spread: // <point spread> [cm]
2524
2525 // get point spread
2526
2527 sscanf(line, "%s %f", token, &ct_PSpread_mean);
2528
2529 if ( verbose >= VERBOSE_MINIMAL )
2530 log( "read_ct_file", "<Point spread>: %f cm\n", ct_PSpread_mean );
2531
2532 break;
2533
2534 case point_std: // s(point spread) [cm]
2535
2536 // get point spread
2537
2538 sscanf(line, "%s %f", token, &ct_PSpread_std);
2539
2540 if ( verbose >= VERBOSE_MINIMAL )
2541 log( "read_ct_file", "s(Point spread): %f cm\n", ct_PSpread_std );
2542
2543 break;
2544
2545 case adjustment_dev: // s(adjustment_dev) [cm]
2546
2547 // get point spread
2548
2549 sscanf(line, "%s %f", token, &ct_Adjustment_std);
2550
2551 if ( verbose >= VERBOSE_MINIMAL )
2552 log( "read_ct_file", "s(Adjustment): %f cm\n", ct_Adjustment_std );
2553
2554 break;
2555
2556 case black_spot: // radius of the black spot in the center [cm]
2557
2558 // get black spot radius
2559
2560 sscanf(line, "%s %f", token, &ct_BlackSpot_rad);
2561
2562 if ( verbose >= VERBOSE_MINIMAL )
2563 log( "read_ct_file", "Radius of the black spots: %f cm\n",
2564 ct_BlackSpot_rad);
2565
2566 break;
2567
2568 case r_mirror: // radius of the mirrors [cm]
2569
2570 // get radius of mirror
2571
2572 sscanf(line, "%s %f", token, &ct_RMirror);
2573
2574 if ( verbose >= VERBOSE_MINIMAL )
2575 log( "read_ct_file", "Radii of the mirrors: %f cm\n", ct_RMirror );
2576
2577 break;
2578
2579 case n_mirrors: // number of mirrors
2580
2581 // get the name of the output_file from the line
2582
2583 sscanf(line, "%s %d", token, &ct_NMirrors);
2584
2585 if ( verbose >= VERBOSE_MINIMAL )
2586 log( "read_ct_file", "Number of mirrors: %d\n", ct_NMirrors );
2587
2588 break;
2589
2590 case camera_width: // camera width [cm]
2591
2592 // get the name of the ct_file from the line
2593
2594 sscanf(line, "%s %f", token, &ct_CameraWidth);
2595
2596 if ( verbose >= VERBOSE_MINIMAL )
2597 log( "read_ct_file", "Camera width: %f cm\n", ct_CameraWidth );
2598
2599 break;
2600
2601 case n_pixels: // number of pixels
2602
2603 // get the name of the output_file from the line
2604
2605 sscanf(line, "%s %d", token, &ct_NPixels);
2606
2607 if ( verbose >= VERBOSE_MINIMAL )
2608 log( "read_ct_file", "Number of pixels: %d\n", ct_NPixels );
2609
2610 break;
2611
2612 case pixel_width: // pixel width [cm]
2613
2614 // get the name of the ct_file from the line
2615
2616 sscanf(line, "%s %f", token, &ct_PixelWidth);
2617
2618 if ( verbose >= VERBOSE_MINIMAL )
2619 log( "read_ct_file", "Pixel width: %f cm\n", ct_PixelWidth );
2620
2621 break;
2622
2623 case define_mirrors: // read table with the parameters of the mirrors
2624
2625 if ( verbose >= VERBOSE_MINIMAL )
2626 log( "read_ct_file", "Table of mirrors data:\n" );
2627
2628 // check whether the number of mirrors was already set
2629
2630 if ( ct_NMirrors == 0 )
2631 error( "read_ct_file", "NMirrors was not set.\n" );
2632
2633 // allocate memory for paths list
2634
2635 if ( verbose >= VERBOSE_MINIMAL )
2636 log( "read_ct_file", "Allocating memory for ct_data\n" );
2637
2638 ct_data = new float*[ct_NMirrors];
2639
2640 for (i=0; i<ct_NMirrors; i++)
2641 ct_data[i] = new float[CT_NDATA];
2642
2643 // read data
2644
2645 // look for binary version of mirror data (faster to load)
2646
2647 sprintf( fmirr, "%s.mirr", get_ct_filename() );
2648
2649 ifstream ctmirr;
2650 ctmirr.open( fmirr );
2651
2652 if ( ctmirr.bad() ) {
2653
2654 // the file does not exist
2655
2656 ctmirr.close(); // close not existing file
2657
2658 ofstream ctmirr_write ( fmirr ); // open file to save the data
2659
2660 if ( verbose >= VERBOSE_MINIMAL )
2661 log( "read_ct_file", "Reading mirrors data...\n" );
2662
2663 // read ASCII data
2664
2665 for (i=0; i<ct_NMirrors; i++) {
2666 // for (j=0; j<CT_NDATA; j++)
2667
2668 ctin.getline(line, LINE_MAX_LENGTH);
2669 sscanf(line, "%f %f %f %f %f %f %f %f %f %f %f %f",
2670 &ct_data[i][0],
2671 &ct_data[i][1],
2672 &ct_data[i][2],
2673 &ct_data[i][3],
2674 &ct_data[i][4],
2675 &ct_data[i][5],
2676 &ct_data[i][6],
2677 &ct_data[i][7],
2678 &ct_data[i][8],
2679 &ct_data[i][9],
2680 &ct_data[i][10],
2681 &ct_data[i][11]);
2682
2683 cout << '[' << i << ']' << flush;
2684 }
2685 cout << endl << flush;
2686
2687 // save binary data and close output file
2688
2689 for (i=0; i<ct_NMirrors; i++)
2690 ctmirr_write.write( (char*)ct_data[i], CT_NDATA*sizeof(float) );
2691
2692 ctmirr_write.close();
2693
2694 } else {
2695
2696 // Excellent!
2697 // read binary data and close input file
2698
2699 for (i=0; i<ct_NMirrors; i++)
2700 ctmirr.read( (char*)ct_data[i], CT_NDATA*sizeof(float) );
2701
2702 ctmirr.close();
2703
2704 }
2705
2706 break;
2707
2708 } // switch ( i )
2709
2710 } // while (! is_end)
2711
2712 // read reflectivity file
2713
2714 read_reflectivity();
2715
2716 // read axis deviations file
2717
2718 read_axisdev();
2719
2720 // read focals file
2721
2722 read_focals();
2723
2724 // end
2725
2726 if ( verbose >= VERBOSE_MINIMAL )
2727 log( "read_ct_file", "done.\n" );
2728 return;
2729}
2730//!@}
2731
2732
2733//!---------------------------------------------------------------------
2734// @name dist_r_P
2735//
2736// @desc distance straight line r - point P
2737//
2738// @var a coord. X of a fixed point of the straight line
2739// @var b coord. Y of a fixed point of the straight line
2740// @var c coord. Z of a fixed point of the straight line
2741// @var l component X of a vector of the straight line
2742// @var m component Y of a vector of the straight line
2743// @var n component Z of a vector of the straight line
2744// @var x coord. X of the point P
2745// @var y coord. Y of the point P
2746// @var z coord. Z of the point P
2747//
2748// @return Returns the distance d(r,P) between the line and the point P
2749//
2750// @date Sat Jun 27 05:58:56 MET DST 1998
2751//----------------------------------------------------------------------
2752// @function
2753
2754//!@{
2755float
2756dist_r_P(float a, float b, float c,
2757 float l, float m, float n,
2758 float x, float y, float z)
2759{
2760 return (
2761 sqrt((SQR((a-x)*m-(b-y)*l) +
2762 SQR((b-y)*n-(c-z)*m) +
2763 SQR((c-z)*l-(a-x)*n))/
2764 (SQR(l)+SQR(m)+SQR(n))
2765 )
2766 );
2767}
2768//!@}
2769
2770
2771//!---------------------------------------------------------------------
2772// @name read_reflectivity
2773//
2774// @desc read reflectivity data for the mirrors
2775//
2776// @date Sat Jun 27 05:58:56 MET DST 1998
2777//----------------------------------------------------------------------
2778// @function
2779
2780//!@{
2781void
2782read_reflectivity(void)
2783{
2784 ifstream reflfile;
2785 ofstream reflfile_write;
2786 int i, errors=0;
2787 char line[LINE_MAX_LENGTH];
2788
2789 // try to open the file
2790 reflfile.open( REFLECTIVITY_FILE );
2791
2792 // if it is wrong or does not exist:
2793
2794 while ( reflfile.bad() ) {
2795
2796 ++errors;
2797
2798 if ( verbose >= VERBOSE_MINIMAL )
2799 log( "read_reflectivity", "Cannot open reflectivity file: %s\n",
2800 REFLECTIVITY_FILE );
2801
2802 if ( errors > 5 )
2803 error( "read_reflectivity", "Exiting.");
2804
2805 // try to re-generate the file
2806
2807 if ( verbose >= VERBOSE_MINIMAL )
2808 log( "read_reflectivity", "Generating file: %s\n",
2809 REFLECTIVITY_FILE );
2810
2811 // try to open the file
2812 reflfile_write.open( REFLECTIVITY_FILE );
2813
2814 // write short comment in the beginning
2815
2816 reflfile_write << "# reflectivity file" << endl
2817 << "# reflectivity for each mirror in the frame" << endl
2818 << "# J C Gonzalez, 1998" << endl
2819 << "#" << endl;
2820
2821 // write data
2822
2823 // write number of data-points
2824
2825 reflfile_write << "# number of datapoints" << endl
2826 << 35 << endl;
2827
2828 // write data-points
2829
2830 reflfile_write << "# datapoints" << endl;
2831
2832 for ( i=0; i<35; ++i )
2833 reflfile_write << 270.0+i*10.0 << ' ' << .9 << endl;
2834
2835 // close it and try to re-open
2836
2837 reflfile_write.close();
2838 reflfile.close();
2839 reflfile.open( REFLECTIVITY_FILE );
2840
2841 }
2842
2843 if ( verbose >= VERBOSE_MINIMAL ) {
2844 if ( errors > 0 )
2845 log( "read_reflectivity", "Reading reflectivity file after %d fail%c\n",
2846 errors, ((errors>1)?'s':' ') );
2847 else
2848 log( "read_reflectivity", "Reading reflectivity file . . .\n" );
2849 }
2850
2851 // we set this to -1, in order to read in the first place the
2852 // number of data-points for the Reflectivity curve
2853 i = -1;
2854
2855 // scan the whole file
2856
2857 while ( ! reflfile.eof() ) {
2858
2859 // get line from the file
2860
2861 reflfile.getline(line, LINE_MAX_LENGTH);
2862
2863 // skip if comment
2864
2865 if ( *line == '#' )
2866 continue;
2867
2868 // get the value
2869
2870 if ( i < 0 ) { // if it is the first number we read
2871
2872 // read the number of datapoints
2873
2874 sscanf(line, "%d", &nReflectivity);
2875
2876 // allocate memory for Reflectivity table
2877
2878 Reflectivity = new float * [2];
2879 Reflectivity[0] = new float[nReflectivity];
2880 Reflectivity[1] = new float[nReflectivity];
2881
2882 } else { // if not, then it is a datapoint
2883
2884 // get the datapoint (wavelength, reflectivity)
2885
2886 sscanf(line, "%f %f", &Reflectivity[0][i], &Reflectivity[1][i]);
2887
2888 }
2889
2890 ++i;
2891
2892 }
2893
2894 // close it
2895
2896 reflfile.close();
2897
2898}
2899//!@}
2900
2901
2902//!---------------------------------------------------------------------
2903// @name read_axisdev
2904//
2905// @desc read axis deviations data for the mirrors
2906//
2907// @date Sat Jun 27 05:58:56 MET DST 1998
2908//----------------------------------------------------------------------
2909// @function
2910
2911//!@{
2912void
2913read_axisdev(void)
2914{
2915 ifstream axisdevfile;
2916 ofstream axisdevfile_write;
2917 int i, errors=0;
2918 char line[LINE_MAX_LENGTH];
2919
2920 // try to open the file
2921 axisdevfile.open( AXISDEVIATION_FILE );
2922
2923 // if it is wrong or does not exist:
2924
2925 while ( axisdevfile.bad() ) {
2926
2927 ++errors;
2928
2929 if ( verbose >= VERBOSE_MINIMAL )
2930 log( "read_axisdev", "Cannot open axis deviation file: %s\n",
2931 AXISDEVIATION_FILE );
2932
2933 if ( errors > 5 )
2934 error( "read_axisdev", "Exiting.");
2935
2936 // try to re-generate the file
2937
2938 if ( verbose >= VERBOSE_MINIMAL )
2939 log( "read_axisdev", "Generating file: %s\n",
2940 AXISDEVIATION_FILE );
2941
2942 // try to open the file
2943 axisdevfile_write.open( AXISDEVIATION_FILE );
2944
2945 // write short comment in the beginning
2946
2947 axisdevfile_write << "# axis deviation file" << endl
2948 << "# axis deviation for each mirror" << endl
2949 << "# J C Gonzalez, 1998" << endl
2950 << "#" << endl;
2951
2952 // write data
2953
2954 for ( i=0; i<ct_NMirrors; ++i )
2955 axisdevfile_write << gennor(0.0,ct_Adjustment_std)
2956 << ' '
2957 << gennor(0.0,ct_Adjustment_std)
2958 << endl;
2959
2960 // close it and try to re-open
2961
2962 axisdevfile_write.close();
2963 axisdevfile.close();
2964 axisdevfile.open( AXISDEVIATION_FILE );
2965
2966 }
2967
2968 if ( verbose >= VERBOSE_MINIMAL ) {
2969 if ( errors > 0 )
2970 log( "read_axisdev", "Reading axis deviation file after %d fail%c\n",
2971 errors, ((errors>1)?'s':' ') );
2972 else
2973 log( "read_axisdev", "Reading axis deviation file . . .\n" );
2974 }
2975
2976 //allocate memory for AxisDeviation table
2977
2978 AxisDeviation = new float * [2];
2979 AxisDeviation[0] = new float[ct_NMirrors];
2980 AxisDeviation[1] = new float[ct_NMirrors];
2981
2982 // scan the whole file
2983
2984 i = 0;
2985
2986 while ( (! axisdevfile.eof()) && (i<ct_NMirrors) ) {
2987
2988 // get line from the file
2989
2990 axisdevfile.getline(line, LINE_MAX_LENGTH);
2991
2992 // skip if comment
2993
2994 if ( *line == '#' )
2995 continue;
2996
2997 // get the value (dx, dy)
2998
2999 sscanf(line, "%f %f", &AxisDeviation[0][i], &AxisDeviation[1][i]);
3000
3001 ++i;
3002
3003 }
3004
3005 // did i do it for all the mirrors ?
3006
3007 if (i < ct_NMirrors) {
3008 error("read_axis_dev", "%s%s (%d < %d)\n",
3009 "Number of pairs in axisdev.dat file is",
3010 "\n\t\tsmaller than number of mirrors", i, ct_NMirrors);
3011 }
3012
3013 // close it
3014
3015 axisdevfile.close();
3016
3017}
3018//!@}
3019
3020
3021//!---------------------------------------------------------------------
3022// @name read_focals
3023//
3024// @desc read focals data for the mirrors
3025//
3026// @date Sat Jun 27 05:58:56 MET DST 1998
3027//----------------------------------------------------------------------
3028// @function
3029
3030//!@{
3031void
3032read_focals(void)
3033{
3034 ifstream focalfile;
3035 ofstream focalfile_write;
3036 int i, errors=0;
3037 char line[LINE_MAX_LENGTH];
3038
3039 /*!@'
3040
3041 Here we read the focals of each mirror, but only for CT1. For
3042 MAGIC we have already the Focals in the definition file. So, we
3043 copy these values to the |ct\_Focals| vector.
3044
3045 */
3046
3047 // >>>>>>>>>> MAGIC <<<<<<<<<<
3048 if ( ct_Type == 1) {
3049
3050 // allocate memory for the focals table
3051 ct_Focal = new float[ct_NMirrors];
3052
3053 for (i=0; i<ct_NMirrors; ++i)
3054 ct_Focal[i] = ct_data[i][CT_FOCAL];
3055
3056 return;
3057
3058 }
3059 // ! >>>>>>>>>> MAGIC <<<<<<<<<<
3060
3061 // try to open the file
3062 focalfile.open( FOCALS_FILE );
3063
3064 // if it is wrong or does not exist:
3065
3066 while ( focalfile.bad() ) {
3067
3068 ++errors;
3069
3070 if ( verbose >= VERBOSE_MINIMAL )
3071 log( "read_focals", "Cannot open focals file: %s\n",
3072 FOCALS_FILE );
3073
3074 if ( errors > 5 )
3075 error( "read_focals", "Exiting.");
3076
3077 // try to re-generate the file
3078
3079 if ( verbose >= VERBOSE_MINIMAL )
3080 log( "read_focals", "Generating file: %s\n",
3081 FOCALS_FILE );
3082
3083 // try to open the file
3084 focalfile_write.open( FOCALS_FILE );
3085
3086 // write short comment in the beginning
3087
3088 focalfile_write << "# focals file" << endl
3089 << "# focals for each mirror in the frame" << endl
3090 << "# J C Gonzalez, 1998" << endl;
3091
3092 // write data
3093
3094 rnormal( NormalRandomNumbers, ct_NMirrors+1 );
3095
3096 for (i=0; i<ct_NMirrors; i++)
3097 focalfile_write << NormalRandomNumbers[i] * ct_Focal_std + ct_Focal_mean
3098 << endl;
3099
3100 // close it and try to re-open
3101
3102 focalfile_write.close();
3103 focalfile.close();
3104 focalfile.open( FOCALS_FILE );
3105
3106 }
3107
3108 // allocate memory for the focals table
3109
3110 ct_Focal = new float[ct_NMirrors];
3111
3112 if ( verbose >= VERBOSE_MINIMAL ) {
3113
3114 if ( errors > 0 )
3115 log( "read_focals", "Reading focals file after %d fail%c\n",
3116 errors, ((errors>1) ? 's' : ' ') );
3117
3118 else
3119
3120 log( "read_focals", "Reading focals file . . .\n" );
3121 }
3122
3123 // scan the whole file
3124
3125 i = 0;
3126
3127 while ( (! focalfile.eof()) && (i<ct_NMirrors) ) {
3128
3129 // get line from the file
3130
3131 focalfile.getline(line, LINE_MAX_LENGTH);
3132
3133 // skip if comment
3134
3135 if ( *line == '#' )
3136 continue;
3137
3138 // get the value
3139
3140 sscanf(line, "%f", &ct_Focal[i]);
3141
3142 ++i;
3143
3144 }
3145
3146 // did i do it for all the mirrors ?
3147
3148 if (i < ct_NMirrors) {
3149 error("read_focals", "%s%s (%d < %d)\n",
3150 "Number of data in focals.dat file is",
3151 "\n\t\tsmaller than number of mirrors", i, ct_NMirrors);
3152 }
3153
3154 // close it
3155
3156 focalfile.close();
3157
3158}
3159//!@}
3160
3161
3162//!---------------------------------------------------------------------
3163// @name rnormal
3164//
3165// @desc returns n(=2k) normaly distributed numbers
3166//
3167// @var *r pointer to a vector where we write the numbers
3168// @var n how many numbers do we generate
3169//
3170// @date Sat Jun 27 05:58:56 MET DST 1998
3171//----------------------------------------------------------------------
3172// @function
3173
3174//!@{
3175void
3176rnormal(double *r, int n)
3177{
3178
3179 double z1, z2;
3180 int i;
3181
3182 for (i=0; i<n; i+=2) {
3183
3184 z1 = RandomNumber;
3185 z2 = RandomNumber;
3186
3187 r[i] = sqrt(-2.0*log(z1)) * cos(2.0*M_PI*z2);
3188 r[i+1] = sqrt(-2.0*log(z1)) * sin(2.0*M_PI*z2);
3189
3190 }
3191
3192}
3193//!@}
3194
3195
3196//!---------------------------------------------------------------------
3197// @name isA
3198//
3199// @desc returns TRUE(FALSE), if the flag is(is not) the given
3200//
3201// @var s1 String to be compared with s2
3202// @var s2 String in the form of a FLAG
3203//
3204// @return TRUE: both strings match; FALSE: otherwise
3205//
3206// @date Wed Jul 8 15:25:39 MET DST 1998
3207//----------------------------------------------------------------------
3208// @function
3209
3210//!@{
3211int
3212isA( char * s1, const char * flag )
3213{
3214 return ( (strncmp((char *)s1, flag, 8)==0) ? TRUE : FALSE );
3215}
3216//!@}
3217
3218
3219//!---------------------------------------------------------------------
3220// @name Curv2Lin
3221//
3222// @desc Curvilinear to Linear (Euclidean) distance
3223//
3224// @var s Curvilinear distance over the parabolic shape
3225//
3226// @return Radial distance from the axis of the paraboloid
3227//
3228// @date Wed Jul 8 15:25:39 MET DST 1998
3229//----------------------------------------------------------------------
3230// @function
3231
3232//!@{
3233float
3234Curv2Lin(float s)
3235{
3236 float x;
3237 short i;
3238
3239 x = s;
3240 for (i = 0; i < 4; i++)
3241 x = (s / 100.) / (1. + (float) 0.000144175317185 * x * x);
3242
3243 return (x*100.);
3244}
3245//!@}
3246
3247
3248//!---------------------------------------------------------------------
3249// @name Lin2Curv
3250//
3251// @desc Linear (Euclidean) to Curvilinear distance
3252//
3253// @var x Radial distance from the axis of the paraboloid
3254//
3255// @return Curvilinear distance over the parabolic shape
3256//
3257// @date Wed Jul 8 15:25:39 MET DST 1998
3258//----------------------------------------------------------------------
3259// @function
3260
3261//!@{
3262float
3263Lin2Curv(float x)
3264{
3265 x /= 100.;
3266 return ((x + (float) 0.000144175317185 * x * x * x)*100.);
3267}
3268//!@}
3269
3270
3271//!---------------------------------------------------------------------
3272// @name get_stdin_files
3273//
3274// @desc get data from stdin and create cerXXXXXX and staXXXXXX
3275//
3276// @var ncerf number of next Cherenkov file cerXXXXXX to be created
3277// @var El Lower limit for requested energy
3278// @var Eu Upper limit for requested energy
3279// @var flag Boolean flag: TRUE: initialize; FALSE: normal
3280//
3281// @return TRUE: succesful; FALSE: failure
3282//
3283// @date Tue Dec 15 15:03:40 MET 1998
3284//----------------------------------------------------------------------
3285// @function
3286
3287//!@{
3288int
3289get_stdin_files(int ncerf, float El, float Eu, int flag)
3290{
3291 static float *buffer; // buffer for read-out of STDIN
3292 static COREventHeader *evth;// Event Header class (from CORSIKA)
3293 static float Elow, Eup; // range for Energy
3294 static int reject;
3295
3296 char *chkbuffer;
3297
3298 char cername[256]; // output filename
3299 ofstream cerfile; // output file (stream)
3300 char staname[256]; // output filename sta
3301 ofstream stafile; // output file sta (stream)
3302 int bytes;
3303 float energy;
3304
3305 /*!@'
3306
3307 The main features of this function, as well as the algorithm,
3308 are taken from |corfilter|.
3309
3310 */
3311
3312 // if it's the first time, allocate buffers
3313
3314 if ( flag == TRUE ) {
3315
3316 if (verbose)
3317 log("get_stdin_files", "Allocating memory for buffer...\n");
3318
3319 buffer = new float [BUFFER_LENGTH];
3320 evth = (COREventHeader*)buffer;
3321
3322 Elow = El;
3323 Eup = Eu;
3324
3325 return ( TRUE );
3326 }
3327
3328 // main loop
3329
3330 while ( ! cin.eof() ) {
3331
3332 // read buffer from STDIN
3333
3334 cin.read( (char*) buffer, BUFFER_LENGTH * sizeof(float) );
3335 bytes = cin.gcount();
3336
3337 if (bytes < (BUFFER_LENGTH * sizeof(float))) {
3338 if (verbose)
3339 logerr("get_stdin_files",
3340 "Can read only %d bytes, instead of required %d.\n",
3341 bytes, BUFFER_LENGTH * sizeof(float));
3342 }
3343
3344 if (bytes == 0)
3345 return ( FALSE );
3346
3347 // check whether EVTH is there
3348
3349 chkbuffer = (char*)buffer;
3350 if (strstr(chkbuffer,EVTH) == 0) {
3351 if (verbose)
3352 logerr("get_stdin_files", "EVTH not found in this block\n");
3353 continue;
3354 } else {
3355 if (verbose)
3356 logerr("get_stdin_files", "EVTH was FOUND in this block!!\n");
3357 }
3358
3359
3360 // show log
3361 if (verbose)
3362 logerr("get_stdin_files", " (E=%8.1f) Reading event %d... ",
3363 (float)evth->Etotal, (int)evth->EvtNumber);
3364
3365 // CRITERIA
3366 energy = evth->Etotal;
3367
3368 if ( (energy<Elow) || (energy>Eup) ) {
3369 if (verbose)
3370 log("REJECTED", "\n");
3371 reject = TRUE;
3372 } else {
3373 if (verbose)
3374 log("ACCEPTED", "\n");
3375 reject = FALSE;
3376 }
3377
3378 if ( !reject ) {
3379
3380 // update no. of cer. file
3381 sprintf(cername, "%s/cer%06d", TMP_STDIN_DIR, ncerf+1);
3382
3383 cerfile.open( cername );
3384
3385 if ( verbose )
3386 log("get_stdin_files", " Writing event %d, will be %s\n",
3387 (int)evth->EvtNumber, cername);
3388
3389 // write it
3390 cerfile.write( (char*)buffer, BUFFER_LENGTH * sizeof(float) );
3391
3392 }
3393
3394 // read blocks till new EVTH is found
3395 // (now must be at the beginning of a block)
3396
3397 if (verbose)
3398 log("get_stdin_files", " Reading photons...\n");
3399
3400 bytes = 0;
3401 while ( TRUE ) {
3402
3403 // read photon
3404 cin.read( (char*)buffer, 7 * sizeof(float) );
3405
3406 // the last photon is reached when the EVTH is read, but
3407 // only when bytes % SIZE_OF_BLOCK == 0
3408 // SIZE_OF_BLOCK = 22932
3409 // with this, we only execute strstr once each 819 photons
3410
3411 // are we in a new block?
3412
3413 if ((bytes % SIZE_OF_BLOCK) == 0) {
3414 // is this the last photon?
3415 if ( strstr(chkbuffer, EVTH) == chkbuffer )
3416 break;
3417 }
3418
3419 // if not, write it to the cer-file
3420 if ( !reject )
3421 cerfile.write( (char*)buffer, 7 * sizeof(float) );
3422
3423 bytes += 7 * sizeof(float);
3424
3425 }
3426
3427 if (verbose)
3428 log("get_stdin_files", " Finished: %d bytes\n", bytes);
3429
3430 if ( !reject ) {
3431 // close cer-file
3432 cerfile.close();
3433 }
3434
3435 // at this point, we must be reading sta-file
3436
3437 // read the rest of the sta-file
3438 if (verbose)
3439 log("get_stdin_files", " Reading sta-data...\n");
3440 cin.read( (char*)buffer + 7*sizeof(float), 5956 - 7*sizeof(float));
3441
3442 if ( !reject ) {
3443
3444 // log
3445 if (verbose)
3446 log("get_stdin_files", " Saving sta-file...\n");
3447
3448 // open file
3449 sprintf(staname, "%s/sta%06d", TMP_STDIN_DIR, ncerf+1);
3450 stafile.open( staname );
3451
3452 // write it
3453 stafile.write( (char*)buffer, 5956 );
3454
3455 // close sta-file
3456 stafile.close();
3457
3458 // get out from here, return to main program
3459 return( TRUE );
3460
3461 }
3462
3463 }
3464
3465 return( TRUE );
3466
3467}
3468//!@}
3469
3470
3471//!---------------------------------------------------------------------
3472// @name get_new_ct_pointing
3473//
3474// @desc returns new random position of CT around a given one
3475//
3476// @var theta Theta (ZA) of the shower
3477// @var phi Phi (AZ) of the shower
3478// @var range Maximum allowed distance between CT and shower axis
3479// @var newtheta Resulting Theta of the CT
3480// @var newphi Resulting Phi of the CT
3481//
3482// @return Angular distance between orig. and new directions
3483//
3484// @date Sat Jun 27 05:58:56 MET DST 1998
3485//----------------------------------------------------------------------
3486// @function
3487
3488//!@{
3489float
3490get_new_ct_pointing(float theta, float phi, float range,
3491 float *newtheta, float *newphi)
3492{
3493 float distance;
3494 float it, ip,Nit;
3495 float sin_theta, cos_theta;
3496 float sin_newtheta, cos_newtheta;
3497 float sin_iphi, cos_iphi;
3498 float iphi;
3499 float range_aux;
3500 int i;
3501
3502 //Random position of the CT. It is the distance
3503 //between shower axis and CT position and its
3504 //distribution follows sin(it).
3505
3506 range_aux=1.0-sin(M_PI/2-range);
3507 Nit = RandomNumber * range_aux;
3508 ip = RandomNumber * 2.0 * M_PI;
3509 it = M_PI/2-asin((1.0-Nit));
3510
3511 sin_theta = sin(theta);
3512 cos_theta = cos(theta);
3513
3514 cos_newtheta = cos_theta*cos(it) - sin_theta*sin(it)*cos(ip);
3515 *newtheta = acos( cos_newtheta );
3516 sin_newtheta = sin( *newtheta );
3517
3518 if(theta!=0){
3519 sin_iphi = sin(it)*sin(ip)/sin_newtheta;
3520 cos_iphi = ((cos(it)-cos_newtheta*cos_theta)/(sin_newtheta*sin_theta));
3521
3522 iphi = atan2( sin_iphi, cos_iphi );
3523 }
3524 else iphi=ip;
3525
3526 *newphi = phi + iphi;
3527
3528 return( it );
3529}
3530//!@}
3531
3532
3533//=------------------------------------------------------------
3534//!@subsection Log of this file.
3535
3536//!@{
3537//
3538// Revision 1.3 2000/01/03 12:41:26 harald
3539// There was a small mistake with +/- signs in the prevois versions. This
3540// line was founded by Jose Carlos. His remark is following:
3541//
3542// Look for the comment I include here in this block: you will see that
3543// now the following expression has a +- sign, where before it was only
3544// + (or -, I don't remember, but only one)
3545//
3546// Revision 1.2 1999/11/01 11:05:54 harald
3547// Small changes to comile the reflector program under linux.
3548// (Different use of NULL on DECalphas-osf1 and on linux)
3549//
3550// Revision 1.1.1.1 1999/10/29 07:00:33 harald
3551// This is the startpoint for the futher development of the Reflector program
3552// of Jose Carlos. For all developments use this CVS-controlled directory.
3553//
3554// Revision 1.18 1999/10/05 11:11:12 gonzalez
3555// Sep.1999
3556//
3557// Revision 1.17 1999/03/24 16:33:01 gonzalez
3558// REFLECTOR 1.1: Release
3559//
3560// Revision 1.16 1999/01/21 16:03:44 gonzalez
3561// Only small modifications
3562//
3563// Revision 1.15 1999/01/19 18:07:16 gonzalez
3564// Bugs in STDIN-STDOUT version corrected.
3565//
3566// Revision 1.14 1999/01/14 17:35:47 gonzalez
3567// Both reading from STDIN (data_from_stdin) and
3568// writing to STDOUT (data_to_STDOUT) working.
3569//
3570// Revision 1.13 1999/01/13 12:41:12 gonzalez
3571// THIS IS A WORKING (last?) VERSION
3572//
3573// Revision 1.12 1998/12/17 16:26:09 gonzalez
3574// *************************************************
3575// *************************************************
3576// WARNING:: Version 1.11 is completly wrong!!
3577// *************************************************
3578// *************************************************
3579//
3580// Revision 1.10 1998/12/15 10:47:22 gonzalez
3581// RELEASE-1.0
3582//
3583// Revision 1.9 1998/11/25 16:30:47 gonzalez
3584// Commit after inclusion of 'Blocking'
3585//
3586//!@}
3587
3588//=EOF
Note: See TracBrowser for help on using the repository browser.