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

Last change on this file was 437, checked in by magicsol, 24 years ago
Few modifications to make the program more robust. They do not afect the algorithm
File size: 96.7 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.11 $
23// $Author: magicsol $
24// $Date: 2000-10-24 14:30:07 $
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//@: number of times a shower is going to be processed
798int nRepeat_Random;
799
800//@: number of times a shower is already processed
801int nRepeated;
802
803//@: maximum random pointing distance
804float Random_Pointing_MaxDist; // [radians]
805
806//!@}
807
808//=---------------------------------------------------------------------
809//!@subsection Main program.
810
811/*!@"
812
813 In the main program we basically follow the items expressed in the
814 first section. Please, have a look into previous sections. I tried
815 to make a code already well commented. Most things should be clear.
816 Nevertheless, I'll try to comment again everything in following
817 releases.
818
819 @"*/
820
821//!@{
822
823//++++++++++++++++++++++++++++++++++++++++
824// MAIN PROGRAM
825//----------------------------------------
826
827int
828main(int argc, char **argv)
829{
830
831 DIR *directory; //@< directory of data
832 struct dirent *de; //@< directory entry
833
834 char pathname[256]; //@< name of directory
835 char cername[256]; //@< name of Cherenkov file
836 char staname[256]; //@< name of Statistics file
837 char outname[256]; //@< name of the output file
838 char parname[256]; //@< name of parameters file
839
840 ifstream cerfile; //@< stream for Cherenkov file
841 ifstream pmfile; //@< stream for pixmap file
842 ofstream outputfile; //@< stream for output file
843 ofstream pbfile; //@< stream for parallel beam file
844
845 COREventHeader evth; //@< Event Header class (from CORSIKA)
846 CORParticle photon; //@< Particle (photon) class (from CORSIKA)
847 CORStatfile stat; //@< Statistics file class (from CORSIKA)
848
849 MCEventHeader mcevth; //@< Event Header class
850 MCCphoton cphoton; //@< Particle (photon) class
851
852 float thetaCT, phiCT, xiCT; //@< parameters of a given shower
853 float coreD, coreX, coreY; //@< core position and distance
854 float u, v, w; //@< cosenos directores
855
856 float r[3]; //@< trayectory
857 float x[3]; //@< position of the photon in the ground
858
859 float rCT[3]; //@< trayectory in the system of the CT
860 float xCT[3]; //@< position of the photon in the ground (CT)
861 float rm[3]; //@< trayectory in the system of a mirror
862 float xmm[3]; //@< intermediate values
863 float xm[3]; //@< position of photon in ground
864 float xcut[3]; //@< location of the cut sphere-trayectory
865 float xcutCT[3]; //@< location of the cut sphere-trayectory (CT)
866
867 float rnor[3], rnorm; //@< normal in that point
868
869 float rrefl[3]; //@< reflected vector, from the mirror
870 float rreflCT[3]; //@< reflected vector, from the CT
871 float xcam[3]; //@< where the photon hits the camera plane
872
873 float calpha; //@< cos(alpha=angle incident/normal)
874 float phi; //@< angle between photon and camera plane
875
876 float a, b, c, t; //@< intermediate variables
877
878 float d; //@< minimum distance trayectory-mirror center
879
880 float wl; //@< wavelength of the Cphoton
881 float reflec; //@< reflectivity for a Cphoton
882
883 int ncer; //@< number of the shower;
884 int np; //@< number of path
885 int i, k; //@< simple counters
886 int i_mirror=-1; //@< number of a given mirror
887
888 int nCphotons; //@< number of a Cphotons written on output file
889
890 float TR; //@< atmospheric transmitance
891 int nbeforeTR, nafterTR; //@< number of Cph. before and after transmission
892 int num_cer_files; //@< number of cer files;
893 int max_num_cer_files; //@< maximum number of cer files to read
894 int f_ev, l_ev; //@< first and last events of a range
895
896 float lE, uE; //@< lower and upper edges of energy range allowed
897
898 float distmirr, distmirr2; //@< distances used in MAGIC reflexion routine
899 float sx, sy;
900 float t1, t2;
901 float pbx, pby, pbt, pbp;
902
903 char line[200];
904
905 int ch, errflg = 0; //@< used by getopt
906
907 /*!@'
908
909 We start in the main program. First we (could) make some
910 presentation, and follows the reading of the parameters file (now
911 from the STDIN), the reading of the CT parameters file, and the
912 creation of the output file, where the processed data will be
913 stored.
914
915 */
916
917 //++
918 // START
919 //--
920
921 // make unbuffered output
922
923 cout.setf ( ios::stdio );
924
925 // parse command line options (see reflector.h)
926
927 parname[0] = '\0';
928
929 optarg = 0;
930 while ( !errflg && ((ch = getopt(argc, argv, COMMAND_LINE_OPTIONS)) != -1) )
931 switch (ch) {
932 case 'f':
933 strcpy(parname, optarg);
934 break;
935 case 'h':
936 usage();
937 break;
938 default :
939 errflg++;
940 }
941
942 // show help if error
943
944 if ( errflg>0 )
945 usage();
946
947 // make some sort of presentation
948 // present();
949
950 // read parameters from the stdin
951
952 if ( strlen(parname) < 1 )
953 readparam(0);
954 else
955 readparam(parname);
956
957 // blocking-mode?
958
959 Block = get_block();
960
961 // reading data from STDIN?
962
963 Data_From_STDIN = get_data_from_stdin();
964
965 // writing data to STDOUT?
966
967 Data_To_STDOUT = get_data_to_stdout();
968
969 // make it verbosely?
970
971 verbose = get_verbose();
972
973 /*!@'
974
975 @#### Parallel beam of light.
976
977 In order to test the behaviour of the program, we can feed it with
978 an artificial data file. This data file will be generated by the
979 program itself and it will consist in a grid of photons, forming
980 all together a parallel beam of light with certain
981 characteristics. These characteristics will be given by the user
982 in the |parameters file|, in the following form.
983
984 parallel\_beam |<theta> <phi> <Xcen> <Ycen> <lenX> <lenY> <nX> <nY> <H>|
985
986 or
987
988 parallel\_beam\_pm |<filename> <scale> <H>|
989
990 In the first form, the first two values give the arrival direction
991 of the beam (zenith angle and azimuth) in degrees, and the last
992 six values specify the geometry of the beam, in cm, and the height
993 of production, also in cm. In the second form, we use a pixmap
994 file in PPM ASCII format, and scale it to create (for 0 degrees of
995 zenith angle) a parallel beam of light.
996
997 */
998
999 //++
1000 // PARALLEL BEAM OF LIGHT
1001 //--
1002
1003 // are we going to use a parallel beam of light?
1004
1005 pb_ParallelBeam = is_parallel_beam();
1006 pb_ParallelBeamPM = is_parallel_beam_pm();
1007
1008 if ( pb_ParallelBeam == TRUE ) {
1009
1010 if ( pb_ParallelBeamPM == FALSE ) {
1011
1012 // get parameters of the parallel beam of light
1013
1014 get_parallel_beam(&pb_Theta, &pb_Phi,
1015 &pb_X, &pb_Y,
1016 &pb_LengthX, &pb_LengthY,
1017 &pb_NX, &pb_NY, &pb_Height);
1018
1019 // generate artificial file(s)
1020
1021 // a. cerfile
1022
1023 sprintf( cername, "%s/cer999999", get_path_name(0) );
1024
1025 pbfile.open( cername, ios::out );
1026
1027 evth.fill ( 1., 1., 1000., 1., 1.e6, 0., 0., 0., pb_Theta, pb_Phi,
1028 pb_X, pb_Y);
1029
1030 evth.write( pbfile );
1031
1032 for ( pbx = (pb_X - pb_LengthX/2.0);
1033 pbx <= (pb_X + pb_LengthX/2.0);
1034 pbx += (pb_LengthX/pb_NX) )
1035
1036 for ( pby = (pb_Y - pb_LengthY/2.0);
1037 pby <= (pb_Y + pb_LengthY/2.0);
1038 pby += (pb_LengthY/pb_NY) ) {
1039
1040 pbt = -atan2(sqrt(SQR(pbx)+SQR(pby)),pb_Height);
1041 pbp = atan2(pby, pbx);
1042 photon.fill(450.0,
1043 pbx, pby,
1044 -sin(pbt) * cos(pbp),
1045 -sin(pbt) * sin(pbp),
1046 0., pb_Height);
1047
1048 photon.write( pbfile );
1049
1050 }
1051
1052 photon.fill( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 );
1053
1054 photon.write( pbfile );
1055
1056 pbfile.close();
1057
1058 // b. statfile
1059
1060 sprintf( staname, "%s/sta999999", get_path_name(0) );
1061
1062 pbfile.open( staname );
1063
1064 stat.write( pbfile );
1065
1066 pbfile.close();
1067
1068 } else {
1069
1070 // get parameters of the parallel beam of light from PIXMAP
1071
1072 strcpy(pb_Filename, get_parallel_beam_pm(&pb_Scale,&pb_Height));
1073
1074 pmfile.open( pb_Filename, ios::out );
1075 pmfile.getline( line, 300 );
1076
1077 do {
1078 pmfile.getline( line, 300 );
1079 } while (line[0] == '#');
1080
1081 sscanf( line, "%f %f", &pb_NX, &pb_NY );
1082 pb_LengthX = pb_NX * pb_Scale;
1083 pb_LengthY = pb_NY * pb_Scale;
1084
1085 // generate artificial file(s)
1086
1087 // a. cerfile
1088
1089 sprintf( cername, "%s/cer999999", get_path_name(0) );
1090
1091 pbfile.open( cername, ios::out );
1092
1093 evth.fill ( 1., 1., 1000., 1., 1.e6, 0., 0., 0., 0., 0., 0., 0.);
1094
1095 evth.write( pbfile );
1096
1097 for ( pbx = (pb_X - pb_LengthX/2.0);
1098 pbx <= (pb_X + pb_LengthX/2.0);
1099 pbx += (pb_LengthX/pb_NX) )
1100
1101 for ( pby = (pb_Y - pb_LengthY/2.0);
1102 pby <= (pb_Y + pb_LengthY/2.0);
1103 pby += (pb_LengthY/pb_NY) ) {
1104
1105 pmfile >> i;
1106
1107 if ( i < 1 )
1108 continue;
1109
1110 pbt = -atan2(sqrt(SQR(pbx)+SQR(pby)),pb_Height);
1111 pbp = atan2(pby, pbx);
1112 photon.fill(450.0,
1113 pbx, pby,
1114 -sin(pbt) * cos(pbp),
1115 -sin(pbt) * sin(pbp),
1116 0., pb_Height);
1117
1118 photon.write( pbfile );
1119
1120 }
1121
1122 photon.fill( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 );
1123
1124 photon.write( pbfile );
1125
1126 pbfile.close();
1127
1128 // b. statfile
1129
1130 sprintf( staname, "%s/sta999999", get_path_name(0) );
1131
1132 pbfile.open( staname );
1133
1134 stat.write( pbfile );
1135
1136 pbfile.close();
1137
1138 pmfile.close();
1139
1140 }
1141
1142 } // YES parallel beam
1143
1144 // read parameters from the stdin
1145
1146 read_ct_file();
1147
1148 // set all random numbers seeds
1149
1150 setall( get_seeds(0), get_seeds(1) );
1151
1152 // get name of the output file
1153
1154 strcpy( outname, get_output_filename() );
1155
1156 // open((re)create) output file
1157
1158 if ( verbose >= VERBOSE_MINIMAL )
1159 log( SIGNATURE, "Openning/creating output file %s\n", outname );
1160
1161 outputfile.open( outname );
1162
1163 if ( outputfile.bad() )
1164 error( SIGNATURE, "Cannot open output file: %s\n", outname );
1165
1166 // save signature to the output file
1167
1168 if ( Data_To_STDOUT )
1169 cout.write( SIGNATURE, sizeof(SIGNATURE) );
1170 else
1171 outputfile.write( SIGNATURE, sizeof(SIGNATURE) );
1172
1173 // get random pointing variables
1174
1175 Random_Pointing = get_random_pointing( &Random_Pointing_MaxDist );
1176
1177 // generate a sort of log information
1178
1179 if ( verbose >= VERBOSE_MINIMAL ) {
1180 log( SIGNATURE, "Random poi.: %f\n", Random_Pointing_MaxDist );
1181 log( SIGNATURE, "Atmospheric model: %s\n", get_atm_mod() );
1182 log( SIGNATURE, "Number of paths: %d\n", get_num_of_paths() );
1183 log( SIGNATURE, "Output filename: %s\n", outname );
1184 log( SIGNATURE, "CT def. filename: %s\n", get_ct_filename() );
1185 }
1186
1187 // get energy range
1188
1189 get_energy_cuts( &lE, &uE );
1190
1191 // get number of times which a shower is going to be used
1192 // at different pointing directions
1193 nRepeat_Random = get_repeat_random();
1194
1195 /*!@'
1196
1197 After reading the parameters file (from STDIN), and the CT
1198 definition file, we start the analysis. Basically, there are three
1199 loops. The outermost loop is a loop in the E@data directories@ we
1200 introduced in the parameters file; the middle loop follows each
1201 |cer*| file in the directory; and the innermost loop gets each
1202 Cherenkov from these files.
1203
1204 */
1205
1206 // for each path of data files
1207
1208 for (np=0; np<get_num_of_paths(); np++) {
1209
1210 strcpy(pathname, get_path_name(np));
1211
1212 // open directory
1213
1214 if ( verbose >= VERBOSE_MINIMAL )
1215 log( SIGNATURE, "Openning directory %s\n", pathname );
1216
1217 directory = opendir(pathname);
1218
1219 if ( directory == 0 )
1220 error( SIGNATURE,
1221 "Cannot open directory %s\n", pathname );
1222
1223 // write start-of-run mark
1224
1225 if ( Data_To_STDOUT )
1226 cout.write( FLAG_START_OF_RUN, SIZE_OF_FLAGS );
1227 else
1228 outputfile.write( FLAG_START_OF_RUN, SIZE_OF_FLAGS );
1229
1230 // trace the number of events (cer* files)
1231
1232 num_cer_files = 0;
1233
1234 // get maximum number, first and last event
1235
1236 max_num_cer_files = get_max_events();
1237
1238 get_range_events(&f_ev, &l_ev);
1239
1240 // if reading from disk, save position of handler in directory
1241
1242 posdir = telldir( directory );
1243
1244 if ( Data_From_STDIN ) {
1245
1246 if ( verbose >= VERBOSE_MINIMAL )
1247 log(SIGNATURE, " * * * READING DATA FROM STDIN * * *\n");
1248
1249 get_stdin_files(0, lE, uE, TRUE);
1250
1251 }
1252
1253 //--- LOOP ON FILES ---
1254 // for each directory entry (files),
1255 // and(or) while still data has to be read
1256 // this condition is:
1257 // while ( ( there are files AND
1258 // maximum number is not reached )
1259 // OR
1260 // ( reading from stdin AND
1261 // maximum number is not reached ) )
1262
1263 nRepeated = 0;
1264
1265 while ( (
1266 ((de = readdir( directory )) != 0) &&
1267 (num_cer_files < max_num_cer_files)
1268 )
1269 ||
1270 (
1271 ( Data_From_STDIN ) &&
1272 (num_cer_files < max_num_cer_files)
1273 ) ) {
1274 /*
1275 // increment the number of times this files is used
1276 nRepeated++;
1277 */
1278
1279 // if Block > 0, then we wait till a file __DOIT is present
1280
1281 if ( Block > 0 ) {
1282
1283 if ( ! WORKING ) {
1284
1285 while ( ! WORKING ) {
1286 // sleep 10 minutes
1287 sleep( 60 );
1288
1289 // look for the file __DOIT
1290 obfile.open("__DOIT");
1291 if ( ! obfile.bad() ) {
1292 WORKING = TRUE;
1293 obfile.close();
1294 }
1295 }
1296
1297 obfile.open("__WORKING");
1298 remove("__DOIT");
1299
1300 rewinddir( directory );
1301 seekdir( directory, posdir );
1302
1303 nfilesBlock = 0;
1304
1305 continue;
1306
1307 } else {
1308
1309 if ( nfilesBlock == Block ) {
1310
1311 posdir = telldir( directory );
1312 obfile.close();
1313 obfile.open("__READY");
1314 obfile << nfilesBlock << " files/block, "
1315 << ncer << " showers.\n";
1316 obfile.close();
1317 remove("__WORKING");
1318 WORKING = FALSE;
1319
1320 continue;
1321 }
1322
1323 }
1324
1325 } // if Block
1326
1327 if ( Data_From_STDIN ) {
1328
1329 if ( verbose >= VERBOSE_MINIMAL )
1330 log(SIGNATURE, "Try to get data from STDIN\n");
1331
1332 if ( ! get_stdin_files(num_cer_files) ) {
1333 if ( verbose >= VERBOSE_MINIMAL )
1334 log(SIGNATURE, "get_stdin_files returned FALSE.\n");
1335 if (cin.eof()) {
1336 num_cer_files = max_num_cer_files + 1;
1337 continue;
1338 }
1339 }
1340
1341 ++num_cer_files;
1342 ncer = num_cer_files;
1343
1344 } else {
1345
1346 // skip removed files
1347
1348 if ( de->d_ino == 0)
1349 continue;
1350
1351 // keep only cer* files
1352
1353 if ( strstr(de->d_name, "cer") != de->d_name )
1354 continue;
1355
1356 if ( Block > 0 ) {
1357 ++nfilesBlock;
1358 obfile << nfilesBlock;
1359 }
1360
1361 // get cer* file number (number of the shower)
1362
1363 // increments the pointer in 3 to calculate the number
1364
1365 ncer = atoi(de->d_name + 3);
1366
1367 if ( (ncer < f_ev) || (ncer > l_ev) )
1368 continue;
1369
1370 // it is a real cer* file, and we want to read it
1371
1372 ++num_cer_files;
1373
1374 if ( verbose >= VERBOSE_NORMAL )
1375 log( SIGNATURE, "cerfile: %s\n", de->d_name );
1376
1377 } // data from disk
1378
1379 // Here we initialize the counters of Cherenkov photons at
1380 // different levels of the analysis.
1381
1382 nCphotons = nbeforeTR = nafterTR = 0;
1383
1384 // full names of the input files cer* and sta*
1385
1386 sprintf(cername, "%s/cer%06d", pathname, ncer);
1387 sprintf(staname, "%s/sta%06d", pathname, ncer);
1388
1389 // try to open cer* file
1390
1391 if ( verbose >= VERBOSE_MAXIMAL )
1392 log( SIGNATURE, "Opening %s\n", cername );
1393
1394 cerfile.open( cername );
1395
1396 if ( cerfile.bad() )
1397 error( SIGNATURE,
1398 "Cannot open input file: %s\n", cername );
1399
1400 /*!@'
1401
1402 Each shower has associated three files: E@ Cherenkov-photons
1403 file @ (|cer*|), E@ Particles file @ (|dat*|) and E@
1404 Statistics file @ (|sta*|). First we read the data from this
1405 last file, which is explained below, and then we read the E@
1406 Event-Header @ block from the |cer*| file.
1407
1408 */
1409
1410 // try to open the statistics file
1411
1412 if ( verbose >= VERBOSE_MAXIMAL )
1413 log( SIGNATURE, "Opening %s\n", staname );
1414
1415 stat.openfile ( staname );
1416
1417 // reading data from this sta* file
1418
1419 if ( verbose >= VERBOSE_MAXIMAL )
1420 log( SIGNATURE, "Reading %s\n", staname );
1421
1422 stat.read();
1423
1424 // read the Event Header from the cer* file
1425
1426 evth.read( cerfile );
1427
1428 if ( verbose >= VERBOSE_MAXIMAL ) {
1429 log( SIGNATURE, "EVENT HEADER:\n");
1430 evth.print();
1431 }
1432
1433 // check energy
1434
1435 if ( (evth.get_energy() < lE) ||
1436 (evth.get_energy() > uE) ) {
1437
1438 // close files
1439 cerfile.close();
1440 stat.closefile();
1441
1442 // skip it
1443 if ( verbose >= VERBOSE_MAXIMAL )
1444 log( SIGNATURE, "Bad energy.\n");
1445
1446 continue;
1447 }
1448
1449 // write start-of-event mark
1450
1451 if ( Data_To_STDOUT )
1452 cout.write( FLAG_START_OF_EVENT, SIZE_OF_FLAGS );
1453 else
1454 outputfile.write( FLAG_START_OF_EVENT, SIZE_OF_FLAGS );
1455
1456 // save information from the cer* to the mc-evth
1457
1458 mcevth.transport( &evth );
1459
1460 // save time interval in the mc-evth
1461
1462 mcevth.put_times ( stat.get_tfirst(), stat.get_tlast() );
1463
1464 // get direction where the CT is pointing to
1465 // (or, better, from when the shower is coming from)
1466
1467 // first, put the values from the shower
1468 thetaCT = evth.get_theta();
1469 phiCT = evth.get_phi();
1470
1471 // then, may be we have to change them
1472
1473 // do we want random pointing (around shower axis) ?
1474 if ( Random_Pointing == TRUE ) {
1475
1476 // we do, then get a random position
1477 xiCT = get_new_ct_pointing(thetaCT, phiCT,
1478 Random_Pointing_MaxDist,
1479 &thetaCT, &phiCT);
1480
1481 xiCT = xiCT; // thi is to avoid "was set but not used" messages
1482
1483 // save deviations in mc-evth
1484 mcevth.put_deviations (thetaCT - evth.get_theta(),
1485 phiCT - evth.get_phi());
1486
1487 if ( verbose >= VERBOSE_MAXIMAL )
1488 log( SIGNATURE, "NEW POINTING: (%f,%f) -> (%f,%f)\n",
1489 evth.get_theta(), evth.get_phi(), thetaCT, phiCT);
1490
1491 } else {
1492
1493 mcevth.put_deviations (0.0, 0.0);
1494
1495 if ( get_fixed_target( &thetaCT, &phiCT ) == FALSE ) {
1496 thetaCT = evth.get_theta();
1497 phiCT = evth.get_phi();
1498 }
1499
1500 }
1501
1502 // save event-header to the output file
1503
1504 if ( Data_To_STDOUT )
1505 cout.write ( (char *)&mcevth, mcevth.mysize() );
1506 else
1507 mcevth.write( outputfile );
1508
1509 if ( verbose >= VERBOSE_MAXIMAL )
1510 mcevth.print();
1511
1512 // get core distance and position
1513
1514 coreD = evth.get_core(&coreX, &coreY);
1515 coreD=coreD; // to avoid "was set but not used" messages
1516
1517 // calculate matrices for changing of coordinate system from-to CT
1518
1519 makeOmega( thetaCT, phiCT );
1520 makeOmegaI( thetaCT, phiCT );
1521
1522 memcpy( OmegaCT, Omega, 9*sizeof(float) );
1523 memcpy( OmegaICT, OmegaI, 9*sizeof(float) );
1524
1525 // read each and every single photon (particle) from cer* file
1526
1527 while ( ! (cerfile.eof() || cerfile.bad()) ) {
1528
1529 // read photon data from the file
1530
1531 photon.read( cerfile );
1532
1533 wl = photon.get_wl();
1534
1535 if ( wl < 1.0 )
1536 break;
1537
1538 // get director cosines in z-axis
1539
1540 w = r[2] = photon.get_w();
1541
1542 /*!@'
1543
1544 @#### Atmospheric attenuation.
1545
1546 The routine |atm()| calls the proper attenuation routine,
1547 and returns the transmission coefficient for that
1548 whavelength, height and zenith angle.
1549
1550 */
1551
1552 //++
1553 // FILTER No. 1: ATMOSPHERE
1554 //--
1555
1556 ++nbeforeTR;
1557
1558 // get transmitance of the atmosphere
1559
1560 TR = atm( photon.get_wl(), photon.get_h(), acos(w) );
1561
1562 // passes the atmosphere?
1563
1564 // if random > TR then goes to the TOP of the loop again
1565
1566 if ( RandomNumber > TR )
1567 continue;
1568
1569 ++nafterTR;
1570
1571 /*!@'
1572
1573 @#### Reflectivity of the mirrors.
1574
1575 We make a 3rd. order interpolation using Lagrange
1576 polynomials, in order to calculate the reflectivity of the
1577 mirror for that wavelength. Further developments will
1578 include also a incidence-angle dependence (which is not very
1579 important).
1580
1581 */
1582
1583 //++
1584 // FILTER No. 2: REFLECTIVITY R(lambda)
1585 //--
1586
1587 // If photon has a larger wl than the largest in Reflectivity
1588 // array then it is skiped.
1589
1590 if (wl>Reflectivity[0][nReflectivity-1]) continue;
1591
1592 // find data point to be used in Lagrange interpolation (-> k)
1593
1594 FindLagrange(Reflectivity,k,wl);
1595
1596 // if random > reflectivity then goes to the TOP of the loop again
1597
1598 reflec = Lagrange(Reflectivity,k,wl);
1599
1600 if ( RandomNumber > reflec )
1601 continue;
1602
1603 /*!@'
1604
1605 @#### Reflexion in the mirrors.
1606
1607 The reflexion is made now is a somehow dirty way: if the
1608 value of the variable |ct\_Type| is 0, the telescope is CT1,
1609 and if the value is 1, the telescope is MAGIC. Due to the
1610 completly different geometries of the dishes and the
1611 mirrors, now we support two fully separated and different
1612 algorithms.
1613
1614 */
1615
1616 //++
1617 // REFLEXION
1618 //--
1619
1620 // get director cosines x,y and save original trayectory
1621
1622 u = r[0] = photon.get_u();
1623 v = r[1] = photon.get_v();
1624 u=u, v=v; // to avoid "was set but not used" messages
1625
1626 // get position in the ground
1627
1628 x[0] = photon.get_x() - coreX;
1629 x[1] = photon.get_y() - coreY;
1630 x[2] = 0.0;
1631
1632 if ( verbose >= VERBOSE_MAXIMAL ) {
1633 cout << '@' << '1'
1634 << ' ' << ncer
1635 << ' ' << t
1636 << ' ' << x[0]
1637 << ' ' << x[1]
1638 << ' ' << x[2]
1639 << ' ' << r[0]
1640 << ' ' << r[1]
1641 << ' ' << r[2]
1642 << ' ' << thetaCT
1643 << ' ' << phiCT
1644 << endl;
1645 }
1646
1647 // change to the system of the CT
1648
1649 applyMxV( OmegaCT, x, xCT );
1650 applyMxV( OmegaCT, r, rCT );
1651
1652 // before moving to the system of the mirror, for MAGIC,
1653 // first we look whether the photon hits a mirror or not
1654
1655 // calculate the intersection of the trayectory of the photon
1656 // with the GLOBAL DISH !!!
1657 // we reproduce the calculation of the coefficients of the
1658 // second order polynomial in z (=xCT[2]), made with
1659 // Mathematica
1660
1661 /*
1662 * In[1]:= parab:=z-(x^2+y^2)/(4F)
1663 * par1=parab /. {x->x0+u/w(z-z0),y->y0+v/w(z-z0)}
1664 *
1665 * Out[1]=
1666 * u (z - z0) 2 v (z - z0) 2
1667 * (x0 + ----------) + (y0 + ----------)
1668 * w w
1669 * z - ---------------------------------------
1670 * 4 F
1671 *
1672 * In[2]:= CoefficientList[ExpandAll[par1*4F*w^2],z]
1673 *
1674 * Out[2]=
1675 * 2 2 2 2
1676 * {-(w x0 ) - w y0 + 2 u w x0 z0 +
1677 *
1678 * 2 2 2 2
1679 * 2 v w y0 z0 - u z0 - v z0 ,
1680 *
1681 * 2 2
1682 * 4 F w - 2 u w x0 - 2 v w y0 + 2 u z0 +
1683 *
1684 * 2 2 2
1685 * 2 v z0, -u - v }
1686 */
1687
1688 // the z coordinate is calculated
1689
1690 a = - SQR(rCT[0]) - SQR(rCT[1]);
1691 b = 4.0*ct_Focal_mean*SQR(rCT[2])
1692 - 2.0*rCT[0]*rCT[2]*xCT[0] - 2.0*rCT[1]*rCT[2]*xCT[1]
1693 + 2.0*SQR(rCT[0])*xCT[2] + 2.0*SQR(rCT[1])*xCT[2];
1694 c = 2*rCT[0]*rCT[2]*x[0]*x[2] + 2*rCT[1]*rCT[2]*x[1]*x[2]
1695 - SQR(rCT[2])*SQR(x[0]) - SQR(rCT[2])*SQR(x[1])
1696 - SQR(rCT[0])*SQR(x[2]) - SQR(rCT[1])*SQR(x[2]);
1697
1698 if ( fabs(a) < 1.e-6 ) {
1699
1700 // only one value
1701
1702 xcut[2] = -c / b;
1703
1704 } else {
1705
1706 d = sqrt( b*b - 4.0*a*c );
1707
1708 // two possible values for z
1709
1710 t1 = (-b+d) / (2.0*a);
1711 t2 = (-b-d) / (2.0*a);
1712
1713 // z must be the minimum of t1 and t2
1714
1715 xcut[2] = (t1 < t2) ? t1 : t2;
1716
1717 }
1718
1719 // xcut[] is NOW the cut between the GLOBAL dish of MAGIC and
1720 // the trayectory of the photon
1721
1722 xcut[0] = xCT[0] + rCT[0]/rCT[2]*(xcut[2]-xCT[2]);
1723 xcut[1] = xCT[1] + rCT[1]/rCT[2]*(xcut[2]-xCT[2]);
1724
1725 if ( ct_Type == 1 ) {
1726
1727 // ++ >>>>>>>>>> MAGIC <<<<<<<<<<
1728
1729 // convert to Curvilinear distance over the parabolic dish
1730
1731 sx = Lin2Curv( xcut[0] );
1732 sy = Lin2Curv( xcut[1] );
1733
1734 // is it outside the dish?
1735
1736 if ((fabs(sx) > 850.0) || (fabs(sy) > 850.0)) {
1737 //cout << "CONDICION 1 !" << endl << flush;
1738 //cout << '1';
1739 continue;
1740 }
1741
1742 // -- >>>>>>>>>> MAGIC <<<<<<<<<<
1743
1744 } // endif (ct_Type == 1)
1745
1746 // calculate the mirror to be used
1747
1748 distmirr = 1000000.;
1749
1750 for (i=0; i<ct_NMirrors && distmirr>=ct_RMirror; ++i) {
1751 distmirr2 = sqrt(SQR(ct_data[i][CT_X] - xcut[0]) +
1752 SQR(ct_data[i][CT_Y] - xcut[1]) +
1753 SQR(ct_data[i][CT_Z] - xcut[2]));
1754 if (distmirr2 < distmirr) {
1755 i_mirror = i;
1756 distmirr = distmirr2;
1757 }
1758 }
1759
1760 // the mirror to use is i_mirror (calculated several lines above)
1761 // check whether the photon is outside the nearest (this) mirror
1762
1763 if ( ct_Type == 0 ) {
1764
1765 // ++ >>>>>>>>>> CT1 <<<<<<<<<<
1766 if (distmirr > ct_RMirror) {
1767 //cout << "CONDICION 2 !" << endl << flush;
1768 //cout << '2';
1769 continue;
1770 }
1771 // -- >>>>>>>>>> CT1 <<<<<<<<<<
1772
1773 } else {
1774
1775 // ++ >>>>>>>>>> MAGIC <<<<<<<<<<
1776 if ((fabs(ct_data[i_mirror][CT_SX] - sx) > ct_RMirror) ||
1777 (fabs(ct_data[i_mirror][CT_SY] - sy) > ct_RMirror)) {
1778 //cout << "CONDICION 2 !" << endl << flush;
1779 //cout << '2';
1780 continue;
1781 }
1782 // -- >>>>>>>>>> MAGIC <<<<<<<<<<
1783
1784 }
1785
1786 // calculate matrices for the mirror
1787
1788 makeOmega (-ct_data[i_mirror][CT_THETAN],
1789 ct_data[i_mirror][CT_PHIN]);
1790 makeOmegaI(-ct_data[i_mirror][CT_THETAN],
1791 ct_data[i_mirror][CT_PHIN]);
1792
1793 // change to the system of the mirror
1794
1795 // first translation...
1796
1797 xmm[0] = xCT[0] - ct_data[i_mirror][CT_X];
1798 xmm[1] = xCT[1] - ct_data[i_mirror][CT_Y];
1799 xmm[2] = xCT[2] - ct_data[i_mirror][CT_Z];
1800
1801 // ...then rotation
1802
1803 applyMxV( Omega, xmm, xm );
1804 applyMxV( Omega, rCT, rm );
1805
1806 // the vector rCT should be normalized, and
1807 // so the vector rm remains normalized as well, but, anyhow...
1808
1809 rnorm = NORM( rm );
1810 rm[0] /= rnorm;
1811 rm[1] /= rnorm;
1812 rm[2] /= rnorm;
1813
1814 // calculate the intersection of the trayectory of the photon
1815 // with the mirror
1816 // we reproduce the calculation of the coefficients of the
1817 // second order polynomial in z (=xm[2]), made with
1818 // Mathematica
1819
1820 /*
1821 * In[1]:= esfera:=x^2+y^2+(z-R)^2-R^2;
1822 * recta:={x->x0+u/w(z-z0),y->y0+v/w(z-z0)}
1823 *
1824 * In[2]:= esfera
1825 *
1826 * 2 2 2 2
1827 * Out[2]= -R + x + y + (-R + z)
1828 *
1829 * In[3]:= recta
1830 *
1831 * u (z - z0) v (z - z0)
1832 * Out[3]= {x -> x0 + ----------, y -> y0 + ----------}
1833 * w w
1834 *
1835 * In[4]:= esf=esfera /. recta
1836 *
1837 * 2 2 u (z - z0) 2 v (z - z0) 2
1838 * Out[4]= -R + (-R + z) + (x0 + ----------) + (y0 + ----------)
1839 * w w
1840 *
1841 * In[5]:= coefs=CoefficientList[ExpandAll[esf],z]
1842 *
1843 * 2 2 2 2
1844 * 2 2 2 u x0 z0 2 v y0 z0 u z0 v z0
1845 * Out[5]= {x0 + y0 - --------- - --------- + ------ + ------,
1846 * w w 2 2
1847 * w w
1848 *
1849 * 2 2 2 2
1850 * 2 u x0 2 v y0 2 u z0 2 v z0 u v
1851 * > -2 R + ------ + ------ - ------- - -------, 1 + -- + --}
1852 * w w 2 2 2 2
1853 * w w w w
1854 * In[6]:= Simplify[ExpandAll[coefs*w^2]]
1855 *
1856 * 2 2 2 2 2 2
1857 * Out[6]= {w (x0 + y0 ) - 2 w (u x0 + v y0) z0 + (u + v ) z0 ,
1858 *
1859 * 2 2 2 2 2
1860 * > -2 (R w - u w x0 + u z0 + v (-(w y0) + v z0)), u + v + w }
1861 *
1862 */
1863
1864 // the z coordinate is calculated, using the coefficients
1865 // shown above
1866
1867 a = SQR(rm[0]) + SQR(rm[1]) + SQR(rm[2]);
1868 if ( ct_Type == 0 ) {
1869 // CT1
1870 b = -2*(2.*ct_Focal_mean*SQR(rm[2])
1871 - rm[0]*rm[2]*xm[0]
1872 + SQR(rm[0])*xm[2]
1873 + rm[1]*(-(rm[2]*xm[1]) + rm[1]*xm[2]));
1874 } else {
1875 // MAGIC
1876 b = -2*(2.*ct_data[i_mirror][CT_FOCAL]*SQR(rm[2])
1877 - rm[0]*rm[2]*xm[0]
1878 + SQR(rm[0])*xm[2]
1879 + rm[1]*(-(rm[2]*xm[1]) + rm[1]*xm[2]));
1880 }
1881 c = (SQR(rm[2])*(SQR(xm[0]) + SQR(xm[1]))
1882 - 2*rm[2]*(rm[0]*xm[0] + rm[1]*xm[1])*xm[2] +
1883 (SQR(rm[0]) + SQR(rm[1]))*SQR(xm[2]));
1884
1885 d = sqrt( b*b - 4.0*a*c );
1886
1887 // two possible values for z
1888
1889 t1 = (-b+d) / (2.0*a);
1890 t2 = (-b-d) / (2.0*a);
1891
1892 // z must be the minimum of t1 and t2
1893
1894 xcut[2] = (t1 < t2) ? t1 : t2;
1895 xcut[0] = xm[0] + rm[0]/rm[2]*(xcut[2]-xm[2]);
1896 xcut[1] = xm[1] + rm[1]/rm[2]*(xcut[2]-xm[2]);
1897
1898 //++
1899 // BLACK SPOTS: If the photon hits the black spot, it's lost
1900 //--
1901
1902 if ( sqrt(SQR(xcut[0]) + SQR(xcut[1])) < ct_BlackSpot_rad ) {
1903 //cout << "CONDITION 3!\n" << flush;
1904 //cout << '3';
1905 continue;
1906 }
1907
1908 // if still we have the photon, continue with the reflexion
1909
1910 // calculate normal vector in this point
1911 // (and normalize, with the sign changed)
1912
1913 rnor[0] = 2.0*xcut[0];
1914 rnor[1] = 2.0*xcut[1];
1915 rnor[2] = 2.0*(xcut[2] - 2.0*ct_Focal[i_mirror]);
1916
1917 rnorm = -NORM( rnor );
1918 rnor[0] /= rnorm;
1919 rnor[1] /= rnorm;
1920 rnor[2] /= rnorm;
1921
1922 // now, both "normal" vector and original trayectory are
1923 // normalized
1924 // just project the original vector in the normal, and
1925 // take it as the "mean" position of the original and
1926 // the "reflected" vector
1927 // from this, we can calculate the "reflected" vector
1928 // calpha = cos(angle(rnor,rm))
1929
1930 calpha = fabs(rnor[0]*rm[0] + rnor[1]*rm[1] + rnor[2]*rm[2]);
1931
1932 // finally!!! we have the reflected trayectory of the photon
1933
1934 rrefl[0] = 2.0*rnor[0]*calpha - rm[0];
1935 rrefl[1] = 2.0*rnor[1]*calpha - rm[1];
1936 rrefl[2] = 2.0*rnor[2]*calpha - rm[2];
1937
1938 rnorm = NORM( rrefl );
1939 rrefl[0] /= rnorm;
1940 rrefl[1] /= rnorm;
1941 rrefl[2] /= rnorm;
1942
1943 // let's go back to the coordinate system of the CT
1944
1945 // first rotation...
1946 applyMxV( OmegaI, xcut, xcutCT);
1947 applyMxV( OmegaI, rrefl, rreflCT);
1948
1949 // ...then translation
1950 xcutCT[0] += ct_data[i_mirror][CT_X];
1951 xcutCT[1] += ct_data[i_mirror][CT_Y];
1952 xcutCT[2] += ct_data[i_mirror][CT_Z];
1953
1954 // calculate intersection of this trayectory and the camera plane
1955 // in the system of the CT, this plane is z = ct_Focal
1956
1957 t = (ct_Focal_mean - xcutCT[2]) / rreflCT[2];
1958
1959 xcam[0] = xcutCT[0] + rreflCT[0]*t;
1960 xcam[1] = xcutCT[1] + rreflCT[1]*t;
1961 xcam[2] = xcutCT[2] + rreflCT[2]*t;
1962
1963 //++
1964 // AXIS DEVIATION: We introduce it here just as a first order
1965 // correction, by modifying the position of the reflected photon.
1966 //--
1967
1968 xcam[0] += AxisDeviation[0][i_mirror];
1969 xcam[1] += AxisDeviation[1][i_mirror];
1970
1971 //++
1972 // SMEARING: We apply the point spread function for the mirrors
1973 //--
1974
1975 // get two N(0;1) random numbers
1976
1977 rnormal( NormalRandomNumbers, 2 );
1978
1979 // modify the Cphoton position in the camera
1980
1981 xcam[0] += NormalRandomNumbers[0] * ct_PSpread_mean;
1982 xcam[1] += NormalRandomNumbers[1] * ct_PSpread_mean;
1983
1984 // check whether the photon goes out of the camera
1985
1986 if ( (SQR(xcam[0])+SQR(xcam[1])) > SQR(ct_CameraWidth) ) {
1987 continue;
1988 }
1989
1990 //++
1991 // ANGLE OF INCIDENCE
1992 //--
1993
1994 // calculate angle of incidence between tray. and camera plane
1995 // the camera plane is
1996 // 0 y + 0 y + z - ct_Focal = 0 => (A,B,C,D) = (0,0,1,-ct_Focal)
1997 // from Table 3.20 "Tasch. der Math."
1998
1999 phi = asin(rreflCT[2]);
2000
2001 //++
2002 // TIMING
2003 //--
2004
2005 // calculate the new time of the photon (in the camera)
2006
2007 // initial time since first interaction
2008
2009 t = photon.get_t();
2010
2011 // substract path from the mirror till the ground, 'cos
2012 // the photon actually hit the mirror!!
2013
2014
2015 t = t + ((( xm[2] > 0. ) ? -1.0 : +1.0) *
2016 sqrt( SQR(xm[0] - xcut[0]) +
2017 SQR(xm[1] - xcut[1]) +
2018 SQR(xm[2] - xcut[2]) ) / Speed_of_Light_air_cmns);
2019
2020 // add path from the mirror till the camera
2021
2022 t = t + sqrt( SQR(xcutCT[0] - xcam[0]) +
2023 SQR(xcutCT[1] - xcam[1]) +
2024 SQR(xcutCT[2] - xcam[2]) ) / Speed_of_Light_air_cmns;
2025
2026 // save data in the photon variable
2027
2028 cphoton.fill( photon.get_id(),
2029 xcam[0], xcam[1],
2030 photon.get_u(), photon.get_v(),
2031 photon.get_h(), t, phi);
2032
2033 // show it
2034
2035 if ( verbose >= VERBOSE_MAXIMAL ) {
2036 cout << '@' << '2'
2037 << ' ' << ncer
2038 << ' ' << xCT[0]
2039 << ' ' << xCT[1]
2040 << ' ' << xCT[2]
2041 << ' ' << rCT[0]
2042 << ' ' << rCT[1]
2043 << ' ' << rCT[2]
2044 << ' ' << xcut[0]
2045 << ' ' << xcut[1]
2046 << ' ' << xcut[2]
2047 << endl;
2048 cout << '@' << '3'
2049 << ' ' << ncer
2050 << ' ' << sx
2051 << ' ' << sy
2052 << ' ' << i_mirror
2053 << ' ' << ct_data[i_mirror][CT_SX]
2054 << ' ' << ct_data[i_mirror][CT_SY]
2055 << ' ' << ct_data[i_mirror][CT_SX] - sx
2056 << ' ' << ct_data[i_mirror][CT_SY] - sy
2057 << endl;
2058 if ( pb_ParallelBeam == FALSE ) {
2059 cout << '@' << '4'
2060 << ' ' << ncer
2061 << ' ' << xcam[0]
2062 << ' ' << xcam[1]
2063 << ' ' << xcam[2]
2064 << ' ' << phi
2065 << ' ' << photon.get_t()-stat.get_tfirst()
2066 << ' ' << t-stat.get_tfirst()
2067 << endl << flush;
2068 } else {
2069 cout << '@' << '4'
2070 << ' ' << ncer
2071 << ' ' << xcam[0]
2072 << ' ' << xcam[1]
2073 << ' ' << xcam[2]
2074 << ' ' << phi
2075 << endl << flush;
2076 }
2077 }
2078
2079 // write the photon data to the output file
2080
2081 if ( Data_To_STDOUT )
2082 cout.write ( (char *)&cphoton, cphoton.mysize() );
2083 else
2084 cphoton.write( outputfile );
2085
2086 // increase number of Cphotons written
2087
2088 ++nCphotons;
2089
2090 } // while still there are photons left
2091
2092 if ( verbose >= VERBOSE_NORMAL )
2093 cout << nCphotons << ' '
2094 << nbeforeTR << ' '
2095 << nafterTR << ' '
2096 << endl << flush;
2097
2098 // write end-of-event mark
2099
2100 if ( Data_To_STDOUT )
2101 cout.write( FLAG_END_OF_EVENT, SIZE_OF_FLAGS );
2102 else
2103 outputfile.write( FLAG_END_OF_EVENT, SIZE_OF_FLAGS );
2104
2105 // close files
2106
2107 cerfile.close();
2108 stat.closefile();
2109
2110 // show how many photons were written
2111
2112 if ( verbose >= VERBOSE_MAXIMAL )
2113 log( SIGNATURE, "%d C-photons written.\n", nCphotons );
2114
2115 // if we are reading data from STDIN, delete the last used files
2116
2117 if ( Data_From_STDIN ) {
2118
2119 unlink( cername );
2120 unlink( staname );
2121
2122 }
2123
2124 } // while there are still data left
2125
2126 // write end-of-run mark
2127
2128 if ( Data_To_STDOUT )
2129 cout.write( FLAG_END_OF_RUN, SIZE_OF_FLAGS );
2130 else
2131 outputfile.write( FLAG_END_OF_RUN, SIZE_OF_FLAGS );
2132
2133 // close directory
2134
2135 closedir( directory );
2136
2137 } // for each data directory
2138
2139 // write end-of-file mark
2140
2141 if ( Data_To_STDOUT )
2142 cout.write( FLAG_END_OF_FILE, SIZE_OF_FLAGS );
2143 else
2144 outputfile.write( FLAG_END_OF_FILE, SIZE_OF_FLAGS );
2145
2146 // close output file
2147
2148 if ( verbose >= VERBOSE_MINIMAL )
2149 log( SIGNATURE, "Closing output file %s\n", outname );
2150
2151 if ( ! Data_To_STDOUT )
2152 outputfile.close();
2153
2154 // clean everything
2155
2156 if ( verbose >= VERBOSE_MINIMAL )
2157 log( SIGNATURE, "Cleanning . . .\n" );
2158
2159 clean();
2160
2161 // program finished
2162
2163 if ( verbose >= VERBOSE_MINIMAL )
2164 log( SIGNATURE, "Done.\n");
2165
2166 return ( 0 );
2167}
2168//!@}
2169
2170// @T \newpage
2171
2172//=---------------------------------------------------------------------
2173// @subsection Functions definition
2174
2175//!---------------------------------------------------------------------
2176// @name present
2177//
2178// @desc Make some presentation
2179//
2180// @date Sat Jun 27 05:58:56 MET DST 1998
2181//----------------------------------------------------------------------
2182// @function
2183
2184//!@{
2185void
2186present(void)
2187{
2188 cout << "##################################################\n"
2189 << SIGNATURE << '\n' << '\n'
2190 << "Simulation of the reflector\n"
2191 << "J C Gonzalez, May 1998\n"
2192 << "##################################################\n\n";
2193}
2194//!@}
2195
2196
2197//!---------------------------------------------------------------------
2198// @name usage
2199//
2200// @desc show help
2201//
2202// @date Tue Dec 15 16:23:30 MET 1998
2203//----------------------------------------------------------------------
2204// @function
2205
2206//!@{
2207void
2208usage(void)
2209{
2210 present();
2211 cout << "\nusage ::\n\n"
2212 << "\t reflector "
2213 << " [ -@ paramfile ] "
2214 << " [ -h ] "
2215 << "\n\n or \n\n"
2216 << "\t reflector < paramfile"
2217 << "\n\n";
2218 exit(0);
2219}
2220//!@}
2221
2222
2223//!---------------------------------------------------------------------
2224// @name clean
2225//
2226// @desc basically, frees memory
2227//
2228// @date Sat Jun 27 05:58:56 MET DST 1998
2229//----------------------------------------------------------------------
2230// @function
2231
2232//!@{
2233void
2234clean(void)
2235{
2236 int i;
2237
2238 // delete focals table
2239
2240 delete [] ct_Focal;
2241 ct_Focal = 0;
2242
2243 // delete reflectivity table
2244
2245 for (i=0; i<2; i++) {
2246 delete [] Reflectivity[i];
2247 }
2248
2249 delete [] Reflectivity;
2250 Reflectivity = 0;
2251
2252 // delete mirrors' data table
2253
2254 for (i=0; i<ct_NMirrors; i++) {
2255 delete [] ct_data[i];
2256 }
2257
2258 delete [] ct_data;
2259 ct_data = 0;
2260}
2261//!@}
2262
2263
2264//!---------------------------------------------------------------------
2265// @name log
2266//
2267// @desc function to send log information
2268//
2269// @var funct Name of the caller function
2270// @var fmt Format to be used (message)
2271// @var ... Other information to be shown
2272//
2273// @date Sat Jun 27 05:58:56 MET DST 1998
2274//----------------------------------------------------------------------
2275// @function
2276
2277//!@{
2278void
2279log(const char *funct, char *fmt, ...)
2280{
2281 va_list args;
2282
2283 // Display the name of the function that called error
2284 printf("[%s]: ", funct);
2285
2286 // Display the remainder of the message
2287 va_start(args, fmt);
2288 vprintf(fmt, args);
2289 va_end(args);
2290}
2291//!@}
2292
2293
2294//!---------------------------------------------------------------------
2295// @name logerr
2296//
2297// @desc function to send log information to stderr
2298//
2299// @var funct Name of the caller function
2300// @var fmt Format to be used (message)
2301// @var ... Other information to be shown
2302//
2303// @date Sat Jun 27 05:58:56 MET DST 1998
2304//----------------------------------------------------------------------
2305// @function
2306
2307//!@{
2308void
2309logerr(const char *funct, char *fmt, ...)
2310{
2311 va_list args;
2312
2313 // Display the name of the function that called error
2314 printf("[%s]: ", funct);
2315
2316 // Display the remainder of the message
2317 va_start(args, fmt);
2318 vfprintf(stderr, fmt, args);
2319 va_end(args);
2320}
2321//!@}
2322
2323
2324//!---------------------------------------------------------------------
2325// @name error
2326//
2327// @desc function to send an error message, and abort the program
2328//
2329// @var funct Name of the caller function
2330// @var fmt Format to be used (message)
2331// @var ... Other information to be shown
2332//
2333// @date Sat Jun 27 05:58:56 MET DST 1998
2334//----------------------------------------------------------------------
2335// @function
2336
2337//!@{
2338void
2339error(const char *funct, char *fmt, ...)
2340{
2341 va_list args;
2342
2343 // Display the name of the function that called error
2344 fprintf(stderr, "ERROR in %s: ", funct);
2345
2346 // Display the remainder of the message
2347 va_start(args, fmt);
2348 vfprintf(stderr, fmt, args);
2349 va_end(args);
2350
2351 abort();
2352}
2353//!@}
2354
2355
2356//!---------------------------------------------------------------------
2357// @name makeOmega
2358//
2359// @desc function to calculate the matrix Omega(theta,phi)
2360//
2361// @var theta Angle theta of the transformation
2362// @var phi Angle phi of the transformation
2363//
2364// @date Sat Jun 27 05:58:56 MET DST 1998
2365//----------------------------------------------------------------------
2366// @function
2367
2368//!@{
2369void
2370makeOmega (float theta, float phi)
2371{
2372 static float ct, st, cp, sp;
2373
2374 // shortcuts for cosine and sine of theta and phi
2375 ct = cos(theta);
2376 st = sin(theta);
2377 cp = cos(phi);
2378 sp = sin(phi);
2379
2380 // save values in the array (see top of file)
2381 Omega[0][0] = cp*ct;
2382 Omega[0][1] = sp*ct;
2383 Omega[0][2] = -st;
2384
2385 Omega[1][0] = -sp;
2386 Omega[1][1] = cp;
2387 Omega[1][2] = 0;
2388
2389 Omega[2][0] = cp*st;
2390 Omega[2][1] = sp*st;
2391 Omega[2][2] = ct;
2392}
2393//!@}
2394
2395
2396//!---------------------------------------------------------------------
2397// @name makeOmegaI
2398//
2399// @desc function to calculate the matrix Omega-1(theta,phi)
2400//
2401// @var theta Angle theta of the transformation
2402// @var phi Angle phi of the transformation
2403//
2404// @date Sat Jun 27 05:58:56 MET DST 1998
2405//----------------------------------------------------------------------
2406// @function
2407
2408//!@{
2409void
2410makeOmegaI(float theta, float phi)
2411{
2412 static float ct, st, cp, sp;
2413
2414 // shortcuts for cosine and sine of theta and phi
2415 ct = cos(theta);
2416 st = sin(theta);
2417 cp = cos(phi);
2418 sp = sin(phi);
2419
2420 // save values in the array (see top of file)
2421 OmegaI[0][0] = cp*ct;
2422 OmegaI[0][1] = -sp;
2423 OmegaI[0][2] = cp*st;
2424
2425 OmegaI[1][0] = sp*ct;
2426 OmegaI[1][1] = cp;
2427 OmegaI[1][2] = sp*st;
2428
2429 OmegaI[2][0] = -st;
2430 OmegaI[2][1] = 0;
2431 OmegaI[2][2] = ct;
2432}
2433//!@}
2434
2435
2436//!---------------------------------------------------------------------
2437// @name applyMxv
2438//
2439// @desc returns the vector v' such that v' = M x v
2440//
2441// @var M matrix of the transformation
2442// @var v vector to be multiplied
2443// @var vi resulting vector
2444//
2445// @date Sat Jun 27 05:58:56 MET DST 1998
2446//----------------------------------------------------------------------
2447// @function
2448
2449//!@{
2450void
2451applyMxV(float M[3][3], float *V, float *Vp)
2452{
2453 Vp[0] = (M[0][0] * V[0] +
2454 M[0][1] * V[1] +
2455 M[0][2] * V[2]);
2456 Vp[1] = (M[1][0] * V[0] +
2457 M[1][1] * V[1] +
2458 M[1][2] * V[2]);
2459 Vp[2] = (M[2][0] * V[0] +
2460 M[2][1] * V[1] +
2461 M[2][2] * V[2]);
2462}
2463//!@}
2464
2465
2466//!---------------------------------------------------------------------
2467// @name read_ct_file
2468//
2469// @desc read CT definition file
2470//
2471// @date Sat Jun 27 05:58:56 MET DST 1998
2472//----------------------------------------------------------------------
2473// @function
2474
2475//!@{
2476void
2477read_ct_file(void)
2478{
2479 char line[LINE_MAX_LENGTH]; // line to get from the ctin
2480 char token[ITEM_MAX_LENGTH]; // a single token
2481 char fmirr[40]; // name of the binary version of mirrors data
2482 int i; // dummy counters
2483
2484 if ( verbose >= VERBOSE_MINIMAL )
2485 log( "read_ct_file", "start.\n" );
2486
2487 ifstream ctin ( get_ct_filename() );
2488
2489 if ( ctin.bad() )
2490 error( "read_ct_file",
2491 "Cannot open CT def. file: %s\n", get_ct_filename() );
2492
2493 // loop till the "end" directive is reached
2494
2495 while (!ctin.eof()) {
2496
2497 // get line from stdin
2498
2499 ctin.getline(line, LINE_MAX_LENGTH);
2500
2501 // look for each item at the beginning of the line
2502
2503 for (i=0; i<=define_mirrors; i++)
2504 if (strstr(line, CT_ITEM_NAMES[i]) == line)
2505 break;
2506
2507 // if it is not a valid line, just ignore it
2508
2509 if (i == define_mirrors+1)
2510 continue;
2511
2512 // case block for each directive
2513
2514 switch ( i ) {
2515
2516 case type: // <type of telescope> (0:CT1 ¦ 1:MAGIC)
2517
2518 // get focal distance
2519
2520 sscanf(line, "%s %d", token, &ct_Type);
2521
2522 if ( verbose >= VERBOSE_MINIMAL )
2523 log( "read_ct_file", "<Type of Telescope>: %s\n",
2524 ((ct_Type==0) ? "CT1" : "MAGIC") );
2525
2526 break;
2527
2528 case focal_distance: // <focal distance> [cm]
2529
2530 // get focal distance
2531
2532 sscanf(line, "%s %f", token, &ct_Focal_mean);
2533
2534 if ( verbose >= VERBOSE_MINIMAL )
2535 log( "read_ct_file", "<Focal distance>: %f cm\n", ct_Focal_mean );
2536
2537 break;
2538
2539 case focal_std: // s(focal distance) [cm]
2540
2541 // get focal distance
2542
2543 sscanf(line, "%s %f", token, &ct_Focal_std);
2544
2545 if ( verbose >= VERBOSE_MINIMAL )
2546 log( "read_ct_file", "s(Focal distance): %f cm\n", ct_Focal_std );
2547
2548 break;
2549
2550 case point_spread: // <point spread> [cm]
2551
2552 // get point spread
2553
2554 sscanf(line, "%s %f", token, &ct_PSpread_mean);
2555
2556 if ( verbose >= VERBOSE_MINIMAL )
2557 log( "read_ct_file", "<Point spread>: %f cm\n", ct_PSpread_mean );
2558
2559 break;
2560
2561 case point_std: // s(point spread) [cm]
2562
2563 // get point spread
2564
2565 sscanf(line, "%s %f", token, &ct_PSpread_std);
2566
2567 if ( verbose >= VERBOSE_MINIMAL )
2568 log( "read_ct_file", "s(Point spread): %f cm\n", ct_PSpread_std );
2569
2570 break;
2571
2572 case adjustment_dev: // s(adjustment_dev) [cm]
2573
2574 // get point spread
2575
2576 sscanf(line, "%s %f", token, &ct_Adjustment_std);
2577
2578 if ( verbose >= VERBOSE_MINIMAL )
2579 log( "read_ct_file", "s(Adjustment): %f cm\n", ct_Adjustment_std );
2580
2581 break;
2582
2583 case black_spot: // radius of the black spot in the center [cm]
2584
2585 // get black spot radius
2586
2587 sscanf(line, "%s %f", token, &ct_BlackSpot_rad);
2588
2589 if ( verbose >= VERBOSE_MINIMAL )
2590 log( "read_ct_file", "Radius of the black spots: %f cm\n",
2591 ct_BlackSpot_rad);
2592
2593 break;
2594
2595 case r_mirror: // radius of the mirrors [cm]
2596
2597 // get radius of mirror
2598
2599 sscanf(line, "%s %f", token, &ct_RMirror);
2600
2601 if ( verbose >= VERBOSE_MINIMAL )
2602 log( "read_ct_file", "Radii of the mirrors: %f cm\n", ct_RMirror );
2603
2604 break;
2605
2606 case n_mirrors: // number of mirrors
2607
2608 // get the name of the output_file from the line
2609
2610 sscanf(line, "%s %d", token, &ct_NMirrors);
2611
2612 if ( verbose >= VERBOSE_MINIMAL )
2613 log( "read_ct_file", "Number of mirrors: %d\n", ct_NMirrors );
2614
2615 break;
2616
2617 case camera_width: // camera width [cm]
2618
2619 // get the name of the ct_file from the line
2620
2621 sscanf(line, "%s %f", token, &ct_CameraWidth);
2622
2623 if ( verbose >= VERBOSE_MINIMAL )
2624 log( "read_ct_file", "Camera width: %f cm\n", ct_CameraWidth );
2625
2626 break;
2627
2628 case n_pixels: // number of pixels
2629
2630 // get the name of the output_file from the line
2631
2632 sscanf(line, "%s %d", token, &ct_NPixels);
2633
2634 if ( verbose >= VERBOSE_MINIMAL )
2635 log( "read_ct_file", "Number of pixels: %d\n", ct_NPixels );
2636
2637 break;
2638
2639 case pixel_width: // pixel width [cm]
2640
2641 // get the name of the ct_file from the line
2642
2643 sscanf(line, "%s %f", token, &ct_PixelWidth);
2644
2645 if ( verbose >= VERBOSE_MINIMAL )
2646 log( "read_ct_file", "Pixel width: %f cm\n", ct_PixelWidth );
2647
2648 break;
2649
2650 case define_mirrors: // read table with the parameters of the mirrors
2651
2652 if ( verbose >= VERBOSE_MINIMAL )
2653 log( "read_ct_file", "Table of mirrors data:\n" );
2654
2655 // check whether the number of mirrors was already set
2656
2657 if ( ct_NMirrors == 0 )
2658 error( "read_ct_file", "NMirrors was not set.\n" );
2659
2660 // allocate memory for paths list
2661
2662 if ( verbose >= VERBOSE_MINIMAL )
2663 log( "read_ct_file", "Allocating memory for ct_data\n" );
2664
2665 ct_data = new float*[ct_NMirrors];
2666
2667 for (i=0; i<ct_NMirrors; i++)
2668 ct_data[i] = new float[CT_NDATA];
2669
2670 // read data
2671
2672 // look for binary version of mirror data (faster to load)
2673
2674 sprintf( fmirr, "%s.mirr", get_ct_filename() );
2675
2676 ifstream ctmirr;
2677 ctmirr.open( fmirr );
2678
2679 if ( ctmirr.bad() ) {
2680
2681 // the file does not exist
2682
2683 ctmirr.close(); // close not existing file
2684
2685 ofstream ctmirr_write ( fmirr ); // open file to save the data
2686
2687 if ( verbose >= VERBOSE_MINIMAL )
2688 log( "read_ct_file", "Reading mirrors data...\n" );
2689
2690 // read ASCII data
2691
2692 for (i=0; i<ct_NMirrors; i++) {
2693 // for (j=0; j<CT_NDATA; j++)
2694
2695 ctin.getline(line, LINE_MAX_LENGTH);
2696 sscanf(line, "%f %f %f %f %f %f %f %f %f %f %f %f",
2697 &ct_data[i][0],
2698 &ct_data[i][1],
2699 &ct_data[i][2],
2700 &ct_data[i][3],
2701 &ct_data[i][4],
2702 &ct_data[i][5],
2703 &ct_data[i][6],
2704 &ct_data[i][7],
2705 &ct_data[i][8],
2706 &ct_data[i][9],
2707 &ct_data[i][10],
2708 &ct_data[i][11]);
2709
2710 cout << '[' << i << ']' << flush;
2711 }
2712 cout << endl << flush;
2713
2714 // save binary data and close output file
2715
2716 for (i=0; i<ct_NMirrors; i++)
2717 ctmirr_write.write( (char*)ct_data[i], CT_NDATA*sizeof(float) );
2718
2719 ctmirr_write.close();
2720
2721 } else {
2722
2723 // Excellent!
2724 // read binary data and close input file
2725
2726 for (i=0; i<ct_NMirrors; i++)
2727 ctmirr.read( (char*)ct_data[i], CT_NDATA*sizeof(float) );
2728
2729 ctmirr.close();
2730
2731 }
2732
2733 break;
2734
2735 } // switch ( i )
2736
2737 } // while (! is_end)
2738
2739 // read reflectivity file
2740
2741 read_reflectivity();
2742
2743 // read axis deviations file
2744
2745 read_axisdev();
2746
2747 // read focals file
2748
2749 read_focals();
2750
2751 // end
2752
2753 if ( verbose >= VERBOSE_MINIMAL )
2754 log( "read_ct_file", "done.\n" );
2755 return;
2756}
2757//!@}
2758
2759
2760//!---------------------------------------------------------------------
2761// @name dist_r_P
2762//
2763// @desc distance straight line r - point P
2764//
2765// @var a coord. X of a fixed point of the straight line
2766// @var b coord. Y of a fixed point of the straight line
2767// @var c coord. Z of a fixed point of the straight line
2768// @var l component X of a vector of the straight line
2769// @var m component Y of a vector of the straight line
2770// @var n component Z of a vector of the straight line
2771// @var x coord. X of the point P
2772// @var y coord. Y of the point P
2773// @var z coord. Z of the point P
2774//
2775// @return Returns the distance d(r,P) between the line and the point P
2776//
2777// @date Sat Jun 27 05:58:56 MET DST 1998
2778//----------------------------------------------------------------------
2779// @function
2780
2781//!@{
2782float
2783dist_r_P(float a, float b, float c,
2784 float l, float m, float n,
2785 float x, float y, float z)
2786{
2787 return (
2788 sqrt((SQR((a-x)*m-(b-y)*l) +
2789 SQR((b-y)*n-(c-z)*m) +
2790 SQR((c-z)*l-(a-x)*n))/
2791 (SQR(l)+SQR(m)+SQR(n))
2792 )
2793 );
2794}
2795//!@}
2796
2797
2798//!---------------------------------------------------------------------
2799// @name read_reflectivity
2800//
2801// @desc read reflectivity data for the mirrors
2802//
2803// @date Sat Jun 27 05:58:56 MET DST 1998
2804//----------------------------------------------------------------------
2805// @function
2806
2807//!@{
2808void
2809read_reflectivity(void)
2810{
2811 ifstream reflfile;
2812 ofstream reflfile_write;
2813 int i, errors=0;
2814 char line[LINE_MAX_LENGTH];
2815
2816 // try to open the file
2817 reflfile.open( REFLECTIVITY_FILE );
2818
2819 // if it is wrong or does not exist:
2820
2821 while ( reflfile.bad() ) {
2822
2823 ++errors;
2824
2825 if ( verbose >= VERBOSE_MINIMAL )
2826 log( "read_reflectivity", "Cannot open reflectivity file: %s\n",
2827 REFLECTIVITY_FILE );
2828
2829 if ( errors > 5 )
2830 error( "read_reflectivity", "Exiting.");
2831
2832 // try to re-generate the file
2833
2834 if ( verbose >= VERBOSE_MINIMAL )
2835 log( "read_reflectivity", "Generating file: %s\n",
2836 REFLECTIVITY_FILE );
2837
2838 // try to open the file
2839 reflfile_write.open( REFLECTIVITY_FILE );
2840
2841 // write short comment in the beginning
2842
2843 reflfile_write << "# reflectivity file" << endl
2844 << "# reflectivity for each mirror in the frame" << endl
2845 << "# J C Gonzalez, 1998" << endl
2846 << "#" << endl;
2847
2848 // write data
2849
2850 // write number of data-points
2851
2852 reflfile_write << "# number of datapoints" << endl
2853 << 35 << endl;
2854
2855 // write data-points
2856
2857 reflfile_write << "# datapoints" << endl;
2858
2859 for ( i=0; i<35; ++i )
2860 reflfile_write << 270.0+i*10.0 << ' ' << .9 << endl;
2861
2862 // close it and try to re-open
2863
2864 reflfile_write.close();
2865 reflfile.close();
2866 reflfile.open( REFLECTIVITY_FILE );
2867
2868 }
2869
2870 if ( verbose >= VERBOSE_MINIMAL ) {
2871 if ( errors > 0 )
2872 log( "read_reflectivity", "Reading reflectivity file after %d fail%c\n",
2873 errors, ((errors>1)?'s':' ') );
2874 else
2875 log( "read_reflectivity", "Reading reflectivity file . . .\n" );
2876 }
2877
2878 // we set this to -1, in order to read in the first place the
2879 // number of data-points for the Reflectivity curve
2880 i = -1;
2881
2882 // scan the whole file
2883
2884 while ( ! reflfile.eof() ) {
2885
2886 // get line from the file
2887
2888 reflfile.getline(line, LINE_MAX_LENGTH);
2889
2890 // skip if comment
2891
2892 if ( *line == '#' )
2893 continue;
2894
2895 // get the value
2896
2897 if ( i < 0 ) { // if it is the first number we read
2898
2899 // read the number of datapoints
2900
2901 sscanf(line, "%d", &nReflectivity);
2902
2903 // allocate memory for Reflectivity table
2904
2905 Reflectivity = new float * [2];
2906 Reflectivity[0] = new float[nReflectivity];
2907 Reflectivity[1] = new float[nReflectivity];
2908
2909 } else { // if not, then it is a datapoint
2910
2911 // get the datapoint (wavelength, reflectivity)
2912
2913 sscanf(line, "%f %f", &Reflectivity[0][i], &Reflectivity[1][i]);
2914
2915 }
2916
2917 ++i;
2918
2919 }
2920
2921 // close it
2922
2923 reflfile.close();
2924
2925}
2926//!@}
2927
2928
2929//!---------------------------------------------------------------------
2930// @name read_axisdev
2931//
2932// @desc read axis deviations data for the mirrors
2933//
2934// @date Sat Jun 27 05:58:56 MET DST 1998
2935//----------------------------------------------------------------------
2936// @function
2937
2938//!@{
2939void
2940read_axisdev(void)
2941{
2942 ifstream axisdevfile;
2943 ofstream axisdevfile_write;
2944 int i, errors=0;
2945 char line[LINE_MAX_LENGTH];
2946
2947 // try to open the file
2948 axisdevfile.open( AXISDEVIATION_FILE );
2949
2950 // if it is wrong or does not exist:
2951
2952 while ( axisdevfile.bad() ) {
2953
2954 ++errors;
2955
2956 if ( verbose >= VERBOSE_MINIMAL )
2957 log( "read_axisdev", "Cannot open axis deviation file: %s\n",
2958 AXISDEVIATION_FILE );
2959
2960 if ( errors > 5 )
2961 error( "read_axisdev", "Exiting.");
2962
2963 // try to re-generate the file
2964
2965 if ( verbose >= VERBOSE_MINIMAL )
2966 log( "read_axisdev", "Generating file: %s\n",
2967 AXISDEVIATION_FILE );
2968
2969 // try to open the file
2970 axisdevfile_write.open( AXISDEVIATION_FILE );
2971
2972 // write short comment in the beginning
2973
2974 axisdevfile_write << "# axis deviation file" << endl
2975 << "# axis deviation for each mirror" << endl
2976 << "# J C Gonzalez, 1998" << endl
2977 << "#" << endl;
2978
2979 // write data
2980
2981 for ( i=0; i<ct_NMirrors; ++i )
2982 axisdevfile_write << gennor(0.0,ct_Adjustment_std)
2983 << ' '
2984 << gennor(0.0,ct_Adjustment_std)
2985 << endl;
2986
2987 // close it and try to re-open
2988
2989 axisdevfile_write.close();
2990 axisdevfile.close();
2991 axisdevfile.open( AXISDEVIATION_FILE );
2992
2993 }
2994
2995 if ( verbose >= VERBOSE_MINIMAL ) {
2996 if ( errors > 0 )
2997 log( "read_axisdev", "Reading axis deviation file after %d fail%c\n",
2998 errors, ((errors>1)?'s':' ') );
2999 else
3000 log( "read_axisdev", "Reading axis deviation file . . .\n" );
3001 }
3002
3003 //allocate memory for AxisDeviation table
3004
3005 AxisDeviation = new float * [2];
3006 AxisDeviation[0] = new float[ct_NMirrors];
3007 AxisDeviation[1] = new float[ct_NMirrors];
3008
3009 // scan the whole file
3010
3011 i = 0;
3012
3013 while ( (! axisdevfile.eof()) && (i<ct_NMirrors) ) {
3014
3015 // get line from the file
3016
3017 axisdevfile.getline(line, LINE_MAX_LENGTH);
3018
3019 // skip if comment
3020
3021 if ( *line == '#' )
3022 continue;
3023
3024 // get the value (dx, dy)
3025
3026 sscanf(line, "%f %f", &AxisDeviation[0][i], &AxisDeviation[1][i]);
3027
3028 ++i;
3029
3030 }
3031
3032 // did i do it for all the mirrors ?
3033
3034 if (i < ct_NMirrors) {
3035 error("read_axis_dev", "%s%s (%d < %d)\n",
3036 "Number of pairs in axisdev.dat file is",
3037 "\n\t\tsmaller than number of mirrors", i, ct_NMirrors);
3038 }
3039
3040 // close it
3041
3042 axisdevfile.close();
3043
3044}
3045//!@}
3046
3047
3048//!---------------------------------------------------------------------
3049// @name read_focals
3050//
3051// @desc read focals data for the mirrors
3052//
3053// @date Sat Jun 27 05:58:56 MET DST 1998
3054//----------------------------------------------------------------------
3055// @function
3056
3057//!@{
3058void
3059read_focals(void)
3060{
3061 ifstream focalfile;
3062 ofstream focalfile_write;
3063 int i, errors=0;
3064 char line[LINE_MAX_LENGTH];
3065
3066 /*!@'
3067
3068 Here we read the focals of each mirror, but only for CT1. For
3069 MAGIC we have already the Focals in the definition file. So, we
3070 copy these values to the |ct\_Focals| vector.
3071
3072 */
3073
3074 // >>>>>>>>>> MAGIC <<<<<<<<<<
3075 if ( ct_Type == 1) {
3076
3077 // allocate memory for the focals table
3078 ct_Focal = new float[ct_NMirrors];
3079
3080 for (i=0; i<ct_NMirrors; ++i)
3081 ct_Focal[i] = ct_data[i][CT_FOCAL];
3082
3083 return;
3084
3085 }
3086 // ! >>>>>>>>>> MAGIC <<<<<<<<<<
3087
3088 // try to open the file
3089 focalfile.open( FOCALS_FILE );
3090
3091 // if it is wrong or does not exist:
3092
3093 while ( focalfile.bad() ) {
3094
3095 ++errors;
3096
3097 if ( verbose >= VERBOSE_MINIMAL )
3098 log( "read_focals", "Cannot open focals file: %s\n",
3099 FOCALS_FILE );
3100
3101 if ( errors > 5 )
3102 error( "read_focals", "Exiting.");
3103
3104 // try to re-generate the file
3105
3106 if ( verbose >= VERBOSE_MINIMAL )
3107 log( "read_focals", "Generating file: %s\n",
3108 FOCALS_FILE );
3109
3110 // try to open the file
3111 focalfile_write.open( FOCALS_FILE );
3112
3113 // write short comment in the beginning
3114
3115 focalfile_write << "# focals file" << endl
3116 << "# focals for each mirror in the frame" << endl
3117 << "# J C Gonzalez, 1998" << endl;
3118
3119 // write data
3120
3121 rnormal( NormalRandomNumbers, ct_NMirrors+1 );
3122
3123 for (i=0; i<ct_NMirrors; i++)
3124 focalfile_write << NormalRandomNumbers[i] * ct_Focal_std + ct_Focal_mean
3125 << endl;
3126
3127 // close it and try to re-open
3128
3129 focalfile_write.close();
3130 focalfile.close();
3131 focalfile.open( FOCALS_FILE );
3132
3133 }
3134
3135 // allocate memory for the focals table
3136
3137 ct_Focal = new float[ct_NMirrors];
3138
3139 if ( verbose >= VERBOSE_MINIMAL ) {
3140
3141 if ( errors > 0 )
3142 log( "read_focals", "Reading focals file after %d fail%c\n",
3143 errors, ((errors>1) ? 's' : ' ') );
3144
3145 else
3146
3147 log( "read_focals", "Reading focals file . . .\n" );
3148 }
3149
3150 // scan the whole file
3151
3152 i = 0;
3153
3154 while ( (! focalfile.eof()) && (i<ct_NMirrors) ) {
3155
3156 // get line from the file
3157
3158 focalfile.getline(line, LINE_MAX_LENGTH);
3159
3160 // skip if comment
3161
3162 if ( *line == '#' )
3163 continue;
3164
3165 // get the value
3166
3167 sscanf(line, "%f", &ct_Focal[i]);
3168
3169 ++i;
3170
3171 }
3172
3173 // did i do it for all the mirrors ?
3174
3175 if (i < ct_NMirrors) {
3176 error("read_focals", "%s%s (%d < %d)\n",
3177 "Number of data in focals.dat file is",
3178 "\n\t\tsmaller than number of mirrors", i, ct_NMirrors);
3179 }
3180
3181 // close it
3182
3183 focalfile.close();
3184
3185}
3186//!@}
3187
3188
3189//!---------------------------------------------------------------------
3190// @name rnormal
3191//
3192// @desc returns n(=2k) normaly distributed numbers
3193//
3194// @var *r pointer to a vector where we write the numbers
3195// @var n how many numbers do we generate
3196//
3197// @date Sat Jun 27 05:58:56 MET DST 1998
3198//----------------------------------------------------------------------
3199// @function
3200
3201//!@{
3202void
3203rnormal(double *r, int n)
3204{
3205
3206 double z1, z2;
3207 int i;
3208
3209 for (i=0; i<n; i+=2) {
3210
3211 z1 = RandomNumber;
3212 z2 = RandomNumber;
3213
3214 r[i] = sqrt(-2.0*log(z1)) * cos(2.0*M_PI*z2);
3215 r[i+1] = sqrt(-2.0*log(z1)) * sin(2.0*M_PI*z2);
3216
3217 }
3218
3219}
3220//!@}
3221
3222
3223//!---------------------------------------------------------------------
3224// @name isA
3225//
3226// @desc returns TRUE(FALSE), if the flag is(is not) the given
3227//
3228// @var s1 String to be compared with s2
3229// @var s2 String in the form of a FLAG
3230//
3231// @return TRUE: both strings match; FALSE: otherwise
3232//
3233// @date Wed Jul 8 15:25:39 MET DST 1998
3234//----------------------------------------------------------------------
3235// @function
3236
3237//!@{
3238int
3239isA( char * s1, const char * flag )
3240{
3241 return ( (strncmp((char *)s1, flag, 8)==0) ? TRUE : FALSE );
3242}
3243//!@}
3244
3245
3246//!---------------------------------------------------------------------
3247// @name Curv2Lin
3248//
3249// @desc Curvilinear to Linear (Euclidean) distance
3250//
3251// @var s Curvilinear distance over the parabolic shape
3252//
3253// @return Radial distance from the axis of the paraboloid
3254//
3255// @date Wed Jul 8 15:25:39 MET DST 1998
3256//----------------------------------------------------------------------
3257// @function
3258
3259//!@{
3260float
3261Curv2Lin(float s)
3262{
3263 float x;
3264 short i;
3265
3266 x = s;
3267 for (i = 0; i < 4; i++)
3268 x = (s / 100.) / (1. + (float) 0.000144175317185 * x * x);
3269
3270 return (x*100.);
3271}
3272//!@}
3273
3274
3275//!---------------------------------------------------------------------
3276// @name Lin2Curv
3277//
3278// @desc Linear (Euclidean) to Curvilinear distance
3279//
3280// @var x Radial distance from the axis of the paraboloid
3281//
3282// @return Curvilinear distance over the parabolic shape
3283//
3284// @date Wed Jul 8 15:25:39 MET DST 1998
3285//----------------------------------------------------------------------
3286// @function
3287
3288//!@{
3289float
3290Lin2Curv(float x)
3291{
3292 x /= 100.;
3293 return ((x + (float) 0.000144175317185 * x * x * x)*100.);
3294}
3295//!@}
3296
3297
3298//!---------------------------------------------------------------------
3299// @name get_stdin_files
3300//
3301// @desc get data from stdin and create cerXXXXXX and staXXXXXX
3302//
3303// @var ncerf number of next Cherenkov file cerXXXXXX to be created
3304// @var El Lower limit for requested energy
3305// @var Eu Upper limit for requested energy
3306// @var flag Boolean flag: TRUE: initialize; FALSE: normal
3307//
3308// @return TRUE: succesful; FALSE: failure
3309//
3310// @date Tue Dec 15 15:03:40 MET 1998
3311//----------------------------------------------------------------------
3312// @function
3313
3314//!@{
3315int
3316get_stdin_files(int ncerf, float El, float Eu, int flag)
3317{
3318 static float *buffer; // buffer for read-out of STDIN
3319 static COREventHeader *evth;// Event Header class (from CORSIKA)
3320 static float Elow, Eup; // range for Energy
3321 static int reject;
3322
3323 char *chkbuffer;
3324
3325 char cername[256]; // output filename
3326 ofstream cerfile; // output file (stream)
3327 char staname[256]; // output filename sta
3328 ofstream stafile; // output file sta (stream)
3329 int bytes;
3330 float energy;
3331
3332 /*!@'
3333
3334 The main features of this function, as well as the algorithm,
3335 are taken from |corfilter|.
3336
3337 */
3338
3339 // if it's the first time, allocate buffers
3340
3341 if ( flag == TRUE ) {
3342
3343 if (verbose)
3344 log("get_stdin_files", "Allocating memory for buffer...\n");
3345
3346 buffer = new float [BUFFER_LENGTH];
3347 evth = (COREventHeader*)buffer;
3348
3349 Elow = El;
3350 Eup = Eu;
3351
3352 return ( TRUE );
3353 }
3354
3355 // main loop
3356
3357 while ( ! cin.eof() ) {
3358
3359 // read buffer from STDIN
3360
3361 cin.read( (char*) buffer, BUFFER_LENGTH * sizeof(float) );
3362 bytes = cin.gcount();
3363
3364 if (bytes < (BUFFER_LENGTH * sizeof(float))) {
3365 if (verbose)
3366 logerr("get_stdin_files",
3367 "Can read only %d bytes, instead of required %d.\n",
3368 bytes, BUFFER_LENGTH * sizeof(float));
3369 }
3370
3371 if (bytes == 0)
3372 return ( FALSE );
3373
3374 // check whether EVTH is there
3375
3376 chkbuffer = (char*)buffer;
3377 if (strstr(chkbuffer,EVTH) == 0) {
3378 if (verbose)
3379 logerr("get_stdin_files", "EVTH not found in this block\n");
3380 continue;
3381 } else {
3382 if (verbose)
3383 logerr("get_stdin_files", "EVTH was FOUND in this block!!\n");
3384 }
3385
3386
3387 // show log
3388 if (verbose)
3389 logerr("get_stdin_files", " (E=%8.1f) Reading event %d... ",
3390 (float)evth->Etotal, (int)evth->EvtNumber);
3391
3392 // CRITERIA
3393 energy = evth->Etotal;
3394
3395 if ( (energy<Elow) || (energy>Eup) ) {
3396 if (verbose)
3397 log("REJECTED", "\n");
3398 reject = TRUE;
3399 } else {
3400 if (verbose)
3401 log("ACCEPTED", "\n");
3402 reject = FALSE;
3403 }
3404
3405 if ( !reject ) {
3406
3407 // update no. of cer. file
3408 sprintf(cername, "%s/cer%06d", TMP_STDIN_DIR, ncerf+1);
3409
3410 cerfile.open( cername );
3411
3412 if ( verbose )
3413 log("get_stdin_files", " Writing event %d, will be %s\n",
3414 (int)evth->EvtNumber, cername);
3415
3416 // write it
3417 cerfile.write( (char*)buffer, BUFFER_LENGTH * sizeof(float) );
3418
3419 }
3420
3421 // read blocks till new EVTH is found
3422 // (now must be at the beginning of a block)
3423
3424 if (verbose)
3425 log("get_stdin_files", " Reading photons...\n");
3426
3427 bytes = 0;
3428 while ( TRUE ) {
3429
3430 // read photon
3431 cin.read( (char*)buffer, 7 * sizeof(float) );
3432
3433 // the last photon is reached when the EVTH is read, but
3434 // only when bytes % SIZE_OF_BLOCK == 0
3435 // SIZE_OF_BLOCK = 22932
3436 // with this, we only execute strstr once each 819 photons
3437
3438 // are we in a new block?
3439
3440 if ((bytes % SIZE_OF_BLOCK) == 0) {
3441 // is this the last photon?
3442 if ( strstr(chkbuffer, EVTH) == chkbuffer )
3443 break;
3444 }
3445
3446 // if not, write it to the cer-file
3447 if ( !reject )
3448 cerfile.write( (char*)buffer, 7 * sizeof(float) );
3449
3450 bytes += 7 * sizeof(float);
3451
3452 }
3453
3454 if (verbose)
3455 log("get_stdin_files", " Finished: %d bytes\n", bytes);
3456
3457 if ( !reject ) {
3458 // close cer-file
3459 cerfile.close();
3460 }
3461
3462 // at this point, we must be reading sta-file
3463
3464 // read the rest of the sta-file
3465 if (verbose)
3466 log("get_stdin_files", " Reading sta-data...\n");
3467 cin.read( (char*)buffer + 7*sizeof(float), 5956 - 7*sizeof(float));
3468
3469 if ( !reject ) {
3470
3471 // log
3472 if (verbose)
3473 log("get_stdin_files", " Saving sta-file...\n");
3474
3475 // open file
3476 sprintf(staname, "%s/sta%06d", TMP_STDIN_DIR, ncerf+1);
3477 stafile.open( staname );
3478
3479 // write it
3480 stafile.write( (char*)buffer, 5956 );
3481
3482 // close sta-file
3483 stafile.close();
3484
3485 // get out from here, return to main program
3486 return( TRUE );
3487
3488 }
3489
3490 }
3491
3492 return( TRUE );
3493
3494}
3495//!@}
3496
3497
3498//!---------------------------------------------------------------------
3499// @name get_new_ct_pointing
3500//
3501// @desc returns new random position of CT around a given one
3502//
3503// @var theta Theta (ZA) of the shower
3504// @var phi Phi (AZ) of the shower
3505// @var range Maximum allowed distance between CT and shower axis
3506// @var newtheta Resulting Theta of the CT
3507// @var newphi Resulting Phi of the CT
3508//
3509// @return Angular distance between orig. and new directions
3510//
3511// @date Sat Jun 27 05:58:56 MET DST 1998
3512//----------------------------------------------------------------------
3513// @function
3514
3515//!@{
3516float
3517get_new_ct_pointing(float theta, float phi, float range,
3518 float *newtheta, float *newphi)
3519{
3520 float distance;
3521 float it, ip;
3522 float sin_theta, cos_theta;
3523 float sin_newtheta, cos_newtheta;
3524 float sin_iphi, cos_iphi;
3525 float iphi;
3526
3527 // for the moment, we only simulate an uniform distribution,
3528 // since our theta distribution in the generation of events is
3529 // already uniform for hadrons, which are the main targets for
3530 // using this option
3531
3532 it = acos(cos(range) + RandomNumber * (1 - cos(range)));
3533 ip = RandomNumber * 2.0 * M_PI;
3534
3535 if ( theta == 0.0 ) {
3536
3537 *newtheta = it;
3538 *newphi = ip;
3539
3540 } else {
3541
3542 sin_theta = sin(theta);
3543 cos_theta = cos(theta);
3544
3545 cos_newtheta = cos_theta*cos(it) + sin_theta*sin(it)*cos(ip);
3546 *newtheta = acos( cos_newtheta );
3547 sin_newtheta = sin( *newtheta );
3548
3549 sin_iphi = sin(it)*sin(ip) / sin_newtheta;
3550 cos_iphi = (( cos(it) - cos_newtheta * cos_theta ) /
3551 ( sin_newtheta * sin_theta ));
3552
3553 iphi = atan2( sin_iphi, cos_iphi );
3554
3555 *newphi = phi + iphi;
3556
3557 }
3558
3559 return( it );
3560}
3561//!@}
3562
3563
3564//=------------------------------------------------------------
3565//!@subsection Log of this file.
3566
3567//!@{
3568//
3569// $Log: not supported by cvs2svn $
3570// Revision 1.10 2000/01/31 20:53:57 harald
3571// A smaller change concerning the random pointing. I got from Jose Carlos
3572// the information:
3573//
3574// >>The point was that the get_new_ct_pointing routine now generates
3575// >>the right coordinates randomly from the original direction with
3576// >>a maximum deviatiation given by the user, but this new directions
3577// >>were not uniformly distributed. Now they are !!
3578//
3579// need to be check soon!!
3580//
3581// Revision 1.9 2000/01/28 09:19:54 harald
3582// A new version from JoseCarlosGonzalez. The old routine for random_pointing
3583// was not correct. This one should be okay!!
3584//
3585// Revision 1.20 2000/01/27 10:47:54 gonzalez
3586// JAN2000-STABLE
3587//
3588// Revision 1.19 1999/11/19 20:52:31 gonzalez
3589// *** empty log message ***
3590//
3591// Revision 1.18 1999/10/05 11:11:12 gonzalez
3592// Sep.1999
3593//
3594// Revision 1.17 1999/03/24 16:33:01 gonzalez
3595// REFLECTOR 1.1: Release
3596//
3597// Revision 1.16 1999/01/21 16:03:44 gonzalez
3598// Only small modifications
3599//
3600// Revision 1.15 1999/01/19 18:07:16 gonzalez
3601// Bugs in STDIN-STDOUT version corrected.
3602//
3603// Revision 1.14 1999/01/14 17:35:47 gonzalez
3604// Both reading from STDIN (data_from_stdin) and
3605// writing to STDOUT (data_to_STDOUT) working.
3606//
3607// Revision 1.13 1999/01/13 12:41:12 gonzalez
3608// THIS IS A WORKING (last?) VERSION
3609//
3610// Revision 1.12 1998/12/17 16:26:09 gonzalez
3611// *************************************************
3612// *************************************************
3613// WARNING:: Version 1.11 is completly wrong!!
3614// *************************************************
3615// *************************************************
3616//
3617// Revision 1.10 1998/12/15 10:47:22 gonzalez
3618// RELEASE-1.0
3619//
3620// Revision 1.9 1998/11/25 16:30:47 gonzalez
3621// Commit after inclusion of 'Blocking'
3622//
3623//!@}
3624
3625//=EOF
Note: See TracBrowser for help on using the repository browser.