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

Last change on this file since 362 was 349, checked in by harald, 25 years ago
A smaller change concerning the random pointing. I got from Jose Carlos the information: >>The point was that the get_new_ct_pointing routine now generates >>the right coordinates randomly from the original direction with >>a maximum deviatiation given by the user, but this new directions >>were not uniformly distributed. Now they are !! need to be check soon!!
File size: 96.1 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.10 $
23// $Author: harald $
24// $Date: 2000-01-31 20:53:57 $
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 // find data point to be used in Lagrange interpolation (-> k)
1588
1589 FindLagrange(Reflectivity,k,wl);
1590
1591 // if random > reflectivity then goes to the TOP of the loop again
1592
1593 reflec = Lagrange(Reflectivity,k,wl);
1594
1595 if ( RandomNumber > reflec )
1596 continue;
1597
1598 /*!@'
1599
1600 @#### Reflexion in the mirrors.
1601
1602 The reflexion is made now is a somehow dirty way: if the
1603 value of the variable |ct\_Type| is 0, the telescope is CT1,
1604 and if the value is 1, the telescope is MAGIC. Due to the
1605 completly different geometries of the dishes and the
1606 mirrors, now we support two fully separated and different
1607 algorithms.
1608
1609 */
1610
1611 //++
1612 // REFLEXION
1613 //--
1614
1615 // get director cosines x,y and save original trayectory
1616
1617 u = r[0] = photon.get_u();
1618 v = r[1] = photon.get_v();
1619 u=u, v=v; // to avoid "was set but not used" messages
1620
1621 // get position in the ground
1622
1623 x[0] = photon.get_x() - coreX;
1624 x[1] = photon.get_y() - coreY;
1625 x[2] = 0.0;
1626
1627 if ( verbose >= VERBOSE_MAXIMAL ) {
1628 cout << '@' << '1'
1629 << ' ' << ncer
1630 << ' ' << t
1631 << ' ' << x[0]
1632 << ' ' << x[1]
1633 << ' ' << x[2]
1634 << ' ' << r[0]
1635 << ' ' << r[1]
1636 << ' ' << r[2]
1637 << ' ' << thetaCT
1638 << ' ' << phiCT
1639 << endl;
1640 }
1641
1642 // change to the system of the CT
1643
1644 applyMxV( OmegaCT, x, xCT );
1645 applyMxV( OmegaCT, r, rCT );
1646
1647 // before moving to the system of the mirror, for MAGIC,
1648 // first we look whether the photon hits a mirror or not
1649
1650 // calculate the intersection of the trayectory of the photon
1651 // with the GLOBAL DISH !!!
1652 // we reproduce the calculation of the coefficients of the
1653 // second order polynomial in z (=xCT[2]), made with
1654 // Mathematica
1655
1656 /*
1657 * In[1]:= parab:=z-(x^2+y^2)/(4F)
1658 * par1=parab /. {x->x0+u/w(z-z0),y->y0+v/w(z-z0)}
1659 *
1660 * Out[1]=
1661 * u (z - z0) 2 v (z - z0) 2
1662 * (x0 + ----------) + (y0 + ----------)
1663 * w w
1664 * z - ---------------------------------------
1665 * 4 F
1666 *
1667 * In[2]:= CoefficientList[ExpandAll[par1*4F*w^2],z]
1668 *
1669 * Out[2]=
1670 * 2 2 2 2
1671 * {-(w x0 ) - w y0 + 2 u w x0 z0 +
1672 *
1673 * 2 2 2 2
1674 * 2 v w y0 z0 - u z0 - v z0 ,
1675 *
1676 * 2 2
1677 * 4 F w - 2 u w x0 - 2 v w y0 + 2 u z0 +
1678 *
1679 * 2 2 2
1680 * 2 v z0, -u - v }
1681 */
1682
1683 // the z coordinate is calculated
1684
1685 a = - SQR(rCT[0]) - SQR(rCT[1]);
1686 b = 4.0*ct_Focal_mean*SQR(rCT[2])
1687 - 2.0*rCT[0]*rCT[2]*xCT[0] - 2.0*rCT[1]*rCT[2]*xCT[1]
1688 + 2.0*SQR(rCT[0])*xCT[2] + 2.0*SQR(rCT[1])*xCT[2];
1689 c = 2*rCT[0]*rCT[2]*x[0]*x[2] + 2*rCT[1]*rCT[2]*x[1]*x[2]
1690 - SQR(rCT[2])*SQR(x[0]) - SQR(rCT[2])*SQR(x[1])
1691 - SQR(rCT[0])*SQR(x[2]) - SQR(rCT[1])*SQR(x[2]);
1692
1693 if ( fabs(a) < 1.e-6 ) {
1694
1695 // only one value
1696
1697 xcut[2] = -c / b;
1698
1699 } else {
1700
1701 d = sqrt( b*b - 4.0*a*c );
1702
1703 // two possible values for z
1704
1705 t1 = (-b+d) / (2.0*a);
1706 t2 = (-b-d) / (2.0*a);
1707
1708 // z must be the minimum of t1 and t2
1709
1710 xcut[2] = (t1 < t2) ? t1 : t2;
1711
1712 }
1713
1714 // xcut[] is NOW the cut between the GLOBAL dish of MAGIC and
1715 // the trayectory of the photon
1716
1717 xcut[0] = xCT[0] + rCT[0]/rCT[2]*(xcut[2]-xCT[2]);
1718 xcut[1] = xCT[1] + rCT[1]/rCT[2]*(xcut[2]-xCT[2]);
1719
1720 if ( ct_Type == 1 ) {
1721
1722 // ++ >>>>>>>>>> MAGIC <<<<<<<<<<
1723
1724 // convert to Curvilinear distance over the parabolic dish
1725
1726 sx = Lin2Curv( xcut[0] );
1727 sy = Lin2Curv( xcut[1] );
1728
1729 // is it outside the dish?
1730
1731 if ((fabs(sx) > 850.0) || (fabs(sy) > 850.0)) {
1732 //cout << "CONDICION 1 !" << endl << flush;
1733 //cout << '1';
1734 continue;
1735 }
1736
1737 // -- >>>>>>>>>> MAGIC <<<<<<<<<<
1738
1739 } // endif (ct_Type == 1)
1740
1741 // calculate the mirror to be used
1742
1743 distmirr = 1000000.;
1744
1745 for (i=0; i<ct_NMirrors && distmirr>=ct_RMirror; ++i) {
1746 distmirr2 = sqrt(SQR(ct_data[i][CT_X] - xcut[0]) +
1747 SQR(ct_data[i][CT_Y] - xcut[1]) +
1748 SQR(ct_data[i][CT_Z] - xcut[2]));
1749 if (distmirr2 < distmirr) {
1750 i_mirror = i;
1751 distmirr = distmirr2;
1752 }
1753 }
1754
1755 // the mirror to use is i_mirror (calculated several lines above)
1756 // check whether the photon is outside the nearest (this) mirror
1757
1758 if ( ct_Type == 0 ) {
1759
1760 // ++ >>>>>>>>>> CT1 <<<<<<<<<<
1761 if (distmirr > ct_RMirror) {
1762 //cout << "CONDICION 2 !" << endl << flush;
1763 //cout << '2';
1764 continue;
1765 }
1766 // -- >>>>>>>>>> CT1 <<<<<<<<<<
1767
1768 } else {
1769
1770 // ++ >>>>>>>>>> MAGIC <<<<<<<<<<
1771 if ((fabs(ct_data[i_mirror][CT_SX] - sx) > ct_RMirror) ||
1772 (fabs(ct_data[i_mirror][CT_SY] - sy) > ct_RMirror)) {
1773 //cout << "CONDICION 2 !" << endl << flush;
1774 //cout << '2';
1775 continue;
1776 }
1777 // -- >>>>>>>>>> MAGIC <<<<<<<<<<
1778
1779 }
1780
1781 // calculate matrices for the mirror
1782
1783 makeOmega (-ct_data[i_mirror][CT_THETAN],
1784 ct_data[i_mirror][CT_PHIN]);
1785 makeOmegaI(-ct_data[i_mirror][CT_THETAN],
1786 ct_data[i_mirror][CT_PHIN]);
1787
1788 // change to the system of the mirror
1789
1790 // first translation...
1791
1792 xmm[0] = xCT[0] - ct_data[i_mirror][CT_X];
1793 xmm[1] = xCT[1] - ct_data[i_mirror][CT_Y];
1794 xmm[2] = xCT[2] - ct_data[i_mirror][CT_Z];
1795
1796 // ...then rotation
1797
1798 applyMxV( Omega, xmm, xm );
1799 applyMxV( Omega, rCT, rm );
1800
1801 // the vector rCT should be normalized, and
1802 // so the vector rm remains normalized as well, but, anyhow...
1803
1804 rnorm = NORM( rm );
1805 rm[0] /= rnorm;
1806 rm[1] /= rnorm;
1807 rm[2] /= rnorm;
1808
1809 // calculate the intersection of the trayectory of the photon
1810 // with the mirror
1811 // we reproduce the calculation of the coefficients of the
1812 // second order polynomial in z (=xm[2]), made with
1813 // Mathematica
1814
1815 /*
1816 * In[1]:= esfera:=x^2+y^2+(z-R)^2-R^2;
1817 * recta:={x->x0+u/w(z-z0),y->y0+v/w(z-z0)}
1818 *
1819 * In[2]:= esfera
1820 *
1821 * 2 2 2 2
1822 * Out[2]= -R + x + y + (-R + z)
1823 *
1824 * In[3]:= recta
1825 *
1826 * u (z - z0) v (z - z0)
1827 * Out[3]= {x -> x0 + ----------, y -> y0 + ----------}
1828 * w w
1829 *
1830 * In[4]:= esf=esfera /. recta
1831 *
1832 * 2 2 u (z - z0) 2 v (z - z0) 2
1833 * Out[4]= -R + (-R + z) + (x0 + ----------) + (y0 + ----------)
1834 * w w
1835 *
1836 * In[5]:= coefs=CoefficientList[ExpandAll[esf],z]
1837 *
1838 * 2 2 2 2
1839 * 2 2 2 u x0 z0 2 v y0 z0 u z0 v z0
1840 * Out[5]= {x0 + y0 - --------- - --------- + ------ + ------,
1841 * w w 2 2
1842 * w w
1843 *
1844 * 2 2 2 2
1845 * 2 u x0 2 v y0 2 u z0 2 v z0 u v
1846 * > -2 R + ------ + ------ - ------- - -------, 1 + -- + --}
1847 * w w 2 2 2 2
1848 * w w w w
1849 * In[6]:= Simplify[ExpandAll[coefs*w^2]]
1850 *
1851 * 2 2 2 2 2 2
1852 * Out[6]= {w (x0 + y0 ) - 2 w (u x0 + v y0) z0 + (u + v ) z0 ,
1853 *
1854 * 2 2 2 2 2
1855 * > -2 (R w - u w x0 + u z0 + v (-(w y0) + v z0)), u + v + w }
1856 *
1857 */
1858
1859 // the z coordinate is calculated, using the coefficients
1860 // shown above
1861
1862 a = SQR(rm[0]) + SQR(rm[1]) + SQR(rm[2]);
1863 if ( ct_Type == 0 ) {
1864 // CT1
1865 b = -2*(2.*ct_Focal_mean*SQR(rm[2])
1866 - rm[0]*rm[2]*xm[0]
1867 + SQR(rm[0])*xm[2]
1868 + rm[1]*(-(rm[2]*xm[1]) + rm[1]*xm[2]));
1869 } else {
1870 // MAGIC
1871 b = -2*(2.*ct_data[i_mirror][CT_FOCAL]*SQR(rm[2])
1872 - rm[0]*rm[2]*xm[0]
1873 + SQR(rm[0])*xm[2]
1874 + rm[1]*(-(rm[2]*xm[1]) + rm[1]*xm[2]));
1875 }
1876 c = (SQR(rm[2])*(SQR(xm[0]) + SQR(xm[1]))
1877 - 2*rm[2]*(rm[0]*xm[0] + rm[1]*xm[1])*xm[2] +
1878 (SQR(rm[0]) + SQR(rm[1]))*SQR(xm[2]));
1879
1880 d = sqrt( b*b - 4.0*a*c );
1881
1882 // two possible values for z
1883
1884 t1 = (-b+d) / (2.0*a);
1885 t2 = (-b-d) / (2.0*a);
1886
1887 // z must be the minimum of t1 and t2
1888
1889 xcut[2] = (t1 < t2) ? t1 : t2;
1890 xcut[0] = xm[0] + rm[0]/rm[2]*(xcut[2]-xm[2]);
1891 xcut[1] = xm[1] + rm[1]/rm[2]*(xcut[2]-xm[2]);
1892
1893 //++
1894 // BLACK SPOTS: If the photon hits the black spot, it's lost
1895 //--
1896
1897 if ( sqrt(SQR(xcut[0]) + SQR(xcut[1])) < ct_BlackSpot_rad ) {
1898 //cout << "CONDITION 3!\n" << flush;
1899 //cout << '3';
1900 continue;
1901 }
1902
1903 // if still we have the photon, continue with the reflexion
1904
1905 // calculate normal vector in this point
1906 // (and normalize, with the sign changed)
1907
1908 rnor[0] = 2.0*xcut[0];
1909 rnor[1] = 2.0*xcut[1];
1910 rnor[2] = 2.0*(xcut[2] - 2.0*ct_Focal[i_mirror]);
1911
1912 rnorm = -NORM( rnor );
1913 rnor[0] /= rnorm;
1914 rnor[1] /= rnorm;
1915 rnor[2] /= rnorm;
1916
1917 // now, both "normal" vector and original trayectory are
1918 // normalized
1919 // just project the original vector in the normal, and
1920 // take it as the "mean" position of the original and
1921 // the "reflected" vector
1922 // from this, we can calculate the "reflected" vector
1923 // calpha = cos(angle(rnor,rm))
1924
1925 calpha = fabs(rnor[0]*rm[0] + rnor[1]*rm[1] + rnor[2]*rm[2]);
1926
1927 // finally!!! we have the reflected trayectory of the photon
1928
1929 rrefl[0] = 2.0*rnor[0]*calpha - rm[0];
1930 rrefl[1] = 2.0*rnor[1]*calpha - rm[1];
1931 rrefl[2] = 2.0*rnor[2]*calpha - rm[2];
1932
1933 rnorm = NORM( rrefl );
1934 rrefl[0] /= rnorm;
1935 rrefl[1] /= rnorm;
1936 rrefl[2] /= rnorm;
1937
1938 // let's go back to the coordinate system of the CT
1939
1940 // first rotation...
1941 applyMxV( OmegaI, xcut, xcutCT);
1942 applyMxV( OmegaI, rrefl, rreflCT);
1943
1944 // ...then translation
1945 xcutCT[0] += ct_data[i_mirror][CT_X];
1946 xcutCT[1] += ct_data[i_mirror][CT_Y];
1947 xcutCT[2] += ct_data[i_mirror][CT_Z];
1948
1949 // calculate intersection of this trayectory and the camera plane
1950 // in the system of the CT, this plane is z = ct_Focal
1951
1952 t = (ct_Focal_mean - xcutCT[2]) / rreflCT[2];
1953
1954 xcam[0] = xcutCT[0] + rreflCT[0]*t;
1955 xcam[1] = xcutCT[1] + rreflCT[1]*t;
1956 xcam[2] = xcutCT[2] + rreflCT[2]*t;
1957
1958 //++
1959 // AXIS DEVIATION: We introduce it here just as a first order
1960 // correction, by modifying the position of the reflected photon.
1961 //--
1962
1963 xcam[0] += AxisDeviation[0][i_mirror];
1964 xcam[1] += AxisDeviation[1][i_mirror];
1965
1966 //++
1967 // SMEARING: We apply the point spread function for the mirrors
1968 //--
1969
1970 // get two N(0;1) random numbers
1971
1972 rnormal( NormalRandomNumbers, 2 );
1973
1974 // modify the Cphoton position in the camera
1975
1976 xcam[0] += NormalRandomNumbers[0] * ct_PSpread_mean;
1977 xcam[1] += NormalRandomNumbers[1] * ct_PSpread_mean;
1978
1979 // check whether the photon goes out of the camera
1980
1981 if ( (SQR(xcam[0])+SQR(xcam[1])) > SQR(ct_CameraWidth) ) {
1982 continue;
1983 }
1984
1985 //++
1986 // ANGLE OF INCIDENCE
1987 //--
1988
1989 // calculate angle of incidence between tray. and camera plane
1990 // the camera plane is
1991 // 0 y + 0 y + z - ct_Focal = 0 => (A,B,C,D) = (0,0,1,-ct_Focal)
1992 // from Table 3.20 "Tasch. der Math."
1993
1994 phi = asin(rreflCT[2]);
1995
1996 //++
1997 // TIMING
1998 //--
1999
2000 // calculate the new time of the photon (in the camera)
2001
2002 // initial time since first interaction
2003
2004 t = photon.get_t();
2005
2006 // substract path from the mirror till the ground, 'cos
2007 // the photon actually hit the mirror!!
2008
2009
2010 t = t + ((( xm[2] > 0. ) ? -1.0 : +1.0) *
2011 sqrt( SQR(xm[0] - xcut[0]) +
2012 SQR(xm[1] - xcut[1]) +
2013 SQR(xm[2] - xcut[2]) ) / Speed_of_Light_air_cmns);
2014
2015 // add path from the mirror till the camera
2016
2017 t = t + sqrt( SQR(xcutCT[0] - xcam[0]) +
2018 SQR(xcutCT[1] - xcam[1]) +
2019 SQR(xcutCT[2] - xcam[2]) ) / Speed_of_Light_air_cmns;
2020
2021 // save data in the photon variable
2022
2023 cphoton.fill( photon.get_id(),
2024 xcam[0], xcam[1],
2025 photon.get_u(), photon.get_v(),
2026 photon.get_h(), t, phi);
2027
2028 // show it
2029
2030 if ( verbose >= VERBOSE_MAXIMAL ) {
2031 cout << '@' << '2'
2032 << ' ' << ncer
2033 << ' ' << xCT[0]
2034 << ' ' << xCT[1]
2035 << ' ' << xCT[2]
2036 << ' ' << rCT[0]
2037 << ' ' << rCT[1]
2038 << ' ' << rCT[2]
2039 << ' ' << xcut[0]
2040 << ' ' << xcut[1]
2041 << ' ' << xcut[2]
2042 << endl;
2043 cout << '@' << '3'
2044 << ' ' << ncer
2045 << ' ' << sx
2046 << ' ' << sy
2047 << ' ' << i_mirror
2048 << ' ' << ct_data[i_mirror][CT_SX]
2049 << ' ' << ct_data[i_mirror][CT_SY]
2050 << ' ' << ct_data[i_mirror][CT_SX] - sx
2051 << ' ' << ct_data[i_mirror][CT_SY] - sy
2052 << endl;
2053 if ( pb_ParallelBeam == FALSE ) {
2054 cout << '@' << '4'
2055 << ' ' << ncer
2056 << ' ' << xcam[0]
2057 << ' ' << xcam[1]
2058 << ' ' << xcam[2]
2059 << ' ' << phi
2060 << ' ' << photon.get_t()-stat.get_tfirst()
2061 << ' ' << t-stat.get_tfirst()
2062 << endl << flush;
2063 } else {
2064 cout << '@' << '4'
2065 << ' ' << ncer
2066 << ' ' << xcam[0]
2067 << ' ' << xcam[1]
2068 << ' ' << xcam[2]
2069 << ' ' << phi
2070 << endl << flush;
2071 }
2072 }
2073
2074 // write the photon data to the output file
2075
2076 if ( Data_To_STDOUT )
2077 cout.write ( (char *)&cphoton, cphoton.mysize() );
2078 else
2079 cphoton.write( outputfile );
2080
2081 // increase number of Cphotons written
2082
2083 ++nCphotons;
2084
2085 } // while still there are photons left
2086
2087 if ( verbose >= VERBOSE_NORMAL )
2088 cout << nCphotons << ' '
2089 << nbeforeTR << ' '
2090 << nafterTR << ' '
2091 << endl << flush;
2092
2093 // write end-of-event mark
2094
2095 if ( Data_To_STDOUT )
2096 cout.write( FLAG_END_OF_EVENT, SIZE_OF_FLAGS );
2097 else
2098 outputfile.write( FLAG_END_OF_EVENT, SIZE_OF_FLAGS );
2099
2100 // close files
2101
2102 cerfile.close();
2103 stat.closefile();
2104
2105 // show how many photons were written
2106
2107 if ( verbose >= VERBOSE_MAXIMAL )
2108 log( SIGNATURE, "%d C-photons written.\n", nCphotons );
2109
2110 // if we are reading data from STDIN, delete the last used files
2111
2112 if ( Data_From_STDIN ) {
2113
2114 unlink( cername );
2115 unlink( staname );
2116
2117 }
2118
2119 } // while there are still data left
2120
2121 // write end-of-run mark
2122
2123 if ( Data_To_STDOUT )
2124 cout.write( FLAG_END_OF_RUN, SIZE_OF_FLAGS );
2125 else
2126 outputfile.write( FLAG_END_OF_RUN, SIZE_OF_FLAGS );
2127
2128 // close directory
2129
2130 closedir( directory );
2131
2132 } // for each data directory
2133
2134 // write end-of-file mark
2135
2136 if ( Data_To_STDOUT )
2137 cout.write( FLAG_END_OF_FILE, SIZE_OF_FLAGS );
2138 else
2139 outputfile.write( FLAG_END_OF_FILE, SIZE_OF_FLAGS );
2140
2141 // close output file
2142
2143 if ( verbose >= VERBOSE_MINIMAL )
2144 log( SIGNATURE, "Closing output file %s\n", outname );
2145
2146 if ( ! Data_To_STDOUT )
2147 outputfile.close();
2148
2149 // clean everything
2150
2151 if ( verbose >= VERBOSE_MINIMAL )
2152 log( SIGNATURE, "Cleanning . . .\n" );
2153
2154 clean();
2155
2156 // program finished
2157
2158 if ( verbose >= VERBOSE_MINIMAL )
2159 log( SIGNATURE, "Done.\n");
2160
2161 return ( 0 );
2162}
2163//!@}
2164
2165// @T \newpage
2166
2167//=---------------------------------------------------------------------
2168// @subsection Functions definition
2169
2170//!---------------------------------------------------------------------
2171// @name present
2172//
2173// @desc Make some presentation
2174//
2175// @date Sat Jun 27 05:58:56 MET DST 1998
2176//----------------------------------------------------------------------
2177// @function
2178
2179//!@{
2180void
2181present(void)
2182{
2183 cout << "##################################################\n"
2184 << SIGNATURE << '\n' << '\n'
2185 << "Simulation of the reflector\n"
2186 << "J C Gonzalez, May 1998\n"
2187 << "##################################################\n\n";
2188}
2189//!@}
2190
2191
2192//!---------------------------------------------------------------------
2193// @name usage
2194//
2195// @desc show help
2196//
2197// @date Tue Dec 15 16:23:30 MET 1998
2198//----------------------------------------------------------------------
2199// @function
2200
2201//!@{
2202void
2203usage(void)
2204{
2205 present();
2206 cout << "\nusage ::\n\n"
2207 << "\t reflector "
2208 << " [ -@ paramfile ] "
2209 << " [ -h ] "
2210 << "\n\n or \n\n"
2211 << "\t reflector < paramfile"
2212 << "\n\n";
2213 exit(0);
2214}
2215//!@}
2216
2217
2218//!---------------------------------------------------------------------
2219// @name clean
2220//
2221// @desc basically, frees memory
2222//
2223// @date Sat Jun 27 05:58:56 MET DST 1998
2224//----------------------------------------------------------------------
2225// @function
2226
2227//!@{
2228void
2229clean(void)
2230{
2231 int i;
2232
2233 // delete focals table
2234
2235 delete [] ct_Focal;
2236 ct_Focal = 0;
2237
2238 // delete reflectivity table
2239
2240 for (i=0; i<2; i++) {
2241 delete [] Reflectivity[i];
2242 }
2243
2244 delete [] Reflectivity;
2245 Reflectivity = 0;
2246
2247 // delete mirrors' data table
2248
2249 for (i=0; i<ct_NMirrors; i++) {
2250 delete [] ct_data[i];
2251 }
2252
2253 delete [] ct_data;
2254 ct_data = 0;
2255}
2256//!@}
2257
2258
2259//!---------------------------------------------------------------------
2260// @name log
2261//
2262// @desc function to send log information
2263//
2264// @var funct Name of the caller function
2265// @var fmt Format to be used (message)
2266// @var ... Other information to be shown
2267//
2268// @date Sat Jun 27 05:58:56 MET DST 1998
2269//----------------------------------------------------------------------
2270// @function
2271
2272//!@{
2273void
2274log(const char *funct, char *fmt, ...)
2275{
2276 va_list args;
2277
2278 // Display the name of the function that called error
2279 printf("[%s]: ", funct);
2280
2281 // Display the remainder of the message
2282 va_start(args, fmt);
2283 vprintf(fmt, args);
2284 va_end(args);
2285}
2286//!@}
2287
2288
2289//!---------------------------------------------------------------------
2290// @name logerr
2291//
2292// @desc function to send log information to stderr
2293//
2294// @var funct Name of the caller function
2295// @var fmt Format to be used (message)
2296// @var ... Other information to be shown
2297//
2298// @date Sat Jun 27 05:58:56 MET DST 1998
2299//----------------------------------------------------------------------
2300// @function
2301
2302//!@{
2303void
2304logerr(const char *funct, char *fmt, ...)
2305{
2306 va_list args;
2307
2308 // Display the name of the function that called error
2309 printf("[%s]: ", funct);
2310
2311 // Display the remainder of the message
2312 va_start(args, fmt);
2313 vfprintf(stderr, fmt, args);
2314 va_end(args);
2315}
2316//!@}
2317
2318
2319//!---------------------------------------------------------------------
2320// @name error
2321//
2322// @desc function to send an error message, and abort the program
2323//
2324// @var funct Name of the caller function
2325// @var fmt Format to be used (message)
2326// @var ... Other information to be shown
2327//
2328// @date Sat Jun 27 05:58:56 MET DST 1998
2329//----------------------------------------------------------------------
2330// @function
2331
2332//!@{
2333void
2334error(const char *funct, char *fmt, ...)
2335{
2336 va_list args;
2337
2338 // Display the name of the function that called error
2339 fprintf(stderr, "ERROR in %s: ", funct);
2340
2341 // Display the remainder of the message
2342 va_start(args, fmt);
2343 vfprintf(stderr, fmt, args);
2344 va_end(args);
2345
2346 abort();
2347}
2348//!@}
2349
2350
2351//!---------------------------------------------------------------------
2352// @name makeOmega
2353//
2354// @desc function to calculate the matrix Omega(theta,phi)
2355//
2356// @var theta Angle theta of the transformation
2357// @var phi Angle phi of the transformation
2358//
2359// @date Sat Jun 27 05:58:56 MET DST 1998
2360//----------------------------------------------------------------------
2361// @function
2362
2363//!@{
2364void
2365makeOmega (float theta, float phi)
2366{
2367 static float ct, st, cp, sp;
2368
2369 // shortcuts for cosine and sine of theta and phi
2370 ct = cos(theta);
2371 st = sin(theta);
2372 cp = cos(phi);
2373 sp = sin(phi);
2374
2375 // save values in the array (see top of file)
2376 Omega[0][0] = cp*ct;
2377 Omega[0][1] = sp*ct;
2378 Omega[0][2] = -st;
2379
2380 Omega[1][0] = -sp;
2381 Omega[1][1] = cp;
2382 Omega[1][2] = 0;
2383
2384 Omega[2][0] = cp*st;
2385 Omega[2][1] = sp*st;
2386 Omega[2][2] = ct;
2387}
2388//!@}
2389
2390
2391//!---------------------------------------------------------------------
2392// @name makeOmegaI
2393//
2394// @desc function to calculate the matrix Omega-1(theta,phi)
2395//
2396// @var theta Angle theta of the transformation
2397// @var phi Angle phi of the transformation
2398//
2399// @date Sat Jun 27 05:58:56 MET DST 1998
2400//----------------------------------------------------------------------
2401// @function
2402
2403//!@{
2404void
2405makeOmegaI(float theta, float phi)
2406{
2407 static float ct, st, cp, sp;
2408
2409 // shortcuts for cosine and sine of theta and phi
2410 ct = cos(theta);
2411 st = sin(theta);
2412 cp = cos(phi);
2413 sp = sin(phi);
2414
2415 // save values in the array (see top of file)
2416 OmegaI[0][0] = cp*ct;
2417 OmegaI[0][1] = -sp;
2418 OmegaI[0][2] = cp*st;
2419
2420 OmegaI[1][0] = sp*ct;
2421 OmegaI[1][1] = cp;
2422 OmegaI[1][2] = sp*st;
2423
2424 OmegaI[2][0] = -st;
2425 OmegaI[2][1] = 0;
2426 OmegaI[2][2] = ct;
2427}
2428//!@}
2429
2430
2431//!---------------------------------------------------------------------
2432// @name applyMxv
2433//
2434// @desc returns the vector v' such that v' = M x v
2435//
2436// @var M matrix of the transformation
2437// @var v vector to be multiplied
2438// @var vi resulting vector
2439//
2440// @date Sat Jun 27 05:58:56 MET DST 1998
2441//----------------------------------------------------------------------
2442// @function
2443
2444//!@{
2445void
2446applyMxV(float M[3][3], float *V, float *Vp)
2447{
2448 Vp[0] = (M[0][0] * V[0] +
2449 M[0][1] * V[1] +
2450 M[0][2] * V[2]);
2451 Vp[1] = (M[1][0] * V[0] +
2452 M[1][1] * V[1] +
2453 M[1][2] * V[2]);
2454 Vp[2] = (M[2][0] * V[0] +
2455 M[2][1] * V[1] +
2456 M[2][2] * V[2]);
2457}
2458//!@}
2459
2460
2461//!---------------------------------------------------------------------
2462// @name read_ct_file
2463//
2464// @desc read CT definition file
2465//
2466// @date Sat Jun 27 05:58:56 MET DST 1998
2467//----------------------------------------------------------------------
2468// @function
2469
2470//!@{
2471void
2472read_ct_file(void)
2473{
2474 char line[LINE_MAX_LENGTH]; // line to get from the ctin
2475 char token[ITEM_MAX_LENGTH]; // a single token
2476 char fmirr[40]; // name of the binary version of mirrors data
2477 int i; // dummy counters
2478
2479 if ( verbose >= VERBOSE_MINIMAL )
2480 log( "read_ct_file", "start.\n" );
2481
2482 ifstream ctin ( get_ct_filename() );
2483
2484 if ( ctin.bad() )
2485 error( "read_ct_file",
2486 "Cannot open CT def. file: %s\n", get_ct_filename() );
2487
2488 // loop till the "end" directive is reached
2489
2490 while (!ctin.eof()) {
2491
2492 // get line from stdin
2493
2494 ctin.getline(line, LINE_MAX_LENGTH);
2495
2496 // look for each item at the beginning of the line
2497
2498 for (i=0; i<=define_mirrors; i++)
2499 if (strstr(line, CT_ITEM_NAMES[i]) == line)
2500 break;
2501
2502 // if it is not a valid line, just ignore it
2503
2504 if (i == define_mirrors+1)
2505 continue;
2506
2507 // case block for each directive
2508
2509 switch ( i ) {
2510
2511 case type: // <type of telescope> (0:CT1 ¦ 1:MAGIC)
2512
2513 // get focal distance
2514
2515 sscanf(line, "%s %d", token, &ct_Type);
2516
2517 if ( verbose >= VERBOSE_MINIMAL )
2518 log( "read_ct_file", "<Type of Telescope>: %s\n",
2519 ((ct_Type==0) ? "CT1" : "MAGIC") );
2520
2521 break;
2522
2523 case focal_distance: // <focal distance> [cm]
2524
2525 // get focal distance
2526
2527 sscanf(line, "%s %f", token, &ct_Focal_mean);
2528
2529 if ( verbose >= VERBOSE_MINIMAL )
2530 log( "read_ct_file", "<Focal distance>: %f cm\n", ct_Focal_mean );
2531
2532 break;
2533
2534 case focal_std: // s(focal distance) [cm]
2535
2536 // get focal distance
2537
2538 sscanf(line, "%s %f", token, &ct_Focal_std);
2539
2540 if ( verbose >= VERBOSE_MINIMAL )
2541 log( "read_ct_file", "s(Focal distance): %f cm\n", ct_Focal_std );
2542
2543 break;
2544
2545 case point_spread: // <point spread> [cm]
2546
2547 // get point spread
2548
2549 sscanf(line, "%s %f", token, &ct_PSpread_mean);
2550
2551 if ( verbose >= VERBOSE_MINIMAL )
2552 log( "read_ct_file", "<Point spread>: %f cm\n", ct_PSpread_mean );
2553
2554 break;
2555
2556 case point_std: // s(point spread) [cm]
2557
2558 // get point spread
2559
2560 sscanf(line, "%s %f", token, &ct_PSpread_std);
2561
2562 if ( verbose >= VERBOSE_MINIMAL )
2563 log( "read_ct_file", "s(Point spread): %f cm\n", ct_PSpread_std );
2564
2565 break;
2566
2567 case adjustment_dev: // s(adjustment_dev) [cm]
2568
2569 // get point spread
2570
2571 sscanf(line, "%s %f", token, &ct_Adjustment_std);
2572
2573 if ( verbose >= VERBOSE_MINIMAL )
2574 log( "read_ct_file", "s(Adjustment): %f cm\n", ct_Adjustment_std );
2575
2576 break;
2577
2578 case black_spot: // radius of the black spot in the center [cm]
2579
2580 // get black spot radius
2581
2582 sscanf(line, "%s %f", token, &ct_BlackSpot_rad);
2583
2584 if ( verbose >= VERBOSE_MINIMAL )
2585 log( "read_ct_file", "Radius of the black spots: %f cm\n",
2586 ct_BlackSpot_rad);
2587
2588 break;
2589
2590 case r_mirror: // radius of the mirrors [cm]
2591
2592 // get radius of mirror
2593
2594 sscanf(line, "%s %f", token, &ct_RMirror);
2595
2596 if ( verbose >= VERBOSE_MINIMAL )
2597 log( "read_ct_file", "Radii of the mirrors: %f cm\n", ct_RMirror );
2598
2599 break;
2600
2601 case n_mirrors: // number of mirrors
2602
2603 // get the name of the output_file from the line
2604
2605 sscanf(line, "%s %d", token, &ct_NMirrors);
2606
2607 if ( verbose >= VERBOSE_MINIMAL )
2608 log( "read_ct_file", "Number of mirrors: %d\n", ct_NMirrors );
2609
2610 break;
2611
2612 case camera_width: // camera width [cm]
2613
2614 // get the name of the ct_file from the line
2615
2616 sscanf(line, "%s %f", token, &ct_CameraWidth);
2617
2618 if ( verbose >= VERBOSE_MINIMAL )
2619 log( "read_ct_file", "Camera width: %f cm\n", ct_CameraWidth );
2620
2621 break;
2622
2623 case n_pixels: // number of pixels
2624
2625 // get the name of the output_file from the line
2626
2627 sscanf(line, "%s %d", token, &ct_NPixels);
2628
2629 if ( verbose >= VERBOSE_MINIMAL )
2630 log( "read_ct_file", "Number of pixels: %d\n", ct_NPixels );
2631
2632 break;
2633
2634 case pixel_width: // pixel width [cm]
2635
2636 // get the name of the ct_file from the line
2637
2638 sscanf(line, "%s %f", token, &ct_PixelWidth);
2639
2640 if ( verbose >= VERBOSE_MINIMAL )
2641 log( "read_ct_file", "Pixel width: %f cm\n", ct_PixelWidth );
2642
2643 break;
2644
2645 case define_mirrors: // read table with the parameters of the mirrors
2646
2647 if ( verbose >= VERBOSE_MINIMAL )
2648 log( "read_ct_file", "Table of mirrors data:\n" );
2649
2650 // check whether the number of mirrors was already set
2651
2652 if ( ct_NMirrors == 0 )
2653 error( "read_ct_file", "NMirrors was not set.\n" );
2654
2655 // allocate memory for paths list
2656
2657 if ( verbose >= VERBOSE_MINIMAL )
2658 log( "read_ct_file", "Allocating memory for ct_data\n" );
2659
2660 ct_data = new float*[ct_NMirrors];
2661
2662 for (i=0; i<ct_NMirrors; i++)
2663 ct_data[i] = new float[CT_NDATA];
2664
2665 // read data
2666
2667 // look for binary version of mirror data (faster to load)
2668
2669 sprintf( fmirr, "%s.mirr", get_ct_filename() );
2670
2671 ifstream ctmirr;
2672 ctmirr.open( fmirr );
2673
2674 if ( ctmirr.bad() ) {
2675
2676 // the file does not exist
2677
2678 ctmirr.close(); // close not existing file
2679
2680 ofstream ctmirr_write ( fmirr ); // open file to save the data
2681
2682 if ( verbose >= VERBOSE_MINIMAL )
2683 log( "read_ct_file", "Reading mirrors data...\n" );
2684
2685 // read ASCII data
2686
2687 for (i=0; i<ct_NMirrors; i++) {
2688 // for (j=0; j<CT_NDATA; j++)
2689
2690 ctin.getline(line, LINE_MAX_LENGTH);
2691 sscanf(line, "%f %f %f %f %f %f %f %f %f %f %f %f",
2692 &ct_data[i][0],
2693 &ct_data[i][1],
2694 &ct_data[i][2],
2695 &ct_data[i][3],
2696 &ct_data[i][4],
2697 &ct_data[i][5],
2698 &ct_data[i][6],
2699 &ct_data[i][7],
2700 &ct_data[i][8],
2701 &ct_data[i][9],
2702 &ct_data[i][10],
2703 &ct_data[i][11]);
2704
2705 cout << '[' << i << ']' << flush;
2706 }
2707 cout << endl << flush;
2708
2709 // save binary data and close output file
2710
2711 for (i=0; i<ct_NMirrors; i++)
2712 ctmirr_write.write( (char*)ct_data[i], CT_NDATA*sizeof(float) );
2713
2714 ctmirr_write.close();
2715
2716 } else {
2717
2718 // Excellent!
2719 // read binary data and close input file
2720
2721 for (i=0; i<ct_NMirrors; i++)
2722 ctmirr.read( (char*)ct_data[i], CT_NDATA*sizeof(float) );
2723
2724 ctmirr.close();
2725
2726 }
2727
2728 break;
2729
2730 } // switch ( i )
2731
2732 } // while (! is_end)
2733
2734 // read reflectivity file
2735
2736 read_reflectivity();
2737
2738 // read axis deviations file
2739
2740 read_axisdev();
2741
2742 // read focals file
2743
2744 read_focals();
2745
2746 // end
2747
2748 if ( verbose >= VERBOSE_MINIMAL )
2749 log( "read_ct_file", "done.\n" );
2750 return;
2751}
2752//!@}
2753
2754
2755//!---------------------------------------------------------------------
2756// @name dist_r_P
2757//
2758// @desc distance straight line r - point P
2759//
2760// @var a coord. X of a fixed point of the straight line
2761// @var b coord. Y of a fixed point of the straight line
2762// @var c coord. Z of a fixed point of the straight line
2763// @var l component X of a vector of the straight line
2764// @var m component Y of a vector of the straight line
2765// @var n component Z of a vector of the straight line
2766// @var x coord. X of the point P
2767// @var y coord. Y of the point P
2768// @var z coord. Z of the point P
2769//
2770// @return Returns the distance d(r,P) between the line and the point P
2771//
2772// @date Sat Jun 27 05:58:56 MET DST 1998
2773//----------------------------------------------------------------------
2774// @function
2775
2776//!@{
2777float
2778dist_r_P(float a, float b, float c,
2779 float l, float m, float n,
2780 float x, float y, float z)
2781{
2782 return (
2783 sqrt((SQR((a-x)*m-(b-y)*l) +
2784 SQR((b-y)*n-(c-z)*m) +
2785 SQR((c-z)*l-(a-x)*n))/
2786 (SQR(l)+SQR(m)+SQR(n))
2787 )
2788 );
2789}
2790//!@}
2791
2792
2793//!---------------------------------------------------------------------
2794// @name read_reflectivity
2795//
2796// @desc read reflectivity data for the mirrors
2797//
2798// @date Sat Jun 27 05:58:56 MET DST 1998
2799//----------------------------------------------------------------------
2800// @function
2801
2802//!@{
2803void
2804read_reflectivity(void)
2805{
2806 ifstream reflfile;
2807 ofstream reflfile_write;
2808 int i, errors=0;
2809 char line[LINE_MAX_LENGTH];
2810
2811 // try to open the file
2812 reflfile.open( REFLECTIVITY_FILE );
2813
2814 // if it is wrong or does not exist:
2815
2816 while ( reflfile.bad() ) {
2817
2818 ++errors;
2819
2820 if ( verbose >= VERBOSE_MINIMAL )
2821 log( "read_reflectivity", "Cannot open reflectivity file: %s\n",
2822 REFLECTIVITY_FILE );
2823
2824 if ( errors > 5 )
2825 error( "read_reflectivity", "Exiting.");
2826
2827 // try to re-generate the file
2828
2829 if ( verbose >= VERBOSE_MINIMAL )
2830 log( "read_reflectivity", "Generating file: %s\n",
2831 REFLECTIVITY_FILE );
2832
2833 // try to open the file
2834 reflfile_write.open( REFLECTIVITY_FILE );
2835
2836 // write short comment in the beginning
2837
2838 reflfile_write << "# reflectivity file" << endl
2839 << "# reflectivity for each mirror in the frame" << endl
2840 << "# J C Gonzalez, 1998" << endl
2841 << "#" << endl;
2842
2843 // write data
2844
2845 // write number of data-points
2846
2847 reflfile_write << "# number of datapoints" << endl
2848 << 35 << endl;
2849
2850 // write data-points
2851
2852 reflfile_write << "# datapoints" << endl;
2853
2854 for ( i=0; i<35; ++i )
2855 reflfile_write << 270.0+i*10.0 << ' ' << .9 << endl;
2856
2857 // close it and try to re-open
2858
2859 reflfile_write.close();
2860 reflfile.close();
2861 reflfile.open( REFLECTIVITY_FILE );
2862
2863 }
2864
2865 if ( verbose >= VERBOSE_MINIMAL ) {
2866 if ( errors > 0 )
2867 log( "read_reflectivity", "Reading reflectivity file after %d fail%c\n",
2868 errors, ((errors>1)?'s':' ') );
2869 else
2870 log( "read_reflectivity", "Reading reflectivity file . . .\n" );
2871 }
2872
2873 // we set this to -1, in order to read in the first place the
2874 // number of data-points for the Reflectivity curve
2875 i = -1;
2876
2877 // scan the whole file
2878
2879 while ( ! reflfile.eof() ) {
2880
2881 // get line from the file
2882
2883 reflfile.getline(line, LINE_MAX_LENGTH);
2884
2885 // skip if comment
2886
2887 if ( *line == '#' )
2888 continue;
2889
2890 // get the value
2891
2892 if ( i < 0 ) { // if it is the first number we read
2893
2894 // read the number of datapoints
2895
2896 sscanf(line, "%d", &nReflectivity);
2897
2898 // allocate memory for Reflectivity table
2899
2900 Reflectivity = new float * [2];
2901 Reflectivity[0] = new float[nReflectivity];
2902 Reflectivity[1] = new float[nReflectivity];
2903
2904 } else { // if not, then it is a datapoint
2905
2906 // get the datapoint (wavelength, reflectivity)
2907
2908 sscanf(line, "%f %f", &Reflectivity[0][i], &Reflectivity[1][i]);
2909
2910 }
2911
2912 ++i;
2913
2914 }
2915
2916 // close it
2917
2918 reflfile.close();
2919
2920}
2921//!@}
2922
2923
2924//!---------------------------------------------------------------------
2925// @name read_axisdev
2926//
2927// @desc read axis deviations data for the mirrors
2928//
2929// @date Sat Jun 27 05:58:56 MET DST 1998
2930//----------------------------------------------------------------------
2931// @function
2932
2933//!@{
2934void
2935read_axisdev(void)
2936{
2937 ifstream axisdevfile;
2938 ofstream axisdevfile_write;
2939 int i, errors=0;
2940 char line[LINE_MAX_LENGTH];
2941
2942 // try to open the file
2943 axisdevfile.open( AXISDEVIATION_FILE );
2944
2945 // if it is wrong or does not exist:
2946
2947 while ( axisdevfile.bad() ) {
2948
2949 ++errors;
2950
2951 if ( verbose >= VERBOSE_MINIMAL )
2952 log( "read_axisdev", "Cannot open axis deviation file: %s\n",
2953 AXISDEVIATION_FILE );
2954
2955 if ( errors > 5 )
2956 error( "read_axisdev", "Exiting.");
2957
2958 // try to re-generate the file
2959
2960 if ( verbose >= VERBOSE_MINIMAL )
2961 log( "read_axisdev", "Generating file: %s\n",
2962 AXISDEVIATION_FILE );
2963
2964 // try to open the file
2965 axisdevfile_write.open( AXISDEVIATION_FILE );
2966
2967 // write short comment in the beginning
2968
2969 axisdevfile_write << "# axis deviation file" << endl
2970 << "# axis deviation for each mirror" << endl
2971 << "# J C Gonzalez, 1998" << endl
2972 << "#" << endl;
2973
2974 // write data
2975
2976 for ( i=0; i<ct_NMirrors; ++i )
2977 axisdevfile_write << gennor(0.0,ct_Adjustment_std)
2978 << ' '
2979 << gennor(0.0,ct_Adjustment_std)
2980 << endl;
2981
2982 // close it and try to re-open
2983
2984 axisdevfile_write.close();
2985 axisdevfile.close();
2986 axisdevfile.open( AXISDEVIATION_FILE );
2987
2988 }
2989
2990 if ( verbose >= VERBOSE_MINIMAL ) {
2991 if ( errors > 0 )
2992 log( "read_axisdev", "Reading axis deviation file after %d fail%c\n",
2993 errors, ((errors>1)?'s':' ') );
2994 else
2995 log( "read_axisdev", "Reading axis deviation file . . .\n" );
2996 }
2997
2998 //allocate memory for AxisDeviation table
2999
3000 AxisDeviation = new float * [2];
3001 AxisDeviation[0] = new float[ct_NMirrors];
3002 AxisDeviation[1] = new float[ct_NMirrors];
3003
3004 // scan the whole file
3005
3006 i = 0;
3007
3008 while ( (! axisdevfile.eof()) && (i<ct_NMirrors) ) {
3009
3010 // get line from the file
3011
3012 axisdevfile.getline(line, LINE_MAX_LENGTH);
3013
3014 // skip if comment
3015
3016 if ( *line == '#' )
3017 continue;
3018
3019 // get the value (dx, dy)
3020
3021 sscanf(line, "%f %f", &AxisDeviation[0][i], &AxisDeviation[1][i]);
3022
3023 ++i;
3024
3025 }
3026
3027 // did i do it for all the mirrors ?
3028
3029 if (i < ct_NMirrors) {
3030 error("read_axis_dev", "%s%s (%d < %d)\n",
3031 "Number of pairs in axisdev.dat file is",
3032 "\n\t\tsmaller than number of mirrors", i, ct_NMirrors);
3033 }
3034
3035 // close it
3036
3037 axisdevfile.close();
3038
3039}
3040//!@}
3041
3042
3043//!---------------------------------------------------------------------
3044// @name read_focals
3045//
3046// @desc read focals data for the mirrors
3047//
3048// @date Sat Jun 27 05:58:56 MET DST 1998
3049//----------------------------------------------------------------------
3050// @function
3051
3052//!@{
3053void
3054read_focals(void)
3055{
3056 ifstream focalfile;
3057 ofstream focalfile_write;
3058 int i, errors=0;
3059 char line[LINE_MAX_LENGTH];
3060
3061 /*!@'
3062
3063 Here we read the focals of each mirror, but only for CT1. For
3064 MAGIC we have already the Focals in the definition file. So, we
3065 copy these values to the |ct\_Focals| vector.
3066
3067 */
3068
3069 // >>>>>>>>>> MAGIC <<<<<<<<<<
3070 if ( ct_Type == 1) {
3071
3072 // allocate memory for the focals table
3073 ct_Focal = new float[ct_NMirrors];
3074
3075 for (i=0; i<ct_NMirrors; ++i)
3076 ct_Focal[i] = ct_data[i][CT_FOCAL];
3077
3078 return;
3079
3080 }
3081 // ! >>>>>>>>>> MAGIC <<<<<<<<<<
3082
3083 // try to open the file
3084 focalfile.open( FOCALS_FILE );
3085
3086 // if it is wrong or does not exist:
3087
3088 while ( focalfile.bad() ) {
3089
3090 ++errors;
3091
3092 if ( verbose >= VERBOSE_MINIMAL )
3093 log( "read_focals", "Cannot open focals file: %s\n",
3094 FOCALS_FILE );
3095
3096 if ( errors > 5 )
3097 error( "read_focals", "Exiting.");
3098
3099 // try to re-generate the file
3100
3101 if ( verbose >= VERBOSE_MINIMAL )
3102 log( "read_focals", "Generating file: %s\n",
3103 FOCALS_FILE );
3104
3105 // try to open the file
3106 focalfile_write.open( FOCALS_FILE );
3107
3108 // write short comment in the beginning
3109
3110 focalfile_write << "# focals file" << endl
3111 << "# focals for each mirror in the frame" << endl
3112 << "# J C Gonzalez, 1998" << endl;
3113
3114 // write data
3115
3116 rnormal( NormalRandomNumbers, ct_NMirrors+1 );
3117
3118 for (i=0; i<ct_NMirrors; i++)
3119 focalfile_write << NormalRandomNumbers[i] * ct_Focal_std + ct_Focal_mean
3120 << endl;
3121
3122 // close it and try to re-open
3123
3124 focalfile_write.close();
3125 focalfile.close();
3126 focalfile.open( FOCALS_FILE );
3127
3128 }
3129
3130 // allocate memory for the focals table
3131
3132 ct_Focal = new float[ct_NMirrors];
3133
3134 if ( verbose >= VERBOSE_MINIMAL ) {
3135
3136 if ( errors > 0 )
3137 log( "read_focals", "Reading focals file after %d fail%c\n",
3138 errors, ((errors>1) ? 's' : ' ') );
3139
3140 else
3141
3142 log( "read_focals", "Reading focals file . . .\n" );
3143 }
3144
3145 // scan the whole file
3146
3147 i = 0;
3148
3149 while ( (! focalfile.eof()) && (i<ct_NMirrors) ) {
3150
3151 // get line from the file
3152
3153 focalfile.getline(line, LINE_MAX_LENGTH);
3154
3155 // skip if comment
3156
3157 if ( *line == '#' )
3158 continue;
3159
3160 // get the value
3161
3162 sscanf(line, "%f", &ct_Focal[i]);
3163
3164 ++i;
3165
3166 }
3167
3168 // did i do it for all the mirrors ?
3169
3170 if (i < ct_NMirrors) {
3171 error("read_focals", "%s%s (%d < %d)\n",
3172 "Number of data in focals.dat file is",
3173 "\n\t\tsmaller than number of mirrors", i, ct_NMirrors);
3174 }
3175
3176 // close it
3177
3178 focalfile.close();
3179
3180}
3181//!@}
3182
3183
3184//!---------------------------------------------------------------------
3185// @name rnormal
3186//
3187// @desc returns n(=2k) normaly distributed numbers
3188//
3189// @var *r pointer to a vector where we write the numbers
3190// @var n how many numbers do we generate
3191//
3192// @date Sat Jun 27 05:58:56 MET DST 1998
3193//----------------------------------------------------------------------
3194// @function
3195
3196//!@{
3197void
3198rnormal(double *r, int n)
3199{
3200
3201 double z1, z2;
3202 int i;
3203
3204 for (i=0; i<n; i+=2) {
3205
3206 z1 = RandomNumber;
3207 z2 = RandomNumber;
3208
3209 r[i] = sqrt(-2.0*log(z1)) * cos(2.0*M_PI*z2);
3210 r[i+1] = sqrt(-2.0*log(z1)) * sin(2.0*M_PI*z2);
3211
3212 }
3213
3214}
3215//!@}
3216
3217
3218//!---------------------------------------------------------------------
3219// @name isA
3220//
3221// @desc returns TRUE(FALSE), if the flag is(is not) the given
3222//
3223// @var s1 String to be compared with s2
3224// @var s2 String in the form of a FLAG
3225//
3226// @return TRUE: both strings match; FALSE: otherwise
3227//
3228// @date Wed Jul 8 15:25:39 MET DST 1998
3229//----------------------------------------------------------------------
3230// @function
3231
3232//!@{
3233int
3234isA( char * s1, const char * flag )
3235{
3236 return ( (strncmp((char *)s1, flag, 8)==0) ? TRUE : FALSE );
3237}
3238//!@}
3239
3240
3241//!---------------------------------------------------------------------
3242// @name Curv2Lin
3243//
3244// @desc Curvilinear to Linear (Euclidean) distance
3245//
3246// @var s Curvilinear distance over the parabolic shape
3247//
3248// @return Radial distance from the axis of the paraboloid
3249//
3250// @date Wed Jul 8 15:25:39 MET DST 1998
3251//----------------------------------------------------------------------
3252// @function
3253
3254//!@{
3255float
3256Curv2Lin(float s)
3257{
3258 float x;
3259 short i;
3260
3261 x = s;
3262 for (i = 0; i < 4; i++)
3263 x = (s / 100.) / (1. + (float) 0.000144175317185 * x * x);
3264
3265 return (x*100.);
3266}
3267//!@}
3268
3269
3270//!---------------------------------------------------------------------
3271// @name Lin2Curv
3272//
3273// @desc Linear (Euclidean) to Curvilinear distance
3274//
3275// @var x Radial distance from the axis of the paraboloid
3276//
3277// @return Curvilinear distance over the parabolic shape
3278//
3279// @date Wed Jul 8 15:25:39 MET DST 1998
3280//----------------------------------------------------------------------
3281// @function
3282
3283//!@{
3284float
3285Lin2Curv(float x)
3286{
3287 x /= 100.;
3288 return ((x + (float) 0.000144175317185 * x * x * x)*100.);
3289}
3290//!@}
3291
3292
3293//!---------------------------------------------------------------------
3294// @name get_stdin_files
3295//
3296// @desc get data from stdin and create cerXXXXXX and staXXXXXX
3297//
3298// @var ncerf number of next Cherenkov file cerXXXXXX to be created
3299// @var El Lower limit for requested energy
3300// @var Eu Upper limit for requested energy
3301// @var flag Boolean flag: TRUE: initialize; FALSE: normal
3302//
3303// @return TRUE: succesful; FALSE: failure
3304//
3305// @date Tue Dec 15 15:03:40 MET 1998
3306//----------------------------------------------------------------------
3307// @function
3308
3309//!@{
3310int
3311get_stdin_files(int ncerf, float El, float Eu, int flag)
3312{
3313 static float *buffer; // buffer for read-out of STDIN
3314 static COREventHeader *evth;// Event Header class (from CORSIKA)
3315 static float Elow, Eup; // range for Energy
3316 static int reject;
3317
3318 char *chkbuffer;
3319
3320 char cername[256]; // output filename
3321 ofstream cerfile; // output file (stream)
3322 char staname[256]; // output filename sta
3323 ofstream stafile; // output file sta (stream)
3324 int bytes;
3325 float energy;
3326
3327 /*!@'
3328
3329 The main features of this function, as well as the algorithm,
3330 are taken from |corfilter|.
3331
3332 */
3333
3334 // if it's the first time, allocate buffers
3335
3336 if ( flag == TRUE ) {
3337
3338 if (verbose)
3339 log("get_stdin_files", "Allocating memory for buffer...\n");
3340
3341 buffer = new float [BUFFER_LENGTH];
3342 evth = (COREventHeader*)buffer;
3343
3344 Elow = El;
3345 Eup = Eu;
3346
3347 return ( TRUE );
3348 }
3349
3350 // main loop
3351
3352 while ( ! cin.eof() ) {
3353
3354 // read buffer from STDIN
3355
3356 cin.read( (char*) buffer, BUFFER_LENGTH * sizeof(float) );
3357 bytes = cin.gcount();
3358
3359 if (bytes < (BUFFER_LENGTH * sizeof(float))) {
3360 if (verbose)
3361 logerr("get_stdin_files",
3362 "Can read only %d bytes, instead of required %d.\n",
3363 bytes, BUFFER_LENGTH * sizeof(float));
3364 }
3365
3366 if (bytes == 0)
3367 return ( FALSE );
3368
3369 // check whether EVTH is there
3370
3371 chkbuffer = (char*)buffer;
3372 if (strstr(chkbuffer,EVTH) == 0) {
3373 if (verbose)
3374 logerr("get_stdin_files", "EVTH not found in this block\n");
3375 continue;
3376 } else {
3377 if (verbose)
3378 logerr("get_stdin_files", "EVTH was FOUND in this block!!\n");
3379 }
3380
3381
3382 // show log
3383 if (verbose)
3384 logerr("get_stdin_files", " (E=%8.1f) Reading event %d... ",
3385 (float)evth->Etotal, (int)evth->EvtNumber);
3386
3387 // CRITERIA
3388 energy = evth->Etotal;
3389
3390 if ( (energy<Elow) || (energy>Eup) ) {
3391 if (verbose)
3392 log("REJECTED", "\n");
3393 reject = TRUE;
3394 } else {
3395 if (verbose)
3396 log("ACCEPTED", "\n");
3397 reject = FALSE;
3398 }
3399
3400 if ( !reject ) {
3401
3402 // update no. of cer. file
3403 sprintf(cername, "%s/cer%06d", TMP_STDIN_DIR, ncerf+1);
3404
3405 cerfile.open( cername );
3406
3407 if ( verbose )
3408 log("get_stdin_files", " Writing event %d, will be %s\n",
3409 (int)evth->EvtNumber, cername);
3410
3411 // write it
3412 cerfile.write( (char*)buffer, BUFFER_LENGTH * sizeof(float) );
3413
3414 }
3415
3416 // read blocks till new EVTH is found
3417 // (now must be at the beginning of a block)
3418
3419 if (verbose)
3420 log("get_stdin_files", " Reading photons...\n");
3421
3422 bytes = 0;
3423 while ( TRUE ) {
3424
3425 // read photon
3426 cin.read( (char*)buffer, 7 * sizeof(float) );
3427
3428 // the last photon is reached when the EVTH is read, but
3429 // only when bytes % SIZE_OF_BLOCK == 0
3430 // SIZE_OF_BLOCK = 22932
3431 // with this, we only execute strstr once each 819 photons
3432
3433 // are we in a new block?
3434
3435 if ((bytes % SIZE_OF_BLOCK) == 0) {
3436 // is this the last photon?
3437 if ( strstr(chkbuffer, EVTH) == chkbuffer )
3438 break;
3439 }
3440
3441 // if not, write it to the cer-file
3442 if ( !reject )
3443 cerfile.write( (char*)buffer, 7 * sizeof(float) );
3444
3445 bytes += 7 * sizeof(float);
3446
3447 }
3448
3449 if (verbose)
3450 log("get_stdin_files", " Finished: %d bytes\n", bytes);
3451
3452 if ( !reject ) {
3453 // close cer-file
3454 cerfile.close();
3455 }
3456
3457 // at this point, we must be reading sta-file
3458
3459 // read the rest of the sta-file
3460 if (verbose)
3461 log("get_stdin_files", " Reading sta-data...\n");
3462 cin.read( (char*)buffer + 7*sizeof(float), 5956 - 7*sizeof(float));
3463
3464 if ( !reject ) {
3465
3466 // log
3467 if (verbose)
3468 log("get_stdin_files", " Saving sta-file...\n");
3469
3470 // open file
3471 sprintf(staname, "%s/sta%06d", TMP_STDIN_DIR, ncerf+1);
3472 stafile.open( staname );
3473
3474 // write it
3475 stafile.write( (char*)buffer, 5956 );
3476
3477 // close sta-file
3478 stafile.close();
3479
3480 // get out from here, return to main program
3481 return( TRUE );
3482
3483 }
3484
3485 }
3486
3487 return( TRUE );
3488
3489}
3490//!@}
3491
3492
3493//!---------------------------------------------------------------------
3494// @name get_new_ct_pointing
3495//
3496// @desc returns new random position of CT around a given one
3497//
3498// @var theta Theta (ZA) of the shower
3499// @var phi Phi (AZ) of the shower
3500// @var range Maximum allowed distance between CT and shower axis
3501// @var newtheta Resulting Theta of the CT
3502// @var newphi Resulting Phi of the CT
3503//
3504// @return Angular distance between orig. and new directions
3505//
3506// @date Sat Jun 27 05:58:56 MET DST 1998
3507//----------------------------------------------------------------------
3508// @function
3509
3510//!@{
3511float
3512get_new_ct_pointing(float theta, float phi, float range,
3513 float *newtheta, float *newphi)
3514{
3515 float distance;
3516 float it, ip;
3517 float sin_theta, cos_theta;
3518 float sin_newtheta, cos_newtheta;
3519 float sin_iphi, cos_iphi;
3520 float iphi;
3521
3522 // for the moment, we only simulate an uniform distribution,
3523 // since our theta distribution in the generation of events is
3524 // already uniform for hadrons, which are the main targets for
3525 // using this option
3526
3527 it = acos(cos(range) + RandomNumber * (1 - cos(range)));
3528 ip = RandomNumber * 2.0 * M_PI;
3529
3530 if ( theta == 0.0 ) {
3531
3532 *newtheta = it;
3533 *newphi = ip;
3534
3535 } else {
3536
3537 sin_theta = sin(theta);
3538 cos_theta = cos(theta);
3539
3540 cos_newtheta = cos_theta*cos(it) + sin_theta*sin(it)*cos(ip);
3541 *newtheta = acos( cos_newtheta );
3542 sin_newtheta = sin( *newtheta );
3543
3544 sin_iphi = sin(it)*sin(ip) / sin_newtheta;
3545 cos_iphi = (( cos(it) - cos_newtheta * cos_theta ) /
3546 ( sin_newtheta * sin_theta ));
3547
3548 iphi = atan2( sin_iphi, cos_iphi );
3549
3550 *newphi = phi + iphi;
3551
3552 }
3553
3554 return( it );
3555}
3556//!@}
3557
3558
3559//=------------------------------------------------------------
3560//!@subsection Log of this file.
3561
3562//!@{
3563//
3564// $Log: not supported by cvs2svn $
3565// Revision 1.9 2000/01/28 09:19:54 harald
3566// A new version from JoseCarlosGonzalez. The old routine for random_pointing
3567// was not correct. This one should be okay!!
3568//
3569// Revision 1.20 2000/01/27 10:47:54 gonzalez
3570// JAN2000-STABLE
3571//
3572// Revision 1.19 1999/11/19 20:52:31 gonzalez
3573// *** empty log message ***
3574//
3575// Revision 1.18 1999/10/05 11:11:12 gonzalez
3576// Sep.1999
3577//
3578// Revision 1.17 1999/03/24 16:33:01 gonzalez
3579// REFLECTOR 1.1: Release
3580//
3581// Revision 1.16 1999/01/21 16:03:44 gonzalez
3582// Only small modifications
3583//
3584// Revision 1.15 1999/01/19 18:07:16 gonzalez
3585// Bugs in STDIN-STDOUT version corrected.
3586//
3587// Revision 1.14 1999/01/14 17:35:47 gonzalez
3588// Both reading from STDIN (data_from_stdin) and
3589// writing to STDOUT (data_to_STDOUT) working.
3590//
3591// Revision 1.13 1999/01/13 12:41:12 gonzalez
3592// THIS IS A WORKING (last?) VERSION
3593//
3594// Revision 1.12 1998/12/17 16:26:09 gonzalez
3595// *************************************************
3596// *************************************************
3597// WARNING:: Version 1.11 is completly wrong!!
3598// *************************************************
3599// *************************************************
3600//
3601// Revision 1.10 1998/12/15 10:47:22 gonzalez
3602// RELEASE-1.0
3603//
3604// Revision 1.9 1998/11/25 16:30:47 gonzalez
3605// Commit after inclusion of 'Blocking'
3606//
3607//!@}
3608
3609//=EOF
Note: See TracBrowser for help on using the repository browser.