Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

camera.cxx@ 10113

Visit:

Last change on this file since 10113 was 308, checked in by harald, 25 years ago

This the starting point for CVS controlled further developments of the camera program. The program was originally written by Jose Carlos. But here you can find a "rootified" version to the program. This means that there is no hbook stuff in it now. Also the output of the program changed to the MagicRawDataFormat. The "rootification" was done by Dirk Petry and Harald Kornmayer. In the following you can see the README file of that version: ================================================== Fri Oct 22 1999 D.P. The MAGIC Monte Carlo System Camera Simulation Programme --------------------------- 1) Description This version is the result of the fusion of H.K.'s root_camera which is described below (section 2) and another version by D.P. which had a few additional useful features. The version compiles under Linux with ROOT 2.22 installed (variable ROOTSYS has to be set). Compile as before simply using "make" in the root_camera directory. All features of H.K.'s root_camera were retained. Additional features of this version are: a) HBOOK is no longer used and all references are removed. b) Instead of HBOOK, the user is given now the possibility of having Diagnostic data in ROOT format as a complement to the ROOT Raw data. This data is written to the file which is determined by the new input parameter "diag_file" in the camera parameter file. All source code file belonging to this part have filenames starting with "MDiag". The user can read the output file using the following commands in an interactive ROOT session: root [0] .L MDiag.so root [1] new TFile("diag.root"); root [2] new TTreeViewer("T"); This brings up a viewer from which all variables of the TTree can be accessed and histogrammed. This example assumes that you have named the file "diag.root", that you are using ROOT version 2.22 or later and that you have the shared object library "MDiag.so" which is produced by the Makefile along with the executable "camera". ! The contents of the so-called diag file is not yet fixed. ! At the moment it is what J.C.G. used to put into the HBOOK ! ntuple. In future versions the moments calculation can be ! removed and the parameter list be modified correspondingly. c) Now concatenated reflector files can be read. This is useful if you have run the reflector with different parameters but you want to continue the analysis with all reflector data going into ONE ROOT outputfile. The previous camera version contained a bug which made reading of two or more concatenated reflector files impossible. d) The reflector output format was changed. It is now version 0.4 . The change solely consists in a shortening of the flag definition in the file include-MC/MCCphoton.hxx ! IF YOU WANT TO READ REFLECTOR FORMAT 0.3, you can easily ! do so by recompiling camera with the previous version of ! include-MC/MCCphoton.hxx. The change was necessary for saving space and better debugging. From now on, this format can be frozen. ! For producing reflector output in the new format, you ! of course have to recompile your reflector with the ! new include-MC/MCCphoton.hxx . e) A first version of the pixelization with the larger outer pixels is implemented. THIS IS NOT YET FULLY TESTED, but first rough tests show that it works at least to a good approximation. The present version implements the camera outline with 18 "gap-pixels" and 595 pixels in total as shown in http://sarastro.ifae.es/internal/home/hardware/camera/numbering.ps This change involved (i) The file pixels.dat is no longer needed. Instead the coordinates are generated by the program itself (takes maybe 1 second). In the file pixel-coords.txt in the same directory as this README, you find a list of the coordinates generated by this new routine. It has the format number i j x y size-factor where i and j are J.C.G.'s so called biaxis hexagonal coordinates (for internal use) and x and y are the coordinates of the pixel centers in the standard camera coordinate system in units of centimeters. The value of "size-factor" determines the linear size of the pixel relative to the central pixels. (ii) The magic.def file has two additional parameters which give the number of central pixels and the number of gap pixels (iii) In camera.h and camera.cxx several changes were necessary, among them the introduction of several new functions The newly suggested outline with asymmetric Winston cones will be implemented in a later version. f) phe files can no longer be read since this contradicts our philosophy that the analysis should be done with other programs like e.g. EVITA and not with "camera" itself. This possibility was removed. g) ROOT is no longer invoked with an interactive interface. In this way, camera can better be run as a batch program and it uses less memory. h) small changes concerning the variable "t_chan" were necessary in order to avoid segmentation faults: The variable is used as an index and it went sometimes outside the limits when camera was reading proton data. This is because the reflector files don't contain the photons in a chronological order and also the timespread can be considerably longer that the foreseen digitisation timespan. Please see the source code of camera.cxx round about line 1090. j) several unused variables were removed, a few warning messages occur when you compile camera.cxx but these can be ignored at the moment. In general the program is of course not finished. It still needs debugging, proper trigger simulation, simulation of the asymmetric version of the outer pixels, proper NSB simulation, adaption of the diag "ntuple" contents to our need and others small improvements. In the directory rfl-files there is now a file in reflector format 0.4 containing a single event produced by the starfiled adder. It has a duration of 30 ns and represents the region around the Crab Nebula. Using the enclosed input parameter file, camera should process this file without problems. 2) The README for the previous version of root_camera README for a preliminary version of the root_camera program. root_camera is based on the program "camera"of Jose Carlos Gonzalez. It was changed in the way that only the pixelisation and the distibution of the phe to the FADCs works in a first version. Using the #undef command most possibilities of the orignal program are switched of. The new parts are signed by - ROOT or __ROOT__ nearly all important codelines for ROOT output are enclosed in structures like #ifdef __ROOT__ code #endif __ROOT__ In same case the new lines are signed by a comment with the word ROOT in it. For timing of the pulse some variable names are changed. (t0, t1, t --> t_ini, t_fin, t_1st, t_chan,...) Look also for this changes. For the new root-file is also a change in readparm-files - __DETAIL_TRIGGER__ This is for the implementation of the current work on trigger studies. Because the class MTrigger is not well documented it isn´t a part of this tar file. Only a dummy File exists. With all files in the archive, the root_camera program should run. A reflector file is in the directory rfl-files ================================================== From now on, use CVS for development!!!!

File size: 78.6 KB

Line
1	//!/////////////////////////////////////////////////////////////////////
2	//
3	// camera
4	//
5	// @file camera.cxx
6	// @title Camera simulation
7	// @subtitle Code for the simulation of the camera phase
8	// @desc Code for the simulation of the camera of CT1 and MAGIC
9	// @author J C Gonzalez
10	// @email gonzalez@mppmu.mpg.de
11	// @date Thu May 7 16:24:22 1998
12	//
13	//----------------------------------------------------------------------
14	//
15	// Created: Thu May 7 16:24:22 1998
16	// Author: Jose Carlos Gonzalez
17	// Purpose: Program for reflector simulation
18	// Notes: See files README for details
19	//
20	//----------------------------------------------------------------------
21	//
22	// $RCSfile: camera.cxx,v $
23	// $Revision: 1.1.1.1 $
24	// $Author: harald $
25	// $Date: 1999-11-05 11:59:31 $
26	//
27	////////////////////////////////////////////////////////////////////////
28	// @tableofcontents @coverpage
29
30	//=-----------------------------------------------------------
31	//!@section Source code of \|camera.cxx\|.
32
33	/*!@"
34
35	In this section we show the (commented) code of the program for the
36	read-out of the output files generated by the simulator of the
37	reflector, \|reflector 0.3\|.
38
39	@"*/
40
41	//=-----------------------------------------------------------
42	//!@subsection Includes and Global variables definition.
43
44	//!@{
45
46	// includes for ROOT
47	// BEWARE: the order matters!
48
49	#include "TROOT.h"
50
51	#include "TFile.h"
52	#include "TTree.h"
53	#include "TBranch.h"
54
55	#include "MDiag.h"
56
57	#include "MTrigger.hxx"
58
59	#include "MRawEvt.h"
60	#include "MMcEvt.h"
61
62	/*!@"
63
64	All the defines are located in the file \|camera.h\|.
65
66	@"*/
67
68	#include "camera.h"
69	//!@}
70
71	/*!@"
72
73	The following set of flags are used in time of compilation. They do
74	not affect directly the behaviour of the program at run-time
75	(though, of course, if you disconnected the option for
76	implementation of the Trigger logic, you will not be able to use any
77	trigger at all. The 'default' values mean default in the sense of
78	what you got from the server when you obtained this program.
79
80	@"*/
81
82	//!@{
83
84	// flag for debugging (default: OFF )
85	#define __DEBUG__
86	#undef __DEBUG__
87
88	// flag for NNT in CT1 camera (default: ON )
89	#undef __CT1_NO_NEIGHBOURS__
90	#define __CT1_NO_NEIGHBOURS__
91
92	// flag for calculation of NSB (default: ON )
93	#undef __NSB__
94	#define __NSB__
95
96	// flag for calculation of QE for pixels (default: ON )
97	#undef __QE__
98	#define __QE__
99
100
101	// flag for implementation of DETAIL_TRIGGER (default: ON )
102	//
103	// This is the new implementation of Trigger studies
104	// It relies on a better simulation of the time stucture
105	// of the PhotoMultiplier. For more details see the
106	// documentation of the --> class MTrigger <--
107	#define __DETAIL_TRIGGER__
108	#undef __DETAIL_TRIGGER__
109
110	// flag for implementation of TRIGGER (default: ON )
111	#undef __TRIGGER__
112	#define __TRIGGER__
113
114	// flag for implementation of Tail-Cut (default: ON )
115	#undef __TAILCUT__
116	#define __TAILCUT__
117
118	// flag for calculation of islands stat. (default: ON )
119	#define __ISLANDS__
120	#undef __ISLANDS__
121
122	// flag for calculation of image parameters (default: ON )
123	#undef __MOMENTS__
124	#define __MOMENTS__
125
126	// flag for ROOT (default: ON )
127	#undef __ROOT__
128	#define __ROOT__
129
130	//!@}
131
132	//=-----------------------------------------------------------
133	//!@subsection Definition of global variables.
134
135	/*!@"
136
137	Now we define some global variables with data about the telescope,
138	such as "focal distance", number of pixels/mirrors,
139	"size of the camera", and so on.
140
141	@"*/
142
143	/*!@"
144
145	Depending on the telescope we are using (CT1 or MAGIC), the
146	information stored in the definition file is different.
147	The variable \|ct_Type\| has the value 0 when we use
148	CT1, and 1 when we use MAGIC.
149
150	@"*/
151
152	//!@{
153	static int ct_Type; //@< Type of telescope: 0:CT1, 1:MAGIC
154	//!@}
155
156	/*!@"
157
158	And this is the information about the whole telescope.
159
160	@"*/
161
162	//!@{
163
164	// parameters of the CT (from the CT definition file)
165
166	////@: Focal distances [cm]
167	//static float *ct_Focal;
168
169	//@: Mean Focal distances [cm]
170	static float ct_Focal_mean;
171
172	//@: STDev. Focal distances [cm]
173	static float ct_Focal_std;
174
175	//@: Mean Point Spread function [cm]
176	static float ct_PSpread_mean;
177
178	//@: STDev. Point Spread function [cm]
179	static float ct_PSpread_std;
180
181	//@: STDev. Adjustmente deviation [cm]
182	static float ct_Adjustment_std;
183
184	//@: Radius of the Black Spot in mirror [cm]
185	static float ct_BlackSpot_rad;
186
187	//@: Radius of one mirror [cm]
188	static float ct_RMirror;
189
190	//@: Camera width [cm]
191	static float ct_CameraWidth;
192
193	//@: Pixel width [cm]
194	static float ct_PixelWidth;
195
196	//@: ct_PixelWidth_corner_2_corner = ct_PixelWidth / cos(60)
197	static float ct_PixelWidth_corner_2_corner;
198
199	//@: ct_PixelWidth_corner_2_corner / 2
200	static float ct_PixelWidth_corner_2_corner_half;
201
202	//@: Number of mirrors
203	static int ct_NMirrors = 0;
204
205	//@: Number of pixels
206	static int ct_NPixels;
207
208	//@: Number of pixels
209	static int ct_NCentralPixels;
210
211	//@: Number of pixels
212	static int ct_NGapPixels;
213
214	//@: ct_Apot = ct_PixelWidth / 2
215	static float ct_Apot;
216
217	//@: ct_2Apot = 2 * ct_Apot = ct_PixelWidth
218	static float ct_2Apot;
219
220	//@: name of the CT definition file to use
221	static char ct_filename[256];
222
223	//@: list of showers to be skipped
224	static int *Skip;
225
226	//@: number of showers to be skipped
227	static int nSkip=0;
228
229	//@: flag: TRUE: data come from STDIN; FALSE: from file
230	static int Data_From_STDIN = FALSE;
231
232	//@: flag: TRUE: write all images to output; FALSE: only triggered showers
233	static int Write_All_Images = FALSE;
234
235	//@: flag: TRUE: write all data to output; FALSE: only triggered showers
236	static int Write_All_Data = FALSE;
237
238	//@: flag: TRUE: selection on the energy
239	static int Select_Energy = TRUE;
240
241	//@: Lower edge of the selected energy range (in GeV)
242	static float Select_Energy_le = 0.0;
243
244	//@: Upper edge of the selected energy range (in GeV)
245	static float Select_Energy_ue = 100000.0;
246
247	//!@}
248
249	/*!@"
250
251	The following double-pointer is a 2-dimensional table with information
252	about each pixel. The routine read_pixels will generate
253	the information for filling it using igen_pixel_coordinates().
254
255	@"*/
256
257	//!@{
258	// Pointer to a tables/Arrays with information about the pixels
259	// and data stored on them with information about the pixels
260
261	//@: table for IJ sys.
262	static float pixels[PIX_ARRAY_SIDE][PIX_ARRAY_SIDE][4];
263
264	//@: coordinates x,y for each pixel
265	static float **pixary;
266
267	//@: indexes of pixels neighbours of a given one
268	static int **pixneig;
269
270	//@: number of neighbours a pixel have
271	static int *npixneig;
272
273	//@: contents of the pixels (ph.e.)
274	static float *fnpix;
275
276	//@: contents of the pixels (ph.e.) after cleanning
277	static float *fnpixclean;
278
279	//@: contents of the sum of all ph.e. in one timeslice of 1 ns
280	static float *fnslicesum ;
281
282
283	//!@}
284
285	/*!@"
286
287	The following double-pointer is a 2-dimensional table with the
288	Quantum Efficiency @$QE@$ of each pixel in the camera, as a function
289	of the wavelength @$\lambda@$. The routine \|read_pixels()\| will read
290	also this information from the file \|qe.dat\|.
291
292	@"*/
293
294	//!@{
295	// Pointer to a table with QE, number of datapoints, and wavelengths
296
297	//@: table of QE
298	static float ***QE;
299
300	//@: number of datapoints for the QE curve
301	static int pointsQE;
302
303	//@: table of QE
304	static float *QElambda;
305	//!@}
306
307	/*!@"
308
309	The following double-pointer is a 2-dimensional table with information
310	about each mirror in the dish. The routine \|read_ct_file()\| will read
311	this information from the CT definition file.
312
313	@"*/
314
315	//!@{
316	// Pointer to a table with the following info.:
317
318	static float **ct_data;
319
320	/*
321	* TYPE=0 (CT1)
322	* i s rho theta x y z thetan phin xn yn zn
323	*
324	* i : number of the mirror
325	* s : arc length [cm]
326	* rho : polar rho of the position of the center of the mirror [cm]
327	* theta : polar angle of the position of the center of the mirror [cm]
328	* x : x coordinate of the center of the mirror [cm]
329	* y : y coordinate of the center of the mirror [cm]
330	* z : z coordinate of the center of the mirror [cm]
331	* thetan : polar theta angle of the direction where the mirror points to
332	* phin : polar phi angle of the direction where the mirror points to
333	* xn : xn coordinate of the normal vector in the center (normalized)
334	* yn : yn coordinate of the normal vector in the center (normalized)
335	* zn : zn coordinate of the normal vector in the center (normalized)
336	*
337	* TYPE=1 (MAGIC)
338	* i f sx sy x y z thetan phin
339	*
340	* i : number of the mirror
341	* f : focal distance of that mirror
342	* sx : curvilinear coordinate of mirror's center in X[cm]
343	* sy : curvilinear coordinate of mirror's center in X[cm]
344	* x : x coordinate of the center of the mirror [cm]
345	* y : y coordinate of the center of the mirror [cm]
346	* z : z coordinate of the center of the mirror [cm]
347	* thetan : polar theta angle of the direction where the mirror points to
348	* phin : polar phi angle of the direction where the mirror points to
349	* xn : xn coordinate of the normal vector in the center (normalized)
350	* yn : yn coordinate of the normal vector in the center (normalized)
351	* zn : zn coordinate of the normal vector in the center (normalized)
352	*/
353	//!@}
354
355	/*!@"
356
357	We define a table into where random numbers will be stored.
358	The routines used for random number generation are provided by
359	\|RANLIB\| (taken from NETLIB, \|www.netlib.org\|), and by
360	the routine \|double drand48(void)\| (prototype defined in
361	\|stdlib.h\|) through the macro \|RandomNumber\| defined in
362	\|camera.h\|.
363
364	@"*/
365
366	//!@{
367	// table of random numbers
368
369	// (unused)
370	//static double RandomNumbers[500];
371	//!@}
372
373	/*!@"
374
375	The following is a variable to count the number of Cphotons
376	in the different steps of the simulation.
377	The definition is as follows:
378	@[
379	\mbox{CountCphotons}[ \mbox{FILTER} ] \equiv
380	\mbox{\it Number of photons after the filter} \mbox{FILTER}
381	@]
382	The filters are defined and can be found in the file \|camera.h\|.
383
384	@"*/
385
386	//!@{
387	// vector to count photons at any given step of the simulation
388
389	static int CountCphotons[10];
390	//!@}
391
392	/*!@"
393
394	The following are the set of parameters calculated for each image.
395	The routines for their calculations are in \|moments.cxx\|.
396
397	@"*/
398
399	//!@{
400	// parameters of the images
401
402	static Moments_Info *moments_ptr;
403	static LenWid_Info *lenwid_ptr;
404
405	static float *maxs;
406	static int *nmaxs;
407	static float length, width, dist, xdist, azw, miss, alpha, *conc;
408	static float phiasym, asymx, asymy;
409	static float charge, smax, maxtrigthr_phe;
410
411	//!@}
412
413	extern char FileName[];
414
415
416	//=-----------------------------------------------------------
417	// @subsection Main program.
418
419	//!@{
420
421	//++++++++++++++++++++++++++++++++++++++++
422	// MAIN PROGRAM
423	//----------------------------------------
424
425	int main(int argc, char **argv)
426	{
427
428	//!@' @#### Definition of variables.
429	//@'
430
431	char inname[256]; //@< input file name
432	char outname[256]; //@< output file name
433	char datname[256]; //@< data (ASCII) output file name
434	char diagname[256]; //@< diagnistic output file (ROOT format)
435	char rootname[256] ; //@< ROOT file name
436
437	char parname[256]; //@< parameters file name
438
439	char sign[20]; //@< initialize sign
440
441	char flag[SIZE_OF_FLAGS + 1]; //@< flags in the .rfl file
442
443	ifstream inputfile; //@< stream for the input file
444	ofstream outputfile; //@< stream for the output file
445	ofstream datafile; //@< stream for the data file
446
447	MCEventHeader mcevth; //@< Event Header class (MC)
448	MCCphoton cphoton; //@< Cherenkov Photon class (MC)
449
450	float thetaCT, phiCT; //@< parameters of a given shower
451	float thetashw, phishw; //@< parameters of a given shower
452	float coreD, coreX, coreY; //@< core position and distance
453	float impactD; //@< impact parameter
454	float l1, m1, n1; //@< auxiliary variables
455	float l2, m2, n2; //@< auxiliary variables
456	float num, den; //@< auxiliary variables
457
458	int nshow=0; //@< partial number of shower in a given run
459	int ntshow=0; //@< total number of showers
460	int ncph=0; //@< partial number of shower in a given run
461	int ntcph=0; //@< total number of showers
462
463	int i, ii, j, k; //@< simple counters
464
465	float t_ini; //@< time of the first Cphoton read in
466	float t; //@< time for a single photon
467	int t_chan ; //@< the bin (channel) in time of a single photon
468
469	int startchan ; //@< the first bin with entries in the time slices
470
471	float cx, cy, ci, cj; //@< coordinates in the XY and IJ systems
472	int ici, icj, iici, iicj; //@< coordinates in the IJ (integers)
473
474	int nPMT; //@< number of pixel
475
476	float wl, last_wl; //@< wavelength of the photon
477	float qe; //@< quantum efficiency
478	float **qeptr;
479
480	int simulateNSB; //@< Will we simulate NSB?
481	float meanNSB; //@< NSB mean value (per pixel)
482	float qThreshold; //@< Threshold value
483	float qTailCut; //@< Tail Cut value
484	int nIslandsCut; //@< Islands Cut value
485	int countIslands; //@< Will we count the islands?
486	int anaPixels;
487
488	float fCorrection; //@< Factor to apply to pixel values (def. 1.)
489
490	float q0; //@< trigger threshold ( intermediate variable )
491	float maxcharge; //@< maximum charge in pixels
492	int noverq0, novq0; //@< number of pixels above threshold
493	int ngrpq0, mxgrp; //@< number of pixels in a group
494
495	int trigger; //@< trigger flag
496	int itrigger; //@< index of pixel fired
497	int ntrigger = 0; //@< number of triggers in the whole file
498	unsigned char triggerBits; //@< byte for trigger condition check (MAGIC)
499	int bit; //@< intermediate variable
500
501	float plateScale_cm2deg; //@< plate scale (deg/cm)
502	float degTriggerZone; //@< trigger area in the camera (radius, in deg.)
503
504	float dtheta, dphi; //@< deviations of CT from shower axis
505
506	int still_in_loop = FALSE;
507
508	char Signature[20];
509
510	float *image_data;
511	int nvar, hidt;
512
513	struct camera cam; // structure holding the camera definition
514
515	#ifdef __DETAIL_TRIGGER__
516
517	MTrigger Trigger ; //@< A instance of the Class MTrigger
518
519	#endif __DETAIL_TRIGGER__
520
521	//!@' @#### Definition of variables for \|getopt()\|.
522	//@'
523
524	int ch, errflg = 0; //@< used by getopt
525
526	/*!@'
527
528	@#### Beginning of the program.
529
530	We start with the main program. First we (could) make some
531	presentation, and follows the reading of the parameters file (now
532	from the \|stdin\|), the reading of the CT parameters file, and the
533	creation of the output file, where the processed data will be
534	stored.
535
536	*/
537
538	//++
539	// START
540	//--
541
542	// make unbuffered output
543
544	cout.setf ( ios::stdio );
545
546	// parse command line options (see reflector.h)
547
548	parname[0] = '\0';
549
550	optarg = NULL;
551	while ( !errflg && ((ch = getopt(argc, argv, COMMAND_LINE_OPTIONS)) != -1) )
552	switch (ch) {
553	case 'f':
554	strcpy(parname, optarg);
555	break;
556	case 'h':
557	usage();
558	break;
559	default :
560	errflg++;
561	}
562
563	// show help if error
564
565	if ( errflg>0 )
566	usage();
567
568	// make some sort of presentation
569
570	present();
571
572	// read parameters file
573
574	if ( strlen(parname) < 1 )
575	readparam(NULL);
576	else
577	readparam(parname);
578
579	// read data from file or from STDIN?
580
581	Data_From_STDIN = get_data_from_stdin();
582
583	// write all images, even those without trigger?
584
585	Write_All_Images = get_write_all_images();
586
587	// write all data (i.e., ph.e.s in pixels)
588
589	Write_All_Data = get_write_all_data();
590
591	// get filenames
592
593	strcpy( inname, get_input_filename() );
594	strcpy( outname, get_output_filename() );
595	strcpy( datname, get_data_filename() );
596	strcpy( diagname, get_diag_filename() );
597	strcpy( rootname, get_root_filename() );
598	strcpy( ct_filename, get_ct_filename() );
599
600	// get different parameters of the simulation
601
602	qThreshold = get_threshold();
603	qTailCut = get_tail_cut();
604	simulateNSB = get_nsb( &meanNSB );
605	countIslands = get_islands_cut( &nIslandsCut );
606
607	// get selections on the parameters
608
609	Select_Energy = get_select_energy( &Select_Energy_le, &Select_Energy_ue);
610
611	// log filenames information
612
613	log(SIGNATURE,
614	"%s:\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s:\t%s\n",
615	"Filenames",
616	"In", inname,
617	"Out", outname,
618	"Data", datname,
619	"Diag", diagname,
620	"ROOT", rootname,
621	"CT", ct_filename);
622
623
624	// log flags information
625
626	log(SIGNATURE,
627	"%s:\n\t%20s: %s\n\t%20s: %s\n\t%20s: %s\n",
628	"Flags",
629	"Data_From_STDIN", ONoff(Data_From_STDIN),
630	"Write_All_Images", ONoff(Write_All_Images),
631	"Write_All_Data", ONoff(Write_All_Data));
632
633	// log parameters information
634
635	log(SIGNATURE,
636	"%s:\n\t%20s: %f\n\t%20s: %f\n\t%20s: %f %s\n\t%20s: %f %s\n",
637	"Parameters",
638	"q0 (Threshold)", qThreshold,
639	"t0 (Tail-cut)", qTailCut,
640	"NSB (phes/pixel)", meanNSB, ONoff(simulateNSB),
641	"i0 (Islands-cut)", nIslandsCut, ONoff(countIslands));
642
643	// log selections
644
645	log(SIGNATURE,
646	"%s:\n\t%20s: %s (%f:%f)\n",
647	"Selections:",
648	"Energy", ONoff(Select_Energy), Select_Energy_le, Select_Energy_ue);
649
650	// set all random numbers seeds
651
652	setall( get_seeds(0), get_seeds(1) );
653
654	// get list of showers to evt. skip
655
656	nSkip = get_nskip_showers();
657
658	if (nSkip > 0) {
659	Skip = new int[ nSkip ];
660	get_skip_showers( Skip );
661
662	log(SIGNATURE, "There are some showers to skip:\n");
663	for (i=0; i<nSkip; ++i)
664	log(SIGNATURE, "\tshower # %d\n", Skip[i]);
665	}
666
667	// read parameters from the ct.def file
668
669	read_ct_file();
670
671	// read pixels data
672
673	read_pixels(&cam);
674
675	// initialise ROOT
676
677	TROOT simple("simple", "MAGIC Telescope Monte Carlo");
678
679	// prepare ROOT tree for the diagnostic data
680
681	TFile *hfile;
682
683	hfile = new TFile( diagname,"RECREATE", "MAGIC Telescope MC diagnostic data");
684
685
686	// Create the ROOT Tree for the diagnostic data
687
688	TTree *tree = new TTree("T","MAGIC Telescope MC diagnostic data");
689	tree->SetAutoSave(100000000);
690
691	Int_t split = 1;
692	Int_t bsize = 64000;
693	MDiagEventobject *event = 0;
694
695	// Create one branch. If splitlevel is set, event is a superbranch
696	// creating a sub branch for each data member of the Eventobject event.
697
698	tree->Branch("event", "MDiagEventobject", &event, bsize, split);
699
700	#ifdef __ROOT__
701
702	MRawEvt *Evt = new MRawEvt() ;
703	MMcEvt *McEvt = new MMcEvt ();
704
705	// initalize the ROOT file
706	//
707	// erzeuge ein Root file
708	//
709
710	TFile outfile ( rootname , "RECREATE" ) ;
711
712	//
713	// create a Tree for the Event data stream
714	//
715
716	TTree EvtTree("EvtTree","Events of Run");
717
718	bsize=128000; split=1;
719
720	EvtTree.Branch("MRawEvt","MRawEvt",
721	&Evt, bsize, split);
722
723	EvtTree.Branch("MMcEvt","MMcEvt",
724	&McEvt, bsize, split);
725
726
727	float flli = 0. ;
728	unsigned short ulli = 0 ;
729
730	#endif __ROOT__
731
732	// for safety and for dimensioning image_data: count the elements in the
733	// diagnostic data branch
734
735	i=0;
736	i++; // "n"
737	i++; // "primary"
738	i++; // "energy"
739	i++; // "cored"
740	i++; // "impact"
741	i++; // "xcore"
742	i++; // "ycore"
743	i++; // "theta"
744	i++; // "phi"
745	i++; // "deviations"
746	i++; // "dtheta"
747	i++; // "dphi"
748	i++; // "trigger"
749	i++; // "ncphs"
750	i++; // "maxpassthr_phe"
751	i++; // "nphes"
752	i++; // "nphes2"
753	i++; // "length"
754	i++; // "width"
755	i++; // "dist"
756	i++; // "xdist"
757	i++; // "azw"
758	i++; // "miss"
759	i++; // "alpha"
760	i++; // "conc2"
761	i++; // "conc3"
762	i++; // "conc4"
763	i++; // "conc5"
764	i++; // "conc6"
765	i++; // "conc7"
766	i++; // "conc8"
767	i++; // "conc9"
768	i++; // "conc10"
769	i++; // "asymx"
770	i++; // "asymy"
771	i++; // "phiasym"
772
773	nvar = i;
774	image_data = new float[nvar];
775
776	// set plate scale (deg/cm) and trigger area (deg)
777
778	plateScale_cm2deg = ( ct_Type == 0 ) ? (0.244/2.1) : 0.030952381;
779
780	if ( ! get_trigger_radius( &degTriggerZone ) )
781	degTriggerZone = ( ct_Type == 0 ) ? (5.0) : (5.0);
782
783	if ( ! get_correction( &fCorrection ) )
784	fCorrection = 1.0;
785
786	// number of pixels for parameters
787
788	anaPixels = get_ana_pixels();
789	anaPixels = (anaPixels == -1) ? ct_NPixels : anaPixels;
790
791	// open input file if we DO read data from a file
792
793	if (! Data_From_STDIN) {
794	log( SIGNATURE, "Openning input \"rfl\" file %s\n", inname );
795	inputfile.open( inname );
796	if ( inputfile.bad() )
797	error( SIGNATURE, "Cannot open input file: %s\n", inname );
798	}
799
800	// get signature, and check it
801	// NOTE: this part repeats further down in the code;
802	// if you change something here you probably want to change it
803	// there as well
804
805	strcpy(Signature, REFL_SIGNATURE);
806
807	strcpy(sign, Signature);
808
809	if ( Data_From_STDIN )
810	cin.read( (char *)sign, strlen(Signature));
811	else
812	inputfile.read( (char *)sign, strlen(Signature));
813
814	if (strcmp(sign, Signature) != 0) {
815	cerr << "ERROR: Signature of .rfl file is not correct\n";
816	cerr << '"' << sign << '"' << '\n';
817	cerr << "should be: " << Signature << '\n';
818	exit(1);
819	}
820
821	if ( Data_From_STDIN )
822	cin.read( (char *)sign, 1);
823	else
824	inputfile.read( (char *)sign, 1);
825
826	// open output file
827
828	log( SIGNATURE, "Openning output \"phe\" file %s\n", outname );
829	outputfile.open( outname );
830
831	if ( outputfile.bad() )
832	error( SIGNATURE, "Cannot open output file: %s\n", outname );
833
834	// open data file
835
836	log( SIGNATURE, "Openning data \"dat\" file %s\n", datname );
837	datafile.open( datname );
838
839	if ( outputfile.bad() )
840	error( SIGNATURE, "Cannot open output file: %s\n", outname );
841
842	// write signature
843
844	outputfile.write( SIGNATURE, sizeof(SIGNATURE) );
845
846	// initializes flag
847
848	strcpy( flag, " \0" );
849
850	// allocate space for PMTs numbers of pixels
851
852	fnpix = new float [ ct_NPixels ];
853	fnpixclean = new float [ ct_NPixels ];
854
855	#ifdef __ROOT__
856
857	fnslicesum = new float [ (2 * SLICES) ] ;
858
859	float slices [ct_NPixels][ (2 * SLICES) ] ;
860	float slices2 [ct_NPixels][ SLICES ] ;
861
862	float trans [ SLICES ] ;
863	#endif __ROOT__
864
865
866	moments_ptr = moments( anaPixels, NULL, NULL, 0.0, 1 );
867
868	//!@' @#### Main loop.
869	//@'
870
871	//begin my version
872
873	// get flag
874
875	if ( Data_From_STDIN )
876	cin.read( flag, SIZE_OF_FLAGS );
877	else
878	inputfile.read ( flag, SIZE_OF_FLAGS );
879
880	// loop over the file
881
882	still_in_loop = TRUE;
883
884	while (
885	((! Data_From_STDIN) && (! inputfile.eof()))
886	\|\|
887	(Data_From_STDIN && still_in_loop)
888	) {
889
890	// reading .rfl files
891	if(!isA( flag, FLAG_START_OF_RUN )){
892	error( SIGNATURE, "Expected start of run flag, but found: %s\n", flag );
893	}
894	else { // found start of run
895	nshow=0;
896	// read next flag
897
898	if ( Data_From_STDIN )
899	cin.read( flag, SIZE_OF_FLAGS );
900	else
901	inputfile.read ( flag, SIZE_OF_FLAGS );
902
903	while( isA( flag, FLAG_START_OF_EVENT )){ // while there is a next event
904	/*!@'
905
906	For the case \|FLAG\_START\_OF\_EVENT\|,
907	we read each Cherenkov photon, and follow these steps:
908
909	@enumerate
910
911	@- Transform XY-coordinates to IJ-coordinates.
912
913	@- With this, we obtain the pixel where the photon hits.
914
915	@- Use the wavelength $\lambda$ and the table of QE, and
916	calculate the estimated (third order interpolated) quantum
917	efficiency for that photon. The photon can be rejected.
918
919	@- If accepted, then add to the pixel.
920
921	@endenumerate
922
923	In principle, we should stop here, and use another program to
924	'smear' the image, to add the Night Sky Background, and to
925	simulate the trigger logic, but we will make this program
926	quick and dirty, and include all here.
927
928	If we are reading PHE files, we jump to the point where the
929	pixelization process already has finished.
930
931	*/
932
933	++nshow;
934	log(SIGNATURE, "Event %d(+%d)\n", nshow, ntshow);
935
936	// get MCEventHeader
937
938	if ( Data_From_STDIN )
939	cin.read( (char*)&mcevth, mcevth.mysize() );
940	else
941	mcevth.read( inputfile );
942
943	// calculate core distance and impact parameter
944
945	coreD = mcevth.get_core(&coreX, &coreY);
946
947	// calculate impact parameter (shortest distance betwee the original
948	// trajectory of the primary (assumed shower-axis) and the
949	// direction where the telescope points to
950	//
951	// we use the following equation, given that the shower core position
952	// is (x1,y1,z1)=(x,y,0),the trajectory is given by (l1,m1,n1),
953	// and the telescope position and orientation are (x2,y2,z2)=(0,0,0)
954	// and (l2,m2,n2)
955	//
956	// \| \|
957	// \| x1-x2 y1-y2 z1-z2 \|
958	// \| \|
959	// + \| l1 m1 n1 \|
960	// - \| \|
961	// \| l2 m2 n2 \|
962	// \| \|
963	// dist = ------------------------------------ ( > 0 )
964	// [ \|l1 m1\|2 \|m1 n1\|2 \|n1 l1\|2 ]1/2
965	// [ \| \| + \| \| + \| \| ]
966	// [ \|l2 m2\| \|m2 n2\| \|n2 l2\| ]
967	//
968	// playing a little bit, we get this reduced for in our case:
969	//
970	//
971	// dist = (- m2 n1 x + m1 n2 x + l2 n1 y - l1 n2 y - l2 m1 z + l1 m2 z) /
972	// [(l2^2 (m1^2 + n1^2) + (m2 n1 - m1 n2)^2 -
973	// 2 l1 l2 (m1 m2 + n1 n2) + l1^2 (m2^2 + n2^2) ] ^(1/2)
974
975	// read the direction of the incoming shower
976
977	thetashw = mcevth.get_theta();
978	phishw = mcevth.get_phi();
979
980	// calculate vector for shower
981
982	l1 = sin(thetashw)*cos(phishw);
983	m1 = sin(thetashw)*sin(phishw);
984	n1 = cos(thetashw);
985
986	// read the deviation of the telescope with respect to the shower
987
988	mcevth.get_deviations ( &thetaCT, &phiCT );
989
990	if ( (thetaCT == 0.) && (phiCT == 0.) ) {
991
992	// CT was looking to the source (both lines are parallel)
993	// therefore, we calculate the impact parameter as the distance
994	// between the CT axis and the core position
995
996	impactD = dist_r_P( 0., 0., 0., l1, m1, n1, coreX, coreY, 0. );
997
998	} else {
999
1000	// the shower comes off-axis
1001
1002	// obtain with this the final direction of the CT
1003
1004	thetaCT += thetashw;
1005	phiCT += phishw;
1006
1007	// calculate vector for telescope
1008
1009	l2 = sin(thetaCT)*cos(phiCT);
1010	m2 = sin(thetaCT)*sin(phiCT);
1011	n2 = cos(thetaCT);
1012
1013	num = (m1n2coreX - m2n1coreX + l2n1coreY - l1n2coreY);
1014	den = (SQR(l1m2 - l2m1) +
1015	SQR(m1n2 - m2n1) +
1016	SQR(n1l2 - n2l1));
1017	den = sqrt(den);
1018
1019	impactD = fabs(num)/den;
1020
1021	// fprintf(stderr, "[%f %f,%f %f] (%f %f %f) (%f %f %f) %f/%f = ",
1022	// thetashw, phishw, thetaCT, phiCT, l1, m1, n1, l2, m2, n2,
1023	// num, den);
1024
1025	}
1026
1027	// clear camera
1028
1029	for ( i=0; i<ct_NPixels; ++i ){
1030
1031	fnpix[i] = 0.0;
1032	#ifdef __ROOT__
1033	for ( ii=0; ii<(2 * SLICES); ii++ ) {
1034	slices [i][ii] = 0 ;
1035	}
1036	#endif __ROOT__
1037	}
1038
1039	ntcph +=ncph;
1040	ncph = 0;
1041
1042	#ifdef __DETAIL_TRIGGER__
1043	//
1044	// clear Trigger
1045	//
1046
1047	Trigger.Reset() ;
1048	#endif __DETAIL_TRIGGER__
1049
1050	//- - - - - - - - - - - - - - - - - - - - - - - - -
1051	// read photons and "map" them into the pixels
1052	//--------------------------------------------------
1053
1054	// initialize CPhoton
1055
1056	cphoton.fill(0., 0., 0., 0., 0., 0., 0., 0.);
1057
1058	// read the photons data
1059
1060	if ( Data_From_STDIN )
1061	cin.read( flag, SIZE_OF_FLAGS );
1062	else
1063	inputfile.read ( flag, SIZE_OF_FLAGS );
1064
1065	// loop over the photons
1066
1067	t_ini = -99999;
1068
1069	while ( !isA( flag, FLAG_END_OF_EVENT ) ) {
1070
1071	memcpy( (char*)&cphoton, flag, SIZE_OF_FLAGS );
1072
1073	if ( Data_From_STDIN )
1074	cin.read( ((char*)&cphoton)+SIZE_OF_FLAGS, cphoton.mysize()-SIZE_OF_FLAGS );
1075	else
1076	inputfile.read( ((char*)&cphoton)+SIZE_OF_FLAGS, cphoton.mysize()-SIZE_OF_FLAGS );
1077
1078	// increase number of photons
1079
1080	ncph++;
1081
1082	t = cphoton.get_t() ;
1083
1084	if(t_ini == -99999){ // this is the first photon we read from this event
1085	t_ini = t; // memorize time
1086	}
1087
1088	// The photons don't come in chronological order!
1089	// Put the first photon at the center of the array by adding the constant SLICES
1090
1091	t_chan = (int) ((t - t_ini )/ WIDTH_TIMESLICE ) + SLICES ;
1092
1093	if (t_chan > (2 * SLICES)){
1094	log(SIGNATURE, "warning, channel number (%d) exceded limit (%d).\n Setting it to %d .\n",
1095	t_chan, (2 * SLICES), (2 * SLICES));
1096	t_chan = (2 * SLICES);
1097	}
1098	else if(t_chan < 0){
1099	log(SIGNATURE, "warning, channel number (%d) below limit (%d).\n Setting it to %d .\n",
1100	t_chan, 0, 0);
1101	t_chan = 0;
1102	}
1103	/*!@'
1104
1105	@#### Pixelization (for the central pixels).
1106
1107	In order to calculate the coordinates, we use the
1108	change of system described in the documentation
1109	of the source code of \|pixel\_coord.cxx\|.
1110	Then, we will use simply the matrix of change
1111	from one system to the other. In our case, this is:
1112
1113	@[
1114	\begin{bmatrix}X\\Y\\\end{bmatrix}
1115	=
1116	\begin{bmatrix}
1117	1 & \cos(60^\circ)\\
1118	0 & \sin(60^\circ)\\
1119	\end{bmatrix}
1120	\begin{bmatrix}I\\J\\\end{bmatrix}
1121	@]
1122
1123	and hence
1124
1125	@[
1126	\begin{bmatrix}I\\J\\\end{bmatrix}
1127	=
1128	\begin{bmatrix}
1129	1 & -\frac{\cos(60^\circ)}{\sin(60^\circ)}\\
1130	0 &\frac{1}{\sin(60^\circ)}\\
1131	\end{bmatrix}
1132	\begin{bmatrix}X\\Y\\\end{bmatrix}
1133	@]
1134
1135	*/
1136
1137	//+++
1138	// Pixelization
1139	//---
1140
1141	// calculate ij-coordinates
1142
1143	// We use a change of coordinate system, using the following
1144	// matrix of change (m^-1) (this is taken from Mathematica output).
1145	/*
1146	* In[1]:= m={{1,cos60},{0,sin60}}; MatrixForm[m]
1147	*
1148	* Out[1]//MatrixForm= 1 cos60
1149	*
1150	* 0 sin60
1151	*
1152	* In[2]:= inv=Inverse[m]; MatrixForm[inv]
1153	*
1154	* Out[2]//MatrixForm= cos60
1155	* -(-----)
1156	* 1 sin60
1157	*
1158	* 1
1159	* -----
1160	* 0 sin60
1161	*
1162	*/
1163
1164	// go to IJ-coordinate system
1165
1166	cx = cphoton.get_x();
1167	cy = cphoton.get_y();
1168
1169	// get wavelength
1170
1171	last_wl = wl;
1172	wl = cphoton.get_wl();
1173
1174	if ( wl < 1.0 )
1175	break;
1176
1177	if ( (wl > 600.0) \|\| (wl < 290.0) )
1178	break;
1179
1180	// check if photon is inside outermost camera radius
1181
1182	if(sqrt(cxcx + cycy) > (cam.dxc[ct_NPixels-1]+1.5*ct_PixelWidth)){
1183
1184	// read next CPhoton
1185	if ( Data_From_STDIN )
1186	cin.read( flag, SIZE_OF_FLAGS );
1187	else
1188	inputfile.read ( flag, SIZE_OF_FLAGS );
1189
1190	// go to beginning of loop, the photon is lost
1191	continue;
1192
1193	}
1194
1195	// cout << "@#1 " << nshow << ' ' << cx << ' ' << cy << endl;
1196
1197	ci = floor( (cx - cy*COS60/SIN60)/ ct_2Apot + 0.5);
1198	cj = floor( (cy/SIN60) / ct_2Apot + 0.5);
1199
1200	ici = (int)(ci);
1201	icj = (int)(cj);
1202
1203	iici = ici+PIX_ARRAY_HALF_SIDE;
1204	iicj = icj+PIX_ARRAY_HALF_SIDE;
1205
1206	// is it inside the array?
1207
1208	if ( (iici > 0) && (iici < PIX_ARRAY_SIDE) &&
1209	(iicj > 0) && (iicj < PIX_ARRAY_SIDE) ) {
1210
1211	// try to put into pixel
1212
1213	// obtain the pixel number for this photon
1214
1215	nPMT = (int)
1216	pixels[ici+PIX_ARRAY_HALF_SIDE][icj+PIX_ARRAY_HALF_SIDE][PIXNUM];
1217
1218	}
1219	else{
1220
1221	nPMT = -1;
1222
1223	}
1224
1225	// check if outside the central camera
1226
1227	if ( (nPMT < 0) \|\| (nPMT >= ct_NCentralPixels) ) {
1228
1229	// check the outer pixels
1230	nPMT = -1;
1231
1232	for(i=ct_NCentralPixels; i<ct_NPixels; i++){
1233	if( bpoint_is_in_pix( cx, cy, i, &cam) ){
1234	nPMT = i;
1235	break;
1236	}
1237	}
1238
1239	if(nPMT==-1){// the photon is in none of the pixels
1240
1241	// read next CPhoton
1242	if ( Data_From_STDIN )
1243	cin.read( flag, SIZE_OF_FLAGS );
1244	else
1245	inputfile.read ( flag, SIZE_OF_FLAGS );
1246
1247	// go to beginning of loop, the photon is lost
1248	continue;
1249	}
1250
1251	}
1252
1253	#ifdef __QE__
1254
1255	//!@' @#### QE simulation.
1256	//@'
1257
1258	//+++
1259	// QE simulation
1260	//---
1261
1262	// find data point to be used in Lagrange interpolation (-> k)
1263
1264	qeptr = (float **)QE[nPMT];
1265
1266	FindLagrange(qeptr,k,wl);
1267
1268	// if random > quantum efficiency, reject it
1269
1270	qe = Lagrange(qeptr,k,wl) / 100.0;
1271
1272	// fprintf(stdout, "%f\n", qe);
1273
1274	if ( RandomNumber > qe ) {
1275
1276	// read next CPhoton
1277	if ( Data_From_STDIN )
1278	cin.read( flag, SIZE_OF_FLAGS );
1279	else
1280	inputfile.read ( flag, SIZE_OF_FLAGS );
1281
1282	// go to beginning of loop
1283	continue;
1284
1285	}
1286
1287	#endif // __QE__
1288
1289	//+++
1290	// Cphoton is accepted
1291	//---
1292
1293	// increase the number of Cphs. in the PMT, i.e.,
1294	// increase in one unit the counter of the photons
1295	// stored in the pixel nPMT
1296
1297	fnpix[nPMT] += 1.0;
1298
1299	#ifdef __ROOT__
1300	fnslicesum[t_chan] += 1.0 ;
1301	slices[nPMT][t_chan] += 1.0 ;
1302	#endif __ROOT__
1303
1304	#ifdef __DETAIL_TRIGGER__
1305	//
1306	// fill the Trigger class with this phe
1307	//
1308	//
1309
1310	Trigger.Fill( nPMT, ( t - t_ini ) ) ;
1311	#endif __DETAIL_TRIGGER__
1312
1313	// read next CPhoton
1314
1315	if ( Data_From_STDIN )
1316	cin.read( flag, SIZE_OF_FLAGS );
1317	else
1318	inputfile.read ( flag, SIZE_OF_FLAGS );
1319
1320	} // end while, i.e. found end of event
1321
1322	log(SIGNATURE, "End of this event: %d cphs(+%d). . .\n",
1323	ncph, ntcph);
1324
1325	// show number of photons
1326
1327	//cout << ncph << " photons read . . . " << endl << flush;
1328
1329	// skip it ?
1330
1331	for ( i=0; i<nSkip; ++i ) {
1332	if (Skip[i] == (nshow+ntshow)) {
1333	i = -1;
1334	break;
1335	}
1336	}
1337
1338	// if after the previous loop, the exit value of i is -1
1339	// then the shower number is in the list of showers to be
1340	// skipped
1341
1342	if (i == -1) {
1343	log(SIGNATURE, "\t\tskipped!\n");
1344	continue;
1345	}
1346
1347	/*!@'
1348
1349	After reading all the Cherenkov photons for a given event,
1350	we have in the table of number of photons for each pixel
1351	only the 'raw' amount of Cherenkov photons @$n_p@$. Now, we
1352	should take this number as the mean value of the
1353	distribution of photons in that pixel @$p@$, following a
1354	Poisson distribution.
1355
1356	@[ n_p \equiv \mu_p @]
1357
1358	and with this number the amount of light coming from the
1359	shower is calculated @$\hat{n}_p@$.
1360
1361	Then, we calculate the amount of Night Sky Background we
1362	must introduce in that pixel @$p@$. We calculate this using
1363	again a Poisson distribution with mean @$\mu_\mathrm{NSB}@$
1364	(defined in the \|camera.h\| file). The value of
1365	@$\mu_\mathrm{NSB}@$ is obtained from measurements. With this
1366	value, the amount of photons @$\hat{n}_\mathrm{NSB}@$ coming
1367	from the Night Sky Background is calculated.
1368
1369	Finally, the amount of photons for that pixels is:
1370	@[ \hat{n}_p^\mathrm{final} = \hat{n}_p + \hat{n}_\mathrm{NSB} @]
1371
1372	*/
1373
1374	// after reading all the photons, our camera is filled
1375
1376	if ( Select_Energy ) {
1377	if (( mcevth.get_energy() < Select_Energy_le ) \|\|
1378	( mcevth.get_energy() > Select_Energy_ue )) {
1379	log(SIGNATURE, "select_energy: shower rejected.\n");
1380	continue;
1381	}
1382	}
1383
1384	#ifdef __NSB__
1385
1386	//!@' @#### NSB (Night Sky Background) simulation.
1387	//@'
1388
1389	//+++
1390	// NSB simulation
1391	//---
1392
1393	// add NSB "noise"
1394	// TO DO: make meanNSB an array and read the contents from a file!
1395
1396	if ( simulateNSB )
1397	for ( i=0; i<ct_NPixels; ++i )
1398	fnpix[i] += (float)ignpoi( meanNSB );
1399
1400	#endif // __NSB__
1401
1402	// if we should apply any kind of correction, do it here.
1403
1404	for ( i=0; i<ct_NPixels; ++i )
1405	fnpix[i] *= fCorrection;
1406
1407	#ifdef __DETAIL_TRIGGER__
1408	// Trigger.Print() ;
1409	cout << Trigger.Diskriminate() << endl << endl ;
1410	#endif __DETAIL_TRIGGER__
1411
1412	#ifdef __ROOT__
1413
1414	//
1415	// Fill the header of this event
1416	//
1417
1418	Evt->FillHeader ( (UShort_t) (ntshow + nshow) , 20 ) ;
1419
1420	// now put out the data of interest
1421	//
1422	// 1. -> look for the first slice with signal
1423	//
1424
1425	for ( i=0; i<(2 * SLICES); i++ )
1426	if ( fnslicesum[i] > 0. )
1427	break ;
1428
1429	startchan = i ;
1430
1431	//
1432	// copy the slices out of the big array
1433	//
1434	// put the first slice with signal to slice 4
1435	//
1436
1437	for (i=0; i<ct_NPixels; i++ )
1438	for ( ii=(startchan-3); ii < (startchan+12); ii++ )
1439	slices2 [i][ii-startchan+3] = slices [i][ii] ;
1440
1441
1442	//
1443	// if a pixes has a signal put it to the MRawEvt
1444	//
1445
1446	for (i=0; i<ct_NPixels; i++ ) {
1447	if ( fnpix[i] > 0 ) {
1448
1449	for ( ii=0; ii < 15; ii++ ) {
1450	trans [ii] = slices2[i][ii] ;
1451	}
1452
1453	Evt->FillPixel ( (UShort_t) i , trans ) ;
1454
1455	}
1456	}
1457
1458	//
1459	//
1460	//
1461
1462	McEvt->Fill( (UShort_t) mcevth.get_primary() ,
1463	mcevth.get_energy(),
1464	mcevth.get_theta(),
1465	mcevth.get_phi(),
1466	mcevth.get_core(),
1467	mcevth.get_coreX(),
1468	mcevth.get_coreY(),
1469	flli,
1470	ulli, ulli, ulli, ulli, ulli ) ;
1471
1472	//
1473	// write it out to the file outfile
1474	//
1475
1476	EvtTree.Fill() ;
1477
1478	// clear all
1479
1480	Evt->Clear() ;
1481	McEvt->Clear() ;
1482
1483	#endif __ROOT__
1484
1485	//++++++++++++++++++++++++++++++++++++++++++++++++++
1486	// at this point we have a camera full of
1487	// ph.e.s
1488	// we should first apply the trigger condition,
1489	// and if there's trigger, then clean the image,
1490	// calculate the islands statistics and the
1491	// other parameters of the image (Hillas' parameters
1492	// and so on).
1493	//--------------------------------------------------
1494
1495	#ifdef __DEBUG__
1496	printf("\n");
1497
1498	for ( ici=0; ici<PIX_ARRAY_SIDE; ++ici ) {
1499
1500	for ( icj=0; icj<PIX_ARRAY_SIDE; ++icj ) {
1501
1502	if ( (int)pixels[ici][icj][PIXNUM] > -1 ) {
1503
1504	if ( fnpix[(int)pixels[ici][icj][PIXNUM]] > 0. ) {
1505
1506	printf ("@@ %4d %4d %10f %10f %4f (%4d %4d)\n", nshow,
1507	(int)pixels[ici][icj][PIXNUM],
1508	pixels[ici][icj][PIXX],
1509	pixels[ici][icj][PIXY],
1510	fnpix[(int)pixels[ici][icj][PIXNUM]], ici, icj);
1511
1512	}
1513
1514	}
1515
1516	}
1517
1518	}
1519
1520	for (i=0; i<ct_NPixels; ++i) {
1521	printf("%d (%d): ", i, npixneig[i]);
1522	for (j=0; j<npixneig[i]; ++i)
1523	printf(" %d", pixneig[i][j]);
1524	printf("\n");
1525	}
1526
1527	#endif // __DEBUG__
1528
1529	#ifdef __TRIGGER__
1530
1531	/*!@'
1532
1533	@#### Trigger logic simulation.
1534
1535	In the following block we look at the pixel contents, looking
1536	for pixels fulfilling the trigger condition. This condition,
1537	in this current version of the program, is the following:
1538
1539	@itemize
1540
1541	@- \|CT1\|: Two neighbour pixels with charge above the threshold
1542	@$q_0@$. For the old CT1 data, however, the trigger condition
1543	was 'any two pixels with charge above the threshold @$q_0@$'.
1544
1545	@- \|MAGIC\|: A 'closed-packet' of four neighbour pixels, each
1546	of them with charge above the threshold @$q_0@$.
1547
1548	@enditemize
1549
1550	In the following figure you can find a sort of description
1551	about the meanning of 'closed-packet'.
1552
1553	@F
1554
1555	\begin{figure}[htbp]
1556	\begin{center}
1557	\includegraphics{closepck.eps}
1558	\caption{Meanning of the expression ``{\it close-packet}''}
1559	\label{fig:closepacket}
1560	\end{center}
1561	\end{figure}
1562
1563	@F
1564
1565	*/
1566
1567	//++
1568	// TRIGGER Condition
1569	//--
1570
1571	//@ If the input parameter "threshold" is 0 we find the maximum
1572	//@ trigger threshold this event can pass
1573
1574	for(k=0; ( qThreshold == 0 ? (k <= iMAX_THRESHOLD_PHE) : (k == 1) ); k++){
1575
1576	// is there trigger?
1577
1578	noverq0 = 0;
1579	q0 = ( qThreshold == 0. ? (float) k : qThreshold );
1580	trigger = FALSE;
1581	mxgrp = 0;
1582	maxcharge = 0.0;
1583
1584	// Warning! NOT all the camera is able to give trigger
1585	// only up to 'degTrigger' degrees
1586
1587	for ( i=0 ; (i<ct_NCentralPixels) && (trigger==FALSE) ; ++i ) {
1588
1589	// calculate absolute maximum
1590
1591	maxcharge = MAX(fnpix[i],maxcharge);
1592
1593	// is this pixel above threshold ?
1594
1595	if ( fnpix[i] <= q0 )
1596	continue;
1597
1598	// it is: increment the number of pixels above threshold
1599
1600	++noverq0;
1601
1602	// if the trigger already fired, just count the pixels
1603	// above threshold
1604
1605	if ( trigger == TRUE )
1606	continue;
1607
1608	// is this pixel inside the trigger zone in the camera ?
1609
1610	if ( (sqrt(SQR(pixary[i][0]) +
1611	SQR(pixary[i][1]))*plateScale_cm2deg) > degTriggerZone)
1612	continue;
1613
1614	// 'ngrpq0' is the number of neighbours of pixel i with q > q0
1615
1616	ngrpq0 = 0;
1617
1618	// look at each pixel in the neighborhood, and count
1619	// those above threshold q0
1620
1621	for ( j=0 ; j<npixneig[i] && pixneig[i][j]>-1 ; ++j )
1622	if ( fnpix[pixneig[i][j]] > q0 )
1623	++ngrpq0;
1624
1625	// check whether we have trigger
1626
1627	if ( ct_Type == 0 ) {
1628
1629	//++ >>>>> CT1 <<<<<
1630
1631	#ifdef __CT1_NO_NEIGHBOURS__
1632
1633	if ( noverq0 > 1 )
1634	trigger = TRUE;
1635
1636	#else
1637
1638	if ( ngrpq0 > 0 )
1639	trigger = TRUE;
1640
1641	#endif
1642
1643	//-- >>>>> CT1 <<<<<
1644
1645	} else {
1646
1647	//++ >>>>> MAGIC <<<<<
1648
1649	// (at least 4 packed with at least q0 phes)
1650
1651	// there are 3 cases
1652	// 1. zero, one or two neighbours have enough charge: no trigger
1653	// 2. five or six neighbours with enough charge: trigger? sure!!
1654	// 3. three or four neighbours with q > q0 : we must look
1655	// for 'closeness'.
1656
1657	switch ( ngrpq0 ) {
1658
1659	case 0:
1660	case 1:
1661	case 2:
1662
1663	trigger = FALSE;
1664	break;
1665
1666	case 3:
1667	case 4:
1668
1669	// if reaches this line, it means three or four neighbours
1670	// around the central pixel
1671
1672	triggerBits = 1;
1673
1674	for ( j=0 ; j<npixneig[i] && pixneig[i][j]>-1; ++j ) {
1675
1676	if ( fnpix[pixneig[i][j]] > q0 ) {
1677
1678	if ( pixary[pixneig[i][j]][0] > pixary[i][0] ) {
1679
1680	if ( nint(pixary[pixneig[i][j]][1]*10.0) >
1681	nint(pixary[i][1]*10.0) )
1682	bit = 2;
1683	else if ( nint(pixary[pixneig[i][j]][1]*10.0) <
1684	nint(pixary[i][1]*10.0) )
1685	bit = 6;
1686	else
1687	bit = 1;
1688
1689	} else {
1690
1691	if ( nint(pixary[pixneig[i][j]][1]*10.0) >
1692	nint(pixary[i][1]*10.0) )
1693	bit = 3;
1694	else if ( nint(pixary[pixneig[i][j]][1]*10.0) <
1695	nint(pixary[i][1]*10.0) )
1696	bit = 5;
1697	else
1698	bit = 4;
1699
1700	}
1701
1702	triggerBits \|= (1<<bit);
1703
1704	}
1705
1706	}
1707
1708	if ( ngrpq0 == 3 ) { // 4-fold trigger
1709
1710	switch ( triggerBits ) {
1711
1712	case 0x0f: // 0 000111 1
1713	case 0x1d: // 0 001110 1
1714	case 0x39: // 0 011100 1
1715	case 0x71: // 0 111000 1
1716	case 0x63: // 0 110001 1
1717	case 0x47: // 0 100011 1
1718
1719	trigger = TRUE;
1720	break;
1721
1722	default:
1723
1724	trigger = FALSE;
1725
1726	}
1727
1728	} else { // 4-fold trigger
1729
1730	switch ( triggerBits ) {
1731
1732	case 0x1f: // 0 001111 1
1733	case 0x3d: // 0 011110 1
1734	case 0x79: // 0 111100 1
1735	case 0x73: // 0 111001 1
1736	case 0x67: // 0 110011 1
1737	case 0x4f: // 0 100111 1
1738
1739	trigger = TRUE;
1740	break;
1741
1742	default:
1743
1744	trigger = FALSE;
1745
1746	}
1747
1748	}
1749
1750	mxgrp = MAX(ngrpq0,mxgrp);
1751
1752	break;
1753
1754	case 5:
1755	case 6:
1756
1757	trigger = TRUE;
1758	break;
1759
1760	default:
1761
1762	trigger = FALSE;
1763	error( SIGNATURE, "Number of neighbours > 6 !!! Exiting.\n\n");
1764	break;
1765
1766	} // switch (ngrpq0)
1767
1768	} // ct_Type
1769
1770	} // for each pixel i
1771
1772	if ( trigger == FALSE ) {
1773	break;
1774	} // end if
1775
1776	} //end for each threshold
1777	maxtrigthr_phe = (float) k-1; // i.e. maxtrigthr_phe < 0. means that
1778	// the event doesn't even pass threshold 0.
1779	// maxtrigthr_phe >= 0 means, the event passes some threshold
1780	// or (in case the input parameter "threshold" was > 0), the event
1781	// passes the threshold given by the input parameter.
1782	if ( maxtrigthr_phe >= 0. ) {
1783	trigger = TRUE;
1784	}
1785
1786
1787	novq0 = noverq0;
1788
1789	if ( trigger == TRUE ) {
1790
1791	itrigger = i;
1792	++ntrigger;
1793
1794	memcpy( fnpixclean, fnpix, sizeof(float) * ct_NPixels );
1795
1796	#ifdef __TAILCUT__
1797
1798	//!@' @#### Tail-cut condition.
1799	//@'
1800
1801	//++
1802	// tail-cut
1803	//--
1804
1805	// Tail-Cut = 0 : No Tail-Cut
1806	// Tail-Cut > 0 : Make Tail-Cut
1807	// Tail-Cut < 0 : Make Tail-Cut with t_0 = Sqrt[ maximum ]
1808
1809	if (qTailCut > 0.0) {
1810
1811	for ( i=0; i<ct_NPixels; ++i )
1812	if ( fnpix[i] < qTailCut )
1813	fnpixclean[i] = 0.0;
1814
1815	} else if (qTailCut < 0.0) {
1816
1817	maxcharge = sqrt(maxcharge);
1818	for ( i=0; i<ct_NPixels; ++i )
1819	if ( fnpix[i] < maxcharge )
1820	fnpixclean[i] = 0.0;
1821
1822	}
1823
1824	#endif // __TAILCUT__
1825
1826	#ifdef __ISLANDS__
1827
1828	//!@' @#### Islands algorithm.
1829	//@'
1830
1831	//++
1832	// islands counting, and cleanning
1833	//--
1834
1835	if ( countIslands )
1836	do_islands( ct_NPixels, fnpixclean, pixneig, npixneig,
1837	countIslands, nIslandsCut);
1838
1839	#endif // __ISLANDS__
1840
1841	#ifdef __MOMENTS__
1842
1843	//!@' @#### Calculation of parameters of the image.
1844	//@'
1845
1846	//++
1847	// moments calculation
1848	//--
1849
1850	// calculate moments and other things
1851
1852	moments_ptr = moments( anaPixels, fnpixclean, pixary,
1853	plateScale_cm2deg, 0 );
1854
1855	charge = moments_ptr->charge ;
1856	smax = moments_ptr->smax ;
1857	maxs = moments_ptr->maxs ;
1858	nmaxs = moments_ptr->nmaxs ;
1859	length = moments_ptr->length ;
1860	width = moments_ptr->width ;
1861	dist = moments_ptr->dist ;
1862	xdist = moments_ptr->xdist ;
1863	azw = moments_ptr->azw ;
1864	miss = moments_ptr->miss ;
1865	alpha = moments_ptr->alpha ;
1866	conc = moments_ptr->conc ;
1867	asymx = moments_ptr->asymx ;
1868	asymx = moments_ptr->asymx ;
1869	phiasym= moments_ptr->phi;
1870
1871	lenwid_ptr = lenwid( anaPixels, fnpixclean, pixary,
1872	plateScale_cm2deg,
1873	ct_PixelWidth_corner_2_corner_half);
1874
1875
1876	// fill the diagnostic Tree
1877
1878	event = new MDiagEventobject();
1879
1880	i=0;
1881	image_data[i] = event->n = hidt/10; i++;
1882	image_data[i] = event->primary = mcevth.get_primary(); i++;
1883	image_data[i] = event->energy = mcevth.get_energy(); i++;
1884	image_data[i] = event->cored = coreD; i++;
1885	image_data[i] = event->impact = impactD; i++;
1886	image_data[i] = event->xcore = coreX; i++;
1887	image_data[i] = event->ycore = coreY; i++;
1888	image_data[i] = event->theta = mcevth.get_theta(); i++;
1889	image_data[i] = event->phi = mcevth.get_phi(); i++;
1890	image_data[i] = event->deviations = mcevth.get_deviations (&dtheta, &dphi); i++;
1891	image_data[i] = event->dtheta = dtheta; i++;
1892	image_data[i] = event->dphi = dphi; i++;
1893	image_data[i] = event->trigger = trigger; i++;
1894	image_data[i] = event->ncphs = ncph; i++;
1895	image_data[i] = event->maxpassthr_phe = maxtrigthr_phe; i++;
1896	image_data[i] = event->nphes = charge; i++;
1897	image_data[i] = event->nphes2 = smax; i++;
1898	image_data[i] = event->length = length; i++;
1899	image_data[i] = event->width = width; i++;
1900	image_data[i] = event->dist = dist; i++;
1901	image_data[i] = event->xdist = xdist; i++;
1902	image_data[i] = event->azw = azw; i++;
1903	image_data[i] = event->miss = miss; i++;
1904	image_data[i] = event->alpha = alpha; i++;
1905	image_data[i] = event->conc2 = conc[0]; i++;
1906	image_data[i] = event->conc3 = conc[1]; i++;
1907	image_data[i] = event->conc4 = conc[2]; i++;
1908	image_data[i] = event->conc5 = conc[3]; i++;
1909	image_data[i] = event->conc6 = conc[4]; i++;
1910	image_data[i] = event->conc7 = conc[5]; i++;
1911	image_data[i] = event->conc8 = conc[6]; i++;
1912	image_data[i] = event->conc9 = conc[7]; i++;
1913	image_data[i] = event->conc10 = conc[8]; i++;
1914	image_data[i] = event->asymx = asymx; i++;
1915	image_data[i] = event->asymy = asymy; i++;
1916	image_data[i] = event->phiasym = phiasym; i++;
1917
1918	// there should be "nvar" variables
1919
1920	if ( i != nvar )
1921	error( SIGNATURE, "Wrong entry number for diagnostic data.\n" );
1922
1923	tree->Fill();
1924	delete event;
1925
1926	// put information in the data file,
1927
1928	datafile << ntrigger;
1929	for(i=0;i<nvar;i++) {
1930	datafile << ' ' << image_data[i];
1931	}
1932
1933
1934	#endif // __MOMENTS__
1935
1936
1937	// revert the fnpixclean matrix into fnpix
1938	// (now we do this, but maybe in a future we want to
1939	// use both fnpix and fnpixclean for different things
1940
1941	memcpy( fnpix, fnpixclean, sizeof(float) * ct_NPixels );
1942
1943	// put this information in the data file,
1944
1945	if ( Write_All_Data ) {
1946	datafile << ' ' << -9999;
1947	for ( i=0; i<ct_NPixels; ++i )
1948	datafile << ' ' << fnpix[i];
1949	}
1950
1951	datafile << endl;
1952
1953	mcevth.set_trigger( TRUE );
1954
1955	log(SIGNATURE, "TRIGGER\n");
1956
1957	} else { // ( trigger == FALSE )
1958
1959	event = new MDiagEventobject();
1960
1961	i=0;
1962	image_data[i] = event->n = hidt/10; i++;
1963	image_data[i] = event->primary = mcevth.get_primary(); i++;
1964	image_data[i] = event->energy = mcevth.get_energy(); i++;
1965	image_data[i] = event->cored = coreD = mcevth.get_core(&coreX, &coreY); i++;
1966	image_data[i] = event->impact = coreD; i++;
1967	image_data[i] = event->xcore = coreX; i++;
1968	image_data[i] = event->ycore = coreY; i++;
1969	image_data[i] = event->theta = mcevth.get_theta(); i++;
1970	image_data[i] = event->phi = mcevth.get_phi(); i++;
1971	image_data[i] = event->deviations = mcevth.get_deviations(&dtheta, &dphi); i++;
1972	image_data[i] = event->dtheta = dtheta; i++;
1973	image_data[i] = event->dphi = dphi; i++;
1974	image_data[i] = event->trigger = trigger; i++;
1975	image_data[i] = event->ncphs = ncph; i++;
1976	image_data[i] = event->maxpassthr_phe = maxtrigthr_phe; i++;
1977	image_data[i] = -1.; i++;
1978	image_data[i] = -1.; i++;
1979	image_data[i] = -1.; i++;
1980	image_data[i] = -1.; i++;
1981	image_data[i] = -1.; i++;
1982	image_data[i] = -1.; i++;
1983	image_data[i] = -1.; i++;
1984	image_data[i] = -1.; i++;
1985	image_data[i] = -1.; i++;
1986	image_data[i] = -1.; i++;
1987	image_data[i] = -1.; i++;
1988	image_data[i] = -1.; i++;
1989	image_data[i] = -1.; i++;
1990	image_data[i] = -1.; i++;
1991	image_data[i] = -1.; i++;
1992	image_data[i] = -1.; i++;
1993	image_data[i] = -1.; i++;
1994	image_data[i] = -1.; i++;
1995	image_data[i] = -1.; i++;
1996	image_data[i] = -1.; i++;
1997	image_data[i] = -1.; i++;
1998
1999	// there should be "nvar" variables
2000
2001	if ( i != nvar )
2002	error( SIGNATURE, "Wrong entry length for Ntuple.\n" );
2003
2004	tree->Fill();
2005	delete event;
2006
2007	// put this information in the data file,
2008
2009	if ( Write_All_Data ) {
2010
2011	datafile << ntrigger;
2012	for ( i=0; i<nvar; ++i )
2013	datafile << ' ' << image_data[i];
2014
2015	datafile << -9999;
2016	for ( i=0; i<ct_NPixels; ++i )
2017	datafile << ' ' << fnpix[i];
2018
2019	datafile << endl;
2020	}
2021
2022	mcevth.set_trigger( FALSE );
2023
2024	} // trigger == FALSE
2025
2026	#endif // __TRIGGER__
2027
2028	//!@' @#### Save data.
2029	//@'
2030
2031	//++++++++++++++++++++++++++++++++++++++++++++++++++
2032	// we now have all information we want
2033	// the only thing we must do now is writing it to
2034	// the output file
2035	//--------------------------------------------------
2036
2037	//++
2038	// save the image to the file
2039	//--
2040
2041	// write MCEventHeader to output file
2042
2043	outputfile.write( (char *)&mcevth, mcevth.mysize() );
2044
2045	#ifdef __TRIGGER__
2046
2047	// save the image
2048
2049	if ( (trigger == TRUE) \|\| (Write_All_Images == TRUE) )
2050	outputfile.write( (char ) fnpix, ct_NPixels sizeof( float ) );
2051
2052	#else
2053
2054	// save the image
2055
2056	outputfile.write( (char ) fnpix, ct_NPixels sizeof( float ) );
2057
2058	#endif // __TRIGGER__
2059
2060	if ( Data_From_STDIN )
2061	cin.read( flag, SIZE_OF_FLAGS );
2062	else
2063	inputfile.read ( flag, SIZE_OF_FLAGS );
2064
2065	} // end while there is a next event
2066
2067	if( !isA( flag, FLAG_END_OF_RUN )){
2068	error( SIGNATURE, "Expected end of run flag, but found: %s\n", flag );
2069	}
2070	else { // found end of run
2071	ntshow += nshow;
2072	log(SIGNATURE, "End of this run with %d events . . .\n", nshow);
2073
2074	if ( Data_From_STDIN )
2075	cin.read( flag, SIZE_OF_FLAGS );
2076	else
2077	inputfile.read ( flag, SIZE_OF_FLAGS );
2078
2079	if( isA( flag, FLAG_END_OF_FILE ) ){ // end of file
2080	log(SIGNATURE, "End of file . . .\n");
2081	still_in_loop = FALSE;
2082
2083	if ((! Data_From_STDIN) && (! inputfile.eof())){
2084
2085	// we have concatenated input files.
2086	// get signature of the next part and check it.
2087	// NOTE: this part repeats further up in the code;
2088	// if you change something here you probably want to change it
2089	// there as well
2090
2091	strcpy(Signature, REFL_SIGNATURE);
2092
2093	strcpy(sign, Signature);
2094
2095	inputfile.read( (char *)sign, strlen(Signature));
2096
2097	if (strcmp(sign, Signature) != 0) {
2098	cerr << "ERROR: Signature of .rfl file is not correct\n";
2099	cerr << '"' << sign << '"' << '\n';
2100	cerr << "should be: " << Signature << '\n';
2101	exit(1);
2102	}
2103
2104	if ( Data_From_STDIN )
2105	cin.read( (char *)sign, 1);
2106	else
2107	inputfile.read( (char *)sign, 1);
2108
2109	}
2110
2111	} // end if found end of file
2112	} // end if found end of run
2113	if ( Data_From_STDIN )
2114	cin.read( flag, SIZE_OF_FLAGS );
2115	else
2116	inputfile.read ( flag, SIZE_OF_FLAGS );
2117	} // end if else found start of run
2118	} // end big while loop
2119
2120	//!@' @#### End of program.
2121	//@'
2122
2123	//end my version
2124
2125	#ifdef __ROOT__
2126	//++
2127	// put the Event to the root file
2128	//--
2129
2130	EvtTree.Write() ;
2131	outfile.Write() ;
2132	outfile.Close() ;
2133
2134	#endif __ROOT__
2135
2136	// close input file
2137
2138	ntcph += ncph;
2139	log( SIGNATURE, "%d event(s), with a total of %d C.photons\n",
2140	ntshow, ntcph );
2141	log( SIGNATURE, "Fraction of triggers: %5.1f%% (%d out of %d)\n",
2142	((float)ntrigger) / ((float)ntshow) * 100.0, ntrigger, ntshow);
2143
2144	// close files
2145
2146	log( SIGNATURE, "Closing files\n" );
2147
2148	inputfile.close();
2149	outputfile.close();
2150	datafile.close();
2151
2152	hfile->Write();
2153
2154	hfile->Close();
2155
2156	#ifdef __DETAIL_TRIGGER__
2157	// Output of Trigger statistics
2158	//
2159
2160	Trigger.PrintStat() ;
2161	#endif __DETAIL_TRIGGER__
2162
2163	// program finished
2164
2165	log( SIGNATURE, "Done.\n");
2166
2167	return( 0 );
2168	}
2169	//!@}
2170
2171	// @T \newpage
2172
2173	//!@subsection Functions definition.
2174
2175	//!-----------------------------------------------------------
2176	// @name present
2177	//
2178	// @desc Make some presentation
2179	//
2180	// @date Sat Jun 27 05:58:56 MET DST 1998
2181	//------------------------------------------------------------
2182	// @function
2183
2184	//!@{
2185	void
2186	present(void)
2187	{
2188	cout << "##################################################\n"
2189	<< SIGNATURE << '\n' << '\n'
2190	<< "Processor of the reflector output\n"
2191	<< "J C Gonzalez, Jun 1998\n"
2192	<< "##################################################\n\n"
2193	<< flush ;
2194	}
2195	//!@}
2196
2197
2198	//!-----------------------------------------------------------
2199	// @name usage
2200	//
2201	// @desc show help
2202	//
2203	// @date Tue Dec 15 16:23:30 MET 1998
2204	//------------------------------------------------------------
2205	// @function
2206
2207	//!@{
2208	void
2209	usage(void)
2210	{
2211	present();
2212	cout << "\nusage ::\n\n"
2213	<< "\t camera "
2214	<< " [ -@ paramfile ] "
2215	<< " [ -h ] "
2216	<< "\n\n or \n\n"
2217	<< "\t camera < paramfile"
2218	<< "\n\n";
2219	exit(0);
2220	}
2221	//!@}
2222
2223
2224	//!-----------------------------------------------------------
2225	// @name log
2226	//
2227	// @desc function to send log information
2228	//
2229	// @var funct Name of the caller function
2230	// @var fmt Format to be used (message)
2231	// @var ... Other information to be shown
2232	//
2233	// @date Sat Jun 27 05:58:56 MET DST 1998
2234	//------------------------------------------------------------
2235	// @function
2236
2237	//!@{
2238	void
2239	log(const char funct, char fmt, ...)
2240	{
2241	va_list args;
2242
2243	// Display the name of the function that called error
2244	printf("[%s]: ", funct);
2245
2246	// Display the remainder of the message
2247	va_start(args, fmt);
2248	vprintf(fmt, args);
2249	va_end(args);
2250	}
2251	//!@}
2252
2253
2254	//!-----------------------------------------------------------
2255	// @name error
2256	//
2257	// @desc function to send an error message, and abort the program
2258	//
2259	// @var funct Name of the caller function
2260	// @var fmt Format to be used (message)
2261	// @var ... Other information to be shown
2262	//
2263	// @date Sat Jun 27 05:58:56 MET DST 1998
2264	//------------------------------------------------------------
2265	// @function
2266
2267	//!@{
2268	void
2269	error(const char funct, char fmt, ...)
2270	{
2271	va_list args;
2272
2273	// Display the name of the function that called error
2274	fprintf(stderr, "ERROR in %s: ", funct);
2275
2276	// Display the remainder of the message
2277	va_start(args, fmt);
2278	vfprintf(stderr, fmt, args);
2279	va_end(args);
2280
2281	perror(funct);
2282
2283	abort();
2284	}
2285	//!@}
2286
2287
2288	//!-----------------------------------------------------------
2289	// @name isA
2290	//
2291	// @desc returns TRUE(FALSE), if the flag is(is not) the given
2292	//
2293	// @var s1 String to be searched
2294	// @var flag Flag to compare with string s1
2295	// @return TRUE: both strings match; FALSE: oth.
2296	//
2297	// @date Wed Jul 8 15:25:39 MET DST 1998
2298	//------------------------------------------------------------
2299	// @function
2300
2301	//!@{
2302	int
2303	isA( char * s1, const char * flag ) {
2304	return ( (strncmp((char *)s1, flag, SIZE_OF_FLAGS)==0) ? 1 : 0 );
2305	}
2306	//!@}
2307
2308
2309	//!-----------------------------------------------------------
2310	// @name read_ct_file
2311	//
2312	// @desc read CT definition file
2313	//
2314	// @date Sat Jun 27 05:58:56 MET DST 1998
2315	//------------------------------------------------------------
2316	// @function
2317
2318	//!@{
2319	void
2320	read_ct_file(void)
2321	{
2322	char line[LINE_MAX_LENGTH]; //@< line to get from the ctin
2323	char token[ITEM_MAX_LENGTH]; //@< a single token
2324	int i, j; //@< dummy counters
2325
2326	log( "read_ct_file", "start.\n" );
2327
2328	ifstream ctin ( ct_filename );
2329
2330	if ( ctin.bad() )
2331	error( "read_ct_file",
2332	"Cannot open CT def. file: %s\n", ct_filename );
2333
2334	// loop till the "end" directive is reached
2335
2336	while (!ctin.eof()) {
2337
2338	// get line from stdin
2339
2340	ctin.getline(line, LINE_MAX_LENGTH);
2341
2342	// look for each item at the beginning of the line
2343
2344	for (i=0; i<=define_mirrors; i++)
2345	if (strstr(line, CT_ITEM_NAMES[i]) == line)
2346	break;
2347
2348	// if it is not a valid line, just ignore it
2349
2350	if (i == define_mirrors+1)
2351	continue;
2352
2353	// case block for each directive
2354
2355	switch ( i ) {
2356
2357	case type: // <type of telescope> (0:CT1 ¦ 1:MAGIC)
2358
2359	// get focal distance
2360
2361	sscanf(line, "%s %d", token, &ct_Type);
2362
2363	log( "read_ct_file", "<Type of Telescope>: %s\n",
2364	((ct_Type==0) ? "CT1" : "MAGIC") );
2365
2366	break;
2367
2368	case focal_distance: // <focal distance> [cm]
2369
2370	// get focal distance
2371
2372	sscanf(line, "%s %f", token, &ct_Focal_mean);
2373
2374	log( "read_ct_file", "<Focal distance>: %f cm\n", ct_Focal_mean );
2375
2376	break;
2377
2378	case focal_std: // s(focal distance) [cm]
2379
2380	// get focal distance
2381
2382	sscanf(line, "%s %f", token, &ct_Focal_std);
2383
2384	log( "read_ct_file", "s(Focal distance): %f cm\n", ct_Focal_std );
2385
2386	break;
2387
2388	case point_spread: // <point spread> [cm]
2389
2390	// get point spread
2391
2392	sscanf(line, "%s %f", token, &ct_PSpread_mean);
2393
2394	log( "read_ct_file", "<Point spread>: %f cm\n", ct_PSpread_mean );
2395
2396	break;
2397
2398	case point_std: // s(point spread) [cm]
2399
2400	// get point spread
2401
2402	sscanf(line, "%s %f", token, &ct_PSpread_std);
2403
2404	log( "read_ct_file", "s(Point spread): %f cm\n", ct_PSpread_std );
2405
2406	break;
2407
2408	case adjustment_dev: // s(adjustment_dev) [cm]
2409
2410	// get point spread
2411
2412	sscanf(line, "%s %f", token, &ct_Adjustment_std);
2413
2414	log( "read_ct_file", "s(Adjustment): %f cm\n", ct_Adjustment_std );
2415
2416	break;
2417
2418	case black_spot: // radius of the black spot in the center [cm]
2419
2420	// get black spot radius
2421
2422	sscanf(line, "%s %f", token, &ct_BlackSpot_rad);
2423
2424	log( "read_ct_file", "Radius of the black spots: %f cm\n",
2425	ct_BlackSpot_rad);
2426
2427	break;
2428
2429	case r_mirror: // radius of the mirrors [cm]
2430
2431	// get radius of mirror
2432
2433	sscanf(line, "%s %f", token, &ct_RMirror);
2434
2435	log( "read_ct_file", "Radii of the mirrors: %f cm\n", ct_RMirror );
2436
2437	break;
2438
2439	case n_mirrors: // number of mirrors
2440
2441	// get the name of the output_file from the line
2442
2443	sscanf(line, "%s %d", token, &ct_NMirrors);
2444
2445	log( "read_ct_file", "Number of mirrors: %d\n", ct_NMirrors );
2446
2447	break;
2448
2449	case camera_width: // camera width [cm]
2450
2451	// get the name of the ct_file from the line
2452
2453	sscanf(line, "%s %f", token, &ct_CameraWidth);
2454
2455	log( "read_ct_file", "Camera width: %f cm\n", ct_CameraWidth );
2456
2457	break;
2458
2459	case n_pixels: // number of pixels
2460
2461	// get the name of the output_file from the line
2462
2463	sscanf(line, "%s %d", token, &ct_NPixels);
2464
2465	log( "read_ct_file", "Number of pixels: %d\n", ct_NPixels );
2466
2467	break;
2468
2469	case n_centralpixels: // number of central pixels
2470
2471	// get the name of the output_file from the line
2472
2473	sscanf(line, "%s %d", token, &ct_NCentralPixels);
2474
2475	log( "read_ct_file", "Number of central pixels: %d\n", ct_NCentralPixels );
2476
2477	break;
2478
2479	case n_gappixels: // number of gap pixels
2480
2481	// get the name of the output_file from the line
2482
2483	sscanf(line, "%s %d", token, &ct_NGapPixels);
2484
2485	log( "read_ct_file", "Number of gap pixels: %d\n", ct_NGapPixels );
2486
2487	break;
2488
2489	case pixel_width: // pixel width [cm]
2490
2491	// get the name of the ct_file from the line
2492
2493	sscanf(line, "%s %f", token, &ct_PixelWidth);
2494
2495	ct_PixelWidth_corner_2_corner = ct_PixelWidth / cos(RAD(30.0));
2496	ct_PixelWidth_corner_2_corner_half =
2497	ct_PixelWidth_corner_2_corner * 0.50;
2498	ct_Apot = ct_PixelWidth / 2;
2499	ct_2Apot = ct_Apot * 2.0;
2500
2501	log( "read_ct_file", "Pixel width: %f cm\n", ct_PixelWidth );
2502
2503	break;
2504
2505	case define_mirrors: // read table with the parameters of the mirrors
2506
2507	log( "read_ct_file", "Table of mirrors data:\n" );
2508
2509	// check whether the number of mirrors was already set
2510
2511	if ( ct_NMirrors == 0 )
2512	error( "read_ct_file", "NMirrors was not set.\n" );
2513
2514	// allocate memory for paths list
2515
2516	log( "read_ct_file", "Allocating memory for ct_data\n" );
2517
2518	ct_data = new float*[ct_NMirrors];
2519
2520	for (i=0; i<ct_NMirrors; i++)
2521	ct_data[i] = new float[CT_NDATA];
2522
2523	// read data
2524
2525	log( "read_ct_file", "Reading mirrors data...\n" );
2526
2527	for (i=0; i<ct_NMirrors; i++)
2528	for (j=0; j<CT_NDATA; j++)
2529	ctin >> ct_data[i][j];
2530
2531	break;
2532
2533	} // switch ( i )
2534
2535	} // end while
2536
2537	// end
2538
2539	log( "read_ct_file", "done.\n" );
2540
2541	return;
2542	}
2543	//!@}
2544
2545
2546	//!-----------------------------------------------------------
2547	// @name read_pixels
2548	//
2549	// @desc read pixels data
2550	//
2551	// @date Fri Mar 12 16:33:34 MET 1999
2552	//------------------------------------------------------------
2553	// @function
2554
2555	//!@{
2556	void
2557	read_pixels(struct camera *pcam)
2558	{
2559	ifstream qefile;
2560	char line[LINE_MAX_LENGTH];
2561	int n, i, j, k;
2562	float qe;
2563
2564	//------------------------------------------------------------
2565	// first, pixels' coordinates
2566
2567	pcam->inumpixels = ct_NPixels;
2568	pcam->inumcentralpixels = ct_NCentralPixels;
2569	pcam->inumgappixels = ct_NGapPixels;
2570	pcam->inumbigpixels = ct_NPixels - ct_NCentralPixels - ct_NGapPixels;
2571	pcam->dpixdiameter_cm = ct_PixelWidth;
2572
2573	// initialize pixel numbers
2574
2575	for ( i=0; i<PIX_ARRAY_SIDE; ++i )
2576	for ( j=0; j<PIX_ARRAY_SIDE; ++j )
2577	pixels[i][j][PIXNUM] = -1;
2578
2579	pixary = new float* [2*ct_NCentralPixels];
2580	for ( i=0; i<2*ct_NCentralPixels; ++i )
2581	pixary[i] = new float[2];
2582
2583	pixneig = new int* [ct_NCentralPixels];
2584	for ( i=0; i<ct_NCentralPixels; ++i ) {
2585	pixneig[i] = new int[6];
2586	for ( j=0; j<6; ++j )
2587	pixneig[i][j] = -1;
2588	}
2589
2590	npixneig = new int[ct_NCentralPixels];
2591	for ( i=0; i<ct_NCentralPixels; ++i )
2592	npixneig[i] = 0;
2593
2594	// generate all coordinates
2595
2596	igen_pixel_coordinates(pcam);
2597
2598	// transfer coordinates to the working arrays for
2599	// the central pixels
2600
2601	for(k=0; k<ct_NCentralPixels; k++){
2602
2603	i = (int) pcam->di[k];
2604	j = (int) pcam->dj[k];
2605
2606	pixels[i+PIX_ARRAY_HALF_SIDE][j+PIX_ARRAY_HALF_SIDE][PIXNUM] = k;
2607	pixels[i+PIX_ARRAY_HALF_SIDE][j+PIX_ARRAY_HALF_SIDE][PIXX] = pcam->dxc[k];
2608	pixels[i+PIX_ARRAY_HALF_SIDE][j+PIX_ARRAY_HALF_SIDE][PIXY] = pcam->dyc[k];
2609
2610	pixary[k][0] = pcam->dxc[k];
2611	pixary[k][1] = pcam->dyc[k];
2612	}
2613
2614	// calculate tables of neighbours
2615
2616	#ifdef __DEBUG__
2617	for ( n=0 ; n<ct_NPixels ; ++n ) {
2618	cout << "Para el pixel " << n << ": ";
2619	for ( i=n+1 ; (i<ct_NPixels)&&(npixneig[n]<6) ; ++i) {
2620	if ( pixels_are_neig(n,i) == TRUE ) {
2621	pixneig[n][npixneig[n]] = i;
2622	pixneig[i][npixneig[i]] = n;
2623	cout << i << ' ';
2624	++npixneig[n];
2625	++npixneig[i];
2626	}
2627	}
2628	cout << endl << flush;
2629	}
2630	#else // ! __DEBUG__
2631	for ( n=0 ; n<ct_NCentralPixels ; ++n )
2632	for ( i=n+1 ; (i<ct_NCentralPixels)&&(npixneig[n]<6) ; ++i)
2633	if ( pixels_are_neig(n,i) == TRUE ) {
2634	pixneig[n][npixneig[n]] = i;
2635	pixneig[i][npixneig[i]] = n;
2636	++npixneig[n];
2637	++npixneig[i];
2638	}
2639	#endif // ! __DEBUG__
2640
2641	#ifdef __DEBUG__
2642	for ( n=0 ; n<ct_NPixels ; ++n ) {
2643	cout << n << ':';
2644	for ( j=0; j<npixneig[n]; ++j)
2645	cout << ' ' << pixneig[n][j];
2646	cout << endl << flush;
2647	}
2648	#endif // __DEBUG__
2649
2650	//------------------------------------------------------------
2651	// second, pixels' QE
2652
2653	// try to open the file
2654
2655	log("read_pixels", "Openning the file \"%s\" . . .\n", QE_FILE);
2656
2657	qefile.open( QE_FILE );
2658
2659	// if it is wrong or does not exist, exit
2660
2661	if ( qefile.bad() )
2662	error( "read_pixels", "Cannot open \"%s\". Exiting.\n", QE_FILE );
2663
2664	// read file
2665
2666	log("read_pixels", "Reading data . . .\n");
2667
2668	i=-1;
2669
2670	while ( ! qefile.eof() ) {
2671
2672	// get line from the file
2673
2674	qefile.getline(line, LINE_MAX_LENGTH);
2675
2676	// skip if comment
2677
2678	if ( *line == '#' )
2679	continue;
2680
2681	// if it is the first valid value, it is the number of QE data points
2682
2683	if ( i < 0 ) {
2684
2685	// get the number of datapoints
2686
2687	sscanf(line, "%d", &pointsQE);
2688
2689	// allocate memory for the table of QEs
2690
2691	QE = new float ** [ct_NPixels];
2692
2693	for ( i=0; i<ct_NPixels; ++i ) {
2694	QE[i] = new float * [2];
2695	QE[i][0] = new float[pointsQE];
2696	QE[i][1] = new float[pointsQE];
2697	}
2698
2699	QElambda = new float [pointsQE];
2700
2701	for ( i=0; i<pointsQE; ++i ) {
2702	qefile.getline(line, LINE_MAX_LENGTH);
2703	sscanf(line, "%f", &QElambda[i]);
2704	}
2705
2706	i=0;
2707
2708	continue;
2709	}
2710
2711	// get the values (num-pixel, num-datapoint, QE-value)
2712
2713	sscanf(line, "%d %d %f", &i, &j, &qe);
2714
2715	if ( ((i-1) < ct_NPixels) && ((i-1) > -1) &&
2716	((j-1) < pointsQE) && ((j-1) > -1) ) {
2717	QE[i-1][0][j-1] = QElambda[j-1];
2718	QE[i-1][1][j-1] = qe;
2719	}
2720
2721	}
2722
2723	// close file
2724
2725	qefile.close();
2726
2727	// end
2728
2729	log("read_pixels", "Done.\n");
2730
2731	}
2732	//!@}
2733
2734
2735	//!-----------------------------------------------------------
2736	// @name pixels_are_neig
2737	//
2738	// @desc check whether two pixels are neighbours
2739	//
2740	// @var pix1 Number of the first pixel
2741	// @var pix2 Number of the second pixel
2742	// @return TRUE: both pixels are neighbours; FALSE: oth.
2743	//
2744	// @date Wed Sep 9 17:58:37 MET DST 1998
2745	//------------------------------------------------------------
2746	// @function
2747
2748	//!@{
2749	int
2750	pixels_are_neig(int pix1, int pix2)
2751	{
2752	if ( sqrt(SQR( pixary[pix1][0] - pixary[pix2][0] ) +
2753	SQR( pixary[pix1][1] - pixary[pix2][1] ) )
2754	> ct_PixelWidth_corner_2_corner )
2755	return ( FALSE );
2756	else
2757	return ( TRUE );
2758	}
2759	//!@}
2760
2761	//!-----------------------------------------------------------
2762	// @name igen_pixel_coordinates
2763	//
2764	// @desc generate the pixel center coordinates
2765	//
2766	// @var *pcam structure camera containing all the
2767	// camera information
2768	// @return total number of pixels
2769	//
2770	// DP
2771	//
2772	// @date Thu Oct 14 10:41:03 CEST 1999
2773	//------------------------------------------------------------
2774	// @function
2775
2776	//!@{
2777	/****** igen_pixel_coordinates() *******************************/
2778
2779	int igen_pixel_coordinates(struct camera *pcam) {
2780	/* generate pixel coordinates, return value is number of pixels */
2781
2782	int i, itot_inside_ring, iN, in, ipixno, iring_no, ipix_in_ring, isegment;
2783	float fsegment_fract;
2784	double dtsize;
2785	double dhsize;
2786	double dpsize;
2787	double dxfirst_pix;
2788	double dyfirst_pix;
2789	double ddxseg1, ddxseg2, ddxseg3, ddxseg4, ddxseg5, ddxseg6;
2790	double ddyseg1, ddyseg2, ddyseg3, ddyseg4, ddyseg5, ddyseg6;
2791
2792
2793	double dstartx, dstarty; /* for the gap pixels and outer pixels */
2794	int j, nrow;
2795
2796	dpsize = pcam->dpixdiameter_cm;
2797	dtsize = dpsize * sqrt(3.) / 2.;
2798	dhsize = dpsize / 2.;
2799
2800	/* Loop over central pixels to generate co-ordinates */
2801
2802	for(ipixno=1; ipixno <= pcam->inumcentralpixels; ipixno++){
2803
2804	/* Initialise variables. The central pixel = ipixno 1 in ring iring_no 0 */
2805
2806	pcam->dpixsizefactor[ipixno] = 1.;
2807
2808	in = 0;
2809
2810	i = 0;
2811	itot_inside_ring = 0;
2812	iring_no = 0;
2813
2814	/* Calculate the number of pixels out to and including the ring containing pixel number */
2815	/* ipixno e.g. for pixel number 17 in ring number 2 itot_inside_ring = 19 */
2816
2817	while (itot_inside_ring == 0){
2818
2819	iN = 3(i(i+1)) + 1;
2820
2821	if (ipixno <= iN){
2822	iring_no = i;
2823	itot_inside_ring = iN;
2824	}
2825
2826	i++;
2827	}
2828
2829
2830	/* Find the number of pixels which make up ring number iring_no e.g. ipix_in_ring = 6 for ring 1 */
2831
2832	ipix_in_ring = 0;
2833	for (i = 0; i < iring_no; ++i){
2834
2835	ipix_in_ring = ipix_in_ring + 6;
2836	}
2837
2838	/* The camera is viewed as 6 radial segments ("pie slices"). Knowing the number of pixels in its */
2839	/* ring calculate which segment the pixel ipixno is in. Then find how far across this segment it is */
2840	/* as a fraction of the number of pixels in this sixth of the ring (ask SMB). */
2841
2842	isegment = 0;
2843	fsegment_fract = 0.;
2844	if (iring_no > 0) {
2845
2846	isegment = (int)((ipixno - itot_inside_ring + ipix_in_ring - 0.5) / iring_no + 1); /* integer division ! numbering starts at 1 */
2847
2848	fsegment_fract = (ipixno - (itot_inside_ring - ipix_in_ring)) - ((isegment-1)*iring_no) - 1 ;
2849
2850	}
2851
2852	/* the first pixel in each ring lies on the positive x axis at a distance dxfirst_pix = iring_no * the */
2853	/* pixel width (flat to flat) dpsize. */
2854
2855	dxfirst_pix = dpsize*iring_no;
2856	dyfirst_pix = 0.;
2857
2858	/* the vector between the first and last pixels in a segment n is (ddxsegn, ddysegn) */
2859
2860	ddxseg1 = - dhsize*iring_no;
2861	ddyseg1 = dtsize*iring_no;
2862	ddxseg2 = -dpsize*iring_no;
2863	ddyseg2 = 0.;
2864	ddxseg3 = ddxseg1;
2865	ddyseg3 = -ddyseg1;
2866	ddxseg4 = -ddxseg1;
2867	ddyseg4 = -ddyseg1;
2868	ddxseg5 = -ddxseg2;
2869	ddyseg5 = 0.;
2870	ddxseg6 = -ddxseg1;
2871	ddyseg6 = ddyseg1;
2872
2873	/* to find the position of pixel ipixno take the position of the first pixel in the ring and move */
2874	/* anti-clockwise around the ring by adding the segment to segment vectors. */
2875
2876	switch (isegment) {
2877
2878	case 0:
2879
2880	pcam->dxc[ipixno-1] = 0.;
2881	pcam->dyc[ipixno-1] = 0.;
2882
2883	case 1:
2884	pcam->dxc[ipixno-1] = dxfirst_pix - dhsize*fsegment_fract;
2885	pcam->dyc[ipixno-1] = dyfirst_pix + dtsize*fsegment_fract;
2886
2887	break;
2888
2889	case 2:
2890
2891	pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 - dpsize*fsegment_fract;
2892	pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + 0.;
2893
2894	break;
2895
2896	case 3:
2897
2898	pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 - dhsize*fsegment_fract;
2899	pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 - dtsize*fsegment_fract;
2900
2901	break;
2902
2903	case 4:
2904
2905	pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + dhsize*fsegment_fract;
2906	pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 - dtsize*fsegment_fract;
2907
2908	break;
2909
2910	case 5:
2911
2912	pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + ddxseg4 + dpsize*fsegment_fract;
2913	pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 + ddyseg4 + 0.;
2914
2915	break;
2916
2917	case 6:
2918
2919	pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + ddxseg4 + ddxseg5 + dhsize*fsegment_fract;
2920	pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 + ddyseg4 + ddyseg5 + dtsize*fsegment_fract;
2921
2922	break;
2923
2924	default:
2925
2926	fprintf(stderr, "ERROR: problem in coordinate generation for pixel %d\n", ipixno);
2927	return(0);
2928
2929	} /* end switch */
2930
2931	} /* end for */
2932
2933	dstartx = pcam->dxc[pcam->inumcentralpixels - 1] + dhsize;
2934	dstarty = pcam->dyc[pcam->inumcentralpixels - 1] + dtsize;
2935
2936	if(pcam->inumgappixels > 0){ /* generate the positions of the gap pixels */
2937
2938	j = pcam->inumcentralpixels;
2939
2940	for(i=0; i<pcam->inumgappixels; i=i+6){
2941	pcam->dxc[j + i ] = dstartx + 2. * (i/6 + 1) * dpsize;
2942	pcam->dyc[j + i ] = dstarty;
2943	pcam->dpixsizefactor[j + i] = 1.;
2944	pcam->dxc[j + i + 1] = pcam->dxc[j + i ] / 2.;
2945	pcam->dyc[j + i + 1] = sqrt(3.) * pcam->dxc[j + i + 1];
2946	pcam->dpixsizefactor[j + i + 1] = 1.;
2947	pcam->dxc[j + i + 2] = - pcam->dxc[j + i + 1];
2948	pcam->dyc[j + i + 2] = pcam->dyc[j + i + 1];
2949	pcam->dpixsizefactor[j + i+ 2] = 1.;
2950	pcam->dxc[j + i + 3] = - pcam->dxc[j + i];
2951	pcam->dyc[j + i + 3] = dstarty;
2952	pcam->dpixsizefactor[j + i+ 3] = 1.;
2953	pcam->dxc[j + i + 4] = pcam->dxc[j + i + 2];
2954	pcam->dyc[j + i + 4] = - pcam->dyc[j + i + 2];
2955	pcam->dpixsizefactor[j + i+ 4] = 1.;
2956	pcam->dxc[j + i + 5] = pcam->dxc[j + i + 1];
2957	pcam->dyc[j + i + 5] = - pcam->dyc[j + i + 1];
2958	pcam->dpixsizefactor[j + i + 5] = 1.;
2959	} /* end for */
2960	} /* end if */
2961
2962	/* generate positions of the outer pixels */
2963
2964	if( pcam->inumbigpixels > 0 ){
2965
2966	j = pcam->inumcentralpixels + pcam->inumgappixels;
2967
2968	for(i=0; i<pcam->inumbigpixels; i++){
2969	pcam->dpixsizefactor[j + i] = 2.;
2970	}
2971
2972	in = 0;
2973
2974	nrow = (int) ceil(dstartx / 2. / dpsize);
2975
2976	while(in < pcam->inumbigpixels){
2977
2978	pcam->dxc[j + in] = dstartx + dpsize;
2979	pcam->dyc[j + in] = dstarty + 2 * dpsize / sqrt(3.);
2980	pcam->dxc[j + in + nrow] = dstartx / 2. - dpsize / 2.;
2981	pcam->dyc[j + in + nrow] = sqrt(3.)/2. * dstartx + 2.5 * dpsize/sqrt(3.);
2982	pcam->dxc[j + in + 3 * nrow - 1] = - pcam->dxc[j + in];
2983	pcam->dyc[j + in + 3 * nrow - 1] = pcam->dyc[j + in];
2984	pcam->dxc[j + in + 3 * nrow] = - pcam->dxc[j + in];
2985	pcam->dyc[j + in + 3 * nrow] = - pcam->dyc[j + in];
2986	pcam->dxc[j + in + 5 * nrow - 1] = pcam->dxc[j + in + nrow];
2987	pcam->dyc[j + in + 5 * nrow - 1] = - pcam->dyc[j + in + nrow];
2988	pcam->dxc[j + in + 6 * nrow - 1] = pcam->dxc[j + in];
2989	pcam->dyc[j + in + 6 * nrow - 1] = - pcam->dyc[j + in];
2990	for(i=1; i<nrow; i++){
2991	pcam->dxc[j + in + i] = pcam->dxc[j + in] - i * dpsize;
2992	pcam->dyc[j + in + i] = pcam->dyc[j + in] + i * dpsize * sqrt(3.);
2993	pcam->dxc[j + in + i + nrow] = pcam->dxc[j + in + nrow] - i * 2 * dpsize;
2994	pcam->dyc[j + in + i + nrow] = pcam->dyc[j + in + nrow];
2995	pcam->dxc[j + in + 3 * nrow - 1 - i] = - pcam->dxc[j + in + i];
2996	pcam->dyc[j + in + 3 * nrow - 1- i] = pcam->dyc[j + in + i];
2997	pcam->dxc[j + in + i + 3 * nrow] = - pcam->dxc[j + in + i];
2998	pcam->dyc[j + in + i + 3 * nrow] = - pcam->dyc[j + in + i];
2999	pcam->dxc[j + in + 5 * nrow - 1 - i] = pcam->dxc[j + in + i + nrow];
3000	pcam->dyc[j + in + 5 * nrow - 1 - i] = - pcam->dyc[j + in + i + nrow];
3001	pcam->dxc[j + in + 6 * nrow - 1 - i] = pcam->dxc[j + in + i];
3002	pcam->dyc[j + in + 6 * nrow - 1 - i] = - pcam->dyc[j + in + i];
3003	}
3004	in = in + 6 * nrow;
3005	dstartx = dstartx + 2. * dpsize;
3006	nrow = nrow + 1;
3007	} /* end while */
3008
3009	} /* end if */
3010
3011	/* generate the ij coordinates */
3012
3013	for(i=0; i<pcam->inumpixels; i++){
3014	pcam->dj[i] = pcam->dyc[i]/SIN60/dpsize;
3015	pcam->di[i] = pcam->dxc[i]/dpsize - pcam->dj[i]*COS60;
3016
3017	// fprintf(stdout, "%d %f %f %f %f %f\n",
3018	// i+1, pcam->di[i], pcam->dj[i], pcam->dxc[i], pcam->dyc[i],
3019	// pcam->dpixsizefactor[i]);
3020
3021	}
3022
3023	return(pcam->inumpixels);
3024
3025	}
3026	//!@}
3027
3028	//!-----------------------------------------------------------
3029	// @name bpoint_is_in_pix
3030	//
3031	// @desc check if a point (x,y) in camera coordinates is inside a given pixel
3032	//
3033	// @var *pcam structure camera containing all the
3034	// camera information
3035	// @var dx, dy point coordinates in centimeters
3036	// @var ipixnum pixel number (starting at 0)
3037	// @return TRUE if the point is inside the pixel, FALSE otherwise
3038	//
3039	// DP
3040	//
3041	// @date Thu Oct 14 16:59:04 CEST 1999
3042	//------------------------------------------------------------
3043	// @function
3044
3045	//!@{
3046
3047	/****** bpoint_is_in_pix() *************************************/
3048
3049	int bpoint_is_in_pix(double dx, double dy, int ipixnum, struct camera *pcam){
3050	/* return TRUE if point (dx, dy) is in pixel number ipixnum, else return FALSE (use camera coordinate system) */
3051	/* the pixel is assumed to be a "closed set" */
3052
3053	double a, b; /* a = length of one of the edges of one pixel, b = half the width of one pixel */
3054	double c, xx, yy; /* auxiliary variable */
3055
3056	b = pcam->dpixdiameter_cm / 2. * pcam->dpixsizefactor[ipixnum];
3057	a = pcam->dpixdiameter_cm / sqrt(3.) * pcam->dpixsizefactor[ipixnum];
3058	c = 1. - 1./sqrt(3.);
3059	if((ipixnum < 0)\|\|(ipixnum >= pcam->inumpixels)){
3060	fprintf(stderr, "Error in bpoint_is_in_pix: invalid pixel number %d\n", ipixnum);
3061	fprintf(stderr, "Exiting.\n");
3062	exit(203);
3063	}
3064	xx = dx - pcam->dxc[ipixnum];
3065	yy = dy - pcam->dyc[ipixnum];
3066
3067	if(((-b <= xx) && (xx <= 0.) && ((-c * xx - a) <= yy) && (yy <= ( c * xx + a))) \|\|
3068	((0. < xx) && (xx <= b ) && (( c * xx - a) <= yy) && (yy <= (-c * xx + a))) ){
3069	return(TRUE); /* point is inside */
3070	}
3071	else{
3072	return(FALSE); /* point is outside */
3073	}
3074	}
3075
3076	//!@}
3077
3078	//------------------------------------------------------------
3079	// @name dist_r_P
3080	//
3081	// @desc distance straight line r - point P
3082	//
3083	// @date Sat Jun 27 05:58:56 MET DST 1998
3084	// @function @code
3085	//------------------------------------------------------------
3086	// dist_r_P
3087	//
3088	// distance straight line r - point P
3089	//------------------------------------------------------------
3090
3091	float
3092	dist_r_P(float a, float b, float c,
3093	float l, float m, float n,
3094	float x, float y, float z)
3095	{
3096	return (
3097	sqrt((SQR((a-x)m-(b-y)l) +
3098	SQR((b-y)n-(c-z)m) +
3099	SQR((c-z)l-(a-x)n))/
3100	(SQR(l)+SQR(m)+SQR(n))
3101	)
3102	);
3103	}
3104	// @endcode
3105
3106
3107	//=------------------------------------------------------------
3108	//!@subsection Log of this file.
3109
3110	//!@{
3111	//
3112	// $Log: not supported by cvs2svn $
3113	// Revision 1.3 1999/10/22 15:01:28 petry
3114	// version sent to H.K. and N.M. on Fri Oct 22 1999
3115	//
3116	// Revision 1.2 1999/10/22 09:44:23 petry
3117	// first synthesized version which compiles and runs without crashing;
3118	//
3119	// Revision 1.1.1.1 1999/10/21 16:35:10 petry
3120	// first synthesised version
3121	//
3122	// Revision 1.13 1999/03/15 14:59:05 gonzalez
3123	// camera-1_1
3124	//
3125	// Revision 1.12 1999/03/02 09:56:10 gonzalez
3126	// * empty log message *
3127	//
3128	//
3129	//!@}
3130
3131	//=EOF

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: branches/start/MagicSoft/Simulation/Detector/Camera/camera.cxx@ 10113

Download in other formats: